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