Jpp  19.0.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMonitorL1dt.cc
Go to the documentation of this file.
6 #include "JDAQ/JDAQEvaluator.hh"
7 #include "JDetector/JDetector.hh"
10 #include "JTrigger/JHitL1.hh"
11 #include "JTrigger/JBuildL2.hh"
12 #include "JTrigger/JPMTSelector.hh"
14 #include "JSupport/JTreeScanner.hh"
15 #include "JSupport/JSupport.hh"
18 #include "Jeep/JParser.hh"
19 #include "JROOT/JManager.hh"
20 #include "JTools/JCombinatorics.hh"
21 
22 #include "TFile.h"
23 #include "TH1D.h"
24 #include "TH2D.h"
25 #include "TMath.h"
26 
27 #include <string>
28 #include <iostream>
29 
30 
31 using namespace JTOOLS;
32 namespace MONITORL1DT {
33 
34  /**
35  * Data structure for hit time and DOM identifier.
36  */
37  class JElement {
38  public:
40  id(0),
41  t (0.0)
42  {}
43 
44  JElement(const int __id,
45  const double __t) :
46  id(__id),
47  t (__t)
48  {}
49 
50  int id;
51  double t;
52  };
53 
54 
55  /**
56  * Sort JElement in time.
57  *
58  * \param first first element
59  * \param second second element
60  * \return true if second later than first; else false
61  */
62  inline bool operator<(const JElement& first, const JElement& second)
63  {
64  return first.t < second.t;
65  }
66 
67 
68  /**
69  * Auxiliary data structure for histogram management.
70  */
71  struct JHistogram {
72  /**
73  * Default constructor.
74  */
76  h2s(NULL),
77  h1l(NULL)
78  {}
79 
80 
81  /**
82  * Constructor.
83  *
84  * \param __h2s 2D histogram for signal
85  * \param __h1l 1D histogram for background
86  */
87  JHistogram(TH2D* __h2s,
88  TH1D* __h1l) :
89  h2s(__h2s),
90  h1l(__h1l)
91  {
92  h2s->Sumw2();
93  h1l->Sumw2();
94  }
95 
96  TH2D* h2s;
97  TH1D* h1l;
98  };
99 
100  // Taken from C. Pieterse nanobeacon code
102  inline bool comparepair(const pair_type& A, const pair_type& B) {
103  return( abs(A.first - A.second) > abs(B.first - B.second) ) ;
104  }
105 }
106 
107 
108 /**
109  * \file
110  *
111  * Auxiliary program to plot L1 time difference correlations between DOMs.
112  *
113  * \author kmelis, lnauta
114  */
115 int main(int argc, char **argv)
116 {
117  using namespace std;
118  using namespace MONITORL1DT;
119  using namespace JPP;
120 
122  string outputFile;
123  string detectorFile;
124  double Timewindow_ns;
125  int binwidth; // since the data is taken on a ns level, having bins of non-integer value is non-sensical
126  double TmaxL1_ns;
127  JROOTClassSelector selector;
128  unsigned int multiplicity;
129  bool correct_time;
130  double livetime_s;
131  int debug;
132 
133  try {
134 
135  JParser<> zap("Program to create L1 hit time difference histograms from raw data.");
136 
137  zap['f'] = make_field(inputFile, "input file");
138  zap['o'] = make_field(outputFile, "output file") = "monitor.root";
139  zap['a'] = make_field(detectorFile, "detector file");
140  zap['t'] = make_field(TmaxL1_ns, "max time between L1 hits [ns]") = 1000.0;
141  zap['T'] = make_field(Timewindow_ns, "time window around t=0 [ns]") = 2400.0;
142  zap['w'] = make_field(binwidth, "binwidth [ns]") = 1;
143  zap['C'] = make_field(selector, "datastream selector") = getROOTClassSelection<JDAQTimesliceTypes_t>();
144  zap['m'] = make_field(multiplicity, "minimal multiplicity of the L1 hits") = 2;
145  zap['c'] = make_field(correct_time, "subtract expected arrival time from delta-t");
146  zap['L'] = make_field(livetime_s, "livetime of the data, set to positive value") = -1.0; // set to positive value if you want to set the livetime manually (i.e. MC)
147  zap['d'] = make_field(debug) = 1;
148 
149  if (zap.read(argc, argv) != 0)
150  return 1;
151  }
152  catch(const exception &error) {
153  FATAL(error.what() << endl);
154  }
155 
156 
158 
159  try {
160  load(detectorFile, detector);
161  }
162  catch(const JException& error) {
163  FATAL(error);
164  }
165 
166  if (detector.empty()) {
167  FATAL("Empty detector." << endl);
168  }
169 
170  const JModuleRouter router(detector);
171 
172  typedef vector<JHitL1> JFrameL1_t;
173 
174  const double ctMin = -1; // full angle over the DOM
175  const JBuildL2<JHitL1> buildL2(JL2Parameters(multiplicity, TmaxL1_ns, ctMin));
176 
177  // Manager histograms for single-DOM output
178  vector<JHistogram> zmap(detector.size());
179 
180  const int nx = detector.size();
181  const double xmin = -0.5;
182  const double xmax = nx - 0.5;
183 
184  const double ymin = -floor(Timewindow_ns) + 0.5;
185  const double ymax = +floor(Timewindow_ns) + 0.5; // changed - to + to get an even number for Nbins, rebin needs numbers that have divisors.
186  const int ny = (int) ((ymax - ymin) / binwidth);
187 
188  // Manage histograms for single-DU output
189  Int_t num_of_floors = getNumberOfFloors(detector);
190  JCombinatorics c(num_of_floors);
191  Int_t npairs = c.getNumberOfPairs();
193 
194  JManager<int, TH2D>* manager;
195  manager = new JManager <int, TH2D>(new TH2D("h%", "", npairs, 0.5, npairs+0.5, ny, ymin, ymax));
196 
197  NOTICE("Running JMonitorL1dt: Monitoring L1 time differences and creating histograms." << endl);
198  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
199 
200  const JModuleAddress& address = router.getAddress(module->getID());
201 
202  STATUS("Booking histograms for module " << module->getID() << endl);
203 
204  const JString title(module->getID());
205  JString titleString1D;
206  JString titleString2D;
207  titleString1D = title + ".1L";
208  titleString2D = title + ".2S";
209 
210  zmap[address.first] = JHistogram(new TH2D((titleString2D).c_str(), NULL, nx, xmin, xmax, ny, ymin, ymax),
211  new TH1D((titleString1D).c_str(), NULL, nx, xmin, xmax));
212 
213  for (JDetector::iterator mod = detector.begin(); mod != detector.end(); ++mod) {
214  zmap[address.first].h2s->GetXaxis()->SetBinLabel(distance(detector.begin(), mod)+1, Form("%i", mod->getID()));
215  zmap[address.first].h1l->GetXaxis()->SetBinLabel(distance(detector.begin(), mod)+1, Form("%i", mod->getID()));
216  }
217  }
218 
219 
220  JTreeScanner<JDAQSummaryslice, JDAQEvaluator> summaryFile(inputFile);
221 
223 
224  int counter = 0;
225 
226  for ( ; in.hasNext() && counter != inputFile.getLimit(); ++counter) {
227 
228  STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
229 
230  const JDAQTimeslice* timeslice = in.next();
231 
232  JFrameL1_t frameL1;
234  vector<bool> DOM_OK(detector.size(), true);
235 
236  for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) {
237  if (router.hasModule(super_frame->getModuleID()) && !super_frame->empty()) {
238 
239  const JModuleAddress& address = router.getAddress(super_frame->getModuleID());
240  const JModule& module = detector.getModule(address);
241 
242  frameL1.clear();
243 
244  buildL2(*super_frame, module, back_inserter(frameL1));
245 
246  for (JFrameL1_t::iterator L1hit = frameL1.begin(); L1hit != frameL1.end(); ++L1hit) {
247  buffer.push_back(JElement(address.first, L1hit->begin()->getT()));
248  }
249  }
250  }
251 
252  for (vector<JHistogram>::iterator h1 = zmap.begin(); h1 != zmap.end(); ++h1) {
253  if (!DOM_OK[distance(zmap.begin(), h1)]) {
254  continue;
255  }
256  for (unsigned int i = 0; i < detector.size(); ++i) {
257  if (DOM_OK[i]) { h1->h1l->Fill(i, getFrameTime() * 1e-9); } // fill the th1 with the frametime (how long did it measure) in seconds, all weight will give total measured time
258  }
259  }
260 
261  // fill histograms with correlations
262  sort(buffer.begin(), buffer.end());
263 
264  for (vector<JElement>::const_iterator p = buffer.begin(); p != buffer.end(); ) {
266 
267  const JModule& module_p = detector.getModule((JModuleAddress)p->id);
268 
269  while (++q != buffer.end() && q->t - p->t <= Timewindow_ns ) {
270  // if (p->id == q->id) { continue; }
271 
272  const JModule& module_q = detector.getModule((JModuleAddress)q->id);
273 
274  double dom_distance = getDistance(module_p.getPosition(), module_q.getPosition());
275  double time_correction = (correct_time ? (dom_distance / getSpeedOfLight()) : 0);
276 
277  zmap[p->id].h2s->Fill(q->id, q->t - p->t - time_correction);
278  zmap[q->id].h2s->Fill(p->id, p->t - q->t + time_correction);
279 
280  if (module_p.getString() == module_q.getString()) {
281  // indices are pairs running from 0-17 between 0-152, while floors are 1-18 and bins are 1-153
282  int xbin = c.getIndex(module_p.getFloor() - 1, module_q.getFloor() - 1) + 1;
283  (*manager)[module_p.getString()]->Fill(xbin, q->t - p->t - time_correction);
284  }
285  }
286  p++; // move to the next L1 hit
287  }
288  buffer.clear();
289  }
290  STATUS(endl);
291 
292  // livetime in case set manually
293  if (livetime_s > 0.0) {
294  for (vector<JHistogram>::iterator i = zmap.begin(); i != zmap.end(); ++i) {
295  TH1D* hl = i->h1l;
296  for (int ibin = 1; ibin <= hl->GetNbinsX(); ++ibin) {
297  hl->SetBinContent(ibin, livetime_s);
298  hl->SetBinError(ibin, 0.0000001);
299  }
300  }
301  }
302 
303  TFile out(outputFile.c_str(), "recreate");
304 
305  for (vector<JHistogram>::iterator i = zmap.begin(); i != zmap.end(); ++i) {
306  i->h2s->Write();
307  i->h1l->Write();
308  }
309 
310  manager->Write(out);
311  out.Write();
312  out.Close();
313 }
const double xmax
Definition: JQuadrature.cc:24
Utility class to parse command line options.
Definition: JParser.hh:1711
General exception.
Definition: JException.hh:24
JCombinatorics::pair_type pair_type
Auxiliary class to convert pair of indices to unique index and back.
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
Wrapper class around STL string class.
Definition: JString.hh:27
int getFloor() const
Get floor number.
Definition: JLocation.hh:145
Data structure for a composite optical module.
Definition: JModule.hh:67
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
bool operator<(const Head &first, const Head &second)
Less than operator.
Definition: JHead.hh:1817
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:89
Auxiliary class to select ROOT class based on class name.
Router for direct addressing of module data in detector data structure.
JHistogram(TH2D *__h2s, TH1D *__h1l)
Constructor.
Definition: JMonitorL1dt.cc:87
Data structure for a pair of indices.
Dynamic ROOT object management.
Auxiliary class for multiplexing object iterators.
string outputFile
Data structure for detector geometry and calibration.
void sort(JComparator_t comparator)
Sort address pairs.
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
Q JDAQEvent livetime_s
Definition: JDataQuality.sh:57
int first
index of module in detector data structure
size_t getNumberOfPairs() const
Get number of pairs.
Template L2 builder.
Definition: JBuildL2.hh:45
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
Detector file.
Definition: JHead.hh:226
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
int getIndex(const int first, const int second) const
Get index of pair of indices.
Data structure for hit time and DOM identifier.
Definition: JMonitorL1dt.cc:37
double getFrameTime()
Get frame time duration.
Definition: JDAQClock.hh:162
#define NOTICE(A)
Definition: JMessage.hh:64
int getNumberOfFloors(const JDetector &detector)
Get number of floors.
Data time slice.
JHistogram()
Default constructor.
Definition: JMonitorL1dt.cc:75
virtual const pointer_type & next() override
Get next element.
Address of module in detector data structure.
virtual bool hasNext() override
Check availability of next element.
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
#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.
Data structure for L2 parameters.
$WORKDIR ev_configure_dqsimulator txt echo process $DQ_SIMULATOR $i $SOURCE_HOST[$index] csh c(setenv ROOTSYS $ROOTSYS &&source $JPP_DIR/setenv.csh $JPP_DIR &&($DQ_SIMULATOR\-u\$NAME\$\-H\$SERVER\$\-M\$LOGGER\$\-d $DEBUG</dev/null > &/dev/null &))'
JElement(const int __id, const double __t)
Definition: JMonitorL1dt.cc:44
const double xmin
Definition: JQuadrature.cc:23
int getString() const
Get string number.
Definition: JLocation.hh:134
const double getSpeedOfLight()
Get speed of light.
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.
then fatal The output file must have the wildcard in the e g root fi eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Definition: JCanberra.sh:48
Utility class to parse command line options.
const JModuleAddress & getAddress(const JObjectID &id) const
Get address of module.
bool hasModule(const JObjectID &id) const
Has module.
bool comparepair(const pair_type &A, const pair_type &B)
do set_variable DETECTOR_TXT $WORKDIR detector
KM3NeT DAQ constants, bit handling, etc.
Auxiliary data structure for histogram management.
Definition: JMonitorL1dt.cc:71
source $JPP_DIR setenv csh $JPP_DIR &dev null eval JShellParser o a A
Basic data structure for L1 hit.
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
Template definition of histogram object interface.
Definition: JHistogram.hh:27