|
|
|
#include <iostream>#include <string.h>#include <fstream>#include <cmath>#include <TFile.h>#include <TH1F.h>#include <TMath.h>#include <TTree.h>#include "Utilities.h"Go to the source code of this file.
Functions | |
| int | main (int argc, char **argv) |
| int main | ( | int | argc, |
| char ** | argv | ||
| ) |
Definition at line 13 of file ReadTree.cc.
References alignmentValidation::c1, CommandLine::check(), gather_cfg::cout, getBin(), CommandLine::getValue(), CommandLine::getVector(), i, EcalCondDB::inf, j, mergeVDriftHistosByStation::name, NETA, CommandLine::parse(), and CommandLine::print().
{
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();
}
1.7.3