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/JMeta.hh"

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

#include "JGizmo/JGizmoToolkit.hh"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.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 norm=false){

  double normalisation = 1;
  if(norm) normalisation = h1->Integral();
  
  for(int i = 1; i <= h1->GetNbinsX(); ++i){

    h2->Fill(run, h1->GetXaxis()->GetBinCenter(i), h1->GetBinContent(i)/normalisation);       
  }
}

int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  vector<string>     inputFile;
  string             outputFile;
  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['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;
  
  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {

    TFile in(i->c_str(), "read");
    
    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();

    }

    in.Close();
  }

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

  TFile out(outputFile.c_str(), "recreate");
  
  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);
  
  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {

    DEBUG("Processing " << *i << endl) ;
    
    TFile in(i->c_str(), "read");

    cout << i->c_str() << endl;
    
    int run = 0;
    
    for (JRootFileReader<JDAQHeader> in(i->c_str()); in.hasNext(); ) {
      
      const JDAQHeader* p = in.next();

      run = p->getRunNumber();

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

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

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

    dir->GetObject("h_Trigger_bit_hit", h1d_thb);
    fill2DHistogram(run, h1d_thb, h_Trigger_bit_hit);

    dir->GetObject("h_Triggered_hits", h1d_th);
    fill2DHistogram(run, h1d_th, h_Triggered_hits);

    dir->GetObject("h_Triggered_hits_3dmuon", h1d_th3D);
    fill2DHistogram(run, h1d_th3D, h_Triggered_hits_3dmuon);

    dir->GetObject("h_Snapshot_hits", h1d_sh);
    fill2DHistogram(run, h1d_sh, h_Snapshot_hits);

    dir->GetObject("h_Triggered_over_Snapshot_hits", h1d_tosh);
    fill2DHistogram(run, h1d_tosh, h_Triggered_over_Snapshot_hits);

    dir->GetObject("h_Number_of_overlays", h1d_noo);
    fill2DHistogram(run, h1d_noo, h_Number_of_overlays);

    dir->GetObject("h_pmt_distribution_triggered_hits", h1d_pdth);
    fill2DHistogram(run, h1d_pdth, h_pmt_distribution_triggered_hits);

    dir->GetObject("h_pmt_distribution_snapshot_hits", h1d_pdsh);
    fill2DHistogram(run, h1d_pdsh, h_pmt_distribution_snapshot_hits);

    dir->GetObject("h_event_duration", h1d_ed);
    fill2DHistogram(run, h1d_ed, h_event_duration);

    in.Close();

  }

  out.Write();
  out.Close();
  
}