Jpp 21.0.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
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/JDeltaRays.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
26namespace {
27
28 const char* const MUON_t = "muon";
29 const char* const SHOWER_t = "shower";
30}
31
32/**
33 * \file
34 *
35 * Auxiliary program to determine L0 and L1 hit probability as a function of
36 * - minimal distance of approach of a muon to a PMT; and
37 * - distance between vertex and PMT for shower at given angle of emission.
38 * \author mdejong
39 */
40int main(int argc, char **argv)
41{
42 using namespace std;
43 using namespace JPP;
44
45 string fileDescriptor;
46 string outputFile;
47 string option;
48 double E;
49 double cd;
50 JAngle3D dir;
51 double TMax_ns;
52 string particle;
53 int debug;
54
55 try {
56
57 JParser<> zap;
58
59 zap['F'] = make_field(fileDescriptor);
60 zap['o'] = make_field(outputFile) = "multipmt.root";
61 zap['O'] = make_field(option) = "KM3NeT", "Antares";
62 zap['E'] = make_field(E, "muon/shower energy [GeV]");
63 zap['c'] = make_field(cd, "cosine emission angle");
64 zap['D'] = make_field(dir, "(theta, phi) of PMT [rad]");
65 zap['T'] = make_field(TMax_ns, "L1 coincidence time window [ns]") = 20.0;
66 zap['P'] = make_field(particle) = MUON_t, SHOWER_t;
67 zap['d'] = make_field(debug) = 0;
68
69 zap(argc, argv);
70 }
71 catch(const exception &error) {
72 FATAL(error.what() << endl);
73 }
74
75
76 const double epsilon = 1.0e-6; // precision angle [rad]
77 const JRange<double> pi(epsilon, PI - epsilon); // constrain angle
78
79
80 // set-up
81
82 vector<JAngle3D> PMT;
83
84 if (option == "KM3NeT") {
85
86 PMT.push_back(JAngle3D(3.142, 0.000)); // 1
87 PMT.push_back(JAngle3D(2.582, 0.524));
88 PMT.push_back(JAngle3D(2.582, 1.571));
89 PMT.push_back(JAngle3D(2.582, 2.618));
90 PMT.push_back(JAngle3D(2.582, 3.665));
91 PMT.push_back(JAngle3D(2.582, 4.712));
92 PMT.push_back(JAngle3D(2.582, 5.760));
93 PMT.push_back(JAngle3D(2.162, 0.000));
94 PMT.push_back(JAngle3D(2.162, 1.047));
95 PMT.push_back(JAngle3D(2.162, 2.094)); // 10
96 PMT.push_back(JAngle3D(2.162, 3.142));
97 PMT.push_back(JAngle3D(2.162, 4.189));
98 PMT.push_back(JAngle3D(2.162, 5.236));
99 PMT.push_back(JAngle3D(1.872, 0.524));
100 PMT.push_back(JAngle3D(1.872, 1.571));
101 PMT.push_back(JAngle3D(1.872, 2.618));
102 PMT.push_back(JAngle3D(1.872, 3.665));
103 PMT.push_back(JAngle3D(1.872, 4.712));
104 PMT.push_back(JAngle3D(1.872, 5.760));
105 PMT.push_back(JAngle3D(1.270, 0.000)); // 20
106 PMT.push_back(JAngle3D(1.270, 1.047));
107 PMT.push_back(JAngle3D(1.270, 2.094));
108 PMT.push_back(JAngle3D(1.270, 3.142));
109 PMT.push_back(JAngle3D(1.270, 4.189));
110 PMT.push_back(JAngle3D(1.270, 5.236));
111 PMT.push_back(JAngle3D(0.980, 0.524));
112 PMT.push_back(JAngle3D(0.980, 1.571));
113 PMT.push_back(JAngle3D(0.980, 2.618));
114 PMT.push_back(JAngle3D(0.980, 3.665));
115 PMT.push_back(JAngle3D(0.980, 4.712)); // 30
116 PMT.push_back(JAngle3D(0.980, 5.760));
117
118 } else if (option == "Antares") {
119
120 PMT.push_back(JAngle3D(2.36, +1.05));
121 PMT.push_back(JAngle3D(2.36, 3.14));
122 PMT.push_back(JAngle3D(2.36, -1.05));
123
124 } else {
125
126 FATAL("Fatal error at detector option " << option << endl);
127 };
128
129
130 // rotate PMTs according specified orientation
131
132 const JRotation3D R(dir);
133
134 for (vector<JAngle3D>::iterator i = PMT.begin(); i != PMT.end(); ++i) {
135 *i = JDirection3D(*i).rotate(R);
136 }
137
138
139 TFile out(outputFile.c_str(), "recreate");
140
141
142 TH1D h00("L00", NULL, 300, 1.0, 151.0);
143
144 TH1D h0m("L0m", NULL, 300, 1.0, 151.0);
145 TH1D h1m("L1m", NULL, 300, 1.0, 151.0);
146
147 TH1D h0s("L0s", NULL, 300, 1.0, 151.0);
148 TH1D h1s("L1s", NULL, 300, 1.0, 151.0);
149
150
151 if (particle == MUON_t) {
152
153 typedef JSplineFunction1S_t JFunction1D_t;
154 typedef JMAPLIST<JPolint1FunctionalMap,
155 JPolint1FunctionalGridMap,
156 JPolint1FunctionalGridMap>::maplist JMapList_t;
157 typedef JPDFTable<JFunction1D_t, JMapList_t> JPDF_t;
158
159 /**
160 * PDF types.
161 */
162 const JPDFType_t type[] = {
163 DIRECT_LIGHT_FROM_MUON,
164 SCATTERED_LIGHT_FROM_MUON,
165 DIRECT_LIGHT_FROM_EMSHOWERS,
166 SCATTERED_LIGHT_FROM_EMSHOWERS,
167 DIRECT_LIGHT_FROM_DELTARAYS,
168 SCATTERED_LIGHT_FROM_DELTARAYS
169 };
170
171 const double TTS = 2.0; // [ns]
172 const int numberOfPoints = 25;
173 const double epsilon = 1.0e-10;
174
175 const int NUMBER_OF_PDFS = sizeof(type)/sizeof(type[0]);
176
177 JPDF_t pdf[NUMBER_OF_PDFS]; // PDF
178
179 const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
180
181 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
182
183 try {
184
185 const string file_name = getFilename(fileDescriptor, type[i]);
186
187 NOTICE("loading PDF from file " << file_name << "... " << flush);
188
189 pdf[i].load(file_name.c_str());
190
191 NOTICE("OK" << endl);
192
193 pdf[i].setExceptionHandler(supervisor);
194 pdf[i].blur(TTS, numberOfPoints, epsilon);
195 }
196 catch(const JException& error) {
197 FATAL(error.what() << endl);
198 }
199 }
200
201 const double Tmin = -15.0; // [ns]
202 const double Tmax = +250.0; // [ns]
203 const double dt = 1.0; // [ns]
204
205 for (int ix = 1; ix <= h00.GetNbinsX(); ++ix) {
206
207 const double R = h00.GetBinCenter(ix);
208
209 double W = 0.0;
210
211 for (vector<JAngle3D>::const_iterator pmt = PMT.begin(); pmt != PMT.end(); ++pmt) {
212
213 const double theta = pi.constrain(pmt->getTheta());
214 const double phi = pi.constrain(fabs(pmt->getPhi()));
215
216 double w = 0.0;
217
218 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
219
220 double yi = pdf[i](R, theta, phi, Tmax).v;
221
222 if (is_bremsstrahlung(type[i])) {
223 yi *= E;
224 } else if (is_deltarays(type[i])) {
225 yi *= JDeltaRays::getEnergyLossFromMuon(E);
226 }
227
228 w += yi;
229 }
230
231 W += 1.0 - TMath::PoissonI(0,w);
232 }
233
234 W /= PMT.size();
235
236 h00.SetBinContent(ix, W);
237 }
238
239 for (int ix = 1; ix <= h0m.GetNbinsX(); ++ix) {
240
241 const double R = h0m.GetBinCenter(ix);
242
243 double V = 0.0;
244
245 for (vector<JAngle3D>::const_iterator pmt = PMT.begin(); pmt != PMT.end(); ++pmt) {
246
247 const double theta = pi.constrain(pmt->getTheta());
248 const double phi = pi.constrain(fabs(pmt->getPhi()));
249
250 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
251
252 double yi = pdf[i](R, theta, phi, Tmax).v;
253
254 if (is_bremsstrahlung(type[i])) {
255 yi *= E;
256 } else if (is_deltarays(type[i])) {
257 yi *= JDeltaRays::getEnergyLossFromMuon(E);
258 }
259
260 V += yi;
261 }
262
263 h0m.SetBinContent(ix, 1.0 - TMath::PoissonI(0,V));
264 }
265 }
266
267
268 const int NUMBER_OF_PMTS = PMT.size();
269
270 double Vi[NUMBER_OF_PMTS]; // integrals
271
272 for (int ix = 1; ix <= h1m.GetNbinsX(); ++ix) {
273
274 const double R = h1m.GetBinCenter(ix);
275
276 double Y = 0.0;
277
278 for (double x = Tmin; x <= Tmax; x += dt) {
279
280 double V = 0.0;
281
282 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
283
284 Vi[pmt] = 0.0;
285
286 const double theta = pi.constrain(PMT[pmt].getTheta());
287 const double phi = pi.constrain(fabs(PMT[pmt].getPhi()));
288
289 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
290
291 double yi[] = {
292 pdf[i](R, theta, phi, x).v,
293 pdf[i](R, theta, phi, x + TMax_ns).v
294 };
295
296 if (is_bremsstrahlung(type[i])) {
297 yi[0] *= E;
298 yi[1] *= E;
299 } else if (is_deltarays(type[i])) {
300 yi[0] *= JDeltaRays::getEnergyLossFromMuon(E);
301 yi[1] *= JDeltaRays::getEnergyLossFromMuon(E);
302 }
303
304 Vi[pmt] += yi[1] - yi[0];
305 }
306
307 V += Vi[pmt];
308 }
309
310 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
311
312 const double theta = pi.constrain(PMT[pmt].getTheta());
313 const double phi = pi.constrain(fabs(PMT[pmt].getPhi()));
314
315 double y = 0.0;
316
317 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
318
319 double yi = pdf[i](R, theta, phi, x).f;
320
321 if (is_bremsstrahlung(type[i])) {
322 yi *= E;
323 } else if (is_deltarays(type[i])) {
324 yi *= JDeltaRays::getEnergyLossFromMuon(E);
325 }
326
327 y += yi;
328 }
329
330 if (y > 0.0) {
331 Y += y * (1.0 - TMath::PoissonI(0, V - Vi[pmt])) * dt;
332 }
333 }
334 }
335
336 h1m.SetBinContent(ix, 1.0 - TMath::PoissonI(0,Y));
337 }
338 }
339
340
341 if (particle == SHOWER_t) {
342
343 typedef JSplineFunction1S_t JFunction1D_t;
344 typedef JMAPLIST<JPolint1FunctionalMap,
345 JPolint1FunctionalMap,
346 JPolint1FunctionalGridMap,
347 JPolint1FunctionalGridMap>::maplist JMapList_t;
348 typedef JPDFTable<JFunction1D_t, JMapList_t> JPDF_t;
349
350 /**
351 * PDF types.
352 */
353 const JPDFType_t type[] = {
354 DIRECT_LIGHT_FROM_EMSHOWER,
355 SCATTERED_LIGHT_FROM_EMSHOWER
356 };
357
358 const double TTS = 2.0; // [ns]
359 const int numberOfPoints = 25;
360 const double epsilon = 1.0e-10;
361
362 const int NUMBER_OF_PDFS = sizeof(type)/sizeof(type[0]);
363
364 JPDF_t pdf[NUMBER_OF_PDFS]; // PDF
365
366 const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
367
368 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
369
370 try {
371
372 const string file_name = getFilename(fileDescriptor, type[i]);
373
374 NOTICE("loading PDF from file " << file_name << "... " << flush);
375
376 pdf[i].load(file_name.c_str());
377
378 NOTICE("OK" << endl);
379
380 pdf[i].setExceptionHandler(supervisor);
381 pdf[i].blur(TTS, numberOfPoints, epsilon);
382 }
383 catch(const JException& error) {
384 FATAL(error.what() << endl);
385 }
386 }
387
388 const double Tmin = -15.0; // [ns]
389 const double Tmax = +250.0; // [ns]
390 const double dt = 1.0; // [ns]
391
392 for (int ix = 1; ix <= h0s.GetNbinsX(); ++ix) {
393
394 const double D = h0s.GetBinCenter(ix);
395
396 double V = 0.0;
397
398 for (vector<JAngle3D>::const_iterator pmt = PMT.begin(); pmt != PMT.end(); ++pmt) {
399
400 const double theta = pi.constrain(pmt->getTheta());
401 const double phi = pi.constrain(fabs(pmt->getPhi()));
402
403 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
404
405 double yi = pdf[i](D, cd, theta, phi, Tmax).V;
406
407 yi *= E;
408
409 V += yi;
410 }
411 }
412
413 h0s.SetBinContent(ix, 1.0 - TMath::PoissonI(0,V));
414 }
415
416
417 const int NUMBER_OF_PMTS = PMT.size();
418
419 double Vi[NUMBER_OF_PMTS]; // integrals
420
421 for (int ix = 1; ix <= h1s.GetNbinsX(); ++ix) {
422
423 const double D = h1s.GetBinCenter(ix);
424
425 double Y = 0.0;
426
427 for (double x = Tmin; x <= Tmax; x += dt) {
428
429 double V = 0.0;
430
431 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
432
433 Vi[pmt] = 0.0;
434
435 const double theta = pi.constrain(PMT[pmt].getTheta());
436 const double phi = pi.constrain(fabs(PMT[pmt].getPhi()));
437
438 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
439
440 double yi[] = {
441 pdf[i](D, cd, theta, phi, x).v,
442 pdf[i](D, cd, theta, phi, x + TMax_ns).v
443 };
444
445 yi[0] *= E;
446 yi[1] *= E;
447
448 Vi[pmt] += yi[1] - yi[0];
449 }
450
451 V += Vi[pmt];
452 }
453
454 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
455
456 const double theta = pi.constrain(PMT[pmt].getTheta());
457 const double phi = pi.constrain(fabs(PMT[pmt].getPhi()));
458
459 double y = 0.0;
460
461 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
462
463 double yi = pdf[i](D, cd, theta, phi, x).f;
464
465 yi *= E;
466
467 y += yi;
468 }
469
470 if (y > 0.0) {
471 Y += y * (1.0 - TMath::PoissonI(0, V - Vi[pmt])) * dt;
472 }
473 }
474 }
475
476 h1s.SetBinContent(ix, 1.0 - TMath::PoissonI(0,Y));
477 }
478 }
479
480
481 out.Write();
482 out.Close();
483}
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:74
main
int main(int argc, char **argv)
Definition JMultiPMT.cc:40
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:2140
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
JLANG::JException
General exception.
Definition JException.hh:25
JLANG::JException::what
virtual const char * what() const override
Get error message.
Definition JException.hh:65
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1697
JPHYSICS::JPDFTable
Multi-dimensional PDF table for arrival time of Cherenkov light.
Definition JPDFTable.hh:44
JTOOLS::JPolintMap
Functional map with polynomial interpolation.
Definition JPolint.hh:1153
JTOOLS::JRange
Range of values.
Definition JRange.hh:42
JTOOLS::JRange::constrain
T constrain(argument_type x) const
Constrain value to range.
Definition JRange.hh:350
JTOOLS::JSplineFunction1D
Template class for spline interpolation in 1D.
Definition JSpline.hh:734
JMATH::zero
static const JZero zero
Function object to assign zero value.
Definition JZero.hh:105
JPHYSICS::JPDFType_t
JPDFType_t
PDF types.
Definition JPDFTypes.hh:24
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:48
JTOOLS::JMAPLIST
Auxiliary class for recursive map list generation.
Definition JMapList.hh:109
Generated by doxygen 1.12.0