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

int main ( int  argc,
char **  argv 
)

Definition at line 19 of file L3Correction.cc.

References alignmentValidation::c1, CommandLine::check(), gather_cfg::cout, CommandLine::getValue(), CommandLine::getVector(), HistoExists(), mps_fire::i, spu::inf(), createfilelist::int, crabWrapper::key, dataset::name, GetRecoTauVFromDQM_MC_cff::next, CommandLine::parse(), CommandLine::print(), and AlCaHLTBitMon_QueryRunRegistry::string.

20 {
21  CommandLine c1;
22  c1.parse(argc,argv);
23 
24  std::string HistoFilename = c1.getValue<string>("HistoFilename");
25  std::string FitterFilename = c1.getValue<string>("FitterFilename");
26  std::string L3ResponseTxtFilename = c1.getValue<string>("L3ResponseTxtFilename");
27  std::string L3CorrectionTxtFilename = c1.getValue<string>("L3CorrectionTxtFilename");
28  std::string L3OutputROOTFilename = c1.getValue<string>("L3OutputROOTFilename");
29  std::string L2CorrectionTxtFilename = c1.getValue<string>("L2CorrectionTxtFilename");
30  std::string L2OutputROOTFilename = c1.getValue<string>("L2OutputROOTFilename");
31  std::vector<double> pt_vec = c1.getVector<double>("RefPtBoundaries");
32  std::vector<double> eta_vec = c1.getVector<double>("EtaBoundaries");
33  if (!c1.check()) 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()) return(0);
66  TIter next(inf->GetListOfKeys());
67  while ((key = (TKey*)next()))
68  HistoNamesList.push_back(key->GetName());
70  sprintf(name,"MeanCaloPt");
71  if (!HistoExists(HistoNamesList,name)) return(0);
72  hCalo = (TH1F*)inf->Get(name);
74  sprintf(name,"MeanRefPt");
75  if (!HistoExists(HistoNamesList,name)) return(0);
76  hRef = (TH1F*)inf->Get(name);
78  sprintf(name,"Correction");
79  if (!HistoExists(HistoNamesList,name)) return(0);
80  hCor = (TH1F*)inf->Get(name);
81  sprintf(name,"Response");
82  if (!HistoExists(HistoNamesList,name)) return(0);
83  hResp = (TH1F*)inf->Get(name);
84  auxi_cor=0;
85  auxi_resp=0;
86  std::cout<<"RefPt bin"<<setw(12)<<"<CaloPt>"<<setw(12)<<"Correction"<<setw(12)<<"Error"<<setw(12)<<"<RefPt>"<<setw(12)<<"Response"<<setw(12)<<"Error"<<std::endl;
87  for(i=0;i<NRefPtBins;i++)
88  {
89  cor = hCor->GetBinContent(i+1);
90  e_cor = hCor->GetBinError(i+1);
91  std::cout<<"["<<pt_vec[i]<<","<<pt_vec[i+1]<<"]"<<setw(12)<<hCalo->GetBinContent(i+1)<<setw(12)<<cor<<setw(12)<<e_cor;
92  resp = hResp->GetBinContent(i+1);
93  e_resp = hResp->GetBinError(i+1);
94  std::cout<<setw(12)<<hRef->GetBinContent(i+1)<<setw(12)<<resp<<setw(12)<<e_resp<<std::endl;
95  if (cor>0 && e_cor>0.000001 && e_cor<0.2)
96  {
97  Correction[auxi_cor] = cor;
98  CorrectionError[auxi_cor] = e_cor;
99  xCaloPt[auxi_cor] = hCalo->GetBinContent(i+1);
100  exCaloPt[auxi_cor] = hCalo->GetBinError(i+1);
101  auxi_cor++;
102  }
103  if (resp>0 && e_resp>0.000001 && e_resp<0.2)
104  {
105  Response[auxi_resp] = resp;
106  ResponseError[auxi_resp] = e_resp;
107  xRefPt[auxi_resp] = hRef->GetBinContent(i+1);
108  exRefPt[auxi_resp] = hRef->GetBinError(i+1);
109  auxi_resp++;
110  }
111  }
112  g_Cor = new TGraphErrors(auxi_cor,xCaloPt,Correction,exCaloPt,CorrectionError);
113  sprintf(name,"CorFit");
114  CorFit = new TF1(name,"[0]+[1]/(pow(log10(x),[2])+[3])",xCaloPt[1],xCaloPt[auxi_cor-1]);
115  CorFit->SetParameter(0,1.);
116  CorFit->SetParameter(1,7.);
117  CorFit->SetParameter(2,4.);
118  CorFit->SetParameter(3,4.);
119  std::cout<<"Fitting correction in the range: "<<xCaloPt[1]<<" "<<xCaloPt[auxi_cor-1]<<std::endl;
120  g_Cor->Fit(CorFit,"RQ");
121  fitter = TVirtualFitter::GetFitter();
122  COV_Cor = new TMatrixD(4,4,fitter->GetCovarianceMatrix());
123  g_Resp = new TGraphErrors(auxi_resp,xRefPt,Response,exRefPt,ResponseError);
124  sprintf(name,"RespFit");
125  RespFit = new TF1(name,"[0]-[1]/(pow(log10(x),[2])+[3])+[4]/x",xRefPt[0],xRefPt[auxi_resp-1]);
126  RespFit->SetParameter(0,1.);
127  RespFit->SetParameter(1,1.);
128  RespFit->SetParameter(2,1.);
129  RespFit->SetParameter(3,1.);
130  RespFit->SetParameter(4,1.);
131  std::cout<<"Fitting response in the range: "<<xRefPt[0]<<" "<<xRefPt[auxi_resp-1]<<std::endl;
132  g_Resp->Fit(RespFit,"RQ");
133  fitter = TVirtualFitter::GetFitter();
134  COV_Resp = new TMatrixD(5,5,fitter->GetCovarianceMatrix());
136  MinCaloPt = 4.;
137  MaxCaloPt = 5000.;
138  MinRefPt = 4.;
139  MaxRefPt = 5000.;
140  CorFit->SetRange(MinCaloPt,MaxCaloPt);
141  RespFit->SetRange(MinRefPt,MaxRefPt);
143  L3CorrectionFile.setf(ios::left);
144  L3CorrectionFile <<setw(12)<<-5.191
145  <<setw(12)<<5.191
146  <<setw(12)<<(int)6
147  <<setw(12)<<MinCaloPt
148  <<setw(12)<<MaxCaloPt
149  <<setw(12)<<CorFit->GetParameter(0)
150  <<setw(12)<<CorFit->GetParameter(1)
151  <<setw(12)<<CorFit->GetParameter(2)
152  <<setw(12)<<CorFit->GetParameter(3);
153  L3CorrectionFile.close();
154 
155  L3ResponseFile.setf(ios::left);
156  L3ResponseFile <<setw(12)<<RespFit->GetParameter(0)
157  <<setw(12)<<RespFit->GetParameter(1)
158  <<setw(12)<<RespFit->GetParameter(2)
159  <<setw(12)<<RespFit->GetParameter(3)
160  <<setw(12)<<RespFit->GetParameter(4);
161  L3ResponseFile.close();
162 
163  outf = new TFile(L3OutputROOTFilename.c_str(),"RECREATE");
164  g_Cor->Write("Correction_vs_CaloPt");
165  COV_Cor->Write("CovMatrix_Correction");
166  g_Resp->Write("Response_vs_RefPt");
167  COV_Resp->Write("CovMatrix_Resp");
168  outf->Close();
169 }
mps_fire.i
i
Definition: mps_fire.py:269
CommandLine::check
bool check()
Definition: Utilities.h:240
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
CommandLine::getVector
std::vector< T > getVector(const std::string &name)
Definition: Utilities.h:144
GetRecoTauVFromDQM_MC_cff.next
next
Definition: GetRecoTauVFromDQM_MC_cff.py:30
CommandLine::print
void print()
Definition: Utilities.h:262
CommandLine::getValue
T getValue(const std::string &name)
Definition: Utilities.h:74
spu::inf
int inf(FILE *, FILE *)
Definition: SherpackUtilities.cc:70
gather_cfg.cout
cout
Definition: gather_cfg.py:143
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.11