JMonitorL1dt.cc

Go to the documentation of this file.

#include "km3net-dataformat/online/JDAQ.hh"
#include "km3net-dataformat/online/JDAQTimeslice.hh"
#include "km3net-dataformat/online/JDAQSummaryslice.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JDAQ/JDAQEvaluator.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JTrigger/JHitL1.hh"
#include "JTrigger/JBuildL2.hh"
#include "JTrigger/JPMTSelector.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JSupport.hh"
#include "JLang/JObjectMultiplexer.hh"
#include "JROOT/JROOTClassSelector.hh"
#include "Jeep/JParser.hh"
#include "JROOT/JManager.hh"
#include "JTools/JCombinatorics.hh"

#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TMath.h"

#include <string>
#include <iostream>


using namespace JTOOLS;
namespace MONITORL1DT {

  /**
   * Data structure for hit time and DOM identifier.
   */
  class JElement {
    public:
      JElement() :
        id(0),
        t (0.0)
    {}

      JElement(const int    __id,
          const double __t) :
        id(__id),
        t (__t)
    {}

      int    id;
      double t;
  };


  /**
   * Sort JElement in time.
   *
   * \param  first        first  element
   * \param  second       second element
   * \return              true if second later than first; else false
   */
  inline bool operator<(const JElement& first, const JElement& second)
  {
    return first.t < second.t;
  }


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


    /**
     * Constructor.
     *
     * \param __h2s   2D histogram for signal
     * \param __h1l   1D histogram for background
     */
    JHistogram(TH2D* __h2s,
        TH1D* __h1l) :
      h2s(__h2s),
      h1l(__h1l)
    {
      h2s->Sumw2();
      h1l->Sumw2();
    }

    TH2D* h2s;
    TH1D* h1l;
  };

  // Taken from C. Pieterse nanobeacon code
  typedef JCombinatorics::pair_type pair_type;
  inline bool comparepair(const pair_type& A, const pair_type& B) {
    return( abs(A.first - A.second) > abs(B.first - B.second) ) ;
  }
}



