Jpp
 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 /**
110  * \file
111  *
112  * Auxiliary program to plot L1 time difference correlations between DOMs.
113  *
114  * \author kmelis, lnauta
115  */
116 int main(int argc, char **argv)
117 {
118  using namespace std;
119  using namespace MONITORL1DT;
120  using namespace JPP;
121 
123  string outputFile;
124  string detectorFile;
125  double Timewindow_ns;
126  int binwidth; // since the data is taken on a ns level, having bins of non-integer value is non-sensical
127  double TmaxL1_ns;
128  JROOTClassSelector selector;
129  unsigned int multiplicity;
130  bool correct_time;
131  double livetime_s;
132  int debug;
133 
134  try {
135 
136  JParser<> zap("Program to create L1 hit time difference histograms from raw data.");
137 
138  zap['f'] = make_field(inputFile, "input file");
139  zap['o'] = make_field(outputFile, "output file") = "monitor.root";
140  zap['a'] = make_field(detectorFile, "detector file");
141  zap['t'] = make_field(TmaxL1_ns, "max time between L1 hits [ns]") = 1000.0;
142  zap['T'] = make_field(Timewindow_ns, "time window around t=0 [ns]") = 2400.0;
143  zap['w'] = make_field(binwidth, "binwidth [ns]") = 1;
144  zap['C'] = make_field(selector, "datastream selector") = getROOTClassSelection<JDAQTimesliceTypes_t>();
145  zap['m'] = make_field(multiplicity, "minimal multiplicity of the L1 hits") = 2;
146  zap['c'] = make_field(correct_time, "subtract expected arrival time from delta-t");
147  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)
148  zap['d'] = make_field(debug) = 1;
149 
150  if (zap.read(argc, argv) != 0)
151  return 1;
152  }
153  catch(const exception &error) {
154  FATAL(error.what() << endl);
155  }
156 
157  cout.tie(&cerr);
158 
160 
161  try {
162  load(detectorFile, detector);
163  }
164  catch(const JException& error) {
165  FATAL(error);
166  }
167 
168  if (detector.empty()) {
169  FATAL("Empty detector." << endl);
170  }
171 
172  const JModuleRouter router(detector);
173 
174  typedef vector<JHitL1> JFrameL1_t;
175 
176  const double ctMin = -1; // full angle over the DOM
177  const JBuildL2<JHitL1> buildL2(JL2Parameters(multiplicity, TmaxL1_ns, ctMin));
178 
179  // Manager histograms for single-DOM output
180  vector<JHistogram> zmap(detector.size());
181 
182  const int nx = detector.size();
183  const double xmin = -0.5;
184  const double xmax = nx - 0.5;
185 
186  const double ymin = -floor(Timewindow_ns) + 0.5;
187  const double ymax = +floor(Timewindow_ns) + 0.5; // changed - to + to get an even number for Nbins, rebin needs numbers that have divisors.
188  const int ny = (int) ((ymax - ymin) / binwidth);
189 
190  // Manage histograms for single-DU output
191  Int_t num_of_floors = getNumberOfFloors(detector);
192  JCombinatorics c(num_of_floors);
193  Int_t npairs = c.getNumberOfPairs();
195 
196  JManager<int, TH2D>* manager;
197  manager = new JManager <int, TH2D>(new TH2D("h%", "", npairs, 0.5, npairs+0.5, ny, ymin, ymax));
198 
199  NOTICE("Running JMonitorL1dt: Monitoring L1 time differences and creating histograms." << endl);
200  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
201 
202  const JModuleAddress& address = router.getAddress(module->getID());
203 
204  STATUS("Booking histograms for module " << module->getID() << endl);
205 
206  const JString title(module->getID());
207  JString titleString1D;
208  JString titleString2D;
209  titleString1D = title + ".1L";
210  titleString2D = title + ".2S";
211 
212  zmap[address.first] = JHistogram(new TH2D((titleString2D).c_str(), NULL, nx, xmin, xmax, ny, ymin, ymax),
213  new TH1D((titleString1D).c_str(), NULL, nx, xmin, xmax));
214 
215  for (JDetector::iterator mod = detector.begin(); mod != detector.end(); ++mod) {
216  zmap[address.first].h2s->GetXaxis()->SetBinLabel(distance(detector.begin(), mod)+1, Form("%i", mod->getID()));
217  zmap[address.first].h1l->GetXaxis()->SetBinLabel(distance(detector.begin(), mod)+1, Form("%i", mod->getID()));
218  }
219  }
220 
221 
222  JTreeScanner<JDAQSummaryslice, JDAQEvaluator> summaryFile(inputFile);
223 
225 
226  int counter = 0;
227 
228  for ( ; in.hasNext() && counter != inputFile.getLimit(); ++counter) {
229 
230  STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
231 
232  const JDAQTimeslice* timeslice = in.next();
233 
234  JFrameL1_t frameL1;
236  vector<bool> DOM_OK(detector.size(), true);
237 
238  for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) {
239  if (router.hasModule(super_frame->getModuleID()) && !super_frame->empty()) {
240 
241  const JModuleAddress& address = router.getAddress(super_frame->getModuleID());
242  const JModule& module = detector.getModule(address);
243 
244  frameL1.clear();
245 
246  buildL2(*super_frame, module, back_inserter(frameL1));
247 
248  for (JFrameL1_t::iterator L1hit = frameL1.begin(); L1hit != frameL1.end(); ++L1hit) {
249  buffer.push_back(JElement(address.first, L1hit->begin()->getT()));
250  }
251  }
252  }
253 
254  for (vector<JHistogram>::iterator h1 = zmap.begin(); h1 != zmap.end(); ++h1) {
255  if (!DOM_OK[distance(zmap.begin(), h1)]) {
256  continue;
257  }
258  for (unsigned int i = 0; i < detector.size(); ++i) {
259  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
260  }
261  }
262 
263  // fill histograms with correlations
264  sort(buffer.begin(), buffer.end());
265 
266  for (vector<JElement>::const_iterator p = buffer.begin(); p != buffer.end(); ) {
268 
269  const JModule& module_p = detector.getModule((JModuleAddress)p->id);
270 
271  while (++q != buffer.end() && q->t - p->t <= Timewindow_ns ) {
272  // if (p->id == q->id) { continue; }
273 
274  const JModule& module_q = detector.getModule((JModuleAddress)q->id);
275 
276  double dom_distance = getDistance(module_p.getPosition(), module_q.getPosition());
277  double time_correction = (correct_time ? (dom_distance / getSpeedOfLight()) : 0);
278 
279  zmap[p->id].h2s->Fill(q->id, q->t - p->t - time_correction);
280  zmap[q->id].h2s->Fill(p->id, p->t - q->t + time_correction);
281 
282  if (module_p.getString() == module_q.getString()) {
283  // indices are pairs running from 0-17 between 0-152, while floors are 1-18 and bins are 1-153
284  int xbin = c.getIndex(module_p.getFloor() - 1, module_q.getFloor() - 1) + 1;
285  (*manager)[module_p.getString()]->Fill(xbin, q->t - p->t - time_correction);
286  }
287  }
288  p++; // move to the next L1 hit
289  }
290  buffer.clear();
291  }
292  STATUS(endl);
293 
294  // livetime in case set manually
295  if (livetime_s > 0.0) {
296  for (vector<JHistogram>::iterator i = zmap.begin(); i != zmap.end(); ++i) {
297  TH1D* hl = i->h1l;
298  for (int ibin = 1; ibin <= hl->GetNbinsX(); ++ibin) {
299  hl->SetBinContent(ibin, livetime_s);
300  hl->SetBinError(ibin, 0.0000001);
301  }
302  }
303  }
304 
305  TFile out(outputFile.c_str(), "recreate");
306 
307  for (vector<JHistogram>::iterator i = zmap.begin(); i != zmap.end(); ++i) {
308  i->h2s->Write();
309  i->h1l->Write();
310  }
311 
312  manager->Write(out);
313  out.Write();
314  out.Close();
315 }
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
JCombinatorics::pair_type pair_type
Auxiliary class to convert pair of indices to unique index and back.
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:57
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:1534
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:80
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
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:69
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.
int first
index of module in detector data structure
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:196
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
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.
int debug
debug level
Definition: JSirene.cc:63
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.
JElement(const int __id, const double __t)
Definition: JMonitorL1dt.cc:44
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.
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
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY source JAcoustics sh $DETECTOR_ID typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:36
unsigned int getNumberOfPairs() const
Get number of pairs.
source $JPP_DIR setenv csh $JPP_DIR eval JShellParser o a A
Basic data structure for L1 hit.
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
Template definition of histogram object interface.
Definition: JHistogram.hh:27
int main(int argc, char *argv[])
Definition: Main.cpp:15