File indexing completed on 2025-04-04 08:05: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 "EventAction.hh"
0037 #include "HistoManager.hh"
0038 #include "PrimaryGeneratorAction.hh"
0039
0040 #include "G4Event.hh"
0041 #include "G4Material.hh"
0042 #include "G4SystemOfUnits.hh"
0043 #include "G4UnitsTable.hh"
0044
0045 #include <iomanip>
0046
0047
0048
0049 Run::Run(DetectorConstruction* detector) : fDetector(detector)
0050 {
0051 fTotEdep[1] = fEleak[1] = fEtotal[1] = joule;
0052
0053 for (G4int i = 0; i < kMaxAbsor; ++i) {
0054 fEdeposit[i] = 0.;
0055 fEmin[i] = joule;
0056 fEmax[i] = 0.;
0057 }
0058 }
0059
0060
0061
0062 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0063 {
0064 fParticle = particle;
0065 fEkin = energy;
0066 }
0067
0068
0069
0070 void Run::CountProcesses(const G4VProcess* process)
0071 {
0072 if (process == nullptr) return;
0073 G4String procName = process->GetProcessName();
0074 std::map<G4String, G4int>::iterator it = fProcCounter.find(procName);
0075 if (it == fProcCounter.end()) {
0076 fProcCounter[procName] = 1;
0077 }
0078 else {
0079 fProcCounter[procName]++;
0080 }
0081 }
0082
0083
0084
0085 void Run::ParticleCount(G4int k, G4String name, G4double Ekin, G4double meanLife)
0086 {
0087 std::map<G4String, ParticleData>::iterator it = fParticleDataMap[k].find(name);
0088 if (it == fParticleDataMap[k].end()) {
0089 (fParticleDataMap[k])[name] = ParticleData(1, Ekin, Ekin, Ekin, meanLife);
0090 }
0091 else {
0092 ParticleData& data = it->second;
0093 data.fCount++;
0094 data.fEmean += Ekin;
0095
0096 G4double emin = data.fEmin;
0097 if (Ekin < emin) data.fEmin = Ekin;
0098 G4double emax = data.fEmax;
0099 if (Ekin > emax) data.fEmax = Ekin;
0100 data.fTmean = meanLife;
0101 }
0102 }
0103
0104
0105
0106 void Run::AddEdep(G4int i, G4double e)
0107 {
0108 if (e > 0.) {
0109 fEdeposit[i] += e;
0110 if (e < fEmin[i]) fEmin[i] = e;
0111 if (e > fEmax[i]) fEmax[i] = e;
0112 }
0113 }
0114
0115
0116
0117 void Run::AddTotEdep(G4double e)
0118 {
0119 if (e > 0.) {
0120 fTotEdep[0] += e;
0121 if (e < fTotEdep[1]) fTotEdep[1] = e;
0122 if (e > fTotEdep[2]) fTotEdep[2] = e;
0123 }
0124 }
0125
0126
0127
0128 void Run::AddEleak(G4double e)
0129 {
0130 if (e > 0.) {
0131 fEleak[0] += e;
0132 if (e < fEleak[1]) fEleak[1] = e;
0133 if (e > fEleak[2]) fEleak[2] = e;
0134 }
0135 }
0136
0137
0138
0139 void Run::AddEtotal(G4double e)
0140 {
0141 if (e > 0.) {
0142 fEtotal[0] += e;
0143 if (e < fEtotal[1]) fEtotal[1] = e;
0144 if (e > fEtotal[2]) fEtotal[2] = e;
0145 }
0146 }
0147
0148
0149
0150 void Run::AddTrackStatus(G4int i)
0151 {
0152 fStatus[i]++;
0153 }
0154
0155
0156
0157 void Run::Merge(const G4Run* run)
0158 {
0159 const Run* localRun = static_cast<const Run*>(run);
0160
0161
0162 fParticle = localRun->fParticle;
0163 fEkin = localRun->fEkin;
0164
0165
0166
0167 G4int nbOfAbsor = fDetector->GetNbOfAbsor();
0168 for (G4int i = 1; i <= nbOfAbsor; ++i) {
0169 fEdeposit[i] += localRun->fEdeposit[i];
0170
0171 G4double min, max;
0172 min = localRun->fEmin[i];
0173 max = localRun->fEmax[i];
0174 if (fEmin[i] > min) fEmin[i] = min;
0175 if (fEmax[i] < max) fEmax[i] = max;
0176 }
0177
0178 for (G4int i = 0; i < 3; ++i)
0179 fStatus[i] += localRun->fStatus[i];
0180
0181
0182 fTotEdep[0] += localRun->fTotEdep[0];
0183 G4double min, max;
0184 min = localRun->fTotEdep[1];
0185 max = localRun->fTotEdep[2];
0186 if (fTotEdep[1] > min) fTotEdep[1] = min;
0187 if (fTotEdep[2] < max) fTotEdep[2] = max;
0188
0189
0190 fEleak[0] += localRun->fEleak[0];
0191 min = localRun->fEleak[1];
0192 max = localRun->fEleak[2];
0193 if (fEleak[1] > min) fEleak[1] = min;
0194 if (fEleak[2] < max) fEleak[2] = max;
0195
0196
0197 fEtotal[0] += localRun->fEtotal[0];
0198 min = localRun->fEtotal[1];
0199 max = localRun->fEtotal[2];
0200 if (fEtotal[1] > min) fEtotal[1] = min;
0201 if (fEtotal[2] < max) fEtotal[2] = max;
0202
0203
0204 std::map<G4String, G4int>::const_iterator itp;
0205 for (itp = localRun->fProcCounter.begin(); itp != localRun->fProcCounter.end(); ++itp) {
0206 G4String procName = itp->first;
0207 G4int localCount = itp->second;
0208 if (fProcCounter.find(procName) == fProcCounter.end()) {
0209 fProcCounter[procName] = localCount;
0210 }
0211 else {
0212 fProcCounter[procName] += localCount;
0213 }
0214 }
0215
0216
0217 for (G4int k = 0; k <= nbOfAbsor; ++k) {
0218 std::map<G4String, ParticleData>::const_iterator itc;
0219 for (itc = localRun->fParticleDataMap[k].begin(); itc != localRun->fParticleDataMap[k].end();
0220 ++itc)
0221 {
0222 G4String name = itc->first;
0223 const ParticleData& localData = itc->second;
0224 if (fParticleDataMap[k].find(name) == fParticleDataMap[k].end()) {
0225 (fParticleDataMap[k])[name] = ParticleData(
0226 localData.fCount, localData.fEmean, localData.fEmin, localData.fEmax, localData.fTmean);
0227 }
0228 else {
0229 ParticleData& data = (fParticleDataMap[k])[name];
0230 data.fCount += localData.fCount;
0231 data.fEmean += localData.fEmean;
0232 G4double emin = localData.fEmin;
0233 if (emin < data.fEmin) data.fEmin = emin;
0234 G4double emax = localData.fEmax;
0235 if (emax > data.fEmax) data.fEmax = emax;
0236 data.fTmean = localData.fTmean;
0237 }
0238 }
0239 }
0240
0241 G4Run::Merge(run);
0242 }
0243
0244
0245
0246 void Run::EndOfRun()
0247 {
0248 G4int prec = 5, wid = prec + 2;
0249 G4int dfprec = G4cout.precision(prec);
0250
0251
0252
0253 G4String partName = fParticle->GetParticleName();
0254 G4int nbOfAbsor = fDetector->GetNbOfAbsor();
0255
0256 G4cout << "\n ======================== run summary =====================\n";
0257 G4cout << "\n The run is " << numberOfEvent << " " << partName << " of "
0258 << G4BestUnit(fEkin, "Energy") << " through " << nbOfAbsor << " absorbers: \n";
0259 for (G4int i = 1; i <= nbOfAbsor; i++) {
0260 G4Material* material = fDetector->GetAbsorMaterial(i);
0261 G4double thickness = fDetector->GetAbsorThickness(i);
0262 G4double density = material->GetDensity();
0263 G4cout << std::setw(5) << i << std::setw(10) << G4BestUnit(thickness, "Length") << " of "
0264 << material->GetName() << " (density: " << G4BestUnit(density, "Volumic Mass") << ")"
0265 << G4endl;
0266 }
0267
0268 if (numberOfEvent == 0) {
0269 G4cout.precision(dfprec);
0270 return;
0271 }
0272
0273 G4cout.precision(3);
0274
0275
0276
0277 G4cout << "\n Process calls frequency :" << G4endl;
0278 G4int index = 0;
0279 std::map<G4String, G4int>::iterator it;
0280 for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
0281 G4String procName = it->first;
0282 G4int count = it->second;
0283 G4String space = " ";
0284 if (++index % 3 == 0) space = "\n";
0285 G4cout << " " << std::setw(20) << procName << "=" << std::setw(7) << count << space;
0286 }
0287 G4cout << G4endl;
0288
0289
0290
0291 for (G4int i = 1; i <= nbOfAbsor; i++) {
0292 fEdeposit[i] /= numberOfEvent;
0293
0294 G4cout << "\n Edep in absorber " << i << " = " << G4BestUnit(fEdeposit[i], "Energy") << "\t("
0295 << G4BestUnit(fEmin[i], "Energy") << "-->" << G4BestUnit(fEmax[i], "Energy") << ")";
0296 }
0297 G4cout << G4endl;
0298
0299 if (nbOfAbsor > 1) {
0300 fTotEdep[0] /= numberOfEvent;
0301 G4cout << "\n Edep in all absorb = " << G4BestUnit(fTotEdep[0], "Energy") << "\t("
0302 << G4BestUnit(fTotEdep[1], "Energy") << "-->" << G4BestUnit(fTotEdep[2], "Energy") << ")"
0303 << G4endl;
0304 }
0305
0306
0307
0308 fEleak[0] /= numberOfEvent;
0309 G4cout << " Energy leakage = " << G4BestUnit(fEleak[0], "Energy") << "\t("
0310 << G4BestUnit(fEleak[1], "Energy") << "-->" << G4BestUnit(fEleak[2], "Energy") << ")"
0311 << G4endl;
0312
0313
0314
0315 fEtotal[0] /= numberOfEvent;
0316 G4cout << " Energy total = " << G4BestUnit(fEtotal[0], "Energy") << "\t("
0317 << G4BestUnit(fEtotal[1], "Energy") << "-->" << G4BestUnit(fEtotal[2], "Energy") << ")"
0318 << G4endl;
0319
0320
0321
0322 for (G4int k = 1; k <= nbOfAbsor; k++) {
0323 G4cout << "\n List of created particles in absorber " << k << ":" << G4endl;
0324
0325 std::map<G4String, ParticleData>::iterator itc;
0326 for (itc = fParticleDataMap[k].begin(); itc != fParticleDataMap[k].end(); itc++) {
0327 G4String name = itc->first;
0328 ParticleData data = itc->second;
0329 G4int count = data.fCount;
0330 G4double eMean = data.fEmean / count;
0331 G4double eMin = data.fEmin;
0332 G4double eMax = data.fEmax;
0333 G4double meanLife = data.fTmean;
0334
0335 G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
0336 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
0337 << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")";
0338 if (meanLife >= 0.)
0339 G4cout << "\tmean life = " << G4BestUnit(meanLife, "Time") << G4endl;
0340 else
0341 G4cout << "\tstable" << G4endl;
0342 }
0343 }
0344
0345
0346 G4cout << "\n List of particles emerging from absorbers :" << G4endl;
0347
0348 std::map<G4String, ParticleData>::iterator itc;
0349 for (itc = fParticleDataMap[0].begin(); itc != fParticleDataMap[0].end(); itc++) {
0350 G4String name = itc->first;
0351 ParticleData data = itc->second;
0352 G4int count = data.fCount;
0353 G4double eMean = data.fEmean / count;
0354 G4double eMin = data.fEmin;
0355 G4double eMax = data.fEmax;
0356
0357
0358 G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
0359 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
0360 << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")" << G4endl;
0361 }
0362
0363
0364
0365 G4double dNofEvents = double(numberOfEvent);
0366 G4double absorbed = 100. * fStatus[0] / dNofEvents;
0367 G4double transmit = 100. * fStatus[1] / dNofEvents;
0368 G4double reflected = 100. * fStatus[2] / dNofEvents;
0369
0370 G4cout.precision(2);
0371 G4cout << "\n Nb of events with primary absorbed = " << absorbed << " %,"
0372 << " transmit = " << transmit << " %,"
0373 << " reflected = " << reflected << " %" << G4endl;
0374
0375
0376
0377 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0378 G4int ih = 10;
0379 G4double binWidth = analysisManager->GetH1Width(ih) * analysisManager->GetH1Unit(ih);
0380 G4double fac = (1. / (numberOfEvent * binWidth)) * (mm / MeV);
0381 analysisManager->ScaleH1(ih, fac);
0382
0383
0384 fProcCounter.clear();
0385 for (G4int k = 0; k <= nbOfAbsor; k++)
0386 fParticleDataMap[k].clear();
0387
0388
0389 G4cout.precision(dfprec);
0390 }
0391
0392