/**
 * \file
 *
 * Auxiliary program to plot L1 time difference correlations between DOMs.
 *
 * \author kmelis, lnauta
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace MONITORL1DT;
  using namespace JPP;

  JMultipleFileScanner<JDAQTimesliceTypes_t> inputFile;
  string             outputFile;
  string             detectorFile;
  double             Timewindow_ns;
  int                binwidth; // since the data is taken on a ns level, having bins of non-integer value is non-sensical
  double             TmaxL1_ns;
  JROOTClassSelector selector;
  unsigned int       multiplicity;
  bool               correct_time;
  double             livetime_s;
  int                debug;

  try {

    JParser<> zap("Program to create L1 hit time difference histograms from raw data.");

    zap['f'] = make_field(inputFile,     "input file");
    zap['o'] = make_field(outputFile,    "output file")                                 = "monitor.root";
    zap['a'] = make_field(detectorFile,  "detector file");
    zap['t'] = make_field(TmaxL1_ns,     "max time between L1 hits [ns]")               = 1000.0;
    zap['T'] = make_field(Timewindow_ns, "time window around t=0 [ns]")                 = 2400.0;
    zap['w'] = make_field(binwidth,      "binwidth [ns]")                               = 1;
    zap['C'] = make_field(selector,      "datastream selector")                         = getROOTClassSelection<JDAQTimesliceTypes_t>();
    zap['m'] = make_field(multiplicity,  "minimal multiplicity of the L1 hits")         = 2;
    zap['c'] = make_field(correct_time,  "subtract expected arrival time from delta-t");
    zap['L'] = make_field(livetime_s,    "livetime of the data, set to positive value") = -1.0; // set to positive value if you want to set the livetime manually (i.e. MC)
    zap['d'] = make_field(debug) = 1;

    if (zap.read(argc, argv) != 0)
      return 1;
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  cout.tie(&cerr);

  JDetector detector;

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

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

  const JModuleRouter router(detector);

  typedef vector<JHitL1> JFrameL1_t;

  const double ctMin = -1; // full angle over the DOM
  const JBuildL2<JHitL1> buildL2(JL2Parameters(multiplicity, TmaxL1_ns, ctMin));

  // Manager histograms for single-DOM output
  vector<JHistogram> zmap(detector.size());  

  const int    nx   = detector.size();
  const double xmin = -0.5;
  const double xmax = nx - 0.5;

  const double ymin = -floor(Timewindow_ns) + 0.5;
  const double ymax = +floor(Timewindow_ns) + 0.5; // changed - to + to get an even number for Nbins, rebin needs numbers that have divisors.
  const int    ny   = (int) ((ymax - ymin) / binwidth);

  // Manage histograms for single-DU output
  Int_t num_of_floors = getNumberOfFloors(detector);
  JCombinatorics c(num_of_floors);
  Int_t npairs = c.getNumberOfPairs();
  c.sort(MONITORL1DT::comparepair);

  JManager<int, TH2D>* manager;
  manager = new JManager <int, TH2D>(new TH2D("h%", "", npairs, 0.5, npairs+0.5, ny, ymin, ymax));

  NOTICE("Running JMonitorL1dt: Monitoring L1 time differences and creating histograms." << endl);
  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {

    const JModuleAddress& address = router.getAddress(module->getID());

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

    const JString title(module->getID());
    JString titleString1D;
    JString titleString2D;
    titleString1D = title + ".1L";
    titleString2D = title + ".2S";

    zmap[address.first] = JHistogram(new TH2D((titleString2D).c_str(), NULL, nx, xmin, xmax, ny, ymin, ymax),
        new TH1D((titleString1D).c_str(), NULL, nx, xmin, xmax));

    for (JDetector::iterator mod = detector.begin(); mod != detector.end(); ++mod) {
      zmap[address.first].h2s->GetXaxis()->SetBinLabel(distance(detector.begin(), mod)+1, Form("%i", mod->getID()));
      zmap[address.first].h1l->GetXaxis()->SetBinLabel(distance(detector.begin(), mod)+1, Form("%i", mod->getID()));
    }
  }


  JTreeScanner<JDAQSummaryslice, JDAQEvaluator> summaryFile(inputFile);

  JLANG::JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(inputFile, selector);

  int counter = 0;

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

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

    const JDAQTimeslice* timeslice = in.next();

    JFrameL1_t frameL1;
    vector<MONITORL1DT::JElement> buffer;
    vector<bool> DOM_OK(detector.size(), true);

    for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) {
      if (router.hasModule(super_frame->getModuleID()) && !super_frame->empty()) {

        const JModuleAddress& address  = router.getAddress(super_frame->getModuleID());
        const JModule&        module   = detector.getModule(address);

        frameL1.clear();

        buildL2(*super_frame, module, back_inserter(frameL1));

        for (JFrameL1_t::iterator L1hit = frameL1.begin(); L1hit != frameL1.end(); ++L1hit) {
          buffer.push_back(JElement(address.first, L1hit->begin()->getT()));
        }
      }
    }

    for (vector<JHistogram>::iterator h1 = zmap.begin(); h1 != zmap.end(); ++h1) {
      if (!DOM_OK[distance(zmap.begin(), h1)]) { 
        continue; 
      }
      for (unsigned int i = 0; i < detector.size(); ++i) {
        if (DOM_OK[i]) { h1->h1l->Fill(i, getFrameTime() * 1e-9); } // fill the th1 with the frametime (how long did it measure) in seconds, all weight will give total measured time
      }
    }

    // fill histograms with correlations
    sort(buffer.begin(), buffer.end());

    for (vector<JElement>::const_iterator p = buffer.begin(); p != buffer.end(); ) {
      vector<JElement>::const_iterator q = p;

      const JModule& module_p = detector.getModule((JModuleAddress)p->id);

      while (++q != buffer.end() && q->t - p->t <= Timewindow_ns ) {
        //        if (p->id == q->id) { continue; }

        const JModule& module_q = detector.getModule((JModuleAddress)q->id);

        double dom_distance    = getDistance(module_p.getPosition(), module_q.getPosition());
        double time_correction = (correct_time ? (dom_distance / getSpeedOfLight()) : 0);

        zmap[p->id].h2s->Fill(q->id, q->t - p->t - time_correction);
        zmap[q->id].h2s->Fill(p->id, p->t - q->t + time_correction);

        if (module_p.getString() == module_q.getString()) {
          // indices are pairs running from 0-17 between 0-152, while floors are 1-18 and bins are 1-153
          int xbin = c.getIndex(module_p.getFloor() - 1, module_q.getFloor() - 1) + 1;
          (*manager)[module_p.getString()]->Fill(xbin, q->t - p->t - time_correction);
        }
      }
      p++;  // move to the next L1 hit
    }
    buffer.clear();
  }
  STATUS(endl);

  // livetime in case set manually
  if (livetime_s > 0.0) {
    for (vector<JHistogram>::iterator i = zmap.begin(); i != zmap.end(); ++i) {
      TH1D* hl = i->h1l;
      for (int ibin = 1; ibin <= hl->GetNbinsX(); ++ibin) {
        hl->SetBinContent(ibin, livetime_s);
        hl->SetBinError(ibin, 0.0000001);
      }
    }
  }

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

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

  manager->Write(out);
  out.Write();
  out.Close();
}