Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4PenelopeCrossSection.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 // 18 Mar 2010   L Pandola    First implementation
 32 // 09 Mar 2012   L. Pandola   Add public method (and machinery) to return
 33 //                            the absolute and the normalized shell cross
 34 //                            sections independently.
 35 //
 36 #include "G4PenelopeCrossSection.hh"
 37 #include "G4SystemOfUnits.hh"
 38 #include "G4PhysicsTable.hh"
 39 #include "G4PhysicsFreeVector.hh"
 40 #include "G4DataVector.hh"
 41 #include "G4Exp.hh"
 42 #include "G4Log.hh"
 43 
 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
 45 G4PenelopeCrossSection::G4PenelopeCrossSection(size_t nPointsE,size_t nShells) :
 46   fSoftCrossSections(nullptr),
 47   fHardCrossSections(nullptr),fShellCrossSections(nullptr),
 48   fShellNormalizedCrossSections(nullptr),
 49   fNumberOfEnergyPoints(nPointsE),fNumberOfShells(nShells)
 50 {
 51   //check the number of points is not zero
 52   if (!fNumberOfEnergyPoints)
 53     {
 54       G4ExceptionDescription ed;
 55       ed << "G4PenelopeCrossSection: invalid number of energy points " << G4endl;
 56       G4Exception("G4PenelopeCrossSection::G4PenelopeCrossSection()",
 57       "em2017",FatalException,ed);
 58     }
 59 
 60   fIsNormalized = false;
 61 
 62   // 1) soft XS table
 63   fSoftCrossSections = new G4PhysicsTable();
 64   //the table contains 3 G4PhysicsFreeVectors,
 65   //(fSoftCrossSections)[0] -->  log XS0 vs. log E
 66   //(fSoftCrossSections)[1] -->  log XS1 vs. log E
 67   //(fSoftCrossSections)[2] -->  log XS2 vs. log E
 68 
 69   //everything is log-log
 70   for (size_t i=0;i<3;i++)
 71     fSoftCrossSections->push_back(new G4PhysicsFreeVector(fNumberOfEnergyPoints));
 72 
 73   //2) hard XS table
 74   fHardCrossSections = new G4PhysicsTable();
 75   //the table contains 3 G4PhysicsFreeVectors,
 76   //(fHardCrossSections)[0] -->  log XH0 vs. log E
 77   //(fHardCrossSections)[1] -->  log XH1 vs. log E
 78   //(fHardCrossSections)[2] -->  log XH2 vs. log E
 79 
 80   //everything is log-log
 81   for (size_t i=0;i<3;i++)
 82     fHardCrossSections->push_back(new G4PhysicsFreeVector(fNumberOfEnergyPoints));
 83 
 84   //3) shell XS table, if it is the case
 85   if (fNumberOfShells)
 86     {
 87       fShellCrossSections = new G4PhysicsTable();
 88       fShellNormalizedCrossSections = new G4PhysicsTable();
 89       //the table has to contain numberofShells G4PhysicsFreeVectors,
 90       //(theTable)[ishell] --> cross section for shell #ishell
 91       for (size_t i=0;i<fNumberOfShells;i++)
 92   {
 93     fShellCrossSections->push_back(new G4PhysicsFreeVector(fNumberOfEnergyPoints));
 94     fShellNormalizedCrossSections->push_back(new G4PhysicsFreeVector(fNumberOfEnergyPoints));
 95   }
 96     }
 97 }
 98 
 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
