5 #include "HepMC/PythiaWrapper6_4.h"
14 fVerbose(iConfig.getUntrackedParameter(
"verbose",0)),
15 token_(consumes<edm::
HepMCProduct>(edm::
InputTag(iConfig.getUntrackedParameter(
"moduleLabel",std::
string(
"generator")),
"unsmeared"))),
16 particleID(iConfig.getUntrackedParameter(
"ParticleID", 0)),
17 motherID(iConfig.getUntrackedParameter(
"MotherID", 0)),
18 chargeconju(iConfig.getUntrackedParameter(
"ChargeConjugation",
true)),
19 ndaughters(iConfig.getUntrackedParameter(
"NumberDaughters", 0)),
21 maxptcut(iConfig.getUntrackedParameter(
"MaxPt", 14000.))
27 defdauID.push_back(0);
29 vector<double> defminptcut;
30 defminptcut.push_back(0.);
32 vector<double> defminetacut;
33 defminetacut.push_back(-10.);
35 vector<double> defmaxetacut;
36 defmaxetacut.push_back(10.);
39 cout <<
"----------------------------------------------------------------------" << endl;
40 cout <<
"--- PythiaDauVFilter" << endl;
41 for (
unsigned int i=0;
i<
dauIDs.size(); ++
i) {
47 cout <<
"----------------------------------------------------------------------" << endl;
68 bool accepted =
false;
73 vector<int> vparticles;
75 HepMC::GenEvent *myGenEvent =
new HepMC::GenEvent(*(evt->GetEvent()));
81 for (HepMC::GenEvent::particle_iterator
p = myGenEvent->particles_begin();
p != myGenEvent->particles_end(); ++
p) {
88 for (HepMC::GenVertex::particles_in_const_iterator des = (*p)->production_vertex()->particles_in_const_begin();
89 des != (*p)->production_vertex()->particles_in_const_end();
92 cout <<
"mother: " << (*des)->pdg_id() <<
" pT: " << (*des)->momentum().perp() <<
" eta: " << (*des)->momentum().eta() << endl;
100 if (0 == OK)
continue;
106 cout <<
"found ID: " << (*p)->pdg_id() <<
" pT: " << (*p)->momentum().perp() <<
" eta: " << (*p)->momentum().eta() << endl;
108 vparticles.push_back((*p)->pdg_id());
109 if ((*p)->end_vertex()) {
110 for (HepMC::GenVertex::particle_iterator des=(*p)->end_vertex()->particles_begin(HepMC::children);
111 des != (*p)->end_vertex()->particles_end(HepMC::children);
115 cout <<
"ID: " << (*des)->pdg_id() <<
" pT: " << (*des)->momentum().perp() <<
" eta: " << (*des)->momentum().eta() << endl;
117 for (
unsigned int i=0;
i<
dauIDs.size(); ++
i) {
118 if ((*des)->pdg_id() !=
dauIDs[
i] )
continue ;
120 cout <<
"i = " <<
i <<
" pT = " << (*des)->momentum().perp() <<
" eta = " << (*des)->momentum().eta() << endl;
122 if ((*des)->momentum().perp() >
minptcut[
i] &&
123 (*des)->momentum().perp() <
maxptcut &&
127 vparticles.push_back((*des)->pdg_id());
129 cout <<
" accepted this particle " << (*des)->pdg_id()
130 <<
" pT = " << (*des)->momentum().perp() <<
" eta = " << (*des)->momentum().eta() << endl;
142 cout <<
" accepted this decay: ";
143 for (
unsigned int iv = 0; iv < vparticles.size(); ++iv)
cout << vparticles[iv] <<
" ";
155 for (HepMC::GenEvent::particle_iterator
p = myGenEvent->particles_begin();
156 p != myGenEvent->particles_end(); ++
p) {
163 for (HepMC::GenVertex::particles_in_const_iterator des = (*p)->production_vertex()->particles_in_const_begin();
164 des != (*p)->production_vertex()->particles_in_const_end();
167 cout <<
"mother: " << (*des)->pdg_id() <<
" pT: " << (*des)->momentum().perp() <<
" eta: " << (*des)->momentum().eta() << endl;
175 if (0 == OK)
continue;
178 cout <<
"found ID: " << (*p)->pdg_id() <<
" pT: " << (*p)->momentum().perp() <<
" eta: " << (*p)->momentum().eta() << endl;
180 vparticles.push_back((*p)->pdg_id());
183 if ((*p)->end_vertex()) {
184 for (HepMC::GenVertex::particle_iterator des=(*p)->end_vertex()->particles_begin(HepMC::children);
185 des != (*p)->end_vertex()->particles_end(HepMC::children);
189 cout <<
"ID: " << (*des)->pdg_id() <<
" pT: " << (*des)->momentum().perp() <<
" eta: " << (*des)->momentum().eta() << endl;
191 for (
unsigned int i=0;
i<
dauIDs.size(); ++
i) {
194 int has_antipart = pydat2.kchg[3-1][pythiaCode-1];
195 if (has_antipart == 0) IDanti =
dauIDs[
i];
196 if ((*des)->pdg_id() != IDanti)
continue ;
198 cout <<
"i = " << i <<
" pT = " << (*des)->momentum().perp() <<
" eta = " << (*des)->momentum().eta() << endl;
200 if ((*des)->momentum().perp() >
minptcut[
i] &&
201 (*des)->momentum().perp() <
maxptcut &&
205 vparticles.push_back((*des)->pdg_id());
207 cout <<
" accepted this particle " << (*des)->pdg_id()
208 <<
" pT = " << (*des)->momentum().perp() <<
" eta = " << (*des)->momentum().eta() << endl;
218 cout <<
" accepted this anti-decay: ";
219 for (
unsigned int iv = 0; iv < vparticles.size(); ++iv)
cout << vparticles[iv] <<
" ";
T getUntrackedParameter(std::string const &, T const &) const
bool getByToken(EDGetToken token, Handle< PROD > &result) const
#define DEFINE_FWK_MODULE(type)
virtual bool filter(edm::Event &, const edm::EventSetup &)
std::vector< double > maxetacut
std::vector< int > dauIDs
std::vector< double > minetacut
Abs< T >::type abs(const T &t)
edm::EDGetTokenT< edm::HepMCProduct > token_
std::pair< int, edm::FunctionWithDict > OK
std::vector< double > minptcut
PythiaDauVFilter(const edm::ParameterSet &)