File indexing completed on 2026-06-03 07:56:51
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 "HistoManager.hh"
0032 #include "PrimaryGeneratorAction.hh"
0033
0034 #include "G4PhysicalConstants.hh"
0035 #include "G4SystemOfUnits.hh"
0036 #include "G4UnitsTable.hh"
0037
0038
0039
0040 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0041 {
0042 fParticle = particle;
0043 fEkin = energy;
0044 }
0045
0046
0047
0048 void Run::ParticleCount(G4String name, G4double Ekin, G4double meanLife)
0049 {
0050 std::map<G4String, ParticleData>::iterator it = fParticleDataMap.find(name);
0051 if (it == fParticleDataMap.end()) {
0052 fParticleDataMap[name] = ParticleData(1, Ekin, Ekin, Ekin, meanLife);
0053 }
0054 else {
0055 ParticleData& data = it->second;
0056 data.fCount++;
0057 data.fEmean += Ekin;
0058
0059 G4double emin = data.fEmin;
0060 if (Ekin < emin) data.fEmin = Ekin;
0061 G4double emax = data.fEmax;
0062 if (Ekin > emax) data.fEmax = Ekin;
0063 data.fTmean = meanLife;
0064 }
0065 }
0066
0067
0068
0069 void Run::SetTimeWindow(G4double t1, G4double t2)
0070 {
0071 fTimeWindow1 = t1;
0072 fTimeWindow2 = t2;
0073 }
0074
0075
0076
0077 void Run::CountInTimeWindow(G4String name, G4bool life1, G4bool life2, G4bool decay)
0078 {
0079 std::map<G4String, ActivityData>::iterator it = fActivityMap.find(name);
0080 if (it == fActivityMap.end()) {
0081 G4int n1(0), n2(0), nd(0);
0082 if (life1) n1 = 1;
0083 if (life2) n2 = 1;
0084 if (decay) nd = 1;
0085 fActivityMap[name] = ActivityData(n1, n2, nd);
0086 }
0087 else {
0088 ActivityData& data = it->second;
0089 if (life1) data.fNlife_t1++;
0090 if (life2) data.fNlife_t2++;
0091 if (decay) data.fNdecay_t1t2++;
0092 }
0093 }
0094
0095
0096
0097 void Run::Balance(G4double Ekin, G4double Pbal)
0098 {
0099 fDecayCount++;
0100 fEkinTot[0] += Ekin;
0101
0102 if (fDecayCount == 1) fEkinTot[1] = fEkinTot[2] = Ekin;
0103 if (Ekin < fEkinTot[1]) fEkinTot[1] = Ekin;
0104 if (Ekin > fEkinTot[2]) fEkinTot[2] = Ekin;
0105
0106 fPbalance[0] += Pbal;
0107
0108 if (fDecayCount == 1) fPbalance[1] = fPbalance[2] = Pbal;
0109 if (Pbal < fPbalance[1]) fPbalance[1] = Pbal;
0110 if (Pbal > fPbalance[2]) fPbalance[2] = Pbal;
0111 }
0112
0113
0114
0115 void Run::EventTiming(G4double time)
0116 {
0117 fTimeCount++;
0118 fEventTime[0] += time;
0119 if (fTimeCount == 1) fEventTime[1] = fEventTime[2] = time;
0120 if (time < fEventTime[1]) fEventTime[1] = time;
0121 if (time > fEventTime[2]) fEventTime[2] = time;
0122 }
0123
0124
0125
0126 void Run::PrimaryTiming(G4double ptime)
0127 {
0128 fPrimaryTime += ptime;
0129 }
0130
0131
0132
0133 void Run::EvisEvent(G4double Evis)
0134 {
0135 fEvisEvent[0] += Evis;
0136 if (fTimeCount == 1) fEvisEvent[1] = fEvisEvent[2] = Evis;
0137 if (Evis < fEvisEvent[1]) fEvisEvent[1] = Evis;
0138 if (Evis > fEvisEvent[2]) fEvisEvent[2] = Evis;
0139 }
0140
0141
0142
0143 void Run::Merge(const G4Run* run)
0144 {
0145 const Run* localRun = static_cast<const Run*>(run);
0146
0147
0148
0149 fParticle = localRun->fParticle;
0150 fEkin = localRun->fEkin;
0151
0152
0153
0154 fDecayCount += localRun->fDecayCount;
0155 fTimeCount += localRun->fTimeCount;
0156 fPrimaryTime += localRun->fPrimaryTime;
0157
0158 fEkinTot[0] += localRun->fEkinTot[0];
0159 fPbalance[0] += localRun->fPbalance[0];
0160 fEventTime[0] += localRun->fEventTime[0];
0161 fEvisEvent[0] += localRun->fEvisEvent[0];
0162
0163 G4double min, max;
0164 min = localRun->fEkinTot[1];
0165 max = localRun->fEkinTot[2];
0166 if (fEkinTot[1] > min) fEkinTot[1] = min;
0167 if (fEkinTot[2] < max) fEkinTot[2] = max;
0168
0169 min = localRun->fPbalance[1];
0170 max = localRun->fPbalance[2];
0171 if (fPbalance[1] > min) fPbalance[1] = min;
0172 if (fPbalance[2] < max) fPbalance[2] = max;
0173
0174 min = localRun->fEventTime[1];
0175 max = localRun->fEventTime[2];
0176 if (fEventTime[1] > min) fEventTime[1] = min;
0177 if (fEventTime[2] < max) fEventTime[2] = max;
0178
0179 min = localRun->fEvisEvent[1];
0180 max = localRun->fEvisEvent[2];
0181 if (fEvisEvent[1] > min) fEvisEvent[1] = min;
0182 if (fEvisEvent[2] < max) fEvisEvent[2] = max;
0183
0184
0185 std::map<G4String, ParticleData>::const_iterator itn;
0186 for (itn = localRun->fParticleDataMap.begin(); itn != localRun->fParticleDataMap.end(); ++itn) {
0187 G4String name = itn->first;
0188 const ParticleData& localData = itn->second;
0189 if (fParticleDataMap.find(name) == fParticleDataMap.end()) {
0190 fParticleDataMap[name] = ParticleData(localData.fCount, localData.fEmean, localData.fEmin,
0191 localData.fEmax, localData.fTmean);
0192 }
0193 else {
0194 ParticleData& data = fParticleDataMap[name];
0195 data.fCount += localData.fCount;
0196 data.fEmean += localData.fEmean;
0197 G4double emin = localData.fEmin;
0198 if (emin < data.fEmin) data.fEmin = emin;
0199 G4double emax = localData.fEmax;
0200 if (emax > data.fEmax) data.fEmax = emax;
0201 data.fTmean = localData.fTmean;
0202 }
0203 }
0204
0205
0206 fTimeWindow1 = localRun->fTimeWindow1;
0207 fTimeWindow2 = localRun->fTimeWindow2;
0208
0209 std::map<G4String, ActivityData>::const_iterator ita;
0210 for (ita = localRun->fActivityMap.begin(); ita != localRun->fActivityMap.end(); ++ita) {
0211 G4String name = ita->first;
0212 const ActivityData& localData = ita->second;
0213 if (fActivityMap.find(name) == fActivityMap.end()) {
0214 fActivityMap[name] =
0215 ActivityData(localData.fNlife_t1, localData.fNlife_t2, localData.fNdecay_t1t2);
0216 }
0217 else {
0218 ActivityData& data = fActivityMap[name];
0219 data.fNlife_t1 += localData.fNlife_t1;
0220 data.fNlife_t2 += localData.fNlife_t2;
0221 data.fNdecay_t1t2 += localData.fNdecay_t1t2;
0222 }
0223 }
0224
0225 G4Run::Merge(run);
0226 }
0227
0228
0229
0230 void Run::EndOfRun()
0231 {
0232 G4int nbEvents = numberOfEvent;
0233 G4String partName = fParticle->GetParticleName();
0234
0235 G4cout << "\n ======================== run summary ======================";
0236 G4cout << "\n The run was " << nbEvents << " " << partName << " of "
0237 << G4BestUnit(fEkin, "Energy");
0238 G4cout << "\n ===========================================================\n";
0239 G4cout << G4endl;
0240 if (nbEvents == 0) {
0241 return;
0242 }
0243
0244 G4int prec = 4, wid = prec + 2;
0245 G4int dfprec = G4cout.precision(prec);
0246
0247
0248
0249 G4cout << " Nb of generated particles: \n" << G4endl;
0250
0251 std::map<G4String, ParticleData>::iterator it;
0252 for (it = fParticleDataMap.begin(); it != fParticleDataMap.end(); it++) {
0253 G4String name = it->first;
0254 ParticleData data = it->second;
0255 G4int count = data.fCount;
0256 G4double eMean = data.fEmean / count;
0257 G4double eMin = data.fEmin;
0258 G4double eMax = data.fEmax;
0259 G4double meanLife = data.fTmean;
0260
0261 G4cout << " " << std::setw(15) << name << ": " << std::setw(7) << count
0262 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
0263 << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")";
0264 if (meanLife > 0.)
0265 G4cout << "\tmean life = " << G4BestUnit(meanLife, "Time") << G4endl;
0266 else if (meanLife < 0.)
0267 G4cout << "\tstable" << G4endl;
0268 else
0269 G4cout << G4endl;
0270 }
0271
0272
0273
0274
0275 if (fDecayCount > 0) {
0276 G4double Ebmean = fEkinTot[0] / fDecayCount;
0277 G4double Pbmean = fPbalance[0] / fDecayCount;
0278
0279 G4cout << "\n Ekin Total (Q single decay): mean = " << std::setw(wid)
0280 << G4BestUnit(Ebmean, "Energy") << "\t( " << G4BestUnit(fEkinTot[1], "Energy") << " --> "
0281 << G4BestUnit(fEkinTot[2], "Energy") << ")" << G4endl;
0282
0283 G4cout << "\n Momentum balance (excluding gamma desexcitation): mean = " << std::setw(wid)
0284 << G4BestUnit(Pbmean, "Energy") << "\t( " << G4BestUnit(fPbalance[1], "Energy")
0285 << " --> " << G4BestUnit(fPbalance[2], "Energy") << ")" << G4endl;
0286 }
0287
0288
0289
0290 if (fTimeCount > 0) {
0291 G4double Tmean = fEventTime[0] / fTimeCount;
0292 G4double halfLife = Tmean * std::log(2.);
0293
0294 G4cout << "\n Total time of life (full chain): mean = " << std::setw(wid)
0295 << G4BestUnit(Tmean, "Time") << " half-life = " << std::setw(wid)
0296 << G4BestUnit(halfLife, "Time") << " ( " << G4BestUnit(fEventTime[1], "Time")
0297 << " --> " << G4BestUnit(fEventTime[2], "Time") << ")" << G4endl;
0298 }
0299
0300
0301
0302 if (fTimeCount > 0) {
0303 G4double Evmean = fEvisEvent[0] / fTimeCount;
0304
0305 G4cout << "\n Total visible energy (full chain) : mean = " << std::setw(wid)
0306 << G4BestUnit(Evmean, "Energy") << " ( " << G4BestUnit(fEvisEvent[1], "Energy")
0307 << " --> " << G4BestUnit(fEvisEvent[2], "Energy") << ")" << G4endl;
0308 }
0309
0310
0311
0312 G4double pTimeMean = fPrimaryTime / nbEvents;
0313 G4double molMass = fParticle->GetAtomicMass() * g / mole;
0314 G4double nAtoms_perUnitOfMass = Avogadro / molMass;
0315 G4double Activity_perUnitOfMass = 0.0;
0316 if (pTimeMean > 0.0) {
0317 Activity_perUnitOfMass = nAtoms_perUnitOfMass / pTimeMean;
0318 }
0319
0320 G4cout << "\n Activity of " << partName << " = " << std::setw(wid)
0321 << Activity_perUnitOfMass * g / becquerel << " Bq/g ("
0322 << Activity_perUnitOfMass * g / curie << " Ci/g) \n"
0323 << G4endl;
0324
0325
0326
0327 if (fTimeWindow2 > 0.) {
0328 G4double dt = fTimeWindow2 - fTimeWindow1;
0329 G4cout << " Activities in time window [t1, t2] = [" << G4BestUnit(fTimeWindow1, "Time")
0330 << ", " << G4BestUnit(fTimeWindow2, "Time")
0331 << "] (delta time = " << G4BestUnit(dt, "Time") << ") : \n"
0332 << G4endl;
0333
0334 std::map<G4String, ActivityData>::iterator ita;
0335 for (ita = fActivityMap.begin(); ita != fActivityMap.end(); ita++) {
0336 G4String name = ita->first;
0337 ActivityData data = ita->second;
0338 G4int n1 = data.fNlife_t1;
0339 G4int n2 = data.fNlife_t2;
0340 G4int ndecay = data.fNdecay_t1t2;
0341 G4double actv = ndecay / dt;
0342
0343 G4cout << " " << std::setw(15) << name << ": "
0344 << " n(t1) = " << std::setw(7) << n1 << "\tn(t2) = " << std::setw(7) << n2
0345 << "\t decays = " << std::setw(7) << ndecay
0346 << " ---> <actv> = " << G4BestUnit(actv, "Activity") << "\n";
0347 }
0348 }
0349 G4cout << G4endl;
0350
0351
0352
0353 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0354 G4double factor = 100. / nbEvents;
0355 analysisManager->ScaleH1(1, factor);
0356 analysisManager->ScaleH1(2, factor);
0357 analysisManager->ScaleH1(3, factor);
0358 analysisManager->ScaleH1(4, factor);
0359 analysisManager->ScaleH1(5, factor);
0360
0361
0362
0363 fParticleDataMap.clear();
0364 fActivityMap.clear();
0365
0366
0367
0368 G4cout.precision(dfprec);
0369 }
0370
0371
