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

Auxiliary program to determine PMT parameters from K40 data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <limits>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TMath.h"
#include "JTools/JRange.hh"
#include "JDAQ/JDAQ.hh"
#include "JDAQ/JDAQTimeslice.hh"
#include "JDAQ/JDAQSummaryslice.hh"
#include "JDAQ/JDAQEvaluator.hh"
#include "JTrigger/JSuperFrame1D.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JHitToolkit.hh"
#include "JTrigger/JTriggerParameters.hh"
#include "JTrigger/JTriggerToolkit.hh"
#include "JTrigger/JSummaryRouter.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "JSupport/JTriggerParametersSupportkit.hh"
#include "JLang/JObjectMultiplexer.hh"
#include "JROOT/JROOTClassSelector.hh"
#include "JCalibrate/JCalibrateK40.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

Auxiliary program to determine PMT parameters from K40 data.

Also consult JMonitor.pdf in documentation.

Author
mdejong

Definition in file JCalibrateK40.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 102 of file JCalibrateK40.cc.

103 {
104  using namespace std;
105  using namespace JPP;
106  using namespace KM3NETDAQ;
107 
108 
109  //-----------------------------------------------------------
110  //parameter interface
111  //-----------------------------------------------------------
112 
113  JMultipleFileScanner<JDAQTimesliceTypes_t> inputFile;
114  JLimit_t& numberOfEvents = inputFile.getLimit();
115  string outputFile;
116  string detectorFile;
117  Double_t TMax_ns;
118  JRange<double> rateRange_Hz;
119  JRange<double> totRange_ns;
120  JRange<double> Tail_ns;
121  JRange<int> multiplicity;
122  double deadTime_us;
123  JROOTClassSelector selector;
124  string background;
125  int debug;
126 
127  try {
128 
129  JParser<> zap("Auxiliary program to determine PMT parameters from K40 data.");
130 
131  zap['f'] = make_field(inputFile, "input file.");
132  zap['o'] = make_field(outputFile, "output file.") = "calibrate_k40.root";
133  zap['n'] = make_field(numberOfEvents) = JLimit::max();
134  zap['a'] = make_field(detectorFile, "detector file.");
135  zap['T'] = make_field(TMax_ns, "time window [ns].") = 20.0;
136  zap['V'] = make_field(rateRange_Hz, "PMT rate range [Hz].") = JRange<double>(0.0, 20.0e3);
137  zap['t'] = make_field(totRange_ns, "PMT time-over-threshold range [ns].") = JRange<double>(0.0, ToTmax_ns);
138  zap['b'] = make_field(background, "background estimation method.") = rates_t, count_t, tails_t, rndms_t;
139  zap['B'] = make_field(Tail_ns, "time window used for background estimation.") = JRange<double>(15.0, 20.0);
140  zap['M'] = make_field(multiplicity, "multiplicity range of hits on DOM.") = JRange<int>(2, 31);
141  zap['D'] = make_field(deadTime_us, "L1 dead time (us)") = 0.0;
142  zap['C'] = make_field(selector, "timeslice selector, e.g. JDAQTimesliceL1.") = getROOTClassSelection<JDAQTimesliceTypes_t>();
143  zap['d'] = make_field(debug, "debug flag.") = 1;
144 
145  zap(argc, argv);
146  }
147  catch(const exception &error) {
148  FATAL(error.what() << endl);
149  }
150 
151  //-----------------------------------------------------------
152  // check the input parameters
153  //-----------------------------------------------------------
154 
155  JTriggerParameters parameters;
156 
157  try {
158  parameters = getTriggerParameters(inputFile);
159  }
160  catch(const JException& error) {
161  FATAL("No trigger parameters from input:" << error.what() << endl);
162  }
163 
164  cout.tie(&cerr);
165 
166  if (selector == JDAQTimesliceL2::Class_Name() ||
167  selector == JDAQTimesliceSN::Class_Name()) {
168  FATAL("Option -C <selector> " << selector << " not compatible with calibration method." << endl);
169  }
170 
171  if ((selector == JDAQTimeslice ::Class_Name() && parameters.writeL1.prescale > 0) ||
172  (selector == JDAQTimesliceL1::Class_Name())) {
173 
174  if (parameters.TMaxLocal_ns < TMax_ns) {
175  FATAL("Option -T <TMax_ns> = " << TMax_ns << " is larger than in the trigger " << parameters.TMaxLocal_ns << endl);
176  }
177 
178  if (background == rndms_t ||
179  background == count_t) {
180  FATAL("Option -C <selector> " << selector << " incompatible with option -b <background> " << background << endl);
181  }
182  }
183 
184  if (!multiplicity.is_valid() || multiplicity.getLowerLimit() < 2) {
185  FATAL("Invalid option -M <multiplicity> " << multiplicity << endl);
186  }
187 
188  if (!totRange_ns.is_valid()) {
189  FATAL("Invalid option -t <totRange_ns> " << totRange_ns << endl);
190  }
191 
192  if (background == tails_t) {
193 
194  if (!Tail_ns.is_valid()) {
195  FATAL("Invalid option -B <Tail_ns> " << Tail_ns << endl);
196  }
197 
198  if (Tail_ns.getUpperLimit() > TMax_ns) {
199  FATAL("Invalid option -B <Tail_ns> " << Tail_ns << "; upper limit larger than option -T <TMax_ns> " << TMax_ns << endl);
200  }
201  }
202 
203  //-----------------------------------------------------------
204  // load detector file
205  //-----------------------------------------------------------
206 
207  JDetector detector;
208 
209  try {
210  load(detectorFile, detector);
211  }
212  catch(const JException& error) {
213  FATAL(error);
214  }
215 
216  if (detector.empty()) {
217  FATAL("Empty detector." << endl);
218  }
219 
220  const JModuleRouter router(detector);
221  JSummaryRouter summaryRouter;
222 
223  //-----------------------------------------------------------
224  // determine time offsets for background method rndms_t
225  //-----------------------------------------------------------
226 
227  vector<double> t0; // time offsets for random background evaluation [ns]
228 
229  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
230  t0.push_back(i * 2 * TMax_ns);
231  }
232 
233  const JMatch_t match(TMax_ns); // time window self-coincidences [ns]
234 
235  //-----------------------------------------------------------
236  // correction factor for rate due to specified time-over-threshold range
237  //-----------------------------------------------------------
238 
239  const JPMTAnalogueSignalProcessor cpu;
240 
241  const int NPE = 1;
242  const double RTU = (cpu.getIntegralOfProbability(cpu.getNPE(totRange_ns.getLowerLimit()), cpu.getNPE(totRange_ns.getUpperLimit()), NPE)
243  /
244  cpu.getIntegralOfProbability(cpu.getNPE(0.0), cpu.getNPE(ToTmax_ns), NPE));
245 
246 
247 
248  //-----------------------------------------------------------
249  // initialise histograms
250  //-----------------------------------------------------------
251 
252  JCombinatorics combinatorics(NUMBER_OF_PMTS);
253 
254  vector<JHistogram> zmap(detector.size());
255 
256  TH1D weights_hist(weights_hist_t, NULL, 3, 0.0, 3.0);
257 
258  weights_hist.GetXaxis()->SetBinLabel(1, W1_overall_t); // [s]
259  weights_hist.GetXaxis()->SetBinLabel(2, WS_t); // [ns]
260  weights_hist.GetXaxis()->SetBinLabel(3, WB_t); // [ns]
261 
262  const int nx = combinatorics.getNumberOfPairs();
263  const double xmin = -0.5;
264  const double xmax = nx - 0.5;
265 
266  const double ymin = -floor(TMax_ns) + 0.5;
267  const double ymax = +floor(TMax_ns) - 0.5;
268  const int ny = (int) ((ymax - ymin) / 1.0);
269 
270 
271  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
272 
273  const JModuleAddress& address = router.getAddress(module->getID());
274 
275  NOTICE("Booking histograms for module " << module->getID() << endl);
276 
277  zmap[address.first] = JHistogram(new TH2D(MAKE_CSTRING(module->getID() << _2S), NULL, nx, xmin, xmax, ny, ymin, ymax),
278  new TH1D(MAKE_CSTRING(module->getID() << _1B), NULL, nx, xmin, xmax),
279  new TH1D(MAKE_CSTRING(module->getID() << _1L), NULL, nx, xmin, xmax));
280  }
281 
282 
283  JTreeScanner<JDAQSummaryslice, JDAQEvaluator> summaryFile(inputFile);
284 
285  typedef JSuperFrame2D<JHitR0> JSuperFrame2D_t;
286  typedef JSuperFrame1D<JHitR0> JSuperFrame1D_t;
287 
288  const double deadTime_ns = deadTime_us * 1e3;
289 
290  JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(inputFile, selector);
291 
292  counter_type counter = 0;
293 
294  for ( ; in.hasNext() && counter != inputFile.getLimit(); ++counter) {
295 
296  STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
297 
298  const JDAQTimeslice* timeslice = in.next();
299 
300  if (background == rates_t) {
301 
302  const Long64_t index = summaryFile.find(*timeslice);
303 
304  if (index == -1) {
305  FATAL("Missing summary data at " << timeslice->getFrameIndex() << endl);
306  }
307 
308  summaryRouter.update(summaryFile.getEntry(index));
309 
310  if (timeslice->getFrameIndex() != summaryRouter.getFrameIndex()) {
311 
312  ERROR("Frame indices do not match "
313  << "[counter=" << counter << "]"
314  << "[timeslice=" << timeslice->getFrameIndex() << "]"
315  << "[summaryslice=" << summaryRouter.getFrameIndex() << "]" << endl);
316 
317  if (!in.hasNext())
318  break;
319  else
320  FATAL("ROOT TTrees not aligned; Run number " << timeslice->getRunNumber() << endl);
321  }
322  }
323 
324 
325  for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
326 
327  if (router.hasModule(frame->getModuleID())) {
328 
329  //-----------------------------------------------------------
330  // determine background rates
331  //-----------------------------------------------------------
332 
333  vector<double> rate_Hz(NUMBER_OF_PMTS, 0.0);
334  vector<bool> veto (NUMBER_OF_PMTS, false);
335 
336  JDAQFrameStatus status;
337 
338  if (background == count_t) {
339 
340  status = frame->getDAQFrameStatus();
341 
342  for (JDAQSuperFrame::const_iterator i = frame->begin(); i != frame->end(); ++i) {
343  rate_Hz[i->getPMT()] += RTU * 1e9 / getFrameTime();
344  }
345 
346  } else if (background == rates_t) {
347 
348  if (!summaryRouter.hasSummaryFrame(frame->getModuleID())) {
349  FATAL("Summary frame not found " << timeslice->getFrameIndex() << ":" << frame->getModuleID() << endl);
350  }
351 
352  const JDAQSummaryFrame& summary_frame = summaryRouter.getSummaryFrame(frame->getModuleID());
353 
354  status = summary_frame.getDAQFrameStatus();
355 
356  for (int i = 0; i != NUMBER_OF_PMTS; ++i){
357  rate_Hz[i] = RTU * summary_frame.getRate(i);
358  }
359  }
360 
361  // Set veto if off-shore high-rate veto OR FIFO full.
362  // Hits of PMTs with veto will not be counted in livetime nor coincidences nor background
363 
364  if (status.testHighRateVeto()) {
365  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
366  if (status.testHighRateVeto(i)) {
367  veto[i] = true;
368  }
369  }
370  }
371 
372  if (status.testFIFOStatus()) {
373  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
374  if (status.testFIFOStatus(i)) {
375  veto[i] = true;
376  }
377  }
378  }
379 
380  //-----------------------------------------------------------
381  // sort the PMT pairs by opening angle in the order largest angle-->smallest angle
382  //-----------------------------------------------------------
383 
384  const JModuleAddress& address = router.getAddress(frame->getModuleID());
385  const JModule& module = detector.getModule(address);
386 
387  combinatorics.configure(module.size());
388 
389  combinatorics.sort(JPairwiseComparator(module));
390 
391  TH2D* h2s = zmap[address.first].h2s;
392  TH1D* h1b = zmap[address.first].h1b;
393  TH1D* h1L = zmap[address.first].h1L;
394 
395  //-----------------------------------------------------------
396  // fill the livetime by PMT pairs
397  //-----------------------------------------------------------
398 
399  for (size_t i = 0; i < combinatorics.getNumberOfPairs(); ++i) {
400 
401  const int pmt1 = combinatorics.getPair(i).first;
402  const int pmt2 = combinatorics.getPair(i).second;
403 
404  if (!veto[pmt1] && rateRange_Hz(rate_Hz[pmt1]) &&
405  !veto[pmt2] && rateRange_Hz(rate_Hz[pmt2]) ) {
406 
407  h1L->Fill(i, getFrameTime()*1e-9);
408  }
409  }
410 
411  //-----------------------------------------------------------
412  // process data
413  //-----------------------------------------------------------
414 
415  JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, module);
416 
417  for (JSuperFrame2D_t::iterator i = buffer.begin(); i != buffer.end(); ++i) {
418 
419  const int pmt = i->getPMTAddress();
420 
421  if ( veto[pmt] || !rateRange_Hz(rate_Hz[pmt]) ) {
422  i->reset();
423  } else {
424  i->join(match);
425  }
426 
427  }
428 
429  JSuperFrame1D_t& data = JSuperFrame1D_t::multiplex(buffer);
430 
431  data.pop_back(); // remove end marker
432 
433  DEBUG("Number of hits " << timeslice->getFrameIndex() << ":" << frame->getModuleID() << ' ' << frame->size() << ' ' << data.size() << endl);
434 
435  if (data.empty()) {
436  continue;
437  }
438 
439  // Signal;
440  // Hits from PMTs that do not comply with rate cuts have been exluded above
441 
442  double t1 = -numeric_limits<double>::max();
443 
444  for (vector<JHitR0>::const_iterator p = data.begin(); p != data.end(); ) {
445 
447 
448  while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}
449 
450  const int N = distance(p,q);
451 
452  if (multiplicity(N)) {
453 
454  if (p->getT() > t1 + deadTime_ns) {
455 
456  for (vector<JHitR0>::const_iterator __p = p; __p != q; ++__p) {
457  for (vector<JHitR0>::const_iterator __q = __p; ++__q != q; ) {
458 
459  if (totRange_ns( __p->getToT()) && totRange_ns(__q->getToT())) {
460  h2s->Fill((double) combinatorics.getIndex(__p->getPMT(),__q->getPMT()),
461  JCombinatorics::getSign(__p->getPMT(),__q->getPMT()) * (__q->getT() - __p->getT()));
462  }
463  }
464  }
465  }
466 
467  t1 = p->getT();
468  }
469 
470  p = q;
471  }
472 
473  // Background;
474  // Note that rate cuts and veto have already been accounted for when data buffer was filled
475 
476  if (background == rndms_t) {
477 
478  for (vector<JHitR0>::iterator i = data.begin(); i != data.end(); ++i) {
479  *i = JHitR0(i->getPMT(), JHit(i->getT() + t0[i->getPMT()], i->getToT()));
480  }
481 
482  sort(data.begin(), data.end());
483 
484  for (vector<JHitR0>::const_iterator p = data.begin(); p != data.end(); ) {
485 
487 
488  while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}
489 
490  const int N = distance(p,q);
491 
492  if (multiplicity(N)) {
493 
494  for (vector<JHitR0>::const_iterator __p = p; __p != q; ++__p) {
495  for (vector<JHitR0>::const_iterator __q = __p; ++__q != q; ) {
496 
497  if (totRange_ns(__p->getToT()) && totRange_ns(__q->getToT())) {
498  h1b->Fill((double) combinatorics.getIndex(p->getPMT(),q->getPMT()), 1.0);
499  }
500  }
501  }
502  }
503 
504  p = q;
505  }
506 
507  } else if (background == rates_t ||
508  background == count_t) {
509 
510  Double_t Rs = 0.0; // total rate [Hz]
511 
512  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
513 
514  const Double_t R1 = rate_Hz[i]; // [Hz]
515 
516  if (!veto[i] && rateRange_Hz(R1)) {
517  Rs += R1;
518  }
519  }
520 
521  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
522  for (int j = i; ++j != NUMBER_OF_PMTS; ) {
523 
524  const Double_t R1 = rate_Hz[i]; // [Hz]
525  const Double_t R2 = rate_Hz[j]; // [Hz]
526 
527  if (!veto[i] && rateRange_Hz(R1) &&
528  !veto[j] && rateRange_Hz(R2)) {
529 
530  // expectation values (counts) within TMax_ns * 2 for the optical module and the two PMTs
531  const double ED = (Rs - rate_Hz[i] - rate_Hz[j]) * 2 * TMax_ns * 1e-9;
532  const double E1 = rate_Hz[i] * 2 * TMax_ns * 1e-9;
533  const double E2 = rate_Hz[j] * 2 * TMax_ns * 1e-9;
534 
535  // expected rate of coincidences = P * # of trials
536  Double_t N = coincidenceP(E1, E2, ED,
537  multiplicity.getLowerLimit(),
538  multiplicity.getUpperLimit()) * getFrameTime() / (2*TMax_ns);
539 
540  h1b->Fill((double) combinatorics.getIndex(i,j), N);
541  }
542  }
543  }
544  }
545  }
546  }
547  }
548  STATUS(endl);
549 
550  if (background == tails_t ) {
551 
552  for (vector<JHistogram>::iterator hist = zmap.begin() ; hist != zmap.end() ; ++hist) {
553 
554  TH2D* h2s = hist->h2s;
555  TH1D* h1b = hist->h1b;
556 
557  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
558  for (int j = i; ++j != NUMBER_OF_PMTS; ) {
559 
560  double Y = 0.0; // sum of counts in tail regions
561  double W = 0.0; // square of error of Y
562  int N = 0; // number of bins in tail regions
563  int k = combinatorics.getIndex(i,j); // PMT pair index (x-axis of h2s)
564 
565  for (int l = 1; l <= h2s->GetYaxis()->GetNbins(); ++l) {
566 
567  const Double_t dt = h2s->GetYaxis()->GetBinCenter(l) ;
568 
569  if (Tail_ns(fabs(dt))) {
570  Y += h2s->GetBinContent(k+1,l);
571  W += h2s->GetBinError(k+1,l) * h2s->GetBinError(k+1,l);
572  N++;
573  }
574  }
575 
576  h1b->SetBinContent(k+1, Y / N * (2.0*TMax_ns) / ((ymax - ymin)/ny) );
577  h1b->SetBinError (k+1, sqrt(W/N) * (2.0*TMax_ns) / ((ymax - ymin)/ny) );
578  }
579  }
580  }
581  }
582 
583  //---------------------------------------------
584  // save normalisation constants and store histograms
585  //---------------------------------------------
586 
587  weights_hist.Fill(W1_overall_t, counter*getFrameTime()*1e-9); // [s]
588  weights_hist.Fill(WS_t, (ymax - ymin)/ny); // [ns]
589  weights_hist.Fill(WB_t, 2.0*TMax_ns); // [ns]
590 
591 
592  TFile out(outputFile.c_str(), "recreate");
593 
594  putObject(&out, JMeta(argc, argv));
595 
596  weights_hist.Write() ;
597 
598  for (vector<JHistogram>::iterator i = zmap.begin(); i != zmap.end(); ++i) {
599  i->h2s->Write();
600  i->h1b->Write();
601  i->h1L->Write();
602  }
603 
604  out.Write();
605  out.Close();
606 }
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
Data structure for all trigger parameters.
double getRate(const int tdc, const double factor=1.0) const
Get count rate.
static const char *const WB_t
Named bin for value of TMax_ns in JCalibrateK40.
#define STATUS(A)
Definition: JMessage.hh:61
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:611
Long64_t counter_type
Type definition for counter.
Definition: JCounter.hh:24
JSuperFrame2D< hit_type > JSuperFrame2D_t
Definition: JDataFilter.cc:93
string outputFile
int getRunNumber() const
Get run number.
int getFrameIndex() const
Get frame index.
JLimit JLimit_t
Type definition of limit.
Definition: JLimit.hh:214
static const char *const WS_t
Named bin for time residual bin width.
Hit data structure.
Definition: JDAQHit.hh:36
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
static const char *const _1B
Name extension for 1D background.
double getFrameTime()
Get frame time duration.
Definition: JDAQClock.hh:162
Data storage class for rate measurements of all PMTs in one module.
#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
#define FATAL(A)
Definition: JMessage.hh:65
double coincidenceP(Double_t E1, Double_t E2, Double_t ED, int M_min, int M_max)
Coincidence probability of two PMTs.
bool testHighRateVeto() const
Test high-rate veto status.
static const char *const weights_hist_t
Histogram naming.
static const char *const _1L
Name extension for 1D live time.
static const char *const W1_overall_t
Named bin for duration of the run.
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
const JDAQFrameStatus & getDAQFrameStatus() const
Get DAQ frame status.
JTriggerParameters getTriggerParameters(const JMultipleFileScanner_t &file_list)
Get trigger parameters.
static const char *const _2S
Name extension for 2D counts.
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
JSuperFrame1D< hit_type > JSuperFrame1D_t
Definition: JDataFilter.cc:92
bool testFIFOStatus() const
Test FIFO status.
bool putObject(TDirectory *dir, const T &object)
Write object to ROOT directory.