CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SignAlgoResolutions.cc
Go to the documentation of this file.
2 // -*- C++ -*-
3 //
4 // Package: METAlgorithms
5 // Class: SignAlgoResolutions
6 //
14 //
15 // Original Author: Kyle Story, Freya Blekman (Cornell University)
16 // Created: Fri Apr 18 11:58:33 CEST 2008
17 // $Id: SignAlgoResolutions.cc,v 1.6 2009/11/23 14:38:07 fblekman Exp $
18 //
19 //
23 
24 #include <math.h>
25 
26 #include <cstdlib>
27 #include <iostream>
28 #include <sstream>
29 #include <string>
30 #include <sys/stat.h>
31 
33  functionmap_(),
34  ptResol_(0),
35  phiResol_(0)
36 {
37  addResolutions(iConfig);
38 }
39 
40 
41 
42 double metsig::SignAlgoResolutions::eval(const resolutionType & type, const resolutionFunc & func, const double & et, const double & phi, const double & eta) const {
43  // derive p from et and eta;
44  double theta = 2*atan(exp(-eta));
45  double p = et / sin(theta); // rough assumption: take e and p equivalent...
46  return eval(type,func,et,phi,eta,p);
47 
48 }
49 
50 double metsig::SignAlgoResolutions::eval(const resolutionType & type, const resolutionFunc & func, const double & et, const double & phi, const double & eta, const double &p) const {
51 
52  functionPars x(4);
53  x[0]=et;
54  x[1]=phi;
55  x[2]=eta;
56  x[3]=p;
57  // std::cout << "getting function of type " << type << " " << func << " " << x[0] << " " << x[1] << " " << x[2] << " " << x[3] << std::endl;
58  return getfunc(type,func,x);
59 
60 }
61 
63  double eta = candidate->eta();
64  double phi = candidate->phi();
65  double et = candidate->energy()*sin(candidate->theta());
66  resolutionType thetype;
67  std::string name;
68  int type = candidate->particleId();
69  switch (type)
70  {
71  case 1:
72  thetype=PFtype1;name="PFChargedHadron"; break;
73  case 2:
74  thetype=PFtype2;name="PFChargedEM"; break;
75  case 3:
76  thetype=PFtype3;name="PFMuon"; break;
77  case 4:
78  thetype=PFtype4;name="PFNeutralEM"; break;
79  case 5:
80  thetype=PFtype5;name="PFNeutralHadron"; break;
81  case 6:
82  thetype=PFtype6;name="PFtype6"; break;
83  case 7:
84  thetype=PFtype7;name="PFtype7"; break;
85  default:
86  thetype=PFtype7;name="PFunknown"; break;
87  }
88 
89  double d_et=0, d_phi=0; //d_phi here is the error on phi component of the et
90  reco::TrackRef trackRef = candidate->trackRef();
91  if(!trackRef.isNull() && type!=2){
92  d_et = trackRef->ptError();
93  d_phi = et*trackRef->phiError();
94  }
95  else{
96  d_et = eval(thetype,ET,et,phi,eta);
97  d_phi = eval(thetype,PHI,et,phi,eta);
98  }
99 
100  metsig::SigInputObj resultingobj(name,et,phi,d_et,d_phi);
101  return resultingobj;
102 }
103 
104 
107 
108  double jpt = jet->pt();
109  double jphi = jet->phi();
110  double jeta = jet->eta();
111  double jdeltapt = 999.;
112  double jdeltapphi = 999.;
113 
114  if(jpt<ptResolThreshold_ && jpt<20.){ //use temporary fix for low pT jets
115  double feta = TMath::Abs(jeta);
116  int ieta = feta<5.? int(feta/0.5) : 9; //bin size = 0.5
117  int ipt = jpt>3. ? int(jpt-3./2) : 0; //bin size =2, starting from ptmin=3GeV
118  jdeltapt = jdpt[ieta][ipt];
119  jdeltapphi = jpt*jdphi[ieta][ipt];
120  }
121  else{
122  TF1* fPtEta = ptResol_->parameterEta("sigma",jeta);
123  TF1* fPhiEta = phiResol_->parameterEta("sigma",jeta);
124  jdeltapt = jpt>ptResolThreshold_ ? jpt*fPtEta->Eval(jpt) : jpt*fPtEta->Eval(ptResolThreshold_);
125  jdeltapphi = jpt>ptResolThreshold_ ? jpt*fPhiEta->Eval(jpt) : jpt*fPhiEta->Eval(ptResolThreshold_);
126  delete fPtEta;
127  delete fPhiEta;
128  }
129 
130  std::string inputtype = "jet";
131  metsig::SigInputObj obj_jet(inputtype,jpt,jphi,jdeltapt,jdeltapphi);
132  //std::cout << "RESOLUTIONS JET: " << jpt << " " << jphi<< " " <<jdeltapt << " " << jdeltapphi << std::endl;
133  return obj_jet;
134 }
135 
136 
138  using namespace std;
139 
140  // Jet Resolutions - for now load from the files. Migrate to EventSetup asap.
142 
143  ptResolThreshold_ = iConfig.getParameter<double>("ptresolthreshold");
144 
145 
146  //get temporary low pT pfjet resolutions
147  for (int ieta=0; ieta<10; ieta++){
148  jdpt[ieta] = iConfig.getParameter<std::vector<double> >(Form("jdpt%d", ieta));
149  jdphi[ieta] = iConfig.getParameter<std::vector<double> >(Form("jdphi%d", ieta));
150  }
151 
152 
153  // for now: do this by hand - this can obviously also be done via ESSource etc.
154  functionPars etparameters(3,0);
155  functionPars phiparameters(1,0);
156  // set the parameters per function:
157  // ECAL, BARREL:
158  std::vector<double> ebet = iConfig.getParameter<std::vector<double> >("EB_EtResPar");
159  std::vector<double> ebphi = iConfig.getParameter<std::vector<double> >("EB_PhiResPar");
160 
161  etparameters[0]=ebet[0];
162  etparameters[1]=ebet[1];
163  etparameters[2]=ebet[2];
164  phiparameters[0]=ebphi[0];
165  addfunction(caloEB,ET,etparameters);
166  addfunction(caloEB,PHI,phiparameters);
167  // ECAL, ENDCAP:
168  std::vector<double> eeet = iConfig.getParameter<std::vector<double> >("EE_EtResPar");
169  std::vector<double> eephi = iConfig.getParameter<std::vector<double> >("EE_PhiResPar");
170 
171  etparameters[0]=eeet[0];
172  etparameters[1]=eeet[1];
173  etparameters[2]=eeet[2];
174  phiparameters[0]=eephi[0];
175  addfunction(caloEE,ET,etparameters);
176  addfunction(caloEE,PHI,phiparameters);
177  // HCAL, BARREL:
178  std::vector<double> hbet = iConfig.getParameter<std::vector<double> >("HB_EtResPar");
179  std::vector<double> hbphi = iConfig.getParameter<std::vector<double> >("HB_PhiResPar");
180 
181  etparameters[0]=hbet[0];
182  etparameters[1]=hbet[1];
183  etparameters[2]=hbet[2];
184  phiparameters[0]=hbphi[0];
185  addfunction(caloHB,ET,etparameters);
186  addfunction(caloHB,PHI,phiparameters);
187  // HCAL, ENDCAP:
188  std::vector<double> heet = iConfig.getParameter<std::vector<double> >("HE_EtResPar");
189  std::vector<double> hephi = iConfig.getParameter<std::vector<double> >("HE_PhiResPar");
190 
191  etparameters[0]=heet[0];
192  etparameters[1]=heet[1];
193  etparameters[2]=heet[2];
194  phiparameters[0]=hephi[0];
195  addfunction(caloHE,ET,etparameters);
196  addfunction(caloHE,PHI,phiparameters);
197  // HCAL, Outer
198  std::vector<double> hoet = iConfig.getParameter<std::vector<double> >("HO_EtResPar");
199  std::vector<double> hophi = iConfig.getParameter<std::vector<double> >("HO_PhiResPar");
200 
201 
202  etparameters[0]=hoet[0];
203  etparameters[1]=hoet[1];
204  etparameters[2]=hoet[2];
205  phiparameters[0]=hophi[0];
206  addfunction(caloHO,ET,etparameters);
207  addfunction(caloHO,PHI,phiparameters);
208  // HCAL, Forward
209  std::vector<double> hfet = iConfig.getParameter<std::vector<double> >("HF_EtResPar");
210  std::vector<double> hfphi = iConfig.getParameter<std::vector<double> >("HF_PhiResPar");
211 
212  etparameters[0]=hfet[0];
213  etparameters[1]=hfet[1];
214  etparameters[2]=hfet[2];
215  phiparameters[0]=hfphi[0];
216  addfunction(caloHF,ET,etparameters);
217  addfunction(caloHF,PHI,phiparameters);
218 
219 
220  // PF objects:
221  // type 1:
222  std::vector<double> pf1et = iConfig.getParameter<std::vector<double> >("PF_EtResType1");
223  std::vector<double> pf1phi = iConfig.getParameter<std::vector<double> >("PF_PhiResType1");
224  etparameters[0]=pf1et[0];
225  etparameters[1]=pf1et[1];
226  etparameters[2]=pf1et[2];
227  phiparameters[0]=pf1phi[0];
228  addfunction(PFtype1,ET,etparameters);
229  addfunction(PFtype1,PHI,phiparameters);
230 
231  // PF objects:
232  // type 2:
233  std::vector<double> pf2et = iConfig.getParameter<std::vector<double> >("PF_EtResType2");
234  std::vector<double> pf2phi = iConfig.getParameter<std::vector<double> >("PF_PhiResType2");
235  etparameters[0]=pf2et[0];
236  etparameters[1]=pf2et[1];
237  etparameters[2]=pf2et[2];
238  phiparameters[0]=pf2phi[0];
239  addfunction(PFtype2,ET,etparameters);
240  addfunction(PFtype2,PHI,phiparameters);
241 
242  // PF objects:
243  // type 3:
244  std::vector<double> pf3et = iConfig.getParameter<std::vector<double> >("PF_EtResType3");
245  std::vector<double> pf3phi = iConfig.getParameter<std::vector<double> >("PF_PhiResType3");
246  etparameters[0]=pf3et[0];
247  etparameters[1]=pf3et[1];
248  etparameters[2]=pf3et[2];
249  phiparameters[0]=pf3phi[0];
250  addfunction(PFtype3,ET,etparameters);
251  addfunction(PFtype3,PHI,phiparameters);
252 
253  // PF objects:
254  // type 4:
255  std::vector<double> pf4et = iConfig.getParameter<std::vector<double> >("PF_EtResType4");
256  std::vector<double> pf4phi = iConfig.getParameter<std::vector<double> >("PF_PhiResType4");
257  etparameters[0]=pf4et[0];
258  etparameters[1]=pf4et[1];
259  etparameters[2]=pf4et[2];
260  //phiparameters[0]=pf4phi[0];
261  addfunction(PFtype4,ET,etparameters);
262  addfunction(PFtype4,PHI,pf4phi); //use the same functional form for photon phi error as for pT, pass whole vector
263 
264  // PF objects:
265  // type 5:
266  std::vector<double> pf5et = iConfig.getParameter<std::vector<double> >("PF_EtResType5");
267  std::vector<double> pf5phi = iConfig.getParameter<std::vector<double> >("PF_PhiResType5");
268  etparameters[0]=pf5et[0];
269  etparameters[1]=pf5et[1];
270  etparameters[2]=pf5et[2];
271  phiparameters[0]=pf5phi[0];
272  addfunction(PFtype5,ET,etparameters);
273  addfunction(PFtype5,PHI,pf5phi);
274 
275  // PF objects:
276  // type 6:
277  std::vector<double> pf6et = iConfig.getParameter<std::vector<double> >("PF_EtResType6");
278  std::vector<double> pf6phi = iConfig.getParameter<std::vector<double> >("PF_PhiResType6");
279  etparameters[0]=pf6et[0];
280  etparameters[1]=pf6et[1];
281  etparameters[2]=pf6et[2];
282  phiparameters[0]=pf6phi[0];
283  addfunction(PFtype6,ET,etparameters);
284  addfunction(PFtype6,PHI,phiparameters);
285 
286 
287  // PF objects:
288  // type 7:
289  std::vector<double> pf7et = iConfig.getParameter<std::vector<double> >("PF_EtResType7");
290  std::vector<double> pf7phi = iConfig.getParameter<std::vector<double> >("PF_PhiResType7");
291  etparameters[0]=pf7et[0];
292  etparameters[1]=pf7et[1];
293  etparameters[2]=pf7et[2];
294  phiparameters[0]=pf7phi[0];
295  addfunction(PFtype7,ET,etparameters);
296  addfunction(PFtype7,PHI,phiparameters);
297 
298  return;
299 }
300 
302 
303  functionCombo mypair(type,func);
304  functionmap_[mypair]=parameters;
305 
306 }
307 
309 
310  double result=0;
311  functionCombo mypair(type,func);
312 
313  if(functionmap_.count(mypair)==0){
314  return result;
315  }
316 
317  functionPars values = (functionmap_.find(mypair))->second;
318  switch ( func ){
319  case metsig::ET :
320  return EtFunction(x,values);
321  case metsig::PHI :
322  return PhiFunction(x,values);
323  case metsig::TRACKP :
324  return PFunction(x,values);
325  case metsig::CONSTPHI :
326  return PhiConstFunction(x,values);
327  }
328 
329  // std::cout << "returning function " << type << " " << func << " " << result << " " << x[0] << std::endl;
330 
331  return result;
332 }
333 
335 {
336  if(par.size()<3)
337  return 0.;
338  if(x.size()<1)
339  return 0.;
340  double et=x[0];
341  if(et<=0.)
342  return 0.;
343  double result = et*sqrt((par[2]*par[2])+(par[1]*par[1]/et)+(par[0]*par[0]/(et*et)));
344  return result;
345 }
346 
347 
349 {
350  double et=x[0];
351  if(et<=0.){
352  return 0.;
353  }
354 
355  //if 1 parameter is C provided, returns C*pT, if three parameters N, S, C are provided, it returns the usual resolution value, as for sigmaPt
356  if(par.size()!=1 && par.size()!=3){//only 1 or 3 parameters supported for phi function
357  return 0.;
358  }
359  else if(par.size()==1){
360  return par[0]*et;
361  }
362  else{
363  return et*sqrt((par[2]*par[2])+(par[1]*par[1]/et)+(par[0]*par[0]/(et*et)));
364  }
365 
366 }
368 {
369  // not currently implemented
370  return 0;
371 }
372 
374 {
375  return par[0];
376 }
377 
378 void
380 
381  using namespace std;
382 
383  // only reinitialize the resolutsion if the pointers are zero
384  if ( ptResol_ == 0 ) {
385  string resolutionsAlgo = iConfig.getParameter<std::string>("resolutionsAlgo");
386  string resolutionsEra = iConfig.getParameter<std::string>("resolutionsEra");
387 
388  string cmssw_base(getenv("CMSSW_BASE"));
389  string cmssw_release_base(getenv("CMSSW_RELEASE_BASE"));
390  string path = cmssw_base + "/src/CondFormats/JetMETObjects/data";
391  struct stat st;
392  if (stat(path.c_str(),&st)!=0) {
393  path = cmssw_release_base + "/src/CondFormats/JetMETObjects/data";
394  }
395  if (stat(path.c_str(),&st)!=0) {
396  cerr<<"ERROR: tried to set path but failed, abort."<<endl;
397  }
398  string era(resolutionsEra);
399  string alg(resolutionsAlgo);
400  string ptFileName = path + "/" + era + "_PtResolution_" +alg+".txt";
401  string phiFileName = path + "/" + era + "_PhiResolution_"+alg+".txt";
402 
403  ptResol_ = new JetResolution(ptFileName,false);
404  phiResol_ = new JetResolution(phiFileName,false);
405  }
406 }
type
Definition: HCALResponse.h:22
T getParameter(std::string const &) const
double PhiFunction(const functionPars &x, const functionPars &par) const
double getfunc(const resolutionType &type, const resolutionFunc &func, std::vector< double > &x) const
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
metsig::SigInputObj evalPF(const reco::PFCandidate *candidate) const
Geom::Theta< T > theta() const
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
double eval(const resolutionType &type, const resolutionFunc &func, const double &et, const double &phi, const double &eta, const double &p) const
T eta() const
Jets made from PFObjects.
Definition: PFJet.h:22
virtual double eta() const
momentum pseudorapidity
int path() const
Definition: HLTadd.h:3
U second(std::pair< T, U > const &p)
virtual double energy() const
energy
bool isNull() const
Checks for null.
Definition: Ref.h:244
T sqrt(T t)
Definition: SSEVec.h:28
tuple result
Definition: query.py:137
virtual double theta() const
momentum polar angle
double EtFunction(const functionPars &x, const functionPars &par) const
std::vector< double > functionPars
virtual double pt() const
transverse momentum
metsig::SigInputObj evalPFJet(const reco::PFJet *jet) const
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:34
list cmssw_release_base
Definition: cmsBenchmark.py:61
void initializeJetResolutions(const edm::ParameterSet &iConfig)
double PhiConstFunction(const functionPars &x, const functionPars &par) const
virtual ParticleType particleId() const
Definition: PFCandidate.h:299
void addfunction(const resolutionType type, const resolutionFunc func, std::vector< double > parameters)
virtual double phi() const
momentum azimuthal angle
std::pair< metsig::resolutionType, metsig::resolutionFunc > functionCombo
list cmssw_base
Definition: launcher.py:19
double PFunction(const functionPars &x, const functionPars &par) const
const double par[8 *NPar][4]
void addResolutions(const edm::ParameterSet &iConfig)
reco::TrackRef trackRef() const
Definition: PFCandidate.h:134
Definition: DDAxes.h:10