File indexing completed on 2024-10-08 23:09:48
0001 #include <iostream>
0002 #include <cstring>
0003 #include <fstream>
0004 #include <cmath>
0005 #include <TFile.h>
0006 #include <TH1F.h>
0007 #include <TMath.h>
0008 #include <TTree.h>
0009 #include "Utilities.h"
0010
0011 using namespace std;
0012
0013 int main(int argc, char **argv) {
0014 CommandLine c1;
0015 c1.parse(argc, argv);
0016
0017 std::string TreeFilename = c1.getValue<string>("TreeFilename");
0018 std::string HistoFilename = c1.getValue<string>("HistoFilename");
0019 bool UseRatioForResponse = c1.getValue<bool>("UseRatioForResponse", true);
0020 bool IsPFJet = c1.getValue<bool>("IsPFJet", false);
0021 int NJETS_MAX = c1.getValue<int>("NJETS_MAX", 2);
0022 double JetPT_MIN = c1.getValue<double>("JetPT_MIN", 1.0);
0023 double DR_MIN = c1.getValue<double>("DR_MIN", 0.25);
0024 std::vector<double> pt_vec = c1.getVector<double>("RefPtBoundaries");
0025 std::vector<double> eta_vec = c1.getVector<double>("EtaBoundaries");
0026 if (!c1.check())
0027 return 0;
0028 c1.print();
0029
0030 const int MAX_NETA = 83;
0031 const int MAX_NPT = 30;
0032 TH1F *hPtRef[MAX_NPT][MAX_NETA];
0033 TH1F *hPtJet[MAX_NPT][MAX_NETA];
0034 TH1F *hPtRefBarrel[MAX_NPT];
0035 TH1F *hPtJetBarrel[MAX_NPT];
0036 TH1F *hResponseBarrel[MAX_NPT];
0037 TH1F *hResponse[MAX_NPT][MAX_NETA];
0038 int NPT = pt_vec.size() - 1;
0039 int NETA = eta_vec.size() - 1;
0040 char name[1024];
0041 bool cut_ptmin, cut_dR, cut_njets;
0042 float ptJet, ptGen, etaJet, emfJet, chfJet, etaGen, phiJet, phiGen, dR, resp;
0043 int rank;
0044 int i, j, ind_pt, ind_eta, responseBins(2200);
0045 double responseLow(-1000.), responseHigh(100.);
0046 unsigned int entry;
0047 if (UseRatioForResponse) {
0048 responseBins = 200;
0049 responseLow = 0.;
0050 responseHigh = 2;
0051 }
0052 std::vector<std::string> HistoNamesList;
0053 TFile *inf = new TFile(TreeFilename.c_str(), "R");
0054 if (inf->IsZombie())
0055 return 0;
0056 TFile *outf = new TFile(HistoFilename.c_str(), "RECREATE");
0057
0058 TTree *tr = (TTree *)inf->Get("mcTruthTree");
0059 TBranch *br_ptJet = (TBranch *)tr->GetBranch("ptJet");
0060 br_ptJet->SetAddress(&ptJet);
0061 TBranch *br_ptGen = (TBranch *)tr->GetBranch("ptGen");
0062 br_ptGen->SetAddress(&ptGen);
0063 TBranch *br_emfJet, *br_chfJet;
0064 if (!IsPFJet) {
0065 br_emfJet = (TBranch *)tr->GetBranch("emfJet");
0066 br_emfJet->SetAddress(&emfJet);
0067 } else {
0068 br_chfJet = (TBranch *)tr->GetBranch("chfJet");
0069 br_chfJet->SetAddress(&chfJet);
0070 }
0071 TBranch *br_etaJet = (TBranch *)tr->GetBranch("etaJet");
0072 br_etaJet->SetAddress(&etaJet);
0073 TBranch *br_etaGen = (TBranch *)tr->GetBranch("etaGen");
0074 br_etaGen->SetAddress(&etaGen);
0075 TBranch *br_phiJet = (TBranch *)tr->GetBranch("phiJet");
0076 br_phiJet->SetAddress(&phiJet);
0077 TBranch *br_phiGen = (TBranch *)tr->GetBranch("phiGen");
0078 br_phiGen->SetAddress(&phiGen);
0079 TBranch *br_dR = (TBranch *)tr->GetBranch("dR");
0080 br_dR->SetAddress(&dR);
0081 TBranch *br_rank = (TBranch *)tr->GetBranch("rank");
0082 br_rank->SetAddress(&rank);
0083 std::cout << "Total entries:" << tr->GetEntries() << std::endl;
0084 if (NETA > 1) {
0085 for (i = 0; i < NPT; i++) {
0086 sprintf(name, "ptRef_RefPt%d_Barrel", i);
0087 hPtRefBarrel[i] = new TH1F(name, name, 2000, 0, 2000);
0088 sprintf(name, "ptJet_RefPt%d_Barrel", i);
0089 hPtJetBarrel[i] = new TH1F(name, name, 2000, 0, 2000);
0090 sprintf(name, "Response_RefPt%d_Barrel", i);
0091 hResponseBarrel[i] = new TH1F(name, name, responseBins, responseLow, responseHigh);
0092 for (j = 0; j < NETA; j++) {
0093 sprintf(name, "ptRef_RefPt%d_Eta%d", i, j);
0094 hPtRef[i][j] = new TH1F(name, name, 2000, 0, 2000);
0095 sprintf(name, "ptCalo_RefPt%d_Eta%d", i, j);
0096 hPtJet[i][j] = new TH1F(name, name, 2000, 0, 2000);
0097 sprintf(name, "Response_RefPt%d_Eta%d", i, j);
0098 hResponse[i][j] = new TH1F(name, name, responseBins, responseLow, responseHigh);
0099 }
0100 }
0101 } else {
0102 for (i = 0; i < NPT; i++) {
0103 sprintf(name, "ptRef_RefPt%d", i);
0104 hPtRefBarrel[i] = new TH1F(name, name, 2000, 0, 2000);
0105 sprintf(name, "ptCalo_RefPt%d", i);
0106 hPtJetBarrel[i] = new TH1F(name, name, 2000, 0, 2000);
0107 sprintf(name, "Response_RefPt%d", i);
0108 hResponseBarrel[i] = new TH1F(name, name, responseBins, responseLow, responseHigh);
0109 }
0110 }
0111 std::cout << "Histograms booked" << std::endl;
0112 for (entry = 0; entry < tr->GetEntries(); entry++) {
0113 if (entry % (tr->GetEntries() / 10) == 0)
0114 std::cout << "Entries: " << entry << std::endl;
0115 tr->GetEntry(entry);
0116 cut_ptmin = (ptJet > JetPT_MIN);
0117 cut_dR = (dR < DR_MIN);
0118 cut_njets = (rank < NJETS_MAX);
0119 if (cut_ptmin && cut_dR && cut_njets) {
0120 ind_pt = getBin(ptGen, pt_vec);
0121 ind_eta = getBin(etaJet, eta_vec);
0122 if (UseRatioForResponse && ptGen > 0)
0123 resp = ptJet / ptGen;
0124 else
0125 resp = ptJet - ptGen;
0126 if (NETA > 1) {
0127 if (fabs(etaJet) < 1.3 && ind_pt >= 0) {
0128 hPtRefBarrel[ind_pt]->Fill(ptGen);
0129 hPtJetBarrel[ind_pt]->Fill(ptJet);
0130 hResponseBarrel[ind_pt]->Fill(resp);
0131 }
0132 if (ind_pt >= 0 && ind_eta >= 0 && NETA > 1) {
0133 hPtRef[ind_pt][ind_eta]->Fill(ptGen);
0134 hPtJet[ind_pt][ind_eta]->Fill(ptJet);
0135 hResponse[ind_pt][ind_eta]->Fill(resp);
0136 }
0137 }
0138 if (NETA == 1)
0139 if (ind_pt >= 0) {
0140 hPtRefBarrel[ind_pt]->Fill(ptGen);
0141 hPtJetBarrel[ind_pt]->Fill(ptJet);
0142 hResponseBarrel[ind_pt]->Fill(resp);
0143 }
0144 }
0145 }
0146 outf->Write();
0147 outf->Close();
0148 return 0;
0149 }
