HiGenCommon/src/HadronDecayGenEvtSelector.cc

0001 #include <iostream>
0002 #include "GeneratorInterface/HiGenCommon/interface/HadronDecayGenEvtSelector.h"
0003 #include "FWCore/Utilities/interface/EDMException.h"
0004 using namespace std;
0005
0006 HadronDecayGenEvtSelector::HadronDecayGenEvtSelector(const edm::ParameterSet& pset) : BaseHiGenEvtSelector(pset) {
0007   hadronId_ = pset.getParameter<vector<int> >("hadrons");
0008   hadronStatus_ = pset.getParameter<vector<int> >("hadronStatus");
0009   hadronEtaMax_ = pset.getParameter<vector<double> >("hadronEtaMax");
0010   hadronEtaMin_ = pset.getParameter<vector<double> >("hadronEtaMin");
0011   hadronPMin_ = pset.getParameter<vector<double> >("hadronPMin");
0012   hadronPtMax_ = pset.getParameter<vector<double> >("hadronPtMax");
0013   hadronPtMin_ = pset.getParameter<vector<double> >("hadronPtMin");
0014
0015   decayId_ = pset.getParameter<int>("decays");
0016   decayStatus_ = pset.getParameter<int>("decayStatus");
0017   decayEtaMax_ = pset.getParameter<double>("decayEtaMax");
0018   decayEtaMin_ = pset.getParameter<double>("decayEtaMin");
0019   decayPMin_ = pset.getParameter<double>("decayPMin");
0020   decayPtMax_ = pset.getParameter<double>("decayPtMax");
0021   decayPtMin_ = pset.getParameter<double>("decayPtMin");
0022   decayNtrig_ = pset.getParameter<int>("decayNtrig");
0023
0024   int id = hadronId_.size();
0025   int st = hadronStatus_.size();
0026   int etamax = hadronEtaMax_.size();
0027   int etamin = hadronEtaMin_.size();
0028   int pmin = hadronPMin_.size();
0029   int ptmax = hadronPtMax_.size();
0030   int ptmin = hadronPtMin_.size();
0031
0032   if (id != st || id != etamax || id != etamin || id != ptmax || id != ptmin || id != pmin) {
0033     throw edm::Exception(edm::errors::LogicError)
0034         << "Hadron selection parameters: " << id << st << etamax << etamin << pmin << ptmax << ptmin << endl;
0035   }
0036 }
0037
0038 //____________________________________________________________________________________________
0039 bool HadronDecayGenEvtSelector::filter(HepMC::GenEvent* evt) {
0040   // loop over HepMC event, and search for  products of interest
0041
0042   HepMC::GenEvent::particle_const_iterator begin = evt->particles_begin();
0043   HepMC::GenEvent::particle_const_iterator end = evt->particles_end();
0044
0045   bool foundHadron = false;
0046   bool foundDecay = false;
0047
0048   int foundtrig = 0;
0049   HepMC::GenEvent::particle_const_iterator it = begin;
0050   while ((!foundHadron || !foundDecay) && it != end) {
0051     for (unsigned i = 0; i < hadronId_.size(); ++i) {
0052       if (selectParticle(*it,
0053                          hadronStatus_[i],
0054                          hadronId_[i],
0055                          hadronEtaMax_[i],
0056                          hadronEtaMin_[i],
0057                          hadronPMin_[i],
0058                          hadronPtMax_[i],
0059                          hadronPtMin_[i]))
0060         foundHadron = true;
0061     }
0062
0063     if (selectParticle(*it, decayStatus_, decayId_, decayEtaMax_, decayEtaMin_, decayPMin_, decayPtMax_, decayPtMin_))
0064       foundtrig++;
0065     if (decayNtrig_ == foundtrig)
0066       foundDecay = true;
0067
0068     ++it;
0069   }
0070
0071   return (foundHadron && foundDecay);
0072 }
0073
0074 //____________________________________________________________________________________________