JCalibrateK40.cc

Go to the documentation of this file.

#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 "km3net-dataformat/online/JDAQ.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/JSummaryRouter.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/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"

namespace {
      
  // Methods for background estimation
  
  static const char* rates_t  =  "rates";      //!< singles rates from summary data
  static const char* count_t  =  "counts";     //!< singles rates from L0 data
  static const char* tails_t  =  "tails";      //!< tails of time residual distribution
  static const char* rndms_t  =  "randoms";    //!< random time offsets (only L0 data)

  /**
   * Maximal time-over-threshold.
   */
  static const double ToTmax_ns = std::numeric_limits<KM3NETDAQ::JDAQHit::JTOT_t>::max();

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


    /**
     * Constructor.
     *
     * \param __h2s        2D histogram for signal
     * \param __h1b        1D histogram for background rates of PMT pairs
     * \param __h1L        1D histogram for livetimes of PMT pairs
     */
    JHistogram(TH2D* __h2s,
               TH1D* __h1b,
               TH1D* __h1L) :
      h2s(__h2s),
      h1b(__h1b),
      h1L(__h1L)
    {
      h2s->Sumw2();
      h1b->Sumw2();
      h1L->Sumw2();
    }

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

    TH2D* h2s;
    TH1D* h1b; //background rates per PMT pair 
    TH1D* h1L; //livetimes per PMT pair
  };
}


