Jpp  18.3.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 
17 
18 #include "JDAQ/JDAQEventIO.hh"
19 #include "JDAQ/JDAQTimesliceIO.hh"
20 #include "JDAQ/JDAQEvaluator.hh"
21 
22 #include "JDetector/JDetector.hh"
25 #include "JTrigger/JBuildL0.hh"
26 #include "JTrigger/JBuildL1.hh"
27 #include "JLang/JSinglePointer.hh"
28 #include "JROOT/JManager.hh"
29 #include "JTools/JWeight.hh"
33 #include "JSupport/JTreeScanner.hh"
35 #include "JSupport/JSupport.hh"
36 #include "JSupport/JMeta.hh"
37 
38 #include "Jeep/JPrint.hh"
39 #include "Jeep/JParser.hh"
40 #include "Jeep/JMessage.hh"
41 
42 namespace {
43 
44  using JTOOLS::JWeight;
45 
46  // status of module
47 
48  enum JStatus_t {
49  DEFAULT = 0, //!< module exists in detector
50  NODATA = -2, //!< module not present in data or no entries
51  OUT_SYNC = -3, //!< module is out-of-sync
52  ERROR = -4, //!< module is out-of-sync, possibly multiple time offsets
53  IN_SYNC = +1 //!< module is in-sync
54  };
55 
56  /**
57  * Get estimated number of background events.
58  *
59  * \param sn weight of signal
60  * \param bg weight of background
61  * \return background
62  */
63  double getBackground(const JWeight& sn, const JWeight& bg)
64  {
65  return std::max(bg.getTotal(), 1.0) * (double) sn.getCount() / (double) bg.getCount();
66  }
67 }
68 
69 /**
70  * \file
71  *
72  * Example program to search for correlations between triggered events and time slice data.\n
73  *
74  * This application can be used to identify optical modules that are out-of-sync with respect to the central White Rabbit clock.\n
75  * To this end, time correlations are searched for between triggered events and local coincidences, referred to as L1 hits.
76  * The local coincidences are independently obtained from time slice data and can thereby cross the time slice border of a triggered event.\n
77  * The MODULE_OUT_OF_SYNC of status the modules that are found to be out-of-sync are set.
78  *
79  * 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
80  * The larger this value, the bigger time offsets can be covered.
81  * To identify the modules that are out-of-sync as fast as possible (without determining the time offset),
82  * 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.
83  *
84  * \author mdejong
85  */
86 int main(int argc, char **argv)
87 {
88  using namespace std;
89  using namespace JPP;
90  using namespace KM3NETDAQ;
91 
92  JSingleFileScanner<> inputFile;
93  JLimit_t& numberOfEvents = inputFile.getLimit();
94  string outputFile;
95  string detectorFile;
96  bool overwriteDetector;
97  double TMaxLocal_ns;
98  int numberOfTimeslices;
99  double binWidth_ns;
100  double deadTime_us;
101  double Pmin;
102  JROOTClassSelector selector;
103  int qaqc;
104  int debug;
105 
106  try {
107 
108  JParser<> zap("Example program to search for correlations between triggered events and timeslice data.");
109 
110  zap['f'] = make_field(inputFile, "input file.");
111  zap['o'] = make_field(outputFile, "optional output file containing ROOT histograms.") = "";
112  zap['a'] = make_field(detectorFile, "detector file.");
113  zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with module (PMT) status.");
114  zap['n'] = make_field(numberOfEvents) = JLimit::max();
115  zap['T'] = make_field(TMaxLocal_ns, "time window for local coincidences (L1),") = 20.0;
116  zap['N'] = make_field(numberOfTimeslices, "time slice difference between triggered event and L1.") = 100;
117  zap['W'] = make_field(binWidth_ns, "bin width for output histograms." ) = 10e3;
118  zap['D'] = make_field(deadTime_us, "L1 dead time (us)") = 200.0;
119  zap['p'] = make_field(Pmin, "minimal probability for background to be signal.") = 1.0e-7;
120  zap['C'] = make_field(selector, "data type selection (default is all).") = "", getROOTClassSelection<JDAQTimesliceTypes_t>();
121  zap['Q'] = make_field(qaqc) = 0;
122  zap['d'] = make_field(debug) = 2;
123 
124  zap(argc, argv);
125  }
126  catch(const exception& error) {
127  FATAL(error.what() << endl);
128  }
129 
130 
132 
133  try {
134  load(detectorFile, detector);
135  }
136  catch(const JException& error) {
137  FATAL(error);
138  }
139 
140  if (detector.empty()) {
141  FATAL("Empty detector " << detectorFile << endl);
142  }
143 
144  const JDAQHitRouter router(detector);
145 
146  const double TMax_ns = max(binWidth_ns, getMaximalTime(detector)); // signal time window
147  const double BOOST = 20.0; // scaling of time window for background
148  const double deadTime_ns = deadTime_us * 1e3;
149 
150  NOTICE("Time window " << FIXED(7,1) << TMax_ns << " [ns]" << endl);
151 
152  typedef double hit_type;
153 
154  JBuildL0<hit_type> buildL0;
157 
158  typedef JManager<int, TH1D> JManager_t;
159 
160  const Double_t xmin = -(numberOfTimeslices + 1) * getFrameTime() - 0.5*binWidth_ns;
161  const Double_t xmax = +(numberOfTimeslices + 1) * getFrameTime() + 0.5*binWidth_ns;
162  const Int_t nx = (Int_t) ((xmax - xmin) / binWidth_ns);
163 
164  JManager_t manager(new TH1D("M_%", NULL, nx, xmin, xmax));
165 
167 
169 
171 
172  if (selector == "") {
173 
174  if ((ps = new JTreeScanner< JAssertConversion<JDAQTimesliceSN, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
175  (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL2, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
176  (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL1, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
177  (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL0, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
178  (ps = new JTreeScanner< JAssertConversion<JDAQTimeslice, JDAQTimeslice> >(inputFile))->getEntries() != 0) {
179  } else {
180  FATAL("No timeslice data." << endl);
181  }
182 
183  NOTICE("Selected class " << ps->getClassName() << endl);
184 
185  } else {
186 
187  ps = zmap[selector];
188 
189  ps->configure(inputFile);
190  }
191 
192  ps->setLimit(inputFile.getLimit());
193 
194  typedef map<int, vector<double> > map_type; // frame index -> event times
195 
196  map_type buffer;
197 
198  for (JTreeScanner<JDAQEvent> in(inputFile.getFilename()); in.hasNext(); ) {
199 
200  STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
201 
202  JDAQEvent* event = in.next();
203 
204  // Average time of triggered hits
205 
206  double t0 = 0.0;
207 
208  for (JDAQEvent::const_iterator<JDAQTriggeredHit> hit = event->begin<JDAQTriggeredHit>(); hit != event->end<JDAQTriggeredHit>(); ++hit) {
209  t0 += getTime(*hit, router.getPMT(*hit));
210  }
211 
212  t0 /= event->size<JDAQTriggeredHit>();
213  t0 += event->getFrameIndex() * getFrameTime();
214 
215  buffer[event->getFrameIndex()].push_back(t0);
216  }
217  STATUS(endl);
218 
219 
220  while (ps->hasNext()) {
221 
222  STATUS("event: " << setw(10) << ps->getCounter() << '\r'); DEBUG(endl);
223 
224  JDAQTimeslice* timeslice = ps->next();
225 
226  map_type::const_iterator p = buffer.lower_bound(timeslice->getFrameIndex() - numberOfTimeslices);
227  map_type::const_iterator q = buffer.upper_bound(timeslice->getFrameIndex() + numberOfTimeslices);
228 
229  int number_of_events = 0;
230 
231  for (map_type::const_iterator i = p; i != q; ++i) {
232  number_of_events += i->second.size();
233  }
234 
235  if (number_of_events == 0) {
236  continue;
237  }
238 
239  for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
240 
241  data.clear();
242 
243  buildL1(*frame, router.getModule(frame->getModuleID()), back_inserter(data));
244 
245  TH1D* h1 = manager[frame->getModuleID()];
246 
247  double t1 = numeric_limits<double>::lowest();
248 
249  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
250 
251  const double t2 = *hit + frame->getFrameIndex() * getFrameTime();
252 
253  if (t2 > t1 + deadTime_ns) {
254 
255  for (map_type::const_iterator i = p; i != q; ++i) {
256  for (map_type::mapped_type::const_iterator t0 = i->second.begin(); t0 != i->second.end(); ++t0) {
257  h1->Fill(t2 - *t0);
258  }
259  }
260 
261  t1 = t2;
262  }
263  }
264  }
265  }
266  STATUS(endl);
267 
268 
269  // fit function
270 
271  TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
272 
273  string option = "L";
274 
275  if (debug < debug_t) {
276  option += "Q";
277  }
278 
279 
280  map<int, JStatus_t> status;
281 
282  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
283 
284  if (module->getFloor() != 0) {
285 
286  status[module->getID()] = DEFAULT;
287 
288  JManager_t::iterator p = manager.find(module->getID());
289 
290  if (p == manager.end() || p->second->GetEntries() == 0) {
291 
292  status[module->getID()] = NODATA;
293 
294  continue;
295  }
296 
297  TH1D* h1 = p->second;
298 
299  // start values
300 
301  Double_t ymax = 0.0;
302  Double_t x0 = 0.0; // [ns]
303  Double_t sigma = 250.0; // [ns]
304  Double_t ymin = 0.0;
305 
306  for (int i = 1; i != h1->GetNbinsX(); ++i) {
307 
308  const Double_t x = h1->GetBinCenter (i);
309  const Double_t y = h1->GetBinContent(i);
310 
311  if (y > ymax) {
312  ymax = y;
313  x0 = x;
314  }
315 
316  ymin += y;
317  }
318 
319  ymin /= h1->GetNbinsX();
320 
321  f1.SetParameter(0, ymax);
322  f1.SetParameter(1, x0);
323  if (binWidth_ns < sigma)
324  f1.SetParameter(2, sigma);
325  else
326  f1.FixParameter(2, binWidth_ns/sqrt(12.0));
327  f1.SetParameter(3, ymin);
328 
329  // fit
330 
331  h1->Fit(&f1, option.c_str(), "same", x0 - 5 * sigma, x0 + 5 * sigma);
332 
333 
334  if (fabs(f1.GetParameter(1)) <= 0.5*getFrameTime() && // position
335  f1.GetParameter(0) >= f1.GetParameter(3)) { // S/N
336 
337  status[module->getID()] = IN_SYNC;
338  }
339 
340 
341  // look for peaks at time offsets equal to a multiple of frame time
342 
343  map<int, JWeight> bg; // background
344  map<int, JWeight> sn; // signal(s)
345 
346  for (Int_t i = 1, ns = -(numberOfTimeslices + 1); i <= h1->GetNbinsX(); ++i) {
347 
348  const Double_t x = h1->GetBinCenter (i);
349  const Double_t y = h1->GetBinContent(i);
350 
351  while (x > (ns + 1) * getFrameTime() - BOOST * TMax_ns) {
352  ++ns;
353  }
354 
355  if (fabs(x - ns * getFrameTime()) < TMax_ns)
356  sn[ns].put(y);
357  else if (fabs(x - ns * getFrameTime()) < BOOST * TMax_ns)
358  bg[ns].put(y);
359  }
360 
361  for (map<int, JWeight>::iterator i = sn.begin(); i != sn.end(); ) {
362 
363  const Double_t y = getBackground(i->second, bg[i->first]);
364  const Double_t P = TMath::PoissonI(i->second.getTotal(), y);
365 
366  STATUS("Module " << setw(10) << module->getID()
367  << " peak at "
368  << setw(4) << i->first << " [frame time] "
369  << "S/N = "
370  << FIXED(7,1) << i->second.getTotal() << " / "
371  << FIXED(7,1) << y << ' '
372  << SCIENTIFIC(12,3) << P << endl);
373 
374  if (i->second.getTotal() > y && P < Pmin)
375  ++i; // keep
376  else
377  sn.erase(i++); // remove
378  }
379 
380  if (!(sn.size() == 1 &&
381  sn.begin()->first == 0)) {
382 
383  ERROR("Module " << setw(8) << module->getID() << " number of peaks " << sn.size() << endl);
384 
385  for (map<int, JWeight>::const_iterator i = sn.begin(); i != sn.end(); ++i) {
386 
387  ERROR("Peak at " << setw(4) << i->first << " [frame time] "
388  << "S/N = "
389  << FIXED(7,1) << i->second.getTotal() << " / "
390  << FIXED(7,1) << getBackground(i->second, bg[i->first]) << endl);
391  }
392 
393  status[module->getID()] = (sn.size() == 1 ? OUT_SYNC : ERROR);
394  }
395  }
396  }
397 
398  if (overwriteDetector) {
399 
400  detector.comment.add(JMeta(argc, argv));
401 
402  for (map<int, JStatus_t>::const_iterator i = status.begin(); i != status.end(); ++i) {
403 
404  if (i->second != IN_SYNC &&
405  i->second != NODATA) {
406 
407  NOTICE("Module " << setw(8) << i->first << " set out-of-sync." << endl);
408 
409  JModule& module = detector.getModule(router.getAddress(i->first));
410 
411  module.getStatus().set(MODULE_OUT_OF_SYNC);
412 
413  for (JModule::iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
414  pmt->getStatus().set(OUT_OF_SYNC);
415  }
416  }
417  }
418 
419  store(detectorFile.c_str(), detector);
420  }
421 
422  if (outputFile != "") {
423  manager.Write(outputFile.c_str());
424  }
425 
426  int nin = 0;
427  int nout = 0;
428 
429  for (map<int, JStatus_t>::const_iterator i = status.begin(); i != status.end(); ++i) {
430  if (i->second == IN_SYNC) {
431  ++nin;
432  }
433  if (i->second != IN_SYNC &&
434  i->second != NODATA) {
435  ++nout;
436  }
437  }
438 
439  NOTICE("Number of modules out/in sync " << nout << '/' << nin << endl);
440 
441  QAQC(nin << ' ' << nout << endl);
442 
443  return 0;
444 }
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:1514
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:48
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
std::map< int, buffer_type > map_type
string -&gt; hits
Definition: JPerth.cc:68
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:1989
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.
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
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.
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 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
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