Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/hadronic/Hadr03/src/Run.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 /// \file Run.cc
 27 /// \brief Implementation of the Run class
 28 //
 29 //
 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 32 
 33 #include "Run.hh"
 34 
 35 #include "DetectorConstruction.hh"
 36 #include "HistoManager.hh"
 37 #include "PrimaryGeneratorAction.hh"
 38 
 39 #include "G4HadronicProcess.hh"
 40 #include "G4HadronicProcessStore.hh"
 41 #include "G4Neutron.hh"
 42 #include "G4ProcessTable.hh"
 43 #include "G4SystemOfUnits.hh"
 44 #include "G4UnitsTable.hh"
 45 
 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 47 
 48 Run::Run(DetectorConstruction* det) : fDetector(det)
 49 {
 50   for (G4int i = 0; i < 3; i++) {
 51     fPbalance[i] = 0.;
 52   }
 53   for (G4int i = 0; i < 3; i++) {
 54     fNbGamma[i] = 0;
 55   }
 56   fPbalance[1] = DBL_MAX;
 57   fNbGamma[1] = 10000;
 58 }
 59 
 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 61 
 62 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
 63 {
 64   fParticle = particle;
 65   fEkin = energy;
 66 }
 67 
 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 69 
 70 void Run::SetTargetXXX(G4bool flag)
 71 {
 72   fTargetXXX = flag;
 73 }
 74 
 75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 76 
 77 void Run::CountProcesses(G4VProcess* process)
 78 {
 79   if (process == nullptr) return;
 80   G4String procName = process->GetProcessName();
 81   std::map<G4String, G4int>::iterator it = fProcCounter.find(procName);
 82   if (it == fProcCounter.end()) {
 83     fProcCounter[procName] = 1;
 84   }
 85   else {
 86     fProcCounter[procName]++;
 87   }
 88 }
 89 
 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 91 
 92 void Run::SumTrack(G4double trackl)
 93 {
 94   fTotalCount++;
 95   fSumTrack += trackl;
 96   fSumTrack2 += trackl * trackl;
 97 }
 98 
 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100 
