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

Last change on this file since 1347 was 1347, checked in by garnier, 13 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.