Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JDeltaRays.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <iomanip>
4
5#include "TROOT.h"
6#include "TFile.h"
7#include "TH1D.h"
8#include "TF1.h"
9#include "TFitResult.h"
10
11#include "JROOT/JMinimizer.hh"
12
15
17
18#include "JTools/JRange.hh"
19
20#include "Jeep/JPrint.hh"
21#include "Jeep/JParser.hh"
22#include "Jeep/JMessage.hh"
23
24
25namespace {
26
27 using namespace JPP;
28
29 static const char* const electron = "electron";
30 static const char* const positron = "positron";
31 static const char* const muon = "muon";
32 static const char* const tau = "tau";
33
34
35 /**
36 * Auxiliary class for 2nd order polynomial form factor.
37 */
38 struct JFormFactor
39 {
40 /**
41 * Constructor.
42 *
43 * \param a second order polynomial coefficient [GeV^-2]
44 * \param b first order polynomial coefficient [GeV^-1]
45 * \param c zeroth order polynomial coefficient [unit]
46 */
47 JFormFactor(const double a,
48 const double b,
49 const double c) :
50 a(a), b(b), c(c)
51 {}
52
53
54 /**
55 * Get form factor for given delta-ray kinetic energy.
56 *
57 * \param T delta-ray kinetic energy [GeV]
58 * \return form factor [unit]
59 */
60 double operator()(const double T) const
61 {
62 return c + T*(b + T*a);
63 }
64
65
66 private:
67 double a; // second order polynomial coefficient [GeV^-2]
68 double b; // first order polynomial coefficient [GeV^-1]
69 double c; // zeroth order polynomial coefficient [unit]
70 };
71
72
73 /**
74 * Get number of delta-rays due to delta-rays per unit track length\n
75 * for an ionising particle with given energy and given mass.
76 *
77 * N.B: For this function a form factor \f$F(T) = \frac{1}{2E^2}T^2 - \frac{\beta^2}{T_{\rm max}} + 1\f$ is assumed, where\n:
78 * - \f$T\f$ corresponds to the delta-ray kinetic energy,\n
79 * - \f$T_{\rm max}\f$ to the maximum delta-ray kinetic energy,\n
80 * - \f$E\f$ to the ionising particle's energy and
81 * - \f$\beta\f$ to the corresponding particle velocity, relative to the speed of light.
82 *
83 * \param E particle energy [GeV]
84 * \param M particle mass [GeV]
85 * \param Tmin minimum delta-ray kinetic energy [GeV]
86 * \param Tmax maximum delta-ray kinetic energy [GeV]
87 * \param Z atomic number [unit]
88 * \param A atomic mass [g/mol]
89 * \return delta-ray count [m^-1]
90 */
91 inline double getCount(const double E,
92 const double M,
93 const double Tmin,
94 const double Tmax,
95 const double Z,
96 const double A)
97 {
98 if (Tmax > Tmin) {
99
100 const double K = 0.307075; // [MeV mol^-1 cm^2]
101
102 const double gamma = E / M;
103 const double beta = sqrt((1.0 - 1.0/gamma) * (1.0 + 1.0/gamma));
104
105 const double a = 0.5/(E*E);
106 const double b = beta*beta/Tmax;
107 const double c = 1.0;
108
109 const double W = 0.5 * K * (Z/A) * (1.0/(beta*beta)); // [MeV g^-1 cm^2]
110
111 const double dT = Tmax - Tmin;
112 const double rT = Tmax / Tmin;
113 const double mT = Tmax * Tmin;
114
115 return W * (a*dT - b*log(rT) + c*dT/mT) * DENSITY_SEA_WATER * 1.0e2; // [m^-1]
116
117 } else {
118
119 return 0.0;
120 }
121 }
122
123
124 /**
125 * Get number of delta-rays due to delta-rays per unit track length\n
126 * for an ionising particle with given energy and given mass.\n\n
127 *
128 * The template parameter corresponds to a class which contains an `operator()(const double)`\n
129 * to compute the form factor corresponding to a given delta-ray kinetic energy.
130 *
131 * \param E particle energy [GeV]
132 * \param M particle mass [GeV]
133 * \param Tmin minimum delta-ray kinetic energy [GeV]
134 * \param Tmax maximum delta-ray kinetic energy [GeV]
135 * \param Z atomic number [unit]
136 * \param A atomic mass [g/mol]
137 * \param F form factor functor
138 * \param N number of points for numeric integration
139 * \return delta-ray count [m^-1]
140 */
141 template<class JFormFactor_t>
142 inline double getCount(const double E,
143 const double M,
144 const double Tmin,
145 const double Tmax,
146 const double Z,
147 const double A,
148 const JFormFactor_t& F,
149 const int N = 1000000)
150 {
151 using namespace JPP;
152
153 if (Tmax > Tmin) {
154
155 const double K = 0.307075; // [MeV mol^-1 cm^2]
156
157 const double gamma = E / M;
158 const double beta = sqrt((1.0 + 1.0/gamma) * (1 - 1.0/gamma));
159
160 const double W = 0.5 * K * (Z/A) * (1.0/(beta*beta)); // [MeV g^-1 cm^2]
161
162 const double xmin = 1.0/Tmax;
163 const double xmax = 1.0/Tmin;
164 const double dx = (xmax - xmin) / ((double) N);
165
166 double count = 0.0;
167
168 for (double x = xmin; x <= xmax; x += dx) {
169
170 const double T = 1.0/x;
171 const double y = W * F(T) * dx; // [MeV g^-1 cm^-2]
172
173 count += y;
174 }
175
176 return count * DENSITY_SEA_WATER * 1.0e2; // [m^-1]
177
178 } else {
179
180 return 0.0;
181 }
182 }
183}
184
185
186/**
187 * \file
188 *
189 * Program to determine the energy loss due to visible delta-rays.\n
190 * Energy loss due to delta-rays and maximal kinetic energy (Tmax) are taken from reference
191 * https://pdg.lbl.gov/2020/reviews/rpp2020-rev-passage-particles-matter.pdf
192 * \author mdejong, bjung
193 */
194int main(int argc, char **argv)
195{
196 using namespace std;
197 using namespace JPP;
198
199 typedef JRange<double> JRange_t;
200
201 string outputFile;
202 JRange_t T_GeV;
203 int numberOfPoints;
204 bool use_numerical;
205 string lepton;
206 double A;
207 double Z;
208 string option;
209 double precision;
210 int debug;
211
212 try {
213
214 JParser<> zap("Program to determine the energy loss due to visible delta-rays.");
215
216 zap['o'] = make_field(outputFile) = "delta-rays.root";
217 zap['n'] = make_field(numberOfPoints, "points for integration") = 1000000;
218 zap['N'] = make_field(use_numerical, "perform numeric integration");
219 zap['T'] = make_field(T_GeV, "kinetic energy range of electron [GeV]") = JRange_t();
220 zap['L'] = make_field(lepton) = muon, tau, positron, electron;
221 zap['A'] = make_field(A, "atomic mass") = 18.0;
222 zap['Z'] = make_field(Z, "atomic number") = 10.0;
223 zap['O'] = make_field(option) = "";
224 zap['e'] = make_field(precision) = 1.0e-6;
225 zap['d'] = make_field(debug) = 1;
226
227 zap(argc, argv);
228 }
229 catch(const exception &error) {
230 FATAL(error.what() << endl);
231 }
232
233 if (option.find('R') == string::npos) { option += 'R'; }
234 if (option.find('S') == string::npos) { option += 'S'; }
235 if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
236
237 double MASS_LEPTON;
238
239 if (lepton == muon)
240 MASS_LEPTON = MASS_MUON;
241 else if (lepton == tau)
242 MASS_LEPTON = MASS_TAU;
243 else if (lepton == positron)
244 MASS_LEPTON = MASS_ELECTRON;
245 else if (lepton == electron)
246 MASS_LEPTON = MASS_ELECTRON;
247 else
248 FATAL("Invalid lepton " << lepton << endl);
249
250 const double Tmin = (T_GeV.is_valid() ?
252
253 TFile out(outputFile.c_str(), "recreate");
254
255 const double error = 1.0e-12;
256
257 const double xmin = log10(MASS_LEPTON);
258 const double xmax = 9.25;
259
260 TH1D h0("h0", NULL, 320, xmin, xmax);
261 TH1D h1("h1", NULL, 320, xmin, xmax);
262
263 for (int i = 1; i <= h0.GetNbinsX(); ++i) {
264
265 const double x = h0.GetBinCenter(i);
266 const double E = pow(10.0,x); // particle energy [GeV]
267
268 double Tmax = (lepton != electron ?
269 getDeltaRayTmax (E, MASS_LEPTON) :
271
272 if (T_GeV.is_valid()) { Tmax = T_GeV.constrain(Tmax); }
273
274 double dEdx = 0.0;
275 double dNdx = 0.0;
276
277 if (use_numerical) {
278
279 const double gamma = E / MASS_LEPTON;
280 const double beta = sqrt((1.0 + 1.0/gamma) * (1 - 1.0/gamma));
281
282 const JFormFactor F(0.5/(E*E), -beta*beta/Tmax, 1.0);
283
284 dEdx = getDeltaRays(E, MASS_LEPTON, Tmin, Tmax, Z, A, F, numberOfPoints) / DENSITY_SEA_WATER * 1e1; // [MeV g^-1 cm^2]
285 dNdx = getCount (E, MASS_LEPTON, Tmin, Tmax, Z, A, F, numberOfPoints) / DENSITY_SEA_WATER * 1e-2; // [g^-1 cm^2]
286
287 } else {
288
289 dEdx = getDeltaRays(E, MASS_LEPTON, Tmin, Tmax, Z, A) / DENSITY_SEA_WATER * 1e1; // [MeV g^-1 cm^2]
290 dNdx = getCount (E, MASS_LEPTON, Tmin, Tmax, Z, A) / DENSITY_SEA_WATER * 1e-2; // [g^-1 cm^2]
291 }
292
293 DEBUG("E [GeV] " << SCIENTIFIC(8,2) << E << endl);
294 DEBUG("Tmin [GeV] " << SCIENTIFIC(8,2) << Tmin << endl);
295 DEBUG("Tmax [GeV] " << SCIENTIFIC(8,2) << Tmax << endl);
296 DEBUG("dE/dx [MeV g^-1 cm^2] " << FIXED(5,4) << dEdx << endl);
297 DEBUG("dE/dx [g^-1 cm^2] " << FIXED(5,2) << dNdx << endl);
298
299 h0.SetBinContent(i, dEdx);
300 h0.SetBinError (i, error);
301 h1.SetBinContent(i, dNdx);
302 }
303
304 if (use_numerical) {
305
306 TF1 f1("f1", "[0] + [1]*x + [2]*x*x + [3]*x*x*x");
307
308 if (option.find('W') == string::npos) {
309 option += "W";
310 }
311
312 f1.SetParameter(0, h0.GetMinimum());
313 f1.SetParameter(1, (h0.GetMaximum() - h0.GetMinimum()) / (h0.GetXaxis()->GetXmax() - h0.GetXaxis()->GetXmin()));
314 f1.SetParameter(2, 0.0);
315 f1.SetParameter(3, 0.0);
316
317 f1.SetLineColor(kRed);
318
319 // fit
320
321 const TFitResultPtr result = h0.Fit(&f1, option.c_str(), "same", xmin, xmax);
322
323 cout << "chi2/NDF " << result->Chi2() << '/' << result->Ndf() << endl;
324
325 cout << "\t// " << lepton << endl;
326 cout << "\t// dE/dX = a + bx + cx^2 + dx^3; // [MeV g^-1 cm^2]; x = log10(E/GeV);" << endl;
327
328 for (int i = 0; i != 4; ++i) {
329 cout << "\tstatic const double " << (char) ('a' + i) << " = " << SCIENTIFIC(10,3) << f1.GetParameter(i) << ";" << endl;
330 }
331
332 double Emin = 1 * MASS_LEPTON;
333 double Emax = 5 * MASS_LEPTON;
334
335 for (double E = 0.5 * (Emin + Emax); ; E = 0.5 * (Emin + Emax)) {
336
337 const double Tmax = getDeltaRayTmax(E, MASS_LEPTON);
338
339 if (fabs(Tmax - Tmin) < precision) {
340 cout << "\tstatic const double Emin = " << FIXED(7,5) << E << "; // [GeV]" << endl;
341 break;
342 }
343
344 if (Tmax > Tmin)
345 Emax = E;
346 else
347 Emin = E;
348 }
349 }
350
351 out.Write();
352 out.Close();
353}
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
string outputFile
int main(int argc, char **argv)
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
Auxiliary methods for PDF calculations.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Physics constants.
I/O formatting auxiliaries.
Auxiliary class to define a range between two values.
int numberOfPoints
Definition JResultPDF.cc:22
Utility class to parse command line options.
Definition JParser.hh:1698
Range of values.
Definition JRange.hh:42
bool is_valid() const
Check validity of range.
Definition JRange.hh:311
T constrain(argument_type x) const
Constrain value to range.
Definition JRange.hh:350
static const uint32_t K[64]
Definition crypt.cc:77
const double a
const double xmax
const double xmin
double getKineticEnergy(const Trk &trk)
Get track kinetic energy.
@ debug_t
debug
Definition JMessage.hh:29
size_t getCount(const array_type< T > &buffer, const JCompare_t &compare)
Count number of unique values.
T pow(const T &x, const double y)
Power .
Definition JMath.hh:97
double getDeltaRayTmin()
Get minimum delta-ray kinetic energy.
static const double DENSITY_SEA_WATER
Fixed environment values.
static const double MASS_MUON
muon mass [GeV]
static const double MASS_ELECTRON
electron mass [GeV]
double getDeltaRayTmax(const double E, const double M)
Get maximum delta-ray kinetic energy for given lepton energy and mass.
static const double MASS_TAU
tau mass [GeV]
double getDeltaRays(const double E, const double M, const double Tmin, const double Tmax, const double Z, const double A)
Get equivalent EM-shower energy due to delta-rays per unit track length for an ionising particle with...
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
return result
Definition JPolint.hh:862
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
Type definition of range.
Definition JHead.hh:43
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488