Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/eRosita/physics/src/G4LowEnergyPolarizedCompton.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 //
 27 //
 28 // ------------------------------------------------------------
 29 //      GEANT 4 class implementation file
 30 //      CERN Geneva Switzerland
 31 //
 32 
 33 // --------- G4LowEnergyPolarizedCompton class -----
 34 //
 35 //           by G.Depaola & F.Longo (21 may 2001)
 36 //
 37 // 21 May 2001 - MGP      Modified to inherit from G4VDiscreteProcess
 38 //                        Applies same algorithm as LowEnergyCompton
 39 //                        if the incoming photon is not polarised
 40 //                        Temporary protection to avoid crash in the case 
 41 //                        of polarisation || incident photon direction
 42 //
 43 // 17 October 2001 - F.Longo - Revised according to a design iteration
 44 //
 45 // 21 February 2002 - F.Longo Revisions with A.Zoglauer and G.Depaola
 46 //                            - better description of parallelism
 47 //                            - system of ref change method improved
 48 // 22 January  2003 - V.Ivanchenko Cut per region
 49 // 24 April    2003 - V.Ivanchenko Cut per region mfpt
 50 //
 51 //
 52 // ************************************************************
 53 //
 54 // Corrections by Rui Curado da Silva (2000)
 55 // New Implementation by G.Depaola & F.Longo
 56 //
 57 // - sampling of phi
 58 // - polarization of scattered photon
 59 //
 60 // --------------------------------------------------------------
 61 
 62 #include "G4LowEnergyPolarizedCompton.hh"
 63 #include "Randomize.hh"
 64 #include "G4PhysicalConstants.hh"
 65 #include "G4SystemOfUnits.hh"
 66 #include "G4ParticleDefinition.hh"
 67 #include "G4Track.hh"
 68 #include "G4Step.hh"
 69 #include "G4ForceCondition.hh"
 70 #include "G4Gamma.hh"
 71 #include "G4Electron.hh"
 72 #include "G4DynamicParticle.hh"
 73 #include "G4VParticleChange.hh"
 74 #include "G4ThreeVector.hh"
 75 #include "G4RDVCrossSectionHandler.hh"
 76 #include "G4RDCrossSectionHandler.hh"
 77 #include "G4RDVEMDataSet.hh"
 78 #include "G4RDCompositeEMDataSet.hh"
 79 #include "G4RDVDataSetAlgorithm.hh"
 80 #include "G4RDLogLogInterpolation.hh"
 81 #include "G4RDVRangeTest.hh"
 82 #include "G4RDRangeTest.hh"
 83 #include "G4MaterialCutsCouple.hh"
 84 
 85 // constructor
 86 
 87 G4LowEnergyPolarizedCompton::G4LowEnergyPolarizedCompton(const G4String& processName)
 88   : G4VDiscreteProcess(processName),
 89     lowEnergyLimit (250*eV),              // initialization
 90     highEnergyLimit(100*GeV),
 91     intrinsicLowEnergyLimit(10*eV),
 92     intrinsicHighEnergyLimit(100*GeV)
 93 {
 94   if (lowEnergyLimit < intrinsicLowEnergyLimit ||
 95       highEnergyLimit > intrinsicHighEnergyLimit)
 96     {
 97       G4Exception("G4LowEnergyPolarizedCompton::G4LowEnergyPolarizedCompton()",
 98                   "OutOfRange", FatalException,
 99                   "Energy outside intrinsic process validity range!");
100     }
101 
102   crossSectionHandler = new G4RDCrossSectionHandler;
103 
104 
105   G4RDVDataSetAlgorithm* scatterInterpolation = new G4RDLogLogInterpolation;
106   G4String scatterFile = "comp/ce-sf-";
107   scatterFunctionData = new G4RDCompositeEMDataSet(scatterInterpolation,1.,1.);
108   scatterFunctionData->LoadData(scatterFile);
109 
110   meanFreePathTable = 0;
111 
112   rangeTest = new G4RDRangeTest;
113 
114   // For Doppler broadening
115   shellData.SetOccupancyData();
116 
117 
118    if (verboseLevel > 0)
119      {
120        G4cout << GetProcessName() << " is created " << G4endl
121         << "Energy range: "
122         << lowEnergyLimit / keV << " keV - "
123         << highEnergyLimit / GeV << " GeV"
124         << G4endl;
125      }
126 }
127 
128 // destructor
129 
130 G4LowEnergyPolarizedCompton::~G4LowEnergyPolarizedCompton()
131 {
132   delete meanFreePathTable;
133   delete crossSectionHandler;
134   delete scatterFunctionData;
135   delete rangeTest;
136 }
137 
138 
139 void G4LowEnergyPolarizedCompton::BuildPhysicsTable(const G4ParticleDefinition& )
140 {
141 
142   crossSectionHandler->Clear();
143   G4String crossSectionFile = "comp/ce-cs-";
144   crossSectionHandler->LoadData(crossSectionFile);
145   delete meanFreePathTable;
146   meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
147 
148   // For Doppler broadening
149   G4String file = "/doppler/shell-doppler";
150   shellData.LoadData(file);
151 
152 }
153 
154 G4VParticleChange* G4LowEnergyPolarizedCompton::PostStepDoIt(const G4Track& aTrack,
155                    const G4Step&  aStep)
156 {
157   // The scattered gamma energy is sampled according to Klein - Nishina formula.
158   // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15).
159   // GEANT4 internal units
160   //
161   // Note : Effects due to binding of atomic electrons are negliged.
162 
163   aParticleChange.Initialize(aTrack);
164 
165   // Dynamic particle quantities
166   const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
167   G4double gammaEnergy0 = incidentPhoton->GetKineticEnergy();
168   G4ThreeVector gammaPolarization0 = incidentPhoton->GetPolarization();
169 
170   //  gammaPolarization0 = gammaPolarization0.unit(); //
171 
172   // Protection: a polarisation parallel to the
173   // direction causes problems;
174   // in that case find a random polarization
175 
176   G4ThreeVector gammaDirection0 = incidentPhoton->GetMomentumDirection();
177   // ---- MGP ---- Next two lines commented out to remove compilation warnings
178   // G4double scalarproduct = gammaPolarization0.dot(gammaDirection0);
179   // G4double angle = gammaPolarization0.angle(gammaDirection0);
180 
181   // Make sure that the polarization vector is perpendicular to the
182   // gamma direction. If not
183 
184   if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
185     { // only for testing now
186       gammaPolarization0 = GetRandomPolarization(gammaDirection0);
187     }
188   else
189     {
190       if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
191   {
192     gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
193   }
194     }
195 
196   // End of Protection
197 
198   // Within energy limit?
199 
200   if(gammaEnergy0 <= lowEnergyLimit)
201     {
202       aParticleChange.ProposeTrackStatus(fStopAndKill);
203       aParticleChange.ProposeEnergy(0.);
204       aParticleChange.ProposeLocalEnergyDeposit(gammaEnergy0);
205       return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
206     }
207 
208   G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
209 
210   // Select randomly one element in the current material
211 
212   const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
213   G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
214 
215   // Sample the energy and the polarization of the scattered photon
216 
217   G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
218 
219   G4double epsilon0 = 1./(1. + 2*E0_m);
220   G4double epsilon0Sq = epsilon0*epsilon0;
221   G4double alpha1   = - std::log(epsilon0);
222   G4double alpha2 = 0.5*(1.- epsilon0Sq);
223 
224   G4double wlGamma = h_Planck*c_light/gammaEnergy0;
225   G4double gammaEnergy1;
226   G4ThreeVector gammaDirection1;
227 
228   do {
229     if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
230       {
231   epsilon   = std::exp(-alpha1*G4UniformRand());  
232   epsilonSq = epsilon*epsilon; 
233       }
234     else 
235       {
236   epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand();
237   epsilon   = std::sqrt(epsilonSq);
238       }
239 
240     onecost = (1.- epsilon)/(epsilon*E0_m);
241     sinThetaSqr   = onecost*(2.-onecost);
242 
243     // Protection
244     if (sinThetaSqr > 1.)
245       {
246   if (verboseLevel>0) G4cout
247     << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
248     << "sin(theta)**2 = "
249     << sinThetaSqr
250     << "; set to 1"
251     << G4endl;
252   sinThetaSqr = 1.;
253       }
254     if (sinThetaSqr < 0.)
255       {
256   if (verboseLevel>0) G4cout
257     << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
258     << "sin(theta)**2 = "
259     << sinThetaSqr
260     << "; set to 0"
261     << G4endl;
262   sinThetaSqr = 0.;
263       }
264     // End protection
265 
266     G4double x =  std::sqrt(onecost/2.) / (wlGamma/cm);;
267     G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
268     greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
269 
270   } while(greject < G4UniformRand()*Z);
271 
272 
273   // ****************************************************
274   //    Phi determination
275   // ****************************************************
276 
277   G4double phi = SetPhi(epsilon,sinThetaSqr);
278 
279   //
280   // scattered gamma angles. ( Z - axis along the parent gamma)
281   //
282 
283   G4double cosTheta = 1. - onecost;
284 
285   // Protection
286 
287   if (cosTheta > 1.)
288     {
289       if (verboseLevel>0) G4cout
290   << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
291   << "cosTheta = "
292   << cosTheta
293   << "; set to 1"
294   << G4endl;
295       cosTheta = 1.;
296     }
297   if (cosTheta < -1.)
298     {
299       if (verboseLevel>0) G4cout 
300   << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
301   << "cosTheta = " 
302   << cosTheta
303   << "; set to -1"
304   << G4endl;
305       cosTheta = -1.;
306     }
307   // End protection      
308   
309   
310   G4double sinTheta = std::sqrt (sinThetaSqr);
311   
312   // Protection
313   if (sinTheta > 1.)
314     {
315       if (verboseLevel>0) G4cout 
316   << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
317   << "sinTheta = " 
318   << sinTheta
319   << "; set to 1"
320   << G4endl;
321       sinTheta = 1.;
322     }
323   if (sinTheta < -1.)
324     {
325       if (verboseLevel>0) G4cout 
326   << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
327   << "sinTheta = " 
328   << sinTheta
329   << "; set to -1" 
330   << G4endl;
331       sinTheta = -1.;
332     }
333   // End protection
334   
335       
336   G4double dirx = sinTheta*std::cos(phi);
337   G4double diry = sinTheta*std::sin(phi);
338   G4double dirz = cosTheta ;
339   
340 
341   // oneCosT , eom
342 
343 
344 
345   // Doppler broadening -  Method based on:
346   // Y. Namito, S. Ban and H. Hirayama, 
347   // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code" 
348   // NIM A 349, pp. 489-494, 1994
349   
350   // Maximum number of sampling iterations
351 
352   G4int maxDopplerIterations = 1000;
353   G4double bindingE = 0.;
354   G4double photonEoriginal = epsilon * gammaEnergy0;
355   G4double photonE = -1.;
356   G4int iteration = 0;
357   G4double eMax = gammaEnergy0;
358 
359   do
360     {
361       iteration++;
362       // Select shell based on shell occupancy
363       G4int shell = shellData.SelectRandomShell(Z);
364       bindingE = shellData.BindingEnergy(Z,shell);
365       
366       eMax = gammaEnergy0 - bindingE;
367      
368       // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
369       G4double pSample = profileData.RandomSelectMomentum(Z,shell);
370       // Rescale from atomic units
371       G4double pDoppler = pSample * fine_structure_const;
372       G4double pDoppler2 = pDoppler * pDoppler;
373       G4double var2 = 1. + onecost * E0_m;
374       G4double var3 = var2*var2 - pDoppler2;
375       G4double var4 = var2 - pDoppler2 * cosTheta;
376       G4double var = var4*var4 - var3 + pDoppler2 * var3;
377       if (var > 0.)
378   {
379     G4double varSqrt = std::sqrt(var);        
380     G4double scale = gammaEnergy0 / var3;  
381           // Random select either root
382     if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;               
383     else photonE = (var4 + varSqrt) * scale;
384   } 
385       else
386   {
387     photonE = -1.;
388   }
389    } while ( iteration <= maxDopplerIterations && 
390        (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
391  
392   // End of recalculation of photon energy with Doppler broadening
393   // Revert to original if maximum number of iterations threshold has been reached
394   if (iteration >= maxDopplerIterations)
395     {
396       photonE = photonEoriginal;
397       bindingE = 0.;
398     }
399 
400   gammaEnergy1 = photonE;
401  
402   // G4cout << "--> PHOTONENERGY1 = " << photonE/keV << G4endl;
403 
404 
405   /// Doppler Broadeing 
406 
407 
408 
409 
410   //
411   // update G4VParticleChange for the scattered photon 
412   //
413 
414   //  gammaEnergy1 = epsilon*gammaEnergy0;
415 
416 
417   // New polarization
418 
419   G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
420               sinThetaSqr,
421               phi,
422               cosTheta);
423 
424   // Set new direction
425   G4ThreeVector tmpDirection1( dirx,diry,dirz );
426   gammaDirection1 = tmpDirection1;
427 
428   // Change reference frame.
429 
430   SystemOfRefChange(gammaDirection0,gammaDirection1,
431         gammaPolarization0,gammaPolarization1);
432 
433   if (gammaEnergy1 > 0.)
434     {
435       aParticleChange.ProposeEnergy( gammaEnergy1 ) ;
436       aParticleChange.ProposeMomentumDirection( gammaDirection1 );
437       aParticleChange.ProposePolarization( gammaPolarization1 );
438     }
439   else
440     {
441       aParticleChange.ProposeEnergy(0.) ;
442       aParticleChange.ProposeTrackStatus(fStopAndKill);
443     }
444 
445   //
446   // kinematic of the scattered electron
447   //
448 
449   G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
450 
451 
452   // Generate the electron only if with large enough range w.r.t. cuts and safety
453 
454   G4double safety = aStep.GetPostStepPoint()->GetSafety();
455 
456 
457   if (rangeTest->Escape(G4Electron::Electron(),couple,ElecKineEnergy,safety))
458     {
459       G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
460       G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
461            gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
462       G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
463       aParticleChange.SetNumberOfSecondaries(1);
464       aParticleChange.AddSecondary(electron);
465       //      aParticleChange.ProposeLocalEnergyDeposit(0.); 
466       aParticleChange.ProposeLocalEnergyDeposit(bindingE);
467     }
468   else
469     {
470       aParticleChange.SetNumberOfSecondaries(0);
471       aParticleChange.ProposeLocalEnergyDeposit(ElecKineEnergy+bindingE);
472     }
473   
474   return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep);
475   
476 }
477 
478 
479 G4double G4LowEnergyPolarizedCompton::SetPhi(G4double energyRate,
480                G4double sinSqrTh)
481 {
482   G4double rand1;
483   G4double rand2;
484   G4double phiProbability;
485   G4double phi;
486   G4double a, b;
487 
488   do
489     {
490       rand1 = G4UniformRand();
491       rand2 = G4UniformRand();
492       phiProbability=0.;
493       phi = twopi*rand1;
494       
495       a = 2*sinSqrTh;
496       b = energyRate + 1/energyRate;
497       
498       phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
499 
500       
501  
502     }
503   while ( rand2 > phiProbability );
504   return phi;
505 }
506 
507 
508 G4ThreeVector G4LowEnergyPolarizedCompton::SetPerpendicularVector(G4ThreeVector& a)
509 {
510   G4double dx = a.x();
511   G4double dy = a.y();
512   G4double dz = a.z();
513   G4double x = dx < 0.0 ? -dx : dx;
514   G4double y = dy < 0.0 ? -dy : dy;
515   G4double z = dz < 0.0 ? -dz : dz;
516   if (x < y) {
517     return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
518   }else{
519     return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
520   }
521 }
522 
523 G4ThreeVector G4LowEnergyPolarizedCompton::GetRandomPolarization(G4ThreeVector& direction0)
524 {
525   G4ThreeVector d0 = direction0.unit();
526   G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
527   G4ThreeVector a0 = a1.unit(); // unit vector
528 
529   G4double rand1 = G4UniformRand();
530   
531   G4double angle = twopi*rand1; // random polar angle
532   G4ThreeVector b0 = d0.cross(a0); // cross product
533   
534   G4ThreeVector c;
535   
536   c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
537   c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
538   c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
539   
540   G4ThreeVector c0 = c.unit();
541 
542   return c0;
543   
544 }
545 
546 
547 G4ThreeVector G4LowEnergyPolarizedCompton::GetPerpendicularPolarization
548 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
549 {
550 
551   // 
552   // The polarization of a photon is always perpendicular to its momentum direction.
553   // Therefore this function removes those vector component of gammaPolarization, which
554   // points in direction of gammaDirection
555   //
556   // Mathematically we search the projection of the vector a on the plane E, where n is the
557   // plains normal vector.
558   // The basic equation can be found in each geometry book (e.g. Bronstein):
559   // p = a - (a o n)/(n o n)*n
560   
561   return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
562 }
563 
564 
565 G4ThreeVector G4LowEnergyPolarizedCompton::SetNewPolarization(G4double epsilon,
566                     G4double sinSqrTh, 
567                     G4double phi,
568                     G4double costheta) 
569 {
570   G4double rand1;
571   G4double rand2;
572   G4double cosPhi = std::cos(phi);
573   G4double sinPhi = std::sin(phi);
574   G4double sinTheta = std::sqrt(sinSqrTh);
575   G4double cosSqrPhi = cosPhi*cosPhi;
576   //  G4double cossqrth = 1.-sinSqrTh;
577   //  G4double sinsqrphi = sinPhi*sinPhi;
578   G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
579  
580 
581   // Determination of Theta 
582   
583   // ---- MGP ---- Commented out the following 3 lines to avoid compilation 
584   // warnings (unused variables)
585   // G4double thetaProbability;
586   G4double theta;
587   // G4double a, b;
588   // G4double cosTheta;
589 
590   /*
591 
592   depaola method
593   
594   do
595   {
596       rand1 = G4UniformRand();
597       rand2 = G4UniformRand();
598       thetaProbability=0.;
599       theta = twopi*rand1;
600       a = 4*normalisation*normalisation;
601       b = (epsilon + 1/epsilon) - 2;
602       thetaProbability = (b + a*std::cos(theta)*std::cos(theta))/(a+b);
603       cosTheta = std::cos(theta);
604     }
605   while ( rand2 > thetaProbability );
606   
607   G4double cosBeta = cosTheta;
608 
609   */
610 
611 
612   // Dan Xu method (IEEE TNS, 52, 1160 (2005))
613 
614   rand1 = G4UniformRand();
615   rand2 = G4UniformRand();
616 
617   if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
618     {
619       if (rand2<0.5)
620   theta = pi/2.0;
621       else
622   theta = 3.0*pi/2.0;
623     }
624   else
625     {
626       if (rand2<0.5)
627   theta = 0;
628       else
629   theta = pi;
630     }
631   G4double cosBeta = std::cos(theta);
632   G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
633   
634   G4ThreeVector gammaPolarization1;
635 
636   G4double xParallel = normalisation*cosBeta;
637   G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
638   G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
639   G4double xPerpendicular = 0.;
640   G4double yPerpendicular = (costheta)*sinBeta/normalisation;
641   G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
642 
643   G4double xTotal = (xParallel + xPerpendicular);
644   G4double yTotal = (yParallel + yPerpendicular);
645   G4double zTotal = (zParallel + zPerpendicular);
646   
647   gammaPolarization1.setX(xTotal);
648   gammaPolarization1.setY(yTotal);
649   gammaPolarization1.setZ(zTotal);
650   
651   return gammaPolarization1;
652 
653 }
654 
655 
656 void G4LowEnergyPolarizedCompton::SystemOfRefChange(G4ThreeVector& direction0,
657                 G4ThreeVector& direction1,
658                 G4ThreeVector& polarization0,
659                 G4ThreeVector& polarization1)
660 {
661   // direction0 is the original photon direction ---> z
662   // polarization0 is the original photon polarization ---> x
663   // need to specify y axis in the real reference frame ---> y 
664   G4ThreeVector Axis_Z0 = direction0.unit();
665   G4ThreeVector Axis_X0 = polarization0.unit();
666   G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
667 
668   G4double direction_x = direction1.getX();
669   G4double direction_y = direction1.getY();
670   G4double direction_z = direction1.getZ();
671   
672   direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
673   G4double polarization_x = polarization1.getX();
674   G4double polarization_y = polarization1.getY();
675   G4double polarization_z = polarization1.getZ();
676 
677   polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
678 
679 }
680 
681 
682 G4bool G4LowEnergyPolarizedCompton::IsApplicable(const G4ParticleDefinition& particle)
683 {
684   return ( &particle == G4Gamma::Gamma() ); 
685 }
686 
687 
688 G4double G4LowEnergyPolarizedCompton::GetMeanFreePath(const G4Track& track,
689                   G4double,
690                   G4ForceCondition*)
691 {
692   const G4DynamicParticle* photon = track.GetDynamicParticle();
693   G4double energy = photon->GetKineticEnergy();
694   const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
695   size_t materialIndex = couple->GetIndex();
696   G4double meanFreePath;
697   if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
698   else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
699   else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
700   return meanFreePath;
701 }
702 
703 
704 
705 
706 
707 
708 
709 
710 
711 
712 
713 
714 
715 
716