Pythia filter with multiple possible mothers

add HINEcalGenEvtSelector, based on pyquen version

(cherry picked from commit c88da8c)

GeneratorInterface/GenFilters/interface/HINEcalGenEvtSelector.h

@@ -0,0 +1,46 @@
+#ifndef HINECALGENEVTSELECTOR_h
+#define HINECALGENEVTSELECTOR_h
+
+
+// system include files
+#include <memory>
+
+// user include files
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/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 HINEcalGenEvtSelector : public edm::EDFilter {
+   public:
+      explicit HINEcalGenEvtSelector(const edm::ParameterSet&);
+      ~HINEcalGenEvtSelector();
+
+
+      virtual bool filter(edm::Event&, const edm::EventSetup&);
+   private:
+      // ----------member data ---------------------------
+
+      edm::EDGetTokenT<edm::HepMCProduct> token_;
+      std::vector<int> partonId_;
+      std::vector<int> partonStatus_;
+      std::vector<double> partonPt_;
+
+      std::vector<int> particleId_;
+      std::vector<int> particleStatus_;
+      std::vector<double> particlePt_;
+
+      double etaMax_;
+};
+#endif

GeneratorInterface/GenFilters/interface/PythiaFilterMultiMother.h

@@ -0,0 +1,54 @@
+#ifndef PYTHIAFILTERMULTIMOTHER_h
+#define PYTHIAFILTERMULTIMOTHER_h
+
+
+// system include files
+#include <memory>
+
+// user include files
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/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 PythiaFilterMultiMother : public edm::EDFilter {
+   public:
+      explicit PythiaFilterMultiMother(const edm::ParameterSet&);
+      ~PythiaFilterMultiMother();
+
+
+      virtual bool filter(edm::Event&, const edm::EventSetup&);
+   private:
+      // ----------member data ---------------------------
+
+       edm::EDGetTokenT<edm::HepMCProduct> token_;
+       int particleID;
+       double minpcut;
+       double maxpcut;
+       double minptcut;
+       double maxptcut;
+       double minetacut;
+       double maxetacut;
+       double minrapcut;
+       double maxrapcut;
+       double minphicut;
+       double maxphicut;
+
+       double rapidity;
+
+       int status;
+       std::vector<int> motherIDs;
+       int processID;
+};
+#endif

GeneratorInterface/GenFilters/src/HINEcalGenEvtSelector.cc

