source: trunk/source/processes/electromagnetic/xrays/test/testG4HeatedKleinNishina.cc@ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 9.4 KB
Line 
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// Unit test for coherent elastic models
28//
29// 18.05.07 V. Grichine
30//
31//
32
33#include "G4ios.hh"
34#include <fstream>
35#include <cmath>
36#include "globals.hh"
37#include "Randomize.hh"
38#include "G4UnitsTable.hh"
39#include <iomanip>
40#include <complex>
41
42
43#include "G4Element.hh"
44#include "G4NistManager.hh"
45
46
47#include "G4ParticleDefinition.hh"
48#include "G4PionPlus.hh"
49#include "G4PionMinus.hh"
50#include "G4MuonPlus.hh"
51#include "G4MuonMinus.hh"
52#include "G4Proton.hh"
53#include "G4Electron.hh"
54#include "G4Gamma.hh"
55
56#include "G4DynamicParticle.hh"
57#include "G4ParticleMomentum.hh"
58
59#include "G4VEmModel.hh"
60#include "G4KleinNishinaCompton.hh"
61#include "G4HeatedKleinNishinaCompton.hh"
62#include "G4PEEffectModel.hh"
63#include "G4PairProductionRelModel.hh"
64#include "G4BetheHeitlerModel.hh"
65
66
67using namespace std;
68
69
70
71int main()
72{
73
74 G4int i, j, k, iMax;
75 G4double x;
76 G4double expXrad=0., g4Xrad;
77
78 std::ofstream writef("angle.dat", std::ios::out ) ;
79 writef.setf( std::ios::scientific, std::ios::floatfield );
80
81
82 G4Element* theElement;
83 G4Material* theMaterial;
84 G4NistManager* man = G4NistManager::Instance();
85 man->SetVerbose(1);
86
87 G4cout << " 1 hydrogen" << G4endl;
88 G4cout << " 2 helium" << G4endl;
89 G4cout << " 4 berillium" << G4endl;
90 G4cout << " 6 carbon" << G4endl;
91 G4cout << " 7 nitrogen" << G4endl;
92 G4cout << " 8 oxigen" << G4endl;
93 G4cout << "13 aluminium" << G4endl;
94 G4cout << "14 silicon" << G4endl;
95 G4cout << "18 argon" << G4endl;
96 G4cout << "26 iron" << G4endl;
97 G4cout << "29 copper" << G4endl;
98 G4cout << "48 cadmium" << G4endl;
99 G4cout << "54 xenon" << G4endl;
100 G4cout << "74 tugnsten" << G4endl;
101 G4cout << "77 iridium" << G4endl;
102 G4cout << "82 lead" << G4endl;
103 G4cout << "92 uranium" << G4endl;
104 G4int choice;
105 // G4cin >> choice;
106 choice = 6;
107
108
109 switch (choice)
110 {
111 case 1:
112
113 theElement = man->FindOrBuildElement("H");
114 theMaterial = man->FindOrBuildMaterial("G4_H");
115 g4Xrad = theMaterial->GetRadlen();
116 break;
117
118 case 2:
119
120 theElement = man->FindOrBuildElement("He");
121 theMaterial = man->FindOrBuildMaterial("G4_He");
122 g4Xrad = theMaterial->GetRadlen();
123 break;
124
125 case 4:
126
127 theElement = man->FindOrBuildElement("Be");
128 theMaterial = man->FindOrBuildMaterial("G4_Be");
129 g4Xrad = theMaterial->GetRadlen();
130 break;
131
132 case 6:
133
134 theElement = man->FindOrBuildElement("C");
135 theMaterial = man->FindOrBuildMaterial("G4_C");
136 g4Xrad = theMaterial->GetRadlen();
137 expXrad = 19.6*cm;
138 break;
139
140 case 7:
141
142 theElement = man->FindOrBuildElement("N");
143 theMaterial = man->FindOrBuildMaterial("G4_N");
144 g4Xrad = theMaterial->GetRadlen();
145 break;
146
147
148 case 8:
149
150 theElement = man->FindOrBuildElement("O");
151 theMaterial = man->FindOrBuildMaterial("G4_O");
152 g4Xrad = theMaterial->GetRadlen();
153 break;
154
155 case 13:
156
157 theElement = man->FindOrBuildElement("Al");
158 theMaterial = man->FindOrBuildMaterial("G4_Al");
159 g4Xrad = theMaterial->GetRadlen();
160 expXrad = 8.9*cm;
161 break;
162
163 case 14:
164
165 theElement = man->FindOrBuildElement("Si");
166 theMaterial = man->FindOrBuildMaterial("G4_Si");
167 g4Xrad = theMaterial->GetRadlen();
168 break;
169
170 case 18:
171
172 theElement = man->FindOrBuildElement("Ar");
173 theMaterial = man->FindOrBuildMaterial("G4_Ar");
174 break;
175
176 case 26:
177
178 theElement = man->FindOrBuildElement("Fe");
179 theMaterial = man->FindOrBuildMaterial("G4_Fe");
180 g4Xrad = theMaterial->GetRadlen();
181 expXrad = 1.76*cm;
182 break;
183
184 case 29:
185
186 theElement = man->FindOrBuildElement("Cu");
187 theMaterial = man->FindOrBuildMaterial("G4_Cu");
188 g4Xrad = theMaterial->GetRadlen();
189 break;
190
191 case 48:
192
193 theElement = man->FindOrBuildElement("Cd");
194 theMaterial = man->FindOrBuildMaterial("G4_Cd");
195 g4Xrad = theMaterial->GetRadlen();
196 break;
197
198
199 case 54:
200
201 theElement = man->FindOrBuildElement("Xe");
202 theMaterial = man->FindOrBuildMaterial("G4_Xe");
203 g4Xrad = theMaterial->GetRadlen();
204 break;
205
206
207 case 74:
208
209 theElement = man->FindOrBuildElement("W");
210 theMaterial = man->FindOrBuildMaterial("G4_W");
211 g4Xrad = theMaterial->GetRadlen();
212 expXrad = 0.35*cm;
213 break;
214
215 case 77:
216
217 theElement = man->FindOrBuildElement("Ir");
218 theMaterial = man->FindOrBuildMaterial("G4_Ir");
219 g4Xrad = theMaterial->GetRadlen();
220 break;
221
222 case 82:
223
224 theElement = man->FindOrBuildElement("Pb");
225 theMaterial = man->FindOrBuildMaterial("G4_Pb");
226 g4Xrad = theMaterial->GetRadlen();
227 expXrad = 0.56*cm;
228 break;
229
230 case 92:
231
232 theElement = man->FindOrBuildElement("U");
233 theMaterial = man->FindOrBuildMaterial("G4_U");
234 g4Xrad = theMaterial->GetRadlen();
235 expXrad = 0.35*cm;
236 break;
237 }
238
239// Particle definition
240
241
242 G4cout << " 0 gamma" << G4endl;
243 G4cout << " 1 electron" << G4endl;
244 G4cout << " 2 proton" << G4endl;
245 G4cout << " 3 pion+" << G4endl;
246 G4cout << " 4 pion-" << G4endl;
247 G4cout << " 4 muon+" << G4endl;
248 G4cout << " 5 muon-" << G4endl;
249
250 // G4cin >> choice;
251 choice = 0;
252
253 G4ParticleDefinition* theParticleDefinition;
254
255
256 switch (choice)
257 {
258 case 0:
259
260 theParticleDefinition = G4Gamma::GammaDefinition();
261 break;
262
263 case 1:
264
265 theParticleDefinition = G4Electron::ElectronDefinition();
266 break;
267
268 case 2:
269
270 theParticleDefinition = G4Proton::ProtonDefinition();
271 break;
272
273 case 3:
274
275 theParticleDefinition = G4PionPlus::PionPlusDefinition();
276 break;
277
278 case 4:
279
280 theParticleDefinition = G4PionMinus::PionMinusDefinition();
281 break;
282
283 case 5:
284
285 theParticleDefinition = G4MuonPlus::MuonPlusDefinition();
286 break;
287
288 case 6:
289
290 theParticleDefinition = G4MuonMinus::MuonMinusDefinition();
291 break;
292
293 }
294
295 G4double energyMscXR, xsc, kinEnergy;
296
297
298 // kinEnergy = 8*GeV; // 25.0*GeV;
299
300 kinEnergy = 25.0*GeV;
301
302 G4VEmModel* comp = new G4KleinNishinaCompton();
303 // G4VEmModel* comp = new G4HeatedKleinNishinaCompton();
304 G4VEmModel* photo = new G4PEEffectModel();
305 G4VEmModel* pair = new G4PairProductionRelModel(theParticleDefinition,"pp");
306 G4VEmModel* bhpair = new G4BetheHeitlerModel(theParticleDefinition,"bhpp");
307
308 G4DynamicParticle* theDynamicParticle = new G4DynamicParticle(theParticleDefinition,
309 G4ParticleMomentum(0.,0.,1.),
310 kinEnergy);
311
312 G4double m1 = theParticleDefinition->GetPDGMass();
313 G4double plab = theDynamicParticle->GetTotalMomentum();
314 G4cout <<"lab momentum, plab = "<<plab/GeV<<" GeV"<<G4endl;
315 G4double plabLowLimit = 20.0*MeV;
316
317 G4int Z = G4int(theElement->GetZ());
318 G4int A = G4int(theElement->GetN()+0.5);
319
320 G4double step = 4.10*mm;
321
322 step = expXrad;
323
324 G4double m2 = man->GetAtomicMassAmu(Z)*GeV;
325 // G4double m2 = man->GetAtomicMass( Z, A);
326 G4cout <<" target mass, m2 = "<<m2/GeV<<" GeV"<<G4endl<<G4endl;
327 G4cout <<"step = "<<step<<" mm; g4Xrad = "<<g4Xrad<<" mm; expXrad = "
328 <<expXrad<<" mm"<<G4endl<<G4endl;
329
330
331 // energyMscXR = 2*MeV;
332 //xsc = pair->ComputeCrossSectionPerAtom(theParticleDefinition, energyMscXR, G4double(Z ) );
333
334
335 iMax = 130;
336
337 writef<<iMax<<G4endl;
338
339
340
341
342 for( i = 0; i < iMax; i++ )
343 {
344 energyMscXR = std::exp(i*0.2)*0.001*keV;
345 // xsc = comp->ComputeCrossSectionPerAtom(theParticleDefinition,energyMscXR,G4double(Z ) );
346 // xsc = photo->ComputeCrossSectionPerAtom(theParticleDefinition, energyMscXR, G4double(Z ) );
347 // xsc = pair->ComputeCrossSectionPerAtom(theParticleDefinition, energyMscXR, G4double(Z ) );
348 xsc = bhpair->ComputeCrossSectionPerAtom(theParticleDefinition, energyMscXR, G4double(Z ) );
349 G4cout<<energyMscXR/MeV<<"\t\t"<<xsc/barn<<G4endl;
350 writef<<energyMscXR/MeV<<"\t\t"<<xsc/barn<<G4endl;
351 }
352
353
354
355
356
357 return 1;
358} // end of main
Note: See TracBrowser for help on using the repository browser.