CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
L2Correction.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <iomanip>
3 #include <string.h>
4 #include <fstream>
5 #include <cmath>
6 #include <TFile.h>
7 #include <TH1F.h>
8 #include <TF1.h>
9 #include <TGraph.h>
10 #include <TGraphErrors.h>
11 #include <TMath.h>
12 #include <TKey.h>
13 #include <TList.h>
14 #include "Utilities.h"
15 
16 using namespace std;
17 
18 int main(int argc, char**argv)
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()) 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) return(0);
41  int i,auxi,etabin,ptbin;
42  char name[100],func[1024];
43  double cor,e_cor;
44  double MinCaloPt[MAX_NETA],MaxCaloPt[MAX_NETA];
45  double correction_x[MAX_NREFPTBINS],correction_y[MAX_NREFPTBINS],correction_ex[MAX_NREFPTBINS],correction_ey[MAX_NREFPTBINS];
46  double cor_rel[MAX_NREFPTBINS],ref_pt,calo_pt,control_pt;
47  double L3_resp[5],L2_cor[10];
48  TKey *key;
49  TFile *inf;
50  TFile *outf;
51  TH1F *hcorrection[MAX_NREFPTBINS];
52  TH1F *h;
53  TF1 *Correction[MAX_NETA];
54  TF1 *L2Correction[MAX_NETA];
55  TF1 *L3Response;
56  std::ifstream L3ResponseFile;
57  std::ofstream L2File;
58  TGraph *g_L2Correction[MAX_NETA];
59  TGraphErrors *g_EtaCorrection[MAX_NETA];
60  double aux_CaloPt[NCaloPtValues] = {10,15,20,30,40,50,75,100,150,200,300,400,500,750,1000,1500,2000,3000};
61  std::vector<std::string> HistoNamesList;
62  inf = new TFile(FitterFilename.c_str(),"r");
63  if (inf->IsZombie()) return(0);
64  TIter next(inf->GetListOfKeys());
65  while ((key = (TKey*)next()))
66  HistoNamesList.push_back(key->GetName());
68  L3ResponseFile.open(L3ResponseTxtFilename.c_str());
69  L3ResponseFile>>L3_resp[0]>>L3_resp[1]>>L3_resp[2]>>L3_resp[3]>>L3_resp[4];
70  L3Response = new TF1("L3Response","[0]-[1]/(pow(log10(x),[2])+[3])+[4]/x",pt_vec[0],pt_vec[NRefPtBins-1]);
71  for(i=0;i<5;i++)
72  L3Response->SetParameter(i,L3_resp[i]);
73  L3ResponseFile.close();
75  for (i=0;i<NRefPtBins;i++)
76  {
77  sprintf(name,"Correction_vs_Eta_RefPt%d",i);
78  if (!HistoExists(HistoNamesList,name)) return(0);
79  hcorrection[i] = (TH1F*)inf->Get(name);
80  }
81  for (etabin=0;etabin<NETA;etabin++)
82  {
83  sprintf(name,"MeanCaloPt_Eta%d",etabin);
84  if (!HistoExists(HistoNamesList,name)) return(0);
85  h = (TH1F*)inf->Get(name);
87  auxi = 0;
88  for (ptbin=0;ptbin<NRefPtBins;ptbin++)
89  {
90  cor = hcorrection[ptbin]->GetBinContent(etabin+1);
91  e_cor = hcorrection[ptbin]->GetBinError(etabin+1);
92  if (cor>0 && e_cor>0.0001 && e_cor<0.3)
93  {
94  correction_x[auxi] = h->GetBinContent(ptbin+1);//average CaloPt for the eta bin
95  correction_ex[auxi] = 0.;
96  correction_y[auxi] = cor;
97  correction_ey[auxi] = e_cor;
98  auxi++;
99  }
100  }
101  sprintf(name,"Correction%d",etabin);
102  if (auxi>1)
103  {
104  MaxCaloPt[etabin]=correction_x[auxi-1];
105  MinCaloPt[etabin]=correction_x[1];
106  if (auxi>10)
107  sprintf(func,"[0]+[1]/(pow(log10(x),[2])+[3])");
108  else
109  sprintf(func,"[0]+[1]*log10(x)+[2]*pow(log10(x),2)");
110  Correction[etabin] = new TF1(name,func,MinCaloPt[etabin],MaxCaloPt[etabin]);
111  }
112  else
113  {
114  std::cout<<name<<": not enough points"<<std::endl;
115  sprintf(func,"[0]");
116  correction_x[0] = 10;
117  correction_x[1] = 100;
118  correction_ex[0] = 0.;
119  correction_ex[1] = 0.;
120  correction_y[0] = 1.;
121  correction_y[1] = 1.;
122  correction_ey[0] = 0.;
123  correction_ey[1] = 0.;
124  auxi = 2;
125  Correction[etabin] = new TF1(name,func,10,100);
126  MaxCaloPt[etabin] = 0;
127  MinCaloPt[etabin] = 0;
128  }
129  g_EtaCorrection[etabin] = new TGraphErrors(auxi,correction_x,correction_y,correction_ex,correction_ey);
130  Correction[etabin]->SetParameter(0,0.);
131  Correction[etabin]->SetParameter(1,0.);
132  Correction[etabin]->SetParameter(2,0.);
133  Correction[etabin]->SetParameter(3,0.);
134  g_EtaCorrection[etabin]->Fit(name,"RQ");
135  std::cout<<name<<" fitted....."<<std::endl;
137  auxi = 0;
138  for(ptbin=0;ptbin<NCaloPtValues;ptbin++)
139  {
140  calo_pt = aux_CaloPt[ptbin];
141  if (calo_pt>=MinCaloPt[etabin] && calo_pt<=MaxCaloPt[etabin])
142  {
143  ref_pt = calo_pt*Correction[etabin]->Eval(calo_pt);
144  control_pt = ref_pt*(L3Response->Eval(ref_pt));
145  cor_rel[auxi] = control_pt/calo_pt;
146  correction_x[auxi] = calo_pt;
147  auxi++;
148  }
149  }
150  sprintf(name,"L2Correction%d",etabin);
151  if (auxi>=2)
152  {
153  sprintf(func,"[0]+[1]*log10(x)+[2]*pow(log10(x),2)");
154  if (auxi==2)
155  sprintf(func,"[0]+[1]*log10(x)");
156  }
157  else
158  {
159  sprintf(func,"[0]");
160  correction_x[0] = 10;
161  correction_x[1] = 100;
162  cor_rel[0] = 1.;
163  cor_rel[1] = 1.;
164  auxi = 2;
165  }
166  g_L2Correction[etabin] = new TGraph(auxi,correction_x,cor_rel);
167  L2Correction[etabin] = new TF1(name,func,correction_x[0],correction_x[auxi-1]);
168  L2Correction[etabin]->SetParameter(0,0.);
169  L2Correction[etabin]->SetParameter(1,0.);
170  L2Correction[etabin]->SetParameter(2,0.);
171  L2Correction[etabin]->SetParameter(3,0.);
172  L2Correction[etabin]->SetParameter(4,0.);
173  L2Correction[etabin]->SetParameter(5,0.);
174  g_L2Correction[etabin]->Fit(name,"RQ");
175  std::cout<<name<<" fitted....."<<std::endl;
176  }//end of eta bin loop
178  L2File.open(L2CorrectionTxtFilename.c_str());
179  L2File.setf(ios::right);
180  for(etabin=0;etabin<NETA;etabin++)
181  {
182  for(i=0;i<6;i++)
183  L2_cor[i] = L2Correction[etabin]->GetParameter(i);
184  L2File << setw(11) << eta_vec[etabin]
185  << setw(11) << eta_vec[etabin+1]
186  << setw(11) << (int)8
187  << setw(12) << MinCaloPt[etabin]
188  << setw(12) << MaxCaloPt[etabin]
189  << setw(13) << L2_cor[0]
190  << setw(13) << L2_cor[1]
191  << setw(13) << L2_cor[2]
192  << setw(13) << L2_cor[3]
193  << setw(13) << L2_cor[4]
194  << setw(13) << L2_cor[5]
195  << "\n";
196  }
197  L2File.close();
198  std::cout<<L2CorrectionTxtFilename<<" written...."<<std::endl;
199  outf = new TFile(L2OutputROOTFilename.c_str(),"RECREATE");
200  for(etabin=0;etabin<NETA;etabin++)
201  {
202  sprintf(name,"Correction_EtaBin%d",etabin);
203  g_EtaCorrection[etabin]->Write(name);
204  sprintf(name,"L2Correction_EtaBin%d",etabin);
205  g_L2Correction[etabin]->Write(name);
206  }
207  outf->Close();
208 }
i
int i
Definition: DBlmapReader.cc:9
NETA
const int NETA
Definition: CaloTowerBoundries.h:4
CommandLine::check
bool check()
Definition: Utilities.h:240
CommandLine::getVector
std::vector< T > getVector(const std::string &name)
Definition: Utilities.h:144
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
EcalCondDB.inf
string inf
Definition: EcalCondDB.py:94
CommandLine::print
void print()
Definition: Utilities.h:262
h
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
CommandLine::getValue
T getValue(const std::string &name)
Definition: Utilities.h:74
relval_steps.key
string key
FastSim: produces sample of signal events, overlayed with premixed minbias events.
Definition: relval_steps.py:736
dir2webdir.argc
tuple argc
Definition: dir2webdir.py:38
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
GetRecoTauVFromDQM_MC_cff.next
tuple next
Definition: GetRecoTauVFromDQM_MC_cff.py:30
CommandLine::parse
bool parse(int argc, char **argv)
Definition: Utilities.h:199
HistoExists
bool HistoExists(std::vector< std::string > LIST, std::string hname)
Definition: Utilities.h:496
Generated for CMSSW Reference Manual by doxygen 1.8.5