File indexing completed on 2025-04-04 08:05:17
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 "G4SystemOfUnits.hh"
0040 #include "G4UnitsTable.hh"
0041
0042
0043
0044 Run::Run(DetectorConstruction* det) : fDetector(det) {}
0045
0046
0047
0048 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0049 {
0050 fParticle = particle;
0051 fEkin = energy;
0052 }
0053
0054
0055
0056 void Run::CountProcesses(const G4VProcess* process)
0057 {
0058 if (process == nullptr) return;
0059 G4String procName = process->GetProcessName();
0060 std::map<G4String, G4int>::iterator it = fProcCounter.find(procName);
0061 if (it == fProcCounter.end()) {
0062 fProcCounter[procName] = 1;
0063 }
0064 else {
0065 fProcCounter[procName]++;
0066 }
0067 }
0068
0069
0070
0071 void Run::ParticleCount(G4String name, G4double Ekin, G4double meanLife)
0072 {
0073 std::map<G4String, ParticleData>::iterator it = fParticleDataMap1.find(name);
0074 if (it == fParticleDataMap1.end()) {
0075 fParticleDataMap1[name] = ParticleData(1, Ekin, Ekin, Ekin, meanLife);
0076 }
0077 else {
0078 ParticleData& data = it->second;
0079 data.fCount++;
0080 data.fEmean += Ekin;
0081
0082 G4double emin = data.fEmin;
0083 if (Ekin < emin) data.fEmin = Ekin;
0084 G4double emax = data.fEmax;
0085 if (Ekin > emax) data.fEmax = Ekin;
0086 data.fTmean = meanLife;
0087 }
0088 }
0089
0090
0091
0092 void Run::SumEnergies(G4double edep, G4double eflow, G4double etot)
0093 {
0094 fEnergyDeposit += edep;
0095 fEnergyDeposit2 += edep * edep;
0096
0097 fEnergyFlow += eflow;
0098 fEnergyFlow2 += eflow * eflow;
0099
0100 fEnergyTotal += etot;
0101 fEnergyTotal2 += etot * etot;
0102 }
0103
0104
0105
0106 void Run::ParticleFlux(G4String name, G4double Ekin)
0107 {
0108 std::map<G4String, ParticleData>::iterator it = fParticleDataMap2.find(name);
0109 if (it == fParticleDataMap2.end()) {
0110 fParticleDataMap2[name] = ParticleData(1, Ekin, Ekin, Ekin, -1 * ns);
0111 }
0112 else {
0113 ParticleData& data = it->second;
0114 data.fCount++;
0115 data.fEmean += Ekin;
0116
0117 G4double emin = data.fEmin;
0118 if (Ekin < emin) data.fEmin = Ekin;
0119 G4double emax = data.fEmax;
0120 if (Ekin > emax) data.fEmax = Ekin;
0121 data.fTmean = -1 * ns;
0122 }
0123 }
0124
0125
0126
0127 void Run::Merge(const G4Run* run)
0128 {
0129 const Run* localRun = static_cast<const Run*>(run);
0130
0131
0132
0133 fParticle = localRun->fParticle;
0134 fEkin = localRun->fEkin;
0135
0136
0137
0138 fEnergyDeposit += localRun->fEnergyDeposit;
0139 fEnergyDeposit2 += localRun->fEnergyDeposit2;
0140 fEnergyFlow += localRun->fEnergyFlow;
0141 fEnergyFlow2 += localRun->fEnergyFlow2;
0142 fEnergyTotal += localRun->fEnergyTotal;
0143 fEnergyTotal2 += localRun->fEnergyTotal2;
0144
0145
0146 std::map<G4String, G4int>::const_iterator itp;
0147 for (itp = localRun->fProcCounter.begin(); itp != localRun->fProcCounter.end(); ++itp) {
0148 G4String procName = itp->first;
0149 G4int localCount = itp->second;
0150 if (fProcCounter.find(procName) == fProcCounter.end()) {
0151 fProcCounter[procName] = localCount;
0152 }
0153 else {
0154 fProcCounter[procName] += localCount;
0155 }
0156 }
0157
0158
0159 std::map<G4String, ParticleData>::const_iterator itc;
0160 for (itc = localRun->fParticleDataMap1.begin(); itc != localRun->fParticleDataMap1.end(); ++itc) {
0161 G4String name = itc->first;
0162 const ParticleData& localData = itc->second;
0163 if (fParticleDataMap1.find(name) == fParticleDataMap1.end()) {
0164 fParticleDataMap1[name] = ParticleData(localData.fCount, localData.fEmean, localData.fEmin,
0165 localData.fEmax, localData.fTmean);
0166 }
0167 else {
0168 ParticleData& data = fParticleDataMap1[name];
0169 data.fCount += localData.fCount;
0170 data.fEmean += localData.fEmean;
0171 G4double emin = localData.fEmin;
0172 if (emin < data.fEmin) data.fEmin = emin;
0173 G4double emax = localData.fEmax;
0174 if (emax > data.fEmax) data.fEmax = emax;
0175 data.fTmean = localData.fTmean;
0176 }
0177 }
0178
0179
0180 std::map<G4String, ParticleData>::const_iterator itn;
0181 for (itn = localRun->fParticleDataMap2.begin(); itn != localRun->fParticleDataMap2.end(); ++itn) {
0182 G4String name = itn->first;
0183 const ParticleData& localData = itn->second;
0184 if (fParticleDataMap2.find(name) == fParticleDataMap2.end()) {
0185 fParticleDataMap2[name] = ParticleData(localData.fCount, localData.fEmean, localData.fEmin,
0186 localData.fEmax, localData.fTmean);
0187 }
0188 else {
0189 ParticleData& data = fParticleDataMap2[name];
0190 data.fCount += localData.fCount;
0191 data.fEmean += localData.fEmean;
0192 G4double emin = localData.fEmin;
0193 if (emin < data.fEmin) data.fEmin = emin;
0194 G4double emax = localData.fEmax;
0195 if (emax > data.fEmax) data.fEmax = emax;
0196 data.fTmean = localData.fTmean;
0197 }
0198 }
0199
0200 G4Run::Merge(run);
0201 }
0202
0203
0204
0205 void Run::EndOfRun()
0206 {
0207 G4int prec = 5, wid = prec + 2;
0208 G4int dfprec = G4cout.precision(prec);
0209
0210
0211
0212 G4Material* material = fDetector->GetMaterial();
0213 G4double density = material->GetDensity();
0214
0215 G4String Particle = fParticle->GetParticleName();
0216 G4cout << "\n The run is " << numberOfEvent << " " << Particle << " of "
0217 << G4BestUnit(fEkin, "Energy") << " through "
0218 << G4BestUnit(fDetector->GetRadius(), "Length") << " of " << material->GetName()
0219 << " (density: " << G4BestUnit(density, "Volumic Mass") << ")" << G4endl;
0220
0221 if (numberOfEvent == 0) {
0222 G4cout.precision(dfprec);
0223 return;
0224 }
0225
0226
0227
0228 G4cout << "\n Process calls frequency :" << G4endl;
0229 G4int index = 0;
0230 std::map<G4String, G4int>::iterator it;
0231 for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
0232 G4String procName = it->first;
0233 G4int count = it->second;
0234 G4String space = " ";
0235 if (++index % 3 == 0) space = "\n";
0236 G4cout << " " << std::setw(20) << procName << "=" << std::setw(7) << count << space;
0237 }
0238 G4cout << G4endl;
0239
0240
0241
0242 G4int TotNbofEvents = numberOfEvent;
0243 fEnergyDeposit /= TotNbofEvents;
0244 fEnergyDeposit2 /= TotNbofEvents;
0245 G4double rmsEdep = fEnergyDeposit2 - fEnergyDeposit * fEnergyDeposit;
0246 if (rmsEdep > 0.)
0247 rmsEdep = std::sqrt(rmsEdep);
0248 else
0249 rmsEdep = 0.;
0250
0251 G4cout << "\n Mean energy deposit per event = " << G4BestUnit(fEnergyDeposit, "Energy")
0252 << "; rms = " << G4BestUnit(rmsEdep, "Energy") << G4endl;
0253
0254
0255
0256 fEnergyFlow /= TotNbofEvents;
0257 fEnergyFlow2 /= TotNbofEvents;
0258 G4double rmsEflow = fEnergyFlow2 - fEnergyFlow * fEnergyFlow;
0259 if (rmsEflow > 0.)
0260 rmsEflow = std::sqrt(rmsEflow);
0261 else
0262 rmsEflow = 0.;
0263
0264 G4cout << " Mean energy leakage per event = " << G4BestUnit(fEnergyFlow, "Energy")
0265 << "; rms = " << G4BestUnit(rmsEflow, "Energy") << G4endl;
0266
0267
0268
0269 fEnergyTotal /= TotNbofEvents;
0270 fEnergyTotal2 /= TotNbofEvents;
0271 G4double rmsEtotal = fEnergyTotal2 - fEnergyTotal * fEnergyTotal;
0272 if (rmsEtotal > 0.)
0273 rmsEtotal = std::sqrt(rmsEtotal);
0274 else
0275 rmsEflow = 0.;
0276
0277 G4cout << "\n Mean energy total per event = " << G4BestUnit(fEnergyTotal, "Energy")
0278 << "; rms = " << G4BestUnit(rmsEtotal, "Energy") << G4endl;
0279
0280
0281
0282 if (fParticleDataMap1.size() > 0) {
0283 G4cout << "\n List of particles at creation :" << G4endl;
0284 std::map<G4String, ParticleData>::iterator itc;
0285 for (itc = fParticleDataMap1.begin(); itc != fParticleDataMap1.end(); itc++) {
0286 G4String name = itc->first;
0287 ParticleData data = itc->second;
0288 G4int count = data.fCount;
0289 G4double eMean = data.fEmean / count;
0290 G4double eMin = data.fEmin;
0291 G4double eMax = data.fEmax;
0292 G4double meanLife = data.fTmean;
0293
0294 G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
0295 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
0296 << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")";
0297 if (meanLife >= 0.)
0298 G4cout << "\tmean life = " << G4BestUnit(meanLife, "Time") << G4endl;
0299 else
0300 G4cout << "\tstable" << G4endl;
0301 }
0302 }
0303
0304
0305
0306 G4cout << "\n List of particles emerging from the absorber :" << G4endl;
0307
0308 std::map<G4String, ParticleData>::iterator itn;
0309 for (itn = fParticleDataMap2.begin(); itn != fParticleDataMap2.end(); itn++) {
0310 G4String name = itn->first;
0311 ParticleData data = itn->second;
0312 G4int count = data.fCount;
0313 G4double eMean = data.fEmean / count;
0314 G4double eMin = data.fEmin;
0315 G4double eMax = data.fEmax;
0316 G4double Eflow = data.fEmean / TotNbofEvents;
0317
0318 G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
0319 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
0320 << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy")
0321 << ") \tEleak/event = " << G4BestUnit(Eflow, "Energy") << G4endl;
0322 }
0323
0324
0325 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0326
0327 G4int ih = 2;
0328 G4double binWidth = analysisManager->GetH1Width(ih);
0329 G4double fac = (1. / (numberOfEvent * binWidth)) * (mm / MeV);
0330 analysisManager->ScaleH1(ih, fac);
0331
0332
0333 fProcCounter.clear();
0334 fParticleDataMap2.clear();
0335
0336
0337 G4cout.precision(dfprec);
0338 }
0339
0340
