GenFilters/plugins/MCDisplacementFilter.cc

0001
0002 // -*- C++ -*-
0003 //
0004 // Package:    MCDisplacementFilter
0005 // Class:      MCDisplacementFilter
0006 //
0007
0008 // system include files
0009 //#include <memory>
0010 //#include <iostream>
0011 #include <vector>
0012
0013 // user include files
0014 #include "FWCore/Framework/interface/Frameworkfwd.h"
0015 #include "FWCore/Framework/interface/global/EDFilter.h"
0016 #include "FWCore/Utilities/interface/InputTag.h"
0017
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/MakerMacros.h"
0020
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0023
0024 using namespace edm;
0025 using namespace std;
0026
0027 //Filter particles based on their minimum and/or maximum displacement on the transverse plane and optionally on their pdgIds
0028 //To run independently of pdgId, do not insert the particleIDs entry in filter declaration
0029
0030 // class decleration
0031 //
0032 namespace edm {
0033   class HepMCProduct;
0034   class ConfigurationDescriptions;
0035 }  // namespace edm
0036
0037 class MCDisplacementFilter : public edm::global::EDFilter<> {
0038 public:
0039   explicit MCDisplacementFilter(const edm::ParameterSet&);
0040   ~MCDisplacementFilter() override = default;
0041
0042   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0043
0044   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0045
0046 private:
0047   // ----------member data ---------------------------
0048   edm::InputTag hepMCProductTag_;
0049   const edm::EDGetTokenT<edm::HepMCProduct> token_;
0050   // To chose on which pdgIds the filter is applied - if ParticleIDs.at(0)==0 runs on all particles
0051   std::vector<int> particleIDs_;
0052
0053   const float theUpperCut_;  // Maximum displacement accepted
0054   const float theLowerCut_;  // Minimum displacement accepted
0055 };
0056
0057 //methods implementation
0058 //
0059 //Class initialization
0060 MCDisplacementFilter::MCDisplacementFilter(const edm::ParameterSet& iConfig)
0061     : hepMCProductTag_(iConfig.getParameter<edm::InputTag>("hepMCProductTag")),
0062       token_(consumes<edm::HepMCProduct>(hepMCProductTag_)),
0063       particleIDs_(iConfig.getParameter<std::vector<int>>("ParticleIDs")),
0064       theUpperCut_(iConfig.getParameter<double>("LengMax")),
0065       theLowerCut_(iConfig.getParameter<double>("LengMin")) {}
0066
0067 //Filter description
0068 void MCDisplacementFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0069   edm::ParameterSetDescription desc;
0070   desc.add<edm::InputTag>("hepMCProductTag", edm::InputTag("generator", "unsmeared"));
0071   desc.add<std::vector<int>>("ParticleIDs", std::vector<int>{0});
0072   desc.add<double>("LengMax", -1.);
0073   desc.add<double>("LengMin", -1.);
0074   descriptions.addDefault(desc);
0075 }
0076
0077 // ------------ method called to skim the data  ------------
0078 bool MCDisplacementFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0079   using namespace edm;
0080
0081   Handle<HepMCProduct> evt;
0082
0083   iEvent.getByToken(token_, evt);
0084
0085   bool pass = false;
0086   bool matchedID = true;
0087
0088   const float theUpperCut2 = theUpperCut_ * theUpperCut_;
0089   const float theLowerCut2 = theLowerCut_ * theLowerCut_;
0090
0091   const HepMC::GenEvent* generated_event = evt->GetEvent();
0092   HepMC::GenEvent::particle_const_iterator p;
0093
0094   for (p = generated_event->particles_begin(); p != generated_event->particles_end(); p++) {
0095     //matchedID might be moved to false to true for a particle in the event, it needs to be resetted everytime
0096     if (particleIDs_.at(0) != 0)
0097       matchedID = false;
0098
0099     //if a list of pdgId is provided, loop only on particles with those pdgId
0100     for (unsigned int idx = 0; idx < particleIDs_.size(); idx++) {
0101       if (abs((*p)->pdg_id()) == abs(particleIDs_.at(idx))) {  //compares absolute values of pdgIds
0102         matchedID = true;
0103         break;
0104       }
0105     }
0106
0107     if (matchedID) {
0108       if (theLowerCut_ <= 0. && theUpperCut_ <= 0.) {
0109         pass = true;
0110         break;
0111       }
0112       if (((*p)->production_vertex() != nullptr) && ((*p)->end_vertex() != nullptr)) {
0113         float dist2 = (((*p)->production_vertex())->position().x() - ((*p)->end_vertex())->position().x()) *
0114                           (((*p)->production_vertex())->position().x() - ((*p)->end_vertex())->position().x()) +
0115                       (((*p)->production_vertex())->position().y() - ((*p)->end_vertex())->position().y()) *
0116                           (((*p)->production_vertex())->position().y() - ((*p)->end_vertex())->position().y());
0117         // lower/upper cuts can be also <= 0 - prompt particle needs to be accepted in that case
0118         if ((theLowerCut_ <= 0. || dist2 >= theLowerCut2) && (theUpperCut_ <= 0. || dist2 < theUpperCut2)) {
0119           pass = true;
0120           break;
0121         }
0122       }
0123       if (((*p)->production_vertex() == nullptr) && (((*p)->end_vertex() != nullptr))) {
0124         // lower/upper cuts can be also 0 - prompt particle needs to be accepted in that case
0125         float distEndVert = (*p)->end_vertex()->position().perp();
0126         if ((theLowerCut_ <= 0. || distEndVert >= theLowerCut_) && (theUpperCut_ <= 0. || distEndVert < theUpperCut_)) {
0127           pass = true;
0128           break;
0129         }
0130       }
0131     }
0132   }
0133
0134   return pass;
0135 }
0136 DEFINE_FWK_MODULE(MCDisplacementFilter);