Merge pull request #33432 from mgratti/from-CMSSW_10_2_X_2021-04-11-0000

[backport] Added new filter, thought for a B -> l HNL X analysis

GeneratorInterface/GenFilters/interface/PythiaFilterMotherSister.h

@@ -0,0 +1,74 @@
+#ifndef PYTHIAFILTERMOTHERSISTER_h
+#define PYTHIAFILTERMOTHERSISTER_h
+// -*- C++ -*-
+//
+// Package:    PythiaFilterMotherSister
+// Class:      PythiaFilterMotherSister
+// 
+/**\class PythiaFilterMotherSister PythiaFilterMotherSister.cc IOMC/PythiaFilterMotherSister/src/PythiaFilterMotherSister.cc
+
+ Description: A filter to identify a particle with given id and kinematic 
+                                                && given mother id (multiple mothers possible)
+                                                && given id and 3d displacement of one among mother's daughters
+
+ Implementation:
+     Inspired by PythiaFilterMultiMother.cc
+*/
+//
+// 
+//         
+//
+//
+
+
+// system include files
+#include <memory>
+
+// user include files
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/global/EDFilter.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+//
+// class decleration
+//
+namespace edm {
+  class HepMCProduct;
+}
+
+class PythiaFilterMotherSister : public edm::global::EDFilter<> {
+   public:
+      explicit PythiaFilterMotherSister(const edm::ParameterSet&);
+      ~PythiaFilterMotherSister() override;
+
+
+      bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
+   private:
+      // ----------member data ---------------------------
+
+       const edm::EDGetTokenT<edm::HepMCProduct> token_;
+       const int particleID;
+       const double minpcut;
+       const double maxpcut;
+       const double minptcut;
+       const double maxptcut;
+       const double minetacut;
+       const double maxetacut;
+       const double minrapcut;
+       const double maxrapcut;
+       const double minphicut;
+       const double maxphicut;
+
+       //const int status; 
+       std::vector<int> motherIDs;
+       const int sisterID;   
+       //const int processID;    
+
+       const double betaBoost;
+       const double maxSisDisplacement;
+};
+#endif

GeneratorInterface/GenFilters/src/PythiaFilterMotherSister.cc

