Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-05 08:09:38

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2020 CERN for the benefit of the Acts project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
0008 
0009 #include "Acts/Material/detail/AverageMaterials.hpp"
0010 
0011 #include "Acts/Material/Material.hpp"
0012 
0013 #include <cmath>
0014 
0015 Acts::MaterialSlab Acts::detail::combineSlabs(const MaterialSlab& slab1,
0016                                               const MaterialSlab& slab2) {
0017   const auto& mat1 = slab1.material();
0018   const auto& mat2 = slab2.material();
0019 
0020   // NOTE 2020-08-26 msmk
0021   // the following computations provide purely geometrical averages of the
0022   // material properties. what we are really interested in are the material
0023   // properties that best describe the material interaction of the averaged
0024   // material. It is unclear if the latter can be computed in a general fashion,
0025   // so we stick to the purely geometric averages for now.
0026 
0027   // NOTE 2021-03-24 corentin
0028   // In the case of the atomic number, the averaging is based on the log
0029   // to properly account for the energy loss in multiple material.
0030 
0031   // use double for (intermediate) computations to avoid precision loss
0032 
0033   // the thickness properties are purely additive
0034   double thickness = static_cast<double>(slab1.thickness()) +
0035                      static_cast<double>(slab2.thickness());
0036 
0037   // if the two materials are the same there is no need for additional
0038   // computation
0039   if (mat1 == mat2) {
0040     return {mat1, static_cast<float>(thickness)};
0041   }
0042 
0043   double thicknessInX0 = static_cast<double>(slab1.thicknessInX0()) +
0044                          static_cast<double>(slab2.thicknessInX0());
0045   double thicknessInL0 = static_cast<double>(slab1.thicknessInL0()) +
0046                          static_cast<double>(slab2.thicknessInL0());
0047 
0048   // radiation/interaction length follows from consistency argument
0049   float x0 = thickness / thicknessInX0;
0050   float l0 = thickness / thicknessInL0;
0051 
0052   // molar amount-of-substance assuming a unit area, i.e. volume = thickness*1*1
0053   double molarAmount1 = static_cast<double>(mat1.molarDensity()) *
0054                         static_cast<double>(slab1.thickness());
0055   double molarAmount2 = static_cast<double>(mat2.molarDensity()) *
0056                         static_cast<double>(slab2.thickness());
0057   double molarAmount = molarAmount1 + molarAmount2;
0058 
0059   // handle vacuum specially
0060   if (!(0.0 < molarAmount)) {
0061     return {Material(), static_cast<float>(thickness)};
0062   }
0063 
0064   // compute average molar density by dividing the total amount-of-substance by
0065   // the total volume for the same unit area, i.e. volume = totalThickness*1*1
0066   float molarDensity = molarAmount / thickness;
0067 
0068   // assume two slabs of materials with N1,N2 atoms/molecules each with atomic
0069   // masses A1,A2 and nuclear charges. We have a total of N = N1 + N2
0070   // atoms/molecules and should have average atomic masses and nuclear charges
0071   // of
0072   //
0073   //     A = (N1*A1 + N2*A2) / (N1+N2) = (N1/N)*A1 + (N2/N)*A2 = W1*A1 + W2*A2
0074   //
0075   // the number of atoms/molecules in a given volume V with molar density rho is
0076   //
0077   //     N = V * rho * Na
0078   //
0079   // where Na is the Avogadro constant. the weighting factors follow as
0080   //
0081   //     Wi = (Vi*rhoi*Na) / (V1*rho1*Na + V2*rho2*Na)
0082   //        = (Vi*rhoi) / (V1*rho1 + V2*rho2)
0083   //
0084   // which can be computed from the molar amount-of-substance above.
0085 
0086   double molarWeight1 = molarAmount1 / molarAmount;
0087   double molarWeight2 = molarAmount2 / molarAmount;
0088   float ar = molarWeight1 * mat1.Ar() + molarWeight2 * mat2.Ar();
0089 
0090   // In the case of the atomic number, its main use is the computation
0091   // of the mean excitation energy approximated in ATL-SOFT-PUB-2008-003 as :
0092   //     I = 16 eV * Z^(0.9)
0093   // This mean excitation energy will then be used to compute energy loss
0094   // which will be proportional to :
0095   //     Eloss ~ ln(1/I)*thickness
0096   // In the case of two successive material :
0097   //     Eloss = El1 + El2
0098   //           ~ ln(Z1)*t1 + ln(Z2)*t2
0099   //           ~ ln(Z)*t
0100   // To respect this the average atomic number thus need to be defined as :
0101   //     ln(Z)*t = ln(Z1)*t1 + ln(Z2)*t2
0102   //           Z = Exp( ln(Z1)*t1/t + ln(Z2)*t2/t )
0103 
0104   double thicknessWeight1 = slab1.thickness() / thickness;
0105   double thicknessWeight2 = slab2.thickness() / thickness;
0106   float z = 0;
0107   if (mat1.Z() != 0 && mat2.Z() != 0) {
0108     z = exp(thicknessWeight1 * log(mat1.Z()) +
0109             thicknessWeight2 * log(mat2.Z()));
0110   } else {
0111     z = thicknessWeight1 * mat1.Z() + thicknessWeight2 * mat2.Z();
0112   }
0113 
0114   return {Material::fromMolarDensity(x0, l0, ar, z, molarDensity),
0115           static_cast<float>(thickness)};
0116 }