Geant4  10.03
G4AdjointhIonisationModel.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: G4AdjointhIonisationModel.cc 66892 2013-01-17 10:57:59Z gunter $
27 //
29 
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "G4AdjointCSManager.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 "G4BraggModel.hh"
41 #include "G4Proton.hh"
42 #include "G4NistManager.hh"
43 
45 //
47  G4VEmAdjointModel("Adjoint_hIonisation")
48 {
49 
50 
51 
52  UseMatrix =true;
53  UseMatrixPerElement = true;
54  ApplyCutInRange = true;
58 
59  //The direct EM Modfel is taken has BetheBloch it is only used for the computation
60  // of the differential cross section.
61  //The Bragg model could be used as an alternative as it offers the same differential cross section
62 
63  theDirectEMModel = new G4BetheBlochModel(projectileDefinition);
64  theBraggDirectEMModel = new G4BraggModel(projectileDefinition);
66 
67  theDirectPrimaryPartDef = projectileDefinition;
69  if (projectileDefinition == G4Proton::Proton()) {
71 
72  }
73 
75 }
77 //
79 {;}
80 
81 
83 //
85  G4bool IsScatProjToProjCase,
86  G4ParticleChange* fParticleChange)
87 {
88  if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
89 
90  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
91 
92  //Elastic inverse scattering
93  //---------------------------------------------------------
94  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
95  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
96 
97  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
98  return;
99  }
100 
101  //Sample secondary energy
102  //-----------------------
103  G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
104  CorrectPostStepWeight(fParticleChange,
105  aTrack.GetWeight(),
106  adjointPrimKinEnergy,
107  projectileKinEnergy,
108  IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
109 
110 
111  //Kinematic:
112  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
113  // him part of its energy
114  //----------------------------------------------------------------------------------------
115 
117  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
118  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
119 
120 
121 
122  //Companion
123  //-----------
125  if (IsScatProjToProjCase) {
127  }
128  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
129  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
130 
131 
132  //Projectile momentum
133  //--------------------
134  G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
135  G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
136  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
137  G4double phi =G4UniformRand()*2.*3.1415926;
138  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
139  projectileMomentum.rotateUz(dir_parallel);
140 
141 
142 
143  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
144  fParticleChange->ProposeTrackStatus(fStopAndKill);
145  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
146  //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
147  }
148  else {
149  fParticleChange->ProposeEnergy(projectileKinEnergy);
150  fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
151  }
152 
153 
154 
155 
156 }
157 
159 //
161  G4bool IsScatProjToProjCase,
162  G4ParticleChange* fParticleChange)
163 {
164 
165  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
167 
168 
169  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
170  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
171 
172  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
173  return;
174  }
175 
176  G4double projectileKinEnergy =0.;
177  G4double eEnergy=0.;
178  G4double newCS=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
179  if (!IsScatProjToProjCase){//1/E^2 distribution
180 
181  eEnergy=adjointPrimKinEnergy;
182  G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
184  if (Emin>=Emax) return;
185  G4double a=1./Emax;
186  G4double b=1./Emin;
187  newCS=newCS*(b-a)/eEnergy;
188 
189  projectileKinEnergy =1./(b- (b-a)*G4UniformRand());
190 
191 
192  }
193  else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
195  if (Emin>=Emax) return;
196  G4double diff1=Emin-adjointPrimKinEnergy;
197  G4double diff2=Emax-adjointPrimKinEnergy;
198 
199  G4double t1=adjointPrimKinEnergy*(1./diff1-1./diff2);
200  G4double t2=adjointPrimKinEnergy*(1./Emin-1./Emax);
201  /*G4double f31=diff1/Emin;
202  G4double f32=diff2/Emax/f31;*/
203  G4double t3=2.*std::log(Emax/Emin);
204  G4double sum_t=t1+t2+t3;
205  newCS=newCS*sum_t/adjointPrimKinEnergy/adjointPrimKinEnergy;
206  G4double t=G4UniformRand()*sum_t;
207  if (t <=t1 ){
208  G4double q= G4UniformRand()*t1/adjointPrimKinEnergy ;
209  projectileKinEnergy =adjointPrimKinEnergy +1./(1./diff1-q);
210 
211  }
212  else if (t <=t2 ) {
213  G4double q= G4UniformRand()*t2/adjointPrimKinEnergy;
214  projectileKinEnergy =1./(1./Emin-q);
215  }
216  else {
217  projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand());
218 
219  }
220  eEnergy=projectileKinEnergy-adjointPrimKinEnergy;
221 
222 
223  }
224 
225 
226 
227  G4double diffCS_perAtom_Used=twopi_mc2_rcl2*mass*adjointPrimKinEnergy/projectileKinEnergy/projectileKinEnergy/eEnergy/eEnergy;
228 
229 
230 
231  //Weight correction
232  //-----------------------
233  //First w_corr is set to the ratio between adjoint total CS and fwd total CS
235 
236  //G4cout<<w_corr<<G4endl;
237  w_corr*=newCS/lastCS;
238  //G4cout<<w_corr<<G4endl;
239  //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
240  //Here we consider the true diffCS as the one obtained by the numerical differentiation over Tcut of the direct CS
241 
242  G4double diffCS = DiffCrossSectionPerAtomPrimToSecond(projectileKinEnergy, eEnergy,1,1);
243  w_corr*=diffCS/diffCS_perAtom_Used;
244  //G4cout<<w_corr<<G4endl;
245 
246  G4double new_weight = aTrack.GetWeight()*w_corr;
247  fParticleChange->SetParentWeightByProcess(false);
248  fParticleChange->SetSecondaryWeightByProcess(false);
249  fParticleChange->ProposeParentWeight(new_weight);
250 
251 
252 
253 
254  //Kinematic:
255  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
256  // him part of its energy
257  //----------------------------------------------------------------------------------------
258 
260  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
261  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
262 
263 
264 
265  //Companion
266  //-----------
268  if (IsScatProjToProjCase) {
270  }
271  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
272  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
273 
274 
275  //Projectile momentum
276  //--------------------
277  G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
278  G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
279  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
280  G4double phi =G4UniformRand()*2.*3.1415926;
281  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
282  projectileMomentum.rotateUz(dir_parallel);
283 
284 
285 
286  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
287  fParticleChange->ProposeTrackStatus(fStopAndKill);
288  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
289  //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
290  }
291  else {
292  fParticleChange->ProposeEnergy(projectileKinEnergy);
293  fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
294  }
295 
296 
297 
298 
299 
300 
301 
302 }
303 
305 //
307  G4double kinEnergyProj,
308  G4double kinEnergyProd,
309  G4double Z,
310  G4double A)
311 {//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
312 
313 
314 
315  G4double dSigmadEprod=0;
316  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
317  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
318 
319 
320  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
321  G4double Tmax=kinEnergyProj;
322 
323  G4double E1=kinEnergyProd;
324  G4double E2=kinEnergyProd*1.000001;
325  G4double dE=(E2-E1);
326  G4double sigma1,sigma2;
327  if (kinEnergyProj >2.*MeV){
328  sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
329  sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
330  }
331  else {
332  sigma1=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
333  sigma2=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
334  }
335 
336 
337  dSigmadEprod=(sigma1-sigma2)/dE;
338  if (dSigmadEprod>1.) {
339  G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
340  G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
341  G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
342 
343  }
344 
345 
346 
347  //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
348  //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
349  //to test the rejection of a secondary
350  //-------------------------
351 
352  //Source code taken from G4BetheBlochModel::SampleSecondaries
353 
354  G4double deltaKinEnergy = kinEnergyProd;
355 
356  //Part of the taken code
357  //----------------------
358 
359 
360 
361  // projectile formfactor - suppresion of high energy
362  // delta-electron production at high energy
363  G4double x = formfact*deltaKinEnergy;
364  if(x > 1.e-6) {
365 
366 
367  G4double totEnergy = kinEnergyProj + mass;
368  G4double etot2 = totEnergy*totEnergy;
369  G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
370  G4double f;
371  G4double f1 = 0.0;
372  f = 1.0 - beta2*deltaKinEnergy/Tmax;
373  if( 0.5 == spin ) {
374  f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
375  f += f1;
376  }
377  G4double x1 = 1.0 + x;
378  G4double gg = 1.0/(x1*x1);
379  if( 0.5 == spin ) {
380  G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
381  gg *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
382  }
383  if(gg > 1.0) {
384  G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g
385  << G4endl;
386  gg=1.;
387  }
388  //G4cout<<"gg"<<gg<<G4endl;
389  dSigmadEprod*=gg;
390  }
391 
392  }
393 
394  return dSigmadEprod;
395 }
396 
397 
398 
400 //
402 {
403  //Slightly modified code taken from G4BetheBlochModel::SetParticle
404  //------------------------------------------------
406  if (theDirectPrimaryPartDef->GetParticleType() == "nucleus" &&
407  pname != "deuteron" && pname != "triton") {
408  isIon = true;
409  }
410 
415  chargeSquare = q*q;
416  ratio = electron_mass_c2/mass;
417  ratio2 = ratio*ratio;
418  one_plus_ratio_2=(1+ratio)*(1+ratio);
419  one_minus_ratio_2=(1-ratio)*(1-ratio);
421  *mass/(0.5*eplus*hbar_Planck*c_squared);
422  magMoment2 = magmom*magmom - 1.0;
423  formfact = 0.0;
425  G4double x = 0.8426*GeV;
426  if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
427  else if(mass > GeV) {
428  x /= G4NistManager::Instance()->GetZ13(mass/proton_mass_c2);
429  // tlimit = 51.2*GeV*A13[iz]*A13[iz];
430  }
431  formfact = 2.0*electron_mass_c2/(x*x);
432  tlimit = 2.0/formfact;
433  }
434 }
435 
437 //
439  G4double primEnergy,
440  G4bool IsScatProjToProjCase)
441 {
442  if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
443  DefineCurrentMaterial(aCouple);
444 
445 
446  G4double Cross=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
447 
448  if (!IsScatProjToProjCase ){
449  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
450  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
451  if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) {
452  Cross*=(1./Emin_proj -1./Emax_proj)/primEnergy;
453  }
454  else Cross=0.;
455 
456 
457 
458 
459 
460 
461  }
462  else {
465  G4double diff1=Emin_proj-primEnergy;
466  G4double diff2=Emax_proj-primEnergy;
467  G4double t1=(1./diff1+1./Emin_proj-1./diff2-1./Emax_proj)/primEnergy;
468  //G4double t2=2.*std::log(diff2*Emin_proj/Emax_proj/diff1)/primEnergy/primEnergy;
469  G4double t2=2.*std::log(Emax_proj/Emin_proj)/primEnergy/primEnergy;
470  Cross*=(t1+t2);
471 
472 
473  }
474  lastCS =Cross;
475  return Cross;
476 }
478 //
480 {
481  G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
482  return Tmax;
483 }
485 //
487 { return PrimAdjEnergy+Tcut;
488 }
490 //
492 { return HighEnergyLimit;
493 }
495 //
497 { G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
498  return Tmin;
499 }
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetPostStepWeightCorrection()
const G4DynamicParticle * GetDynamicParticle() const
G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex, G4double prim_energy, G4bool IsScatProjToProjCase)
G4ParticleDefinition * theDirectPrimaryPartDef
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
void ProposeParentWeight(G4double finalWeight)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
void RapidSampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
static G4AdjointElectron * AdjointElectron()
static constexpr double g
Definition: G4SIunits.hh:183
G4Material * currentMaterial
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
static G4NistManager * Instance()
const G4String & GetParticleName() const
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
G4double GetTotalMomentum() const
G4double GetZ13(G4double Z) 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)
G4double GetElectronDensity() const
Definition: G4Material.hh:217
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
static constexpr double eplus
Definition: G4SIunits.hh:199
static G4Proton * Proton()
Definition: G4Proton.cc:93
const G4String & GetParticleType() const
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:321
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
G4bool UseOnlyOneMatrixForAllElements
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
G4double GetPDGMass() const
static G4AdjointProton * AdjointProton()
virtual void SampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
void ProposeEnergy(G4double finalEnergy)
static constexpr double GeV
Definition: G4SIunits.hh:217
void AddSecondary(G4Track *aSecondary)
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
G4double GetWeight() const
static const G4double Emin
G4double GetPDGSpin() const
static const G4double Emax
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
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)
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
G4double GetPDGCharge() const
static G4AdjointCSManager * GetAdjointCSManager()
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)
virtual G4double DiffCrossSectionPerAtomPrimToSecond(G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
G4AdjointhIonisationModel(G4ParticleDefinition *projectileDefinition)
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)