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