ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4AdjointIonIonisationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4AdjointIonIonisationModel.cc
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 //
28 #include "G4AdjointCSManager.hh"
29 
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "G4Integrator.hh"
33 #include "G4TrackStatus.hh"
34 #include "G4ParticleChange.hh"
35 #include "G4AdjointElectron.hh"
36 #include "G4AdjointProton.hh"
37 #include "G4AdjointInterpolator.hh"
38 #include "G4BetheBlochModel.hh"
39 #include "G4BraggIonModel.hh"
40 #include "G4Proton.hh"
41 #include "G4GenericIon.hh"
42 #include "G4NistManager.hh"
43 
45 //
47  G4VEmAdjointModel("Adjoint_IonIonisation")
48 {
49 
50 
51  UseMatrix =true;
52  UseMatrixPerElement = true;
53  ApplyCutInRange = true;
57  use_only_bragg = false; // for the Ion ionisation using the parametrised table model the cross sections and the sample of secondaries is done
58  // as in the BraggIonModel, Therefore the use of this flag;
59 
60  //The direct EM Model is taken has BetheBloch it is only used for the computation
61  // of the differential cross section.
62  //The Bragg model could be used as an alternative as it offers the same differential cross section
63 
69  /* theDirectPrimaryPartDef =fwd_ion;
70  theAdjEquivOfDirectPrimPartDef =adj_ion;
71 
72  DefineProjectileProperty();*/
73 
74 
75 
76 
77 }
79 //
81 {;}
83 //
85  G4bool IsScatProjToProjCase,
86  G4ParticleChange* fParticleChange)
87 {
88  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
89 
90  //Elastic inverse scattering
91  //---------------------------------------------------------
92  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
93  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
94 
95  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
96  return;
97  }
98 
99  //Sample secondary energy
100  //-----------------------
101  G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
102  CorrectPostStepWeight(fParticleChange,
103  aTrack.GetWeight(),
104  adjointPrimKinEnergy,
105  projectileKinEnergy,
106  IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
107 
108 
109  //Kinematic:
110  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
111  // him part of its energy
112  //----------------------------------------------------------------------------------------
113 
115  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
116  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
117 
118 
119 
120  //Companion
121  //-----------
123  if (IsScatProjToProjCase) {
125  }
126  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
127  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
128 
129 
130  //Projectile momentum
131  //--------------------
132  G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
133  G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
134  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
135  G4double phi =G4UniformRand()*2.*3.1415926;
136  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
137  projectileMomentum.rotateUz(dir_parallel);
138 
139 
140 
141  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
142  fParticleChange->ProposeTrackStatus(fStopAndKill);
143  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
144  //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
145  }
146  else {
147  fParticleChange->ProposeEnergy(projectileKinEnergy);
148  fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
149  }
150 
151 
152 
153 
154 }
155 
157 //
159  G4double kinEnergyProj,
160  G4double kinEnergyProd,
161  G4double Z,
162  G4double A)
163 {//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
164 
165 
166 
167  G4double dSigmadEprod=0;
168  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
169  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
170 
171  G4double kinEnergyProjScaled = massRatio*kinEnergyProj;
172 
173 
174  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
175  G4double Tmax=kinEnergyProj;
176 
177  G4double E1=kinEnergyProd;
178  G4double E2=kinEnergyProd*1.000001;
179  G4double dE=(E2-E1);
180  G4double sigma1,sigma2;
182  if (kinEnergyProjScaled >2.*MeV && !use_only_bragg) theDirectEMModel = theBetheBlochDirectEMModel; //Bethe Bloch Model
183  sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
184  sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
185 
186  dSigmadEprod=(sigma1-sigma2)/dE;
187 
188  //G4double chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPrimaryPartDef,currentMaterial,E);
189 
190 
191 
192  if (dSigmadEprod>1.) {
193  G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
194  G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
195  G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
196 
197  }
198 
199 
200 
201 
202 
203 
205  //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
206  //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
207  //to test the rejection of a secondary
208  //-------------------------
209 
210  //Source code taken from G4BetheBlochModel::SampleSecondaries
211 
212  G4double deltaKinEnergy = kinEnergyProd;
213 
214  //Part of the taken code
215  //----------------------
216 
217 
218 
219  // projectile formfactor - suppresion of high energy
220  // delta-electron production at high energy
221 
222 
223  G4double x = formfact*deltaKinEnergy;
224  if(x > 1.e-6) {
225  G4double totEnergy = kinEnergyProj + mass;
226  G4double etot2 = totEnergy*totEnergy;
227  G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
228  G4double f;
229  G4double f1 = 0.0;
230  f = 1.0 - beta2*deltaKinEnergy/Tmax;
231  if( 0.5 == spin ) {
232  f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
233  f += f1;
234  }
235  G4double x1 = 1.0 + x;
236  G4double gg = 1.0/(x1*x1);
237  if( 0.5 == spin ) {
238  G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
239  gg *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
240  }
241  if(gg > 1.0) {
242  G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: gg= " << gg
243  << G4endl;
244  gg=1.;
245  }
246  //G4cout<<"gg"<<gg<<G4endl;
247  dSigmadEprod*=gg;
248  }
249  }
250 
251  }
252 
253  return dSigmadEprod;
254 }
256 //
258 { theDirectPrimaryPartDef =fwd_ion;
260 
262 }
264 //
266  G4double adjointPrimKinEnergy, G4double projectileKinEnergy,G4bool )
267 {
268  //It is needed because the direct cross section used to compute the differential cross section is not the one used in
269  // the direct model where the GenericIon stuff is considered with correction of effective charge. In the direct model the samnepl of secondaries does
270  // not reflect the integral cross section. The integral fwd cross section that we used to compute the differential CS
271  // 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
272  // weight correction is needed at the end.
273 
274 
275  G4double new_weight=old_weight;
276 
277  //the correction of CS due to the problem explained above
278  G4double kinEnergyProjScaled = massRatio*projectileKinEnergy;
280  if (kinEnergyProjScaled >2.*MeV && !use_only_bragg) theDirectEMModel = theBetheBlochDirectEMModel; //Bethe Bloch Model
282  G4double chargeSqRatio =1.;
283  if (chargeSquare>1.) chargeSqRatio = theDirectEMModel->GetChargeSquareRatio(theDirectPrimaryPartDef,currentMaterial,projectileKinEnergy);
284  G4double CorrectFwdCS = chargeSqRatio*theDirectEMModel->ComputeCrossSectionPerAtom(G4GenericIon::GenericIon(),kinEnergyProjScaled,1,1 ,currentTcutForDirectSecond,1.e20);
285  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,
286 
287 
288  //additional CS crorrection needed for cross section biasing in general.
289  //May be wrong for ions!!! Most of the time not used!
290  G4double w_corr =1./CS_biasing_factor;
292  new_weight*=w_corr;
293 
294  new_weight*=projectileKinEnergy/adjointPrimKinEnergy;
295 
296  fParticleChange->SetParentWeightByProcess(false);
297  fParticleChange->SetSecondaryWeightByProcess(false);
298  fParticleChange->ProposeParentWeight(new_weight);
299 }
300 
301 
303 //
305 {
306  //Slightly modified code taken from G4BetheBlochModel::SetParticle
307  //------------------------------------------------
309  if (theDirectPrimaryPartDef->GetParticleType() == "nucleus" &&
310  pname != "deuteron" && pname != "triton") {
311  isIon = true;
312  }
313 
319  chargeSquare = q*q;
321  ratio2 = ratio*ratio;
322  one_plus_ratio_2=(1+ratio)*(1+ratio);
323  one_minus_ratio_2=(1-ratio)*(1-ratio);
326  magMoment2 = magmom*magmom - 1.0;
327  formfact = 0.0;
329  G4double x = 0.8426*GeV;
330  if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
331  else if(mass > GeV) {
333  // tlimit = 51.2*GeV*A13[iz]*A13[iz];
334  }
335  formfact = 2.0*electron_mass_c2/(x*x);
336  tlimit = 2.0/formfact;
337  }
338 }
339 
340 
342 //
344 {
345  G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
346  return Tmax;
347 }
349 //
351 { return PrimAdjEnergy+Tcut;
352 }
354 //
356 { return HighEnergyLimit;
357 }
359 //
361 { G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
362  return Tmin;
363 }