JMergeRunAnalyzer.cc

Go to the documentation of this file.

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

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TKey.h"
#include "TString.h"
#include "TMath.h"
#include "TTimeStamp.h"
 
#include "km3net-dataformat/online/JDAQHeader.hh"
 
#include "JROOT/JRootFileReader.hh"
#include "JROOT/JRootFileWriter.hh"

#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSingleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JMeta.hh"

#include "JTrigger/JTriggerInterface.hh"
#include "JTrigger/JTriggerBits.hh"

#include "JGizmo/JGizmoToolkit.hh"

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

#include "JRunAnalyzer/JRunInfo.hh"

/*
 * \author adomi
 * Auxiliary program to merge JRunAnalyzer histograms.
 */

/* 
 * Fill TH2 from TH1 of each run.
 */
inline void fill2DHistogram(int run, TH1D* h1, TH2D* h2, bool w_mc, double daq_livetime = 1, double mc_livetime = 1){

  if(h1 != NULL && h2 != NULL){

    double w = 1;
      
    for(int i = 1; i <= h1->GetNbinsX(); ++i){
    
      if(w_mc){
        w = daq_livetime / mc_livetime;
      }
  
      h2->Fill(run, h1->GetXaxis()->GetBinCenter(i), h1->GetBinContent(i) * w);       
    }
  }
}


