File indexing completed on 2025-02-23 09:22:18
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033 #include "Run.hh"
0034
0035 #include "DetectorConstruction.hh"
0036 #include "HistoManager.hh"
0037 #include "PrimaryGeneratorAction.hh"
0038
0039 #include "G4EmCalculator.hh"
0040 #include "G4SystemOfUnits.hh"
0041 #include "G4Track.hh"
0042 #include "G4UnitsTable.hh"
0043 #include "G4VPhysicalVolume.hh"
0044
0045 #include <iomanip>
0046
0047
0048
0049 Run::Run(DetectorConstruction* det, HistoManager* histoMgr)
0050 : G4Run(), fDetector(det), fHistoMgr(histoMgr)
0051 {
0052 fSumR = 0;
0053 fNevt = fNelec = fNposit = fNgam = fNstep = fNgamPh = fNgamTar = fNeTar = fNePh = fNstepTarget =
0054 0;
0055
0056 fAnalysisManager = G4AnalysisManager::Instance();
0057
0058 fHistoId = fHistoMgr->GetHistoIdentifiers();
0059 fNHisto = fHistoId.size();
0060
0061 fCheckVolume = fDetector->GetCheckVolume();
0062 fGasVolume = fDetector->GetGasVolume();
0063 fPhantom = fDetector->GetPhantom();
0064 fTarget1 = fDetector->GetTarget1();
0065 fTarget2 = fDetector->GetTarget2();
0066
0067 fNBinsR = fDetector->GetNumberDivR();
0068 fNBinsZ = fDetector->GetNumberDivZ();
0069
0070 fScoreZ = fDetector->GetScoreZ();
0071 fAbsorberZ = fDetector->GetAbsorberZ();
0072 fAbsorberR = fDetector->GetAbsorberR();
0073 fMaxEnergy = fDetector->GetMaxEnergy();
0074
0075 fNBinsE = fDetector->GetNumberDivE();
0076 fMaxEnergy = fDetector->GetMaxEnergy();
0077
0078 fStepZ = fAbsorberZ / (G4double)fNBinsZ;
0079 fStepR = fAbsorberR / (G4double)fNBinsR;
0080 fStepE = fMaxEnergy / (G4double)fNBinsE;
0081 fScoreBin = (G4int)(fScoreZ / fStepZ + 0.5);
0082
0083 fVerbose = fDetector->GetVerbose();
0084
0085 fGamma = G4Gamma::Gamma();
0086 fElectron = G4Electron::Electron();
0087 fPositron = G4Positron::Positron();
0088
0089 G4cout << " " << fNBinsR << " bins R stepR= " << fStepR / mm << " mm " << G4endl;
0090 G4cout << " " << fNBinsZ << " bins Z stepZ= " << fStepZ / mm << " mm " << G4endl;
0091 G4cout << " " << fNBinsE << " bins E stepE= " << fStepE / MeV << " MeV " << G4endl;
0092 G4cout << " " << fScoreBin << "-th bin in Z is used for R distribution" << G4endl;
0093
0094 fVolumeR.clear();
0095 fEdep.clear();
0096 fGammaE.clear();
0097
0098 fVolumeR.resize(fNBinsR, 0.0);
0099 fEdep.resize(fNBinsR, 0.0);
0100 fGammaE.resize(fNBinsE, 0.0);
0101
0102 G4double r1 = 0.0;
0103 G4double r2 = fStepR;
0104 for (G4int i = 0; i < fNBinsR; ++i) {
0105 fVolumeR[i] = cm * cm / (CLHEP::pi * (r2 * r2 - r1 * r1));
0106 r1 = r2;
0107 r2 += fStepR;
0108 }
0109
0110 if (fAnalysisManager->GetFileName() != "") fHistoMgr->Update(det, true);
0111 }
0112
0113
0114
0115 Run::~Run() {}
0116
0117
0118
0119 void Run::Merge(const G4Run* run)
0120 {
0121 const Run* localRun = static_cast<const Run*>(run);
0122
0123 fNevt += localRun->fNevt;
0124
0125 fNelec += localRun->fNelec;
0126 fNgam += localRun->fNgam;
0127 fNposit += localRun->fNposit;
0128
0129 fNgamTar += localRun->fNgamTar;
0130 fNgamPh += localRun->fNgamPh;
0131 fNeTar += localRun->fNeTar;
0132 fNePh += localRun->fNePh;
0133
0134 fNstep += localRun->fNstep;
0135 fNstepTarget += localRun->fNstepTarget;
0136
0137 fSumR += localRun->fSumR;
0138
0139 for (int i = 0; i < (int)fEdep.size(); i++)
0140 fEdep[i] += localRun->fEdep[i];
0141 for (int i = 0; i < (int)fGammaE.size(); i++)
0142 fGammaE[i] += localRun->fGammaE[i];
0143
0144 G4Run::Merge(run);
0145 }
0146
0147
0148
0149 void Run::EndOfRun()
0150 {
0151 G4cout << "Histo: End of run actions are started" << G4endl;
0152
0153
0154 G4cout << "========================================================" << G4endl;
0155 G4double x = (G4double)fNevt;
0156 if (fNevt > 0) {
0157 x = 1.0 / x;
0158 }
0159 G4double xe = x * (G4double)fNelec;
0160 G4double xg = x * (G4double)fNgam;
0161 G4double xp = x * (G4double)fNposit;
0162 G4double xs = x * (G4double)fNstep;
0163 G4double xph = x * (G4double)fNgamPh;
0164 G4double xes = x * (G4double)fNstepTarget;
0165 G4double xgt = x * (G4double)fNgamTar;
0166 G4double xet = x * (G4double)fNeTar;
0167 G4double xphe = x * (G4double)fNePh;
0168
0169 G4cout << "Number of events " << std::setprecision(8) << fNevt
0170 << G4endl;
0171 G4cout << std::setprecision(4) << "Average number of e- " << xe << G4endl;
0172 G4cout << std::setprecision(4) << "Average number of gamma " << xg << G4endl;
0173 G4cout << std::setprecision(4) << "Average number of e+ " << xp << G4endl;
0174 G4cout << std::setprecision(4) << "Average number of steps in the phantom " << xs << G4endl;
0175 G4cout << std::setprecision(4) << "Average number of e- steps in the target " << xes
0176 << G4endl;
0177 G4cout << std::setprecision(4) << "Average number of g produced in the target " << xgt
0178 << G4endl;
0179 G4cout << std::setprecision(4) << "Average number of e- produced in the target " << xet
0180 << G4endl;
0181 G4cout << std::setprecision(4) << "Average number of g produced in the phantom " << xph
0182 << G4endl;
0183 G4cout << std::setprecision(4) << "Average number of e- produced in the phantom " << xphe
0184 << G4endl;
0185 G4cout << std::setprecision(4) << "Total fGamma fluence in front of the phantom "
0186 << x * fSumR / MeV << " MeV " << G4endl;
0187 G4cout << "========================================================" << G4endl;
0188 G4cout << G4endl;
0189 G4cout << G4endl;
0190
0191 G4double sum = 0.0;
0192 for (G4int i = 0; i < fNBinsR; i++) {
0193 fEdep[i] *= x;
0194 sum += fEdep[i];
0195 }
0196
0197 if (fAnalysisManager) {
0198 if (fAnalysisManager->IsActive()) {
0199
0200 for (G4int i = 0; i < fNHisto; i++) {
0201 fAnalysisManager->GetH1(fHistoId[i])->scale(x);
0202 }
0203 G4double nr = fEdep[0];
0204 if (nr > 0.0) {
0205 nr = 1. / nr;
0206 }
0207 fAnalysisManager->GetH1(fHistoId[0])->scale(nr);
0208
0209 nr = sum * fStepR;
0210 if (nr > 0.0) {
0211 nr = 1. / nr;
0212 }
0213 fAnalysisManager->GetH1(fHistoId[1])->scale(nr);
0214
0215 fAnalysisManager->GetH1(fHistoId[3])
0216 ->scale(1000.0 * cm3 / (CLHEP::pi * fAbsorberR * fAbsorberR * fStepZ));
0217 fAnalysisManager->GetH1(fHistoId[4])->scale(1000.0 * cm3 * fVolumeR[0] / fStepZ);
0218
0219
0220 if (!fAnalysisManager->Write()) {
0221 G4Exception("Histo::Save()", "hist01", FatalException, "Cannot write ROOT file.");
0222 }
0223 G4cout << "### Histo::Save: Histograms are saved" << G4endl;
0224 if (fAnalysisManager->CloseFile() && fVerbose) {
0225 G4cout << " File is closed" << G4endl;
0226 }
0227 }
0228 }
0229 }
0230
0231
0232
0233 void Run::AddTargetPhoton(const G4DynamicParticle* p)
0234 {
0235 ++fNgamTar;
0236 if (fAnalysisManager) {
0237 fAnalysisManager->FillH1(fHistoId[5], p->GetKineticEnergy() / MeV, 1.0);
0238 }
0239 }
0240
0241
0242
0243 void Run::AddPhantomPhoton(const G4DynamicParticle*)
0244 {
0245 ++fNgamPh;
0246 }
0247
0248
0249
0250 void Run::AddTargetElectron(const G4DynamicParticle* p)
0251 {
0252 ++fNeTar;
0253 if (fAnalysisManager) {
0254 fAnalysisManager->FillH1(fHistoId[8], p->GetKineticEnergy() / MeV, 1.0);
0255 }
0256 }
0257
0258
0259
0260 void Run::AddPhantomElectron(const G4DynamicParticle* p)
0261 {
0262 ++fNePh;
0263 if (fAnalysisManager) {
0264 fAnalysisManager->FillH1(fHistoId[7], p->GetKineticEnergy() / MeV, 1.0);
0265 }
0266 }
0267
0268
0269
0270 void Run::ScoreNewTrack(const G4Track* aTrack)
0271 {
0272
0273 const G4ParticleDefinition* particle = aTrack->GetParticleDefinition();
0274 G4int pid = aTrack->GetParentID();
0275 G4VPhysicalVolume* pv = aTrack->GetVolume();
0276 const G4DynamicParticle* dp = aTrack->GetDynamicParticle();
0277
0278
0279 if (0 == pid) {
0280 ++fNevt;
0281 if (fVerbose) {
0282 G4ThreeVector pos = aTrack->GetVertexPosition();
0283 G4ThreeVector dir = aTrack->GetMomentumDirection();
0284 G4cout << "TrackingAction: Primary " << particle->GetParticleName()
0285 << " Ekin(MeV)= " << aTrack->GetKineticEnergy() / MeV << "; pos= " << pos
0286 << "; dir= " << dir << G4endl;
0287 }
0288
0289 }
0290 else if (0 < pid && particle == fElectron) {
0291 if (fVerbose) {
0292 G4cout << "TrackingAction: Secondary electron " << G4endl;
0293 }
0294 AddElectron();
0295 if (pv == fPhantom) {
0296 AddPhantomElectron(dp);
0297 }
0298 else if (pv == fTarget1 || pv == fTarget2) {
0299 AddTargetElectron(dp);
0300 }
0301
0302
0303 }
0304 else if (0 < pid && particle == fPositron) {
0305 if (fVerbose) {
0306 G4cout << "TrackingAction: Secondary positron " << G4endl;
0307 }
0308 AddPositron();
0309
0310
0311 }
0312 else if (0 < pid && particle == fGamma) {
0313 if (fVerbose) {
0314 G4cout << "TrackingAction: Secondary gamma; parentID= " << pid
0315 << " E= " << aTrack->GetKineticEnergy() << G4endl;
0316 }
0317 AddPhoton();
0318 if (pv == fPhantom) {
0319 AddPhantomPhoton(dp);
0320 }
0321 else if (pv == fTarget1 || pv == fTarget2) {
0322 AddTargetPhoton(dp);
0323 }
0324 }
0325 }
0326
0327
0328
0329 void Run::AddPhantomGamma(G4double e, G4double r)
0330 {
0331 e /= MeV;
0332 fSumR += e;
0333 G4int bin = (G4int)(e / fStepE);
0334 if (bin >= fNBinsE) {
0335 bin = fNBinsE - 1;
0336 }
0337 fGammaE[bin] += e;
0338 G4int bin1 = (G4int)(r / fStepR);
0339 if (bin1 >= fNBinsR) {
0340 bin1 = fNBinsR - 1;
0341 }
0342 if (fAnalysisManager) {
0343 G4AnalysisManager::Instance()->FillH1(fHistoId[6], e, 1.0);
0344 G4AnalysisManager::Instance()->FillH1(fHistoId[9], r, e * fVolumeR[bin1]);
0345 }
0346 }
0347
0348
0349
0350 void Run::AddPhantomStep(G4double edep, G4double r1, G4double z1, G4double r2, G4double z2,
0351 G4double r0, G4double z0)
0352 {
0353 ++fNstep;
0354 G4int nzbin = (G4int)(z0 / fStepZ);
0355 if (fVerbose) {
0356 G4cout << "Histo: edep(MeV)= " << edep / MeV << " at binz= " << nzbin << " r1= " << r1
0357 << " z1= " << z1 << " r2= " << r2 << " z2= " << z2 << " r0= " << r0 << " z0= " << z0
0358 << G4endl;
0359 }
0360 if (nzbin == fScoreBin) {
0361 G4int bin = (G4int)(r0 / fStepR);
0362 if (bin >= fNBinsR) {
0363 bin = fNBinsR - 1;
0364 }
0365 double w = edep * fVolumeR[bin];
0366 fEdep[bin] += w;
0367
0368 if (fAnalysisManager) {
0369 G4AnalysisManager::Instance()->FillH1(fHistoId[0], r0, w);
0370 G4AnalysisManager::Instance()->FillH1(fHistoId[1], r0, w);
0371 G4AnalysisManager::Instance()->FillH1(fHistoId[2], r0, w);
0372 }
0373 }
0374 G4int bin1 = (G4int)(z1 / fStepZ);
0375 if (bin1 >= fNBinsZ) {
0376 bin1 = fNBinsZ - 1;
0377 }
0378 G4int bin2 = (G4int)(z2 / fStepZ);
0379 if (bin2 >= fNBinsZ) {
0380 bin2 = fNBinsZ - 1;
0381 }
0382 if (bin1 == bin2) {
0383 if (fAnalysisManager) {
0384 G4AnalysisManager::Instance()->FillH1(fHistoId[3], z0, edep);
0385 if (r1 < fStepR) {
0386 G4double w = edep;
0387 if (r2 > fStepR) {
0388 w *= (fStepR - r1) / (r2 - r1);
0389 }
0390 G4AnalysisManager::Instance()->FillH1(fHistoId[4], z0, w);
0391 }
0392 }
0393 }
0394 else {
0395 G4int bin;
0396
0397 if (bin2 < bin1) {
0398 bin = bin2;
0399 G4double z = z2;
0400 bin2 = bin1;
0401 z2 = z1;
0402 bin1 = bin;
0403 z1 = z;
0404 }
0405 G4double zz1 = z1;
0406 G4double zz2 = (bin1 + 1) * fStepZ;
0407 G4double rr1 = r1;
0408 G4double dz = z2 - z1;
0409 G4double dr = r2 - r1;
0410 G4double rr2 = r1 + dr * (zz2 - zz1) / dz;
0411 for (bin = bin1; bin <= bin2; bin++) {
0412 if (fAnalysisManager) {
0413 G4double de = edep * (zz2 - zz1) / dz;
0414 G4double zf = (zz1 + zz2) * 0.5;
0415 {
0416 G4AnalysisManager::Instance()->FillH1(fHistoId[3], zf, de);
0417 }
0418 if (rr1 < fStepR) {
0419 G4double w = de;
0420 if (rr2 > fStepR) w *= (fStepR - rr1) / (rr2 - rr1);
0421 {
0422 G4AnalysisManager::Instance()->FillH1(fHistoId[4], zf, w);
0423 }
0424 }
0425 }
0426 zz1 = zz2;
0427 zz2 = std::min(z2, zz1 + fStepZ);
0428 rr1 = rr2;
0429 rr2 = rr1 + dr * (zz2 - zz1) / dz;
0430 }
0431 }
0432 }
0433
0434
