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

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