Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:21:35

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       resp = 0.;
0123       if (UseRatioForResponse && ptGen > 0)
0124         resp = ptJet / ptGen;
0125       else
0126         resp = ptJet - ptGen;
0127       if (NETA > 1) {
0128         if (fabs(etaJet) < 1.3 && ind_pt >= 0) {
0129           hPtRefBarrel[ind_pt]->Fill(ptGen);
0130           hPtJetBarrel[ind_pt]->Fill(ptJet);
0131           hResponseBarrel[ind_pt]->Fill(resp);
0132         }
0133         if (ind_pt >= 0 && ind_eta >= 0 && NETA > 1) {
0134           hPtRef[ind_pt][ind_eta]->Fill(ptGen);
0135           hPtJet[ind_pt][ind_eta]->Fill(ptJet);
0136           hResponse[ind_pt][ind_eta]->Fill(resp);
0137         }
0138       }
0139       if (NETA == 1)
0140         if (ind_pt >= 0) {
0141           hPtRefBarrel[ind_pt]->Fill(ptGen);
0142           hPtJetBarrel[ind_pt]->Fill(ptJet);
0143           hResponseBarrel[ind_pt]->Fill(resp);
0144         }
0145     }
0146   }
0147   outf->Write();
0148   outf->Close();
0149   return 0;
0150 }