LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
RunAction.cc
Go to the documentation of this file.
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 //
27 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
28 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
29 
30 #include "RunAction.hh"
31 #include "DetectorConstruction.hh"
33 
34 #include "G4Run.hh"
35 #include "G4ProcessManager.hh"
36 #include "G4UnitsTable.hh"
37 #include "G4EmCalculator.hh"
38 #include "G4Electron.hh"
39 #include "G4SystemOfUnits.hh"
40 
41 #include <vector>
42 
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44 
46 :detector(det), primary(kin)
47 { }
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 
52 { }
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55 
56 void RunAction::BeginOfRunAction(const G4Run*)
57 {
58  //set precision for printing
59  G4int prec = G4cout.precision(6);
60 
61  // get particle
62  G4ParticleDefinition* particle = primary->GetParticleGun()
63  ->GetParticleDefinition();
64  G4String partName = particle->GetParticleName();
65  G4double charge = particle->GetPDGCharge();
66  G4double energy = primary->GetParticleGun()->GetParticleEnergy();
67 
68  // get material
69  G4Material* material = detector->GetMaterial();
70  G4String matName = material->GetName();
71  G4double density = material->GetDensity();
72  G4double radl = material->GetRadlen();
73 
74  G4cout << "\n " << partName << " ("
75  << G4BestUnit(energy,"Energy") << ") in "
76  << material->GetName() << " (density: "
77  << G4BestUnit(density,"Volumic Mass") << "; radiation length: "
78  << G4BestUnit(radl, "Length") << ")" << G4endl;
79 
80  // get cuts
81  GetCuts();
82  if (charge != 0.) {
83  G4cout << "\n Range cuts : \t gamma "
84  << std::setw(8) << G4BestUnit(rangeCut[0],"Length")
85  << "\t e- " << std::setw(8) << G4BestUnit(rangeCut[1],"Length");
86  G4cout << "\n Energy cuts : \t gamma "
87  << std::setw(8) << G4BestUnit(energyCut[0],"Energy")
88  << "\t e- " << std::setw(8) << G4BestUnit(energyCut[1],"Energy")
89  << G4endl;
90  }
91 
92  // get processList and extract EM processes (but not MultipleScattering)
93  G4ProcessVector* plist = particle->GetProcessManager()->GetProcessList();
94  G4String procName;
95  G4double cut;
96  std::vector<G4String> emName;
97  std::vector<G4double> enerCut;
98  size_t length = plist->size();
99  for (size_t j=0; j<length; j++) {
100  procName = (*plist)[j]->GetProcessName();
101  cut = energyCut[1];
102  if ((procName == "eBrem")||(procName == "muBrems")) cut = energyCut[0];
103  if (((*plist)[j]->GetProcessType() == fElectromagnetic) &&
104  (procName != "msc")) {
105  emName.push_back(procName);
106  enerCut.push_back(cut);
107  }
108  }
109 
110  // print list of processes
111  G4cout << "\n processes : ";
112  for (size_t j=0; j<emName.size();j++)
113  G4cout << "\t" << std::setw(13) << emName[j] << "\t";
114  G4cout << "\t" << std::setw(13) <<"total";
115 
116  //instanciate EmCalculator
117  G4EmCalculator emCal;
118  // emCal.SetVerbose(2);
119 
120  //compute cross section per atom (only for single material)
121  if (material->GetNumberOfElements() == 1) {
122  G4double Z = material->GetZ();
123  G4double A = material->GetA();
124 
125  std::vector<G4double> sigma0;
126  G4double sig, sigtot = 0.;
127 
128  for (size_t j=0; j<emName.size();j++) {
129  sig = emCal.ComputeCrossSectionPerAtom
130  (energy,particle,emName[j],Z,A,enerCut[j]);
131  sigtot += sig;
132  sigma0.push_back(sig);
133  }
134  sigma0.push_back(sigtot);
135 
136  G4cout << "\n \n cross section per atom : ";
137  for (size_t j=0; j<sigma0.size();j++) {
138  G4cout << "\t" << std::setw(13) << G4BestUnit(sigma0[j], "Surface");
139  }
140  G4cout << G4endl;
141  }
142 
143  //get cross section per volume
144  std::vector<G4double> sigma1;
145  std::vector<G4double> sigma2;
146  G4double Sig, Sigtot = 0.;
147 
148  for (size_t j=0; j<emName.size();j++) {
149  Sig = emCal.GetCrossSectionPerVolume(energy,particle,emName[j],material);
150  if (Sig == 0.) Sig = emCal.ComputeCrossSectionPerVolume
151  (energy,particle,emName[j],material,enerCut[j]);
152  Sigtot += Sig;
153  sigma1.push_back(Sig);
154  sigma2.push_back(Sig/density);
155  }
156  sigma1.push_back(Sigtot);
157  sigma2.push_back(Sigtot/density);
158 
159  //print cross sections
160  G4cout << "\n \n cross section per volume : ";
161  for (size_t j=0; j<sigma1.size();j++) {
162  G4cout << "\t" << std::setw(13) << sigma1[j]*cm << " cm^-1";
163  }
164 
165  G4cout << "\n cross section per mass : ";
166  for (size_t j=0; j<sigma2.size();j++) {
167  G4cout << "\t" << std::setw(13) << G4BestUnit(sigma2[j], "Surface/Mass");
168  }
169 
170  //print mean free path
171 
172  G4double lambda;
173 
174  G4cout << "\n \n mean free path : ";
175  for (size_t j=0; j<sigma1.size();j++) {
176  lambda = DBL_MAX;
177  if (sigma1[j] > 0.) lambda = 1/sigma1[j];
178  G4cout << "\t" << std::setw(13) << G4BestUnit( lambda, "Length");
179  }
180 
181  //mean free path (g/cm2)
182  G4cout << "\n (g/cm2) : ";
183  for (size_t j=0; j<sigma2.size();j++) {
184  lambda = DBL_MAX;
185  if (sigma2[j] > 0.) lambda = 1/sigma2[j];
186  G4cout << "\t" << std::setw(13) << G4BestUnit( lambda, "Mass/Surface");
187  }
188  G4cout << G4endl;
189 
190  if (charge == 0.) {
191  G4cout.precision(prec);
192  G4cout << "\n-------------------------------------------------------------\n"
193  << G4endl;
194  return;
195  }
196 
197  //get stopping power
198  std::vector<G4double> dedx1;
199  std::vector<G4double> dedx2;
200  G4double dedx, dedxtot = 0.;
201 
202  for (size_t j=0; j<emName.size();j++) {
203  dedx = emCal.ComputeDEDX(energy,particle,emName[j],material,enerCut[j]);
204  dedx1.push_back(dedx);
205  dedx2.push_back(dedx/density);
206  }
207  dedxtot = emCal.GetDEDX(energy,particle,material);
208  dedx1.push_back(dedxtot);
209  dedx2.push_back(dedxtot/density);
210 
211  //print stopping power
212  G4cout << "\n \n restricted dE/dx : ";
213  for (size_t j=0; j<sigma1.size();j++) {
214  G4cout << "\t" << std::setw(13) << G4BestUnit(dedx1[j],"Energy/Length");
215  }
216 
217  G4cout << "\n (MeV/g/cm2) : ";
218  for (size_t j=0; j<sigma2.size();j++) {
219  G4cout << "\t" << std::setw(13) << G4BestUnit(dedx2[j],"Energy*Surface/Mass");
220  }
221 
222  //get range from restricted dedx
223  G4double range1 = emCal.GetRangeFromRestricteDEDX(energy,particle,material);
224  G4double range2 = range1*density;
225 
226  //get range from full dedx
227  G4double Range1 = emCal.GetCSDARange(energy,particle,material);
228  G4double Range2 = Range1*density;
229 
230  //print range
231  G4cout << "\n \n range from restrict dE/dx: "
232  << "\t" << std::setw(8) << G4BestUnit(range1,"Length")
233  << " (" << std::setw(8) << G4BestUnit(range2,"Mass/Surface") << ")";
234 
235  G4cout << "\n range from full dE/dx : "
236  << "\t" << std::setw(8) << G4BestUnit(Range1,"Length")
237  << " (" << std::setw(8) << G4BestUnit(Range2,"Mass/Surface") << ")";
238 
239  //get transport mean free path (for multiple scattering)
240  G4double MSmfp1 = emCal.GetMeanFreePath(energy,particle,"msc",material);
241  G4double MSmfp2 = MSmfp1*density;
242 
243  //print transport mean free path
244  G4cout << "\n \n transport mean free path : "
245  << "\t" << std::setw(8) << G4BestUnit(MSmfp1,"Length")
246  << " (" << std::setw(8) << G4BestUnit(MSmfp2,"Mass/Surface") << ")";
247 
248  if (particle == G4Electron::Electron()) CriticalEnergy();
249 
250  G4cout << "\n-------------------------------------------------------------\n";
251  G4cout << G4endl;
252 
253  // reset default precision
254  G4cout.precision(prec);
255 }
256 
257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
258 
259 void RunAction::EndOfRunAction(const G4Run* )
260 { }
261 
262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
263 
264 #include "G4ProductionCutsTable.hh"
265 
266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
267 
269 {
270  G4ProductionCutsTable* theCoupleTable =
271  G4ProductionCutsTable::GetProductionCutsTable();
272 
273  size_t numOfCouples = theCoupleTable->GetTableSize();
274  const G4MaterialCutsCouple* couple = 0;
275  G4int index = 0;
276  for (size_t i=0; i<numOfCouples; i++) {
277  couple = theCoupleTable->GetMaterialCutsCouple(i);
278  if (couple->GetMaterial() == detector->GetMaterial()) {index = i; break;}
279  }
280 
281  rangeCut[0] =
282  (*(theCoupleTable->GetRangeCutsVector(idxG4GammaCut)))[index];
283  rangeCut[1] =
284  (*(theCoupleTable->GetRangeCutsVector(idxG4ElectronCut)))[index];
285  rangeCut[2] =
286  (*(theCoupleTable->GetRangeCutsVector(idxG4PositronCut)))[index];
287 
288  energyCut[0] =
289  (*(theCoupleTable->GetEnergyCutsVector(idxG4GammaCut)))[index];
290  energyCut[1] =
291  (*(theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut)))[index];
292  energyCut[2] =
293  (*(theCoupleTable->GetEnergyCutsVector(idxG4PositronCut)))[index];
294 
295 }
296 
297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
298 
300 {
301  // compute e- critical energy (Rossi definition) and Moliere radius.
302  // Review of Particle Physics - Eur. Phys. J. C3 (1998) page 147
303  //
304  G4EmCalculator emCal;
305 
306  const G4Material* material = detector->GetMaterial();
307  const G4double radl = material->GetRadlen();
308  G4double ekin = 5*MeV;
309  G4double deioni;
310  G4double err = 1., errmax = 0.001;
311  G4int iter = 0 , itermax = 10;
312  while (err > errmax && iter < itermax) {
313  iter++;
314  deioni = radl*
315  emCal.ComputeDEDX(ekin,G4Electron::Electron(),"eIoni",material);
316  err = std::abs(deioni - ekin)/ekin;
317  ekin = deioni;
318  }
319  G4cout << "\n \n critical energy (Rossi) : "
320  << "\t" << std::setw(8) << G4BestUnit(ekin,"Energy");
321 
322  //Pdg formula (only for single material)
323  G4double pdga[2] = { 610*MeV, 710*MeV };
324  G4double pdgb[2] = { 1.24, 0.92 };
325  G4double EcPdg = 0.;
326 
327  if (material->GetNumberOfElements() == 1) {
328  G4int istat = 0;
329  if (material->GetState() == kStateGas) istat = 1;
330  G4double Zeff = material->GetZ() + pdgb[istat];
331  EcPdg = pdga[istat]/Zeff;
332  G4cout << "\t\t\t (from Pdg formula : "
333  << std::setw(8) << G4BestUnit(EcPdg,"Energy") << ")";
334  }
335 
336  const G4double Es = 21.2052*MeV;
337  G4double rMolier1 = Es/ekin, rMolier2 = rMolier1*radl;
338  G4cout << "\n Moliere radius : "
339  << "\t" << std::setw(8) << rMolier1 << " X0 "
340  << "= " << std::setw(8) << G4BestUnit(rMolier2,"Length");
341 
342  if (material->GetNumberOfElements() == 1) {
343  G4double rMPdg = radl*Es/EcPdg;
344  G4cout << "\t (from Pdg formula : "
345  << std::setw(8) << G4BestUnit(rMPdg,"Length") << ")";
346  }
347 }
348 
349 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EndOfRunAction(const G4Run *)
Definition: RunAction.cc:259
DetectorConstruction * detector
Definition: RunAction.hh:58
void GetCuts()
Definition: RunAction.cc:268
constexpr auto abs(T v)
Returns the absolute value of the argument.
G4ParticleGun * GetParticleGun()
G4double energyCut[3]
Definition: RunAction.hh:61
void BeginOfRunAction(const G4Run *)
Definition: RunAction.cc:56
Float_t Z
Definition: plot.C:37
double energy
Definition: plottest35.C:25
void CriticalEnergy()
Definition: RunAction.cc:299
G4double rangeCut[3]
Definition: RunAction.hh:60
PrimaryGeneratorAction * primary
Definition: RunAction.hh:59
~RunAction()
Definition: RunAction.cc:51
RunAction(DetectorConstruction *, PrimaryGeneratorAction *)
Definition: RunAction.cc:45