Egamma/plugins/HLTPMMassFilter.cc

0001 /** \class HLTPMMassFilter
0002  *
0003  * Original Author: Jeremy Werner
0004  * Institution: Princeton University, USA
0005  * Contact: Jeremy.Werner@cern.ch
0006  * Date: February 21, 2007
0007  */
0008 #include "HLTPMMassFilter.h"
0009
0010 #include <cstdlib>
0011 #include <cmath>
0012
0013 //
0014 // constructors and destructor
0015 //
0016 HLTPMMassFilter::HLTPMMassFilter(const edm::ParameterSet& iConfig) : HLTFilter(iConfig), magFieldToken_(esConsumes()) {
0017   candTag_ = iConfig.getParameter<edm::InputTag>("candTag");
0018   beamSpot_ = iConfig.getParameter<edm::InputTag>("beamSpot");
0019   l1EGTag_ = iConfig.getParameter<edm::InputTag>("l1EGCand");
0020
0021   lowerMassCut_ = iConfig.getParameter<double>("lowerMassCut");
0022   upperMassCut_ = iConfig.getParameter<double>("upperMassCut");
0023   nZcandcut_ = iConfig.getParameter<int>("nZcandcut");
0024   reqOppCharge_ = iConfig.getUntrackedParameter<bool>("reqOppCharge", false);
0025   isElectron1_ = iConfig.getUntrackedParameter<bool>("isElectron1", true);
0026   isElectron2_ = iConfig.getUntrackedParameter<bool>("isElectron2", true);
0027
0028   candToken_ = consumes<trigger::TriggerFilterObjectWithRefs>(candTag_);
0029   beamSpotToken_ = consumes<reco::BeamSpot>(beamSpot_);
0030 }
0031
0032 HLTPMMassFilter::~HLTPMMassFilter() = default;
0033
0034 void HLTPMMassFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0035   edm::ParameterSetDescription desc;
0036   makeHLTFilterDescription(desc);
0037   desc.add<edm::InputTag>("candTag", edm::InputTag("hltL1NonIsoDoublePhotonEt5UpsHcalIsolFilter"));
0038   desc.add<edm::InputTag>("beamSpot", edm::InputTag("hltOfflineBeamSpot"));
0039   desc.add<double>("lowerMassCut", 8.0);
0040   desc.add<double>("upperMassCut", 11.0);
0041   desc.add<int>("nZcandcut", 1);
0042   desc.addUntracked<bool>("reqOppCharge", true);
0043   desc.addUntracked<bool>("isElectron1", false);
0044   desc.addUntracked<bool>("isElectron2", false);
0045   desc.add<edm::InputTag>("l1EGCand", edm::InputTag("hltL1IsoRecoEcalCandidate"));
0046   descriptions.add("hltPMMassFilter", desc);
0047 }
0048
0049 // ------------ method called to produce the data  ------------
0050 bool HLTPMMassFilter::hltFilter(edm::Event& iEvent,
0051                                 const edm::EventSetup& iSetup,
0052                                 trigger::TriggerFilterObjectWithRefs& filterproduct) const {
0053   using namespace std;
0054   using namespace edm;
0055   using namespace reco;
0056   using namespace trigger;
0057
0058   // The filter object
0059   if (saveTags()) {
0060     filterproduct.addCollectionTag(l1EGTag_);
0061   }
0062
0063   auto const& theMagField = iSetup.getData(magFieldToken_);
0064
0065   edm::Handle<trigger::TriggerFilterObjectWithRefs> PrevFilterOutput;
0066   iEvent.getByToken(candToken_, PrevFilterOutput);
0067
0068   // beam spot
0069   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0070   iEvent.getByToken(beamSpotToken_, recoBeamSpotHandle);
0071   // gets its position
0072   const GlobalPoint vertexPos(
0073       recoBeamSpotHandle->position().x(), recoBeamSpotHandle->position().y(), recoBeamSpotHandle->position().z());
0074
0075   int n = 0;
0076
0077   // REMOVED USAGE OF STATIC ARRAYS
0078   // double px[66];
0079   // double py[66];
0080   // double pz[66];
0081   // double energy[66];
0082   std::vector<TLorentzVector> pEleCh1;
0083   std::vector<TLorentzVector> pEleCh2;
0084   std::vector<double> charge;
0085
0086   if (isElectron1_ && isElectron2_) {
0087     Ref<ElectronCollection> refele;
0088
0089     vector<Ref<ElectronCollection> > electrons;
0090     PrevFilterOutput->getObjects(TriggerElectron, electrons);
0091
0092     for (auto& electron : electrons) {
0093       refele = electron;
0094
0095       TLorentzVector pThisEle(refele->px(), refele->py(), refele->pz(), refele->energy());
0096       pEleCh1.push_back(pThisEle);
0097       charge.push_back(refele->charge());
0098     }
0099
0100     std::vector<bool> save_cand(electrons.size(), true);
0101     for (unsigned int jj = 0; jj < electrons.size(); jj++) {
0102       for (unsigned int ii = jj + 1; ii < electrons.size(); ii++) {
0103         if (reqOppCharge_ && charge[jj] * charge[ii] > 0)
0104           continue;
0105         if (isGoodPair(pEleCh1[jj], pEleCh1[ii])) {
0106           n++;
0107           for (auto const idx : {jj, ii}) {
0108             if (save_cand[idx])
0109               filterproduct.addObject(TriggerElectron, electrons[idx]);
0110             save_cand[idx] = false;
0111           }
0112         }
0113       }
0114     }
0115
0116   } else {
0117     Ref<RecoEcalCandidateCollection> refsc;
0118
0119     vector<Ref<RecoEcalCandidateCollection> > scs;
0120     PrevFilterOutput->getObjects(TriggerCluster, scs);
0121     if (scs.empty())
0122       PrevFilterOutput->getObjects(TriggerPhoton, scs);  //we dont know if its type trigger cluster or trigger photon
0123
0124     for (auto& i : scs) {
0125       refsc = i;
0126       const reco::SuperClusterRef sc = refsc->superCluster();
0127       TLorentzVector pscPos = approxMomAtVtx(theMagField, vertexPos, sc, 1);
0128       pEleCh1.push_back(pscPos);
0129
0130       TLorentzVector pscEle = approxMomAtVtx(theMagField, vertexPos, sc, -1);
0131       pEleCh2.push_back(pscEle);
0132     }
0133
0134     std::vector<bool> save_cand(scs.size(), true);
0135     for (unsigned int jj = 0; jj < scs.size(); jj++) {
0136       for (unsigned int ii = jj + 1; ii < scs.size(); ii++) {
0137         if (isGoodPair(pEleCh1[jj], pEleCh2[ii]) or isGoodPair(pEleCh2[jj], pEleCh1[ii])) {
0138           n++;
0139           for (auto const idx : {jj, ii}) {
0140             if (save_cand[idx])
0141               filterproduct.addObject(TriggerCluster, scs[idx]);
0142             save_cand[idx] = false;
0143           }
0144         }
0145       }
0146     }
0147   }
0148
0149   // filter decision
0150   bool accept(n >= nZcandcut_);
0151
0152   return accept;
0153 }
0154
0155 bool HLTPMMassFilter::isGoodPair(TLorentzVector const& v1, TLorentzVector const& v2) const {
0156   if (std::abs(v1.E() - v2.E()) < 0.00001)
0157     return false;
0158
0159   auto const mass = (v1 + v2).M();
0160   return (mass >= lowerMassCut_ and mass <= upperMassCut_);
0161 }
0162
0163 TLorentzVector HLTPMMassFilter::approxMomAtVtx(const MagneticField& magField,
0164                                                const GlobalPoint& xvert,
0165                                                const reco::SuperClusterRef sc,
0166                                                int charge) const {
0167   GlobalPoint xsc(sc->position().x(), sc->position().y(), sc->position().z());
0168   float energy = sc->energy();
0169   auto theFTS = trackingTools::ftsFromVertexToPoint(magField, xsc, xvert, energy, charge);
0170   float theApproxMomMod = theFTS.momentum().x() * theFTS.momentum().x() +
0171                           theFTS.momentum().y() * theFTS.momentum().y() + theFTS.momentum().z() * theFTS.momentum().z();
0172   TLorentzVector theApproxMom(
0173       theFTS.momentum().x(), theFTS.momentum().y(), theFTS.momentum().z(), sqrt(theApproxMomMod + 2.61121E-7));
0174   return theApproxMom;
0175 }
0176
0177 // declare this class as a framework plugin
0178 #include "FWCore/Framework/interface/MakerMacros.h"
0179 DEFINE_FWK_MODULE(HLTPMMassFilter);