CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
GenSpecificAlgo.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: METAlgorithms
4 // Class: GenSpecificAlgo
5 //
6 // Original Authors: R. Cavanaugh (taken from F.Ratnikov, UMd)
7 // Created: June 6, 2006
8 //
9 //
10 //____________________________________________________________________________||
12 #include "TMath.h"
13 
14 #include <set>
15 
16 //____________________________________________________________________________||
17 reco::GenMET GenSpecificAlgo::addInfo(edm::Handle<edm::View<reco::Candidate> > particles, CommonMETData *met, double globalThreshold, bool onlyFiducial,bool applyFiducialThresholdForFractions, bool usePt)
18 {
19  fillCommonMETData(met, particles, globalThreshold, onlyFiducial, usePt);
20 
21  SpecificGenMETData specific = mkSpecificGenMETData(particles, globalThreshold, onlyFiducial, applyFiducialThresholdForFractions, usePt);
22 
23  const LorentzVector p4( met->mex, met->mey, met->mez, met->met );
24  const Point vtx( 0.0, 0.0, 0.0 );
25 
26  reco::GenMET genMet(specific, met->sumet, p4, vtx );
27  return genMet;
28 }
29 
30 //____________________________________________________________________________||
32 {
33  double sum_et = 0.0;
34  double sum_ex = 0.0;
35  double sum_ey = 0.0;
36  double sum_ez = 0.0;
37 
38  for( edm::View<reco::Candidate>::const_iterator iParticle = (particles.product())->begin(); iParticle != (particles.product())->end(); ++iParticle)
39  {
40  if( (onlyFiducial && TMath::Abs(iParticle->eta()) >= 5.0)) continue;
41 
42  if( iParticle->et() <= globalThreshold ) continue;
43 
44  if(usePt)
45  {
46  double phi = iParticle->phi();
47  double et = iParticle->pt();
48  sum_ez += iParticle->pz();
49  sum_et += et;
50  sum_ex += et*cos(phi);
51  sum_ey += et*sin(phi);
52  }
53  else
54  {
55  double phi = iParticle->phi();
56  double theta = iParticle->theta();
57  double e = iParticle->energy();
58  double et = e*sin(theta);
59  sum_ez += e*cos(theta);
60  sum_et += et;
61  sum_ex += et*cos(phi);
62  sum_ey += et*sin(phi);
63  }
64  }
65 
66  met->mex = -sum_ex;
67  met->mey = -sum_ey;
68  met->mez = -sum_ez;
69  met->met = sqrt( sum_ex*sum_ex + sum_ey*sum_ey );
70  met->sumet = sum_et;
71 }
72 
73 //____________________________________________________________________________||
74 SpecificGenMETData GenSpecificAlgo::mkSpecificGenMETData(edm::Handle<edm::View<reco::Candidate> >& particles,double globalThreshold, bool onlyFiducial,bool applyFiducialThresholdForFractions, bool usePt)
75 {
76  const static int neutralEMpdgId[] = { 22 /* photon */ };
77  const static std::set<int> neutralEMpdgIdSet(neutralEMpdgId, neutralEMpdgId + sizeof(neutralEMpdgId)/sizeof(int));
78 
79  const static int chargedEMpdgId[] = { 11 /* e */ };
80  const static std::set<int> chargedEMpdgIdSet(chargedEMpdgId, chargedEMpdgId + sizeof(chargedEMpdgId)/sizeof(int));
81 
82  const static int muonpdgId[] = { 13 /* muon */ };
83  const static std::set<int> muonpdgIdSet(muonpdgId, muonpdgId + sizeof(muonpdgId)/sizeof(int));
84 
85  const static int neutralHADpdgId[] = {
86  130 /* K_long */,
87  310 /* K_short */,
88  3122 /* Lambda */,
89  2112 /* n */,
90  3222 /* Neutral Cascade */
91  };
92  const static std::set<int> neutralHADpdgIdSet(neutralHADpdgId, neutralHADpdgId + sizeof(neutralHADpdgId)/sizeof(int));
93 
94  const static int chargedHADpdgId[] = {
95  211 /* pi */,
96  321 /* K+/K- */,
97  2212 /* p */,
98  3312 /* Cascade - */,
99  3112 /* Sigma - */,
100  3322 /* Sigma + */,
101  3334 /* Omega - */
102  };
103  const static std::set<int> chargedHADpdgIdSet(chargedHADpdgId, chargedHADpdgId + sizeof(chargedHADpdgId)/sizeof(int));
104 
105  const static int invisiblepdgId[] = {
106  12 /* e_nu */,
107  14 /* mu_nu */,
108  16 /* tau_nu */,
109  1000022 /* Neutral ~Chi_0 */,
110  1000012 /* LH ~e_nu */,
111  1000014 /* LH ~mu_nu */,
112  1000016 /* LH ~tau_nu */,
113  2000012 /* RH ~e_nu */,
114  2000014 /* RH ~mu_nu */,
115  2000016 /* RH ~tau_nu */,
116  39 /* G */,
117  1000039 /* ~G */,
118  5100039 /* KK G */,
119  4000012 /* excited e_nu */,
120  4000014 /* excited mu_nu */,
121  4000016 /* excited tau_nu */,
122  9900012 /* Maj e_nu */,
123  9900014 /* Maj mu_nu */,
124  9900016 /* Maj tau_nu */,
125  };
126  const static std::set<int> invisiblepdgIdSet(invisiblepdgId, invisiblepdgId + sizeof(invisiblepdgId)/sizeof(int));
127 
129  double Et_unclassified = 0.0;
130 
131  for(edm::View<reco::Candidate>::const_iterator iParticle = (particles.product())->begin(); iParticle != (particles.product())->end(); ++iParticle)
132  {
133  if(applyFiducialThresholdForFractions) if( onlyFiducial && (TMath::Abs(iParticle->eta()) >= 5.0) ) continue;
134  if(applyFiducialThresholdForFractions) if( iParticle->et() <= globalThreshold ) continue;
135 
136  int pdgId = TMath::Abs( iParticle->pdgId() ) ;
137  double pt = (usePt) ? iParticle->pt() : iParticle->et();
138  if(neutralEMpdgIdSet.count(pdgId)) specific.NeutralEMEtFraction += pt;
139  else if(chargedEMpdgIdSet.count(pdgId)) specific.ChargedEMEtFraction += pt;
140  else if(muonpdgIdSet.count(pdgId)) specific.MuonEtFraction += pt;
141  else if(neutralHADpdgIdSet.count(pdgId)) specific.NeutralHadEtFraction += pt;
142  else if(chargedHADpdgIdSet.count(pdgId)) specific.ChargedHadEtFraction += pt;
143  else if(invisiblepdgIdSet.count(pdgId)) specific.InvisibleEtFraction += pt;
144  else Et_unclassified += pt;
145  }
146 
147  double Et_Total = specific.NeutralEMEtFraction + specific.NeutralHadEtFraction + specific.ChargedEMEtFraction +
148  specific.ChargedHadEtFraction + specific.MuonEtFraction + specific.InvisibleEtFraction + Et_unclassified;
149 
150  if(Et_Total)
151  {
152  specific.NeutralEMEtFraction /= Et_Total;
153  specific.NeutralHadEtFraction /= Et_Total;
154  specific.ChargedEMEtFraction /= Et_Total;
155  specific.ChargedHadEtFraction /= Et_Total;
156  specific.MuonEtFraction /= Et_Total;
157  specific.InvisibleEtFraction /= Et_Total;
158  }
159 
160  return specific;
161 }
162 
163 //____________________________________________________________________________||
dictionary specific
reco::GenMET addInfo(edm::Handle< edm::View< reco::Candidate > > particles, CommonMETData *met, double globalThreshold=0, bool onlyFiducial=false, bool applyFiducialThresholdForFractions=false, bool usePt=false)
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
math::XYZPoint Point
Geom::Theta< T > theta() const
T sqrt(T t)
Definition: SSEVec.h:48
double p4[4]
Definition: TauolaWrapper.h:92
void fillCommonMETData(CommonMETData *met, edm::Handle< edm::View< reco::Candidate > > &particles, double globalThreshold, bool onlyFiducial, bool usePt)
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
T Abs(T a)
Definition: MathUtil.h:49
Structure containing data common to all types of MET.
Definition: CommonMETData.h:12
#define end
Definition: vmac.h:37
TLorentzVector genMet(const HepMC::GenEvent *all, double etamin=-9999., double etamax=9999.)
#define begin
Definition: vmac.h:30
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:85
SpecificGenMETData mkSpecificGenMETData(edm::Handle< edm::View< reco::Candidate > > &particles, double globalThreshold, bool onlyFiducial, bool applyFiducialThresholdForFractions, bool usePt)
math::XYZTLorentzVector LorentzVector