int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  vector<string>     inputFile;
  JFileRecorder<JTYPELIST<run_info, JMeta, TH1D, TH2D>::typelist>  outputFile;
  bool               weight_mc;
  int                debug;
 
  try { 
 
    JParser<> zap("Auxiliary program to merge JRunAnalyzer histograms.");
     
    zap['f'] = make_field(inputFile,      "input file (output from JRunAnalyzer).");
    zap['o'] = make_field(outputFile,     "output file.")   = "merge-jra.root";
    zap['w'] = make_field(weight_mc,      "weight mc events");
    zap['d'] = make_field(debug,          "debug.")         = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  int min_run = numeric_limits<int>::max(), max_run = 0;

  size_t min_date = numeric_limits<int>::max(), max_date = 0;

  run_info runinfo;

  double total_mc_livetime = 0;
  int run = 0;

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

    {
      JSingleFileScanner<run_info> in(*i);
      if(in.hasNext()){
        runinfo = *in.next();
        if(in.hasNext()){
          FATAL("MULTIPLE RUN INFO");
        }
        
        int run_id = runinfo.run_id;

        if(run_id != run){
          
          run = run_id;

          total_mc_livetime = runinfo.mc_livetime;
        
        } else {

        total_mc_livetime += runinfo.mc_livetime;
        }
      
      } else {
        FATAL("NO RUN INFO");
      }
    }

    for (JRootFileReader<JDAQHeader> in(i->c_str()); in.hasNext(); ) {

      const JDAQHeader* p = in.next();

      if(p->getTimesliceStart().getUTCseconds() < min_date) min_date = p->getTimesliceStart().getUTCseconds();
      if(p->getTimesliceStart().getUTCseconds() > max_date) max_date = p->getTimesliceStart().getUTCseconds();
      if(p->getRunNumber() < min_run) min_run = p->getRunNumber();
      if(p->getRunNumber() > max_run) max_run = p->getRunNumber();   
    }
  }
  
  time_t tmin = min_date, tmax = max_date;
  
  TTimeStamp date_min(tmin, 0);
  TTimeStamp date_max(tmax, 0);
  
  cout << "Minimum RUN: " << min_run << ", "; date_min.Print("s");
  cout << "Maximum RUN: " << max_run << ", "; date_max.Print("s");
  
  int xbins = (max_run - min_run) + 1;
  string title = string("Data taken from: ") + date_min.AsString("s") + " - " + date_max.AsString("s");
  
  TH2D* h_pmt_rate_distribution = new TH2D("h_pmt_rate_distribution", string(title + "; RUN; PMT Rate Distribution [kHz]").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, JDAQRate::getN(), JDAQRate::getData(1.0e-4));  
  TH2D* h_Trigger_bit_hit = new TH2D("h_Trigger_bit_hit", string(title + "; RUN; Hits trigger bit").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, NUMBER_OF_TRIGGER_BITS, -0.5, NUMBER_OF_TRIGGER_BITS - 0.5); 
  TH2D* h_Triggered_hits = new TH2D("h_Triggered_hits", string(title + "; RUN; Number of triggered hits").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5,  50, 0, 4);
  setLogarithmicY(h_Triggered_hits);
  TH2D* h_Triggered_hits_3dmuon = new TH2D("h_Triggered_hits_3dmuon", string(title + "; RUN; Number of triggered hits - JTRIGGER3DMUON").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 50, 0, 3);
  setLogarithmicY(h_Triggered_hits_3dmuon);
  TH2D* h_Snapshot_hits = new TH2D("h_Snapshot_hits", string(title + "; RUN; Number of snapshot hits").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5,  50, 0, 4);
  setLogarithmicY(h_Snapshot_hits);
  TH2D* h_Triggered_over_Snapshot_hits = new TH2D("h_Triggered_over_Snapshot_hits", string(title + "; RUN; Triggered/Snapshot hits").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 100, 0, 0.5);
  TH2D* h_Number_of_overlays = new TH2D("h_Number_of_overlays", string(title + ";RUN;Number of Overlays").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 1000, -0.5, 1000 - 0.5);
  TH2D* h_pmt_distribution_triggered_hits = new TH2D("h_pmt_distribution_triggered_hits", string(title + ";RUN;Normalised PMT Distrib triggered hits (upper PMTs=[0-11])").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, NUMBER_OF_PMTS, -0.5 , NUMBER_OF_PMTS - 0.5);
  TH2D* h_pmt_distribution_snapshot_hits = new TH2D("h_pmt_distribution_snapshot_hits", string(title + ";RUN; Normalised PMT Distrib snapshot hits (upper PMTs=[0-11])").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, NUMBER_OF_PMTS, -0.5 , NUMBER_OF_PMTS - 0.5);
  TH2D* h_event_duration = new TH2D("h_event_duration", string(title + ";RUN;Event Duration [ns]").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 60, 1, 6);
  setLogarithmicY(h_event_duration);
  TH2D* h_tres = new TH2D("h_tres", string(title + ";RUN;Time residuals [ns]").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 100, -50, 150);
  TH2D* h_likelihood = new TH2D ("h_likelihood", string(title + ";RUN;Likelihood").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 100, -1000, 1000);
  TH2D* h_beta0 = new TH2D("h_beta0", string(title + ";RUN;beta0").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 20, 0, 1);
  TH2D* h_energy = new TH2D ("h_energy", string(title + ";RUN;Energy [GeV];").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 65, 0, 9);
  setLogarithmicY(h_energy);
  TH2D* h_zenith = new TH2D("h_zenith", string(title + ";RUN;cosZenith").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 20, -1, 1);
  TH2D* h_azimuth = new TH2D("h_azimuth", string(title + ";RUN;cosAzimuth").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 20, -1, 1);
  TH2D* h_radial_position = new TH2D ("h_radial_position", string(title + ";RUN;Radial Position [m]").c_str(),  xbins, (double) min_run-0.5, (double) max_run+0.5, 60, 0, 4.2);
  setLogarithmicY(h_radial_position);
  TH2D* h_z_position = new TH2D ("h_z_position", string(title + ";RUN;Z Position [m]").c_str(), xbins, (double) min_run-0.5, (double) max_run+0.5, 50, 0, 3.2);
  setLogarithmicY(h_z_position);

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

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

    {
      JSingleFileScanner<run_info> in(*i);
      if(in.hasNext()){ 
        runinfo = *in.next();
        if(in.hasNext()){
          FATAL("MULTIPLE RUN INFO");
        }
      } else {
        FATAL("NO RUN INFO");
      } 
    }

    TFile in(i->c_str(), "read");

    cout << i->c_str() << endl;
    
    int run = runinfo.run_id;
    
    TH1D* h1d_prd  = new TH1D();
    TH1D* h1d_thb  = new TH1D();
    TH1D* h1d_th   = new TH1D();
    TH1D* h1d_th3D = new TH1D();
    TH1D* h1d_sh   = new TH1D();
    TH1D* h1d_tosh = new TH1D();
    TH1D* h1d_noo  = new TH1D();
    TH1D* h1d_pdth = new TH1D();
    TH1D* h1d_pdsh = new TH1D();
    TH1D* h1d_ed   = new TH1D();
    TH1D* h1d_tr   = new TH1D();
    TH1D* h1d_lk   = new TH1D();
    TH1D* h1d_b0   = new TH1D();
    TH1D* h1d_en   = new TH1D();
    TH1D* h1d_zn   = new TH1D();
    TH1D* h1d_az   = new TH1D();
    TH1D* h1d_rp   = new TH1D();
    TH1D* h1d_zp   = new TH1D();

    in.GetDirectory("Detector")->GetObject("h_pmt_rate_distribution", h1d_prd);
    fill2DHistogram(run, h1d_prd, h_pmt_rate_distribution, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    TDirectory* dir = in.GetDirectory("JDAQEvent");

    dir->GetObject("h_Trigger_bit_hit", h1d_thb);
    fill2DHistogram(run, h1d_thb, h_Trigger_bit_hit, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    dir->GetObject("h_Triggered_hits", h1d_th);
    fill2DHistogram(run, h1d_th, h_Triggered_hits, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    dir->GetObject("h_Triggered_hits_3dmuon", h1d_th3D);
    fill2DHistogram(run, h1d_th3D, h_Triggered_hits_3dmuon, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    dir->GetObject("h_Snapshot_hits", h1d_sh);
    fill2DHistogram(run, h1d_sh, h_Snapshot_hits, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    dir->GetObject("h_Triggered_over_Snapshot_hits", h1d_tosh);
    fill2DHistogram(run, h1d_tosh, h_Triggered_over_Snapshot_hits, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    dir->GetObject("h_Number_of_overlays", h1d_noo);
    fill2DHistogram(run, h1d_noo, h_Number_of_overlays, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    dir->GetObject("h_pmt_distribution_triggered_hits", h1d_pdth);
    fill2DHistogram(run, h1d_pdth, h_pmt_distribution_triggered_hits, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    dir->GetObject("h_pmt_distribution_snapshot_hits", h1d_pdsh);
    fill2DHistogram(run, h1d_pdsh, h_pmt_distribution_snapshot_hits, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    dir->GetObject("h_event_duration", h1d_ed);
    fill2DHistogram(run, h1d_ed, h_event_duration, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    TDirectory* dir_r = in.GetDirectory("Reco");

    dir_r->GetObject("h_tres", h1d_tr);
    fill2DHistogram(run, h1d_tr, h_tres, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    dir_r->GetObject("h_likelihood", h1d_lk);
    fill2DHistogram(run, h1d_lk, h_likelihood, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    dir_r->GetObject("h_beta0", h1d_b0);
    fill2DHistogram(run, h1d_b0, h_beta0, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    dir_r->GetObject("h_energy", h1d_en);
    fill2DHistogram(run, h1d_en, h_energy, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    dir_r->GetObject("h_zenith", h1d_zn);
    fill2DHistogram(run, h1d_zn, h_zenith, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    dir_r->GetObject("h_azimuth", h1d_az);
    fill2DHistogram(run, h1d_az, h_azimuth, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    dir_r->GetObject("h_radial_position", h1d_rp);
    fill2DHistogram(run, h1d_rp, h_radial_position, weight_mc, runinfo.daq_livetime, total_mc_livetime);

    dir_r->GetObject("h_z_position", h1d_zp);
    fill2DHistogram(run, h1d_zp, h_z_position, weight_mc, runinfo.daq_livetime, total_mc_livetime);
  }

  outputFile.open();

  outputFile.put(*h_event_duration);
  outputFile.put(*h_pmt_distribution_snapshot_hits);
  outputFile.put(*h_pmt_distribution_triggered_hits);
  outputFile.put(*h_Number_of_overlays);
  outputFile.put(*h_Triggered_over_Snapshot_hits);
  outputFile.put(*h_Snapshot_hits);
  outputFile.put(*h_Triggered_hits_3dmuon);
  outputFile.put(*h_Triggered_hits);
  outputFile.put(*h_Trigger_bit_hit);
  outputFile.put(*h_tres);
  outputFile.put(*h_likelihood);
  outputFile.put(*h_beta0);
  outputFile.put(*h_energy);
  outputFile.put(*h_zenith);
  outputFile.put(*h_azimuth);
  outputFile.put(*h_radial_position);
  outputFile.put(*h_z_position);

  JMultipleFileScanner<run_info> io(inputFile);

  io >> outputFile;

  outputFile.close();
}