Geant4  10.02.p01
G4AdjointIonIonisationModel.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4AdjointIonIonisationModel.cc 66892 2013-01-17 10:57:59Z gunter $
27 //
29 #include "G4AdjointCSManager.hh"
30 
31 #include "G4PhysicalConstants.hh"
32 #include "G4SystemOfUnits.hh"
33 #include "G4Integrator.hh"
34 #include "G4TrackStatus.hh"
35 #include "G4ParticleChange.hh"
36 #include "G4AdjointElectron.hh"
37 #include "G4AdjointProton.hh"
38 #include "G4AdjointInterpolator.hh"
39 #include "G4BetheBlochModel.hh"
40 #include "G4BraggIonModel.hh"
41 #include "G4Proton.hh"
42 #include "G4GenericIon.hh"
43 #include "G4NistManager.hh"
44 
46 //
48  G4VEmAdjointModel("Adjoint_IonIonisation")
49 {
50 
51 
52  UseMatrix =true;
53  UseMatrixPerElement = true;
54  ApplyCutInRange = true;
58  use_only_bragg = false; // for the Ion ionisation using the parametrised table model the cross sections and the sample of secondaries is done
59  // as in the BraggIonModel, Therefore the use of this flag;
60 
61  //The direct EM Model is taken has BetheBloch it is only used for the computation
62  // of the differential cross section.
63  //The Bragg model could be used as an alternative as it offers the same differential cross section
64 
70  /* theDirectPrimaryPartDef =fwd_ion;
71  theAdjEquivOfDirectPrimPartDef =adj_ion;
72 
73  DefineProjectileProperty();*/
74 
75 
76 
77 
78 }
80 //
82 {;}
84 //
86  G4bool IsScatProjToProjCase,
87  G4ParticleChange* fParticleChange)
88 {
89  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
90 
91  //Elastic inverse scattering
92  //---------------------------------------------------------
93  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
94  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
95 
96  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
97  return;
98  }
99 
100  //Sample secondary energy
101  //-----------------------
102  G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
103  CorrectPostStepWeight(fParticleChange,
104  aTrack.GetWeight(),
105  adjointPrimKinEnergy,
106  projectileKinEnergy,
107  IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
108 
109 
110  //Kinematic:
111  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
112  // him part of its energy
113  //----------------------------------------------------------------------------------------
114 
116  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
117  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
118 
119 
120 
121  //Companion
122  //-----------
124  if (IsScatProjToProjCase) {
126  }
127  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
128  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
129 
130 
131  //Projectile momentum
132  //--------------------
133  G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
134  G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
135  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
136  G4double phi =G4UniformRand()*2.*3.1415926;
137  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
138  projectileMomentum.rotateUz(dir_parallel);
139 
140 
141 
142  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
143  fParticleChange->ProposeTrackStatus(fStopAndKill);
144  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
145  //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
146  }
147  else {
148  fParticleChange->ProposeEnergy(projectileKinEnergy);
149  fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
150  }
151 
152 
153 
154 
155 }
156 
158 //
160  G4double kinEnergyProj,
161  G4double kinEnergyProd,
162  G4double Z,
163  G4double A)
164 {//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
165 
166 
167 
168  G4double dSigmadEprod=0;
169  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
170  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
171 
172  G4double kinEnergyProjScaled = massRatio*kinEnergyProj;
173 
174 
175  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
176  G4double Tmax=kinEnergyProj;
177 
178  G4double E1=kinEnergyProd;
179  G4double E2=kinEnergyProd*1.000001;
180  G4double dE=(E2-E1);
181  G4double sigma1,sigma2;
183  if (kinEnergyProjScaled >2.*MeV && !use_only_bragg) theDirectEMModel = theBetheBlochDirectEMModel; //Bethe Bloch Model
184  sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
185  sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
186 
187  dSigmadEprod=(sigma1-sigma2)/dE;
188 
189  //G4double chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPrimaryPartDef,currentMaterial,E);
190 
191 
192 
193  if (dSigmadEprod>1.) {
194  G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
195  G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
196  G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
197 
198  }
199 
200 
201 
202 
203 
204 
206  //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
207  //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
208  //to test the rejection of a secondary
209  //-------------------------
210 
211  //Source code taken from G4BetheBlochModel::SampleSecondaries
212 
213  G4double deltaKinEnergy = kinEnergyProd;
214 
215  //Part of the taken code
216  //----------------------
217 
218 
219 
220  // projectile formfactor - suppresion of high energy
221  // delta-electron production at high energy
222 
223 
224  G4double x = formfact*deltaKinEnergy;
225  if(x > 1.e-6) {
226  G4double totEnergy = kinEnergyProj + mass;
227  G4double etot2 = totEnergy*totEnergy;
228  G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
229  G4double f;
230  G4double f1 = 0.0;
231  f = 1.0 - beta2*deltaKinEnergy/Tmax;
232  if( 0.5 == spin ) {
233  f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
234  f += f1;
235  }
236  G4double x1 = 1.0 + x;
237  G4double gg = 1.0/(x1*x1);
238  if( 0.5 == spin ) {
239  G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
240  gg *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
241  }
242  if(gg > 1.0) {
243  G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: gg= " << gg
244  << G4endl;
245  gg=1.;
246  }
247  //G4cout<<"gg"<<gg<<G4endl;
248  dSigmadEprod*=gg;
249  }
250  }
251 
252  }
253 
254  return dSigmadEprod;
255 }
257 //
259 { theDirectPrimaryPartDef =fwd_ion;
261 
263 }
265 //
267  G4double adjointPrimKinEnergy, G4double projectileKinEnergy,G4bool )
268 {
269  //It is needed because the direct cross section used to compute the differential cross section is not the one used in
270  // the direct model where the GenericIon stuff is considered with correction of effective charge. In the direct model the samnepl of secondaries does
271  // not reflect the integral cross section. The integral fwd cross section that we used to compute the differential CS
272  // 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
273  // weight correction is needed at the end.
274 
275 
276  G4double new_weight=old_weight;
277 
278  //the correction of CS due to the problem explained above
279  G4double kinEnergyProjScaled = massRatio*projectileKinEnergy;
281  if (kinEnergyProjScaled >2.*MeV && !use_only_bragg) theDirectEMModel = theBetheBlochDirectEMModel; //Bethe Bloch Model
283  G4double chargeSqRatio =1.;
284  if (chargeSquare>1.) chargeSqRatio = theDirectEMModel->GetChargeSquareRatio(theDirectPrimaryPartDef,currentMaterial,projectileKinEnergy);
285  G4double CorrectFwdCS = chargeSqRatio*theDirectEMModel->ComputeCrossSectionPerAtom(G4GenericIon::GenericIon(),kinEnergyProjScaled,1,1 ,currentTcutForDirectSecond,1.e20);
286  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,
287 
288 
289  //additional CS crorrection needed for cross section biasing in general.
290  //May be wrong for ions!!! Most of the time not used!
291  G4double w_corr =1./CS_biasing_factor;
293  new_weight*=w_corr;
294 
295  new_weight*=projectileKinEnergy/adjointPrimKinEnergy;
296 
297  fParticleChange->SetParentWeightByProcess(false);
298  fParticleChange->SetSecondaryWeightByProcess(false);
299  fParticleChange->ProposeParentWeight(new_weight);
300 }
301 
302 
304 //
306 {
307  //Slightly modified code taken from G4BetheBlochModel::SetParticle
308  //------------------------------------------------
310  if (theDirectPrimaryPartDef->GetParticleType() == "nucleus" &&
311  pname != "deuteron" && pname != "triton") {
312  isIon = true;
313  }
314 
320  chargeSquare = q*q;
321  ratio = electron_mass_c2/mass;
322  ratio2 = ratio*ratio;
323  one_plus_ratio_2=(1+ratio)*(1+ratio);
324  one_minus_ratio_2=(1-ratio)*(1-ratio);
326  *mass/(0.5*eplus*hbar_Planck*c_squared);
327  magMoment2 = magmom*magmom - 1.0;
328  formfact = 0.0;
330  G4double x = 0.8426*GeV;
331  if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
332  else if(mass > GeV) {
333  x /= G4NistManager::Instance()->GetZ13(mass/proton_mass_c2);
334  // tlimit = 51.2*GeV*A13[iz]*A13[iz];
335  }
336  formfact = 2.0*electron_mass_c2/(x*x);
337  tlimit = 2.0/formfact;
338  }
339 }
340 
341 
343 //
345 {
346  G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
347  return Tmax;
348 }
350 //
352 { return PrimAdjEnergy+Tcut;
353 }
355 //
357 { return HighEnergyLimit;
358 }
360 //
362 { G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
363  return Tmin;
364 }
static const double MeV
Definition: G4SIunits.hh:211
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetPostStepWeightCorrection()
const G4DynamicParticle * GetDynamicParticle() const
G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex, G4double prim_energy, G4bool IsScatProjToProjCase)
G4ParticleDefinition * theDirectPrimaryPartDef
virtual G4double DiffCrossSectionPerAtomPrimToSecond(G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
void ProposeParentWeight(G4double finalWeight)
G4double GetZ13(G4double Z)
static G4AdjointElectron * AdjointElectron()
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)
G4Material * currentMaterial
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
static G4NistManager * Instance()
const G4String & GetParticleName() const
G4double GetTotalMomentum() const
void SetSecondaryWeightByProcess(G4bool)
void SetParentWeightByProcess(G4bool)
static const G4double dE
G4VEmModel * theDirectEMModel
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
void SetIon(G4ParticleDefinition *adj_ion, G4ParticleDefinition *fwd_ion)
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
static const double GeV
Definition: G4SIunits.hh:214
const G4String & GetParticleType() const
static const G4double f1
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:313
G4bool UseOnlyOneMatrixForAllElements
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:345
static G4GenericIon * GenericIon()
Definition: G4GenericIon.cc:93
virtual void SampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
G4double GetPDGMass() const
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
const G4double x[NPOINTSGL]
void ProposeEnergy(G4double finalEnergy)
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
void AddSecondary(G4Track *aSecondary)
G4double GetWeight() const
G4double GetPDGSpin() const
#define G4endl
Definition: G4ios.hh:61
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
G4double currentTcutForDirectSecond
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4double GetPDGMagneticMoment() const
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
static const double eplus
Definition: G4SIunits.hh:196
G4double GetPDGCharge() const
static G4AdjointCSManager * GetAdjointCSManager()