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