JCalibrateK40.cc File Reference

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

#include <string>
#include <iostream>
#include <iomanip>
#include <limits>
#include <array>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "km3net-dataformat/online/JDAQ.hh"
#include "km3net-dataformat/online/JDAQHeader.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JDAQ/JDAQEvaluator.hh"
#include "JTrigger/JSuperFrame1D.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JMatchL0.hh"
#include "JTrigger/JHitToolkit.hh"
#include "JTrigger/JTriggerParameters.hh"
#include "JTrigger/JTriggerToolkit.hh"
#include "JTrigger/JPreprocessor.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"
#include "JSupport/JSingleFileScanner.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSummaryFileRouter.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "JSupport/JTriggerParametersSupportkit.hh"
#include "JTools/JRange.hh"
#include "JTools/JAbstractHistogram.hh"
#include "JROOT/JROOTClassSelector.hh"
#include "JROOT/JRootToolkit.hh"
#include "JROOT/JRootFileWriter.hh"
#include "JLang/JObjectMultiplexer.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.

By default, the combinatorial background is estimated from the singles rate.
In case L0 data are taken (and no summary data are available), the singles rate can be determined from the amount of L0 data.
One can also estimate the background from the tails in the time residual distribution.
In that case, option -B can be used to specify the minimal and maximal time offset in ns.
For example -B "10 20".
Note that the minimal and maximal time should be larger that the width of the time residual distribution and less than the time window of the L1 coincidence, respectively.
The time window is applied symmetrically around a time offset of zero.
Finally, if L0 data are available, one can estimate the background using the same procedure but with a random (read wrong) time calibration.

The option -t can be used to specify three time-over-threshold values.
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.
Furthermore, the time-over-threshold of one of the two hits should be greater or equals the second value.
If you change the (default) values, you may also want to change the corresponding values in JEditPMTParameters.cc.

The option -M can be used to specify a range of multiplicities.
For example -M "2 2" will strictly select two-fold coincidences. Note that such a selection can bias the result.

Also consult documentation.

Author

mdejong

Definition in file JCalibrateK40.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 201 of file JCalibrateK40.cc.

