METAlgorithms/src/GenSpecificAlgo.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    METAlgorithms
0004 // Class:      GenSpecificAlgo
0005 //
0006 // Original Authors:  R. Cavanaugh (taken from F.Ratnikov, UMd)
0007 //          Created:  June 6, 2006
0008 //
0009 //
0010 //____________________________________________________________________________||
0011 #include "RecoMET/METAlgorithms/interface/GenSpecificAlgo.h"
0012 #include "TMath.h"
0013
0014 #include <set>
0015
0016 //____________________________________________________________________________||
0017 reco::GenMET GenSpecificAlgo::addInfo(edm::Handle<edm::View<reco::Candidate> > particles,
0018                                       CommonMETData* met,
0019                                       double globalThreshold,
0020                                       bool onlyFiducial,
0021                                       bool applyFiducialThresholdForFractions,
0022                                       bool usePt) {
0023   fillCommonMETData(met, particles, globalThreshold, onlyFiducial, usePt);
0024
0025   SpecificGenMETData specific =
0026       mkSpecificGenMETData(particles, globalThreshold, onlyFiducial, applyFiducialThresholdForFractions, usePt);
0027
0028   const LorentzVector p4(met->mex, met->mey, met->mez, met->met);
0029   const Point vtx(0.0, 0.0, 0.0);
0030
0031   reco::GenMET genMet(specific, met->sumet, p4, vtx);
0032   return genMet;
0033 }
0034
0035 //____________________________________________________________________________||
0036 void GenSpecificAlgo::fillCommonMETData(CommonMETData* met,
0037                                         edm::Handle<edm::View<reco::Candidate> >& particles,
0038                                         double globalThreshold,
0039                                         bool onlyFiducial,
0040                                         bool usePt) {
0041   double sum_et = 0.0;
0042   double sum_ex = 0.0;
0043   double sum_ey = 0.0;
0044   double sum_ez = 0.0;
0045
0046   for (edm::View<reco::Candidate>::const_iterator iParticle = (particles.product())->begin();
0047        iParticle != (particles.product())->end();
0048        ++iParticle) {
0049     if ((onlyFiducial && TMath::Abs(iParticle->eta()) >= 5.0))
0050       continue;
0051
0052     if (iParticle->et() <= globalThreshold)
0053       continue;
0054
0055     if (usePt) {
0056       double phi = iParticle->phi();
0057       double et = iParticle->pt();
0058       sum_ez += iParticle->pz();
0059       sum_et += et;
0060       sum_ex += et * cos(phi);
0061       sum_ey += et * sin(phi);
0062     } else {
0063       double phi = iParticle->phi();
0064       double theta = iParticle->theta();
0065       double e = iParticle->energy();
0066       double et = e * sin(theta);
0067       sum_ez += e * cos(theta);
0068       sum_et += et;
0069       sum_ex += et * cos(phi);
0070       sum_ey += et * sin(phi);
0071     }
0072   }
0073
0074   met->mex = -sum_ex;
0075   met->mey = -sum_ey;
0076   met->mez = -sum_ez;
0077   met->met = sqrt(sum_ex * sum_ex + sum_ey * sum_ey);
0078   met->sumet = sum_et;
0079 }
0080
0081 //____________________________________________________________________________||
0082 SpecificGenMETData GenSpecificAlgo::mkSpecificGenMETData(edm::Handle<edm::View<reco::Candidate> >& particles,
0083                                                          double globalThreshold,
0084                                                          bool onlyFiducial,
0085                                                          bool applyFiducialThresholdForFractions,
0086                                                          bool usePt) {
0087   const static int neutralEMpdgId[] = {22 /* photon */};
0088   const static std::set<int> neutralEMpdgIdSet(neutralEMpdgId, neutralEMpdgId + sizeof(neutralEMpdgId) / sizeof(int));
0089
0090   const static int chargedEMpdgId[] = {11 /* e */};
0091   const static std::set<int> chargedEMpdgIdSet(chargedEMpdgId, chargedEMpdgId + sizeof(chargedEMpdgId) / sizeof(int));
0092
0093   const static int muonpdgId[] = {13 /* muon */};
0094   const static std::set<int> muonpdgIdSet(muonpdgId, muonpdgId + sizeof(muonpdgId) / sizeof(int));
0095
0096   const static int neutralHADpdgId[] = {
0097       130 /* K_long          */,
0098       310 /* K_short         */,
0099       3122 /* Lambda          */,
0100       2112 /* n               */,
0101       3222 /* Neutral Cascade */
0102   };
0103   const static std::set<int> neutralHADpdgIdSet(neutralHADpdgId,
0104                                                 neutralHADpdgId + sizeof(neutralHADpdgId) / sizeof(int));
0105
0106   const static int chargedHADpdgId[] = {
0107       211 /* pi        */,
0108       321 /* K+/K-     */,
0109       2212 /* p         */,
0110       3312 /* Cascade - */,
0111       3112 /* Sigma -   */,
0112       3322 /* Sigma +   */,
0113       3334 /* Omega -   */
0114   };
0115   const static std::set<int> chargedHADpdgIdSet(chargedHADpdgId,
0116                                                 chargedHADpdgId + sizeof(chargedHADpdgId) / sizeof(int));
0117
0118   const static int invisiblepdgId[] = {
0119       12 /* e_nu            */,      14 /* mu_nu           */,      16 /* tau_nu          */,
0120       1000022 /* Neutral ~Chi_0  */, 1000012 /* LH ~e_nu        */, 1000014 /* LH ~mu_nu       */,
0121       1000016 /* LH ~tau_nu      */, 2000012 /* RH ~e_nu        */, 2000014 /* RH ~mu_nu       */,
0122       2000016 /* RH ~tau_nu      */, 39 /* G               */,      1000039 /* ~G              */,
0123       5100039 /* KK G            */, 4000012 /* excited e_nu    */, 4000014 /* excited mu_nu   */,
0124       4000016 /* excited tau_nu  */, 9900012 /* Maj e_nu        */, 9900014 /* Maj mu_nu       */,
0125       9900016 /* Maj tau_nu      */,
0126   };
0127   const static std::set<int> invisiblepdgIdSet(invisiblepdgId, invisiblepdgId + sizeof(invisiblepdgId) / sizeof(int));
0128
0129   SpecificGenMETData specific;
0130   double Et_unclassified = 0.0;
0131
0132   for (edm::View<reco::Candidate>::const_iterator iParticle = (particles.product())->begin();
0133        iParticle != (particles.product())->end();
0134        ++iParticle) {
0135     if (applyFiducialThresholdForFractions)
0136       if (onlyFiducial && (TMath::Abs(iParticle->eta()) >= 5.0))
0137         continue;
0138     if (applyFiducialThresholdForFractions)
0139       if (iParticle->et() <= globalThreshold)
0140         continue;
0141
0142     int pdgId = TMath::Abs(iParticle->pdgId());
0143     double pt = (usePt) ? iParticle->pt() : iParticle->et();
0144     if (neutralEMpdgIdSet.count(pdgId))
0145       specific.NeutralEMEtFraction += pt;
0146     else if (chargedEMpdgIdSet.count(pdgId))
0147       specific.ChargedEMEtFraction += pt;
0148     else if (muonpdgIdSet.count(pdgId))
0149       specific.MuonEtFraction += pt;
0150     else if (neutralHADpdgIdSet.count(pdgId))
0151       specific.NeutralHadEtFraction += pt;
0152     else if (chargedHADpdgIdSet.count(pdgId))
0153       specific.ChargedHadEtFraction += pt;
0154     else if (invisiblepdgIdSet.count(pdgId))
0155       specific.InvisibleEtFraction += pt;
0156     else
0157       Et_unclassified += pt;
0158   }
0159
0160   double Et_Total = specific.NeutralEMEtFraction + specific.NeutralHadEtFraction + specific.ChargedEMEtFraction +
0161                     specific.ChargedHadEtFraction + specific.MuonEtFraction + specific.InvisibleEtFraction +
0162                     Et_unclassified;
0163
0164   if (Et_Total) {
0165     specific.NeutralEMEtFraction /= Et_Total;
0166     specific.NeutralHadEtFraction /= Et_Total;
0167     specific.ChargedEMEtFraction /= Et_Total;
0168     specific.ChargedHadEtFraction /= Et_Total;
0169     specific.MuonEtFraction /= Et_Total;
0170     specific.InvisibleEtFraction /= Et_Total;
0171   }
0172
0173   return specific;
0174 }
0175
0176 //____________________________________________________________________________||