MCJet/bin/L2Correction.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 <TGraph.h>
0010 #include <TGraphErrors.h>
0011 #include <TMath.h>
0012 #include <TKey.h>
0013 #include <TList.h>
0014 #include "Utilities.h"
0015
0016 using namespace std;
0017
0018 int main(int argc, char **argv) {
0019   CommandLine c1;
0020   c1.parse(argc, argv);
0021
0022   std::string HistoFilename = c1.getValue<string>("HistoFilename");
0023   std::string FitterFilename = c1.getValue<string>("FitterFilename");
0024   std::string L3ResponseTxtFilename = c1.getValue<string>("L3ResponseTxtFilename");
0025   std::string L3CorrectionTxtFilename = c1.getValue<string>("L3CorrectionTxtFilename");
0026   std::string L3OutputROOTFilename = c1.getValue<string>("L3OutputROOTFilename");
0027   std::string L2CorrectionTxtFilename = c1.getValue<string>("L2CorrectionTxtFilename");
0028   std::string L2OutputROOTFilename = c1.getValue<string>("L2OutputROOTFilename");
0029   std::vector<double> pt_vec = c1.getVector<double>("RefPtBoundaries");
0030   std::vector<double> eta_vec = c1.getVector<double>("EtaBoundaries");
0031   if (!c1.check())
0032     return 0;
0033   c1.print();
0034   /////////////////////////////////////////////////////////////////////////
0035   const int MAX_NETA = 83;
0036   const int MAX_NREFPTBINS = 30;
0037   const int NCaloPtValues = 18;
0038   int NRefPtBins = pt_vec.size() - 1;
0039   int NETA = eta_vec.size() - 1;
0040   if (NETA < 2)
0041     return (0);
0042   int i, auxi, etabin, ptbin;
0043   char name[100], func[1024];
0044   double cor, e_cor;
0045   double MinCaloPt[MAX_NETA], MaxCaloPt[MAX_NETA];
0046   double correction_x[MAX_NREFPTBINS], correction_y[MAX_NREFPTBINS], correction_ex[MAX_NREFPTBINS],
0047       correction_ey[MAX_NREFPTBINS];
0048   double cor_rel[MAX_NREFPTBINS], ref_pt, calo_pt, control_pt;
0049   double L3_resp[5], L2_cor[10];
0050   TKey *key;
0051   TFile *inf;
0052   TFile *outf;
0053   TH1F *hcorrection[MAX_NREFPTBINS];
0054   TH1F *h;
0055   TF1 *Correction[MAX_NETA];
0056   TF1 *L2Correction[MAX_NETA];
0057   TF1 *L3Response;
0058   std::ifstream L3ResponseFile;
0059   std::ofstream L2File;
0060   TGraph *g_L2Correction[MAX_NETA];
0061   TGraphErrors *g_EtaCorrection[MAX_NETA];
0062   double aux_CaloPt[NCaloPtValues] = {
0063       10, 15, 20, 30, 40, 50, 75, 100, 150, 200, 300, 400, 500, 750, 1000, 1500, 2000, 3000};
0064   std::vector<std::string> HistoNamesList;
0065   inf = new TFile(FitterFilename.c_str(), "r");
0066   if (inf->IsZombie())
0067     return (0);
0068   TIter next(inf->GetListOfKeys());
0069   while ((key = (TKey *)next()))
0070     HistoNamesList.push_back(key->GetName());
0071   //////////////////// Reading the L3 Response ////////////////////////////////
0072   L3ResponseFile.open(L3ResponseTxtFilename.c_str());
0073   L3ResponseFile >> L3_resp[0] >> L3_resp[1] >> L3_resp[2] >> L3_resp[3] >> L3_resp[4];
0074   L3Response = new TF1("L3Response", "[0]-[1]/(pow(log10(x),[2])+[3])+[4]/x", pt_vec[0], pt_vec[NRefPtBins - 1]);
0075   for (i = 0; i < 5; i++)
0076     L3Response->SetParameter(i, L3_resp[i]);
0077   L3ResponseFile.close();
0078   /////////////////// Calculating the L2 Correction //////////////////////////////
0079   for (i = 0; i < NRefPtBins; i++) {
0080     sprintf(name, "Correction_vs_Eta_RefPt%d", i);
0081     if (!HistoExists(HistoNamesList, name))
0082       return (0);
0083     hcorrection[i] = (TH1F *)inf->Get(name);
0084   }
0085   for (etabin = 0; etabin < NETA; etabin++) {
0086     sprintf(name, "MeanCaloPt_Eta%d", etabin);
0087     if (!HistoExists(HistoNamesList, name))
0088       return (0);
0089     h = (TH1F *)inf->Get(name);
0090     ///////////// Absolute correction calculation for every eta bin  //////////
0091     auxi = 0;
0092     for (ptbin = 0; ptbin < NRefPtBins; ptbin++) {
0093       cor = hcorrection[ptbin]->GetBinContent(etabin + 1);
0094       e_cor = hcorrection[ptbin]->GetBinError(etabin + 1);
0095       if (cor > 0 && e_cor > 0.0001 && e_cor < 0.3) {
0096         correction_x[auxi] = h->GetBinContent(ptbin + 1);  //average CaloPt for the eta bin
0097         correction_ex[auxi] = 0.;
0098         correction_y[auxi] = cor;
0099         correction_ey[auxi] = e_cor;
0100         auxi++;
0101       }
0102     }
0103     sprintf(name, "Correction%d", etabin);
0104     if (auxi > 1) {
0105       MaxCaloPt[etabin] = correction_x[auxi - 1];
0106       MinCaloPt[etabin] = correction_x[1];
0107       if (auxi > 10)
0108         sprintf(func, "[0]+[1]/(pow(log10(x),[2])+[3])");
0109       else
0110         sprintf(func, "[0]+[1]*log10(x)+[2]*pow(log10(x),2)");
0111       Correction[etabin] = new TF1(name, func, MinCaloPt[etabin], MaxCaloPt[etabin]);
0112     } else {
0113       std::cout << name << ": not enough points" << std::endl;
0114       sprintf(func, "[0]");
0115       correction_x[0] = 10;
0116       correction_x[1] = 100;
0117       correction_ex[0] = 0.;
0118       correction_ex[1] = 0.;
0119       correction_y[0] = 1.;
0120       correction_y[1] = 1.;
0121       correction_ey[0] = 0.;
0122       correction_ey[1] = 0.;
0123       auxi = 2;
0124       Correction[etabin] = new TF1(name, func, 10, 100);
0125       MaxCaloPt[etabin] = 0;
0126       MinCaloPt[etabin] = 0;
0127     }
0128     g_EtaCorrection[etabin] = new TGraphErrors(auxi, correction_x, correction_y, correction_ex, correction_ey);
0129     Correction[etabin]->SetParameter(0, 0.);
0130     Correction[etabin]->SetParameter(1, 0.);
0131     Correction[etabin]->SetParameter(2, 0.);
0132     Correction[etabin]->SetParameter(3, 0.);
0133     g_EtaCorrection[etabin]->Fit(name, "RQ");
0134     std::cout << name << " fitted....." << std::endl;
0135     ///////////// L2 Relative correction calculation for every eta bin /////////
0136     auxi = 0;
0137     for (ptbin = 0; ptbin < NCaloPtValues; ptbin++) {
0138       calo_pt = aux_CaloPt[ptbin];
0139       if (calo_pt >= MinCaloPt[etabin] && calo_pt <= MaxCaloPt[etabin]) {
0140         ref_pt = calo_pt * Correction[etabin]->Eval(calo_pt);
0141         control_pt = ref_pt * (L3Response->Eval(ref_pt));
0142         cor_rel[auxi] = control_pt / calo_pt;
0143         correction_x[auxi] = calo_pt;
0144         auxi++;
0145       }
0146     }
0147     sprintf(name, "L2Correction%d", etabin);
0148     if (auxi >= 2) {
0149       sprintf(func, "[0]+[1]*log10(x)+[2]*pow(log10(x),2)");
0150       if (auxi == 2)
0151         sprintf(func, "[0]+[1]*log10(x)");
0152     } else {
0153       sprintf(func, "[0]");
0154       correction_x[0] = 10;
0155       correction_x[1] = 100;
0156       cor_rel[0] = 1.;
0157       cor_rel[1] = 1.;
0158       auxi = 2;
0159     }
0160     g_L2Correction[etabin] = new TGraph(auxi, correction_x, cor_rel);
0161     L2Correction[etabin] = new TF1(name, func, correction_x[0], correction_x[auxi - 1]);
0162     L2Correction[etabin]->SetParameter(0, 0.);
0163     L2Correction[etabin]->SetParameter(1, 0.);
0164     L2Correction[etabin]->SetParameter(2, 0.);
0165     L2Correction[etabin]->SetParameter(3, 0.);
0166     L2Correction[etabin]->SetParameter(4, 0.);
0167     L2Correction[etabin]->SetParameter(5, 0.);
0168     g_L2Correction[etabin]->Fit(name, "RQ");
0169     std::cout << name << " fitted....." << std::endl;
0170   }  //end of eta bin loop
0171   //////////////////////// Writing //////////////////////////////
0172   L2File.open(L2CorrectionTxtFilename.c_str());
0173   L2File.setf(ios::right);
0174   for (etabin = 0; etabin < NETA; etabin++) {
0175     for (i = 0; i < 6; i++)
0176       L2_cor[i] = L2Correction[etabin]->GetParameter(i);
0177     L2File << setw(11) << eta_vec[etabin] << setw(11) << eta_vec[etabin + 1] << setw(11) << (int)8 << setw(12)
0178            << MinCaloPt[etabin] << setw(12) << MaxCaloPt[etabin] << setw(13) << L2_cor[0] << setw(13) << L2_cor[1]
0179            << setw(13) << L2_cor[2] << setw(13) << L2_cor[3] << setw(13) << L2_cor[4] << setw(13) << L2_cor[5] << "\n";
0180   }
0181   L2File.close();
0182   std::cout << L2CorrectionTxtFilename << " written...." << std::endl;
0183   outf = new TFile(L2OutputROOTFilename.c_str(), "RECREATE");
0184   for (etabin = 0; etabin < NETA; etabin++) {
0185     sprintf(name, "Correction_EtaBin%d", etabin);
0186     g_EtaCorrection[etabin]->Write(name);
0187     sprintf(name, "L2Correction_EtaBin%d", etabin);
0188     g_L2Correction[etabin]->Write(name);
0189   }
0190   outf->Close();
0191 }