GenFilters/plugins/PythiaFilterZJetWithOutBg.cc

0001 /** \class PythiaFilterZJetWithOutBg
0002  *
0003  *  PythiaFilterZJetWithOutBg filter implements generator-level preselections
0004  *  for photon+jet like events to be used in jet energy calibration.
0005  *  Ported from fortran code written by V.Konoplianikov.
0006  *
0007  * \author A.Ulyanov, ITEP
0008  *
0009  ************************************************************/
0010
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "DataFormats/Math/interface/LorentzVector.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/Frameworkfwd.h"
0015 #include "FWCore/Framework/interface/global/EDFilter.h"
0016 #include "FWCore/Framework/interface/MakerMacros.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/Utilities/interface/EDGetToken.h"
0019 #include "FWCore/Utilities/interface/InputTag.h"
0020 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0021
0022 #include <cmath>
0023 #include <cstdlib>
0024 #include <string>
0025 #include <vector>
0026
0027 class PythiaFilterZJetWithOutBg : public edm::global::EDFilter<> {
0028 public:
0029   explicit PythiaFilterZJetWithOutBg(const edm::ParameterSet&);
0030
0031   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0032
0033 private:
0034   const edm::EDGetTokenT<edm::HepMCProduct> token_;
0035   const double etaMuMax;
0036   const double ptMuMin;
0037   const double ptZMin;
0038   const double ptZMax;
0039   const double m_z;
0040   const double dm_z;
0041 };
0042
0043 PythiaFilterZJetWithOutBg::PythiaFilterZJetWithOutBg(const edm::ParameterSet& iConfig)
0044     : token_(consumes<edm::HepMCProduct>(
0045           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0046       etaMuMax(iConfig.getUntrackedParameter<double>("MaxMuonEta", 2.5)),
0047       ptMuMin(iConfig.getUntrackedParameter<double>("MinMuonPt", 3.5)),
0048       ptZMin(iConfig.getUntrackedParameter<double>("MinZPt")),
0049       ptZMax(iConfig.getUntrackedParameter<double>("MaxZPt")),
0050       m_z(91.19),
0051       dm_z(10.) {}
0052
0053 bool PythiaFilterZJetWithOutBg::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0054   bool accepted = false;
0055   edm::Handle<edm::HepMCProduct> evt;
0056   iEvent.getByToken(token_, evt);
0057
0058   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0059
0060   std::vector<const HepMC::GenParticle*> mu;
0061
0062   for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0063        ++p) {
0064     if (std::abs((*p)->pdg_id()) == 13 && (*p)->status() == 1)
0065       mu.push_back(*p);
0066     if (mu.size() > 1)
0067       break;
0068   }
0069
0070   //    std::cout<<" Number of muons "<<mu.size()<<" "<<mu[0]->pdg_id()<<" "<<mu[1]->pdg_id()<<std::endl;
0071
0072   if (mu.size() != 2)
0073     return false;
0074
0075   if (mu[0]->pdg_id() * (mu[1]->pdg_id()) > 0) {
0076     return false;
0077   }
0078
0079   //      std::cout<<" Muons sign accept "<<mu[0]->momentum().perp()<<" "<<mu[1]->momentum().perp()<<std::endl;
0080
0081   if (mu[0]->momentum().perp() < ptMuMin || mu[1]->momentum().perp() < ptMuMin)
0082     return false;
0083
0084   //    std::cout<<" Muons pt accept "<<std::fabs(mu[0]->momentum().eta())<<" "<<std::fabs(mu[1]->momentum().eta())<<std::endl;
0085
0086   if (std::fabs(mu[0]->momentum().eta()) > etaMuMax)
0087     return false;
0088   if (std::fabs(mu[1]->momentum().eta()) > etaMuMax)
0089     return false;
0090
0091   double mmup = mu[0]->generatedMass();
0092   double mmum = mu[1]->generatedMass();
0093   double pxZ = mu[0]->momentum().x() + mu[1]->momentum().x();
0094   double pyZ = mu[0]->momentum().y() + mu[1]->momentum().y();
0095   double pzZ = mu[0]->momentum().z() + mu[1]->momentum().z();
0096
0097   double pmup2 = mu[0]->momentum().x() * mu[0]->momentum().x() + mu[0]->momentum().y() * mu[0]->momentum().y() +
0098                  mu[0]->momentum().z() * mu[0]->momentum().z();
0099   double pmum2 = mu[1]->momentum().x() * mu[1]->momentum().x() + mu[1]->momentum().y() * mu[1]->momentum().y() +
0100                  mu[1]->momentum().z() * mu[1]->momentum().z();
0101   double emup = std::sqrt(pmup2 + mmup * mmup);
0102   double emum = std::sqrt(pmum2 + mmum * mmum);
0103
0104   double massZ = std::sqrt((emup + emum) * (emup + emum) - pxZ * pxZ - pyZ * pyZ - pzZ * pzZ);
0105
0106   //    std::cout<<" Muons eta accept "<<massZ<<std::endl;
0107
0108   if (std::fabs(massZ - m_z) > dm_z)
0109     return false;
0110
0111   //    double ptZ= (mu[0]->momentum() + mu[1]->momentum()).perp();
0112   //    std::cout<<" MassZ accept "<<ptZ<<std::endl;
0113
0114   math::XYZTLorentzVector tot_mom(mu[0]->momentum());
0115   math::XYZTLorentzVector mom2(mu[1]->momentum());
0116   tot_mom += mom2;
0117   //    double ptZ= (mu[0]->momentum() + mu[1]->momentum()).perp();
0118   double ptZ = tot_mom.pt();
0119
0120   if (ptZ > ptZMin && ptZ < ptZMax)
0121     accepted = true;
0122
0123   return accepted;
0124 }
0125
0126 DEFINE_FWK_MODULE(PythiaFilterZJetWithOutBg);