CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Functions
ReadTree.cc File Reference
#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)
 

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 13 of file ReadTree.cc.

References alignmentValidation::c1, CommandLine::check(), gather_cfg::cout, PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi::dR, HLT_25ns14e33_v1_cff::etaJet, getBin(), CommandLine::getValue(), CommandLine::getVector(), i, EcalCondDB::inf, j, mergeVDriftHistosByStation::name, NETA, CommandLine::parse(), CommandLine::print(), and AlCaHLTBitMon_QueryRunRegistry::string.

14 {
16  c1.parse(argc,argv);
17 
18  std::string TreeFilename = c1.getValue<string> ("TreeFilename" );
19  std::string HistoFilename = c1.getValue<string> ("HistoFilename" );
20  bool UseRatioForResponse = c1.getValue<bool> ("UseRatioForResponse", true);
21  bool IsPFJet = c1.getValue<bool> ("IsPFJet", false);
22  int NJETS_MAX = c1.getValue<int> ("NJETS_MAX", 2);
23  double JetPT_MIN = c1.getValue<double> ("JetPT_MIN", 1.0);
24  double DR_MIN = c1.getValue<double> ("DR_MIN", 0.25);
25  std::vector<double> pt_vec = c1.getVector<double>("RefPtBoundaries" );
26  std::vector<double> eta_vec = c1.getVector<double>("EtaBoundaries" );
27  if (!c1.check()) return 0;
28  c1.print();
30  const int MAX_NETA = 83;
31  const int MAX_NPT = 30;
32  TH1F *hPtRef[MAX_NPT][MAX_NETA];
33  TH1F *hPtJet[MAX_NPT][MAX_NETA];
34  TH1F *hPtRefBarrel[MAX_NPT];
35  TH1F *hPtJetBarrel[MAX_NPT];
36  TH1F *hResponseBarrel[MAX_NPT];
37  TH1F *hResponse[MAX_NPT][MAX_NETA];
38  int NPT = pt_vec.size()-1;
39  int NETA = eta_vec.size()-1;
40  char name[1024];
41  bool cut_ptmin,cut_dR,cut_njets;
42  float ptJet,ptGen,etaJet,emfJet,chfJet,etaGen,phiJet,phiGen,dR,resp;
43  int rank;
44  int i,j,ind_pt,ind_eta,responseBins(2200);
45  double responseLow(-1000.),responseHigh(100.);
46  unsigned int entry;
47  if (UseRatioForResponse)
48  {
49  responseBins = 200;
50  responseLow = 0.;
51  responseHigh = 2;
52  }
53  std::vector<std::string> HistoNamesList;
54  TFile *inf = new TFile(TreeFilename.c_str(),"R");
55  if (inf->IsZombie()) return 0;
56  TFile *outf = new TFile(HistoFilename.c_str(),"RECREATE");
57 
58  TTree *tr = (TTree*)inf->Get("mcTruthTree");
59  TBranch *br_ptJet = (TBranch*)tr->GetBranch("ptJet");
60  br_ptJet->SetAddress(&ptJet);
61  TBranch *br_ptGen = (TBranch*)tr->GetBranch("ptGen");
62  br_ptGen->SetAddress(&ptGen);
63  TBranch *br_emfJet,*br_chfJet;
64  if (!IsPFJet)
65  {
66  br_emfJet = (TBranch*)tr->GetBranch("emfJet");
67  br_emfJet->SetAddress(&emfJet);
68  }
69  else
70  {
71  br_chfJet = (TBranch*)tr->GetBranch("chfJet");
72  br_chfJet->SetAddress(&chfJet);
73  }
74  TBranch *br_etaJet = (TBranch*)tr->GetBranch("etaJet");
75  br_etaJet->SetAddress(&etaJet);
76  TBranch *br_etaGen = (TBranch*)tr->GetBranch("etaGen");
77  br_etaGen->SetAddress(&etaGen);
78  TBranch *br_phiJet = (TBranch*)tr->GetBranch("phiJet");
79  br_phiJet->SetAddress(&phiJet);
80  TBranch *br_phiGen = (TBranch*)tr->GetBranch("phiGen");
81  br_phiGen->SetAddress(&phiGen);
82  TBranch *br_dR = (TBranch*)tr->GetBranch("dR");
83  br_dR->SetAddress(&dR);
84  TBranch *br_rank = (TBranch*)tr->GetBranch("rank");
85  br_rank->SetAddress(&rank);
86  std::cout<<"Total entries:"<<tr->GetEntries()<<std::endl;
87  if (NETA>1)
88  {
89  for (i=0;i<NPT;i++)
90  {
91  sprintf(name,"ptRef_RefPt%d_Barrel",i);
92  hPtRefBarrel[i] = new TH1F(name,name,2000,0,2000);
93  sprintf(name,"ptJet_RefPt%d_Barrel",i);
94  hPtJetBarrel[i] = new TH1F(name,name,2000,0,2000);
95  sprintf(name,"Response_RefPt%d_Barrel",i);
96  hResponseBarrel[i] = new TH1F(name,name,responseBins,responseLow,responseHigh);
97  for (j=0;j<NETA;j++)
98  {
99  sprintf(name,"ptRef_RefPt%d_Eta%d",i,j);
100  hPtRef[i][j] = new TH1F(name,name,2000,0,2000);
101  sprintf(name,"ptCalo_RefPt%d_Eta%d",i,j);
102  hPtJet[i][j] = new TH1F(name,name,2000,0,2000);
103  sprintf(name,"Response_RefPt%d_Eta%d",i,j);
104  hResponse[i][j] = new TH1F(name,name,responseBins,responseLow,responseHigh);
105  }
106  }
107  }
108  else
109  {
110  for (i=0;i<NPT;i++)
111  {
112  sprintf(name,"ptRef_RefPt%d",i);
113  hPtRefBarrel[i] = new TH1F(name,name,2000,0,2000);
114  sprintf(name,"ptCalo_RefPt%d",i);
115  hPtJetBarrel[i] = new TH1F(name,name,2000,0,2000);
116  sprintf(name,"Response_RefPt%d",i);
117  hResponseBarrel[i] = new TH1F(name,name,responseBins,responseLow,responseHigh);
118  }
119  }
120  std::cout<<"Histograms booked"<<std::endl;
121  for(entry=0;entry<tr->GetEntries();entry++)
122  {
123  if (entry % (tr->GetEntries()/10) == 0)
124  std::cout<<"Entries: "<<entry<<std::endl;
125  tr->GetEntry(entry);
126  cut_ptmin = (ptJet>JetPT_MIN);
127  cut_dR = (dR<DR_MIN);
128  cut_njets = (rank<NJETS_MAX);
129  if (cut_ptmin && cut_dR && cut_njets)
130  {
131  ind_pt = getBin(ptGen,pt_vec);
132  ind_eta = getBin(etaJet,eta_vec);
133  resp = 0.;
134  if (UseRatioForResponse && ptGen>0)
135  resp = ptJet/ptGen;
136  else
137  resp = ptJet - ptGen;
138  if (NETA>1)
139  {
140  if (fabs(etaJet)<1.3 && ind_pt>=0)
141  {
142  hPtRefBarrel[ind_pt]->Fill(ptGen);
143  hPtJetBarrel[ind_pt]->Fill(ptJet);
144  hResponseBarrel[ind_pt]->Fill(resp);
145  }
146  if (ind_pt>=0 && ind_eta>=0 && NETA>1)
147  {
148  hPtRef[ind_pt][ind_eta]->Fill(ptGen);
149  hPtJet[ind_pt][ind_eta]->Fill(ptJet);
150  hResponse[ind_pt][ind_eta]->Fill(resp);
151  }
152  }
153  if (NETA==1)
154  if (ind_pt>=0)
155  {
156  hPtRefBarrel[ind_pt]->Fill(ptGen);
157  hPtJetBarrel[ind_pt]->Fill(ptJet);
158  hResponseBarrel[ind_pt]->Fill(resp);
159  }
160  }
161  }
162  outf->Write();
163  outf->Close();
164  return 0;
165 }
int i
Definition: DBlmapReader.cc:9
const int NETA
string inf
Definition: EcalCondDB.py:94
int j
Definition: DBlmapReader.cc:9
T getValue(const std::string &name)
Definition: Utilities.h:74
int getBin(double x, std::vector< double > boundaries)
Definition: Utilities.h:512
tuple argc
Definition: dir2webdir.py:38
tuple cout
Definition: gather_cfg.py:121
bool parse(int argc, char **argv)
Definition: Utilities.h:199