Geant4  10.03
G4AdjointComptonModel.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: G4AdjointComptonModel.cc 100666 2016-10-31 10:27:00Z gcosmo $
27 //
28 #include "G4AdjointComptonModel.hh"
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 "G4AdjointGamma.hh"
38 #include "G4Gamma.hh"
39 #include "G4KleinNishinaCompton.hh"
40 
41 
43 //
45  G4VEmAdjointModel("AdjointCompton")
46 
47 { SetApplyCutInRange(false);
48  SetUseMatrix(false);
55  theDirectEMModel=new G4KleinNishinaCompton(G4Gamma::Gamma(),"ComptonDirectModel");
56  G4direct_CS = 0.;
57 }
59 //
61 {;}
63 //
65  G4bool IsScatProjToProjCase,
66  G4ParticleChange* fParticleChange)
67 {
68  if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
69 
70  //A recall of the compton scattering law is
71  //Egamma2=Egamma1/(1+(Egamma1/E0_electron)(1.-cos_th))
72  //Therefore Egamma2_max= Egamma2(cos_th=1) = Egamma1
73  //Therefore Egamma2_min= Egamma2(cos_th=-1) = Egamma1/(1+2.(Egamma1/E0_electron))
74 
75 
76  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
77  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
78  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
79  return;
80  }
81 
82 
83  //Sample secondary energy
84  //-----------------------
85  G4double gammaE1;
86  gammaE1 = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy,
87  IsScatProjToProjCase);
88 
89 
90  //gammaE2
91  //-----------
92 
93  G4double gammaE2 = adjointPrimKinEnergy;
94  if (!IsScatProjToProjCase) gammaE2 = gammaE1 - adjointPrimKinEnergy;
95 
96 
97 
98 
99 
100 
101  //Cos th
102  //-------
103 // G4cout<<"Compton scattering "<<gammaE1<<'\t'<<gammaE2<<G4endl;
104  G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2);
105  if (!IsScatProjToProjCase) {
106  G4double p_elec=theAdjointPrimary->GetTotalMomentum();
107  cos_th = (gammaE1 - gammaE2*cos_th)/p_elec;
108  }
109  G4double sin_th = 0.;
110  if (std::abs(cos_th)>1){
111  //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl;
112  if (cos_th>0) {
113  cos_th=1.;
114  }
115  else cos_th=-1.;
116  sin_th=0.;
117  }
118  else sin_th = std::sqrt(1.-cos_th*cos_th);
119 
120 
121 
122 
123  //gamma0 momentum
124  //--------------------
125 
126 
127  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
128  G4double phi =G4UniformRand()*2.*3.1415926;
129  G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th);
130  gammaMomentum1.rotateUz(dir_parallel);
131 // G4cout<<gamma0Energy<<'\t'<<gamma0Momentum<<G4endl;
132 
133 
134  //It is important to correct the weight of particles before adding the secondary
135  //------------------------------------------------------------------------------
136  CorrectPostStepWeight(fParticleChange,
137  aTrack.GetWeight(),
138  adjointPrimKinEnergy,
139  gammaE1,
140  IsScatProjToProjCase);
141 
142  if (!IsScatProjToProjCase){ //kill the primary and add a secondary
143  fParticleChange->ProposeTrackStatus(fStopAndKill);
144  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
145  //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl;
146  }
147  else {
148  fParticleChange->ProposeEnergy(gammaE1);
149  fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
150  }
151 
152 
153 }
155 //
157  G4bool IsScatProjToProjCase,
158  G4ParticleChange* fParticleChange)
159 {
160 
161  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
163 
164 
165  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
166 
167 
168  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
169  return;
170  }
171 
172 
173 
174  G4double diffCSUsed=0.1*currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
175  G4double gammaE1=0.;
176  G4double gammaE2=0.;
177  if (!IsScatProjToProjCase){
178 
179  G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
180  G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);;
181  if (Emin>=Emax) return;
182  G4double f1=(Emin-adjointPrimKinEnergy)/Emin;
183  G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1;
184  gammaE1=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand()));;
185  gammaE2=gammaE1-adjointPrimKinEnergy;
186  diffCSUsed= diffCSUsed*(1.+2.*std::log(1.+electron_mass_c2/adjointPrimKinEnergy))*adjointPrimKinEnergy/gammaE1/gammaE2;
187 
188 
189  }
190  else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
192  if (Emin>=Emax) return;
193  gammaE2 =adjointPrimKinEnergy;
194  gammaE1=Emin*std::pow(Emax/Emin,G4UniformRand());
195  diffCSUsed= diffCSUsed/gammaE1;
196  }
197 
198 
199 
200 
201  //Weight correction
202  //-----------------------
203  //First w_corr is set to the ratio between adjoint total CS and fwd total CS
207  }
208  //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the
209  //one consistent with the direct model
210 
211 
212  G4double diffCS = DiffCrossSectionPerAtomPrimToScatPrim(gammaE1, gammaE2,1,0.);
213  if (diffCS >0) diffCS /=G4direct_CS; // here we have the normalised diffCS
214  //An we remultiply by the lambda of the forward process
215  diffCS*=theDirectEMProcess->GetLambda(gammaE1,currentCouple);
216  //diffCS*=theDirectEMModel->CrossSectionPerVolume(currentMaterial,G4Gamma::Gamma(),gammaE1,0.,2.*gammaE1);
217  //G4cout<<"diffCS/diffCSUsed "<<diffCS/diffCSUsed<<'\t'<<gammaE1<<'\t'<<gammaE2<<G4endl;
218 
219  w_corr*=diffCS/diffCSUsed;
220 
221  G4double new_weight = aTrack.GetWeight()*w_corr;
222  fParticleChange->SetParentWeightByProcess(false);
223  fParticleChange->SetSecondaryWeightByProcess(false);
224  fParticleChange->ProposeParentWeight(new_weight);
225 
226 
227 
228  //Cos th
229  //-------
230 
231  G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2);
232  if (!IsScatProjToProjCase) {
233  G4double p_elec=theAdjointPrimary->GetTotalMomentum();
234  cos_th = (gammaE1 - gammaE2*cos_th)/p_elec;
235  }
236  G4double sin_th = 0.;
237  if (std::abs(cos_th)>1){
238  //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl;
239  if (cos_th>0) {
240  cos_th=1.;
241  }
242  else cos_th=-1.;
243  sin_th=0.;
244  }
245  else sin_th = std::sqrt(1.-cos_th*cos_th);
246 
247 
248 
249 
250  //gamma0 momentum
251  //--------------------
252 
253 
254  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
255  G4double phi =G4UniformRand()*2.*3.1415926;
256  G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th);
257  gammaMomentum1.rotateUz(dir_parallel);
258 
259 
260 
261 
262  if (!IsScatProjToProjCase){ //kill the primary and add a secondary
263  fParticleChange->ProposeTrackStatus(fStopAndKill);
264  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
265  //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl;
266  }
267  else {
268  fParticleChange->ProposeEnergy(gammaE1);
269  fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
270  }
271 
272 
273 
274 }
275 
276 
278 //
279 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
281  G4double gamEnergy0,
282  G4double kinEnergyElec,
283  G4double Z,
284  G4double A)
285 {
286  G4double gamEnergy1 = gamEnergy0 - kinEnergyElec;
287  G4double dSigmadEprod=0.;
288  if (gamEnergy1>0.) dSigmadEprod=DiffCrossSectionPerAtomPrimToScatPrim(gamEnergy0,gamEnergy1,Z,A);
289  return dSigmadEprod;
290 }
291 
292 
294 //
296  G4double gamEnergy0,
297  G4double gamEnergy1,
298  G4double Z,
299  G4double )
300 { //Based on Klein Nishina formula
301  // In the forward case (see G4KleinNishinaCompton) the cross section is parametrised
302  // but the secondaries are sampled from the
303  // Klein Nishida differential cross section
304  // The used diffrential cross section here is therefore the cross section multiplied by the normalised
305  //differential Klein Nishida cross section
306 
307 
308  //Klein Nishida Cross Section
309  //-----------------------------
310  G4double epsilon = gamEnergy0 / electron_mass_c2 ;
311  G4double one_plus_two_epsi =1.+2.*epsilon;
312  G4double gamEnergy1_max = gamEnergy0;
313  G4double gamEnergy1_min = gamEnergy0/one_plus_two_epsi;
314  if (gamEnergy1 >gamEnergy1_max || gamEnergy1<gamEnergy1_min) {
315  /*G4cout<<"the differential CS is null"<<G4endl;
316  G4cout<<gamEnergy0<<G4endl;
317  G4cout<<gamEnergy1<<G4endl;
318  G4cout<<gamEnergy1_min<<G4endl;*/
319  return 0.;
320  }
321 
322 
323  G4double epsi2 = epsilon *epsilon ;
324  G4double one_plus_two_epsi_2=one_plus_two_epsi*one_plus_two_epsi;
325 
326 
327  G4double CS=std::log(one_plus_two_epsi)*(1.- 2.*(1.+epsilon)/epsi2);
328  CS+=4./epsilon +0.5*(1.-1./one_plus_two_epsi_2);
329  CS/=epsilon;
330  //Note that the pi*re2*Z factor is neglected because it is suppresed when computing dCS_dE1/CS;
331  // in the differential cross section
332 
333 
334  //Klein Nishida Differential Cross Section
335  //-----------------------------------------
336  G4double epsilon1 = gamEnergy1 / electron_mass_c2 ;
337  G4double v= epsilon1/epsilon;
338  G4double term1 =1.+ 1./epsilon -1/epsilon1;
339  G4double dCS_dE1= 1./v +v + term1*term1 -1.;
340  dCS_dE1 *=1./epsilon/gamEnergy0;
341 
342 
343  //Normalised to the CS used in G4
344  //-------------------------------
345 
347  gamEnergy0,
348  Z, 0., 0.,0.);
349 
350  dCS_dE1 *= G4direct_CS/CS;
351 /* G4cout<<"the differential CS is not null"<<G4endl;
352  G4cout<<gamEnergy0<<G4endl;
353  G4cout<<gamEnergy1<<G4endl;*/
354 
355  return dCS_dE1;
356 
357 
358 }
359 
361 //
363 { G4double inv_e_max = 1./PrimAdjEnergy - 2./electron_mass_c2;
364  G4double e_max = HighEnergyLimit;
365  if (inv_e_max > 0. ) e_max=std::min(1./inv_e_max,HighEnergyLimit);
366  return e_max;
367 }
369 //
371 { G4double half_e=PrimAdjEnergy/2.;
372  G4double term=std::sqrt(half_e*(electron_mass_c2+half_e));
373  G4double emin=half_e+term;
374  return emin;
375 }
377 //
379  G4double primEnergy,
380  G4bool IsScatProjToProjCase)
381 {
382  if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
383  DefineCurrentMaterial(aCouple);
384 
385 
386  float Cross=0.;
387  float Emax_proj =0.;
388  float Emin_proj =0.;
389  if (!IsScatProjToProjCase ){
390  Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
391  Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
392  if (Emax_proj>Emin_proj ){
393  Cross= 0.1*std::log((Emax_proj-float (primEnergy))*Emin_proj/Emax_proj/(Emin_proj-primEnergy))
394  *(1.+2.*std::log(float(1.+electron_mass_c2/primEnergy)));
395  }
396  }
397  else {
398  Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
399  Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,0.);
400  if (Emax_proj>Emin_proj) {
401  Cross = 0.1*std::log(Emax_proj/Emin_proj);
402  //+0.5*primEnergy*primEnergy(1./(Emin_proj*Emin_proj) - 1./(Emax_proj*Emax_proj)); neglected at the moment
403  }
404 
405 
406  }
407 
408  Cross*=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
409  lastCS=Cross;
410  return double(Cross);
411 }
413 //
415  G4double primEnergy,
416  G4bool IsScatProjToProjCase)
417 { return AdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
418  //return G4VEmAdjointModel::GetAdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
419 }
static G4AdjointGamma * AdjointGamma()
G4double GetLambda(G4double &kinEnergy, const G4MaterialCutsCouple *couple)
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
virtual void SampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetPostStepWeightCorrection()
const G4DynamicParticle * GetDynamicParticle() const
G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex, G4double prim_energy, G4bool IsScatProjToProjCase)
G4ParticleDefinition * theDirectPrimaryPartDef
void ProposeParentWeight(G4double finalWeight)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
void RapidSampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
static G4AdjointElectron * AdjointElectron()
virtual G4double GetAdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
G4Material * currentMaterial
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
G4double GetTotalMomentum() const
void SetUseMatrixPerElement(G4bool aBool)
void SetSecondaryWeightByProcess(G4bool)
void SetParentWeightByProcess(G4bool)
G4VEmModel * theDirectEMModel
#define G4UniformRand()
Definition: Randomize.hh:97
double A(double temperature)
G4double GetElectronDensity() const
Definition: G4Material.hh:217
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
void SetUseMatrix(G4bool aBool)
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:321
G4MaterialCutsCouple * currentCouple
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
G4bool correct_weight_for_post_step_in_model
void SetUseOnlyOneMatrixForAllElements(G4bool aBool)
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
G4double additional_weight_correction_factor_for_post_step_outside_model
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
void ProposeEnergy(G4double finalEnergy)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void AddSecondary(G4Track *aSecondary)
G4double GetWeight() const
static const G4double Emin
static const G4double Emax
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
G4double currentTcutForDirectSecond
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
void SetApplyCutInRange(G4bool aBool)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
virtual G4double DiffCrossSectionPerAtomPrimToSecond(G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
static G4AdjointCSManager * GetAdjointCSManager()
double epsilon(double density, double temperature)
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)
virtual G4double DiffCrossSectionPerAtomPrimToScatPrim(G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)