{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  JMultipleFileScanner_t inputFile;
  counter_type       numberOfEvents;
  string             outputFile;
  string             detectorFile;
  Double_t           TMax_ns;
  JRange<double>     rateRange_Hz;
  array <double, 3>  totRange_ns;
  JRange<double>     Tail_ns;
  JRange<int>        multiplicity;
  double             deadTime_us;
  JROOTClassSelector selector;
  string             background;
  JPreprocessor      option;
  int                debug;
    
  try { 

    JParser<> zap("Auxiliary program to determine PMT parameters from K40 data.");
    
    zap['f'] = make_field(inputFile,      "input file.");
    zap['o'] = make_field(outputFile,     "output file.")                                = "calibrate_k40.root";
    zap['n'] = make_field(numberOfEvents)                                                = JLimit::max();
    zap['a'] = make_field(detectorFile,   "detector file.");
    zap['T'] = make_field(TMax_ns,        "time window [ns].")                           = 20.0;
    zap['V'] = make_field(rateRange_Hz,   "PMT rate range [Hz].")                        = JRange<double>(0.0, 20.0e3);
    zap['t'] = make_field(totRange_ns,    "PMT time-over-threshold range [ns].")         = { 4.0, 7.0, ToTmax_ns };
    zap['b'] = make_field(background,     "background estimation method.")               = rates_t, count_t, tails_t, rndms_t;
    zap['B'] = make_field(Tail_ns,        "time window used for background estimation.") = JRange<double>(15.0, 20.0);
    zap['M'] = make_field(multiplicity,   "multiplicity range of hits on DOM.")          = JRange<int>(2, 31);
    zap['D'] = make_field(deadTime_us,    "L1 dead time (us)")                           = 0.0;
    zap['C'] = make_field(selector,       "timeslice selector, e.g. JDAQTimesliceL1.")   = getROOTClassSelection<JDAQTimesliceTypes_t>();
    zap['O'] = make_field(option,         "hit pre-processing option.")                  = JPreprocessor::getOptions(JPreprocessor::filter_t);
    zap['d'] = make_field(debug,          "debug flag.")                                 = 1;
     
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  //-----------------------------------------------------------
  // check the input parameters
  //-----------------------------------------------------------


  if (selector == JDAQTimesliceL2::Class_Name() ||
      selector == JDAQTimesliceSN::Class_Name()) {
    FATAL("Option -C <selector> " << selector << " not compatible with calibration method." << endl);
  }

  if (selector == JDAQTimeslice  ::Class_Name() ||
      selector == JDAQTimesliceL1::Class_Name()) {
  
    JTriggerParameters parameters;

    try {
      parameters = getTriggerParameters(inputFile);
    }
    catch(const JException& error) {
      FATAL("No trigger parameters from input:" << error.what() << endl);
    }

    if ((selector == JDAQTimeslice  ::Class_Name() && parameters.writeL1.prescale > 0) ||
        (selector == JDAQTimesliceL1::Class_Name())) {

      if (parameters.TMaxLocal_ns < TMax_ns) {
        FATAL("Option -T <TMax_ns> = " << TMax_ns << " is larger than in the trigger " << parameters.TMaxLocal_ns << endl);
      }

      if (background == rndms_t ||
          background == count_t) {
        FATAL("Option -C <selector> " << selector << " incompatible with option -b <background> " << background << endl);
      }
    }
  }

  if (!multiplicity.is_valid() || multiplicity.getLowerLimit() < 2 || multiplicity.getUpperLimit() > NUMBER_OF_PMTS) {
    FATAL("Invalid option -M <multiplicity> " << multiplicity << endl);
  }

  if (totRange_ns[0] > totRange_ns[1] ||
      totRange_ns[1] > totRange_ns[2]) {
    FATAL("Invalid option -t <totRange_ns> " << JEEPZ() << totRange_ns << endl);
  }

  if (background == tails_t) {

    if (!Tail_ns.is_valid()) {
      FATAL("Invalid option -B <Tail_ns> " << Tail_ns << endl);
    }

    if (Tail_ns.getUpperLimit() > TMax_ns) {
      FATAL("Invalid option -B <Tail_ns> " << Tail_ns << "; upper limit larger than option -T <TMax_ns> " << TMax_ns << endl);
    }
  }

  //-----------------------------------------------------------
  // load detector file
  //-----------------------------------------------------------

  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  if (detector.empty()) {
    FATAL("Empty detector." << endl);
  }

  const JModuleRouter router(detector);

  //-----------------------------------------------------------
  // determine time offsets for background method rndms_t
  //-----------------------------------------------------------

  vector<double> t0;  // time offsets for random background evaluation [ns]

  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
    t0.push_back(i * 2 * TMax_ns);
  }

  //-----------------------------------------------------------
  // correction factor for rate due to specified time-over-threshold range
  //-----------------------------------------------------------

  const JPMTAnalogueSignalProcessor cpu;
  
  const int    NPE = 1;
  const double RTU = (cpu.getIntegralOfChargeProbability(cpu.getNPE(totRange_ns[0]), cpu.getNPE(totRange_ns[2]), NPE)
                      /
                      cpu.getIntegralOfChargeProbability(cpu.getNPE(0.0),            cpu.getNPE(ToTmax_ns),      NPE));
  

  //-----------------------------------------------------------
  // initialise histograms
  //-----------------------------------------------------------

  vector<JHistogram> zmap(detector.size());  

  const Double_t ymin = -floor(TMax_ns) + 0.5;
  const Double_t ymax = +floor(TMax_ns) - 0.5;
  const Int_t    ny   = (Int_t) ((ymax - ymin) / 1.0);

  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
  
    const JModuleAddress& address = router.getAddress(module->getID());

    if (!module->empty()) {

      NOTICE("Booking histograms for module " << module->getID() << endl);

      zmap[address.first] = JHistogram(*module, JAbstractHistogram<Double_t>(ny, ymin, ymax));
    }
  }


  typedef JHitR0        hit_type;
  typedef JSuperFrame2D<hit_type>   JSuperFrame2D_t;
  typedef JSuperFrame1D<hit_type>   JSuperFrame1D_t;

  const JMatchL0<hit_type> match(TMax_ns);     // time window self-coincidences [ns]

  const double deadTime_ns  =  deadTime_us * 1e3;


  TFile out(outputFile.c_str(), "recreate");

  putObject(out, JMeta(argc, argv));

  counter_type counter = 0;

  for (JMultipleFileScanner_t::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {

    STATUS("Processing: " << *i << endl);

    JSummaryFileRouter summary(*i);

    JSingleFileScanner<JDAQTimesliceTypes_t>                reader(*i);
    JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(reader, selector);

    for (JDAQHeader header; in.hasNext() && counter != numberOfEvents; ++counter) {

      STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);

      const JDAQTimeslice* timeslice = in.next();

      if (header.getDetectorID() != timeslice->getDetectorID() ||
          header.getRunNumber () != timeslice->getRunNumber ()) {

        header = timeslice->getDAQHeader();

        putObject(out, header);
      }

      if (background == rates_t) {

        summary.update(timeslice->getDAQHeader());

        if (timeslice->getFrameIndex() != summary.getSummaryslice()->getFrameIndex()) {

          ERROR("Frame indices do not match at " 
                << "[counter = "      << counter                                    << "]" 
                << "[timeslice = "    << timeslice->getFrameIndex()                 << "]"
                << "[summaryslice = " << summary.getSummaryslice()->getFrameIndex() << "]" << endl);

          continue;
        }
      }

      for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {

        if (router.hasModule(frame->getModuleID())) {

          const JModule& module = router.getModule(frame->getModuleID());

          //-----------------------------------------------------------
          // determine background rates and veto per PMT
          //-----------------------------------------------------------

          vector<double> rate_Hz(NUMBER_OF_PMTS, 0.0);
          vector<bool>   veto   (NUMBER_OF_PMTS, false);

          if        (background == count_t) {

            for (JDAQSuperFrame::const_iterator i = frame->begin(); i != frame->end(); ++i) {
              rate_Hz[i->getPMT()] += RTU * 1e9 / getFrameTime();
            }

          } else if (background == tails_t) {

            // see below

          } else if (background == rndms_t) {

            // see below

          } else if (background == rates_t) {

            if (summary.hasSummaryFrame(frame->getModuleID())) {

              const JDAQSummaryFrame& sum = summary.getSummaryFrame(frame->getModuleID());

              for (int i = 0; i != NUMBER_OF_PMTS; ++i){
                rate_Hz[i] = sum.getRate (i) *RTU;
                veto   [i] = sum.getValue(i) == 0;
              }

            } else {

              FATAL("Summary frame of module " << frame->getModuleID() << " not found at frame index " << timeslice->getFrameIndex() << endl);
            }
          }

          const JDAQFrameStatus status = frame->getDAQFrameStatus();
        
          // Set veto according DAQ or PMT status and rate;
          // Hits of PMTs with veto will not be counted in livetime nor coincidences nor background

          for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
            if (!getDAQStatus(status, module, i) ||
                !getPMTStatus(status, module, i) ||
                !rateRange_Hz(rate_Hz[i])) {
              veto[i] = true;
            }
          }

          const JHistogram&     histogram     = zmap[router.getAddress(frame->getModuleID()).first];
          const JCombinatorics& combinatorics = histogram.getCombinatorics();

          TH2D* h2s = histogram.h2s;
          TH1D* h1b = histogram.h1b;
          TH1D* h1L = histogram.h1L;

          //-----------------------------------------------------------
          // process data
          //-----------------------------------------------------------

          JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, module);

          buffer.preprocess(option, match);

          for (JSuperFrame2D_t::iterator i = buffer.begin(); i != buffer.end(); ++i) {
            if (veto[i->getPMTAddress()]) {
              i->reset();
            }
          }

          JSuperFrame1D_t& data = JSuperFrame1D_t::multiplex(buffer);

          DEBUG("Number of hits " << timeslice->getFrameIndex() << ":" << frame->getModuleID() << ' ' << frame->size() << ' ' << data.size() << endl);

          // Signal;
          // Hits from PMTs that do not comply with rate cuts have been exluded above

          size_t numberOfSignalEvents = 0;

          double t1 = numeric_limits<double>::lowest();
          
          for (vector<hit_type>::const_iterator p = data.begin(); p != data.end(); ) {
            
            vector<hit_type>::const_iterator q = p;
            
            while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}

            const int N = distance(p,q);

            if (multiplicity(N)) {
            
              numberOfSignalEvents += 1;

              if (p->getT() > t1 + deadTime_ns) {
            
                for (vector<hit_type>::const_iterator __p = p; __p != q; ++__p) {
                  for (vector<hit_type>::const_iterator __q = __p; ++__q != q; ) {

                    if ((__p->getToT() >= totRange_ns[0])    &&
                        (__q->getToT() >= totRange_ns[0])    &&
                        (__p->getToT() >= totRange_ns[1] ||
                         __q->getToT() >= totRange_ns[1])    &&
                        (__p->getToT() <= totRange_ns[2])    &&
                        (__q->getToT() <= totRange_ns[2])) {

                      h2s->Fill((double) combinatorics.getIndex(__p->getPMT(),__q->getPMT()), 
                                JCombinatorics::getSign(__p->getPMT(),__q->getPMT()) * (__q->getT() - __p->getT()));
                    }
                  }
                }

                t1 = p->getT();
              }
            }

            p = q;
          }

          // Background;
          // Note that rate cuts and veto have already been accounted for when data buffer was filled

          if        (background == rndms_t) {

            for (vector<hit_type>::iterator i = data.begin(); i != data.end(); ++i) {
              *i = hit_type(i->getPMT(), JHit(i->getT() + t0[i->getPMT()], i->getToT()));
            }
            
            sort(data.begin(), data.end());

            double t1 = numeric_limits<double>::lowest();

            for (vector<hit_type>::const_iterator p = data.begin(); p != data.end(); ) {
            
              vector<hit_type>::const_iterator q = p;
               
              while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}

              const int N = distance(p,q);

              if (multiplicity(N)) {

                if (p->getT() > t1 + deadTime_ns) {

                  for (vector<hit_type>::const_iterator __p = p; __p != q; ++__p) {
                    for (vector<hit_type>::const_iterator __q = __p; ++__q != q; ) {
                  
                      if ((__p->getToT() >= totRange_ns[0])    &&
                          (__q->getToT() >= totRange_ns[0])    &&
                          (__p->getToT() >= totRange_ns[1] ||
                           __q->getToT() >= totRange_ns[1])    &&
                          (__p->getToT() <= totRange_ns[2])    &&
                          (__q->getToT() <= totRange_ns[2])) {

                        h1b->Fill((double) combinatorics.getIndex(p->getPMT(),q->getPMT()), 1.0);
                      }
                    }
                  }

                  t1 = p->getT();
                }
              }

              p = q;
            }

          } else if (background == rates_t ||
                     background == count_t) {

            double Rs = 0.0;                                       // total rate [Hz]
            
            for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
              if (!veto[i]) {
                Rs += rate_Hz[i];                                  // [Hz]
              }
            }
            
            for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
              for (int j = i; ++j != NUMBER_OF_PMTS; ) {
                
                if (!veto[i] && !veto[j]) {

                  const double R1 = rate_Hz[i];                    // [Hz]
                  const double R2 = rate_Hz[j];                    // [Hz]

                  // evaluate expected counts within a time window of 2 * TMax_ns for the two PMTs and the other PMTs inside the optical module;
                  // expected rate due to random coincidences then is product of the probability for the specific observation and the number of trials
                  // 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
                  // corresponding to the given multiplicity range 

                  const double N = getP((R1)           * 2 * TMax_ns * 1e-9,
                                        (R2)           * 2 * TMax_ns * 1e-9,
                                        (Rs - R1 - R2) * 2 * TMax_ns * 1e-9,
                                        multiplicity.getLowerLimit(), 
                                        multiplicity.getUpperLimit()) * getFrameTime() / (2*TMax_ns);

                  h1b->Fill((double) combinatorics.getIndex(i,j), N);
                }
              }
            }
          }

          //-----------------------------------------------------------
          // fill the livetime by PMT pairs
          //-----------------------------------------------------------
        
          const double livetime_s = getFrameTime()*1e-9  *  exp(-deadTime_ns * numberOfSignalEvents/getFrameTime());

          for (size_t i = 0; i < combinatorics.getNumberOfPairs(); ++i) {

            const int pmt1 = combinatorics.getPair(i).first;
            const int pmt2 = combinatorics.getPair(i).second;

            if (!veto[pmt1] && !veto[pmt2]) {
              h1L->Fill(i, livetime_s);
            }
          }
        }
      }
    }
  }
  STATUS(endl);

  if (background == tails_t ) {

    const double R = (2.0*TMax_ns) / ((ymax - ymin)/ny);

    for (vector<JHistogram>::iterator hr = zmap.begin() ; hr != zmap.end() ; ++hr) {

      if (hr->is_valid()) {

        TH2D* h2s = hr->h2s;
        TH1D* h1b = hr->h1b;

        for (int ix = 1; ix <= h1b->GetXaxis()->GetNbins(); ++ix) {

          double Y = 0.0;                          // sum of counts  in tail regions
          double W = 0.0;                          // square of error of Y
          int    N = 0;                            // number of bins in tail regions

          for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {

            const double y = h2s->GetYaxis()->GetBinCenter(iy) ;

            if (Tail_ns(fabs(y))) { 
              Y += h2s->GetBinContent(ix,iy);
              W += h2s->GetBinError  (ix,iy) * h2s->GetBinError(ix,iy);
              N += 1;
            }
          }

          h1b->SetBinContent(ix,      Y/N  * R);
          h1b->SetBinError  (ix, sqrt(W/N) * R);
        }
      }
    }
  }

  //---------------------------------------------
  // save normalisation constants and store histograms
  //---------------------------------------------

  TH1D weights_hist(weights_hist_t, NULL, 3, 0.0, 3.0); 

  weights_hist.GetXaxis()->SetBinLabel(1, W1_overall_t);         // [s]
  weights_hist.GetXaxis()->SetBinLabel(2, WS_t);                 // [ns]
  weights_hist.GetXaxis()->SetBinLabel(3, WB_t);                 // [ns]

  weights_hist.Fill(W1_overall_t, counter*getFrameTime()*1e-9);  // [s]
  weights_hist.Fill(WS_t,         (ymax - ymin)/ny);             // [ns]
  weights_hist.Fill(WB_t,         2.0*TMax_ns);                  // [ns]

  out << weights_hist;

  for (vector<JHistogram>::iterator i = zmap.begin(); i != zmap.end(); ++i) {
    if (i->is_valid()) {
      out << *(i->h2s);
      out << *(i->h1b);
      out << *(i->h1L);
    }
  }

  out.Write();
  out.Close();
}