101 void Run::CountNuclearChannel(G4String name, G4double Q)
102 {
103   std::map<G4String, NuclChannel>::iterator it = fNuclChannelMap.find(name);
104   if (it == fNuclChannelMap.end()) {
105     fNuclChannelMap[name] = NuclChannel(1, Q);
106   }
107   else {
108     NuclChannel& data = it->second;
109     data.fCount++;
110     data.fQ += Q;
111   }
112 }
113 
114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
115 
116 void Run::ParticleCount(G4String name, G4double Ekin)
117 {
118   std::map<G4String, ParticleData>::iterator it = fParticleDataMap.find(name);
119   if (it == fParticleDataMap.end()) {
120     fParticleDataMap[name] = ParticleData(1, Ekin, Ekin, Ekin);
121   }
122   else {
123     ParticleData& data = it->second;
124     data.fCount++;
125     data.fEmean += Ekin;
126     // update min max
127     G4double emin = data.fEmin;
128     if (Ekin < emin) data.fEmin = Ekin;
129     G4double emax = data.fEmax;
130     if (Ekin > emax) data.fEmax = Ekin;
131   }
132 }
133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
134 
135 void Run::Balance(G4double Pbal)
136 {
137   fPbalance[0] += Pbal;
138   // update min max
139   if (fTotalCount == 1) fPbalance[1] = fPbalance[2] = Pbal;
140   if (Pbal < fPbalance[1]) fPbalance[1] = Pbal;
141   if (Pbal > fPbalance[2]) fPbalance[2] = Pbal;
142 }
143 
144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
145 
146 void Run::CountGamma(G4int nGamma)
147 {
148   fGammaCount++;
149   fNbGamma[0] += nGamma;
150   // update min max
151   if (fGammaCount == 1) fNbGamma[1] = fNbGamma[2] = nGamma;
152   if (nGamma < fNbGamma[1]) fNbGamma[1] = nGamma;
153   if (nGamma > fNbGamma[2]) fNbGamma[2] = nGamma;
154 }
155 
156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
157 
158 void Run::Merge(const G4Run* run)
159 {
160   const Run* localRun = static_cast<const Run*>(run);
161 
162   // primary particle info
163   //
164   fParticle = localRun->fParticle;
165   fEkin = localRun->fEkin;
166 
167   // accumulate sums
168   //
169   fTotalCount += localRun->fTotalCount;
170   fGammaCount += localRun->fGammaCount;
171   fSumTrack += localRun->fSumTrack;
172   fSumTrack2 += localRun->fSumTrack2;
173 
174   fPbalance[0] += localRun->fPbalance[0];
175   G4double min, max;
176   min = localRun->fPbalance[1];
177   max = localRun->fPbalance[2];
178   if (fPbalance[1] > min) fPbalance[1] = min;
179   if (fPbalance[2] < max) fPbalance[2] = max;
180 
181   fNbGamma[0] += localRun->fNbGamma[0];
182   G4int nbmin, nbmax;
183   nbmin = localRun->fNbGamma[1];
184   nbmax = localRun->fNbGamma[2];
185   if (fNbGamma[1] > nbmin) fNbGamma[1] = nbmin;
186   if (fNbGamma[2] < nbmax) fNbGamma[2] = nbmax;
187 
188   // map: processes count
189   std::map<G4String, G4int>::const_iterator itp;
190   for (itp = localRun->fProcCounter.begin(); itp != localRun->fProcCounter.end(); ++itp) {
191     G4String procName = itp->first;
192     G4int localCount = itp->second;
193     if (fProcCounter.find(procName) == fProcCounter.end()) {
194       fProcCounter[procName] = localCount;
195     }
196     else {
197       fProcCounter[procName] += localCount;
198     }
199   }
200 
201   // map: nuclear channels
202   std::map<G4String, NuclChannel>::const_iterator itc;
203   for (itc = localRun->fNuclChannelMap.begin(); itc != localRun->fNuclChannelMap.end(); ++itc) {
204     G4String name = itc->first;
205     const NuclChannel& localData = itc->second;
206     if (fNuclChannelMap.find(name) == fNuclChannelMap.end()) {
207       fNuclChannelMap[name] = NuclChannel(localData.fCount, localData.fQ);
208     }
209     else {
210       NuclChannel& data = fNuclChannelMap[name];
211       data.fCount += localData.fCount;
212       data.fQ += localData.fQ;
213     }
214   }
215 
216   // map: particles count
217   std::map<G4String, ParticleData>::const_iterator itn;
218   for (itn = localRun->fParticleDataMap.begin(); itn != localRun->fParticleDataMap.end(); ++itn) {
219     G4String name = itn->first;
220     const ParticleData& localData = itn->second;
221     if (fParticleDataMap.find(name) == fParticleDataMap.end()) {
222       fParticleDataMap[name] =
223         ParticleData(localData.fCount, localData.fEmean, localData.fEmin, localData.fEmax);
224     }
225     else {
226       ParticleData& data = fParticleDataMap[name];
227       data.fCount += localData.fCount;
228       data.fEmean += localData.fEmean;
229       G4double emin = localData.fEmin;
230       if (emin < data.fEmin) data.fEmin = emin;
231       G4double emax = localData.fEmax;
232       if (emax > data.fEmax) data.fEmax = emax;
233     }
234   }
235 
236   G4Run::Merge(run);
237 }
238 
239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
240 
241 void Run::EndOfRun(G4bool print)
242 {
243   G4int prec = 5, wid = prec + 2;
244   G4int dfprec = G4cout.precision(prec);
245 
246   // run condition
247   //
248   const G4Material* material = fDetector->GetMaterial();
249   G4double density = material->GetDensity();
250 
251   G4String Particle = fParticle->GetParticleName();
252   G4cout << "\n The run is " << numberOfEvent << " " << Particle << " of "
253          << G4BestUnit(fEkin, "Energy") << " through " << G4BestUnit(fDetector->GetSize(), "Length")
254          << " of " << material->GetName() << " (density: " << G4BestUnit(density, "Volumic Mass")
255          << ")" << G4endl;
256 
257   if (numberOfEvent == 0) {
258     G4cout.precision(dfprec);
259     return;
260   }
261 
262   // frequency of processes
263   //
264   G4cout << "\n Process calls frequency:" << G4endl;
265   G4int survive = 0;
266   std::map<G4String, G4int>::iterator it;
267   for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
268     G4String procName = it->first;
269     G4int count = it->second;
270     G4cout << "\t" << procName << "= " << count;
271     if (procName == "Transportation") survive = count;
272   }
273   G4cout << G4endl;
274 
275   if (survive > 0) {
276     G4cout << "\n Nb of incident particles surviving after "
277            << G4BestUnit(fDetector->GetSize(), "Length") << " of " << material->GetName() << " : "
278            << survive << G4endl;
279   }
280 
281   if (fTotalCount == 0) fTotalCount = 1;  // force printing anyway
282 
283   // compute mean free path and related quantities
284   //
285   G4double MeanFreePath = fSumTrack / fTotalCount;
286   G4double MeanTrack2 = fSumTrack2 / fTotalCount;
287   G4double rms = std::sqrt(std::fabs(MeanTrack2 - MeanFreePath * MeanFreePath));
288   G4double CrossSection = 0.0;
289   if (MeanFreePath > 0.0) {
290     CrossSection = 1. / MeanFreePath;
291   }
292   G4double massicMFP = MeanFreePath * density;
293   G4double massicCS = 0.0;
294   if (massicMFP > 0.0) {
295     massicCS = 1. / massicMFP;
296   }
297 
298   G4cout << "\n\n MeanFreePath:\t" << G4BestUnit(MeanFreePath, "Length") << " +- "
299          << G4BestUnit(rms, "Length") << "\tmassic: " << G4BestUnit(massicMFP, "Mass/Surface")
300          << "\n CrossSection:\t" << CrossSection * cm << " cm^-1 "
301          << "\t\tmassic: " << G4BestUnit(massicCS, "Surface/Mass") << G4endl;
302 
303   // cross section per atom (only for single material)
304   //
305   if (material->GetNumberOfElements() == 1) {
306     G4double nbAtoms = material->GetTotNbOfAtomsPerVolume();
307     G4double crossSection = CrossSection / nbAtoms;
308     G4cout << " crossSection per atom:\t" << G4BestUnit(crossSection, "Surface") << G4endl;
309   }
310   // check cross section from G4HadronicProcessStore
311   //
312   G4cout << "\n Verification: "
313          << "crossSections from G4HadronicProcessStore" << G4endl;
314 
315   G4ProcessTable* processTable = G4ProcessTable::GetProcessTable();
316   G4HadronicProcessStore* store = G4HadronicProcessStore::Instance();
317   G4double sumc1 = 0.0, sumc2 = 0.0;
318   const G4Element* element =
319     (material->GetNumberOfElements() == 1) ? material->GetElement(0) : nullptr;
320   for (it = fProcCounter.begin(); it != fProcCounter.end(); ++it) {
321     G4String procName = it->first;
322     const G4VProcess* process = processTable->FindProcess(procName, fParticle);
323     PrintXS(process, material, element, store, density, sumc1, sumc2);
324   }
325   if (sumc1 > 0.0) {
326     G4cout << "\n"
327            << std::setw(20) << "total"
328            << " = " << G4BestUnit(sumc1, "Surface/Mass") << "\t";
329     if (sumc2 > 0.0) {
330       G4cout << G4BestUnit(sumc2, "Surface");
331     }
332     G4cout << G4endl;
333   }
334   else {
335     G4cout << " not available" << G4endl;
336   }
337 
338   // nuclear channel count
339   //
340   G4cout << "\n List of nuclear reactions: \n" << G4endl;
341   std::map<G4String, NuclChannel>::iterator ic;
342   for (ic = fNuclChannelMap.begin(); ic != fNuclChannelMap.end(); ic++) {
343     G4String name = ic->first;
344     NuclChannel data = ic->second;
345     G4int count = data.fCount;
346     G4double Q = data.fQ / count;
347     if (print)
348       G4cout << "  " << std::setw(60) << name << ": " << std::setw(7) << count
349              << "   Q = " << std::setw(wid) << G4BestUnit(Q, "Energy") << G4endl;
350   }
351 
352   // Gamma count
353   //
354   if (print && (fGammaCount > 0)) {
355     G4cout << "\n"
356            << std::setw(58) << "number of gamma or e- (ic): N = " << fNbGamma[1] << " --> "
357            << fNbGamma[2] << G4endl;
358   }
359 
360   if (print && fTargetXXX) {
361     G4cout << "\n   --> NOTE: XXXX because neutronHP is unable to return target nucleus" << G4endl;
362   }
363 
364   // particles count
365   //
366   G4cout << "\n List of generated particles:" << G4endl;
367 
368   std::map<G4String, ParticleData>::iterator itn;
369   for (itn = fParticleDataMap.begin(); itn != fParticleDataMap.end(); itn++) {
370     G4String name = itn->first;
371     ParticleData data = itn->second;
372     G4int count = data.fCount;
373     G4double eMean = data.fEmean / count;
374     G4double eMin = data.fEmin;
375     G4double eMax = data.fEmax;
376     if (print)
377       G4cout << "  " << std::setw(13) << name << ": " << std::setw(7) << count
378              << "  Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
379              << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")"
380              << G4endl;
381   }
382 
383   // energy momentum balance
384   //
385   if (fTotalCount > 1) {
386     G4double Pbmean = fPbalance[0] / fTotalCount;
387     G4cout << "\n   Momentum balance: Pmean = " << std::setw(wid) << G4BestUnit(Pbmean, "Energy")
388            << "\t( " << G4BestUnit(fPbalance[1], "Energy") << " --> "
389            << G4BestUnit(fPbalance[2], "Energy") << ") \n"
390            << G4endl;
391   }
392 
393   // normalize histograms
394   ////G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
395   ////G4double factor = 1./numberOfEvent;
396   ////analysisManager->ScaleH1(3,factor);
397 
398   // remove all contents in fProcCounter, fCount
399   fProcCounter.clear();
400   fNuclChannelMap.clear();
401   fParticleDataMap.clear();
402 
403   // restore default format
404   G4cout.precision(dfprec);
405 }
406 
407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
408 
409 void Run::PrintXS(const G4VProcess* proc, const G4Material* mat, const G4Element* elm,
410                   G4HadronicProcessStore* store, G4double density, G4double& sum1, G4double& sum2)
411 {
412   if (nullptr == proc) {
413     return;
414   }
415   G4double xs1 = store->GetCrossSectionPerVolume(fParticle, fEkin, proc, mat);
416   G4double massSigma = xs1 / density;
417   sum1 += massSigma;
418   if (nullptr != elm) {
419     G4double xs2 = store->GetCrossSectionPerAtom(fParticle, fEkin, proc, elm, mat);
420     sum2 += xs2;
421     G4cout << "\n"
422            << std::setw(20) << proc->GetProcessName() << " = "
423            << G4BestUnit(massSigma, "Surface/Mass") << "\t" << G4BestUnit(xs2, "Surface");
424   }
425   else {
426     G4cout << "\n"
427            << std::setw(20) << proc->GetProcessName() << " = "
428            << G4BestUnit(massSigma, "Surface/Mass");
429   }
430 }
431 
432 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
433