Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4PenelopeBremsstrahlungModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 //
 27 // Author: Luciano Pandola
 28 //
 29 // History:
 30 // --------
 31 // 23 Nov 2010   L Pandola    First complete implementation, Penelope v2008
 32 // 24 May 2011   L. Pandola   Renamed (make default Penelope)
 33 // 13 Mar 2012   L. Pandola   Updated the interface for the angular generator
 34 // 18 Jul 2012   L. Pandola   Migrate to the new interface of the angular generator, which
 35 //                            now provides the G4ThreeVector and takes care of rotation
 36 // 02 Oct 2013   L. Pandola   Migrated to MT
 37 // 17 Oct 2013   L. Pandola   Partially revert the MT migration: the angular generator is
 38 //                             kept as thread-local, and created/managed by the workers.
 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........oooOO0OOooo........oooOO0OOooo....
 61 namespace { G4Mutex  PenelopeBremsstrahlungModelMutex = G4MUTEX_INITIALIZER; }
 62 
 63 G4PenelopeBremsstrahlungModel::G4PenelopeBremsstrahlungModel(const G4ParticleDefinition* part,
 64                    const G4String& nam)
 65   :G4VEmModel(nam),fParticleChange(nullptr),fParticle(nullptr),
 66    fPenelopeFSHelper(nullptr),fPenelopeAngular(nullptr),fEnergyGrid(nullptr),
 67    fXSTableElectron(nullptr),fXSTablePositron(nullptr),
 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::GetOscillatorManager();
 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 openings, sampling of atoms
 88   // 4 = entering in methods
 89 
 90   // Atomic deexcitation model activated by default
 91   SetDeexcitationFlag(true);
 92 
 93 }
 94 
 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 96 
 97 G4PenelopeBremsstrahlungModel::~G4PenelopeBremsstrahlungModel()
 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........oooOO0OOooo........oooOO0OOooo....
