Geant4  10.02.p02
G4PenelopeBremsstrahlungFS.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: G4PenelopeBremsstrahlungFS.cc 97613 2016-06-06 12:24:51Z gcosmo $
27 //
28 // Author: Luciano Pandola
29 //
30 // History:
31 // --------
32 // 23 Nov 2010 L Pandola First complete implementation
33 // 02 May 2011 L.Pandola Remove dependency on CLHEP::HepMatrix
34 // 24 May 2011 L.Pandola Renamed (make v2008 as default Penelope)
35 // 03 Oct 2013 L.Pandola Migration to MT
36 // 30 Oct 2013 L.Pandola Use G4Cache to avoid new/delete of the
37 // data vector on the fly in SampleGammaEnergy()
38 //
39 
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 #include "G4PhysicsFreeVector.hh"
43 #include "G4PhysicsTable.hh"
44 #include "G4Material.hh"
45 #include "Randomize.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4Exp.hh"
49 
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 
53  theReducedXSTable(0),theEffectiveZSq(0),theSamplingTable(0),
54  thePBcut(0),fVerbosity(verbosity)
55 {
56  fCache.Put(0);
57  G4double tempvector[nBinsX] =
58  {1.0e-12,0.025e0,0.05e0,0.075e0,0.1e0,0.15e0,0.2e0,0.25e0,
59  0.3e0,0.35e0,0.4e0,0.45e0,0.5e0,0.55e0,0.6e0,0.65e0,0.7e0,
60  0.75e0,0.8e0,0.85e0,0.9e0,0.925e0,0.95e0,0.97e0,0.99e0,
61  0.995e0,0.999e0,0.9995e0,0.9999e0,0.99995e0,0.99999e0,1.0e0};
62 
63  for (size_t ix=0;ix<nBinsX;ix++)
64  theXGrid[ix] = tempvector[ix];
65 
66  for (size_t i=0;i<nBinsE;i++)
67  theEGrid[i] = 0.;
68 
69  theElementData = new std::map<G4int,G4DataVector*>;
70 }
71 
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73 
75 {
76  ClearTables();
77 
78  //The G4Physics*Vector pointers contained in the fCache are automatically deleted by
79  //the G4Allocator, so there is no need to take care of them manually
80 
81  //Clear manually theElementData
82  if (theElementData)
83  {
84  std::map<G4int,G4DataVector*>::iterator i;
85  for (i=theElementData->begin(); i != theElementData->end(); i++)
86  delete i->second;
87  delete theElementData;
88  theElementData = 0;
89  }
90 
91 }
92 
93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
94 
95 
97 {
98  //Just to check
99  if (!isMaster)
100  G4Exception("G4PenelopeBremsstrahlungFS::ClearTables()",
101  "em0100",FatalException,"Worker thread in this method");
102 
103  std::map< std::pair<const G4Material*,G4double> ,G4PhysicsTable*>::iterator j;
104 
105  if (theReducedXSTable)
106  {
107  for (j=theReducedXSTable->begin(); j != theReducedXSTable->end(); j++)
108  {
109  G4PhysicsTable* tab = j->second;
110  //tab->clearAndDestroy();
111  delete tab;
112  }
113  delete theReducedXSTable;
114  theReducedXSTable = 0;
115  }
116 
117  if (theSamplingTable)
118  {
119  for (j=theSamplingTable->begin(); j != theSamplingTable->end(); j++)
120  {
121  G4PhysicsTable* tab = j->second;
122  // tab->clearAndDestroy();
123  delete tab;
124  }
125  delete theSamplingTable;
126  theSamplingTable = 0;
127  }
128  if (thePBcut)
129  {
130  /*
131  std::map< std::pair<const G4Material*,G4double> ,G4PhysicsFreeVector*>::iterator kk;
132  for (kk=thePBcut->begin(); kk != thePBcut->end(); kk++)
133  delete kk->second;
134  */
135  delete thePBcut;
136  thePBcut = 0;
137  }
138 
139 
140  if (theEffectiveZSq)
141  {
142  delete theEffectiveZSq;
143  theEffectiveZSq = 0;
144  }
145 
146  return;
147 }
148 
149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150 
152 {
153  if (!theEffectiveZSq)
154  {
156  ed << "The container for the <Z^2> values is not initialized" << G4endl;
157  G4Exception("G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()",
158  "em2007",FatalException,ed);
159  return 0;
160  }
161  //found in the table: return it
162  if (theEffectiveZSq->count(material))
163  return theEffectiveZSq->find(material)->second;
164  else
165  {
167  ed << "The value of <Z^2> is not properly set for material " <<
168  material->GetName() << G4endl;
169  //requires running of BuildScaledXSTable()
170  G4Exception("G4PenelopeBremsstrahlungFS::GetEffectiveZSquared()",
171  "em2008",FatalException,ed);
172  }
173  return 0;
174 }
175 
176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
177 
179  G4double cut,G4bool isMaster)
180 {
181  //Corresponds to subroutines EBRaW and EBRaR of PENELOPE
182  /*
183  This method generates the table of the scaled energy-loss cross section from
184  bremsstrahlung emission for the given material. Original data are read from
185  file. The table is normalized according to the Berger-Seltzer cross section.
186  */
187 
188  //Just to check
189  if (!isMaster)
190  G4Exception("G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
191  "em0100",FatalException,"Worker thread in this method");
192 
193  if (fVerbosity > 2)
194  {
195  G4cout << "Entering in G4PenelopeBremsstrahlungFS::BuildScaledXSTable for " <<
196  material->GetName() << G4endl;
197  G4cout << "Threshold = " << cut/keV << " keV, isMaster= " << isMaster <<
198  G4endl;
199  }
200 
201  //This method should be accessed by the master only
202  if (!theSamplingTable)
204  new std::map< std::pair<const G4Material*,G4double> , G4PhysicsTable*>;
205  if (!thePBcut)
206  thePBcut =
207  new std::map< std::pair<const G4Material*,G4double> , G4PhysicsFreeVector* >;
208 
209  //check if the container exists (if not, create it)
210  if (!theReducedXSTable)
211  theReducedXSTable = new std::map< std::pair<const G4Material*,G4double> ,
212  G4PhysicsTable*>;
213  if (!theEffectiveZSq)
214  theEffectiveZSq = new std::map<const G4Material*,G4double>;
215 
216 
217 
218  //*********************************************************************
219  //Determine the equivalent atomic number <Z^2>
220  //*********************************************************************
221  std::vector<G4double> *StechiometricFactors = new std::vector<G4double>;
222  G4int nElements = material->GetNumberOfElements();
223  const G4ElementVector* elementVector = material->GetElementVector();
224  const G4double* fractionVector = material->GetFractionVector();
225  for (G4int i=0;i<nElements;i++)
226  {
227  G4double fraction = fractionVector[i];
228  G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
229  StechiometricFactors->push_back(fraction/atomicWeigth);
230  }
231  //Find max
232  G4double MaxStechiometricFactor = 0.;
233  for (G4int i=0;i<nElements;i++)
234  {
235  if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
236  MaxStechiometricFactor = (*StechiometricFactors)[i];
237  }
238  //Normalize
239  for (G4int i=0;i<nElements;i++)
240  (*StechiometricFactors)[i] /= MaxStechiometricFactor;
241 
242  G4double sumz2 = 0;
243  G4double sums = 0;
244  for (G4int i=0;i<nElements;i++)
245  {
246  G4double Z = (*elementVector)[i]->GetZ();
247  sumz2 += (*StechiometricFactors)[i]*Z*Z;
248  sums += (*StechiometricFactors)[i];
249  }
250  G4double ZBR2 = sumz2/sums;
251 
252  theEffectiveZSq->insert(std::make_pair(material,ZBR2));
253 
254 
255  //*********************************************************************
256  // loop on elements and read data files
257  //*********************************************************************
258  G4DataVector* tempData = new G4DataVector(nBinsE);
259  G4DataVector* tempMatrix = new G4DataVector(nBinsE*nBinsX,0.);
260 
261  for (G4int iel=0;iel<nElements;iel++)
262  {
263  G4double Z = (*elementVector)[iel]->GetZ();
264  G4int iZ = (G4int) Z;
265  G4double wgt = (*StechiometricFactors)[iel]*Z*Z/ZBR2;
266  //
267 
268  //the element is not already loaded
269  if (!theElementData->count(iZ))
270  {
271  ReadDataFile(iZ);
272  if (!theElementData->count(iZ))
273  {
275  ed << "Error in G4PenelopeBremsstrahlungFS::BuildScaledXSTable" << G4endl;
276  ed << "Unable to retrieve data for element " << iZ << G4endl;
277  G4Exception("G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
278  "em2009",FatalException,ed);
279  }
280  }
281 
282  G4DataVector* atomData = theElementData->find(iZ)->second;
283 
284  for (size_t ie=0;ie<nBinsE;ie++)
285  {
286  (*tempData)[ie] += wgt*(*atomData)[ie*(nBinsX+1)+nBinsX]; //last column contains total XS
287  for (size_t ix=0;ix<nBinsX;ix++)
288  (*tempMatrix)[ie*nBinsX+ix] += wgt*(*atomData)[ie*(nBinsX+1)+ix];
289  }
290  }
291 
292  //*********************************************************************
293  // the total energy loss spectrum is re-normalized to reproduce the total
294  // scaled cross section of Berger and Seltzer
295  //*********************************************************************
296  for (size_t ie=0;ie<nBinsE;ie++)
297  {
298  //for each energy, calculate integral of dSigma/dx over dx
299  G4double* tempData2 = new G4double[nBinsX];
300  for (size_t ix=0;ix<nBinsX;ix++)
301  tempData2[ix] = (*tempMatrix)[ie*nBinsX+ix];
302  G4double rsum = GetMomentumIntegral(tempData2,1.0,0);
303  delete[] tempData2;
304  G4double fact = millibarn*(theEGrid[ie]+electron_mass_c2)*(1./fine_structure_const)/
305  (classic_electr_radius*classic_electr_radius*(theEGrid[ie]+2.0*electron_mass_c2));
306  G4double fnorm = (*tempData)[ie]/(rsum*fact);
307  G4double TST = 100.*std::fabs(fnorm-1.0);
308  if (TST > 1.0)
309  {
311  ed << "G4PenelopeBremsstrahlungFS. Corrupted data files?" << G4endl;
312  G4cout << "TST= " << TST << "; fnorm = " << fnorm << G4endl;
313  G4cout << "rsum = " << rsum << G4endl;
314  G4cout << "fact = " << fact << G4endl;
315  G4cout << ie << " " << theEGrid[ie]/keV << " " << (*tempData)[ie]/barn << G4endl;
316  G4Exception("G4PenelopeBremsstrahlungFS::BuildScaledXSTable()",
317  "em2010",FatalException,ed);
318  }
319  for (size_t ix=0;ix<nBinsX;ix++)
320  (*tempMatrix)[ie*nBinsX+ix] *= fnorm;
321  }
322 
323  //*********************************************************************
324  // create and fill the tables
325  //*********************************************************************
326  G4PhysicsTable* thePhysicsTable = new G4PhysicsTable();
327  // the table will contain 32 G4PhysicsFreeVectors with different
328  // values of x. Each of the G4PhysicsFreeVectors has a profile of
329  // log(XS) vs. log(E)
330 
331  //reserve space of the vectors. Everything is log-log
332  //I add one extra "fake" point at low energy, since the Penelope
333  //table starts at 1 keV
334  for (size_t i=0;i<nBinsX;i++)
335  thePhysicsTable->push_back(new G4PhysicsFreeVector(nBinsE+1));
336 
337  for (size_t ix=0;ix<nBinsX;ix++)
338  {
339  G4PhysicsFreeVector* theVec =
340  (G4PhysicsFreeVector*) ((*thePhysicsTable)[ix]);
341  for (size_t ie=0;ie<nBinsE;ie++)
342  {
343  G4double logene = std::log(theEGrid[ie]);
344  G4double aValue = (*tempMatrix)[ie*nBinsX+ix];
345  if (aValue < 1e-20*millibarn) //protection against log(0)
346  aValue = 1e-20*millibarn;
347  theVec->PutValue(ie+1,logene,std::log(aValue));
348  }
349  //Add fake point at 1 eV using an extrapolation with the derivative
350  //at the first valid point (Penelope approach)
351  G4double derivative = ((*theVec)[2]-(*theVec)[1])/(theVec->Energy(2) - theVec->Energy(1));
352  G4double log1eV = std::log(1*eV);
353  G4double val1eV = (*theVec)[1]+derivative*(log1eV-theVec->Energy(1));
354  //fake point at very low energy
355  theVec->PutValue(0,log1eV,val1eV);
356  }
357  std::pair<const G4Material*,G4double> theKey = std::make_pair(material,cut);
358  theReducedXSTable->insert(std::make_pair(theKey,thePhysicsTable));
359 
360  delete StechiometricFactors;
361  delete tempData;
362  delete tempMatrix;
363 
364  //Do here also the initialization of the energy sampling
365  if (!(theSamplingTable->count(theKey)))
366  InitializeEnergySampling(material,cut);
367 
368  return;
369 }
370 
371 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
372 
374 {
375 
376  char* path = getenv("G4LEDATA");
377  if (!path)
378  {
379  G4String excep = "G4PenelopeBremsstrahlungFS - G4LEDATA environment variable not set!";
380  G4Exception("G4PenelopeBremsstrahlungFS::ReadDataFile()",
381  "em0006",FatalException,excep);
382  return;
383  }
384  /*
385  Read the cross section file
386  */
387  std::ostringstream ost;
388  if (Z>9)
389  ost << path << "/penelope/bremsstrahlung/pdebr" << Z << ".p08";
390  else
391  ost << path << "/penelope/bremsstrahlung/pdebr0" << Z << ".p08";
392  std::ifstream file(ost.str().c_str());
393  if (!file.is_open())
394  {
395  G4String excep = "G4PenelopeBremsstrahlungFS - data file " +
396  G4String(ost.str()) + " not found!";
397  G4Exception("G4PenelopeBremsstrahlungFS::ReadDataFile()",
398  "em0003",FatalException,excep);
399  return;
400  }
401 
402  G4int readZ =0;
403  file >> readZ;
404 
405  //check the right file is opened.
406  if (readZ != Z)
407  {
409  ed << "Corrupted data file for Z=" << Z << G4endl;
410  G4Exception("G4PenelopeBremsstrahlungFS::ReadDataFile()",
411  "em0005",FatalException,ed);
412  return;
413  }
414 
415  G4DataVector* theMatrix = new G4DataVector(nBinsE*(nBinsX+1),0.); //initialized with zeros
416 
417  for (size_t ie=0;ie<nBinsE;ie++)
418  {
419  G4double myDouble = 0;
420  file >> myDouble; //energy (eV)
421  if (!theEGrid[ie]) //fill only the first time
422  theEGrid[ie] = myDouble*eV;
423  //
424  for (size_t ix=0;ix<nBinsX;ix++)
425  {
426  file >> myDouble;
427  (*theMatrix)[ie*(nBinsX+1)+ix] = myDouble*millibarn;
428  }
429  file >> myDouble; //total cross section
430  (*theMatrix)[ie*(nBinsX+1)+nBinsX] = myDouble*millibarn;
431  }
432 
433  if (theElementData)
434  theElementData->insert(std::make_pair(Z,theMatrix));
435  else
436  delete theMatrix;
437  file.close();
438  return;
439 }
440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
441 
443  G4double xup,G4int momOrder) const
444 //x is always the gridX
445 {
446  //Corresponds to the function RLMOM of Penelope
447  //This method performs the calculation of the integral of (x^momOrder)*y over the interval
448  //from x[0] to xup, obtained by linear interpolation on a table of y.
449  //The independent variable is assumed to take positive values only.
450  //
451  size_t size = nBinsX;
452  const G4double eps = 1e-35;
453 
454  //Check that the call is valid
455  if (momOrder<-1 || size<2 || theXGrid[0]<0)
456  {
457  G4Exception("G4PenelopeBremsstrahlungFS::GetMomentumIntegral()",
458  "em2011",FatalException,"Invalid call");
459  }
460 
461  for (size_t i=1;i<size;i++)
462  {
463  if (theXGrid[i]<0 || theXGrid[i]<theXGrid[i-1])
464  {
466  ed << "Invalid call for bin " << i << G4endl;
467  G4Exception("G4PenelopeBremsstrahlungFS::GetMomentumIntegral()",
468  "em2012",FatalException,ed);
469  }
470  }
471 
472  //Compute the integral
473  G4double result = 0;
474  if (xup < theXGrid[0])
475  return result;
476  bool loopAgain = true;
477  G4double xt = std::min(xup,theXGrid[size-1]);
478  G4double xtc = 0;
479  for (size_t i=0;i<size-1;i++)
480  {
481  G4double x1 = std::max(theXGrid[i],eps);
482  G4double y1 = y[i];
483  G4double x2 = std::max(theXGrid[i+1],eps);
484  G4double y2 = y[i+1];
485  if (xt < x2)
486  {
487  xtc = xt;
488  loopAgain = false;
489  }
490  else
491  xtc = x2;
492  G4double dx = x2-x1;
493  G4double dy = y2-y1;
494  G4double ds = 0;
495  if (std::fabs(dx)>1e-14*std::fabs(dy))
496  {
497  G4double b=dy/dx;
498  G4double a=y1-b*x1;
499  if (momOrder == -1)
500  ds = a*std::log(xtc/x1)+b*(xtc-x1);
501  else if (momOrder == 0) //speed it up, not using pow()
502  ds = a*(xtc-x1) + 0.5*b*(xtc*xtc-x1*x1);
503  else
504  ds = a*(std::pow(xtc,momOrder+1)-std::pow(x1,momOrder+1))/((G4double) (momOrder + 1))
505  + b*(std::pow(xtc,momOrder+2)-std::pow(x1,momOrder+2))/((G4double) (momOrder + 2));
506  }
507  else
508  ds = 0.5*(y1+y2)*(xtc-x1)*std::pow(xtc,momOrder);
509  result += ds;
510  if (!loopAgain)
511  return result;
512  }
513  return result;
514 }
515 
516 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
517 
519  const G4double cut) const
520 {
521  //check if it already contains the entry
522  std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
523 
524  if (!(theReducedXSTable->count(theKey)))
525  {
526  G4Exception("G4PenelopeBremsstrahlungFS::GetScaledXSTable()",
527  "em2013",FatalException,"Unable to retrieve the cross section table");
528  }
529 
530  return theReducedXSTable->find(theKey)->second;
531 }
532 
533 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
534 
536  G4double cut)
537 {
538  if (fVerbosity > 2)
539  G4cout << "Entering in G4PenelopeBremsstrahlungFS::InitializeEnergySampling() for " <<
540  material->GetName() << G4endl;
541 
542  //This method should be accessed by the master only
543  std::pair<const G4Material*,G4double> theKey = std::make_pair(material,cut);
544 
545  G4PhysicsTable* thePhysicsTable = new G4PhysicsTable();
546  // the table will contain 57 G4PhysicsFreeVectors with different
547  // values of E.
549 
550  //I reserve space of the vectors.
551  for (size_t i=0;i<nBinsE;i++)
552  thePhysicsTable->push_back(new G4PhysicsFreeVector(nBinsX));
553 
554  //Retrieve the table. Must already exist at this point, because this
555  //method is invoked by GetScaledXSTable()
556  if (!(theReducedXSTable->count(theKey)))
557  G4Exception("G4PenelopeBremsstrahlungFS::InitializeEnergySampling()",
558  "em2013",FatalException,"Unable to retrieve the cross section table");
559  G4PhysicsTable* theTableReduced = theReducedXSTable->find(theKey)->second;
560 
561  for (size_t ie=0;ie<nBinsE;ie++)
562  {
563  G4PhysicsFreeVector* theVec =
564  (G4PhysicsFreeVector*) ((*thePhysicsTable)[ie]);
565  //Fill the table
566  G4double value = 0; //first value
567  theVec->PutValue(0,theXGrid[0],value);
568  for (size_t ix=1;ix<nBinsX;ix++)
569  {
570  //Here calculate the cumulative distribution
571  // int_{0}^{x} dSigma(x',E)/dx' (1/x') dx'
572  G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTableReduced)[ix-1];
573  G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTableReduced)[ix];
574 
575  G4double x1=std::max(theXGrid[ix-1],1.0e-35);
576  //Remember: the table theReducedXSTable has a fake first point in energy
577  //so, it contains one more bin than nBinsE.
578  G4double y1=G4Exp((*v1)[ie+1]);
579  G4double x2=std::max(theXGrid[ix],1.0e-35);
580  G4double y2=G4Exp((*v2)[ie+1]);
581  G4double B = (y2-y1)/(x2-x1);
582  G4double A = y1-B*x1;
583  G4double dS = A*std::log(x2/x1)+B*(x2-x1);
584  value += dS;
585  theVec->PutValue(ix,theXGrid[ix],value);
586  }
587 
588  //fill the PB vector
589  G4double xc = cut/theEGrid[ie];
590  //Fill a temp data vector
591  G4double* tempData = new G4double[nBinsX];
592  for (size_t ix=0;ix<nBinsX;ix++)
593  {
594  G4PhysicsFreeVector* vv = (G4PhysicsFreeVector*) (*theTableReduced)[ix];
595  tempData[ix] = G4Exp((*vv)[ie+1]);
596  }
597  G4double pbval = (xc<=1) ?
598  GetMomentumIntegral(tempData,xc,-1) :
599  GetMomentumIntegral(tempData,1.0,-1);
600  thePBvec->PutValue(ie,theEGrid[ie],pbval);
601  delete[] tempData;
602  }
603 
604  theSamplingTable->insert(std::make_pair(theKey,thePhysicsTable));
605  thePBcut->insert(std::make_pair(theKey,thePBvec));
606  return;
607 }
608 
609 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
610 
612  const G4double cut) const
613 {
614  std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
615  if (!(theSamplingTable->count(theKey)) || !(thePBcut->count(theKey)))
616  {
618  ed << "Unable to retrieve the SamplingTable: " <<
619  theSamplingTable->count(theKey) << " " <<
620  thePBcut->count(theKey) << G4endl;
621  G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
622  "em2014",FatalException,ed);
623  }
624 
625 
626  const G4PhysicsTable* theTableInte = theSamplingTable->find(theKey)->second;
627  const G4PhysicsTable* theTableRed = theReducedXSTable->find(theKey)->second;
628 
629  //Find the energy bin using bi-partition
630  size_t eBin = 0;
631  G4bool firstOrLastBin = false;
632 
633  if (energy < theEGrid[0]) //below first bin
634  {
635  eBin = 0;
636  firstOrLastBin = true;
637  }
638  else if (energy > theEGrid[nBinsE-1]) //after last bin
639  {
640  eBin = nBinsE-1;
641  firstOrLastBin = true;
642  }
643  else
644  {
645  size_t i=0;
646  size_t j=nBinsE-1;
647  while ((j-i)>1)
648  {
649  size_t k = (i+j)/2;
650  if (energy > theEGrid[k])
651  i = k;
652  else
653  j = k;
654  }
655  eBin = i;
656  }
657 
658  //Get the appropriate physics vector
659  const G4PhysicsFreeVector* theVec1 = (G4PhysicsFreeVector*) (*theTableInte)[eBin];
660 
661  //Use a "temporary" vector which contains the linear interpolation of the x spectra
662  //in energy. The temporary vector is thread-local, so that there is no conflict.
663  //This is achieved via G4Cache. The theTempVect is allocated only once per thread
664  //(member variable), but it is overwritten at every call of this method
665  //(because the interpolation factors change!)
666  G4PhysicsFreeVector* theTempVec = fCache.Get();
667  if (!theTempVec) //First time this thread gets the cache
668  {
669  theTempVec = new G4PhysicsFreeVector(nBinsX);
670  //The G4Physics*Vector pointers are automatically deleted by the G4Allocator,
671  //so there is no need to take care of it manually
672  fCache.Put(theTempVec);
673  if (fVerbosity > 4)
674  G4cout << "Creating new instance of G4PhysicsFreeVector() on the worker" << G4endl;
675  }
676 
677  //theTempVect is allocated only once (member variable), but it is overwritten at
678  //every call of this method (because the interpolation factors change!)
679  if (!firstOrLastBin)
680  {
681  const G4PhysicsFreeVector* theVec2 = (G4PhysicsFreeVector*) (*theTableInte)[eBin+1];
682  for (size_t iloop=0;iloop<nBinsX;iloop++)
683  {
684  G4double val = (*theVec1)[iloop]+(((*theVec2)[iloop]-(*theVec1)[iloop]))*
685  (energy-theEGrid[eBin])/(theEGrid[eBin+1]-theEGrid[eBin]);
686  theTempVec->PutValue(iloop,theXGrid[iloop],val);
687  }
688  }
689  else //first or last bin, no interpolation
690  {
691  for (size_t iloop=0;iloop<nBinsX;iloop++)
692  theTempVec->PutValue(iloop,theXGrid[iloop],(*theVec1)[iloop]);
693  }
694 
695  //Start the game
696  G4double pbcut = (*(thePBcut->find(theKey)->second))[eBin];
697 
698  if (!firstOrLastBin) //linear interpolation on pbcut as well
699  {
700  pbcut = (*(thePBcut->find(theKey)->second))[eBin] +
701  ((*(thePBcut->find(theKey)->second))[eBin+1]-(*(thePBcut->find(theKey)->second))[eBin])*
702  (energy-theEGrid[eBin])/(theEGrid[eBin+1]-theEGrid[eBin]);
703  }
704 
705  G4double pCumulative = (*theTempVec)[nBinsX-1]; //last value
706 
707  G4double eGamma = 0;
708  G4int nIterations = 0;
709  do
710  {
711  G4double pt = pbcut + G4UniformRand()*(pCumulative - pbcut);
712  nIterations++;
713 
714  //find where it is
715  size_t ibin = 0;
716  if (pt < (*theTempVec)[0])
717  ibin = 0;
718  else if (pt > (*theTempVec)[nBinsX-1])
719  {
720  //We observed problems due to numerical rounding here (STT).
721  //delta here is a tiny positive number
722  G4double delta = pt-(*theTempVec)[nBinsX-1];
723  if (delta < pt*1e-10) // very small! Numerical rounding only
724  {
725  ibin = nBinsX-2;
727  ed << "Found that (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt <<
728  " , (*theTempVec)[nBinsX-1] = " << (*theTempVec)[nBinsX-1] << " and delta = " <<
729  (pt-(*theTempVec)[nBinsX-1]) << G4endl;
730  ed << "Possible symptom of problem with numerical precision" << G4endl;
731  G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
732  "em2015",JustWarning,ed);
733  }
734  else //real problem
735  {
737  ed << "Crash at (pt > (*theTempVec)[nBinsX-1]) with pt = " << pt <<
738  " , (*theTempVec)[nBinsX-1]=" << (*theTempVec)[nBinsX-1] << " and nBinsX = " <<
739  nBinsX << G4endl;
740  ed << "Material: " << mat->GetName() << ", energy = " << energy/keV << " keV" <<
741  G4endl;
742  G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()",
743  "em2015",FatalException,ed);
744  }
745  }
746  else
747  {
748  size_t i=0;
749  size_t j=nBinsX-1;
750  while ((j-i)>1)
751  {
752  size_t k = (i+j)/2;
753  if (pt > (*theTempVec)[k])
754  i = k;
755  else
756  j = k;
757  }
758  ibin = i;
759  }
760 
761  G4double w1 = theXGrid[ibin];
762  G4double w2 = theXGrid[ibin+1];
763 
764  const G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTableRed)[ibin];
765  const G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTableRed)[ibin+1];
766  //Remember: the table theReducedXSTable has a fake first point in energy
767  //so, it contains one more bin than nBinsE.
768  G4double pdf1 = G4Exp((*v1)[eBin+1]);
769  G4double pdf2 = G4Exp((*v2)[eBin+1]);
770  G4double deltaW = w2-w1;
771  G4double dpdfb = pdf2-pdf1;
772  G4double B = dpdfb/deltaW;
773  G4double A = pdf1-B*w1;
774  //I already made an interpolation in energy, so I can use the actual value for the
775  //calculation of the wbcut, instead of the grid values (except for the last bin)
776  G4double wbcut = (cut < energy) ? cut/energy : 1.0;
777  if (firstOrLastBin) //this is an particular case: no interpolation available
778  wbcut = (cut < theEGrid[eBin]) ? cut/theEGrid[eBin] : 1.0;
779 
780  if (w1 < wbcut)
781  w1 = wbcut;
782  if (w2 < w1)
783  {
784  //This configuration can happen if initially wbcut > w2 > w1. Due to the previous
785  //statement, (w1 = wbcut), it becomes wbcut = w1 > w2. In this case, it is not a
786  //real problem. It becomes a problem if w2 < w1 before the w1 = wbcut statement. Issue
787  //a warning only in this specific case.
788  if (w2 > wbcut)
789  {
791  ed << "Warning in G4PenelopeBremsstrahlungFS::SampleX()" << G4endl;
792  ed << "Conflicting end-point values: w1=" << w1 << "; w2 = " << w2 << G4endl;
793  ed << "wbcut = " << wbcut << " energy= " << energy/keV << " keV" << G4endl;
794  ed << "cut = " << cut/keV << " keV" << G4endl;
795  G4Exception("G4PenelopeBremsstrahlungFS::SampleGammaEnergy()","em2015",
796  JustWarning,ed);
797  }
798  return w1*energy;
799  }
800 
801  G4double pmax = std::max(A+B*w1,A+B*w2);
802  G4bool loopAgain = false;
803  do
804  {
805  loopAgain = false;
806  eGamma = w1* std::pow((w2/w1),G4UniformRand());
807  if (G4UniformRand()*pmax > (A+B*eGamma))
808  loopAgain = true;
809  }while(loopAgain);
810  eGamma *= energy;
811  if (nIterations > 100) //protection against infinite loops
812  return eGamma;
813  }while(eGamma < cut); //repeat if sampled sub-cut!
814 
815  return eGamma;
816 }
std::map< G4int, G4DataVector * > * theElementData
G4double GetEffectiveZSquared(const G4Material *mat) const
Master and workers (do not touch tables) All of them are const.
G4double GetMomentumIntegral(G4double *y, G4double up, G4int momOrder) const
std::vector< G4Element * > G4ElementVector
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
const G4double fact
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
const G4String & GetName() const
Definition: G4Material.hh:178
void push_back(G4PhysicsVector *)
value_type & Get() const
Definition: G4Cache.hh:282
std::map< std::pair< const G4Material *, G4double >, G4PhysicsTable * > * theSamplingTable
G4double SampleGammaEnergy(G4double energy, const G4Material *, const G4double cut) const
std::map< const G4Material *, G4double > * theEffectiveZSq
static const G4double eps
G4double a
Definition: TRTMaterials.hh:39
double B(double temperature)
void InitializeEnergySampling(const G4Material *, G4double cut)
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
void BuildScaledXSTable(const G4Material *material, G4double cut, G4bool isMaster)
void ClearTables(G4bool isMaster=true)
Reserved for the master model: they build and handle tables.
bool G4bool
Definition: G4Types.hh:79
std::ostream & tab(std::ostream &)
Definition: CCalutils.cc:89
G4double Energy(size_t index) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static const double eV
Definition: G4SIunits.hh:212
std::map< std::pair< const G4Material *, G4double >, G4PhysicsFreeVector * > * thePBcut
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4PenelopeBremsstrahlungFS(G4int verbosity=0)
Only master models are supposed to create instances.
G4double energy(const ThreeVector &p, const G4double m)
static const double g
Definition: G4SIunits.hh:180
G4Cache< G4PhysicsFreeVector * > fCache
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static const double millibarn
Definition: G4SIunits.hh:105
std::map< std::pair< const G4Material *, G4double >, G4PhysicsTable * > * theReducedXSTable
static const double mole
Definition: G4SIunits.hh:283
#define G4endl
Definition: G4ios.hh:61
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
static const double barn
Definition: G4SIunits.hh:104
static const double keV
Definition: G4SIunits.hh:213
double G4double
Definition: G4Types.hh:76
const G4double * GetFractionVector() const
Definition: G4Material.hh:194
const G4PhysicsTable * GetScaledXSTable(const G4Material *, const G4double cut) const
void Put(const value_type &val) const
Definition: G4Cache.hh:286