Jpp  master_rocky
the software that should make you happy
JCalibrateK40.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <limits>
5 #include <array>
6 
7 #include "TROOT.h"
8 #include "TFile.h"
9 #include "TH1D.h"
10 #include "TH2D.h"
11 
14 
15 #include "JDAQ/JDAQTimesliceIO.hh"
17 #include "JDAQ/JDAQEvaluator.hh"
20 #include "JTrigger/JHitR0.hh"
21 #include "JTrigger/JMatchL0.hh"
22 #include "JTrigger/JHitToolkit.hh"
27 #include "JDetector/JDetector.hh"
34 #include "JSupport/JSupport.hh"
35 #include "JSupport/JMeta.hh"
37 #include "JTools/JRange.hh"
39 
41 #include "JROOT/JRootToolkit.hh"
42 #include "JROOT/JRootFileWriter.hh"
45 
46 #include "Jeep/JPrint.hh"
47 #include "Jeep/JParser.hh"
48 #include "Jeep/JMessage.hh"
49 
50 namespace {
51 
52  using JLANG::JObjectID;
55  using JDETECTOR::JModule;
56 
57  // Methods for background estimation
58 
59  static const char* const rates_t = "rates"; //!< singles rates from summary data
60  static const char* const count_t = "counts"; //!< singles rates from L0 data
61  static const char* const tails_t = "tails"; //!< tails of time residual distribution
62  static const char* const rndms_t = "randoms"; //!< random time offsets (only L0 data)
63 
64  /**
65  * Maximal time-over-threshold.
66  */
67  static const double ToTmax_ns = 250;
68 
69  /**
70  * Auxiliary data structure for histogram management.
71  */
72  struct JHistogram :
73  public JObjectID,
74  public JCombinatorics
75  {
76  /**
77  * Default constructor.
78  */
79  JHistogram() :
80  h2s(NULL),
81  h1b(NULL),
82  h1L(NULL)
83  {}
84 
85 
86  /**
87  * Constructor.
88  *
89  * \param module module
90  * \param axis axis
91  */
92  JHistogram(const JModule& module, const JAbstractHistogram<Double_t>& axis) :
93  JObjectID(module.getID())
94  {
95  using namespace JPP;
96 
97  // sort the PMT pairs by opening angle in the order largest angle -> smallest angle
98 
99  this->configure(module.size());
100 
101  this->sort(JPairwiseComparator(module));
102 
103  const Int_t nx = this->getNumberOfPairs();
104  const Double_t xmin = -0.5;
105  const Double_t xmax = nx - 0.5;
106 
107  h2s = new TH2D(MAKE_CSTRING(module.getID() << _2S), NULL, nx, xmin, xmax, axis.getNumberOfBins(), axis.getLowerLimit(), axis.getUpperLimit());
108  h1b = new TH1D(MAKE_CSTRING(module.getID() << _1B), NULL, nx, xmin, xmax);
109  h1L = new TH1D(MAKE_CSTRING(module.getID() << _1L), NULL, nx, xmin, xmax);
110 
111  h2s->Sumw2();
112  h1b->Sumw2();
113  h1L->Sumw2();
114  }
115 
116  /**
117  * Check validity.
118  *
119  * \return true if valid; else false
120  */
121  bool is_valid() const
122  {
123  return (h2s != NULL &&
124  h1b != NULL &&
125  h1L != NULL);
126  }
127 
128  TH2D* h2s; //!< time distribution per PMT pair
129  TH1D* h1b; //!< background rates per PMT pair
130  TH1D* h1L; //!< livetimes per PMT pair
131  };
132 
133 
134  /**
135  * Get coincidence probability of two PMTs within one module due to random background.
136  *
137  * \param E1 expected counts of 1st PMT
138  * \param E2 expected counts of 2nd PMT
139  * \param ED expected counts of module, excluding PMT 1 and PMT 2
140  * \param M_min multiplicity lower limit
141  * \param M_max multiplicity upper limit
142  * \return probability
143  */
144  inline double getP(const double E1, const double E2, const double ED, const int M_min, const int M_max)
145  {
146  // 1st, always calculate the probability for twofold multiplicity
147 
148  double P = E1 * E2 * exp(-(E1 + E2 + ED));
149 
150  // 2nd, calculate the probability of multiplicity at lower limit
151 
152  int MD = 1;
153 
154  for ( ; MD + 2 <= M_min; ++MD) {
155  P *= ED/MD;
156  }
157 
158  // 3rd, add the probabilities of multiplicities up to upper limit
159 
160  double P_sum = P;
161 
162  for ( ; MD + 2 <= M_max; ++MD) {
163  P *= ED/MD;
164  P_sum += P;
165  }
166 
167  return P_sum;
168  }
169 }
170 
171 
172 /**
173  * \file
174  *
175  * Auxiliary program to determine PMT parameters from K40 data.
176  *
177  * By default, the combinatorial background is estimated from the singles rate.\n
178  * In case L0 data are taken (and no summary data are available),
179  * the singles rate can be determined from the amount of L0 data.\n
180  * One can also estimate the background from the tails in the time residual distribution.\n
181  * In that case, option -B can be used to specify the minimal and maximal time offset in ns.\n
182  * For example -B "10 20".\n
183  * Note that the minimal and maximal time should be larger that the width of
184  * the time residual distribution and less than the time window of the L1 coincidence, respectively.\n
185  * The time window is applied symmetrically around a time offset of zero.\n
186  * Finally, if L0 data are available, one can estimate the background using
187  * the same procedure but with a random (read wrong) time calibration.
188  *
189  * The option -t can be used to specify three time-over-threshold values.\n
190  * The time-over-threshold of each hit in the coincidence should be greater or equal the first value and less or equal the last value.\n
191  * Furthermore, the time-over-threshold of one of the two hits should be greater or equals the second value.\n
192  * If you change the (default) values, you may also want to change the corresponding values in JEditPMTParameters.cc.
193  *
194  * The option -M can be used to specify a range of multiplicities.\n
195  * For example -M "2 2" will strictly select two-fold coincidences.
196  * Note that such a selection can bias the result.
197  *
198  * Also consult <a href="https://common.pages.km3net.de/jpp/JMonitor.PDF">documentation</a>.
199  * \author mdejong
200  */
201 int main(int argc, char **argv)
202 {
203  using namespace std;
204  using namespace JPP;
205  using namespace KM3NETDAQ;
206 
207  JMultipleFileScanner_t inputFile;
208  counter_type numberOfEvents;
209  string outputFile;
210  string detectorFile;
211  Double_t TMax_ns;
212  JRange<double> rateRange_Hz;
213  array <double, 3> totRange_ns;
214  JRange<double> Tail_ns;
215  JRange<int> multiplicity;
216  double deadTime_us;
217  JROOTClassSelector selector;
218  string background;
219  JPreprocessor option;
220  int debug;
221 
222  try {
223 
224  JParser<> zap("Auxiliary program to determine PMT parameters from K40 data.");
225 
226  zap['f'] = make_field(inputFile, "input file.");
227  zap['o'] = make_field(outputFile, "output file.") = "calibrate_k40.root";
228  zap['n'] = make_field(numberOfEvents) = JLimit::max();
229  zap['a'] = make_field(detectorFile, "detector file.");
230  zap['T'] = make_field(TMax_ns, "time window [ns].") = 20.0;
231  zap['V'] = make_field(rateRange_Hz, "PMT rate range [Hz].") = JRange<double>(0.0, 20.0e3);
232  zap['t'] = make_field(totRange_ns, "PMT time-over-threshold range [ns].") = { 4.0, 7.0, ToTmax_ns };
233  zap['b'] = make_field(background, "background estimation method.") = rates_t, count_t, tails_t, rndms_t;
234  zap['B'] = make_field(Tail_ns, "time window used for background estimation.") = JRange<double>(15.0, 20.0);
235  zap['M'] = make_field(multiplicity, "multiplicity range of hits on DOM.") = JRange<int>(2, 31);
236  zap['D'] = make_field(deadTime_us, "L1 dead time (us)") = 0.0;
237  zap['C'] = make_field(selector, "timeslice selector, e.g. JDAQTimesliceL1.") = getROOTClassSelection<JDAQTimesliceTypes_t>();
238  zap['O'] = make_field(option, "hit pre-processing option.") = JPreprocessor::getOptions(JPreprocessor::filter_t);
239  zap['d'] = make_field(debug, "debug flag.") = 1;
240 
241  zap(argc, argv);
242  }
243  catch(const exception &error) {
244  FATAL(error.what() << endl);
245  }
246 
247  //-----------------------------------------------------------
248  // check the input parameters
249  //-----------------------------------------------------------
250 
251 
252  if (selector == JDAQTimesliceL2::Class_Name() ||
253  selector == JDAQTimesliceSN::Class_Name()) {
254  FATAL("Option -C <selector> " << selector << " not compatible with calibration method." << endl);
255  }
256 
257  if (selector == JDAQTimeslice ::Class_Name() ||
258  selector == JDAQTimesliceL1::Class_Name()) {
259 
260  JTriggerParameters parameters;
261 
262  try {
263  parameters = getTriggerParameters(inputFile);
264  }
265  catch(const JException& error) {
266  FATAL("No trigger parameters from input:" << error.what() << endl);
267  }
268 
269  if ((selector == JDAQTimeslice ::Class_Name() && parameters.writeL1.prescale > 0) ||
270  (selector == JDAQTimesliceL1::Class_Name())) {
271 
272  if (parameters.TMaxLocal_ns < TMax_ns) {
273  FATAL("Option -T <TMax_ns> = " << TMax_ns << " is larger than in the trigger " << parameters.TMaxLocal_ns << endl);
274  }
275 
276  if (background == rndms_t ||
277  background == count_t) {
278  FATAL("Option -C <selector> " << selector << " incompatible with option -b <background> " << background << endl);
279  }
280  }
281  }
282 
283  if (!multiplicity.is_valid() || multiplicity.getLowerLimit() < 2 || multiplicity.getUpperLimit() > NUMBER_OF_PMTS) {
284  FATAL("Invalid option -M <multiplicity> " << multiplicity << endl);
285  }
286 
287  if (totRange_ns[0] > totRange_ns[1] ||
288  totRange_ns[1] > totRange_ns[2]) {
289  FATAL("Invalid option -t <totRange_ns> " << JEEPZ() << totRange_ns << endl);
290  }
291 
292  if (background == tails_t) {
293 
294  if (!Tail_ns.is_valid()) {
295  FATAL("Invalid option -B <Tail_ns> " << Tail_ns << endl);
296  }
297 
298  if (Tail_ns.getUpperLimit() > TMax_ns) {
299  FATAL("Invalid option -B <Tail_ns> " << Tail_ns << "; upper limit larger than option -T <TMax_ns> " << TMax_ns << endl);
300  }
301  }
302 
303  //-----------------------------------------------------------
304  // load detector file
305  //-----------------------------------------------------------
306 
308 
309  try {
310  load(detectorFile, detector);
311  }
312  catch(const JException& error) {
313  FATAL(error);
314  }
315 
316  if (detector.empty()) {
317  FATAL("Empty detector." << endl);
318  }
319 
320  const JModuleRouter router(detector);
321 
322  //-----------------------------------------------------------
323  // determine time offsets for background method rndms_t
324  //-----------------------------------------------------------
325 
326  vector<double> t0; // time offsets for random background evaluation [ns]
327 
328  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
329  t0.push_back(i * 2 * TMax_ns);
330  }
331 
332  //-----------------------------------------------------------
333  // correction factor for rate due to specified time-over-threshold range
334  //-----------------------------------------------------------
335 
336  const JPMTAnalogueSignalProcessor cpu;
337 
338  const int NPE = 1;
339  const double RTU = (cpu.getIntegralOfChargeProbability(cpu.getNPE(totRange_ns[0]), cpu.getNPE(totRange_ns[2]), NPE)
340  /
341  cpu.getIntegralOfChargeProbability(cpu.getNPE(0.0), cpu.getNPE(ToTmax_ns), NPE));
342 
343 
344  //-----------------------------------------------------------
345  // initialise histograms
346  //-----------------------------------------------------------
347 
348  vector<JHistogram> zmap(detector.size());
349 
350  const Double_t ymin = -floor(TMax_ns) + 0.5;
351  const Double_t ymax = +floor(TMax_ns) - 0.5;
352  const Int_t ny = (Int_t) ((ymax - ymin) / 1.0);
353 
354  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
355 
356  const JModuleAddress& address = router.getAddress(module->getID());
357 
358  if (!module->empty()) {
359 
360  NOTICE("Booking histograms for module " << module->getID() << endl);
361 
362  zmap[address.first] = JHistogram(*module, JAbstractHistogram<Double_t>(ny, ymin, ymax));
363  }
364  }
365 
366 
367  typedef JHitR0 hit_type;
368  typedef JSuperFrame2D<hit_type> JSuperFrame2D_t;
369  typedef JSuperFrame1D<hit_type> JSuperFrame1D_t;
370 
371  const JMatchL0<hit_type> match(TMax_ns); // time window self-coincidences [ns]
372 
373  const double deadTime_ns = deadTime_us * 1e3;
374 
375 
376  TFile out(outputFile.c_str(), "recreate");
377 
378  putObject(out, JMeta(argc, argv));
379 
380  counter_type counter = 0;
381 
382  for (JMultipleFileScanner_t::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
383 
384  STATUS("Processing: " << *i << endl);
385 
386  JSummaryFileRouter summary(*i);
387 
390 
391  for (JDAQHeader header; in.hasNext() && counter != numberOfEvents; ++counter) {
392 
393  STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
394 
395  const JDAQTimeslice* timeslice = in.next();
396 
397  if (header.getDetectorID() != timeslice->getDetectorID() ||
398  header.getRunNumber () != timeslice->getRunNumber ()) {
399 
400  header = timeslice->getDAQHeader();
401 
402  putObject(out, header);
403  }
404 
405  if (background == rates_t) {
406 
407  summary.update(timeslice->getDAQHeader());
408 
409  if (timeslice->getFrameIndex() != summary.getSummaryslice()->getFrameIndex()) {
410 
411  ERROR("Frame indices do not match at "
412  << "[counter = " << counter << "]"
413  << "[timeslice = " << timeslice->getFrameIndex() << "]"
414  << "[summaryslice = " << summary.getSummaryslice()->getFrameIndex() << "]" << endl);
415 
416  continue;
417  }
418  }
419 
420  for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
421 
422  if (router.hasModule(frame->getModuleID())) {
423 
424  const JModule& module = router.getModule(frame->getModuleID());
425 
426  //-----------------------------------------------------------
427  // determine background rates and veto per PMT
428  //-----------------------------------------------------------
429 
430  vector<double> rate_Hz(NUMBER_OF_PMTS, 0.0);
431  vector<bool> veto (NUMBER_OF_PMTS, false);
432 
433  if (background == count_t) {
434 
435  for (JDAQSuperFrame::const_iterator i = frame->begin(); i != frame->end(); ++i) {
436  rate_Hz[i->getPMT()] += RTU * 1e9 / getFrameTime();
437  }
438 
439  } else if (background == tails_t) {
440 
441  // see below
442 
443  } else if (background == rndms_t) {
444 
445  // see below
446 
447  } else if (background == rates_t) {
448 
449  if (summary.hasSummaryFrame(frame->getModuleID())) {
450 
451  const JDAQSummaryFrame& sum = summary.getSummaryFrame(frame->getModuleID());
452 
453  for (int i = 0; i != NUMBER_OF_PMTS; ++i){
454  rate_Hz[i] = sum.getRate (i) *RTU;
455  veto [i] = sum.getValue(i) == 0;
456  }
457 
458  } else {
459 
460  FATAL("Summary frame of module " << frame->getModuleID() << " not found at frame index " << timeslice->getFrameIndex() << endl);
461  }
462  }
463 
464  const JDAQFrameStatus status = frame->getDAQFrameStatus();
465 
466  // Set veto according DAQ or PMT status and rate;
467  // Hits of PMTs with veto will not be counted in livetime nor coincidences nor background
468 
469  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
470  if (!getDAQStatus(status, module, i) ||
471  !getPMTStatus(status, module, i) ||
472  !rateRange_Hz(rate_Hz[i])) {
473  veto[i] = true;
474  }
475  }
476 
477  const JHistogram& histogram = zmap[router.getAddress(frame->getModuleID()).first];
478  const JCombinatorics& combinatorics = histogram.getCombinatorics();
479 
480  TH2D* h2s = histogram.h2s;
481  TH1D* h1b = histogram.h1b;
482  TH1D* h1L = histogram.h1L;
483 
484  //-----------------------------------------------------------
485  // process data
486  //-----------------------------------------------------------
487 
488  JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, module);
489 
490  buffer.preprocess(option, match);
491 
492  for (JSuperFrame2D_t::iterator i = buffer.begin(); i != buffer.end(); ++i) {
493  if (veto[i->getPMTAddress()]) {
494  i->reset();
495  }
496  }
497 
498  JSuperFrame1D_t& data = JSuperFrame1D_t::multiplex(buffer);
499 
500  DEBUG("Number of hits " << timeslice->getFrameIndex() << ":" << frame->getModuleID() << ' ' << frame->size() << ' ' << data.size() << endl);
501 
502  // Signal;
503  // Hits from PMTs that do not comply with rate cuts have been exluded above
504 
505  size_t numberOfSignalEvents = 0;
506 
507  double t1 = numeric_limits<double>::lowest();
508 
509  for (vector<hit_type>::const_iterator p = data.begin(); p != data.end(); ) {
510 
512 
513  while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}
514 
515  const int N = distance(p,q);
516 
517  if (multiplicity(N)) {
518 
519  numberOfSignalEvents += 1;
520 
521  if (p->getT() > t1 + deadTime_ns) {
522 
523  for (vector<hit_type>::const_iterator __p = p; __p != q; ++__p) {
524  for (vector<hit_type>::const_iterator __q = __p; ++__q != q; ) {
525 
526  if ((__p->getToT() >= totRange_ns[0]) &&
527  (__q->getToT() >= totRange_ns[0]) &&
528  (__p->getToT() >= totRange_ns[1] ||
529  __q->getToT() >= totRange_ns[1]) &&
530  (__p->getToT() <= totRange_ns[2]) &&
531  (__q->getToT() <= totRange_ns[2])) {
532 
533  h2s->Fill((double) combinatorics.getIndex(__p->getPMT(),__q->getPMT()),
534  JCombinatorics::getSign(__p->getPMT(),__q->getPMT()) * (__q->getT() - __p->getT()));
535  }
536  }
537  }
538 
539  t1 = p->getT();
540  }
541  }
542 
543  p = q;
544  }
545 
546  // Background;
547  // Note that rate cuts and veto have already been accounted for when data buffer was filled
548 
549  if (background == rndms_t) {
550 
551  for (vector<hit_type>::iterator i = data.begin(); i != data.end(); ++i) {
552  *i = hit_type(i->getPMT(), JHit(i->getT() + t0[i->getPMT()], i->getToT()));
553  }
554 
555  sort(data.begin(), data.end());
556 
557  double t1 = numeric_limits<double>::lowest();
558 
559  for (vector<hit_type>::const_iterator p = data.begin(); p != data.end(); ) {
560 
562 
563  while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}
564 
565  const int N = distance(p,q);
566 
567  if (multiplicity(N)) {
568 
569  if (p->getT() > t1 + deadTime_ns) {
570 
571  for (vector<hit_type>::const_iterator __p = p; __p != q; ++__p) {
572  for (vector<hit_type>::const_iterator __q = __p; ++__q != q; ) {
573 
574  if ((__p->getToT() >= totRange_ns[0]) &&
575  (__q->getToT() >= totRange_ns[0]) &&
576  (__p->getToT() >= totRange_ns[1] ||
577  __q->getToT() >= totRange_ns[1]) &&
578  (__p->getToT() <= totRange_ns[2]) &&
579  (__q->getToT() <= totRange_ns[2])) {
580 
581  h1b->Fill((double) combinatorics.getIndex(p->getPMT(),q->getPMT()), 1.0);
582  }
583  }
584  }
585 
586  t1 = p->getT();
587  }
588  }
589 
590  p = q;
591  }
592 
593  } else if (background == rates_t ||
594  background == count_t) {
595 
596  double Rs = 0.0; // total rate [Hz]
597 
598  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
599  if (!veto[i]) {
600  Rs += rate_Hz[i]; // [Hz]
601  }
602  }
603 
604  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
605  for (int j = i; ++j != NUMBER_OF_PMTS; ) {
606 
607  if (!veto[i] && !veto[j]) {
608 
609  const double R1 = rate_Hz[i]; // [Hz]
610  const double R2 = rate_Hz[j]; // [Hz]
611 
612  // evaluate expected counts within a time window of 2 * TMax_ns for the two PMTs and the other PMTs inside the optical module;
613  // expected rate due to random coincidences then is product of the probability for the specific observation and the number of trials
614  // the observation is one hit in PMT 1, one hit in PMT 2 and a number of hits in the other PMTs inside the same module
615  // corresponding to the given multiplicity range
616 
617  const double N = getP((R1) * 2 * TMax_ns * 1e-9,
618  (R2) * 2 * TMax_ns * 1e-9,
619  (Rs - R1 - R2) * 2 * TMax_ns * 1e-9,
620  multiplicity.getLowerLimit(),
621  multiplicity.getUpperLimit()) * getFrameTime() / (2*TMax_ns);
622 
623  h1b->Fill((double) combinatorics.getIndex(i,j), N);
624  }
625  }
626  }
627  }
628 
629  //-----------------------------------------------------------
630  // fill the livetime by PMT pairs
631  //-----------------------------------------------------------
632 
633  const double livetime_s = getFrameTime()*1e-9 * exp(-deadTime_ns * numberOfSignalEvents/getFrameTime());
634 
635  for (size_t i = 0; i < combinatorics.getNumberOfPairs(); ++i) {
636 
637  const int pmt1 = combinatorics.getPair(i).first;
638  const int pmt2 = combinatorics.getPair(i).second;
639 
640  if (!veto[pmt1] && !veto[pmt2]) {
641  h1L->Fill(i, livetime_s);
642  }
643  }
644  }
645  }
646  }
647  }
648  STATUS(endl);
649 
650  if (background == tails_t ) {
651 
652  const double R = (2.0*TMax_ns) / ((ymax - ymin)/ny);
653 
654  for (vector<JHistogram>::iterator hr = zmap.begin() ; hr != zmap.end() ; ++hr) {
655 
656  if (hr->is_valid()) {
657 
658  TH2D* h2s = hr->h2s;
659  TH1D* h1b = hr->h1b;
660 
661  for (int ix = 1; ix <= h1b->GetXaxis()->GetNbins(); ++ix) {
662 
663  double Y = 0.0; // sum of counts in tail regions
664  double W = 0.0; // square of error of Y
665  int N = 0; // number of bins in tail regions
666 
667  for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
668 
669  const double y = h2s->GetYaxis()->GetBinCenter(iy) ;
670 
671  if (Tail_ns(fabs(y))) {
672  Y += h2s->GetBinContent(ix,iy);
673  W += h2s->GetBinError (ix,iy) * h2s->GetBinError(ix,iy);
674  N += 1;
675  }
676  }
677 
678  h1b->SetBinContent(ix, Y/N * R);
679  h1b->SetBinError (ix, sqrt(W/N) * R);
680  }
681  }
682  }
683  }
684 
685  //---------------------------------------------
686  // save normalisation constants and store histograms
687  //---------------------------------------------
688 
689  TH1D weights_hist(weights_hist_t, NULL, 3, 0.0, 3.0);
690 
691  weights_hist.GetXaxis()->SetBinLabel(1, W1_overall_t); // [s]
692  weights_hist.GetXaxis()->SetBinLabel(2, WS_t); // [ns]
693  weights_hist.GetXaxis()->SetBinLabel(3, WB_t); // [ns]
694 
695  weights_hist.Fill(W1_overall_t, counter*getFrameTime()*1e-9); // [s]
696  weights_hist.Fill(WS_t, (ymax - ymin)/ny); // [ns]
697  weights_hist.Fill(WB_t, 2.0*TMax_ns); // [ns]
698 
699  out << weights_hist;
700 
701  for (vector<JHistogram>::iterator i = zmap.begin(); i != zmap.end(); ++i) {
702  if (i->is_valid()) {
703  out << *(i->h2s);
704  out << *(i->h1b);
705  out << *(i->h1L);
706  }
707  }
708 
709  out.Write();
710  out.Close();
711 }
int main(int argc, char **argv)
string outputFile
KM3NeT DAQ constants, bit handling, etc.
Data structure for detector geometry and calibration.
Basic data structure for L0 hit.
Tools for handling different hit types.
Match operator for consecutive hits.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
#define STATUS(A)
Definition: JMessage.hh:63
#define ERROR(A)
Definition: JMessage.hh:66
#define NOTICE(A)
Definition: JMessage.hh:64
#define FATAL(A)
Definition: JMessage.hh:67
int debug
debug level
Definition: JSirene.cc:69
ROOT I/O of application specific meta data.
Direct access to module in detector data structure.
Scanning of objects from multiple files according a format that follows from the extension of each fi...
PMT analogue signal processor.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2142
Auxiliaries for pre-processing of hits.
I/O formatting auxiliaries.
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:72
Auxiliary class to define a range between two values.
Scanning of objects from a single file according a format that follows from the extension of each fil...
ROOT TTree parameter settings of various packages.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Detector data structure.
Definition: JDetector.hh:96
Address of module in detector data structure.
int first
index of module in detector data structure
Router for direct addressing of module data in detector data structure.
const JModuleAddress & getAddress(const JObjectID &id) const
Get address of module.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
bool hasModule(const JObjectID &id) const
Has module.
Data structure for a composite optical module.
Definition: JModule.hh:75
General exception.
Definition: JException.hh:24
virtual const char * what() const override
Get error message.
Definition: JException.hh:64
Auxiliary class for object identification.
Definition: JObjectID.hh:25
int getID() const
Get identifier.
Definition: JObjectID.hh:50
Auxiliary class for multiplexing object iterators.
virtual bool hasNext() override
Check availability of next element.
virtual const pointer_type & next() override
Get next element.
Utility class to parse command line options.
Definition: JParser.hh:1698
Object reading from a list of files.
File router for fast addressing of summary data.
void update(const JDAQHeader &header)
Update router.
bool hasSummaryFrame(const JDAQModuleIdentifier &module) const
Has summary frame.
const JDAQSummaryFrame & getSummaryFrame() const
Get default summary frame.
const JDAQSummaryslice * getSummaryslice() const
Get summary slice.
Auxiliary class to convert pair of indices to unique index and back.
int getIndex(const int first, const int second) const
Get index of pair of indices.
size_t getNumberOfPairs() const
Get number of pairs.
const pair_type & getPair(const int index) const
Get pair of indices for given index.
Template definition of histogram object interface.
Definition: JHistogram.hh:28
bool is_valid() const
Check validity of range.
Definition: JRange.hh:311
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
Reduced data structure for L0 hit.
Definition: JHitR0.hh:27
L0 match criterion.
Definition: JMatchL0.hh:29
1-dimensional frame with time calibrated data from one optical module.
2-dimensional frame with time calibrated data from one optical module.
int getDetectorID() const
Get detector identifier.
int getRunNumber() const
Get run number.
int getFrameIndex() const
Get frame index.
const JDAQFrameStatus & getDAQFrameStatus() const
Get DAQ frame status.
const JDAQHeader & getDAQHeader() const
Get DAQ header.
Definition: JDAQHeader.hh:49
Hit data structure.
Definition: JDAQHit.hh:35
Data storage class for rate measurements of all PMTs in one module.
double getRate(const int tdc, const double factor=1.0) const
Get count rate.
JRate_t getValue(const int tdc) const
Get value.
#define R1(x)
const double xmax
Definition: JQuadrature.cc:24
const double xmin
Definition: JQuadrature.cc:23
static const char *const _2S
Name extension for 2D counts.
static const char *const WB_t
Named bin for value of TMax_ns in JCalibrateK40.
static const char *const weights_hist_t
Histogram naming.
static const char *const _1B
Name extension for 1D background.
static const char *const WS_t
Named bin for time residual bin width.
static const char *const W1_overall_t
Named bin for duration of the run.
static const char *const _1L
Name extension for 1D live time.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
double getP(const double expval, bool hit)
Get Poisson probability to observe a hit or not for given expectation value for the number of hits.
Definition: JFitToolkit.hh:41
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
Long64_t counter_type
Type definition for counter.
bool is_valid(const json &js)
Check validity of JSon data.
JTriggerParameters getTriggerParameters(const JMultipleFileScanner_t &file_list)
Get trigger parameters.
void configure(const T &value, const JAbstractCollection< JAbscissa_t > &bounds, JBool< false > option)
Configuration of value.
int j
Definition: JPolint.hh:792
bool getDAQStatus(const JDAQFrameStatus &frame, const JStatus &status)
Test status of DAQ.
bool getPMTStatus(const JStatus &status)
Test status of PMT.
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
double getFrameTime()
Get frame time duration.
Definition: JDAQClock.hh:162
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
Definition: JSTDTypes.hh:14
Detector file.
Definition: JHead.hh:227
Acoustics hit.
Auxiliary class to sort pairs of PMT addresses within optical module.
virtual double getNPE(const double tot_ns) const override
Get number of photo-electrons.
double getIntegralOfChargeProbability(const double xmin, const double xmax, const int NPE) const
Get integral of probability.
Auxiliary data structure for streaming of STL containers.
Auxiliary class to select ROOT class based on class name.
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:72
Auxiliary base class for list of file names.
Simple data structure for histogram binning.
int getNumberOfBins() const
Get number of bins.
Auxiliary class for specifying the way of pre-processing of hits.