G4LivermorePolarizedComptonModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // Authors: G.Depaola & F.Longo
00029 //
00030 // History:
00031 // --------
00032 // 02 May 2009   S Incerti as V. Ivanchenko proposed in G4LivermoreComptonModel.cc
00033 //
00034 // Cleanup initialisation and generation of secondaries:
00035 //                  - apply internal high-energy limit only in constructor 
00036 //                  - do not apply low-energy limit (default is 0)
00037 //                  - remove GetMeanFreePath method and table
00038 //                  - added protection against numerical problem in energy sampling 
00039 //                  - use G4ElementSelector
00040 
00041 #include "G4LivermorePolarizedComptonModel.hh"
00042 #include "G4PhysicalConstants.hh"
00043 #include "G4SystemOfUnits.hh"
00044 
00045 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00046 
00047 using namespace std;
00048 
00049 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00050 
00051 G4LivermorePolarizedComptonModel::G4LivermorePolarizedComptonModel(const G4ParticleDefinition*,
00052                                              const G4String& nam)
00053   :G4VEmModel(nam),fParticleChange(0),isInitialised(false),
00054    meanFreePathTable(0),scatterFunctionData(0),crossSectionHandler(0)
00055 {
00056   lowEnergyLimit = 250 * eV; 
00057   highEnergyLimit = 100 * GeV;
00058   //SetLowEnergyLimit(lowEnergyLimit);
00059   SetHighEnergyLimit(highEnergyLimit);
00060   
00061   verboseLevel= 0;
00062   // Verbosity scale:
00063   // 0 = nothing 
00064   // 1 = warning for energy non-conservation 
00065   // 2 = details of energy budget
00066   // 3 = calculation of cross sections, file openings, sampling of atoms
00067   // 4 = entering in methods
00068 
00069   if( verboseLevel>0 ) { 
00070   G4cout << "Livermore Polarized Compton is constructed " << G4endl
00071          << "Energy range: "
00072          << lowEnergyLimit / eV << " eV - "
00073          << highEnergyLimit / GeV << " GeV"
00074          << G4endl;
00075   }
00076 }
00077 
00078 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00079 
00080 G4LivermorePolarizedComptonModel::~G4LivermorePolarizedComptonModel()
00081 {  
00082   if (meanFreePathTable)   delete meanFreePathTable;
00083   if (crossSectionHandler) delete crossSectionHandler;
00084   if (scatterFunctionData) delete scatterFunctionData;
00085 }
00086 
00087 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00088 
00089 void G4LivermorePolarizedComptonModel::Initialise(const G4ParticleDefinition* particle,
00090                                        const G4DataVector& cuts)
00091 {
00092   if (verboseLevel > 3)
00093     G4cout << "Calling G4LivermorePolarizedComptonModel::Initialise()" << G4endl;
00094 
00095   if (crossSectionHandler)
00096   {
00097     crossSectionHandler->Clear();
00098     delete crossSectionHandler;
00099   }
00100 
00101   // Reading of data files - all materials are read
00102   
00103   crossSectionHandler = new G4CrossSectionHandler;
00104   crossSectionHandler->Clear();
00105   G4String crossSectionFile = "comp/ce-cs-";
00106   crossSectionHandler->LoadData(crossSectionFile);
00107 
00108   meanFreePathTable = 0;
00109   meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
00110 
00111   G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
00112   G4String scatterFile = "comp/ce-sf-";
00113   scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
00114   scatterFunctionData->LoadData(scatterFile);
00115 
00116   // For Doppler broadening
00117   shellData.SetOccupancyData();
00118   G4String file = "/doppler/shell-doppler";
00119   shellData.LoadData(file);
00120 
00121   if (verboseLevel > 2) 
00122     G4cout << "Loaded cross section files for Livermore Polarized Compton model" << G4endl;
00123 
00124   InitialiseElementSelectors(particle,cuts);
00125 
00126   if(  verboseLevel>0 ) { 
00127     G4cout << "Livermore Polarized Compton model is initialized " << G4endl
00128          << "Energy range: "
00129          << LowEnergyLimit() / eV << " eV - "
00130          << HighEnergyLimit() / GeV << " GeV"
00131          << G4endl;
00132   }
00133   
00134   //
00135     
00136   if(isInitialised) return;
00137   fParticleChange = GetParticleChangeForGamma();
00138   isInitialised = true;
00139 }
00140 
00141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00142 
00143 G4double G4LivermorePolarizedComptonModel::ComputeCrossSectionPerAtom(
00144                                        const G4ParticleDefinition*,
00145                                              G4double GammaEnergy,
00146                                              G4double Z, G4double,
00147                                              G4double, G4double)
00148 {
00149   if (verboseLevel > 3)
00150     G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" << G4endl;
00151 
00152   if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0.0;
00153 
00154   G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
00155   return cs;
00156 }
00157 
00158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00159 
00160 void G4LivermorePolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00161                                               const G4MaterialCutsCouple* couple,
00162                                               const G4DynamicParticle* aDynamicGamma,
00163                                               G4double,
00164                                               G4double)
00165 {
00166   // The scattered gamma energy is sampled according to Klein - Nishina formula.
00167   // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15).
00168   // GEANT4 internal units
00169   //
00170   // Note : Effects due to binding of atomic electrons are negliged.
00171 
00172   if (verboseLevel > 3)
00173     G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" << G4endl;
00174 
00175   G4double gammaEnergy0 = aDynamicGamma->GetKineticEnergy();
00176   G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
00177 
00178   // Protection: a polarisation parallel to the
00179   // direction causes problems;
00180   // in that case find a random polarization
00181 
00182   G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
00183 
00184   // Make sure that the polarization vector is perpendicular to the
00185   // gamma direction. If not
00186 
00187   if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
00188     { // only for testing now
00189       gammaPolarization0 = GetRandomPolarization(gammaDirection0);
00190     }
00191   else
00192     {
00193       if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
00194         {
00195           gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
00196         }
00197     }
00198 
00199   // End of Protection
00200 
00201   // Within energy limit?
00202 
00203   if(gammaEnergy0 <= lowEnergyLimit)
00204     {
00205       fParticleChange->ProposeTrackStatus(fStopAndKill);
00206       fParticleChange->SetProposedKineticEnergy(0.);
00207       fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0);
00208       return;
00209     }
00210 
00211   G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
00212 
00213   // Select randomly one element in the current material
00214   //G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
00215   const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
00216   const G4Element* elm = SelectRandomAtom(couple,particle,gammaEnergy0);
00217   G4int Z = (G4int)elm->GetZ();
00218 
00219   // Sample the energy and the polarization of the scattered photon
00220 
00221   G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
00222 
00223   G4double epsilon0Local = 1./(1. + 2*E0_m);
00224   G4double epsilon0Sq = epsilon0Local*epsilon0Local;
00225   G4double alpha1   = - std::log(epsilon0Local);
00226   G4double alpha2 = 0.5*(1.- epsilon0Sq);
00227 
00228   G4double wlGamma = h_Planck*c_light/gammaEnergy0;
00229   G4double gammaEnergy1;
00230   G4ThreeVector gammaDirection1;
00231 
00232   do {
00233     if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
00234       {
00235         epsilon   = std::exp(-alpha1*G4UniformRand());  
00236         epsilonSq = epsilon*epsilon; 
00237       }
00238     else 
00239       {
00240         epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand();
00241         epsilon   = std::sqrt(epsilonSq);
00242       }
00243 
00244     onecost = (1.- epsilon)/(epsilon*E0_m);
00245     sinThetaSqr   = onecost*(2.-onecost);
00246 
00247     // Protection
00248     if (sinThetaSqr > 1.)
00249       {
00250         G4cout
00251           << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
00252           << "sin(theta)**2 = "
00253           << sinThetaSqr
00254           << "; set to 1"
00255           << G4endl;
00256         sinThetaSqr = 1.;
00257       }
00258     if (sinThetaSqr < 0.)
00259       {
00260         G4cout
00261           << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
00262           << "sin(theta)**2 = "
00263           << sinThetaSqr
00264           << "; set to 0"
00265           << G4endl;
00266         sinThetaSqr = 0.;
00267       }
00268     // End protection
00269 
00270     G4double x =  std::sqrt(onecost/2.) / (wlGamma/cm);;
00271     G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
00272     greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
00273 
00274   } while(greject < G4UniformRand()*Z);
00275 
00276 
00277   // ****************************************************
00278   //            Phi determination
00279   // ****************************************************
00280 
00281   G4double phi = SetPhi(epsilon,sinThetaSqr);
00282 
00283   //
00284   // scattered gamma angles. ( Z - axis along the parent gamma)
00285   //
00286 
00287   G4double cosTheta = 1. - onecost;
00288 
00289   // Protection
00290 
00291   if (cosTheta > 1.)
00292     {
00293       G4cout
00294         << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
00295         << "cosTheta = "
00296         << cosTheta
00297         << "; set to 1"
00298         << G4endl;
00299       cosTheta = 1.;
00300     }
00301   if (cosTheta < -1.)
00302     {
00303       G4cout 
00304         << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
00305         << "cosTheta = " 
00306         << cosTheta
00307         << "; set to -1"
00308         << G4endl;
00309       cosTheta = -1.;
00310     }
00311   // End protection      
00312   
00313   
00314   G4double sinTheta = std::sqrt (sinThetaSqr);
00315   
00316   // Protection
00317   if (sinTheta > 1.)
00318     {
00319       G4cout 
00320         << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
00321         << "sinTheta = " 
00322         << sinTheta
00323         << "; set to 1"
00324         << G4endl;
00325       sinTheta = 1.;
00326     }
00327   if (sinTheta < -1.)
00328     {
00329       G4cout 
00330         << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
00331         << "sinTheta = " 
00332         << sinTheta
00333         << "; set to -1" 
00334         << G4endl;
00335       sinTheta = -1.;
00336     }
00337   // End protection
00338   
00339       
00340   G4double dirx = sinTheta*std::cos(phi);
00341   G4double diry = sinTheta*std::sin(phi);
00342   G4double dirz = cosTheta ;
00343   
00344 
00345   // oneCosT , eom
00346 
00347   // Doppler broadening -  Method based on:
00348   // Y. Namito, S. Ban and H. Hirayama, 
00349   // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code" 
00350   // NIM A 349, pp. 489-494, 1994
00351   
00352   // Maximum number of sampling iterations
00353 
00354   G4int maxDopplerIterations = 1000;
00355   G4double bindingE = 0.;
00356   G4double photonEoriginal = epsilon * gammaEnergy0;
00357   G4double photonE = -1.;
00358   G4int iteration = 0;
00359   G4double eMax = gammaEnergy0;
00360 
00361   do
00362     {
00363       iteration++;
00364       // Select shell based on shell occupancy
00365       G4int shell = shellData.SelectRandomShell(Z);
00366       bindingE = shellData.BindingEnergy(Z,shell);
00367       
00368       eMax = gammaEnergy0 - bindingE;
00369      
00370       // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
00371       G4double pSample = profileData.RandomSelectMomentum(Z,shell);
00372       // Rescale from atomic units
00373       G4double pDoppler = pSample * fine_structure_const;
00374       G4double pDoppler2 = pDoppler * pDoppler;
00375       G4double var2 = 1. + onecost * E0_m;
00376       G4double var3 = var2*var2 - pDoppler2;
00377       G4double var4 = var2 - pDoppler2 * cosTheta;
00378       G4double var = var4*var4 - var3 + pDoppler2 * var3;
00379       if (var > 0.)
00380         {
00381           G4double varSqrt = std::sqrt(var);        
00382           G4double scale = gammaEnergy0 / var3;  
00383           // Random select either root
00384           if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;               
00385           else photonE = (var4 + varSqrt) * scale;
00386         } 
00387       else
00388         {
00389           photonE = -1.;
00390         }
00391    } while ( iteration <= maxDopplerIterations && 
00392              (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
00393  
00394   // End of recalculation of photon energy with Doppler broadening
00395   // Revert to original if maximum number of iterations threshold has been reached
00396   if (iteration >= maxDopplerIterations)
00397     {
00398       photonE = photonEoriginal;
00399       bindingE = 0.;
00400     }
00401 
00402   gammaEnergy1 = photonE;
00403  
00404   //
00405   // update G4VParticleChange for the scattered photon 
00406   //
00407 
00408   //  gammaEnergy1 = epsilon*gammaEnergy0;
00409 
00410 
00411   // New polarization
00412 
00413   G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
00414                                                         sinThetaSqr,
00415                                                         phi,
00416                                                         cosTheta);
00417 
00418   // Set new direction
00419   G4ThreeVector tmpDirection1( dirx,diry,dirz );
00420   gammaDirection1 = tmpDirection1;
00421 
00422   // Change reference frame.
00423 
00424   SystemOfRefChange(gammaDirection0,gammaDirection1,
00425                     gammaPolarization0,gammaPolarization1);
00426 
00427   if (gammaEnergy1 > 0.)
00428     {
00429       fParticleChange->SetProposedKineticEnergy( gammaEnergy1 ) ;
00430       fParticleChange->ProposeMomentumDirection( gammaDirection1 );
00431       fParticleChange->ProposePolarization( gammaPolarization1 );
00432     }
00433   else
00434     {
00435       gammaEnergy1 = 0.;
00436       fParticleChange->SetProposedKineticEnergy(0.) ;
00437       fParticleChange->ProposeTrackStatus(fStopAndKill);
00438     }
00439 
00440   //
00441   // kinematic of the scattered electron
00442   //
00443 
00444   G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
00445 
00446   // SI -protection against negative final energy: no e- is created
00447   // like in G4LivermoreComptonModel.cc
00448   if(ElecKineEnergy < 0.0) {
00449     fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0 - gammaEnergy1);
00450     return;
00451   }
00452  
00453   // SI - Removed range test
00454   
00455   G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
00456 
00457   G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
00458                                    gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
00459 
00460   fParticleChange->ProposeLocalEnergyDeposit(bindingE);
00461   
00462   G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
00463   fvect->push_back(dp);
00464 
00465 }
00466 
00467 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00468 
00469 G4double G4LivermorePolarizedComptonModel::SetPhi(G4double energyRate,
00470                                              G4double sinSqrTh)
00471 {
00472   G4double rand1;
00473   G4double rand2;
00474   G4double phiProbability;
00475   G4double phi;
00476   G4double a, b;
00477 
00478   do
00479     {
00480       rand1 = G4UniformRand();
00481       rand2 = G4UniformRand();
00482       phiProbability=0.;
00483       phi = twopi*rand1;
00484       
00485       a = 2*sinSqrTh;
00486       b = energyRate + 1/energyRate;
00487       
00488       phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
00489 
00490       
00491  
00492     }
00493   while ( rand2 > phiProbability );
00494   return phi;
00495 }
00496 
00497 
00498 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00499 
00500 G4ThreeVector G4LivermorePolarizedComptonModel::SetPerpendicularVector(G4ThreeVector& a)
00501 {
00502   G4double dx = a.x();
00503   G4double dy = a.y();
00504   G4double dz = a.z();
00505   G4double x = dx < 0.0 ? -dx : dx;
00506   G4double y = dy < 0.0 ? -dy : dy;
00507   G4double z = dz < 0.0 ? -dz : dz;
00508   if (x < y) {
00509     return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
00510   }else{
00511     return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
00512   }
00513 }
00514 
00515 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00516 
00517 G4ThreeVector G4LivermorePolarizedComptonModel::GetRandomPolarization(G4ThreeVector& direction0)
00518 {
00519   G4ThreeVector d0 = direction0.unit();
00520   G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
00521   G4ThreeVector a0 = a1.unit(); // unit vector
00522 
00523   G4double rand1 = G4UniformRand();
00524   
00525   G4double angle = twopi*rand1; // random polar angle
00526   G4ThreeVector b0 = d0.cross(a0); // cross product
00527   
00528   G4ThreeVector c;
00529   
00530   c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
00531   c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
00532   c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
00533   
00534   G4ThreeVector c0 = c.unit();
00535 
00536   return c0;
00537   
00538 }
00539 
00540 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00541 
00542 G4ThreeVector G4LivermorePolarizedComptonModel::GetPerpendicularPolarization
00543 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
00544 {
00545 
00546   // 
00547   // The polarization of a photon is always perpendicular to its momentum direction.
00548   // Therefore this function removes those vector component of gammaPolarization, which
00549   // points in direction of gammaDirection
00550   //
00551   // Mathematically we search the projection of the vector a on the plane E, where n is the
00552   // plains normal vector.
00553   // The basic equation can be found in each geometry book (e.g. Bronstein):
00554   // p = a - (a o n)/(n o n)*n
00555   
00556   return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
00557 }
00558 
00559 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00560 
00561 G4ThreeVector G4LivermorePolarizedComptonModel::SetNewPolarization(G4double epsilon,
00562                                                               G4double sinSqrTh, 
00563                                                               G4double phi,
00564                                                               G4double costheta) 
00565 {
00566   G4double rand1;
00567   G4double rand2;
00568   G4double cosPhi = std::cos(phi);
00569   G4double sinPhi = std::sin(phi);
00570   G4double sinTheta = std::sqrt(sinSqrTh);
00571   G4double cosSqrPhi = cosPhi*cosPhi;
00572   //  G4double cossqrth = 1.-sinSqrTh;
00573   //  G4double sinsqrphi = sinPhi*sinPhi;
00574   G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
00575  
00576 
00577   // Determination of Theta 
00578   
00579   // ---- MGP ---- Commented out the following 3 lines to avoid compilation 
00580   // warnings (unused variables)
00581   // G4double thetaProbability;
00582   G4double theta;
00583   // G4double a, b;
00584   // G4double cosTheta;
00585 
00586   /*
00587 
00588   depaola method
00589   
00590   do
00591   {
00592       rand1 = G4UniformRand();
00593       rand2 = G4UniformRand();
00594       thetaProbability=0.;
00595       theta = twopi*rand1;
00596       a = 4*normalisation*normalisation;
00597       b = (epsilon + 1/epsilon) - 2;
00598       thetaProbability = (b + a*std::cos(theta)*std::cos(theta))/(a+b);
00599       cosTheta = std::cos(theta);
00600     }
00601   while ( rand2 > thetaProbability );
00602   
00603   G4double cosBeta = cosTheta;
00604 
00605   */
00606 
00607 
00608   // Dan Xu method (IEEE TNS, 52, 1160 (2005))
00609 
00610   rand1 = G4UniformRand();
00611   rand2 = G4UniformRand();
00612 
00613   if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
00614     {
00615       if (rand2<0.5)
00616         theta = pi/2.0;
00617       else
00618         theta = 3.0*pi/2.0;
00619     }
00620   else
00621     {
00622       if (rand2<0.5)
00623         theta = 0;
00624       else
00625         theta = pi;
00626     }
00627   G4double cosBeta = std::cos(theta);
00628   G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
00629   
00630   G4ThreeVector gammaPolarization1;
00631 
00632   G4double xParallel = normalisation*cosBeta;
00633   G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
00634   G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
00635   G4double xPerpendicular = 0.;
00636   G4double yPerpendicular = (costheta)*sinBeta/normalisation;
00637   G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
00638 
00639   G4double xTotal = (xParallel + xPerpendicular);
00640   G4double yTotal = (yParallel + yPerpendicular);
00641   G4double zTotal = (zParallel + zPerpendicular);
00642   
00643   gammaPolarization1.setX(xTotal);
00644   gammaPolarization1.setY(yTotal);
00645   gammaPolarization1.setZ(zTotal);
00646   
00647   return gammaPolarization1;
00648 
00649 }
00650 
00651 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00652 
00653 void G4LivermorePolarizedComptonModel::SystemOfRefChange(G4ThreeVector& direction0,
00654                                                     G4ThreeVector& direction1,
00655                                                     G4ThreeVector& polarization0,
00656                                                     G4ThreeVector& polarization1)
00657 {
00658   // direction0 is the original photon direction ---> z
00659   // polarization0 is the original photon polarization ---> x
00660   // need to specify y axis in the real reference frame ---> y 
00661   G4ThreeVector Axis_Z0 = direction0.unit();
00662   G4ThreeVector Axis_X0 = polarization0.unit();
00663   G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
00664 
00665   G4double direction_x = direction1.getX();
00666   G4double direction_y = direction1.getY();
00667   G4double direction_z = direction1.getZ();
00668   
00669   direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
00670   G4double polarization_x = polarization1.getX();
00671   G4double polarization_y = polarization1.getY();
00672   G4double polarization_z = polarization1.getZ();
00673 
00674   polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
00675 
00676 }
00677 
00678 

Generated on Mon May 27 17:48:48 2013 for Geant4 by  doxygen 1.4.7