Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:13:42 

0001 // -*- C++ -*-
0002 //
0003 // Package:    PythiaMomDauFilter
0004 // Class:      PythiaMomDauFilter
0005 //
0006 /**\class PythiaMomDauFilter PythiaMomDauFilter.cc
0007 
0008  Description: Filter events using MotherId and ChildrenIds infos
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Daniele Pedrini
0015 //         Created:  Oct 27 2015
0016 // Fixed : Ta-Wei Wang, Dec 11 2015
0017 // $Id: PythiaMomDauFilter.h,v 1.1 2015/10/27  pedrini Exp $
0018 //
0019 
0020 #include "DataFormats/Common/interface/Handle.h"
0021 #include "FWCore/Concurrency/interface/SharedResourceNames.h"
0022 #include "FWCore/Framework/interface/Event.h"
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/one/EDFilter.h"
0025 #include "FWCore/Framework/interface/MakerMacros.h"
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027 #include "FWCore/Utilities/interface/EDGetToken.h"
0028 #include "FWCore/Utilities/interface/InputTag.h"
0029 #include "GeneratorInterface/GenFilters/plugins/MCFilterZboostHelper.h"
0030 #include "HepMC/PythiaWrapper6_4.h"
0031 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0032 
0033 #include <string>
0034 #include <vector>
0035 
0036 #define PYCOMP pycomp_
0037 extern "C" {
0038 int PYCOMP(int& ip);
0039 }
0040 
0041 // Must be a "one" because it uses a Pythia 6 FORTRAN function
0042 class PythiaMomDauFilter : public edm::one::EDFilter<edm::one::SharedResources> {
0043 public:
0044   explicit PythiaMomDauFilter(const edm::ParameterSet&);
0045 
0046   bool filter(edm::Event&, const edm::EventSetup&) override;
0047 
0048 private:
0049   // ----------member data ---------------------------
0050 
0051   edm::EDGetTokenT<edm::HepMCProduct> label_;
0052   std::vector<int> dauIDs;
0053   std::vector<int> desIDs;
0054   int particleID;
0055   int daughterID;
0056   bool chargeconju;
0057   int ndaughters;
0058   int ndescendants;
0059   double minptcut;
0060   double maxptcut;
0061   double minetacut;
0062   double maxetacut;
0063   double mom_minptcut;
0064   double mom_maxptcut;
0065   double mom_minetacut;
0066   double mom_maxetacut;
0067   double betaBoost;
0068 };
0069 
0070 using namespace std;
0071 
0072 PythiaMomDauFilter::PythiaMomDauFilter(const edm::ParameterSet& iConfig)
0073     : label_(consumes<edm::HepMCProduct>(
0074           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0075       particleID(iConfig.getUntrackedParameter<int>("ParticleID", 0)),
0076       daughterID(iConfig.getUntrackedParameter<int>("DaughterID", 0)),
0077       chargeconju(iConfig.getUntrackedParameter<bool>("ChargeConjugation", true)),
0078       ndaughters(iConfig.getUntrackedParameter<int>("NumberDaughters", 0)),
0079       ndescendants(iConfig.getUntrackedParameter<int>("NumberDescendants", 0)),
0080       minptcut(iConfig.getUntrackedParameter<double>("MinPt", 0.)),
0081       maxptcut(iConfig.getUntrackedParameter<double>("MaxPt", 14000.)),
0082       minetacut(iConfig.getUntrackedParameter<double>("MinEta", -10.)),
0083       maxetacut(iConfig.getUntrackedParameter<double>("MaxEta", 10.)),
0084       mom_minptcut(iConfig.getUntrackedParameter<double>("MomMinPt", 0.)),
0085       mom_maxptcut(iConfig.getUntrackedParameter<double>("MomMaxPt", 14000.)),
0086       mom_minetacut(iConfig.getUntrackedParameter<double>("MomMinEta", -10.)),
0087       mom_maxetacut(iConfig.getUntrackedParameter<double>("MomMaxEta", 10.)),
0088       betaBoost(iConfig.getUntrackedParameter("BetaBoost", 0.)) {
0089   usesResource(edm::SharedResourceNames::kPythia6);
0090 
0091   //now do what ever initialization is needed
0092   vector<int> defdauID;
0093   defdauID.push_back(0);
0094   vector<int> defdesID;
0095   defdesID.push_back(0);
0096   dauIDs = iConfig.getUntrackedParameter<vector<int> >("DaughterIDs", defdauID);
0097   desIDs = iConfig.getUntrackedParameter<vector<int> >("DescendantsIDs", defdesID);
0098 }
0099 
0100 bool PythiaMomDauFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0101   bool accepted = false;
0102   bool mom_accepted = false;
0103   edm::Handle<edm::HepMCProduct> evt;
0104   iEvent.getByToken(label_, evt);
0105 
0106   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0107   int ndauac = 0;
0108   int ndau = 0;
0109   int ndesac = 0;
0110   int ndes = 0;
0111 
0112   for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0113        ++p) {
0114     if ((*p)->pdg_id() != particleID)
0115       continue;
0116     HepMC::FourVector mom = MCFilterZboostHelper::zboost((*p)->momentum(), betaBoost);
0117     if ((*p)->momentum().perp() > mom_minptcut && (*p)->momentum().perp() < mom_maxptcut && mom.eta() > mom_minetacut &&
0118         mom.eta() < mom_maxetacut) {
0119       mom_accepted = true;
0120       ndauac = 0;
0121       ndau = 0;
0122       ndesac = 0;
0123       ndes = 0;
0124       if ((*p)->end_vertex()) {
0125         for (HepMC::GenVertex::particle_iterator dau = (*p)->end_vertex()->particles_begin(HepMC::children);
0126              dau != (*p)->end_vertex()->particles_end(HepMC::children);
0127              ++dau) {
0128           ++ndau;
0129           for (unsigned int i = 0; i < dauIDs.size(); ++i) {
0130             if ((*dau)->pdg_id() != dauIDs[i])
0131               continue;
0132             ++ndauac;
0133             break;
0134           }
0135           if ((*dau)->pdg_id() == daughterID) {
0136             for (HepMC::GenVertex::particle_iterator des = (*dau)->end_vertex()->particles_begin(HepMC::children);
0137                  des != (*des)->end_vertex()->particles_end(HepMC::children);
0138                  ++des) {
0139               ++ndes;
0140               for (unsigned int i = 0; i < desIDs.size(); ++i) {
0141                 if ((*des)->pdg_id() != desIDs[i])
0142                   continue;
0143                 HepMC::FourVector dau_i = MCFilterZboostHelper::zboost((*des)->momentum(), betaBoost);
0144                 if ((*des)->momentum().perp() > minptcut && (*des)->momentum().perp() < maxptcut &&
0145                     dau_i.eta() > minetacut && dau_i.eta() < maxetacut) {
0146                   ++ndesac;
0147                   break;
0148                 }
0149               }
0150             }
0151           }
0152         }
0153       }
0154     }
0155     if (ndau == ndaughters && ndauac == ndaughters && mom_accepted && ndes == ndescendants && ndesac == ndescendants) {
0156       accepted = true;
0157       break;
0158     }
0159   }
0160 
0161   if (!accepted && chargeconju) {
0162     mom_accepted = false;
0163     for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0164          ++p) {
0165       if ((*p)->pdg_id() != -particleID)
0166         continue;
0167       HepMC::FourVector mom = MCFilterZboostHelper::zboost((*p)->momentum(), betaBoost);
0168       if ((*p)->momentum().perp() > mom_minptcut && (*p)->momentum().perp() < mom_maxptcut &&
0169           mom.eta() > mom_minetacut && mom.eta() < mom_maxetacut) {
0170         mom_accepted = true;
0171         ndauac = 0;
0172         ndau = 0;
0173         ndesac = 0;
0174         ndes = 0;
0175         if ((*p)->end_vertex()) {
0176           for (HepMC::GenVertex::particle_iterator dau = (*p)->end_vertex()->particles_begin(HepMC::children);
0177                dau != (*p)->end_vertex()->particles_end(HepMC::children);
0178                ++dau) {
0179             ++ndau;
0180             for (unsigned int i = 0; i < dauIDs.size(); ++i) {
0181               int IDanti = -dauIDs[i];
0182               int pythiaCode = PYCOMP(dauIDs[i]);
0183               int has_antipart = pydat2.kchg[3 - 1][pythiaCode - 1];
0184               if (has_antipart == 0)
0185                 IDanti = dauIDs[i];
0186               if ((*dau)->pdg_id() != IDanti)
0187                 continue;
0188               ++ndauac;
0189               break;
0190             }
0191             int daughterIDanti = -daughterID;
0192             int pythiaCode = PYCOMP(daughterID);
0193             int has_antipart = pydat2.kchg[3 - 1][pythiaCode - 1];
0194             if (has_antipart == 0)
0195               daughterIDanti = daughterID;
0196             if ((*dau)->pdg_id() == daughterIDanti) {
0197               for (HepMC::GenVertex::particle_iterator des = (*dau)->end_vertex()->particles_begin(HepMC::children);
0198                    des != (*des)->end_vertex()->particles_end(HepMC::children);
0199                    ++des) {
0200                 ++ndes;
0201                 for (unsigned int i = 0; i < desIDs.size(); ++i) {
0202                   int IDanti = -desIDs[i];
0203                   int pythiaCode = PYCOMP(desIDs[i]);
0204                   int has_antipart = pydat2.kchg[3 - 1][pythiaCode - 1];
0205                   if (has_antipart == 0)
0206                     IDanti = desIDs[i];
0207                   if ((*des)->pdg_id() != IDanti)
0208                     continue;
0209                   HepMC::FourVector dau_i = MCFilterZboostHelper::zboost((*des)->momentum(), betaBoost);
0210                   if ((*des)->momentum().perp() > minptcut && (*des)->momentum().perp() < maxptcut &&
0211                       dau_i.eta() > minetacut && dau_i.eta() < maxetacut) {
0212                     ++ndesac;
0213                     break;
0214                   }
0215                 }
0216               }
0217             }
0218           }
0219         }
0220       }
0221       if (ndau == ndaughters && ndauac == ndaughters && mom_accepted && ndes == ndescendants &&
0222           ndesac == ndescendants) {
0223         accepted = true;
0224         break;
0225       }
0226     }
0227   }
0228   return accepted;
0229 }
0230 
0231 DEFINE_FWK_MODULE(PythiaMomDauFilter);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!