Jpp  master_rocky
the software that should make you happy
JMultiPMT.cc
Go to the documentation of this file.
1 
2 #include <iostream>
3 #include <iomanip>
4 #include <string>
5 #include <vector>
6 
7 #include "TROOT.h"
8 #include "TFile.h"
9 #include "TH1D.h"
10 #include "TMath.h"
11 
12 #include "JTools/JFunction1D_t.hh"
13 #include "JTools/JFunctionalMap_t.hh"
14 #include "JTools/JRange.hh"
15 #include "JPhysics/JPDFTransformer.hh"
16 #include "JPhysics/JPDFTable.hh"
17 #include "JPhysics/JPDFTypes.hh"
18 #include "JPhysics/JPDFToolkit.hh"
19 #include "JGeometry3D/JAngle3D.hh"
20 #include "JGeometry3D/JDirection3D.hh"
21 #include "JGeometry3D/JRotation3D.hh"
22 #include "Jeep/JParser.hh"
23 #include "Jeep/JMessage.hh"
24 
25 
26 /**
27  * \file
28  *
29  * Auxiliary program to determine L0 and L1 hit probability as a function of
30  * - minimal distance of approach of a muon to a PMT; and
31  * - distance between vertex and PMT for shower at given angle of emission.
32  * \author mdejong
33  */
34 int main(int argc, char **argv)
35 {
36  using namespace std;
37  using namespace JPP;
38 
39  string fileDescriptor;
40  string outputFile;
41  string option;
42  double E;
43  double cd;
44  JAngle3D dir;
45  double TMax_ns;
46  int debug;
47 
48  try {
49 
50  JParser<> zap;
51 
52  zap['f'] = make_field(fileDescriptor);
53  zap['o'] = make_field(outputFile) = "multipmt.root";
54  zap['O'] = make_field(option) = "KM3NeT", "Antares";
55  zap['E'] = make_field(E, "muon/shower energy [GeV]");
56  zap['c'] = make_field(cd, "cosine emission angle");
57  zap['D'] = make_field(dir, "(theta, phi) of PMT [rad]");
58  zap['T'] = make_field(TMax_ns, "L1 coincidence time window [ns]") = 20.0;
59  zap['d'] = make_field(debug) = 0;
60 
61  zap['O'] = JPARSER::not_initialised();
62 
63  zap(argc, argv);
64  }
65  catch(const exception &error) {
66  FATAL(error.what() << endl);
67  }
68 
69 
70  const double epsilon = 1.0e-6; // precision angle [rad]
71  const JRange<double> pi(epsilon, PI - epsilon); // constrain angle
72 
73 
74  // set-up
75 
76  vector<JAngle3D> PMT;
77 
78  if (option == "KM3NeT") {
79 
80  PMT.push_back(JAngle3D(3.142, 0.000)); // 1
81  PMT.push_back(JAngle3D(2.582, 0.524));
82  PMT.push_back(JAngle3D(2.582, 1.571));
83  PMT.push_back(JAngle3D(2.582, 2.618));
84  PMT.push_back(JAngle3D(2.582, 3.665));
85  PMT.push_back(JAngle3D(2.582, 4.712));
86  PMT.push_back(JAngle3D(2.582, 5.760));
87  PMT.push_back(JAngle3D(2.162, 0.000));
88  PMT.push_back(JAngle3D(2.162, 1.047));
89  PMT.push_back(JAngle3D(2.162, 2.094)); // 10
90  PMT.push_back(JAngle3D(2.162, 3.142));
91  PMT.push_back(JAngle3D(2.162, 4.189));
92  PMT.push_back(JAngle3D(2.162, 5.236));
93  PMT.push_back(JAngle3D(1.872, 0.524));
94  PMT.push_back(JAngle3D(1.872, 1.571));
95  PMT.push_back(JAngle3D(1.872, 2.618));
96  PMT.push_back(JAngle3D(1.872, 3.665));
97  PMT.push_back(JAngle3D(1.872, 4.712));
98  PMT.push_back(JAngle3D(1.872, 5.760));
99  PMT.push_back(JAngle3D(1.270, 0.000)); // 20
100  PMT.push_back(JAngle3D(1.270, 1.047));
101  PMT.push_back(JAngle3D(1.270, 2.094));
102  PMT.push_back(JAngle3D(1.270, 3.142));
103  PMT.push_back(JAngle3D(1.270, 4.189));
104  PMT.push_back(JAngle3D(1.270, 5.236));
105  PMT.push_back(JAngle3D(0.980, 0.524));
106  PMT.push_back(JAngle3D(0.980, 1.571));
107  PMT.push_back(JAngle3D(0.980, 2.618));
108  PMT.push_back(JAngle3D(0.980, 3.665));
109  PMT.push_back(JAngle3D(0.980, 4.712)); // 30
110  PMT.push_back(JAngle3D(0.980, 5.760));
111 
112  } else if (option == "Antares") {
113 
114  PMT.push_back(JAngle3D(2.36, +1.05));
115  PMT.push_back(JAngle3D(2.36, 3.14));
116  PMT.push_back(JAngle3D(2.36, -1.05));
117 
118  } else {
119 
120  FATAL("Fatal error at detector option " << option << endl);
121  };
122 
123 
124  // rotate PMTs according specified orientation
125 
126  const JRotation3D R(dir);
127 
128  for (vector<JAngle3D>::iterator i = PMT.begin(); i != PMT.end(); ++i) {
129  *i = JDirection3D(*i).rotate(R);
130  }
131 
132 
133  TFile out(outputFile.c_str(), "recreate");
134 
135 
136  TH1D h0m("L0m", NULL, 300, 1.0, 151.0);
137  TH1D h1m("L1m", NULL, 300, 1.0, 151.0);
138 
139  TH1D h0s("L0s", NULL, 300, 1.0, 151.0);
140  TH1D h1s("L1s", NULL, 300, 1.0, 151.0);
141 
142 
143  {
144  typedef JSplineFunction1S_t JFunction1D_t;
145  typedef JMAPLIST<JPolint1FunctionalMap,
146  JPolint1FunctionalGridMap,
147  JPolint1FunctionalGridMap>::maplist JMapList_t;
148  typedef JPDFTable<JFunction1D_t, JMapList_t> JPDF_t;
149 
150  /**
151  * PDF types.
152  */
153  const JPDFType_t type[] = { DIRECT_LIGHT_FROM_MUON,
154  SCATTERED_LIGHT_FROM_MUON,
155  DIRECT_LIGHT_FROM_EMSHOWERS,
156  SCATTERED_LIGHT_FROM_EMSHOWERS,
157  DIRECT_LIGHT_FROM_DELTARAYS,
158  SCATTERED_LIGHT_FROM_DELTARAYS };
159 
160  const double TTS = 2.0; // [ns]
161  const int numberOfPoints = 25;
162  const double epsilon = 1.0e-10;
163 
164  const int NUMBER_OF_PDFS = sizeof(type)/sizeof(type[0]);
165 
166  JPDF_t pdf[NUMBER_OF_PDFS]; // PDF
167 
168  const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
169 
170  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
171 
172  try {
173 
174  const string file_name = getFilename(fileDescriptor, type[i]);
175 
176  NOTICE("loading PDF from file " << file_name << "... " << flush);
177 
178  pdf[i].load(file_name.c_str());
179 
180  NOTICE("OK" << endl);
181 
182  pdf[i].setExceptionHandler(supervisor);
183  pdf[i].blur(TTS, numberOfPoints, epsilon);
184  }
185  catch(const JException& error) {
186  FATAL(error.what() << endl);
187  }
188  }
189 
190  const double Tmin = -15.0; // [ns]
191  const double Tmax = +250.0; // [ns]
192  const double dt = 1.0; // [ns]
193 
194  for (int ix = 1; ix <= h0m.GetNbinsX(); ++ix) {
195 
196  const double R = h0m.GetBinCenter(ix);
197 
198  double V = 0.0;
199 
200  for (vector<JAngle3D>::const_iterator pmt = PMT.begin(); pmt != PMT.end(); ++pmt) {
201 
202  const double theta = pi.constrain(pmt->getTheta());
203  const double phi = pi.constrain(fabs(pmt->getPhi()));
204 
205  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
206 
207  double yi = pdf[i](R, theta, phi, Tmax).V;
208 
209  if (is_bremsstrahlung(type[i])) {
210  yi *= E;
211  } else if (is_deltarays(type[i])) {
212  yi *= getDeltaRaysFromMuon(E);
213  }
214 
215  V += yi;
216  }
217 
218  h0m.SetBinContent(ix, 1.0 - TMath::PoissonI(0,V));
219  }
220  }
221 
222  const int NUMBER_OF_PMTS = PMT.size();
223 
224  double Vi[NUMBER_OF_PMTS]; // integrals
225 
226  for (int ix = 1; ix <= h1m.GetNbinsX(); ++ix) {
227 
228  const double R = h1m.GetBinCenter(ix);
229 
230  double Y = 0.0;
231 
232  for (double x = Tmin; x <= Tmax; x += dt) {
233 
234  double V = 0.0;
235 
236  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
237 
238  Vi[pmt] = 0.0;
239 
240  const double theta = pi.constrain(PMT[pmt].getTheta());
241  const double phi = pi.constrain(fabs(PMT[pmt].getPhi()));
242 
243  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
244 
245  double yi[] = {
246  pdf[i](R, theta, phi, x).v,
247  pdf[i](R, theta, phi, x + TMax_ns).v
248  };
249 
250  if (is_bremsstrahlung(type[i])) {
251  yi[0] *= E;
252  yi[1] *= E;
253  } else if (is_deltarays(type[i])) {
254  yi[0] *= getDeltaRaysFromMuon(E);
255  yi[1] *= getDeltaRaysFromMuon(E);
256  }
257 
258  Vi[pmt] += yi[1] - yi[0];
259  }
260 
261  V += Vi[pmt];
262  }
263 
264  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
265 
266  const double theta = pi.constrain(PMT[pmt].getTheta());
267  const double phi = pi.constrain(fabs(PMT[pmt].getPhi()));
268 
269  double y = 0.0;
270 
271  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
272 
273  double yi = pdf[i](R, theta, phi, x).f;
274 
275  if (is_bremsstrahlung(type[i])) {
276  yi *= E;
277  } else if (is_deltarays(type[i])) {
278  yi *= getDeltaRaysFromMuon(E);
279  }
280 
281  y += yi;
282  }
283 
284  if (y > 0.0) {
285  Y += y * (1.0 - TMath::PoissonI(0, V - Vi[pmt])) * dt;
286  }
287  }
288  }
289 
290  h1m.SetBinContent(ix, 1.0 - TMath::PoissonI(0,Y));
291  }
292  }
293 
294 
295  {
296  typedef JSplineFunction1S_t JFunction1D_t;
297  typedef JMAPLIST<JPolint1FunctionalMap,
298  JPolint1FunctionalMap,
299  JPolint1FunctionalGridMap,
300  JPolint1FunctionalGridMap>::maplist JMapList_t;
301  typedef JPDFTable<JFunction1D_t, JMapList_t> JPDF_t;
302 
303  /**
304  * PDF types.
305  */
306  const JPDFType_t type[] = { DIRECT_LIGHT_FROM_EMSHOWER,
307  SCATTERED_LIGHT_FROM_EMSHOWER };
308 
309  const double TTS = 2.0; // [ns]
310  const int numberOfPoints = 25;
311  const double epsilon = 1.0e-10;
312 
313  const int NUMBER_OF_PDFS = sizeof(type)/sizeof(type[0]);
314 
315  JPDF_t pdf[NUMBER_OF_PDFS]; // PDF
316 
317  const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
318 
319  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
320 
321  try {
322 
323  const string file_name = getFilename(fileDescriptor, type[i]);
324 
325  NOTICE("loading PDF from file " << file_name << "... " << flush);
326 
327  pdf[i].load(file_name.c_str());
328 
329  NOTICE("OK" << endl);
330 
331  pdf[i].setExceptionHandler(supervisor);
332  pdf[i].blur(TTS, numberOfPoints, epsilon);
333  }
334  catch(const JException& error) {
335  FATAL(error.what() << endl);
336  }
337  }
338 
339  const double Tmin = -15.0; // [ns]
340  const double Tmax = +250.0; // [ns]
341  const double dt = 1.0; // [ns]
342 
343  for (int ix = 1; ix <= h0s.GetNbinsX(); ++ix) {
344 
345  const double D = h0s.GetBinCenter(ix);
346 
347  double V = 0.0;
348 
349  for (vector<JAngle3D>::const_iterator pmt = PMT.begin(); pmt != PMT.end(); ++pmt) {
350 
351  const double theta = pi.constrain(pmt->getTheta());
352  const double phi = pi.constrain(fabs(pmt->getPhi()));
353 
354  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
355 
356  double yi = pdf[i](D, cd, theta, phi, Tmax).V;
357 
358  yi *= E;
359 
360  V += yi;
361  }
362  }
363 
364  h0s.SetBinContent(ix, 1.0 - TMath::PoissonI(0,V));
365  }
366 
367 
368  const int NUMBER_OF_PMTS = PMT.size();
369 
370  double Vi[NUMBER_OF_PMTS]; // integrals
371 
372  for (int ix = 1; ix <= h1s.GetNbinsX(); ++ix) {
373 
374  const double D = h1s.GetBinCenter(ix);
375 
376  double Y = 0.0;
377 
378  for (double x = Tmin; x <= Tmax; x += dt) {
379 
380  double V = 0.0;
381 
382  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
383 
384  Vi[pmt] = 0.0;
385 
386  const double theta = pi.constrain(PMT[pmt].getTheta());
387  const double phi = pi.constrain(fabs(PMT[pmt].getPhi()));
388 
389  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
390 
391  double yi[] = {
392  pdf[i](D, cd, theta, phi, x).v,
393  pdf[i](D, cd, theta, phi, x + TMax_ns).v
394  };
395 
396  yi[0] *= E;
397  yi[1] *= E;
398 
399  Vi[pmt] += yi[1] - yi[0];
400  }
401 
402  V += Vi[pmt];
403  }
404 
405  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
406 
407  const double theta = pi.constrain(PMT[pmt].getTheta());
408  const double phi = pi.constrain(fabs(PMT[pmt].getPhi()));
409 
410  double y = 0.0;
411 
412  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
413 
414  double yi = pdf[i](D, cd, theta, phi, x).f;
415 
416  yi *= E;
417 
418  y += yi;
419  }
420 
421  if (y > 0.0) {
422  Y += y * (1.0 - TMath::PoissonI(0, V - Vi[pmt])) * dt;
423  }
424  }
425  }
426 
427  h1s.SetBinContent(ix, 1.0 - TMath::PoissonI(0,Y));
428  }
429  }
430 
431 
432  out.Write();
433  out.Close();
434 }
PMT
JDAQPMTIdentifier PMT
Command line options.
Definition: JDAQTimesliceSelector.cc:35
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JFunction1D_t.hh
JFunctionalMap_t.hh
Various implementations of functional maps.
JMessage.hh
General purpose messaging.
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
debug
int debug
debug level
Definition: JSirene.cc:69
main
int main(int argc, char **argv)
Definition: JMultiPMT.cc:34
JPDFToolkit.hh
Auxiliary methods for PDF calculations.
JPDFTypes.hh
Numbering scheme for PDF types.
JParser.hh
Utility class to parse command line options.
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2142
JRange.hh
Auxiliary class to define a range between two values.
numberOfPoints
int numberOfPoints
Definition: JResultPDF.cc:22
JGEOMETRY3D::JAngle3D
Data structure for angles in three dimensions.
Definition: JAngle3D.hh:35
JGEOMETRY3D::JDirection3D
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:35
JGEOMETRY3D::JDirection3D::rotate
JDirection3D & rotate(const JRotation3D &R)
Rotate.
Definition: JDirection3D.hh:179
JGEOMETRY3D::JRotation3D
Rotation matrix.
Definition: JRotation3D.hh:114
JLANG::JException
General exception.
Definition: JException.hh:24
JLANG::JException::what
virtual const char * what() const override
Get error message.
Definition: JException.hh:64
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1698
JPHYSICS::JPDFTable
Multi-dimensional PDF table for arrival time of Cherenkov light.
Definition: JPDFTable.hh:44
JTOOLS::JRange::constrain
T constrain(argument_type x) const
Constrain value to range.
Definition: JRange.hh:350
GAUSS_LEGENDRE::epsilon
const double epsilon
Definition: JQuadrature.cc:21
JEEP::getFilename
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
JMATH::zero
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
JMATH::PI
static const double PI
Mathematical constants.
Definition: JMath/JConstants.hh:20
JPHYSICS::getDeltaRaysFromMuon
double getDeltaRaysFromMuon(const double E, const JRange< double > T_GeV=JRange< double >(DELTARAY_TMIN, DELTARAY_TMAX))
Equivalent EM-shower energy due to delta-rays per unit muon track length.
Definition: JPDFToolkit.hh:260
JPHYSICS::is_deltarays
bool is_deltarays(const int pdf)
Test if given PDF type corresponds to Cherenkov light from delta-rays.
Definition: JPDFTypes.hh:151
JPHYSICS::is_bremsstrahlung
bool is_bremsstrahlung(const int pdf)
Test if given PDF type corresponds to Cherenkov light from Bremsstrahlung.
Definition: JPDFTypes.hh:137
JPHYSICS::JPDFType_t
JPDFType_t
PDF types.
Definition: JPDFTypes.hh:24
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWER
@ SCATTERED_LIGHT_FROM_EMSHOWER
scattered light from EM shower
Definition: JPDFTypes.hh:38
JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS
@ SCATTERED_LIGHT_FROM_DELTARAYS
scattered light from delta-rays
Definition: JPDFTypes.hh:33
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWERS
@ DIRECT_LIGHT_FROM_EMSHOWERS
direct light from EM showers
Definition: JPDFTypes.hh:29
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWERS
@ SCATTERED_LIGHT_FROM_EMSHOWERS
scattered light from EM showers
Definition: JPDFTypes.hh:30
JPHYSICS::SCATTERED_LIGHT_FROM_MUON
@ SCATTERED_LIGHT_FROM_MUON
scattered light from muon
Definition: JPDFTypes.hh:27
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER
@ DIRECT_LIGHT_FROM_EMSHOWER
direct light from EM shower
Definition: JPDFTypes.hh:37
JPHYSICS::DIRECT_LIGHT_FROM_DELTARAYS
@ DIRECT_LIGHT_FROM_DELTARAYS
direct light from delta-rays
Definition: JPDFTypes.hh:32
JPHYSICS::DIRECT_LIGHT_FROM_MUON
@ DIRECT_LIGHT_FROM_MUON
direct light from muon
Definition: JPDFTypes.hh:26
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:43
JTOOLS::v
data_type v[N+1][M+1]
Definition: JPolint.hh:866
KM3NETDAQ::NUMBER_OF_PMTS
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
std
Definition: JSTDTypes.hh:14
JPARSER::not_initialised
Empty structure for specification of parser element that is not initialised (i.e. does require input)...
Definition: JParser.hh:74
JTOOLS::JMAPLIST
Auxiliary class for recursive map list generation.
Definition: JMapList.hh:109
JTOOLS::JPolint1FunctionalGridMap
Type definition of a 1st degree polynomial interpolation based on a JGridMap implementation.
Definition: JFunctionalMap_t.hh:93
JTOOLS::JPolint1FunctionalMap
Type definition of a 1st degree polynomial interpolation based on a JMap implementation.
Definition: JFunctionalMap_t.hh:57
JTOOLS::JSplineFunction1S_t
Type definition of a spline interpolation method based on a JCollection with JResultPDF result type.
Definition: JFunction1D_t.hh:53
Generated by doxygen 1.9.1