Core/src/GenericDauHepMCFilter.cc

0001
0002 #include "GeneratorInterface/Core/interface/GenericDauHepMCFilter.h"
0003
0004 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0005 #include <iostream>
0006
0007 using namespace edm;
0008 using namespace std;
0009
0010 GenericDauHepMCFilter::GenericDauHepMCFilter(const edm::ParameterSet& iConfig)
0011     : particleID(iConfig.getParameter<int>("ParticleID")),
0012       chargeconju(iConfig.getParameter<bool>("ChargeConjugation")),
0013       ndaughters(iConfig.getParameter<int>("NumberDaughters")),
0014       dauIDs(iConfig.getParameter<std::vector<int> >("DaughterIDs")),
0015       minptcut(iConfig.getParameter<double>("MinPt")),
0016       maxptcut(iConfig.getParameter<double>("MaxPt")),
0017       minetacut(iConfig.getParameter<double>("MinEta")),
0018       maxetacut(iConfig.getParameter<double>("MaxEta")) {}
0019
0020 GenericDauHepMCFilter::~GenericDauHepMCFilter() {}
0021
0022 //
0023 // member functions
0024 //
0025
0026 // ------------ method called to produce the data  ------------
0027 bool GenericDauHepMCFilter::filter(const HepMC::GenEvent* evt) {
0028   for (HepMC::GenEvent::particle_const_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p) {
0029     if ((*p)->pdg_id() != particleID)
0030       continue;
0031
0032     int ndauac = 0;
0033     int ndau = 0;
0034     if ((*p)->end_vertex()) {
0035       for (HepMC::GenVertex::particle_iterator des = (*p)->end_vertex()->particles_begin(HepMC::children);
0036            des != (*p)->end_vertex()->particles_end(HepMC::children);
0037            ++des) {
0038         ++ndau;
0039         //       cout << "   -> Daughter = " << (*des)->pdg_id() << endl;
0040         for (unsigned int i = 0; i < dauIDs.size(); ++i) {
0041           if ((*des)->pdg_id() != dauIDs[i])
0042             continue;
0043           if ((*des)->momentum().perp() > minptcut && (*des)->momentum().perp() < maxptcut &&
0044               (*des)->momentum().eta() > minetacut && (*des)->momentum().eta() < maxetacut) {
0045             ++ndauac;
0046             //           cout << "     -> pT = " << (*des)->momentum().perp() << endl;
0047             break;
0048           }
0049         }
0050       }
0051     }
0052     if (ndau == ndaughters && ndauac == ndaughters) {
0053       return true;
0054     }
0055   }
0056
0057   if (chargeconju) {
0058     for (HepMC::GenEvent::particle_const_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p) {
0059       if ((*p)->pdg_id() != -particleID)
0060         continue;
0061       int ndauac = 0;
0062       int ndau = 0;
0063       if ((*p)->end_vertex()) {
0064         for (HepMC::GenVertex::particle_iterator des = (*p)->end_vertex()->particles_begin(HepMC::children);
0065              des != (*p)->end_vertex()->particles_end(HepMC::children);
0066              ++des) {
0067           ++ndau;
0068           for (unsigned int i = 0; i < dauIDs.size(); ++i) {
0069             bool has_antipart = !(dauIDs[i] == 22 || dauIDs[i] == 23);
0070             int IDanti = has_antipart ? -dauIDs[i] : dauIDs[i];
0071             if ((*des)->pdg_id() != IDanti)
0072               continue;
0073             if ((*des)->momentum().perp() > minptcut && (*des)->momentum().perp() < maxptcut &&
0074                 (*des)->momentum().eta() > minetacut && (*des)->momentum().eta() < maxetacut) {
0075               ++ndauac;
0076               break;
0077             }
0078           }
0079         }
0080       }
0081       if (ndau == ndaughters && ndauac == ndaughters) {
0082         return true;
0083       }
0084     }
0085   }
0086
0087   return false;
0088 }