Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 69844 2013-05-16 09:19:33Z 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 //
44 G4AdjointComptonModel::G4AdjointComptonModel():
45  G4VEmAdjointModel("AdjointCompton")
46 
47 { SetApplyCutInRange(false);
48  SetUseMatrix(false);
49  SetUseMatrixPerElement(true);
50  SetUseOnlyOneMatrixForAllElements(true);
51  theAdjEquivOfDirectPrimPartDef =G4AdjointGamma::AdjointGamma();
52  theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
53  theDirectPrimaryPartDef=G4Gamma::Gamma();
54  second_part_of_same_type=false;
55  theDirectEMModel=new G4KleinNishinaCompton(G4Gamma::Gamma(),"ComptonDirectModel");
56  G4direct_CS = 0.;
57 }
59 //
60 G4AdjointComptonModel::~G4AdjointComptonModel()
61 {;}
63 //
64 void G4AdjointComptonModel::SampleSecondaries(const G4Track& aTrack,
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 //
156 void G4AdjointComptonModel::RapidSampleSecondaries(const G4Track& aTrack,
157  G4bool IsScatProjToProjCase,
158  G4ParticleChange* fParticleChange)
159 {
160 
161  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
162  DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
163 
164 
165  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
166 
167 
168  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
169  return;
170  }
171 
172 
173 
174  G4double diffCSUsed=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);
191  G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
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
205  G4double 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  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
278 G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToSecond(
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 //
293 G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim(
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 while
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 
344  G4direct_CS = theDirectEMModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),
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 //
360 G4double G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
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 //
368 G4double G4AdjointComptonModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
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 //
376 G4double G4AdjointComptonModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
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= 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 = 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 //
412 G4double G4AdjointComptonModel::GetAdjointCrossSection(const G4MaterialCutsCouple* aCouple,
413  G4double primEnergy,
414  G4bool IsScatProjToProjCase)
415 { return AdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
416  //return G4VEmAdjointModel::GetAdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
417 }
Generated on Sat May 25 2013 14:33:29 for Geant4 by   doxygen 1.8.4