| 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 // 24 May 2011 L. Pandola Renamed (make de
33 // 13 Mar 2012 L. Pandola Updated the inte
34 // 18 Jul 2012 L. Pandola Migrate to the n
35 // now provides the
36 // 02 Oct 2013 L. Pandola Migrated to MT
37 // 17 Oct 2013 L. Pandola Partially revert
38 // kept as thread-
39 //
40
41 #include "G4PenelopeBremsstrahlungModel.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4PenelopeBremsstrahlungFS.hh"
45 #include "G4PenelopeBremsstrahlungAngular.hh"
46 #include "G4ParticleDefinition.hh"
47 #include "G4MaterialCutsCouple.hh"
48 #include "G4ProductionCutsTable.hh"
49 #include "G4DynamicParticle.hh"
50 #include "G4Gamma.hh"
51 #include "G4Electron.hh"
52 #include "G4Positron.hh"
53 #include "G4PenelopeOscillatorManager.hh"
54 #include "G4PenelopeCrossSection.hh"
55 #include "G4PhysicsFreeVector.hh"
56 #include "G4PhysicsLogVector.hh"
57 #include "G4PhysicsTable.hh"
58 #include "G4Exp.hh"
59
60 //....oooOO0OOooo........oooOO0OOooo........oo
61 namespace { G4Mutex PenelopeBremsstrahlungMod
62
63 G4PenelopeBremsstrahlungModel::G4PenelopeBrems
64 const G4String& nam)
65 :G4VEmModel(nam),fParticleChange(nullptr),fP
66 fPenelopeFSHelper(nullptr),fPenelopeAngular
67 fXSTableElectron(nullptr),fXSTablePositron(
68 fIsInitialised(false),fLocalTable(false)
69
70 {
71 fIntrinsicLowEnergyLimit = 100.0*eV;
72 fIntrinsicHighEnergyLimit = 100.0*GeV;
73 nBins = 200;
74
75 if (part)
76 SetParticle(part);
77
78 SetHighEnergyLimit(fIntrinsicHighEnergyLimit
79 //
80 fOscManager = G4PenelopeOscillatorManager::G
81 //
82 fVerboseLevel= 0;
83 // Verbosity scale:
84 // 0 = nothing
85 // 1 = warning for energy non-conservation
86 // 2 = details of energy budget
87 // 3 = calculation of cross sections, file o
88 // 4 = entering in methods
89
90 // Atomic deexcitation model activated by de
91 SetDeexcitationFlag(true);
92
93 }
94
95 //....oooOO0OOooo........oooOO0OOooo........oo
96
97 G4PenelopeBremsstrahlungModel::~G4PenelopeBrem
98 {
99 if (IsMaster() || fLocalTable)
100 {
101 ClearTables();
102 if (fPenelopeFSHelper)
103 delete fPenelopeFSHelper;
104 }
105 // This is thread-local at the moment
106 if (fPenelopeAngular)
107 delete fPenelopeAngular;
108 }
109
110 //....oooOO0OOooo........oooOO0OOooo........oo
111
112 void G4PenelopeBremsstrahlungModel::Initialise
113 c
114 {
115 if (fVerboseLevel > 3)
116 G4cout << "Calling G4PenelopeBremsstrahlun
117
118 SetParticle(part);
119
120 if (IsMaster() && part == fParticle)
121 {
122 if (!fPenelopeFSHelper)
123 fPenelopeFSHelper = new G4PenelopeBremsstrah
124 if (!fPenelopeAngular)
125 fPenelopeAngular = new G4PenelopeBremsstrahl
126 //Clear and re-build the tables
127 ClearTables();
128
129 //forces the cleaning of tables, in this
130 if (fPenelopeAngular)
131 fPenelopeAngular->Initialize();
132
133 //Set the number of bins for the tables.
134 nBins = (std::size_t) (20*std::log10(Hig
135 nBins = std::max(nBins,(std::size_t)100)
136 fEnergyGrid = new G4PhysicsLogVector(Low
137 HighEner
138 nBins-1)
139
140 fXSTableElectron = new
141 std::map< std::pair<const G4Material*,G4doub
142 fXSTablePositron = new
143 std::map< std::pair<const G4Material*,G4doub
144
145 G4ProductionCutsTable* theCoupleTable =
146 G4ProductionCutsTable::GetProductionCutsTabl
147
148 //Build tables for all materials
149 for (G4int i=0;i<(G4int)theCoupleTable->
150 {
151 const G4Material* theMat =
152 theCoupleTable->GetMaterialCutsCouple(i)
153 //Forces the building of the helper tables
154 fPenelopeFSHelper->BuildScaledXSTable(theM
155 fPenelopeAngular->PrepareTables(theMat,IsM
156 BuildXSTable(theMat,theCuts.at(i));
157
158 }
159
160 if (fVerboseLevel > 2) {
161 G4cout << "Penelope Bremsstrahlung model v20
162 << "Energy range: "
163 << LowEnergyLimit() / keV << " keV -
164 << HighEnergyLimit() / GeV << " GeV."
165 << G4endl;
166 }
167 }
168
169 if(fIsInitialised) return;
170 fParticleChange = GetParticleChangeForLoss()
171 fIsInitialised = true;
172 }
173
174 //....oooOO0OOooo........oooOO0OOooo........oo
175
176 void G4PenelopeBremsstrahlungModel::Initialise
177 G4VEmModel *masterModel)
178 {
179 if (fVerboseLevel > 3)
180 G4cout << "Calling G4PenelopeBremsstrahlu
181 //
182 //Check that particle matches: one might hav
183 //for e+ and e-).
184 //
185 if (part == fParticle)
186 {
187 //Get the const table pointers from the
188 const G4PenelopeBremsstrahlungModel* the
189 static_cast<G4PenelopeBremsstrahlungModel*>
190
191 //Copy pointers to the data tables
192 fEnergyGrid = theModel->fEnergyGrid;
193 fXSTableElectron = theModel->fXSTableEle
194 fXSTablePositron = theModel->fXSTablePos
195 fPenelopeFSHelper = theModel->fPenelopeF
196
197 //created in each thread and initialized
198 if (!fPenelopeAngular)
199 fPenelopeAngular = new G4PenelopeBremsstrahl
200 //forces the cleaning of tables, in this
201 if (fPenelopeAngular)
202 fPenelopeAngular->Initialize();
203
204 G4ProductionCutsTable* theCoupleTable =
205 G4ProductionCutsTable::GetProductionCutsTabl
206 //Build tables for all materials
207 for (G4int i=0;i<(G4int)theCoupleTable->
208 {
209 const G4Material* theMat =
210 theCoupleTable->GetMaterialCutsCouple(i)
211 fPenelopeAngular->PrepareTables(theMat,IsM
212 }
213
214 //copy the data
215 nBins = theModel->nBins;
216
217 //Same verbosity for all workers, as the
218 fVerboseLevel = theModel->fVerboseLevel;
219 }
220 return;
221 }
222
223 //....oooOO0OOooo........oooOO0OOooo........oo
224
225 G4double G4PenelopeBremsstrahlungModel::CrossS
226 con
227 G4d
228 G4d
229 G4d
230 {
231 //
232 if (fVerboseLevel > 3)
233 G4cout << "Calling CrossSectionPerVolume()
234
235 SetupForMaterial(theParticle, material, ener
236 G4double crossPerMolecule = 0.;
237
238 G4PenelopeCrossSection* theXS = GetCrossSect
239
240 if (theXS)
241 crossPerMolecule = theXS->GetHardCrossSect
242
243 G4double atomDensity = material->GetTotNbOfA
244 G4double atPerMol = fOscManager->GetAtomsPe
245
246 if (fVerboseLevel > 3)
247 G4cout << "Material " << material->GetName
248 "atoms per molecule" << G4endl;
249
250 G4double moleculeDensity = 0.;
251 if (atPerMol)
252 moleculeDensity = atomDensity/atPerMol;
253
254 G4double crossPerVolume = crossPerMolecule*m
255
256 if (fVerboseLevel > 2)
257 {
258 G4cout << "G4PenelopeBremsstrahlungModel "
259 G4cout << "Mean free path for gamma emissi
260 energy/keV << " keV = " << (1./crossPerV
261 }
262
263 return crossPerVolume;
264 }
265
266 //....oooOO0OOooo........oooOO0OOooo........oo
267
268 //This is a dummy method. Never inkoved by the
269 //a warning if one tries to get Cross Sections
270 //G4EmCalculator.
271 G4double G4PenelopeBremsstrahlungModel::Comput
272 G4double,
273 G4double,
274 G4double,
275 G4double,
276 G4double)
277 {
278 G4cout << "*** G4PenelopeBremsstrahlungModel
279 G4cout << "Penelope Bremsstrahlung model v20
280 G4cout << "so the result is always zero. For
281 G4cout << "GetCrossSectionPerVolume() or Get
282 return 0;
283 }
284
285 //....oooOO0OOooo........oooOO0OOooo........oo
286
287 G4double G4PenelopeBremsstrahlungModel::Comput
288 const G4ParticleDefinition*
289 G4double kineticEnergy,
290 G4double cutEnergy)
291 {
292 if (fVerboseLevel > 3)
293 G4cout << "Calling ComputeDEDX() of G4Pene
294
295 G4PenelopeCrossSection* theXS = GetCrossSect
296
297 G4double sPowerPerMolecule = 0.0;
298 if (theXS)
299 sPowerPerMolecule = theXS->GetSoftStopping
300
301 G4double atomDensity = material->GetTotNbOfA
302 G4double atPerMol = fOscManager->GetAtomsPe
303
304 G4double moleculeDensity = 0.;
305 if (atPerMol)
306 moleculeDensity = atomDensity/atPerMol;
307
308 G4double sPowerPerVolume = sPowerPerMolecule
309
310 if (fVerboseLevel > 2)
311 {
312 G4cout << "G4PenelopeBremsstrahlungModel
313 G4cout << "Stopping power < " << cutEner
314 kineticEnergy/keV << " keV = " <<
315 sPowerPerVolume/(keV/mm) << " keV/mm"
316 }
317 return sPowerPerVolume;
318 }
319
320 //....oooOO0OOooo........oooOO0OOooo........oo
321
322 void G4PenelopeBremsstrahlungModel::SampleSeco
323 const G4MaterialCutsCouple*
324 const G4DynamicParticle* aDy
325 G4double cutG,
326 G4double)
327 {
328 if (fVerboseLevel > 3)
329 G4cout << "Calling SampleSecondaries() of
330
331 G4double kineticEnergy = aDynamicParticle->G
332 const G4Material* material = couple->GetMate
333
334 if (kineticEnergy <= fIntrinsicLowEnergyLimi
335 {
336 fParticleChange->SetProposedKineticEnergy
337 fParticleChange->ProposeLocalEnergyDeposi
338 return ;
339 }
340
341 G4ParticleMomentum particleDirection0 = aDyn
342 //This is the momentum
343 G4ThreeVector initialMomentum = aDynamicPar
344
345 //Not enough energy to produce a secondary!
346 if (kineticEnergy < cutG)
347 return;
348
349 if (fVerboseLevel > 3)
350 G4cout << "Going to sample gamma energy fo
351 "energy = " << kineticEnergy/keV << ", c
352
353 //Sample gamma's energy according to the sp
354 G4double gammaEnergy =
355 fPenelopeFSHelper->SampleGammaEnergy(kinet
356
357 if (fVerboseLevel > 3)
358 G4cout << "Sampled gamma energy: " << gamm
359
360 //Now sample the direction for the Gamma. No
361 //Z is unused here, I plug 0. The informatio
362 G4ThreeVector gammaDirection1 =
363 fPenelopeAngular->SampleDirection(aDynami
364
365 if (fVerboseLevel > 3)
366 G4cout << "Sampled cosTheta for e-: " << g
367
368 G4double residualPrimaryEnergy = kineticEner
369 if (residualPrimaryEnergy < 0)
370 {
371 //Ok we have a problem, all energy goes
372 gammaEnergy += residualPrimaryEnergy;
373 residualPrimaryEnergy = 0.0;
374 }
375
376 //Produce final state according to momentum
377 G4ThreeVector particleDirection1 = initialMo
378 particleDirection1 = particleDirection1.unit
379
380 //Update the primary particle
381 if (residualPrimaryEnergy > 0.)
382 {
383 fParticleChange->ProposeMomentumDirectio
384 fParticleChange->SetProposedKineticEnerg
385 }
386 else
387 {
388 fParticleChange->SetProposedKineticEnerg
389 }
390
391 //Now produce the photon
392 G4DynamicParticle* theGamma = new G4DynamicP
393
394
395 fvect->push_back(theGamma);
396
397 if (fVerboseLevel > 1)
398 {
399 G4cout << "-----------------------------
400 G4cout << "Energy balance from G4Penelop
401 G4cout << "Incoming primary energy: " <<
402 G4cout << "-----------------------------
403 G4cout << "Outgoing primary energy: " <<
404 G4cout << "Bremsstrahlung photon " << ga
405 G4cout << "Total final state: " << (resi
406 << " keV" << G4endl;
407 G4cout << "-----------------------------
408 }
409
410 if (fVerboseLevel > 0)
411 {
412 G4double energyDiff = std::fabs(residual
413 if (energyDiff > 0.05*keV)
414 G4cout << "Warning from G4PenelopeBrem
415 <<
416 (residualPrimaryEnergy+gammaEnergy)/
417 " keV (final) vs. " <<
418 kineticEnergy/keV << " keV (initial)
419 }
420 return;
421 }
422
423 //....oooOO0OOooo........oooOO0OOooo........oo
424
425 void G4PenelopeBremsstrahlungModel::ClearTable
426 {
427 if (!IsMaster() && !fLocalTable)
428 //Should not be here!
429 G4Exception("G4PenelopeBremsstrahlungModel
430 "em0100",FatalException,"Worker thread in
431
432 if (fXSTableElectron)
433 {
434 for (auto& item : (*fXSTableElectron))
435 delete item.second;
436 delete fXSTableElectron;
437 fXSTableElectron = nullptr;
438 }
439 if (fXSTablePositron)
440 {
441 for (auto& item : (*fXSTablePositron))
442 delete item.second;
443 delete fXSTablePositron;
444 fXSTablePositron = nullptr;
445 }
446 /*
447 if (fEnergyGrid)
448 delete fEnergyGrid;
449 */
450 if (fPenelopeFSHelper)
451 fPenelopeFSHelper->ClearTables(IsMaster())
452
453 if (fVerboseLevel > 2)
454 G4cout << "G4PenelopeBremsstrahlungModel:
455
456 return;
457 }
458
459 //....oooOO0OOooo........oooOO0OOooo........oo
460
461 G4double G4PenelopeBremsstrahlungModel::MinEne
462 const G4MaterialCutsCouple*
463 {
464 return fIntrinsicLowEnergyLimit;
465 }
466
467 //....oooOO0OOooo........oooOO0OOooo........oo
468
469 void G4PenelopeBremsstrahlungModel::BuildXSTab
470 {
471 if (!IsMaster() && !fLocalTable)
472 //Should not be here!
473 G4Exception("G4PenelopeBremsstrahlungModel
474 "em0100",FatalException,"Worker thread in
475
476 //The key of the map
477 std::pair<const G4Material*,G4double> theKey
478
479 //The table already exists
480 if (fXSTableElectron->count(theKey) && fXSTa
481 return;
482
483 //
484 //This method fills the G4PenelopeCrossSecti
485 //and for the given material/cut couple.
486 //Equivalent of subroutines EBRaT and PINaT
487 //
488 if (fVerboseLevel > 2)
489 {
490 G4cout << "G4PenelopeBremsstrahlungModel
491 G4cout << "for e+/e- in " << mat->GetNam
492 cut/keV << " keV " << G4endl;
493 }
494
495 //Tables have been already created (checked
496 if (fEnergyGrid->GetVectorLength() != nBins)
497 {
498 G4ExceptionDescription ed;
499 ed << "Energy Grid looks not initialized
500 ed << nBins << " " << fEnergyGrid->GetVe
501 G4Exception("G4PenelopeBremsstrahlungMod
502 "em2016",FatalException,ed);
503 }
504
505 G4PenelopeCrossSection* XSEntry = new G4Pene
506 G4PenelopeCrossSection* XSEntry2 = new G4Pen
507 const G4PhysicsTable* table = fPenelopeFSHel
508
509 //loop on the energy grid
510 for (std::size_t bin=0;bin<nBins;++bin)
511 {
512 G4double energy = fEnergyGrid->GetLowEd
513 G4double XH0=0, XH1=0, XH2=0;
514 G4double XS0=0, XS1=0, XS2=0;
515
516 //Global xs factor
517 G4double fact = fPenelopeFSHelper->GetE
518 ((energy+electron_mass_c2)*(energy+electron
519 (energy*(energy+2.0*electron_mass_c2)));
520
521 G4double restrictedCut = cut/energy;
522
523 //Now I need the dSigma/dX profile - in
524 //the 32-point x grid. Interpolation is
525 std::size_t nBinsX = fPenelopeFSHelper-
526 G4double* tempData = new G4double[nBins
527 G4double logene = G4Log(energy);
528 for (std::size_t ix=0;ix<nBinsX;++ix)
529 {
530 //find dSigma/dx for the given E. X belon
531 G4double val = (*table)[ix]->Value(logene
532 tempData[ix] = G4Exp(val); //back to the
533 }
534
535 G4double XH0A = 0.;
536 if (restrictedCut <= 1) //calculate onl
537 XH0A = fPenelopeFSHelper->GetMomentumIntegr
538 fPenelopeFSHelper->GetMomentumIntegral(te
539 G4double XS1A = fPenelopeFSHelper->GetM
540 restrictedCut,0);
541 G4double XS2A = fPenelopeFSHelper->GetM
542 restrictedCut,1);
543 G4double XH1A=0, XH2A=0;
544 if (restrictedCut <=1)
545 {
546 XH1A = fPenelopeFSHelper->GetMomentumInte
547 XS1A;
548 XH2A = fPenelopeFSHelper->GetMomentumInte
549 XS2A;
550 }
551 delete[] tempData;
552
553 XH0 = XH0A*fact;
554 XS1 = XS1A*fact*energy;
555 XH1 = XH1A*fact*energy;
556 XS2 = XS2A*fact*energy*energy;
557 XH2 = XH2A*fact*energy*energy;
558
559 XSEntry->AddCrossSectionPoint(bin,energ
560
561 //take care of positrons
562 G4double posCorrection = GetPositronXSC
563 XSEntry2->AddCrossSectionPoint(bin,ener
564 XH1*posCorrection,
565 XH2*posCorrection,
566 XS0,
567 XS1*posCorrection,
568 XS2*posCorrection);
569 }
570
571 //Insert in the appropriate table
572 fXSTableElectron->insert(std::make_pair(theK
573 fXSTablePositron->insert(std::make_pair(theK
574
575 return;
576 }
577
578
579 //....oooOO0OOooo........oooOO0OOooo........oo
580
581 G4PenelopeCrossSection*
582 G4PenelopeBremsstrahlungModel::GetCrossSection
583 const G4Material* mat,
584 G4double cut)
585 {
586 if (part != G4Electron::Electron() && part !
587 {
588 G4ExceptionDescription ed;
589 ed << "Invalid particle: " << part->GetP
590 G4Exception("G4PenelopeBremsstrahlungMod
591 "em0001",FatalException,ed);
592 return nullptr;
593 }
594
595 if (part == G4Electron::Electron())
596 {
597 //Either Initialize() was not called, or
598 //not invoked
599 if (!fXSTableElectron)
600 {
601 //create a **thread-local** version of the
602 //Unit Tests
603 G4String excep = "The Cross Section
604 G4Exception("G4PenelopeBremsstrahlun
605 "em2013",JustWarning,excep);
606 fLocalTable = true;
607 fXSTableElectron = new
608 std::map< std::pair<const G4Material*,G4
609 if (!fEnergyGrid)
610 fEnergyGrid = new G4PhysicsLogVector(Low
611 HighEnergyLimit(),
612 nBins-1); //one hidden bin is adde
613 if (!fPenelopeFSHelper)
614 fPenelopeFSHelper = new G4PenelopeBremss
615 }
616 std::pair<const G4Material*,G4double> th
617 if (fXSTableElectron->count(theKey)) //t
618 return fXSTableElectron->find(theKey)-
619 else
620 {
621 //If we are here, it means that Initialize
622 //not filled up. This can happen in a Unit
623 if (fVerboseLevel > 0)
624 {
625 //G4Exception (warning) is issued only
626 G4ExceptionDescription ed;
627 ed << "Unable to find e- table for " <
628 << cut/keV << " keV " << G4endl;
629 ed << "This can happen only in Unit Te
630 G4Exception("G4PenelopeBremsstrahlungM
631 "em2009",JustWarning,ed);
632 }
633 //protect file reading via autolock
634 G4AutoLock lock(&PenelopeBremsstrahlungMod
635 fPenelopeFSHelper->BuildScaledXSTabl
636 BuildXSTable(mat,cut);
637 lock.unlock();
638 //now it should be ok
639 return fXSTableElectron->find(theKey)->sec
640 }
641 }
642 if (part == G4Positron::Positron())
643 {
644 //Either Initialize() was not called, or
645 //not invoked
646 if (!fXSTablePositron)
647 {
648 G4String excep = "The Cross Section Table
649 G4Exception("G4PenelopeBremsstrahlun
650 "em2013",JustWarning,excep);
651 fLocalTable = true;
652 fXSTablePositron = new
653 std::map< std::pair<const G4Material*,G4
654 if (!fEnergyGrid)
655 fEnergyGrid = new G4PhysicsLogVector(Low
656 HighEnergyLimit(),
657 nBins-1); //one hidden bin is adde
658 if (!fPenelopeFSHelper)
659 fPenelopeFSHelper = new G4Penelope
660 }
661 std::pair<const G4Material*,G4double> th
662 if (fXSTablePositron->count(theKey)) //t
663 return fXSTablePositron->find(theKey)-
664 else
665 {
666 //If we are here, it means that Initialize
667 //not filled up. This can happen in a Unit
668 if (fVerboseLevel > 0)
669 {
670 //Issue a G4Exception (warning) only i
671 G4ExceptionDescription ed;
672 ed << "Unable to find e+ table for " <
673 << cut/keV << " keV " << G4endl;
674 ed << "This can happen only in Unit Te
675 G4Exception("G4PenelopeBremsstrahlungM
676 "em2009",JustWarning,ed);
677 }
678 //protect file reading via autolock
679 G4AutoLock lock(&PenelopeBremsstrahlungMod
680 fPenelopeFSHelper->BuildScaledXSTabl
681 BuildXSTable(mat,cut);
682 lock.unlock();
683 //now it should be ok
684 return fXSTablePositron->find(theKey)->sec
685 }
686 }
687 return nullptr;
688 }
689
690 //....oooOO0OOooo........oooOO0OOooo........oo
691
692 G4double G4PenelopeBremsstrahlungModel::GetPos
693 G4double energy)
694 {
695 //The electron-to-positron correction factor
696 //radiative stopping powers for positrons an
697 //by Kim et al. (1986) (cf. Berger and Seltz
698 //analytical approximation which reproduces
699 //accuracy
700 G4double t=G4Log(1.0+1e6*energy/
701 (electron_mass_c2*fPenelopeFSHelper-
702 G4double corr = 1.0-G4Exp(-t*(1.2359e-1-t*(6
703 (3.1516e-2-t*(7.7446e-3-t*(1.0595
704 (7.0568e-5-t*
705 1.8080e-6)))))));
706 return corr;
707 }
708
709 //....oooOO0OOooo........oooOO0OOooo........oo
710
711 void G4PenelopeBremsstrahlungModel::SetParticl
712 {
713 if(!fParticle) {
714 fParticle = p;
715 }
716 }
717