MCJet/bin/L3Correction.cc

0001 #include <iostream>
0002 #include <iomanip>
0003 #include <cstring>
0004 #include <fstream>
0005 #include <cmath>
0006 #include <TFile.h>
0007 #include <TH1F.h>
0008 #include <TF1.h>
0009 #include <TGraphErrors.h>
0010 #include <TMath.h>
0011 #include <TKey.h>
0012 #include <TList.h>
0013 #include <TVirtualFitter.h>
0014 #include <TMatrixD.h>
0015 #include "Utilities.h"
0016
0017 using namespace std;
0018
0019 int main(int argc, char **argv) {
0020   CommandLine c1;
0021   c1.parse(argc, argv);
0022
0023   std::string HistoFilename = c1.getValue<string>("HistoFilename");
0024   std::string FitterFilename = c1.getValue<string>("FitterFilename");
0025   std::string L3ResponseTxtFilename = c1.getValue<string>("L3ResponseTxtFilename");
0026   std::string L3CorrectionTxtFilename = c1.getValue<string>("L3CorrectionTxtFilename");
0027   std::string L3OutputROOTFilename = c1.getValue<string>("L3OutputROOTFilename");
0028   std::string L2CorrectionTxtFilename = c1.getValue<string>("L2CorrectionTxtFilename");
0029   std::string L2OutputROOTFilename = c1.getValue<string>("L2OutputROOTFilename");
0030   std::vector<double> pt_vec = c1.getVector<double>("RefPtBoundaries");
0031   std::vector<double> eta_vec = c1.getVector<double>("EtaBoundaries");
0032   if (!c1.check())
0033     return 0;
0034   c1.print();
0035   /////////////////////////////////////////////////////////////////////////
0036   const int MAX_NREFPTBINS = 30;
0037   int NRefPtBins = pt_vec.size() - 1;
0038   char name[100];
0039   int i, auxi_cor, auxi_resp;
0040   double Correction[MAX_NREFPTBINS] = {};
0041   double CorrectionError[MAX_NREFPTBINS] = {};
0042   double Response[MAX_NREFPTBINS] = {};
0043   double ResponseError[MAX_NREFPTBINS] = {};
0044   double xCaloPt[MAX_NREFPTBINS] = {};
0045   double exCaloPt[MAX_NREFPTBINS] = {};
0046   double xRefPt[MAX_NREFPTBINS] = {};
0047   double exRefPt[MAX_NREFPTBINS] = {};
0048   double MinCaloPt, MaxCaloPt, MinRefPt, MaxRefPt;
0049   double cor, e_cor, resp, e_resp;
0050   TFile *inf;
0051   TFile *outf;
0052   TH1F *hCor, *hResp, *hRef, *hCalo;
0053   TF1 *CorFit;
0054   TF1 *RespFit;
0055   TGraphErrors *g_Cor;
0056   TGraphErrors *g_Resp;
0057   TVirtualFitter *fitter;
0058   TMatrixD *COV_Cor, *COV_Resp;
0059   std::ofstream L3CorrectionFile, L3ResponseFile;
0060   L3CorrectionFile.open(L3CorrectionTxtFilename.c_str());
0061   L3ResponseFile.open(L3ResponseTxtFilename.c_str());
0062   TKey *key;
0063   std::vector<std::string> HistoNamesList;
0064   inf = new TFile(FitterFilename.c_str(), "r");
0065   if (inf->IsZombie())
0066     return (0);
0067   TIter next(inf->GetListOfKeys());
0068   while ((key = (TKey *)next()))
0069     HistoNamesList.push_back(key->GetName());
0070   ///////////// CaloPt /////////////////////////////////////
0071   sprintf(name, "MeanCaloPt");
0072   if (!HistoExists(HistoNamesList, name))
0073     return (0);
0074   hCalo = (TH1F *)inf->Get(name);
0075   ///////////// RefPt /////////////////////////////////////
0076   sprintf(name, "MeanRefPt");
0077   if (!HistoExists(HistoNamesList, name))
0078     return (0);
0079   hRef = (TH1F *)inf->Get(name);
0080   ///////////////////////// Correction //////////////////////////////////////////////
0081   sprintf(name, "Correction");
0082   if (!HistoExists(HistoNamesList, name))
0083     return (0);
0084   hCor = (TH1F *)inf->Get(name);
0085   sprintf(name, "Response");
0086   if (!HistoExists(HistoNamesList, name))
0087     return (0);
0088   hResp = (TH1F *)inf->Get(name);
0089   auxi_cor = 0;
0090   auxi_resp = 0;
0091   std::cout << "RefPt bin" << setw(12) << "<CaloPt>" << setw(12) << "Correction" << setw(12) << "Error" << setw(12)
0092             << "<RefPt>" << setw(12) << "Response" << setw(12) << "Error" << std::endl;
0093   for (i = 0; i < NRefPtBins; i++) {
0094     cor = hCor->GetBinContent(i + 1);
0095     e_cor = hCor->GetBinError(i + 1);
0096     std::cout << "[" << pt_vec[i] << "," << pt_vec[i + 1] << "]" << setw(12) << hCalo->GetBinContent(i + 1) << setw(12)
0097               << cor << setw(12) << e_cor;
0098     resp = hResp->GetBinContent(i + 1);
0099     e_resp = hResp->GetBinError(i + 1);
0100     std::cout << setw(12) << hRef->GetBinContent(i + 1) << setw(12) << resp << setw(12) << e_resp << std::endl;
0101     if (cor > 0 && e_cor > 0.000001 && e_cor < 0.2) {
0102       Correction[auxi_cor] = cor;
0103       CorrectionError[auxi_cor] = e_cor;
0104       xCaloPt[auxi_cor] = hCalo->GetBinContent(i + 1);
0105       exCaloPt[auxi_cor] = hCalo->GetBinError(i + 1);
0106       auxi_cor++;
0107     }
0108     if (resp > 0 && e_resp > 0.000001 && e_resp < 0.2) {
0109       Response[auxi_resp] = resp;
0110       ResponseError[auxi_resp] = e_resp;
0111       xRefPt[auxi_resp] = hRef->GetBinContent(i + 1);
0112       exRefPt[auxi_resp] = hRef->GetBinError(i + 1);
0113       auxi_resp++;
0114     }
0115   }
0116   g_Cor = new TGraphErrors(auxi_cor, xCaloPt, Correction, exCaloPt, CorrectionError);
0117   sprintf(name, "CorFit");
0118   CorFit = new TF1(name, "[0]+[1]/(pow(log10(x),[2])+[3])", xCaloPt[1], xCaloPt[auxi_cor - 1]);
0119   CorFit->SetParameter(0, 1.);
0120   CorFit->SetParameter(1, 7.);
0121   CorFit->SetParameter(2, 4.);
0122   CorFit->SetParameter(3, 4.);
0123   std::cout << "Fitting correction in the range: " << xCaloPt[1] << " " << xCaloPt[auxi_cor - 1] << std::endl;
0124   g_Cor->Fit(CorFit, "RQ");
0125   fitter = TVirtualFitter::GetFitter();
0126   COV_Cor = new TMatrixD(4, 4, fitter->GetCovarianceMatrix());
0127   g_Resp = new TGraphErrors(auxi_resp, xRefPt, Response, exRefPt, ResponseError);
0128   sprintf(name, "RespFit");
0129   RespFit = new TF1(name, "[0]-[1]/(pow(log10(x),[2])+[3])+[4]/x", xRefPt[0], xRefPt[auxi_resp - 1]);
0130   RespFit->SetParameter(0, 1.);
0131   RespFit->SetParameter(1, 1.);
0132   RespFit->SetParameter(2, 1.);
0133   RespFit->SetParameter(3, 1.);
0134   RespFit->SetParameter(4, 1.);
0135   std::cout << "Fitting response in the range: " << xRefPt[0] << " " << xRefPt[auxi_resp - 1] << std::endl;
0136   g_Resp->Fit(RespFit, "RQ");
0137   fitter = TVirtualFitter::GetFitter();
0138   COV_Resp = new TMatrixD(5, 5, fitter->GetCovarianceMatrix());
0139   /////////////////////////////////////////////////////////
0140   MinCaloPt = 4.;
0141   MaxCaloPt = 5000.;
0142   MinRefPt = 4.;
0143   MaxRefPt = 5000.;
0144   CorFit->SetRange(MinCaloPt, MaxCaloPt);
0145   RespFit->SetRange(MinRefPt, MaxRefPt);
0146   //////////////////////////// Writing /////////////////////////////////////////////
0147   L3CorrectionFile.setf(ios::left);
0148   L3CorrectionFile << setw(12) << -5.191 << setw(12) << 5.191 << setw(12) << (int)6 << setw(12) << MinCaloPt << setw(12)
0149                    << MaxCaloPt << setw(12) << CorFit->GetParameter(0) << setw(12) << CorFit->GetParameter(1)
0150                    << setw(12) << CorFit->GetParameter(2) << setw(12) << CorFit->GetParameter(3);
0151   L3CorrectionFile.close();
0152
0153   L3ResponseFile.setf(ios::left);
0154   L3ResponseFile << setw(12) << RespFit->GetParameter(0) << setw(12) << RespFit->GetParameter(1) << setw(12)
0155                  << RespFit->GetParameter(2) << setw(12) << RespFit->GetParameter(3) << setw(12)
0156                  << RespFit->GetParameter(4);
0157   L3ResponseFile.close();
0158
0159   outf = new TFile(L3OutputROOTFilename.c_str(), "RECREATE");
0160   g_Cor->Write("Correction_vs_CaloPt");
0161   COV_Cor->Write("CovMatrix_Correction");
0162   g_Resp->Write("Response_vs_RefPt");
0163   COV_Resp->Write("CovMatrix_Resp");
0164   outf->Close();
0165 }