| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // Author: Luciano Pandola
28 //
29 // History:
30 // --------
31 // 23 Nov 2010 L Pandola First complete i
32 // 02 May 2011 L.Pandola Remove dependenc
33 // 24 May 2011 L.Pandola Renamed (make v2
34 // 03 Oct 2013 L.Pandola Migration to MT
35 // 30 Oct 2013 L.Pandola Use G4Cache to a
36 // data vector on
37 //
38
39 //....oooOO0OOooo........oooOO0OOooo........oo
40 #include "G4PenelopeBremsstrahlungFS.hh"
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........oo
51
52 G4PenelopeBremsstrahlungFS::G4PenelopeBremsstr
53 fReducedXSTable(nullptr),fEffectiveZSq(nullp
54 fPBcut(nullptr),fVerbosity(verbosity)
55 {
56 fCache.Put(0);
57 G4double tempvector[fNBinsX] =
58 {1.0e-12,0.025e0,0.05e0,0.075e0,0.1e0,0.15
59 0.3e0,0.35e0,0.4e0,0.45e0,0.5e0,0.55e0,0.6
60 0.75e0,0.8e0,0.85e0,0.9e0,0.925e0,0.95e0,0
61 0.995e0,0.999e0,0.9995e0,0.9999e0,0.99995e
62
63 for (std::size_t ix=0;ix<fNBinsX;++ix)
64 theXGrid[ix] = tempvector[ix];
65
66 for (std::size_t i=0;i<fNBinsE;++i)
67 theEGrid[i] = 0.;
68
69 fElementData = new std::map<G4int,G4DataVect
70 }
71
72 //....oooOO0OOooo........oooOO0OOooo........oo
73
74 G4PenelopeBremsstrahlungFS::~G4PenelopeBremsst
75 {
76 ClearTables();
77
78 //The G4Physics*Vector pointers contained in
79 //the G4AutoDelete so there is no need to ta
80
81 //Clear manually fElementData
82 if (fElementData)
83 {
84 for (auto& item : (*fElementData))
85 delete item.second;
86 delete fElementData;
87 fElementData = nullptr;
88 }
89 }
90
91 //....oooOO0OOooo........oooOO0OOooo........oo
92
93
94 void G4PenelopeBremsstrahlungFS::ClearTables(G
95 {
96 //Just to check
97 if (!isMaster)
98 G4Exception("G4PenelopeBremsstrahlungFS::C
99 "em0100",FatalException,"Worker thread in
100
101 if (fReducedXSTable)
102 {
103 for (auto& item : (*fReducedXSTable))
104 {
105 G4PhysicsTable* tab = item.second;
106 tab->clearAndDestroy();
107 delete tab;
108 }
109 fReducedXSTable->clear();
110 delete fReducedXSTable;
111 fReducedXSTable = nullptr;
112 }
113
114 if (fSamplingTable)
115 {
116 for (auto& item : (*fSamplingTable))
117 {
118 G4PhysicsTable* tab = item.second;
119 tab->clearAndDestroy();
120 delete tab;
121 }
122 fSamplingTable->clear();
123 delete fSamplingTable;
124 fSamplingTable = nullptr;
125 }
126 if (fPBcut)
127 {
128 /*
129 std::map< std::pair<const G4Material*,G4doub
130 for (kk=fPBcut->begin(); kk != fPBcut->end()
131 delete kk->second;
132 */
133 delete fPBcut;
134 fPBcut = nullptr;
135 }
136
137 if (fEffectiveZSq)
138 {
139 delete fEffectiveZSq;
140 fEffectiveZSq = nullptr;
141 }
142
143 return;
144 }
145
146 //....oooOO0OOooo........oooOO0OOooo........oo
147
148 G4double G4PenelopeBremsstrahlungFS::GetEffect
149 {
150 if (!fEffectiveZSq)
151 {
152 G4ExceptionDescription ed;
153 ed << "The container for the <Z^2> value
154 G4Exception("G4PenelopeBremsstrahlungFS:
155 "em2007",FatalException,ed);
156 return 0;
157 }
158 //found in the table: return it
159 if (fEffectiveZSq->count(material))
160 return fEffectiveZSq->find(material)->seco
161 else
162 {
163 G4ExceptionDescription ed;
164 ed << "The value of <Z^2> is not proper
165 material->GetName() << G4endl;
166 //requires running of BuildScaledXSTable
167 G4Exception("G4PenelopeBremsstrahlungFS:
168 "em2008",FatalException,ed);
169 }
170 return 0;
171 }
172
173 //....oooOO0OOooo........oooOO0OOooo........oo
174
175 void G4PenelopeBremsstrahlungFS::BuildScaledXS
176 G4double cut,G4bool isMaster)
177 {
178 //Corresponds to subroutines EBRaW and EBRaR
179 /*
180 This method generates the table of the sca
181 bremsstrahlung emission for the given mate
182 file. The table is normalized according to
183 */
184
185 //Just to check
186 if (!isMaster)
187 G4Exception("G4PenelopeBremsstrahlungFS::B
188 "em0100",FatalException,"Worker thread in
189
190 if (fVerbosity > 2)
191 {
192 G4cout << "Entering in G4PenelopeBremsst
193 material->GetName() << G4endl;
194 G4cout << "Threshold = " << cut/keV << "
195 G4endl;
196 }
197
198 //This method should be accessed by the mast
199 if (!fSamplingTable)
200 fSamplingTable =
201 new std::map< std::pair<const G4Material
202 if (!fPBcut)
203 fPBcut =
204 new std::map< std::pair<const G4Material
205
206 //check if the container exists (if not, cre
207 if (!fReducedXSTable)
208 fReducedXSTable = new std::map< std::pair<
209 G4PhysicsTable*>;
210 if (!fEffectiveZSq)
211 fEffectiveZSq = new std::map<const G4Mater
212
213 //******************************************
214 //Determine the equivalent atomic number <Z^
215 //******************************************
216 std::vector<G4double> *StechiometricFactors
217 std::size_t nElements = material->GetNumberO
218 const G4ElementVector* elementVector = mater
219 const G4double* fractionVector = material->G
220 for (std::size_t i=0;i<nElements;i++)
221 {
222 G4double fraction = fractionVector[i];
223 G4double atomicWeigth = (*elementVector)
224 StechiometricFactors->push_back(fraction
225 }
226 //Find max
227 G4double MaxStechiometricFactor = 0.;
228 for (std::size_t i=0;i<nElements;i++)
229 {
230 if ((*StechiometricFactors)[i] > MaxStec
231 MaxStechiometricFactor = (*Stechiometr
232 }
233 //Normalize
234 for (std::size_t i=0;i<nElements;i++)
235 (*StechiometricFactors)[i] /= MaxStechiom
236
237 G4double sumz2 = 0;
238 G4double sums = 0;
239 for (std::size_t i=0;i<nElements;i++)
240 {
241 G4double Z = (*elementVector)[i]->GetZ()
242 sumz2 += (*StechiometricFactors)[i]*Z*Z;
243 sums += (*StechiometricFactors)[i];
244 }
245 G4double ZBR2 = sumz2/sums;
246
247 fEffectiveZSq->insert(std::make_pair(materia
248
249 //******************************************
250 // loop on elements and read data files
251 //******************************************
252 G4DataVector* tempData = new G4DataVector(fN
253 G4DataVector* tempMatrix = new G4DataVector(
254
255 for (std::size_t iel=0;iel<nElements;iel++)
256 {
257 G4double Z = (*elementVector)[iel]->GetZ
258 G4int iZ = (G4int) Z;
259 G4double wgt = (*StechiometricFactors)[i
260
261 //the element is not already loaded
262 if (!fElementData->count(iZ))
263 {
264 ReadDataFile(iZ);
265 if (!fElementData->count(iZ))
266 {
267 G4ExceptionDescription ed;
268 ed << "Error in G4PenelopeBremsstrahlu
269 ed << "Unable to retrieve data for ele
270 G4Exception("G4PenelopeBremsstrahlungF
271 "em2009",FatalException,ed);
272 }
273 }
274
275 G4DataVector* atomData = fElementData->f
276
277 for (std::size_t ie=0;ie<fNBinsE;++ie)
278 {
279 (*tempData)[ie] += wgt*(*atomData)[ie*(fNB
280 for (std::size_t ix=0;ix<fNBinsX;++ix)
281 (*tempMatrix)[ie*fNBinsX+ix] += wgt*(*at
282 }
283 }
284
285 //******************************************
286 // the total energy loss spectrum is re-norm
287 // scaled cross section of Berger and Seltze
288 //******************************************
289 for (std::size_t ie=0;ie<fNBinsE;++ie)
290 {
291 //for each energy, calculate integral of
292 G4double* tempData2 = new G4double[fNBin
293 for (std::size_t ix=0;ix<fNBinsX;++ix)
294 tempData2[ix] = (*tempMatrix)[ie*fNBinsX+ix]
295 G4double rsum = GetMomentumIntegral(temp
296 delete[] tempData2;
297 G4double fact = millibarn*(theEGrid[ie]+
298 (classic_electr_radius*classic_electr_radius
299 G4double fnorm = (*tempData)[ie]/(rsum*f
300 G4double TST = 100.*std::fabs(fnorm-1.0)
301 if (TST > 1.0)
302 {
303 G4ExceptionDescription ed;
304 ed << "G4PenelopeBremsstrahlungFS. Corrupt
305 G4cout << "TST= " << TST << "; fnorm = " <
306 G4cout << "rsum = " << rsum << G4endl;
307 G4cout << "fact = " << fact << G4endl;
308 G4cout << ie << " " << theEGrid[ie]/keV <<
309 G4Exception("G4PenelopeBremsstrahlungFS::B
310 "em2010",FatalException,ed);
311 }
312 for (std::size_t ix=0;ix<fNBinsX;++ix)
313 (*tempMatrix)[ie*fNBinsX+ix] *= fnorm;
314 }
315
316 //******************************************
317 // create and fill the tables
318 //******************************************
319 G4PhysicsTable* thePhysicsTable = new G4Phys
320 // the table will contain 32 G4PhysicsFreeVe
321 // values of x. Each of the G4PhysicsFreeVec
322 // log(XS) vs. log(E)
323
324 //reserve space of the vectors. Everything i
325 //I add one extra "fake" point at low energy
326 //table starts at 1 keV
327 for (std::size_t i=0;i<fNBinsX;++i)
328 thePhysicsTable->push_back(new G4PhysicsFr
329
330 for (std::size_t ix=0;ix<fNBinsX;++ix)
331 {
332 G4PhysicsFreeVector* theVec =
333 (G4PhysicsFreeVector*) ((*thePhysicsTable)[i
334 for (std::size_t ie=0;ie<fNBinsE;++ie)
335 {
336 G4double logene = G4Log(theEGrid[ie]);
337 G4double aValue = (*tempMatrix)[ie*fNBinsX
338 if (aValue < 1e-20*millibarn) //protection
339 aValue = 1e-20*millibarn;
340 theVec->PutValues(ie+1,logene,G4Log(aValue
341 }
342 //Add fake point at 1 eV using an extrap
343 //at the first valid point (Penelope app
344 G4double derivative = ((*theVec)[2]-(*th
345 G4double log1eV = G4Log(1*eV);
346 G4double val1eV = (*theVec)[1]+derivativ
347 //fake point at very low energy
348 theVec->PutValues(0,log1eV,val1eV);
349 }
350 std::pair<const G4Material*,G4double> theKey
351 fReducedXSTable->insert(std::make_pair(theKe
352
353 delete StechiometricFactors;
354 delete tempData;
355 delete tempMatrix;
356
357 //Do here also the initialization of the ene
358 if (!(fSamplingTable->count(theKey)))
359 InitializeEnergySampling(material,cut);
360
361 return;
362 }
363
364 //....oooOO0OOooo........oooOO0OOooo........oo
365
366 void G4PenelopeBremsstrahlungFS::ReadDataFile(
367 {
368 const char* path = G4FindDataDir("G4LEDATA")
369 if (!path)
370 {
371 G4String excep = "G4PenelopeBremsstrahlu
372 G4Exception("G4PenelopeBremsstrahlungFS:
373 "em0006",FatalException,excep);
374 return;
375 }
376 /*
377 Read the cross section file
378 */
379 std::ostringstream ost;
380 if (Z>9)
381 ost << path << "/penelope/bremsstrahlung/p
382 else
383 ost << path << "/penelope/bremsstrahlung/p
384 std::ifstream file(ost.str().c_str());
385 if (!file.is_open())
386 {
387 G4String excep = "G4PenelopeBremsstrahlu
388 G4String(ost.str()) + " not found!";
389 G4Exception("G4PenelopeBremsstrahlungFS:
390 "em0003",FatalException,excep);
391 return;
392 }
393
394 G4int readZ =0;
395 file >> readZ;
396
397 //check the right file is opened.
398 if (readZ != Z)
399 {
400 G4ExceptionDescription ed;
401 ed << "Corrupted data file for Z=" << Z
402 G4Exception("G4PenelopeBremsstrahlungFS:
403 "em0005",FatalException,ed);
404 return;
405 }
406
407 G4DataVector* theMatrix = new G4DataVector(f
408
409 for (std::size_t ie=0;ie<fNBinsE;++ie)
410 {
411 G4double myDouble = 0;
412 file >> myDouble; //energy (eV)
413 if (!theEGrid[ie]) //fill only the first
414 theEGrid[ie] = myDouble*eV;
415 //
416 for (std::size_t ix=0;ix<fNBinsX;++ix)
417 {
418 file >> myDouble;
419 (*theMatrix)[ie*(fNBinsX+1)+ix] = myDouble
420 }
421 file >> myDouble; //total cross section
422 (*theMatrix)[ie*(fNBinsX+1)+fNBinsX] = m
423 }
424
425 if (fElementData)
426 fElementData->insert(std::make_pair(Z,theM
427 else
428 delete theMatrix;
429 file.close();
430 return;
431 }
432
433 //....oooOO0OOooo........oooOO0OOooo........oo
434
435 G4double G4PenelopeBremsstrahlungFS::GetMoment
436 G4double xup,G4int momOrder)
437 //x is always the gridX
438 {
439 //Corresponds to the function RLMOM of Penel
440 //This method performs the calculation of th
441 //from x[0] to xup, obtained by linear inter
442 //The independent variable is assumed to tak
443 //
444 std::size_t size = fNBinsX;
445 const G4double eps = 1e-35;
446
447 //Check that the call is valid
448 if (momOrder<-1 || size<2 || theXGrid[0]<0)
449 {
450 G4Exception("G4PenelopeBremsstrahlungFS:
451 "em2011",FatalException,"Invalid call");
452 }
453
454 for (std::size_t i=1;i<size;++i)
455 {
456 if (theXGrid[i]<0 || theXGrid[i]<theXGri
457 {
458 G4ExceptionDescription ed;
459 ed << "Invalid call for bin " << i << G4en
460 G4Exception("G4PenelopeBremsstrahlungFS::G
461 "em2012",FatalException,ed);
462 }
463 }
464
465 //Compute the integral
466 G4double result = 0;
467 if (xup < theXGrid[0])
468 return result;
469 G4bool loopAgain = true;
470 G4double xt = std::min(xup,theXGrid[size-1])
471 G4double xtc = 0;
472 for (std::size_t i=0;i<size-1;++i)
473 {
474 G4double x1 = std::max(theXGrid[i],eps);
475 G4double y1 = y[i];
476 G4double x2 = std::max(theXGrid[i+1],eps
477 G4double y2 = y[i+1];
478 if (xt < x2)
479 {
480 xtc = xt;
481 loopAgain = false;
482 }
483 else
484 xtc = x2;
485 G4double dx = x2-x1;
486 G4double dy = y2-y1;
487 G4double ds = 0;
488 if (std::fabs(dx)>1e-14*std::fabs(dy))
489 {
490 G4double b=dy/dx;
491 G4double a=y1-b*x1;
492 if (momOrder == -1)
493 ds = a*G4Log(xtc/x1)+b*(xtc-x1);
494 else if (momOrder == 0) //speed it up, not
495 ds = a*(xtc-x1) + 0.5*b*(xtc*xtc-x1*x1);
496 else
497 ds = a*(std::pow(xtc,momOrder+1)-std::po
498 + b*(std::pow(xtc,momOrder+2)-std::pow
499 }
500 else
501 ds = 0.5*(y1+y2)*(xtc-x1)*std::pow(xtc,momOr
502 result += ds;
503 if (!loopAgain)
504 return result;
505 }
506 return result;
507 }
508
509 //....oooOO0OOooo........oooOO0OOooo........oo
510
511 const G4PhysicsTable* G4PenelopeBremsstrahlung
512 const G4double cut) const
513 {
514 //check if it already contains the entry
515 std::pair<const G4Material*,G4double> theKey
516
517 if (!(fReducedXSTable->count(theKey)))
518 {
519 G4Exception("G4PenelopeBremsstrahlungFS:
520 "em2013",FatalException,"Unable to retri
521 }
522
523 return fReducedXSTable->find(theKey)->second
524 }
525
526 //....oooOO0OOooo........oooOO0OOooo........oo
527
528 void G4PenelopeBremsstrahlungFS::InitializeEne
529 G4double cut)
530 {
531 if (fVerbosity > 2)
532 G4cout << "Entering in G4PenelopeBremsstra
533 material->GetName() << G4endl;
534
535 //This method should be accessed by the mast
536 std::pair<const G4Material*,G4double> theKey
537
538 G4PhysicsTable* thePhysicsTable = new G4Phys
539 // the table will contain 57 G4PhysicsFreeVe
540 // values of E.
541 G4PhysicsFreeVector* thePBvec = new G4Physic
542
543 //I reserve space of the vectors.
544 for (std::size_t i=0;i<fNBinsE;++i)
545 thePhysicsTable->push_back(new G4PhysicsFr
546
547 //Retrieve the table. Must already exist at
548 //method is invoked by GetScaledXSTable()
549 if (!(fReducedXSTable->count(theKey)))
550 G4Exception("G4PenelopeBremsstrahlungFS::I
551 "em2013",FatalException,"Unable to retriev
552 G4PhysicsTable* theTableReduced = fReducedXS
553
554 for (std::size_t ie=0;ie<fNBinsE;++ie)
555 {
556 G4PhysicsFreeVector* theVec =
557 (G4PhysicsFreeVector*) ((*thePhysicsTable)[i
558 //Fill the table
559 G4double value = 0; //first value
560 theVec->PutValues(0,theXGrid[0],value);
561 for (std::size_t ix=1;ix<fNBinsX;++ix)
562 {
563 //Here calculate the cumulative distributi
564 // int_{0}^{x} dSigma(x',E)/dx' (1/x') dx'
565 G4PhysicsFreeVector* v1 = (G4PhysicsFreeVe
566 G4PhysicsFreeVector* v2 = (G4PhysicsFreeVe
567
568 G4double x1=std::max(theXGrid[ix-1],1.0e-3
569 //Remember: the table fReducedXSTable has
570 //so, it contains one more bin than fNBins
571 G4double y1=G4Exp((*v1)[ie+1]);
572 G4double x2=std::max(theXGrid[ix],1.0e-35)
573 G4double y2=G4Exp((*v2)[ie+1]);
574 G4double B = (y2-y1)/(x2-x1);
575 G4double A = y1-B*x1;
576 G4double dS = A*G4Log(x2/x1)+B*(x2-x1);
577 value += dS;
578 theVec->PutValues(ix,theXGrid[ix],value);
579 }
580 //fill the PB vector
581 G4double xc = cut/theEGrid[ie];
582 //Fill a temp data vector
583 G4double* tempData = new G4double[fNBins
584 for (std::size_t ix=0;ix<fNBinsX;++ix)
585 {
586 G4PhysicsFreeVector* vv = (G4PhysicsFreeVe
587 tempData[ix] = G4Exp((*vv)[ie+1]);
588 }
589 G4double pbval = (xc<=1) ?
590 GetMomentumIntegral(tempData,xc,-1) :
591 GetMomentumIntegral(tempData,1.0,-1);
592 thePBvec->PutValues(ie,theEGrid[ie],pbva
593 delete[] tempData;
594 }
595
596 fSamplingTable->insert(std::make_pair(theKey
597 fPBcut->insert(std::make_pair(theKey,thePBve
598 return;
599 }
600
601 //....oooOO0OOooo........oooOO0OOooo........oo
602
603 G4double G4PenelopeBremsstrahlungFS::SampleGam
604 const G4double cut) const
605 {
606 std::pair<const G4Material*,G4double> theKey
607 if (!(fSamplingTable->count(theKey)) || !(fP
608 {
609 G4ExceptionDescription ed;
610 ed << "Unable to retrieve the SamplingTa
611 fSamplingTable->count(theKey) << " " <<
612 fPBcut->count(theKey) << G4endl;
613 G4Exception("G4PenelopeBremsstrahlungFS:
614 "em2014",FatalException,ed);
615 }
616 const G4PhysicsTable* theTableInte = fSampli
617 const G4PhysicsTable* theTableRed = fReduced
618
619 //Find the energy bin using bi-partition
620 std::size_t eBin = 0;
621 G4bool firstOrLastBin = false;
622
623 if (energy < theEGrid[0]) //below first bin
624 {
625 eBin = 0;
626 firstOrLastBin = true;
627 }
628 else if (energy > theEGrid[fNBinsE-1]) //aft
629 {
630 eBin = fNBinsE-1;
631 firstOrLastBin = true;
632 }
633 else
634 {
635 std::size_t i=0;
636 std::size_t j=fNBinsE-1;
637 while ((j-i)>1)
638 {
639 std::size_t k = (i+j)/2;
640 if (energy > theEGrid[k])
641 i = k;
642 else
643 j = k;
644 }
645 eBin = i;
646 }
647
648 //Get the appropriate physics vector
649 const G4PhysicsFreeVector* theVec1 = (G4Phys
650
651 //Use a "temporary" vector which contains th
652 //in energy. The temporary vector is thread-
653 //This is achieved via G4Cache. The theTempV
654 //(member variable), but it is overwritten a
655 //(because the interpolation factors change!
656 G4PhysicsFreeVector* theTempVec = fCache.Get
657 if (!theTempVec) //First time this thread ge
658 {
659 theTempVec = new G4PhysicsFreeVector(fNB
660 fCache.Put(theTempVec);
661 // The G4AutoDelete takes care here to c
662 G4AutoDelete::Register(theTempVec);
663 if (fVerbosity > 4)
664 G4cout << "Creating new instance of G4Physic
665 }
666
667 //theTempVect is allocated only once (member
668 //every call of this method (because the int
669 if (!firstOrLastBin)
670 {
671 const G4PhysicsFreeVector* theVec2 = (G4
672 for (std::size_t iloop=0;iloop<fNBinsX;+
673 {
674 G4double val = (*theVec1)[iloop]+(((*theVe
675 (energy-theEGrid[eBin])/(theEGrid[eBin+1
676 theTempVec->PutValues(iloop,theXGrid[iloop
677 }
678 }
679 else //first or last bin, no interpolation
680 {
681 for (std::size_t iloop=0;iloop<fNBinsX;+
682 theTempVec->PutValues(iloop,theXGrid[iloop],
683 }
684
685 //Start the game
686 G4double pbcut = (*(fPBcut->find(theKey)->se
687
688 if (!firstOrLastBin) //linear interpolation
689 {
690 pbcut = (*(fPBcut->find(theKey)->second)
691 ((*(fPBcut->find(theKey)->second))[eBin+1]-(
692 (energy-theEGrid[eBin])/(theEGrid[eBin+1]-th
693 }
694
695 G4double pCumulative = (*theTempVec)[fNBinsX
696
697 G4double eGamma = 0;
698 G4int nIterations = 0;
699 do
700 {
701 G4double pt = pbcut + G4UniformRand()*(p
702 nIterations++;
703
704 //find where it is
705 std::size_t ibin = 0;
706 if (pt < (*theTempVec)[0])
707 ibin = 0;
708 else if (pt > (*theTempVec)[fNBinsX-1])
709 {
710 //We observed problems due to numerical ro
711 //delta here is a tiny positive number
712 G4double delta = pt-(*theTempVec)[fNBinsX-
713 if (delta < pt*1e-10) // very small! Numer
714 {
715 ibin = fNBinsX-2;
716 G4ExceptionDescription ed;
717 ed << "Found that (pt > (*theTempVec)[
718 " , (*theTempVec)[fNBinsX-1] = " << (*theT
719 (pt-(*theTempVec)[fNBinsX-1]) << G4endl;
720 ed << "Possible symptom of problem wit
721 G4Exception("G4PenelopeBremsstrahlungF
722 "em2015",JustWarning,ed);
723 }
724 else //real problem
725 {
726 G4ExceptionDescription ed;
727 ed << "Crash at (pt > (*theTempVec)[fN
728 " , (*theTempVec)[fNBinsX-1]=" << (*theTem
729 fNBinsX << G4endl;
730 ed << "Material: " << mat->GetName() <
731 G4endl;
732 G4Exception("G4PenelopeBremsstrahlungF
733 "em2015",FatalException,ed);
734 }
735 }
736 else
737 {
738 std::size_t i=0;
739 std::size_t j=fNBinsX-1;
740 while ((j-i)>1)
741 {
742 std::size_t k = (i+j)/2;
743 if (pt > (*theTempVec)[k])
744 i = k;
745 else
746 j = k;
747 }
748 ibin = i;
749 }
750
751 G4double w1 = theXGrid[ibin];
752 G4double w2 = theXGrid[ibin+1];
753
754 const G4PhysicsFreeVector* v1 = (G4Physi
755 const G4PhysicsFreeVector* v2 = (G4Physi
756 //Remember: the table fReducedXSTable ha
757 //so, it contains one more bin than fNBi
758 G4double pdf1 = G4Exp((*v1)[eBin+1]);
759 G4double pdf2 = G4Exp((*v2)[eBin+1]);
760 G4double deltaW = w2-w1;
761 G4double dpdfb = pdf2-pdf1;
762 G4double B = dpdfb/deltaW;
763 G4double A = pdf1-B*w1;
764 //I already made an interpolation in ene
765 //calculation of the wbcut, instead of t
766 G4double wbcut = (cut < energy) ? cut/e
767 if (firstOrLastBin) //this is an particu
768 wbcut = (cut < theEGrid[eBin]) ? cut/theEGr
769
770 if (w1 < wbcut)
771 w1 = wbcut;
772 if (w2 < w1)
773 {
774 //This configuration can happen if initial
775 //statement, (w1 = wbcut), it becomes wbcu
776 //real problem. It becomes a problem if w2
777 //a warning only in this specific case.
778 if (w2 > wbcut)
779 {
780 G4ExceptionDescription ed;
781 ed << "Warning in G4PenelopeBremsstrah
782 ed << "Conflicting end-point values: w
783 ed << "wbcut = " << wbcut << " energy=
784 ed << "cut = " << cut/keV << " keV" <<
785 G4Exception("G4PenelopeBremsstrahlungF
786 JustWarning,ed);
787 }
788 return w1*energy;
789 }
790
791 G4double pmax = std::max(A+B*w1,A+B*w2);
792 G4bool loopAgain = false;
793 do
794 {
795 loopAgain = false;
796 eGamma = w1* std::pow((w2/w1),G4UniformRan
797 if (G4UniformRand()*pmax > (A+B*eGamma))
798 loopAgain = true;
799 }while(loopAgain);
800 eGamma *= energy;
801 if (nIterations > 100) //protection agai
802 return eGamma;
803 }while(eGamma < cut); //repeat if sampled
804
805 return eGamma;
806 }
807