111 
112 void G4PenelopeBremsstrahlungModel::Initialise(const G4ParticleDefinition* part,
113                                              const G4DataVector& theCuts)
114 {
115   if (fVerboseLevel > 3)
116     G4cout << "Calling G4PenelopeBremsstrahlungModel::Initialise()" << G4endl;
117 
118   SetParticle(part);
119 
120   if (IsMaster() && part == fParticle)
121     {
122       if (!fPenelopeFSHelper)
123   fPenelopeFSHelper = new G4PenelopeBremsstrahlungFS(fVerboseLevel);
124       if (!fPenelopeAngular)
125   fPenelopeAngular = new G4PenelopeBremsstrahlungAngular();
126       //Clear and re-build the tables
127       ClearTables();
128 
129       //forces the cleaning of tables, in this specific case
130       if (fPenelopeAngular)
131   fPenelopeAngular->Initialize();
132 
133       //Set the number of bins for the tables. 20 points per decade
134       nBins = (std::size_t) (20*std::log10(HighEnergyLimit()/LowEnergyLimit()));
135       nBins = std::max(nBins,(std::size_t)100);
136       fEnergyGrid = new G4PhysicsLogVector(LowEnergyLimit(),
137                                       HighEnergyLimit(),
138                                       nBins-1); //one hidden bin is added
139 
140       fXSTableElectron = new
141   std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
142       fXSTablePositron = new
143   std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
144 
145       G4ProductionCutsTable* theCoupleTable =
146   G4ProductionCutsTable::GetProductionCutsTable();
147 
148       //Build tables for all materials
149       for (G4int i=0;i<(G4int)theCoupleTable->GetTableSize();++i)
150   {
151     const G4Material* theMat =
152       theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
153     //Forces the building of the helper tables
154     fPenelopeFSHelper->BuildScaledXSTable(theMat,theCuts.at(i),IsMaster());
155     fPenelopeAngular->PrepareTables(theMat,IsMaster());
156     BuildXSTable(theMat,theCuts.at(i));
157 
158   }
159 
160       if (fVerboseLevel > 2) {
161   G4cout << "Penelope Bremsstrahlung model v2008 is initialized " << G4endl
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........oooOO0OOooo........oooOO0OOooo....
175 
176 void G4PenelopeBremsstrahlungModel::InitialiseLocal(const G4ParticleDefinition* part,
177                 G4VEmModel *masterModel)
178 {
179   if (fVerboseLevel > 3)
180     G4cout << "Calling  G4PenelopeBremsstrahlungModel::InitialiseLocal()" << G4endl;
181   //
182   //Check that particle matches: one might have multiple master models (e.g.
183   //for e+ and e-).
184   //
185   if (part == fParticle)
186     {
187       //Get the const table pointers from the master to the workers
188       const G4PenelopeBremsstrahlungModel* theModel =
189   static_cast<G4PenelopeBremsstrahlungModel*> (masterModel);
190 
191       //Copy pointers to the data tables
192       fEnergyGrid = theModel->fEnergyGrid;
193       fXSTableElectron = theModel->fXSTableElectron;
194       fXSTablePositron = theModel->fXSTablePositron;
195       fPenelopeFSHelper = theModel->fPenelopeFSHelper;
196 
197       //created in each thread and initialized.
198       if (!fPenelopeAngular)
199   fPenelopeAngular = new G4PenelopeBremsstrahlungAngular();
200       //forces the cleaning of tables, in this specific case
201       if (fPenelopeAngular)
202   fPenelopeAngular->Initialize();
203 
204       G4ProductionCutsTable* theCoupleTable =
205   G4ProductionCutsTable::GetProductionCutsTable();
206       //Build tables for all materials
207       for (G4int i=0;i<(G4int)theCoupleTable->GetTableSize();++i)
208   {
209     const G4Material* theMat =
210       theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
211     fPenelopeAngular->PrepareTables(theMat,IsMaster());
212   }
213 
214       //copy the data
215       nBins = theModel->nBins;
216 
217       //Same verbosity for all workers, as the master
218       fVerboseLevel = theModel->fVerboseLevel;
219     }
220   return;
221 }
222 
223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
224 
225 G4double G4PenelopeBremsstrahlungModel::CrossSectionPerVolume(const G4Material* material,
226                                            const G4ParticleDefinition* theParticle,
227                                            G4double energy,
228                                            G4double cutEnergy,
229                                            G4double)
230 {
231   //
232   if (fVerboseLevel > 3)
233     G4cout << "Calling CrossSectionPerVolume() of G4PenelopeBremsstrahlungModel" << G4endl;
234 
235   SetupForMaterial(theParticle, material, energy);
236   G4double crossPerMolecule = 0.;
237 
238   G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(theParticle,material,
239                                                                 cutEnergy);
240   if (theXS)
241     crossPerMolecule = theXS->GetHardCrossSection(energy);
242 
243   G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
244   G4double atPerMol =  fOscManager->GetAtomsPerMolecule(material);
245 
246   if (fVerboseLevel > 3)
247     G4cout << "Material " << material->GetName() << " has " << atPerMol <<
248       "atoms per molecule" << G4endl;
249 
250   G4double moleculeDensity = 0.;
251   if (atPerMol)
252     moleculeDensity = atomDensity/atPerMol;
253 
254   G4double crossPerVolume = crossPerMolecule*moleculeDensity;
255 
256   if (fVerboseLevel > 2)
257   {
258     G4cout << "G4PenelopeBremsstrahlungModel " << G4endl;
259     G4cout << "Mean free path for gamma emission > " << cutEnergy/keV << " keV at " <<
260       energy/keV << " keV = " << (1./crossPerVolume)/mm << " mm" << G4endl;
261   }
262 
263   return crossPerVolume;
264 }
265 
266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
267 
268 //This is a dummy method. Never inkoved by the tracking, it just issues
269 //a warning if one tries to get Cross Sections per Atom via the
270 //G4EmCalculator.
271 G4double G4PenelopeBremsstrahlungModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
272                      G4double,
273                      G4double,
274                      G4double,
275                      G4double,
276                      G4double)
277 {
278   G4cout << "*** G4PenelopeBremsstrahlungModel -- WARNING ***" << G4endl;
279   G4cout << "Penelope Bremsstrahlung model v2008 does not calculate cross section _per atom_ " << G4endl;
280   G4cout << "so the result is always zero. For physics values, please invoke " << G4endl;
281   G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl;
282   return 0;
283 }
284 
285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
286 
287 G4double G4PenelopeBremsstrahlungModel::ComputeDEDXPerVolume(const G4Material* material,
288                    const G4ParticleDefinition* theParticle,
289                    G4double kineticEnergy,
290                    G4double cutEnergy)
291 {
292   if (fVerboseLevel > 3)
293     G4cout << "Calling ComputeDEDX() of G4PenelopeBremsstrahlungModel" << G4endl;
294 
295   G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(theParticle,material,
296                                                                 cutEnergy);
297   G4double sPowerPerMolecule = 0.0;
298   if (theXS)
299     sPowerPerMolecule = theXS->GetSoftStoppingPower(kineticEnergy);
300 
301   G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
302   G4double atPerMol =  fOscManager->GetAtomsPerMolecule(material);
303 
304   G4double moleculeDensity = 0.;
305   if (atPerMol)
306     moleculeDensity = atomDensity/atPerMol;
307 
308   G4double sPowerPerVolume = sPowerPerMolecule*moleculeDensity;
309 
310   if (fVerboseLevel > 2)
311     {
312       G4cout << "G4PenelopeBremsstrahlungModel " << G4endl;
313       G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
314         kineticEnergy/keV << " keV = " <<
315         sPowerPerVolume/(keV/mm) << " keV/mm" << G4endl;
316     }
317   return sPowerPerVolume;
318 }
319 
320 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
321 
322 void G4PenelopeBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>*fvect,
323                   const G4MaterialCutsCouple* couple,
324                   const G4DynamicParticle* aDynamicParticle,
325                   G4double cutG,
326                   G4double)
327 {
328   if (fVerboseLevel > 3)
329     G4cout << "Calling SampleSecondaries() of G4PenelopeBremsstrahlungModel" << G4endl;
330 
331   G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
332   const G4Material* material = couple->GetMaterial();
333 
334   if (kineticEnergy <= fIntrinsicLowEnergyLimit)
335    {
336      fParticleChange->SetProposedKineticEnergy(0.);
337      fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
338      return ;
339    }
340 
341   G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
342   //This is the momentum
343   G4ThreeVector initialMomentum =  aDynamicParticle->GetMomentum();
344 
345   //Not enough energy to produce a secondary! Return with nothing happened
346   if (kineticEnergy < cutG)
347     return;
348 
349   if (fVerboseLevel > 3)
350     G4cout << "Going to sample gamma energy for: " <<material->GetName() << " " <<
351       "energy = " << kineticEnergy/keV << ", cut = " << cutG/keV << G4endl;
352 
353    //Sample gamma's energy according to the spectrum
354   G4double gammaEnergy =
355     fPenelopeFSHelper->SampleGammaEnergy(kineticEnergy,material,cutG);
356 
357   if (fVerboseLevel > 3)
358     G4cout << "Sampled gamma energy: " << gammaEnergy/keV << " keV" << G4endl;
359 
360   //Now sample the direction for the Gamma. Notice that the rotation is already done
361   //Z is unused here, I plug 0. The information is in the material pointer
362    G4ThreeVector gammaDirection1 =
363      fPenelopeAngular->SampleDirection(aDynamicParticle,gammaEnergy,0,material);
364 
365   if (fVerboseLevel > 3)
366     G4cout << "Sampled cosTheta for e-: " << gammaDirection1.cosTheta() << G4endl;
367 
368   G4double residualPrimaryEnergy = kineticEnergy-gammaEnergy;
369   if (residualPrimaryEnergy < 0)
370     {
371       //Ok we have a problem, all energy goes with the gamma
372       gammaEnergy += residualPrimaryEnergy;
373       residualPrimaryEnergy = 0.0;
374     }
375 
376   //Produce final state according to momentum conservation
377   G4ThreeVector particleDirection1 = initialMomentum - gammaEnergy*gammaDirection1;
378   particleDirection1 = particleDirection1.unit(); //normalize
379 
380   //Update the primary particle
381   if (residualPrimaryEnergy > 0.)
382     {
383       fParticleChange->ProposeMomentumDirection(particleDirection1);
384       fParticleChange->SetProposedKineticEnergy(residualPrimaryEnergy);
385     }
386   else
387     {
388       fParticleChange->SetProposedKineticEnergy(0.);
389     }
390 
391   //Now produce the photon
392   G4DynamicParticle* theGamma = new G4DynamicParticle(G4Gamma::Gamma(),
393                                                       gammaDirection1,
394                                                       gammaEnergy);
395   fvect->push_back(theGamma);
396 
397   if (fVerboseLevel > 1)
398     {
399       G4cout << "-----------------------------------------------------------" << G4endl;
400       G4cout << "Energy balance from G4PenelopeBremsstrahlung" << G4endl;
401       G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
402       G4cout << "-----------------------------------------------------------" << G4endl;
403       G4cout << "Outgoing primary energy: " << residualPrimaryEnergy/keV << " keV" << G4endl;
404       G4cout << "Bremsstrahlung photon " << gammaEnergy/keV << " keV" << G4endl;
405       G4cout << "Total final state: " << (residualPrimaryEnergy+gammaEnergy)/keV
406              << " keV" << G4endl;
407       G4cout << "-----------------------------------------------------------" << G4endl;
408     }
409 
410   if (fVerboseLevel > 0)
411     {
412       G4double energyDiff = std::fabs(residualPrimaryEnergy+gammaEnergy-kineticEnergy);
413       if (energyDiff > 0.05*keV)
414         G4cout << "Warning from G4PenelopeBremsstrahlung: problem with energy conservation: "
415          <<
416           (residualPrimaryEnergy+gammaEnergy)/keV <<
417           " keV (final) vs. " <<
418           kineticEnergy/keV << " keV (initial)" << G4endl;
419     }
420   return;
421 }
422 
423 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
424 
425 void G4PenelopeBremsstrahlungModel::ClearTables()
426 {
427   if (!IsMaster() && !fLocalTable)
428     //Should not be here!
429     G4Exception("G4PenelopeBremsstrahlungModel::ClearTables()",
430     "em0100",FatalException,"Worker thread in this method");
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: cleared tables" << G4endl;
455 
456  return;
457 }
458 
459 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
460 
461 G4double G4PenelopeBremsstrahlungModel::MinEnergyCut(const G4ParticleDefinition*,
462                    const G4MaterialCutsCouple*)
463 {
464   return fIntrinsicLowEnergyLimit;
465 }
466 
467 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
468 
469 void G4PenelopeBremsstrahlungModel::BuildXSTable(const G4Material* mat,G4double cut)
470 {
471   if (!IsMaster() && !fLocalTable)
472     //Should not be here!
473     G4Exception("G4PenelopeBremsstrahlungModel::BuildXSTable()",
474     "em0100",FatalException,"Worker thread in this method");
475 
476   //The key of the map
477   std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
478 
479   //The table already exists
480   if (fXSTableElectron->count(theKey) && fXSTablePositron->count(theKey))
481     return;
482 
483   //
484   //This method fills the G4PenelopeCrossSection containers for electrons or positrons
485   //and for the given material/cut couple.
486   //Equivalent of subroutines EBRaT and PINaT of Penelope
487   //
488   if (fVerboseLevel > 2)
489     {
490       G4cout << "G4PenelopeBremsstrahlungModel: going to build cross section table " << G4endl;
491       G4cout << "for e+/e- in " << mat->GetName() << " for Ecut(gamma)= " <<
492   cut/keV << " keV " << G4endl;
493     }
494 
495   //Tables have been already created (checked by GetCrossSectionTableForCouple)
496   if (fEnergyGrid->GetVectorLength() != nBins)
497     {
498       G4ExceptionDescription ed;
499       ed << "Energy Grid looks not initialized" << G4endl;
500       ed << nBins << " " << fEnergyGrid->GetVectorLength() << G4endl;
501       G4Exception("G4PenelopeBremsstrahlungModel::BuildXSTable()",
502       "em2016",FatalException,ed);
503     }
504 
505   G4PenelopeCrossSection* XSEntry = new G4PenelopeCrossSection(nBins);
506   G4PenelopeCrossSection* XSEntry2 = new G4PenelopeCrossSection(nBins);
507   const G4PhysicsTable* table = fPenelopeFSHelper->GetScaledXSTable(mat,cut);
508 
509   //loop on the energy grid
510   for (std::size_t bin=0;bin<nBins;++bin)
511     {
512        G4double energy = fEnergyGrid->GetLowEdgeEnergy(bin);
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->GetEffectiveZSquared(mat)*
518    ((energy+electron_mass_c2)*(energy+electron_mass_c2)/
519     (energy*(energy+2.0*electron_mass_c2)));
520 
521        G4double restrictedCut = cut/energy;
522 
523        //Now I need the dSigma/dX profile - interpolated on energy - for
524        //the 32-point x grid. Interpolation is log-log
525        std::size_t nBinsX = fPenelopeFSHelper->GetNBinsX();
526        G4double* tempData = new G4double[nBinsX];
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 belongs to the 32-point grid.
531      G4double val = (*table)[ix]->Value(logene);
532      tempData[ix] = G4Exp(val); //back to the real value!
533    }
534 
535        G4double XH0A = 0.;
536        if (restrictedCut <= 1) //calculate only if we are above threshold!
537    XH0A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,-1) -
538      fPenelopeFSHelper->GetMomentumIntegral(tempData,restrictedCut,-1);
539        G4double XS1A = fPenelopeFSHelper->GetMomentumIntegral(tempData,
540                     restrictedCut,0);
541        G4double XS2A = fPenelopeFSHelper->GetMomentumIntegral(tempData,
542                     restrictedCut,1);
543        G4double XH1A=0, XH2A=0;
544        if (restrictedCut <=1)
545    {
546      XH1A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,0) -
547        XS1A;
548      XH2A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,1) -
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,energy,XH0,XH1,XH2,XS0,XS1,XS2);
560 
561        //take care of positrons
562        G4double posCorrection = GetPositronXSCorrection(mat,energy);
563        XSEntry2->AddCrossSectionPoint(bin,energy,XH0*posCorrection,
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(theKey,XSEntry));
573   fXSTablePositron->insert(std::make_pair(theKey,XSEntry2));
574 
575   return;
576 }
577 
578 
579 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
580 
581 G4PenelopeCrossSection*
582 G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple(const G4ParticleDefinition* part,
583                    const G4Material* mat,
584                    G4double cut)
585 {
586   if (part != G4Electron::Electron() && part != G4Positron::Positron())
587     {
588       G4ExceptionDescription ed;
589       ed << "Invalid particle: " << part->GetParticleName() << G4endl;
590       G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
591       "em0001",FatalException,ed);
592       return nullptr;
593     }
594 
595   if (part == G4Electron::Electron())
596     {
597       //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
598       //not invoked
599       if (!fXSTableElectron)
600         {
601     //create a **thread-local** version of the table. Used only for G4EmCalculator and
602     //Unit Tests
603           G4String excep = "The Cross Section Table for e- was not initialized correctly!";
604           G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
605           "em2013",JustWarning,excep);
606     fLocalTable = true;
607     fXSTableElectron = new
608       std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
609     if (!fEnergyGrid)
610       fEnergyGrid = new G4PhysicsLogVector(LowEnergyLimit(),
611             HighEnergyLimit(),
612             nBins-1); //one hidden bin is added
613     if (!fPenelopeFSHelper)
614       fPenelopeFSHelper = new G4PenelopeBremsstrahlungFS(fVerboseLevel);
615         }
616       std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
617       if (fXSTableElectron->count(theKey)) //table already built
618         return fXSTableElectron->find(theKey)->second;
619       else
620   {
621     //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
622     //not filled up. This can happen in a UnitTest or via G4EmCalculator
623     if (fVerboseLevel > 0)
624       {
625         //G4Exception (warning) is issued only in verbose mode
626         G4ExceptionDescription ed;
627         ed << "Unable to find e- table for " << mat->GetName() << " at Ecut(gamma)= "
628      << cut/keV << " keV " << G4endl;
629         ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
630         G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
631         "em2009",JustWarning,ed);
632       }
633     //protect file reading via autolock
634     G4AutoLock lock(&PenelopeBremsstrahlungModelMutex);
635           fPenelopeFSHelper->BuildScaledXSTable(mat,cut,true); //pretend to be a master
636     BuildXSTable(mat,cut);
637     lock.unlock();
638     //now it should be ok
639     return fXSTableElectron->find(theKey)->second;
640   }
641     }
642   if (part == G4Positron::Positron())
643     {
644       //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
645       //not invoked
646       if (!fXSTablePositron)
647         {
648     G4String excep = "The Cross Section Table for e+ was not initialized correctly!";
649           G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
650           "em2013",JustWarning,excep);
651     fLocalTable = true;
652     fXSTablePositron = new
653       std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
654     if (!fEnergyGrid)
655       fEnergyGrid = new G4PhysicsLogVector(LowEnergyLimit(),
656             HighEnergyLimit(),
657             nBins-1); //one hidden bin is added
658     if (!fPenelopeFSHelper)
659             fPenelopeFSHelper = new G4PenelopeBremsstrahlungFS(fVerboseLevel);
660         }
661       std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
662       if (fXSTablePositron->count(theKey)) //table already built
663         return fXSTablePositron->find(theKey)->second;
664       else
665         {
666     //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
667     //not filled up. This can happen in a UnitTest or via G4EmCalculator
668     if (fVerboseLevel > 0)
669       {
670         //Issue a G4Exception (warning) only in verbose mode
671         G4ExceptionDescription ed;
672         ed << "Unable to find e+ table for " << mat->GetName() << " at Ecut(gamma)= "
673      << cut/keV << " keV " << G4endl;
674         ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
675         G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
676         "em2009",JustWarning,ed);
677       }
678     //protect file reading via autolock
679     G4AutoLock lock(&PenelopeBremsstrahlungModelMutex);
680           fPenelopeFSHelper->BuildScaledXSTable(mat,cut,true); //pretend to be a master
681     BuildXSTable(mat,cut);
682     lock.unlock();
683     //now it should be ok
684     return fXSTablePositron->find(theKey)->second;
685         }
686     }
687   return nullptr;
688 }
689 
690 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
691 
692 G4double G4PenelopeBremsstrahlungModel::GetPositronXSCorrection(const G4Material* mat,
693                   G4double energy)
694 {
695   //The electron-to-positron correction factor is set equal to the ratio of the
696   //radiative stopping powers for positrons and electrons, which has been calculated
697   //by Kim et al. (1986) (cf. Berger and Seltzer, 1982). Here, it is used an
698   //analytical approximation which reproduces the tabulated values with 0.5%
699   //accuracy
700   G4double t=G4Log(1.0+1e6*energy/
701           (electron_mass_c2*fPenelopeFSHelper->GetEffectiveZSquared(mat)));
702   G4double corr = 1.0-G4Exp(-t*(1.2359e-1-t*(6.1274e-2-t*
703              (3.1516e-2-t*(7.7446e-3-t*(1.0595e-3-t*
704                       (7.0568e-5-t*
705                        1.8080e-6)))))));
706   return corr;
707 }
708 
709 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
710 
711 void G4PenelopeBremsstrahlungModel::SetParticle(const G4ParticleDefinition* p)
712 {
713   if(!fParticle) {
714     fParticle = p;
715   }
716 }
717