@@ -0,0 +1,85 @@
+
+#include "GeneratorInterface/GenFilters/interface/HINEcalGenEvtSelector.h"
+
+#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+#include <iostream>
+
+using namespace edm;
+using namespace std;
+
+bool selectParticle(HepMC::GenParticle* par, int status, int pdg /*Absolute*/, double ptMin, double etaMax){
+  return (par->status() == status && abs(par->pdg_id()) == pdg && par->momentum().perp() > ptMin && fabs(par->momentum().eta()) < etaMax);
+}
+
+HINEcalGenEvtSelector::HINEcalGenEvtSelector(const edm::ParameterSet& iConfig) :
+  token_(consumes<edm::HepMCProduct>(iConfig.getUntrackedParameter("moduleLabel",std::string("generator"))))
+{
+  //now do what ever initialization is needed
+  partonId_ = iConfig.getParameter<vector<int> >("partons");
+  partonStatus_ = iConfig.getParameter<vector<int> >("partonStatus");
+  partonPt_ = iConfig.getParameter<vector<double> >("partonPt");
+
+  particleId_ = iConfig.getParameter<vector<int> >("particles");
+  particleStatus_ = iConfig.getParameter<vector<int> >("particleStatus");
+  particlePt_ = iConfig.getParameter<vector<double> >("particlePt");
+
+  etaMax_ = iConfig.getParameter<double>("etaMax");
+
+  int id = partonId_.size();
+  int st = partonStatus_.size();
+  int pt = partonPt_.size();
+
+  if(partonId_.size() != partonStatus_.size() || partonId_.size() != partonPt_.size()){
+    throw edm::Exception(edm::errors::LogicError)<<id<<st<<pt<<endl;
+  }
+
+  id = particleId_.size();
+  st = particleStatus_.size();
+  pt = particlePt_.size();
+
+  if(particleId_.size() != particleStatus_.size() || particleId_.size() != particlePt_.size()){
+    throw edm::Exception(edm::errors::LogicError)<<id<<st<<pt<<endl;
+  }
+}
+
+
+HINEcalGenEvtSelector::~HINEcalGenEvtSelector()
+{
+
+   // 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 HINEcalGenEvtSelector::filter(edm::Event& iEvent, const edm::EventSetup& iSetup)
+{
+   using namespace edm;
+   Handle<HepMCProduct> evt;
+   iEvent.getByToken(token_, evt);
+
+   const HepMC::GenEvent * myGenEvent = evt->GetEvent();
+   HepMC::GenEvent::particle_const_iterator begin = myGenEvent->particles_begin();
+   HepMC::GenEvent::particle_const_iterator end = myGenEvent->particles_end();
+
+   bool foundParticle = false;
+   bool foundParton = false;
+
+   HepMC::GenEvent::particle_const_iterator it = begin;
+   while((!foundParton || !foundParticle) && it != end){
+     for(unsigned i = 0; i < partonId_.size(); ++i){
+       if(selectParticle(*it, partonStatus_[i], partonId_[i], partonPt_[i], etaMax_)) foundParton = true;
+     }
+     for(unsigned i = 0; i < particleId_.size(); ++i){
+       if(selectParticle(*it, particleStatus_[i], particleId_[i], particlePt_[i], etaMax_)) foundParticle = true;
+     }
+     ++it;
+   }
+
+   return (foundParton && foundParticle);
+}

GeneratorInterface/GenFilters/src/PythiaFilterMultiMother.cc

@@ -0,0 +1,124 @@
+
+#include "GeneratorInterface/GenFilters/interface/PythiaFilterMultiMother.h"
+
+#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+#include <iostream>
+
+using namespace edm;
+using namespace std;
+
+
+PythiaFilterMultiMother::PythiaFilterMultiMother(const edm::ParameterSet& iConfig) :
+token_(consumes<edm::HepMCProduct>(iConfig.getUntrackedParameter("moduleLabel",std::string("generator")))),
+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)),
+status(iConfig.getUntrackedParameter("Status", 0)),
+motherIDs(iConfig.getUntrackedParameter("MotherIDs", std::vector<int>{0})),
+processID(iConfig.getUntrackedParameter("ProcessID", 0))
+{
+   //now do what ever initialization is needed
+
+}
+
+
+PythiaFilterMultiMother::~PythiaFilterMultiMother()
+{
+
+   // 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 PythiaFilterMultiMother::filter(edm::Event& iEvent, const edm::EventSetup& iSetup)
+{
+   using namespace edm;
+   bool accepted = false;
+   Handle<HepMCProduct> evt;
+   iEvent.getByToken(token_, evt);
+
+   const HepMC::GenEvent * myGenEvent = evt->GetEvent();
+
+   if(processID == 0 || processID == myGenEvent->signal_process_id()) {
+
+     for ( HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();
+	   p != myGenEvent->particles_end(); ++p ) {
+
+       rapidity = 0.5*log( ((*p)->momentum().e()+(*p)->momentum().pz()) / ((*p)->momentum().e()-(*p)->momentum().pz()) );
+
+       if ( abs((*p)->pdg_id()) == particleID
+	    && (*p)->momentum().rho() > minpcut
+	    && (*p)->momentum().rho() < maxpcut
+	    && (*p)->momentum().perp() > minptcut
+	    && (*p)->momentum().perp() < maxptcut
+	    && (*p)->momentum().eta() > minetacut
+	    && (*p)->momentum().eta() < maxetacut
+	    && rapidity > minrapcut
+	    && rapidity < maxrapcut
+	    && (*p)->momentum().phi() > minphicut
+	    && (*p)->momentum().phi() < maxphicut ) {
+
+
+	 for(std::vector<int>::const_iterator motherID = motherIDs.begin(); motherID != motherIDs.end(); ++motherID) {
+	   if (status == 0 && *motherID == 0){
+	     accepted = true;
+	   }
+	   if (status != 0 && *motherID == 0){
+	     if ((*p)->status() == status)
+	       accepted = true;
+	   }
+
+	   HepMC::GenParticle* mother = (*((*p)->production_vertex()->particles_in_const_begin()));
+
+	   if (status == 0 && *motherID != 0){
+	     if (abs(mother->pdg_id()) == abs(*motherID)) {
+	       accepted = true;
+	     }
+	   }
+	   if (status != 0 && *motherID != 0){
+
+	     if ((*p)->status() == status && abs(mother->pdg_id()) == abs(*motherID)){
+	       accepted = true;
+
+	     }
+	   }
+	 }
+
+	 /*
+	   if (status == 0 && motherID != 0){
+	   if (abs(((*p)->mother())->pdg_id()) == abs(motherID)) {
+	   accepted = true;
+	   }
+	   }
+	   if (status != 0 && motherID != 0){
+
+	   if ((*p)->status() == status && abs(((*p)->mother())->pdg_id()) == abs(motherID)){
+	   accepted = true;
+
+	   }
+	   }
+	 */
+
+       }
+     }
+
+   } else { accepted = true; }
+
+   if (accepted){
+   return true; } else {return false;}
+
+}

GeneratorInterface/GenFilters/src/SealModule.cc

@@ -1,6 +1,8 @@
 #include "FWCore/Framework/interface/InputSourceMacros.h"
 #include "FWCore/Framework/interface/MakerMacros.h"
 #include "GeneratorInterface/GenFilters/interface/PythiaFilter.h"
+#include "GeneratorInterface/GenFilters/interface/PythiaFilterMultiMother.h"
+#include "GeneratorInterface/GenFilters/interface/HINEcalGenEvtSelector.h"
 #include "GeneratorInterface/GenFilters/interface/PythiaDauFilter.h"
 #include "GeneratorInterface/GenFilters/interface/PythiaFilterGammaJet.h"
 #include "GeneratorInterface/GenFilters/interface/PythiaFilterGammaGamma.h"
@@ -47,6 +49,8 @@
   using cms::BHFilter;
   DEFINE_FWK_MODULE(LQGenFilter);
   DEFINE_FWK_MODULE(PythiaFilter);
+  DEFINE_FWK_MODULE(PythiaFilterMultiMother);
+  DEFINE_FWK_MODULE(HINEcalGenEvtSelector);
   DEFINE_FWK_MODULE(PythiaDauFilter);
   DEFINE_FWK_MODULE(PythiaFilterGammaJet);
   DEFINE_FWK_MODULE(PythiaFilterGammaGamma);