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