ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4AdjointhIonisationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4AdjointhIonisationModel.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 //
27 #include "G4AdjointhIonisationModel.hh"
28 
29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4AdjointCSManager.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 "G4BraggModel.hh"
40 #include "G4Proton.hh"
41 #include "G4NistManager.hh"
42 
44 //
45 G4AdjointhIonisationModel::G4AdjointhIonisationModel(G4ParticleDefinition* projectileDefinition):
46  G4VEmAdjointModel("Adjoint_hIonisation")
47 {
48 
49 
50 
51  UseMatrix =true;
52  UseMatrixPerElement = true;
53  ApplyCutInRange = true;
54  UseOnlyOneMatrixForAllElements = true;
55  CS_biasing_factor =1.;
56  second_part_of_same_type =false;
57 
58  //The direct EM Modfel is taken has BetheBloch it is only used for the computation
59  // of the differential cross section.
60  //The Bragg model could be used as an alternative as it offers the same differential cross section
61 
62  theDirectEMModel = new G4BetheBlochModel(projectileDefinition);
63  theBraggDirectEMModel = new G4BraggModel(projectileDefinition);
64  theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
65 
66  theDirectPrimaryPartDef = projectileDefinition;
67  theAdjEquivOfDirectPrimPartDef = 0;
68  if (projectileDefinition == G4Proton::Proton()) {
69  theAdjEquivOfDirectPrimPartDef = G4AdjointProton::AdjointProton();
70 
71  }
72 
73  DefineProjectileProperty();
74 }
76 //
77 G4AdjointhIonisationModel::~G4AdjointhIonisationModel()
78 {;}
79 
80 
82 //
83 void G4AdjointhIonisationModel::SampleSecondaries(const G4Track& aTrack,
84  G4bool IsScatProjToProjCase,
85  G4ParticleChange* fParticleChange)
86 {
87  if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,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 
115  G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
116  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
117  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
118 
119 
120 
121  //Companion
122  //-----------
123  G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
124  if (IsScatProjToProjCase) {
125  companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
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 //
159 void G4AdjointhIonisationModel::RapidSampleSecondaries(const G4Track& aTrack,
160  G4bool IsScatProjToProjCase,
161  G4ParticleChange* fParticleChange)
162 {
163 
164  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
165  DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
166 
167 
168  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
169  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
170 
171  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
172  return;
173  }
174 
175  G4double projectileKinEnergy =0.;
176  G4double eEnergy=0.;
177  G4double newCS=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
178  if (!IsScatProjToProjCase){//1/E^2 distribution
179 
180  eEnergy=adjointPrimKinEnergy;
181  G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
182  G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);
183  if (Emin>=Emax) return;
184  G4double a=1./Emax;
185  G4double b=1./Emin;
186  newCS=newCS*(b-a)/eEnergy;
187 
188  projectileKinEnergy =1./(b- (b-a)*G4UniformRand());
189 
190 
191  }
192  else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
193  G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
194  if (Emin>=Emax) return;
195  G4double diff1=Emin-adjointPrimKinEnergy;
196  G4double diff2=Emax-adjointPrimKinEnergy;
197 
198  G4double t1=adjointPrimKinEnergy*(1./diff1-1./diff2);
199  G4double t2=adjointPrimKinEnergy*(1./Emin-1./Emax);
200  /*G4double f31=diff1/Emin;
201  G4double f32=diff2/Emax/f31;*/
202  G4double t3=2.*std::log(Emax/Emin);
203  G4double sum_t=t1+t2+t3;
204  newCS=newCS*sum_t/adjointPrimKinEnergy/adjointPrimKinEnergy;
205  G4double t=G4UniformRand()*sum_t;
206  if (t <=t1 ){
207  G4double q= G4UniformRand()*t1/adjointPrimKinEnergy ;
208  projectileKinEnergy =adjointPrimKinEnergy +1./(1./diff1-q);
209 
210  }
211  else if (t <=t2 ) {
212  G4double q= G4UniformRand()*t2/adjointPrimKinEnergy;
213  projectileKinEnergy =1./(1./Emin-q);
214  }
215  else {
216  projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand());
217 
218  }
219  eEnergy=projectileKinEnergy-adjointPrimKinEnergy;
220 
221 
222  }
223 
224 
225 
226  G4double diffCS_perAtom_Used=twopi_mc2_rcl2*mass*adjointPrimKinEnergy/projectileKinEnergy/projectileKinEnergy/eEnergy/eEnergy;
227 
228 
229 
230  //Weight correction
231  //-----------------------
232  //First w_corr is set to the ratio between adjoint total CS and fwd total CS
233  G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
234 
235  //G4cout<<w_corr<<G4endl;
236  w_corr*=newCS/lastCS;
237  //G4cout<<w_corr<<G4endl;
238  //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
239  //Here we consider the true diffCS as the one obtained by the numerical differentiation over Tcut of the direct CS
240 
241  G4double diffCS = DiffCrossSectionPerAtomPrimToSecond(projectileKinEnergy, eEnergy,1,1);
242  w_corr*=diffCS/diffCS_perAtom_Used;
243  //G4cout<<w_corr<<G4endl;
244 
245  G4double new_weight = aTrack.GetWeight()*w_corr;
246  fParticleChange->SetParentWeightByProcess(false);
247  fParticleChange->SetSecondaryWeightByProcess(false);
248  fParticleChange->ProposeParentWeight(new_weight);
249 
250 
251 
252 
253  //Kinematic:
254  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
255  // him part of its energy
256  //----------------------------------------------------------------------------------------
257 
258  G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
259  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
260  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
261 
262 
263 
264  //Companion
265  //-----------
266  G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
267  if (IsScatProjToProjCase) {
268  companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
269  }
270  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
271  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
272 
273 
274  //Projectile momentum
275  //--------------------
276  G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
277  G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
278  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
279  G4double phi =G4UniformRand()*2.*3.1415926;
280  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
281  projectileMomentum.rotateUz(dir_parallel);
282 
283 
284 
285  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
286  fParticleChange->ProposeTrackStatus(fStopAndKill);
287  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
288  //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
289  }
290  else {
291  fParticleChange->ProposeEnergy(projectileKinEnergy);
292  fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
293  }
294 
295 
296 
297 
298 
299 
300 
301 }
302 
304 //
305 G4double G4AdjointhIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
306  G4double kinEnergyProj,
307  G4double kinEnergyProd,
308  G4double Z,
309  G4double A)
310 {//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
311 
312 
313 
314  G4double dSigmadEprod=0;
315  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
316  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
317 
318 
319  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
320  G4double Tmax=kinEnergyProj;
321 
322  G4double E1=kinEnergyProd;
323  G4double E2=kinEnergyProd*1.000001;
324  G4double dE=(E2-E1);
325  G4double sigma1,sigma2;
326  if (kinEnergyProj >2.*MeV){
327  sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
328  sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
329  }
330  else {
331  sigma1=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
332  sigma2=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
333  }
334 
335 
336  dSigmadEprod=(sigma1-sigma2)/dE;
337  if (dSigmadEprod>1.) {
338  G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
339  G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
340  G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
341 
342  }
343 
344 
345 
346  //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
347  //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
348  //to test the rejection of a secondary
349  //-------------------------
350 
351  //Source code taken from G4BetheBlochModel::SampleSecondaries
352 
353  G4double deltaKinEnergy = kinEnergyProd;
354 
355  //Part of the taken code
356  //----------------------
357 
358 
359 
360  // projectile formfactor - suppresion of high energy
361  // delta-electron production at high energy
362  G4double x = formfact*deltaKinEnergy;
363  if(x > 1.e-6) {
364 
365 
366  G4double totEnergy = kinEnergyProj + mass;
367  G4double etot2 = totEnergy*totEnergy;
368  G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
369  G4double f;
370  G4double f1 = 0.0;
371  f = 1.0 - beta2*deltaKinEnergy/Tmax;
372  if( 0.5 == spin ) {
373  f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
374  f += f1;
375  }
376  G4double x1 = 1.0 + x;
377  G4double gg = 1.0/(x1*x1);
378  if( 0.5 == spin ) {
379  G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
380  gg *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
381  }
382  if(gg > 1.0) {
383  G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g
384  << G4endl;
385  gg=1.;
386  }
387  //G4cout<<"gg"<<gg<<G4endl;
388  dSigmadEprod*=gg;
389  }
390 
391  }
392 
393  return dSigmadEprod;
394 }
395 
396 
397 
399 //
400 void G4AdjointhIonisationModel::DefineProjectileProperty()
401 {
402  //Slightly modified code taken from G4BetheBlochModel::SetParticle
403  //------------------------------------------------
404  G4String pname = theDirectPrimaryPartDef->GetParticleName();
405  if (theDirectPrimaryPartDef->GetParticleType() == "nucleus" &&
406  pname != "deuteron" && pname != "triton") {
407  isIon = true;
408  }
409 
410  mass = theDirectPrimaryPartDef->GetPDGMass();
411  mass_ratio_projectile = proton_mass_c2/theDirectPrimaryPartDef->GetPDGMass();;
412  spin = theDirectPrimaryPartDef->GetPDGSpin();
413  G4double q = theDirectPrimaryPartDef->GetPDGCharge()/eplus;
414  chargeSquare = q*q;
415  ratio = electron_mass_c2/mass;
416  ratio2 = ratio*ratio;
417  one_plus_ratio_2=(1+ratio)*(1+ratio);
418  one_minus_ratio_2=(1-ratio)*(1-ratio);
419  G4double magmom = theDirectPrimaryPartDef->GetPDGMagneticMoment()
420  *mass/(0.5*eplus*hbar_Planck*c_squared);
421  magMoment2 = magmom*magmom - 1.0;
422  formfact = 0.0;
423  if(theDirectPrimaryPartDef->GetLeptonNumber() == 0) {
424  G4double x = 0.8426*GeV;
425  if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
426  else if(mass > GeV) {
427  x /= G4NistManager::Instance()->GetZ13(mass/proton_mass_c2);
428  // tlimit = 51.2*GeV*A13[iz]*A13[iz];
429  }
430  formfact = 2.0*electron_mass_c2/(x*x);
431  tlimit = 2.0/formfact;
432  }
433 }
434 
436 //
437 G4double G4AdjointhIonisationModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
438  G4double primEnergy,
439  G4bool IsScatProjToProjCase)
440 {
441  if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
442  DefineCurrentMaterial(aCouple);
443 
444 
445  G4double Cross=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
446 
447  if (!IsScatProjToProjCase ){
448  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
449  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
450  if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) {
451  Cross*=(1./Emin_proj -1./Emax_proj)/primEnergy;
452  }
453  else Cross=0.;
454 
455 
456 
457 
458 
459 
460  }
461  else {
462  G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
463  G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond);
464  G4double diff1=Emin_proj-primEnergy;
465  G4double diff2=Emax_proj-primEnergy;
466  G4double t1=(1./diff1+1./Emin_proj-1./diff2-1./Emax_proj)/primEnergy;
467  //G4double t2=2.*std::log(diff2*Emin_proj/Emax_proj/diff1)/primEnergy/primEnergy;
468  G4double t2=2.*std::log(Emax_proj/Emin_proj)/primEnergy/primEnergy;
469  Cross*=(t1+t2);
470 
471 
472  }
473  lastCS =Cross;
474  return Cross;
475 }
477 //
478 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
479 {
480  G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
481  return Tmax;
482 }
484 //
485 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
486 { return PrimAdjEnergy+Tcut;
487 }
489 //
490 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
491 { return HighEnergyLimit;
492 }
494 //
495 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
496 { G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
497  return Tmin;
498 }