ResponseFitter.cc

Go to the documentation of this file.

#include <iostream>
#include <cstring>
#include <fstream>
#include <cmath>
#include <TFile.h>
#include <TH1F.h>
#include <TMath.h>
#include <TKey.h>
#include <TList.h>
#include "Utilities.h"

using namespace std;

int main(int argc, char **argv) {
  CommandLine c1;
  c1.parse(argc, argv);

  std::string HistoFilename = c1.getValue<string>("HistoFilename");
  std::string FitterFilename = c1.getValue<string>("FitterFilename");
  std::string L3ResponseTxtFilename = c1.getValue<string>("L3ResponseTxtFilename");
  std::string L3CorrectionTxtFilename = c1.getValue<string>("L3CorrectionTxtFilename");
  std::string L3OutputROOTFilename = c1.getValue<string>("L3OutputROOTFilename");
  std::string L2CorrectionTxtFilename = c1.getValue<string>("L2CorrectionTxtFilename");
  std::string L2OutputROOTFilename = c1.getValue<string>("L2OutputROOTFilename");
  bool UseRatioForResponse = c1.getValue<bool>("UseRatioForResponse");
  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_NREFPTBINS = 30;
  int NRefPtBins = pt_vec.size() - 1;
  int NETA = eta_vec.size() - 1;
  char name[100];
  int j, etabin;
  double e, mR, eR, sR, seR, mRBarrel, eRBarrel, sRBarrel, seRBarrel, r, c;
  double mCaloPt, eCaloPt, sCaloPt, mRefPt, eRefPt, sRefPt;
  double mRefPtEtaBin, eRefPtEtaBin, sRefPtEtaBin, mCaloPtEtaBin, eCaloPtEtaBin, sCaloPtEtaBin;
  double EtaBoundaries[MAX_NETA], RefPtBoundaries[MAX_NREFPTBINS];
  std::vector<std::string> HistoNamesList;
  for (j = 0; j <= NRefPtBins; j++)
    RefPtBoundaries[j] = pt_vec[j];
  for (j = 0; j <= NETA; j++)
    EtaBoundaries[j] = eta_vec[j];
  TFile *inf;                                   //Input file containing the response histograms.
  TFile *outf;                                  //Output file containing the fitter results.
  TH1F *BarrelResponse;                         //Histogram with the barrel response in RefPt bins.
  TH1F *BarrelCorrection;                       //Histogram with the barrel correction in RefPt bins.
  TH1F *MeanRefPt_Barrel;                       //Histogram with the barrel average RefPt in RefPt bins.
  TH1F *MeanCaloPt_Barrel;                      //Histogram with the barrel average CaloPt in RefPt bins.
  TH1F *MeanRefPt_EtaBin[MAX_NETA];             //Histograms with the average RefPt in Eta & RefPt bins.
  TH1F *MeanCaloPt_EtaBin[MAX_NETA];            //Histograms with the average CaloPt in Eta & RefPt bins.
  TH1F *ResponseVsEta_RefPt[MAX_NREFPTBINS];    //Histograms with the average response vs Eta in RefPt bins.
  TH1F *CorrectionVsEta_RefPt[MAX_NREFPTBINS];  //Histograms with the average correction vs Eta in RefPt bins.
  TH1F *h;                                      //Auxilary histogram;
  TKey *key;
  inf = new TFile(HistoFilename.c_str(), "r");
  if (inf->IsZombie())
    return (0);
  TIter next(inf->GetListOfKeys());
  while ((key = (TKey *)next()))
    HistoNamesList.push_back(key->GetName());
  outf = new TFile(FitterFilename.c_str(), "RECREATE");
  TDirectory *dir_Response =
      (TDirectory *)outf->mkdir("FittedHistograms");  //Directory in output file to store the fitted histograms.
  BarrelResponse = new TH1F("Response", "Response", NRefPtBins, RefPtBoundaries);
  BarrelCorrection = new TH1F("Correction", "Correction", NRefPtBins, RefPtBoundaries);
  MeanRefPt_Barrel = new TH1F("MeanRefPt", "MeanRefPt", NRefPtBins, RefPtBoundaries);
  MeanCaloPt_Barrel = new TH1F("MeanCaloPt", "MeanCaloPt", NRefPtBins, RefPtBoundaries);
  if (NETA > 1)  //multiple eta bins: used for L2+L3 correction calculation
  {
    std::cout << "************* Fitting Response Histograms in multiple Eta bins. ************" << std::endl;
    for (etabin = 0; etabin < NETA; etabin++) {
      sprintf(name, "MeanRefPt_Eta%d", etabin);
      MeanRefPt_EtaBin[etabin] = new TH1F(name, name, NRefPtBins, RefPtBoundaries);
      sprintf(name, "MeanCaloPt_Eta%d", etabin);
      MeanCaloPt_EtaBin[etabin] = new TH1F(name, name, NRefPtBins, RefPtBoundaries);
    }
    for (j = 0; j < NRefPtBins; j++)  //loop over RefPt bins
    {
      std::cout << "RefJetPt Bin: [" << RefPtBoundaries[j] << "," << RefPtBoundaries[j + 1] << "] GeV" << std::endl;
      sprintf(name, "ptRef_RefPt%d_Barrel", j);
      if (!HistoExists(HistoNamesList, name))
        return (0);
      h = (TH1F *)inf->Get(name);
      GetMEAN(h, mRefPt, eRefPt, sRefPt);
      sprintf(name, "ptCalo_RefPt%d_Barrel", j);
      if (!HistoExists(HistoNamesList, name))
        return (0);
      h = (TH1F *)inf->Get(name);
      GetMEAN(h, mCaloPt, eCaloPt, sCaloPt);
      sprintf(name, "Response_RefPt%d_Barrel", j);
      if (!HistoExists(HistoNamesList, name))
        return (0);
      h = (TH1F *)inf->Get(name);
      GetMPV(name, h, dir_Response, mRBarrel, eRBarrel, sRBarrel, seRBarrel);
      MeanRefPt_Barrel->SetBinContent(j + 1, mRefPt);
      MeanRefPt_Barrel->SetBinError(j + 1, eRefPt);
      MeanCaloPt_Barrel->SetBinContent(j + 1, mCaloPt);
      MeanCaloPt_Barrel->SetBinError(j + 1, eCaloPt);
      CalculateResponse(UseRatioForResponse, mRefPt, eRefPt, mRBarrel, eRBarrel, r, e);
      BarrelResponse->SetBinContent(j + 1, r);
      BarrelResponse->SetBinError(j + 1, e);
      CalculateCorrection(UseRatioForResponse, mRefPt, eRefPt, mRBarrel, eRBarrel, c, e);
      BarrelCorrection->SetBinContent(j + 1, c);
      BarrelCorrection->SetBinError(j + 1, e);
      sprintf(name, "Response_vs_Eta_RefPt%d", j);
      ResponseVsEta_RefPt[j] = new TH1F(name, name, NETA, EtaBoundaries);
      sprintf(name, "Correction_vs_Eta_RefPt%d", j);
      CorrectionVsEta_RefPt[j] = new TH1F(name, name, NETA, EtaBoundaries);
      for (etabin = 0; etabin < NETA; etabin++)  //loop over eta bins
      {
        sprintf(name, "Response_RefPt%d_Eta%d", j, etabin);
        if (!HistoExists(HistoNamesList, name))
          return (0);
        h = (TH1F *)inf->Get(name);
        GetMPV(name, h, dir_Response, mR, eR, sR, seR);
        sprintf(name, "ptRef_RefPt%d_Eta%d", j, etabin);
        if (!HistoExists(HistoNamesList, name))
          return (0);
        h = (TH1F *)inf->Get(name);
        GetMEAN(h, mRefPtEtaBin, eRefPtEtaBin, sRefPtEtaBin);
        sprintf(name, "ptCalo_RefPt%d_Eta%d", j, etabin);
        if (!HistoExists(HistoNamesList, name))
          return (0);
        h = (TH1F *)inf->Get(name);
        GetMEAN(h, mCaloPtEtaBin, eCaloPtEtaBin, sCaloPtEtaBin);
        MeanRefPt_EtaBin[etabin]->SetBinContent(j + 1, mRefPtEtaBin);
        MeanRefPt_EtaBin[etabin]->SetBinError(j + 1, eRefPtEtaBin);
        MeanCaloPt_EtaBin[etabin]->SetBinContent(j + 1, mCaloPtEtaBin);
        MeanCaloPt_EtaBin[etabin]->SetBinError(j + 1, eCaloPtEtaBin);
        CalculateResponse(UseRatioForResponse, mRefPtEtaBin, eRefPtEtaBin, mR, eR, r, e);
        ResponseVsEta_RefPt[j]->SetBinContent(etabin + 1, r);
        ResponseVsEta_RefPt[j]->SetBinError(etabin + 1, e);
        CalculateCorrection(UseRatioForResponse, mRefPtEtaBin, eRefPtEtaBin, mR, eR, c, e);
        CorrectionVsEta_RefPt[j]->SetBinContent(etabin + 1, c);
        CorrectionVsEta_RefPt[j]->SetBinError(etabin + 1, e);
      }   //end of EtaBin loop
    }     // end of Pt loop
  } else  //single eta bin: used for L3 correction calculation
  {
    std::cout << "************* Fitting Response Histograms in single eta bin. ************" << std::endl;
    for (j = 0; j < NRefPtBins; j++)  //loop over Pt bins
    {
      std::cout << "RefJetPt Bin: [" << RefPtBoundaries[j] << "," << RefPtBoundaries[j + 1] << "] GeV" << std::endl;
      sprintf(name, "ptRef_RefPt%d", j);
      if (!HistoExists(HistoNamesList, name))
        return (0);
      h = (TH1F *)inf->Get(name);
      GetMEAN(h, mRefPt, eRefPt, sRefPt);
      sprintf(name, "ptCalo_RefPt%d", j);
      if (!HistoExists(HistoNamesList, name))
        return (0);
      h = (TH1F *)inf->Get(name);
      GetMEAN(h, mCaloPt, eCaloPt, sCaloPt);
      sprintf(name, "Response_RefPt%d", j);
      if (!HistoExists(HistoNamesList, name))
        return (0);
      h = (TH1F *)inf->Get(name);
      GetMPV(name, h, dir_Response, mRBarrel, eRBarrel, sRBarrel, seRBarrel);
      MeanRefPt_Barrel->SetBinContent(j + 1, mRefPt);
      MeanRefPt_Barrel->SetBinError(j + 1, eRefPt);
      MeanCaloPt_Barrel->SetBinContent(j + 1, mCaloPt);
      MeanCaloPt_Barrel->SetBinError(j + 1, eCaloPt);
      CalculateResponse(UseRatioForResponse, mRefPt, eRefPt, mRBarrel, eRBarrel, r, e);
      BarrelResponse->SetBinContent(j + 1, r);
      BarrelResponse->SetBinError(j + 1, e);
      CalculateCorrection(UseRatioForResponse, mRefPt, eRefPt, mRBarrel, eRBarrel, c, e);
      BarrelCorrection->SetBinContent(j + 1, c);
      BarrelCorrection->SetBinError(j + 1, e);
    }  // end of Pt loop
  }
  outf->cd();
  MeanRefPt_Barrel->Write();
  MeanCaloPt_Barrel->Write();
  BarrelResponse->Write();
  BarrelCorrection->Write();
  if (NETA > 1) {
    for (etabin = 0; etabin < NETA; etabin++) {
      MeanRefPt_EtaBin[etabin]->Write();
      MeanCaloPt_EtaBin[etabin]->Write();
    }
    for (j = 0; j < NRefPtBins; j++) {
      ResponseVsEta_RefPt[j]->Write();
      CorrectionVsEta_RefPt[j]->Write();
    }
  }
  outf->Close();
}