GenFilters/plugins/PythiaFilter.cc

0001
0002 #include "GeneratorInterface/GenFilters/plugins/PythiaFilter.h"
0003 #include "GeneratorInterface/GenFilters/plugins/MCFilterZboostHelper.h"
0004
0005 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0006 #include <iostream>
0007
0008 using namespace edm;
0009 using namespace std;
0010
0011 PythiaFilter::PythiaFilter(const edm::ParameterSet& iConfig)
0012     : token_(consumes<edm::HepMCProduct>(
0013           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0014       particleID(iConfig.getUntrackedParameter("ParticleID", 0)),
0015       minpcut(iConfig.getUntrackedParameter("MinP", 0.)),
0016       maxpcut(iConfig.getUntrackedParameter("MaxP", 10000.)),
0017       minptcut(iConfig.getUntrackedParameter("MinPt", 0.)),
0018       maxptcut(iConfig.getUntrackedParameter("MaxPt", 10000.)),
0019       minetacut(iConfig.getUntrackedParameter("MinEta", -10.)),
0020       maxetacut(iConfig.getUntrackedParameter("MaxEta", 10.)),
0021       minrapcut(iConfig.getUntrackedParameter("MinRapidity", -20.)),
0022       maxrapcut(iConfig.getUntrackedParameter("MaxRapidity", 20.)),
0023       minphicut(iConfig.getUntrackedParameter("MinPhi", -3.5)),
0024       maxphicut(iConfig.getUntrackedParameter("MaxPhi", 3.5)),
0025       status(iConfig.getUntrackedParameter("Status", 0)),
0026       motherID(iConfig.getUntrackedParameter("MotherID", 0)),
0027       processID(iConfig.getUntrackedParameter("ProcessID", 0)),
0028       betaBoost(iConfig.getUntrackedParameter("BetaBoost", 0.)) {
0029   //now do what ever initialization is needed
0030 }
0031
0032 PythiaFilter::~PythiaFilter() {
0033   // do anything here that needs to be done at desctruction time
0034   // (e.g. close files, deallocate resources etc.)
0035 }
0036
0037 //
0038 // member functions
0039 //
0040
0041 // ------------ method called to produce the data  ------------
0042 bool PythiaFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0043   using namespace edm;
0044   bool accepted = false;
0045   Handle<HepMCProduct> evt;
0046   iEvent.getByToken(token_, evt);
0047
0048   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0049
0050   if (processID == 0 || processID == myGenEvent->signal_process_id()) {
0051     for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0052          ++p) {
0053       HepMC::FourVector mom = MCFilterZboostHelper::zboost((*p)->momentum(), betaBoost);
0054       double rapidity = 0.5 * log((mom.e() + mom.pz()) / (mom.e() - mom.pz()));
0055
0056       if (abs((*p)->pdg_id()) == particleID && mom.rho() > minpcut && mom.rho() < maxpcut &&
0057           (*p)->momentum().perp() > minptcut && (*p)->momentum().perp() < maxptcut && mom.eta() > minetacut &&
0058           mom.eta() < maxetacut && rapidity > minrapcut && rapidity < maxrapcut && (*p)->momentum().phi() > minphicut &&
0059           (*p)->momentum().phi() < maxphicut) {
0060         if (status == 0 && motherID == 0) {
0061           accepted = true;
0062         }
0063         if (status != 0 && motherID == 0) {
0064           if ((*p)->status() == status)
0065             accepted = true;
0066         }
0067
0068         HepMC::GenParticle* mother = (*((*p)->production_vertex()->particles_in_const_begin()));
0069
0070         if (status == 0 && motherID != 0) {
0071           if (abs(mother->pdg_id()) == abs(motherID)) {
0072             accepted = true;
0073           }
0074         }
0075         if (status != 0 && motherID != 0) {
0076           if ((*p)->status() == status && abs(mother->pdg_id()) == abs(motherID)) {
0077             accepted = true;
0078           }
0079         }
0080
0081         /*
0082        if (status == 0 && motherID != 0){
0083        if (abs(((*p)->mother())->pdg_id()) == abs(motherID)) {
0084        accepted = true;
0085        }
0086        }
0087        if (status != 0 && motherID != 0){
0088
0089        if ((*p)->status() == status && abs(((*p)->mother())->pdg_id()) == abs(motherID)){
0090        accepted = true;
0091
0092        }
0093        }
0094      */
0095       }
0096     }
0097
0098   } else {
0099     accepted = true;
0100   }
0101
0102   if (accepted) {
0103     return true;
0104   } else {
0105     return false;
0106   }
0107 }