Jpp  15.0.0-rc.2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JTurbotIter.cc
Go to the documentation of this file.
1 
2 
3 #include <string>
4 #include <iostream>
5 #include <iomanip>
6 #include <limits>
7 #include <map>
8 #include <vector>
9 
10 #include "TROOT.h"
11 #include "TFile.h"
12 #include "TH1D.h"
13 #include "TF1.h"
14 #include "TMath.h"
15 
17 #include "JDAQ/JDAQEventIO.hh"
18 #include "JDAQ/JDAQTimesliceIO.hh"
19 #include "JDAQ/JDAQEvaluator.hh"
20 #include "JDetector/JDetector.hh"
24 #include "JTrigger/JBuildL0.hh"
25 #include "JTrigger/JBuildL1.hh"
26 #include "JLang/JSinglePointer.hh"
27 #include "JROOT/JManager.hh"
28 #include "JTools/JWeight.hh"
32 #include "JSupport/JTreeScanner.hh"
34 #include "JSupport/JSupport.hh"
35 #include "JSupport/JMeta.hh"
36 
37 #include "Jeep/JPrint.hh"
38 #include "Jeep/JParser.hh"
39 #include "Jeep/JMessage.hh"
40 
41 namespace {
42 
43  // status of module
44  enum JStatus_t {
45  DEFAULT = 0, //!< module exists in detector
46  NODATA = -2, //!< module not present in data or no entries
47  OUT_SYNC = -3, //!< module is out-of-sync
48  ERROR = -4, //!< module is out-of-sync, possibly multiple time offsets
49  IN_SYNC = +1 //!< module is in-sync
50  };
51 }
52 
53 
54 /**
55  * \file
56  *
57  * Example program to search for correlations between triggered events and timeslice data.
58  * \author mdejong
59  */
60 int main(int argc, char **argv)
61 {
62  using namespace std;
63  using namespace JPP;
64  using namespace KM3NETDAQ;
65 
66  JSingleFileScanner<> inputFile;
67  JLimit_t& numberOfEvents = inputFile.getLimit();
68  string outputFile;
69  string detectorFile;
70  double TMaxLocal_ns;
71  string pmtFile;
72  int numberOfTimeslices;
73  double binWidth_ns;
74  double Pmin;
75  JROOTClassSelector selector;
76  int qaqc;
77  int debug;
78 
79  string peaks;
80 
81  try {
82 
83  JParser<> zap("Example program to search for correlations between triggered events and timeslice data.");
84 
85  zap['f'] = make_field(inputFile);
86  zap['o'] = make_field(outputFile) = "";
87  zap['a'] = make_field(detectorFile);
88  zap['n'] = make_field(numberOfEvents) = JLimit::max();
89  zap['T'] = make_field(TMaxLocal_ns) = 20.0;
90  zap['P'] = make_field(pmtFile) = "";
91  zap['N'] = make_field(numberOfTimeslices) = 100;
92  zap['W'] = make_field(binWidth_ns) = 10e3;
93  zap['p'] = make_field(Pmin) = 1.0e-6;
94  zap['C'] = make_field(selector) = "", getROOTClassSelection<JDAQTimesliceTypes_t>();
95  zap['Q'] = make_field(qaqc) = 0;
96  zap['d'] = make_field(debug) = 2;
97 
98  zap['O'] = make_field(peaks) = ""; //Make a small summary file, with the list of the peaks
99 
100  zap.read(argc, argv);
101  }
102  catch(const exception& error) {
103  FATAL(error.what() << endl);
104  }
105 
106  cout.tie(&cerr);
107 
109 
110  try {
111  load(detectorFile, detector);
112  }
113  catch(const JException& error) {
114  FATAL(error);
115  }
116 
117  const JDAQHitRouter router(detector);
118 
119  const double TMax_ns = max(binWidth_ns, getMaximalTime(detector));
120 
121  typedef double hit_type;
122 
123  JBuildL0<hit_type> buildL0;
125  vector <hit_type> data;
126 
127  typedef JManager<int, TH1D> JManager_t;
128 
129  const Double_t xmin = -(numberOfTimeslices + 1) * getFrameTime() - 0.5*binWidth_ns;
130  const Double_t xmax = +(numberOfTimeslices + 1) * getFrameTime() + 0.5*binWidth_ns;
131  const Int_t nx = (Int_t) ((xmax - xmin) / binWidth_ns);
132 
133  JManager_t manager(new TH1D("M_%", NULL, nx, xmin, xmax));
134 
136 
138 
140 
141  if (selector == "") {
142 
143  if ((ps = new JTreeScanner< JAssertConversion<JDAQTimesliceSN, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
144  (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL2, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
145  (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL1, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
146  (ps = new JTreeScanner< JAssertConversion<JDAQTimeslice, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
147  (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL0, JDAQTimeslice> >(inputFile))->getEntries() != 0) {
148  } else {
149  FATAL("No timeslice data." << endl);
150  }
151 
152  NOTICE("Selected class " << ps->getClassName() << endl);
153 
154  } else {
155 
156  ps = zmap[selector];
157 
158  ps->configure(inputFile);
159  }
160 
161  /* Updated method :
162  * OOS DOMs can have an impact on the histogram of other DOMs
163  * especialy OOS DOMs on adjacent In Sync or OOS DOMs (presence of ghost peaks)
164  * To remove these ghost peaks, the idea is to iterate the method, removing
165  * at each loop the most OOS DOM, and recomputing the histograms without
166  * them (removing events where they are involved)
167  * Once there is no more OOS DOMs in the remaining DOMs, end of precedure.
168  * Steps of procedure (from J. Hofestaedt, code updated by R. Le Breton)
169  * 1) Check if any DOM shows a peak away from 0 (OOS)
170  * If no OOS -> finish procedure
171  * 2) If at least one OOS, then identify the DOM that is most OOS,
172  * ie. highest peak (other than delta time = 0)
173  * 3) Remove the most OOS DOM and calculate the histograms again but do not
174  * consider the OOS DOM (in later iterations more than one DOM is OOS and not
175  * considered for histogram filling).
176  * 4) start with 1) again
177  */
178 
179  ps->setLimit(inputFile.getLimit());
180 
181  typedef map<int, vector<double> > map_type; // frame index -> event times
182 
183  bool OOS = true; //Boolean: to continue or not the procedure
184  map_type buffer;
185 
186  map_type peaksPerDoms;
187 
188  map<int, JStatus_t> status2;
189 
190  //OOS DOMs IDs : if at least one OOS DOM, one ID is added into oosDomsId at each iteration
191  map<int, bool> oosDomsId;
192 
193  //Default value for all modules : false
194  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
195  oosDomsId[module->getID()] = false;
196  }
197 
198  do{
199  for (JTreeScanner<JDAQEvent> in(inputFile.getFilename()); in.hasNext(); ) {
200 
201  STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
202 
203  JDAQEvent* event = in.next();
204 
205  // Average time of triggered hits
206 
207  double t0 = 0.0;
208  bool test_OOS = false; //If an OOS DOM is involved in an event, remove this event
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  if(oosDomsId[hit->getModuleID()]) {
214  test_OOS = true;
215  break;
216  }
217  }
218 
219  if(test_OOS) continue; //If true, event not taken into account
220 
221  t0 /= event->size<JDAQTriggeredHit>();
222  t0 += event->getFrameIndex() * getFrameTime();
223 
224  buffer[event->getFrameIndex()].push_back(t0);
225  }
226  STATUS(endl);
227 
228  //Events : start from the beginning
229  ps->rewind();
230 
231  while (ps->hasNext()) {
232 
233  STATUS("event: " << setw(10) << ps->getCounter() << '\r'); DEBUG(endl);
234 
235  JDAQTimeslice* timeslice = ps->next();
236 
237  map_type::const_iterator p = buffer.lower_bound(timeslice->getFrameIndex() - numberOfTimeslices);
238  map_type::const_iterator q = buffer.upper_bound(timeslice->getFrameIndex() + numberOfTimeslices);
239 
240  int number_of_events = 0;
241 
242  for (map_type::const_iterator i = p; i != q; ++i) {
243  number_of_events += i->second.size();
244  }
245 
246  if (number_of_events == 0) {
247  continue;
248  }
249 
250  for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
251 
252  data.clear();
253 
254  buildL1(*frame, router.getModule(frame->getModuleID()), back_inserter(data));
255 
256  TH1D* h1 = manager[frame->getModuleID()];
257 
258  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
259 
260  const double t1 = *hit + frame->getFrameIndex() * getFrameTime();
261 
262  for (map_type::const_iterator i = p; i != q; ++i) {
263  for (map_type::mapped_type::const_iterator j = i->second.begin(); j != i->second.end(); ++j) {
264 
265  const double t0 = *j;
266 
267  h1->Fill(t1 - t0);
268  }
269  }
270  }
271  }
272  }
273 
274  buffer.clear(); //Clear the buffer for the next loop, if any
275 
276  STATUS(endl);
277 
278  Double_t most_OOS_peak = 0;
279  Int_t most_OOS_ID = 0;
280 
281  TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
282 
283  string option = "L";
284 
285  if (debug < debug_t) {
286  option += "Q";
287  }
288 
289  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
290 
291  bool status = false;
292 
293  status2[module->getID()] = DEFAULT;
294 
295  JManager_t::iterator p = manager.find(module->getID());
296 
297  if (p == manager.end() || p->second->GetEntries() == 0) {
298 
299  status2[module->getID()] = NODATA;
300 
301  continue;
302  }
303 
304  TH1D* h1 = p->second;
305 
306  // start values
307 
308  Double_t ymax = 0.0;
309  Double_t x0 = 0.0; // [ns]
310  Double_t sigma = 250.0; // [ns]
311  Double_t ymin = 0.0;
312 
313  for (int i = 1; i != h1->GetNbinsX(); ++i) {
314 
315  const Double_t x = h1->GetBinCenter (i);
316  const Double_t y = h1->GetBinContent(i);
317 
318  if (y > ymax) {
319  ymax = y;
320  x0 = x;
321  }
322 
323  ymin += y;
324  }
325 
326  ymin /= h1->GetNbinsX();
327 
328  f1.SetParameter(0, ymax);
329  f1.SetParameter(1, x0);
330  if (binWidth_ns < sigma)
331  f1.SetParameter(2, sigma);
332  else
333  f1.FixParameter(2, binWidth_ns/sqrt(12.0));
334  f1.SetParameter(3, ymin);
335 
336  // fit
337 
338  h1->Fit(&f1, option.c_str(), "same", x0 - 5 * sigma, x0 + 5 * sigma);
339 
340  status = (fabs(f1.GetParameter(1)) <= 0.5*getFrameTime() && // position
341  f1.GetParameter(0) >= f1.GetParameter(3)); // S/N
342 
343  if (status) status2[module->getID()] = IN_SYNC;
344 
345 
346  // look for peaks at offsets equal to a multiple of frame time
347  JWeight bg; // background
348  map<int, JWeight> sn; // signal(s)
349 
350  for (Int_t i = 1, ns = -(numberOfTimeslices + 1); i <= h1->GetNbinsX(); ++i) {
351 
352  const Double_t x = h1->GetBinCenter (i);
353  const Double_t y = h1->GetBinContent(i);
354 
355  while (x > (ns + 1) * getFrameTime() - TMax_ns) { ++ns; }
356 
357  if (fabs(x - ns * getFrameTime()) < TMax_ns)
358  sn[ns].put(y);
359  else
360  bg .put(y);
361  }
362 
363  DEBUG("Module " << setw(8) << module->getID() << endl);
364 
365  for (map<int, JWeight>::iterator i = sn.begin(); i != sn.end(); ) {
366 
367  const Double_t y = bg.getTotal() * i->second.getCount() / bg.getCount();
368  const Double_t P = TMath::PoissonI(i->second.getTotal(), y);
369 
370  DEBUG("Peak at " << setw(4) << i->first << " [frame time] "
371  << "S/N = "
372  << FIXED(7,1) << i->second.getTotal() << " / "
373  << FIXED(7,1) << y << ' '
374  << SCIENTIFIC(12,3) << P << endl);
375 
376  if (P < Pmin && i->second.getTotal() > y) //Avoid edge effect (second bolean)
377  ++i; // keep
378  else
379  sn.erase(i++); // remove
380  }
381 
382  if (!(sn.size() == 1 &&
383  sn.begin()->first == 0)) {
384 
385  WARNING("Module " << setw(8) << module->getID() << " number of peaks " << sn.size() << endl);
386 
387  for (map<int, JWeight>::const_iterator i = sn.begin(); i != sn.end(); ++i) {
388 
389  const Double_t noise = bg.getTotal() * i->second.getCount() / bg.getCount();
390 
391  WARNING("Peak at " << setw(4) << i->first << " [frame time] "
392  << "S/N = "
393  << FIXED(7,1) << i->second.getTotal() << " / "
394  << FIXED(7,1) << noise << endl);
395 
396  peaksPerDoms[module->getID()].push_back(sn.size());
397  peaksPerDoms[module->getID()].push_back(i->first);
398  peaksPerDoms[module->getID()].push_back(i->second.getTotal());
399  peaksPerDoms[module->getID()].push_back(noise);
400  }
401 
402  status2[module->getID()] = (sn.size() == 1 ? OUT_SYNC : ERROR);
403  }
404 
405 
406  //if DOM is OOS
407  if (!status) {
408  //if OOS peak greater than the last and not in OOS DOMs vector
409  if ((f1.GetParameter(0) > most_OOS_peak) && !(oosDomsId[module->getID()])){
410 
411  most_OOS_peak = f1.GetParameter(0);
412  most_OOS_ID = module->getID();
413  }
414  }
415  }
416 
417  if (most_OOS_ID != 0) {
418  oosDomsId[most_OOS_ID] = true;//exclude the DOM to compute the histograms
419 
420  //Reseting all histograms
421  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
422  manager[module->getID()]->Reset("ICESM");
423  }
424  peaksPerDoms.clear();
425  }
426  else{
427  OOS = false;
428  }
429  }while(OOS);
430 
431  //To manage non working and OOS DOMs
432  NOTICE("Managing non-working and OOS DOMs." << endl);
433  if (pmtFile != "") {
434 
436 
437  try {
438  parameters.load(pmtFile.c_str());
439  }
440  catch(const JException& error) {}
441 
442  for (map<int, JStatus_t>::const_iterator i = status2.begin(); i != status2.end(); ++i) {
443 
444  if (i->second != IN_SYNC) {
445 
446  NOTICE("Module " << setw(8) << i->first << " set QEs of all PMTs to zero." << endl);
447 
448  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
449  parameters[JPMTIdentifier(i->first, pmt)].QE = 0.0;
450  }
451  }
452  }
453 
454  ofstream out(pmtFile.c_str());
455 
456  parameters.comment.add(JMeta(argc, argv));
457 
458  out << parameters << endl;
459 
460  out.close();
461  }
462 
463  if (outputFile != "") {
464  manager.Write(outputFile.c_str());
465  }
466 
467  //print small summary file of the peaks
468  if (peaks != ""){
469  ofstream stream(peaks);
470 
471  stream << "#DOM NB_PEAKS TIMESLICE_SHIFT SIGNAL NOISE\n";
472 
473  for (map_type::const_iterator dom = peaksPerDoms.begin(); dom != peaksPerDoms.end(); ++dom) {
474  stream << dom->first << " ";
475  int i = 1;
476  for (map_type::mapped_type::const_iterator data = dom->second.begin(); data != dom->second.end(); ++data) {
477  stream << *data;
478  ((i%4 == 0) ? stream << "\n" : stream << " ");
479  i++;
480  }
481  }
482  stream.close();
483  }
484 
485  //QAQC
486  int nin = 0;
487  int nout = 0;
488 
489  for (map<int, JStatus_t>::const_iterator i = status2.begin(); i != status2.end(); ++i) {
490  if (i->second == IN_SYNC) {
491  ++nin;
492  }
493  if (i->second != IN_SYNC &&
494  i->second != NODATA) {
495  ++nout;
496  }
497  }
498 
499  NOTICE("Number of modules out/in sync " << nout << '/' << nin << endl);
500 
501  QAQC(nin << ' ' << nout << endl);
502 
503  return 0;
504 
505 }
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
Direct access to PMT data in detector data structure for DAQ hits.
void clear()
Clear data.
#define WARNING(A)
Definition: JMessage.hh:65
debug
Definition: JMessage.hh:29
JStatus_t
Definition: JTurbot.cc:46
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
double getTotal() const
Get total weight.
Definition: JWeight.hh:79
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:89
Template const_iterator.
Definition: JDAQEvent.hh:68
Auxiliary class to select ROOT class based on class name.
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:71
Auxialiary class to assert type conversion.
Definition: JConversion.hh:66
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.
Definition: JTreeScanner.hh:91
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.
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...
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
void put(const double w)
Put weight.
Definition: JWeight.hh:55
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:196
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
double getFrameTime()
Get frame time duration.
Definition: JDAQClock.hh:162
ROOT I/O of application specific meta data.
long long int getCount() const
Get total count.
Definition: JWeight.hh:68
#define NOTICE(A)
Definition: JMessage.hh:64
#define ERROR(A)
Definition: JMessage.hh:66
Data time slice.
Auxiliary class for map of PMT parameters.
Weight calculator.
Definition: JWeight.hh:23
int debug
debug level
Definition: JSirene.cc:63
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
void load(const char *file_name)
Load from input file.
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
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
Auxiliary data structure for L1 build parameters.
Definition: JBuildL1.hh:37
int j
Definition: JPolint.hh:666
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
JComment & add(const std::string &comment)
Add comment.
Definition: JComment.hh:100
Template L0 hit builder.
Definition: JBuildL0.hh:35
do set_variable DETECTOR_TXT $WORKDIR detector
KM3NeT DAQ constants, bit handling, etc.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
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 source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:41
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:484
int qaqc
QA/QC file descriptor.
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:62