ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TabulateEnergyLoss.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TabulateEnergyLoss.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019 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 
11 #include <cstddef>
12 #include <cstdlib>
13 #include <iomanip>
14 #include <iostream>
15 
20 #include "Acts/Utilities/Units.hpp"
22 
23 using namespace Acts::UnitLiterals;
24 
25 static constexpr int width = 11;
26 static constexpr int precision = 3;
27 static constexpr char separator = ' ';
28 
29 static void printHeader(std::ostream& os, const Acts::MaterialProperties& slab,
30  Acts::PdgParticle pdg, float mass, float charge) {
31  os << "# material: " << slab << '\n';
32  os << "# particle pdg id: " << pdg << '\n';
33  os << "# particle mass: " << mass / 1_MeV << "MeV\n";
34  os << "# particle charge: " << charge / 1_e << "e\n";
35  os << "# particle momentum is given in GeV\n";
36  os << "# tabulated energy loss is given in MeV\n";
37  os << "# delta is the total energy loss\n";
38  os << "# delta_ion is the energy loss due to ionisation and excitation\n";
39  os << "# delta_rad is the energy loss due to radiative effects\n";
40  os << "# sigma is the width of the enery loss distribution\n";
41  // column names
42  os << std::left;
43  os << std::setw(width) << "momentum" << separator;
44  os << std::setw(width) << "beta" << separator;
45  os << std::setw(width) << "beta_gamma" << separator;
46  os << std::setw(width) << "delta" << separator;
47  os << std::setw(width) << "delta_ion" << separator;
48  os << std::setw(width) << "delta_rad" << separator;
49  os << std::setw(width) << "sigma" << '\n';
50 }
51 
52 static void printLine(std::ostream& os, float mass, float momentum, float delta,
53  float deltaIon, float deltaRad, float sigma) {
54  const auto energy = std::sqrt(mass * mass + momentum * momentum);
55  const auto beta = momentum / energy;
56  const auto betaGamma = momentum / mass;
57  os << std::right << std::fixed << std::setprecision(precision);
58  os << std::setw(width) << momentum / 1_GeV << separator;
59  os << std::setw(width) << beta << separator;
60  os << std::setw(width) << betaGamma << separator;
61  os << std::setw(width) << delta / 1_MeV << separator;
62  os << std::setw(width) << deltaIon / 1_MeV << separator;
63  os << std::setw(width) << deltaRad / 1_MeV << separator;
64  os << std::setw(width) << sigma / 1_MeV << '\n';
65 }
66 
67 int main(int argc, char const* argv[]) {
68  // handle input arguments
69  if (argc != 6) {
70  std::cerr << "usage: " << argv[0] << " thickness pdg p_min p_max n\n";
71  std::cerr << "\n";
72  std::cerr << "tabulate energy loss over a configurable momentum range.\n";
73  std::cerr << "\n";
74  std::cerr << "parameters:\n";
75  std::cerr << " thickness: material thickness in mm\n";
76  std::cerr << " pdg: PDG particle type identifier\n";
77  std::cerr << " p_{min/max}: momentum range in GeV\n";
78  std::cerr << " n: number of points in the momentum range\n";
79  return EXIT_FAILURE;
80  }
81  const auto thickness = atof(argv[1]) * 1_mm;
82  const auto pdg = static_cast<Acts::PdgParticle>(atoi(argv[2]));
83  const auto mass = ActsFatras::findMass(pdg);
84  const auto charge = ActsFatras::findCharge(pdg);
85  const auto pmin = atof(argv[3]) * 1_GeV;
86  const auto pmax = atof(argv[4]) * 1_GeV;
87  const auto deltap = (pmax - pmin) / atoi(argv[5]);
88 
89  // use fixed material (beryllium) for now
90  // TODO make material configurable by command line
91  const Acts::Material material(35.28_cm, 42.10_cm, 9.012, 4, 1.848_g / 1_cm3);
92  const Acts::MaterialProperties slab(material, thickness);
93 
94  printHeader(std::cout, slab, pdg, mass, charge);
95  // scan momentum
96  for (auto p = pmin; p < pmax; p += deltap) {
97  const auto qOverP = charge / p;
98 
99  // TODO make mean/mode configurable by command line
100  const auto delta = computeEnergyLossMean(slab, pdg, mass, qOverP, charge);
101  const auto deltaIon =
102  Acts::computeEnergyLossBethe(slab, pdg, mass, qOverP, charge);
103  const auto deltaRad =
104  computeEnergyLossRadiative(slab, pdg, mass, qOverP, charge);
105  const auto sigma =
107 
108  printLine(std::cout, mass, p, delta, deltaIon, deltaRad, sigma);
109  }
110 
111  return EXIT_SUCCESS;
112 }