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