EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AverageMaterials.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file AverageMaterials.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2020 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
10 
12  const MaterialSlab& slab2) {
13  const auto& mat1 = slab1.material();
14  const auto& mat2 = slab2.material();
15 
16  // NOTE 2020-08-26 msmk
17  // the following computations provide purely geometrical averages of the
18  // material properties. what we are really interested in are the material
19  // properties that best describe the material interaction of the averaged
20  // material. It is unclear if the latter can be computed in a general fashion,
21  // so we stick to the purely geometric averages for now.
22 
23  // use double for (intermediate) computations to avoid precision loss
24 
25  // the thickness properties are purely additive
26  double thickness = static_cast<double>(slab1.thickness()) +
27  static_cast<double>(slab2.thickness());
28  double thicknessInX0 = static_cast<double>(slab1.thicknessInX0()) +
29  static_cast<double>(slab2.thicknessInX0());
30  double thicknessInL0 = static_cast<double>(slab1.thicknessInL0()) +
31  static_cast<double>(slab2.thicknessInL0());
32 
33  // radiation/interaction length follows from consistency argument
34  float x0 = thickness / thicknessInX0;
35  float l0 = thickness / thicknessInL0;
36 
37  // molar amount-of-substance assuming a unit area, i.e. volume = thickness*1*1
38  double molarAmount1 = static_cast<double>(mat1.molarDensity()) *
39  static_cast<double>(slab1.thickness());
40  double molarAmount2 = static_cast<double>(mat2.molarDensity()) *
41  static_cast<double>(slab2.thickness());
42  double molarAmount = molarAmount1 + molarAmount2;
43 
44  // handle vacuum specially
45  if (not(0.0 < molarAmount)) {
46  return {Material(), static_cast<float>(thickness)};
47  }
48 
49  // compute average molar density by divding the total amount-of-substance by
50  // the total volume for the same unit area, i.e. volume = totalThickness*1*1
51  float molarDensity = molarAmount / thickness;
52 
53  // assume two slabs of materials with N1,N2 atoms/molecules each with atomic
54  // masses A1,A2 and nuclear charges. We have a total of N = N1 + N2
55  // atoms/molecules and should have average atomic masses and nuclear charges
56  // of
57  //
58  // A = (N1*A1 + N2*A2) / (N1+N2) = (N1/N)*A1 + (N2/N)*A2 = W1*A1 + W2*A2
59  // Z = (N1*Z1 + N2*Z2) / (N1+N2) = (N1/N)*Z1 + (N2/N)*Z2 = W1*Z1 * W2*Z2
60  //
61  // the number of atoms/molecues in a given volume V with molar density rho is
62  //
63  // N = V * rho * Na
64  //
65  // where Na is the Avogadro constant. the weighting factors follow as
66  //
67  // Wi = (Vi*rhoi*Na) / (V1*rho1*Na + V2*rho2*Na)
68  // = (Vi*rhoi) / (V1*rho1 + V2*rho2)
69  //
70  // which can be computed from the molar amount-of-substance above.
71  double weight1 = molarAmount1 / molarAmount;
72  double weight2 = molarAmount2 / molarAmount;
73  float ar = weight1 * mat1.Ar() + weight2 * mat2.Ar();
74  float z = weight1 * mat1.Z() + weight2 * mat2.Z();
75 
76  return {Material::fromMolarDensity(x0, l0, ar, z, molarDensity),
77  static_cast<float>(thickness)};
78 }