G4LivermorePolarizedGammaConversionModel.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: G4LivermorePolarizedGammaConversionModel.hh,v 1.2 2010-11-23 16:42:15 flongo Exp $
00027 // GEANT4 tag $Name: not supported by cvs2svn $
00028 //
00029 // Authors: G.Depaola & F.Longo
00030 //
00031 //
00032 
00033 #include "G4LivermorePolarizedGammaConversionModel.hh"
00034 #include "G4PhysicalConstants.hh"
00035 #include "G4SystemOfUnits.hh"
00036 
00037 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00038 
00039 using namespace std;
00040 
00041 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00042 
00043 G4LivermorePolarizedGammaConversionModel::G4LivermorePolarizedGammaConversionModel(
00044    const G4ParticleDefinition*, const G4String& nam)
00045   :G4VEmModel(nam),fParticleChange(0),
00046    isInitialised(false),meanFreePathTable(0),crossSectionHandler(0)
00047 {
00048   lowEnergyLimit = 2*electron_mass_c2;
00049   highEnergyLimit = 100 * GeV;
00050   SetLowEnergyLimit(lowEnergyLimit);
00051   SetHighEnergyLimit(highEnergyLimit);
00052   smallEnergy = 2.*MeV;
00053 
00054   Phi=0.;
00055   Psi=0.;
00056   
00057   verboseLevel= 0;
00058   // Verbosity scale:
00059   // 0 = nothing 
00060   // 1 = warning for energy non-conservation 
00061   // 2 = details of energy budget
00062   // 3 = calculation of cross sections, file openings, samping of atoms
00063   // 4 = entering in methods
00064 
00065   if(verboseLevel > 0) {
00066     G4cout << "Livermore Polarized GammaConversion is constructed " << G4endl
00067            << "Energy range: "
00068            << lowEnergyLimit / keV << " keV - "
00069            << highEnergyLimit / GeV << " GeV"
00070            << G4endl;
00071   }
00072 
00073   crossSectionHandler = new G4CrossSectionHandler();
00074   crossSectionHandler->Initialise(0,lowEnergyLimit,highEnergyLimit,400);
00075 }
00076 
00077 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00078 
00079 G4LivermorePolarizedGammaConversionModel::~G4LivermorePolarizedGammaConversionModel()
00080 {  
00081   delete crossSectionHandler;
00082 }
00083 
00084 
00085 
00086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00087 
00088 void G4LivermorePolarizedGammaConversionModel::Initialise(const G4ParticleDefinition* particle,
00089                                        const G4DataVector& cuts)
00090 {
00091   if (verboseLevel > 3)
00092     G4cout << "Calling G4LivermorePolarizedGammaConversionModel::Initialise()" << G4endl;
00093 
00094   if (crossSectionHandler)
00095   {
00096     crossSectionHandler->Clear();
00097     delete crossSectionHandler;
00098   }
00099 
00100   // Energy limits
00101   /*
00102   // V.Ivanchenko: this was meanless check  
00103   if (LowEnergyLimit() < lowEnergyLimit)
00104     {
00105       G4cout << "G4LivermorePolarizedGammaConversionModel: low energy limit increased from " << 
00106         LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
00107       //      SetLowEnergyLimit(lowEnergyLimit);
00108     }
00109   */
00110   if (HighEnergyLimit() > highEnergyLimit)
00111     {
00112       G4cout << "G4LivermorePolarizedGammaConversionModel: high energy limit decreased from " << 
00113         HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
00114       // V.Ivanchenko: this is forbidden 
00115      // SetHighEnergyLimit(highEnergyLimit);
00116     }
00117 
00118   // Reading of data files - all materials are read
00119   
00120   crossSectionHandler = new G4CrossSectionHandler;
00121   crossSectionHandler->Clear();
00122   G4String crossSectionFile = "pair/pp-cs-";
00123   crossSectionHandler->LoadData(crossSectionFile);
00124 
00125   //
00126   if (verboseLevel > 2) {
00127     G4cout << "Loaded cross section files for Livermore Polarized GammaConversion model" 
00128            << G4endl;
00129   }
00130   InitialiseElementSelectors(particle,cuts);
00131 
00132   if(verboseLevel > 0) {
00133     G4cout << "Livermore Polarized GammaConversion model is initialized " << G4endl
00134            << "Energy range: "
00135            << LowEnergyLimit() / keV << " keV - "
00136            << HighEnergyLimit() / GeV << " GeV"
00137            << G4endl;
00138   }
00139 
00140   //    
00141   if(!isInitialised) { 
00142     isInitialised = true;
00143     fParticleChange = GetParticleChangeForGamma();
00144   }
00145 }
00146 
00147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00148 
00149 G4double G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom(
00150                                        const G4ParticleDefinition*,
00151                                        G4double GammaEnergy,
00152                                        G4double Z, G4double,
00153                                        G4double, G4double)
00154 {
00155   if (verboseLevel > 3) {
00156     G4cout << "G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom()" 
00157            << G4endl;
00158   }
00159   if(Z < 0.9 || GammaEnergy <= lowEnergyLimit) { return 0.0; }
00160   G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
00161   return cs;
00162 }
00163 
00164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00165 
00166 void 
00167 G4LivermorePolarizedGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00168                                                             const G4MaterialCutsCouple* couple,
00169                                                             const G4DynamicParticle* aDynamicGamma,
00170                                                             G4double,
00171                                                             G4double)
00172 {
00173 
00174   // Fluorescence generated according to:
00175   // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
00176   // subshell ionisation by a particle or due to deexcitation or decay of radionuclides",
00177   // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
00178 
00179   if (verboseLevel > 3)
00180     G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedGammaConversionModel" << G4endl;
00181 
00182   G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
00183   // Within energy limit?
00184 
00185   if(photonEnergy <= lowEnergyLimit)
00186     {
00187       fParticleChange->ProposeTrackStatus(fStopAndKill);
00188       fParticleChange->SetProposedKineticEnergy(0.);
00189       return;
00190     }
00191 
00192 
00193   G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
00194   G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
00195 
00196   // Make sure that the polarization vector is perpendicular to the
00197   // gamma direction. If not
00198 
00199   if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
00200     { // only for testing now
00201       gammaPolarization0 = GetRandomPolarization(gammaDirection0);
00202     }
00203   else
00204     {
00205       if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
00206         {
00207           gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
00208         }
00209     }
00210 
00211   // End of Protection
00212 
00213 
00214   G4double epsilon ;
00215   G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
00216 
00217   // Do it fast if photon energy < 2. MeV
00218 
00219   if (photonEnergy < smallEnergy )
00220     {
00221       epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand();
00222     }
00223   else
00224     {
00225 
00226  // Select randomly one element in the current material
00227 
00228       //     G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
00229       //const G4Element* element = crossSectionHandler->SelectRandomElement(couple,photonEnergy);
00230 
00231       const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
00232       const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
00233 
00234       /*
00235       if (element == 0)
00236         {
00237           G4cout << "G4LivermorePolarizedGammaConversionModel::PostStepDoIt - element = 0" << G4endl;
00238         }
00239       */
00240       
00241       G4IonisParamElm* ionisation = element->GetIonisation();
00242       
00243       /*
00244       if (ionisation == 0)
00245         {
00246           G4cout << "G4LivermorePolarizedGammaConversionModel::PostStepDoIt - ionisation = 0" << G4endl;
00247         }
00248       */
00249 
00250 
00251       // Extract Coulomb factor for this Element
00252 
00253       G4double fZ = 8. * (ionisation->GetlogZ3());
00254       if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
00255 
00256       // Limits of the screening variable
00257       G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ;
00258       G4double screenMax = exp ((42.24 - fZ)/8.368) - 0.952 ;
00259       G4double screenMin = std::min(4.*screenFactor,screenMax) ;
00260 
00261       // Limits of the energy sampling
00262       G4double epsilon1 = 0.5 - 0.5 * sqrt(1. - screenMin / screenMax) ;
00263       G4double epsilonMin = std::max(epsilon0Local,epsilon1);
00264       G4double epsilonRange = 0.5 - epsilonMin ;
00265 
00266       // Sample the energy rate of the created electron (or positron)
00267       G4double screen;
00268       G4double gReject ;
00269 
00270       G4double f10 = ScreenFunction1(screenMin) - fZ;
00271       G4double f20 = ScreenFunction2(screenMin) - fZ;
00272       G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
00273       G4double normF2 = std::max(1.5 * f20,0.);
00274 
00275       do {
00276         if (normF1 / (normF1 + normF2) > G4UniformRand() )
00277           {
00278             epsilon = 0.5 - epsilonRange * pow(G4UniformRand(), 0.3333) ;
00279             screen = screenFactor / (epsilon * (1. - epsilon));
00280             gReject = (ScreenFunction1(screen) - fZ) / f10 ;
00281           }
00282         else
00283           {
00284             epsilon = epsilonMin + epsilonRange * G4UniformRand();
00285             screen = screenFactor / (epsilon * (1 - epsilon));
00286             gReject = (ScreenFunction2(screen) - fZ) / f20 ;
00287 
00288 
00289           }
00290       } while ( gReject < G4UniformRand() );
00291 
00292     }   //  End of epsilon sampling
00293   
00294   // Fix charges randomly
00295   
00296   G4double electronTotEnergy;
00297   G4double positronTotEnergy;
00298 
00299 
00300   if (G4int(2*G4UniformRand()))  
00301     {
00302       electronTotEnergy = (1. - epsilon) * photonEnergy;
00303       positronTotEnergy = epsilon * photonEnergy;
00304     }
00305   else
00306     {
00307       positronTotEnergy = (1. - epsilon) * photonEnergy;
00308       electronTotEnergy = epsilon * photonEnergy;
00309     }
00310 
00311   // Scattered electron (positron) angles. ( Z - axis along the parent photon)
00312   // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
00313   // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
00314 
00315 /*
00316   G4double u;
00317   const G4double a1 = 0.625;
00318   G4double a2 = 3. * a1;
00319 
00320   if (0.25 > G4UniformRand())
00321     {
00322       u = - log(G4UniformRand() * G4UniformRand()) / a1 ;
00323     }
00324   else
00325     {
00326       u = - log(G4UniformRand() * G4UniformRand()) / a2 ;
00327     }
00328 */
00329 
00330   G4double Ene = electronTotEnergy/electron_mass_c2; // Normalized energy
00331 
00332   G4double cosTheta = 0.;
00333   G4double sinTheta = 0.;
00334 
00335   SetTheta(&cosTheta,&sinTheta,Ene);
00336 
00337   //  G4double theta = u * electron_mass_c2 / photonEnergy ;
00338   //  G4double phi  = twopi * G4UniformRand() ;
00339 
00340   G4double phi,psi=0.;
00341 
00342   //corrected e+ e- angular angular distribution //preliminary!
00343 
00344   //  if(photonEnergy>50*MeV)
00345   // {
00346   phi = SetPhi(photonEnergy);
00347   psi = SetPsi(photonEnergy,phi);
00348   //  }
00349   //else
00350   // {
00351   //psi = G4UniformRand()*2.*pi;
00352   //phi = pi; // coplanar
00353   // }
00354 
00355   Psi = psi;
00356   Phi = phi;
00357   //G4cout << "PHI " << phi << G4endl;
00358   //G4cout << "PSI " << psi << G4endl;
00359 
00360   G4double phie = psi; //azimuthal angle for the electron
00361 
00362   G4double dirX = sinTheta*cos(phie);
00363   G4double dirY = sinTheta*sin(phie);
00364   G4double dirZ = cosTheta;
00365   G4ThreeVector electronDirection(dirX,dirY,dirZ);
00366   // Kinematics of the created pair:
00367   // the electron and positron are assumed to have a symetric angular
00368   // distribution with respect to the Z axis along the parent photon
00369 
00370   //G4double localEnergyDeposit = 0. ;
00371 
00372   G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
00373 
00374   SystemOfRefChange(gammaDirection0,electronDirection,
00375                     gammaPolarization0);
00376 
00377   G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
00378                                                         electronDirection,
00379                                                         electronKineEnergy);
00380 
00381   // The e+ is always created (even with kinetic energy = 0) for further annihilation
00382 
00383   Ene = positronTotEnergy/electron_mass_c2; // Normalized energy
00384 
00385   cosTheta = 0.;
00386   sinTheta = 0.;
00387 
00388   SetTheta(&cosTheta,&sinTheta,Ene);
00389   G4double phip = phie+phi; //azimuthal angle for the positron
00390 
00391   dirX = sinTheta*cos(phip);
00392   dirY = sinTheta*sin(phip);
00393   dirZ = cosTheta;
00394   G4ThreeVector positronDirection(dirX,dirY,dirZ);
00395 
00396   G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
00397   SystemOfRefChange(gammaDirection0,positronDirection,
00398                     gammaPolarization0);
00399 
00400   // Create G4DynamicParticle object for the particle2
00401   G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
00402                                                        positronDirection, positronKineEnergy);
00403 
00404 
00405   fvect->push_back(particle1);
00406   fvect->push_back(particle2);
00407 
00408 
00409 
00410   // Kill the incident photon
00411 
00412 
00413 
00414   // Create lists of pointers to DynamicParticles (photons and electrons)
00415   // (Is the electron vector necessary? To be checked)
00416   //  std::vector<G4DynamicParticle*>* photonVector = 0;
00417   //std::vector<G4DynamicParticle*> electronVector;
00418 
00419   fParticleChange->ProposeMomentumDirection( 0., 0., 0. );
00420   fParticleChange->SetProposedKineticEnergy(0.);
00421   fParticleChange->ProposeTrackStatus(fStopAndKill);
00422 
00423 }
00424 
00425 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00426 
00427 G4double G4LivermorePolarizedGammaConversionModel::ScreenFunction1(G4double screenVariable)
00428 {
00429   // Compute the value of the screening function 3*phi1 - phi2
00430 
00431   G4double value;
00432 
00433   if (screenVariable > 1.)
00434     value = 42.24 - 8.368 * log(screenVariable + 0.952);
00435   else
00436     value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
00437 
00438   return value;
00439 }
00440 
00441 
00442 
00443 G4double G4LivermorePolarizedGammaConversionModel::ScreenFunction2(G4double screenVariable)
00444 {
00445   // Compute the value of the screening function 1.5*phi1 - 0.5*phi2
00446 
00447   G4double value;
00448 
00449   if (screenVariable > 1.)
00450     value = 42.24 - 8.368 * log(screenVariable + 0.952);
00451   else
00452     value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
00453 
00454   return value;
00455 }
00456 
00457 
00458 void G4LivermorePolarizedGammaConversionModel::SetTheta(G4double* p_cosTheta, G4double* p_sinTheta, G4double Energy)
00459 {
00460 
00461   // to avoid computational errors since Theta could be very small
00462   // Energy in Normalized Units (!)
00463 
00464   G4double Momentum = sqrt(Energy*Energy -1);
00465   G4double Rand = G4UniformRand();
00466 
00467   *p_cosTheta = (Energy*((2*Rand)- 1) + Momentum)/((Momentum*(2*Rand-1))+Energy);
00468   *p_sinTheta = (2*sqrt(Rand*(1-Rand)))/(Momentum*(2*Rand-1)+Energy);
00469 }
00470 
00471 
00472 
00473 G4double G4LivermorePolarizedGammaConversionModel::SetPhi(G4double Energy)
00474 {
00475 
00476 
00477   G4double value = 0.;
00478   G4double Ene = Energy/MeV;
00479 
00480   G4double pl[4];
00481 
00482 
00483   G4double pt[2];
00484   G4double xi = 0;
00485   G4double xe = 0.;
00486   G4double n1=0.;
00487   G4double n2=0.;
00488 
00489 
00490   if (Ene>=50.)
00491     {
00492       const G4double ay0=5.6, by0=18.6, aa0=2.9, ba0 = 8.16E-3;
00493       const G4double aw = 0.0151, bw = 10.7, cw = -410.;
00494 
00495       const G4double axc = 3.1455, bxc = -1.11, cxc = 310.;
00496 
00497       pl[0] = Fln(ay0,by0,Ene);
00498       pl[1] = aa0 + ba0*(Ene);
00499       pl[2] = Poli(aw,bw,cw,Ene);
00500       pl[3] = Poli(axc,bxc,cxc,Ene);
00501 
00502       const G4double abf = 3.1216, bbf = 2.68;
00503       pt[0] = -1.4;
00504       pt[1] = abf + bbf/Ene;
00505 
00506 
00507 
00508       //G4cout << "PL > 50. "<< pl[0] << " " << pl[1] << " " << pl[2] << " " <<pl[3] << " " << G4endl;
00509 
00510       xi = 3.0;
00511       xe = Encu(pl,pt,xi);
00512       //G4cout << "ENCU "<< xe << G4endl;
00513       n1 = Fintlor(pl,pi) - Fintlor(pl,xe);
00514       n2 = Finttan(pt,xe) - Finttan(pt,0.);
00515     }
00516   else
00517     {
00518       const G4double ay0=0.144, by0=0.11;
00519       const G4double aa0=2.7, ba0 = 2.74;
00520       const G4double aw = 0.21, bw = 10.8, cw = -58.;
00521       const G4double axc = 3.17, bxc = -0.87, cxc = -6.;
00522 
00523       pl[0] = Fln(ay0, by0, Ene);
00524       pl[1] = Fln(aa0, ba0, Ene);
00525       pl[2] = Poli(aw,bw,cw,Ene);
00526       pl[3] = Poli(axc,bxc,cxc,Ene);
00527 
00528       //G4cout << "PL < 50."<< pl[0] << " " << pl[1] << " " << pl[2] << " " <<pl[3] << " " << G4endl;
00529       //G4cout << "ENCU "<< xe << G4endl;
00530       n1 = Fintlor(pl,pi) - Fintlor(pl,xe);
00531 
00532     }
00533 
00534 
00535   G4double n=0.;
00536   n = n1+n2;
00537 
00538   G4double c1 = 0.;
00539   c1 = Glor(pl, xe);
00540 
00541 /*
00542   G4double xm = 0.;
00543   xm = Flor(pl,pl[3])*Glor(pl,pl[3]);
00544 */
00545 
00546   G4double r1,r2,r3;
00547   G4double xco=0.;
00548 
00549   if (Ene>=50.)
00550     {
00551       r1= G4UniformRand();
00552       if( r1>=n2/n)
00553         {
00554           do
00555             {
00556               r2 = G4UniformRand();
00557               value = Finvlor(pl,xe,r2);
00558               xco = Glor(pl,value)/c1;
00559               r3 = G4UniformRand();
00560             } while(r3>=xco);
00561         }
00562       else
00563         {
00564           value = Finvtan(pt,n,r1);
00565         }
00566     }
00567   else
00568     {
00569       do
00570         {
00571           r2 = G4UniformRand();
00572           value = Finvlor(pl,xe,r2);
00573           xco = Glor(pl,value)/c1;
00574           r3 = G4UniformRand();
00575         } while(r3>=xco);
00576     }
00577 
00578   //  G4cout << "PHI = " <<value <<  G4endl;
00579   return value;
00580 }
00581 G4double G4LivermorePolarizedGammaConversionModel::SetPsi(G4double Energy, G4double PhiLocal)
00582 {
00583 
00584   G4double value = 0.;
00585   G4double Ene = Energy/MeV;
00586 
00587   G4double p0l[4];
00588   G4double ppml[4];
00589   G4double p0t[2];
00590   G4double ppmt[2];
00591 
00592   G4double xi = 0.;
00593   G4double xe0 = 0.;
00594   G4double xepm = 0.;
00595 
00596   if (Ene>=50.)
00597     {
00598       const G4double ay00 = 3.4, by00 = 9.8, aa00 = 1.34, ba00 = 5.3;
00599       const G4double aw0 = 0.014, bw0 = 9.7, cw0 = -2.E4;
00600       const G4double axc0 = 3.1423, bxc0 = -2.35, cxc0 = 0.;
00601       const G4double ay0p = 1.53, by0p = 3.2, aap = 0.67, bap = 8.5E-3;
00602       const G4double awp = 6.9E-3, bwp = 12.6, cwp = -3.8E4;
00603       const G4double axcp = 2.8E-3,bxcp = -3.133;
00604       const G4double abf0 = 3.1213, bbf0 = 2.61;
00605       const G4double abfpm = 3.1231, bbfpm = 2.84;
00606 
00607       p0l[0] = Fln(ay00, by00, Ene);
00608       p0l[1] = Fln(aa00, ba00, Ene);
00609       p0l[2] = Poli(aw0, bw0, cw0, Ene);
00610       p0l[3] = Poli(axc0, bxc0, cxc0, Ene);
00611 
00612       ppml[0] = Fln(ay0p, by0p, Ene);
00613       ppml[1] = aap + bap*(Ene);
00614       ppml[2] = Poli(awp, bwp, cwp, Ene);
00615       ppml[3] = Fln(axcp,bxcp,Ene);
00616 
00617       p0t[0] = -0.81;
00618       p0t[1] = abf0 + bbf0/Ene;
00619       ppmt[0] = -0.6;
00620       ppmt[1] = abfpm + bbfpm/Ene;
00621 
00622       //G4cout << "P0L > 50"<< p0l[0] << " " << p0l[1] << " " << p0l[2] << " " <<p0l[3] << " " << G4endl;
00623       //G4cout << "PPML > 50"<< ppml[0] << " " << ppml[1] << " " << ppml[2] << " " <<ppml[3] << " " << G4endl;
00624 
00625       xi = 3.0;
00626       xe0 = Encu(p0l, p0t, xi);
00627       //G4cout << "ENCU1 "<< xe0 << G4endl;
00628       xepm = Encu(ppml, ppmt, xi);
00629       //G4cout << "ENCU2 "<< xepm << G4endl;
00630     }
00631   else
00632     {
00633       const G4double ay00 = 2.82, by00 = 6.35;
00634       const G4double aa00 = -1.75, ba00 = 0.25;
00635 
00636       const G4double aw0 = 0.028, bw0 = 5., cw0 = -50.;
00637       const G4double axc0 = 3.14213, bxc0 = -2.3, cxc0 = 5.7;
00638       const G4double ay0p = 1.56, by0p = 3.6;
00639       const G4double aap = 0.86, bap = 8.3E-3;
00640       const G4double awp = 0.022, bwp = 7.4, cwp = -51.;
00641       const G4double xcp = 3.1486;
00642 
00643       p0l[0] = Fln(ay00, by00, Ene);
00644       p0l[1] = aa00+pow(Ene, ba00);
00645       p0l[2] = Poli(aw0, bw0, cw0, Ene);
00646       p0l[3] = Poli(axc0, bxc0, cxc0, Ene);
00647       ppml[0] = Fln(ay0p, by0p, Ene);
00648       ppml[1] = aap + bap*(Ene);
00649       ppml[2] = Poli(awp, bwp, cwp, Ene);
00650       ppml[3] = xcp;
00651 
00652     }
00653 
00654   G4double a,b=0.;
00655 
00656   if (Ene>=50.)
00657     {
00658       if (PhiLocal>xepm)
00659         {
00660           b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml,PhiLocal));
00661         }
00662       else
00663         {
00664           b = Ftan(ppmt,PhiLocal);
00665         }
00666       if (PhiLocal>xe0)
00667         {
00668           a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,PhiLocal));
00669         }
00670       else
00671         {
00672           a = Ftan(p0t,PhiLocal);
00673         }
00674     }
00675   else
00676     {
00677       b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml,PhiLocal));
00678       a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,PhiLocal));
00679     }
00680   G4double nr =0.;
00681 
00682   if (b>a)
00683     {
00684       nr = 1./b;
00685     }
00686   else
00687     {
00688       nr = 1./a;
00689     }
00690 
00691   G4double r1,r2=0.;
00692   G4double r3 =-1.;
00693   do
00694     {
00695       r1 = G4UniformRand();
00696       r2 = G4UniformRand();
00697       value = r2*pi;
00698       r3 = nr*(a*cos(value)*cos(value) + b*sin(value)*sin(value));
00699     }while(r1>r3);
00700 
00701   return value;
00702 }
00703 
00704 
00705 G4double G4LivermorePolarizedGammaConversionModel::Poli
00706 (G4double a, G4double b, G4double c, G4double x)
00707 {
00708   G4double value=0.;
00709   if(x>0.)
00710     {
00711       value =(a + b/x + c/(x*x*x));
00712     }
00713   else
00714     {
00715       //G4cout << "ERROR in Poli! " << G4endl;
00716     }
00717   return value;
00718 }
00719 G4double G4LivermorePolarizedGammaConversionModel::Fln
00720 (G4double a, G4double b, G4double x)
00721 {
00722   G4double value=0.;
00723   if(x>0.)
00724     {
00725       value =(a*log(x)-b);
00726     }
00727   else
00728     {
00729       //G4cout << "ERROR in Fln! " << G4endl;
00730     }
00731   return value;
00732 }
00733 
00734 
00735 G4double G4LivermorePolarizedGammaConversionModel::Encu
00736 (G4double* p_p1, G4double* p_p2, G4double x0)
00737 {
00738   G4int i=0;
00739   G4double fx = 1.;
00740   G4double x = x0;
00741   const G4double xmax = 3.0;
00742 
00743   for(i=0; i<100; ++i)
00744     {
00745       fx = (Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p_p2, x))/
00746         (Fdlor(p_p1,x) - Fdtan(p_p2,x));
00747       x -= fx;
00748       if(x > xmax) { return xmax; }
00749       //      x -= (Flor(p_p1, x)*Glor(p_p1,x) - Ftan(p_p2, x))/
00750       //  (Fdlor(p_p1,x) - Fdtan(p_p2,x));
00751       // fx = Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p_p2, x);
00752       // G4cout << std::fabs(fx) << " " << i << " " << x << "dentro ENCU " << G4endl;
00753       if(std::fabs(fx) <= x*1.0e-6) { break; }
00754     } 
00755 
00756   if(x < 0.0) { x = 0.0; }
00757   return x;
00758 }
00759 
00760 
00761 G4double G4LivermorePolarizedGammaConversionModel::Flor(G4double* p_p1, G4double x)
00762 {
00763   G4double value =0.;
00764   // G4double y0 = p_p1[0];
00765   // G4double A = p_p1[1];
00766   G4double w = p_p1[2];
00767   G4double xc = p_p1[3];
00768 
00769   value = 1./(pi*(w*w + 4.*(x-xc)*(x-xc)));
00770   return value;
00771 }
00772 
00773 
00774 G4double G4LivermorePolarizedGammaConversionModel::Glor(G4double* p_p1, G4double x)
00775 {
00776   G4double value =0.;
00777   G4double y0 = p_p1[0];
00778   G4double A = p_p1[1];
00779   G4double w = p_p1[2];
00780   G4double xc = p_p1[3];
00781 
00782   value = (y0 *pi*(w*w +  4.*(x-xc)*(x-xc)) + 2.*A*w);
00783   return value;
00784 }
00785 
00786 
00787 G4double G4LivermorePolarizedGammaConversionModel::Fdlor(G4double* p_p1, G4double x)
00788 {
00789   G4double value =0.;
00790   //G4double y0 = p_p1[0];
00791   G4double A = p_p1[1];
00792   G4double w = p_p1[2];
00793   G4double xc = p_p1[3];
00794 
00795   value = (-16.*A*w*(x-xc))/
00796     (pi*(w*w+4.*(x-xc)*(x-xc))*(w*w+4.*(x-xc)*(x-xc)));
00797   return value;
00798 }
00799 
00800 
00801 G4double G4LivermorePolarizedGammaConversionModel::Fintlor(G4double* p_p1, G4double x)
00802 {
00803   G4double value =0.;
00804   G4double y0 = p_p1[0];
00805   G4double A = p_p1[1];
00806   G4double w = p_p1[2];
00807   G4double xc = p_p1[3];
00808 
00809   value = y0*x + A*atan( 2*(x-xc)/w) / pi;
00810   return value;
00811 }
00812 
00813 
00814 G4double G4LivermorePolarizedGammaConversionModel::Finvlor(G4double* p_p1, G4double x, G4double r)
00815 {
00816   G4double value = 0.;
00817   G4double nor = 0.;
00818   //G4double y0 = p_p1[0];
00819   //  G4double A = p_p1[1];
00820   G4double w = p_p1[2];
00821   G4double xc = p_p1[3];
00822 
00823   nor = atan(2.*(pi-xc)/w)/(2.*pi*w) - atan(2.*(x-xc)/w)/(2.*pi*w);
00824   value = xc - (w/2.)*tan(-2.*r*nor*pi*w+atan(2*(xc-x)/w));
00825 
00826   return value;
00827 }
00828 
00829 
00830 G4double G4LivermorePolarizedGammaConversionModel::Ftan(G4double* p_p1, G4double x)
00831 {
00832   G4double value =0.;
00833   G4double a = p_p1[0];
00834   G4double b = p_p1[1];
00835 
00836   value = a /(x-b);
00837   return value;
00838 }
00839 
00840 
00841 G4double G4LivermorePolarizedGammaConversionModel::Fdtan(G4double* p_p1, G4double x)
00842 {
00843   G4double value =0.;
00844   G4double a = p_p1[0];
00845   G4double b = p_p1[1];
00846 
00847   value = -1.*a / ((x-b)*(x-b));
00848   return value;
00849 }
00850 
00851 
00852 G4double G4LivermorePolarizedGammaConversionModel::Finttan(G4double* p_p1, G4double x)
00853 {
00854   G4double value =0.;
00855   G4double a = p_p1[0];
00856   G4double b = p_p1[1];
00857 
00858 
00859   value = a*log(b-x);
00860   return value;
00861 }
00862 
00863 
00864 G4double G4LivermorePolarizedGammaConversionModel::Finvtan(G4double* p_p1, G4double cnor, G4double r)
00865 {
00866   G4double value =0.;
00867   G4double a = p_p1[0];
00868   G4double b = p_p1[1];
00869 
00870   value = b*(1-exp(r*cnor/a));
00871 
00872   return value;
00873 }
00874 
00875 
00876 
00877 
00878 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00879 
00880 G4ThreeVector G4LivermorePolarizedGammaConversionModel::SetPerpendicularVector(G4ThreeVector& a)
00881 {
00882   G4double dx = a.x();
00883   G4double dy = a.y();
00884   G4double dz = a.z();
00885   G4double x = dx < 0.0 ? -dx : dx;
00886   G4double y = dy < 0.0 ? -dy : dy;
00887   G4double z = dz < 0.0 ? -dz : dz;
00888   if (x < y) {
00889     return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
00890   }else{
00891     return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
00892   }
00893 }
00894 
00895 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00896 
00897 G4ThreeVector G4LivermorePolarizedGammaConversionModel::GetRandomPolarization(G4ThreeVector& direction0)
00898 {
00899   G4ThreeVector d0 = direction0.unit();
00900   G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
00901   G4ThreeVector a0 = a1.unit(); // unit vector
00902 
00903   G4double rand1 = G4UniformRand();
00904   
00905   G4double angle = twopi*rand1; // random polar angle
00906   G4ThreeVector b0 = d0.cross(a0); // cross product
00907   
00908   G4ThreeVector c;
00909   
00910   c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
00911   c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
00912   c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
00913   
00914   G4ThreeVector c0 = c.unit();
00915 
00916   return c0;
00917   
00918 }
00919 
00920 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00921 
00922 G4ThreeVector G4LivermorePolarizedGammaConversionModel::GetPerpendicularPolarization
00923 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
00924 {
00925 
00926   // 
00927   // The polarization of a photon is always perpendicular to its momentum direction.
00928   // Therefore this function removes those vector component of gammaPolarization, which
00929   // points in direction of gammaDirection
00930   //
00931   // Mathematically we search the projection of the vector a on the plane E, where n is the
00932   // plains normal vector.
00933   // The basic equation can be found in each geometry book (e.g. Bronstein):
00934   // p = a - (a o n)/(n o n)*n
00935   
00936   return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
00937 }
00938 
00939 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00940 
00941 
00942 void G4LivermorePolarizedGammaConversionModel::SystemOfRefChange
00943     (G4ThreeVector& direction0,G4ThreeVector& direction1,
00944      G4ThreeVector& polarization0)
00945 {
00946   // direction0 is the original photon direction ---> z
00947   // polarization0 is the original photon polarization ---> x
00948   // need to specify y axis in the real reference frame ---> y 
00949   G4ThreeVector Axis_Z0 = direction0.unit();
00950   G4ThreeVector Axis_X0 = polarization0.unit();
00951   G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
00952   
00953   G4double direction_x = direction1.getX();
00954   G4double direction_y = direction1.getY();
00955   G4double direction_z = direction1.getZ();
00956   
00957   direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 +  direction_z*Axis_Z0).unit();
00958   
00959 }
00960 
00961 
00962 
00963 

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