Jpp
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 
10 #include "JTools/JConstants.hh"
11 #include "JTools/JRange.hh"
12 
13 #include "Jeep/JPrint.hh"
14 #include "Jeep/JParser.hh"
15 #include "Jeep/JMessage.hh"
16 
17 
18 namespace {
19 
20  using namespace JPP;
21 
22  static const double K = 0.307075; // [MeV g^-1 cm^2]
23  static const double Z = 10.0; // atomic number [unit]
24  static const double A = 18.0; // atomic mass [g/mol]
25 
26  static const char* electron = "electron";
27  static const char* positron = "positron";
28  static const char* muon = "muon";
29  static const char* tau = "tau";
30 }
31 
32 
33 /**
34  * \file
35  *
36  * Program to determine the energy loss due to visible delta-rays.
37  * \author mdejong
38  */
39 int main(int argc, char **argv)
40 {
41  using namespace std;
42 
43  typedef JTOOLS::JRange<double> JRange_t;
44 
45  string outputFile;
46  JRange_t T_GeV;
47  int numberOfPoints;
48  string lepton;
49  string option;
50  double precision;
51  int debug;
52 
53  try {
54 
55  JParser<> zap("Program to determine the energy loss due to visible delta-rays.");
56 
57  zap['o'] = make_field(outputFile) = "delta-rays.root";
58  zap['T'] = make_field(T_GeV, "kinetic energy range of electron [GeV]") = JRange_t();
59  zap['n'] = make_field(numberOfPoints, "points for integration") = 1000000;
60  zap['L'] = make_field(lepton) = muon, tau, positron, electron;
61  zap['O'] = make_field(option) = "";
62  zap['e'] = make_field(precision) = 1.0e-6;
63  zap['d'] = make_field(debug) = 1;
64 
65  zap(argc, argv);
66  }
67  catch(const exception &error) {
68  FATAL(error.what() << endl);
69  }
70 
71  using namespace JPP;
72 
73  double MASS_LEPTON;
74 
75  if (lepton == muon)
76  MASS_LEPTON = MASS_MUON;
77  else if (lepton == tau)
78  MASS_LEPTON = MASS_TAU;
79  else if (lepton == positron)
80  MASS_LEPTON = MASS_ELECTRON;
81  else if (lepton == electron)
82  MASS_LEPTON = MASS_ELECTRON;
83  else
84  FATAL("Invalid lepton " << lepton << endl);
85 
86 
87  const double RATIO = MASS_ELECTRON / MASS_LEPTON;
88 
89 
90  // Minimal energy of electron for Cherenkov threshold [GeV]
91 
92  const double EMIN = MASS_ELECTRON * INDEX_OF_REFRACTION_WATER / sqrt(INDEX_OF_REFRACTION_WATER*INDEX_OF_REFRACTION_WATER - 1.0);
93 
94  // Minimal kinetic energy of electron [GeV]
95 
96  double Tmin = sqrt(EMIN*EMIN - MASS_ELECTRON*MASS_ELECTRON);
97 
98  NOTICE("Tmin [GeV] " << FIXED(8,6) << Tmin << ' ' << FIXED(8,6) << T_GeV.getLowerLimit() << endl);
99 
100  if (T_GeV.is_valid()) {
101  if (Tmin < T_GeV.getLowerLimit()) {
102  Tmin = T_GeV.getLowerLimit();
103  }
104  }
105 
106 
107  TFile out(outputFile.c_str(), "recreate");
108 
109  TH1D h0("h0", NULL, 80, log10(MASS_LEPTON), +9.25);
110 
111  for (int i = 1; i <= h0.GetNbinsX(); ++i) {
112 
113  const double x = h0.GetBinCenter(i);
114  const double E = pow(10.0,x); // muon energy [GeV]
115 
116  const double gamma = E / MASS_LEPTON;
117  const double beta = sqrt(1.0 - 1.0/(gamma*gamma));
118 
119 
120  // Maximal kinetic energy of electron [GeV]
121 
122  double Tmax = (lepton == electron ?
123  (0.5*sqrt(E*E - MASS_ELECTRON*MASS_ELECTRON)) :
124  (2.0*MASS_ELECTRON*beta*beta*gamma*gamma) / (1.0 + 2.0*gamma*RATIO + RATIO*RATIO));
125 
126  if (T_GeV.is_valid()) {
127  if (Tmax > T_GeV.getUpperLimit()) {
128  Tmax = T_GeV.getUpperLimit();
129  }
130  }
131 
132 
133  double value = 0.0; // energy of delta-rays
134 
135  if (Tmax > Tmin) {
136 
137  const double a = 1.0; // Form factor formula
138  const double b = -beta*beta/Tmax; //
139  const double c = 0.5/(E*E); // F(T) = a + b*T + c*T^2
140 
141  const double xmin = log(Tmin); // T * 1/T^2 * dT = d log(T)
142  const double xmax = log(Tmax);
143  const double dx = (xmax - xmin) / numberOfPoints;
144 
145  const double W = dx * 0.5 * K * (Z/A) * (1.0/(beta*beta));
146 
147  for (double x = xmin; x <= xmax; x += dx) {
148 
149  const double T = exp(x);
150 
151  const double F = a + T*(b + T*(c)); // Form factor formula
152  const double y = W * F;
153 
154  value += y;
155  }
156  }
157 
158  DEBUG("E [GeV] " << SCIENTIFIC(8,2) << E << endl);
159  DEBUG("Tmin [GeV] " << SCIENTIFIC(8,2) << Tmin << endl);
160  DEBUG("Tmax [GeV] " << SCIENTIFIC(8,2) << Tmax << endl);
161  DEBUG("dE/dx [MeV g^-1 cm^2] " << FIXED(5,4) << value << endl);
162 
163  h0.SetBinContent(i, value);
164  }
165 
166 
167  TF1 f1("f1", "[0] + [1]*x + [2]*x*x + [3]*x*x*x");
168 
169  if (option.find('W') == string::npos) {
170  option += "W";
171  }
172 
173  f1.SetParameter(0, h0.GetMinimum());
174  f1.SetParameter(1, (h0.GetMaximum() - h0.GetMinimum()) / (h0.GetXaxis()->GetXmax() - h0.GetXaxis()->GetXmin()));
175  f1.SetParameter(2, 0.0);
176  f1.SetParameter(3, 0.0);
177 
178  f1.SetLineColor(kRed);
179 
180  // fit
181 
182  h0.Fit(&f1, option.c_str(), "same");
183 
184 
185  cout << "\t// " << lepton << endl;
186  cout << "\t// dE/dX = a + bx + cx^2 + dx^3; // [MeV g^-1 cm^2]; x = log10(E/GeV);" << endl;
187 
188  for (int i = 0; i != 4; ++i) {
189  cout << "\tstatic const double " << (char) ('a' + i) << " = " << SCIENTIFIC(10,3) << f1.GetParameter(i) << ";" << endl;
190  }
191 
192  {
193  double Tmin = sqrt(EMIN*EMIN - MASS_ELECTRON*MASS_ELECTRON);
194 
195  double Emin = 1 * MASS_LEPTON;
196  double Emax = 5 * MASS_LEPTON;
197 
198  for (double E = 0.5 * (Emin + Emax); ; E = 0.5 * (Emin + Emax)) {
199 
200  const double gamma = E / MASS_LEPTON;
201  const double beta = sqrt(1.0 - 1.0/(gamma*gamma));
202 
203  const double Tmax = (lepton == electron ?
204  (0.5*sqrt(E*E - MASS_ELECTRON*MASS_ELECTRON)) :
205  (2.0*MASS_ELECTRON*beta*beta*gamma*gamma) / (1.0 + 2.0*gamma*RATIO + RATIO*RATIO));
206 
207  if (fabs(Tmax - Tmin) < precision) {
208  cout << "\tstatic const double Emin = " << FIXED(7,5) << E << "; // [GeV]" << endl;
209  break;
210  }
211 
212  if (Tmax > Tmin)
213  Emax = E;
214  else
215  Emin = E;
216  }
217  }
218 
219  out.Write();
220  out.Close();
221 }
JTOOLS::MASS_TAU
static const double MASS_TAU
tau mass [GeV]
Definition: JConstants.hh:60
main
int main(int argc, char **argv)
Definition: JDeltaRays.cc:39
FIXED
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:481
JTOOLS::MASS_MUON
static const double MASS_MUON
muon mass [GeV]
Definition: JConstants.hh:59
JMessage.hh
JPrint.hh
numberOfPoints
int numberOfPoints
Definition: JResultPDF.cc:22
JTOOLS::JRange< double >
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1493
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
JRange.hh
debug
int debug
debug level
Definition: JSirene.cc:59
JConstants.hh
K
static const uint32_t K[64]
Definition: crypt.cc:77
JTOOLS::INDEX_OF_REFRACTION_WATER
static const double INDEX_OF_REFRACTION_WATER
average index of refraction of water
Definition: JConstants.hh:37
SCIENTIFIC
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:518
JParser.hh
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
JTOOLS::MASS_ELECTRON
static const double MASS_ELECTRON
electron mass [GeV]
Definition: JConstants.hh:58
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std
Definition: jaanetDictionary.h:36
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
Generated by   doxygen 1.8.16