| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4AdjointComptonModel.cc 75591 2013-11-04 12:33:11Z gcosmo $
26 // 27 //
27 <<
28 #include "G4AdjointComptonModel.hh" 28 #include "G4AdjointComptonModel.hh"
29 <<
30 #include "G4AdjointCSManager.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"
31 #include "G4AdjointElectron.hh" 36 #include "G4AdjointElectron.hh"
32 #include "G4AdjointGamma.hh" 37 #include "G4AdjointGamma.hh"
33 #include "G4Gamma.hh" 38 #include "G4Gamma.hh"
34 #include "G4KleinNishinaCompton.hh" 39 #include "G4KleinNishinaCompton.hh"
35 #include "G4ParticleChange.hh" << 40
36 #include "G4PhysicalConstants.hh" <<
37 #include "G4SystemOfUnits.hh" <<
38 #include "G4TrackStatus.hh" <<
39 #include "G4VEmProcess.hh" <<
40 41
41 ////////////////////////////////////////////// 42 ////////////////////////////////////////////////////////////////////////////////
42 G4AdjointComptonModel::G4AdjointComptonModel() << 43 //
43 : G4VEmAdjointModel("AdjointCompton") << 44 G4AdjointComptonModel::G4AdjointComptonModel():
44 { << 45 G4VEmAdjointModel("AdjointCompton")
45 SetApplyCutInRange(false); << 46
>> 47 { SetApplyCutInRange(false);
46 SetUseMatrix(false); 48 SetUseMatrix(false);
47 SetUseMatrixPerElement(true); 49 SetUseMatrixPerElement(true);
48 SetUseOnlyOneMatrixForAllElements(true); 50 SetUseOnlyOneMatrixForAllElements(true);
49 fAdjEquivDirectPrimPart = G4AdjointGamma:: << 51 theAdjEquivOfDirectPrimPartDef =G4AdjointGamma::AdjointGamma();
50 fAdjEquivDirectSecondPart = G4AdjointElectro << 52 theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
51 fDirectPrimaryPart = G4Gamma::Gamma() << 53 theDirectPrimaryPartDef=G4Gamma::Gamma();
52 fSecondPartSameType = false; << 54 second_part_of_same_type=false;
53 fDirectModel = << 55 theDirectEMModel=new G4KleinNishinaCompton(G4Gamma::Gamma(),"ComptonDirectModel");
54 new G4KleinNishinaCompton(G4Gamma::Gamma() << 56 G4direct_CS = 0.;
55 } 57 }
56 <<
57 ////////////////////////////////////////////// 58 ////////////////////////////////////////////////////////////////////////////////
58 G4AdjointComptonModel::~G4AdjointComptonModel( << 59 //
59 << 60 G4AdjointComptonModel::~G4AdjointComptonModel()
>> 61 {;}
60 ////////////////////////////////////////////// 62 ////////////////////////////////////////////////////////////////////////////////
>> 63 //
61 void G4AdjointComptonModel::SampleSecondaries( 64 void G4AdjointComptonModel::SampleSecondaries(const G4Track& aTrack,
62 << 65 G4bool IsScatProjToProjCase,
63 << 66 G4ParticleChange* fParticleChange)
64 { << 67 {
65 if(!fUseMatrix) << 68 if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
66 return RapidSampleSecondaries(aTrack, isSc << 69
67 << 70 //A recall of the compton scattering law is
68 // A recall of the compton scattering law: << 71 //Egamma2=Egamma1/(1+(Egamma1/E0_electron)(1.-cos_th))
69 // Egamma2=Egamma1/(1+(Egamma1/E0_electron)( << 72 //Therefore Egamma2_max= Egamma2(cos_th=1) = Egamma1
70 // Therefore Egamma2_max= Egamma2(cos_th=1) << 73 //Therefore Egamma2_min= Egamma2(cos_th=-1) = Egamma1/(1+2.(Egamma1/E0_electron))
71 // and Egamma2_min= Egamma2(cos_th=-1) << 74
72 // Egamma1/(1+2.(Egamma1/E0_elec << 75
73 const G4DynamicParticle* theAdjointPrimary = << 76 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
74 G4double adjointPrimKinEnergy = theAdjointPr 77 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
75 if(adjointPrimKinEnergy > GetHighEnergyLimit << 78 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
76 { << 79 return;
77 return; <<
78 } <<
79 <<
80 // Sample secondary energy <<
81 G4double gammaE1; <<
82 gammaE1 = <<
83 SampleAdjSecEnergyFromCSMatrix(adjointPrim <<
84 <<
85 // gammaE2 <<
86 G4double gammaE2 = adjointPrimKinEnergy; <<
87 if(!isScatProjToProj) <<
88 gammaE2 = gammaE1 - adjointPrimKinEnergy; <<
89 <<
90 // Cos th <<
91 G4double cos_th = 1. + electron_mass_c2 * (1 <<
92 if(!isScatProjToProj) <<
93 { <<
94 cos_th = <<
95 (gammaE1 - gammaE2 * cos_th) / theAdjoin <<
96 } <<
97 G4double sin_th = 0.; <<
98 if(std::abs(cos_th) > 1.) <<
99 { <<
100 if(cos_th > 0.) <<
101 { <<
102 cos_th = 1.; <<
103 } <<
104 else <<
105 cos_th = -1.; <<
106 sin_th = 0.; <<
107 } <<
108 else <<
109 sin_th = std::sqrt(1. - cos_th * cos_th); <<
110 <<
111 // gamma0 momentum <<
112 G4ThreeVector dir_parallel = theAdjointPrima <<
113 G4double phi = G4UniformRand() <<
114 G4ThreeVector gammaMomentum1 = <<
115 gammaE1 * <<
116 G4ThreeVector(std::cos(phi) * sin_th, std: <<
117 gammaMomentum1.rotateUz(dir_parallel); <<
118 <<
119 // correct the weight of particles before ad <<
120 CorrectPostStepWeight(fParticleChange, aTrac <<
121 adjointPrimKinEnergy, <<
122 <<
123 if(!isScatProjToProj) <<
124 { // kill the primary and add a secondary <<
125 fParticleChange->ProposeTrackStatus(fStopA <<
126 fParticleChange->AddSecondary( <<
127 new G4DynamicParticle(fAdjEquivDirectPri <<
128 } <<
129 else <<
130 { <<
131 fParticleChange->ProposeEnergy(gammaE1); <<
132 fParticleChange->ProposeMomentumDirection( <<
133 } 80 }
134 } << 81
135 << 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 }
136 ////////////////////////////////////////////// 154 ////////////////////////////////////////////////////////////////////////////////
137 void G4AdjointComptonModel::RapidSampleSeconda << 155 //
138 const G4Track& aTrack, G4bool isScatProjToPr << 156 void G4AdjointComptonModel::RapidSampleSecondaries(const G4Track& aTrack,
139 G4ParticleChange* fParticleChange) << 157 G4bool IsScatProjToProjCase,
140 { << 158 G4ParticleChange* fParticleChange)
141 const G4DynamicParticle* theAdjointPrimary = << 159 {
142 DefineCurrentMaterial(aTrack.GetMaterialCuts << 160
143 << 161 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
144 G4double adjointPrimKinEnergy = theAdjointPr << 162 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
145 << 163
146 if(adjointPrimKinEnergy > GetHighEnergyLimit << 164
147 { << 165 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
148 return; << 166
149 } << 167
150 << 168 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
151 G4double diffCSUsed = << 169 return;
152 0.1 * fCurrentMaterial->GetElectronDensity << 170 }
153 G4double gammaE1 = 0.; << 171
154 G4double gammaE2 = 0.; << 172
155 if(!isScatProjToProj) << 173
156 { << 174 G4double diffCSUsed=0.1*currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
157 G4double Emax = GetSecondAdjEnergyMaxForPr << 175 G4double gammaE1=0.;
158 G4double Emin = GetSecondAdjEnergyMinForPr << 176 G4double gammaE2=0.;
159 if(Emin >= Emax) << 177 if (!IsScatProjToProjCase){
160 return; << 178
161 G4double f1 = (Emin - adjointPrimKinEnergy << 179 G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
162 G4double f2 = (Emax - adjointPrimKinEnergy << 180 G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);;
163 gammaE1 = adjointPrimKinEnergy / (1. - f1 << 181 if (Emin>=Emax) return;
164 gammaE2 = gammaE1 - adjointPrimKinEnergy; << 182 G4double f1=(Emin-adjointPrimKinEnergy)/Emin;
165 diffCSUsed = << 183 G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1;
166 diffCSUsed * << 184 gammaE1=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand()));;
167 (1. + 2. * std::log(1. + electron_mass_c << 185 gammaE2=gammaE1-adjointPrimKinEnergy;
168 adjointPrimKinEnergy / gammaE1 / gammaE2 << 186 diffCSUsed= diffCSUsed*(1.+2.*std::log(1.+electron_mass_c2/adjointPrimKinEnergy))*adjointPrimKinEnergy/gammaE1/gammaE2;
169 } << 187
170 else << 188
171 { << 189 }
172 G4double Emax = << 190 else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
173 GetSecondAdjEnergyMaxForScatProjToProj(a << 191 G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
174 G4double Emin = << 192 if (Emin>=Emax) return;
175 GetSecondAdjEnergyMinForScatProjToProj(a << 193 gammaE2 =adjointPrimKinEnergy;
176 if(Emin >= Emax) << 194 gammaE1=Emin*std::pow(Emax/Emin,G4UniformRand());
177 return; << 195 diffCSUsed= diffCSUsed/gammaE1;
178 gammaE2 = adjointPrimKinEnergy; << 196
179 gammaE1 = Emin * std::pow(Emax / Emin, << 197 }
180 diffCSUsed = diffCSUsed / gammaE1; << 198
181 } << 199
182 << 200
183 // Weight correction << 201
184 // First w_corr is set to the ratio between << 202 //Weight correction
185 G4double w_corr = fOutsideWeightFactor; << 203 //-----------------------
186 if(fInModelWeightCorr) << 204 //First w_corr is set to the ratio between adjoint total CS and fwd total CS
187 { << 205 G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
188 w_corr = << 206
189 G4AdjointCSManager::GetAdjointCSManager( << 207 //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the
190 } << 208 //one consistent with the direct model
191 // Then another correction is needed because << 209
192 // been used rather than the one consistent << 210
193 << 211 G4double diffCS = DiffCrossSectionPerAtomPrimToScatPrim(gammaE1, gammaE2,1,0.);
194 G4double diffCS = << 212 if (diffCS >0) diffCS /=G4direct_CS; // here we have the normalised diffCS
195 DiffCrossSectionPerAtomPrimToScatPrim(gamm << 213 diffCS*=theDirectEMProcess->GetLambda(gammaE1,currentCouple);
196 if(diffCS > 0.) << 214 //diffCS*=theDirectEMModel->CrossSectionPerVolume(currentMaterial,G4Gamma::Gamma(),gammaE1,0.,2.*gammaE1);
197 diffCS /= fDirectCS; // here we have the << 215 //G4cout<<"diffCS/diffCSUsed "<<diffCS/diffCSUsed<<'\t'<<gammaE1<<'\t'<<gammaE2<<G4endl;
198 // And we remultiply by the lambda of the fo << 216
199 diffCS *= fDirectProcess->GetCrossSection(ga << 217 w_corr*=diffCS/diffCSUsed;
200 << 218
201 w_corr *= diffCS / diffCSUsed; << 219 G4double new_weight = aTrack.GetWeight()*w_corr;
202 << 220 fParticleChange->SetParentWeightByProcess(false);
203 G4double new_weight = aTrack.GetWeight() * w << 221 fParticleChange->SetSecondaryWeightByProcess(false);
204 fParticleChange->SetParentWeightByProcess(fa << 222 fParticleChange->ProposeParentWeight(new_weight);
205 fParticleChange->SetSecondaryWeightByProcess << 223
206 fParticleChange->ProposeParentWeight(new_wei << 224
207 << 225
208 G4double cos_th = 1. + electron_mass_c2 * (1 << 226 //Cos th
209 if(!isScatProjToProj) << 227 //-------
210 { << 228
211 G4double p_elec = theAdjointPrimary->GetTo << 229 G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2);
212 cos_th = (gammaE1 - gammaE2 * cos << 230 if (!IsScatProjToProjCase) {
213 } << 231 G4double p_elec=theAdjointPrimary->GetTotalMomentum();
214 G4double sin_th = 0.; << 232 cos_th = (gammaE1 - gammaE2*cos_th)/p_elec;
215 if(std::abs(cos_th) > 1.) << 233 }
216 { << 234 G4double sin_th = 0.;
217 if(cos_th > 0.) << 235 if (std::abs(cos_th)>1){
218 { << 236 //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl;
219 cos_th = 1.; << 237 if (cos_th>0) {
220 } << 238 cos_th=1.;
221 else << 239 }
222 cos_th = -1.; << 240 else cos_th=-1.;
223 } << 241 sin_th=0.;
224 else << 242 }
225 sin_th = std::sqrt(1. - cos_th * cos_th); << 243 else sin_th = std::sqrt(1.-cos_th*cos_th);
226 << 244
227 // gamma0 momentum << 245
228 G4ThreeVector dir_parallel = theAdjointPrima << 246
229 G4double phi = G4UniformRand() << 247
230 G4ThreeVector gammaMomentum1 = << 248 //gamma0 momentum
231 gammaE1 * << 249 //--------------------
232 G4ThreeVector(std::cos(phi) * sin_th, std: << 250
233 gammaMomentum1.rotateUz(dir_parallel); << 251
234 << 252 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
235 if(!isScatProjToProj) << 253 G4double phi =G4UniformRand()*2.*3.1415926;
236 { // kill the primary and add a secondary << 254 G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th);
237 fParticleChange->ProposeTrackStatus(fStopA << 255 gammaMomentum1.rotateUz(dir_parallel);
238 fParticleChange->AddSecondary( << 256
239 new G4DynamicParticle(fAdjEquivDirectPri << 257
240 } << 258
241 else << 259
242 { << 260 if (!IsScatProjToProjCase){ //kill the primary and add a secondary
243 fParticleChange->ProposeEnergy(gammaE1); << 261 fParticleChange->ProposeTrackStatus(fStopAndKill);
244 fParticleChange->ProposeMomentumDirection( << 262 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
245 } << 263 //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl;
246 } << 264 }
>> 265 else {
>> 266 fParticleChange->ProposeEnergy(gammaE1);
>> 267 fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
>> 268 }
>> 269
>> 270
>> 271
>> 272 }
247 273
>> 274
248 ////////////////////////////////////////////// 275 ////////////////////////////////////////////////////////////////////////////////
249 // The implementation here is correct for ener << 276 //
250 // photoelectric and compton scattering the me << 277 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
251 G4double G4AdjointComptonModel::DiffCrossSecti 278 G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToSecond(
252 G4double gamEnergy0, G4double kinEnergyElec, << 279 G4double gamEnergy0,
253 { << 280 G4double kinEnergyElec,
254 G4double gamEnergy1 = gamEnergy0 - kinEner << 281 G4double Z,
255 G4double dSigmadEprod = 0.; << 282 G4double A)
256 if(gamEnergy1 > 0.) << 283 {
257 dSigmadEprod = << 284 G4double gamEnergy1 = gamEnergy0 - kinEnergyElec;
258 DiffCrossSectionPerAtomPrimToScatPrim(ga << 285 G4double dSigmadEprod=0.;
259 return dSigmadEprod; << 286 if (gamEnergy1>0.) dSigmadEprod=DiffCrossSectionPerAtomPrimToScatPrim(gamEnergy0,gamEnergy1,Z,A);
>> 287 return dSigmadEprod;
260 } 288 }
261 289
>> 290
262 ////////////////////////////////////////////// 291 ////////////////////////////////////////////////////////////////////////////////
>> 292 //
263 G4double G4AdjointComptonModel::DiffCrossSecti 293 G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim(
264 G4double gamEnergy0, G4double gamEnergy1, G4 << 294 G4double gamEnergy0,
265 { << 295 G4double gamEnergy1,
266 // Based on Klein Nishina formula << 296 G4double Z,
267 // In the forward case (see G4KleinNishinaCo << 297 G4double )
268 // parametrised but the secondaries are samp << 298 { //Based on Klein Nishina formula
269 // differential cross section. The different << 299 // In the forward case (see G4KleinNishinaModel) the cross section is parametrised while
270 // is therefore the cross section multiplied << 300 // the secondaries are sampled from the
271 // differential Klein Nishina cross section << 301 // Klein Nishida differential cross section
272 << 302 // The used diffrential cross section here is therefore the cross section multiplied by the normalised
273 // Klein Nishina Cross Section << 303 //differential Klein Nishida cross section
274 G4double epsilon = gamEnergy0 / el << 304
275 G4double one_plus_two_epsi = 1. + 2. * epsil << 305
276 if(gamEnergy1 > gamEnergy0 || gamEnergy1 < g << 306 //Klein Nishida Cross Section
277 { << 307 //-----------------------------
278 return 0.; << 308 G4double epsilon = gamEnergy0 / electron_mass_c2 ;
279 } << 309 G4double one_plus_two_epsi =1.+2.*epsilon;
280 << 310 G4double gamEnergy1_max = gamEnergy0;
281 G4double CS = std::log(one_plus_two_epsi) * << 311 G4double gamEnergy1_min = gamEnergy0/one_plus_two_epsi;
282 (1. - 2. * (1. + epsilon) / (e << 312 if (gamEnergy1 >gamEnergy1_max || gamEnergy1<gamEnergy1_min) {
283 CS += << 313 /*G4cout<<"the differential CS is null"<<G4endl;
284 4. / epsilon + 0.5 * (1. - 1. / (one_plus_ << 314 G4cout<<gamEnergy0<<G4endl;
285 CS /= epsilon; << 315 G4cout<<gamEnergy1<<G4endl;
286 // Note that the pi*re2*Z factor is neglecte << 316 G4cout<<gamEnergy1_min<<G4endl;*/
287 // computing dCS_dE1/CS in the differential << 317 return 0.;
288 << 318 }
289 // Klein Nishina Differential Cross Section << 319
290 G4double epsilon1 = gamEnergy1 / electron_ma << 320
291 G4double v = epsilon1 / epsilon; << 321 G4double epsi2 = epsilon *epsilon ;
292 G4double term1 = 1. + 1. / epsilon - 1. / << 322 G4double one_plus_two_epsi_2=one_plus_two_epsi*one_plus_two_epsi;
293 G4double dCS_dE1 = 1. / v + v + term1 * ter << 323
294 dCS_dE1 *= 1. / epsilon / gamEnergy0; << 324
295 << 325 G4double CS=std::log(one_plus_two_epsi)*(1.- 2.*(1.+epsilon)/epsi2);
296 // Normalised to the CS used in G4 << 326 CS+=4./epsilon +0.5*(1.-1./one_plus_two_epsi_2);
297 fDirectCS = fDirectModel->ComputeCrossSectio << 327 CS/=epsilon;
298 G4Gamma::Gamma(), gamEnergy0, Z, 0., 0., 0 << 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;
299 354
300 dCS_dE1 *= fDirectCS / CS; <<
301 355
302 return dCS_dE1; <<
303 } 356 }
304 << 357
305 ////////////////////////////////////////////// 358 ////////////////////////////////////////////////////////////////////////////////
306 G4double G4AdjointComptonModel::GetSecondAdjEn << 359 //
307 G4double primAdjEnergy) << 360 G4double G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
308 { << 361 { G4double inv_e_max = 1./PrimAdjEnergy - 2./electron_mass_c2;
309 G4double inv_e_max = 1. / primAdjEnergy - 2. << 362 G4double e_max = HighEnergyLimit;
310 G4double e_max = GetHighEnergyLimit(); << 363 if (inv_e_max > 0. ) e_max=std::min(1./inv_e_max,HighEnergyLimit);
311 if(inv_e_max > 0.) << 364 return e_max;
312 e_max = std::min(1. / inv_e_max, e_max); <<
313 return e_max; <<
314 } 365 }
315 <<
316 ////////////////////////////////////////////// 366 ////////////////////////////////////////////////////////////////////////////////
317 G4double G4AdjointComptonModel::GetSecondAdjEn << 367 //
318 G4double primAdjEnergy) << 368 G4double G4AdjointComptonModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
319 { << 369 { G4double half_e=PrimAdjEnergy/2.;
320 G4double half_e = primAdjEnergy / 2.; << 370 G4double term=std::sqrt(half_e*(electron_mass_c2+half_e));
321 return half_e + std::sqrt(half_e * (electron << 371 G4double emin=half_e+term;
>> 372 return emin;
322 } 373 }
323 <<
324 ////////////////////////////////////////////// 374 ////////////////////////////////////////////////////////////////////////////////
325 G4double G4AdjointComptonModel::AdjointCrossSe << 375 //
326 const G4MaterialCutsCouple* aCouple, G4doubl << 376 G4double G4AdjointComptonModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
327 G4bool isScatProjToProj) << 377 G4double primEnergy,
328 { << 378 G4bool IsScatProjToProjCase)
329 if(fUseMatrix) << 379 {
330 return G4VEmAdjointModel::AdjointCrossSect << 380 if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
331 <<
332 DefineCurrentMaterial(aCouple); 381 DefineCurrentMaterial(aCouple);
333 << 382
334 G4float Cross = 0.; << 383
335 G4float Emax_proj = 0.; << 384 float Cross=0.;
336 G4float Emin_proj = 0.; << 385 float Emax_proj =0.;
337 if(!isScatProjToProj) << 386 float Emin_proj =0.;
338 { << 387 if (!IsScatProjToProjCase ){
339 Emax_proj = GetSecondAdjEnergyMaxForProdTo << 388 Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
340 Emin_proj = GetSecondAdjEnergyMinForProdTo << 389 Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
341 if(Emax_proj > Emin_proj) << 390 if (Emax_proj>Emin_proj ){
342 { << 391 Cross= 0.1*std::log((Emax_proj-float (primEnergy))*Emin_proj/Emax_proj/(Emin_proj-primEnergy))
343 Cross = 0.1 * << 392 *(1.+2.*std::log(float(1.+electron_mass_c2/primEnergy)));
344 std::log((Emax_proj - G4float(pr << 393 }
345 Emax_proj / (Emin_proj << 394 }
346 (1. + 2. * std::log(G4float(1. + << 395 else {
347 } << 396 Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
348 } << 397 Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,0.);
349 else << 398 if (Emax_proj>Emin_proj) {
350 { << 399 Cross = 0.1*std::log(Emax_proj/Emin_proj);
351 Emax_proj = GetSecondAdjEnergyMaxForScatPr << 400 //+0.5*primEnergy*primEnergy(1./(Emin_proj*Emin_proj) - 1./(Emax_proj*Emax_proj)); neglected at the moment
352 Emin_proj = GetSecondAdjEnergyMinForScatPr << 401 }
353 if(Emax_proj > Emin_proj) << 402
354 { << 403
355 Cross = 0.1 * std::log(Emax_proj / Emin_ << 404 }
356 } << 405
357 } << 406 Cross*=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
358 << 407 lastCS=Cross;
359 Cross *= fCurrentMaterial->GetElectronDensit << 408 return double(Cross);
360 fLastCS = Cross; << 409 }
361 return double(Cross); << 410 ////////////////////////////////////////////////////////////////////////////////
>> 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);
362 } 417 }
363 418