JDomAnalyser.cc

Go to the documentation of this file.

// c++ standard library
#include <iostream>
#include <iomanip>

#include "JSupport/JSupport.hh"
#include "JSupport/JSupportToolkit.hh"
#include "JSupport/JMeta.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JAutoTreeScanner.hh"

#include "JROOT/JROOTClassSelector.hh"
#include "JLang/JObjectMultiplexer.hh"

#include "JROOT/JManager.hh"

#include "Jeep/JParser.hh"
#include "km3net-dataformat/online/JDAQPMTIdentifier.hh"
#include "km3net-dataformat/online/JDAQ.hh"
#include "JDAQ/JDAQTimesliceIO.hh"

#include "TH1D.h"
#include "TMath.h"

using namespace KM3NETDAQ;
using namespace JLANG;   
using namespace JPP;
using namespace std;

/*
 * Rebins a histogram with constant bin width in a linear scale, so that the bins will have constant width in log10 scale. 
 *
 * \param h The histogram.
 */
inline void BinLogX (TH1* h){
  
  TAxis *axis      = h->GetXaxis();
  int bins         = axis->GetNbins();
  Axis_t min       = axis->GetXmin();
  Axis_t max       = axis->GetXmax();
  Axis_t width     = (max - min) / bins;
  Axis_t *new_bins = new Axis_t[bins + 1];

  for (int i = 0; i <= bins; i++) {
    new_bins[i] = TMath::Power(10, min + i * width);
  }
  
  axis->Set(bins, new_bins);
  delete new_bins;
}


int main(int argc, char **argv) {

  JMultipleFileScanner<JDAQTimesliceTypes_t> inputFile;
  JLimit_t&                                  nSlices = inputFile.getLimit();
  JROOTClassSelector                         selector;
  string                                     outputFile;
  string                                     detectorFile;

  try {

    JParser<> zap;
    
    zap['f'] = make_field(inputFile  , "Path to input file ");
    zap['o'] = make_field(outputFile , "Path to output file") = "out.root"; 
    zap['C'] = make_field(selector   , "timeslice selector, e.g. JDAQTimesliceL1.")  = getROOTClassSelection<JDAQTimesliceTypes_t>();
    zap['s'] = make_field(nSlices    , "number of slices"   ) = JLimit::max(); 

    zap(argc,argv);
  }
  catch(const exception &error) {
    ERROR(error.what() << endl);
  }

  const double factor = 1./(1.0e-6 * getFrameTime()); // [kHz]
  
  TH1D* h1 = new TH1D("rate_%", NULL, 100, -2, 2);
  TH1D* h2 = new TH1D("time_%", NULL, 1, -0.5, 0.5);

  JAutoTreeScanner<JDAQTimeslice, JDAQEvaluator> map = JType<JDAQTimesliceTypes_t>();                  
  JTreeScannerInterface<JDAQTimeslice>*            s = map[selector];        
  s->configure(inputFile);   

  if (s->hasNext()){
    JFrameIndexRange range (s->begin()->getFrameIndex() , s->rbegin()->getFrameIndex());
    h2->SetBins(range.second - range.first + 1, range.first - 0.5, range.second + 0.5);
  }

  BinLogX(h1);

  JManager<JDAQPMTIdentifier,TH1D> r (h1);
  JManager<JDAQPMTIdentifier,TH1D> t (h2);

  JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> scanner(inputFile, selector);

  counter_type counter = 0;

  for ( ; scanner.hasNext() && counter != inputFile.getLimit(); ++counter) {
             
    JDAQTimeslice* slice = scanner.next();

    for(JDAQTimeslice::const_iterator frame = slice->begin() ; frame != slice->end() ; ++frame) {
      
      vector <int>  pmt_hit_count (NUMBER_OF_PMTS , 0) ;

      for (JDAQSuperFrame::const_iterator hit = frame->begin() ; hit != frame->end() ; ++hit){
        pmt_hit_count[hit->getPMT()]++;
      }
        
      for (int pmt = 0 ; pmt != NUMBER_OF_PMTS ; ++pmt) {
        r[JDAQPMTIdentifier(frame->getModuleID(),pmt)]->Fill(pmt_hit_count[pmt]*factor);
        t[JDAQPMTIdentifier(frame->getModuleID(),pmt)]->Fill(slice->getFrameIndex(), pmt_hit_count[pmt]*factor);                  
      }      
    }
  }

  for (typename JManager < JDAQPMTIdentifier , TH1D >::const_iterator i = r.begin() ; i != r.end() ; ++i){

    for(int j=1 ; j < i -> second -> GetNbinsX() ; ++j){
      
      double width = i->second->GetXaxis()->GetBinWidth(j);
      i -> second->SetBinContent(j ,i->second->GetBinContent(j)/width);
      i -> second->SetBinError  (j ,i->second->GetBinError  (j)/width);
    }
  } 
  
  TFile out(outputFile.c_str(), "recreate");
  putObject(&out, JMeta(argc, argv));
  r.Write(out);
  t.Write(out);
}