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 //
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 
33 #include "RunAction.hh"
34 #include "DetectorConstruction.hh"
35 #include "PrimaryGeneratorAction.hh"
36 
37 #include "G4Run.hh"
38 #include "G4ProcessManager.hh"
39 #include "G4LossTableManager.hh"
40 #include "G4UnitsTable.hh"
41 #include "G4EmCalculator.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4Electron.hh"
45 #include "G4Positron.hh"
46 
47 #include <vector>
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 
51 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin)
52 :G4UserRunAction(),fDetector(det), fPrimary(kin)
53 { }
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 
57 RunAction::~RunAction()
58 { }
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
62 void RunAction::BeginOfRunAction(const G4Run*)
63 {
64  //set precision for printing
65  G4int prec = G4cout.precision(6);
66 
67  //instanciate EmCalculator
68  G4EmCalculator emCal;
69  // emCal.SetVerbose(2);
70 
71  // get particle
72  G4ParticleDefinition* particle = fPrimary->GetParticleGun()
73  ->GetParticleDefinition();
74  G4String partName = particle->GetParticleName();
75  G4double charge = particle->GetPDGCharge();
76  G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
77 
78  // get material
79  const G4Material* material = fDetector->GetMaterial();
80  G4String matName = material->GetName();
81  G4double density = material->GetDensity();
82  G4double radl = material->GetRadlen();
83 
84  G4cout << "\n " << partName << " ("
85  << G4BestUnit(energy,"Energy") << ") in "
86  << material->GetName() << " (density: "
87  << G4BestUnit(density,"Volumic Mass") << "; radiation length: "
88  << G4BestUnit(radl, "Length") << ")" << G4endl;
89 
90  // get cuts
91  GetCuts();
92  if (charge != 0.) {
93  G4cout << "\n Range cuts : \t gamma "
94  << std::setw(8) << G4BestUnit(fRangeCut[0],"Length")
95  << "\t e- " << std::setw(8) << G4BestUnit(fRangeCut[1],"Length");
96  G4cout << "\n Energy cuts : \t gamma "
97  << std::setw(8) << G4BestUnit(fEnergyCut[0],"Energy")
98  << "\t e- " << std::setw(8) << G4BestUnit(fEnergyCut[1],"Energy")
99  << G4endl;
100  }
101 
102  // max energy transfert
103  if (charge != 0.) {
104  G4double Mass_c2 = particle->GetPDGMass();
105  G4double moverM = electron_mass_c2/Mass_c2;
106  G4double gamM1 = energy/Mass_c2, gam = gamM1 + 1., gamP1 = gam + 1.;
107  G4double Tmax = energy;
108  if(particle == G4Electron::Electron()) {
109  Tmax *= 0.5;
110  } else if(particle != G4Positron::Positron()) {
111  Tmax = (2*electron_mass_c2*gamM1*gamP1)/(1.+2*gam*moverM+moverM*moverM);
112  }
113  G4double range = emCal.GetCSDARange(Tmax,G4Electron::Electron(),material);
114 
115  G4cout << "\n Max_energy _transferable : " << G4BestUnit(Tmax,"Energy")
116  << " (" << G4BestUnit(range,"Length") << ")" << G4endl;
117  }
118 
119  // get processList and extract EM processes (but not MultipleScattering)
120  G4ProcessVector* plist = particle->GetProcessManager()->GetProcessList();
121  G4String procName;
122  G4double cut;
123  std::vector<G4String> emName;
124  std::vector<G4double> enerCut;
125  size_t length = plist->size();
126  for (size_t j=0; j<length; j++) {
127  procName = (*plist)[j]->GetProcessName();
128  cut = fEnergyCut[1];
129  if ((procName == "eBrem")||(procName == "muBrems")) cut = fEnergyCut[0];
130  if (((*plist)[j]->GetProcessType() == fElectromagnetic) &&
131  (procName != "msc")) {
132  emName.push_back(procName);
133  enerCut.push_back(cut);
134  }
135  }
136 
137  // write html documentation, if requested
138  char* htmlDocName = std::getenv("G4PhysListName"); // file name
139  char* htmlDocDir = std::getenv("G4PhysListDocDir"); // directory
140  if (htmlDocName && htmlDocDir) {
141  G4LossTableManager::Instance()->DumpHtml();
142  }
143 
144  // print list of processes
145  G4cout << "\n processes : ";
146  for (size_t j=0; j<emName.size();j++)
147  G4cout << "\t" << std::setw(13) << emName[j] << "\t";
148  G4cout << "\t" << std::setw(13) <<"total";
149 
150  //compute cross section per atom (only for single material)
151  if (material->GetNumberOfElements() == 1) {
152  G4double Z = material->GetZ();
153  G4double A = material->GetA();
154 
155  std::vector<G4double> sigma0;
156  G4double sig, sigtot = 0.;
157 
158  for (size_t j=0; j<emName.size();j++) {
159  sig = emCal.ComputeCrossSectionPerAtom
160  (energy,particle,emName[j],Z,A,enerCut[j]);
161  sigtot += sig;
162  sigma0.push_back(sig);
163  }
164  sigma0.push_back(sigtot);
165 
166  G4cout << "\n \n cross section per atom : ";
167  for (size_t j=0; j<sigma0.size();j++) {
168  G4cout << "\t" << std::setw(13) << G4BestUnit(sigma0[j], "Surface");
169  }
170  G4cout << G4endl;
171  }
172 
173  //get cross section per volume
174  std::vector<G4double> sigma0;
175  std::vector<G4double> sigma1;
176  std::vector<G4double> sigma2;
177  G4double Sig, SigtotComp = 0., Sigtot = 0.;
178 
179  for (size_t j=0; j<emName.size();j++) {
180  Sig = emCal.ComputeCrossSectionPerVolume
181  (energy,particle,emName[j],material,enerCut[j]);
182  SigtotComp += Sig;
183  sigma0.push_back(Sig);
184  Sig = emCal.GetCrossSectionPerVolume(energy,particle,emName[j],material);
185  Sigtot += Sig;
186  sigma1.push_back(Sig);
187  sigma2.push_back(Sig/density);
188  }
189  sigma0.push_back(SigtotComp);
190  sigma1.push_back(Sigtot);
191  sigma2.push_back(Sigtot/density);
192 
193  //print cross sections
194  G4cout << "\n \n compCrossSectionPerVolume : ";
195  for (size_t j=0; j<sigma0.size();j++) {
196  G4cout << "\t" << std::setw(13) << sigma0[j]*cm << " cm^-1";
197  }
198  G4cout << "\n cross section per volume : ";
199  for (size_t j=0; j<sigma1.size();j++) {
200  G4cout << "\t" << std::setw(13) << sigma1[j]*cm << " cm^-1";
201  }
202 
203  G4cout << "\n cross section per mass : ";
204  for (size_t j=0; j<sigma2.size();j++) {
205  G4cout << "\t" << std::setw(13)
206  << G4BestUnit(sigma2[j], "Surface/Mass");
207  }
208 
209  //print mean free path
210 
211  G4double lambda;
212 
213  G4cout << "\n \n mean free path : ";
214  for (size_t j=0; j<sigma1.size();j++) {
215  lambda = DBL_MAX;
216  if (sigma1[j] > 0.) lambda = 1/sigma1[j];
217  G4cout << "\t" << std::setw(13) << G4BestUnit( lambda, "Length");
218  }
219 
220  //mean free path (g/cm2)
221  G4cout << "\n (g/cm2) : ";
222  for (size_t j=0; j<sigma2.size();j++) {
223  lambda = DBL_MAX;
224  if (sigma2[j] > 0.) lambda = 1/sigma2[j];
225  G4cout << "\t" << std::setw(13) << G4BestUnit( lambda, "Mass/Surface");
226  }
227  G4cout << G4endl;
228 
229  if (charge == 0.) {
230  G4cout.precision(prec);
231  G4cout << "\n-----------------------------------------------------------\n"
232  << G4endl;
233  return;
234  }
235 
236  //get stopping power
237  std::vector<G4double> dedx1;
238  std::vector<G4double> dedx2;
239  G4double dedx, dedxtot = 0.;
240  size_t nproc = emName.size();
241 
242  for (size_t j=0; j<nproc; j++) {
243  dedx = emCal.ComputeDEDX(energy,particle,emName[j],material,enerCut[j]);
244  dedxtot += dedx;
245  dedx1.push_back(dedx);
246  dedx2.push_back(dedx/density);
247  }
248  dedx1.push_back(dedxtot);
249  dedx2.push_back(dedxtot/density);
250 
251  //print stopping power
252  G4cout << "\n \n restricted dE/dx : ";
253  for (size_t j=0; j<=nproc; j++) {
254  G4cout << "\t" << std::setw(13)
255  << G4BestUnit(dedx1[j],"Energy/Length");
256  }
257 
258  G4cout << "\n (MeV/g/cm2) : ";
259  for (size_t j=0; j<=nproc; j++) {
260  G4cout << "\t" << std::setw(13)
261  << G4BestUnit(dedx2[j],"Energy*Surface/Mass");
262  }
263  dedxtot = 0.;
264 
265  for (size_t j=0; j<nproc; j++) {
266  dedx = emCal.ComputeDEDX(energy,particle,emName[j],material,energy);
267  dedxtot += dedx;
268  dedx1[j] = dedx;
269  dedx2[j] = dedx/density;
270  }
271  dedx1[nproc] = dedxtot;
272  dedx2[nproc] = dedxtot/density;
273 
274  //print stopping power
275  G4cout << "\n \n unrestricted dE/dx : ";
276  for (size_t j=0; j<=nproc; j++) {
277  G4cout << "\t" << std::setw(13) << G4BestUnit(dedx1[j],"Energy/Length");
278  }
279 
280  G4cout << "\n (MeV/g/cm2) : ";
281  for (size_t j=0; j<=nproc; j++) {
282  G4cout << "\t" << std::setw(13)
283  << G4BestUnit(dedx2[j],"Energy*Surface/Mass");
284  }
285 
286  //get range from restricted dedx
287  G4double range1 = emCal.GetRangeFromRestricteDEDX(energy,particle,material);
288  G4double range2 = range1*density;
289 
290  //print range
291  G4cout << "\n \n range from restrict dE/dx: "
292  << "\t" << std::setw(8) << G4BestUnit(range1,"Length")
293  << " (" << std::setw(8) << G4BestUnit(range2,"Mass/Surface") << ")";
294 
295  //get range from full dedx
296  G4double EmaxTable = G4EmParameters::Instance()->MaxEnergyForCSDARange();
297  if(energy < EmaxTable) {
298  G4double Range1 = emCal.GetCSDARange(energy,particle,material);
299  G4double Range2 = Range1*density;
300 
301  G4cout << "\n range from full dE/dx : "
302  << "\t" << std::setw(8) << G4BestUnit(Range1,"Length")
303  << " (" << std::setw(8) << G4BestUnit(Range2,"Mass/Surface") << ")";
304  }
305 
306  //get transport mean free path (for multiple scattering)
307  G4double MSmfp1 = emCal.GetMeanFreePath(energy,particle,"msc",material);
308  G4double MSmfp2 = MSmfp1*density;
309 
310  //print transport mean free path
311  G4cout << "\n \n transport mean free path : "
312  << "\t" << std::setw(8) << G4BestUnit(MSmfp1,"Length")
313  << " (" << std::setw(8) << G4BestUnit(MSmfp2,"Mass/Surface") << ")";
314 
315  if (particle == G4Electron::Electron()) CriticalEnergy();
316 
317  G4cout << "\n-------------------------------------------------------------\n";
318  G4cout << G4endl;
319 
320  // reset default precision
321  G4cout.precision(prec);
322 }
323 
324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
325 
326 void RunAction::EndOfRunAction(const G4Run* )
327 { }
328 
329 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
330 
331 #include "G4ProductionCutsTable.hh"
332 
333 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
334 
335 void RunAction::GetCuts()
336 {
337  G4ProductionCutsTable* theCoupleTable =
338  G4ProductionCutsTable::GetProductionCutsTable();
339 
340  size_t numOfCouples = theCoupleTable->GetTableSize();
341  const G4MaterialCutsCouple* couple = 0;
342  G4int index = 0;
343  for (size_t i=0; i<numOfCouples; i++) {
344  couple = theCoupleTable->GetMaterialCutsCouple(i);
345  if (couple->GetMaterial() == fDetector->GetMaterial()) {index = i; break;}
346  }
347 
348  fRangeCut[0] =
349  (*(theCoupleTable->GetRangeCutsVector(idxG4GammaCut)))[index];
350  fRangeCut[1] =
351  (*(theCoupleTable->GetRangeCutsVector(idxG4ElectronCut)))[index];
352  fRangeCut[2] =
353  (*(theCoupleTable->GetRangeCutsVector(idxG4PositronCut)))[index];
354 
355  fEnergyCut[0] =
356  (*(theCoupleTable->GetEnergyCutsVector(idxG4GammaCut)))[index];
357  fEnergyCut[1] =
358  (*(theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut)))[index];
359  fEnergyCut[2] =
360  (*(theCoupleTable->GetEnergyCutsVector(idxG4PositronCut)))[index];
361 
362 }
363 
364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
365 
366 void RunAction::CriticalEnergy()
367 {
368  // compute e- critical energy (Rossi definition) and Moliere radius.
369  // Review of Particle Physics - Eur. Phys. J. C3 (1998) page 147
370  //
371  G4EmCalculator emCal;
372 
373  const G4Material* material = fDetector->GetMaterial();
374  const G4double radl = material->GetRadlen();
375  G4double ekin = 5*MeV;
376  G4double deioni;
377  G4double err = 1., errmax = 0.001;
378  G4int iter = 0 , itermax = 10;
379  while (err > errmax && iter < itermax) {
380  iter++;
381  deioni = radl*
382  emCal.ComputeDEDX(ekin,G4Electron::Electron(),"eIoni",material);
383  err = std::abs(deioni - ekin)/ekin;
384  ekin = deioni;
385  }
386  G4cout << "\n \n critical energy (Rossi) : "
387  << "\t" << std::setw(8) << G4BestUnit(ekin,"Energy");
388 
389  //Pdg formula (only for single material)
390  G4double pdga[2] = { 610*MeV, 710*MeV };
391  G4double pdgb[2] = { 1.24, 0.92 };
392  G4double EcPdg = 0.;
393 
394  if (material->GetNumberOfElements() == 1) {
395  G4int istat = 0;
396  if (material->GetState() == kStateGas) istat = 1;
397  G4double Zeff = material->GetZ() + pdgb[istat];
398  EcPdg = pdga[istat]/Zeff;
399  G4cout << "\t\t\t (from Pdg formula : "
400  << std::setw(8) << G4BestUnit(EcPdg,"Energy") << ")";
401  }
402 
403  const G4double Es = 21.2052*MeV;
404  G4double rMolier1 = Es/ekin, rMolier2 = rMolier1*radl;
405  G4cout << "\n Moliere radius : "
406  << "\t" << std::setw(8) << rMolier1 << " X0 "
407  << "= " << std::setw(8) << G4BestUnit(rMolier2,"Length");
408 
409  if (material->GetNumberOfElements() == 1) {
410  G4double rMPdg = radl*Es/EcPdg;
411  G4cout << "\t (from Pdg formula : "
412  << std::setw(8) << G4BestUnit(rMPdg,"Length") << ")";
413  }
414 }
415 
416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......