Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 #include "GeneratorInterface/GenFilters/plugins/PythiaDauVFilter.h"
0002 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0003 #include <iostream>
0004 #include <memory>
0005 
0006 #include <vector>
0007 
0008 using namespace edm;
0009 using namespace std;
0010 using namespace Pythia8;
0011 
0012 PythiaDauVFilter::PythiaDauVFilter(const edm::ParameterSet& iConfig)
0013     : fVerbose(iConfig.getUntrackedParameter("verbose", 0)),
0014       token_(consumes<edm::HepMCProduct>(
0015           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0016       particleID(iConfig.getUntrackedParameter("ParticleID", 0)),
0017       motherID(iConfig.getUntrackedParameter("MotherID", 0)),
0018       chargeconju(iConfig.getUntrackedParameter("ChargeConjugation", true)),
0019       ndaughters(iConfig.getUntrackedParameter("NumberDaughters", 0)),
0020       maxptcut(iConfig.getUntrackedParameter("MaxPt", 14000.)) {
0021   //now do what ever initialization is needed
0022   vector<int> defdauID;
0023   defdauID.push_back(0);
0024   dauIDs = iConfig.getUntrackedParameter<vector<int> >("DaughterIDs", defdauID);
0025   vector<double> defminptcut;
0026   defminptcut.push_back(0.);
0027   minptcut = iConfig.getUntrackedParameter<vector<double> >("MinPt", defminptcut);
0028   vector<double> defminetacut;
0029   defminetacut.push_back(-10.);
0030   minetacut = iConfig.getUntrackedParameter<vector<double> >("MinEta", defminetacut);
0031   vector<double> defmaxetacut;
0032   defmaxetacut.push_back(10.);
0033   maxetacut = iConfig.getUntrackedParameter<vector<double> >("MaxEta", defmaxetacut);
0034 
0035   edm::LogInfo("PythiaDauVFilter") << "----------------------------------------------------------------------" << endl;
0036   edm::LogInfo("PythiaDauVFilter") << "--- PythiaDauVFilter" << endl;
0037   for (unsigned int i = 0; i < dauIDs.size(); ++i) {
0038     edm::LogInfo("PythiaDauVFilter") << "ID: " << dauIDs[i] << " pT > " << minptcut[i] << " " << minetacut[i]
0039                                      << " eta < " << maxetacut[i] << endl;
0040   }
0041   edm::LogInfo("PythiaDauVFilter") << "maxptcut   = " << maxptcut << endl;
0042   edm::LogInfo("PythiaDauVFilter") << "particleID = " << particleID << endl;
0043   edm::LogInfo("PythiaDauVFilter") << "motherID   = " << motherID << endl;
0044 
0045   // create pythia8 instance to access particle data
0046   edm::LogInfo("PythiaDauVFilter") << "Creating pythia8 instance for particle properties" << endl;
0047   if (!fLookupGen.get())
0048     fLookupGen = std::make_unique<Pythia>();
0049 }
0050 
0051 PythiaDauVFilter::~PythiaDauVFilter() {
0052   // do anything here that needs to be done at desctruction time
0053   // (e.g. close files, deallocate resources etc.)
0054 }
0055 
0056 //
0057 // member functions
0058 //
0059 
0060 // ------------ method called to produce the data  ------------
0061 bool PythiaDauVFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0062   using namespace edm;
0063   bool accepted = false;
0064   Handle<HepMCProduct> evt;
0065   iEvent.getByToken(token_, evt);
0066 
0067   int OK(1);
0068   vector<int> vparticles;
0069 
0070   HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
0071 
0072   if (fVerbose > 5) {
0073     edm::LogInfo("PythiaDauVFilter") << "looking for " << particleID << endl;
0074   }
0075 
0076   for (HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0077     if ((*p)->pdg_id() != particleID)
0078       continue;
0079 
0080     // -- Check for mother of this particle
0081     if (0 != motherID) {
0082       OK = 0;
0083       for (HepMC::GenVertex::particles_in_const_iterator des = (*p)->production_vertex()->particles_in_const_begin();
0084            des != (*p)->production_vertex()->particles_in_const_end();
0085            ++des) {
0086         if (fVerbose > 10) {
0087           edm::LogInfo("PythiaDauVFilter") << "mother: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0088                                            << " eta: " << (*des)->momentum().eta() << endl;
0089         }
0090         if (abs(motherID) == abs((*des)->pdg_id())) {
0091           OK = 1;
0092           break;
0093         }
0094       }
0095     }
0096     if (0 == OK)
0097       continue;
0098 
0099     // -- check for daugthers
0100     int ndauac = 0;
0101     int ndau = 0;
0102     if (fVerbose > 5) {
0103       edm::LogInfo("PythiaDauVFilter") << "found ID: " << (*p)->pdg_id() << " pT: " << (*p)->momentum().perp()
0104                                        << " eta: " << (*p)->momentum().eta() << endl;
0105     }
0106     vparticles.push_back((*p)->pdg_id());
0107     if ((*p)->end_vertex()) {
0108       for (HepMC::GenVertex::particle_iterator des = (*p)->end_vertex()->particles_begin(HepMC::children);
0109            des != (*p)->end_vertex()->particles_end(HepMC::children);
0110            ++des) {
0111         ++ndau;
0112         if (fVerbose > 5) {
0113           edm::LogInfo("PythiaDauVFilter") << "ID: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0114                                            << " eta: " << (*des)->momentum().eta() << endl;
0115         }
0116         for (unsigned int i = 0; i < dauIDs.size(); ++i) {
0117           if ((*des)->pdg_id() != dauIDs[i])
0118             continue;
0119           if (fVerbose > 5) {
0120             edm::LogInfo("PythiaDauVFilter") << "i = " << i << " pT = " << (*des)->momentum().perp()
0121                                              << " eta = " << (*des)->momentum().eta() << endl;
0122           }
0123           if ((*des)->momentum().perp() > minptcut[i] && (*des)->momentum().perp() < maxptcut &&
0124               (*des)->momentum().eta() > minetacut[i] && (*des)->momentum().eta() < maxetacut[i]) {
0125             ++ndauac;
0126             vparticles.push_back((*des)->pdg_id());
0127             if (fVerbose > 2) {
0128               edm::LogInfo("PythiaDauVFilter")
0129                   << "  accepted this particle " << (*des)->pdg_id() << " pT = " << (*des)->momentum().perp()
0130                   << " eta = " << (*des)->momentum().eta() << endl;
0131             }
0132             break;
0133           }
0134         }
0135       }
0136     }
0137 
0138     // -- allow photons
0139     if (ndau >= ndaughters && ndauac == ndaughters) {
0140       accepted = true;
0141       if (fVerbose > 0) {
0142         edm::LogInfo("PythiaDauVFilter") << "  accepted this decay: ";
0143         for (unsigned int iv = 0; iv < vparticles.size(); ++iv)
0144           edm::LogInfo("PythiaDauVFilter") << vparticles[iv] << " ";
0145         edm::LogInfo("PythiaDauVFilter") << " from mother = " << motherID << endl;
0146       }
0147       break;
0148     }
0149   }
0150 
0151   if (!accepted && chargeconju) {
0152     OK = 1;
0153 
0154     for (HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0155       if ((*p)->pdg_id() != -particleID)
0156         continue;
0157 
0158       // -- Check for mother of this particle
0159       if (0 != motherID) {
0160         OK = 0;
0161         for (HepMC::GenVertex::particles_in_const_iterator des = (*p)->production_vertex()->particles_in_const_begin();
0162              des != (*p)->production_vertex()->particles_in_const_end();
0163              ++des) {
0164           if (fVerbose > 10) {
0165             edm::LogInfo("PythiaDauVFilter") << "mother: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0166                                              << " eta: " << (*des)->momentum().eta() << endl;
0167           }
0168           if (abs(motherID) == abs((*des)->pdg_id())) {
0169             OK = 1;
0170             break;
0171           }
0172         }
0173       }
0174       if (0 == OK)
0175         continue;
0176 
0177       if (fVerbose > 5) {
0178         edm::LogInfo("PythiaDauVFilter") << "found ID: " << (*p)->pdg_id() << " pT: " << (*p)->momentum().perp()
0179                                          << " eta: " << (*p)->momentum().eta() << endl;
0180       }
0181       vparticles.push_back((*p)->pdg_id());
0182       int ndauac = 0;
0183       int ndau = 0;
0184       if ((*p)->end_vertex()) {
0185         for (HepMC::GenVertex::particle_iterator des = (*p)->end_vertex()->particles_begin(HepMC::children);
0186              des != (*p)->end_vertex()->particles_end(HepMC::children);
0187              ++des) {
0188           ++ndau;
0189           if (fVerbose > 5) {
0190             edm::LogInfo("PythiaDauVFilter") << "ID: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0191                                              << " eta: " << (*des)->momentum().eta() << endl;
0192           }
0193           for (unsigned int i = 0; i < dauIDs.size(); ++i) {
0194             int IDanti = -dauIDs[i];
0195             if (!(fLookupGen->particleData.isParticle(IDanti)))
0196               IDanti = dauIDs[i];
0197             if ((*des)->pdg_id() != IDanti)
0198               continue;
0199             if (fVerbose > 5) {
0200               edm::LogInfo("PythiaDauVFilter") << "i = " << i << " pT = " << (*des)->momentum().perp()
0201                                                << " eta = " << (*des)->momentum().eta() << endl;
0202             }
0203             if ((*des)->momentum().perp() > minptcut[i] && (*des)->momentum().perp() < maxptcut &&
0204                 (*des)->momentum().eta() > minetacut[i] && (*des)->momentum().eta() < maxetacut[i]) {
0205               ++ndauac;
0206               vparticles.push_back((*des)->pdg_id());
0207               if (fVerbose > 2) {
0208                 edm::LogInfo("PythiaDauVFilter")
0209                     << "  accepted this particle " << (*des)->pdg_id() << " pT = " << (*des)->momentum().perp()
0210                     << " eta = " << (*des)->momentum().eta() << endl;
0211               }
0212               break;
0213             }
0214           }
0215         }
0216       }
0217       if (ndau >= ndaughters && ndauac == ndaughters) {
0218         accepted = true;
0219         if (fVerbose > 0) {
0220           edm::LogInfo("PythiaDauVFilter") << "  accepted this anti-decay: ";
0221           for (unsigned int iv = 0; iv < vparticles.size(); ++iv)
0222             edm::LogInfo("PythiaDauVFilter") << vparticles[iv] << " ";
0223           edm::LogInfo("PythiaDauVFilter") << " from mother = " << motherID << endl;
0224         }
0225         break;
0226       }
0227     }
0228   }
0229 
0230   delete myGenEvent;
0231   return accepted;
0232 }
0233 
0234 //define this as a plug-in
0235 DEFINE_FWK_MODULE(PythiaDauVFilter);
[ 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!