source: trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedGammaConversionModel.cc@ 1314

Last change on this file since 1314 was 1197, checked in by garnier, 16 years ago

nvx fichiers dans CVS
File size: 26.0 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#include "G4LivermorePolarizedGammaConversionModel.hh"
28
29//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
30
31using namespace std;
32
33//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
34
35G4LivermorePolarizedGammaConversionModel::G4LivermorePolarizedGammaConversionModel(const G4ParticleDefinition*,
36 const G4String& nam)
37 :G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),crossSectionHandler(0)
38{
39 lowEnergyLimit = 1.0220000 * MeV;
40 highEnergyLimit = 100 * GeV;
41 SetLowEnergyLimit(lowEnergyLimit);
42 SetHighEnergyLimit(highEnergyLimit);
43 smallEnergy = 2.*MeV;
44
45
46 verboseLevel= 0;
47 // Verbosity scale:
48 // 0 = nothing
49 // 1 = warning for energy non-conservation
50 // 2 = details of energy budget
51 // 3 = calculation of cross sections, file openings, samping of atoms
52 // 4 = entering in methods
53
54 G4cout << "Livermore Polarized GammaConversion is constructed " << G4endl
55 << "Energy range: "
56 << lowEnergyLimit / keV << " keV - "
57 << highEnergyLimit / GeV << " GeV"
58 << G4endl;
59
60 crossSectionHandler = new G4CrossSectionHandler();
61 crossSectionHandler->Initialise(0,1.0220*MeV,100.*GeV,400);
62
63
64}
65
66//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67
68G4LivermorePolarizedGammaConversionModel::~G4LivermorePolarizedGammaConversionModel()
69{
70 // if (meanFreePathTable) delete meanFreePathTable;
71 if (crossSectionHandler) delete crossSectionHandler;
72}
73
74
75
76//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77
78void G4LivermorePolarizedGammaConversionModel::Initialise(const G4ParticleDefinition* particle,
79 const G4DataVector& cuts)
80{
81 if (verboseLevel > 3)
82 G4cout << "Calling G4LivermorePolarizedGammaConversionModel::Initialise()" << G4endl;
83
84 if (crossSectionHandler)
85 {
86 crossSectionHandler->Clear();
87 delete crossSectionHandler;
88 }
89
90
91 // Energy limits
92
93 if (LowEnergyLimit() < lowEnergyLimit)
94 {
95 G4cout << "G4LivermorePolarizedGammaConversionModel: low energy limit increased from " <<
96 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
97 SetLowEnergyLimit(lowEnergyLimit);
98 }
99
100 if (HighEnergyLimit() > highEnergyLimit)
101 {
102 G4cout << "G4LivermorePolarizedGammaConversionModel: high energy limit decreased from " <<
103 HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
104 SetHighEnergyLimit(highEnergyLimit);
105 }
106
107 // Reading of data files - all materials are read
108
109 crossSectionHandler = new G4CrossSectionHandler;
110 crossSectionHandler->Clear();
111 G4String crossSectionFile = "pair/pp-cs-";
112 crossSectionHandler->LoadData(crossSectionFile);
113
114 // meanFreePathTable = 0;
115 //meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
116
117
118 //
119 if (verboseLevel > 2)
120 G4cout << "Loaded cross section files for Livermore Polarized GammaConversion model" << G4endl;
121
122 InitialiseElementSelectors(particle,cuts);
123
124 G4cout << "Livermore Polarized GammaConversion model is initialized " << G4endl
125 << "Energy range: "
126 << LowEnergyLimit() / keV << " keV - "
127 << HighEnergyLimit() / GeV << " GeV"
128 << G4endl;
129
130 //
131
132 if(isInitialised) return;
133
134 if(pParticleChange)
135 fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
136 else
137 fParticleChange = new G4ParticleChangeForGamma();
138
139 isInitialised = true;
140}
141
142//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
143
144G4double G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom(
145 const G4ParticleDefinition*,
146 G4double GammaEnergy,
147 G4double Z, G4double,
148 G4double, G4double)
149{
150 if (verboseLevel > 3)
151 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedGammaConversionModel" << G4endl;
152
153 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
154 return cs;
155}
156
157//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
158
159void G4LivermorePolarizedGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
160 const G4MaterialCutsCouple* couple,
161 const G4DynamicParticle* aDynamicGamma,
162 G4double,
163 G4double)
164{
165
166 // Fluorescence generated according to:
167 // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
168 // subshell ionisation by a particle or due to deexcitation or decay of radionuclides",
169 // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
170
171 if (verboseLevel > 3)
172 G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedGammaConversionModel" << G4endl;
173
174 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
175 // Within energy limit?
176
177 if(photonEnergy <= lowEnergyLimit)
178 {
179 fParticleChange->ProposeTrackStatus(fStopAndKill);
180 fParticleChange->SetProposedKineticEnergy(0.);
181 return;
182 }
183
184
185 G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
186 G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
187
188 // Make sure that the polarization vector is perpendicular to the
189 // gamma direction. If not
190
191 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
192 { // only for testing now
193 gammaPolarization0 = GetRandomPolarization(gammaDirection0);
194 }
195 else
196 {
197 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
198 {
199 gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
200 }
201 }
202
203 // End of Protection
204
205
206 G4double epsilon ;
207 G4double epsilon0 = electron_mass_c2 / photonEnergy ;
208
209 // Do it fast if photon energy < 2. MeV
210
211 if (photonEnergy < smallEnergy )
212 {
213 epsilon = epsilon0 + (0.5 - epsilon0) * G4UniformRand();
214 }
215 else
216 {
217
218 // Select randomly one element in the current material
219
220 // G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
221 const G4Element* element = crossSectionHandler->SelectRandomElement(couple,photonEnergy);
222
223 if (element == 0)
224 {
225 G4cout << "G4LivermorePolarizedGammaConversionModel::PostStepDoIt - element = 0" << G4endl;
226 }
227 G4IonisParamElm* ionisation = element->GetIonisation();
228 if (ionisation == 0)
229 {
230 G4cout << "G4LivermorePolarizedGammaConversionModel::PostStepDoIt - ionisation = 0" << G4endl;
231 }
232
233
234
235 // Extract Coulomb factor for this Element
236
237 G4double fZ = 8. * (ionisation->GetlogZ3());
238 if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
239
240 // Limits of the screening variable
241 G4double screenFactor = 136. * epsilon0 / (element->GetIonisation()->GetZ3()) ;
242 G4double screenMax = exp ((42.24 - fZ)/8.368) - 0.952 ;
243 G4double screenMin = std::min(4.*screenFactor,screenMax) ;
244
245 // Limits of the energy sampling
246 G4double epsilon1 = 0.5 - 0.5 * sqrt(1. - screenMin / screenMax) ;
247 G4double epsilonMin = std::max(epsilon0,epsilon1);
248 G4double epsilonRange = 0.5 - epsilonMin ;
249
250 // Sample the energy rate of the created electron (or positron)
251 G4double screen;
252 G4double gReject ;
253
254 G4double f10 = ScreenFunction1(screenMin) - fZ;
255 G4double f20 = ScreenFunction2(screenMin) - fZ;
256 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
257 G4double normF2 = std::max(1.5 * f20,0.);
258
259 do {
260 if (normF1 / (normF1 + normF2) > G4UniformRand() )
261 {
262 epsilon = 0.5 - epsilonRange * pow(G4UniformRand(), 0.3333) ;
263 screen = screenFactor / (epsilon * (1. - epsilon));
264 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
265 }
266 else
267 {
268 epsilon = epsilonMin + epsilonRange * G4UniformRand();
269 screen = screenFactor / (epsilon * (1 - epsilon));
270 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
271
272
273 }
274 } while ( gReject < G4UniformRand() );
275
276 } // End of epsilon sampling
277
278 // Fix charges randomly
279
280 G4double electronTotEnergy;
281 G4double positronTotEnergy;
282
283
284 if (CLHEP::RandBit::shootBit())
285 {
286 electronTotEnergy = (1. - epsilon) * photonEnergy;
287 positronTotEnergy = epsilon * photonEnergy;
288 }
289 else
290 {
291 positronTotEnergy = (1. - epsilon) * photonEnergy;
292 electronTotEnergy = epsilon * photonEnergy;
293 }
294
295 // Scattered electron (positron) angles. ( Z - axis along the parent photon)
296 // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
297 // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
298
299 G4double u;
300 const G4double a1 = 0.625;
301 G4double a2 = 3. * a1;
302
303 if (0.25 > G4UniformRand())
304 {
305 u = - log(G4UniformRand() * G4UniformRand()) / a1 ;
306 }
307 else
308 {
309 u = - log(G4UniformRand() * G4UniformRand()) / a2 ;
310 }
311
312 G4double Ene = electronTotEnergy/electron_mass_c2; // Normalized energy
313
314 G4double cosTheta = 0.;
315 G4double sinTheta = 0.;
316
317 SetTheta(&cosTheta,&sinTheta,Ene);
318
319 // G4double theta = u * electron_mass_c2 / photonEnergy ;
320 // G4double phi = twopi * G4UniformRand() ;
321
322 G4double phi,psi=0.;
323
324 //corrected e+ e- angular angular distribution //preliminary!
325
326 // if(photonEnergy>50*MeV)
327 // {
328 phi = SetPhi(photonEnergy);
329 psi = SetPsi(photonEnergy,phi);
330 // }
331 //else
332 // {
333 //psi = G4UniformRand()*2.*pi;
334 //phi = pi; // coplanar
335 // }
336
337 Psi = psi;
338 Phi = phi;
339 //G4cout << "PHI " << phi << G4endl;
340 //G4cout << "PSI " << psi << G4endl;
341
342 G4double phie = psi; //azimuthal angle for the electron
343
344 G4double dirX = sinTheta*cos(phie);
345 G4double dirY = sinTheta*sin(phie);
346 G4double dirZ = cosTheta;
347 G4ThreeVector electronDirection(dirX,dirY,dirZ);
348 // Kinematics of the created pair:
349 // the electron and positron are assumed to have a symetric angular
350 // distribution with respect to the Z axis along the parent photon
351
352 //G4double localEnergyDeposit = 0. ;
353
354 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
355
356 SystemOfRefChange(gammaDirection0,electronDirection,
357 gammaPolarization0);
358
359 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
360 electronDirection,
361 electronKineEnergy);
362
363 // The e+ is always created (even with kinetic energy = 0) for further annihilation
364
365 Ene = positronTotEnergy/electron_mass_c2; // Normalized energy
366
367 cosTheta = 0.;
368 sinTheta = 0.;
369
370 SetTheta(&cosTheta,&sinTheta,Ene);
371 G4double phip = phie+phi; //azimuthal angle for the positron
372
373 dirX = sinTheta*cos(phip);
374 dirY = sinTheta*sin(phip);
375 dirZ = cosTheta;
376 G4ThreeVector positronDirection(dirX,dirY,dirZ);
377
378 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
379 SystemOfRefChange(gammaDirection0,positronDirection,
380 gammaPolarization0);
381
382 // Create G4DynamicParticle object for the particle2
383 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
384 positronDirection, positronKineEnergy);
385
386
387 fvect->push_back(particle1);
388 fvect->push_back(particle2);
389
390
391
392 // Kill the incident photon
393
394
395
396 // Create lists of pointers to DynamicParticles (photons and electrons)
397 // (Is the electron vector necessary? To be checked)
398 // std::vector<G4DynamicParticle*>* photonVector = 0;
399 //std::vector<G4DynamicParticle*> electronVector;
400
401 fParticleChange->ProposeMomentumDirection( 0., 0., 0. );
402 fParticleChange->SetProposedKineticEnergy(0.);
403 fParticleChange->ProposeTrackStatus(fStopAndKill);
404
405}
406
407//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
408
409G4double G4LivermorePolarizedGammaConversionModel::ScreenFunction1(G4double screenVariable)
410{
411 // Compute the value of the screening function 3*phi1 - phi2
412
413 G4double value;
414
415 if (screenVariable > 1.)
416 value = 42.24 - 8.368 * log(screenVariable + 0.952);
417 else
418 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
419
420 return value;
421}
422
423
424
425G4double G4LivermorePolarizedGammaConversionModel::ScreenFunction2(G4double screenVariable)
426{
427 // Compute the value of the screening function 1.5*phi1 - 0.5*phi2
428
429 G4double value;
430
431 if (screenVariable > 1.)
432 value = 42.24 - 8.368 * log(screenVariable + 0.952);
433 else
434 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
435
436 return value;
437}
438
439
440void G4LivermorePolarizedGammaConversionModel::SetTheta(G4double* p_cosTheta, G4double* p_sinTheta, G4double Energy)
441{
442
443 // to avoid computational errors since Theta could be very small
444 // Energy in Normalized Units (!)
445
446 G4double Momentum = sqrt(Energy*Energy -1);
447 G4double Rand = G4UniformRand();
448
449 *p_cosTheta = (Energy*((2*Rand)- 1) + Momentum)/((Momentum*(2*Rand-1))+Energy);
450 *p_sinTheta = (2*sqrt(Rand*(1-Rand)))/(Momentum*(2*Rand-1)+Energy);
451}
452
453
454
455G4double G4LivermorePolarizedGammaConversionModel::SetPhi(G4double Energy)
456{
457
458
459 G4double value = 0.;
460 G4double Ene = Energy/MeV;
461
462 G4double pl[4];
463
464
465 G4double pt[2];
466 G4double xi = 0;
467 G4double xe = 0.;
468 G4double n1=0.;
469 G4double n2=0.;
470
471
472 if (Ene>=50.)
473 {
474 const G4double ay0=5.6, by0=18.6, aa0=2.9, ba0 = 8.16E-3;
475 const G4double aw = 0.0151, bw = 10.7, cw = -410.;
476
477 const G4double axc = 3.1455, bxc = -1.11, cxc = 310.;
478
479 pl[0] = Fln(ay0,by0,Ene);
480 pl[1] = aa0 + ba0*(Ene);
481 pl[2] = Poli(aw,bw,cw,Ene);
482 pl[3] = Poli(axc,bxc,cxc,Ene);
483
484 const G4double abf = 3.1216, bbf = 2.68;
485 pt[0] = -1.4;
486 pt[1] = abf + bbf/Ene;
487
488
489
490 //G4cout << "PL > 50. "<< pl[0] << " " << pl[1] << " " << pl[2] << " " <<pl[3] << " " << G4endl;
491
492 xi = 3.0;
493 xe = Encu(pl,pt,xi);
494 //G4cout << "ENCU "<< xe << G4endl;
495 n1 = Fintlor(pl,pi) - Fintlor(pl,xe);
496 n2 = Finttan(pt,xe) - Finttan(pt,0.);
497 }
498 else
499 {
500 const G4double ay0=0.144, by0=0.11;
501 const G4double aa0=2.7, ba0 = 2.74;
502 const G4double aw = 0.21, bw = 10.8, cw = -58.;
503 const G4double axc = 3.17, bxc = -0.87, cxc = -6.;
504
505 pl[0] = Fln(ay0, by0, Ene);
506 pl[1] = Fln(aa0, ba0, Ene);
507 pl[2] = Poli(aw,bw,cw,Ene);
508 pl[3] = Poli(axc,bxc,cxc,Ene);
509
510 //G4cout << "PL < 50."<< pl[0] << " " << pl[1] << " " << pl[2] << " " <<pl[3] << " " << G4endl;
511 //G4cout << "ENCU "<< xe << G4endl;
512 n1 = Fintlor(pl,pi) - Fintlor(pl,xe);
513
514 }
515
516
517 G4double n=0.;
518 n = n1+n2;
519
520 G4double c1 = 0.;
521 c1 = Glor(pl, xe);
522
523 G4double xm = 0.;
524 xm = Flor(pl,pl[3])*Glor(pl,pl[3]);
525
526 G4double r1,r2,r3;
527 G4double xco=0.;
528
529 if (Ene>=50.)
530 {
531 r1= G4UniformRand();
532 if( r1>=n2/n)
533 {
534 do
535 {
536 r2 = G4UniformRand();
537 value = Finvlor(pl,xe,r2);
538 xco = Glor(pl,value)/c1;
539 r3 = G4UniformRand();
540 } while(r3>=xco);
541 }
542 else
543 {
544 value = Finvtan(pt,n,r1);
545 }
546 }
547 else
548 {
549 do
550 {
551 r2 = G4UniformRand();
552 value = Finvlor(pl,xe,r2);
553 xco = Glor(pl,value)/c1;
554 r3 = G4UniformRand();
555 } while(r3>=xco);
556 }
557
558 // G4cout << "PHI = " <<value << G4endl;
559 return value;
560}
561G4double G4LivermorePolarizedGammaConversionModel::SetPsi(G4double Energy, G4double Phi)
562{
563
564 G4double value = 0.;
565 G4double Ene = Energy/MeV;
566
567 G4double p0l[4];
568 G4double ppml[4];
569 G4double p0t[2];
570 G4double ppmt[2];
571
572 G4double xi = 0.;
573 G4double xe0 = 0.;
574 G4double xepm = 0.;
575
576 if (Ene>=50.)
577 {
578 const G4double ay00 = 3.4, by00 = 9.8, aa00 = 1.34, ba00 = 5.3;
579 const G4double aw0 = 0.014, bw0 = 9.7, cw0 = -2.E4;
580 const G4double axc0 = 3.1423, bxc0 = -2.35, cxc0 = 0.;
581 const G4double ay0p = 1.53, by0p = 3.2, aap = 0.67, bap = 8.5E-3;
582 const G4double awp = 6.9E-3, bwp = 12.6, cwp = -3.8E4;
583 const G4double axcp = 2.8E-3,bxcp = -3.133;
584 const G4double abf0 = 3.1213, bbf0 = 2.61;
585 const G4double abfpm = 3.1231, bbfpm = 2.84;
586
587 p0l[0] = Fln(ay00, by00, Ene);
588 p0l[1] = Fln(aa00, ba00, Ene);
589 p0l[2] = Poli(aw0, bw0, cw0, Ene);
590 p0l[3] = Poli(axc0, bxc0, cxc0, Ene);
591
592 ppml[0] = Fln(ay0p, by0p, Ene);
593 ppml[1] = aap + bap*(Ene);
594 ppml[2] = Poli(awp, bwp, cwp, Ene);
595 ppml[3] = Fln(axcp,bxcp,Ene);
596
597 p0t[0] = -0.81;
598 p0t[1] = abf0 + bbf0/Ene;
599 ppmt[0] = -0.6;
600 ppmt[1] = abfpm + bbfpm/Ene;
601
602 //G4cout << "P0L > 50"<< p0l[0] << " " << p0l[1] << " " << p0l[2] << " " <<p0l[3] << " " << G4endl;
603 //G4cout << "PPML > 50"<< ppml[0] << " " << ppml[1] << " " << ppml[2] << " " <<ppml[3] << " " << G4endl;
604
605 xi = 3.0;
606 xe0 = Encu(p0l, p0t, xi);
607 //G4cout << "ENCU1 "<< xe0 << G4endl;
608 xepm = Encu(ppml, ppmt, xi);
609
610
611 //G4cout << "ENCU2 "<< xepm << G4endl;
612
613 }
614 else
615 {
616 const G4double ay00 = 2.82, by00 = 6.35;
617 const G4double aa00 = -1.75, ba00 = 0.25;
618
619 const G4double aw0 = 0.028, bw0 = 5., cw0 = -50.;
620 const G4double axc0 = 3.14213, bxc0 = -2.3, cxc0 = 5.7;
621 const G4double ay0p = 1.56, by0p = 3.6;
622 const G4double aap = 0.86, bap = 8.3E-3;
623 const G4double awp = 0.022, bwp = 7.4, cwp = -51.;
624 const G4double xcp = 3.1486;
625
626 p0l[0] = Fln(ay00, by00, Ene);
627 p0l[1] = aa00+pow(Ene, ba00);
628 p0l[2] = Poli(aw0, bw0, cw0, Ene);
629 p0l[3] = Poli(axc0, bxc0, cxc0, Ene);
630 ppml[0] = Fln(ay0p, by0p, Ene);
631 ppml[1] = aap + bap*(Ene);
632 ppml[2] = Poli(awp, bwp, cwp, Ene);
633 ppml[3] = xcp;
634
635 }
636
637
638
639 G4double a,b=0.;
640
641 if (Ene>=50.)
642 {
643 if (Phi>xepm)
644 {
645 b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml,Phi));
646 }
647 else
648 {
649 b = Ftan(ppmt,Phi);
650 }
651 if (Phi>xe0)
652 {
653 a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,Phi));
654 }
655 else
656 {
657 a = Ftan(p0t,Phi);
658 }
659 }
660 else
661 {
662 b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml,Phi));
663 a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,Phi));
664 }
665 G4double nr =0.;
666
667 if (b>a)
668 {
669 nr = 1./b;
670 }
671 else
672 {
673 nr = 1./a;
674 }
675
676 G4double r1,r2=0.;
677 G4double r3 =-1.;
678 do
679 {
680 r1 = G4UniformRand();
681 r2 = G4UniformRand();
682 value = r2*pi;
683 r3 = nr*(a*cos(value)*cos(value) + b*sin(value)*sin(value));
684 }while(r1>r3);
685
686 return value;
687}
688
689
690G4double G4LivermorePolarizedGammaConversionModel::Poli
691(G4double a, G4double b, G4double c, G4double x)
692{
693 G4double value=0.;
694 if(x>0.)
695 {
696 value =(a + b/x + c/(x*x*x));
697 }
698 else
699 {
700 //G4cout << "ERROR in Poli! " << G4endl;
701 }
702 return value;
703}
704G4double G4LivermorePolarizedGammaConversionModel::Fln
705(G4double a, G4double b, G4double x)
706{
707 G4double value=0.;
708 if(x>0.)
709 {
710 value =(a*log(x)-b);
711 }
712 else
713 {
714 //G4cout << "ERROR in Fln! " << G4endl;
715 }
716 return value;
717}
718
719
720G4double G4LivermorePolarizedGammaConversionModel::Encu
721(G4double* p_p1, G4double* p_p2, G4double x)
722{
723 G4double value=0.;
724 G4int i=0;
725 G4double fx = 1.;
726
727 do
728 {
729 x -= (Flor(p_p1, x)*Glor(p_p1,x) - Ftan(p_p2, x))/
730 (Fdlor(p_p1,x) - Fdtan(p_p2,x));
731 fx = Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p_p2, x);
732 i += 1;
733 //G4cout << abs(fx) << " " << i << " " << x << "dentro ENCU " << G4endl;
734 } while( (i<100) && (abs(fx) > 1e-6)) ;
735
736 if (i>100||x>pi) x = 3.0;
737 value = x;
738
739 if (value<0.) value=0.;
740
741 return value;
742}
743
744
745
746
747G4double G4LivermorePolarizedGammaConversionModel::Flor(G4double* p_p1, G4double x)
748{
749 G4double value =0.;
750 // G4double y0 = p_p1[0];
751 // G4double A = p_p1[1];
752 G4double w = p_p1[2];
753 G4double xc = p_p1[3];
754
755 value = 1./(pi*(w*w + 4.*(x-xc)*(x-xc)));
756 return value;
757}
758
759G4double G4LivermorePolarizedGammaConversionModel::Glor(G4double* p_p1, G4double x)
760{
761 G4double value =0.;
762 G4double y0 = p_p1[0];
763 G4double A = p_p1[1];
764 G4double w = p_p1[2];
765 G4double xc = p_p1[3];
766
767 value = (y0 *pi*(w*w + 4.*(x-xc)*(x-xc)) + 2.*A*w);
768 return value;
769}
770
771G4double G4LivermorePolarizedGammaConversionModel::Fdlor(G4double* p_p1, G4double x)
772{
773 G4double value =0.;
774 //G4double y0 = p_p1[0];
775 G4double A = p_p1[1];
776 G4double w = p_p1[2];
777 G4double xc = p_p1[3];
778
779 value = (-16.*A*w*(x-xc))/
780 (pi*(w*w+4.*(x-xc)*(x-xc))*(w*w+4.*(x-xc)*(x-xc)));
781 return value;
782}
783
784
785G4double G4LivermorePolarizedGammaConversionModel::Fintlor(G4double* p_p1, G4double x)
786{
787 G4double value =0.;
788 G4double y0 = p_p1[0];
789 G4double A = p_p1[1];
790 G4double w = p_p1[2];
791 G4double xc = p_p1[3];
792
793 value = y0*x + A*atan( 2*(x-xc)/w) / pi;
794 return value;
795}
796
797
798G4double G4LivermorePolarizedGammaConversionModel::Finvlor(G4double* p_p1, G4double x, G4double r)
799{
800 G4double value = 0.;
801 G4double nor = 0.;
802 //G4double y0 = p_p1[0];
803 // G4double A = p_p1[1];
804 G4double w = p_p1[2];
805 G4double xc = p_p1[3];
806
807 nor = atan(2.*(pi-xc)/w)/(2.*pi*w) - atan(2.*(x-xc)/w)/(2.*pi*w);
808 value = xc - (w/2.)*tan(-2.*r*nor*pi*w+atan(2*(xc-x)/w));
809
810 return value;
811}
812
813G4double G4LivermorePolarizedGammaConversionModel::Ftan(G4double* p_p1, G4double x)
814{
815 G4double value =0.;
816 G4double a = p_p1[0];
817 G4double b = p_p1[1];
818
819 value = a /(x-b);
820 return value;
821}
822
823
824G4double G4LivermorePolarizedGammaConversionModel::Fdtan(G4double* p_p1, G4double x)
825{
826 G4double value =0.;
827 G4double a = p_p1[0];
828 G4double b = p_p1[1];
829
830 value = -1.*a / ((x-b)*(x-b));
831 return value;
832}
833
834
835G4double G4LivermorePolarizedGammaConversionModel::Finttan(G4double* p_p1, G4double x)
836{
837 G4double value =0.;
838 G4double a = p_p1[0];
839 G4double b = p_p1[1];
840
841
842 value = a*log(b-x);
843 return value;
844}
845
846G4double G4LivermorePolarizedGammaConversionModel::Finvtan(G4double* p_p1, G4double cnor, G4double r)
847{
848 G4double value =0.;
849 G4double a = p_p1[0];
850 G4double b = p_p1[1];
851
852 value = b*(1-exp(r*cnor/a));
853
854 return value;
855}
856
857
858
859
860//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
861
862G4ThreeVector G4LivermorePolarizedGammaConversionModel::SetPerpendicularVector(G4ThreeVector& a)
863{
864 G4double dx = a.x();
865 G4double dy = a.y();
866 G4double dz = a.z();
867 G4double x = dx < 0.0 ? -dx : dx;
868 G4double y = dy < 0.0 ? -dy : dy;
869 G4double z = dz < 0.0 ? -dz : dz;
870 if (x < y) {
871 return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
872 }else{
873 return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
874 }
875}
876
877//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
878
879G4ThreeVector G4LivermorePolarizedGammaConversionModel::GetRandomPolarization(G4ThreeVector& direction0)
880{
881 G4ThreeVector d0 = direction0.unit();
882 G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
883 G4ThreeVector a0 = a1.unit(); // unit vector
884
885 G4double rand1 = G4UniformRand();
886
887 G4double angle = twopi*rand1; // random polar angle
888 G4ThreeVector b0 = d0.cross(a0); // cross product
889
890 G4ThreeVector c;
891
892 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
893 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
894 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
895
896 G4ThreeVector c0 = c.unit();
897
898 return c0;
899
900}
901
902//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
903
904G4ThreeVector G4LivermorePolarizedGammaConversionModel::GetPerpendicularPolarization
905(const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
906{
907
908 //
909 // The polarization of a photon is always perpendicular to its momentum direction.
910 // Therefore this function removes those vector component of gammaPolarization, which
911 // points in direction of gammaDirection
912 //
913 // Mathematically we search the projection of the vector a on the plane E, where n is the
914 // plains normal vector.
915 // The basic equation can be found in each geometry book (e.g. Bronstein):
916 // p = a - (a o n)/(n o n)*n
917
918 return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
919}
920
921//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
922
923
924void G4LivermorePolarizedGammaConversionModel::SystemOfRefChange
925 (G4ThreeVector& direction0,G4ThreeVector& direction1,
926 G4ThreeVector& polarization0)
927{
928 // direction0 is the original photon direction ---> z
929 // polarization0 is the original photon polarization ---> x
930 // need to specify y axis in the real reference frame ---> y
931 G4ThreeVector Axis_Z0 = direction0.unit();
932 G4ThreeVector Axis_X0 = polarization0.unit();
933 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
934
935 G4double direction_x = direction1.getX();
936 G4double direction_y = direction1.getY();
937 G4double direction_z = direction1.getZ();
938
939 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
940
941}
942
943
944
945
Note: See TracBrowser for help on using the repository browser.