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