MCJet/bin/ResponseFitter.cc

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 <TKey.h>
0009 #include <TList.h>
0010 #include "Utilities.h"
0011
0012 using namespace std;
0013
0014 int main(int argc, char **argv) {
0015   CommandLine c1;
0016   c1.parse(argc, argv);
0017
0018   std::string HistoFilename = c1.getValue<string>("HistoFilename");
0019   std::string FitterFilename = c1.getValue<string>("FitterFilename");
0020   std::string L3ResponseTxtFilename = c1.getValue<string>("L3ResponseTxtFilename");
0021   std::string L3CorrectionTxtFilename = c1.getValue<string>("L3CorrectionTxtFilename");
0022   std::string L3OutputROOTFilename = c1.getValue<string>("L3OutputROOTFilename");
0023   std::string L2CorrectionTxtFilename = c1.getValue<string>("L2CorrectionTxtFilename");
0024   std::string L2OutputROOTFilename = c1.getValue<string>("L2OutputROOTFilename");
0025   bool UseRatioForResponse = c1.getValue<bool>("UseRatioForResponse");
0026   std::vector<double> pt_vec = c1.getVector<double>("RefPtBoundaries");
0027   std::vector<double> eta_vec = c1.getVector<double>("EtaBoundaries");
0028   if (!c1.check())
0029     return 0;
0030   c1.print();
0031   /////////////////////////////////////////////////////////////////////////
0032   const int MAX_NETA = 83;
0033   const int MAX_NREFPTBINS = 30;
0034   int NRefPtBins = pt_vec.size() - 1;
0035   int NETA = eta_vec.size() - 1;
0036   char name[100];
0037   int j, etabin;
0038   double e, mR, eR, sR, seR, mRBarrel, eRBarrel, sRBarrel, seRBarrel, r, c;
0039   double mCaloPt, eCaloPt, sCaloPt, mRefPt, eRefPt, sRefPt;
0040   double mRefPtEtaBin, eRefPtEtaBin, sRefPtEtaBin, mCaloPtEtaBin, eCaloPtEtaBin, sCaloPtEtaBin;
0041   double EtaBoundaries[MAX_NETA], RefPtBoundaries[MAX_NREFPTBINS];
0042   std::vector<std::string> HistoNamesList;
0043   for (j = 0; j <= NRefPtBins; j++)
0044     RefPtBoundaries[j] = pt_vec[j];
0045   for (j = 0; j <= NETA; j++)
0046     EtaBoundaries[j] = eta_vec[j];
0047   TFile *inf;                                   //Input file containing the response histograms.
0048   TFile *outf;                                  //Output file containing the fitter results.
0049   TH1F *BarrelResponse;                         //Histogram with the barrel response in RefPt bins.
0050   TH1F *BarrelCorrection;                       //Histogram with the barrel correction in RefPt bins.
0051   TH1F *MeanRefPt_Barrel;                       //Histogram with the barrel average RefPt in RefPt bins.
0052   TH1F *MeanCaloPt_Barrel;                      //Histogram with the barrel average CaloPt in RefPt bins.
0053   TH1F *MeanRefPt_EtaBin[MAX_NETA];             //Histograms with the average RefPt in Eta & RefPt bins.
0054   TH1F *MeanCaloPt_EtaBin[MAX_NETA];            //Histograms with the average CaloPt in Eta & RefPt bins.
0055   TH1F *ResponseVsEta_RefPt[MAX_NREFPTBINS];    //Histograms with the average response vs Eta in RefPt bins.
0056   TH1F *CorrectionVsEta_RefPt[MAX_NREFPTBINS];  //Histograms with the average correction vs Eta in RefPt bins.
0057   TH1F *h;                                      //Auxilary histogram;
0058   TKey *key;
0059   //////////////////////////////////////////////////////////////////////////
0060   inf = new TFile(HistoFilename.c_str(), "r");
0061   if (inf->IsZombie())
0062     return (0);
0063   TIter next(inf->GetListOfKeys());
0064   while ((key = (TKey *)next()))
0065     HistoNamesList.push_back(key->GetName());
0066   outf = new TFile(FitterFilename.c_str(), "RECREATE");
0067   TDirectory *dir_Response =
0068       (TDirectory *)outf->mkdir("FittedHistograms");  //Directory in output file to store the fitted histograms.
0069   BarrelResponse = new TH1F("Response", "Response", NRefPtBins, RefPtBoundaries);
0070   BarrelCorrection = new TH1F("Correction", "Correction", NRefPtBins, RefPtBoundaries);
0071   MeanRefPt_Barrel = new TH1F("MeanRefPt", "MeanRefPt", NRefPtBins, RefPtBoundaries);
0072   MeanCaloPt_Barrel = new TH1F("MeanCaloPt", "MeanCaloPt", NRefPtBins, RefPtBoundaries);
0073   if (NETA > 1)  //multiple eta bins: used for L2+L3 correction calculation
0074   {
0075     std::cout << "************* Fitting Response Histograms in multiple Eta bins. ************" << std::endl;
0076     for (etabin = 0; etabin < NETA; etabin++) {
0077       sprintf(name, "MeanRefPt_Eta%d", etabin);
0078       MeanRefPt_EtaBin[etabin] = new TH1F(name, name, NRefPtBins, RefPtBoundaries);
0079       sprintf(name, "MeanCaloPt_Eta%d", etabin);
0080       MeanCaloPt_EtaBin[etabin] = new TH1F(name, name, NRefPtBins, RefPtBoundaries);
0081     }
0082     for (j = 0; j < NRefPtBins; j++)  //loop over RefPt bins
0083     {
0084       std::cout << "RefJetPt Bin: [" << RefPtBoundaries[j] << "," << RefPtBoundaries[j + 1] << "] GeV" << std::endl;
0085       sprintf(name, "ptRef_RefPt%d_Barrel", j);
0086       if (!HistoExists(HistoNamesList, name))
0087         return (0);
0088       h = (TH1F *)inf->Get(name);
0089       GetMEAN(h, mRefPt, eRefPt, sRefPt);
0090       sprintf(name, "ptCalo_RefPt%d_Barrel", j);
0091       if (!HistoExists(HistoNamesList, name))
0092         return (0);
0093       h = (TH1F *)inf->Get(name);
0094       GetMEAN(h, mCaloPt, eCaloPt, sCaloPt);
0095       sprintf(name, "Response_RefPt%d_Barrel", j);
0096       if (!HistoExists(HistoNamesList, name))
0097         return (0);
0098       h = (TH1F *)inf->Get(name);
0099       GetMPV(name, h, dir_Response, mRBarrel, eRBarrel, sRBarrel, seRBarrel);
0100       ///////////////// RefPt in barrel ///////////////////////////////////
0101       MeanRefPt_Barrel->SetBinContent(j + 1, mRefPt);
0102       MeanRefPt_Barrel->SetBinError(j + 1, eRefPt);
0103       ///////////////// CaloPt in barrel ///////////////////////////////////
0104       MeanCaloPt_Barrel->SetBinContent(j + 1, mCaloPt);
0105       MeanCaloPt_Barrel->SetBinError(j + 1, eCaloPt);
0106       ////////////////// Absolute response in barrel //////////////////////////
0107       CalculateResponse(UseRatioForResponse, mRefPt, eRefPt, mRBarrel, eRBarrel, r, e);
0108       BarrelResponse->SetBinContent(j + 1, r);
0109       BarrelResponse->SetBinError(j + 1, e);
0110       ////////////////// Absolute correction in barrel //////////////////////////
0111       CalculateCorrection(UseRatioForResponse, mRefPt, eRefPt, mRBarrel, eRBarrel, c, e);
0112       BarrelCorrection->SetBinContent(j + 1, c);
0113       BarrelCorrection->SetBinError(j + 1, e);
0114       ////////////////// Eta bins /////////////////////////////////////
0115       sprintf(name, "Response_vs_Eta_RefPt%d", j);
0116       ResponseVsEta_RefPt[j] = new TH1F(name, name, NETA, EtaBoundaries);
0117       sprintf(name, "Correction_vs_Eta_RefPt%d", j);
0118       CorrectionVsEta_RefPt[j] = new TH1F(name, name, NETA, EtaBoundaries);
0119       for (etabin = 0; etabin < NETA; etabin++)  //loop over eta bins
0120       {
0121         ///////////////////////////////////////////////////////////////
0122         sprintf(name, "Response_RefPt%d_Eta%d", j, etabin);
0123         if (!HistoExists(HistoNamesList, name))
0124           return (0);
0125         h = (TH1F *)inf->Get(name);
0126         GetMPV(name, h, dir_Response, mR, eR, sR, seR);
0127         sprintf(name, "ptRef_RefPt%d_Eta%d", j, etabin);
0128         if (!HistoExists(HistoNamesList, name))
0129           return (0);
0130         h = (TH1F *)inf->Get(name);
0131         GetMEAN(h, mRefPtEtaBin, eRefPtEtaBin, sRefPtEtaBin);
0132         sprintf(name, "ptCalo_RefPt%d_Eta%d", j, etabin);
0133         if (!HistoExists(HistoNamesList, name))
0134           return (0);
0135         h = (TH1F *)inf->Get(name);
0136         GetMEAN(h, mCaloPtEtaBin, eCaloPtEtaBin, sCaloPtEtaBin);
0137         ///////////////// RefPt in etabin ///////////////////////////////////
0138         MeanRefPt_EtaBin[etabin]->SetBinContent(j + 1, mRefPtEtaBin);
0139         MeanRefPt_EtaBin[etabin]->SetBinError(j + 1, eRefPtEtaBin);
0140         ///////////////// CaloPt in etabin ///////////////////////////////////
0141         MeanCaloPt_EtaBin[etabin]->SetBinContent(j + 1, mCaloPtEtaBin);
0142         MeanCaloPt_EtaBin[etabin]->SetBinError(j + 1, eCaloPtEtaBin);
0143         ////////////////// Absolute response in etabin //////////////////////////
0144         CalculateResponse(UseRatioForResponse, mRefPtEtaBin, eRefPtEtaBin, mR, eR, r, e);
0145         ResponseVsEta_RefPt[j]->SetBinContent(etabin + 1, r);
0146         ResponseVsEta_RefPt[j]->SetBinError(etabin + 1, e);
0147         ////////////////// Absolute correction in etabin //////////////////////////
0148         CalculateCorrection(UseRatioForResponse, mRefPtEtaBin, eRefPtEtaBin, mR, eR, c, e);
0149         CorrectionVsEta_RefPt[j]->SetBinContent(etabin + 1, c);
0150         CorrectionVsEta_RefPt[j]->SetBinError(etabin + 1, e);
0151       }  //end of EtaBin loop
0152     }  // end of Pt loop
0153   } else  //single eta bin: used for L3 correction calculation
0154   {
0155     std::cout << "************* Fitting Response Histograms in single eta bin. ************" << std::endl;
0156     for (j = 0; j < NRefPtBins; j++)  //loop over Pt bins
0157     {
0158       std::cout << "RefJetPt Bin: [" << RefPtBoundaries[j] << "," << RefPtBoundaries[j + 1] << "] GeV" << std::endl;
0159       sprintf(name, "ptRef_RefPt%d", j);
0160       if (!HistoExists(HistoNamesList, name))
0161         return (0);
0162       h = (TH1F *)inf->Get(name);
0163       GetMEAN(h, mRefPt, eRefPt, sRefPt);
0164       sprintf(name, "ptCalo_RefPt%d", j);
0165       if (!HistoExists(HistoNamesList, name))
0166         return (0);
0167       h = (TH1F *)inf->Get(name);
0168       GetMEAN(h, mCaloPt, eCaloPt, sCaloPt);
0169       sprintf(name, "Response_RefPt%d", j);
0170       if (!HistoExists(HistoNamesList, name))
0171         return (0);
0172       h = (TH1F *)inf->Get(name);
0173       GetMPV(name, h, dir_Response, mRBarrel, eRBarrel, sRBarrel, seRBarrel);
0174       ///////////////// RefPt in barrel ///////////////////////////////////
0175       MeanRefPt_Barrel->SetBinContent(j + 1, mRefPt);
0176       MeanRefPt_Barrel->SetBinError(j + 1, eRefPt);
0177       ///////////////// CaloPt in barrel ///////////////////////////////////
0178       MeanCaloPt_Barrel->SetBinContent(j + 1, mCaloPt);
0179       MeanCaloPt_Barrel->SetBinError(j + 1, eCaloPt);
0180       ////////////////// Absolute response in barrel //////////////////////////
0181       CalculateResponse(UseRatioForResponse, mRefPt, eRefPt, mRBarrel, eRBarrel, r, e);
0182       BarrelResponse->SetBinContent(j + 1, r);
0183       BarrelResponse->SetBinError(j + 1, e);
0184       ////////////////// Absolute correction in barrel //////////////////////////
0185       CalculateCorrection(UseRatioForResponse, mRefPt, eRefPt, mRBarrel, eRBarrel, c, e);
0186       BarrelCorrection->SetBinContent(j + 1, c);
0187       BarrelCorrection->SetBinError(j + 1, e);
0188     }  // end of Pt loop
0189   }
0190   ////////////////////// writing ///////////////////////////////
0191   outf->cd();
0192   MeanRefPt_Barrel->Write();
0193   MeanCaloPt_Barrel->Write();
0194   BarrelResponse->Write();
0195   BarrelCorrection->Write();
0196   if (NETA > 1) {
0197     for (etabin = 0; etabin < NETA; etabin++) {
0198       MeanRefPt_EtaBin[etabin]->Write();
0199       MeanCaloPt_EtaBin[etabin]->Write();
0200     }
0201     for (j = 0; j < NRefPtBins; j++) {
0202       ResponseVsEta_RefPt[j]->Write();
0203       CorrectionVsEta_RefPt[j]->Write();
0204     }
0205   }
0206   outf->Close();
0207 }