Jpp  17.0.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JPlotCDF.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 #include <map>
8 
9 #include "TROOT.h"
10 #include "TFile.h"
11 #include "TH1D.h"
12 #include "TProfile.h"
13 #include "TRandom3.h"
14 
15 #include "JTools/JFunction1D_t.hh"
18 #include "JPhysics/JCDFTable.hh"
19 #include "JPhysics/JPDFTypes.hh"
20 #include "JPhysics/JPDFToolkit.hh"
21 #include "JGeometry3D/JAngle3D.hh"
22 #include "JROOT/JRootToolkit.hh"
23 #include "JROOT/JManager.hh"
24 #include "Jeep/JTimer.hh"
25 #include "Jeep/JParser.hh"
26 #include "Jeep/JMessage.hh"
27 
28 
29 /**
30  * \file
31  *
32  * Program to verify generation of arrival times of Cherenkov photons from a muon using tabulated CDF.
33  * \author mdejong
34  */
35 int main(int argc, char **argv)
36 {
37  using namespace std;
38  using namespace JPP;
39 
40  typedef JAbstractHistogram<double> JHistogram_t;
41 
42  vector<string> inputFile;
43  string outputFile;
44  double E;
45  double R;
46  JAngle3D dir;
47  int numberOfEvents;
48  JHistogram_t histogram;
49  int debug;
50 
51  try {
52 
53  JParser<> zap("Program to verify generation of arrival times of Cherenkov photons from a muon using tabulated CDF.");
54 
55  zap['f'] = make_field(inputFile);
56  zap['o'] = make_field(outputFile) = "";
57  zap['E'] = make_field(E, "muon energy [GeV]") = 1.0;
58  zap['R'] = make_field(R, "distance of approach [m]");
59  zap['D'] = make_field(dir, "(theta, phi) of PMT [rad]");
60  zap['n'] = make_field(numberOfEvents) = 0;
61  zap['H'] = make_field(histogram, "histogram binning") = JHistogram_t();
62  zap['d'] = make_field(debug) = 2;
63 
64  zap(argc, argv);
65  }
66  catch(const exception& error) {
67  FATAL(error.what() << endl);
68  }
69 
70 
71  typedef JHermiteSplineFunction1D_t JFunction1D_t;
72  //typedef JSplineFunction1D_t JFunction1D_t;
73  //typedef JPolint1Function1D_t JFunction1D_t;
74  typedef JMAPLIST<JPolint1FunctionalMap,
75  JPolint1FunctionalGridMap,
76  JPolint1FunctionalGridMap>::maplist JMapList_t;
77  typedef JCDFTable<JFunction1D_t, JMapList_t> JCDF_t;
78 
79  const int N = inputFile.size();
80 
81  int type[N];
82  JCDF_t cdf [N];
83 
84  try {
85 
86  for (int i = 0; i != N; ++i) {
87 
88  NOTICE("loading input from file " << inputFile[i] << "... " << flush);
89 
90  type[i] = getPDFType(inputFile[i]);
91 
92  cdf [i].load(inputFile[i].c_str());
93 
94  NOTICE("OK" << endl);
95  }
96  }
97  catch(const JException& error) {
98  FATAL(error.what() << endl);
99  }
100 
101 
102  if (outputFile == "") {
103 
104  for ( ; ; ) {
105 
106  double x;
107 
108  cout << "> " << flush;
109  cin >> x;
110 
111  cout << " --> ";
112 
113  for (int i = 0; i != N; ++i) {
114 
115  try {
116 
117  double npe = cdf[i].getNPE (R, dir.getTheta(), dir.getPhi());
118  double t = cdf[i].getTime(R, dir.getTheta(), dir.getPhi(), x);
119 
120  if (is_bremsstrahlung(type[i])) {
121  npe *= E;
122  } else if (is_deltarays(type[i])) {
123  npe *= getDeltaRaysFromMuon(E);
124  }
125 
126  cout << ' ' << FIXED(6,2) << t << ' ' << FIXED(5,2) << npe;
127  }
128  catch(const exception& error) {
129  ERROR(error.what() << endl);
130  }
131  }
132 
133  cout << endl;
134  }
135 
136  return 0;
137  }
138 
139 
140  int function = 0;
141 
142  if (inputFile.size() == 1 && getPDFType(inputFile[0]) == DIRECT_LIGHT_FROM_MUON) {
143  function = 1;
144  }
145 
146  //const double z0 = -100.0; // [m]
147  //const double t0 = (-z0 + R * getTanThetaC()) / C; // time [ns]
148  const double t0 = 0.0; // time [ns]
149 
150  if (!histogram.is_valid()) {
151 
152  if (function == 1) {
153 
154  histogram = JHistogram_t(t0 - 20.0, t0 + 50.0);
155 
156  histogram.setBinWidth(0.5);
157 
158  } else {
159 
160  histogram = JHistogram_t(t0 - 20.0, t0 + 500.0);
161 
162  histogram.setBinWidth(2.0);
163  }
164  }
165 
166  TH1D h0("h0", NULL, histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit());
167 
168  h0.Sumw2();
169 
170  JManager<int, TH1D> H1(new TH1D("h1[%]", NULL, 100000, 0.0, +1.0));
171 
172  for (int i = 0; i != N; ++i) {
173 
174  for (int j = 1; j <= H1->GetNbinsX(); ++j) {
175 
176  try {
177 
178  const double x = H1->GetBinCenter(j);
179  const double t = cdf[i].getTime(R, dir.getTheta(), dir.getPhi(), x);
180 
181  H1[i]->SetBinContent(j, t);
182  }
183  catch(const exception& error) {
184  ERROR(error.what() << endl);
185  }
186  }
187  }
188 
189  if (numberOfEvents > 0) {
190 
191  JTimer timer;
192 
193  timer.reset();
194 
195  for (int counter = 0; counter != numberOfEvents; ++counter) {
196 
197  if (counter%1000== 0) {
198  STATUS(setw(10) << counter << '\r'); DEBUG(endl);
199  }
200 
201  timer.start();
202 
203  for (int i = 0; i != N; ++i) {
204 
205  try {
206 
207  double npe = cdf[i].getNPE(R, dir.getTheta(), dir.getPhi());
208 
209  if (is_bremsstrahlung(type[i])) {
210  npe *= E;
211  } else if (is_deltarays(type[i])) {
212  npe *= getDeltaRaysFromMuon(E);
213  }
214 
215  for (int j = gRandom->Poisson(npe); j != 0; --j) {
216 
217  const double x = gRandom->Rndm();
218  const double t = cdf[i].getTime(R, dir.getTheta(), dir.getPhi(), x);
219 
220  h0.Fill(t + t0);
221  }
222  }
223  catch(const exception& error) {
224  NOTICE(error.what() << endl);
225  }
226  }
227 
228  timer.stop();
229  }
230  STATUS(endl);
231 
232  const double w = 1.0 / (double) numberOfEvents;
233 
234  if (debug > JEEP::status_t) {
235  timer.print(cout, w);
236  }
237 
238  // normalise histogram (dP/dt)
239 
240  for (int j = 1; j <= h0.GetNbinsX(); ++j) {
241 
242  const double y = h0.GetBinContent(j);
243  const double z = h0.GetBinError (j);
244  const double dx = h0.GetBinWidth (j);
245 
246  h0.SetBinContent(j, y / (dx*numberOfEvents));
247  h0.SetBinError (j, z / (dx*numberOfEvents));
248  }
249  }
250 
251  TFile out(outputFile.c_str(), "recreate");
252 
253  out << h0 << H1;
254 
255  out.Write();
256  out.Close();
257 }
Utility class to parse command line options.
Definition: JParser.hh:1500
data_type w[N+1][M+1]
Definition: JPolint.hh:757
int main(int argc, char *argv[])
Definition: Main.cc:15
Auxiliary methods for PDF calculations.
then usage E
Definition: JMuonPostfit.sh:35
#define STATUS(A)
Definition: JMessage.hh:63
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
Dynamic ROOT object management.
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
status
Definition: JMessage.hh:30
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
#define NOTICE(A)
Definition: JMessage.hh:64
#define ERROR(A)
Definition: JMessage.hh:66
int debug
debug level
Definition: JSirene.cc:66
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:43
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
then for FUNCTION in pdf npe cdf
Definition: JPlotNPE-PDG.sh:73
int getPDFType(const std::string &file_name)
Get PDF type.
Definition: JPDFTypes.hh:77
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:682
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62