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