Jpp  18.0.0-rc.2
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 the minimal distance of approach of a muon to a PMT. More...

#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TRandom3.h"
#include "TMath.h"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JPhysics/JPDFTransformer.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFTypes.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 the minimal distance of approach of a muon to a PMT.

Author
mdejong

Definition in file JMultiPMT.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

PDF types.

Definition at line 32 of file JMultiPMT.cc.

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  typedef JSplineFunction1S_t JFunction1D_t;
67  typedef
68  JMapList<JPolint1FunctionalMap,
69  JMapList<JPolint1FunctionalGridMap,
70  JMapList<JPolint1FunctionalGridMap> > > JMapList_t;
71  typedef JPDFTable<JFunction1D_t, JMapList_t> JPDF_t;
72 
73  /**
74  * PDF types.
75  */
76  const JPDFType_t pdf_t[] = { DIRECT_LIGHT_FROM_MUON,
77  SCATTERED_LIGHT_FROM_MUON,
78  DIRECT_LIGHT_FROM_EMSHOWERS,
79  SCATTERED_LIGHT_FROM_EMSHOWERS };
80 
81  const double TTS = 2.0; // [ns]
82  const int numberOfPoints = 25;
83  const double epsilon = 1.0e-10;
84 
85  const int NUMBER_OF_PDFS = sizeof(pdf_t)/sizeof(pdf_t[0]);
86 
87  JPDF_t pdf[NUMBER_OF_PDFS]; // PDF
88 
89  const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
90 
91  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
92 
93  try {
94 
95  const string file_name = getFilename(fileDescriptor, pdf_t[i]);
96 
97  NOTICE("loading PDF from file " << file_name << "... " << flush);
98 
99  pdf[i].load(file_name.c_str());
100 
101  NOTICE("OK" << endl);
102 
103  pdf[i].setExceptionHandler(supervisor);
104  pdf[i].blur(TTS, numberOfPoints, epsilon);
105  }
106  catch(const JException& error) {
107  FATAL(error.what() << endl);
108  }
109  }
110 
111 
112  // set-up
113 
114  vector<JAngle3D> PMT;
115 
116  if (option == "KM3NeT") {
117 
118  PMT.push_back(JAngle3D(3.142, 0.000)); // 1
119  PMT.push_back(JAngle3D(2.582, 0.524));
120  PMT.push_back(JAngle3D(2.582, 1.571));
121  PMT.push_back(JAngle3D(2.582, 2.618));
122  PMT.push_back(JAngle3D(2.582, 3.665));
123  PMT.push_back(JAngle3D(2.582, 4.712));
124  PMT.push_back(JAngle3D(2.582, 5.760));
125  PMT.push_back(JAngle3D(2.162, 0.000));
126  PMT.push_back(JAngle3D(2.162, 1.047));
127  PMT.push_back(JAngle3D(2.162, 2.094)); // 10
128  PMT.push_back(JAngle3D(2.162, 3.142));
129  PMT.push_back(JAngle3D(2.162, 4.189));
130  PMT.push_back(JAngle3D(2.162, 5.236));
131  PMT.push_back(JAngle3D(1.872, 0.524));
132  PMT.push_back(JAngle3D(1.872, 1.571));
133  PMT.push_back(JAngle3D(1.872, 2.618));
134  PMT.push_back(JAngle3D(1.872, 3.665));
135  PMT.push_back(JAngle3D(1.872, 4.712));
136  PMT.push_back(JAngle3D(1.872, 5.760));
137  PMT.push_back(JAngle3D(1.270, 0.000)); // 20
138  PMT.push_back(JAngle3D(1.270, 1.047));
139  PMT.push_back(JAngle3D(1.270, 2.094));
140  PMT.push_back(JAngle3D(1.270, 3.142));
141  PMT.push_back(JAngle3D(1.270, 4.189));
142  PMT.push_back(JAngle3D(1.270, 5.236));
143  PMT.push_back(JAngle3D(0.980, 0.524));
144  PMT.push_back(JAngle3D(0.980, 1.571));
145  PMT.push_back(JAngle3D(0.980, 2.618));
146  PMT.push_back(JAngle3D(0.980, 3.665));
147  PMT.push_back(JAngle3D(0.980, 4.712)); // 30
148  PMT.push_back(JAngle3D(0.980, 5.760));
149 
150  } else if (option == "Antares") {
151 
152  PMT.push_back(JAngle3D(2.36, +1.05));
153  PMT.push_back(JAngle3D(2.36, 3.14));
154  PMT.push_back(JAngle3D(2.36, -1.05));
155 
156  } else {
157 
158  FATAL("Fatal error at detector option " << option << endl);
159  };
160 
161 
162  // rotate PMTs according specified orientation
163 
164  const JRotation3D R(dir);
165 
166  for (vector<JAngle3D>::iterator i = PMT.begin(); i != PMT.end(); ++i) {
167  *i = JDirection3D(*i).rotate(R);
168  }
169 
170 
171  TFile out(outputFile.c_str(), "recreate");
172 
173 
174  TH1D h0("L0", NULL, 300, 1.0, 151.0);
175  TH1D h1("L1", NULL, 300, 1.0, 151.0);
176 
177 
178  for (int i = 1; i <= h0.GetNbinsX(); ++i) {
179 
180  const double R = h0.GetBinCenter(i);
181 
182  double V = 0.0;
183 
184  for (vector<JAngle3D>::const_iterator pmt = PMT.begin(); pmt != PMT.end(); ++pmt) {
185 
186  V += (pdf[0](R, pmt->getTheta(), fabs(pmt->getPhi()), 0.0).V +
187  pdf[1](R, pmt->getTheta(), fabs(pmt->getPhi()), 0.0).V +
188  pdf[2](R, pmt->getTheta(), fabs(pmt->getPhi()), 0.0).V * E_GeV +
189  pdf[3](R, pmt->getTheta(), fabs(pmt->getPhi()), 0.0).V * E_GeV);
190  }
191 
192  h0.SetBinContent(i, 1.0 - TMath::PoissonI(0,V));
193  }
194 
195 
196  const int NUMBER_OF_PMTS = PMT.size();
197 
198  double Vi[NUMBER_OF_PMTS]; // integrals
199 
200  const double Tmin = -15.0; // [ns]
201  const double Tmax = +250.0; // [ns]
202  const double dt = 1.0; // [ns]
203 
204 
205  for (int i = 1; i <= h1.GetNbinsX(); ++i) {
206 
207  const double R = h1.GetBinCenter(i);
208 
209  double Y = 0.0;
210 
211  for (double x = Tmin; x <= Tmax; x += dt) {
212 
213  double V = 0.0;
214 
215  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
216 
217  Vi[pmt] = ((pdf[0](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x + TMax_ns).v +
218  pdf[1](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x + TMax_ns).v +
219  pdf[2](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x + TMax_ns).v * E_GeV +
220  pdf[3](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x + TMax_ns).v * E_GeV)
221 
222  -
223 
224  (pdf[0](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x).v +
225  pdf[1](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x).v +
226  pdf[2](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x).v * E_GeV +
227  pdf[3](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x).v * E_GeV));
228 
229  if (Vi[pmt] < 0.0) {
230  Vi[pmt] = 0.0;
231  }
232 
233  V += Vi[pmt];
234  }
235 
236  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
237 
238  const double y = (pdf[0](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x).f +
239  pdf[1](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x).f +
240  pdf[2](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x).f * E_GeV +
241  pdf[3](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x).f * E_GeV);
242 
243  if (y > 0.0) {
244  Y += y * (1.0 - TMath::PoissonI(0, V - Vi[pmt])) * dt;
245  }
246  }
247  }
248 
249  h1.SetBinContent(i, 1.0 - TMath::PoissonI(0,Y));
250  }
251 
252  out.Write();
253  out.Close();
254 }
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1514
JDATABASE::PBS::PMT
static const JPBS_t PMT(3, 4, 2, 3)
PBS of photo-multiplier tube (PMT)
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::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:1989
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
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
JPARSER::not_initialised
Empty structure for specification of parser element that is not initialised (i.e. does require input)...
Definition: JParser.hh:89
R
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
numberOfPoints
int numberOfPoints
Definition: JResultPDF.cc:22
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:777
KM3NETDAQ::NUMBER_OF_PMTS
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
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