Jpp  18.6.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JMultiPMT.cc File Reference

Auxiliary program to determine L0 and L1 hit probability as a function of. More...

#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TMath.h"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JRange.hh"
#include "JPhysics/JPDFTransformer.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JGeometry3D/JAngle3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JRotation3D.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Auxiliary program to determine L0 and L1 hit probability as a function of.

Definition in file JMultiPMT.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

PDF types.

PDF types.

Definition at line 34 of file JMultiPMT.cc.

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 }
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1711
E
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
JDATABASE::PBS::PMT
static const JPBS_t PMT(3, 4, 2, 3)
PBS of photo-multiplier tube (PMT)
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWER
scattered light from EM shower
Definition: JPDFTypes.hh:38
f
o $QUALITY_ROOT d $DEBUG!CHECK_EXIT_CODE JPlot1D f
Definition: JDataQuality.sh:76
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWERS
direct light from EM showers
Definition: JPDFTypes.hh:29
JMATH::zero
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
Y
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
Definition: detector-XY:fit1d.sh:33
V
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JPHYSICS::DIRECT_LIGHT_FROM_MUON
direct light from muon
Definition: JPDFTypes.hh:26
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
type
then set_variable PMT_FILE set_variable DAQ_FILE set_variable OUTPUT_FILE set_variable DETECTOR else fatal Wrong number of arguments fi JPrintTree f $DAQ_FILE type
Definition: JRandomEventRunByRun.sh:46
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::SCATTERED_LIGHT_FROM_MUON
scattered light from muon
Definition: JPDFTypes.hh:27
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS
scattered light from delta-rays
Definition: JPDFTypes.hh:33
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER
direct light from EM shower
Definition: JPDFTypes.hh:37
JMATH::PI
static const double PI
Mathematical constants.
Definition: JMath/JConstants.hh:20
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWERS
scattered light from EM showers
Definition: JPDFTypes.hh:30
JPHYSICS::JPDFType_t
JPDFType_t
PDF types.
Definition: JPDFTypes.hh:24
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
JPHYSICS::DIRECT_LIGHT_FROM_DELTARAYS
direct light from delta-rays
Definition: JPDFTypes.hh:32
JPARSER::not_initialised
Empty structure for specification of parser element that is not initialised (i.e. does require input)...
Definition: JParser.hh:90
R
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
numberOfPoints
int numberOfPoints
Definition: JResultPDF.cc:22
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
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
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
D
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
GAUSS_LEGENDRE::epsilon
const double epsilon
Definition: JQuadrature.cc:21
debug
int debug
debug level
Definition: archive-put-wiki-detectors.sh:92
Generated by   doxygen 1.8.5