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