source: trunk/source/processes/electromagnetic/adjoint/src/G4AdjointComptonModel.cc @ 1201

Last change on this file since 1201 was 1196, checked in by garnier, 15 years ago

update CVS release candidate geant4.9.3.01
File size: 14.0 KB
Line 
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,v 1.5 2009/11/20 10:31:20 ldesorgh Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28//
29#include "G4AdjointComptonModel.hh"
30#include "G4AdjointCSManager.hh"
31
32
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
42////////////////////////////////////////////////////////////////////////////////
43//
44G4AdjointComptonModel::G4AdjointComptonModel():
45 G4VEmAdjointModel("AdjointCompton")
46
47{ SetApplyCutInRange(false);
48  SetUseMatrix(true);
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 
57}
58////////////////////////////////////////////////////////////////////////////////
59//
60G4AdjointComptonModel::~G4AdjointComptonModel()
61{;}
62////////////////////////////////////////////////////////////////////////////////
63//
64void 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} 
154////////////////////////////////////////////////////////////////////////////////
155//
156void 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 G4double diffCSUsed=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*pow(f2,G4UniformRand()));;
184        gammaE2=gammaE1-adjointPrimKinEnergy;
185        diffCSUsed= diffCSUsed*(1.+2.*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*pow(Emax/Emin,G4UniformRand());
194        diffCSUsed= diffCSUsed/gammaE1;
195       
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
204 G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
205
206 //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the
207 //one consistent with the direct model
208 
209 
210 G4double diffCS = DiffCrossSectionPerAtomPrimToScatPrim(gammaE1, gammaE2,1,0.);
211 if (diffCS >0)  diffCS /=G4direct_CS;  // here we have the normalised diffCS
212 diffCS*=theDirectEMProcess->GetLambda(gammaE1,currentCouple);
213 //diffCS*=theDirectEMModel->CrossSectionPerVolume(currentMaterial,G4Gamma::Gamma(),gammaE1,0.,2.*gammaE1);
214 //G4cout<<"diffCS/diffCSUsed "<<diffCS/diffCSUsed<<'\t'<<gammaE1<<'\t'<<gammaE2<<G4endl;                                 
215 
216 w_corr*=diffCS/diffCSUsed;
217           
218 G4double new_weight = aTrack.GetWeight()*w_corr;
219 fParticleChange->SetParentWeightByProcess(false);
220 fParticleChange->SetSecondaryWeightByProcess(false);
221 fParticleChange->ProposeParentWeight(new_weight); 
222 
223 
224 
225 //Cos th
226 //-------
227
228 G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2);
229 if (!IsScatProjToProjCase) {
230        G4double p_elec=theAdjointPrimary->GetTotalMomentum();
231        cos_th = (gammaE1 - gammaE2*cos_th)/p_elec;
232 }
233 G4double sin_th = 0.;
234 if (std::abs(cos_th)>1){
235        //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl;
236        if (cos_th>0) {
237                cos_th=1.;
238        }
239        else    cos_th=-1.;
240        sin_th=0.;
241 }
242 else  sin_th = std::sqrt(1.-cos_th*cos_th);
243
244 
245 
246 
247 //gamma0 momentum
248 //--------------------
249
250 
251 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
252 G4double phi =G4UniformRand()*2.*3.1415926;
253 G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th);
254 gammaMomentum1.rotateUz(dir_parallel);
255 
256 
257
258 
259 if (!IsScatProjToProjCase){ //kill the primary and add a secondary
260        fParticleChange->ProposeTrackStatus(fStopAndKill);
261        fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
262        //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl;
263 }
264 else {
265        fParticleChange->ProposeEnergy(gammaE1);
266        fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
267 }
268 
269 
270 
271} 
272
273                       
274////////////////////////////////////////////////////////////////////////////////
275//
276//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine 
277G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToSecond(
278                                      G4double gamEnergy0, 
279                                      G4double kinEnergyElec,
280                                      G4double Z, 
281                                      G4double A)
282{ 
283  G4double gamEnergy1 =  gamEnergy0 - kinEnergyElec;
284  G4double dSigmadEprod=0.;
285  if (gamEnergy1>0.) dSigmadEprod=DiffCrossSectionPerAtomPrimToScatPrim(gamEnergy0,gamEnergy1,Z,A);
286  return dSigmadEprod;   
287}
288
289
290////////////////////////////////////////////////////////////////////////////////
291//
292G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim(
293                                      G4double gamEnergy0, 
294                                      G4double gamEnergy1,
295                                      G4double Z, 
296                                      G4double )
297{ //Based on Klein Nishina formula
298 // In the forward case (see G4KleinNishinaModel)  the  cross section is parametrised while
299 // the secondaries are sampled from the
300 // Klein Nishida differential cross section
301 // The used diffrential cross section here is therefore the cross section multiplied by the normalised
302 //differential Klein Nishida cross section
303 
304 
305 //Klein Nishida Cross Section
306 //-----------------------------
307 G4double epsilon = gamEnergy0 / electron_mass_c2 ;
308 G4double one_plus_two_epsi =1.+2.*epsilon;
309 G4double gamEnergy1_max = gamEnergy0;
310 G4double gamEnergy1_min = gamEnergy0/one_plus_two_epsi;
311 if (gamEnergy1 >gamEnergy1_max ||  gamEnergy1<gamEnergy1_min) {
312        /*G4cout<<"the differential CS is null"<<G4endl;
313        G4cout<<gamEnergy0<<G4endl;
314        G4cout<<gamEnergy1<<G4endl;
315        G4cout<<gamEnergy1_min<<G4endl;*/
316        return 0.;
317 }
318       
319 
320 G4double epsi2 = epsilon *epsilon ;
321 G4double one_plus_two_epsi_2=one_plus_two_epsi*one_plus_two_epsi;
322 
323 
324 G4double CS=std::log(one_plus_two_epsi)*(1.- 2.*(1.+epsilon)/epsi2);
325 CS+=4./epsilon +0.5*(1.-1./one_plus_two_epsi_2);
326 CS/=epsilon;
327 //Note that the pi*re2*Z factor is neglected because it is suppresed when computing dCS_dE1/CS;
328 // in the differential cross section
329 
330 
331 //Klein Nishida Differential Cross Section
332 //-----------------------------------------
333 G4double epsilon1 = gamEnergy1 / electron_mass_c2 ;
334 G4double v= epsilon1/epsilon;
335 G4double term1 =1.+ 1./epsilon -1/epsilon1;
336 G4double dCS_dE1= 1./v +v + term1*term1 -1.;
337 dCS_dE1 *=1./epsilon/gamEnergy0;
338 
339 
340 //Normalised to the CS used in G4
341 //-------------------------------
342 
343 G4direct_CS = theDirectEMModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),
344                                             gamEnergy0,
345                                             Z, 0., 0.,0.);
346 
347 dCS_dE1 *= G4direct_CS/CS;
348/* G4cout<<"the differential CS is not null"<<G4endl;
349 G4cout<<gamEnergy0<<G4endl;
350 G4cout<<gamEnergy1<<G4endl;*/
351 
352 return dCS_dE1;
353
354
355}
356                                     
357////////////////////////////////////////////////////////////////////////////////
358//
359G4double G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
360{ G4double inv_e_max =  1./PrimAdjEnergy - 2./electron_mass_c2;
361  G4double e_max = HighEnergyLimit;
362  if (inv_e_max > 0. ) e_max=std::min(1./inv_e_max,HighEnergyLimit);
363  return  e_max; 
364}
365////////////////////////////////////////////////////////////////////////////////
366//
367G4double G4AdjointComptonModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
368{ G4double half_e=PrimAdjEnergy/2.;
369  G4double term=std::sqrt(half_e*(electron_mass_c2+half_e));
370  G4double emin=half_e+term;
371  return  emin; 
372}
373////////////////////////////////////////////////////////////////////////////////
374//
375G4double G4AdjointComptonModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
376                                             G4double primEnergy,
377                                             G4bool IsScatProjToProjCase)
378{ 
379  if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
380  DefineCurrentMaterial(aCouple);
381 
382 
383  G4double Cross=0.;
384  G4double Emax_proj =0.;
385  G4double Emin_proj =0.;
386  if (!IsScatProjToProjCase ){
387        Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
388        Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
389        if (Emax_proj>Emin_proj ){
390                 Cross= log((Emax_proj-primEnergy)*Emin_proj/Emax_proj/(Emin_proj-primEnergy))
391                                                                *(1.+2.*log(1.+electron_mass_c2/primEnergy));
392        }       
393  }
394  else {
395        Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
396        Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,0.);
397        if (Emax_proj>Emin_proj) {
398                Cross = log(Emax_proj/Emin_proj);
399                //+0.5*primEnergy*primEnergy(1./(Emin_proj*Emin_proj) - 1./(Emax_proj*Emax_proj)); neglected at the moment
400        }
401       
402       
403  }
404 
405  Cross*=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
406  lastCS=Cross;
407  return Cross; 
408}
Note: See TracBrowser for help on using the repository browser.