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