JHalibut.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <map>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TF1.h"

#include "JDAQ/JDAQTimesliceIO.hh"
#include "km3net-dataformat/online/JDAQClock.hh"
#include "JDAQ/JDAQEvaluator.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"

#include "JTrigger/JHitR0.hh"
#include "JTrigger/JMatchL0.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JSuperFrame1D.hh"
#include "JTrigger/JSuperFrame2D.hh"

#include "JROOT/JManager.hh"
#include "JCalibrate/JCalibrateK40.hh"

#include "JROOT/JROOTClassSelector.hh"
#include "JTools/JQuantile.hh"
#include "JTools/JRange.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JAutoTreeScanner.hh"
#include "JSupport/JSupport.hh"

#include "JLang/JObjectMultiplexer.hh"
#include "JMath/JMathToolkit.hh"

#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 * Example program to determine N-fold coincidence rates.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  typedef JRange<int>    JRange_t;

  JMultipleFileScanner<> inputFile;
  JLimit_t&              numberOfEvents = inputFile.getLimit();
  string                 outputFile;
  string                 detectorFile;
  double                 TMax_ns;
  JROOTClassSelector     selector;
  JRange_t               M;
  string                 summaryFile;
  string                 option;
  int                    debug;

  try { 

    JParser<> zap("Example program to determine N-fold coincidence rates.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "halibut.root";
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['a'] = make_field(detectorFile);
    zap['T'] = make_field(TMax_ns)             = 20.0;
    zap['C'] = make_field(selector)            = getROOTClassSelection<JDAQTimesliceTypes_t>();
    zap['M'] = make_field(M)                   = JRange_t(2, 31);
    zap['s'] = make_field(summaryFile)         = "halibut.txt";
    zap['O'] = make_field(option)              = "";
    zap['d'] = make_field(debug)               = 1;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  JDetector detector;

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

  const JModuleRouter router(detector);

  map<int, double>    TMax_s;             // histogram upper limit

  TMax_s[2] = 5.0e-2;
  TMax_s[3] = 0.2e+0;
  TMax_s[4] = 1.0e+0;
  TMax_s[5] = 1.0e+1;
  TMax_s[6] = 1.0e+2;

  typedef JSuperFrame2D<JHitR0>   JSuperFrame2D_t;
  typedef JSuperFrame1D<JHitR0>   JSuperFrame1D_t;
      
  const JMatchL0<JHitR0>      match(TMax_ns);     // time window self-coincidences [ns]

  typedef JManager<int, TH1D>  JManager_t;

  JManager_t H1(new TH1D("H1[%]", NULL, 100, -TMax_ns, +TMax_ns));
  JManager_t T1(new TH1D("T1[%]", NULL, 100, 0.0, TMax_s[M.getLowerLimit()]));

  H1->Sumw2();
  T1->Sumw2();

  JTreeScanner<>::debug = debug;

  JAutoTreeScanner<JDAQTimeslice, JDAQEvaluator> zmap = JType<JDAQTimesliceTypes_t>();

  JTreeScannerInterface<JDAQTimeslice>* ps = zmap[selector];

  counter_type counter = 0;

  // process file-by-file to speed up JTreeScanner::configure();

  for (JMultipleFileScanner<>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {

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

    ps->configure(*i);

    vector<double> t0(detector.size(), 0.0);

    for ( ; ps->hasNext() && counter != inputFile.getLimit(); ++counter) {

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

      const JDAQTimeslice* timeslice = ps->next();

      for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
        
        JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, router.getModule(frame->getModuleID()));
        
        for (JSuperFrame2D_t::iterator i = buffer.begin(); i != buffer.end(); ++i) {
          i->join(match);
        }

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

        if (data.size() > 1) {
          
          TH1D* h1 = H1[frame->getModuleID()];
          TH1D* t1 = T1[frame->getModuleID()];

          for (vector<JHitR0>::const_iterator p = data.begin(); p != data.end(); ) {
            
            vector<JHitR0>::const_iterator q = p;
            
            while (++q != data.end() && q->getT() - p->getT() <= TMax_ns ) {}
            
            const int N = distance(p,q);
            
            if (M(N)) {
              
              const int    i  = router.getIndex(frame->getModuleID());
              const double ts = getTimeOfRTS(frame->getFrameIndex()) + p->getT();
              
              t1->Fill((ts - t0[i]) * 1.0e-9);      // [s]
            
              t0[i] = ts;
              
              const double W = factorial(N, 2);
              
              for (vector<JHitR0>::const_iterator __p = p; __p != q; ++__p) {
                for (vector<JHitR0>::const_iterator __q = __p; ++__q != q; ) {
                  h1->Fill(JCombinatorics::getSign(__p->getPMT(),__q->getPMT()) * (__q->getT() - __p->getT()), 1.0/W);
                }
              }
            }
            
            p = q;
          }
        }
      }
    }
    STATUS(endl);
  }

  if (counter != 0) {

    const double V = (H1->GetXaxis()->GetXmax() - H1->GetXaxis()->GetXmin()) / (double) H1->GetXaxis()->GetNbins();  // [ns]
    const double W = counter * getFrameTime() * 1.0e-9;                                                              // [s]

    for (JManager_t::iterator i = H1.begin(); i != H1.end(); ++i) {
      i->second->Scale(1.0/(V*W));
    }

    for (JManager_t::iterator i = T1.begin(); i != T1.end(); ++i) {
      i->second->Scale(1.0/i->second->GetMaximum());
    }
  }

  if (summaryFile != "") {

    const double V = (H1->GetXaxis()->GetXmax() - H1->GetXaxis()->GetXmin()) / (double) H1->GetXaxis()->GetNbins();  // [ns]

    const int number_of_strings =  getNumberOfStrings(detector);
    const int number_of_floors  =  getNumberOfFloors (detector);
    const int PRECISION         = (M.getLowerLimit() > 2 ? 4 : 3);

    ofstream out(summaryFile.c_str());

    out << "Multiplicity " << M                                      << endl;
    out << "-------------------------------------------------------" << endl;
    out << " location  |   Gauss  |   S - B  |   Total  |   slope  " << endl;
    out << "           |   [Hz]   |   [Hz]   |   [Hz]   |   [Hz]   " << endl;
    out << "-------------------------------------------------------" << endl;

    JQuantile Q[4];

    for (int string = 1; string <= number_of_strings; ++string) {
      for (int floor = number_of_floors; floor >= 1; --floor) {

        const int id = detector.getModule(JLocation(string,floor)).getID();

        out << "  " << setw(3) << string << ' ' << setw(2) << floor << "  ";
          
        TH1D* h1 = (H1.find(id) != H1.end() ? H1[id] : NULL);
        TH1D* t1 = (T1.find(id) != T1.end() ? T1[id] : NULL);

        if (h1 != NULL) {

          TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2*TMath::Pi())*[2]) + [3]");
          
          f1.SetParameter(0, h1->GetMaximum());
          f1.SetParameter(1, 0.0);
          f1.SetParameter(2, h1->GetRMS() * 0.25);
          f1.SetParameter(3, h1->GetMinimum());
          
          h1->Fit(&f1, option.c_str(), "same");
          
          out << " | " << FIXED(8,PRECISION) <<  f1.GetParameter(0);
          out << " | " << FIXED(8,PRECISION) << (h1->GetSumOfWeights() - f1.GetParameter(3) * h1->GetNbinsX()) * V;
          out << " | " << FIXED(8,PRECISION) <<  h1->GetSumOfWeights() * V;
          
          Q[0].put( f1.GetParameter(0));
          Q[1].put((h1->GetSumOfWeights() - f1.GetParameter(3) * h1->GetNbinsX()) * V);
          Q[2].put( h1->GetSumOfWeights() * V);
        }

        if (t1 != NULL) {
          
          TF1 f1("f1", "[0]*exp(-[1]*x)");
          
          f1.SetParameter(0, t1->GetMaximum());
          f1.SetParameter(1, 1.0 / t1->GetRMS());
          
          t1->Fit(&f1, option.c_str(), "same");
          
          out << " | " << FIXED(8,PRECISION) <<  f1.GetParameter(1);
          
          Q[3].put(f1.GetParameter(1));
        }

        out << endl;
      }

      out << endl;
    }

    if (Q[0].getCount() != 0) {

      out << "-------------------------------------------------------" << endl;
      out << setw(10) << left << " average";
      
      for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
        out << " | " << FIXED(8,PRECISION) << Q[i].getMean();
      }

      out << endl;
    }

    out.close();
  }

  if (outputFile != "") {

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

    H1.Write(out);
    T1.Write(out);

    out.Close();
  }
}