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