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