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