| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /*
27 Authors:
28
29 Updated 15 November 2019
30
31 Updates:
32 1. Change reading method for cross section d
33 2. Add warning not to use with polarized pho
34
35 M. Omer and R. Hajima on 17 October 2016
36 contact:
37 omer.mohamed@jaea.go.jp and hajima.ryoichi@q
38 Publication Information:
39 1- M. Omer, R. Hajima, Including Delbrück s
40 Nucl. Instrum. Methods Phys. Res. Sect. B, v
41 https://doi.org/10.1016/j.nimb.2017.05.028
42 2- M. Omer, R. Hajima, Geant4 physics proces
43 JAEA Technical Report 2018-007, 2018.
44 https://doi.org/10.11484/jaea-data-code-2018
45 */
46 #include "G4JAEAElasticScatteringModel.hh"
47 #include "G4AutoLock.hh"
48 #include "G4SystemOfUnits.hh"
49
50 using namespace std;
51 namespace { G4Mutex G4JAEAElasticScatteringMod
52
53 //....oooOO0OOooo........oooOO0OOooo........oo
54 //....oooOO0OOooo........oooOO0OOooo........oo
55
56 G4PhysicsFreeVector* G4JAEAElasticScatteringMo
57 G4DataVector* G4JAEAElasticScatteringModel::ES
58
59 G4JAEAElasticScatteringModel::G4JAEAElasticSca
60 :G4VEmModel("G4JAEAElasticScatteringModel"),
61 {
62 fParticleChange = nullptr;
63 //Low energy limit for G4JAEAElasticScatteri
64 lowEnergyLimit = 10 * keV;
65
66 verboseLevel= 0;
67 // Verbosity scale for debugging purposes:
68 // 0 = nothing
69 // 1 = calculation of cross sections, file o
70 // 2 = entering in methods
71
72 if(verboseLevel > 0)
73 {
74 G4cout << "G4JAEAElasticScatteringModel
75 }
76 }
77
78 //....oooOO0OOooo........oooOO0OOooo........oo
79
80 G4JAEAElasticScatteringModel::~G4JAEAElasticSc
81 {
82 if(IsMaster()) {
83 for(G4int i=0; i<=maxZ; ++i) {
84 if(dataCS[i]) {
85 delete dataCS[i];
86 dataCS[i] = nullptr;
87 }
88 if (ES_Data[i]){
89 delete ES_Data[i];
90 ES_Data[i] = nullptr;
91 }
92 }
93 }
94 }
95
96 //....oooOO0OOooo........oooOO0OOooo........oo
97
98 void G4JAEAElasticScatteringModel::Initialise(
99 const G4DataVector& cuts)
100 {
101 if (verboseLevel > 1)
102 {
103 G4cout << "Calling Initialise() of G4JAE
104 << "Energy range: "
105 << LowEnergyLimit() / eV << " eV - "
106 << HighEnergyLimit() / GeV << " GeV"
107 << G4endl;
108 }
109
110 if(IsMaster()) {
111 // Initialise element selector
112 InitialiseElementSelectors(particle, cuts)
113
114 // Access to elements
115 const char* path = G4FindDataDir("G4LEDATA
116 G4ProductionCutsTable* theCoupleTable =
117 G4ProductionCutsTable::GetProductionCuts
118 G4int numOfCouples = (G4int)theCoupleTable
119
120 for(G4int i=0; i<numOfCouples; ++i)
121 {
122 const G4MaterialCutsCouple* couple =
123 theCoupleTable->GetMaterialCutsCouple(i);
124 const G4Material* material = couple->GetMate
125 const G4ElementVector* theElementVector = ma
126 std::size_t nelm = material->GetNumberOfElem
127
128 for (std::size_t j=0; j<nelm; ++j)
129 {
130 G4int Z = G4lrint((*theElementVector)[j]
131 if(Z < 1) { Z = 1; }
132 else if(Z > maxZ) { Z = maxZ; }
133 if( (!dataCS[Z]) ) { ReadData(Z, path);
134 }
135 }
136 }
137
138 if(isInitialised) { return; }
139 fParticleChange = GetParticleChangeForGamma(
140 isInitialised = true;
141
142 }
143
144 //....oooOO0OOooo........oooOO0OOooo........oo
145
146 void G4JAEAElasticScatteringModel::InitialiseL
147 G4VEmModel* masterModel)
148 {
149 SetElementSelectors(masterModel->GetElementS
150 }
151
152 //....oooOO0OOooo........oooOO0OOooo........oo
153
154 void G4JAEAElasticScatteringModel::ReadData(st
155 {
156 if (verboseLevel > 1)
157 {
158 G4cout << "Calling ReadData() of G4JAEAE
159 << G4endl;
160 }
161
162 if(dataCS[Z]) { return; }
163
164 const char* datadir = path;
165
166 if(!datadir)
167 {
168 datadir = G4FindDataDir("G4LEDATA");
169 if(!datadir)
170 {
171 G4Exception("G4JAEAElasticScatteringModel:
172 FatalException,
173 "Environment variable G4LEDATA not d
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;
195 }
196 else
197 {
198 if(verboseLevel > 3) {
199 G4cout << "File " << ostCS.str()
200 << " is opened by G4JAEAElasticScatte
201 }
202 }
203 if (!ES_Data[Z])
204 ES_Data[Z] = new G4DataVector();
205
206 G4float buffer_var;
207 while (ES_Data_Buffer.read(reinterpret_cast<
208 {
209 ES_Data[Z]->push_back(buffer_var);
210 }
211
212 /*
213 Writing the total cross section data to a
214 This provides an interpolation of the Ener
215 */
216
217 dataCS[Z] = new G4PhysicsFreeVector(300,0.01
218 //Note that the total cross section and ener
219 for (G4int i=0;i<300;++i)
220 dataCS[Z]->PutValues(i,10.*i*1e-3,ES_Data[
221
222 // Activation of spline interpolation
223 dataCS[Z] ->FillSecondDerivatives();
224
225 }
226
227 //....oooOO0OOooo........oooOO0OOooo........oo
228
229 G4double G4JAEAElasticScatteringModel::Compute
230 const G4ParticleDefinition*,
231 G4double GammaEnergy,
232 G4double Z, G4double,
233 G4double, G4double)
234 {
235 if (verboseLevel > 2)
236 {
237 G4cout << "G4JAEAElasticScatteringModel:
238 << G4endl;
239 }
240
241 if(GammaEnergy < lowEnergyLimit) { return 0.
242
243 G4double xs = 0.0;
244
245 G4int intZ = G4lrint(Z);
246
247 if(intZ < 1 || intZ > maxZ) { return xs; }
248
249 G4PhysicsFreeVector* pv = dataCS[intZ];
250
251 // if element was not initialised
252 // do initialisation safely for MT mode
253 if(!pv) {
254 InitialiseForElement(0, intZ);
255 pv = dataCS[intZ];
256 if(!pv) { return xs; }
257 }
258
259 G4int n = G4int(pv->GetVectorLength() - 1);
260
261 G4double e = GammaEnergy;
262 if(e >= pv->Energy(n)) {
263 xs = (*pv)[n];
264 } else if(e >= pv->Energy(0)) {
265 xs = pv->Value(e);
266 }
267
268 if(verboseLevel > 0)
269 {
270 G4cout << "****** DEBUG: tcs value for
271 << e << G4endl;
272 G4cout << " cs (Geant4 internal unit)
273 G4cout << " -> first E*E*cs value i
274 << G4endl;
275 G4cout << " -> last E*E*cs value i
276 << G4endl;
277 G4cout << "***************************
278 << G4endl;
279 }
280
281 return (xs);
282 }
283
284 //....oooOO0OOooo........oooOO0OOooo........oo
285
286 void G4JAEAElasticScatteringModel::SampleSecon
287 std::vector<G4DynamicParticle
288 const G4MaterialCutsCouple* c
289 const G4DynamicParticle* aDyn
290 G4double, G4double)
291 {
292 if (verboseLevel > 2) {
293 G4cout << "Calling SampleSecondaries() of
294 << G4endl;
295 }
296
297 G4double photonEnergy0 = aDynamicGamma->GetK
298
299 // Absorption of low-energy gamma
300 if (photonEnergy0 <= lowEnergyLimit)
301 {
302 fParticleChange->ProposeTrackStatus(fSto
303 fParticleChange->SetProposedKineticEnerg
304 fParticleChange->ProposeLocalEnergyDepos
305 return;
306 }
307
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
315 G4double polarization_magnitude=Xi1*Xi1+Xi2*
316 if ((polarization_magnitude)>0 || (Xi1*Xi1>0
317 {
318 G4cout<<"WARNING: G4JAEAElasticScatterin
319 <<G4endl;
320 G4cout<<"The event is ignored."<<G4endl;
321 return;
322 }
323
324 // Select randomly one element in the curren
325 const G4ParticleDefinition* particle = aDyn
326 const G4Element* elm = SelectRandomAtom(coup
327 G4int Z = G4lrint(elm->GetZ());
328
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
346 //Create the cummulative distribution functi
347 for (G4int i =0;i<=180;++i)
348 pdf[i]=distribution[i]/normdist;
349 cdf[0]=0;
350 G4double cdfsum =0;
351 for (G4int i=0; i<=180;++i)
352 {
353 cdfsum=cdfsum+pdf[i];
354 cdf[i]=cdfsum;
355 }
356 //Sampling the polar angle by inverse transf
357 G4double r = G4UniformRand();
358 G4double *cdfptr=lower_bound(cdf,cdf+181,r);
359 G4int cdfindex = (G4int)(cdfptr-cdf-1);
360 G4double cdfinv = (r-cdf[cdfindex])/(cdf[cdf
361 G4double theta = (cdfindex+cdfinv)/180.;
362 //polar is now ready
363 theta = theta*CLHEP::pi;
364
365
366 /* Alternative sampling using CLHEP function
367 CLHEP::RandGeneral GenDistTheta(distribut
368 G4double theta = CLHEP::pi*GenDistTheta.s
369 theta =theta*CLHEP::pi; //polar is now
370 */
371
372 //Azimuth is uniformally distributed
373 G4double phi = CLHEP::twopi*G4UniformRand()
374
375 G4ThreeVector finaldirection(sin(theta)*cos(
376 finaldirection.rotateUz(aDynamicGamma->GetMo
377 //Sampling the Final State
378 fParticleChange->ProposeMomentumDirection(fi
379 fParticleChange->SetProposedKineticEnergy(ph
380 }
381
382 //....oooOO0OOooo........oooOO0OOooo........oo
383
384 void
385 G4JAEAElasticScatteringModel::InitialiseForEle
386 G4int Z)
387 {
388 G4AutoLock l(&G4JAEAElasticScatteringModelMu
389 if(!dataCS[Z]) { ReadData(Z); }
390 l.unlock();
391 }
392
393 //....oooOO0OOooo........oooOO0OOooo........oo
394