| 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 /* 26 /*
27 Authors: << 27 Authors:
28 28
29 Updated 15 November 2019 << 29 Updated 15 Novebmer 2019
30 30
31 Updates: << 31 Updates:
32 1. Change reading method for cross section d << 32 1. Change reading method for cross section data.
33 2. Add warning not to use with polarized pho << 33 2. Add warning not to use with polarized photons.
34 << 34
35 M. Omer and R. Hajima on 17 October 2016 << 35 M. Omer and R. Hajima on 17 October 2016
36 contact: << 36 contact:
37 omer.mohamed@jaea.go.jp and hajima.ryoichi@q << 37 omer.mohamed@jaea.go.jp and hajima.ryoichi@qst.go.jp
38 Publication Information: << 38 Publication Information:
39 1- M. Omer, R. Hajima, Including Delbrück s << 39 1- M. Omer, R. Hajima, Including Delbrück scattering in Geant4,
40 Nucl. Instrum. Methods Phys. Res. Sect. B, v << 40 Nucl. Instrum. Methods Phys. Res. Sect. B, vol. 405, 2017, pp. 43-49.,
41 https://doi.org/10.1016/j.nimb.2017.05.028 << 41 https://doi.org/10.1016/j.nimb.2017.05.028
42 2- M. Omer, R. Hajima, Geant4 physics proces << 42 2- M. Omer, R. Hajima, Geant4 physics process for elastic scattering of gamma-rays,
43 JAEA Technical Report 2018-007, 2018. << 43 JAEA Technical Report 2018-007, 2018.
44 https://doi.org/10.11484/jaea-data-code-2018 << 44 https://doi.org/10.11484/jaea-data-code-2018-007
45 */ 45 */
>> 46 // on base of G4LivermoreRayleighModel
>> 47
>> 48
46 #include "G4JAEAElasticScatteringModel.hh" 49 #include "G4JAEAElasticScatteringModel.hh"
47 #include "G4AutoLock.hh" <<
48 #include "G4SystemOfUnits.hh" 50 #include "G4SystemOfUnits.hh"
49 51
50 using namespace std; 52 using namespace std;
51 namespace { G4Mutex G4JAEAElasticScatteringMod <<
52 53
53 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54 //....oooOO0OOooo........oooOO0OOooo........oo 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55 56
56 G4PhysicsFreeVector* G4JAEAElasticScatteringMo << 57 G4LPhysicsFreeVector* G4JAEAElasticScatteringModel::dataCS[]={nullptr} ;
57 G4DataVector* G4JAEAElasticScatteringModel::ES 58 G4DataVector* G4JAEAElasticScatteringModel::ES_Data[]={nullptr};
58 59
59 G4JAEAElasticScatteringModel::G4JAEAElasticSca 60 G4JAEAElasticScatteringModel::G4JAEAElasticScatteringModel()
60 :G4VEmModel("G4JAEAElasticScatteringModel"), 61 :G4VEmModel("G4JAEAElasticScatteringModel"),isInitialised(false)
61 { 62 {
62 fParticleChange = nullptr; << 63 fParticleChange = 0;
63 //Low energy limit for G4JAEAElasticScatteri << 64 //Low energy limit for G4JAEAElasticScatteringModel process.
64 lowEnergyLimit = 10 * keV; 65 lowEnergyLimit = 10 * keV;
65 66
66 verboseLevel= 0; 67 verboseLevel= 0;
67 // Verbosity scale for debugging purposes: 68 // Verbosity scale for debugging purposes:
68 // 0 = nothing 69 // 0 = nothing
69 // 1 = calculation of cross sections, file o 70 // 1 = calculation of cross sections, file openings...
70 // 2 = entering in methods 71 // 2 = entering in methods
71 72
72 if(verboseLevel > 0) 73 if(verboseLevel > 0)
73 { << 74 {
74 G4cout << "G4JAEAElasticScatteringModel << 75 G4cout << "G4JAEAElasticScatteringModel is constructed " << G4endl;
75 } << 76 }
>> 77
76 } 78 }
77 79
78 //....oooOO0OOooo........oooOO0OOooo........oo 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79 81
80 G4JAEAElasticScatteringModel::~G4JAEAElasticSc 82 G4JAEAElasticScatteringModel::~G4JAEAElasticScatteringModel()
81 { 83 {
82 if(IsMaster()) { 84 if(IsMaster()) {
83 for(G4int i=0; i<=maxZ; ++i) { 85 for(G4int i=0; i<=maxZ; ++i) {
84 if(dataCS[i]) { 86 if(dataCS[i]) {
85 delete dataCS[i]; 87 delete dataCS[i];
86 dataCS[i] = nullptr; 88 dataCS[i] = nullptr;
87 } 89 }
88 if (ES_Data[i]){ 90 if (ES_Data[i]){
89 delete ES_Data[i]; 91 delete ES_Data[i];
90 ES_Data[i] = nullptr; 92 ES_Data[i] = nullptr;
91 } 93 }
92 } 94 }
93 } 95 }
94 } 96 }
95 97
96 //....oooOO0OOooo........oooOO0OOooo........oo 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97 99
98 void G4JAEAElasticScatteringModel::Initialise( 100 void G4JAEAElasticScatteringModel::Initialise(const G4ParticleDefinition* particle,
99 const G4DataVector& cuts) << 101 const G4DataVector& cuts)
100 { 102 {
101 if (verboseLevel > 1) 103 if (verboseLevel > 1)
102 { << 104 {
103 G4cout << "Calling Initialise() of G4JAE << 105 G4cout << "Calling Initialise() of G4JAEAElasticScatteringModel." << G4endl
104 << "Energy range: " << 106 << "Energy range: "
105 << LowEnergyLimit() / eV << " eV - " << 107 << LowEnergyLimit() / eV << " eV - "
106 << HighEnergyLimit() / GeV << " GeV" << 108 << HighEnergyLimit() / GeV << " GeV"
107 << G4endl; << 109 << G4endl;
108 } << 110 }
109 111
110 if(IsMaster()) { 112 if(IsMaster()) {
>> 113
111 // Initialise element selector 114 // Initialise element selector
112 InitialiseElementSelectors(particle, cuts) 115 InitialiseElementSelectors(particle, cuts);
113 116
114 // Access to elements 117 // Access to elements
115 const char* path = G4FindDataDir("G4LEDATA << 118 char* path = std::getenv("G4LEDATA");
116 G4ProductionCutsTable* theCoupleTable = 119 G4ProductionCutsTable* theCoupleTable =
117 G4ProductionCutsTable::GetProductionCuts 120 G4ProductionCutsTable::GetProductionCutsTable();
118 G4int numOfCouples = (G4int)theCoupleTable << 121 G4int numOfCouples = theCoupleTable->GetTableSize();
119 122
120 for(G4int i=0; i<numOfCouples; ++i) 123 for(G4int i=0; i<numOfCouples; ++i)
121 { 124 {
122 const G4MaterialCutsCouple* couple = 125 const G4MaterialCutsCouple* couple =
123 theCoupleTable->GetMaterialCutsCouple(i); 126 theCoupleTable->GetMaterialCutsCouple(i);
124 const G4Material* material = couple->GetMate 127 const G4Material* material = couple->GetMaterial();
125 const G4ElementVector* theElementVector = ma 128 const G4ElementVector* theElementVector = material->GetElementVector();
126 std::size_t nelm = material->GetNumberOfElem << 129 G4int nelm = material->GetNumberOfElements();
127 130
128 for (std::size_t j=0; j<nelm; ++j) << 131 for (G4int j=0; j<nelm; ++j)
129 { 132 {
130 G4int Z = G4lrint((*theElementVector)[j] 133 G4int Z = G4lrint((*theElementVector)[j]->GetZ());
131 if(Z < 1) { Z = 1; } 134 if(Z < 1) { Z = 1; }
132 else if(Z > maxZ) { Z = maxZ; } 135 else if(Z > maxZ) { Z = maxZ; }
133 if( (!dataCS[Z]) ) { ReadData(Z, path); 136 if( (!dataCS[Z]) ) { ReadData(Z, path); }
134 } 137 }
135 } 138 }
136 } 139 }
137 140
138 if(isInitialised) { return; } 141 if(isInitialised) { return; }
139 fParticleChange = GetParticleChangeForGamma( 142 fParticleChange = GetParticleChangeForGamma();
140 isInitialised = true; 143 isInitialised = true;
141 144
142 } 145 }
143 146
144 //....oooOO0OOooo........oooOO0OOooo........oo 147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
145 148
146 void G4JAEAElasticScatteringModel::InitialiseL 149 void G4JAEAElasticScatteringModel::InitialiseLocal(const G4ParticleDefinition*,
147 G4VEmModel* masterModel) << 150 G4VEmModel* masterModel)
148 { 151 {
149 SetElementSelectors(masterModel->GetElementS 152 SetElementSelectors(masterModel->GetElementSelectors());
150 } 153 }
151 154
152 //....oooOO0OOooo........oooOO0OOooo........oo 155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
153 156
154 void G4JAEAElasticScatteringModel::ReadData(st << 157 void G4JAEAElasticScatteringModel::ReadData(size_t Z, const char* path)
155 { 158 {
156 if (verboseLevel > 1) << 159
157 { << 160 if (verboseLevel > 1)
158 G4cout << "Calling ReadData() of G4JAEAE << 161 {
159 << G4endl; << 162 G4cout << "Calling ReadData() of G4JAEAElasticScatteringModel"
160 } << 163 << G4endl;
>> 164 }
161 165
162 if(dataCS[Z]) { return; } 166 if(dataCS[Z]) { return; }
163 167
164 const char* datadir = path; 168 const char* datadir = path;
165 169
166 if(!datadir) 170 if(!datadir)
167 { << 171 {
168 datadir = G4FindDataDir("G4LEDATA"); << 172 datadir = std::getenv("G4LEDATA");
169 if(!datadir) << 173 if(!datadir)
170 { << 174 {
171 G4Exception("G4JAEAElasticScatteringModel: << 175 G4Exception("G4JAEAElasticScatteringModel::ReadData()","em0006",
172 FatalException, << 176 FatalException,
173 "Environment variable G4LEDATA not d << 177 "Environment variable G4LEDATA not defined");
174 return; <<
175 } <<
176 } <<
177 <<
178 /* Reading all data in the form of 183 * 300 <<
179 The first row is the energy, and the seco <<
180 Rows from the 3rd to the 183rd are the di <<
181 resolution of 1 degree. <<
182 */ <<
183 std::ostringstream ostCS; <<
184 ostCS << datadir << "/JAEAESData/amp_Z_" << <<
185 std::ifstream ES_Data_Buffer(ostCS.str().c_s <<
186 if( !ES_Data_Buffer.is_open() ) <<
187 { <<
188 G4ExceptionDescription ed; <<
189 ed << "G4JAEAElasticScattertingModel dat <<
190 << "> is not opened!" << G4endl; <<
191 G4Exception("G4JAEAElasticScatteringMode <<
192 ed, <<
193 "G4LEDATA version should be G4EMLOW7.11 <<
194 return; 178 return;
195 } 179 }
196 else << 180 }
197 { << 181
198 if(verboseLevel > 3) { << 182 /* Reading all data in the form of 183 * 300 array.
199 G4cout << "File " << ostCS.str() << 183 The first row is the energy, and the second row is the total cross section.
200 << " is opened by G4JAEAElasticScatte << 184 Rows from the 3rd to the 183rd are the differential cross section with an angular resolution of 1 degree.
201 } << 185 */
202 } << 186
203 if (!ES_Data[Z]) << 187
204 ES_Data[Z] = new G4DataVector(); << 188 std::ostringstream ostCS;
>> 189 ostCS << datadir << "/JAEAESData/amp_Z_" << Z ;
>> 190 std::ifstream ES_Data_Buffer(ostCS.str().c_str(),ios::binary);
>> 191 if( !ES_Data_Buffer.is_open() )
>> 192 {
>> 193 G4ExceptionDescription ed;
>> 194 ed << "G4JAEAElasticScattertingModel data file <" << ostCS.str().c_str()
>> 195 << "> is not opened!" << G4endl;
>> 196 G4Exception("G4JAEAElasticScatteringModel::ReadData()","em0003",FatalException,
>> 197 ed,
>> 198 "G4LEDATA version should be G4EMLOW7.11 or later. Elastic Scattering Data are not loaded");
>> 199 return;
>> 200 }
>> 201 else
>> 202 {
>> 203 if(verboseLevel > 3) {
>> 204 G4cout << "File " << ostCS.str()
>> 205 << " is opened by G4JAEAElasticScatteringModel" << G4endl;
>> 206 }
>> 207 }
>> 208 if (!ES_Data[Z])
>> 209 ES_Data[Z] = new G4DataVector();
205 210
206 G4float buffer_var; <<
207 while (ES_Data_Buffer.read(reinterpret_cast< <<
208 { <<
209 ES_Data[Z]->push_back(buffer_var); <<
210 } <<
211 211
212 /* << 212 G4float buffer_var;
213 Writing the total cross section data to a << 213 while (ES_Data_Buffer.read(reinterpret_cast<char*>(&buffer_var),sizeof(float)))
214 This provides an interpolation of the Ener << 214 {
215 */ << 215 ES_Data[Z]->push_back(buffer_var);
216 << 216 }
217 dataCS[Z] = new G4PhysicsFreeVector(300,0.01 << 217
218 //Note that the total cross section and ener << 218
219 for (G4int i=0;i<300;++i) << 219
220 dataCS[Z]->PutValues(i,10.*i*1e-3,ES_Data[ << 220 /*
221 << 221 Writing the total cross section data to a G4LPhysicsFreeVector.
222 // Activation of spline interpolation << 222 This provides an interpolation of the Energy-Total Cross Section data.
223 dataCS[Z] ->FillSecondDerivatives(); << 223 */
224 << 224
>> 225 dataCS[Z] = new G4LPhysicsFreeVector(300,0.01,3.);
>> 226 //Note that the total cross section and energy are converted to the internal units.
>> 227 for (G4int i=0;i<300;++i)
>> 228 dataCS[Z]->PutValue(i,10.*i*1e-3,ES_Data[Z]->at(i)*1e-22);
>> 229
>> 230 // Activation of spline interpolation
>> 231 dataCS[Z] ->SetSpline(true);
>> 232
>> 233
>> 234
225 } 235 }
226 236
227 //....oooOO0OOooo........oooOO0OOooo........oo 237 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
228 238
229 G4double G4JAEAElasticScatteringModel::Compute 239 G4double G4JAEAElasticScatteringModel::ComputeCrossSectionPerAtom(
230 const G4ParticleDefinition*, << 240 const G4ParticleDefinition*,
231 G4double GammaEnergy, << 241 G4double GammaEnergy,
232 G4double Z, G4double, << 242 G4double Z, G4double,
233 G4double, G4double) << 243 G4double, G4double)
234 { 244 {
>> 245
235 if (verboseLevel > 2) 246 if (verboseLevel > 2)
236 { << 247 {
237 G4cout << "G4JAEAElasticScatteringModel: << 248 G4cout << "G4JAEAElasticScatteringModel::ComputeCrossSectionPerAtom()"
238 << G4endl; << 249 << G4endl;
239 } << 250 }
240 251
241 if(GammaEnergy < lowEnergyLimit) { return 0. 252 if(GammaEnergy < lowEnergyLimit) { return 0.0; }
242 253
243 G4double xs = 0.0; 254 G4double xs = 0.0;
244 255
245 G4int intZ = G4lrint(Z); 256 G4int intZ = G4lrint(Z);
246 257
247 if(intZ < 1 || intZ > maxZ) { return xs; } 258 if(intZ < 1 || intZ > maxZ) { return xs; }
248 259
249 G4PhysicsFreeVector* pv = dataCS[intZ]; << 260 G4LPhysicsFreeVector* pv = dataCS[intZ];
250 261
251 // if element was not initialised 262 // if element was not initialised
252 // do initialisation safely for MT mode 263 // do initialisation safely for MT mode
253 if(!pv) { 264 if(!pv) {
254 InitialiseForElement(0, intZ); 265 InitialiseForElement(0, intZ);
255 pv = dataCS[intZ]; 266 pv = dataCS[intZ];
256 if(!pv) { return xs; } 267 if(!pv) { return xs; }
257 } 268 }
258 269
259 G4int n = G4int(pv->GetVectorLength() - 1); << 270 G4int n = pv->GetVectorLength() - 1;
260 271
261 G4double e = GammaEnergy; 272 G4double e = GammaEnergy;
262 if(e >= pv->Energy(n)) { 273 if(e >= pv->Energy(n)) {
263 xs = (*pv)[n]; 274 xs = (*pv)[n];
264 } else if(e >= pv->Energy(0)) { 275 } else if(e >= pv->Energy(0)) {
265 xs = pv->Value(e); 276 xs = pv->Value(e);
266 } 277 }
267 278
268 if(verboseLevel > 0) 279 if(verboseLevel > 0)
269 { << 280 {
270 G4cout << "****** DEBUG: tcs value for << 281 G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)="
271 << e << G4endl; << 282 << e << G4endl;
272 G4cout << " cs (Geant4 internal unit) << 283 G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
273 G4cout << " -> first E*E*cs value i << 284 G4cout << " -> first E*E*cs value in CS data file (iu) =" << (*pv)[0]
274 << G4endl; << 285 << G4endl;
275 G4cout << " -> last E*E*cs value i << 286 G4cout << " -> last E*E*cs value in CS data file (iu) =" << (*pv)[n]
276 << G4endl; << 287 << G4endl;
277 G4cout << "*************************** << 288 G4cout << "*********************************************************"
278 << G4endl; << 289 << G4endl;
279 } << 290 }
280 291
281 return (xs); 292 return (xs);
282 } 293 }
283 294
284 //....oooOO0OOooo........oooOO0OOooo........oo 295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
285 296
286 void G4JAEAElasticScatteringModel::SampleSecon 297 void G4JAEAElasticScatteringModel::SampleSecondaries(
287 std::vector<G4DynamicParticle << 298 std::vector<G4DynamicParticle*>*,
288 const G4MaterialCutsCouple* c << 299 const G4MaterialCutsCouple* couple,
289 const G4DynamicParticle* aDyn << 300 const G4DynamicParticle* aDynamicGamma,
290 G4double, G4double) << 301 G4double, G4double)
291 { 302 {
292 if (verboseLevel > 2) { 303 if (verboseLevel > 2) {
293 G4cout << "Calling SampleSecondaries() of 304 G4cout << "Calling SampleSecondaries() of G4JAEAElasticScatteringModel."
294 << G4endl; 305 << G4endl;
295 } 306 }
296 307
297 G4double photonEnergy0 = aDynamicGamma->GetK 308 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
298 309
299 // Absorption of low-energy gamma 310 // Absorption of low-energy gamma
300 if (photonEnergy0 <= lowEnergyLimit) 311 if (photonEnergy0 <= lowEnergyLimit)
301 { 312 {
302 fParticleChange->ProposeTrackStatus(fSto 313 fParticleChange->ProposeTrackStatus(fStopAndKill);
303 fParticleChange->SetProposedKineticEnerg 314 fParticleChange->SetProposedKineticEnergy(0.);
304 fParticleChange->ProposeLocalEnergyDepos 315 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
305 return; 316 return;
306 } 317 }
307 318
308 //Warning if the incoming photon has polariz <<
309 G4double Xi1=0, Xi2=0, Xi3=0; <<
310 G4ThreeVector gammaPolarization0 = aDynamicG <<
311 Xi1=gammaPolarization0.x(); <<
312 Xi2=gammaPolarization0.y(); <<
313 Xi3=gammaPolarization0.z(); <<
314 319
315 G4double polarization_magnitude=Xi1*Xi1+Xi2* << 320 //Warning if the incoming photon has polarization
316 if ((polarization_magnitude)>0 || (Xi1*Xi1>0 << 321
317 { << 322 G4double Xi1=0, Xi2=0, Xi3=0;
318 G4cout<<"WARNING: G4JAEAElasticScatterin << 323 G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
319 <<G4endl; << 324 Xi1=gammaPolarization0.x();
320 G4cout<<"The event is ignored."<<G4endl; << 325 Xi2=gammaPolarization0.y();
321 return; << 326 Xi3=gammaPolarization0.z();
>> 327
>> 328 G4double polarization_magnitude=Xi1*Xi1+Xi2*Xi2+Xi3*Xi3;
>> 329 if ((polarization_magnitude)>0 || (Xi1*Xi1>0) || (Xi2*Xi2>0) || (Xi3*Xi3>0))
>> 330 {
>> 331 G4cout<<"WARNING: G4JAEAElasticScatteringModel is only compatible with non-polarized photons."<<G4endl;
>> 332 G4cout<<"The event is ignored."<<G4endl;
>> 333 return;
322 } 334 }
323 335
324 // Select randomly one element in the curren 336 // Select randomly one element in the current material
325 const G4ParticleDefinition* particle = aDyn 337 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
326 const G4Element* elm = SelectRandomAtom(coup 338 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
327 G4int Z = G4lrint(elm->GetZ()); 339 G4int Z = G4lrint(elm->GetZ());
328 340
329 G4int energyindex=round(100*photonEnergy0)-1 <<
330 /* <<
331 Getting the normalized probablity distrbut <<
332 normalization factor to create the probabi <<
333 */ <<
334 G4double a1=0, a2=0, a3=0,a4=0; <<
335 G4double normdist=0; <<
336 for (G4int i=0;i<=180;++i) <<
337 { <<
338 a1=ES_Data[Z]->at(4*i+300+181*4*(energyi <<
339 a2=ES_Data[Z]->at(4*i+1+300+181*4*(energ <<
340 a3=ES_Data[Z]->at(4*i+2+300+181*4*(energ <<
341 a4=ES_Data[Z]->at(4*i+3+300+181*4*(energ <<
342 distribution[i]=a1*a1+a2*a2+a3*a3+a4*a4; <<
343 normdist += distribution[i]; <<
344 } <<
345 341
346 //Create the cummulative distribution functi << 342 G4int energyindex=round(100*photonEnergy0)-1;
347 for (G4int i =0;i<=180;++i) << 343 /*
348 pdf[i]=distribution[i]/normdist; << 344 Getting the normalized probablity distrbution function and
349 cdf[0]=0; << 345 normalization factor to create the probability distribution function
350 G4double cdfsum =0; << 346 */
351 for (G4int i=0; i<=180;++i) << 347 G4double a1=0, a2=0, a3=0,a4=0;
352 { << 348 G4double normdist=0;
353 cdfsum=cdfsum+pdf[i]; << 349 for (G4int i=0;i<=180;++i)
354 cdf[i]=cdfsum; << 350 {
355 } << 351 a1=ES_Data[Z]->at(4*i+300+181*4*(energyindex));
356 //Sampling the polar angle by inverse transf << 352 a2=ES_Data[Z]->at(4*i+1+300+181*4*(energyindex));
357 G4double r = G4UniformRand(); << 353 a3=ES_Data[Z]->at(4*i+2+300+181*4*(energyindex));
358 G4double *cdfptr=lower_bound(cdf,cdf+181,r); << 354 a4=ES_Data[Z]->at(4*i+3+300+181*4*(energyindex));
359 G4int cdfindex = (G4int)(cdfptr-cdf-1); << 355 distribution[i]=a1*a1+a2*a2+a3*a3+a4*a4;
360 G4double cdfinv = (r-cdf[cdfindex])/(cdf[cdf << 356 normdist += distribution[i];
361 G4double theta = (cdfindex+cdfinv)/180.; << 357 }
362 //polar is now ready << 358
363 theta = theta*CLHEP::pi; << 359
364 << 360 //Create the cummulative distribution function (cdf)
365 << 361 for (G4int i =0;i<=180;++i)
366 /* Alternative sampling using CLHEP function << 362 pdf[i]=distribution[i]/normdist;
367 CLHEP::RandGeneral GenDistTheta(distribut << 363 cdf[0]=0;
368 G4double theta = CLHEP::pi*GenDistTheta.s << 364 G4double cdfsum =0;
369 theta =theta*CLHEP::pi; //polar is now << 365 for (G4int i=0; i<=180;++i)
370 */ << 366 {
371 << 367 cdfsum=cdfsum+pdf[i];
372 //Azimuth is uniformally distributed << 368 cdf[i]=cdfsum;
373 G4double phi = CLHEP::twopi*G4UniformRand() << 369 }
374 << 370 //Sampling the polar angle by inverse transform uing cdf.
375 G4ThreeVector finaldirection(sin(theta)*cos( << 371 G4double r = G4UniformRand();
376 finaldirection.rotateUz(aDynamicGamma->GetMo << 372 G4double *cdfptr=lower_bound(cdf,cdf+181,r);
377 //Sampling the Final State << 373 G4int cdfindex = (G4int)(cdfptr-cdf-1);
378 fParticleChange->ProposeMomentumDirection(fi << 374 G4double cdfinv = (r-cdf[cdfindex])/(cdf[cdfindex+1]-cdf[cdfindex]);
379 fParticleChange->SetProposedKineticEnergy(ph << 375 G4double theta = (cdfindex+cdfinv)/180.;
>> 376 //polar is now ready
>> 377 theta = theta*CLHEP::pi;
>> 378
>> 379
>> 380 /* Alternative sampling using CLHEP functions
>> 381 CLHEP::RandGeneral GenDistTheta(distribution,181);
>> 382 G4double theta = CLHEP::pi*GenDistTheta.shoot();
>> 383 theta =theta*CLHEP::pi; //polar is now ready
>> 384 */
>> 385
>> 386 //Azimuth is uniformally distributed
>> 387 G4double phi = CLHEP::twopi*G4UniformRand();
>> 388
>> 389 G4ThreeVector finaldirection(sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta));
>> 390 finaldirection.rotateUz(aDynamicGamma->GetMomentumDirection());
>> 391 //Sampling the Final State
>> 392 fParticleChange->ProposeMomentumDirection(finaldirection);
>> 393 fParticleChange->SetProposedKineticEnergy(photonEnergy0);
380 } 394 }
381 395
382 //....oooOO0OOooo........oooOO0OOooo........oo 396 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
383 397
>> 398 #include "G4AutoLock.hh"
>> 399 namespace { G4Mutex G4JAEAElasticScatteringModelMutex = G4MUTEX_INITIALIZER; }
>> 400
384 void 401 void
385 G4JAEAElasticScatteringModel::InitialiseForEle 402 G4JAEAElasticScatteringModel::InitialiseForElement(const G4ParticleDefinition*,
386 G4int Z) << 403 G4int Z)
387 { 404 {
388 G4AutoLock l(&G4JAEAElasticScatteringModelMu 405 G4AutoLock l(&G4JAEAElasticScatteringModelMutex);
389 if(!dataCS[Z]) { ReadData(Z); } 406 if(!dataCS[Z]) { ReadData(Z); }
390 l.unlock(); 407 l.unlock();
391 } 408 }
392 409
393 //....oooOO0OOooo........oooOO0OOooo........oo 410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
394 411