@@ -0,0 +1,116 @@
+
+#include "GeneratorInterface/GenFilters/interface/PythiaFilterMotherSister.h"
+#include "GeneratorInterface/GenFilters/interface/MCFilterZboostHelper.h"
+
+#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+#include <iostream>
+
+using namespace edm;
+using namespace std;
+
+
+PythiaFilterMotherSister::PythiaFilterMotherSister(const edm::ParameterSet& iConfig) :
+token_(consumes<edm::HepMCProduct>(edm::InputTag(iConfig.getUntrackedParameter("moduleLabel",std::string("generator")),"unsmeared"))),
+particleID(iConfig.getUntrackedParameter("ParticleID", 0)),
+minpcut(iConfig.getUntrackedParameter("MinP", 0.)),
+maxpcut(iConfig.getUntrackedParameter("MaxP", 10000.)),
+minptcut(iConfig.getUntrackedParameter("MinPt", 0.)),
+maxptcut(iConfig.getUntrackedParameter("MaxPt", 10000.)),
+minetacut(iConfig.getUntrackedParameter("MinEta", -10.)),
+maxetacut(iConfig.getUntrackedParameter("MaxEta", 10.)),
+minrapcut(iConfig.getUntrackedParameter("MinRapidity", -20.)),
+maxrapcut(iConfig.getUntrackedParameter("MaxRapidity", 20.)),
+minphicut(iConfig.getUntrackedParameter("MinPhi", -3.5)),
+maxphicut(iConfig.getUntrackedParameter("MaxPhi", 3.5)),
+motherIDs(iConfig.getUntrackedParameter("MotherIDs", std::vector<int>{0})),
+sisterID(iConfig.getUntrackedParameter("SisterID", 0)),
+betaBoost(iConfig.getUntrackedParameter("BetaBoost",0.)),
+maxSisDisplacement(iConfig.getUntrackedParameter("MaxSisterDisplacement", -1.))
+{
+   //now do what ever initialization is needed
+
+}
+
+
+PythiaFilterMotherSister::~PythiaFilterMotherSister()
+{
+
+   // do anything here that needs to be done at desctruction time
+   // (e.g. close files, deallocate resources etc.)
+
+}
+
+
+//
+// member functions
+//
+
+// ------------ method called to produce the data  ------------
+bool PythiaFilterMotherSister::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const
+{
+   using namespace edm;
+   Handle<HepMCProduct> evt;
+   iEvent.getByToken(token_, evt);
+
+   const HepMC::GenEvent * myGenEvent = evt->GetEvent();
+
+   for ( HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();
+         p != myGenEvent->particles_end(); ++p ) {
+     HepMC::FourVector mom = MCFilterZboostHelper::zboost((*p)->momentum(),betaBoost);
+     double rapidity = 0.5*log( (mom.e()+mom.pz()) / (mom.e()-mom.pz()) );
+
+     if ( abs((*p)->pdg_id()) == particleID 
+          && mom.rho() > minpcut 
+          && mom.rho() < maxpcut
+          && (*p)->momentum().perp() > minptcut 
+          && (*p)->momentum().perp() < maxptcut
+          && mom.eta() > minetacut
+          && mom.eta() < maxetacut 
+          && rapidity > minrapcut
+          && rapidity < maxrapcut 
+          && (*p)->momentum().phi() > minphicut
+          && (*p)->momentum().phi() < maxphicut ) 
+     {
+
+
+       HepMC::GenParticle* mother = (*((*p)->production_vertex()->particles_in_const_begin()));
+
+       // check various possible mothers
+       for(auto motherID : motherIDs){
+
+         if(abs(mother->pdg_id()) == abs(motherID)){
+
+           // loop over its daughters
+           for ( HepMC::GenVertex::particle_iterator dau = mother->end_vertex()->particles_begin(HepMC::children); 
+                                                     dau != mother->end_vertex()->particles_end(HepMC::children);
+                                                     ++dau ) {
+             // find the daugther you're interested in
+             if(abs((*dau)->pdg_id()) == abs(sisterID)) {
+
+               // calculate displacement of the sister particle, from production to decay
+               HepMC::GenVertex* v1   = (*dau)->production_vertex();
+               HepMC::GenVertex* v2   = (*dau)->end_vertex();
+
+               double lx12 = v1->position().x() - v2->position().x();
+               double ly12 = v1->position().y() - v2->position().y();
+               double lz12 = v1->position().z() - v2->position().z();
+               double lxyz12 = sqrt( lx12*lx12 + ly12*ly12 + lz12*lz12 );
+
+               if(maxSisDisplacement!= -1){
+                 if(lxyz12 < maxSisDisplacement){
+                   return true; 
+                 }
+               } else {
+                  return true; 
+               }
+
+             } 
+           } 
+         } 
+       } 
+     } 
+   } 
+
+   return false;
+
+}

GeneratorInterface/GenFilters/src/SealModule.cc

@@ -3,6 +3,7 @@
 #include "GeneratorInterface/GenFilters/interface/PythiaFilter.h"
 #include "GeneratorInterface/GenFilters/interface/PythiaFilterHT.h"
 #include "GeneratorInterface/GenFilters/interface/PythiaFilterMultiMother.h"
+#include "GeneratorInterface/GenFilters/interface/PythiaFilterMotherSister.h"
 #include "GeneratorInterface/GenFilters/interface/PythiaDauFilter.h"
 #include "GeneratorInterface/GenFilters/interface/PythiaProbeFilter.h"
 #include "GeneratorInterface/GenFilters/interface/PythiaFilterGammaJet.h"
@@ -56,6 +57,7 @@
   DEFINE_FWK_MODULE(PythiaFilter);
   DEFINE_FWK_MODULE(PythiaFilterHT);
   DEFINE_FWK_MODULE(PythiaFilterMultiMother);
+  DEFINE_FWK_MODULE(PythiaFilterMotherSister);
   DEFINE_FWK_MODULE(PythiaDauFilter);
   DEFINE_FWK_MODULE(PythiaProbeFilter);
   DEFINE_FWK_MODULE(PythiaFilterGammaJet);