Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JTurbotIter.cc File Reference

Example program to search for correlations between triggered events and timeslice data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <limits>
#include <map>
#include <vector>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TF1.h"
#include "TMath.h"
#include "JDAQ/JDAQ.hh"
#include "JDAQ/JDAQEvent.hh"
#include "JDAQ/JDAQTimeslice.hh"
#include "JDAQ/JDAQEvaluator.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDAQHitRouter.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JBuildL1.hh"
#include "JLang/JSinglePointer.hh"
#include "JGizmo/JManager.hh"
#include "JTools/JWeight.hh"
#include "JROOT/JROOTClassSelector.hh"
#include "JSupport/JSingleFileScanner.hh"
#include "JSupport/JTreeScannerInterface.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JAutoTreeScanner.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Example program to search for correlations between triggered events and timeslice data.

Author
mdejong

Definition in file JTurbotIter.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 48 of file JTurbotIter.cc.

49 {
50  using namespace std;
51  using namespace JPP;
52  using namespace KM3NETDAQ;
53 
54  JSingleFileScanner<> inputFile;
55  JLimit_t& numberOfEvents = inputFile.getLimit();
56  string outputFile;
57  string detectorFile;
58  double TMaxLocal_ns;
59  string pmtFile;
60  int numberOfTimeslices;
61  double binWidth_ns;
62  multimap<int, int> mask;
63  double Pmin;
64  JROOTClassSelector selector;
65  int debug;
66  string peaks;
67  int multiPeaks;
68 
69  try {
70 
71  JParser<> zap("Example program to search for correlations between triggered events and timeslice data.");
72 
73  zap['f'] = make_field(inputFile);
74  zap['o'] = make_field(outputFile) = "";
75  zap['a'] = make_field(detectorFile);
76  zap['n'] = make_field(numberOfEvents) = JLimit::max();
77  zap['T'] = make_field(TMaxLocal_ns) = 20.0;
78  zap['P'] = make_field(pmtFile) = "";
79  zap['N'] = make_field(numberOfTimeslices) = 100;
80  zap['W'] = make_field(binWidth_ns) = 10e3;
81  zap['M'] = make_field(mask) = JPARSER::initialised();
82  zap['p'] = make_field(Pmin) = 1.0e-7;
83  zap['C'] = make_field(selector) = "", getROOTClassSelection<JDAQTimesliceTypes_t>();
84  zap['d'] = make_field(debug) = 2;
85  zap['O'] = make_field(peaks) = ""; //Make a small summary file, with the list of the peaks
86  zap['m'] = make_field(multiPeaks) = 1; //Enables final multi-peaks analysis
87 
88  if (zap.read(argc, argv) != 0)
89  return 1;
90  }
91  catch(const exception& error) {
92  FATAL(error.what() << endl);
93  }
94 
95  cout.tie(&cerr);
96 
97 
98  map<int, int> MASK; // mask PMTs with permament high-rate veto.
99 
100  const int WR = 0x80000000; // White Rabbit
101 
102  //MASK[808969848] = 0x00000020; // floor 03, pmt 05
103  //MASK[809544061] = 0x00000080; // floor 05, pmt 07
104  //MASK[808432835] = 0x00004000; // floor 18, pmt 14
105 
106  for (multimap<int, int>::const_iterator i = mask.begin(); i != mask.end(); ++i) {
107  JBit(i->second).set(MASK[i->first]);
108  }
109 
110  for (map<int, int>::const_iterator i = MASK.begin(); i != MASK.end(); ++i) {
111  DEBUG("Mask " << setw(8) << i->first << " 0x" << hex << i->second << dec << endl);
112  }
113 
114  JDetector detector;
115 
116  try {
117  load(detectorFile, detector);
118  }
119  catch(const JException& error) {
120  FATAL(error);
121  }
122 
123  const JDAQHitRouter router(detector);
124 
125  const double TMax_ns = max(binWidth_ns, getMaximalTime(detector));
126 
127  typedef double hit_type;
128 
129  JBuildL0<hit_type> buildL0;
130  JBuildL1<hit_type> buildL1(TMaxLocal_ns, true);
131  vector <hit_type> data;
132 
133  typedef JManager<int, TH1D> JManager_t;
134 
135  const Double_t xmin = -(numberOfTimeslices + 1) * getFrameTime() - 0.5*binWidth_ns;
136  const Double_t xmax = +(numberOfTimeslices + 1) * getFrameTime() + 0.5*binWidth_ns;
137  const Int_t nx = (Int_t) ((xmax - xmin) / binWidth_ns);
138 
139  JManager_t manager(new TH1D("M_%", NULL, nx, xmin, xmax));
140 
142 
143  JSinglePointer< JTreeScannerInterface<JDAQTimeslice> > ps;
144 
145  JAutoTreeScanner<JDAQTimeslice> zmap = JType<JDAQTimesliceTypes_t>();
146 
147  if (selector == "") {
148 
149  if ((ps = new JTreeScanner< JAssertConversion<JDAQTimesliceSN, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
150  (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL2, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
151  (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL1, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
152  (ps = new JTreeScanner< JAssertConversion<JDAQTimeslice, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
153  (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL0, JDAQTimeslice> >(inputFile))->getEntries() != 0) {
154  } else {
155  FATAL("No timeslice data." << endl);
156  }
157 
158  NOTICE("Selected class " << ps->getClassName() << endl);
159 
160  } else {
161 
162  ps = zmap[selector];
163 
164  ps->configure(inputFile);
165  }
166 
167  //Vector with OOS DOMs IDs
168  vector<int> oosDomsId;
169 
170  /* Updated method :
171  * OOS DOMs can have an impact on the histogram of other DOMs
172  * especialy OOS DOMs on adjacent In Sync or OOS DOMs (presence of ghost peaks)
173  * To remove these ghost peaks, the idea is to iterate the method, removing
174  * at each loop the most OOS DOM, and recomputing the histograms without
175  * them (removing events where there are involved)
176  * Once there is no more OOS DOMs in the remaining DOMs, end of precedure.
177  * Steps of procedure (from J. Hofestaedt, code updated by R. Le Breton)
178  * 1) Check if any DOM shows a peak away from 0 (OOS)
179  * If no OOS -> finish procedure
180  * 2) If at least one OOS, then identify the DOM that is most OOS,
181  * ie. highest peak (other than delta time = 0)
182  * 3) Remove the most OOS DOM and calculate the histograms again but do not
183  * consider the OOS DOM (in later iterations more than one DOM is OOS and not
184  * considered for histogram filling).
185  * 4) start with 1) again
186  */
187 
188  typedef map<int, vector<double> > map_type; // frame index -> event times
189 
190  bool OOS = true; //Boolean: to continue or not the procedure
191  JPMTParametersMap parameters;
192  map_type buffer;
193 
194  vector<vector<double> > peaksPerDoms;
195 
196  if (pmtFile != "") {
197  try {
198  parameters.load(pmtFile.c_str());
199  }
200  catch(const JException& error) {}
201  }
202 
203  do{
204  for (JTreeScanner<JDAQEvent> in(inputFile.getFilename(), inputFile.getLimit()); in.hasNext(); ) {
205 
206  STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
207 
208  JDAQEvent* event = in.next();
209 
210  // Average time of triggered hits
211 
212  double t0 = 0.0;
213  bool test_OOS = false; //If an OOS DOM is involved in an event, remove this event
214 
215  for (JDAQEvent::const_iterator<JDAQTriggeredHit> hit = event->begin<JDAQTriggeredHit>(); hit != event->end<JDAQTriggeredHit>(); ++hit) {
216  t0 += getTime(*hit, router.getPMT(*hit));
217 
218  if(std::find(oosDomsId.begin(), oosDomsId.end(), hit->getModuleIdentifier()) != oosDomsId.end()) {
219  test_OOS = true;
220  continue;
221  }
222  }
223 
224  if(test_OOS) continue; //If true, event not taken into account
225 
226  t0 /= event->size<JDAQTriggeredHit>();
227  t0 += event->getFrameIndex() * getFrameTime();
228 
229  buffer[event->getFrameIndex()].push_back(t0);
230  }
231  STATUS(endl);
232 
233  //Events : start from the beginning
234  ps->rewind();
235 
236  while (ps->hasNext()) {
237 
238  STATUS("event: " << setw(10) << ps->getCounter() << '\r'); DEBUG(endl);
239 
240  JDAQTimeslice* timeslice = ps->next();
241 
242  map_type::const_iterator p = buffer.lower_bound(timeslice->getFrameIndex() - numberOfTimeslices);
243  map_type::const_iterator q = buffer.upper_bound(timeslice->getFrameIndex() + numberOfTimeslices);
244 
245  int number_of_events = 0;
246 
247  for (map_type::const_iterator i = p; i != q; ++i) {
248  number_of_events += i->second.size();
249  }
250 
251  if (number_of_events == 0) {
252  continue;
253  }
254 
255  for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
256 
257  // Test that none of the PMTs has high-rate veto.
258 
259  if ((frame->getStatus() & ~MASK[frame->getModuleID()] & ~WR) == 0) {
260 
261  data.clear();
262 
263  buildL1(*frame, router.getModule(frame->getModuleID()), back_inserter(data));
264 
265  TH1D* h1 = manager[frame->getModuleID()];
266 
267  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
268 
269  const double t1 = *hit + frame->getFrameIndex() * getFrameTime();
270 
271  for (map_type::const_iterator i = p; i != q; ++i) {
272  for (map_type::mapped_type::const_iterator j = i->second.begin(); j != i->second.end(); ++j) {
273 
274  const double t0 = *j;
275 
276  h1->Fill(t1 - t0);
277  }
278  }
279  }
280  }
281  }
282  }
283 
284  buffer.clear(); //Clear the buffer for the next loop, if any
285 
286  STATUS(endl);
287 
288  if (pmtFile != "") {
289 
290  Double_t most_OOS_peak = 0;
291  Int_t most_OOS_ID = 0;
292 
293  TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
294 
295  string option = "L";
296 
297  if (debug < debug_t) {
298  option += "Q";
299  }
300 
301  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
302 
303  bool status = false;
304 
305  JManager_t::iterator p = manager.find(module->getID());
306 
307  if (p != manager.end()) {
308 
309  TH1D* h1 = p->second;
310 
311  if (h1->GetEntries() != 0) {
312 
313  // start values
314 
315  Double_t ymax = 0.0;
316  Double_t x0 = 0.0; // [ns]
317  Double_t sigma = 250.0; // [ns]
318  Double_t ymin = 0.0;
319 
320  for (int i = 1; i != h1->GetNbinsX(); ++i) {
321 
322  const Double_t x = h1->GetBinCenter (i);
323  const Double_t y = h1->GetBinContent(i);
324 
325  if (y > ymax) {
326  ymax = y;
327  x0 = x;
328  }
329 
330  ymin += y;
331  }
332 
333  ymin /= h1->GetNbinsX();
334 
335  f1.SetParameter(0, ymax);
336  f1.SetParameter(1, x0);
337  if (binWidth_ns < sigma)
338  f1.SetParameter(2, sigma);
339  else
340  f1.FixParameter(2, binWidth_ns/sqrt(12.0));
341  f1.SetParameter(3, ymin);
342 
343  // fit
344 
345  h1->Fit(&f1, option.c_str(), "same", x0 - 5 * sigma, x0 + 5 * sigma);
346 
347  status = (fabs(f1.GetParameter(1)) <= 0.5*getFrameTime() &&
348  f1.GetParameter(0) >= f1.GetParameter(3));
349  }
350 
351  if (h1->GetEntries() != 0) {
352 
353  // look for peaks at offsets equal to a multiple of frame time
354  JWeight bg; // background
355  map<int, JWeight> ps; // signal(s)
356 
357  for (Int_t i = 1, ns = -(numberOfTimeslices + 1); i <= h1->GetNbinsX(); ++i) {
358 
359  const Double_t x = h1->GetBinCenter (i);
360  const Double_t y = h1->GetBinContent(i);
361 
362  while (x > (ns + 1) * getFrameTime() - TMax_ns) { ++ns; }
363 
364  if (fabs(x - ns * getFrameTime()) < TMax_ns)
365  ps[ns].put(y);
366  else
367  bg .put(y);
368  }
369 
370  for (map<int, JWeight>::iterator i = ps.begin(); i != ps.end(); ) {
371 
372  const Double_t noise = bg.getTotal() * i->second.getCount() / bg.getCount();
373  const Double_t P = TMath::PoissonI(i->second.getTotal(), noise);
374 
375  if (P < Pmin && i->second.getTotal() > noise) //Avoid edge effect (second bolean)
376  ++i; // keep
377  else
378  ps.erase(i++); // remove
379  }
380 
381  if (ps.size() != 1) {
382  ERROR("Module " << setw(8) << module->getID() << " number of peaks " << ps.size() << " != 1" << endl);
383  }
384 
385  for (map<int, JWeight>::const_iterator i = ps.begin(); i != ps.end(); ++i) {
386 
387  vector<double> subPeaksPerDoms;
388  const Double_t noise = bg.getTotal() * i->second.getCount() / bg.getCount();
389 
390  if (ps.size() != 1) {
391  WARNING("Peak at " << setw(4) << i->first << " [frame time] S/N = "
392  << FIXED(7,1) << i->second.getTotal() << " / "
393  << FIXED(7,1) << noise << endl);
394  }
395 
396  subPeaksPerDoms.push_back(module->getID());
397  subPeaksPerDoms.push_back(ps.size());
398  subPeaksPerDoms.push_back(i->first);
399  subPeaksPerDoms.push_back(i->second.getTotal());
400  subPeaksPerDoms.push_back(noise);
401  peaksPerDoms.push_back(subPeaksPerDoms);
402  }
403  }
404 
405  //if DOM is OOS
406  if (!status) {
407  //if OOS peak greater than the last and not in OOS DOMs vector
408  if ((f1.GetParameter(0) > most_OOS_peak) && std::find(oosDomsId.begin(), oosDomsId.end(), module->getID()) == oosDomsId.end()){
409 
410  most_OOS_peak = f1.GetParameter(0);
411  most_OOS_ID = module->getID();
412  }
413  }
414  }
415  }
416 
417  if (most_OOS_ID != 0) {
418  oosDomsId.push_back(most_OOS_ID);//exclude the DOM to compute the histograms
419 
420  NOTICE("Module " << most_OOS_ID << " set QEs of all PMTs to zero." << endl);
421 
422  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
423  parameters[JPMTIdentifier(most_OOS_ID, pmt)].QE = 0.0;
424  }
425  manager.clear();
426  peaksPerDoms.clear();
427  }
428  else{
429  OOS = false;
430  }
431  }
432  }while(OOS);
433 
434  //If multipeaks analysis is activated (multiPeaks = 1)
435  //Make ratio : sum_OOS_peaks / sum_all_peaks
436  //If ratio > 0.05, DOM is declared OOS
437 
438  if (multiPeaks == 1) {
439 
440  NOTICE("Multi-peaks analysis activated" << endl);
441 
442  vector<int> workingDOMs;
443 
444  //If DOM is working and not in the list of OOS DOMs
445  for(int i = 0; i != int(peaksPerDoms.size()); i++) {
446  if((std::find(workingDOMs.begin(), workingDOMs.end(), peaksPerDoms[i][0]) == workingDOMs.end()) &&
447  (std::find(oosDomsId.begin(), oosDomsId.end(), peaksPerDoms[i][0]) == oosDomsId.end())) {
448  workingDOMs.push_back(peaksPerDoms[i][0]);
449  }
450  }
451 
452  for(int i = 0; i != int(workingDOMs.size()); i++) {
453 
454  double ratio = 0.;
455  double sumAllPeaks = 0.;
456 
457  for(int j = 0; j != int(peaksPerDoms.size()); j++) {
458  if (peaksPerDoms[j][0] == workingDOMs[i]){
459  if (int(peaksPerDoms[j][2]) != 0) {
460  ratio += peaksPerDoms[j][3] - peaksPerDoms[j][4];
461  }
462  sumAllPeaks += peaksPerDoms[j][3] - peaksPerDoms[j][4];
463  }
464  }
465 
466  ratio /= sumAllPeaks;
467 
468  if (ratio > 0.05){
469  NOTICE("Module " << workingDOMs[i] << " set QEs of all PMTs to zero." << endl);
470  NOTICE("----> Ratio = " << 100*ratio << "%" << endl);
471  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
472  parameters[JPMTIdentifier(workingDOMs[i], pmt)].QE = 0.0;
473  }
474  }
475 
476  }
477  }
478 
479  //To manage non working DOMs
480  NOTICE("Managing non-working DOMs." << endl);
481  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
482  JManager_t::iterator p = manager.find(module->getID());
483  if (p == manager.end()){
484  NOTICE("Module " << module->getID() << " set QEs of all PMTs to zero." << endl);
485  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
486  parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
487  }
488  }
489  }
490 
491  ofstream out(pmtFile.c_str());
492 
493  parameters.comment.add(JMeta(argc, argv));
494 
495  out << parameters << endl;
496 
497  out.close();
498 
499  if (outputFile != "") {
500  manager.Write(outputFile.c_str());
501  }
502 
503  //print small summary file of the peaks
504  if (peaks != ""){
505  ofstream stream(peaks);
506 
507  stream << "#DOM NB_PEAKS TIMESLICE_SHIFT SIGNAL NOISE\n";
508 
509  for(int i = 0; i != int(peaksPerDoms.size()); i++) {
510  for(int j = 0; j != int(peaksPerDoms[i].size()); j++){
511  if (j == 0 || j == 1 || j == 2 || j == 3) stream << int(peaksPerDoms[i][j]) << " " ;
512  else stream << peaksPerDoms[i][j] << " " ;
513  }
514  stream << "\n";
515  }
516  stream.close();
517  }
518 
519 }
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:71
Utility class to parse command line options.
Definition: JParser.hh:1410
void clear()
Clear data.
#define WARNING(A)
Definition: JMessage.hh:63
debug
Definition: JMessage.hh:27
void set(int &mask) const
Set bit in given bit mask.
Definition: JDAQ.hh:77
Auxiliary class to manage set of histograms.
#define STATUS(A)
Definition: JMessage.hh:61
Template const_iterator.
Definition: JDAQEvent.hh:69
Structure to store the ToT mean and standard deviation of the hits produced by a nanobeacon in a sour...
Empty structure for specification of parser element that is initialised (i.e.
Definition: JParser.hh:64
Auxiliary class for a type holder.
Definition: JType.hh:19
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:461
double getTime(const Hit &hit)
Get true time of hit.
string outputFile
int getFrameIndex() const
Get frame index.
JLimit JLimit_t
Type definition of limit.
Definition: JLimit.hh:214
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
double getFrameTime()
Get frame time duration.
Definition: JDAQClock.hh:162
#define NOTICE(A)
Definition: JMessage.hh:62
#define ERROR(A)
Definition: JMessage.hh:64
Data time slice.
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
int debug
debug level
Definition: JSirene.cc:59
double getMaximalTime(const JDetector &detector)
Get maximal time between modules in detector following causality.
#define FATAL(A)
Definition: JMessage.hh:65
Auxiliary data structure for single bit.
Definition: JDAQ.hh:36
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:72
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
JTriggerCounter_t next()
Increment trigger counter.
double hit_type
Definition: JDataFilter.cc:89
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60