Jpp  17.3.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JPlotPDF.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <fstream>
5 #include <iomanip>
6 #include <vector>
7 
8 #include "TROOT.h"
9 #include "TFile.h"
10 #include "TH1D.h"
11 
12 #include "JTools/JFunction1D_t.hh"
16 #include "JPhysics/JPDFTable.hh"
17 #include "JPhysics/JPDFTypes.hh"
18 #include "JPhysics/JPDFToolkit.hh"
19 #include "JGeometry3D/JAngle3D.hh"
20 #include "Jeep/JPrint.hh"
21 #include "Jeep/JParser.hh"
22 #include "Jeep/JMessage.hh"
23 
24 
25 /**
26  * \file
27  *
28  * Program to plot PDF of Cherenkov light from muon using interpolation tables.
29  * \author mdejong
30  */
31 int main(int argc, char **argv)
32 {
33  using namespace std;
34  using namespace JPP;
35 
36  typedef JAbstractHistogram<double> JHistogram_t;
37 
38  vector<string> inputFile;
39  string outputFile;
40  double E;
41  double R;
42  JAngle3D dir;
43  double TTS_ns;
44  JHistogram_t histogram;
45  int debug;
46 
47  try {
48 
49  JParser<> zap("Program to plot PDF of Cherenkov light from muon using interpolation tables.");
50 
51  zap['f'] = make_field(inputFile);
52  zap['o'] = make_field(outputFile) = "";
53  zap['E'] = make_field(E, "muon energy [GeV]") = 1.0;
54  zap['R'] = make_field(R, "distance of approach [m]");
55  zap['D'] = make_field(dir, "(theta, phi) of PMT [rad]");
56  zap['T'] = make_field(TTS_ns, "PMT time smearing [ns]") = 0.0; // [ns]
57  zap['H'] = make_field(histogram, "histogram binning") = JHistogram_t();
58  zap['d'] = make_field(debug) = 0;
59 
60  zap(argc, argv);
61  }
62  catch(const exception &error) {
63  FATAL(error.what() << endl);
64  }
65 
66 
67  typedef JSplineFunction1S_t JFunction1D_t;
68  //typedef JPolint1Function1D_t JFunction1D_t;
69  typedef JMAPLIST<JPolint1FunctionalMap,
70  JPolint1FunctionalGridMap,
71  JPolint1FunctionalGridMap>::maplist JMapList_t;
72  typedef JPDFTable<JFunction1D_t, JMapList_t> JPDF_t;
73 
74  const int N = inputFile.size();
75 
76  int type[N];
77  JPDF_t pdf [N];
78 
79  try {
80 
81  for (int i = 0; i != N; ++i) {
82 
83  NOTICE("loading input from file " << inputFile[i] << "... " << flush);
84 
85  type[i] = getPDFType(inputFile[i]);
86 
87  pdf [i].load(inputFile[i].c_str());
88 
89  pdf [i].setExceptionHandler(new JFunction1D_t::JDefaultResult(JMATH::zero));
90 
91  if (TTS_ns > 0.0) {
92  pdf[i].blur(TTS_ns);
93  }
94 
95  NOTICE("OK" << endl);
96  }
97  }
98  catch(const JException& error) {
99  FATAL(error.what() << endl);
100  }
101 
102 
103  if (outputFile == "") {
104 
105  for (double dt; cin >> dt; ) {
106 
107  for (int i = 0; i != N; ++i) {
108 
109  JFunction1D_t::result_type y = pdf[i](R, dir.getTheta(), dir.getPhi(), dt);
110 
111  if (is_bremsstrahlung(type[i])) {
112  y *= E;
113  } else if (is_deltarays(type[i])) {
114  y *= getDeltaRaysFromMuon(E);
115  }
116 
117  cout << setw(2) << type[i] << ' '
118  << SCIENTIFIC(7,1) << E << ' '
119  << FIXED(5,1) << R << ' '
120  << FIXED(5,2) << dir.getTheta() << ' '
121  << FIXED(5,2) << dir.getPhi() << ' '
122  << FIXED(5,1) << dt << ' '
123  << SCIENTIFIC(9,3) << get_value(y) << endl;
124  }
125  }
126 
127  return 0;
128  }
129 
130 
131  TFile out(outputFile.c_str(), "recreate");
132 
133  int function = 0;
134 
135  if (inputFile.size() == 1 && getPDFType(inputFile[0]) == DIRECT_LIGHT_FROM_MUON) {
136  function = 1;
137  }
138 
139  //const double z0 = -100.0; // [m]
140  //const double t0 = (-z0 + R * getTanThetaC()) / C; // time [ns]
141  const double t0 = 0.0; // time [ns]
142 
143  if (!histogram.is_valid()) {
144 
145  if (function == 1) {
146 
147  histogram = JHistogram_t(t0 - 20.0, t0 + 50.0);
148 
149  histogram.setBinWidth(0.1);
150 
151  } else {
152 
153  histogram = JHistogram_t(t0 - 20.0, t0 + 500.0);
154 
155  histogram.setBinWidth(0.5);
156  }
157  }
158 
159  TH1D h0("h0", NULL, histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit());
160  TH1D h1("h1", NULL, histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit());
161  TH1D h2("h2", NULL, histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit());
162 
163  for (int i = 1; i <= h0.GetNbinsX(); ++i) {
164 
165  const double dt = h0.GetBinCenter(i) - t0;
166 
167  JFunction1D_t::result_type Y = JMATH::zero;
168 
169  for (int j = 0; j != N; ++j) {
170 
171  JFunction1D_t::result_type y = pdf[j](R, dir.getTheta(), dir.getPhi(), dt);
172 
173  if (is_bremsstrahlung(type[j])) {
174  y *= E;
175  } else if (is_deltarays(type[j])) {
176  y *= getDeltaRaysFromMuon(E);
177  }
178 
179  Y += y;
180  }
181 
182  h0.SetBinContent(i, get_value (Y));
183  h1.SetBinContent(i, get_derivative(Y));
184  h2.SetBinContent(i, get_integral (Y));
185  }
186 
187  out.Write();
188  out.Close();
189 }
Utility class to parse command line options.
Definition: JParser.hh:1517
int main(int argc, char *argv[])
Definition: Main.cc:15
Auxiliary methods for PDF calculations.
then usage E
Definition: JMuonPostfit.sh:35
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
string outputFile
direct light from muon
Definition: JPDFTypes.hh:26
Various implementations of functional maps.
Numbering scheme for PDF types.
bool is_bremsstrahlung(const int pdf)
Test if given PDF type corresponds to Cherenkov light from Bremsstrahlung.
Definition: JPDFTypes.hh:137
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1993
#define NOTICE(A)
Definition: JMessage.hh:64
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
int getPDFType(const std::string &file_name)
Get PDF type.
Definition: JPDFTypes.hh:77
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
JResultEvaluator< JResult_t >::result_type get_derivative(const JResult_t &value)
Helper method to convert function value to derivative value.
Definition: JResult.hh:1011
Utility class to parse command line options.
bool is_deltarays(const int pdf)
Test if given PDF type corresponds to Cherenkov light from delta-rays.
Definition: JPDFTypes.hh:151
double getDeltaRaysFromMuon(const double E)
Equivalent EM-shower energy due to delta-rays per unit muon track length.
Definition: JPDFToolkit.hh:67
int j
Definition: JPolint.hh:703
JResultEvaluator< JResult_t >::result_type get_integral(const JResult_t &value)
Helper method to convert function value to partial integral value.
Definition: JResult.hh:1024
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:484
JResultEvaluator< JResult_t >::result_type get_value(const JResult_t &value)
Helper method to recursively evaluate a to function value.
Definition: JResult.hh:998
int debug
debug level