File indexing completed on 2025-04-04 08:05:14
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 "G4HadronicProcess.hh"
0040 #include "G4HadronicProcessStore.hh"
0041 #include "G4Neutron.hh"
0042 #include "G4ProcessTable.hh"
0043 #include "G4SystemOfUnits.hh"
0044 #include "G4UnitsTable.hh"
0045
0046
0047
0048 Run::Run(DetectorConstruction* det) : fDetector(det)
0049 {
0050 for (G4int i = 0; i < 3; i++) {
0051 fPbalance[i] = 0.;
0052 }
0053 for (G4int i = 0; i < 3; i++) {
0054 fNbGamma[i] = 0;
0055 }
0056 fPbalance[1] = DBL_MAX;
0057 fNbGamma[1] = 10000;
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::SetTargetXXX(G4bool flag)
0071 {
0072 fTargetXXX = flag;
0073 }
0074
0075
0076
0077 void Run::CountProcesses(G4VProcess* process)
0078 {
0079 if (process == nullptr) return;
0080 G4String procName = process->GetProcessName();
0081 std::map<G4String, G4int>::iterator it = fProcCounter.find(procName);
0082 if (it == fProcCounter.end()) {
0083 fProcCounter[procName] = 1;
0084 }
0085 else {
0086 fProcCounter[procName]++;
0087 }
0088 }
0089
0090
0091
0092 void Run::SumTrack(G4double trackl)
0093 {
0094 fTotalCount++;
0095 fSumTrack += trackl;
0096 fSumTrack2 += trackl * trackl;
0097 }
0098
0099
0100
0101 void Run::CountNuclearChannel(G4String name, G4double Q)
0102 {
0103 std::map<G4String, NuclChannel>::iterator it = fNuclChannelMap.find(name);
0104 if (it == fNuclChannelMap.end()) {
0105 fNuclChannelMap[name] = NuclChannel(1, Q);
0106 }
0107 else {
0108 NuclChannel& data = it->second;
0109 data.fCount++;
0110 data.fQ += Q;
0111 }
0112 }
0113
0114
0115
0116 void Run::ParticleCount(G4String name, G4double Ekin)
0117 {
0118 std::map<G4String, ParticleData>::iterator it = fParticleDataMap.find(name);
0119 if (it == fParticleDataMap.end()) {
0120 fParticleDataMap[name] = ParticleData(1, Ekin, Ekin, Ekin);
0121 }
0122 else {
0123 ParticleData& data = it->second;
0124 data.fCount++;
0125 data.fEmean += Ekin;
0126
0127 G4double emin = data.fEmin;
0128 if (Ekin < emin) data.fEmin = Ekin;
0129 G4double emax = data.fEmax;
0130 if (Ekin > emax) data.fEmax = Ekin;
0131 }
0132 }
0133
0134
0135 void Run::Balance(G4double Pbal)
0136 {
0137 fPbalance[0] += Pbal;
0138
0139 if (fTotalCount == 1) fPbalance[1] = fPbalance[2] = Pbal;
0140 if (Pbal < fPbalance[1]) fPbalance[1] = Pbal;
0141 if (Pbal > fPbalance[2]) fPbalance[2] = Pbal;
0142 }
0143
0144
0145
0146 void Run::CountGamma(G4int nGamma)
0147 {
0148 fGammaCount++;
0149 fNbGamma[0] += nGamma;
0150
0151 if (fGammaCount == 1) fNbGamma[1] = fNbGamma[2] = nGamma;
0152 if (nGamma < fNbGamma[1]) fNbGamma[1] = nGamma;
0153 if (nGamma > fNbGamma[2]) fNbGamma[2] = nGamma;
0154 }
0155
0156
0157
0158 void Run::Merge(const G4Run* run)
0159 {
0160 const Run* localRun = static_cast<const Run*>(run);
0161
0162
0163
0164 fParticle = localRun->fParticle;
0165 fEkin = localRun->fEkin;
0166
0167
0168
0169 fTotalCount += localRun->fTotalCount;
0170 fGammaCount += localRun->fGammaCount;
0171 fSumTrack += localRun->fSumTrack;
0172 fSumTrack2 += localRun->fSumTrack2;
0173
0174 fPbalance[0] += localRun->fPbalance[0];
0175 G4double min, max;
0176 min = localRun->fPbalance[1];
0177 max = localRun->fPbalance[2];
0178 if (fPbalance[1] > min) fPbalance[1] = min;
0179 if (fPbalance[2] < max) fPbalance[2] = max;
0180
0181 fNbGamma[0] += localRun->fNbGamma[0];
0182 G4int nbmin, nbmax;
0183 nbmin = localRun->fNbGamma[1];
0184 nbmax = localRun->fNbGamma[2];
0185 if (fNbGamma[1] > nbmin) fNbGamma[1] = nbmin;
0186 if (fNbGamma[2] < nbmax) fNbGamma[2] = nbmax;
0187
0188
0189 std::map<G4String, G4int>::const_iterator itp;
0190 for (itp = localRun->fProcCounter.begin(); itp != localRun->fProcCounter.end(); ++itp) {
0191 G4String procName = itp->first;
0192 G4int localCount = itp->second;
0193 if (fProcCounter.find(procName) == fProcCounter.end()) {
0194 fProcCounter[procName] = localCount;
0195 }
0196 else {
0197 fProcCounter[procName] += localCount;
0198 }
0199 }
0200
0201
0202 std::map<G4String, NuclChannel>::const_iterator itc;
0203 for (itc = localRun->fNuclChannelMap.begin(); itc != localRun->fNuclChannelMap.end(); ++itc) {
0204 G4String name = itc->first;
0205 const NuclChannel& localData = itc->second;
0206 if (fNuclChannelMap.find(name) == fNuclChannelMap.end()) {
0207 fNuclChannelMap[name] = NuclChannel(localData.fCount, localData.fQ);
0208 }
0209 else {
0210 NuclChannel& data = fNuclChannelMap[name];
0211 data.fCount += localData.fCount;
0212 data.fQ += localData.fQ;
0213 }
0214 }
0215
0216
0217 std::map<G4String, ParticleData>::const_iterator itn;
0218 for (itn = localRun->fParticleDataMap.begin(); itn != localRun->fParticleDataMap.end(); ++itn) {
0219 G4String name = itn->first;
0220 const ParticleData& localData = itn->second;
0221 if (fParticleDataMap.find(name) == fParticleDataMap.end()) {
0222 fParticleDataMap[name] =
0223 ParticleData(localData.fCount, localData.fEmean, localData.fEmin, localData.fEmax);
0224 }
0225 else {
0226 ParticleData& data = fParticleDataMap[name];
0227 data.fCount += localData.fCount;
0228 data.fEmean += localData.fEmean;
0229 G4double emin = localData.fEmin;
0230 if (emin < data.fEmin) data.fEmin = emin;
0231 G4double emax = localData.fEmax;
0232 if (emax > data.fEmax) data.fEmax = emax;
0233 }
0234 }
0235
0236 G4Run::Merge(run);
0237 }
0238
0239
0240
0241 void Run::EndOfRun(G4bool print)
0242 {
0243 G4int prec = 5, wid = prec + 2;
0244 G4int dfprec = G4cout.precision(prec);
0245
0246
0247
0248 const G4Material* material = fDetector->GetMaterial();
0249 G4double density = material->GetDensity();
0250
0251 G4String Particle = fParticle->GetParticleName();
0252 G4cout << "\n The run is " << numberOfEvent << " " << Particle << " of "
0253 << G4BestUnit(fEkin, "Energy") << " through " << G4BestUnit(fDetector->GetSize(), "Length")
0254 << " of " << material->GetName() << " (density: " << G4BestUnit(density, "Volumic Mass")
0255 << ")" << G4endl;
0256
0257 if (numberOfEvent == 0) {
0258 G4cout.precision(dfprec);
0259 return;
0260 }
0261
0262
0263
0264 G4cout << "\n Process calls frequency:" << G4endl;
0265 G4int survive = 0;
0266 std::map<G4String, G4int>::iterator it;
0267 for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
0268 G4String procName = it->first;
0269 G4int count = it->second;
0270 G4cout << "\t" << procName << "= " << count;
0271 if (procName == "Transportation") survive = count;
0272 }
0273 G4cout << G4endl;
0274
0275 if (survive > 0) {
0276 G4cout << "\n Nb of incident particles surviving after "
0277 << G4BestUnit(fDetector->GetSize(), "Length") << " of " << material->GetName() << " : "
0278 << survive << G4endl;
0279 }
0280
0281 if (fTotalCount == 0) fTotalCount = 1;
0282
0283
0284
0285 G4double MeanFreePath = fSumTrack / fTotalCount;
0286 G4double MeanTrack2 = fSumTrack2 / fTotalCount;
0287 G4double rms = std::sqrt(std::fabs(MeanTrack2 - MeanFreePath * MeanFreePath));
0288 G4double CrossSection = 0.0;
0289 if (MeanFreePath > 0.0) {
0290 CrossSection = 1. / MeanFreePath;
0291 }
0292 G4double massicMFP = MeanFreePath * density;
0293 G4double massicCS = 0.0;
0294 if (massicMFP > 0.0) {
0295 massicCS = 1. / massicMFP;
0296 }
0297
0298 G4cout << "\n\n MeanFreePath:\t" << G4BestUnit(MeanFreePath, "Length") << " +- "
0299 << G4BestUnit(rms, "Length") << "\tmassic: " << G4BestUnit(massicMFP, "Mass/Surface")
0300 << "\n CrossSection:\t" << CrossSection * cm << " cm^-1 "
0301 << "\t\tmassic: " << G4BestUnit(massicCS, "Surface/Mass") << G4endl;
0302
0303
0304
0305 if (material->GetNumberOfElements() == 1) {
0306 G4double nbAtoms = material->GetTotNbOfAtomsPerVolume();
0307 G4double crossSection = CrossSection / nbAtoms;
0308 G4cout << " crossSection per atom:\t" << G4BestUnit(crossSection, "Surface") << G4endl;
0309 }
0310
0311
0312 G4cout << "\n Verification: "
0313 << "crossSections from G4HadronicProcessStore" << G4endl;
0314
0315 G4ProcessTable* processTable = G4ProcessTable::GetProcessTable();
0316 G4HadronicProcessStore* store = G4HadronicProcessStore::Instance();
0317 G4double sumc1 = 0.0, sumc2 = 0.0;
0318 const G4Element* element =
0319 (material->GetNumberOfElements() == 1) ? material->GetElement(0) : nullptr;
0320 for (it = fProcCounter.begin(); it != fProcCounter.end(); ++it) {
0321 G4String procName = it->first;
0322 const G4VProcess* process = processTable->FindProcess(procName, fParticle);
0323 PrintXS(process, material, element, store, density, sumc1, sumc2);
0324 }
0325 if (sumc1 > 0.0) {
0326 G4cout << "\n"
0327 << std::setw(20) << "total"
0328 << " = " << G4BestUnit(sumc1, "Surface/Mass") << "\t";
0329 if (sumc2 > 0.0) {
0330 G4cout << G4BestUnit(sumc2, "Surface");
0331 }
0332 G4cout << G4endl;
0333 }
0334 else {
0335 G4cout << " not available" << G4endl;
0336 }
0337
0338
0339
0340 G4cout << "\n List of nuclear reactions: \n" << G4endl;
0341 std::map<G4String, NuclChannel>::iterator ic;
0342 for (ic = fNuclChannelMap.begin(); ic != fNuclChannelMap.end(); ic++) {
0343 G4String name = ic->first;
0344 NuclChannel data = ic->second;
0345 G4int count = data.fCount;
0346 G4double Q = data.fQ / count;
0347 if (print)
0348 G4cout << " " << std::setw(60) << name << ": " << std::setw(7) << count
0349 << " Q = " << std::setw(wid) << G4BestUnit(Q, "Energy") << G4endl;
0350 }
0351
0352
0353
0354 if (print && (fGammaCount > 0)) {
0355 G4cout << "\n"
0356 << std::setw(58) << "number of gamma or e- (ic): N = " << fNbGamma[1] << " --> "
0357 << fNbGamma[2] << G4endl;
0358 }
0359
0360 if (print && fTargetXXX) {
0361 G4cout << "\n --> NOTE: XXXX because neutronHP is unable to return target nucleus" << G4endl;
0362 }
0363
0364
0365
0366 G4cout << "\n List of generated particles:" << G4endl;
0367
0368 std::map<G4String, ParticleData>::iterator itn;
0369 for (itn = fParticleDataMap.begin(); itn != fParticleDataMap.end(); itn++) {
0370 G4String name = itn->first;
0371 ParticleData data = itn->second;
0372 G4int count = data.fCount;
0373 G4double eMean = data.fEmean / count;
0374 G4double eMin = data.fEmin;
0375 G4double eMax = data.fEmax;
0376 if (print)
0377 G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
0378 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
0379 << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")"
0380 << G4endl;
0381 }
0382
0383
0384
0385 if (fTotalCount > 1) {
0386 G4double Pbmean = fPbalance[0] / fTotalCount;
0387 G4cout << "\n Momentum balance: Pmean = " << std::setw(wid) << G4BestUnit(Pbmean, "Energy")
0388 << "\t( " << G4BestUnit(fPbalance[1], "Energy") << " --> "
0389 << G4BestUnit(fPbalance[2], "Energy") << ") \n"
0390 << G4endl;
0391 }
0392
0393
0394
0395
0396
0397
0398
0399 fProcCounter.clear();
0400 fNuclChannelMap.clear();
0401 fParticleDataMap.clear();
0402
0403
0404 G4cout.precision(dfprec);
0405 }
0406
0407
0408
0409 void Run::PrintXS(const G4VProcess* proc, const G4Material* mat, const G4Element* elm,
0410 G4HadronicProcessStore* store, G4double density, G4double& sum1, G4double& sum2)
0411 {
0412 if (nullptr == proc) {
0413 return;
0414 }
0415 G4double xs1 = store->GetCrossSectionPerVolume(fParticle, fEkin, proc, mat);
0416 G4double massSigma = xs1 / density;
0417 sum1 += massSigma;
0418 if (nullptr != elm) {
0419 G4double xs2 = store->GetCrossSectionPerAtom(fParticle, fEkin, proc, elm, mat);
0420 sum2 += xs2;
0421 G4cout << "\n"
0422 << std::setw(20) << proc->GetProcessName() << " = "
0423 << G4BestUnit(massSigma, "Surface/Mass") << "\t" << G4BestUnit(xs2, "Surface");
0424 }
0425 else {
0426 G4cout << "\n"
0427 << std::setw(20) << proc->GetProcessName() << " = "
0428 << G4BestUnit(massSigma, "Surface/Mass");
0429 }
0430 }
0431
0432
