Geant4: G4AdjointhIonisationModel.cc Source 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: G4AdjointhIonisationModel.cc 69844 2013-05-16 09:19:33Z gcosmo $
00027 //
00028 #include "G4AdjointhIonisationModel.hh"
00029 
00030 #include "G4PhysicalConstants.hh"
00031 #include "G4SystemOfUnits.hh"
00032 #include "G4AdjointCSManager.hh"
00033 #include "G4Integrator.hh"
00034 #include "G4TrackStatus.hh"
00035 #include "G4ParticleChange.hh"
00036 #include "G4AdjointElectron.hh"
00037 #include "G4AdjointProton.hh"
00038 #include "G4AdjointInterpolator.hh"
00039 #include "G4BetheBlochModel.hh"
00040 #include "G4BraggModel.hh"
00041 #include "G4Proton.hh"
00042 #include "G4NistManager.hh"
00043 
00045 //
00046 G4AdjointhIonisationModel::G4AdjointhIonisationModel(G4ParticleDefinition* projectileDefinition):
00047   G4VEmAdjointModel("Adjoint_hIonisation")
00048 { 
00049 
00050 
00051 
00052   UseMatrix =true;
00053   UseMatrixPerElement = true;
00054   ApplyCutInRange = true;
00055   UseOnlyOneMatrixForAllElements = true; 
00056   CS_biasing_factor =1.;
00057   second_part_of_same_type =false;
00058   
00059   //The direct EM Modfel is taken has BetheBloch it is only used for the computation 
00060   // of the differential cross section.
00061   //The Bragg model could be used  as an alternative as  it offers the same differential cross section
00062   
00063   theDirectEMModel = new G4BetheBlochModel(projectileDefinition);
00064   theBraggDirectEMModel = new G4BraggModel(projectileDefinition); 
00065   theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
00066   
00067   theDirectPrimaryPartDef = projectileDefinition;
00068   theAdjEquivOfDirectPrimPartDef = 0;
00069   if (projectileDefinition == G4Proton::Proton()) {
00070         theAdjEquivOfDirectPrimPartDef = G4AdjointProton::AdjointProton();
00071         
00072   }
00073   
00074   DefineProjectileProperty();
00075 }
00077 //
00078 G4AdjointhIonisationModel::~G4AdjointhIonisationModel()
00079 {;}
00080 
00081 
00083 //
00084 void G4AdjointhIonisationModel::SampleSecondaries(const G4Track& aTrack,
00085                                 G4bool IsScatProjToProjCase,
00086                                 G4ParticleChange* fParticleChange)
00087 { 
00088  if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
00089  
00090  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
00091   
00092  //Elastic inverse scattering 
00093  //---------------------------------------------------------
00094  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
00095  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
00096 
00097  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
00098         return;
00099  }
00100  
00101  //Sample secondary energy
00102  //-----------------------
00103  G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
00104  CorrectPostStepWeight(fParticleChange, 
00105                               aTrack.GetWeight(), 
00106                               adjointPrimKinEnergy,
00107                               projectileKinEnergy,
00108                               IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
00109 
00110                  
00111  //Kinematic: 
00112  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest  and gives
00113  // him part of its  energy
00114  //----------------------------------------------------------------------------------------
00115  
00116  G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
00117  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
00118  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
00119  
00120  
00121  
00122  //Companion
00123  //-----------
00124  G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
00125  if (IsScatProjToProjCase) {
00126         companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
00127  } 
00128  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
00129  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;    
00130  
00131  
00132  //Projectile momentum
00133  //--------------------
00134  G4double  P_parallel = (adjointPrimP*adjointPrimP +  projectileP2 - companionP2)/(2.*adjointPrimP);
00135  G4double  P_perp = std::sqrt( projectileP2 -  P_parallel*P_parallel);
00136  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
00137  G4double phi =G4UniformRand()*2.*3.1415926;
00138  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
00139  projectileMomentum.rotateUz(dir_parallel);
00140  
00141  
00142  
00143  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
00144         fParticleChange->ProposeTrackStatus(fStopAndKill);
00145         fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
00146         //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
00147  }
00148  else {
00149         fParticleChange->ProposeEnergy(projectileKinEnergy);
00150         fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
00151  }      
00152         
00153 
00154   
00155 
00156 }
00157 
00159 //
00160 void G4AdjointhIonisationModel::RapidSampleSecondaries(const G4Track& aTrack,
00161                        G4bool IsScatProjToProjCase,
00162                        G4ParticleChange* fParticleChange)
00163 { 
00164 
00165  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
00166  DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
00167  
00168  
00169  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
00170  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
00171  
00172  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
00173         return;
00174  }
00175   
00176  G4double projectileKinEnergy =0.;
00177  G4double eEnergy=0.;
00178  G4double newCS=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
00179  if (!IsScatProjToProjCase){//1/E^2 distribution
00180         
00181         eEnergy=adjointPrimKinEnergy;
00182         G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
00183         G4double Emin=  GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);
00184         if (Emin>=Emax) return;
00185         G4double a=1./Emax;
00186         G4double b=1./Emin;
00187         newCS=newCS*(b-a)/eEnergy;
00188         
00189         projectileKinEnergy =1./(b- (b-a)*G4UniformRand()); 
00190         
00191         
00192  }
00193  else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
00194         G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
00195         if (Emin>=Emax) return;
00196         G4double diff1=Emin-adjointPrimKinEnergy;
00197         G4double diff2=Emax-adjointPrimKinEnergy;
00198         
00199         G4double t1=adjointPrimKinEnergy*(1./diff1-1./diff2);
00200         G4double t2=adjointPrimKinEnergy*(1./Emin-1./Emax);
00201         /*G4double f31=diff1/Emin;
00202         G4double f32=diff2/Emax/f31;*/
00203         G4double t3=2.*std::log(Emax/Emin);
00204         G4double sum_t=t1+t2+t3;
00205         newCS=newCS*sum_t/adjointPrimKinEnergy/adjointPrimKinEnergy;
00206         G4double t=G4UniformRand()*sum_t;
00207         if (t <=t1 ){
00208                 G4double q= G4UniformRand()*t1/adjointPrimKinEnergy ;
00209                 projectileKinEnergy =adjointPrimKinEnergy +1./(1./diff1-q);
00210                 
00211         }
00212         else if (t <=t2 )  {
00213                 G4double q= G4UniformRand()*t2/adjointPrimKinEnergy;
00214                 projectileKinEnergy =1./(1./Emin-q);
00215         }
00216         else {  
00217                 projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand());
00218                 
00219         }
00220         eEnergy=projectileKinEnergy-adjointPrimKinEnergy;
00221         
00222         
00223  }
00224 
00225  
00226 
00227  G4double diffCS_perAtom_Used=twopi_mc2_rcl2*mass*adjointPrimKinEnergy/projectileKinEnergy/projectileKinEnergy/eEnergy/eEnergy; 
00228   
00229   
00230                                                         
00231  //Weight correction
00232  //-----------------------
00233  //First w_corr is set to the ratio between adjoint total CS and fwd total CS
00234  G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
00235  
00236  //G4cout<<w_corr<<G4endl;
00237  w_corr*=newCS/lastCS;
00238  //G4cout<<w_corr<<G4endl;
00239  //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the one consistent with the direct model
00240  //Here we consider the true  diffCS as the one obtained by the numerical differentiation over Tcut of the direct CS
00241  
00242  G4double diffCS = DiffCrossSectionPerAtomPrimToSecond(projectileKinEnergy, eEnergy,1,1);
00243  w_corr*=diffCS/diffCS_perAtom_Used;
00244  //G4cout<<w_corr<<G4endl;
00245            
00246  G4double new_weight = aTrack.GetWeight()*w_corr;
00247  fParticleChange->SetParentWeightByProcess(false);
00248  fParticleChange->SetSecondaryWeightByProcess(false);
00249  fParticleChange->ProposeParentWeight(new_weight);
00250  
00251  
00252  
00253  
00254  //Kinematic: 
00255  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest  and gives
00256  // him part of its  energy
00257  //----------------------------------------------------------------------------------------
00258  
00259  G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
00260  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
00261  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
00262  
00263  
00264  
00265  //Companion
00266  //-----------
00267  G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
00268  if (IsScatProjToProjCase) {
00269         companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
00270  } 
00271  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
00272  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;    
00273  
00274  
00275  //Projectile momentum
00276  //--------------------
00277  G4double  P_parallel = (adjointPrimP*adjointPrimP +  projectileP2 - companionP2)/(2.*adjointPrimP);
00278  G4double  P_perp = std::sqrt( projectileP2 -  P_parallel*P_parallel);
00279  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
00280  G4double phi =G4UniformRand()*2.*3.1415926;
00281  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
00282  projectileMomentum.rotateUz(dir_parallel);
00283  
00284  
00285  
00286  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
00287         fParticleChange->ProposeTrackStatus(fStopAndKill);
00288         fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
00289         //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
00290  }
00291  else {
00292         fParticleChange->ProposeEnergy(projectileKinEnergy);
00293         fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
00294  }      
00295         
00296 
00297  
00298  
00299 
00300 
00301 
00302 } 
00303                 
00305 //
00306 G4double G4AdjointhIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
00307                                       G4double kinEnergyProj, 
00308                                       G4double kinEnergyProd,
00309                                       G4double Z, 
00310                                       G4double A)
00311 {//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
00312   
00313 
00314  
00315  G4double dSigmadEprod=0;
00316  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
00317  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
00318  
00319  
00320  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile 
00321         G4double Tmax=kinEnergyProj;
00322         
00323         G4double E1=kinEnergyProd;
00324         G4double E2=kinEnergyProd*1.000001;
00325         G4double dE=(E2-E1);
00326         G4double sigma1,sigma2;
00327         if (kinEnergyProj >2.*MeV){
00328             sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
00329             sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
00330         }
00331         else {
00332             sigma1=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
00333             sigma2=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
00334         }
00335         
00336         
00337         dSigmadEprod=(sigma1-sigma2)/dE;
00338         if (dSigmadEprod>1.) {
00339                 G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
00340                 G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
00341                 G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
00342                 
00343         }
00344         
00345          
00346         
00347          //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
00348          //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
00349          //to test the rejection of a secondary
00350          //-------------------------
00351          
00352          //Source code taken from    G4BetheBlochModel::SampleSecondaries
00353          
00354          G4double deltaKinEnergy = kinEnergyProd;
00355          
00356          //Part of the taken code
00357          //----------------------
00358          
00359          
00360          
00361          // projectile formfactor - suppresion of high energy
00362          // delta-electron production at high energy
00363          G4double x = formfact*deltaKinEnergy;
00364           if(x > 1.e-6) {
00365                 
00366                 
00367                 G4double totEnergy     = kinEnergyProj + mass;
00368                 G4double etot2         = totEnergy*totEnergy;
00369                 G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
00370                 G4double f;
00371                 G4double f1 = 0.0;
00372                 f = 1.0 - beta2*deltaKinEnergy/Tmax;
00373                 if( 0.5 == spin ) {
00374                         f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
00375                         f += f1;
00376                 }
00377                 G4double x1 = 1.0 + x;
00378                 G4double gg  = 1.0/(x1*x1);
00379                 if( 0.5 == spin ) {
00380                         G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
00381                         gg *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
00382                 }
00383                 if(gg > 1.0) {
00384                         G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g
00385                                 << G4endl;
00386                         gg=1.;
00387                 }
00388                 //G4cout<<"gg"<<gg<<G4endl;
00389                 dSigmadEprod*=gg;               
00390          }
00391          
00392   }
00393 
00394  return dSigmadEprod;   
00395 }
00396 
00397 
00398 
00400 //
00401 void G4AdjointhIonisationModel::DefineProjectileProperty()
00402 {   
00403     //Slightly modified code taken from G4BetheBlochModel::SetParticle
00404     //------------------------------------------------
00405     G4String pname = theDirectPrimaryPartDef->GetParticleName();
00406     if (theDirectPrimaryPartDef->GetParticleType() == "nucleus" &&
00407         pname != "deuteron" && pname != "triton") {
00408       isIon = true;
00409     }
00410     
00411     mass = theDirectPrimaryPartDef->GetPDGMass();
00412     mass_ratio_projectile = proton_mass_c2/theDirectPrimaryPartDef->GetPDGMass();;
00413     spin = theDirectPrimaryPartDef->GetPDGSpin();
00414     G4double q = theDirectPrimaryPartDef->GetPDGCharge()/eplus;
00415     chargeSquare = q*q;
00416     ratio = electron_mass_c2/mass;
00417     ratio2 = ratio*ratio;
00418     one_plus_ratio_2=(1+ratio)*(1+ratio);
00419     one_minus_ratio_2=(1-ratio)*(1-ratio);
00420     G4double magmom = theDirectPrimaryPartDef->GetPDGMagneticMoment()
00421       *mass/(0.5*eplus*hbar_Planck*c_squared);
00422     magMoment2 = magmom*magmom - 1.0;
00423     formfact = 0.0;
00424     if(theDirectPrimaryPartDef->GetLeptonNumber() == 0) {
00425       G4double x = 0.8426*GeV;
00426       if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
00427       else if(mass > GeV) {
00428         x /= G4NistManager::Instance()->GetZ13(mass/proton_mass_c2);
00429         //      tlimit = 51.2*GeV*A13[iz]*A13[iz];
00430       }
00431       formfact = 2.0*electron_mass_c2/(x*x);
00432       tlimit   = 2.0/formfact;
00433    }
00434 }
00435 
00437 //
00438 G4double G4AdjointhIonisationModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
00439                                              G4double primEnergy,
00440                                              G4bool IsScatProjToProjCase)
00441 { 
00442   if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
00443   DefineCurrentMaterial(aCouple);
00444         
00445   
00446   G4double Cross=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
00447   
00448   if (!IsScatProjToProjCase ){
00449         G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
00450         G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
00451         if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) {
00452                 Cross*=(1./Emin_proj -1./Emax_proj)/primEnergy;
00453         } 
00454         else Cross=0.;
00455                 
00456         
00457 
00458         
00459         
00460         
00461   }
00462   else {
00463         G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
00464         G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond);
00465         G4double diff1=Emin_proj-primEnergy;
00466         G4double diff2=Emax_proj-primEnergy;
00467         G4double t1=(1./diff1+1./Emin_proj-1./diff2-1./Emax_proj)/primEnergy;
00468         //G4double t2=2.*std::log(diff2*Emin_proj/Emax_proj/diff1)/primEnergy/primEnergy;
00469         G4double t2=2.*std::log(Emax_proj/Emin_proj)/primEnergy/primEnergy;
00470         Cross*=(t1+t2);
00471 
00472         
00473   }
00474   lastCS =Cross;
00475   return Cross; 
00476 }
00478 //
00479 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
00480 { 
00481   G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
00482   return Tmax;
00483 }
00485 //
00486 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
00487 { return PrimAdjEnergy+Tcut;
00488 }
00490 //                              
00491 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
00492 { return HighEnergyLimit;
00493 }
00495 //
00496 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
00497 {  G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
00498    return Tmin;
00499 }