HiGenCommon/src/PartonHadronDecayGenEvtSelector.cc

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