CMS 3D CMS Logo
Functions
L2Correction.cc File Reference
#include <iostream>
#include <iomanip>
#include <cstring>
#include <fstream>
#include <cmath>
#include <TFile.h>
#include <TH1F.h>
#include <TF1.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TMath.h>
#include <TKey.h>
#include <TList.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 18 of file L2Correction.cc.

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

18  {
19  CommandLine c1;
20  c1.parse(argc, argv);
21 
22  std::string HistoFilename = c1.getValue<string>("HistoFilename");
23  std::string FitterFilename = c1.getValue<string>("FitterFilename");
24  std::string L3ResponseTxtFilename = c1.getValue<string>("L3ResponseTxtFilename");
25  std::string L3CorrectionTxtFilename = c1.getValue<string>("L3CorrectionTxtFilename");
26  std::string L3OutputROOTFilename = c1.getValue<string>("L3OutputROOTFilename");
27  std::string L2CorrectionTxtFilename = c1.getValue<string>("L2CorrectionTxtFilename");
28  std::string L2OutputROOTFilename = c1.getValue<string>("L2OutputROOTFilename");
29  std::vector<double> pt_vec = c1.getVector<double>("RefPtBoundaries");
30  std::vector<double> eta_vec = c1.getVector<double>("EtaBoundaries");
31  if (!c1.check())
32  return 0;
33  c1.print();
35  const int MAX_NETA = 83;
36  const int MAX_NREFPTBINS = 30;
37  const int NCaloPtValues = 18;
38  int NRefPtBins = pt_vec.size() - 1;
39  int NETA = eta_vec.size() - 1;
40  if (NETA < 2)
41  return (0);
42  int i, auxi, etabin, ptbin;
43  char name[100], func[1024];
44  double cor, e_cor;
45  double MinCaloPt[MAX_NETA], MaxCaloPt[MAX_NETA];
46  double correction_x[MAX_NREFPTBINS], correction_y[MAX_NREFPTBINS], correction_ex[MAX_NREFPTBINS],
47  correction_ey[MAX_NREFPTBINS];
48  double cor_rel[MAX_NREFPTBINS], ref_pt, calo_pt, control_pt;
49  double L3_resp[5], L2_cor[10];
50  TKey *key;
51  TFile *inf;
52  TFile *outf;
53  TH1F *hcorrection[MAX_NREFPTBINS];
54  TH1F *h;
55  TF1 *Correction[MAX_NETA];
56  TF1 *L2Correction[MAX_NETA];
57  TF1 *L3Response;
58  std::ifstream L3ResponseFile;
59  std::ofstream L2File;
60  TGraph *g_L2Correction[MAX_NETA];
61  TGraphErrors *g_EtaCorrection[MAX_NETA];
62  double aux_CaloPt[NCaloPtValues] = {
63  10, 15, 20, 30, 40, 50, 75, 100, 150, 200, 300, 400, 500, 750, 1000, 1500, 2000, 3000};
64  std::vector<std::string> HistoNamesList;
65  inf = new TFile(FitterFilename.c_str(), "r");
66  if (inf->IsZombie())
67  return (0);
68  TIter next(inf->GetListOfKeys());
69  while ((key = (TKey *)next()))
70  HistoNamesList.push_back(key->GetName());
72  L3ResponseFile.open(L3ResponseTxtFilename.c_str());
73  L3ResponseFile >> L3_resp[0] >> L3_resp[1] >> L3_resp[2] >> L3_resp[3] >> L3_resp[4];
74  L3Response = new TF1("L3Response", "[0]-[1]/(pow(log10(x),[2])+[3])+[4]/x", pt_vec[0], pt_vec[NRefPtBins - 1]);
75  for (i = 0; i < 5; i++)
76  L3Response->SetParameter(i, L3_resp[i]);
77  L3ResponseFile.close();
79  for (i = 0; i < NRefPtBins; i++) {
80  sprintf(name, "Correction_vs_Eta_RefPt%d", i);
81  if (!HistoExists(HistoNamesList, name))
82  return (0);
83  hcorrection[i] = (TH1F *)inf->Get(name);
84  }
85  for (etabin = 0; etabin < NETA; etabin++) {
86  sprintf(name, "MeanCaloPt_Eta%d", etabin);
87  if (!HistoExists(HistoNamesList, name))
88  return (0);
89  h = (TH1F *)inf->Get(name);
91  auxi = 0;
92  for (ptbin = 0; ptbin < NRefPtBins; ptbin++) {
93  cor = hcorrection[ptbin]->GetBinContent(etabin + 1);
94  e_cor = hcorrection[ptbin]->GetBinError(etabin + 1);
95  if (cor > 0 && e_cor > 0.0001 && e_cor < 0.3) {
96  correction_x[auxi] = h->GetBinContent(ptbin + 1); //average CaloPt for the eta bin
97  correction_ex[auxi] = 0.;
98  correction_y[auxi] = cor;
99  correction_ey[auxi] = e_cor;
100  auxi++;
101  }
102  }
103  sprintf(name, "Correction%d", etabin);
104  if (auxi > 1) {
105  MaxCaloPt[etabin] = correction_x[auxi - 1];
106  MinCaloPt[etabin] = correction_x[1];
107  if (auxi > 10)
108  sprintf(func, "[0]+[1]/(pow(log10(x),[2])+[3])");
109  else
110  sprintf(func, "[0]+[1]*log10(x)+[2]*pow(log10(x),2)");
111  Correction[etabin] = new TF1(name, func, MinCaloPt[etabin], MaxCaloPt[etabin]);
112  } else {
113  std::cout << name << ": not enough points" << std::endl;
114  sprintf(func, "[0]");
115  correction_x[0] = 10;
116  correction_x[1] = 100;
117  correction_ex[0] = 0.;
118  correction_ex[1] = 0.;
119  correction_y[0] = 1.;
120  correction_y[1] = 1.;
121  correction_ey[0] = 0.;
122  correction_ey[1] = 0.;
123  auxi = 2;
124  Correction[etabin] = new TF1(name, func, 10, 100);
125  MaxCaloPt[etabin] = 0;
126  MinCaloPt[etabin] = 0;
127  }
128  g_EtaCorrection[etabin] = new TGraphErrors(auxi, correction_x, correction_y, correction_ex, correction_ey);
129  Correction[etabin]->SetParameter(0, 0.);
130  Correction[etabin]->SetParameter(1, 0.);
131  Correction[etabin]->SetParameter(2, 0.);
132  Correction[etabin]->SetParameter(3, 0.);
133  g_EtaCorrection[etabin]->Fit(name, "RQ");
134  std::cout << name << " fitted....." << std::endl;
136  auxi = 0;
137  for (ptbin = 0; ptbin < NCaloPtValues; ptbin++) {
138  calo_pt = aux_CaloPt[ptbin];
139  if (calo_pt >= MinCaloPt[etabin] && calo_pt <= MaxCaloPt[etabin]) {
140  ref_pt = calo_pt * Correction[etabin]->Eval(calo_pt);
141  control_pt = ref_pt * (L3Response->Eval(ref_pt));
142  cor_rel[auxi] = control_pt / calo_pt;
143  correction_x[auxi] = calo_pt;
144  auxi++;
145  }
146  }
147  sprintf(name, "L2Correction%d", etabin);
148  if (auxi >= 2) {
149  sprintf(func, "[0]+[1]*log10(x)+[2]*pow(log10(x),2)");
150  if (auxi == 2)
151  sprintf(func, "[0]+[1]*log10(x)");
152  } else {
153  sprintf(func, "[0]");
154  correction_x[0] = 10;
155  correction_x[1] = 100;
156  cor_rel[0] = 1.;
157  cor_rel[1] = 1.;
158  auxi = 2;
159  }
160  g_L2Correction[etabin] = new TGraph(auxi, correction_x, cor_rel);
161  L2Correction[etabin] = new TF1(name, func, correction_x[0], correction_x[auxi - 1]);
162  L2Correction[etabin]->SetParameter(0, 0.);
163  L2Correction[etabin]->SetParameter(1, 0.);
164  L2Correction[etabin]->SetParameter(2, 0.);
165  L2Correction[etabin]->SetParameter(3, 0.);
166  L2Correction[etabin]->SetParameter(4, 0.);
167  L2Correction[etabin]->SetParameter(5, 0.);
168  g_L2Correction[etabin]->Fit(name, "RQ");
169  std::cout << name << " fitted....." << std::endl;
170  } //end of eta bin loop
172  L2File.open(L2CorrectionTxtFilename.c_str());
173  L2File.setf(ios::right);
174  for (etabin = 0; etabin < NETA; etabin++) {
175  for (i = 0; i < 6; i++)
176  L2_cor[i] = L2Correction[etabin]->GetParameter(i);
177  L2File << setw(11) << eta_vec[etabin] << setw(11) << eta_vec[etabin + 1] << setw(11) << (int)8 << setw(12)
178  << MinCaloPt[etabin] << setw(12) << MaxCaloPt[etabin] << setw(13) << L2_cor[0] << setw(13) << L2_cor[1]
179  << setw(13) << L2_cor[2] << setw(13) << L2_cor[3] << setw(13) << L2_cor[4] << setw(13) << L2_cor[5] << "\n";
180  }
181  L2File.close();
182  std::cout << L2CorrectionTxtFilename << " written...." << std::endl;
183  outf = new TFile(L2OutputROOTFilename.c_str(), "RECREATE");
184  for (etabin = 0; etabin < NETA; etabin++) {
185  sprintf(name, "Correction_EtaBin%d", etabin);
186  g_EtaCorrection[etabin]->Write(name);
187  sprintf(name, "L2Correction_EtaBin%d", etabin);
188  g_L2Correction[etabin]->Write(name);
189  }
190  outf->Close();
191 }
EcalMonitorTask_cff.func
func
Definition: EcalMonitorTask_cff.py:10
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
h
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
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