CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Functions
L3Correction.cc File Reference
#include <iostream>
#include <iomanip>
#include <string.h>
#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(), i, EcalCondDB::inf, combine::key, mergeVDriftHistosByStation::name, GetRecoTauVFromDQM_MC_cff::next, CommandLine::parse(), and CommandLine::print().

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],CorrectionError[MAX_NREFPTBINS],Response[MAX_NREFPTBINS],ResponseError[MAX_NREFPTBINS];
41  double xCaloPt[MAX_NREFPTBINS],exCaloPt[MAX_NREFPTBINS],xRefPt[MAX_NREFPTBINS],exRefPt[MAX_NREFPTBINS];
42  double MinCaloPt,MaxCaloPt,MinRefPt,MaxRefPt;
43  double cor,e_cor,resp,e_resp;
44  TFile *inf;
45  TFile *outf;
46  TH1F *hCor,*hResp,*hRef,*hCalo;
47  TF1 *CorFit;
48  TF1 *RespFit;
49  TGraphErrors *g_Cor;
50  TGraphErrors *g_Resp;
51  TVirtualFitter *fitter;
52  TMatrixD *COV_Cor, *COV_Resp;
53  ofstream L3CorrectionFile,L3ResponseFile;
54  L3CorrectionFile.open(L3CorrectionTxtFilename.c_str());
55  L3ResponseFile.open(L3ResponseTxtFilename.c_str());
56  TKey *key;
57  std::vector<std::string> HistoNamesList;
58  inf = new TFile(FitterFilename.c_str(),"r");
59  if (inf->IsZombie()) return(0);
60  TIter next(inf->GetListOfKeys());
61  while ((key = (TKey*)next()))
62  HistoNamesList.push_back(key->GetName());
64  sprintf(name,"MeanCaloPt");
65  if (!HistoExists(HistoNamesList,name)) return(0);
66  hCalo = (TH1F*)inf->Get(name);
68  sprintf(name,"MeanRefPt");
69  if (!HistoExists(HistoNamesList,name)) return(0);
70  hRef = (TH1F*)inf->Get(name);
72  sprintf(name,"Correction");
73  if (!HistoExists(HistoNamesList,name)) return(0);
74  hCor = (TH1F*)inf->Get(name);
75  sprintf(name,"Response");
76  if (!HistoExists(HistoNamesList,name)) return(0);
77  hResp = (TH1F*)inf->Get(name);
78  auxi_cor=0;
79  auxi_resp=0;
80  std::cout<<"RefPt bin"<<setw(12)<<"<CaloPt>"<<setw(12)<<"Correction"<<setw(12)<<"Error"<<setw(12)<<"<RefPt>"<<setw(12)<<"Response"<<setw(12)<<"Error"<<std::endl;
81  for(i=0;i<NRefPtBins;i++)
82  {
83  cor = hCor->GetBinContent(i+1);
84  e_cor = hCor->GetBinError(i+1);
85  std::cout<<"["<<pt_vec[i]<<","<<pt_vec[i+1]<<"]"<<setw(12)<<hCalo->GetBinContent(i+1)<<setw(12)<<cor<<setw(12)<<e_cor;
86  resp = hResp->GetBinContent(i+1);
87  e_resp = hResp->GetBinError(i+1);
88  std::cout<<setw(12)<<hRef->GetBinContent(i+1)<<setw(12)<<resp<<setw(12)<<e_resp<<std::endl;
89  if (cor>0 && e_cor>0.000001 && e_cor<0.2)
90  {
91  Correction[auxi_cor] = cor;
92  CorrectionError[auxi_cor] = e_cor;
93  xCaloPt[auxi_cor] = hCalo->GetBinContent(i+1);
94  exCaloPt[auxi_cor] = hCalo->GetBinError(i+1);
95  auxi_cor++;
96  }
97  if (resp>0 && e_resp>0.000001 && e_resp<0.2)
98  {
99  Response[auxi_resp] = resp;
100  ResponseError[auxi_resp] = e_resp;
101  xRefPt[auxi_resp] = hRef->GetBinContent(i+1);
102  exRefPt[auxi_resp] = hRef->GetBinError(i+1);
103  auxi_resp++;
104  }
105  }
106  g_Cor = new TGraphErrors(auxi_cor,xCaloPt,Correction,exCaloPt,CorrectionError);
107  sprintf(name,"CorFit");
108  CorFit = new TF1(name,"[0]+[1]/(pow(log10(x),[2])+[3])",xCaloPt[1],xCaloPt[auxi_cor-1]);
109  CorFit->SetParameter(0,1.);
110  CorFit->SetParameter(1,7.);
111  CorFit->SetParameter(2,4.);
112  CorFit->SetParameter(3,4.);
113  std::cout<<"Fitting correction in the range: "<<xCaloPt[1]<<" "<<xCaloPt[auxi_cor-1]<<std::endl;
114  g_Cor->Fit(CorFit,"RQ");
115  fitter = TVirtualFitter::GetFitter();
116  COV_Cor = new TMatrixD(4,4,fitter->GetCovarianceMatrix());
117  g_Resp = new TGraphErrors(auxi_resp,xRefPt,Response,exRefPt,ResponseError);
118  sprintf(name,"RespFit");
119  RespFit = new TF1(name,"[0]-[1]/(pow(log10(x),[2])+[3])+[4]/x",xRefPt[0],xRefPt[auxi_resp-1]);
120  RespFit->SetParameter(0,1.);
121  RespFit->SetParameter(1,1.);
122  RespFit->SetParameter(2,1.);
123  RespFit->SetParameter(3,1.);
124  RespFit->SetParameter(4,1.);
125  std::cout<<"Fitting response in the range: "<<xRefPt[0]<<" "<<xRefPt[auxi_resp-1]<<std::endl;
126  g_Resp->Fit(RespFit,"RQ");
127  fitter = TVirtualFitter::GetFitter();
128  COV_Resp = new TMatrixD(5,5,fitter->GetCovarianceMatrix());
130  MinCaloPt = 4.;
131  MaxCaloPt = 5000.;
132  MinRefPt = 4.;
133  MaxRefPt = 5000.;
134  CorFit->SetRange(MinCaloPt,MaxCaloPt);
135  RespFit->SetRange(MinRefPt,MaxRefPt);
137  L3CorrectionFile.setf(ios::left);
138  L3CorrectionFile <<setw(12)<<-5.191
139  <<setw(12)<<5.191
140  <<setw(12)<<(int)6
141  <<setw(12)<<MinCaloPt
142  <<setw(12)<<MaxCaloPt
143  <<setw(12)<<CorFit->GetParameter(0)
144  <<setw(12)<<CorFit->GetParameter(1)
145  <<setw(12)<<CorFit->GetParameter(2)
146  <<setw(12)<<CorFit->GetParameter(3);
147  L3CorrectionFile.close();
148 
149  L3ResponseFile.setf(ios::left);
150  L3ResponseFile <<setw(12)<<RespFit->GetParameter(0)
151  <<setw(12)<<RespFit->GetParameter(1)
152  <<setw(12)<<RespFit->GetParameter(2)
153  <<setw(12)<<RespFit->GetParameter(3)
154  <<setw(12)<<RespFit->GetParameter(4);
155  L3ResponseFile.close();
156 
157  outf = new TFile(L3OutputROOTFilename.c_str(),"RECREATE");
158  g_Cor->Write("Correction_vs_CaloPt");
159  COV_Cor->Write("CovMatrix_Correction");
160  g_Resp->Write("Response_vs_RefPt");
161  COV_Resp->Write("CovMatrix_Resp");
162  outf->Close();
163 }
i
int i
Definition: DBlmapReader.cc:9
CommandLine::check
bool check()
Definition: Utilities.h:240
CommandLine::getVector
std::vector< T > getVector(const std::string &name)
Definition: Utilities.h:144
EcalCondDB.inf
string inf
Definition: EcalCondDB.py:94
CommandLine::print
void print()
Definition: Utilities.h:262
CommandLine::getValue
T getValue(const std::string &name)
Definition: Utilities.h:74
dir2webdir.argc
tuple argc
Definition: dir2webdir.py:41
combine.key
list key
Definition: combine.py:13
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