Jpp 20.0.0-27-g39925593c-D
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
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 */
34int 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[] = {
154 DIRECT_LIGHT_FROM_MUON,
155 SCATTERED_LIGHT_FROM_MUON,
156 DIRECT_LIGHT_FROM_EMSHOWERS,
157 SCATTERED_LIGHT_FROM_EMSHOWERS,
158 DIRECT_LIGHT_FROM_DELTARAYS,
159 SCATTERED_LIGHT_FROM_DELTARAYS
160 };
161
162 const double TTS = 2.0; // [ns]
163 const int numberOfPoints = 25;
164 const double epsilon = 1.0e-10;
165
166 const int NUMBER_OF_PDFS = sizeof(type)/sizeof(type[0]);
167
168 JPDF_t pdf[NUMBER_OF_PDFS]; // PDF
169
170 const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
171
172 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
173
174 try {
175
176 const string file_name = getFilename(fileDescriptor, type[i]);
177
178 NOTICE("loading PDF from file " << file_name << "... " << flush);
179
180 pdf[i].load(file_name.c_str());
181
182 NOTICE("OK" << endl);
183
184 pdf[i].setExceptionHandler(supervisor);
185 pdf[i].blur(TTS, numberOfPoints, epsilon);
186 }
187 catch(const JException& error) {
188 FATAL(error.what() << endl);
189 }
190 }
191
192 const double Tmin = -15.0; // [ns]
193 const double Tmax = +250.0; // [ns]
194 const double dt = 1.0; // [ns]
195
196 for (int ix = 1; ix <= h0m.GetNbinsX(); ++ix) {
197
198 const double R = h0m.GetBinCenter(ix);
199
200 double V = 0.0;
201
202 for (vector<JAngle3D>::const_iterator pmt = PMT.begin(); pmt != PMT.end(); ++pmt) {
203
204 const double theta = pi.constrain(pmt->getTheta());
205 const double phi = pi.constrain(fabs(pmt->getPhi()));
206
207 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
208
209 double yi = pdf[i](R, theta, phi, Tmax).V;
210
211 if (is_bremsstrahlung(type[i])) {
212 yi *= E;
213 } else if (is_deltarays(type[i])) {
214 yi *= JDeltaRays::getEnergyLossFromMuon(E);
215 }
216
217 V += yi;
218 }
219
220 h0m.SetBinContent(ix, 1.0 - TMath::PoissonI(0,V));
221 }
222 }
223
224 const int NUMBER_OF_PMTS = PMT.size();
225
226 double Vi[NUMBER_OF_PMTS]; // integrals
227
228 for (int ix = 1; ix <= h1m.GetNbinsX(); ++ix) {
229
230 const double R = h1m.GetBinCenter(ix);
231
232 double Y = 0.0;
233
234 for (double x = Tmin; x <= Tmax; x += dt) {
235
236 double V = 0.0;
237
238 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
239
240 Vi[pmt] = 0.0;
241
242 const double theta = pi.constrain(PMT[pmt].getTheta());
243 const double phi = pi.constrain(fabs(PMT[pmt].getPhi()));
244
245 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
246
247 double yi[] = {
248 pdf[i](R, theta, phi, x).v,
249 pdf[i](R, theta, phi, x + TMax_ns).v
250 };
251
252 if (is_bremsstrahlung(type[i])) {
253 yi[0] *= E;
254 yi[1] *= E;
255 } else if (is_deltarays(type[i])) {
256 yi[0] *= JDeltaRays::getEnergyLossFromMuon(E);
257 yi[1] *= JDeltaRays::getEnergyLossFromMuon(E);
258 }
259
260 Vi[pmt] += yi[1] - yi[0];
261 }
262
263 V += Vi[pmt];
264 }
265
266 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
267
268 const double theta = pi.constrain(PMT[pmt].getTheta());
269 const double phi = pi.constrain(fabs(PMT[pmt].getPhi()));
270
271 double y = 0.0;
272
273 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
274
275 double yi = pdf[i](R, theta, phi, x).f;
276
277 if (is_bremsstrahlung(type[i])) {
278 yi *= E;
279 } else if (is_deltarays(type[i])) {
280 yi *= JDeltaRays::getEnergyLossFromMuon(E);
281 }
282
283 y += yi;
284 }
285
286 if (y > 0.0) {
287 Y += y * (1.0 - TMath::PoissonI(0, V - Vi[pmt])) * dt;
288 }
289 }
290 }
291
292 h1m.SetBinContent(ix, 1.0 - TMath::PoissonI(0,Y));
293 }
294 }
295
296
297 {
298 typedef JSplineFunction1S_t JFunction1D_t;
299 typedef JMAPLIST<JPolint1FunctionalMap,
300 JPolint1FunctionalMap,
301 JPolint1FunctionalGridMap,
302 JPolint1FunctionalGridMap>::maplist JMapList_t;
303 typedef JPDFTable<JFunction1D_t, JMapList_t> JPDF_t;
304
305 /**
306 * PDF types.
307 */
308 const JPDFType_t type[] = {
309 DIRECT_LIGHT_FROM_EMSHOWER,
310 SCATTERED_LIGHT_FROM_EMSHOWER
311 };
312
313 const double TTS = 2.0; // [ns]
314 const int numberOfPoints = 25;
315 const double epsilon = 1.0e-10;
316
317 const int NUMBER_OF_PDFS = sizeof(type)/sizeof(type[0]);
318
319 JPDF_t pdf[NUMBER_OF_PDFS]; // PDF
320
321 const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
322
323 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
324
325 try {
326
327 const string file_name = getFilename(fileDescriptor, type[i]);
328
329 NOTICE("loading PDF from file " << file_name << "... " << flush);
330
331 pdf[i].load(file_name.c_str());
332
333 NOTICE("OK" << endl);
334
335 pdf[i].setExceptionHandler(supervisor);
336 pdf[i].blur(TTS, numberOfPoints, epsilon);
337 }
338 catch(const JException& error) {
339 FATAL(error.what() << endl);
340 }
341 }
342
343 const double Tmin = -15.0; // [ns]
344 const double Tmax = +250.0; // [ns]
345 const double dt = 1.0; // [ns]
346
347 for (int ix = 1; ix <= h0s.GetNbinsX(); ++ix) {
348
349 const double D = h0s.GetBinCenter(ix);
350
351 double V = 0.0;
352
353 for (vector<JAngle3D>::const_iterator pmt = PMT.begin(); pmt != PMT.end(); ++pmt) {
354
355 const double theta = pi.constrain(pmt->getTheta());
356 const double phi = pi.constrain(fabs(pmt->getPhi()));
357
358 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
359
360 double yi = pdf[i](D, cd, theta, phi, Tmax).V;
361
362 yi *= E;
363
364 V += yi;
365 }
366 }
367
368 h0s.SetBinContent(ix, 1.0 - TMath::PoissonI(0,V));
369 }
370
371
372 const int NUMBER_OF_PMTS = PMT.size();
373
374 double Vi[NUMBER_OF_PMTS]; // integrals
375
376 for (int ix = 1; ix <= h1s.GetNbinsX(); ++ix) {
377
378 const double D = h1s.GetBinCenter(ix);
379
380 double Y = 0.0;
381
382 for (double x = Tmin; x <= Tmax; x += dt) {
383
384 double V = 0.0;
385
386 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
387
388 Vi[pmt] = 0.0;
389
390 const double theta = pi.constrain(PMT[pmt].getTheta());
391 const double phi = pi.constrain(fabs(PMT[pmt].getPhi()));
392
393 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
394
395 double yi[] = {
396 pdf[i](D, cd, theta, phi, x).v,
397 pdf[i](D, cd, theta, phi, x + TMax_ns).v
398 };
399
400 yi[0] *= E;
401 yi[1] *= E;
402
403 Vi[pmt] += yi[1] - yi[0];
404 }
405
406 V += Vi[pmt];
407 }
408
409 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
410
411 const double theta = pi.constrain(PMT[pmt].getTheta());
412 const double phi = pi.constrain(fabs(PMT[pmt].getPhi()));
413
414 double y = 0.0;
415
416 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
417
418 double yi = pdf[i](D, cd, theta, phi, x).f;
419
420 yi *= E;
421
422 y += yi;
423 }
424
425 if (y > 0.0) {
426 Y += y * (1.0 - TMath::PoissonI(0, V - Vi[pmt])) * dt;
427 }
428 }
429 }
430
431 h1s.SetBinContent(ix, 1.0 - TMath::PoissonI(0,Y));
432 }
433 }
434
435
436 out.Write();
437 out.Close();
438}
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:72
main
int main(int argc, char **argv)
Definition JMultiPMT.cc:34
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
JLANG::JException
General exception.
Definition JException.hh:24
JLANG::JException::what
virtual const char * what() const override
Get error message.
Definition JException.hh:64
JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition JOscProbInterpolator.hh:71
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JMATH::zero
static const JZero zero
Function object to assign zero value.
Definition JZero.hh:105
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45
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.10.0