100 G4PenelopeCrossSection::~G4PenelopeCrossSection()
101 {
102   //clean up tables
103   if (fShellCrossSections)
104     {
105       fShellCrossSections->clearAndDestroy();
106       delete fShellCrossSections;
107     }
108   if (fShellNormalizedCrossSections)
109     {
110       fShellNormalizedCrossSections->clearAndDestroy();
111       delete fShellNormalizedCrossSections;
112     }
113   if (fSoftCrossSections)
114     {
115       fSoftCrossSections->clearAndDestroy();
116       delete fSoftCrossSections;
117     }
118   if (fHardCrossSections)
119     {
120       fHardCrossSections->clearAndDestroy();
121       delete fHardCrossSections;
122     }
123 }
124 
125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
126 void G4PenelopeCrossSection::AddCrossSectionPoint(size_t binNumber,G4double energy,
127               G4double XH0,
128               G4double XH1, G4double XH2,
129               G4double XS0, G4double XS1,
130               G4double XS2)
131 {
132   if (!fSoftCrossSections || !fHardCrossSections)
133     {
134       G4cout << "Something wrong in G4PenelopeCrossSection::AddCrossSectionPoint" <<
135   G4endl;
136       G4cout << "Trying to fill un-initialized tables" << G4endl;
137       return;
138     }
139 
140   //fill vectors
141   G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*fSoftCrossSections)[0];
142 
143   if (binNumber >= fNumberOfEnergyPoints)
144     {
145       G4cout << "Something wrong in G4PenelopeCrossSection::AddCrossSectionPoint" <<
146   G4endl;
147       G4cout << "Trying to register more points than originally declared" << G4endl;
148       return;
149     }
150    G4double logEne = G4Log(energy);
151 
152    //XS0
153    G4double val = G4Log(std::max(XS0,1e-42*cm2)); //avoid log(0)
154    theVector->PutValues(binNumber,logEne,val);
155 
156    //XS1
157    theVector = (G4PhysicsFreeVector*) (*fSoftCrossSections)[1];
158    val =  G4Log(std::max(XS1,1e-42*eV*cm2)); //avoid log(0)
159    theVector->PutValues(binNumber,logEne,val);
160 
161    //XS2
162    theVector = (G4PhysicsFreeVector*) (*fSoftCrossSections)[2];
163    val =  G4Log(std::max(XS2,1e-42*eV*eV*cm2)); //avoid log(0)
164    theVector->PutValues(binNumber,logEne,val);
165 
166    //XH0
167    theVector = (G4PhysicsFreeVector*) (*fHardCrossSections)[0];
168    val =  G4Log(std::max(XH0,1e-42*cm2)); //avoid log(0)
169    theVector->PutValues(binNumber,logEne,val);
170 
171    //XH1
172    theVector = (G4PhysicsFreeVector*) (*fHardCrossSections)[1];
173    val =  G4Log(std::max(XH1,1e-42*eV*cm2)); //avoid log(0)
174    theVector->PutValues(binNumber,logEne,val);
175 
176     //XH2
177    theVector = (G4PhysicsFreeVector*) (*fHardCrossSections)[2];
178    val =  G4Log(std::max(XH2,1e-42*eV*eV*cm2)); //avoid log(0)
179    theVector->PutValues(binNumber,logEne,val);
180 
181    return;
182 }
183 
184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
185 
186 void G4PenelopeCrossSection::AddShellCrossSectionPoint(size_t binNumber,
187                    size_t shellID,
188                    G4double energy,
189                    G4double xs)
190 {
191   if (!fShellCrossSections)
192     {
193       G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
194   G4endl;
195       G4cout << "Trying to fill un-initialized table" << G4endl;
196       return;
197     }
198 
199   if (shellID >= fNumberOfShells)
200     {
201       G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
202   G4endl;
203       G4cout << "Trying to fill shell #" << shellID << " while the maximum is "
204        <<  fNumberOfShells-1 << G4endl;
205       return;
206     }
207 
208   //fill vector
209   G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*fShellCrossSections)[shellID];
210 
211   if (binNumber >= fNumberOfEnergyPoints)
212     {
213       G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
214   G4endl;
215       G4cout << "Trying to register more points than originally declared" << G4endl;
216       return;
217     }
218    G4double logEne = G4Log(energy);
219    G4double val = G4Log(std::max(xs,1e-42*cm2)); //avoid log(0)
220    theVector->PutValues(binNumber,logEne,val);
221 
222    return;
223 }
224 
225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
226 
227 G4double G4PenelopeCrossSection::GetTotalCrossSection(G4double energy) const
228 {
229   G4double result = 0;
230   //take here XS0 + XH0
231   if (!fSoftCrossSections || !fHardCrossSections)
232     {
233       G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
234   G4endl;
235       G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
236       return result;
237     }
238 
239   // 1) soft part
240   G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*fSoftCrossSections)[0];
241   if (theVector->GetVectorLength() < fNumberOfEnergyPoints)
242     {
243       G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
244   G4endl;
245       G4cout << "Soft cross section table looks not filled" << G4endl;
246       return result;
247     }
248   G4double logene = G4Log(energy);
249   G4double logXS = theVector->Value(logene);
250   G4double softXS = G4Exp(logXS);
251 
252    // 2) hard part
253   theVector = (G4PhysicsFreeVector*) (*fHardCrossSections)[0];
254   if (theVector->GetVectorLength() < fNumberOfEnergyPoints)
255     {
256       G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
257   G4endl;
258       G4cout << "Hard cross section table looks not filled" << G4endl;
259       return result;
260     }
261   logXS = theVector->Value(logene);
262   G4double hardXS = G4Exp(logXS);
263 
264   result = hardXS + softXS;
265   return result;
266 }
267 
268 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
269 
270 G4double G4PenelopeCrossSection::GetHardCrossSection(G4double energy) const
271 {
272   G4double result = 0;
273   //take here XH0
274   if (!fHardCrossSections)
275     {
276       G4cout << "Something wrong in G4PenelopeCrossSection::GetHardCrossSection" <<
277   G4endl;
278       G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
279       return result;
280     }
281 
282   G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*fHardCrossSections)[0];
283   if (theVector->GetVectorLength() < fNumberOfEnergyPoints)
284     {
285       G4cout << "Something wrong in G4PenelopeCrossSection::GetHardCrossSection" <<
286   G4endl;
287       G4cout << "Hard cross section table looks not filled" << G4endl;
288       return result;
289     }
290   G4double logene = G4Log(energy);
291   G4double logXS = theVector->Value(logene);
292   result = G4Exp(logXS);
293 
294   return result;
295 }
296 
297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
298 
299 G4double G4PenelopeCrossSection::GetSoftStoppingPower(G4double energy) const
300 {
301   G4double result = 0;
302   //take here XH0
303   if (!fSoftCrossSections)
304     {
305       G4cout << "Something wrong in G4PenelopeCrossSection::GetSoftStoppingPower" <<
306   G4endl;
307       G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
308       return result;
309     }
310 
311   G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*fSoftCrossSections)[1];
312   if (theVector->GetVectorLength() < fNumberOfEnergyPoints)
313     {
314       G4cout << "Something wrong in G4PenelopeCrossSection::GetSoftStoppingPower" <<
315   G4endl;
316       G4cout << "Soft cross section table looks not filled" << G4endl;
317       return result;
318     }
319   G4double logene = G4Log(energy);
320   G4double logXS = theVector->Value(logene);
321   result = G4Exp(logXS);
322 
323   return result;
324 }
325 
326 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
327 
328 G4double G4PenelopeCrossSection::GetShellCrossSection(size_t shellID,G4double energy) const
329 {
330   G4double result = 0;
331   if (!fShellCrossSections)
332     {
333       G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
334   G4endl;
335       G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
336       return result;
337     }
338   if (shellID >= fNumberOfShells)
339     {
340       G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
341   G4endl;
342       G4cout << "Trying to retrieve shell #" << shellID << " while the maximum is "
343        <<  fNumberOfShells-1 << G4endl;
344       return result;
345     }
346 
347   G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*fShellCrossSections)[shellID];
348 
349   if (theVector->GetVectorLength() < fNumberOfEnergyPoints)
350     {
351       G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
352   G4endl;
353       G4cout << "Shell cross section table looks not filled" << G4endl;
354       return result;
355     }
356   G4double logene = G4Log(energy);
357   G4double logXS = theVector->Value(logene);
358   result = G4Exp(logXS);
359 
360   return result;
361 }
362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
363 
364 G4double G4PenelopeCrossSection::GetNormalizedShellCrossSection(size_t shellID,G4double energy) const
365 {
366   G4double result = 0;
367   if (!fShellNormalizedCrossSections)
368     {
369       G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
370   G4endl;
371       G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
372       return result;
373     }
374 
375   if (!fIsNormalized)
376     {
377       G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" << G4endl;
378       G4cout << "The table of normalized cross section is not initialized" << G4endl;
379     }
380 
381   if (shellID >= fNumberOfShells)
382     {
383       G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
384   G4endl;
385       G4cout << "Trying to retrieve shell #" << shellID << " while the maximum is "
386        <<  fNumberOfShells-1 << G4endl;
387       return result;
388     }
389 
390   const G4PhysicsFreeVector* theVector =
391     (G4PhysicsFreeVector*) (*fShellNormalizedCrossSections)[shellID];
392 
393   if (theVector->GetVectorLength() < fNumberOfEnergyPoints)
394     {
395       G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
396   G4endl;
397       G4cout << "Shell cross section table looks not filled" << G4endl;
398       return result;
399     }
400   G4double logene = G4Log(energy);
401   G4double logXS = theVector->Value(logene);
402   result = G4Exp(logXS);
403 
404   return result;
405 }
406 
407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
408 
409 void G4PenelopeCrossSection::NormalizeShellCrossSections()
410 {
411   if (fIsNormalized) //already done!
412     {
413       G4cout << "G4PenelopeCrossSection::NormalizeShellCrossSections()" << G4endl;
414       G4cout << "already invoked. Ignore it" << G4endl;
415       return;
416     }
417 
418   if (!fShellNormalizedCrossSections)
419     {
420       G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
421   G4endl;
422       G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
423       return;
424     }
425 
426   for (size_t i=0;i<fNumberOfEnergyPoints;i++) //loop on energy
427     {
428       //energy grid is the same for all shells
429 
430       //Recalculate manually the XS factor, to avoid problems with
431       //underflows
432       G4double normFactor = 0.;
433       for (size_t shellID=0;shellID<fNumberOfShells;shellID++)
434   {
435     G4PhysicsFreeVector* theVec =
436       (G4PhysicsFreeVector*) (*fShellCrossSections)[shellID];
437 
438     normFactor += G4Exp((*theVec)[i]);
439   }
440       G4double logNormFactor = G4Log(normFactor);
441       //Normalize
442       for (size_t shellID=0;shellID<fNumberOfShells;shellID++)
443   {
444    G4PhysicsFreeVector* theVec =
445       (G4PhysicsFreeVector*) (*fShellNormalizedCrossSections)[shellID];
446    G4PhysicsFreeVector* theFullVec =
447      (G4PhysicsFreeVector*) (*fShellCrossSections)[shellID];
448    G4double previousValue = (*theFullVec)[i]; //log(XS)
449    G4double logEnergy = theFullVec->GetLowEdgeEnergy(i);
450    //log(XS/normFactor) = log(XS) - log(normFactor)
451    theVec->PutValues(i,logEnergy,previousValue-logNormFactor);
452   }
453     }
454   fIsNormalized = true;
455 
456   return;
457 }
458