CMS 3D CMS Logo
ResponseFitter.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <cstring>
3 #include <fstream>
4 #include <cmath>
5 #include <TFile.h>
6 #include <TH1F.h>
7 #include <TMath.h>
8 #include <TKey.h>
9 #include <TList.h>
10 #include "Utilities.h"
11 
12 using namespace std;
13 
14 int main(int argc, char **argv) {
15  CommandLine c1;
16  c1.parse(argc, argv);
17 
18  std::string HistoFilename = c1.getValue<string>("HistoFilename");
19  std::string FitterFilename = c1.getValue<string>("FitterFilename");
20  std::string L3ResponseTxtFilename = c1.getValue<string>("L3ResponseTxtFilename");
21  std::string L3CorrectionTxtFilename = c1.getValue<string>("L3CorrectionTxtFilename");
22  std::string L3OutputROOTFilename = c1.getValue<string>("L3OutputROOTFilename");
23  std::string L2CorrectionTxtFilename = c1.getValue<string>("L2CorrectionTxtFilename");
24  std::string L2OutputROOTFilename = c1.getValue<string>("L2OutputROOTFilename");
25  bool UseRatioForResponse = c1.getValue<bool>("UseRatioForResponse");
26  std::vector<double> pt_vec = c1.getVector<double>("RefPtBoundaries");
27  std::vector<double> eta_vec = c1.getVector<double>("EtaBoundaries");
28  if (!c1.check())
29  return 0;
30  c1.print();
32  const int MAX_NETA = 83;
33  const int MAX_NREFPTBINS = 30;
34  int NRefPtBins = pt_vec.size() - 1;
35  int NETA = eta_vec.size() - 1;
36  char name[100];
37  int j, etabin;
38  double e, mR, eR, sR, seR, mRBarrel, eRBarrel, sRBarrel, seRBarrel, r, c;
39  double mCaloPt, eCaloPt, sCaloPt, mRefPt, eRefPt, sRefPt;
40  double mRefPtEtaBin, eRefPtEtaBin, sRefPtEtaBin, mCaloPtEtaBin, eCaloPtEtaBin, sCaloPtEtaBin;
41  double EtaBoundaries[MAX_NETA], RefPtBoundaries[MAX_NREFPTBINS];
42  std::vector<std::string> HistoNamesList;
43  for (j = 0; j <= NRefPtBins; j++)
44  RefPtBoundaries[j] = pt_vec[j];
45  for (j = 0; j <= NETA; j++)
46  EtaBoundaries[j] = eta_vec[j];
47  TFile *inf; //Input file containing the response histograms.
48  TFile *outf; //Output file containing the fitter results.
49  TH1F *BarrelResponse; //Histogram with the barrel response in RefPt bins.
50  TH1F *BarrelCorrection; //Histogram with the barrel correction in RefPt bins.
51  TH1F *MeanRefPt_Barrel; //Histogram with the barrel average RefPt in RefPt bins.
52  TH1F *MeanCaloPt_Barrel; //Histogram with the barrel average CaloPt in RefPt bins.
53  TH1F *MeanRefPt_EtaBin[MAX_NETA]; //Histograms with the average RefPt in Eta & RefPt bins.
54  TH1F *MeanCaloPt_EtaBin[MAX_NETA]; //Histograms with the average CaloPt in Eta & RefPt bins.
55  TH1F *ResponseVsEta_RefPt[MAX_NREFPTBINS]; //Histograms with the average response vs Eta in RefPt bins.
56  TH1F *CorrectionVsEta_RefPt[MAX_NREFPTBINS]; //Histograms with the average correction vs Eta in RefPt bins.
57  TH1F *h; //Auxilary histogram;
58  TKey *key;
60  inf = new TFile(HistoFilename.c_str(), "r");
61  if (inf->IsZombie())
62  return (0);
63  TIter next(inf->GetListOfKeys());
64  while ((key = (TKey *)next()))
65  HistoNamesList.push_back(key->GetName());
66  outf = new TFile(FitterFilename.c_str(), "RECREATE");
67  TDirectory *dir_Response =
68  (TDirectory *)outf->mkdir("FittedHistograms"); //Directory in output file to store the fitted histograms.
69  BarrelResponse = new TH1F("Response", "Response", NRefPtBins, RefPtBoundaries);
70  BarrelCorrection = new TH1F("Correction", "Correction", NRefPtBins, RefPtBoundaries);
71  MeanRefPt_Barrel = new TH1F("MeanRefPt", "MeanRefPt", NRefPtBins, RefPtBoundaries);
72  MeanCaloPt_Barrel = new TH1F("MeanCaloPt", "MeanCaloPt", NRefPtBins, RefPtBoundaries);
73  if (NETA > 1) //multiple eta bins: used for L2+L3 correction calculation
74  {
75  std::cout << "************* Fitting Response Histograms in multiple Eta bins. ************" << std::endl;
76  for (etabin = 0; etabin < NETA; etabin++) {
77  sprintf(name, "MeanRefPt_Eta%d", etabin);
78  MeanRefPt_EtaBin[etabin] = new TH1F(name, name, NRefPtBins, RefPtBoundaries);
79  sprintf(name, "MeanCaloPt_Eta%d", etabin);
80  MeanCaloPt_EtaBin[etabin] = new TH1F(name, name, NRefPtBins, RefPtBoundaries);
81  }
82  for (j = 0; j < NRefPtBins; j++) //loop over RefPt bins
83  {
84  std::cout << "RefJetPt Bin: [" << RefPtBoundaries[j] << "," << RefPtBoundaries[j + 1] << "] GeV" << std::endl;
85  sprintf(name, "ptRef_RefPt%d_Barrel", j);
86  if (!HistoExists(HistoNamesList, name))
87  return (0);
88  h = (TH1F *)inf->Get(name);
89  GetMEAN(h, mRefPt, eRefPt, sRefPt);
90  sprintf(name, "ptCalo_RefPt%d_Barrel", j);
91  if (!HistoExists(HistoNamesList, name))
92  return (0);
93  h = (TH1F *)inf->Get(name);
94  GetMEAN(h, mCaloPt, eCaloPt, sCaloPt);
95  sprintf(name, "Response_RefPt%d_Barrel", j);
96  if (!HistoExists(HistoNamesList, name))
97  return (0);
98  h = (TH1F *)inf->Get(name);
99  GetMPV(name, h, dir_Response, mRBarrel, eRBarrel, sRBarrel, seRBarrel);
101  MeanRefPt_Barrel->SetBinContent(j + 1, mRefPt);
102  MeanRefPt_Barrel->SetBinError(j + 1, eRefPt);
104  MeanCaloPt_Barrel->SetBinContent(j + 1, mCaloPt);
105  MeanCaloPt_Barrel->SetBinError(j + 1, eCaloPt);
107  CalculateResponse(UseRatioForResponse, mRefPt, eRefPt, mRBarrel, eRBarrel, r, e);
108  BarrelResponse->SetBinContent(j + 1, r);
109  BarrelResponse->SetBinError(j + 1, e);
111  CalculateCorrection(UseRatioForResponse, mRefPt, eRefPt, mRBarrel, eRBarrel, c, e);
112  BarrelCorrection->SetBinContent(j + 1, c);
113  BarrelCorrection->SetBinError(j + 1, e);
115  sprintf(name, "Response_vs_Eta_RefPt%d", j);
116  ResponseVsEta_RefPt[j] = new TH1F(name, name, NETA, EtaBoundaries);
117  sprintf(name, "Correction_vs_Eta_RefPt%d", j);
118  CorrectionVsEta_RefPt[j] = new TH1F(name, name, NETA, EtaBoundaries);
119  for (etabin = 0; etabin < NETA; etabin++) //loop over eta bins
120  {
122  sprintf(name, "Response_RefPt%d_Eta%d", j, etabin);
123  if (!HistoExists(HistoNamesList, name))
124  return (0);
125  h = (TH1F *)inf->Get(name);
126  GetMPV(name, h, dir_Response, mR, eR, sR, seR);
127  sprintf(name, "ptRef_RefPt%d_Eta%d", j, etabin);
128  if (!HistoExists(HistoNamesList, name))
129  return (0);
130  h = (TH1F *)inf->Get(name);
131  GetMEAN(h, mRefPtEtaBin, eRefPtEtaBin, sRefPtEtaBin);
132  sprintf(name, "ptCalo_RefPt%d_Eta%d", j, etabin);
133  if (!HistoExists(HistoNamesList, name))
134  return (0);
135  h = (TH1F *)inf->Get(name);
136  GetMEAN(h, mCaloPtEtaBin, eCaloPtEtaBin, sCaloPtEtaBin);
138  MeanRefPt_EtaBin[etabin]->SetBinContent(j + 1, mRefPtEtaBin);
139  MeanRefPt_EtaBin[etabin]->SetBinError(j + 1, eRefPtEtaBin);
141  MeanCaloPt_EtaBin[etabin]->SetBinContent(j + 1, mCaloPtEtaBin);
142  MeanCaloPt_EtaBin[etabin]->SetBinError(j + 1, eCaloPtEtaBin);
144  CalculateResponse(UseRatioForResponse, mRefPtEtaBin, eRefPtEtaBin, mR, eR, r, e);
145  ResponseVsEta_RefPt[j]->SetBinContent(etabin + 1, r);
146  ResponseVsEta_RefPt[j]->SetBinError(etabin + 1, e);
148  CalculateCorrection(UseRatioForResponse, mRefPtEtaBin, eRefPtEtaBin, mR, eR, c, e);
149  CorrectionVsEta_RefPt[j]->SetBinContent(etabin + 1, c);
150  CorrectionVsEta_RefPt[j]->SetBinError(etabin + 1, e);
151  } //end of EtaBin loop
152  } // end of Pt loop
153  } else //single eta bin: used for L3 correction calculation
154  {
155  std::cout << "************* Fitting Response Histograms in single eta bin. ************" << std::endl;
156  for (j = 0; j < NRefPtBins; j++) //loop over Pt bins
157  {
158  std::cout << "RefJetPt Bin: [" << RefPtBoundaries[j] << "," << RefPtBoundaries[j + 1] << "] GeV" << std::endl;
159  sprintf(name, "ptRef_RefPt%d", j);
160  if (!HistoExists(HistoNamesList, name))
161  return (0);
162  h = (TH1F *)inf->Get(name);
163  GetMEAN(h, mRefPt, eRefPt, sRefPt);
164  sprintf(name, "ptCalo_RefPt%d", j);
165  if (!HistoExists(HistoNamesList, name))
166  return (0);
167  h = (TH1F *)inf->Get(name);
168  GetMEAN(h, mCaloPt, eCaloPt, sCaloPt);
169  sprintf(name, "Response_RefPt%d", j);
170  if (!HistoExists(HistoNamesList, name))
171  return (0);
172  h = (TH1F *)inf->Get(name);
173  GetMPV(name, h, dir_Response, mRBarrel, eRBarrel, sRBarrel, seRBarrel);
175  MeanRefPt_Barrel->SetBinContent(j + 1, mRefPt);
176  MeanRefPt_Barrel->SetBinError(j + 1, eRefPt);
178  MeanCaloPt_Barrel->SetBinContent(j + 1, mCaloPt);
179  MeanCaloPt_Barrel->SetBinError(j + 1, eCaloPt);
181  CalculateResponse(UseRatioForResponse, mRefPt, eRefPt, mRBarrel, eRBarrel, r, e);
182  BarrelResponse->SetBinContent(j + 1, r);
183  BarrelResponse->SetBinError(j + 1, e);
185  CalculateCorrection(UseRatioForResponse, mRefPt, eRefPt, mRBarrel, eRBarrel, c, e);
186  BarrelCorrection->SetBinContent(j + 1, c);
187  BarrelCorrection->SetBinError(j + 1, e);
188  } // end of Pt loop
189  }
191  outf->cd();
192  MeanRefPt_Barrel->Write();
193  MeanCaloPt_Barrel->Write();
194  BarrelResponse->Write();
195  BarrelCorrection->Write();
196  if (NETA > 1) {
197  for (etabin = 0; etabin < NETA; etabin++) {
198  MeanRefPt_EtaBin[etabin]->Write();
199  MeanCaloPt_EtaBin[etabin]->Write();
200  }
201  for (j = 0; j < NRefPtBins; j++) {
202  ResponseVsEta_RefPt[j]->Write();
203  CorrectionVsEta_RefPt[j]->Write();
204  }
205  }
206  outf->Close();
207 }
cmsBatch.argv
argv
Definition: cmsBatch.py:279
main
int main(int argc, char **argv)
Definition: ResponseFitter.cc:14
dir2webdir.argc
argc
Definition: dir2webdir.py:39
h
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
Definition: L1TUtmAlgorithmRcd.h:4
gather_cfg.cout
cout
Definition: gather_cfg.py:144
GetMEAN
void GetMEAN(TH1F *histo, double &peak, double &error, double &sigma)
Definition: Utilities.h:397
CalculateCorrection
void CalculateCorrection(bool UseRatioForResponse, double x, double ex, double y, double ey, double &c, double &e)
Definition: Utilities.h:425
h
CommandLine
Definition: Utilities.h:24
HistoExists
bool HistoExists(std::vector< std::string > LIST, std::string hname)
Definition: Utilities.h:440
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
CalculateResponse
void CalculateResponse(bool UseRatioForResponse, double x, double ex, double y, double ey, double &r, double &e)
Definition: Utilities.h:410
GetMPV
void GetMPV(char name[100], TH1F *histo, TDirectory *Dir, double &peak, double &error, double &sigma, double &err_sigma)
Definition: Utilities.h:337
alignmentValidation.c1
c1
do drawing
Definition: alignmentValidation.py:1025
alignCSCRings.r
r
Definition: alignCSCRings.py:93
std
Definition: JetResolutionObject.h:76
dqmiodatasetharvest.inf
inf
Definition: dqmiodatasetharvest.py:38
Utilities.h
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
c
auto & c
Definition: CAHitNtupletGeneratorKernelsImpl.h:46
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
crabWrapper.key
key
Definition: crabWrapper.py:19
NETA
const int NETA
Definition: CaloTowerBoundries.h:9
GetRecoTauVFromDQM_MC_cff.next
next
Definition: GetRecoTauVFromDQM_MC_cff.py:31
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16