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