ReadTree.cc

Go to the documentation of this file.

#include <iostream>
#include <cstring>
#include <fstream>
#include <cmath>
#include <TFile.h>
#include <TH1F.h>
#include <TMath.h>
#include <TTree.h>
#include "Utilities.h"

using namespace std;

int main(int argc, char**argv)
{
  CommandLine c1;
  c1.parse(argc,argv);

  std::string TreeFilename      = c1.getValue<string> ("TreeFilename"             );
  std::string HistoFilename     = c1.getValue<string> ("HistoFilename"            );
  bool UseRatioForResponse = c1.getValue<bool>   ("UseRatioForResponse", true);
  bool IsPFJet             = c1.getValue<bool>   ("IsPFJet",            false);
  int NJETS_MAX            = c1.getValue<int>    ("NJETS_MAX",              2);
  double JetPT_MIN         = c1.getValue<double> ("JetPT_MIN",            1.0); 
  double DR_MIN            = c1.getValue<double> ("DR_MIN",              0.25);
  std::vector<double> pt_vec    = c1.getVector<double>("RefPtBoundaries"          );
  std::vector<double> eta_vec   = c1.getVector<double>("EtaBoundaries"            );
  if (!c1.check()) return 0; 
  c1.print();
  const int MAX_NETA = 83; 
  const int MAX_NPT  = 30; 
  TH1F *hPtRef[MAX_NPT][MAX_NETA];
  TH1F *hPtJet[MAX_NPT][MAX_NETA];
  TH1F *hPtRefBarrel[MAX_NPT];
  TH1F *hPtJetBarrel[MAX_NPT];
  TH1F *hResponseBarrel[MAX_NPT];
  TH1F *hResponse[MAX_NPT][MAX_NETA];
  int NPT  = pt_vec.size()-1;
  int NETA = eta_vec.size()-1;
  char name[1024];
  bool cut_ptmin,cut_dR,cut_njets;
  float ptJet,ptGen,etaJet,emfJet,chfJet,etaGen,phiJet,phiGen,dR,resp;
  int rank;
  int i,j,ind_pt,ind_eta,responseBins(2200);
  double responseLow(-1000.),responseHigh(100.);
  unsigned int entry;
  if (UseRatioForResponse)
    {
      responseBins = 200;
      responseLow  = 0.;
      responseHigh = 2;
    } 
  std::vector<std::string> HistoNamesList; 
  TFile *inf = new TFile(TreeFilename.c_str(),"R");
  if (inf->IsZombie()) return 0;
  TFile *outf = new TFile(HistoFilename.c_str(),"RECREATE");
  
  TTree *tr = (TTree*)inf->Get("mcTruthTree");
  TBranch *br_ptJet = (TBranch*)tr->GetBranch("ptJet");
  br_ptJet->SetAddress(&ptJet);
  TBranch *br_ptGen = (TBranch*)tr->GetBranch("ptGen");
  br_ptGen->SetAddress(&ptGen);
  TBranch *br_emfJet,*br_chfJet;
  if (!IsPFJet)
    {
      br_emfJet = (TBranch*)tr->GetBranch("emfJet");
      br_emfJet->SetAddress(&emfJet);
    }
  else
    {
      br_chfJet = (TBranch*)tr->GetBranch("chfJet");
      br_chfJet->SetAddress(&chfJet);
    }
  TBranch *br_etaJet = (TBranch*)tr->GetBranch("etaJet");
  br_etaJet->SetAddress(&etaJet);
  TBranch *br_etaGen = (TBranch*)tr->GetBranch("etaGen");
  br_etaGen->SetAddress(&etaGen);
  TBranch *br_phiJet = (TBranch*)tr->GetBranch("phiJet");
  br_phiJet->SetAddress(&phiJet);
  TBranch *br_phiGen = (TBranch*)tr->GetBranch("phiGen");
  br_phiGen->SetAddress(&phiGen);
  TBranch *br_dR = (TBranch*)tr->GetBranch("dR");
  br_dR->SetAddress(&dR);
  TBranch *br_rank = (TBranch*)tr->GetBranch("rank");
  br_rank->SetAddress(&rank);
  std::cout<<"Total entries:"<<tr->GetEntries()<<std::endl;
  if (NETA>1)
    {
      for (i=0;i<NPT;i++)
       {
         sprintf(name,"ptRef_RefPt%d_Barrel",i);
         hPtRefBarrel[i] = new TH1F(name,name,2000,0,2000);
         sprintf(name,"ptJet_RefPt%d_Barrel",i);
         hPtJetBarrel[i] = new TH1F(name,name,2000,0,2000);
         sprintf(name,"Response_RefPt%d_Barrel",i);
         hResponseBarrel[i] = new TH1F(name,name,responseBins,responseLow,responseHigh);     
         for (j=0;j<NETA;j++)
           {
         sprintf(name,"ptRef_RefPt%d_Eta%d",i,j);
             hPtRef[i][j] = new TH1F(name,name,2000,0,2000);
             sprintf(name,"ptCalo_RefPt%d_Eta%d",i,j);
             hPtJet[i][j] = new TH1F(name,name,2000,0,2000);
         sprintf(name,"Response_RefPt%d_Eta%d",i,j);
             hResponse[i][j] = new TH1F(name,name,responseBins,responseLow,responseHigh);
           }
       }
    }
  else
    {
      for (i=0;i<NPT;i++)
       {
         sprintf(name,"ptRef_RefPt%d",i);
         hPtRefBarrel[i] = new TH1F(name,name,2000,0,2000);
         sprintf(name,"ptCalo_RefPt%d",i);
         hPtJetBarrel[i] = new TH1F(name,name,2000,0,2000);
         sprintf(name,"Response_RefPt%d",i);
         hResponseBarrel[i] = new TH1F(name,name,responseBins,responseLow,responseHigh);     
       }
    }
  std::cout<<"Histograms booked"<<std::endl;
  for(entry=0;entry<tr->GetEntries();entry++)
    {
      if (entry % (tr->GetEntries()/10) == 0)
        std::cout<<"Entries: "<<entry<<std::endl;
      tr->GetEntry(entry);
      cut_ptmin = (ptJet>JetPT_MIN);
      cut_dR    = (dR<DR_MIN);
      cut_njets = (rank<NJETS_MAX);
      if (cut_ptmin && cut_dR && cut_njets)
        {
      ind_pt = getBin(ptGen,pt_vec);
      ind_eta = getBin(etaJet,eta_vec);   
          resp = 0.;
          if (UseRatioForResponse && ptGen>0)
            resp = ptJet/ptGen;
          else
            resp = ptJet - ptGen; 
          if (NETA>1)
            { 
              if (fabs(etaJet)<1.3 && ind_pt>=0)
                {
                  hPtRefBarrel[ind_pt]->Fill(ptGen);
                  hPtJetBarrel[ind_pt]->Fill(ptJet);
                  hResponseBarrel[ind_pt]->Fill(resp);
                }  
              if (ind_pt>=0 && ind_eta>=0 && NETA>1)
                { 
              hPtRef[ind_pt][ind_eta]->Fill(ptGen);
                  hPtJet[ind_pt][ind_eta]->Fill(ptJet);
              hResponse[ind_pt][ind_eta]->Fill(resp);
                }
            }
          if (NETA==1)
            if (ind_pt>=0)
              {
                hPtRefBarrel[ind_pt]->Fill(ptGen);
                hPtJetBarrel[ind_pt]->Fill(ptJet);
                hResponseBarrel[ind_pt]->Fill(resp);
              }  
        }
    }
  outf->Write();
  outf->Close();
  return 0;
}