Jpp  18.5.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JDemoPDF.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <fstream>
5 #include <iomanip>
6 #include <assert.h>
7 
8 #include "TROOT.h"
9 #include "TFile.h"
10 #include "TKey.h"
11 #include "TH1D.h"
12 #include "TProfile.h"
13 
14 #include "JTools/JFunction1D_t.hh"
16 #include "JTools/JQuantiles.hh"
17 #include "JPhysics/JPDFTable.hh"
18 #include "JPhysics/JPDFTypes.hh"
19 #include "JGeometry3D/JAngle3D.hh"
20 #include "Jeep/JTimer.hh"
21 #include "Jeep/JParser.hh"
22 #include "Jeep/JMessage.hh"
23 
24 
25 /**
26  * \file
27  * Demonstration program to plot RMS of arrival time of first hit as a function of
28  * the minimal distance of approach of a muon to the PMT.
29  * \author mdejong
30  */
31 int main(int argc, char **argv)
32 {
33  using namespace std;
34  using namespace JPP;
35 
36  string inputFile;
37  string outputFile;
38  double E_GeV;
39  double R_Hz;
40  JAngle3D dir;
41  int debug;
42 
43  try {
44 
45  JParser<> zap("Demonstration program to plot RMS of arrival time of first hit as a function of the minimal distance of approach of a muon to the PMT.");
46 
47  zap['f'] = make_field(inputFile);
48  zap['o'] = make_field(outputFile) = "demo.root";
49  zap['E'] = make_field(E_GeV, "muon energy [GeV]");
50  zap['R'] = make_field(R_Hz, "background rate [Hz]");
51  zap['D'] = make_field(dir, "(theta, phi) of PMT [rad]");
52  zap['d'] = make_field(debug) = 0;
53 
54  zap(argc, argv);
55  }
56  catch(const exception &error) {
57  FATAL(error.what() << endl);
58  }
59 
60 
61  typedef JSplineFunction1S_t JFunction1D_t;
62  typedef JMAPLIST<JPolint1FunctionalMap,
63  JPolint1FunctionalGridMap,
64  JPolint1FunctionalGridMap>::maplist JMapList_t;
65  typedef JPDFTable<JFunction1D_t, JMapList_t> JPDF_t;
66 
67  JFunction1D_t::JSupervisor supervisor(new JFunction1D_t::JDefaultResult(zero));
68 
69  const JPDFType_t pdf_t[] = { DIRECT_LIGHT_FROM_MUON,
73 
74  const int N = sizeof(pdf_t) / sizeof(pdf_t[0]);
75 
76  JPDF_t pdf[N];
77 
78  for (int i = 0; i != N; ++i) {
79 
80  try {
81 
82  const string file_name = getFilename(inputFile, pdf_t[i]);
83 
84  NOTICE("loading input from file " << file_name << "... " << flush);
85 
86  pdf[i].load(file_name.c_str());
87 
88  NOTICE("OK" << endl);
89  }
90  catch(const JException& error) {
91  FATAL(error.what() << endl);
92  }
93 
94  pdf[i].setExceptionHandler(supervisor);
95  }
96 
97 
98  vector<double> X; // transverse distance [m]
99 
100  if (X.empty()) {
101  X.push_back( 5.0);
102  X.push_back( 15.0);
103  X.push_back( 25.0);
104  X.push_back( 35.0);
105  X.push_back( 45.0);
106  X.push_back( 55.0);
107  X.push_back( 65.0);
108  X.push_back( 75.0);
109  X.push_back( 85.0);
110  X.push_back( 95.0);
111  }
112 
113 
114  TFile out(outputFile.c_str(), "recreate");
115 
116  vector<TH1D*> buffer;
117 
118  for (vector<double>::const_iterator i = X.begin(); i != X.end(); ++i) {
119 
120  ostringstream os;
121 
122  os << "t[" << *i << "]";
123 
124  buffer.push_back(new TH1D(os.str().c_str(), NULL, 5000, -250.0, +250.0));
125  }
126 
128 
129  {
132 
133  x.push_back(*i - 0.5*(*j-*i));
134 
135  for (; j != X.end(); ++i, ++j)
136  x.push_back(0.5*(*i+*j));
137 
138  --i;
139  --j;
140 
141  x.push_back(*j + 0.5*(*j-*i));
142  }
143 
144  TH1D rms("rms", NULL, X.size(), &x[0]);
145 
146 
147  JFunction1D_t::result_type p[4];
148 
149  JTimer timer;
150 
151  timer.reset();
152 
153  const double Tmin = -250.0; // [ns]
154  const double Tmax = +250.0; // [ns]
155 
156  int n = 0;
157 
158  for (size_t i = 0; i != X.size(); ++i) {
159 
160  TH1D* h = buffer[i];
161  const double R = X[i];
162 
163  JSplineFunction1S_t zsp;
164 
165  for (int j = 1; j <= h->GetNbinsX(); ++j) {
166 
167  const double t1 = h->GetBinCenter(j);
168 
169  timer.start();
170 
171  ++n;
172 
173  for (int k = 0; k != 4; ++k) {
174  p[k] = pdf[k](R, dir.getTheta(), dir.getPhi(), t1);
175  }
176 
177  // integral [Tmin,t1]
178  const double v =
179  p[0].v +
180  p[1].v +
181  p[2].v * E_GeV +
182  p[3].v * E_GeV +
183  R_Hz * 1e-9 * (t1 - Tmin);
184 
185  // integral [Tmin,Tmax]
186  const double V =
187  p[0].V +
188  p[1].V +
189  p[2].V * E_GeV +
190  p[3].V * E_GeV +
191  R_Hz * 1e-9 * (Tmax - Tmin);
192 
193  // function value
194  const double y =
195  p[0].f +
196  p[1].f +
197  p[2].f * E_GeV +
198  p[3].f * E_GeV +
199  R_Hz * 1e-9;
200 
201  const double W = exp(-v) * y / (1.0 - exp(-V));
202 
203  timer.stop();
204 
205  zsp[t1] = W;
206 
207  h->SetBinContent(j,W);
208  }
209 
210  zsp.compile();
211 
212  JQuantiles result(zsp);
213 
214  const double sig = result.getFWHM() * 0.5 / sqrt(2.0*log(2.0));
215 
216  rms.Fill(R, sig);
217 
218  NOTICE("integral " << R << ' ' << result.getIntegral() << endl);
219  }
220 
221  const float w = 1.0 / (float) n;
222 
223  if (debug > notice_t) {
224  timer.print(cout, w);
225  }
226 
227  out.Write();
228  out.Close();
229 }
Utility class to parse command line options.
Definition: JParser.hh:1514
data_type w[N+1][M+1]
Definition: JPolint.hh:867
then fatal No hydrophone data file $HYDROPHONE_TXT fi sort gr k
int main(int argc, char *argv[])
Definition: Main.cc:15
then set_variable DIR else fatal Wrong number of arguments fi for INPUT_FILE in ls rt $DIR stage * log
notice
Definition: JMessage.hh:32
direct light from EM showers
Definition: JPDFTypes.hh:29
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
string outputFile
direct light from muon
Definition: JPDFTypes.hh:26
Various implementations of functional maps.
Numbering scheme for PDF types.
const int n
Definition: JPolint.hh:786
scattered light from muon
Definition: JPDFTypes.hh:27
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
#define NOTICE(A)
Definition: JMessage.hh:64
scattered light from EM showers
Definition: JPDFTypes.hh:30
JPDFType_t
PDF types.
Definition: JPDFTypes.hh:24
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:40
Utility class to parse command line options.
no fit printf nominal n $STRING awk v X
std::string getFilename(const std::string &file_name)
Get file name part, i.e. part after last JEEP::PATHNAME_SEPARATOR if any.
Definition: JeepToolkit.hh:128
int j
Definition: JPolint.hh:792
data_type v[N+1][M+1]
Definition: JPolint.hh:866
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
int debug
debug level