Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PenelopeBremsstrahlungAngular.cc
Go to the documentation of this file.
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 // $Id: G4PenelopeBremsstrahlungAngular.cc 76220 2013-11-08 10:15:00Z gcosmo $
27 //
28 // --------------------------------------------------------------
29 //
30 // File name: G4PenelopeBremsstrahlungAngular
31 //
32 // Author: Luciano Pandola
33 //
34 // Creation date: November 2010
35 //
36 // History:
37 // -----------
38 // 23 Nov 2010 L. Pandola 1st implementation
39 // 24 May 2011 L. Pandola Renamed (make v2008 as default Penelope)
40 // 13 Mar 2012 L. Pandola Made a derived class of G4VEmAngularDistribution
41 // and update the interface accordingly
42 // 18 Jul 2012 L. Pandola Migrated to the new basic interface of G4VEmAngularDistribution
43 // Now returns a G4ThreeVector and takes care of the rotation
44 // 03 Oct 2013 L. Pandola Migrated to MT: only the master model handles tables
45 // 17 Oct 2013 L. Pandola Partially revert MT migration. The angular generator is kept as
46 // thread-local, and each worker has full access to it.
47 //
48 //----------------------------------------------------------------
49 
51 
52 #include "globals.hh"
53 #include "G4PhysicalConstants.hh"
54 #include "G4SystemOfUnits.hh"
55 #include "G4PhysicsFreeVector.hh"
56 #include "G4PhysicsTable.hh"
57 #include "G4Material.hh"
58 #include "Randomize.hh"
59 
61  G4VEmAngularDistribution("Penelope"), theEffectiveZSq(0),
62  theLorentzTables1(0),theLorentzTables2(0)
63 
64 {
65  dataRead = false;
66  verbosityLevel = 0;
67 }
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
70 
72 {
73  ClearTables();
74 }
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77 
79 {
80  ClearTables();
81 }
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84 
85 void G4PenelopeBremsstrahlungAngular::ClearTables()
86 {
87  std::map<G4double,G4PhysicsTable*>::iterator j;
88 
89  if (theLorentzTables1)
90  {
91  for (j=theLorentzTables1->begin(); j != theLorentzTables1->end(); j++)
92  {
93  G4PhysicsTable* tab = j->second;
94  //tab->clearAndDestroy();
95  delete tab;
96  }
97  delete theLorentzTables1;
98  theLorentzTables1 = 0;
99  }
100 
101  if (theLorentzTables2)
102  {
103  for (j=theLorentzTables2->begin(); j != theLorentzTables2->end(); j++)
104  {
105  G4PhysicsTable* tab = j->second;
106  //tab->clearAndDestroy();
107  delete tab;
108  }
109  delete theLorentzTables2;
110  theLorentzTables2 = 0;
111  }
112  if (theEffectiveZSq)
113  {
114  delete theEffectiveZSq;
115  theEffectiveZSq = 0;
116  }
117 }
118 
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120 
121 void G4PenelopeBremsstrahlungAngular::ReadDataFile()
122 {
123  //Read information from DataBase file
124  char* path = getenv("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<NumberofKPoints;k++)
147  for (i=0;i<NumberofZPoints;i++)
148  for (j=0;j<NumberofEPoints;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  QQ1[i][j][k]=a1;
158  QQ2[i][j][k]=a2;
159  }
160  else
161  {
163  ed << "Corrupted data file " << pathFile << "?" << G4endl;
164  G4Exception("G4PenelopeBremsstrahlungAngular::ReadDataFile()",
165  "em0005",FatalException,ed);
166  }
167  }
168  file.close();
169  dataRead = true;
170 }
171 
172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
173 
175 {
176  //Unused at the moment: the G4PenelopeBremsstrahlungAngular is thread-local, so each worker
177  //builds its own version of the tables.
178  /*
179  if (!isMaster)
180  //Should not be here!
181  G4Exception("G4PenelopeBremsstrahlungAngular::PrepareTables()",
182  "em0100",FatalException,"Worker thread in this method");
183  */
184 
185  //Check if data file has already been read
186  if (!dataRead)
187  {
188  ReadDataFile();
189  if (!dataRead)
190  G4Exception("G4PenelopeBremsstrahlungAngular::PrepareInterpolationTables()",
191  "em2001",FatalException,"Unable to build interpolation table");
192  }
193 
194  if (!theLorentzTables1)
195  theLorentzTables1 = new std::map<G4double,G4PhysicsTable*>;
196  if (!theLorentzTables2)
197  theLorentzTables2 = new std::map<G4double,G4PhysicsTable*>;
198 
199  G4double Zmat = CalculateEffectiveZ(material);
200 
201  const G4int reducedEnergyGrid=21;
202  //Support arrays.
203  G4double betas[NumberofEPoints]; //betas for interpolation
204  //tables for interpolation
205  G4double Q1[NumberofEPoints][NumberofKPoints];
206  G4double Q2[NumberofEPoints][NumberofKPoints];
207  //expanded tables for interpolation
208  G4double Q1E[NumberofEPoints][reducedEnergyGrid];
209  G4double Q2E[NumberofEPoints][reducedEnergyGrid];
210  G4double pZ[NumberofZPoints] = {2.0,8.0,13.0,47.0,79.0,92.0};
211 
212  G4int i=0,j=0,k=0; // i=index for Z, j=index for E, k=index for K
213  //Interpolation in Z
214  for (i=0;i<NumberofEPoints;i++)
215  {
216  for (j=0;j<NumberofKPoints;j++)
217  {
218  G4PhysicsFreeVector* QQ1vector = new G4PhysicsFreeVector(NumberofZPoints);
219  G4PhysicsFreeVector* QQ2vector = new G4PhysicsFreeVector(NumberofZPoints);
220 
221  //fill vectors
222  for (k=0;k<NumberofZPoints;k++)
223  {
224  QQ1vector->PutValue(k,pZ[k],std::log(QQ1[k][i][j]));
225  QQ2vector->PutValue(k,pZ[k],QQ2[k][i][j]);
226  }
227 
228  QQ1vector->SetSpline(true);
229  QQ2vector->SetSpline(true);
230 
231  Q1[i][j]= std::exp(QQ1vector->Value(Zmat));
232  Q2[i][j]=QQ2vector->Value(Zmat);
233  delete QQ1vector;
234  delete QQ2vector;
235  }
236  }
237  G4double pE[NumberofEPoints] = {1.0e-03*MeV,5.0e-03*MeV,1.0e-02*MeV,5.0e-02*MeV,
238  1.0e-01*MeV,5.0e-01*MeV};
239  G4double pK[NumberofKPoints] = {0.0,0.6,0.8,0.95};
240  G4double ppK[reducedEnergyGrid];
241 
242  for(i=0;i<reducedEnergyGrid;i++)
243  ppK[i]=((G4double) i) * 0.05;
244 
245 
246  for(i=0;i<NumberofEPoints;i++)
247  betas[i]=std::sqrt(pE[i]*(pE[i]+2*electron_mass_c2))/(pE[i]+electron_mass_c2);
248 
249 
250  for (i=0;i<NumberofEPoints;i++)
251  {
252  for (j=0;j<NumberofKPoints;j++)
253  Q1[i][j]=Q1[i][j]/Zmat;
254  }
255 
256  //Expanded table of distribution parameters
257  for (i=0;i<NumberofEPoints;i++)
258  {
259  G4PhysicsFreeVector* Q1vector = new G4PhysicsFreeVector(NumberofKPoints);
260  G4PhysicsFreeVector* Q2vector = new G4PhysicsFreeVector(NumberofKPoints);
261 
262  for (j=0;j<NumberofKPoints;j++)
263  {
264  Q1vector->PutValue(j,pK[j],std::log(Q1[i][j])); //logarithmic
265  Q2vector->PutValue(j,pK[j],Q2[i][j]);
266  }
267 
268  for (j=0;j<reducedEnergyGrid;j++)
269  {
270  Q1E[i][j]=Q1vector->Value(ppK[j]);
271  Q2E[i][j]=Q2vector->Value(ppK[j]);
272  }
273  delete Q1vector;
274  delete Q2vector;
275  }
276  //
277  //TABLES to be stored
278  //
279  G4PhysicsTable* theTable1 = new G4PhysicsTable();
280  G4PhysicsTable* theTable2 = new G4PhysicsTable();
281  // the table will contain reducedEnergyGrid G4PhysicsFreeVectors with different
282  // values of k,
283  // Each of the G4PhysicsFreeVectors has a profile of
284  // y vs. E
285  //
286  //reserve space of the vectors.
287  for (j=0;j<reducedEnergyGrid;j++)
288  {
289  G4PhysicsFreeVector* thevec = new G4PhysicsFreeVector(NumberofEPoints);
290  theTable1->push_back(thevec);
291  G4PhysicsFreeVector* thevec2 = new G4PhysicsFreeVector(NumberofEPoints);
292  theTable2->push_back(thevec2);
293  }
294 
295  for (j=0;j<reducedEnergyGrid;j++)
296  {
297  G4PhysicsFreeVector* thevec = (G4PhysicsFreeVector*) (*theTable1)[j];
298  G4PhysicsFreeVector* thevec2 = (G4PhysicsFreeVector*) (*theTable2)[j];
299  for (i=0;i<NumberofEPoints;i++)
300  {
301  thevec->PutValue(i,betas[i],Q1E[i][j]);
302  thevec2->PutValue(i,betas[i],Q2E[i][j]);
303  }
304  thevec->SetSpline(true);
305  thevec2->SetSpline(true);
306  }
307 
308  if (theLorentzTables1 && theLorentzTables2)
309  {
310  theLorentzTables1->insert(std::make_pair(Zmat,theTable1));
311  theLorentzTables2->insert(std::make_pair(Zmat,theTable2));
312  }
313  else
314  {
316  ed << "Unable to create tables of Lorentz coefficients for " << G4endl;
317  ed << "<Z>= " << Zmat << " in G4PenelopeBremsstrahlungAngular" << G4endl;
318  delete theTable1;
319  delete theTable2;
320  G4Exception("G4PenelopeBremsstrahlungAngular::PrepareInterpolationTables()",
321  "em2005",FatalException,ed);
322  }
323 
324  return;
325 }
326 
327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
328 
330  G4double eGamma,
331  G4int,
332  const G4Material* material)
333 {
334  if (!material)
335  {
336  G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()",
337  "em2040",FatalException,"The pointer to G4Material* is NULL");
338  return fLocalDirection;
339  }
340 
341  //Retrieve the effective Z
342  G4double Zmat = 0;
343 
344  if (!theEffectiveZSq)
345  {
346  G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()",
347  "em2040",FatalException,"EffectiveZ table not available");
348  return fLocalDirection;
349  }
350 
351  //found in the table: return it
352  if (theEffectiveZSq->count(material))
353  Zmat = theEffectiveZSq->find(material)->second;
354  else
355  {
356  G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()",
357  "em2040",FatalException,"Material not found in the effectiveZ table");
358  return fLocalDirection;
359  }
360 
361  if (verbosityLevel > 0)
362  {
363  G4cout << "Effective <Z> for material : " << material->GetName() <<
364  " = " << Zmat << G4endl;
365  }
366 
367  G4double ePrimary = dp->GetKineticEnergy();
368 
369  G4double beta = std::sqrt(ePrimary*(ePrimary+2*electron_mass_c2))/
370  (ePrimary+electron_mass_c2);
371  G4double cdt = 0;
372  G4double sinTheta = 0;
373  G4double phi = 0;
374 
375  //Use a pure dipole distribution for energy above 500 keV
376  if (ePrimary > 500*keV)
377  {
378  cdt = 2.0*G4UniformRand() - 1.0;
379  if (G4UniformRand() > 0.75)
380  {
381  if (cdt<0)
382  cdt = -1.0*std::pow(-cdt,1./3.);
383  else
384  cdt = std::pow(cdt,1./3.);
385  }
386  cdt = (cdt+beta)/(1.0+beta*cdt);
387  //Get primary kinematics
388  sinTheta = std::sqrt(1. - cdt*cdt);
389  phi = twopi * G4UniformRand();
390  fLocalDirection.set(sinTheta* std::cos(phi),
391  sinTheta* std::sin(phi),
392  cdt);
393  //rotate
395  //return
396  return fLocalDirection;
397  }
398 
399  if (!(theLorentzTables1->count(Zmat)) || !(theLorentzTables2->count(Zmat)))
400  {
402  ed << "Unable to retrieve Lorentz tables for Z= " << Zmat << G4endl;
403  G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()",
404  "em2006",FatalException,ed);
405  }
406 
407  //retrieve actual tables
408  const G4PhysicsTable* theTable1 = theLorentzTables1->find(Zmat)->second;
409  const G4PhysicsTable* theTable2 = theLorentzTables2->find(Zmat)->second;
410 
411  G4double RK=20.0*eGamma/ePrimary;
412  G4int ik=std::min((G4int) RK,19);
413 
414  G4double P10=0,P11=0,P1=0;
415  G4double P20=0,P21=0,P2=0;
416 
417  //First coefficient
418  const G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTable1)[ik];
419  const G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTable1)[ik+1];
420  P10 = v1->Value(beta);
421  P11 = v2->Value(beta);
422  P1=P10+(RK-(G4double) ik)*(P11-P10);
423 
424  //Second coefficient
425  const G4PhysicsFreeVector* v3 = (G4PhysicsFreeVector*) (*theTable2)[ik];
426  const G4PhysicsFreeVector* v4 = (G4PhysicsFreeVector*) (*theTable2)[ik+1];
427  P20=v3->Value(beta);
428  P21=v4->Value(beta);
429  P2=P20+(RK-(G4double) ik)*(P21-P20);
430 
431  //Sampling from the Lorenz-trasformed dipole distributions
432  P1=std::min(std::exp(P1)/beta,1.0);
433  G4double betap = std::min(std::max(beta*(1.0+P2/beta),0.0),0.9999);
434 
435  G4double testf=0;
436 
437  if (G4UniformRand() < P1)
438  {
439  do{
440  cdt = 2.0*G4UniformRand()-1.0;
441  testf=2.0*G4UniformRand()-(1.0+cdt*cdt);
442  }while(testf>0);
443  }
444  else
445  {
446  do{
447  cdt = 2.0*G4UniformRand()-1.0;
448  testf=G4UniformRand()-(1.0-cdt*cdt);
449  }while(testf>0);
450  }
451  cdt = (cdt+betap)/(1.0+betap*cdt);
452 
453  //Get primary kinematics
454  sinTheta = std::sqrt(1. - cdt*cdt);
455  phi = twopi * G4UniformRand();
456  fLocalDirection.set(sinTheta* std::cos(phi),
457  sinTheta* std::sin(phi),
458  cdt);
459  //rotate
461  //return
462  return fLocalDirection;
463 
464 }
465 
466 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
467 
469  const G4double ,
470  const G4int )
471 {
472  G4cout << "WARNING: G4PenelopeBremsstrahlungAngular() does NOT support PolarAngle()" << G4endl;
473  G4cout << "Please use the alternative interface SampleDirection()" << G4endl;
474  G4Exception("G4PenelopeBremsstrahlungAngular::PolarAngle()",
475  "em0005",FatalException,"Unsupported interface");
476  return 0;
477 }
478 
479 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
480 
481 G4double G4PenelopeBremsstrahlungAngular::CalculateEffectiveZ(const G4Material* material)
482 {
483  if (!theEffectiveZSq)
484  theEffectiveZSq = new std::map<const G4Material*,G4double>;
485 
486  //Already exists: return it
487  if (theEffectiveZSq->count(material))
488  return theEffectiveZSq->find(material)->second;
489 
490  //Helper for the calculation
491  std::vector<G4double> *StechiometricFactors = new std::vector<G4double>;
492  G4int nElements = material->GetNumberOfElements();
493  const G4ElementVector* elementVector = material->GetElementVector();
494  const G4double* fractionVector = material->GetFractionVector();
495  for (G4int i=0;i<nElements;i++)
496  {
497  G4double fraction = fractionVector[i];
498  G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
499  StechiometricFactors->push_back(fraction/atomicWeigth);
500  }
501  //Find max
502  G4double MaxStechiometricFactor = 0.;
503  for (G4int i=0;i<nElements;i++)
504  {
505  if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
506  MaxStechiometricFactor = (*StechiometricFactors)[i];
507  }
508  //Normalize
509  for (G4int i=0;i<nElements;i++)
510  (*StechiometricFactors)[i] /= MaxStechiometricFactor;
511 
512  G4double sumz2 = 0;
513  G4double sums = 0;
514  for (G4int i=0;i<nElements;i++)
515  {
516  G4double Z = (*elementVector)[i]->GetZ();
517  sumz2 += (*StechiometricFactors)[i]*Z*Z;
518  sums += (*StechiometricFactors)[i];
519  }
520  delete StechiometricFactors;
521 
522  G4double ZBR = std::sqrt(sumz2/sums);
523  theEffectiveZSq->insert(std::make_pair(material,ZBR));
524 
525  return ZBR;
526 }
void set(double x, double y, double z)
std::vector< G4Element * > G4ElementVector
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
const G4String & GetName() const
Definition: G4Material.hh:176
void push_back(G4PhysicsVector *)
void PrepareTables(const G4Material *material, G4bool isMaster)
Reserved for Master Model.
G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double out_energy, G4int Z, const G4Material *mat=0)
Samples the direction of the outgoing photon (in global coordinates).
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:188
int G4int
Definition: G4Types.hh:78
void SetSpline(G4bool)
string material
Definition: eplot.py:19
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
std::ostream & tab(std::ostream &)
Definition: CCalutils.cc:89
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
float electron_mass_c2
Definition: hepunit.py:274
G4double Value(G4double theEnergy, size_t &lastidx) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
#define G4endl
Definition: G4ios.hh:61
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
double G4double
Definition: G4Types.hh:76
G4double PolarAngle(const G4double initial_energy, const G4double final_energy, const G4int Z)
const G4double * GetFractionVector() const
Definition: G4Material.hh:192