Jpp  test_elongated_shower_pde
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* 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.\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 energy of electron for Cherenkov threshold [GeV]
93 
95 
96  // Minimal kinetic energy of electron [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 
113  for (int i = 1; i <= h0.GetNbinsX(); ++i) {
114 
115  const double x = h0.GetBinCenter(i);
116  const double E = pow(10.0,x); // muon energy [GeV]
117 
118  const double gamma = E / MASS_LEPTON;
119  const double beta = sqrt(1.0 - 1.0/(gamma*gamma));
120 
121 
122  // Maximal kinetic energy of electron [GeV]
123 
124  double Tmax = (lepton == electron ?
125  (0.5*sqrt(E*E - MASS_ELECTRON*MASS_ELECTRON)) :
126  (2.0*MASS_ELECTRON*beta*beta*gamma*gamma) / (1.0 + 2.0*gamma*RATIO + RATIO*RATIO));
127 
128  if (T_GeV.is_valid()) {
129  if (Tmax > T_GeV.getUpperLimit()) {
130  Tmax = T_GeV.getUpperLimit();
131  }
132  }
133 
134 
135  double value = 0.0; // energy of delta-rays
136 
137  if (Tmax > Tmin) {
138 
139  const double a = 1.0; // Form factor formula
140  const double b = -beta*beta/Tmax; //
141  const double c = 0.5/(E*E); // F(T) = a + b*T + c*T^2
142 
143  const double xmin = log(Tmin); // T * 1/T^2 * dT = d log(T)
144  const double xmax = log(Tmax);
145  const double dx = (xmax - xmin) / numberOfPoints;
146 
147  const double W = dx * 0.5 * K * (Z/A) * (1.0/(beta*beta));
148 
149  for (double x = xmin; x <= xmax; x += dx) {
150 
151  const double T = exp(x);
152 
153  const double F = a + T*(b + T*(c)); // Form factor formula
154  const double y = W * F;
155 
156  value += y;
157  }
158  }
159 
160  DEBUG("E [GeV] " << SCIENTIFIC(8,2) << E << endl);
161  DEBUG("Tmin [GeV] " << SCIENTIFIC(8,2) << Tmin << endl);
162  DEBUG("Tmax [GeV] " << SCIENTIFIC(8,2) << Tmax << endl);
163  DEBUG("dE/dx [MeV g^-1 cm^2] " << FIXED(5,4) << value << endl);
164 
165  h0.SetBinContent(i, value);
166  }
167 
168 
169  TF1 f1("f1", "[0] + [1]*x + [2]*x*x + [3]*x*x*x");
170 
171  if (option.find('W') == string::npos) {
172  option += "W";
173  }
174 
175  f1.SetParameter(0, h0.GetMinimum());
176  f1.SetParameter(1, (h0.GetMaximum() - h0.GetMinimum()) / (h0.GetXaxis()->GetXmax() - h0.GetXaxis()->GetXmin()));
177  f1.SetParameter(2, 0.0);
178  f1.SetParameter(3, 0.0);
179 
180  f1.SetLineColor(kRed);
181 
182  // fit
183 
184  h0.Fit(&f1, option.c_str(), "same");
185 
186 
187  cout << "\t// " << lepton << endl;
188  cout << "\t// dE/dX = a + bx + cx^2 + dx^3; // [MeV g^-1 cm^2]; x = log10(E/GeV);" << endl;
189 
190  for (int i = 0; i != 4; ++i) {
191  cout << "\tstatic const double " << (char) ('a' + i) << " = " << SCIENTIFIC(10,3) << f1.GetParameter(i) << ";" << endl;
192  }
193 
194  {
195  double Tmin = sqrt(EMIN*EMIN - MASS_ELECTRON*MASS_ELECTRON);
196 
197  double Emin = 1 * MASS_LEPTON;
198  double Emax = 5 * MASS_LEPTON;
199 
200  for (double E = 0.5 * (Emin + Emax); ; E = 0.5 * (Emin + Emax)) {
201 
202  const double gamma = E / MASS_LEPTON;
203  const double beta = sqrt(1.0 - 1.0/(gamma*gamma));
204 
205  const double Tmax = (lepton == electron ?
206  (0.5*sqrt(E*E - MASS_ELECTRON*MASS_ELECTRON)) :
207  (2.0*MASS_ELECTRON*beta*beta*gamma*gamma) / (1.0 + 2.0*gamma*RATIO + RATIO*RATIO));
208 
209  if (fabs(Tmax - Tmin) < precision) {
210  cout << "\tstatic const double Emin = " << FIXED(7,5) << E << "; // [GeV]" << endl;
211  break;
212  }
213 
214  if (Tmax > Tmin)
215  Emax = E;
216  else
217  Emin = E;
218  }
219  }
220 
221  out.Write();
222  out.Close();
223 }
Utility class to parse command line options.
Definition: JParser.hh:1500
int main(int argc, char *argv[])
Definition: Main.cc:15
static const uint32_t K[64]
Definition: crypt.cc:77
then usage E
Definition: JMuonPostfit.sh:35
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
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 JAcoustics.sh -- typeset -A TRIPODS get_tripods $WORKDIR/tripod.txt TRIPODS XMEAN
Type definition of range.
Definition: JHead.hh:39
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
set_variable E_E log10(E_{fit}/E_{#mu})"
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:98
Physics constants.
int debug
debug level
Definition: JSirene.cc:68
then awk F
static const double MASS_ELECTRON
electron mass [GeV]
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
Auxiliary class to define a range between two values.
then JCalibrateToT a
Definition: JTuneHV.sh:116
Utility class to parse command line options.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
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
$WORKDIR ev_configure_domsimulator txt echo process $DOM_SIMULATOR $i $SOURCE_HOST[$index] csh c(setenv ROOTSYS $ROOTSYS &&source $JPP_DIR/setenv.csh $JPP_DIR &&($DOM_SIMULATOR\-u\$NAME\$\-H\$SERVER\$\-M\$LOGGER\$\-d $DEBUG</dev/null > &/dev/null &))'
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:484
then for LEPTON in muon tau positron electron
Definition: JDeltaRays.sh:41