CMS 3D CMS Logo
Functions
L3Correction.cc File Reference
#include <iostream>
#include <iomanip>
#include <cstring>
#include <fstream>
#include <cmath>
#include <TFile.h>
#include <TH1F.h>
#include <TF1.h>
#include <TGraphErrors.h>
#include <TMath.h>
#include <TKey.h>
#include <TList.h>
#include <TVirtualFitter.h>
#include <TMatrixD.h>
#include "Utilities.h"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 19 of file L3Correction.cc.

References dir2webdir::argc, GCPpyPlots::argv, alignmentValidation::c1, gather_cfg::cout, HistoExists(), mps_fire::i, dqmiodatasetharvest::inf, createfilelist::int, submitPVResolutionJobs::key, Skims_PA_cff::name, GetRecoTauVFromDQM_MC_cff::next, and AlCaHLTBitMon_QueryRunRegistry::string.

19  {
20  CommandLine c1;
21  c1.parse(argc, argv);
22 
23  std::string HistoFilename = c1.getValue<string>("HistoFilename");
24  std::string FitterFilename = c1.getValue<string>("FitterFilename");
25  std::string L3ResponseTxtFilename = c1.getValue<string>("L3ResponseTxtFilename");
26  std::string L3CorrectionTxtFilename = c1.getValue<string>("L3CorrectionTxtFilename");
27  std::string L3OutputROOTFilename = c1.getValue<string>("L3OutputROOTFilename");
28  std::string L2CorrectionTxtFilename = c1.getValue<string>("L2CorrectionTxtFilename");
29  std::string L2OutputROOTFilename = c1.getValue<string>("L2OutputROOTFilename");
30  std::vector<double> pt_vec = c1.getVector<double>("RefPtBoundaries");
31  std::vector<double> eta_vec = c1.getVector<double>("EtaBoundaries");
32  if (!c1.check())
33  return 0;
34  c1.print();
36  const int MAX_NREFPTBINS = 30;
37  int NRefPtBins = pt_vec.size() - 1;
38  char name[100];
39  int i, auxi_cor, auxi_resp;
40  double Correction[MAX_NREFPTBINS] = {};
41  double CorrectionError[MAX_NREFPTBINS] = {};
42  double Response[MAX_NREFPTBINS] = {};
43  double ResponseError[MAX_NREFPTBINS] = {};
44  double xCaloPt[MAX_NREFPTBINS] = {};
45  double exCaloPt[MAX_NREFPTBINS] = {};
46  double xRefPt[MAX_NREFPTBINS] = {};
47  double exRefPt[MAX_NREFPTBINS] = {};
48  double MinCaloPt, MaxCaloPt, MinRefPt, MaxRefPt;
49  double cor, e_cor, resp, e_resp;
50  TFile *inf;
51  TFile *outf;
52  TH1F *hCor, *hResp, *hRef, *hCalo;
53  TF1 *CorFit;
54  TF1 *RespFit;
55  TGraphErrors *g_Cor;
56  TGraphErrors *g_Resp;
57  TVirtualFitter *fitter;
58  TMatrixD *COV_Cor, *COV_Resp;
59  std::ofstream L3CorrectionFile, L3ResponseFile;
60  L3CorrectionFile.open(L3CorrectionTxtFilename.c_str());
61  L3ResponseFile.open(L3ResponseTxtFilename.c_str());
62  TKey *key;
63  std::vector<std::string> HistoNamesList;
64  inf = new TFile(FitterFilename.c_str(), "r");
65  if (inf->IsZombie())
66  return (0);
67  TIter next(inf->GetListOfKeys());
68  while ((key = (TKey *)next()))
69  HistoNamesList.push_back(key->GetName());
71  sprintf(name, "MeanCaloPt");
72  if (!HistoExists(HistoNamesList, name))
73  return (0);
74  hCalo = (TH1F *)inf->Get(name);
76  sprintf(name, "MeanRefPt");
77  if (!HistoExists(HistoNamesList, name))
78  return (0);
79  hRef = (TH1F *)inf->Get(name);
81  sprintf(name, "Correction");
82  if (!HistoExists(HistoNamesList, name))
83  return (0);
84  hCor = (TH1F *)inf->Get(name);
85  sprintf(name, "Response");
86  if (!HistoExists(HistoNamesList, name))
87  return (0);
88  hResp = (TH1F *)inf->Get(name);
89  auxi_cor = 0;
90  auxi_resp = 0;
91  std::cout << "RefPt bin" << setw(12) << "<CaloPt>" << setw(12) << "Correction" << setw(12) << "Error" << setw(12)
92  << "<RefPt>" << setw(12) << "Response" << setw(12) << "Error" << std::endl;
93  for (i = 0; i < NRefPtBins; i++) {
94  cor = hCor->GetBinContent(i + 1);
95  e_cor = hCor->GetBinError(i + 1);
96  std::cout << "[" << pt_vec[i] << "," << pt_vec[i + 1] << "]" << setw(12) << hCalo->GetBinContent(i + 1) << setw(12)
97  << cor << setw(12) << e_cor;
98  resp = hResp->GetBinContent(i + 1);
99  e_resp = hResp->GetBinError(i + 1);
100  std::cout << setw(12) << hRef->GetBinContent(i + 1) << setw(12) << resp << setw(12) << e_resp << std::endl;
101  if (cor > 0 && e_cor > 0.000001 && e_cor < 0.2) {
102  Correction[auxi_cor] = cor;
103  CorrectionError[auxi_cor] = e_cor;
104  xCaloPt[auxi_cor] = hCalo->GetBinContent(i + 1);
105  exCaloPt[auxi_cor] = hCalo->GetBinError(i + 1);
106  auxi_cor++;
107  }
108  if (resp > 0 && e_resp > 0.000001 && e_resp < 0.2) {
109  Response[auxi_resp] = resp;
110  ResponseError[auxi_resp] = e_resp;
111  xRefPt[auxi_resp] = hRef->GetBinContent(i + 1);
112  exRefPt[auxi_resp] = hRef->GetBinError(i + 1);
113  auxi_resp++;
114  }
115  }
116  g_Cor = new TGraphErrors(auxi_cor, xCaloPt, Correction, exCaloPt, CorrectionError);
117  sprintf(name, "CorFit");
118  CorFit = new TF1(name, "[0]+[1]/(pow(log10(x),[2])+[3])", xCaloPt[1], xCaloPt[auxi_cor - 1]);
119  CorFit->SetParameter(0, 1.);
120  CorFit->SetParameter(1, 7.);
121  CorFit->SetParameter(2, 4.);
122  CorFit->SetParameter(3, 4.);
123  std::cout << "Fitting correction in the range: " << xCaloPt[1] << " " << xCaloPt[auxi_cor - 1] << std::endl;
124  g_Cor->Fit(CorFit, "RQ");
125  fitter = TVirtualFitter::GetFitter();
126  COV_Cor = new TMatrixD(4, 4, fitter->GetCovarianceMatrix());
127  g_Resp = new TGraphErrors(auxi_resp, xRefPt, Response, exRefPt, ResponseError);
128  sprintf(name, "RespFit");
129  RespFit = new TF1(name, "[0]-[1]/(pow(log10(x),[2])+[3])+[4]/x", xRefPt[0], xRefPt[auxi_resp - 1]);
130  RespFit->SetParameter(0, 1.);
131  RespFit->SetParameter(1, 1.);
132  RespFit->SetParameter(2, 1.);
133  RespFit->SetParameter(3, 1.);
134  RespFit->SetParameter(4, 1.);
135  std::cout << "Fitting response in the range: " << xRefPt[0] << " " << xRefPt[auxi_resp - 1] << std::endl;
136  g_Resp->Fit(RespFit, "RQ");
137  fitter = TVirtualFitter::GetFitter();
138  COV_Resp = new TMatrixD(5, 5, fitter->GetCovarianceMatrix());
140  MinCaloPt = 4.;
141  MaxCaloPt = 5000.;
142  MinRefPt = 4.;
143  MaxRefPt = 5000.;
144  CorFit->SetRange(MinCaloPt, MaxCaloPt);
145  RespFit->SetRange(MinRefPt, MaxRefPt);
147  L3CorrectionFile.setf(ios::left);
148  L3CorrectionFile << setw(12) << -5.191 << setw(12) << 5.191 << setw(12) << (int)6 << setw(12) << MinCaloPt << setw(12)
149  << MaxCaloPt << setw(12) << CorFit->GetParameter(0) << setw(12) << CorFit->GetParameter(1)
150  << setw(12) << CorFit->GetParameter(2) << setw(12) << CorFit->GetParameter(3);
151  L3CorrectionFile.close();
152 
153  L3ResponseFile.setf(ios::left);
154  L3ResponseFile << setw(12) << RespFit->GetParameter(0) << setw(12) << RespFit->GetParameter(1) << setw(12)
155  << RespFit->GetParameter(2) << setw(12) << RespFit->GetParameter(3) << setw(12)
156  << RespFit->GetParameter(4);
157  L3ResponseFile.close();
158 
159  outf = new TFile(L3OutputROOTFilename.c_str(), "RECREATE");
160  g_Cor->Write("Correction_vs_CaloPt");
161  COV_Cor->Write("CovMatrix_Correction");
162  g_Resp->Write("Response_vs_RefPt");
163  COV_Resp->Write("CovMatrix_Resp");
164  outf->Close();
165 }
mps_fire.i
i
Definition: mps_fire.py:429
GetRecoTauVFromDQM_MC_cff.next
next
Definition: GetRecoTauVFromDQM_MC_cff.py:31
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
submitPVResolutionJobs.key
key
prepare the HTCondor submission files and eventually submit them
Definition: submitPVResolutionJobs.py:222
gather_cfg.cout
cout
Definition: gather_cfg.py:144
HistoExists
bool HistoExists(std::vector< std::string > LIST, std::string hname)
Definition: Utilities.h:441
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
Generated for CMSSW Reference Manual by doxygen 1.8.14