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

Last change on this file since 1346 was 1199, checked in by garnier, 16 years ago

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