Geant4: G4AdjointIonIonisationModel.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: G4AdjointIonIonisationModel.cc 69844 2013-05-16 09:19:33Z gcosmo $
00027 //
00028 #include "G4AdjointIonIonisationModel.hh"
00029 #include "G4AdjointCSManager.hh"
00030 
00031 #include "G4PhysicalConstants.hh"
00032 #include "G4SystemOfUnits.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 "G4BraggIonModel.hh"
00041 #include "G4Proton.hh"
00042 #include "G4GenericIon.hh"
00043 #include "G4NistManager.hh"
00044 
00046 //
00047 G4AdjointIonIonisationModel::G4AdjointIonIonisationModel():
00048   G4VEmAdjointModel("Adjoint_IonIonisation")
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   use_only_bragg = false; // for the Ion ionisation using the parametrised table model the cross sections and the sample of secondaries is done
00059                                 // as in the BraggIonModel, Therefore the use of this flag;  
00060   
00061   //The direct EM Model is taken has BetheBloch it is only used for the computation 
00062   // of the differential cross section.
00063   //The Bragg model could be used  as an alternative as  it offers the same differential cross section
00064   
00065   theBetheBlochDirectEMModel = new G4BetheBlochModel(G4GenericIon::GenericIon());
00066   theBraggIonDirectEMModel = new G4BraggIonModel(G4GenericIon::GenericIon()); 
00067   theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
00068   theDirectPrimaryPartDef =0;
00069   theAdjEquivOfDirectPrimPartDef =0;
00070  /* theDirectPrimaryPartDef =fwd_ion;
00071   theAdjEquivOfDirectPrimPartDef =adj_ion;
00072  
00073   DefineProjectileProperty();*/
00074 
00075 
00076 
00077 
00078 }
00080 //
00081 G4AdjointIonIonisationModel::~G4AdjointIonIonisationModel()
00082 {;}
00084 //
00085 void G4AdjointIonIonisationModel::SampleSecondaries(const G4Track& aTrack,
00086                                 G4bool IsScatProjToProjCase,
00087                                 G4ParticleChange* fParticleChange)
00088 { 
00089  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
00090   
00091  //Elastic inverse scattering 
00092  //---------------------------------------------------------
00093  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
00094  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
00095 
00096  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
00097         return;
00098  }
00099  
00100  //Sample secondary energy
00101  //-----------------------
00102  G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
00103  CorrectPostStepWeight(fParticleChange, 
00104                               aTrack.GetWeight(), 
00105                               adjointPrimKinEnergy,
00106                               projectileKinEnergy,
00107                               IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
00108 
00109                  
00110  //Kinematic: 
00111  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest  and gives
00112  // him part of its  energy
00113  //----------------------------------------------------------------------------------------
00114  
00115  G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
00116  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
00117  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
00118  
00119  
00120  
00121  //Companion
00122  //-----------
00123  G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
00124  if (IsScatProjToProjCase) {
00125         companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
00126  } 
00127  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
00128  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;    
00129  
00130  
00131  //Projectile momentum
00132  //--------------------
00133  G4double  P_parallel = (adjointPrimP*adjointPrimP +  projectileP2 - companionP2)/(2.*adjointPrimP);
00134  G4double  P_perp = std::sqrt( projectileP2 -  P_parallel*P_parallel);
00135  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
00136  G4double phi =G4UniformRand()*2.*3.1415926;
00137  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
00138  projectileMomentum.rotateUz(dir_parallel);
00139  
00140  
00141  
00142  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
00143         fParticleChange->ProposeTrackStatus(fStopAndKill);
00144         fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
00145         //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
00146  }
00147  else {
00148         fParticleChange->ProposeEnergy(projectileKinEnergy);
00149         fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
00150  }      
00151         
00152 
00153   
00154 
00155 }
00156                 
00158 //
00159 G4double G4AdjointIonIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
00160                                       G4double kinEnergyProj, 
00161                                       G4double kinEnergyProd,
00162                                       G4double Z, 
00163                                       G4double A)
00164 {//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
00165   
00166 
00167  
00168  G4double dSigmadEprod=0;
00169  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
00170  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
00171  
00172  G4double kinEnergyProjScaled = massRatio*kinEnergyProj;
00173  
00174  
00175  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile 
00176         G4double Tmax=kinEnergyProj;
00177         
00178         G4double E1=kinEnergyProd;
00179         G4double E2=kinEnergyProd*1.000001;
00180         G4double dE=(E2-E1);
00181         G4double sigma1,sigma2;
00182         theDirectEMModel =theBraggIonDirectEMModel;
00183         if (kinEnergyProjScaled >2.*MeV && !use_only_bragg) theDirectEMModel = theBetheBlochDirectEMModel; //Bethe Bloch Model
00184         sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
00185         sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
00186         
00187         dSigmadEprod=(sigma1-sigma2)/dE;
00188         
00189         //G4double chargeSqRatio =  currentModel->GetChargeSquareRatio(theDirectPrimaryPartDef,currentMaterial,E);
00190         
00191         
00192         
00193         if (dSigmadEprod>1.) {
00194                 G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
00195                 G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
00196                 G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
00197                 
00198         }
00199         
00200          
00201         
00202          
00203          
00204          
00205          if (theDirectEMModel == theBetheBlochDirectEMModel ){
00206          //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
00207          //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
00208          //to test the rejection of a secondary
00209          //-------------------------
00210          
00211          //Source code taken from    G4BetheBlochModel::SampleSecondaries
00212          
00213                 G4double deltaKinEnergy = kinEnergyProd;
00214          
00215          //Part of the taken code
00216          //----------------------
00217          
00218          
00219          
00220          // projectile formfactor - suppresion of high energy
00221          // delta-electron production at high energy
00222          
00223          
00224                 G4double x = formfact*deltaKinEnergy;
00225                 if(x > 1.e-6) {
00226                         G4double totEnergy     = kinEnergyProj + mass;
00227                         G4double etot2         = totEnergy*totEnergy;
00228                         G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
00229                         G4double f;
00230                         G4double f1 = 0.0;
00231                         f = 1.0 - beta2*deltaKinEnergy/Tmax;
00232                         if( 0.5 == spin ) {
00233                                 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
00234                                 f += f1;
00235                         }
00236                         G4double x1 = 1.0 + x;
00237                         G4double gg  = 1.0/(x1*x1);
00238                         if( 0.5 == spin ) {
00239                                 G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
00240                                 gg *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
00241                         }
00242                         if(gg > 1.0) {
00243                                 G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: gg= " << gg
00244                                 << G4endl;
00245                                 gg=1.;
00246                         }
00247                         //G4cout<<"gg"<<gg<<G4endl;
00248                         dSigmadEprod*=gg;               
00249                 }
00250         }       
00251          
00252   }
00253 
00254  return dSigmadEprod;   
00255 }
00257 //
00258 void G4AdjointIonIonisationModel::SetIon(G4ParticleDefinition* adj_ion, G4ParticleDefinition* fwd_ion)
00259 { theDirectPrimaryPartDef =fwd_ion;
00260   theAdjEquivOfDirectPrimPartDef =adj_ion;
00261  
00262   DefineProjectileProperty();
00263 }
00265 //
00266 void G4AdjointIonIonisationModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange, G4double old_weight,
00267                                                         G4double adjointPrimKinEnergy, G4double projectileKinEnergy,G4bool )
00268 {
00269  //It is needed because the direct cross section used to compute the differential cross section is not the one used in 
00270  // the direct model where the GenericIon stuff is considered with correction of effective charge.  In the direct model the samnepl of secondaries does
00271  // not reflect the integral cross section. The integral fwd cross section that we used to compute the differential CS
00272  // match the sample of secondaries in the forward case despite the fact that its is not the same total CS than in the FWD case. For this reasion an extra
00273  // weight correction is needed at the end.
00274  
00275 
00276  G4double new_weight=old_weight;
00277  
00278  //the correction of CS due to the problem explained above
00279  G4double kinEnergyProjScaled = massRatio*projectileKinEnergy;
00280  theDirectEMModel =theBraggIonDirectEMModel;
00281  if (kinEnergyProjScaled >2.*MeV && !use_only_bragg) theDirectEMModel = theBetheBlochDirectEMModel; //Bethe Bloch Model
00282  G4double UsedFwdCS=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,projectileKinEnergy,1,1 ,currentTcutForDirectSecond,1.e20);
00283  G4double chargeSqRatio =1.;
00284  if (chargeSquare>1.) chargeSqRatio =  theDirectEMModel->GetChargeSquareRatio(theDirectPrimaryPartDef,currentMaterial,projectileKinEnergy);
00285  G4double  CorrectFwdCS = chargeSqRatio*theDirectEMModel->ComputeCrossSectionPerAtom(G4GenericIon::GenericIon(),kinEnergyProjScaled,1,1 ,currentTcutForDirectSecond,1.e20);
00286  if (UsedFwdCS >0)  new_weight*= CorrectFwdCS/UsedFwdCS;//May be some check is needed if UsedFwdCS ==0 probably that then we should avoid a secondary to be produced, 
00287  
00288  
00289  //additional CS crorrection  needed for cross section biasing in general. 
00290  //May be wrong for ions!!! Most of the time not used!
00291   G4double w_corr =1./CS_biasing_factor;
00292   w_corr*=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
00293   new_weight*=w_corr;
00294  
00295   new_weight*=projectileKinEnergy/adjointPrimKinEnergy;
00296  
00297  fParticleChange->SetParentWeightByProcess(false);
00298  fParticleChange->SetSecondaryWeightByProcess(false);
00299  fParticleChange->ProposeParentWeight(new_weight);
00300 }
00301 
00302 
00304 //
00305 void G4AdjointIonIonisationModel::DefineProjectileProperty()
00306 {   
00307     //Slightly modified code taken from G4BetheBlochModel::SetParticle
00308     //------------------------------------------------
00309     G4String pname = theDirectPrimaryPartDef->GetParticleName();
00310     if (theDirectPrimaryPartDef->GetParticleType() == "nucleus" &&
00311         pname != "deuteron" && pname != "triton") {
00312       isIon = true;
00313     }
00314     
00315     mass = theDirectPrimaryPartDef->GetPDGMass();
00316     massRatio= G4GenericIon::GenericIon()->GetPDGMass()/mass;
00317     mass_ratio_projectile = massRatio;
00318     spin = theDirectPrimaryPartDef->GetPDGSpin();
00319     G4double q = theDirectPrimaryPartDef->GetPDGCharge()/eplus;
00320     chargeSquare = q*q;
00321     ratio = electron_mass_c2/mass;
00322     ratio2 = ratio*ratio;
00323     one_plus_ratio_2=(1+ratio)*(1+ratio);
00324     one_minus_ratio_2=(1-ratio)*(1-ratio);
00325     G4double magmom = theDirectPrimaryPartDef->GetPDGMagneticMoment()
00326       *mass/(0.5*eplus*hbar_Planck*c_squared);
00327     magMoment2 = magmom*magmom - 1.0;
00328     formfact = 0.0;
00329     if(theDirectPrimaryPartDef->GetLeptonNumber() == 0) {
00330       G4double x = 0.8426*GeV;
00331       if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
00332       else if(mass > GeV) {
00333         x /= G4NistManager::Instance()->GetZ13(mass/proton_mass_c2);
00334         //      tlimit = 51.2*GeV*A13[iz]*A13[iz];
00335       }
00336       formfact = 2.0*electron_mass_c2/(x*x);
00337       tlimit   = 2.0/formfact;
00338    }
00339 }
00340 
00341 
00343 //
00344 G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
00345 { 
00346   G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
00347   return Tmax;
00348 }
00350 //
00351 G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
00352 { return PrimAdjEnergy+Tcut;
00353 }
00355 //                              
00356 G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
00357 { return HighEnergyLimit;
00358 }
00360 //
00361 G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
00362 {  G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
00363    return Tmin;
00364 }