Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:04:34 

0001 #include "GeneratorInterface/GenFilters/plugins/PythiaDauVFilterMatchID.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 PythiaDauVFilterMatchID::PythiaDauVFilterMatchID(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("PythiaDauVFilterMatchID")
0036       << "----------------------------------------------------------------------" << endl;
0037   edm::LogInfo("PythiaDauVFilterMatchID") << "--- PythiaDauVFilterMatchID" << endl;
0038   for (unsigned int i = 0; i < dauIDs.size(); ++i) {
0039     edm::LogInfo("PythiaDauVFilterMatchID")
0040         << "ID: " << dauIDs[i] << " pT > " << minptcut[i] << " " << minetacut[i] << " eta < " << maxetacut[i] << endl;
0041   }
0042   edm::LogInfo("PythiaDauVFilterMatchID") << "maxptcut   = " << maxptcut << endl;
0043   edm::LogInfo("PythiaDauVFilterMatchID") << "particleID = " << particleID << endl;
0044   edm::LogInfo("PythiaDauVFilterMatchID") << "motherID   = " << motherID << endl;
0045 
0046   // create pythia8 instance to access particle data
0047   edm::LogInfo("PythiaDauVFilterMatchID") << "Creating pythia8 instance for particle properties" << endl;
0048   if (!fLookupGen.get())
0049     fLookupGen = std::make_unique<Pythia>();
0050 }
0051 
0052 PythiaDauVFilterMatchID::~PythiaDauVFilterMatchID() {
0053   // do anything here that needs to be done at desctruction time
0054   // (e.g. close files, deallocate resources etc.)
0055 }
0056 
0057 //
0058 // member functions
0059 //
0060 
0061 // ------------ method called to produce the data  ------------
0062 bool PythiaDauVFilterMatchID::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0063   using namespace edm;
0064   bool accepted = false;
0065   Handle<HepMCProduct> evt;
0066   iEvent.getByToken(token_, evt);
0067 
0068   int OK(1);
0069   vector<int> vparticles;
0070 
0071   HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
0072 
0073   if (fVerbose > 5) {
0074     edm::LogInfo("PythiaDauVFilterMatchID") << "looking for " << particleID << endl;
0075   }
0076 
0077   for (HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0078     if ((*p)->pdg_id() != particleID)
0079       continue;
0080 
0081     // -- Check for mother of this particle
0082     if (0 != motherID) {
0083       OK = 0;
0084       for (HepMC::GenVertex::particles_in_const_iterator des = (*p)->production_vertex()->particles_in_const_begin();
0085            des != (*p)->production_vertex()->particles_in_const_end();
0086            ++des) {
0087         if (fVerbose > 10) {
0088           edm::LogInfo("PythiaDauVFilterMatchID")
0089               << "mother: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0090               << " eta: " << (*des)->momentum().eta() << endl;
0091         }
0092         if (abs(motherID) == abs((*des)->pdg_id())) {
0093           OK = 1;
0094           break;
0095         }
0096       }
0097     }
0098     if (0 == OK)
0099       continue;
0100 
0101     //generate targets
0102     std::vector<decayTarget> targets;
0103     for (unsigned int i = 0; i < dauIDs.size(); i++) {
0104       decayTarget target;
0105       target.pdgID = dauIDs[i];
0106       target.minPt = minptcut[i];
0107       target.maxPt = maxptcut;
0108       target.minEta = minetacut[i];
0109       target.maxEta = maxetacut[i];
0110       targets.push_back(target);
0111     }
0112     if (fVerbose > 10) {
0113       edm::LogInfo("PythiaDauVFilterMatchID") << "created target: ";
0114       for (unsigned int i = 0; i < targets.size(); i++) {
0115         edm::LogInfo("PythiaDauVFilterMatchID") << targets[i].pdgID << " ";
0116       }
0117       edm::LogInfo("PythiaDauVFilterMatchID") << endl;
0118     }
0119 
0120     // -- check for daugthers
0121     int ndau = 0;
0122     if (fVerbose > 5) {
0123       edm::LogInfo("PythiaDauVFilterMatchID") << "found ID: " << (*p)->pdg_id() << " pT: " << (*p)->momentum().perp()
0124                                               << " eta: " << (*p)->momentum().eta() << endl;
0125     }
0126     vparticles.push_back((*p)->pdg_id());
0127     if ((*p)->end_vertex()) {
0128       for (HepMC::GenVertex::particle_iterator des = (*p)->end_vertex()->particles_begin(HepMC::children);
0129            des != (*p)->end_vertex()->particles_end(HepMC::children);
0130            ++des) {
0131         if (TMath::Abs((*des)->pdg_id()) == 22) {
0132           continue;
0133         }
0134         ++ndau;
0135         if (fVerbose > 5) {
0136           edm::LogInfo("PythiaDauVFilterMatchID") << "ID: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0137                                                   << " eta: " << (*des)->momentum().eta() << endl;
0138         }
0139         {  // protect matchedIdx
0140           int matchedIdx = -1;
0141           for (unsigned int i = 0; i < targets.size(); i++) {
0142             if ((*des)->pdg_id() != targets[i].pdgID) {
0143               continue;
0144             }
0145             if (fVerbose > 5) {
0146               edm::LogInfo("PythiaDauVFilterMatchID") << "i = " << i << " pT = " << (*des)->momentum().perp()
0147                                                       << " eta = " << (*des)->momentum().eta() << endl;
0148             }
0149 
0150             if ((*des)->momentum().perp() > targets[i].minPt && (*des)->momentum().perp() < targets[i].maxPt &&
0151                 (*des)->momentum().eta() > targets[i].minEta && (*des)->momentum().eta() < targets[i].maxEta) {
0152               vparticles.push_back((*des)->pdg_id());
0153               if (fVerbose > 2) {
0154                 edm::LogInfo("PythiaDauVFilterMatchID")
0155                     << "  accepted this particle " << (*des)->pdg_id() << " pT = " << (*des)->momentum().perp()
0156                     << " eta = " << (*des)->momentum().eta() << endl;
0157               }
0158               matchedIdx = i;
0159               break;
0160             }
0161           }
0162           if (matchedIdx > -1) {
0163             targets.erase(targets.begin() + matchedIdx);
0164           }
0165           if (fVerbose > 10) {
0166             edm::LogInfo("PythiaDauVFilterMatchID") << "remaining targets: ";
0167             for (unsigned int i = 0; i < targets.size(); i++) {
0168               edm::LogInfo("PythiaDauVFilterMatchID") << targets[i].pdgID << " ";
0169             }
0170             edm::LogInfo("PythiaDauVFilterMatchID") << endl;
0171           }
0172         }
0173       }
0174     }
0175 
0176     if (ndau == ndaughters && targets.empty()) {
0177       accepted = true;
0178       if (fVerbose > 0) {
0179         edm::LogInfo("PythiaDauVFilterMatchID") << "  accepted this decay: ";
0180         for (unsigned int iv = 0; iv < vparticles.size(); ++iv)
0181           edm::LogInfo("PythiaDauVFilterMatchID") << vparticles[iv] << " ";
0182         edm::LogInfo("PythiaDauVFilterMatchID") << " from mother = " << motherID << endl;
0183       }
0184       break;
0185     }
0186   }
0187 
0188   int anti_particleID = -particleID;
0189   if (!(fLookupGen->particleData.isParticle(anti_particleID))) {
0190     anti_particleID = 0;
0191     if (fVerbose > 5)
0192       edm::LogInfo("PythiaDauVFilterMatchID")
0193           << "Particle " << particleID << " is its own anti-particle, skipping further testing " << endl;
0194   }
0195   if (!accepted && chargeconju && anti_particleID) {
0196     OK = 1;
0197 
0198     for (HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0199       if ((*p)->pdg_id() != anti_particleID)
0200         continue;
0201 
0202       // -- Check for mother of this particle
0203       if (0 != motherID) {
0204         OK = 0;
0205         for (HepMC::GenVertex::particles_in_const_iterator des = (*p)->production_vertex()->particles_in_const_begin();
0206              des != (*p)->production_vertex()->particles_in_const_end();
0207              ++des) {
0208           if (fVerbose > 10) {
0209             edm::LogInfo("PythiaDauVFilterMatchID")
0210                 << "mother: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0211                 << " eta: " << (*des)->momentum().eta() << endl;
0212           }
0213           if (abs(motherID) == abs((*des)->pdg_id())) {
0214             OK = 1;
0215             break;
0216           }
0217         }
0218       }
0219       if (0 == OK)
0220         continue;
0221 
0222       //generate anti targets
0223       std::vector<decayTarget> targets;
0224       for (unsigned int i = 0; i < dauIDs.size(); i++) {
0225         decayTarget target;
0226         int IDanti = -dauIDs[i];
0227         if (!(fLookupGen->particleData.isParticle(IDanti)))
0228           IDanti = dauIDs[i];
0229         target.pdgID = IDanti;
0230         target.minPt = minptcut[i];
0231         target.maxPt = maxptcut;
0232         target.minEta = minetacut[i];
0233         target.maxEta = maxetacut[i];
0234         targets.push_back(target);
0235       }
0236       if (fVerbose > 10) {
0237         edm::LogInfo("PythiaDauVFilterMatchID") << "created anti target: ";
0238         for (unsigned int i = 0; i < targets.size(); i++) {
0239           edm::LogInfo("PythiaDauVFilterMatchID") << targets[i].pdgID << " ";
0240         }
0241         edm::LogInfo("PythiaDauVFilterMatchID") << endl;
0242       }
0243 
0244       if (fVerbose > 5) {
0245         edm::LogInfo("PythiaDauVFilterMatchID") << "found ID: " << (*p)->pdg_id() << " pT: " << (*p)->momentum().perp()
0246                                                 << " eta: " << (*p)->momentum().eta() << endl;
0247       }
0248       vparticles.push_back((*p)->pdg_id());
0249       int ndau = 0;
0250       if ((*p)->end_vertex()) {
0251         for (HepMC::GenVertex::particle_iterator des = (*p)->end_vertex()->particles_begin(HepMC::children);
0252              des != (*p)->end_vertex()->particles_end(HepMC::children);
0253              ++des) {
0254           if (TMath::Abs((*des)->pdg_id()) == 22) {
0255             continue;
0256           }
0257           ++ndau;
0258           if (fVerbose > 5) {
0259             edm::LogInfo("PythiaDauVFilterMatchID")
0260                 << "ID: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0261                 << " eta: " << (*des)->momentum().eta() << endl;
0262           }
0263           {  // protect matchedIdx
0264             int matchedIdx = -1;
0265             for (unsigned int i = 0; i < targets.size(); i++) {
0266               if ((*des)->pdg_id() != targets[i].pdgID) {
0267                 continue;
0268               }
0269               if (fVerbose > 5) {
0270                 edm::LogInfo("PythiaDauVFilterMatchID") << "i = " << i << " pT = " << (*des)->momentum().perp()
0271                                                         << " eta = " << (*des)->momentum().eta() << endl;
0272               }
0273 
0274               if ((*des)->momentum().perp() > targets[i].minPt && (*des)->momentum().perp() < targets[i].maxPt &&
0275                   (*des)->momentum().eta() > targets[i].minEta && (*des)->momentum().eta() < targets[i].maxEta) {
0276                 vparticles.push_back((*des)->pdg_id());
0277                 if (fVerbose > 2) {
0278                   edm::LogInfo("PythiaDauVFilterMatchID")
0279                       << "  accepted this particle " << (*des)->pdg_id() << " pT = " << (*des)->momentum().perp()
0280                       << " eta = " << (*des)->momentum().eta() << endl;
0281                 }
0282                 matchedIdx = i;
0283                 break;
0284               }
0285             }
0286             if (matchedIdx > -1) {
0287               targets.erase(targets.begin() + matchedIdx);
0288             }
0289             if (fVerbose > 10) {
0290               edm::LogInfo("PythiaDauVFilterMatchID") << "remaining targets: ";
0291               for (unsigned int i = 0; i < targets.size(); i++) {
0292                 edm::LogInfo("PythiaDauVFilterMatchID") << targets[i].pdgID << " ";
0293               }
0294               edm::LogInfo("PythiaDauVFilterMatchID") << endl;
0295             }
0296           }
0297         }
0298       }
0299       if (ndau == ndaughters && targets.empty()) {
0300         accepted = true;
0301         if (fVerbose > 0) {
0302           edm::LogInfo("PythiaDauVFilterMatchID") << "  accepted this decay: ";
0303           for (unsigned int iv = 0; iv < vparticles.size(); ++iv)
0304             edm::LogInfo("PythiaDauVFilterMatchID") << vparticles[iv] << " ";
0305           edm::LogInfo("PythiaDauVFilterMatchID") << " from mother = " << motherID << endl;
0306         }
0307         break;
0308       }
0309     }
0310   }
0311 
0312   delete myGenEvent;
0313   return accepted;
0314 }
[ 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!