ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PenelopeCrossSection.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PenelopeCrossSection.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 //
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 //
37 #include "G4SystemOfUnits.hh"
38 #include "G4PhysicsTable.hh"
39 #include "G4PhysicsFreeVector.hh"
40 #include "G4DataVector.hh"
41 #include "G4Exp.hh"
42 
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
44 G4PenelopeCrossSection::G4PenelopeCrossSection(size_t nPointsE,size_t nShells) :
45  numberOfEnergyPoints(nPointsE),numberOfShells(nShells),softCrossSections(0),
46  hardCrossSections(0),shellCrossSections(0),shellNormalizedCrossSections(0)
47 {
48  //check the number of points is not zero
50  {
52  ed << "G4PenelopeCrossSection: invalid number of energy points " << G4endl;
53  G4Exception("G4PenelopeCrossSection::G4PenelopeCrossSection()",
54  "em2017",FatalException,ed);
55  }
56 
57  isNormalized = false;
58 
59  // 1) soft XS table
61  //the table contains 3 G4PhysicsFreeVectors,
62  //(softCrossSections)[0] --> log XS0 vs. log E
63  //(softCrossSections)[1] --> log XS1 vs. log E
64  //(softCrossSections)[2] --> log XS2 vs. log E
65 
66  //everything is log-log
67  for (size_t i=0;i<3;i++)
69 
70  //2) hard XS table
72  //the table contains 3 G4PhysicsFreeVectors,
73  //(hardCrossSections)[0] --> log XH0 vs. log E
74  //(hardCrossSections)[1] --> log XH1 vs. log E
75  //(hardCrossSections)[2] --> log XH2 vs. log E
76 
77  //everything is log-log
78  for (size_t i=0;i<3;i++)
80 
81  //3) shell XS table, if it is the case
82  if (numberOfShells)
83  {
86  //the table has to contain numberofShells G4PhysicsFreeVectors,
87  //(theTable)[ishell] --> cross section for shell #ishell
88  for (size_t i=0;i<numberOfShells;i++)
89  {
92  }
93  }
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
98 {
99  //clean up tables
100  if (shellCrossSections)
101  {
102  //shellCrossSections->clearAndDestroy();
103  delete shellCrossSections;
104  }
106  {
107  //shellNormalizedCrossSections->clearAndDestroy();
109  }
110  if (softCrossSections)
111  {
112  //softCrossSections->clearAndDestroy();
113  delete softCrossSections;
114  }
115  if (hardCrossSections)
116  {
117  //hardCrossSections->clearAndDestroy();
118  delete hardCrossSections;
119  }
120 }
121 
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
124  G4double XH0,
125  G4double XH1, G4double XH2,
126  G4double XS0, G4double XS1,
127  G4double XS2)
128 {
130  {
131  G4cout << "Something wrong in G4PenelopeCrossSection::AddCrossSectionPoint" <<
132  G4endl;
133  G4cout << "Trying to fill un-initialized tables" << G4endl;
134  return;
135  }
136 
137  //fill vectors
139 
140  if (binNumber >= numberOfEnergyPoints)
141  {
142  G4cout << "Something wrong in G4PenelopeCrossSection::AddCrossSectionPoint" <<
143  G4endl;
144  G4cout << "Trying to register more points than originally declared" << G4endl;
145  return;
146  }
147  G4double logEne = std::log(energy);
148 
149  //XS0
150  G4double val = std::log(std::max(XS0,1e-42*cm2)); //avoid log(0)
151  theVector->PutValue(binNumber,logEne,val);
152 
153  //XS1
154  theVector = (G4PhysicsFreeVector*) (*softCrossSections)[1];
155  val = std::log(std::max(XS1,1e-42*eV*cm2)); //avoid log(0)
156  theVector->PutValue(binNumber,logEne,val);
157 
158  //XS2
159  theVector = (G4PhysicsFreeVector*) (*softCrossSections)[2];
160  val = std::log(std::max(XS2,1e-42*eV*eV*cm2)); //avoid log(0)
161  theVector->PutValue(binNumber,logEne,val);
162 
163  //XH0
164  theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
165  val = std::log(std::max(XH0,1e-42*cm2)); //avoid log(0)
166  theVector->PutValue(binNumber,logEne,val);
167 
168  //XH1
169  theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[1];
170  val = std::log(std::max(XH1,1e-42*eV*cm2)); //avoid log(0)
171  theVector->PutValue(binNumber,logEne,val);
172 
173  //XH2
174  theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[2];
175  val = std::log(std::max(XH2,1e-42*eV*eV*cm2)); //avoid log(0)
176  theVector->PutValue(binNumber,logEne,val);
177 
178  return;
179 }
180 
181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
182 
184  size_t shellID,
186  G4double xs)
187 {
188  if (!shellCrossSections)
189  {
190  G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
191  G4endl;
192  G4cout << "Trying to fill un-initialized table" << G4endl;
193  return;
194  }
195 
196  if (shellID >= numberOfShells)
197  {
198  G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
199  G4endl;
200  G4cout << "Trying to fill shell #" << shellID << " while the maximum is "
201  << numberOfShells-1 << G4endl;
202  return;
203  }
204 
205  //fill vector
206  G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
207 
208  if (binNumber >= numberOfEnergyPoints)
209  {
210  G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
211  G4endl;
212  G4cout << "Trying to register more points than originally declared" << G4endl;
213  return;
214  }
215  G4double logEne = std::log(energy);
216  G4double val = std::log(std::max(xs,1e-42*cm2)); //avoid log(0)
217  theVector->PutValue(binNumber,logEne,val);
218 
219  return;
220 }
221 
222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
223 
225 {
226  G4double result = 0;
227  //take here XS0 + XH0
229  {
230  G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
231  G4endl;
232  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
233  return result;
234  }
235 
236  // 1) soft part
238  if (theVector->GetVectorLength() < numberOfEnergyPoints)
239  {
240  G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
241  G4endl;
242  G4cout << "Soft cross section table looks not filled" << G4endl;
243  return result;
244  }
245  G4double logene = std::log(energy);
246  G4double logXS = theVector->Value(logene);
247  G4double softXS = G4Exp(logXS);
248 
249  // 2) hard part
250  theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
251  if (theVector->GetVectorLength() < numberOfEnergyPoints)
252  {
253  G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
254  G4endl;
255  G4cout << "Hard cross section table looks not filled" << G4endl;
256  return result;
257  }
258  logXS = theVector->Value(logene);
259  G4double hardXS = G4Exp(logXS);
260 
261  result = hardXS + softXS;
262  return result;
263 
264 }
265 
266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
267 
269 {
270  G4double result = 0;
271  //take here XH0
272  if (!hardCrossSections)
273  {
274  G4cout << "Something wrong in G4PenelopeCrossSection::GetHardCrossSection" <<
275  G4endl;
276  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
277  return result;
278  }
279 
281  if (theVector->GetVectorLength() < numberOfEnergyPoints)
282  {
283  G4cout << "Something wrong in G4PenelopeCrossSection::GetHardCrossSection" <<
284  G4endl;
285  G4cout << "Hard cross section table looks not filled" << G4endl;
286  return result;
287  }
288  G4double logene = std::log(energy);
289  G4double logXS = theVector->Value(logene);
290  result = G4Exp(logXS);
291 
292  return result;
293 }
294 
295 
296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
297 
299 {
300  G4double result = 0;
301  //take here XH0
302  if (!softCrossSections)
303  {
304  G4cout << "Something wrong in G4PenelopeCrossSection::GetSoftStoppingPower" <<
305  G4endl;
306  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
307  return result;
308  }
309 
311  if (theVector->GetVectorLength() < numberOfEnergyPoints)
312  {
313  G4cout << "Something wrong in G4PenelopeCrossSection::GetSoftStoppingPower" <<
314  G4endl;
315  G4cout << "Soft cross section table looks not filled" << G4endl;
316  return result;
317  }
318  G4double logene = std::log(energy);
319  G4double logXS = theVector->Value(logene);
320  result = G4Exp(logXS);
321 
322  return result;
323 }
324 
325 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
326 
328 {
329  G4double result = 0;
330  if (!shellCrossSections)
331  {
332  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
333  G4endl;
334  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
335  return result;
336  }
337  if (shellID >= numberOfShells)
338  {
339  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
340  G4endl;
341  G4cout << "Trying to retrieve shell #" << shellID << " while the maximum is "
342  << numberOfShells-1 << G4endl;
343  return result;
344  }
345 
346  G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
347 
348  if (theVector->GetVectorLength() < numberOfEnergyPoints)
349  {
350  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
351  G4endl;
352  G4cout << "Shell cross section table looks not filled" << G4endl;
353  return result;
354  }
355  G4double logene = std::log(energy);
356  G4double logXS = theVector->Value(logene);
357  result = G4Exp(logXS);
358 
359  return result;
360 }
361 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
362 
364 {
365  G4double result = 0;
367  {
368  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
369  G4endl;
370  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
371  return result;
372  }
373 
374  if (!isNormalized)
375  {
376  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" << G4endl;
377  G4cout << "The table of normalized cross section is not initialized" << G4endl;
378  }
379 
380 
381  if (shellID >= numberOfShells)
382  {
383  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
384  G4endl;
385  G4cout << "Trying to retrieve shell #" << shellID << " while the maximum is "
386  << numberOfShells-1 << G4endl;
387  return result;
388  }
389 
390  const G4PhysicsFreeVector* theVector =
392 
393  if (theVector->GetVectorLength() < numberOfEnergyPoints)
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 = std::log(energy);
401  G4double logXS = theVector->Value(logene);
402  result = G4Exp(logXS);
403 
404  return result;
405 }
406 
407 
408 
409 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
410 
412 {
413  if (isNormalized) //already done!
414  {
415  G4cout << "G4PenelopeCrossSection::NormalizeShellCrossSections()" << G4endl;
416  G4cout << "already invoked. Ignore it" << G4endl;
417  return;
418  }
419 
421  {
422  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
423  G4endl;
424  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
425  return;
426  }
427 
428  for (size_t i=0;i<numberOfEnergyPoints;i++) //loop on energy
429  {
430  //energy grid is the same for all shells
431 
432  //Recalculate manually the XS factor, to avoid problems with
433  //underflows
434  G4double normFactor = 0.;
435  for (size_t shellID=0;shellID<numberOfShells;shellID++)
436  {
437  G4PhysicsFreeVector* theVec =
439 
440  normFactor += G4Exp((*theVec)[i]);
441  }
442  G4double logNormFactor = std::log(normFactor);
443  //Normalize
444  for (size_t shellID=0;shellID<numberOfShells;shellID++)
445  {
446  G4PhysicsFreeVector* theVec =
448  G4PhysicsFreeVector* theFullVec =
450  G4double previousValue = (*theFullVec)[i]; //log(XS)
451  G4double logEnergy = theFullVec->GetLowEdgeEnergy(i);
452  //log(XS/normFactor) = log(XS) - log(normFactor)
453  theVec->PutValue(i,logEnergy,previousValue-logNormFactor);
454  }
455  }
456 
457  isNormalized = true;
458 
459 
460  /*
461  //TESTING
462  for (size_t shellID=0;shellID<numberOfShells;shellID++)
463  {
464  G4cout << "SHELL " << shellID << G4endl;
465  G4PhysicsFreeVector* theVec =
466  (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
467  for (size_t i=0;i<numberOfEnergyPoints;i++) //loop on energy
468  {
469  G4double logene = theVec->GetLowEdgeEnergy(i);
470  G4cout << G4Exp(logene)/MeV << " " << G4Exp((*theVec)[i]) << G4endl;
471  }
472  }
473  */
474 
475  return;
476 }