source: trunk/source/processes/electromagnetic/xrays/test/testG4MscRadiation.cc @ 1244

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

nvx fichiers dans CVS
File size: 12.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
55#include "G4DynamicParticle.hh"
56#include "G4ParticleMomentum.hh"
57
58
59#include "G4MscRadiation.hh"
60
61
62using namespace std;
63
64
65G4double SimpleMF(G4double x)
66{
67  G4double order = 6.*x;
68  return 1. - std::exp(-order);
69}
70
71
72G4double SimpleMPsi(G4double x)
73{
74  G4double order = 4.*x;
75  return 1. - std::exp(-order);
76}
77
78
79G4double SimpleMG(G4double x)
80{ 
81  return 3.*SimpleMPsi(x) - 2.*SimpleMF(x);
82}
83
84
85
86G4double MediumMF(G4double x)
87{
88  G4double order = 6.*x;
89  order *= 1. + (3. - pi)*x;
90
91  return 1. - std::exp(-order);
92}
93
94
95G4double MediumMPsi(G4double x)
96{
97  G4double order = 4.*x + 8.*x*x;
98  return 1. - std::exp(-order);
99}
100
101
102G4double MediumMG(G4double x)
103{ 
104  return 3.*MediumMPsi(x) - 2.*MediumMF(x);
105}
106
107
108G4double ComplexMF(G4double x)
109{
110  G4double order = 6.*x;
111  order *= 1. + (3. - pi)*x;
112  order -= x*x*x/(0.623+0.796*x+0.658*x*x);
113
114  return 1. - std::exp(-order);
115}
116
117
118G4double ComplexMPsi(G4double x)
119{
120  G4double order = 4.*x + 8.*x*x/(1.+ 3.96*x + 4.97*x*x - 0.05*x*x*x + 7.5*x*x*x*x);
121  return 1. - std::exp(-order);
122}
123
124
125
126G4double ComplexMG(G4double x)
127{ 
128  return 3.*ComplexMPsi(x) - 2.*ComplexMF(x);
129}
130
131
132
133int main()
134{
135
136  G4int i, j, k, iMax;
137  G4double x;
138  G4double expXrad=0., g4Xrad;
139
140  std::ofstream writef("angle.dat", std::ios::out ) ;
141  writef.setf( std::ios::scientific, std::ios::floatfield );
142
143  // iMax = 200;
144  iMax = 0;
145  // writef<<iMax<<G4endl;
146  for( i = 1; i <= iMax; i++ )
147  {
148    x = 0.01*i;
149    // G4cout<<x<<"\t"<<SimpleMF(x)<<"\t"<<MediumMF(x)<<"\t"<<ComplexMF(x)<<G4endl;
150    // writef<<x<<"\t"<<SimpleMF(x)<<"\t"<<MediumMF(x)<<"\t"<<ComplexMF(x)<<G4endl;
151 
152    // G4cout<<x<<"\t"<<SimpleMPsi(x)<<"\t"<<MediumMPsi(x)<<"\t"<<ComplexMPsi(x)<<G4endl;
153    // writef<<x<<"\t"<<SimpleMPsi(x)<<"\t"<<MediumMPsi(x)<<"\t"<<ComplexMPsi(x)<<G4endl;
154
155    // G4cout<<x<<"\t"<<SimpleMG(x)<<"\t"<<MediumMG(x)<<"\t"<<ComplexMG(x)<<G4endl;
156    // writef<<x<<"\t"<<SimpleMG(x)<<"\t"<<MediumMG(x)<<"\t"<<ComplexMG(x)<<G4endl;
157  }
158
159  G4Element*     theElement;
160  G4Material*    theMaterial;
161  G4NistManager* man = G4NistManager::Instance();
162  man->SetVerbose(1);
163
164  G4cout << " 1 hydrogen" << G4endl;
165  G4cout << " 2 helium" << G4endl;
166  G4cout << " 4 berillium" << G4endl;
167  G4cout << " 6 carbon" << G4endl;
168  G4cout << " 7 nitrogen" << G4endl;
169  G4cout << " 8 oxigen" << G4endl;
170  G4cout << "13 aluminium" << G4endl;
171  G4cout << "14 silicon" << G4endl;
172  G4cout << "18 argon" << G4endl;
173  G4cout << "26 iron" << G4endl;
174  G4cout << "29 copper" << G4endl;
175  G4cout << "48 cadmium" << G4endl;
176  G4cout << "74 tugnsten" << G4endl;
177  G4cout << "77 iridium" << G4endl;
178  G4cout << "82 lead" << G4endl;
179  G4cout << "92 uranium" << G4endl;
180  G4int choice;
181  // G4cin >> choice;
182  choice = 13;
183
184
185  switch (choice)
186  {
187    case 1:
188
189      theElement  = man->FindOrBuildElement("H");
190      theMaterial = man->FindOrBuildMaterial("G4_H");
191      g4Xrad = theMaterial->GetRadlen();
192      break;
193
194    case 2:
195
196      theElement  = man->FindOrBuildElement("He");
197      theMaterial = man->FindOrBuildMaterial("G4_He");
198      g4Xrad = theMaterial->GetRadlen();
199      break;
200
201    case 4:
202
203      theElement  = man->FindOrBuildElement("Be");
204      theMaterial = man->FindOrBuildMaterial("G4_Be");
205      g4Xrad = theMaterial->GetRadlen();
206      break;
207
208    case 6:
209
210      theElement  = man->FindOrBuildElement("C");
211      theMaterial = man->FindOrBuildMaterial("G4_C");
212      g4Xrad = theMaterial->GetRadlen();
213      expXrad = 19.6*cm;
214      break;
215
216    case 7:
217
218      theElement  = man->FindOrBuildElement("N");
219      theMaterial = man->FindOrBuildMaterial("G4_N");
220      g4Xrad = theMaterial->GetRadlen();
221      break;
222
223
224    case 8:
225
226      theElement  = man->FindOrBuildElement("O");
227      theMaterial = man->FindOrBuildMaterial("G4_O");
228      g4Xrad = theMaterial->GetRadlen();
229      break;
230
231    case 13:
232
233      theElement  = man->FindOrBuildElement("Al");
234      theMaterial = man->FindOrBuildMaterial("G4_Al");
235      g4Xrad = theMaterial->GetRadlen();
236      expXrad = 8.9*cm;
237      break;
238
239    case 14:
240
241      theElement  = man->FindOrBuildElement("Si");
242      theMaterial = man->FindOrBuildMaterial("G4_Si");
243      g4Xrad = theMaterial->GetRadlen();
244      break;
245
246    case 18:
247
248      theElement  = man->FindOrBuildElement("Ar");
249      theMaterial = man->FindOrBuildMaterial("G4_Ar");
250      break;
251
252    case 26:
253
254      theElement  = man->FindOrBuildElement("Fe");
255      theMaterial = man->FindOrBuildMaterial("G4_Fe");
256      g4Xrad = theMaterial->GetRadlen();
257      expXrad = 1.76*cm;
258      break;
259
260    case 29:
261
262      theElement  = man->FindOrBuildElement("Cu");
263      theMaterial = man->FindOrBuildMaterial("G4_Cu");
264      g4Xrad = theMaterial->GetRadlen();
265      break;
266
267    case 48:
268
269      theElement  = man->FindOrBuildElement("Cd");
270      theMaterial = man->FindOrBuildMaterial("G4_Cd");
271      g4Xrad = theMaterial->GetRadlen();
272      break;
273
274
275    case 74:
276
277      theElement  = man->FindOrBuildElement("W");
278      theMaterial = man->FindOrBuildMaterial("G4_W");
279      g4Xrad = theMaterial->GetRadlen();
280      expXrad = 0.35*cm;
281      break;
282
283    case 77:
284
285      theElement  = man->FindOrBuildElement("Ir");
286      theMaterial = man->FindOrBuildMaterial("G4_Ir");
287      g4Xrad = theMaterial->GetRadlen();
288      break;
289
290    case 82:
291
292      theElement  = man->FindOrBuildElement("Pb");
293      theMaterial = man->FindOrBuildMaterial("G4_Pb");
294      g4Xrad = theMaterial->GetRadlen();
295      expXrad = 0.56*cm;
296      break;
297
298    case 92:
299
300      theElement  = man->FindOrBuildElement("U");
301      theMaterial = man->FindOrBuildMaterial("G4_U");
302      g4Xrad = theMaterial->GetRadlen();
303      expXrad = 0.35*cm;
304      break;
305  }
306
307// Particle definition
308
309  G4cout << " 1 electron" << G4endl;
310  G4cout << " 2 proton" << G4endl;
311  G4cout << " 3 pion+" << G4endl;
312  G4cout << " 4 pion-" << G4endl;
313  G4cout << " 4 muon+" << G4endl;
314  G4cout << " 5 muon-" << G4endl;
315
316  //  G4cin >> choice;
317  choice = 1;
318
319  G4ParticleDefinition* theParticleDefinition;
320
321
322  switch (choice)
323  {
324    case 1:
325
326      theParticleDefinition = G4Electron::ElectronDefinition(); 
327      break;
328
329    case 2:
330
331      theParticleDefinition = G4Proton::ProtonDefinition();
332      break;
333
334    case 3:
335
336      theParticleDefinition = G4PionPlus::PionPlusDefinition(); 
337      break;
338
339    case 4:
340
341      theParticleDefinition = G4PionMinus::PionMinusDefinition();
342      break;
343 
344    case 5:
345     
346      theParticleDefinition = G4MuonPlus::MuonPlusDefinition(); 
347      break;
348
349    case 6:
350     
351      theParticleDefinition = G4MuonMinus::MuonMinusDefinition();
352      break;
353
354  }
355
356  G4double energyMscXR, kinEnergy;
357
358
359  // kinEnergy = 8*GeV;   // 25.0*GeV;
360
361  kinEnergy = 25.0*GeV;
362
363
364  G4DynamicParticle*  theDynamicParticle = new G4DynamicParticle(theParticleDefinition,
365                                              G4ParticleMomentum(0.,0.,1.),
366                                              kinEnergy);
367
368  G4double m1 = theParticleDefinition->GetPDGMass();
369  G4double plab = theDynamicParticle->GetTotalMomentum();
370  G4cout <<"lab momentum, plab = "<<plab/GeV<<" GeV"<<G4endl;
371  G4double plabLowLimit = 20.0*MeV;
372
373  G4int Z   = G4int(theElement->GetZ());
374  G4int A    = G4int(theElement->GetN()+0.5);
375
376  G4double step = 4.10*mm;
377
378  step = expXrad;
379
380  G4double m2 = man->GetAtomicMassAmu(Z)*GeV;
381  // G4double m2 = man->GetAtomicMass( Z, A);
382  G4cout <<" target mass, m2 = "<<m2/GeV<<" GeV"<<G4endl<<G4endl;
383  G4cout <<"step = "<<step<<" mm; g4Xrad = "<<g4Xrad<<" mm; expXrad = "
384         <<expXrad<<"  mm"<<G4endl<<G4endl;
385
386
387  // G4MscRadiation* mscRad = new G4MscRadiation(theMaterial, step);
388  G4MscRadiation* mscRad = new G4MscRadiation(theMaterial, step, theDynamicParticle);
389
390  mscRad->SetVerboseLevel(0);
391
392  G4double numberMscXR, numberMGYXR, numberE146XR, absorption, lincofXR, radLength; 
393
394  G4double sRe, sIm, mcRe, mcIm, tmRe, tmIm, tmc, sf;
395
396  G4complex ms, mc, tm;
397
398  G4double s, phi, xi, G, psi, tmef;
399
400  iMax = 50;
401
402  writef<<iMax<<G4endl;
403
404  for( i = 0; i < iMax; i++ )
405  {
406    energyMscXR = std::exp(i*0.2)*0.1*MeV;
407    lincofXR = mscRad->GetPlateLinearPhotoAbs(energyMscXR);
408    absorption = (1. - std::exp(-lincofXR*step))/lincofXR;
409 
410    /*   
411   
412    numberMscXR  = mscRad->CalculateMscDiffdNdx(theDynamicParticle,energyMscXR);
413    numberMGYXR  = mscRad->CalculateMscMigdalDiffdNdx(theDynamicParticle,energyMscXR);
414    numberE146XR = mscRad->CalculateMscE146DiffdNdx(energyMscXR);
415
416    radLength = mscRad->GetRadLength();
417
418    numberMscXR  *= energyMscXR; // *expXrad;
419    numberMGYXR  *= energyMscXR; // *expXrad;
420    numberE146XR *= energyMscXR; // *expXrad;
421
422    // numberMscXR *= expXrad;
423    // numberMscXR *= expXrad;
424
425    numberMscXR  *= absorption;
426    numberMGYXR  *= absorption;
427    numberE146XR *= absorption;
428
429    G4cout <<"effStep = "<<absorption<<" mm; eXR  = "
430           <<energyMscXR/MeV<<" MeV; MscXR =  "
431           <<numberMscXR<<" "<<"; MGYXR =  "
432           <<numberMGYXR<<" "<<"; E146XR =  "
433           <<numberE146XR<<" "<<G4endl;
434    writef <<energyMscXR/MeV<<"\t"<<numberMscXR<<"\t"<<numberMGYXR<<"\t"<<numberE146XR<<G4endl;
435
436   
437   
438
439    mscRad->CalculateCorrectionTMGY(energyMscXR);
440    tm = mscRad->GetCorrectionTMGY();
441    ms = mscRad->CalculateMigdalS(energyMscXR);
442    mc = mscRad->CalculateCorrectionMsc(energyMscXR);
443    tmc = real(mc/tm);
444   
445
446    // sf = mscRad->SupressionFunction(kinEnergy, energyMscXR);
447
448     mscRad->CalcLPMFunctions(kinEnergy, energyMscXR);
449
450     s = mscRad->GetMigdalS();     
451     phi = mscRad->GetMigdalPhi();     
452     xi = mscRad->GetMigdalXi();     
453     psi = mscRad->GetMigdalPsi();     
454     G = mscRad->GetMigdalG();
455     tmef = mscRad->GetTMEffect();   
456         
457    G4cout <<energyMscXR/MeV<<"\t"<<real(tm)<<"\t"<<imag(tm)
458           <<"\t"<<real(ms)<<"\t"<<imag(ms)
459           <<"\t"<<real(mc)<<"\t"<<imag(mc)<<"\t"<<sf<<G4endl;
460    writef<<energyMscXR/MeV<<"\t"<<real(tm)<<"\t"<<real(mc)<<"\t"<<tmc<<"\t"<<sf<<G4endl;
461   
462   
463
464    // G4cout <<energyMscXR/MeV<<"\t"<<sf<<G4endl;
465    // G4cout <<energyMscXR/MeV<<"\t"<<s<<G4endl;
466    // writef <<energyMscXR/MeV<<"\t"<<s<<G4endl;
467    // writef <<energyMscXR/MeV<<"\t"<<real(mc)<<G4endl;
468    // writef <<energyMscXR/MeV<<"\t"<<s<<G4endl;
469
470
471    G4cout <<energyMscXR/MeV<<"\t"<<phi<<"\t"<<xi<<"\t"<<psi<<"\t"<<G<<"\t"<<tmef<<G4endl;
472    writef <<energyMscXR/MeV<<"\t"<<phi<<"\t"<<xi<<"\t"<<psi<<"\t"<<G<<"\t"<<tmef<<G4endl;
473
474*/ 
475
476
477
478  }
479
480
481  G4double eTkin, lambda;
482  iMax = 50;
483 
484  // writef<<iMax<<G4endl;
485
486  for( i = 0; i < iMax; i++ )
487  {
488
489    eTkin = std::exp(i*0.5)*0.001*MeV;
490
491    G4DynamicParticle*  aDynamicParticle = new G4DynamicParticle(theParticleDefinition,
492                                              G4ParticleMomentum(0.,0.,1.),
493                                              eTkin);
494
495    G4MscRadiation* aMscRad = new G4MscRadiation(theMaterial, step, aDynamicParticle);
496
497    lambda = aMscRad->CalculateMscE146Lambda();
498
499    G4cout <<eTkin/GeV<<"\t"<<lambda<<G4endl;
500    writef <<eTkin/GeV<<"\t"<<lambda<<G4endl;
501
502
503    delete aMscRad;
504    delete aDynamicParticle;
505  }
506
507
508
509
510  return 1;
511} // end of main
Note: See TracBrowser for help on using the repository browser.