/**
 * \file
 *
 * 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 -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 sample.
 *
 * 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;


  //-----------------------------------------------------------
  //parameter interface
  //-----------------------------------------------------------

  JMultipleFileScanner<JDAQTimesliceTypes_t> inputFile;
  JLimit_t& numberOfEvents = inputFile.getLimit();
  string             outputFile;
  string             detectorFile;
  Double_t           TMax_ns;
  JRange<double>     rateRange_Hz;
  JRange<double>     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].")         = JRange<double>(4.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();
    zap['d'] = make_field(debug,          "debug flag.")                                 = 1;
     
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

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

  cout.tie(&cerr);

  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) {
    FATAL("Invalid option -M <multiplicity> " << multiplicity << endl);
  }

  if (!totRange_ns.is_valid()) {
    FATAL("Invalid option -t <totRange_ns> " << 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);
  JSummaryRouter      summaryRouter;

  //-----------------------------------------------------------
  // 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.getLowerLimit()), cpu.getNPE(totRange_ns.getUpperLimit()), NPE) 
                      /
                      cpu.getIntegralOfChargeProbability(cpu.getNPE(0.0),                         cpu.getNPE(ToTmax_ns),                   NPE));
  


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

  JCombinatorics combinatorics(NUMBER_OF_PMTS);

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

  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]

  const int    nx   = combinatorics.getNumberOfPairs();
  const double xmin = -0.5;
  const double xmax = nx - 0.5;
  
  const double ymin = -floor(TMax_ns) + 0.5;
  const double ymax = +floor(TMax_ns) - 0.5;
  const int    ny   = (int) ((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(new TH2D(MAKE_CSTRING(module->getID() << _2S), NULL, nx, xmin, xmax, ny, ymin, ymax),
                                       new TH1D(MAKE_CSTRING(module->getID() << _1B), NULL, nx, xmin, xmax),
                                       new TH1D(MAKE_CSTRING(module->getID() << _1L), NULL, nx, xmin, xmax));
    }
  }


  JTreeScanner<JDAQSummaryslice, JDAQEvaluator> summaryFile(inputFile);

  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;

  JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(inputFile, selector);
  
  counter_type counter = 0;

  for ( ; in.hasNext() && counter != inputFile.getLimit(); ++counter) {

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

    const JDAQTimeslice* timeslice = in.next();

    if (background == rates_t) {

      const Long64_t index = summaryFile.find(*timeslice);
      
      if (index == -1) {
        FATAL("Missing summary data at " << timeslice->getFrameIndex() << endl);
      }

      summaryRouter.update(summaryFile.getEntry(index));

      if (timeslice->getFrameIndex() != summaryRouter.getFrameIndex()) {

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

        if (!in.hasNext())
          break;
        else
          FATAL("ROOT TTrees not aligned; Run number " << timeslice->getRunNumber() << endl);
      }
    }


    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
        //-----------------------------------------------------------

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

        JDAQFrameStatus status;
        
        if        (background == count_t) {
        
          status = frame->getDAQFrameStatus();

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

        } else if (background == tails_t) {

          status = frame->getDAQFrameStatus();

        } else if (background == rndms_t) {

          status = frame->getDAQFrameStatus();

        } else if (background == rates_t) {

          if (!summaryRouter.hasSummaryFrame(frame->getModuleID())) {
            FATAL("Summary frame not found " << timeslice->getFrameIndex() << ":" << frame->getModuleID() << endl);
          }

          const JDAQSummaryFrame& summary_frame = summaryRouter.getSummaryFrame(frame->getModuleID());

          status = summary_frame.getDAQFrameStatus();

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

        // Set veto according DAQ or PMT status
        // 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)) {
            veto[i] = true;
          }
        }

        //-----------------------------------------------------------
        // sort the PMT pairs by opening angle in the order largest angle-->smallest angle
        //-----------------------------------------------------------

        const JModuleAddress& address = router.getAddress(frame->getModuleID());

        combinatorics.configure(module.size());

        combinatorics.sort(JPairwiseComparator(module));

        TH2D* h2s = zmap[address.first].h2s;
        TH1D* h1b = zmap[address.first].h1b;
        TH1D* h1L = zmap[address.first].h1L;

        //-----------------------------------------------------------
        // fill the livetime by PMT pairs
        //-----------------------------------------------------------
        
        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] && rateRange_Hz(rate_Hz[pmt1]) &&
              !veto[pmt2] && rateRange_Hz(rate_Hz[pmt2]) ) {
            
            h1L->Fill(i, getFrameTime()*1e-9);
          }
        }

        //-----------------------------------------------------------
        // 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) {

          const int pmt = i->getPMTAddress();
          
          if ( veto[pmt] || !rateRange_Hz(rate_Hz[pmt]) ) {
            i->reset();
          }
        }

        JSuperFrame1D_t& data = JSuperFrame1D_t::multiplex(buffer);
        
        data.pop_back();   // remove end marker

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

        if (data.empty()) { 
          continue;
        }

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

        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 (totRange_ns( __p->getToT()) && totRange_ns(__q->getToT())) {
                    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());

          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)) {

              for (vector<hit_type>::const_iterator __p = p; __p != q; ++__p) {
                for (vector<hit_type>::const_iterator __q = __p; ++__q != q; ) {
                  
                  if (totRange_ns(__p->getToT()) && totRange_ns(__q->getToT())) {
                    h1b->Fill((double) combinatorics.getIndex(p->getPMT(),q->getPMT()), 1.0);
                  }
                }
              }
            }

            p = q;
          }

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

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

                // expectation values (counts) within TMax_ns * 2 for the optical module and the two PMTs
                const double ED = (Rs - rate_Hz[i] - rate_Hz[j]) * 2 * TMax_ns * 1e-9;
                const double E1 = rate_Hz[i] * 2 * TMax_ns * 1e-9;
                const double E2 = rate_Hz[j] * 2 * TMax_ns * 1e-9;

                // expected rate of coincidences = P * # of trials
                Double_t N = coincidenceP(E1, E2, ED, 
                                          multiplicity.getLowerLimit(), 
                                          multiplicity.getUpperLimit()) * getFrameTime() / (2*TMax_ns);

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

  if (background == tails_t ) {

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

      if (hist->is_valid()) {

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

        for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
          for (int j = i; ++j != NUMBER_OF_PMTS; ) {

            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
            int    k = combinatorics.getIndex(i,j);  // PMT pair index (x-axis of h2s)

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

              const Double_t dt = h2s->GetYaxis()->GetBinCenter(l) ;

              if (Tail_ns(fabs(dt))) { 
                Y += h2s->GetBinContent(k+1,l);
                W += h2s->GetBinError(k+1,l) * h2s->GetBinError(k+1,l);
                N++;
              }
            }

            h1b->SetBinContent(k+1, Y / N     * (2.0*TMax_ns) / ((ymax - ymin)/ny) );
            h1b->SetBinError  (k+1, sqrt(W/N) * (2.0*TMax_ns) / ((ymax - ymin)/ny) );
          }
        }
      }
    }
  }

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

  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]


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

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

  weights_hist.Write() ;

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

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