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

Last change on this file since 1036 was 966, checked in by garnier, 17 years ago

fichier ajoutes
File size: 9.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#include "G4AdjointComptonModel.hh"
27#include "G4AdjointCSManager.hh"
28
29
30#include "G4Integrator.hh"
31#include "G4TrackStatus.hh"
32#include "G4ParticleChange.hh"
33#include "G4AdjointElectron.hh"
34#include "G4AdjointGamma.hh"
35#include "G4Gamma.hh"
36
37
38////////////////////////////////////////////////////////////////////////////////
39//
40G4AdjointComptonModel::G4AdjointComptonModel():
41 G4VEmAdjointModel("AdjointCompton")
42
43{ SetApplyCutInRange(false);
44 SetUseMatrix(true);
45 SetUseMatrixPerElement(true);
46 SetIsIonisation(false);
47 SetUseOnlyOneMatrixForAllElements(true);
48 theAdjEquivOfDirectPrimPartDef =G4AdjointGamma::AdjointGamma();
49 theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
50 theDirectPrimaryPartDef=G4Gamma::Gamma();
51 second_part_of_same_type=false;
52
53}
54////////////////////////////////////////////////////////////////////////////////
55//
56G4AdjointComptonModel::~G4AdjointComptonModel()
57{;}
58////////////////////////////////////////////////////////////////////////////////
59//
60void G4AdjointComptonModel::SampleSecondaries(const G4Track& aTrack,
61 G4bool IsScatProjToProjCase,
62 G4ParticleChange* fParticleChange)
63{
64
65
66 //A recall of the compton scattering law is
67 //Egamma2=Egamma1/(1+(Egamma1/E0_electron)(1.-cos_th))
68 //Therefore Egamma2_max= Egamma2(cos_th=1) = Egamma1
69 //Therefore Egamma2_min= Egamma2(cos_th=-1) = Egamma1/(1+2.(Egamma1/E0_electron))
70
71
72 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
73 //DefineCurrentMaterial(aTrack->GetMaterialCutsCouple());
74 size_t ind= 0;
75
76
77
78
79 //Elastic inverse scattering //not correct in all the cases
80 //---------------------------------------------------------
81 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
82
83 //G4cout<<adjointPrimKinEnergy<<std::endl;
84 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
85 return;
86 }
87
88 //Sample secondary energy
89 //-----------------------
90 G4double gammaE1;
91
92 gammaE1 = SampleAdjSecEnergyFromCSMatrix(ind,
93 adjointPrimKinEnergy,
94 IsScatProjToProjCase);
95
96
97 //gammaE2
98 //-----------
99
100 G4double gammaE2 = adjointPrimKinEnergy;
101 if (!IsScatProjToProjCase) gammaE2 = gammaE1 - adjointPrimKinEnergy;
102
103
104
105
106
107
108 //Cos th
109 //-------
110// G4cout<<"Compton scattering "<<gammaE1<<'\t'<<gammaE2<<std::endl;
111 G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2);
112 if (!IsScatProjToProjCase) {
113 G4double p_elec=theAdjointPrimary->GetTotalMomentum();
114 cos_th = (gammaE1 - gammaE2*cos_th)/p_elec;
115 }
116 G4double sin_th = 0.;
117 if (std::abs(cos_th)>1){
118 //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<std::endl;
119 if (cos_th>0) {
120 cos_th=1.;
121 }
122 else cos_th=-1.;
123 sin_th=0.;
124 }
125 else sin_th = std::sqrt(1.-cos_th*cos_th);
126
127
128
129
130 //gamma0 momentum
131 //--------------------
132
133
134 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
135 G4double phi =G4UniformRand()*2.*3.1415926;
136 G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th);
137 gammaMomentum1.rotateUz(dir_parallel);
138// G4cout<<gamma0Energy<<'\t'<<gamma0Momentum<<std::endl;
139
140
141 //It is important to correct the weight of particles before adding the secondary
142 //------------------------------------------------------------------------------
143 CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), adjointPrimKinEnergy,gammaE1);
144
145 if (!IsScatProjToProjCase && CorrectWeightMode){ //kill the primary and add a secondary
146 fParticleChange->ProposeTrackStatus(fStopAndKill);
147 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
148 //G4cout<<"gamma0Momentum "<<gamma0Momentum<<std::endl;
149 }
150 else {
151 fParticleChange->ProposeEnergy(gammaE1);
152 fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
153 }
154
155
156}
157////////////////////////////////////////////////////////////////////////////////
158//
159//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
160G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToSecond(
161 G4double gamEnergy0,
162 G4double kinEnergyElec,
163 G4double Z,
164 G4double A)
165{
166 G4double gamEnergy1 = gamEnergy0 - kinEnergyElec;
167 G4double dSigmadEprod=0.;
168 if (gamEnergy1>0.) dSigmadEprod=DiffCrossSectionPerAtomPrimToScatPrim(gamEnergy0,gamEnergy1,Z,A);
169 return dSigmadEprod;
170}
171
172
173////////////////////////////////////////////////////////////////////////////////
174//
175G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim(
176 G4double gamEnergy0,
177 G4double gamEnergy1,
178 G4double Z,
179 G4double )
180{ //Based on Klein Nishina formula
181 //* In the forward case (see G4KleinNishinaModel) the cross section is parametrised while the secondaries are sampled from the
182 // Klein Nishida differential cross section
183 // The used diffrential cross section here is therefore the cross section multiplied by the normalidsed differential Klein Nishida cross section
184
185
186 //Klein Nishida Cross Section
187 //-----------------------------
188 G4double epsilon = gamEnergy0 / electron_mass_c2 ;
189 G4double one_plus_two_epsi =1.+2.*epsilon;
190 G4double gamEnergy1_max = gamEnergy0;
191 G4double gamEnergy1_min = gamEnergy0/one_plus_two_epsi;
192 if (gamEnergy1 >gamEnergy1_max || gamEnergy1<gamEnergy1_min) {
193 /*G4cout<<"the differential CS is null"<<std::endl;
194 G4cout<<gamEnergy0<<std::endl;
195 G4cout<<gamEnergy1<<std::endl;
196 G4cout<<gamEnergy1_min<<std::endl;*/
197 return 0.;
198 }
199
200
201 G4double epsi2 = epsilon *epsilon ;
202 G4double one_plus_two_epsi_2=one_plus_two_epsi*one_plus_two_epsi;
203
204
205 G4double CS=std::log(one_plus_two_epsi)*(1.- 2.*(1.+epsilon)/epsi2);
206 CS+=4./epsilon +0.5*(1.-1./one_plus_two_epsi_2);
207 CS/=epsilon;
208 //Note that the pi*re2*Z factor is neglected because it is suppresed when computing dCS_dE1/CS;
209 // in the differential cross section
210
211
212 //Klein Nishida Differential Cross Section
213 //-----------------------------------------
214 G4double epsilon1 = gamEnergy1 / electron_mass_c2 ;
215 G4double v= epsilon1/epsilon;
216 G4double term1 =1.+ 1./epsilon -1/epsilon1;
217 G4double dCS_dE1= 1./v +v + term1*term1 -1.;
218 dCS_dE1 *=1./epsilon/gamEnergy0;
219
220
221 //Normalised to the CS used in G4
222 //-------------------------------
223
224 G4double G4direct_CS = theDirectEMModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),
225 gamEnergy0,
226 Z, 0., 0.,0.);
227
228 dCS_dE1 *= G4direct_CS/CS;
229/* G4cout<<"the differential CS is not null"<<std::endl;
230 G4cout<<gamEnergy0<<std::endl;
231 G4cout<<gamEnergy1<<std::endl;*/
232
233 return dCS_dE1;
234
235
236}
237////////////////////////////////////////////////////////////////////////////////
238//
239G4double G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
240{ G4double inv_e_max = 1./PrimAdjEnergy - 2./electron_mass_c2;
241 G4double e_max = HighEnergyLimit;
242 if (inv_e_max > 0. ) e_max=std::min(1./inv_e_max,HighEnergyLimit);
243 return e_max;
244}
245////////////////////////////////////////////////////////////////////////////////
246//
247G4double G4AdjointComptonModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
248{ G4double half_e=PrimAdjEnergy/2.;
249 G4double term=std::sqrt(half_e*(electron_mass_c2+half_e));
250 G4double emin=half_e+term;
251 return emin;
252}
Note: See TracBrowser for help on using the repository browser.