Jpp  18.4.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 "JPhysics/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* const electron = "electron";
27  static const char* const positron = "positron";
28  static const char* const muon = "muon";
29  static const char* const tau = "tau";
30 }
31 
32 
33 /**
34  * \file
35  *
36  * Program to determine the energy loss due to visible delta-rays.\n
37  * Energy loss due to delta-rays and maximal kinetic energy (Tmax) are taken from reference
38  * https://pdg.lbl.gov/2020/reviews/rpp2020-rev-passage-particles-matter.pdf
39  * \author mdejong
40  */
41 int main(int argc, char **argv)
42 {
43  using namespace std;
44 
46 
47  string outputFile;
48  JRange_t T_GeV;
49  int numberOfPoints;
50  string lepton;
51  string option;
52  double precision;
53  int debug;
54 
55  try {
56 
57  JParser<> zap("Program to determine the energy loss due to visible delta-rays.");
58 
59  zap['o'] = make_field(outputFile) = "delta-rays.root";
60  zap['T'] = make_field(T_GeV, "kinetic energy range of electron [GeV]") = JRange_t();
61  zap['n'] = make_field(numberOfPoints, "points for integration") = 1000000;
62  zap['L'] = make_field(lepton) = muon, tau, positron, electron;
63  zap['O'] = make_field(option) = "";
64  zap['e'] = make_field(precision) = 1.0e-6;
65  zap['d'] = make_field(debug) = 1;
66 
67  zap(argc, argv);
68  }
69  catch(const exception &error) {
70  FATAL(error.what() << endl);
71  }
72 
73  using namespace JPP;
74 
75  double MASS_LEPTON;
76 
77  if (lepton == muon)
78  MASS_LEPTON = MASS_MUON;
79  else if (lepton == tau)
80  MASS_LEPTON = MASS_TAU;
81  else if (lepton == positron)
82  MASS_LEPTON = MASS_ELECTRON;
83  else if (lepton == electron)
84  MASS_LEPTON = MASS_ELECTRON;
85  else
86  FATAL("Invalid lepton " << lepton << endl);
87 
88 
89  const double RATIO = MASS_ELECTRON / MASS_LEPTON;
90 
91 
92  // Minimal total energy of electron for Cherenkov threshold [GeV]
93 
95 
96  // Minimal kinetic energy of electron for Cherenkov threshold [GeV]
97 
98  double Tmin = sqrt(EMIN*EMIN - MASS_ELECTRON*MASS_ELECTRON);
99 
100  NOTICE("Tmin [GeV] " << FIXED(8,6) << Tmin << ' ' << FIXED(8,6) << T_GeV.getLowerLimit() << endl);
101 
102  if (T_GeV.is_valid()) {
103  if (Tmin < T_GeV.getLowerLimit()) {
104  Tmin = T_GeV.getLowerLimit();
105  }
106  }
107 
108 
109  TFile out(outputFile.c_str(), "recreate");
110 
111  TH1D h0("h0", NULL, 80, log10(MASS_LEPTON), +9.25);
112  TH1D h1("h1", NULL, 80, log10(MASS_LEPTON), +9.25);
113 
114  for (int i = 1; i <= h0.GetNbinsX(); ++i) {
115 
116  const double x = h0.GetBinCenter(i);
117  const double E = pow(10.0,x); // particle energy [GeV]
118 
119  const double gamma = E / MASS_LEPTON;
120  const double beta = sqrt(1.0 - 1.0/(gamma*gamma));
121 
122 
123  // Maximal kinetic energy of electron [GeV]
124 
125  double Tmax = (lepton == electron ?
126  (0.5*sqrt(E*E - MASS_ELECTRON*MASS_ELECTRON)) :
127  (2.0*MASS_ELECTRON*beta*beta*gamma*gamma) / (1.0 + 2.0*gamma*RATIO + RATIO*RATIO));
128 
129  if (T_GeV.is_valid()) {
130  if (Tmax > T_GeV.getUpperLimit()) {
131  Tmax = T_GeV.getUpperLimit();
132  }
133  }
134 
135 
136  double weight = 0.0;
137  double count = 0.0;
138 
139  if (Tmax > Tmin) {
140 
141  const double a = 1.0; // Form factor formula
142  const double b = -beta*beta/Tmax; //
143  const double c = 0.5/(E*E); // F(T) = a + b*T + c*T^2
144 
145  const double W = 0.5 * K * (Z/A) * (1.0/(beta*beta));
146 
147  {
148  const double xmin = log(Tmin); // T * 1/T^2 * dT = d log(T)
149  const double xmax = log(Tmax);
150  const double dx = (xmax - xmin) / numberOfPoints;
151 
152  for (double x = xmin; x <= xmax; x += dx) {
153 
154  const double T = exp(x);
155 
156  const double F = a + T*(b + T*(c)); // Form factor formula
157  const double y = W * F * dx;
158 
159  weight += y;
160  }
161  }
162  {
163  const double xmin = 1.0/Tmax; // 1/T^2 * dT = d -1/T
164  const double xmax = 1.0/Tmin;
165  const double dx = (xmax - xmin) / numberOfPoints;
166 
167  for (double x = xmin; x <= xmax; x += dx) {
168 
169  const double T = 1.0/x;
170 
171  const double F = a + T*(b + T*(c)); // Form factor formula
172  const double y = W * F * dx;
173 
174  count += y;
175  }
176  }
177  }
178 
179  DEBUG("E [GeV] " << SCIENTIFIC(8,2) << E << endl);
180  DEBUG("Tmin [GeV] " << SCIENTIFIC(8,2) << Tmin << endl);
181  DEBUG("Tmax [GeV] " << SCIENTIFIC(8,2) << Tmax << endl);
182  DEBUG("dE/dx [MeV g^-1 cm^2] " << FIXED(5,4) << weight << endl);
183  DEBUG("dE/dx [g^-1 cm^2] " << FIXED(5,2) << count << endl);
184 
185  h0.SetBinContent(i, weight);
186  h1.SetBinContent(i, count);
187  }
188 
189 
190  TF1 f1("f1", "[0] + [1]*x + [2]*x*x + [3]*x*x*x");
191 
192  if (option.find('W') == string::npos) {
193  option += "W";
194  }
195 
196  f1.SetParameter(0, h0.GetMinimum());
197  f1.SetParameter(1, (h0.GetMaximum() - h0.GetMinimum()) / (h0.GetXaxis()->GetXmax() - h0.GetXaxis()->GetXmin()));
198  f1.SetParameter(2, 0.0);
199  f1.SetParameter(3, 0.0);
200 
201  f1.SetLineColor(kRed);
202 
203  // fit
204 
205  h0.Fit(&f1, option.c_str(), "same");
206 
207 
208  cout << "\t// " << lepton << endl;
209  cout << "\t// dE/dX = a + bx + cx^2 + dx^3; // [MeV g^-1 cm^2]; x = log10(E/GeV);" << endl;
210 
211  for (int i = 0; i != 4; ++i) {
212  cout << "\tstatic const double " << (char) ('a' + i) << " = " << SCIENTIFIC(10,3) << f1.GetParameter(i) << ";" << endl;
213  }
214 
215  {
216  double Tmin = sqrt(EMIN*EMIN - MASS_ELECTRON*MASS_ELECTRON);
217 
218  double Emin = 1 * MASS_LEPTON;
219  double Emax = 5 * MASS_LEPTON;
220 
221  for (double E = 0.5 * (Emin + Emax); ; E = 0.5 * (Emin + Emax)) {
222 
223  const double gamma = E / MASS_LEPTON;
224  const double beta = sqrt(1.0 - 1.0/(gamma*gamma));
225 
226  const double Tmax = (lepton == electron ?
227  (0.5*sqrt(E*E - MASS_ELECTRON*MASS_ELECTRON)) :
228  (2.0*MASS_ELECTRON*beta*beta*gamma*gamma) / (1.0 + 2.0*gamma*RATIO + RATIO*RATIO));
229 
230  if (fabs(Tmax - Tmin) < precision) {
231  cout << "\tstatic const double Emin = " << FIXED(7,5) << E << "; // [GeV]" << endl;
232  break;
233  }
234 
235  if (Tmax > Tmin)
236  Emax = E;
237  else
238  Emin = E;
239  }
240  }
241 
242  out.Write();
243  out.Close();
244 }
const double xmax
Definition: JQuadrature.cc:24
Utility class to parse command line options.
Definition: JParser.hh:1514
int main(int argc, char *argv[])
Definition: Main.cc:15
static const uint32_t K[64]
Definition: crypt.cc:77
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
then set_variable DIR else fatal Wrong number of arguments fi for INPUT_FILE in ls rt $DIR stage * log
static const double MASS_MUON
muon mass [GeV]
static const double INDEX_OF_REFRACTION_WATER
Average index of refraction of water corresponding to the group velocity.
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
string outputFile
const JPolynome f1(1.0, 2.0, 3.0)
Function.
Type definition of range.
Definition: JHead.hh:41
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
set_variable E_E log10(E_{fit}/E_{#mu})"
then JCalibrateToT a
Definition: JTuneHV.sh:113
static const double MASS_TAU
tau mass [GeV]
do set_variable OUTPUT_DIRECTORY $WORKDIR T
#define NOTICE(A)
Definition: JMessage.hh:64
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:97
Physics constants.
then awk F
static const double MASS_ELECTRON
electron mass [GeV]
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
$WORKDIR ev_configure_dqsimulator txt echo process $DQ_SIMULATOR $i $SOURCE_HOST[$index] csh c(setenv ROOTSYS $ROOTSYS &&source $JPP_DIR/setenv.csh $JPP_DIR &&($DQ_SIMULATOR\-u\$NAME\$\-H\$SERVER\$\-M\$LOGGER\$\-d $DEBUG</dev/null > &/dev/null &))'
const double xmin
Definition: JQuadrature.cc:23
Auxiliary class to define a range between two values.
Utility class to parse command line options.
int numberOfPoints
Definition: JResultPDF.cc:22
do set_variable MODULE getModule a $WORKDIR detector_a datx L $STRING JEditDetector a $WORKDIR detector_a datx M $MODULE setz o $WORKDIR detector_a datx JEditDetector a $WORKDIR detector_b datx M $MODULE setz o $WORKDIR detector_b datx done echo Output stored at $WORKDIR detector_a datx and $WORKDIR tripod_a txt JDrawDetector2D a $WORKDIR detector_a datx a $WORKDIR detector_b datx L BL o detector $FORMAT $BATCH JDrawDetector2D T $WORKDIR tripod_a txt T $WORKDIR tripod_b txt L BL o tripod $FORMAT $BATCH JCompareDetector a $WORKDIR detector_a datx b $WORKDIR detector_b datx o $WORKDIR abc root &dev null for KEY in X Y Z
source $JPP_DIR setenv csh $JPP_DIR &dev null eval JShellParser o a A
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:486
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable STRING $argv[2] set_array QUANTILES set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable MODULE `getModule -a $DETECTOR -L "$STRING 0"` source JAcousticsToolkit.sh typeset -A TRIPODS get_tripods $WORKDIR/tripod.txt TRIPODS XMEAN
then for LEPTON in muon tau positron electron
Definition: JDeltaRays.sh:41
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62