Jpp  17.1.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JPlotNPE.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 "TH2D.h"
11 
12 #include "JTools/JFunction1D_t.hh"
16 #include "JPhysics/JPDFTable.hh"
17 #include "JPhysics/JPDFTypes.hh"
18 #include "JPhysics/JNPETable.hh"
19 #include "JPhysics/JGeanz.hh"
20 #include "JGeometry3D/JAngle3D.hh"
21 #include "Jeep/JPrint.hh"
22 #include "Jeep/JParser.hh"
23 #include "Jeep/JMessage.hh"
24 
25 
26 /**
27  * \file
28  *
29  * Program to plot PDF of Cherenkov light from EM-shower using interpolation tables.
30  * \author lquinn
31  */
32 int main(int argc, char **argv)
33 {
34  using namespace std;
35  using namespace JPP;
36 
37  typedef JAbstractHistogram<double> JHistogram_t;
38 
39  vector<string> inputFile;
40  string outputFile;
41  JAngle3D dir;
42  JHistogram_t x;
43  JHistogram_t y;
44  double E;
45  int numberOfPoints;
46  int debug;
47 
48  try {
49 
50  JParser<> zap("Program to plot PDF of Cherenkov light from EM-shower using interpolation tables.");
51 
52  zap['f'] = make_field(inputFile);
53  zap['o'] = make_field(outputFile) = "";
54  zap['D'] = make_field(dir, "(theta, phi) of PMT [rad]");
55  zap['x'] = make_field(x, "histogram x-binning") = JHistogram_t();
56  zap['y'] = make_field(y, "histogram y-binning") = JHistogram_t();
57  zap['E'] = make_field(E, "Energy [GeV]") = 0.0;
58  zap['N'] = make_field(numberOfPoints) = 10;
59  zap['d'] = make_field(debug) = 0;
60 
61  zap(argc, argv);
62  }
63  catch(const exception &error) {
64  FATAL(error.what() << endl);
65  }
66 
67 
69 
70  if (E > 0.0) {
71  for (int i = 0; i != numberOfPoints; ++i) {
72  Z.push_back(geanz.getLength(E, (i + 0.5) / (double) numberOfPoints));
73  }
74  }
75 
76 
77  typedef JPolint0Function1D_t JFunction1D_t;
78  typedef
79  JMapList<JPolint2FunctionalMap,
80  JMapList<JPolint2FunctionalMap,
81  JMapList<JPolint1FunctionalGridMap,
82  JMapList<JPolint1FunctionalGridMap> > > > JMapList_t;
83  typedef JPDFTable<JFunction1D_t, JMapList_t> JPDF_t;
84  typedef JNPETable<double, double, JMapList_t> JNPE_t;
85 
86  JDistance<double>::precision = 1.0e-10;
87 
88  const int N = inputFile.size();
89 
90  JPDF_t pdf[N];
91  JNPE_t npe[N];
92 
93  try {
94 
95  for (int i = 0; i != N; ++i) {
96 
97  NOTICE("loading input from file " << inputFile[i] << "... " << flush);
98 
99  pdf[i].load(inputFile[i].c_str());
100 
101  pdf[i].setExceptionHandler(new JFunction1D_t::JDefaultResult(JMATH::zero));
102 
103  npe[i] = JNPE_t(pdf[i]);
104 
105  NOTICE("OK" << endl);
106  }
107  }
108  catch(const JException& error) {
109  FATAL(error.what() << endl);
110  }
111 
112 
113  TFile out(outputFile.c_str(), "recreate");
114 
115  if (!x.is_valid()) { x = JHistogram_t(150, 0.0, 150.0); }
116  if (!y.is_valid()) { y = JHistogram_t(200, -1.0, +1.0); }
117 
118  TH2D h0("h0", "PDF Projection; D [m]; cos #theta_{0}",
119  x.getNumberOfBins(), x.getLowerLimit(), x.getUpperLimit(),
120  y.getNumberOfBins(), y.getLowerLimit(), y.getUpperLimit());
121 
122  for (int ix = 1; ix <= h0.GetNbinsX(); ++ix) {
123  for (int iy = 1; iy <= h0.GetNbinsY(); ++iy) {
124 
125  const double cd = h0.GetYaxis()->GetBinCenter(iy);
126  const double D = h0.GetXaxis()->GetBinCenter(ix);
127 
128  double Y = 0.0;
129 
130  if (!Z.empty()) {
131 
132  const double W = 1.0 / (double) Z.size();
133 
134  for (vector<double>::const_iterator z = Z.begin(); z != Z.end(); ++z) {
135 
136  const double __D = sqrt(D*D - 2.0*(D*cd)*(*z) + (*z)*(*z));
137  const double __cd = (D * cd - (*z)) / __D;
138 
139  for (int i = 0; i != N; ++i) {
140  try {
141  Y += W * npe[i](__D, __cd, dir.getTheta(), dir.getPhi());
142  }
143  catch(const exception& error) {}
144  }
145  }
146 
147  } else {
148 
149  for (int i = 0; i != N; ++i) {
150  Y += npe[i](D, cd, dir.getTheta(), dir.getPhi());
151  }
152  }
153 
154  h0.SetBinContent(ix, iy, Y);
155  }
156  }
157 
158  out.Write();
159  out.Close();
160 }
Utility class to parse command line options.
Definition: JParser.hh:1517
int main(int argc, char *argv[])
Definition: Main.cc:15
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
string outputFile
Various implementations of functional maps.
Numbering scheme for PDF types.
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1993
static const JGeanz geanz(1.85, 0.62, 0.54)
Function object for longitudinal EM-shower profile.
#define NOTICE(A)
Definition: JMessage.hh:64
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
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
Longitudinal emission profile EM-shower.
double getLength(const double E, const double P, const double eps=1.0e-3) const
Get shower length for a given integrated probability.
Definition: JGeanz.hh:122
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
int debug
debug level