GenFilters/plugins/ZgammaMassFilter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    ZgammaMassFilter
0004 // Class:      ZgammaMassFilter
0005 //
0006 /*
0007
0008  Description: filter events based on the Pythia particle information
0009
0010  Implementation: inherits from generic EDFilter
0011
0012 */
0013 //
0014 // Original Author:  Alexey Ferapontov
0015 //         Created:  Thu July 26 11:57:54 CDT 2012
0016 // $Id: ZgammaMassFilter.h,v 1.1 2012/08/10 12:46:29 lenzip Exp $
0017 //
0018 //
0019
0020 #include "DataFormats/Common/interface/Handle.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/Frameworkfwd.h"
0023 #include "FWCore/Framework/interface/global/EDFilter.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0026 #include "FWCore/Utilities/interface/EDGetToken.h"
0027 #include "FWCore/Utilities/interface/InputTag.h"
0028 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0029
0030 #include "TLorentzVector.h"
0031
0032 #include <algorithm>
0033 #include <cmath>
0034 #include <cstdlib>
0035
0036 class ZgammaMassFilter : public edm::global::EDFilter<> {
0037 public:
0038   explicit ZgammaMassFilter(const edm::ParameterSet&);
0039
0040   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0041
0042 private:
0043   const edm::EDGetTokenT<edm::HepMCProduct> token_;
0044
0045   const double minPhotonPt;
0046   const double minLeptonPt;
0047
0048   const double minPhotonEta;
0049   const double minLeptonEta;
0050
0051   const double maxPhotonEta;
0052   const double maxLeptonEta;
0053
0054   const double minDileptonMass;
0055   const double minZgMass;
0056 };
0057
0058 namespace {
0059   // order std::vector of TLorentzVector elements
0060   class orderByPt {
0061   public:
0062     bool operator()(TLorentzVector const& a, TLorentzVector const& b) {
0063       if (a.Pt() == b.Pt()) {
0064         return a.Pt() < b.Pt();
0065       } else {
0066         return a.Pt() > b.Pt();
0067       }
0068     }
0069   };
0070 }  // namespace
0071
0072 using namespace edm;
0073 using namespace std;
0074
0075 ZgammaMassFilter::ZgammaMassFilter(const edm::ParameterSet& iConfig)
0076     : token_(consumes<edm::HepMCProduct>(iConfig.getParameter<InputTag>("HepMCProduct"))),
0077       minPhotonPt(iConfig.getParameter<double>("minPhotonPt")),
0078       minLeptonPt(iConfig.getParameter<double>("minLeptonPt")),
0079       minPhotonEta(iConfig.getParameter<double>("minPhotonEta")),
0080       minLeptonEta(iConfig.getParameter<double>("minLeptonEta")),
0081       maxPhotonEta(iConfig.getParameter<double>("maxPhotonEta")),
0082       maxLeptonEta(iConfig.getParameter<double>("maxLeptonEta")),
0083       minDileptonMass(iConfig.getParameter<double>("minDileptonMass")),
0084       minZgMass(iConfig.getParameter<double>("minZgMass")) {}
0085
0086 bool ZgammaMassFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0087   bool accepted = false;
0088   Handle<HepMCProduct> evt;
0089   iEvent.getByToken(token_, evt);
0090   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0091
0092   vector<TLorentzVector> Lepton;
0093   Lepton.clear();
0094   vector<TLorentzVector> Photon;
0095   Photon.clear();
0096   vector<float> Charge;
0097   Charge.clear();
0098
0099   for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0100        ++p) {
0101     if ((*p)->status() == 1 && (abs((*p)->pdg_id()) == 11 || abs((*p)->pdg_id()) == 13 || abs((*p)->pdg_id()) == 15)) {
0102       TLorentzVector LeptP((*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz(), (*p)->momentum().e());
0103       if (LeptP.Pt() > minLeptonPt) {
0104         Lepton.push_back(LeptP);
0105       }  // if pt
0106     }    // if lepton
0107
0108     if (abs((*p)->pdg_id()) == 22 && (*p)->status() == 1) {
0109       TLorentzVector PhotP((*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz(), (*p)->momentum().e());
0110       if (PhotP.Pt() > minPhotonPt) {
0111         Photon.push_back(PhotP);
0112       }  // if pt
0113     }    // if photon
0114
0115   }  // loop over particles
0116
0117   // std::cout << "\n" << "Photon size: " << Photon.size() << std::endl;
0118   // for (unsigned int u=0; u<Photon.size(); u++){
0119   //   std::cout << "BEF photon PT: " << Photon[u].Pt() << std::endl;
0120   // }
0121   // std::cout << "\n" << "Lepton size: " << Lepton.size() << std::endl;
0122   // for (unsigned int u=0; u<Lepton.size(); u++){
0123   //   std::cout << "BEF lepton PT: " << Lepton[u].Pt() << std::endl;
0124   // }
0125
0126   // order Lepton and Photon according to Pt
0127   std::stable_sort(Photon.begin(), Photon.end(), orderByPt());
0128   std::stable_sort(Lepton.begin(), Lepton.end(), orderByPt());
0129
0130   //  std::cout << "\n" << std::endl;
0131   //  std::cout << "\n" << "Photon size: " << Photon.size() << std::endl;
0132   //  for (unsigned int u=0; u<Photon.size(); u++){
0133   //    std::cout << "AFT photon PT: " << Photon[u].Pt() << std::endl;
0134   //  }
0135   //  std::cout << "\n" << "Lepton size: " << Lepton.size() << std::endl;
0136   //  for (unsigned int u=0; u<Lepton.size(); u++){
0137   //    std::cout << "AFT lepton PT: " << Lepton[u].Pt() << std::endl;
0138   //  }
0139   //  std::cout << "\n" << std::endl;
0140
0141   if (!Photon.empty() && Lepton.size() > 1 && Photon[0].Pt() > minPhotonPt && Lepton[0].Pt() > minLeptonPt &&
0142       Lepton[1].Pt() > minLeptonPt && Photon[0].Eta() > minPhotonEta && Lepton[0].Eta() > minLeptonEta &&
0143       Lepton[1].Eta() > minLeptonEta && Photon[0].Eta() < maxPhotonEta && Lepton[0].Eta() < maxLeptonEta &&
0144       Lepton[1].Eta() < maxLeptonEta && (Lepton[0] + Lepton[1]).M() > minDileptonMass &&
0145       (Lepton[0] + Lepton[1] + Photon[0]).M() >
0146           minZgMass) {  // satisfy molteplicity, kinematics, and ll llg minimum mass
0147     accepted = true;
0148   }
0149
0150   //  std::cout << "++ returning: " << accepted << "\n" << std::endl;
0151
0152   return accepted;
0153 }
0154
0155 DEFINE_FWK_MODULE(ZgammaMassFilter);