JCalibrateK40.cc

Go to the documentation of this file.

#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 "JTrigger/JTriggerToolkit.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"

namespace {

  using JLANG::JObjectID;
  using JTOOLS::JCombinatorics;
  using JTOOLS::JAbstractHistogram;
  using JDETECTOR::JModule;
      
  // Methods for background estimation
  
  static const char* const rates_t  =  "rates";      //!< singles rates from summary data
  static const char* const count_t  =  "counts";     //!< singles rates from L0 data
  static const char* const tails_t  =  "tails";      //!< tails of time residual distribution
  static const char* const rndms_t  =  "randoms";    //!< random time offsets (only L0 data)

  /**
   * Maximal time-over-threshold.
   */
  static const double ToTmax_ns = 250;

  /**
   * Auxiliary data structure for histogram management.
   */
  struct JHistogram : 
    public JObjectID,
    public JCombinatorics
  {
    /**
     * Default constructor.
     */
    JHistogram() :
      h2s(NULL),
      h1b(NULL),
      h1L(NULL)
    {}


    /**
     * Constructor.
     *
     * \param  module      module
     * \param  axis        axis
     */
    JHistogram(const JModule& module, const JAbstractHistogram<Double_t>& axis) :
      JObjectID(module.getID())
    {
      using namespace JPP;

      // sort the PMT pairs by opening angle in the order largest angle -> smallest angle
    
      this->configure(module.size());

      this->sort(JPairwiseComparator(module));

      const Int_t    nx   = this->getNumberOfPairs();
      const Double_t xmin = -0.5;
      const Double_t xmax = nx - 0.5;

      h2s = new TH2D(MAKE_CSTRING(module.getID() << _2S), NULL, nx, xmin, xmax, axis.getNumberOfBins(), axis.getLowerLimit(), axis.getUpperLimit());
      h1b = new TH1D(MAKE_CSTRING(module.getID() << _1B), NULL, nx, xmin, xmax);
      h1L = new TH1D(MAKE_CSTRING(module.getID() << _1L), NULL, nx, xmin, xmax);

      h2s->Sumw2();
      h1b->Sumw2();
      h1L->Sumw2();
    }

    /**
     * Check validity.
     *
     * \return             true if valid; else false
     */
    bool is_valid() const
    {
      return (h2s != NULL &&
              h1b != NULL &&
              h1L != NULL);
    }

    TH2D* h2s;           //!< time distribution per PMT pair 
    TH1D* h1b;           //!< background rates  per PMT pair 
    TH1D* h1L;           //!< livetimes         per PMT pair
  };


  /**
   * Get coincidence probability of two PMTs within one module due to random background.
   *
   * \param  E1          expected counts of 1st PMT
   * \param  E2          expected counts of 2nd PMT
   * \param  ED          expected counts of module, excluding PMT 1 and PMT 2
   * \param  M_min       multiplicity lower limit
   * \param  M_max       multiplicity upper limit
   * \return             probability
   */
  inline double getP(const double E1, const double E2, const double ED, const int M_min, const int M_max)
  {
    // 1st, always calculate the probability for twofold multiplicity

    double P = E1 * E2 * exp(-(E1 + E2 + ED));
      
    // 2nd, calculate the probability of multiplicity at lower limit

    int MD = 1;

    for ( ; MD + 2 <= M_min; ++MD) {
      P *= ED/MD;
    }
  
    // 3rd, add the probabilities of multiplicities up to upper limit

    double P_sum = P;

    for ( ; MD + 2 <= M_max; ++MD) {
      P     *= ED/MD;
      P_sum += P;
    }
  
    return P_sum;
  } 
}


/**
 * \file
 *
 * Auxiliary program to determine PMT parameters from K40 data. 
 *
 * By default, the combinatorial background is estimated from the singles rate.\n
 * In case L0 data are taken (and no summary data are available),
 * the singles rate can be determined from the amount of L0 data.\n
 * One can also estimate the background from the tails in the time residual distribution.\n
 * In that case, option -B can be used to specify the minimal and maximal time offset in ns.\n
 * For example -B "10 20".\n
 * 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.\n
 * The time window is applied symmetrically around a time offset of zero.\n
 * 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.\n
 * 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
 * Furthermore, the time-over-threshold of one of the two hits should be greater or equals the second value.\n
 * 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.\n
 * For example -M "2 2" will strictly select two-fold coincidences.
 * Note that such a selection can bias the result.
 *
 * Also consult <a href="https://common.pages.km3net.de/jpp/JMonitor.PDF">documentation</a>.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  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();
}