| 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 // --------------------------------------------------------------
28 //
29 // File name: G4PenelopeBremsstrahlungAngular
30 //
31 // Author: Luciano Pandola
32 //
33 // Creation date: November 2010
34 //
35 // History:
36 // -----------
37 // 23 Nov 2010 L. Pandola 1st implementation
38 // 24 May 2011 L. Pandola Renamed (make v2008 as default Penelope)
39 // 13 Mar 2012 L. Pandola Made a derived class of G4VEmAngularDistribution
40 // and update the interface accordingly
41 // 18 Jul 2012 L. Pandola Migrated to the new basic interface of G4VEmAngularDistribution
42 // Now returns a G4ThreeVector and takes care of the rotation
43 // 03 Oct 2013 L. Pandola Migrated to MT: only the master model handles tables
44 // 17 Oct 2013 L. Pandola Partially revert MT migration. The angular generator is kept as
45 // thread-local, and each worker has full access to it.
46 //
47 //----------------------------------------------------------------
48
49 #include "G4PenelopeBremsstrahlungAngular.hh"
50
51 #include "globals.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "G4PhysicsFreeVector.hh"
55 #include "G4PhysicsTable.hh"
56 #include "G4Material.hh"
57 #include "Randomize.hh"
58 #include "G4Exp.hh"
59
60 G4PenelopeBremsstrahlungAngular::G4PenelopeBremsstrahlungAngular() :
61 G4VEmAngularDistribution("Penelope"), fEffectiveZSq(nullptr),
62 fLorentzTables1(nullptr),fLorentzTables2(nullptr)
63
64 {
65 fDataRead = false;
66 fVerbosityLevel = 0;
67 }
68
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
70
71 G4PenelopeBremsstrahlungAngular::~G4PenelopeBremsstrahlungAngular()
72 {
73 ClearTables();
74 }
75
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77
78 void G4PenelopeBremsstrahlungAngular::Initialize()
79 {
80 ClearTables();
81 }
82
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84
85 void G4PenelopeBremsstrahlungAngular::ClearTables()
86 {
87 if (fLorentzTables1)
88 {
89 for (auto j=fLorentzTables1->begin(); j != fLorentzTables1->end(); ++j)
90 {
91 G4PhysicsTable* tab = j->second;
92 tab->clearAndDestroy();
93 delete tab;
94 }
95 fLorentzTables1->clear();
96 delete fLorentzTables1;
97 fLorentzTables1 = nullptr;
98 }
99
100 if (fLorentzTables2)
101 {
102 for (auto j=fLorentzTables2->begin(); j != fLorentzTables2->end(); ++j)
103 {
104 G4PhysicsTable* tab = j->second;
105 tab->clearAndDestroy();
106 delete tab;
107 }
108 fLorentzTables2->clear();
109 delete fLorentzTables2;
110 fLorentzTables2 = nullptr;
111 }
112 if (fEffectiveZSq)
113 {
114 delete fEffectiveZSq;
115 fEffectiveZSq = nullptr;
116 }
117 }
118
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120
121 void G4PenelopeBremsstrahlungAngular::ReadDataFile()
122 {
123 //Read information from DataBase file
124 const char* path = G4FindDataDir("G4LEDATA");
125 if (!path)
126 {
127 G4String excep =
128 "G4PenelopeBremsstrahlungAngular - G4LEDATA environment variable not set!";
129 G4Exception("G4PenelopeBremsstrahlungAngular::ReadDataFile()",
130 "em0006",FatalException,excep);
131 return;
132 }
133 G4String pathString(path);
134 G4String pathFile = pathString + "/penelope/bremsstrahlung/pdbrang.p08";
135 std::ifstream file(pathFile);
136
137 if (!file.is_open())
138 {
139 G4String excep = "G4PenelopeBremsstrahlungAngular - data file " + pathFile + " not found!";
140 G4Exception("G4PenelopeBremsstrahlungAngular::ReadDataFile()",
141 "em0003",FatalException,excep);
142 return;
143 }
144 G4int i=0,j=0,k=0; // i=index for Z, j=index for E, k=index for K
145
146 for (k=0;k<fNumberofKPoints;k++)
147 for (i=0;i<fNumberofZPoints;i++)
148 for (j=0;j<fNumberofEPoints;j++)
149 {
150 G4double a1,a2;
151 G4int ik1,iz1,ie1;
152 G4double zr,er,kr;
153 file >> iz1 >> ie1 >> ik1 >> zr >> er >> kr >> a1 >> a2;
154 //check the data are correct
155 if ((iz1-1 == i) && (ik1-1 == k) && (ie1-1 == j))
156 {
157 fQQ1[i][j][k]=a1;
158 fQQ2[i][j][k]=a2;
159 }
160 else
161 {
162 G4ExceptionDescription ed;
163 ed << "Corrupted data file " << pathFile << "?" << G4endl;
164 G4Exception("G4PenelopeBremsstrahlungAngular::ReadDataFile()",
165 "em0005",FatalException,ed);
166 }
167 }
168 file.close();
169 fDataRead = true;
170 }
171
172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
173
174 void G4PenelopeBremsstrahlungAngular::PrepareTables(const G4Material* material,G4bool /*isMaster*/ )
175 {
176 //Unused at the moment: the G4PenelopeBremsstrahlungAngular is thread-local, so each worker
177 //builds its own version of the tables.
178
179 //Check if data file has already been read
180 if (!fDataRead)
181 {
182 ReadDataFile();
183 if (!fDataRead)
184 G4Exception("G4PenelopeBremsstrahlungAngular::PrepareInterpolationTables()",
185 "em2001",FatalException,"Unable to build interpolation table");
186 }
187
188 if (!fLorentzTables1)
189 fLorentzTables1 = new std::map<G4double,G4PhysicsTable*>;
190 if (!fLorentzTables2)
191 fLorentzTables2 = new std::map<G4double,G4PhysicsTable*>;
192
193 G4double Zmat = CalculateEffectiveZ(material);
194
195 const G4int reducedEnergyGrid=21;
196 //Support arrays.
197 G4double betas[fNumberofEPoints]; //betas for interpolation
198 //tables for interpolation
199 G4double Q1[fNumberofEPoints][fNumberofKPoints];
200 G4double Q2[fNumberofEPoints][fNumberofKPoints];
201 //expanded tables for interpolation
202 G4double Q1E[fNumberofEPoints][reducedEnergyGrid];
203 G4double Q2E[fNumberofEPoints][reducedEnergyGrid];
204 G4double pZ[fNumberofZPoints] = {2.0,8.0,13.0,47.0,79.0,92.0};
205
206 G4int i=0,j=0,k=0; // i=index for Z, j=index for E, k=index for K
207 //Interpolation in Z
208 for (i=0;i<fNumberofEPoints;i++)
209 {
210 for (j=0;j<fNumberofKPoints;j++)
211 {
212 G4PhysicsFreeVector* QQ1vector =
213 new G4PhysicsFreeVector(fNumberofZPoints,/*spline =*/true);
214 G4PhysicsFreeVector* QQ2vector =
215 new G4PhysicsFreeVector(fNumberofZPoints,/*spline =*/true);
216
217 //fill vectors
218 for (k=0;k<fNumberofZPoints;k++)
219 {
220 QQ1vector->PutValues(k,pZ[k],G4Log(fQQ1[k][i][j]));
221 QQ2vector->PutValues(k,pZ[k],fQQ2[k][i][j]);
222 }
223 //Filled: now calculate derivatives
224 QQ1vector->FillSecondDerivatives();
225 QQ2vector->FillSecondDerivatives();
226
227 Q1[i][j]= G4Exp(QQ1vector->Value(Zmat));
228 Q2[i][j]=QQ2vector->Value(Zmat);
229 delete QQ1vector;
230 delete QQ2vector;
231 }
232 }
233 G4double pE[fNumberofEPoints] = {1.0e-03*MeV,5.0e-03*MeV,1.0e-02*MeV,5.0e-02*MeV,
234 1.0e-01*MeV,5.0e-01*MeV};
235 G4double pK[fNumberofKPoints] = {0.0,0.6,0.8,0.95};
236 G4double ppK[reducedEnergyGrid];
237
238 for(i=0;i<reducedEnergyGrid;i++)
239 ppK[i]=((G4double) i) * 0.05;
240
241
242 for(i=0;i<fNumberofEPoints;i++)
243 betas[i]=std::sqrt(pE[i]*(pE[i]+2*electron_mass_c2))/(pE[i]+electron_mass_c2);
244
245
246 for (i=0;i<fNumberofEPoints;i++)
247 {
248 for (j=0;j<fNumberofKPoints;j++)
249 Q1[i][j]=Q1[i][j]/Zmat;
250 }
251
252 //Expanded table of distribution parameters
253 for (i=0;i<fNumberofEPoints;i++)
254 {
255 G4PhysicsFreeVector* Q1vector = new G4PhysicsFreeVector(fNumberofKPoints);
256 G4PhysicsFreeVector* Q2vector = new G4PhysicsFreeVector(fNumberofKPoints);
257
258 for (j=0;j<fNumberofKPoints;j++)
259 {
260 Q1vector->PutValues(j,pK[j],G4Log(Q1[i][j])); //logarithmic
261 Q2vector->PutValues(j,pK[j],Q2[i][j]);
262 }
263
264 for (j=0;j<reducedEnergyGrid;j++)
265 {
266 Q1E[i][j]=Q1vector->Value(ppK[j]);
267 Q2E[i][j]=Q2vector->Value(ppK[j]);
268 }
269 delete Q1vector;
270 delete Q2vector;
271 }
272 //
273 //TABLES to be stored
274 //
275 G4PhysicsTable* theTable1 = new G4PhysicsTable();
276 G4PhysicsTable* theTable2 = new G4PhysicsTable();
277 // the table will contain reducedEnergyGrid G4PhysicsFreeVectors with different
278 // values of k,
279 // Each of the G4PhysicsFreeVectors has a profile of
280 // y vs. E
281 //
282 //reserve space of the vectors.
283 for (j=0;j<reducedEnergyGrid;j++)
284 {
285 G4PhysicsFreeVector* thevec = new G4PhysicsFreeVector(fNumberofEPoints,/*spline=*/true);
286 theTable1->push_back(thevec);
287 G4PhysicsFreeVector* thevec2 = new G4PhysicsFreeVector(fNumberofEPoints,/*spline=*/true);
288 theTable2->push_back(thevec2);
289 }
290
291 for (j=0;j<reducedEnergyGrid;j++)
292 {
293 G4PhysicsFreeVector* thevec = (G4PhysicsFreeVector*) (*theTable1)[j];
294 G4PhysicsFreeVector* thevec2 = (G4PhysicsFreeVector*) (*theTable2)[j];
295 for (i=0;i<fNumberofEPoints;i++)
296 {
297 thevec->PutValues(i,betas[i],Q1E[i][j]);
298 thevec2->PutValues(i,betas[i],Q2E[i][j]);
299 }
300 //Vectors are filled: calculate derivatives
301 thevec->FillSecondDerivatives();
302 thevec2->FillSecondDerivatives();
303 }
304
305 if (fLorentzTables1 && fLorentzTables2)
306 {
307 fLorentzTables1->insert(std::make_pair(Zmat,theTable1));
308 fLorentzTables2->insert(std::make_pair(Zmat,theTable2));
309 }
310 else
311 {
312 G4ExceptionDescription ed;
313 ed << "Unable to create tables of Lorentz coefficients for " << G4endl;
314 ed << "<Z>= " << Zmat << " in G4PenelopeBremsstrahlungAngular" << G4endl;
315 delete theTable1;
316 delete theTable2;
317 G4Exception("G4PenelopeBremsstrahlungAngular::PrepareInterpolationTables()",
318 "em2005",FatalException,ed);
319 }
320
321 return;
322 }
323
324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
325
326 G4ThreeVector& G4PenelopeBremsstrahlungAngular::SampleDirection(const G4DynamicParticle* dp,
327 G4double eGamma,
328 G4int,
329 const G4Material* material)
330 {
331 if (!material)
332 {
333 G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()",
334 "em2040",FatalException,"The pointer to G4Material* is nullptr");
335 return fLocalDirection;
336 }
337
338 //Retrieve the effective Z
339 G4double Zmat = 0;
340
341 //The model might be initialized incorrectly, if the angular generator is not used with the
342 //G4PenelopeBremsstrahungModel: make sure it works also with other models.
343 if (!fEffectiveZSq)
344 {
345 G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()",
346 "em2040",JustWarning,"EffectiveZSq table does not exist: create it");
347 PrepareTables(material,false);
348 //return fLocalDirection;
349 }
350
351 //found in the table: return it
352 if (fEffectiveZSq->count(material))
353 Zmat = fEffectiveZSq->find(material)->second;
354 else //this can happen in unit tests or when the AngModel is coupled with bremsstrahlunh
355 //models other than G4PenelopeBremsstrahungModel
356 {
357 G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()",
358 "em2040",JustWarning,"Material not found in the effectiveZ table");
359 PrepareTables(material,false);
360 Zmat = fEffectiveZSq->find(material)->second;
361 // return fLocalDirection;
362 }
363
364 if (fVerbosityLevel > 0)
365 {
366 G4cout << "Effective <Z> for material : " << material->GetName() <<
367 " = " << Zmat << G4endl;
368 }
369
370 G4double ePrimary = dp->GetKineticEnergy();
371
372 G4double beta = std::sqrt(ePrimary*(ePrimary+2*electron_mass_c2))/
373 (ePrimary+electron_mass_c2);
374 G4double cdt = 0;
375 G4double sinTheta = 0;
376 G4double phi = 0;
377
378 //Use a pure dipole distribution for energy above 500 keV
379 if (ePrimary > 500*keV)
380 {
381 cdt = 2.0*G4UniformRand() - 1.0;
382 if (G4UniformRand() > 0.75)
383 {
384 if (cdt<0)
385 cdt = -1.0*std::pow(-cdt,1./3.);
386 else
387 cdt = std::pow(cdt,1./3.);
388 }
389 cdt = (cdt+beta)/(1.0+beta*cdt);
390 //Get primary kinematics
391 sinTheta = std::sqrt(1. - cdt*cdt);
392 phi = twopi * G4UniformRand();
393 fLocalDirection.set(sinTheta* std::cos(phi),
394 sinTheta* std::sin(phi),
395 cdt);
396 //rotate
397 fLocalDirection.rotateUz(dp->GetMomentumDirection());
398 //return
399 return fLocalDirection;
400 }
401
402 if (!(fLorentzTables1->count(Zmat)) || !(fLorentzTables2->count(Zmat)))
403 {
404 G4ExceptionDescription ed;
405 ed << "Unable to retrieve Lorentz tables for Z= " << Zmat << G4endl;
406 G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()",
407 "em2006",FatalException,ed);
408 }
409
410 //retrieve actual tables
411 const G4PhysicsTable* theTable1 = fLorentzTables1->find(Zmat)->second;
412 const G4PhysicsTable* theTable2 = fLorentzTables2->find(Zmat)->second;
413
414 G4double RK=20.0*eGamma/ePrimary;
415 G4int ik=std::min((G4int) RK,19);
416
417 G4double P10=0,P11=0,P1=0;
418 G4double P20=0,P21=0,P2=0;
419
420 //First coefficient
421 const G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTable1)[ik];
422 const G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTable1)[ik+1];
423 P10 = v1->Value(beta);
424 P11 = v2->Value(beta);
425 P1=P10+(RK-(G4double) ik)*(P11-P10);
426
427 //Second coefficient
428 const G4PhysicsFreeVector* v3 = (G4PhysicsFreeVector*) (*theTable2)[ik];
429 const G4PhysicsFreeVector* v4 = (G4PhysicsFreeVector*) (*theTable2)[ik+1];
430 P20=v3->Value(beta);
431 P21=v4->Value(beta);
432 P2=P20+(RK-(G4double) ik)*(P21-P20);
433
434 //Sampling from the Lorenz-trasformed dipole distributions
435 P1=std::min(G4Exp(P1)/beta,1.0);
436 G4double betap = std::min(std::max(beta*(1.0+P2/beta),0.0),0.9999);
437
438 G4double testf=0;
439
440 if (G4UniformRand() < P1)
441 {
442 do{
443 cdt = 2.0*G4UniformRand()-1.0;
444 testf=2.0*G4UniformRand()-(1.0+cdt*cdt);
445 }while(testf>0);
446 }
447 else
448 {
449 do{
450 cdt = 2.0*G4UniformRand()-1.0;
451 testf=G4UniformRand()-(1.0-cdt*cdt);
452 }while(testf>0);
453 }
454 cdt = (cdt+betap)/(1.0+betap*cdt);
455
456 //Get primary kinematics
457 sinTheta = std::sqrt(1. - cdt*cdt);
458 phi = twopi * G4UniformRand();
459 fLocalDirection.set(sinTheta* std::cos(phi),
460 sinTheta* std::sin(phi),
461 cdt);
462 //rotate
463 fLocalDirection.rotateUz(dp->GetMomentumDirection());
464 //return
465 return fLocalDirection;
466
467 }
468
469 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
470
471 G4double G4PenelopeBremsstrahlungAngular::CalculateEffectiveZ(const G4Material* material)
472 {
473 if (!fEffectiveZSq)
474 fEffectiveZSq = new std::map<const G4Material*,G4double>;
475
476 //Already exists: return it
477 if (fEffectiveZSq->count(material))
478 return fEffectiveZSq->find(material)->second;
479
480 //Helper for the calculation
481 std::vector<G4double> *StechiometricFactors = new std::vector<G4double>;
482 G4int nElements = (G4int)material->GetNumberOfElements();
483 const G4ElementVector* elementVector = material->GetElementVector();
484 const G4double* fractionVector = material->GetFractionVector();
485 for (G4int i=0;i<nElements;++i)
486 {
487 G4double fraction = fractionVector[i];
488 G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
489 StechiometricFactors->push_back(fraction/atomicWeigth);
490 }
491 //Find max
492 G4double MaxStechiometricFactor = 0.;
493 for (G4int i=0;i<nElements;++i)
494 {
495 if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
496 MaxStechiometricFactor = (*StechiometricFactors)[i];
497 }
498 //Normalize
499 for (G4int i=0;i<nElements;++i)
500 (*StechiometricFactors)[i] /= MaxStechiometricFactor;
501
502 G4double sumz2 = 0;
503 G4double sums = 0;
504 for (G4int i=0;i<nElements;++i)
505 {
506 G4double Z = (*elementVector)[i]->GetZ();
507 sumz2 += (*StechiometricFactors)[i]*Z*Z;
508 sums += (*StechiometricFactors)[i];
509 }
510 delete StechiometricFactors;
511
512 G4double ZBR = std::sqrt(sumz2/sums);
513 fEffectiveZSq->insert(std::make_pair(material,ZBR));
514
515 return ZBR;
516 }
517