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
JMonopole.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <iomanip>
5 #include <limits>
6 #include <vector>
7 #include <set>
8 
9 #include "TROOT.h"
10 #include "TFile.h"
11 #include "TH1D.h"
12 
16 
17 #include "JDetector/JDetector.hh"
20 
23 #include "JGeometry3D/JOmega3D.hh"
24 
25 #include "JTools/JQuantile.hh"
26 
28 #include "JSupport/JSupport.hh"
29 
30 #include "Jeep/JParser.hh"
31 #include "Jeep/JMessage.hh"
32 
33 
34 namespace {
35 
36  using JTOOLS::JQuantile;
39 
40 
41  /**
42  * Auxiliary data structure to store rate as a function of index of module in detector.
43  */
44  struct JRates :
45  public std::vector<JQuantile>
46  {
47  /**
48  * Constructor.
49  *
50  * \param router module router
51  * \param summary summary data
52  */
53  JRates(const JModuleRouter& router, const JDAQSummaryslice& summary) :
54  std::vector<JQuantile>(router.getReference().size())
55  {
56  const double factor = 1.0e-3; // [kHz]
57 
58  for (JDAQSummaryslice::const_iterator frame = summary.begin(); frame != summary.end(); ++frame) {
59 
60  if (router.hasModule(frame->getModuleID())) {
61 
62  const int index = router.getIndex(frame->getModuleID()); // index of module in detector
63 
64  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
65  (*this)[index].put(frame->getRate(pmt, factor)); // rate [kHz]
66  }
67  }
68  }
69  }
70 
71 
72  /**
73  * Constructor.
74  *
75  * \param router module router
76  * \param summary summary data
77  */
78  JRates(const JModuleRouter& router, const ExtendedSummary_TimeSlice& summary) :
79  std::vector<JQuantile>(router.getReference().size())
80  {
81  ANTARES::setClock(); // Antares DAQ clock
82 
83  const double factor = 2.0e6 / getFrameTime(); // 2 ARS / PMT [kHz]
84 
85  for (ExtendedSummary_TimeSlice::const_iterator frame = summary.begin(); frame != summary.end(); ++frame) {
86 
87  if (router.hasModule(frame->lcm_id())) {
88 
89  const int index = router.getIndex(frame->lcm_id()); // index of module in detector
90 
91  const int count = frame->numberOfItemsOrg() << 4; // decompress data
92 
93  (*this)[index].put(count * factor); // rate [kHz]
94  }
95  }
96  }
97  };
98 
99 
100  /**
101  * Main method.
102  *
103  * \param argc number of command line arguments
104  * \param argv list of command line arguments
105  * \return status
106  */
107  template<class T>
108  bool do_main(int argc, char **argv)
109  {
110  using namespace std;
111  using namespace JPP;
112  using namespace KM3NETDAQ;
113 
114  JMultipleFileScanner<T> inputFile;
115  JLimit_t& numberOfEvents = inputFile.getLimit();
116  string outputFile;
117  string detectorFile;
118  double roadWidth_m;
119  double gridAngle_deg;
120  int numberOfPMTs;
121  int debug;
122 
123  try {
124 
125  JParser<> zap("Example program to analyse summary data.");
126 
127  zap['f'] = make_field(inputFile);
128  zap['o'] = make_field(outputFile) = "monopole.root";
129  zap['a'] = make_field(detectorFile);
130  zap['n'] = make_field(numberOfEvents) = JLimit::max();
131  zap['R'] = make_field(roadWidth_m);
132  zap['G'] = make_field(gridAngle_deg);
133  zap['N'] = make_field(numberOfPMTs);
134  zap['d'] = make_field(debug) = 2;
135 
136  zap(argc, argv);
137  }
138  catch(const exception& error) {
139  FATAL(error.what() << endl);
140  }
141 
142 
143  if (roadWidth_m <= 0.0) { FATAL("Invalid road width [m] " << roadWidth_m << endl); }
144  if (gridAngle_deg <= 0.0) { FATAL("Invalid grid angle [deg] " << gridAngle_deg << endl); }
145  if (gridAngle_deg >= 90.0) { FATAL("Invalid grid angle [deg] " << gridAngle_deg << endl); }
146 
147 
149 
150  try {
151  load(detectorFile, detector);
152  }
153  catch(const JException& error) {
154  FATAL(error);
155  }
156 
157  DEBUG("Number of modules in detector " << detector.size() << endl);
158 
159  const JModuleRouter router(detector); // look-up table of module in detector
160  const JOmega3D omega(gridAngle_deg * PI/180.0); // set of assumed monopole directions
161 
162  typedef vector< set<int> > buffer1D; // index of module in detector -> indices of modules in detector
163  typedef vector< buffer1D > buffer2D; // index of direction -> index of module in detector -> indices of modules in detector
164 
165  buffer2D M2(omega.size(), buffer1D(detector.size())); // look-up table of modules within road width around given module for each direction
166 
167  for (size_t i = 0; i != omega.size(); ++i) {
168 
169  const JRotation3D R(omega[i]); // rotates z-axis to assumed monopole direction
170 
171  for (size_t m1 = 0; m1 != detector.size(); ++m1) { // index of first module in detector
172  for (size_t m2 = 0; m2 != m1; ++m2) { // index of second module in detector
173 
174  JPosition3D p1 = detector[m1].getPosition();
175  JPosition3D p2 = detector[m2].getPosition();
176 
177  p1.rotate(R);
178  p2.rotate(R);
179 
180  const double x = p2.getX() - p1.getX();
181  const double y = p2.getY() - p1.getY();
182 
183  if (sqrt(x*x + y*y) <= roadWidth_m) { // application of road width
184  M2[i][m1].insert(m2);
185  }
186  }
187  }
188  }
189 
190 
191  TFile out(outputFile.c_str(), "recreate");
192 
193  TH1D h1("h1", NULL, 1000, 0.0, 1.0e3); // average rate [kHz]
194 
195 
196  while (inputFile.hasNext()) {
197 
198  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
199 
200  const T* summary = inputFile.next();
201 
202  const JRates buffer(router, *summary); // rate as a function of index of module in detector
203 
204  for (size_t i = 0; i != omega.size(); ++i) { // scan directions
205 
206  JQuantile Q; // rate per direction
207 
208  for (size_t m1 = 0; m1 != buffer.size(); ++m1) { // index of module in detector
209 
210  Q += buffer[m1];
211 
212  const set<int>& M1 = M2[i][m1]; // indices of modules in detector
213 
214  for (set<int>::const_iterator m2 = M1.begin(); m2 != M1.end(); ++m2) {
215  Q += buffer[*m2];
216  }
217  }
218 
219  if (Q.getCount() >= numberOfPMTs) {
220  h1.Fill(Q.getMean()); // rate per PMT
221  }
222  }
223  }
224  STATUS(endl);
225 
226  out.Write();
227  out.Close();
228 
229  return (inputFile.getCounter() != 0 ? true : false);
230  }
231 }
232 
233 
234 /**
235  * \file
236  *
237  * Example program to analyse summary data.
238  * \author mdejong
239  */
240 int main(int argc, char **argv)
241 {
242  if (do_main<KM3NETDAQ::JDAQSummaryslice>(argc, argv)) { return 0; }
243  if (do_main<ExtendedSummary_TimeSlice> (argc, argv)) { return 0; }
244 }
void setClock()
Set clock.
Definition: JDAQClock.hh:297
Utility class to parse command line options.
Definition: JParser.hh:1514
General exception.
Definition: JException.hh:23
Q(UTCMax_s-UTCMin_s)-livetime_s
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
TPaveText * p1
Auxiliary data structure for running average, standard deviation and quantiles.
Definition: JQuantile.hh:43
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:89
const int getIndex(const JObjectID &id) const
Get index of module.
Router for direct addressing of module data in detector data structure.
Rotation matrix.
Definition: JRotation3D.hh:111
Direction set covering (part of) solid angle.
Definition: JOmega3D.hh:66
string outputFile
Data structure for detector geometry and calibration.
double getMean() const
Get mean value.
Definition: JQuantile.hh:252
ExtendedSummary time slices.
Definition: TimeSlice.hh:642
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Detector file.
Definition: JHead.hh:226
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
double getFrameTime()
Get frame time duration.
Definition: JDAQClock.hh:162
do set_variable OUTPUT_DIRECTORY $WORKDIR T
static const double PI
Mathematical constants.
double getY() const
Get y position.
Definition: JVector3D.hh:104
long long int getCount() const
Get total count.
Definition: JQuantile.hh:186
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Direct access to module in detector data structure.
p2
Definition: module-Z:fit.sh:74
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
bool hasModule(const JObjectID &id) const
Has module.
double getX() const
Get x position.
Definition: JVector3D.hh:94
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84
do set_variable DETECTOR_TXT $WORKDIR detector
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62