GeoModel - GeoModelCore
GeoMaterial.cxx
1 /*
2  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "GeoModelKernel/GeoMaterial.h"
6 #include <stdexcept>
7 #include <cmath>
8 #include <numeric>
9 #include <cfloat>
10 #include <iostream>
11 #include <algorithm>
12 #include <iterator>
13 
14 // These constants are the ionization potentials, indexed by atomic
15 // number. They have been obtained from
16 //
17 // http://physics.nist.gov/PhysRefData/XrayMassCoef/tab1.html
18 //
19 // Which is the reference of the particle data book. Approximate
20 // semiempirical formulae are not accurate enough, so we take the
21 // tabular values. These are in electron-volts.
22 
23 
24 const double
25  GeoMaterial::s_ionizationPotential[93] = {
26  0.0 *
27  GeoModelKernelUnits::eV,
28  19.2 *
29  GeoModelKernelUnits::eV,
30  41.8 *
31  GeoModelKernelUnits::eV,
32  40.0 *
33  GeoModelKernelUnits::eV,
34  63.7 *
35  GeoModelKernelUnits::eV,
36  76.0 *
37  GeoModelKernelUnits::eV,
38  78.0 *
39  GeoModelKernelUnits::eV,
40  82.0 *
41  GeoModelKernelUnits::eV,
42  95.0 *
43  GeoModelKernelUnits::eV,
44  115.0 *
45  GeoModelKernelUnits::eV,
46  137.0 *
47  GeoModelKernelUnits::eV,
48  149.0 *
49  GeoModelKernelUnits::eV,
50  156.0 *
51  GeoModelKernelUnits::eV,
52  166.0 *
53  GeoModelKernelUnits::eV,
54  173.0 *
55  GeoModelKernelUnits::eV,
56  173.0 *
57  GeoModelKernelUnits::eV,
58  180.0 *
59  GeoModelKernelUnits::eV,
60  174.0 *
61  GeoModelKernelUnits::eV,
62  188.0 *
63  GeoModelKernelUnits::eV,
64  190.0 *
65  GeoModelKernelUnits::eV,
66  191.0 *
67  GeoModelKernelUnits::eV,
68  216.0 *
69  GeoModelKernelUnits::eV,
70  233.0 *
71  GeoModelKernelUnits::eV,
72  245.0 *
73  GeoModelKernelUnits::eV,
74  257.0 *
75  GeoModelKernelUnits::eV,
76  272.0 *
77  GeoModelKernelUnits::eV,
78  286.0 *
79  GeoModelKernelUnits::eV,
80  297.0 *
81  GeoModelKernelUnits::eV,
82  311.0 *
83  GeoModelKernelUnits::eV,
84  322.0 *
85  GeoModelKernelUnits::eV,
86  330.0 *
87  GeoModelKernelUnits::eV,
88  334.0 *
89  GeoModelKernelUnits::eV,
90  350.0 *
91  GeoModelKernelUnits::eV,
92  347.0 *
93  GeoModelKernelUnits::eV,
94  348.0 *
95  GeoModelKernelUnits::eV,
96  343.0 *
97  GeoModelKernelUnits::eV,
98  352.0 *
99  GeoModelKernelUnits::eV,
100  363.0 *
101  GeoModelKernelUnits::eV,
102  366.0 *
103  GeoModelKernelUnits::eV,
104  379.0 *
105  GeoModelKernelUnits::eV,
106  393.0 *
107  GeoModelKernelUnits::eV,
108  417.0 *
109  GeoModelKernelUnits::eV,
110  424.0 *
111  GeoModelKernelUnits::eV,
112  428.0 *
113  GeoModelKernelUnits::eV,
114  441.0 *
115  GeoModelKernelUnits::eV,
116  449.0 *
117  GeoModelKernelUnits::eV,
118  470.0 *
119  GeoModelKernelUnits::eV,
120  470.0 *
121  GeoModelKernelUnits::eV,
122  469.0 *
123  GeoModelKernelUnits::eV,
124  488.0 *
125  GeoModelKernelUnits::eV,
126  488.0 *
127  GeoModelKernelUnits::eV,
128  487.0 *
129  GeoModelKernelUnits::eV,
130  485.0 *
131  GeoModelKernelUnits::eV,
132  491.0 *
133  GeoModelKernelUnits::eV,
134  482.0 *
135  GeoModelKernelUnits::eV,
136  488.0 *
137  GeoModelKernelUnits::eV,
138  491.0 *
139  GeoModelKernelUnits::eV,
140  501.0 *
141  GeoModelKernelUnits::eV,
142  523.0 *
143  GeoModelKernelUnits::eV,
144  535.0 *
145  GeoModelKernelUnits::eV,
146  546.0 *
147  GeoModelKernelUnits::eV,
148  560.0 *
149  GeoModelKernelUnits::eV,
150  574.0 *
151  GeoModelKernelUnits::eV,
152  580.0 *
153  GeoModelKernelUnits::eV,
154  591.0 *
155  GeoModelKernelUnits::eV,
156  614.0 *
157  GeoModelKernelUnits::eV,
158  628.0 *
159  GeoModelKernelUnits::eV,
160  650.0 *
161  GeoModelKernelUnits::eV,
162  658.0 *
163  GeoModelKernelUnits::eV,
164  674.0 *
165  GeoModelKernelUnits::eV,
166  684.0 *
167  GeoModelKernelUnits::eV,
168  694.0 *
169  GeoModelKernelUnits::eV,
170  705.0 *
171  GeoModelKernelUnits::eV,
172  718.0 *
173  GeoModelKernelUnits::eV,
174  727.0 *
175  GeoModelKernelUnits::eV,
176  736.0 *
177  GeoModelKernelUnits::eV,
178  746.0 *
179  GeoModelKernelUnits::eV,
180  757.0 *
181  GeoModelKernelUnits::eV,
182  790.0 *
183  GeoModelKernelUnits::eV,
184  790.0 *
185  GeoModelKernelUnits::eV,
186  800.0 *
187  GeoModelKernelUnits::eV,
188  810.0 *
189  GeoModelKernelUnits::eV,
190  823.0 *
191  GeoModelKernelUnits::eV,
192  823.0 *
193  GeoModelKernelUnits::eV,
194  830.0 *
195  GeoModelKernelUnits::eV,
196  825.0 *
197  GeoModelKernelUnits::eV,
198  794.0 *
199  GeoModelKernelUnits::eV,
200  827.0 *
201  GeoModelKernelUnits::eV,
202  826.0 *
203  GeoModelKernelUnits::eV,
204  841.0 *
205  GeoModelKernelUnits::eV,
206  847.0 *
207  GeoModelKernelUnits::eV,
208  878.0 *
209  GeoModelKernelUnits::eV,
210  890.0 *
211  GeoModelKernelUnits::eV
212 };
213 
214 unsigned int GeoMaterial::s_lastID = 0;
215 
216 GeoMaterial::GeoMaterial (const std::string &Name, double Density)
217  : m_name(Name)
218  , m_density(Density)
219  , m_iD(s_lastID++)
220  , m_radLength(0)
221  , m_intLength(0)
222  , m_dedDxConst(0)
223  , m_deDxI0(0)
224  , m_locked(false)
225 {
226 }
227 
228 GeoMaterial::~GeoMaterial()
229 {
230  for (size_t i = 0; i < m_element.size (); i++)
231  {
232  m_element[i]->unref ();
233  }
234 }
235 
236 void GeoMaterial::add (const GeoElement* element, double fraction)
237 {
238  // You can only add materials until you call "lock"...
239  if (!m_locked)
240  {
241  std::vector <const GeoElement *>::iterator e = std::find(m_element.begin(),m_element.end(),element);
242  if (e==m_element.end()) {
243  m_element.push_back (element);
244  m_fraction.push_back (fraction);
245  element->ref ();
246  }
247  else {
248  int n = e-m_element.begin();
249  m_fraction[n]+=fraction;
250  }
251  }
252  else
253  {
254  throw std::out_of_range ("Element added after material locked");
255  }
256 }
257 
258 void GeoMaterial::add (const GeoMaterial* material, double fraction)
259 {
260  if (!m_locked)
261  {
262  for (size_t e = 0; e < material->getNumElements (); e++)
263  {
264  add(material->m_element[e],fraction * material->m_fraction[e]);
265  }
266  }
267  else
268  {
269  throw std::out_of_range ("Material added after material locked");
270  }
271 }
272 
273 void GeoMaterial::lock ()
274 {
275  if(m_locked) return;
276 
277  m_locked = true;
278 
279  // -------------------------------------------//
280  // Now compute some quantities: //
281  // Source of these calculations is: //
282  // //
283  // For DeDx constants and X0: PDG. //
284  // For Lambda0: G4. //
285  // //
286  // For energy loss //
287  // //
288  const double C0 = 0.00307 * GeoModelKernelUnits::cm3 / GeoModelKernelUnits::gram; //
289  // //
290  // For nuclear absorption length. //
291  const double lambda0 = 35 * GeoModelKernelUnits::gram / GeoModelKernelUnits::cm2; //
292  // //
293  //--------------------------------------------//
294 
295  if (getNumElements () == 0)
296  {
297  throw std::out_of_range ("Attempt to lock a material with no elements");
298  return;
299  }
300 
301 
302  //--------------------------------------------//
303  // //
304  // -------------------------------------------//
305 
306  double dEDxConstant = 0, dEDxI0 = 0, NILinv = 0.0, radInv = 0.0;
307 
308  { // ===============Renormalization================================
309  double wSum=std::accumulate(m_fraction.begin(),m_fraction.end(),0.0);
310  if (fabs(wSum-1.0)>FLT_EPSILON) // 'FLT_EPSILON' defined in <cfloat>
311  {
312  if(std::getenv("GEOMODELKERNEL_VERBOSE")) {
313  std::cerr << "Warning in material "
314  << m_name
315  << ". Mass fractions sum to "
316  << wSum << "; renormalizing to 1.0" << std::endl;
317  }
318  }
319  double inv_wSum = 1. / wSum;
320  for (size_t e=0;e<getNumElements();e++) {m_fraction[e]*=inv_wSum;}
321  } // ==============================================================
322 
323  const double inv_lambda0 = 1. / lambda0;
324  for (size_t e = 0; e < getNumElements (); e++)
325  {
326  double w = getFraction (e); // Weight fraction.
327  double Z = m_element[e]->getZ (); // Atomic number
328  double A = m_element[e]->getA (); // Atomic mass.
329  double N = m_element[e]->getN (); // Number of nucleons.
330  double dovera = m_density ? m_density / A : 0; // don't crash if both are 0
331  double n = m_fraction[e] * GeoModelKernelUnits::Avogadro * dovera; // Number density.
332  int iZ = (int) (m_element[e]->getZ () + 0.5) - 1; // Atomic number index
333 
334  dEDxConstant += w * C0 * dovera * Z;
335  // Make sure we don't overflow the table:
336  // the `Ether' special `element' has Z set to 500.
337  if (iZ >= 0 && iZ < (std::end(s_ionizationPotential)-std::begin(s_ionizationPotential)))
338  dEDxI0 += w * s_ionizationPotential[iZ];
339  NILinv += n * pow (N, 2.0 / 3.0) * GeoModelKernelUnits::amu * inv_lambda0;
340 
341  double nAtomsPerVolume = A ? GeoModelKernelUnits::Avogadro*m_density*m_fraction[e]/A : 0.;
342  radInv += (nAtomsPerVolume*m_element[e]->getRadTsai());
343  }
344  m_dedDxConst = dEDxConstant;
345  m_deDxI0 = dEDxI0 ;
346  m_intLength = NILinv ? 1.0 / NILinv : 0;
347  m_radLength = radInv ? 1.0 / radInv : 0;
348 }
349 
350 double GeoMaterial::getDeDxConstant () const
351 {
352  if (!m_locked)
353  throw std::out_of_range ("Material accessed before lock");
354  return m_dedDxConst;
355 }
356 
357 double GeoMaterial::getDeDxI0 () const
358 {
359  if (!m_locked)
360  throw std::out_of_range ("Material accessed before lock");
361  return m_deDxI0;
362 }
363 
364 double GeoMaterial::getDeDxMin () const
365 {
366  //------------------------------------------------------------//
367  static const double ConstToMin = 11.528; //
368  //------------------------------------------------------------//
369 
370  if (!m_locked)
371  throw std::out_of_range ("Material accessed before lock");
372 
373  // -----------------------------------------------------------//
374  // See: Paul Avery's notes on track fitting, CBX 92-39. //
375  // Good for typical materials //
376  // -----------------------------------------------------------//
377 
378  return m_dedDxConst * ConstToMin;
379 
380 }
381 
382 double GeoMaterial::getRadLength () const
383 {
384  if (!m_locked)
385  throw std::out_of_range ("Material accessed before lock");
386  return m_radLength;
387 }
388 
389 double GeoMaterial::getIntLength () const
390 {
391  if (!m_locked)
392  throw std::out_of_range ("Material accessed before lock");
393  return m_intLength;
394 }
395 
396 unsigned int GeoMaterial::getNumElements () const
397 {
398  if (!m_locked)
399  throw std::out_of_range ("Material accessed before lock");
400  return m_element.size ();
401 }
402 
403 const GeoElement* GeoMaterial::getElement (unsigned int i) const
404 {
405  if (!m_locked)
406  throw std::out_of_range ("Material accessed before lock");
407  return m_element[i];
408 }
409 
410 double GeoMaterial::getFraction (int i) const
411 {
412  if (!m_locked)
413  throw std::out_of_range ("Material accessed before lock");
414  return m_fraction[i];
415 }
GeoElement
This class represents chemical elements. Elements are allocated on the heap and deleted automatically...
Definition: GeoElement.h:20
GeoFullPhysVol::add
void add(GeoGraphNode *graphNode)
Adds a Graph Node to the Geometry Graph.
Definition: GeoFullPhysVol.cxx:26
GeoMaterial
This class represents materials. A material is a collection of elements. In addition it has a density...
Definition: GeoMaterial.h:33
Generated on Thu May 28 2020 05:11:08 for GeoModel - GeoModelCore by   doxygen 1.8.11