Geant4  9.6.p02
 All Classes 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$
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 //
45 //----------------------------------------------------------------
46 
47 #include "G4PenelopeBremsstrahlungAngular.hh"
48 
49 #include "globals.hh"
50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4PhysicsFreeVector.hh"
53 #include "G4PhysicsTable.hh"
54 #include "G4Material.hh"
55 #include "Randomize.hh"
56 
57 G4PenelopeBremsstrahlungAngular::G4PenelopeBremsstrahlungAngular() :
58  G4VEmAngularDistribution("Penelope"), theEffectiveZSq(0),
59  theLorentzTables1(0),theLorentzTables2(0)
60 
61 {
62  dataRead = false;
63  verbosityLevel = 0;
64 }
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67 
68 G4PenelopeBremsstrahlungAngular::~G4PenelopeBremsstrahlungAngular()
69 {
70  ClearTables();
71 }
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74 
75 void G4PenelopeBremsstrahlungAngular::Initialize()
76 {
77  ClearTables();
78 }
79 
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 
82 void G4PenelopeBremsstrahlungAngular::ClearTables()
83 {
84  std::map<G4double,G4PhysicsTable*>::iterator j;
85 
86  if (theLorentzTables1)
87  {
88  for (j=theLorentzTables1->begin(); j != theLorentzTables1->end(); j++)
89  {
90  G4PhysicsTable* tab = j->second;
91  tab->clearAndDestroy();
92  delete tab;
93  }
94  delete theLorentzTables1;
95  theLorentzTables1 = 0;
96  }
97 
98  if (theLorentzTables2)
99  {
100  for (j=theLorentzTables2->begin(); j != theLorentzTables2->end(); j++)
101  {
102  G4PhysicsTable* tab = j->second;
103  tab->clearAndDestroy();
104  delete tab;
105  }
106  delete theLorentzTables2;
107  theLorentzTables2 = 0;
108  }
109  if (theEffectiveZSq)
110  {
111  delete theEffectiveZSq;
112  theEffectiveZSq = 0;
113  }
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
117 
118 void G4PenelopeBremsstrahlungAngular::ReadDataFile()
119 {
120  //Read information from DataBase file
121  char* path = getenv("G4LEDATA");
122  if (!path)
123  {
124  G4String excep =
125  "G4PenelopeBremsstrahlungAngular - G4LEDATA environment variable not set!";
126  G4Exception("G4PenelopeBremsstrahlungAngular::ReadDataFile()",
127  "em0006",FatalException,excep);
128  return;
129  }
130  G4String pathString(path);
131  G4String pathFile = pathString + "/penelope/bremsstrahlung/pdbrang.p08";
132  std::ifstream file(pathFile);
133 
134  if (!file.is_open())
135  {
136  G4String excep = "G4PenelopeBremsstrahlungAngular - data file " + pathFile + " not found!";
137  G4Exception("G4PenelopeBremsstrahlungAngular::ReadDataFile()",
138  "em0003",FatalException,excep);
139  return;
140  }
141  G4int i=0,j=0,k=0; // i=index for Z, j=index for E, k=index for K
142 
143  for (k=0;k<NumberofKPoints;k++)
144  for (i=0;i<NumberofZPoints;i++)
145  for (j=0;j<NumberofEPoints;j++)
146  {
147  G4double a1,a2;
148  G4int ik1,iz1,ie1;
149  G4double zr,er,kr;
150  file >> iz1 >> ie1 >> ik1 >> zr >> er >> kr >> a1 >> a2;
151  //check the data are correct
152  if ((iz1-1 == i) && (ik1-1 == k) && (ie1-1 == j))
153  {
154  QQ1[i][j][k]=a1;
155  QQ2[i][j][k]=a2;
156  }
157  else
158  {
159  G4ExceptionDescription ed;
160  ed << "Corrupted data file " << pathFile << "?" << G4endl;
161  G4Exception("G4PenelopeBremsstrahlungAngular::ReadDataFile()",
162  "em0005",FatalException,ed);
163  }
164  }
165  file.close();
166  dataRead = true;
167 }
168 
169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
170 
171 void G4PenelopeBremsstrahlungAngular::PrepareInterpolationTables(G4double Zmat)
172 {
173  //Check if data file has already been read
174  if (!dataRead)
175  {
176  ReadDataFile();
177  if (!dataRead)
178  G4Exception("G4PenelopeBremsstrahlungAngular::PrepareInterpolationTables()",
179  "em2001",FatalException,"Unable to build interpolation table");
180  }
181 
182  const G4int reducedEnergyGrid=21;
183  //Support arrays.
184  G4double betas[NumberofEPoints]; //betas for interpolation
185  //tables for interpolation
186  G4double Q1[NumberofEPoints][NumberofKPoints];
187  G4double Q2[NumberofEPoints][NumberofKPoints];
188  //expanded tables for interpolation
189  G4double Q1E[NumberofEPoints][reducedEnergyGrid];
190  G4double Q2E[NumberofEPoints][reducedEnergyGrid];
191  G4double pZ[NumberofZPoints] = {2.0,8.0,13.0,47.0,79.0,92.0};
192 
193  G4int i=0,j=0,k=0; // i=index for Z, j=index for E, k=index for K
194  //Interpolation in Z
195  for (i=0;i<NumberofEPoints;i++)
196  {
197  for (j=0;j<NumberofKPoints;j++)
198  {
199  G4PhysicsFreeVector* QQ1vector = new G4PhysicsFreeVector(NumberofZPoints);
200  QQ1vector->SetSpline(true);
201  G4PhysicsFreeVector* QQ2vector = new G4PhysicsFreeVector(NumberofZPoints);
202  QQ2vector->SetSpline(true);
203 
204  //fill vectors
205  for (k=0;k<NumberofZPoints;k++)
206  {
207  QQ1vector->PutValue(k,pZ[k],std::log(QQ1[k][i][j]));
208  QQ2vector->PutValue(k,pZ[k],QQ2[k][i][j]);
209  }
210 
211  Q1[i][j]= std::exp(QQ1vector->Value(Zmat));
212  Q2[i][j]=QQ2vector->Value(Zmat);
213  delete QQ1vector;
214  delete QQ2vector;
215  }
216  }
217  G4double pE[NumberofEPoints] = {1.0e-03*MeV,5.0e-03*MeV,1.0e-02*MeV,5.0e-02*MeV,
218  1.0e-01*MeV,5.0e-01*MeV};
219  G4double pK[NumberofKPoints] = {0.0,0.6,0.8,0.95};
220  G4double ppK[reducedEnergyGrid];
221 
222  for(i=0;i<reducedEnergyGrid;i++)
223  ppK[i]=((G4double) i) * 0.05;
224 
225 
226  for(i=0;i<NumberofEPoints;i++)
227  betas[i]=std::sqrt(pE[i]*(pE[i]+2*electron_mass_c2))/(pE[i]+electron_mass_c2);
228 
229 
230  for (i=0;i<NumberofEPoints;i++)
231  {
232  for (j=0;j<NumberofKPoints;j++)
233  Q1[i][j]=Q1[i][j]/Zmat;
234  }
235 
236  //Expanded table of distribution parameters
237  for (i=0;i<NumberofEPoints;i++)
238  {
239  G4PhysicsFreeVector* Q1vector = new G4PhysicsFreeVector(NumberofKPoints);
240  G4PhysicsFreeVector* Q2vector = new G4PhysicsFreeVector(NumberofKPoints);
241 
242  for (j=0;j<NumberofKPoints;j++)
243  {
244  Q1vector->PutValue(j,pK[j],std::log(Q1[i][j])); //logarithmic
245  Q2vector->PutValue(j,pK[j],Q2[i][j]);
246  }
247 
248  for (j=0;j<reducedEnergyGrid;j++)
249  {
250  Q1E[i][j]=Q1vector->Value(ppK[j]);
251  Q2E[i][j]=Q2vector->Value(ppK[j]);
252  }
253  delete Q1vector;
254  delete Q2vector;
255  }
256  //
257  //TABLES to be stored
258  //
259  G4PhysicsTable* theTable1 = new G4PhysicsTable();
260  G4PhysicsTable* theTable2 = new G4PhysicsTable();
261  // the table will contain reducedEnergyGrid G4PhysicsFreeVectors with different
262  // values of k,
263  // Each of the G4PhysicsFreeVectors has a profile of
264  // y vs. E
265  //
266  //reserve space of the vectors.
267  for (j=0;j<reducedEnergyGrid;j++)
268  {
269  G4PhysicsFreeVector* thevec = new G4PhysicsFreeVector(NumberofEPoints);
270  thevec->SetSpline(true);
271  theTable1->push_back(thevec);
272  G4PhysicsFreeVector* thevec2 = new G4PhysicsFreeVector(NumberofEPoints);
273  thevec2->SetSpline(true);
274  theTable2->push_back(thevec2);
275  }
276 
277  for (j=0;j<reducedEnergyGrid;j++)
278  {
279  G4PhysicsFreeVector* thevec = (G4PhysicsFreeVector*) (*theTable1)[j];
280  G4PhysicsFreeVector* thevec2 = (G4PhysicsFreeVector*) (*theTable2)[j];
281  for (i=0;i<NumberofEPoints;i++)
282  {
283  thevec->PutValue(i,betas[i],Q1E[i][j]);
284  thevec2->PutValue(i,betas[i],Q2E[i][j]);
285  }
286  }
287 
288  if (theLorentzTables1 && theLorentzTables2)
289  {
290  theLorentzTables1->insert(std::make_pair(Zmat,theTable1));
291  theLorentzTables2->insert(std::make_pair(Zmat,theTable2));
292  }
293  else
294  {
295  G4ExceptionDescription ed;
296  ed << "Unable to create tables of Lorentz coefficients for " << G4endl;
297  ed << "<Z>= " << Zmat << " in G4PenelopeBremsstrahlungAngular" << G4endl;
298  delete theTable1;
299  delete theTable2;
300  G4Exception("G4PenelopeBremsstrahlungAngular::PrepareInterpolationTables()",
301  "em2005",FatalException,ed);
302  }
303 
304  return;
305 }
306 
307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
308 
309 G4ThreeVector& G4PenelopeBremsstrahlungAngular::SampleDirection(const G4DynamicParticle* dp,
310  G4double eGamma,
311  G4int,
312  const G4Material* material)
313 {
314  if (!material)
315  {
316  G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()",
317  "em2040",FatalException,"The pointer to G4Material* is NULL");
318  return fLocalDirection;
319  }
320 
321  G4double Zmat = GetEffectiveZ(material);
322  if (verbosityLevel > 0)
323  {
324  G4cout << "Effective <Z> for material : " << material->GetName() <<
325  " = " << Zmat << G4endl;
326  }
327 
328  G4double ePrimary = dp->GetKineticEnergy();
329 
330  G4double beta = std::sqrt(ePrimary*(ePrimary+2*electron_mass_c2))/
331  (ePrimary+electron_mass_c2);
332  G4double cdt = 0;
333  G4double sinTheta = 0;
334  G4double phi = 0;
335 
336  //Use a pure dipole distribution for energy above 500 keV
337  if (ePrimary > 500*keV)
338  {
339  cdt = 2.0*G4UniformRand() - 1.0;
340  if (G4UniformRand() > 0.75)
341  {
342  if (cdt<0)
343  cdt = -1.0*std::pow(-cdt,1./3.);
344  else
345  cdt = std::pow(cdt,1./3.);
346  }
347  cdt = (cdt+beta)/(1.0+beta*cdt);
348  //Get primary kinematics
349  sinTheta = std::sqrt(1. - cdt*cdt);
350  phi = twopi * G4UniformRand();
351  fLocalDirection.set(sinTheta* std::cos(phi),
352  sinTheta* std::sin(phi),
353  cdt);
354  //rotate
355  fLocalDirection.rotateUz(dp->GetMomentumDirection());
356  //return
357  return fLocalDirection;
358  }
359 
360  //Else, retrieve tables and go through the full thing
361  if (!theLorentzTables1)
362  theLorentzTables1 = new std::map<G4double,G4PhysicsTable*>;
363  if (!theLorentzTables2)
364  theLorentzTables2 = new std::map<G4double,G4PhysicsTable*>;
365 
366  //Check if tables exist for the given Zmat
367  if (!(theLorentzTables1->count(Zmat)))
368  PrepareInterpolationTables(Zmat);
369 
370  if (!(theLorentzTables1->count(Zmat)) || !(theLorentzTables2->count(Zmat)))
371  {
372  G4ExceptionDescription ed;
373  ed << "Unable to retrieve Lorentz tables for Z= " << Zmat << G4endl;
374  G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()",
375  "em2006",FatalException,ed);
376  }
377 
378  //retrieve actual tables
379  G4PhysicsTable* theTable1 = theLorentzTables1->find(Zmat)->second;
380  G4PhysicsTable* theTable2 = theLorentzTables2->find(Zmat)->second;
381 
382  G4double RK=20.0*eGamma/ePrimary;
383  G4int ik=std::min((G4int) RK,19);
384 
385  G4double P10=0,P11=0,P1=0;
386  G4double P20=0,P21=0,P2=0;
387 
388  //First coefficient
389  G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTable1)[ik];
390  G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTable1)[ik+1];
391  P10 = v1->Value(beta);
392  P11 = v2->Value(beta);
393  P1=P10+(RK-(G4double) ik)*(P11-P10);
394 
395  //Second coefficient
396  G4PhysicsFreeVector* v3 = (G4PhysicsFreeVector*) (*theTable2)[ik];
397  G4PhysicsFreeVector* v4 = (G4PhysicsFreeVector*) (*theTable2)[ik+1];
398  P20=v3->Value(beta);
399  P21=v4->Value(beta);
400  P2=P20+(RK-(G4double) ik)*(P21-P20);
401 
402  //Sampling from the Lorenz-trasformed dipole distributions
403  P1=std::min(std::exp(P1)/beta,1.0);
404  G4double betap = std::min(std::max(beta*(1.0+P2/beta),0.0),0.9999);
405 
406  G4double testf=0;
407 
408  if (G4UniformRand() < P1)
409  {
410  do{
411  cdt = 2.0*G4UniformRand()-1.0;
412  testf=2.0*G4UniformRand()-(1.0+cdt*cdt);
413  }while(testf>0);
414  }
415  else
416  {
417  do{
418  cdt = 2.0*G4UniformRand()-1.0;
419  testf=G4UniformRand()-(1.0-cdt*cdt);
420  }while(testf>0);
421  }
422  cdt = (cdt+betap)/(1.0+betap*cdt);
423 
424  //Get primary kinematics
425  sinTheta = std::sqrt(1. - cdt*cdt);
426  phi = twopi * G4UniformRand();
427  fLocalDirection.set(sinTheta* std::cos(phi),
428  sinTheta* std::sin(phi),
429  cdt);
430  //rotate
431  fLocalDirection.rotateUz(dp->GetMomentumDirection());
432  //return
433  return fLocalDirection;
434 
435 }
436 
437 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
438 
439 G4double G4PenelopeBremsstrahlungAngular::PolarAngle(const G4double ,
440  const G4double ,
441  const G4int )
442 {
443  G4cout << "WARNING: G4PenelopeBremsstrahlungAngular() does NOT support PolarAngle()" << G4endl;
444  G4cout << "Please use the alternative interface SampleDirection()" << G4endl;
445  G4Exception("G4PenelopeBremsstrahlungAngular::PolarAngle()",
446  "em0005",FatalException,"Unsupported interface");
447  return 0;
448 }
449 
450 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
451 
452 G4double G4PenelopeBremsstrahlungAngular::GetEffectiveZ(const G4Material* material)
453 {
454  if (!theEffectiveZSq)
455  theEffectiveZSq = new std::map<const G4Material*,G4double>;
456 
457  //found in the table: return it
458  if (theEffectiveZSq->count(material))
459  return theEffectiveZSq->find(material)->second;
460 
461  //not found: calculate and return
462  //Helper for the calculation
463  std::vector<G4double> *StechiometricFactors = new std::vector<G4double>;
464  G4int nElements = material->GetNumberOfElements();
465  const G4ElementVector* elementVector = material->GetElementVector();
466  const G4double* fractionVector = material->GetFractionVector();
467  for (G4int i=0;i<nElements;i++)
468  {
469  G4double fraction = fractionVector[i];
470  G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
471  StechiometricFactors->push_back(fraction/atomicWeigth);
472  }
473  //Find max
474  G4double MaxStechiometricFactor = 0.;
475  for (G4int i=0;i<nElements;i++)
476  {
477  if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
478  MaxStechiometricFactor = (*StechiometricFactors)[i];
479  }
480  //Normalize
481  for (G4int i=0;i<nElements;i++)
482  (*StechiometricFactors)[i] /= MaxStechiometricFactor;
483 
484  G4double sumz2 = 0;
485  G4double sums = 0;
486  for (G4int i=0;i<nElements;i++)
487  {
488  G4double Z = (*elementVector)[i]->GetZ();
489  sumz2 += (*StechiometricFactors)[i]*Z*Z;
490  sums += (*StechiometricFactors)[i];
491  }
492  delete StechiometricFactors;
493 
494  G4double ZBR = std::sqrt(sumz2/sums);
495  theEffectiveZSq->insert(std::make_pair(material,ZBR));
496 
497  return ZBR;
498 }
Generated on Sat May 25 2013 14:33:34 for Geant4 by   doxygen 1.8.4