physics/src/G4LowEnergyPolarizedCompton.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 //
0028 // ------------------------------------------------------------
0029 //      GEANT 4 class implementation file
0030 //      CERN Geneva Switzerland
0031 //
0032
0033 // --------- G4LowEnergyPolarizedCompton class -----
0034 //
0035 //           by G.Depaola & F.Longo (21 may 2001)
0036 //
0037 // 21 May 2001 - MGP      Modified to inherit from G4VDiscreteProcess
0038 //                        Applies same algorithm as LowEnergyCompton
0039 //                        if the incoming photon is not polarised
0040 //                        Temporary protection to avoid crash in the case
0041 //                        of polarisation || incident photon direction
0042 //
0043 // 17 October 2001 - F.Longo - Revised according to a design iteration
0044 //
0045 // 21 February 2002 - F.Longo Revisions with A.Zoglauer and G.Depaola
0046 //                            - better description of parallelism
0047 //                            - system of ref change method improved
0048 // 22 January  2003 - V.Ivanchenko Cut per region
0049 // 24 April    2003 - V.Ivanchenko Cut per region mfpt
0050 //
0051 //
0052 // ************************************************************
0053 //
0054 // Corrections by Rui Curado da Silva (2000)
0055 // New Implementation by G.Depaola & F.Longo
0056 //
0057 // - sampling of phi
0058 // - polarization of scattered photon
0059 //
0060 // --------------------------------------------------------------
0061
0062 #include "G4LowEnergyPolarizedCompton.hh"
0063 #include "Randomize.hh"
0064 #include "G4PhysicalConstants.hh"
0065 #include "G4SystemOfUnits.hh"
0066 #include "G4ParticleDefinition.hh"
0067 #include "G4Track.hh"
0068 #include "G4Step.hh"
0069 #include "G4ForceCondition.hh"
0070 #include "G4Gamma.hh"
0071 #include "G4Electron.hh"
0072 #include "G4DynamicParticle.hh"
0073 #include "G4VParticleChange.hh"
0074 #include "G4ThreeVector.hh"
0075 #include "G4RDVCrossSectionHandler.hh"
0076 #include "G4RDCrossSectionHandler.hh"
0077 #include "G4RDVEMDataSet.hh"
0078 #include "G4RDCompositeEMDataSet.hh"
0079 #include "G4RDVDataSetAlgorithm.hh"
0080 #include "G4RDLogLogInterpolation.hh"
0081 #include "G4RDVRangeTest.hh"
0082 #include "G4RDRangeTest.hh"
0083 #include "G4MaterialCutsCouple.hh"
0084
0085 // constructor
0086
0087 G4LowEnergyPolarizedCompton::G4LowEnergyPolarizedCompton(const G4String& processName)
0088   : G4VDiscreteProcess(processName),
0089     lowEnergyLimit (250*eV),              // initialization
0090     highEnergyLimit(100*GeV),
0091     intrinsicLowEnergyLimit(10*eV),
0092     intrinsicHighEnergyLimit(100*GeV)
0093 {
0094   if (lowEnergyLimit < intrinsicLowEnergyLimit ||
0095       highEnergyLimit > intrinsicHighEnergyLimit)
0096     {
0097       G4Exception("G4LowEnergyPolarizedCompton::G4LowEnergyPolarizedCompton()",
0098                   "OutOfRange", FatalException,
0099                   "Energy outside intrinsic process validity range!");
0100     }
0101
0102   crossSectionHandler = new G4RDCrossSectionHandler;
0103
0104
0105   G4RDVDataSetAlgorithm* scatterInterpolation = new G4RDLogLogInterpolation;
0106   G4String scatterFile = "comp/ce-sf-";
0107   scatterFunctionData = new G4RDCompositeEMDataSet(scatterInterpolation,1.,1.);
0108   scatterFunctionData->LoadData(scatterFile);
0109
0110   meanFreePathTable = 0;
0111
0112   rangeTest = new G4RDRangeTest;
0113
0114   // For Doppler broadening
0115   shellData.SetOccupancyData();
0116
0117
0118    if (verboseLevel > 0)
0119      {
0120        G4cout << GetProcessName() << " is created " << G4endl
0121           << "Energy range: "
0122           << lowEnergyLimit / keV << " keV - "
0123           << highEnergyLimit / GeV << " GeV"
0124           << G4endl;
0125      }
0126 }
0127
0128 // destructor
0129
0130 G4LowEnergyPolarizedCompton::~G4LowEnergyPolarizedCompton()
0131 {
0132   delete meanFreePathTable;
0133   delete crossSectionHandler;
0134   delete scatterFunctionData;
0135   delete rangeTest;
0136 }
0137
0138
0139 void G4LowEnergyPolarizedCompton::BuildPhysicsTable(const G4ParticleDefinition& )
0140 {
0141
0142   crossSectionHandler->Clear();
0143   G4String crossSectionFile = "comp/ce-cs-";
0144   crossSectionHandler->LoadData(crossSectionFile);
0145   delete meanFreePathTable;
0146   meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
0147
0148   // For Doppler broadening
0149   G4String file = "/doppler/shell-doppler";
0150   shellData.LoadData(file);
0151
0152 }
0153
0154 G4VParticleChange* G4LowEnergyPolarizedCompton::PostStepDoIt(const G4Track& aTrack,
0155                                  const G4Step&  aStep)
0156 {
0157   // The scattered gamma energy is sampled according to Klein - Nishina formula.
0158   // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15).
0159   // GEANT4 internal units
0160   //
0161   // Note : Effects due to binding of atomic electrons are negliged.
0162
0163   aParticleChange.Initialize(aTrack);
0164
0165   // Dynamic particle quantities
0166   const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
0167   G4double gammaEnergy0 = incidentPhoton->GetKineticEnergy();
0168   G4ThreeVector gammaPolarization0 = incidentPhoton->GetPolarization();
0169
0170   //  gammaPolarization0 = gammaPolarization0.unit(); //
0171
0172   // Protection: a polarisation parallel to the
0173   // direction causes problems;
0174   // in that case find a random polarization
0175
0176   G4ThreeVector gammaDirection0 = incidentPhoton->GetMomentumDirection();
0177   // ---- MGP ---- Next two lines commented out to remove compilation warnings
0178   // G4double scalarproduct = gammaPolarization0.dot(gammaDirection0);
0179   // G4double angle = gammaPolarization0.angle(gammaDirection0);
0180
0181   // Make sure that the polarization vector is perpendicular to the
0182   // gamma direction. If not
0183
0184   if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
0185     { // only for testing now
0186       gammaPolarization0 = GetRandomPolarization(gammaDirection0);
0187     }
0188   else
0189     {
0190       if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
0191     {
0192       gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
0193     }
0194     }
0195
0196   // End of Protection
0197
0198   // Within energy limit?
0199
0200   if(gammaEnergy0 <= lowEnergyLimit)
0201     {
0202       aParticleChange.ProposeTrackStatus(fStopAndKill);
0203       aParticleChange.ProposeEnergy(0.);
0204       aParticleChange.ProposeLocalEnergyDeposit(gammaEnergy0);
0205       return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
0206     }
0207
0208   G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
0209
0210   // Select randomly one element in the current material
0211
0212   const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
0213   G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
0214
0215   // Sample the energy and the polarization of the scattered photon
0216
0217   G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
0218
0219   G4double epsilon0 = 1./(1. + 2*E0_m);
0220   G4double epsilon0Sq = epsilon0*epsilon0;
0221   G4double alpha1   = - std::log(epsilon0);
0222   G4double alpha2 = 0.5*(1.- epsilon0Sq);
0223
0224   G4double wlGamma = h_Planck*c_light/gammaEnergy0;
0225   G4double gammaEnergy1;
0226   G4ThreeVector gammaDirection1;
0227
0228   do {
0229     if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
0230       {
0231     epsilon   = std::exp(-alpha1*G4UniformRand());
0232     epsilonSq = epsilon*epsilon;
0233       }
0234     else
0235       {
0236     epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand();
0237     epsilon   = std::sqrt(epsilonSq);
0238       }
0239
0240     onecost = (1.- epsilon)/(epsilon*E0_m);
0241     sinThetaSqr   = onecost*(2.-onecost);
0242
0243     // Protection
0244     if (sinThetaSqr > 1.)
0245       {
0246     if (verboseLevel>0) G4cout
0247       << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
0248       << "sin(theta)**2 = "
0249       << sinThetaSqr
0250       << "; set to 1"
0251       << G4endl;
0252     sinThetaSqr = 1.;
0253       }
0254     if (sinThetaSqr < 0.)
0255       {
0256     if (verboseLevel>0) G4cout
0257       << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
0258       << "sin(theta)**2 = "
0259       << sinThetaSqr
0260       << "; set to 0"
0261       << G4endl;
0262     sinThetaSqr = 0.;
0263       }
0264     // End protection
0265
0266     G4double x =  std::sqrt(onecost/2.) / (wlGamma/cm);;
0267     G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
0268     greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
0269
0270   } while(greject < G4UniformRand()*Z);
0271
0272
0273   // ****************************************************
0274   //        Phi determination
0275   // ****************************************************
0276
0277   G4double phi = SetPhi(epsilon,sinThetaSqr);
0278
0279   //
0280   // scattered gamma angles. ( Z - axis along the parent gamma)
0281   //
0282
0283   G4double cosTheta = 1. - onecost;
0284
0285   // Protection
0286
0287   if (cosTheta > 1.)
0288     {
0289       if (verboseLevel>0) G4cout
0290     << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
0291     << "cosTheta = "
0292     << cosTheta
0293     << "; set to 1"
0294     << G4endl;
0295       cosTheta = 1.;
0296     }
0297   if (cosTheta < -1.)
0298     {
0299       if (verboseLevel>0) G4cout
0300     << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
0301     << "cosTheta = "
0302     << cosTheta
0303     << "; set to -1"
0304     << G4endl;
0305       cosTheta = -1.;
0306     }
0307   // End protection
0308
0309
0310   G4double sinTheta = std::sqrt (sinThetaSqr);
0311
0312   // Protection
0313   if (sinTheta > 1.)
0314     {
0315       if (verboseLevel>0) G4cout
0316     << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
0317     << "sinTheta = "
0318     << sinTheta
0319     << "; set to 1"
0320     << G4endl;
0321       sinTheta = 1.;
0322     }
0323   if (sinTheta < -1.)
0324     {
0325       if (verboseLevel>0) G4cout
0326     << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
0327     << "sinTheta = "
0328     << sinTheta
0329     << "; set to -1"
0330     << G4endl;
0331       sinTheta = -1.;
0332     }
0333   // End protection
0334
0335
0336   G4double dirx = sinTheta*std::cos(phi);
0337   G4double diry = sinTheta*std::sin(phi);
0338   G4double dirz = cosTheta ;
0339
0340
0341   // oneCosT , eom
0342
0343
0344
0345   // Doppler broadening -  Method based on:
0346   // Y. Namito, S. Ban and H. Hirayama,
0347   // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code"
0348   // NIM A 349, pp. 489-494, 1994
0349
0350   // Maximum number of sampling iterations
0351
0352   G4int maxDopplerIterations = 1000;
0353   G4double bindingE = 0.;
0354   G4double photonEoriginal = epsilon * gammaEnergy0;
0355   G4double photonE = -1.;
0356   G4int iteration = 0;
0357   G4double eMax = gammaEnergy0;
0358
0359   do
0360     {
0361       iteration++;
0362       // Select shell based on shell occupancy
0363       G4int shell = shellData.SelectRandomShell(Z);
0364       bindingE = shellData.BindingEnergy(Z,shell);
0365
0366       eMax = gammaEnergy0 - bindingE;
0367
0368       // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
0369       G4double pSample = profileData.RandomSelectMomentum(Z,shell);
0370       // Rescale from atomic units
0371       G4double pDoppler = pSample * fine_structure_const;
0372       G4double pDoppler2 = pDoppler * pDoppler;
0373       G4double var2 = 1. + onecost * E0_m;
0374       G4double var3 = var2*var2 - pDoppler2;
0375       G4double var4 = var2 - pDoppler2 * cosTheta;
0376       G4double var = var4*var4 - var3 + pDoppler2 * var3;
0377       if (var > 0.)
0378     {
0379       G4double varSqrt = std::sqrt(var);
0380       G4double scale = gammaEnergy0 / var3;
0381           // Random select either root
0382       if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
0383       else photonE = (var4 + varSqrt) * scale;
0384     }
0385       else
0386     {
0387       photonE = -1.;
0388     }
0389    } while ( iteration <= maxDopplerIterations &&
0390          (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
0391
0392   // End of recalculation of photon energy with Doppler broadening
0393   // Revert to original if maximum number of iterations threshold has been reached
0394   if (iteration >= maxDopplerIterations)
0395     {
0396       photonE = photonEoriginal;
0397       bindingE = 0.;
0398     }
0399
0400   gammaEnergy1 = photonE;
0401
0402   // G4cout << "--> PHOTONENERGY1 = " << photonE/keV << G4endl;
0403
0404
0405   /// Doppler Broadeing
0406
0407
0408
0409
0410   //
0411   // update G4VParticleChange for the scattered photon
0412   //
0413
0414   //  gammaEnergy1 = epsilon*gammaEnergy0;
0415
0416
0417   // New polarization
0418
0419   G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
0420                             sinThetaSqr,
0421                             phi,
0422                             cosTheta);
0423
0424   // Set new direction
0425   G4ThreeVector tmpDirection1( dirx,diry,dirz );
0426   gammaDirection1 = tmpDirection1;
0427
0428   // Change reference frame.
0429
0430   SystemOfRefChange(gammaDirection0,gammaDirection1,
0431             gammaPolarization0,gammaPolarization1);
0432
0433   if (gammaEnergy1 > 0.)
0434     {
0435       aParticleChange.ProposeEnergy( gammaEnergy1 ) ;
0436       aParticleChange.ProposeMomentumDirection( gammaDirection1 );
0437       aParticleChange.ProposePolarization( gammaPolarization1 );
0438     }
0439   else
0440     {
0441       aParticleChange.ProposeEnergy(0.) ;
0442       aParticleChange.ProposeTrackStatus(fStopAndKill);
0443     }
0444
0445   //
0446   // kinematic of the scattered electron
0447   //
0448
0449   G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
0450
0451
0452   // Generate the electron only if with large enough range w.r.t. cuts and safety
0453
0454   G4double safety = aStep.GetPostStepPoint()->GetSafety();
0455
0456
0457   if (rangeTest->Escape(G4Electron::Electron(),couple,ElecKineEnergy,safety))
0458     {
0459       G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
0460       G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
0461                    gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
0462       G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
0463       aParticleChange.SetNumberOfSecondaries(1);
0464       aParticleChange.AddSecondary(electron);
0465       //      aParticleChange.ProposeLocalEnergyDeposit(0.);
0466       aParticleChange.ProposeLocalEnergyDeposit(bindingE);
0467     }
0468   else
0469     {
0470       aParticleChange.SetNumberOfSecondaries(0);
0471       aParticleChange.ProposeLocalEnergyDeposit(ElecKineEnergy+bindingE);
0472     }
0473
0474   return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep);
0475
0476 }
0477
0478
0479 G4double G4LowEnergyPolarizedCompton::SetPhi(G4double energyRate,
0480                          G4double sinSqrTh)
0481 {
0482   G4double rand1;
0483   G4double rand2;
0484   G4double phiProbability;
0485   G4double phi;
0486   G4double a, b;
0487
0488   do
0489     {
0490       rand1 = G4UniformRand();
0491       rand2 = G4UniformRand();
0492       phiProbability=0.;
0493       phi = twopi*rand1;
0494
0495       a = 2*sinSqrTh;
0496       b = energyRate + 1/energyRate;
0497
0498       phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
0499
0500
0501
0502     }
0503   while ( rand2 > phiProbability );
0504   return phi;
0505 }
0506
0507
0508 G4ThreeVector G4LowEnergyPolarizedCompton::SetPerpendicularVector(G4ThreeVector& a)
0509 {
0510   G4double dx = a.x();
0511   G4double dy = a.y();
0512   G4double dz = a.z();
0513   G4double x = dx < 0.0 ? -dx : dx;
0514   G4double y = dy < 0.0 ? -dy : dy;
0515   G4double z = dz < 0.0 ? -dz : dz;
0516   if (x < y) {
0517     return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
0518   }else{
0519     return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
0520   }
0521 }
0522
0523 G4ThreeVector G4LowEnergyPolarizedCompton::GetRandomPolarization(G4ThreeVector& direction0)
0524 {
0525   G4ThreeVector d0 = direction0.unit();
0526   G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
0527   G4ThreeVector a0 = a1.unit(); // unit vector
0528
0529   G4double rand1 = G4UniformRand();
0530
0531   G4double angle = twopi*rand1; // random polar angle
0532   G4ThreeVector b0 = d0.cross(a0); // cross product
0533
0534   G4ThreeVector c;
0535
0536   c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
0537   c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
0538   c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
0539
0540   G4ThreeVector c0 = c.unit();
0541
0542   return c0;
0543
0544 }
0545
0546
0547 G4ThreeVector G4LowEnergyPolarizedCompton::GetPerpendicularPolarization
0548 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
0549 {
0550
0551   //
0552   // The polarization of a photon is always perpendicular to its momentum direction.
0553   // Therefore this function removes those vector component of gammaPolarization, which
0554   // points in direction of gammaDirection
0555   //
0556   // Mathematically we search the projection of the vector a on the plane E, where n is the
0557   // plains normal vector.
0558   // The basic equation can be found in each geometry book (e.g. Bronstein):
0559   // p = a - (a o n)/(n o n)*n
0560
0561   return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
0562 }
0563
0564
0565 G4ThreeVector G4LowEnergyPolarizedCompton::SetNewPolarization(G4double epsilon,
0566                                   G4double sinSqrTh,
0567                                   G4double phi,
0568                                   G4double costheta)
0569 {
0570   G4double rand1;
0571   G4double rand2;
0572   G4double cosPhi = std::cos(phi);
0573   G4double sinPhi = std::sin(phi);
0574   G4double sinTheta = std::sqrt(sinSqrTh);
0575   G4double cosSqrPhi = cosPhi*cosPhi;
0576   //  G4double cossqrth = 1.-sinSqrTh;
0577   //  G4double sinsqrphi = sinPhi*sinPhi;
0578   G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
0579
0580
0581   // Determination of Theta
0582
0583   // ---- MGP ---- Commented out the following 3 lines to avoid compilation
0584   // warnings (unused variables)
0585   // G4double thetaProbability;
0586   G4double theta;
0587   // G4double a, b;
0588   // G4double cosTheta;
0589
0590   /*
0591
0592   depaola method
0593
0594   do
0595   {
0596       rand1 = G4UniformRand();
0597       rand2 = G4UniformRand();
0598       thetaProbability=0.;
0599       theta = twopi*rand1;
0600       a = 4*normalisation*normalisation;
0601       b = (epsilon + 1/epsilon) - 2;
0602       thetaProbability = (b + a*std::cos(theta)*std::cos(theta))/(a+b);
0603       cosTheta = std::cos(theta);
0604     }
0605   while ( rand2 > thetaProbability );
0606
0607   G4double cosBeta = cosTheta;
0608
0609   */
0610
0611
0612   // Dan Xu method (IEEE TNS, 52, 1160 (2005))
0613
0614   rand1 = G4UniformRand();
0615   rand2 = G4UniformRand();
0616
0617   if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
0618     {
0619       if (rand2<0.5)
0620     theta = pi/2.0;
0621       else
0622     theta = 3.0*pi/2.0;
0623     }
0624   else
0625     {
0626       if (rand2<0.5)
0627     theta = 0;
0628       else
0629     theta = pi;
0630     }
0631   G4double cosBeta = std::cos(theta);
0632   G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
0633
0634   G4ThreeVector gammaPolarization1;
0635
0636   G4double xParallel = normalisation*cosBeta;
0637   G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
0638   G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
0639   G4double xPerpendicular = 0.;
0640   G4double yPerpendicular = (costheta)*sinBeta/normalisation;
0641   G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
0642
0643   G4double xTotal = (xParallel + xPerpendicular);
0644   G4double yTotal = (yParallel + yPerpendicular);
0645   G4double zTotal = (zParallel + zPerpendicular);
0646
0647   gammaPolarization1.setX(xTotal);
0648   gammaPolarization1.setY(yTotal);
0649   gammaPolarization1.setZ(zTotal);
0650
0651   return gammaPolarization1;
0652
0653 }
0654
0655
0656 void G4LowEnergyPolarizedCompton::SystemOfRefChange(G4ThreeVector& direction0,
0657                             G4ThreeVector& direction1,
0658                             G4ThreeVector& polarization0,
0659                             G4ThreeVector& polarization1)
0660 {
0661   // direction0 is the original photon direction ---> z
0662   // polarization0 is the original photon polarization ---> x
0663   // need to specify y axis in the real reference frame ---> y
0664   G4ThreeVector Axis_Z0 = direction0.unit();
0665   G4ThreeVector Axis_X0 = polarization0.unit();
0666   G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
0667
0668   G4double direction_x = direction1.getX();
0669   G4double direction_y = direction1.getY();
0670   G4double direction_z = direction1.getZ();
0671
0672   direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
0673   G4double polarization_x = polarization1.getX();
0674   G4double polarization_y = polarization1.getY();
0675   G4double polarization_z = polarization1.getZ();
0676
0677   polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
0678
0679 }
0680
0681
0682 G4bool G4LowEnergyPolarizedCompton::IsApplicable(const G4ParticleDefinition& particle)
0683 {
0684   return ( &particle == G4Gamma::Gamma() );
0685 }
0686
0687
0688 G4double G4LowEnergyPolarizedCompton::GetMeanFreePath(const G4Track& track,
0689                               G4double,
0690                               G4ForceCondition*)
0691 {
0692   const G4DynamicParticle* photon = track.GetDynamicParticle();
0693   G4double energy = photon->GetKineticEnergy();
0694   const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
0695   size_t materialIndex = couple->GetIndex();
0696   G4double meanFreePath;
0697   if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
0698   else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
0699   else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
0700   return meanFreePath;
0701 }
0702
0703
0704
0705
0706
0707
0708
0709
0710
0711
0712
0713
0714
0715