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 "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 "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 -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 documentation.

Author

mdejong

Definition in file JCalibrateK40.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 140 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;
  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);

  //-----------------------------------------------------------
  // 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());  

  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));
    }
  }


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

          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 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] = RTU * sum.getRate(i);
              }

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

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

          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

          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 (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());

            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())) {
                        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_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);
                }
              }
            }
          }

          //-----------------------------------------------------------
          // 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] && rateRange_Hz(rate_Hz[pmt1]) &&
                !veto[pmt2] && rateRange_Hz(rate_Hz[pmt2]) ) {
            
              h1L->Fill(i, livetime_s);
            }
          }
        }
      }
    }
  }
  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
  //---------------------------------------------

  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();
}