Jpp - 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 
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 }
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1500
JDATABASE::PBS::PMT
static const JPBS_t PMT(3, 4, 2, 3)
PBS of photo-multiplier tube (PMT)
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWERS
direct light from EM showers
Definition: JPDFTypes.hh:32
h1
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:69
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:32
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JPHYSICS::DIRECT_LIGHT_FROM_MUON
direct light from muon
Definition: JPDFTypes.hh:29
f
then JPizza f
Definition: JPizza.sh:46
JPHYSICS::SCATTERED_LIGHT_FROM_MUON
scattered light from muon
Definition: JPDFTypes.hh:30
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWERS
scattered light from EM showers
Definition: JPDFTypes.hh:33
debug
int debug
debug level
Definition: JSirene.cc:63
JPHYSICS::JPDFType_t
JPDFType_t
PDF types.
Definition: JPDFTypes.hh:27
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:72
R
then echo Variable JPP_DIR undefined exit fi source $JPP_DIR setenv sh $JPP_DIR set_variable NORTH set_variable EAST set_variable SOUTH set_variable WEST set_variable WORKDIR tmp set_variable R set_variable CT set_variable YMAX set_variable YMIN if do_usage *then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:35
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:88
JTOOLS::v
data_type v[N+1][M+1]
Definition: JPolint.hh:740
KM3NETDAQ::NUMBER_OF_PMTS
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
Generated by   doxygen 1.8.5