JDepthDependence.cc

Go to the documentation of this file.

#include <iomanip>
#include <iostream>
#include <sstream>
#include <string>

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

#include "TApplication.h"
#include "TLegend.h"
#include "TGraphErrors.h"

#include "JDAQ/JDAQ.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "JLang/JString.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"

using namespace JPP;  

double getCount (TH1D* hptr, int muon_threshold) {
  double count = 0;
  for (int i = hptr->GetXaxis()->FindBin(muon_threshold); i <= hptr->GetNbinsX(); i++) {
    count += hptr->GetBinContent(i);
  }
  return count;
}

double getLiveTime (TH1D* hptr, const JModule& module) {
  string label = getModuleLabel(module);
  double livetime = hptr->GetBinContent(hptr->GetXaxis()->FindBin(label.c_str()));
  // trick for retrocompatibility with old JMM
  if (livetime == 0) {
    std::replace(label.begin(), label.end() - 1, '0', ' ');
    livetime = hptr->GetBinContent(hptr->GetXaxis()->FindBin(label.c_str()));
  }
  return livetime;
}

/**
 * \file
 *
 * Auxiliary program to analyse coincidence multiplicity;
 * \author mlincetto
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;  
  
  //-----------------------------------------------------------
  // parameter interface
  //-----------------------------------------------------------
  
  const int nFiles = 2;
  vector<string>     inputFile(nFiles,"");
  string             outputFile;
  string             detectorFile;
  int                line;
  int                debug;
  int                muonThreshold;
  int                minLiveTime_s;
  bool               dataLiveTime;
  bool               interactive;

  try { 

    JParser<> zap("Example program to plot the muon depth-dependence relationship from JMonitorMultiplicity output");
    
    zap['f'] = make_field(inputFile[0], "JMM data input file");
    zap['g'] = make_field(inputFile[1], "JMM MC input file") = "";
    zap['o'] = make_field(outputFile, "Name for the .pdf and .root output files")          = "depthdependence";
    zap['a'] = make_field(detectorFile);
    zap['d'] = make_field(debug)                = 1;
    zap['l'] = make_field(line)                 = 1;
    zap['M'] = make_field(muonThreshold, "Minimum multiplicity")        = 8;
    zap['T'] = make_field(minLiveTime_s, "Minimum DOM livetime [s] to be eligible for plotting")        = 3600;
    zap['I'] = make_field(interactive, "Launch TApplication for interactive display.");
    zap['L'] = make_field(dataLiveTime);

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

  cout.tie(&cerr);


  //----------------------
  // loading detector file
  //----------------------
  
  JDetector detector;

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

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

  const JModuleRouter router(detector);

  double utm_z = detector.getUTMZ();

  NOTICE("Detector base UTM z [m]: " << utm_z << endl);

  //-----------------------------------------------------------
  // load input files
  //-----------------------------------------------------------

  vector<TFile*> in;
  TFile* fptr;

  for (unsigned i = 0; i < inputFile.size() && inputFile[i] != ""; i++) {
    
    DEBUG("Loading input file " << inputFile[i] << endl);

    fptr  = TFile::Open(inputFile[i].c_str(), "exist");

    if (fptr == NULL || !fptr->IsOpen()) {
      FATAL("File: " << inputFile[i] << " not opened." << endl);
    }

    in.push_back(fptr);
  }
  
  //-----------------------------------------------------------
  // recover rates from histograms
  //-----------------------------------------------------------

  vector<TH1D*> LiveTime;

  for (unsigned i = 0; i < in.size(); i++) {
    DEBUG("Loading livetime histogram from " << inputFile[i] << endl);
    TH1D* hptr = (TH1D*)in[i]->Get("LiveTime");
    
    if (hptr == NULL) {
      FATAL("Missing live time histogram.");
    }

    LiveTime.push_back(hptr);
  }

  vector<double> depth;
  vector<double> data_rate;
  vector<double> data_err;
  vector<double> mc_rate;
  vector<double> mc_err;

  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
    
    if (module->getString() != line) { continue; }

    STATUS("Loading data histogram for from input file: " << getModuleLabel(*module) << "\t ID = " << module->getID()  << endl);

    TH1D* data_histogram = (TH1D*)in[0]->Get(TString(getModuleLabel(*module)) + "_P");

    if (data_histogram != NULL) {

      double data_count = getCount(data_histogram, muonThreshold);
      double data_livetime = getLiveTime(LiveTime[0], *module);

      if (data_livetime > minLiveTime_s) {
        
        double rate = data_count / data_livetime;
        double err = sqrt(data_count) / data_livetime;

        NOTICE(getModuleLabel(*module) << ": z_rel [m] = " << module->getZ() << ", count [#] = " << data_count << ", livetime [s] = " << data_livetime << ", rate [Hz] = "  << rate << endl);

        data_rate.push_back(rate);
        data_err.push_back(err);
        depth.push_back((-utm_z) - module->getZ());
        
        if (in.size() > 1) {

            LiveTime[1]  = (TH1D*)in[1]->Get("LiveTime");
            TH1D*  mc_histogram = (TH1D*)in[1]->Get(TString(getModuleLabel(*module)) + "_P");
            double mc_count     = getCount(mc_histogram, muonThreshold);
            double mc_livetime  = (dataLiveTime ? data_livetime : getLiveTime(LiveTime[1], *module));

            double rate = mc_count / mc_livetime;
            double err  = sqrt(mc_count) / mc_livetime;
            
            NOTICE("\t MC livetime [s] = " << mc_livetime << ", MC rate [Hz] = " << rate << endl);
            
            mc_rate.push_back(rate);
            mc_err.push_back(err);
          }
      }
    }
    
  }

  //-----------------------------------------------------------
  // time to graph!
  //-----------------------------------------------------------

  TApplication* theApp = NULL;

  if (interactive) {
    theApp = new TApplication("Plotter", &argc, argv);
  }

  TCanvas *c1 = new TCanvas("c1","KM3NeT preliminary");

  TGraph* gr_data = new TGraphErrors(data_rate.size(), &depth[0], &data_rate[0], 0, &data_err[0]);
  gr_data->SetMarkerColor(kRed);
  gr_data->SetLineColor(kRed);
  TGraph* gr_mc;
  TMultiGraph* mgr = new TMultiGraph();
  mgr->Add(gr_data);

  TLegend* leg = new TLegend(0.7, 0.7, 0.9, 0.9);
  leg->AddEntry(gr_data, "Data", "lp");


  if (in.size() > 1) {
    DEBUG("Loading MC plots" << endl);
    gr_mc = new TGraphErrors(mc_rate.size(), &depth[0], &mc_rate[0], 0, &mc_err[0]); 
    mgr->Add(gr_mc);
    DEBUG("MC plot loaded" << endl);
    leg->AddEntry(gr_mc, "MC", "lp");

  }

  mgr->SetTitle(TString("KM3NeT Preliminary; Depth [m]; Inclusive ") + Form("%d", muonThreshold) + TString("-fold coincidence rate [Hz]"));
  
  mgr->Draw("AL");
  mgr->GetYaxis()->SetTitleOffset(1.5);

  leg->SetFillStyle(0);
  leg->Draw();
  
  c1->Print(TString(outputFile.c_str()) + ".pdf");

  if (interactive) {
    theApp->Run();
  }

  //-----------------------------------------------------------
  // writing output
  //-----------------------------------------------------------

  // writing of ROOT file should be included
  
  return 0;

}