Jpp 19.3.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/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 */
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[] = { 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:72
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
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
Range of values.
Definition JRange.hh:42
JTOOLS::JRange::constrain
T constrain(argument_type x) const
Constrain value to range.
Definition JRange.hh:350
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:43
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.12.0