Jpp
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 "TRandom3.h"
11 #include "TMath.h"
12 
13 #include "JTools/JFunction1D_t.hh"
14 #include "JTools/JFunctionalMap_t.hh"
15 #include "JPhysics/JPDFTransformer.hh"
16 #include "JPhysics/JPDFTable.hh"
17 #include "JPhysics/JPDFTypes.hh"
18 #include "JGeometry3D/JAngle3D.hh"
19 #include "JGeometry3D/JDirection3D.hh"
20 #include "JGeometry3D/JRotation3D.hh"
21 #include "Jeep/JParser.hh"
22 #include "Jeep/JMessage.hh"
23 
24 
25 /**
26  * \file
27  *
28  * Auxiliary program to determine L0 and L1 hit probability as a function of
29  * the minimal distance of approach of a muon to a PMT.
30  * \author mdejong
31  */
32 int main(int argc, char **argv)
33 {
34  using namespace std;
35  using namespace JPP;
36 
37  string fileDescriptor;
38  string outputFile;
39  string option;
40  double E_GeV;
41  JAngle3D dir;
42  double TMax_ns;
43  int debug;
44 
45  try {
46 
47  JParser<> zap;
48 
49  zap['f'] = make_field(fileDescriptor);
50  zap['o'] = make_field(outputFile) = "multipmt.root";
51  zap['O'] = make_field(option) = "KM3NeT", "Antares";
52  zap['E'] = make_field(E_GeV);
53  zap['D'] = make_field(dir);
54  zap['T'] = make_field(TMax_ns) = 20.0; // L1 coincidence time window [ns]
55  zap['d'] = make_field(debug) = 0;
56 
57  zap['O'] = JPARSER::not_initialised();
58 
59  zap(argc, argv);
60  }
61  catch(const exception &error) {
62  FATAL(error.what() << endl);
63  }
64 
65 
66 
67  typedef JSplineFunction1S_t JFunction1D_t;
68  typedef
69  JMapList<JPolint1FunctionalMap,
70  JMapList<JPolint1FunctionalGridMap,
71  JMapList<JPolint1FunctionalGridMap> > > JMapList_t;
72  typedef JPDFTable<JFunction1D_t, JMapList_t> JPDF_t;
73 
74  /**
75  * PDF types.
76  */
77  const JPDFType_t pdf_t[] = { DIRECT_LIGHT_FROM_MUON,
78  SCATTERED_LIGHT_FROM_MUON,
79  DIRECT_LIGHT_FROM_EMSHOWERS,
80  SCATTERED_LIGHT_FROM_EMSHOWERS };
81 
82  const double TTS = 2.0; // [ns]
83  const int numberOfPoints = 25;
84  const double epsilon = 1.0e-10;
85 
86  const int NUMBER_OF_PDFS = sizeof(pdf_t)/sizeof(pdf_t[0]);
87 
88  JPDF_t pdf[NUMBER_OF_PDFS]; // PDF
89 
90  const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
91 
92  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
93 
94  try {
95 
96  const string file_name = getFilename(fileDescriptor, pdf_t[i]);
97 
98  NOTICE("loading PDF from file " << file_name << "... " << flush);
99 
100  pdf[i].load(file_name.c_str());
101 
102  NOTICE("OK" << endl);
103 
104  pdf[i].setExceptionHandler(supervisor);
105  pdf[i].blur(TTS, numberOfPoints, epsilon);
106  }
107  catch(const JException& error) {
108  FATAL(error.what() << endl);
109  }
110  }
111 
112 
113  // set-up
114 
115  vector<JAngle3D> PMT;
116 
117  if (option == "KM3NeT") {
118 
119  PMT.push_back(JAngle3D(3.142, 0.000)); // 1
120  PMT.push_back(JAngle3D(2.582, 0.524));
121  PMT.push_back(JAngle3D(2.582, 1.571));
122  PMT.push_back(JAngle3D(2.582, 2.618));
123  PMT.push_back(JAngle3D(2.582, 3.665));
124  PMT.push_back(JAngle3D(2.582, 4.712));
125  PMT.push_back(JAngle3D(2.582, 5.760));
126  PMT.push_back(JAngle3D(2.162, 0.000));
127  PMT.push_back(JAngle3D(2.162, 1.047));
128  PMT.push_back(JAngle3D(2.162, 2.094)); // 10
129  PMT.push_back(JAngle3D(2.162, 3.142));
130  PMT.push_back(JAngle3D(2.162, 4.189));
131  PMT.push_back(JAngle3D(2.162, 5.236));
132  PMT.push_back(JAngle3D(1.872, 0.524));
133  PMT.push_back(JAngle3D(1.872, 1.571));
134  PMT.push_back(JAngle3D(1.872, 2.618));
135  PMT.push_back(JAngle3D(1.872, 3.665));
136  PMT.push_back(JAngle3D(1.872, 4.712));
137  PMT.push_back(JAngle3D(1.872, 5.760));
138  PMT.push_back(JAngle3D(1.270, 0.000)); // 20
139  PMT.push_back(JAngle3D(1.270, 1.047));
140  PMT.push_back(JAngle3D(1.270, 2.094));
141  PMT.push_back(JAngle3D(1.270, 3.142));
142  PMT.push_back(JAngle3D(1.270, 4.189));
143  PMT.push_back(JAngle3D(1.270, 5.236));
144  PMT.push_back(JAngle3D(0.980, 0.524));
145  PMT.push_back(JAngle3D(0.980, 1.571));
146  PMT.push_back(JAngle3D(0.980, 2.618));
147  PMT.push_back(JAngle3D(0.980, 3.665));
148  PMT.push_back(JAngle3D(0.980, 4.712)); // 30
149  PMT.push_back(JAngle3D(0.980, 5.760));
150 
151  } else if (option == "Antares") {
152 
153  PMT.push_back(JAngle3D(2.36, +1.05));
154  PMT.push_back(JAngle3D(2.36, 3.14));
155  PMT.push_back(JAngle3D(2.36, -1.05));
156 
157  } else {
158 
159  FATAL("Fatal error at detector option " << option << endl);
160  };
161 
162 
163  // rotate PMTs according specified orientation
164 
165  const JRotation3D R(dir);
166 
167  for (vector<JAngle3D>::iterator i = PMT.begin(); i != PMT.end(); ++i) {
168  *i = JDirection3D(*i).rotate(R);
169  }
170 
171 
172  TFile out(outputFile.c_str(), "recreate");
173 
174 
175  TH1D h0("L0", NULL, 300, 1.0, 151.0);
176  TH1D h1("L1", NULL, 300, 1.0, 151.0);
177 
178 
179  for (int i = 1; i <= h0.GetNbinsX(); ++i) {
180 
181  const double R = h0.GetBinCenter(i);
182 
183  double V = 0.0;
184 
185  for (vector<JAngle3D>::const_iterator pmt = PMT.begin(); pmt != PMT.end(); ++pmt) {
186 
187  V += (pdf[0](R, pmt->getTheta(), abs(pmt->getPhi()), 0.0).V +
188  pdf[1](R, pmt->getTheta(), abs(pmt->getPhi()), 0.0).V +
189  pdf[2](R, pmt->getTheta(), abs(pmt->getPhi()), 0.0).V * E_GeV +
190  pdf[3](R, pmt->getTheta(), abs(pmt->getPhi()), 0.0).V * E_GeV);
191  }
192 
193  h0.SetBinContent(i, 1.0 - TMath::PoissonI(0,V));
194  }
195 
196 
197  const int NUMBER_OF_PMTS = PMT.size();
198 
199  double Vi[NUMBER_OF_PMTS]; // integrals
200 
201  const double Tmin = -15.0; // [ns]
202  const double Tmax = +250.0; // [ns]
203  const double dt = 1.0; // [ns]
204 
205 
206  for (int i = 1; i <= h1.GetNbinsX(); ++i) {
207 
208  const double R = h1.GetBinCenter(i);
209 
210  double Y = 0.0;
211 
212  for (double x = Tmin; x <= Tmax; x += dt) {
213 
214  double V = 0.0;
215 
216  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
217 
218  Vi[pmt] = ((pdf[0](R, PMT[pmt].getTheta(), abs(PMT[pmt].getPhi()), x + TMax_ns).v +
219  pdf[1](R, PMT[pmt].getTheta(), abs(PMT[pmt].getPhi()), x + TMax_ns).v +
220  pdf[2](R, PMT[pmt].getTheta(), abs(PMT[pmt].getPhi()), x + TMax_ns).v * E_GeV +
221  pdf[3](R, PMT[pmt].getTheta(), abs(PMT[pmt].getPhi()), x + TMax_ns).v * E_GeV)
222 
223  -
224 
225  (pdf[0](R, PMT[pmt].getTheta(), abs(PMT[pmt].getPhi()), x).v +
226  pdf[1](R, PMT[pmt].getTheta(), abs(PMT[pmt].getPhi()), x).v +
227  pdf[2](R, PMT[pmt].getTheta(), abs(PMT[pmt].getPhi()), x).v * E_GeV +
228  pdf[3](R, PMT[pmt].getTheta(), abs(PMT[pmt].getPhi()), x).v * E_GeV));
229 
230  if (Vi[pmt] < 0.0) {
231  Vi[pmt] = 0.0;
232  }
233 
234  V += Vi[pmt];
235  }
236 
237  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
238 
239  const double y = (pdf[0](R, PMT[pmt].getTheta(), abs(PMT[pmt].getPhi()), x).f +
240  pdf[1](R, PMT[pmt].getTheta(), abs(PMT[pmt].getPhi()), x).f +
241  pdf[2](R, PMT[pmt].getTheta(), abs(PMT[pmt].getPhi()), x).f * E_GeV +
242  pdf[3](R, PMT[pmt].getTheta(), abs(PMT[pmt].getPhi()), x).f * E_GeV);
243 
244  if (y > 0.0) {
245  Y += y * (1.0 - TMath::PoissonI(0, V - Vi[pmt])) * dt;
246  }
247  }
248  }
249 
250  h1.SetBinContent(i, 1.0 - TMath::PoissonI(0,Y));
251  }
252 
253  out.Write();
254  out.Close();
255 }
JMessage.hh
JDirection3D.hh
numberOfPoints
int numberOfPoints
Definition: JResultPDF.cc:22
std::vector< JAngle3D >
KM3NETDAQ::NUMBER_OF_PMTS
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1493
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
PMT
JDAQPMTIdentifier PMT
Command line options.
Definition: JDAQTimesliceSelector.cc:35
JPDFTransformer.hh
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
JFunction1D_t.hh
debug
int debug
debug level
Definition: JSirene.cc:59
main
int main(int argc, char **argv)
Definition: JMultiPMT.cc:32
JFunctionalMap_t.hh
JAngle3D.hh
JEEP::getFilename
std::string getFilename(const std::string &file_name)
Get file name part, i.e.
Definition: JeepToolkit.hh:86
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWERS
scattered light from EM showers
Definition: JPDFTypes.hh:33
JParser.hh
JTOOLS::v
data_type v[N+1][M+1]
Definition: JPolint.hh:707
JPHYSICS::JPDFType_t
JPDFType_t
PDF types.
Definition: JPDFTypes.hh:27
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
std
Definition: jaanetDictionary.h:36
JRotation3D.hh
JPHYSICS::DIRECT_LIGHT_FROM_MUON
direct light from muon
Definition: JPDFTypes.hh:29
JMATH::zero
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:94
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWERS
direct light from EM showers
Definition: JPDFTypes.hh:32
JPDFTypes.hh
JPDFTable.hh
JPARSER::not_initialised
Empty structure for specification of parser element that is not initialised (i.e.
Definition: JParser.hh:69
JPHYSICS::SCATTERED_LIGHT_FROM_MUON
scattered light from muon
Definition: JPDFTypes.hh:30
Generated by   doxygen 1.8.16