Jpp  18.6.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JTurbot.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <limits>
5 #include <map>
6 #include <vector>
7 
8 #include "TROOT.h"
9 #include "TFile.h"
10 #include "TH1D.h"
11 #include "TF1.h"
12 #include "TMath.h"
13 
14 #include "JROOT/JMinimizer.hh"
15 
19 
20 #include "JDAQ/JDAQEventIO.hh"
21 #include "JDAQ/JDAQTimesliceIO.hh"
22 #include "JDAQ/JDAQEvaluator.hh"
23 
24 #include "JDetector/JDetector.hh"
27 #include "JTrigger/JBuildL0.hh"
28 #include "JTrigger/JBuildL1.hh"
29 #include "JLang/JSinglePointer.hh"
30 #include "JROOT/JManager.hh"
31 #include "JTools/JWeight.hh"
35 #include "JSupport/JTreeScanner.hh"
37 #include "JSupport/JSupport.hh"
38 #include "JSupport/JMeta.hh"
39 
40 #include "Jeep/JPrint.hh"
41 #include "Jeep/JParser.hh"
42 #include "Jeep/JMessage.hh"
43 
44 namespace {
45 
46  using JTOOLS::JWeight;
47 
48  // status of module
49 
50  enum JStatus_t {
51  DEFAULT = 0, //!< module exists in detector
52  NODATA = -2, //!< module not present in data or no entries
53  OUT_SYNC = -3, //!< module is out-of-sync
54  ERROR = -4, //!< module is out-of-sync, possibly multiple time offsets
55  IN_SYNC = +1 //!< module is in-sync
56  };
57 
58  /**
59  * Get estimated number of background events.
60  *
61  * \param sn weight of signal
62  * \param bg weight of background
63  * \return background
64  */
65  double getBackground(const JWeight& sn, const JWeight& bg)
66  {
67  return std::max(bg.getTotal(), 1.0) * (double) sn.getCount() / (double) bg.getCount();
68  }
69 }
70 
71 /**
72  * \file
73  *
74  * Example program to search for correlations between triggered events and time slice data.\n
75  *
76  * This application can be used to identify optical modules that are out-of-sync with respect to the central White Rabbit clock.\n
77  * To this end, time correlations are searched for between triggered events and local coincidences, referred to as L1 hits.
78  * The local coincidences are independently obtained from time slice data and can thereby cross the time slice border of a triggered event.\n
79  * The MODULE_OUT_OF_SYNC of status the modules that are found to be out-of-sync are set.
80  *
81  * The number of time slices to look for L1 hits away from the triggered event can be set via command line option <tt>-N <number of time slices></tt>.\n
82  * The larger this value, the bigger time offsets can be covered.
83  * To identify the modules that are out-of-sync as fast as possible (without determining the time offset),
84  * this value can be set to 0 with the caveat that a out-of-sync occurring during the data taking run may then be missed.
85  *
86  * \author mdejong
87  */
88 int main(int argc, char **argv)
89 {
90  using namespace std;
91  using namespace JPP;
92  using namespace KM3NETDAQ;
93 
94  JSingleFileScanner<> inputFile;
95  JLimit_t& numberOfEvents = inputFile.getLimit();
96  string outputFile;
97  string detectorFile;
98  bool overwriteDetector;
99  double TMaxLocal_ns;
100  int numberOfTimeslices;
101  double binWidth_ns;
102  double deadTime_us;
103  double Pmin;
104  JROOTClassSelector selector;
105  int qaqc;
106  int debug;
107 
108  try {
109 
110  JParser<> zap("Example program to search for correlations between triggered events and timeslice data.");
111 
112  zap['f'] = make_field(inputFile, "input file.");
113  zap['o'] = make_field(outputFile, "optional output file containing ROOT histograms.") = "";
114  zap['a'] = make_field(detectorFile, "detector file.");
115  zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with module (PMT) status.");
116  zap['n'] = make_field(numberOfEvents) = JLimit::max();
117  zap['T'] = make_field(TMaxLocal_ns, "time window for local coincidences (L1),") = 20.0;
118  zap['N'] = make_field(numberOfTimeslices, "time slice difference between triggered event and L1.") = 100;
119  zap['W'] = make_field(binWidth_ns, "bin width for output histograms." ) = 10e3;
120  zap['D'] = make_field(deadTime_us, "L1 dead time (us)") = 200.0;
121  zap['p'] = make_field(Pmin, "minimal probability for background to be signal.") = 1.0e-7;
122  zap['C'] = make_field(selector, "data type selection (default is all).") = "", getROOTClassSelection<JDAQTimesliceTypes_t>();
123  zap['Q'] = make_field(qaqc) = 0;
124  zap['d'] = make_field(debug) = 2;
125 
126  zap(argc, argv);
127  }
128  catch(const exception& error) {
129  FATAL(error.what() << endl);
130  }
131 
132 
134 
135  try {
136  load(detectorFile, detector);
137  }
138  catch(const JException& error) {
139  FATAL(error);
140  }
141 
142  if (detector.empty()) {
143  FATAL("Empty detector " << detectorFile << endl);
144  }
145 
146  const JDAQHitRouter router(detector);
147 
148  const double TMax_ns = max(binWidth_ns, getMaximalTime(detector)); // signal time window
149  const double BOOST = 20.0; // scaling of time window for background
150  const double deadTime_ns = deadTime_us * 1e3;
151 
152  NOTICE("Time window " << FIXED(7,1) << TMax_ns << " [ns]" << endl);
153 
154  typedef double hit_type;
155 
156  JBuildL0<hit_type> buildL0;
159 
160  typedef JManager<int, TH1D> JManager_t;
161 
162  const Double_t xmin = -(numberOfTimeslices + 1) * getFrameTime() - 0.5*binWidth_ns;
163  const Double_t xmax = +(numberOfTimeslices + 1) * getFrameTime() + 0.5*binWidth_ns;
164  const Int_t nx = (Int_t) ((xmax - xmin) / binWidth_ns);
165 
166  JManager_t manager(new TH1D("M_%", NULL, nx, xmin, xmax));
167 
169 
171 
173 
174  if (selector == "") {
175 
176  if ((ps = new JTreeScanner< JAssertConversion<JDAQTimesliceSN, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
177  (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL2, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
178  (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL1, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
179  (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL0, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
180  (ps = new JTreeScanner< JAssertConversion<JDAQTimeslice, JDAQTimeslice> >(inputFile))->getEntries() != 0) {
181  } else {
182  FATAL("No timeslice data." << endl);
183  }
184 
185  NOTICE("Selected class " << ps->getClassName() << endl);
186 
187  } else {
188 
189  ps = zmap[selector];
190 
191  ps->configure(inputFile);
192  }
193 
194  ps->setLimit(inputFile.getLimit());
195 
196  typedef map<int, vector<double> > map_type; // frame index -> event times
197 
198  map_type buffer;
199 
200  for (JTreeScanner<JDAQEvent> in(inputFile.getFilename()); in.hasNext(); ) {
201 
202  STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
203 
204  JDAQEvent* event = in.next();
205 
206  // Average time of triggered hits
207 
208  double t0 = 0.0;
209 
210  for (JDAQEvent::const_iterator<JDAQTriggeredHit> hit = event->begin<JDAQTriggeredHit>(); hit != event->end<JDAQTriggeredHit>(); ++hit) {
211  t0 += getTime(*hit, router.getPMT(*hit));
212  }
213 
214  t0 /= event->size<JDAQTriggeredHit>();
215  t0 += event->getFrameIndex() * getFrameTime();
216 
217  buffer[event->getFrameIndex()].push_back(t0);
218  }
219  STATUS(endl);
220 
221 
222  while (ps->hasNext()) {
223 
224  STATUS("event: " << setw(10) << ps->getCounter() << '\r'); DEBUG(endl);
225 
226  JDAQTimeslice* timeslice = ps->next();
227 
228  map_type::const_iterator p = buffer.lower_bound(timeslice->getFrameIndex() - numberOfTimeslices);
229  map_type::const_iterator q = buffer.upper_bound(timeslice->getFrameIndex() + numberOfTimeslices);
230 
231  int number_of_events = 0;
232 
233  for (map_type::const_iterator i = p; i != q; ++i) {
234  number_of_events += i->second.size();
235  }
236 
237  if (number_of_events == 0) {
238  continue;
239  }
240 
241  for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
242 
243  data.clear();
244 
245  buildL1(*frame, router.getModule(frame->getModuleID()), back_inserter(data));
246 
247  TH1D* h1 = manager[frame->getModuleID()];
248 
249  double t1 = numeric_limits<double>::lowest();
250 
251  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
252 
253  const double t2 = *hit + frame->getFrameIndex() * getFrameTime();
254 
255  if (t2 > t1 + deadTime_ns) {
256 
257  for (map_type::const_iterator i = p; i != q; ++i) {
258  for (map_type::mapped_type::const_iterator t0 = i->second.begin(); t0 != i->second.end(); ++t0) {
259  h1->Fill(t2 - *t0);
260  }
261  }
262 
263  t1 = t2;
264  }
265  }
266  }
267  }
268  STATUS(endl);
269 
270 
271  // fit function
272 
273  TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
274 
275  string option = "L";
276 
277  if (debug < debug_t) {
278  option += "Q";
279  }
280 
281 
282  map<int, JStatus_t> status;
283 
284  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
285 
286  if (module->getFloor() != 0) {
287 
288  status[module->getID()] = DEFAULT;
289 
290  JManager_t::iterator p = manager.find(module->getID());
291 
292  if (p == manager.end() || p->second->GetEntries() == 0) {
293 
294  status[module->getID()] = NODATA;
295 
296  continue;
297  }
298 
299  TH1D* h1 = p->second;
300 
301  // start values
302 
303  Double_t ymax = 0.0;
304  Double_t x0 = 0.0; // [ns]
305  Double_t sigma = 250.0; // [ns]
306  Double_t ymin = 0.0;
307 
308  for (int i = 1; i != h1->GetNbinsX(); ++i) {
309 
310  const Double_t x = h1->GetBinCenter (i);
311  const Double_t y = h1->GetBinContent(i);
312 
313  if (y > ymax) {
314  ymax = y;
315  x0 = x;
316  }
317 
318  ymin += y;
319  }
320 
321  ymin /= h1->GetNbinsX();
322 
323  f1.SetParameter(0, ymax);
324  f1.SetParameter(1, x0);
325  if (binWidth_ns < sigma)
326  f1.SetParameter(2, sigma);
327  else
328  f1.FixParameter(2, binWidth_ns/sqrt(12.0));
329  f1.SetParameter(3, ymin);
330 
331  for (Int_t i = 0; i != f1.GetNpar(); ++i) {
332  f1.SetParError(i, 0.0);
333  }
334 
335  // fit
336 
337  h1->Fit(&f1, option.c_str(), "same", x0 - 5 * sigma, x0 + 5 * sigma);
338 
339 
340  if (fabs(f1.GetParameter(1)) <= 0.5*getFrameTime() && // position
341  f1.GetParameter(0) >= f1.GetParameter(3)) { // S/N
342 
343  status[module->getID()] = IN_SYNC;
344  }
345 
346 
347  // look for peaks at time offsets equal to a multiple of frame time
348 
349  map<int, JWeight> bg; // background
350  map<int, JWeight> sn; // signal(s)
351 
352  for (Int_t i = 1, ns = -(numberOfTimeslices + 1); i <= h1->GetNbinsX(); ++i) {
353 
354  const Double_t x = h1->GetBinCenter (i);
355  const Double_t y = h1->GetBinContent(i);
356 
357  while (x > (ns + 1) * getFrameTime() - BOOST * TMax_ns) {
358  ++ns;
359  }
360 
361  if (fabs(x - ns * getFrameTime()) < TMax_ns)
362  sn[ns].put(y);
363  else if (fabs(x - ns * getFrameTime()) < BOOST * TMax_ns)
364  bg[ns].put(y);
365  }
366 
367  for (map<int, JWeight>::iterator i = sn.begin(); i != sn.end(); ) {
368 
369  const Double_t y = getBackground(i->second, bg[i->first]);
370  const Double_t P = TMath::PoissonI(i->second.getTotal(), y);
371 
372  STATUS("Module " << setw(10) << module->getID()
373  << " peak at "
374  << setw(4) << i->first << " [frame time] "
375  << "S/N = "
376  << FIXED(7,1) << i->second.getTotal() << " / "
377  << FIXED(7,1) << y << ' '
378  << SCIENTIFIC(12,3) << P << endl);
379 
380  if (i->second.getTotal() > y && P < Pmin)
381  ++i; // keep
382  else
383  sn.erase(i++); // remove
384  }
385 
386  if (!(sn.size() == 1 &&
387  sn.begin()->first == 0)) {
388 
389  ERROR("Module " << setw(8) << module->getID() << " number of peaks " << sn.size() << endl);
390 
391  for (map<int, JWeight>::const_iterator i = sn.begin(); i != sn.end(); ++i) {
392 
393  ERROR("Peak at " << setw(4) << i->first << " [frame time] "
394  << "S/N = "
395  << FIXED(7,1) << i->second.getTotal() << " / "
396  << FIXED(7,1) << getBackground(i->second, bg[i->first]) << endl);
397  }
398 
399  status[module->getID()] = (sn.size() == 1 ? OUT_SYNC : ERROR);
400  }
401  }
402  }
403 
404  if (overwriteDetector) {
405 
406  detector.comment.add(JMeta(argc, argv));
407 
408  for (map<int, JStatus_t>::const_iterator i = status.begin(); i != status.end(); ++i) {
409 
410  if (i->second != IN_SYNC &&
411  i->second != NODATA) {
412 
413  NOTICE("Module " << setw(8) << i->first << " set out-of-sync." << endl);
414 
415  JModule& module = detector.getModule(router.getAddress(i->first));
416 
417  module.getStatus().set(MODULE_OUT_OF_SYNC);
418 
419  for (JModule::iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
420  pmt->getStatus().set(OUT_OF_SYNC);
421  }
422  }
423  }
424 
425  store(detectorFile.c_str(), detector);
426  }
427 
428  if (outputFile != "") {
429  manager.Write(outputFile.c_str());
430  }
431 
432  int nin = 0;
433  int nout = 0;
434 
435  for (map<int, JStatus_t>::const_iterator i = status.begin(); i != status.end(); ++i) {
436  if (i->second == IN_SYNC) {
437  ++nin;
438  }
439  if (i->second != IN_SYNC &&
440  i->second != NODATA) {
441  ++nout;
442  }
443  }
444 
445  NOTICE("Number of modules out/in sync " << nout << '/' << nin << endl);
446 
447  QAQC(nin << ' ' << nout << endl);
448 
449  return 0;
450 }
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
const double xmax
Definition: JQuadrature.cc:24
Utility class to parse command line options.
Definition: JParser.hh:1711
General exception.
Definition: JException.hh:24
Direct access to PMT data in detector data structure for DAQ hits.
debug
Definition: JMessage.hh:29
JStatus_t
Definition: JTurbot.cc:50
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
Data structure for a composite optical module.
Definition: JModule.hh:67
Auxiliary data structure to unify weights of acoustics data according to the number of pings per emit...
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:89
Template const_iterator.
Definition: JDAQEvent.hh:66
Auxiliary class to select ROOT class based on class name.
std::vector< size_t > ns
Auxialiary class to assert type conversion.
Definition: JConversion.hh:43
const JModule & getModule(const JDAQKeyHit &hit) const
Get module parameters.
Dynamic ROOT object management.
Auxiliary class for a type holder.
Definition: JType.hh:19
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
double getTime(const Hit &hit)
Get true time of hit.
string outputFile
Template definition for direct access of elements in ROOT TChain.
Data structure for detector geometry and calibration.
Simple wrapper around JModuleRouter class for direct addressing of PMT data in detector data structur...
int getFrameIndex() const
Get frame index.
const double sigma[]
Definition: JQuadrature.cc:74
The template JSinglePointer class can be used to hold a pointer to an object.
Scanning of objects from a single file according a format that follows from the extension of each fil...
const JPolynome f1(1.0, 2.0, 3.0)
Function.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
int getStatus() const
Get status.
Definition: JStatus.hh:63
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
I/O formatting auxiliaries.
Detector file.
Definition: JHead.hh:226
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
double getFrameTime()
Get frame time duration.
Definition: JDAQClock.hh:162
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
ROOT I/O of application specific meta data.
#define NOTICE(A)
Definition: JMessage.hh:64
#define ERROR(A)
Definition: JMessage.hh:66
Data time slice.
Weight calculator.
Definition: JWeight.hh:23
double getMaximalTime(const double R_Hz)
Get maximal time for given rate.
General purpose messaging.
Auxiliary class to select JTreeScanner based on ROOT class name.
Template L1 hit builder.
Definition: JBuildL1.hh:85
#define FATAL(A)
Definition: JMessage.hh:67
const double xmin
Definition: JQuadrature.cc:23
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
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.
#define QAQC(A)
QA/QC output macro.
Definition: JMessage.hh:100
static const int OUT_OF_SYNC
Enable (disable) synchronous signal from this PMT if this status bit is 0 (1);.
Definition: pmt_status.hh:17
const JModuleAddress & getAddress(const JObjectID &id) const
Get address of module.
static const int MODULE_OUT_OF_SYNC
Enable (disable) synchronous signal from this module if this status bit is 0 (1);.
Auxiliary data structure for L1 build parameters.
Definition: JBuildL1.hh:37
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84
std::map< int, range_type > map_type
Template L0 hit builder.
Definition: JBuildL0.hh:35
do set_variable DETECTOR_TXT $WORKDIR detector
KM3NeT DAQ constants, bit handling, etc.
const JPMT & getPMT(const JDAQKeyHit &hit) const
Get PMT parameters.
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:486
int qaqc
QA/QC file descriptor.
int debug
debug level
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:62
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62