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

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

geant4 tag 9.4
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  fParticleChange = GetParticleChangeForGamma();
135   
136  isInitialised = true;
137}
138
139//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
140
141G4double G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom(
142                                       const G4ParticleDefinition*,
143                                             G4double GammaEnergy,
144                                             G4double Z, G4double,
145                                             G4double, G4double)
146{
147  if (verboseLevel > 3)
148    G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedGammaConversionModel" << G4endl;
149
150  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
151  return cs;
152}
153
154//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
155
156void G4LivermorePolarizedGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
157                                                               const G4MaterialCutsCouple* couple,
158                                                               const G4DynamicParticle* aDynamicGamma,
159                                                               G4double,
160                                                               G4double)
161{
162
163  // Fluorescence generated according to:
164  // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
165  // subshell ionisation by a particle or due to deexcitation or decay of radionuclides",
166  // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
167
168  if (verboseLevel > 3)
169    G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedGammaConversionModel" << G4endl;
170
171  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
172  // Within energy limit?
173
174  if(photonEnergy <= lowEnergyLimit)
175    {
176      fParticleChange->ProposeTrackStatus(fStopAndKill);
177      fParticleChange->SetProposedKineticEnergy(0.);
178      return;
179    }
180
181
182  G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
183  G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
184
185  // Make sure that the polarization vector is perpendicular to the
186  // gamma direction. If not
187
188  if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
189    { // only for testing now
190      gammaPolarization0 = GetRandomPolarization(gammaDirection0);
191    }
192  else
193    {
194      if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
195        {
196          gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
197        }
198    }
199
200  // End of Protection
201
202
203  G4double epsilon ;
204  G4double epsilon0 = electron_mass_c2 / photonEnergy ;
205
206  // Do it fast if photon energy < 2. MeV
207
208  if (photonEnergy < smallEnergy )
209    {
210      epsilon = epsilon0 + (0.5 - epsilon0) * G4UniformRand();
211    }
212  else
213    {
214
215 // Select randomly one element in the current material
216
217      //     G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
218      //const G4Element* element = crossSectionHandler->SelectRandomElement(couple,photonEnergy);
219
220      const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
221      const G4Element* element = SelectRandomAtom(couple,particle,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.