Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4PenelopeBremsstrahlungModel.cc
이 파일의 문서화 페이지로 가기
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 {
74  nBins = 200;
75 
76  if (part)
77  SetParticle(part);
78 
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)
127  if (!fPenelopeAngular)
129  //Clear and re-build the tables
130  ClearTables();
131 
132  //forces the cleaning of tables, in this specific case
133  if (fPenelopeAngular)
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);
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());
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;
201 
202  //fPenelopeAngular = theModel->fPenelopeAngular;
203 
204  //created in each thread and initialized.
205  if (!fPenelopeAngular)
207  //forces the cleaning of tables, in this specific case
208  if (fPenelopeAngular)
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();
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,
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 
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)
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 {
477 }
478 
479 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
480 
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 
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  {
527  G4double XH0=0, XH1=0, XH2=0;
528  G4double XS0=0, XS1=0, XS2=0;
529 
530  //Global xs factor
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);
554  restrictedCut,0);
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 
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)
625  HighEnergyLimit(),
626  nBins-1); //one hidden bin is added
627  if (!fPenelopeFSHelper)
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)
671  HighEnergyLimit(),
672  nBins-1); //one hidden bin is added
673  if (!fPenelopeFSHelper)
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 
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/
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 
729 {
730  if(!fParticle) {
731  fParticle = p;
732  }
733 }
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void PrepareTables(const G4Material *material, G4bool isMaster)
Reserved for Master Model.
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:447
G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double out_energy, G4int Z, const G4Material *mat=0)
Samples the direction of the outgoing photon (in global coordinates).
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
G4double GetLowEdgeEnergy(size_t binNumber) const
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:707
static constexpr double keV
Definition: G4SIunits.hh:216
static constexpr double mm
Definition: G4SIunits.hh:115
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *theParticle, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy=DBL_MAX)
#define G4endl
Definition: G4ios.hh:61
const G4ThreeVector & GetMomentumDirection() const
const char * p
Definition: xmltok.h:285
G4double GetAtomsPerMolecule(const G4Material *)
Returns the total number of atoms per molecule.
G4double GetEffectiveZSquared(const G4Material *mat) const
void BuildScaledXSTable(const G4Material *material, G4double cut, G4bool isMaster)
const G4String & GetParticleName() const
void SetParticle(const G4ParticleDefinition *)
G4PenelopeBremsstrahlungModel(const G4ParticleDefinition *p=0, const G4String &processName="PenBrem")
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:616
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *)
const G4String & GetName() const
Definition: G4Material.hh:179
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:116
G4PenelopeBremsstrahlungFS * fPenelopeFSHelper
static G4PenelopeOscillatorManager * GetOscillatorManager()
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:763
double G4double
Definition: G4Types.hh:76
TString part[npart]
Definition: Style.C:32
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:88
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
void AddCrossSectionPoint(size_t binNumber, G4double energy, G4double XH0, G4double XH1, G4double XH2, G4double XS0, G4double XS1, G4double XS2)
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *)
G4PenelopeBremsstrahlungAngular * fPenelopeAngular
static constexpr double electron_mass_c2
G4double GetMomentumIntegral(G4double *y, G4double up, G4int momOrder) const
double energy
Definition: plottest35.C:25
G4PenelopeOscillatorManager * oscManager
double cosTheta() const
static G4Positron * Positron()
Definition: G4Positron.cc:94
static constexpr double eV
Definition: G4SIunits.hh:215
G4double GetPositronXSCorrection(const G4Material *, G4double energy)
Float_t mat
void ClearTables(G4bool isMaster=true)
Reserved for the master model: they build and handle tables.
static G4Electron * Electron()
Definition: G4Electron.cc:94
float bin[41]
Definition: plottest35.C:14
Hep3Vector unit() const
G4double SampleGammaEnergy(G4double energy, const G4Material *, const G4double cut) const
G4ThreeVector GetMomentum() const
static G4ProductionCutsTable * GetProductionCutsTable()
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *theParticle, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
std::map< std::pair< const G4Material *, G4double >, G4PenelopeCrossSection * > * XSTablePositron
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void BuildXSTable(const G4Material *material, G4double cut)
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
G4bool IsMaster() const
Definition: G4VEmModel.hh:700
const G4PhysicsTable * GetScaledXSTable(const G4Material *, const G4double cut) const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
const G4Material * GetMaterial() const
std::map< std::pair< const G4Material *, G4double >, G4PenelopeCrossSection * > * XSTableElectron
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:210
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double GetHardCrossSection(G4double energy) const
Returns hard cross section at the given energy.
G4double GetSoftStoppingPower(G4double energy) const
Returns the total stopping power due to soft collisions.
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:609
static constexpr double GeV
Definition: G4SIunits.hh:217
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
size_t GetVectorLength() const
G4PenelopeCrossSection * GetCrossSectionTableForCouple(const G4ParticleDefinition *, const G4Material *, G4double cut)
std::mutex G4Mutex
Definition: G4Threading.hh:84