Geant4  10.02.p01
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 91870 2015-08-07 15:21:40Z 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 
202  //Weight correction
203  //-----------------------
204  //First w_corr is set to the ratio between adjoint total CS and fwd total CS
206 
207  //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the
208  //one consistent with the direct model
209 
210 
211  G4double diffCS = DiffCrossSectionPerAtomPrimToScatPrim(gammaE1, gammaE2,1,0.);
212  if (diffCS >0) diffCS /=G4direct_CS; // here we have the normalised diffCS
213  diffCS*=theDirectEMProcess->GetLambda(gammaE1,currentCouple);
214  //diffCS*=theDirectEMModel->CrossSectionPerVolume(currentMaterial,G4Gamma::Gamma(),gammaE1,0.,2.*gammaE1);
215  //G4cout<<"diffCS/diffCSUsed "<<diffCS/diffCSUsed<<'\t'<<gammaE1<<'\t'<<gammaE2<<G4endl;
216 
217  w_corr*=diffCS/diffCSUsed;
218 
219  G4double new_weight = aTrack.GetWeight()*w_corr;
220  fParticleChange->SetParentWeightByProcess(false);
221  fParticleChange->SetSecondaryWeightByProcess(false);
222  fParticleChange->ProposeParentWeight(new_weight);
223 
224 
225 
226  //Cos th
227  //-------
228 
229  G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2);
230  if (!IsScatProjToProjCase) {
231  G4double p_elec=theAdjointPrimary->GetTotalMomentum();
232  cos_th = (gammaE1 - gammaE2*cos_th)/p_elec;
233  }
234  G4double sin_th = 0.;
235  if (std::abs(cos_th)>1){
236  //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl;
237  if (cos_th>0) {
238  cos_th=1.;
239  }
240  else cos_th=-1.;
241  sin_th=0.;
242  }
243  else sin_th = std::sqrt(1.-cos_th*cos_th);
244 
245 
246 
247 
248  //gamma0 momentum
249  //--------------------
250 
251 
252  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
253  G4double phi =G4UniformRand()*2.*3.1415926;
254  G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th);
255  gammaMomentum1.rotateUz(dir_parallel);
256 
257 
258 
259 
260  if (!IsScatProjToProjCase){ //kill the primary and add a secondary
261  fParticleChange->ProposeTrackStatus(fStopAndKill);
262  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
263  //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl;
264  }
265  else {
266  fParticleChange->ProposeEnergy(gammaE1);
267  fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
268  }
269 
270 
271 
272 }
273 
274 
276 //
277 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
279  G4double gamEnergy0,
280  G4double kinEnergyElec,
281  G4double Z,
282  G4double A)
283 {
284  G4double gamEnergy1 = gamEnergy0 - kinEnergyElec;
285  G4double dSigmadEprod=0.;
286  if (gamEnergy1>0.) dSigmadEprod=DiffCrossSectionPerAtomPrimToScatPrim(gamEnergy0,gamEnergy1,Z,A);
287  return dSigmadEprod;
288 }
289 
290 
292 //
294  G4double gamEnergy0,
295  G4double gamEnergy1,
296  G4double Z,
297  G4double )
298 { //Based on Klein Nishina formula
299  // In the forward case (see G4KleinNishinaModel) the cross section is parametrised
300  // the secondaries are sampled from the
301  // Klein Nishida differential cross section
302  // The used diffrential cross section here is therefore the cross section multiplied by the normalised
303  //differential Klein Nishida cross section
304 
305 
306  //Klein Nishida Cross Section
307  //-----------------------------
308  G4double epsilon = gamEnergy0 / electron_mass_c2 ;
309  G4double one_plus_two_epsi =1.+2.*epsilon;
310  G4double gamEnergy1_max = gamEnergy0;
311  G4double gamEnergy1_min = gamEnergy0/one_plus_two_epsi;
312  if (gamEnergy1 >gamEnergy1_max || gamEnergy1<gamEnergy1_min) {
313  /*G4cout<<"the differential CS is null"<<G4endl;
314  G4cout<<gamEnergy0<<G4endl;
315  G4cout<<gamEnergy1<<G4endl;
316  G4cout<<gamEnergy1_min<<G4endl;*/
317  return 0.;
318  }
319 
320 
321  G4double epsi2 = epsilon *epsilon ;
322  G4double one_plus_two_epsi_2=one_plus_two_epsi*one_plus_two_epsi;
323 
324 
325  G4double CS=std::log(one_plus_two_epsi)*(1.- 2.*(1.+epsilon)/epsi2);
326  CS+=4./epsilon +0.5*(1.-1./one_plus_two_epsi_2);
327  CS/=epsilon;
328  //Note that the pi*re2*Z factor is neglected because it is suppresed when computing dCS_dE1/CS;
329  // in the differential cross section
330 
331 
332  //Klein Nishida Differential Cross Section
333  //-----------------------------------------
334  G4double epsilon1 = gamEnergy1 / electron_mass_c2 ;
335  G4double v= epsilon1/epsilon;
336  G4double term1 =1.+ 1./epsilon -1/epsilon1;
337  G4double dCS_dE1= 1./v +v + term1*term1 -1.;
338  dCS_dE1 *=1./epsilon/gamEnergy0;
339 
340 
341  //Normalised to the CS used in G4
342  //-------------------------------
343 
345  gamEnergy0,
346  Z, 0., 0.,0.);
347 
348  dCS_dE1 *= G4direct_CS/CS;
349 /* G4cout<<"the differential CS is not null"<<G4endl;
350  G4cout<<gamEnergy0<<G4endl;
351  G4cout<<gamEnergy1<<G4endl;*/
352 
353  return dCS_dE1;
354 
355 
356 }
357 
359 //
361 { G4double inv_e_max = 1./PrimAdjEnergy - 2./electron_mass_c2;
362  G4double e_max = HighEnergyLimit;
363  if (inv_e_max > 0. ) e_max=std::min(1./inv_e_max,HighEnergyLimit);
364  return e_max;
365 }
367 //
369 { G4double half_e=PrimAdjEnergy/2.;
370  G4double term=std::sqrt(half_e*(electron_mass_c2+half_e));
371  G4double emin=half_e+term;
372  return emin;
373 }
375 //
377  G4double primEnergy,
378  G4bool IsScatProjToProjCase)
379 {
380  if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
381  DefineCurrentMaterial(aCouple);
382 
383 
384  float Cross=0.;
385  float Emax_proj =0.;
386  float Emin_proj =0.;
387  if (!IsScatProjToProjCase ){
388  Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
389  Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
390  if (Emax_proj>Emin_proj ){
391  Cross= 0.1*std::log((Emax_proj-float (primEnergy))*Emin_proj/Emax_proj/(Emin_proj-primEnergy))
392  *(1.+2.*std::log(float(1.+electron_mass_c2/primEnergy)));
393  }
394  }
395  else {
396  Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
397  Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,0.);
398  if (Emax_proj>Emin_proj) {
399  Cross = 0.1*std::log(Emax_proj/Emin_proj);
400  //+0.5*primEnergy*primEnergy(1./(Emin_proj*Emin_proj) - 1./(Emax_proj*Emax_proj)); neglected at the moment
401  }
402 
403 
404  }
405 
406  Cross*=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
407  lastCS=Cross;
408  return double(Cross);
409 }
411 //
413  G4double primEnergy,
414  G4bool IsScatProjToProjCase)
415 { return AdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
416  //return G4VEmAdjointModel::GetAdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
417 }
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)
static const G4double f2
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
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
G4MaterialCutsCouple * currentCouple
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
void SetUseOnlyOneMatrixForAllElements(G4bool aBool)
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
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.)