| Geant4 Cross Reference |
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 /*
27 Authors:
28
29 Updated 15 November 2019
30
31 Updates:
32 1. Change reading method for cross section data.
33 2. Add warning not to use with polarized photons.
34
35 M. Omer and R. Hajima on 17 October 2016
36 contact:
37 omer.mohamed@jaea.go.jp and hajima.ryoichi@qst.go.jp
38 Publication Information:
39 1- M. Omer, R. Hajima, Including Delbrück scattering in Geant4,
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
42 2- M. Omer, R. Hajima, Geant4 physics process for elastic scattering of gamma-rays,
43 JAEA Technical Report 2018-007, 2018.
44 https://doi.org/10.11484/jaea-data-code-2018-007
45 */
46 #include "G4JAEAElasticScatteringModel.hh"
47 #include "G4AutoLock.hh"
48 #include "G4SystemOfUnits.hh"
49
50 using namespace std;
51 namespace { G4Mutex G4JAEAElasticScatteringModelMutex = G4MUTEX_INITIALIZER; }
52
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55
56 G4PhysicsFreeVector* G4JAEAElasticScatteringModel::dataCS[]={nullptr} ;
57 G4DataVector* G4JAEAElasticScatteringModel::ES_Data[]={nullptr};
58
59 G4JAEAElasticScatteringModel::G4JAEAElasticScatteringModel()
60 :G4VEmModel("G4JAEAElasticScatteringModel"),isInitialised(false)
61 {
62 fParticleChange = nullptr;
63 //Low energy limit for G4JAEAElasticScatteringModel process.
64 lowEnergyLimit = 10 * keV;
65
66 verboseLevel= 0;
67 // Verbosity scale for debugging purposes:
68 // 0 = nothing
69 // 1 = calculation of cross sections, file openings...
70 // 2 = entering in methods
71
72 if(verboseLevel > 0)
73 {
74 G4cout << "G4JAEAElasticScatteringModel is constructed " << G4endl;
75 }
76 }
77
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79
80 G4JAEAElasticScatteringModel::~G4JAEAElasticScatteringModel()
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........oooOO0OOooo........oooOO0OOooo....
97
98 void G4JAEAElasticScatteringModel::Initialise(const G4ParticleDefinition* particle,
99 const G4DataVector& cuts)
100 {
101 if (verboseLevel > 1)
102 {
103 G4cout << "Calling Initialise() of G4JAEAElasticScatteringModel." << G4endl
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::GetProductionCutsTable();
118 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
119
120 for(G4int i=0; i<numOfCouples; ++i)
121 {
122 const G4MaterialCutsCouple* couple =
123 theCoupleTable->GetMaterialCutsCouple(i);
124 const G4Material* material = couple->GetMaterial();
125 const G4ElementVector* theElementVector = material->GetElementVector();
126 std::size_t nelm = material->GetNumberOfElements();
127
128 for (std::size_t j=0; j<nelm; ++j)
129 {
130 G4int Z = G4lrint((*theElementVector)[j]->GetZ());
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........oooOO0OOooo........oooOO0OOooo....
145
146 void G4JAEAElasticScatteringModel::InitialiseLocal(const G4ParticleDefinition*,
147 G4VEmModel* masterModel)
148 {
149 SetElementSelectors(masterModel->GetElementSelectors());
150 }
151
152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
153
154 void G4JAEAElasticScatteringModel::ReadData(std::size_t Z, const char* path)
155 {
156 if (verboseLevel > 1)
157 {
158 G4cout << "Calling ReadData() of G4JAEAElasticScatteringModel"
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::ReadData()","em0006",
172 FatalException,
173 "Environment variable G4LEDATA not defined");
174 return;
175 }
176 }
177
178 /* Reading all data in the form of 183 * 300 array.
179 The first row is the energy, and the second row is the total cross section.
180 Rows from the 3rd to the 183rd are the differential cross section with an angular
181 resolution of 1 degree.
182 */
183 std::ostringstream ostCS;
184 ostCS << datadir << "/JAEAESData/amp_Z_" << Z ;
185 std::ifstream ES_Data_Buffer(ostCS.str().c_str(),ios::binary);
186 if( !ES_Data_Buffer.is_open() )
187 {
188 G4ExceptionDescription ed;
189 ed << "G4JAEAElasticScattertingModel data file <" << ostCS.str().c_str()
190 << "> is not opened!" << G4endl;
191 G4Exception("G4JAEAElasticScatteringModel::ReadData()","em0003",FatalException,
192 ed,
193 "G4LEDATA version should be G4EMLOW7.11 or later. Elastic Scattering Data are not loaded");
194 return;
195 }
196 else
197 {
198 if(verboseLevel > 3) {
199 G4cout << "File " << ostCS.str()
200 << " is opened by G4JAEAElasticScatteringModel" << G4endl;
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<char*>(&buffer_var),sizeof(float)))
208 {
209 ES_Data[Z]->push_back(buffer_var);
210 }
211
212 /*
213 Writing the total cross section data to a G4PhysicsFreeVector.
214 This provides an interpolation of the Energy-Total Cross Section data.
215 */
216
217 dataCS[Z] = new G4PhysicsFreeVector(300,0.01,3.,/*spine=*/true);
218 //Note that the total cross section and energy are converted to the internal units.
219 for (G4int i=0;i<300;++i)
220 dataCS[Z]->PutValues(i,10.*i*1e-3,ES_Data[Z]->at(i)*1e-22);
221
222 // Activation of spline interpolation
223 dataCS[Z] ->FillSecondDerivatives();
224
225 }
226
227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
228
229 G4double G4JAEAElasticScatteringModel::ComputeCrossSectionPerAtom(
230 const G4ParticleDefinition*,
231 G4double GammaEnergy,
232 G4double Z, G4double,
233 G4double, G4double)
234 {
235 if (verboseLevel > 2)
236 {
237 G4cout << "G4JAEAElasticScatteringModel::ComputeCrossSectionPerAtom()"
238 << G4endl;
239 }
240
241 if(GammaEnergy < lowEnergyLimit) { return 0.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 Z=" << Z << " at energy (MeV)="
271 << e << G4endl;
272 G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
273 G4cout << " -> first E*E*cs value in CS data file (iu) =" << (*pv)[0]
274 << G4endl;
275 G4cout << " -> last E*E*cs value in CS data file (iu) =" << (*pv)[n]
276 << G4endl;
277 G4cout << "*********************************************************"
278 << G4endl;
279 }
280
281 return (xs);
282 }
283
284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
285
286 void G4JAEAElasticScatteringModel::SampleSecondaries(
287 std::vector<G4DynamicParticle*>*,
288 const G4MaterialCutsCouple* couple,
289 const G4DynamicParticle* aDynamicGamma,
290 G4double, G4double)
291 {
292 if (verboseLevel > 2) {
293 G4cout << "Calling SampleSecondaries() of G4JAEAElasticScatteringModel."
294 << G4endl;
295 }
296
297 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
298
299 // Absorption of low-energy gamma
300 if (photonEnergy0 <= lowEnergyLimit)
301 {
302 fParticleChange->ProposeTrackStatus(fStopAndKill);
303 fParticleChange->SetProposedKineticEnergy(0.);
304 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
305 return;
306 }
307
308 //Warning if the incoming photon has polarization
309 G4double Xi1=0, Xi2=0, Xi3=0;
310 G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
311 Xi1=gammaPolarization0.x();
312 Xi2=gammaPolarization0.y();
313 Xi3=gammaPolarization0.z();
314
315 G4double polarization_magnitude=Xi1*Xi1+Xi2*Xi2+Xi3*Xi3;
316 if ((polarization_magnitude)>0 || (Xi1*Xi1>0) || (Xi2*Xi2>0) || (Xi3*Xi3>0))
317 {
318 G4cout<<"WARNING: G4JAEAElasticScatteringModel is only compatible with non-polarized photons."
319 <<G4endl;
320 G4cout<<"The event is ignored."<<G4endl;
321 return;
322 }
323
324 // Select randomly one element in the current material
325 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
326 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
327 G4int Z = G4lrint(elm->GetZ());
328
329 G4int energyindex=round(100*photonEnergy0)-1;
330 /*
331 Getting the normalized probablity distrbution function and
332 normalization factor to create the probability distribution function
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*(energyindex));
339 a2=ES_Data[Z]->at(4*i+1+300+181*4*(energyindex));
340 a3=ES_Data[Z]->at(4*i+2+300+181*4*(energyindex));
341 a4=ES_Data[Z]->at(4*i+3+300+181*4*(energyindex));
342 distribution[i]=a1*a1+a2*a2+a3*a3+a4*a4;
343 normdist += distribution[i];
344 }
345
346 //Create the cummulative distribution function (cdf)
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 transform uing cdf.
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[cdfindex+1]-cdf[cdfindex]);
361 G4double theta = (cdfindex+cdfinv)/180.;
362 //polar is now ready
363 theta = theta*CLHEP::pi;
364
365
366 /* Alternative sampling using CLHEP functions
367 CLHEP::RandGeneral GenDistTheta(distribution,181);
368 G4double theta = CLHEP::pi*GenDistTheta.shoot();
369 theta =theta*CLHEP::pi; //polar is now ready
370 */
371
372 //Azimuth is uniformally distributed
373 G4double phi = CLHEP::twopi*G4UniformRand();
374
375 G4ThreeVector finaldirection(sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta));
376 finaldirection.rotateUz(aDynamicGamma->GetMomentumDirection());
377 //Sampling the Final State
378 fParticleChange->ProposeMomentumDirection(finaldirection);
379 fParticleChange->SetProposedKineticEnergy(photonEnergy0);
380 }
381
382 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
383
384 void
385 G4JAEAElasticScatteringModel::InitialiseForElement(const G4ParticleDefinition*,
386 G4int Z)
387 {
388 G4AutoLock l(&G4JAEAElasticScatteringModelMutex);
389 if(!dataCS[Z]) { ReadData(Z); }
390 l.unlock();
391 }
392
393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
394