LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
DirectAccess.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 //
28 //
29 //
30 //
31 // ------------------------------------------------------------
32 //
33 // To print cross sections per atom and mean free path for simple material
34 //
35 #include "G4Material.hh"
36 
37 #include "G4PEEffectFluoModel.hh"
38 #include "G4KleinNishinaCompton.hh"
39 #include "G4BetheHeitlerModel.hh"
40 
41 #include "G4eeToTwoGammaModel.hh"
42 
43 #include "G4MollerBhabhaModel.hh"
44 #include "G4SeltzerBergerModel.hh"
45 
46 #include "G4BetheBlochModel.hh"
47 #include "G4BraggModel.hh"
48 
49 #include "G4MuBetheBlochModel.hh"
50 #include "G4MuBremsstrahlungModel.hh"
51 #include "G4MuPairProductionModel.hh"
52 
53 #include "globals.hh"
54 #include "G4UnitsTable.hh"
55 #include "G4SystemOfUnits.hh"
56 
57 #include "G4Gamma.hh"
58 #include "G4Positron.hh"
59 #include "G4Electron.hh"
60 #include "G4Proton.hh"
61 #include "G4MuonPlus.hh"
62 
63 #include "G4DataVector.hh"
64 #include "G4NistManager.hh"
65 #include "G4ParticleTable.hh"
66 
67 int main() {
68 
69  G4UnitDefinition::BuildUnitsTable();
70 
71  G4ParticleDefinition* gamma = G4Gamma::Gamma();
72  G4ParticleDefinition* posit = G4Positron::Positron();
73  G4ParticleDefinition* elec = G4Electron::Electron();
74  G4ParticleDefinition* prot = G4Proton::Proton();
75  G4ParticleDefinition* muon = G4MuonPlus::MuonPlus();
76  G4ParticleTable* partTable = G4ParticleTable::GetParticleTable();
77  partTable->SetReadiness();
78 
79  G4DataVector cuts;
80  cuts.push_back(1*keV);
81 
82  // define materials
83  //
84  G4Material* material =
85  G4NistManager::Instance()->FindOrBuildMaterial("G4_Fe");
86 
87  G4cout << *(G4Material::GetMaterialTable()) << G4endl;
88 
89  G4MaterialCutsCouple* couple = new G4MaterialCutsCouple(material);
90  couple->SetIndex(0);
91 
92  // work only for simple materials
93  G4double Z = material->GetZ();
94  G4double A = material->GetA();
95 
96  // initialise gamma processes (models)
97  //
98  G4VEmModel* phot = new G4PEEffectFluoModel();
99  G4VEmModel* comp = new G4KleinNishinaCompton();
100  G4VEmModel* conv = new G4BetheHeitlerModel();
101  phot->Initialise(gamma, cuts);
102  comp->Initialise(gamma, cuts);
103  conv->Initialise(gamma, cuts);
104 
105  // valid pointer to a couple is needed for this model
106  phot->SetCurrentCouple(couple);
107 
108  // compute CrossSection per atom and MeanFreePath
109  //
110  G4double Emin = 1.01*MeV, Emax = 2.01*MeV, dE = 100*keV;
111 
112  G4cout << "\n #### Gamma : CrossSectionPerAtom and MeanFreePath for "
113  << material->GetName() << G4endl;
114  G4cout << "\n Energy \t PhotoElec \t Compton \t Conversion \t";
115  G4cout << "\t PhotoElec \t Compton \t Conversion" << G4endl;
116 
117  for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
118  G4cout << "\n " << G4BestUnit (Energy, "Energy")
119  << "\t"
120  << G4BestUnit (phot->ComputeCrossSectionPerAtom(gamma,Energy,Z),"Surface")
121  << "\t"
122  << G4BestUnit (comp->ComputeCrossSectionPerAtom(gamma,Energy,Z),"Surface")
123  << "\t"
124  << G4BestUnit (conv->ComputeCrossSectionPerAtom(gamma,Energy,Z),"Surface")
125  << "\t \t"
126  << G4BestUnit (phot->ComputeMeanFreePath(gamma,Energy,material),"Length")
127  << "\t"
128  << G4BestUnit (comp->ComputeMeanFreePath(gamma,Energy,material),"Length")
129  << "\t"
130  << G4BestUnit (conv->ComputeMeanFreePath(gamma,Energy,material),"Length");
131  }
132 
133  G4cout << G4endl;
134 
135  // initialise positron annihilation (model)
136  //
137  G4VEmModel* anni = new G4eeToTwoGammaModel();
138  anni->Initialise(posit, cuts);
139 
140  // compute CrossSection per atom and MeanFreePath
141  //
142  Emin = 1.01*MeV; Emax = 2.01*MeV; dE = 100*keV;
143 
144  G4cout << "\n #### e+ annihilation : CrossSectionPerAtom and MeanFreePath"
145  << " for " << material->GetName() << G4endl;
146  G4cout << "\n Energy \t e+ annihil \t";
147  G4cout << "\t e+ annihil" << G4endl;
148 
149  for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
150  G4cout << "\n " << G4BestUnit (Energy, "Energy")
151  << "\t"
152  << G4BestUnit (anni->ComputeCrossSectionPerAtom(posit,Energy,Z),"Surface")
153  << "\t \t"
154  << G4BestUnit (anni->ComputeMeanFreePath(posit,Energy,material),"Length");
155  }
156 
157  G4cout << G4endl;
158 
159  // initialise electron processes (models)
160  //
161  G4VEmModel* ioni = new G4MollerBhabhaModel();
162  G4VEmModel* brem = new G4SeltzerBergerModel();
163  ioni->Initialise(elec, cuts);
164  brem->Initialise(elec, cuts);
165 
166  // compute CrossSection per atom and MeanFreePath
167  //
168  Emin = 1.01*MeV; Emax = 101.01*MeV; dE = 10*MeV;
169  G4double Ecut = 100*keV;
170 
171  G4cout << "\n ####electron: CrossSection, MeanFreePath and StoppingPower"
172  << " for " << material->GetName()
173  << ";\tEnergy cut = " << G4BestUnit (Ecut, "Energy") << G4endl;
174 
175  G4cout << "\n Energy \t ionization \t bremsstra \t";
176  G4cout << "\t ionization \t bremsstra \t";
177  G4cout << "\t ionization \t bremsstra" << G4endl;
178 
179  for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
180  G4cout << "\n " << G4BestUnit (Energy, "Energy")
181  << "\t"
182  << G4BestUnit (ioni->ComputeCrossSectionPerAtom(elec,Energy,Z,A,Ecut),
183  "Surface")
184  << "\t"
185  << G4BestUnit (brem->ComputeCrossSectionPerAtom(elec,Energy,Z,A,Ecut),
186  "Surface")
187  << "\t \t"
188  << G4BestUnit (ioni->ComputeMeanFreePath(elec,Energy,material,Ecut),
189  "Length")
190  << "\t"
191  << G4BestUnit (brem->ComputeMeanFreePath(elec,Energy,material,Ecut),
192  "Length")
193  << "\t \t"
194  << G4BestUnit (ioni->ComputeDEDXPerVolume(material,elec,Energy,Ecut),
195  "Energy/Length")
196  << "\t"
197  << G4BestUnit (brem->ComputeDEDXPerVolume(material,elec,Energy,Ecut),
198  "Energy/Length");
199  }
200 
201  G4cout << G4endl;
202 
203  // initialise proton processes (models)
204  //
205  ioni = new G4BetheBlochModel();
206  ioni->Initialise(prot, cuts);
207 
208  // compute CrossSection per atom and MeanFreePath
209  //
210  Emin = 1.01*MeV; Emax = 102.01*MeV; dE = 10*MeV;
211  Ecut = 100*keV;
212 
213  G4cout << "\n #### proton : CrossSection, MeanFreePath and StoppingPower"
214  << " for " << material->GetName()
215  << ";\tEnergy cut = " << G4BestUnit (Ecut, "Energy") << G4endl;
216 
217  G4cout << "\n Energy \t ionization \t";
218  G4cout << "\t ionization \t";
219  G4cout << "\t ionization" << G4endl;
220 
221  for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
222  G4cout << "\n " << G4BestUnit (Energy, "Energy")
223  << "\t"
224  << G4BestUnit (ioni->ComputeCrossSectionPerAtom(prot,Energy,Z,A,Ecut),
225  "Surface")
226  << "\t \t"
227  << G4BestUnit (ioni->ComputeMeanFreePath(prot,Energy,material,Ecut),
228  "Length")
229  << "\t \t"
230  << G4BestUnit (ioni->ComputeDEDXPerVolume(material,prot,Energy,Ecut),
231  "Energy/Length");
232  }
233 
234  G4cout << G4endl;
235 
236  // low energy : Bragg Model
237  ioni = new G4BraggModel(prot);
238  ioni->Initialise(prot, cuts);
239 
240  // compute CrossSection per atom and MeanFreePath
241  //
242  Emin = 1.1*keV; Emax = 2.01*MeV; dE = 300*keV;
243  Ecut = 10*keV;
244 
245  G4cout << "\n #### proton : low energy model (Bragg) "
246  << ";\tEnergy cut = " << G4BestUnit (Ecut, "Energy") << G4endl;
247 
248  G4cout << "\n Energy \t ionization \t";
249  G4cout << "\t ionization \t";
250  G4cout << "\t ionization" << G4endl;
251 
252  for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
253  G4cout << "\n " << G4BestUnit (Energy, "Energy")
254  << "\t"
255  << G4BestUnit (ioni->ComputeCrossSectionPerAtom(prot,Energy,Z,A,Ecut),
256  "Surface")
257  << "\t \t"
258  << G4BestUnit (ioni->ComputeMeanFreePath(prot,Energy,material,Ecut),
259  "Length")
260  << "\t \t"
261  << G4BestUnit (ioni->ComputeDEDXPerVolume(material,prot,Energy,Ecut),
262  "Energy/Length");
263  }
264 
265  G4cout << G4endl;
266 
267  // initialise muon processes (models)
268  //
269  ioni = new G4MuBetheBlochModel();
270  brem = new G4MuBremsstrahlungModel();
271  G4VEmModel* pair = new G4MuPairProductionModel();
272  ioni->Initialise(muon, cuts);
273  brem->Initialise(muon, cuts);
274  pair->Initialise(muon, cuts);
275 
276  // compute CrossSection per atom and MeanFreePath
277  //
278  Emin = 1.01*GeV; Emax = 101.01*GeV; dE = 10*GeV;
279  Ecut = 10*MeV;
280 
281  G4cout << "\n ####muon: CrossSection and MeanFreePath for "
282  << material->GetName()
283  << ";\tEnergy cut = " << G4BestUnit (Ecut, "Energy") << G4endl;
284 
285  G4cout << "\n Energy \t ionization \t bremsstra \t pair_prod \t";
286  G4cout << "\t ionization \t bremsstra \t pair_prod" << G4endl;
287 
288  for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
289  G4cout << "\n " << G4BestUnit (Energy, "Energy")
290  << "\t"
291  << G4BestUnit (ioni->ComputeCrossSectionPerAtom(muon,Energy,Z,A,Ecut),
292  "Surface")
293  << "\t"
294  << G4BestUnit (brem->ComputeCrossSectionPerAtom(muon,Energy,Z,A,Ecut),
295  "Surface")
296  << "\t"
297  << G4BestUnit (pair->ComputeCrossSectionPerAtom(muon,Energy,Z,A,Ecut),
298  "Surface")
299  << "\t \t"
300  << G4BestUnit (ioni->ComputeMeanFreePath(muon,Energy,material,Ecut),
301  "Length")
302  << "\t"
303  << G4BestUnit (brem->ComputeMeanFreePath(muon,Energy,material,Ecut),
304  "Length")
305  << "\t"
306  << G4BestUnit (pair->ComputeMeanFreePath(muon,Energy,material,Ecut),
307  "Length");
308  }
309 
310  G4cout << G4endl;
311 
312  G4cout << "\n ####muon: StoppingPower for "
313  << material->GetName()
314  << ";\tEnergy cut = " << G4BestUnit (Ecut, "Energy") << G4endl;
315 
316  G4cout << "\n Energy \t ionization \t bremsstra \t pair_prod \t" << G4endl;
317 
318  for (G4double Energy = Emin; Energy <= Emax; Energy += dE) {
319  G4cout << "\n " << G4BestUnit (Energy, "Energy")
320  << "\t"
321  << G4BestUnit (ioni->ComputeDEDXPerVolume(material,muon,Energy,Ecut),
322  "Energy/Length")
323  << "\t"
324  << G4BestUnit (brem->ComputeDEDXPerVolume(material,muon,Energy,Ecut),
325  "Energy/Length")
326  << "\t"
327  << G4BestUnit (pair->ComputeDEDXPerVolume(material,muon,Energy,Ecut),
328  "Energy/Length");
329  }
330 
331  G4cout << G4endl;
332  return EXIT_SUCCESS;
333 }
334 
335 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
main
int main()
Definition: DirectAccess.cc:67
Z
Float_t Z
Definition: plot.C:37
phot
General LArSoft Utilities.
Definition: BeamFlashTrackMatchTaggerAlg.h:26