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