ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RunAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RunAction.cc
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"
32 #include "PrimaryGeneratorAction.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 
45 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin)
46 :detector(det), primary(kin)
47 { }
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 
51 RunAction::~RunAction()
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 
268 void RunAction::GetCuts()
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 
299 void RunAction::CriticalEnergy()
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......