GenFilters/plugins/MCZll.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    MCZll
0004 // Class:      MCZll
0005 //
0006 /*
0007
0008  Description: filter events based on the Pythia ProcessID and the Pt_hat
0009
0010  Implementation: inherits from generic EDFilter
0011
0012 */
0013 //
0014 // Original Author:  Paolo Meridiani
0015 //
0016 //
0017
0018 #include "DataFormats/Common/interface/Handle.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/Frameworkfwd.h"
0021 #include "FWCore/Framework/interface/global/EDFilter.h"
0022 #include "FWCore/Framework/interface/MakerMacros.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "FWCore/Utilities/interface/EDGetToken.h"
0026 #include "FWCore/Utilities/interface/InputTag.h"
0027 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0028
0029 #include <cmath>
0030 #include <cstdlib>
0031 #include <memory>
0032 #include <sstream>
0033 #include <string>
0034 #include <utility>
0035
0036 class MCZll : public edm::global::EDFilter<> {
0037 public:
0038   explicit MCZll(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   int leptonFlavour_;
0045   double leptonPtMin_;
0046   double leptonPtMax_;
0047   double leptonEtaMin_;
0048   double leptonEtaMax_;
0049   std::pair<double, double> zMassRange_;
0050   bool filter_;
0051 };
0052
0053 using namespace edm;
0054 using namespace std;
0055
0056 MCZll::MCZll(const edm::ParameterSet& iConfig)
0057     : token_(consumes<edm::HepMCProduct>(
0058           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))) {
0059   leptonFlavour_ = iConfig.getUntrackedParameter<int>("leptonFlavour", 11);
0060   leptonPtMin_ = iConfig.getUntrackedParameter<double>("leptonPtMin", 5.);
0061   leptonPtMax_ = iConfig.getUntrackedParameter<double>("leptonPtMax", 99999.);
0062   leptonEtaMin_ = iConfig.getUntrackedParameter<double>("leptonEtaMin", 0.);
0063   leptonEtaMax_ = iConfig.getUntrackedParameter<double>("leptonEtaMax", 2.7);
0064   zMassRange_.first = iConfig.getUntrackedParameter<double>("zMassMin", 60.);
0065   zMassRange_.second = iConfig.getUntrackedParameter<double>("zMassMax", 120.);
0066   filter_ = iConfig.getUntrackedParameter<bool>("filter", true);
0067   std::ostringstream str;
0068   str << "=========================================================\n"
0069       << "Filter MCZll being constructed with parameters: "
0070       << "\nleptonFlavour " << leptonFlavour_ << "\nleptonPtMin " << leptonPtMin_ << "\nleptonPtMax " << leptonPtMax_
0071       << "\nleptonEtaMin " << leptonEtaMin_ << "\nleptonEtaMax " << leptonEtaMax_ << "\nzMassMin " << zMassRange_.first
0072       << "\nzMassMax " << zMassRange_.second << "\n=========================================================";
0073   edm::LogVerbatim("MCZllInfo") << str.str();
0074   if (filter_)
0075     produces<HepMCProduct>();
0076 }
0077
0078 bool MCZll::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0079   std::unique_ptr<HepMCProduct> bare_product(new HepMCProduct());
0080
0081   bool accepted = false;
0082   Handle<HepMCProduct> evt;
0083   iEvent.getByToken(token_, evt);
0084   HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
0085   HepMC::GenEvent* zEvent = new HepMC::GenEvent();
0086
0087   if (myGenEvent->signal_process_id() != 1) {
0088     delete myGenEvent;
0089     delete zEvent;
0090     return false;
0091   }
0092
0093   //found a prompt Z
0094
0095   for (HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0096     if (!accepted && ((*p)->pdg_id() == 23) && (*p)->status() == 3) {
0097       accepted = true;
0098       HepMC::GenVertex* zVertex = new HepMC::GenVertex();
0099       HepMC::GenParticle* myZ = new HepMC::GenParticle(*(*p));
0100       zVertex->add_particle_in(myZ);
0101       //      std::cout << (*p)->momentum().invariantMass() << std::endl;
0102       if ((*p)->momentum().m() < zMassRange_.first || (*p)->momentum().m() > zMassRange_.second)
0103         accepted = false;
0104       std::vector<HepMC::GenParticle*> children;
0105       HepMC::GenVertex* outVertex = (*p)->end_vertex();
0106       for (HepMC::GenVertex::particles_out_const_iterator iter = outVertex->particles_out_const_begin();
0107            iter != outVertex->particles_out_const_end();
0108            iter++)
0109         children.push_back(*iter);
0110       std::vector<HepMC::GenParticle*>::const_iterator aDaughter;
0111       for (aDaughter = children.begin(); aDaughter != children.end(); aDaughter++) {
0112         HepMC::GenParticle* myDa = new HepMC::GenParticle(*(*aDaughter));
0113         zVertex->add_particle_out(myDa);
0114         if ((*aDaughter)->status() == 2)
0115           continue;
0116         //        (*aDaughter)->print();
0117
0118         if (!(abs((*aDaughter)->pdg_id()) == abs(leptonFlavour_)))
0119           accepted = false;
0120         //      std::cout << (*aDaughter)->momentum().perp() << " " << (*aDaughter)->momentum().eta() << std::endl;
0121         if ((*aDaughter)->momentum().perp() < leptonPtMin_)
0122           accepted = false;
0123         if ((*aDaughter)->momentum().perp() > leptonPtMax_)
0124           accepted = false;
0125         if (fabs((*aDaughter)->momentum().eta()) > leptonEtaMax_)
0126           accepted = false;
0127         if (fabs((*aDaughter)->momentum().eta()) < leptonEtaMin_)
0128           accepted = false;
0129       }
0130       zEvent->add_vertex(zVertex);
0131       if (accepted)
0132         break;
0133     }
0134   }
0135
0136   if (accepted) {
0137     if (zEvent)
0138       bare_product->addHepMCData(zEvent);
0139     if (filter_)
0140       iEvent.put(std::move(bare_product));
0141     LogDebug("MCZll") << "Event " << iEvent.id().event() << " accepted" << std::endl;
0142     delete myGenEvent;
0143     return true;
0144   }
0145
0146   delete myGenEvent;
0147   delete zEvent;
0148   return false;
0149 }
0150
0151 DEFINE_FWK_MODULE(MCZll);