Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PenelopeBremsstrahlungModel.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: G4PenelopeBremsstrahlungModel.cc 99415 2016-09-21 09:05:43Z gcosmo $
27 //
28 // Author: Luciano Pandola
29 //
30 // History:
31 // --------
32 // 23 Nov 2010 L Pandola First complete implementation, Penelope v2008
33 // 24 May 2011 L. Pandola Renamed (make default Penelope)
34 // 13 Mar 2012 L. Pandola Updated the interface for the angular generator
35 // 18 Jul 2012 L. Pandola Migrate to the new interface of the angular generator, which
36 // now provides the G4ThreeVector and takes care of rotation
37 // 02 Oct 2013 L. Pandola Migrated to MT
38 // 17 Oct 2013 L. Pandola Partially revert the MT migration: the angular generator is
39 // kept as thread-local, and created/managed by the workers.
40 //
41 
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
47 #include "G4ParticleDefinition.hh"
48 #include "G4MaterialCutsCouple.hh"
49 #include "G4ProductionCutsTable.hh"
50 #include "G4DynamicParticle.hh"
51 #include "G4Gamma.hh"
52 #include "G4Electron.hh"
53 #include "G4Positron.hh"
56 #include "G4PhysicsFreeVector.hh"
57 #include "G4PhysicsLogVector.hh"
58 #include "G4PhysicsTable.hh"
59 #include "G4Exp.hh"
60 
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
62 namespace { G4Mutex PenelopeBremsstrahlungModelMutex = G4MUTEX_INITIALIZER; }
63 
65  const G4String& nam)
66  :G4VEmModel(nam),fParticleChange(0),fParticle(0),
67  isInitialised(false),energyGrid(0),
68  XSTableElectron(0),XSTablePositron(0),fPenelopeFSHelper(0),
69  fPenelopeAngular(0),fLocalTable(false)
70 
71 {
72  fIntrinsicLowEnergyLimit = 100.0*eV;
73  fIntrinsicHighEnergyLimit = 100.0*GeV;
74  nBins = 200;
75 
76  if (part)
77  SetParticle(part);
78 
79  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
80  //
82  //
83  verboseLevel= 0;
84 
85  // Verbosity scale:
86  // 0 = nothing
87  // 1 = warning for energy non-conservation
88  // 2 = details of energy budget
89  // 3 = calculation of cross sections, file openings, sampling of atoms
90  // 4 = entering in methods
91 
92  // Atomic deexcitation model activated by default
93  SetDeexcitationFlag(true);
94 
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98 
100 {
101  if (IsMaster() || fLocalTable)
102  {
103  ClearTables();
104  if (fPenelopeFSHelper)
105  delete fPenelopeFSHelper;
106  }
107  // This is thread-local at the moment
108  if (fPenelopeAngular)
109  delete fPenelopeAngular;
110 }
111 
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
113 
115  const G4DataVector& theCuts)
116 {
117  if (verboseLevel > 3)
118  G4cout << "Calling G4PenelopeBremsstrahlungModel::Initialise()" << G4endl;
119 
120  SetParticle(part);
121 
122  if (IsMaster() && part == fParticle)
123  {
124 
125  if (!fPenelopeFSHelper)
126  fPenelopeFSHelper = new G4PenelopeBremsstrahlungFS(verboseLevel);
127  if (!fPenelopeAngular)
128  fPenelopeAngular = new G4PenelopeBremsstrahlungAngular();
129  //Clear and re-build the tables
130  ClearTables();
131 
132  //forces the cleaning of tables, in this specific case
133  if (fPenelopeAngular)
134  fPenelopeAngular->Initialize();
135 
136  //Set the number of bins for the tables. 20 points per decade
137  nBins = (size_t) (20*std::log10(HighEnergyLimit()/LowEnergyLimit()));
138  nBins = std::max(nBins,(size_t)100);
139  energyGrid = new G4PhysicsLogVector(LowEnergyLimit(),
140  HighEnergyLimit(),
141  nBins-1); //one hidden bin is added
142 
143 
144  XSTableElectron = new
145  std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
146  XSTablePositron = new
147  std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
148 
149  G4ProductionCutsTable* theCoupleTable =
151 
152  //Build tables for all materials
153  for (size_t i=0;i<theCoupleTable->GetTableSize();i++)
154  {
155  const G4Material* theMat =
156  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
157  //Forces the building of the helper tables
158  fPenelopeFSHelper->BuildScaledXSTable(theMat,theCuts.at(i),IsMaster());
159  fPenelopeAngular->PrepareTables(theMat,IsMaster());
160  BuildXSTable(theMat,theCuts.at(i));
161 
162  }
163 
164 
165  if (verboseLevel > 2) {
166  G4cout << "Penelope Bremsstrahlung model v2008 is initialized " << G4endl
167  << "Energy range: "
168  << LowEnergyLimit() / keV << " keV - "
169  << HighEnergyLimit() / GeV << " GeV."
170  << G4endl;
171  }
172  }
173 
174  if(isInitialised) return;
176  isInitialised = true;
177 }
178 
179 
181  G4VEmModel *masterModel)
182 {
183  if (verboseLevel > 3)
184  G4cout << "Calling G4PenelopeBremsstrahlungModel::InitialiseLocal()" << G4endl;
185 
186  //
187  //Check that particle matches: one might have multiple master models (e.g.
188  //for e+ and e-).
189  //
190  if (part == fParticle)
191  {
192  //Get the const table pointers from the master to the workers
193  const G4PenelopeBremsstrahlungModel* theModel =
194  static_cast<G4PenelopeBremsstrahlungModel*> (masterModel);
195 
196  //Copy pointers to the data tables
197  energyGrid = theModel->energyGrid;
198  XSTableElectron = theModel->XSTableElectron;
199  XSTablePositron = theModel->XSTablePositron;
200  fPenelopeFSHelper = theModel->fPenelopeFSHelper;
201 
202  //fPenelopeAngular = theModel->fPenelopeAngular;
203 
204  //created in each thread and initialized.
205  if (!fPenelopeAngular)
206  fPenelopeAngular = new G4PenelopeBremsstrahlungAngular();
207  //forces the cleaning of tables, in this specific case
208  if (fPenelopeAngular)
209  fPenelopeAngular->Initialize();
210 
211  G4ProductionCutsTable* theCoupleTable =
213  //Build tables for all materials
214  for (size_t i=0;i<theCoupleTable->GetTableSize();i++)
215  {
216  const G4Material* theMat =
217  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
218  fPenelopeAngular->PrepareTables(theMat,IsMaster());
219  }
220 
221 
222  //copy the data
223  nBins = theModel->nBins;
224 
225  //Same verbosity for all workers, as the master
226  verboseLevel = theModel->verboseLevel;
227  }
228 
229  return;
230 }
231 
232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
233 
235  const G4ParticleDefinition* theParticle,
237  G4double cutEnergy,
238  G4double)
239 {
240  //
241  if (verboseLevel > 3)
242  G4cout << "Calling CrossSectionPerVolume() of G4PenelopeBremsstrahlungModel" << G4endl;
243 
244  SetupForMaterial(theParticle, material, energy);
245 
246  G4double crossPerMolecule = 0.;
247 
248  G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(theParticle,material,
249  cutEnergy);
250 
251  if (theXS)
252  crossPerMolecule = theXS->GetHardCrossSection(energy);
253 
254  G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
255  G4double atPerMol = oscManager->GetAtomsPerMolecule(material);
256 
257  if (verboseLevel > 3)
258  G4cout << "Material " << material->GetName() << " has " << atPerMol <<
259  "atoms per molecule" << G4endl;
260 
261  G4double moleculeDensity = 0.;
262  if (atPerMol)
263  moleculeDensity = atomDensity/atPerMol;
264 
265  G4double crossPerVolume = crossPerMolecule*moleculeDensity;
266 
267  if (verboseLevel > 2)
268  {
269  G4cout << "G4PenelopeBremsstrahlungModel " << G4endl;
270  G4cout << "Mean free path for gamma emission > " << cutEnergy/keV << " keV at " <<
271  energy/keV << " keV = " << (1./crossPerVolume)/mm << " mm" << G4endl;
272  }
273 
274  return crossPerVolume;
275 }
276 
277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
278 
279 //This is a dummy method. Never inkoved by the tracking, it just issues
280 //a warning if one tries to get Cross Sections per Atom via the
281 //G4EmCalculator.
283  G4double,
284  G4double,
285  G4double,
286  G4double,
287  G4double)
288 {
289  G4cout << "*** G4PenelopeBremsstrahlungModel -- WARNING ***" << G4endl;
290  G4cout << "Penelope Bremsstrahlung model v2008 does not calculate cross section _per atom_ " << G4endl;
291  G4cout << "so the result is always zero. For physics values, please invoke " << G4endl;
292  G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl;
293  return 0;
294 }
295 
296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
297 
299  const G4ParticleDefinition* theParticle,
300  G4double kineticEnergy,
301  G4double cutEnergy)
302 {
303  if (verboseLevel > 3)
304  G4cout << "Calling ComputeDEDX() of G4PenelopeBremsstrahlungModel" << G4endl;
305 
306  G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(theParticle,material,
307  cutEnergy);
308  G4double sPowerPerMolecule = 0.0;
309  if (theXS)
310  sPowerPerMolecule = theXS->GetSoftStoppingPower(kineticEnergy);
311 
312 
313  G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
314  G4double atPerMol = oscManager->GetAtomsPerMolecule(material);
315 
316  G4double moleculeDensity = 0.;
317  if (atPerMol)
318  moleculeDensity = atomDensity/atPerMol;
319 
320  G4double sPowerPerVolume = sPowerPerMolecule*moleculeDensity;
321 
322  if (verboseLevel > 2)
323  {
324  G4cout << "G4PenelopeBremsstrahlungModel " << G4endl;
325  G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
326  kineticEnergy/keV << " keV = " <<
327  sPowerPerVolume/(keV/mm) << " keV/mm" << G4endl;
328  }
329  return sPowerPerVolume;
330 }
331 
332 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
333 
334 void G4PenelopeBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>*fvect,
335  const G4MaterialCutsCouple* couple,
336  const G4DynamicParticle* aDynamicParticle,
337  G4double cutG,
338  G4double)
339 {
340  if (verboseLevel > 3)
341  G4cout << "Calling SampleSecondaries() of G4PenelopeBremsstrahlungModel" << G4endl;
342 
343  G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
344  const G4Material* material = couple->GetMaterial();
345 
346  if (kineticEnergy <= fIntrinsicLowEnergyLimit)
347  {
350  return ;
351  }
352 
353  G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
354  //This is the momentum
355  G4ThreeVector initialMomentum = aDynamicParticle->GetMomentum();
356 
357  //Not enough energy to produce a secondary! Return with nothing happened
358  if (kineticEnergy < cutG)
359  return;
360 
361  if (verboseLevel > 3)
362  G4cout << "Going to sample gamma energy for: " <<material->GetName() << " " <<
363  "energy = " << kineticEnergy/keV << ", cut = " << cutG/keV << G4endl;
364 
365  //Sample gamma's energy according to the spectrum
366  G4double gammaEnergy =
367  fPenelopeFSHelper->SampleGammaEnergy(kineticEnergy,material,cutG);
368 
369  if (verboseLevel > 3)
370  G4cout << "Sampled gamma energy: " << gammaEnergy/keV << " keV" << G4endl;
371 
372  //Now sample the direction for the Gamma. Notice that the rotation is already done
373  //Z is unused here, I plug 0. The information is in the material pointer
374  G4ThreeVector gammaDirection1 =
375  fPenelopeAngular->SampleDirection(aDynamicParticle,gammaEnergy,0,material);
376 
377  if (verboseLevel > 3)
378  G4cout << "Sampled cosTheta for e-: " << gammaDirection1.cosTheta() << G4endl;
379 
380  G4double residualPrimaryEnergy = kineticEnergy-gammaEnergy;
381  if (residualPrimaryEnergy < 0)
382  {
383  //Ok we have a problem, all energy goes with the gamma
384  gammaEnergy += residualPrimaryEnergy;
385  residualPrimaryEnergy = 0.0;
386  }
387 
388  //Produce final state according to momentum conservation
389  G4ThreeVector particleDirection1 = initialMomentum - gammaEnergy*gammaDirection1;
390  particleDirection1 = particleDirection1.unit(); //normalize
391 
392  //Update the primary particle
393  if (residualPrimaryEnergy > 0.)
394  {
395  fParticleChange->ProposeMomentumDirection(particleDirection1);
396  fParticleChange->SetProposedKineticEnergy(residualPrimaryEnergy);
397  }
398  else
399  {
401  }
402 
403  //Now produce the photon
405  gammaDirection1,
406  gammaEnergy);
407  fvect->push_back(theGamma);
408 
409  if (verboseLevel > 1)
410  {
411  G4cout << "-----------------------------------------------------------" << G4endl;
412  G4cout << "Energy balance from G4PenelopeBremsstrahlung" << G4endl;
413  G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
414  G4cout << "-----------------------------------------------------------" << G4endl;
415  G4cout << "Outgoing primary energy: " << residualPrimaryEnergy/keV << " keV" << G4endl;
416  G4cout << "Bremsstrahlung photon " << gammaEnergy/keV << " keV" << G4endl;
417  G4cout << "Total final state: " << (residualPrimaryEnergy+gammaEnergy)/keV
418  << " keV" << G4endl;
419  G4cout << "-----------------------------------------------------------" << G4endl;
420  }
421 
422  if (verboseLevel > 0)
423  {
424  G4double energyDiff = std::fabs(residualPrimaryEnergy+gammaEnergy-kineticEnergy);
425  if (energyDiff > 0.05*keV)
426  G4cout << "Warning from G4PenelopeBremsstrahlung: problem with energy conservation: "
427  <<
428  (residualPrimaryEnergy+gammaEnergy)/keV <<
429  " keV (final) vs. " <<
430  kineticEnergy/keV << " keV (initial)" << G4endl;
431  }
432  return;
433 }
434 
435 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
436 
437 void G4PenelopeBremsstrahlungModel::ClearTables()
438 {
439  if (!IsMaster() && !fLocalTable)
440  //Should not be here!
441  G4Exception("G4PenelopeBremsstrahlungModel::ClearTables()",
442  "em0100",FatalException,"Worker thread in this method");
443 
444  if (XSTableElectron)
445  {
446  for (auto& item : (*XSTableElectron))
447  delete item.second;
448  delete XSTableElectron;
449  XSTableElectron = nullptr;
450  }
451  if (XSTablePositron)
452  {
453  for (auto& item : (*XSTablePositron))
454  delete item.second;
455  delete XSTablePositron;
456  XSTablePositron = nullptr;
457  }
458  /*
459  if (energyGrid)
460  delete energyGrid;
461  */
462  if (fPenelopeFSHelper)
463  fPenelopeFSHelper->ClearTables(IsMaster());
464 
465  if (verboseLevel > 2)
466  G4cout << "G4PenelopeBremsstrahlungModel: cleared tables" << G4endl;
467 
468  return;
469 }
470 
471 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
472 
474  const G4MaterialCutsCouple*)
475 {
476  return fIntrinsicLowEnergyLimit;
477 }
478 
479 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
480 
481 void G4PenelopeBremsstrahlungModel::BuildXSTable(const G4Material* mat,G4double cut)
482 {
483  if (!IsMaster() && !fLocalTable)
484  //Should not be here!
485  G4Exception("G4PenelopeBremsstrahlungModel::BuildXSTable()",
486  "em0100",FatalException,"Worker thread in this method");
487 
488  //The key of the map
489  std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
490 
491  //The table already exists
492  if (XSTableElectron->count(theKey) && XSTablePositron->count(theKey))
493  return;
494 
495  //
496  //This method fills the G4PenelopeCrossSection containers for electrons or positrons
497  //and for the given material/cut couple.
498  //Equivalent of subroutines EBRaT and PINaT of Penelope
499  //
500  if (verboseLevel > 2)
501  {
502  G4cout << "G4PenelopeBremsstrahlungModel: going to build cross section table " << G4endl;
503  G4cout << "for e+/e- in " << mat->GetName() << " for Ecut(gamma)= " <<
504  cut/keV << " keV " << G4endl;
505  }
506 
507  //Tables have been already created (checked by GetCrossSectionTableForCouple)
508  if (energyGrid->GetVectorLength() != nBins)
509  {
511  ed << "Energy Grid looks not initialized" << G4endl;
512  ed << nBins << " " << energyGrid->GetVectorLength() << G4endl;
513  G4Exception("G4PenelopeBremsstrahlungModel::BuildXSTable()",
514  "em2016",FatalException,ed);
515  }
516 
517  G4PenelopeCrossSection* XSEntry = new G4PenelopeCrossSection(nBins);
518  G4PenelopeCrossSection* XSEntry2 = new G4PenelopeCrossSection(nBins);
519 
520  const G4PhysicsTable* table = fPenelopeFSHelper->GetScaledXSTable(mat,cut);
521 
522 
523  //loop on the energy grid
524  for (size_t bin=0;bin<nBins;bin++)
525  {
526  G4double energy = energyGrid->GetLowEdgeEnergy(bin);
527  G4double XH0=0, XH1=0, XH2=0;
528  G4double XS0=0, XS1=0, XS2=0;
529 
530  //Global xs factor
531  G4double fact = fPenelopeFSHelper->GetEffectiveZSquared(mat)*
532  ((energy+electron_mass_c2)*(energy+electron_mass_c2)/
533  (energy*(energy+2.0*electron_mass_c2)));
534 
535  G4double restrictedCut = cut/energy;
536 
537  //Now I need the dSigma/dX profile - interpolated on energy - for
538  //the 32-point x grid. Interpolation is log-log
539  size_t nBinsX = fPenelopeFSHelper->GetNBinsX();
540  G4double* tempData = new G4double[nBinsX];
541  G4double logene = std::log(energy);
542  for (size_t ix=0;ix<nBinsX;ix++)
543  {
544  //find dSigma/dx for the given E. X belongs to the 32-point grid.
545  G4double val = (*table)[ix]->Value(logene);
546  tempData[ix] = G4Exp(val); //back to the real value!
547  }
548 
549  G4double XH0A = 0.;
550  if (restrictedCut <= 1) //calculate only if we are above threshold!
551  XH0A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,-1) -
552  fPenelopeFSHelper->GetMomentumIntegral(tempData,restrictedCut,-1);
553  G4double XS1A = fPenelopeFSHelper->GetMomentumIntegral(tempData,
554  restrictedCut,0);
555  G4double XS2A = fPenelopeFSHelper->GetMomentumIntegral(tempData,
556  restrictedCut,1);
557  G4double XH1A=0, XH2A=0;
558  if (restrictedCut <=1)
559  {
560  XH1A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,0) -
561  XS1A;
562  XH2A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,1) -
563  XS2A;
564  }
565  delete[] tempData;
566 
567  XH0 = XH0A*fact;
568  XS1 = XS1A*fact*energy;
569  XH1 = XH1A*fact*energy;
570  XS2 = XS2A*fact*energy*energy;
571  XH2 = XH2A*fact*energy*energy;
572 
573  XSEntry->AddCrossSectionPoint(bin,energy,XH0,XH1,XH2,XS0,XS1,XS2);
574 
575  //take care of positrons
576  G4double posCorrection = GetPositronXSCorrection(mat,energy);
577  XSEntry2->AddCrossSectionPoint(bin,energy,XH0*posCorrection,
578  XH1*posCorrection,
579  XH2*posCorrection,
580  XS0,
581  XS1*posCorrection,
582  XS2*posCorrection);
583  }
584 
585  //Insert in the appropriate table
586  XSTableElectron->insert(std::make_pair(theKey,XSEntry));
587  XSTablePositron->insert(std::make_pair(theKey,XSEntry2));
588 
589  return;
590 }
591 
592 
593 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
594 
596 G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple(const G4ParticleDefinition* part,
597  const G4Material* mat,
598  G4double cut)
599 {
600  if (part != G4Electron::Electron() && part != G4Positron::Positron())
601  {
603  ed << "Invalid particle: " << part->GetParticleName() << G4endl;
604  G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
605  "em0001",FatalException,ed);
606  return nullptr;
607  }
608 
609  if (part == G4Electron::Electron())
610  {
611  //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
612  //not invoked
613  if (!XSTableElectron)
614  {
615  //create a **thread-local** version of the table. Used only for G4EmCalculator and
616  //Unit Tests
617  G4String excep = "The Cross Section Table for e- was not initialized correctly!";
618  G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
619  "em2013",JustWarning,excep);
620  fLocalTable = true;
621  XSTableElectron = new
622  std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
623  if (!energyGrid)
624  energyGrid = new G4PhysicsLogVector(LowEnergyLimit(),
625  HighEnergyLimit(),
626  nBins-1); //one hidden bin is added
627  if (!fPenelopeFSHelper)
628  fPenelopeFSHelper = new G4PenelopeBremsstrahlungFS(verboseLevel);
629  }
630  std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
631  if (XSTableElectron->count(theKey)) //table already built
632  return XSTableElectron->find(theKey)->second;
633  else
634  {
635  //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
636  //not filled up. This can happen in a UnitTest or via G4EmCalculator
637  if (verboseLevel > 0)
638  {
639  //G4Exception (warning) is issued only in verbose mode
641  ed << "Unable to find e- table for " << mat->GetName() << " at Ecut(gamma)= "
642  << cut/keV << " keV " << G4endl;
643  ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
644  G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
645  "em2009",JustWarning,ed);
646  }
647  //protect file reading via autolock
648  G4AutoLock lock(&PenelopeBremsstrahlungModelMutex);
649  fPenelopeFSHelper->BuildScaledXSTable(mat,cut,true); //pretend to be a master
650  BuildXSTable(mat,cut);
651  lock.unlock();
652  //now it should be ok
653  return XSTableElectron->find(theKey)->second;
654  }
655 
656  }
657  if (part == G4Positron::Positron())
658  {
659  //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
660  //not invoked
661  if (!XSTablePositron)
662  {
663  G4String excep = "The Cross Section Table for e+ was not initialized correctly!";
664  G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
665  "em2013",JustWarning,excep);
666  fLocalTable = true;
667  XSTablePositron = new
668  std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
669  if (!energyGrid)
670  energyGrid = new G4PhysicsLogVector(LowEnergyLimit(),
671  HighEnergyLimit(),
672  nBins-1); //one hidden bin is added
673  if (!fPenelopeFSHelper)
674  fPenelopeFSHelper = new G4PenelopeBremsstrahlungFS(verboseLevel);
675  }
676  std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
677  if (XSTablePositron->count(theKey)) //table already built
678  return XSTablePositron->find(theKey)->second;
679  else
680  {
681  //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
682  //not filled up. This can happen in a UnitTest or via G4EmCalculator
683  if (verboseLevel > 0)
684  {
685  //Issue a G4Exception (warning) only in verbose mode
687  ed << "Unable to find e+ table for " << mat->GetName() << " at Ecut(gamma)= "
688  << cut/keV << " keV " << G4endl;
689  ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
690  G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
691  "em2009",JustWarning,ed);
692  }
693  //protect file reading via autolock
694  G4AutoLock lock(&PenelopeBremsstrahlungModelMutex);
695  fPenelopeFSHelper->BuildScaledXSTable(mat,cut,true); //pretend to be a master
696  BuildXSTable(mat,cut);
697  lock.unlock();
698  //now it should be ok
699  return XSTablePositron->find(theKey)->second;
700  }
701  }
702  return nullptr;
703 }
704 
705 
706 
707 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
708 
709 G4double G4PenelopeBremsstrahlungModel::GetPositronXSCorrection(const G4Material* mat,
710  G4double energy)
711 {
712  //The electron-to-positron correction factor is set equal to the ratio of the
713  //radiative stopping powers for positrons and electrons, which has been calculated
714  //by Kim et al. (1986) (cf. Berger and Seltzer, 1982). Here, it is used an
715  //analytical approximation which reproduces the tabulated values with 0.5%
716  //accuracy
717  G4double t=std::log(1.0+1e6*energy/
718  (electron_mass_c2*fPenelopeFSHelper->GetEffectiveZSquared(mat)));
719  G4double corr = 1.0-G4Exp(-t*(1.2359e-1-t*(6.1274e-2-t*
720  (3.1516e-2-t*(7.7446e-3-t*(1.0595e-3-t*
721  (7.0568e-5-t*
722  1.8080e-6)))))));
723  return corr;
724 }
725 
726 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
727 
728 void G4PenelopeBremsstrahlungModel::SetParticle(const G4ParticleDefinition* p)
729 {
730  if(!fParticle) {
731  fParticle = p;
732  }
733 }
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
G4double GetEffectiveZSquared(const G4Material *mat) const
G4double GetMomentumIntegral(G4double *y, G4double up, G4int momOrder) const
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:118
static constexpr double mm
Definition: G4SIunits.hh:115
tuple bin
Definition: plottest35.py:22
G4double GetSoftStoppingPower(G4double energy) const
Returns the total stopping power due to soft collisions.
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
const char * p
Definition: xmltok.h:285
const G4String & GetName() const
Definition: G4Material.hh:178
double cosTheta() const
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:411
G4double SampleGammaEnergy(G4double energy, const G4Material *, const G4double cut) const
G4double GetHardCrossSection(G4double energy) const
Returns hard cross section at the given energy.
G4bool IsMaster() const
Definition: G4VEmModel.hh:725
void PrepareTables(const G4Material *material, G4bool isMaster)
Reserved for Master Model.
size_t GetVectorLength() const
G4double GetLowEdgeEnergy(size_t binNumber) const
G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double out_energy, G4int Z, const G4Material *mat=0)
Samples the direction of the outgoing photon (in global coordinates).
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void AddCrossSectionPoint(size_t binNumber, G4double energy, G4double XH0, G4double XH1, G4double XH2, G4double XS0, G4double XS1, G4double XS2)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
string material
Definition: eplot.py:19
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *)
G4GLOB_DLL std::ostream G4cout
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *theParticle, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
void BuildScaledXSTable(const G4Material *material, G4double cut, G4bool isMaster)
void ClearTables(G4bool isMaster=true)
Reserved for the master model: they build and handle tables.
G4PenelopeBremsstrahlungModel(const G4ParticleDefinition *p=0, const G4String &processName="PenBrem")
const G4ThreeVector & GetMomentumDirection() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *theParticle, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy=DBL_MAX)
static constexpr double eV
Definition: G4SIunits.hh:215
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
float electron_mass_c2
Definition: hepunit.py:274
static G4PenelopeOscillatorManager * GetOscillatorManager()
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
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
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
static G4ProductionCutsTable * GetProductionCutsTable()
G4int G4Mutex
Definition: G4Threading.hh:173
static G4Positron * Positron()
Definition: G4Positron.cc:94
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *)
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
Hep3Vector unit() const
static constexpr double GeV
Definition: G4SIunits.hh:217
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:788
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4double GetAtomsPerMolecule(const G4Material *)
Returns the total number of atoms per molecule.
const G4PhysicsTable * GetScaledXSTable(const G4Material *, const G4double cut) const
static constexpr double keV
Definition: G4SIunits.hh:216
G4ThreeVector GetMomentum() const
const G4Material * GetMaterial() const