Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 
0002 #include "GeneratorInterface/GenFilters/plugins/PythiaProbeFilter.h"
0003 
0004 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0005 #include <iostream>
0006 #include <memory>
0007 
0008 using namespace edm;
0009 using namespace std;
0010 using namespace Pythia8;
0011 
0012 PythiaProbeFilter::PythiaProbeFilter(const edm::ParameterSet& iConfig)
0013     : token_(consumes<edm::HepMCProduct>(
0014           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0015       particleID(iConfig.getUntrackedParameter("ParticleID", 0)),
0016       MomID(iConfig.getUntrackedParameter("MomID", 0)),
0017       GrandMomID(iConfig.getUntrackedParameter("GrandMomID", 0)),
0018       chargeconju(iConfig.getUntrackedParameter("ChargeConjugation", true)),
0019       nsisters(iConfig.getUntrackedParameter("NumberOfSisters", 0)),
0020       naunts(iConfig.getUntrackedParameter("NumberOfAunts", 0)),
0021       minptcut(iConfig.getUntrackedParameter("MinPt", 0.)),
0022       maxptcut(iConfig.getUntrackedParameter("MaxPt", 14000.)),
0023       minetacut(iConfig.getUntrackedParameter("MinEta", -10.)),
0024       maxetacut(iConfig.getUntrackedParameter("MaxEta", 10.)),
0025       countQEDCorPhotons(iConfig.getUntrackedParameter("countQEDCorPhotons", false)) {
0026   //now do what ever initialization is needed
0027   vector<int> defID;
0028   defID.push_back(0);
0029   exclsisIDs = iConfig.getUntrackedParameter<vector<int> >("SisterIDs", defID);
0030   exclauntIDs = iConfig.getUntrackedParameter<vector<int> >("AuntIDs", defID);
0031   identicalParticle = false;
0032   for (unsigned int ilist = 0; ilist < exclsisIDs.size(); ++ilist) {
0033     if (fabs(exclsisIDs[ilist]) == fabs(particleID))
0034       identicalParticle = true;
0035   }
0036   // create pythia8 instance to access particle data
0037   edm::LogInfo("PythiaProbeFilter::PythiaProbeFilter") << "Creating pythia8 instance for particle properties" << endl;
0038   if (!fLookupGen.get())
0039     fLookupGen = std::make_unique<Pythia>();
0040 }
0041 
0042 PythiaProbeFilter::~PythiaProbeFilter() {
0043   // do anything here that needs to be done at desctruction time
0044   // (e.g. close files, deallocate resources etc.)
0045 }
0046 
0047 bool PythiaProbeFilter::AlreadyExcludedCheck(std::vector<unsigned int> excludedList, unsigned int current_part) const {
0048   bool result = false;
0049   for (unsigned int checkNow : excludedList) {
0050     if (current_part != checkNow)
0051       continue;
0052     result = true;
0053     break;
0054   }
0055   return result;
0056 }
0057 //
0058 // member functions
0059 //
0060 // ------------ method called to produce the data  ------------
0061 bool PythiaProbeFilter::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   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0068 
0069   //access particles
0070   for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0071        ++p) {
0072     //select tag particle
0073     if (fabs((*p)->pdg_id()) == fabs(particleID))
0074       if ((*p)->pdg_id() != particleID && !chargeconju)
0075         continue;
0076 
0077     if (fabs((*p)->pdg_id()) != fabs(particleID) && chargeconju)
0078       continue;
0079     //kinematic properties of tag
0080     if ((*p)->momentum().perp() < minptcut || (*p)->momentum().perp() > maxptcut)
0081       continue;
0082     if ((*p)->momentum().eta() < minetacut || (*p)->momentum().eta() > maxetacut)
0083       continue;
0084     //remove probe side particles
0085     bool excludeTagParticle = false;
0086     if (naunts == 0) {
0087       if ((*p)->production_vertex()) {
0088         for (HepMC::GenVertex::particle_iterator anc = (*p)->production_vertex()->particles_begin(HepMC::parents);
0089              anc != (*p)->production_vertex()->particles_end(HepMC::parents);
0090              ++anc) {
0091           if (fabs((*anc)->pdg_id()) != fabs(MomID) && chargeconju)
0092             continue;
0093           else if ((*anc)->pdg_id() != MomID && !chargeconju)
0094             continue;
0095           int nsis = 0;
0096           int exclsis = 0;
0097           std::vector<unsigned int> checklistSis;
0098           if ((*anc)->end_vertex()) {
0099             for (HepMC::GenVertex::particle_iterator sis = (*anc)->end_vertex()->particles_begin(HepMC::children);
0100                  sis != (*anc)->end_vertex()->particles_end(HepMC::children);
0101                  ++sis) {
0102               //identify the tag particle in the decay
0103               if ((*p)->pdg_id() == (*sis)->pdg_id() && (identicalParticle || !chargeconju))
0104                 continue;
0105               if (fabs((*p)->pdg_id()) == fabs((*sis)->pdg_id()) && !identicalParticle && chargeconju)
0106                 continue;
0107               //remove QED photons
0108               if ((*sis)->pdg_id() == 22 && !countQEDCorPhotons)
0109                 continue;
0110               nsis++;
0111               //check if this sis is excluded already
0112               for (unsigned int ilist = 0; ilist < exclsisIDs.size(); ++ilist) {
0113                 if (AlreadyExcludedCheck(checklistSis, ilist)) {
0114                   continue;
0115                 }
0116                 if (fabs(exclsisIDs[ilist]) == fabs((*sis)->pdg_id()) && chargeconju) {
0117                   exclsis++;
0118                   checklistSis.push_back(ilist);
0119                 }
0120                 if (exclsisIDs[ilist] == (*sis)->pdg_id() && !chargeconju) {
0121                   exclsis++;
0122                   checklistSis.push_back(ilist);
0123                 }
0124               }
0125             }
0126           }
0127           if (nsis == exclsis && nsis == nsisters) {
0128             excludeTagParticle = true;
0129             break;
0130           }
0131         }
0132       }
0133     } else if (naunts > 0) {
0134       //now take into account that we have up 2 generations in the decay
0135       if ((*p)->production_vertex()) {
0136         for (HepMC::GenVertex::particle_iterator anc = (*p)->production_vertex()->particles_begin(HepMC::parents);
0137              anc != (*p)->production_vertex()->particles_end(HepMC::parents);
0138              ++anc) {
0139           if (fabs((*anc)->pdg_id()) != fabs(MomID) && chargeconju)
0140             continue;
0141           else if ((*anc)->pdg_id() != MomID && !chargeconju)
0142             continue;
0143           int nsis = 0;
0144           int exclsis = 0;
0145           std::vector<unsigned int> checklistSis;
0146           int naunt = 0;
0147           int exclaunt = 0;
0148           std::vector<unsigned int> checklistAunt;
0149           if ((*anc)->end_vertex()) {
0150             for (HepMC::GenVertex::particle_iterator sis = (*anc)->end_vertex()->particles_begin(HepMC::children);
0151                  sis != (*anc)->end_vertex()->particles_end(HepMC::children);
0152                  ++sis) {
0153               //identify the particle under study in the decay
0154               if ((*p)->pdg_id() == (*sis)->pdg_id() && (identicalParticle || !chargeconju))
0155                 continue;
0156               if (fabs((*p)->pdg_id()) == fabs((*sis)->pdg_id()) && !identicalParticle && chargeconju)
0157                 continue;
0158               //remove QED photons
0159               if ((*sis)->pdg_id() == 22 && !countQEDCorPhotons)
0160                 continue;
0161               nsis++;
0162               for (unsigned int ilist = 0; ilist < exclsisIDs.size(); ++ilist) {
0163                 if (AlreadyExcludedCheck(checklistSis, ilist))
0164                   continue;
0165                 if (fabs(exclsisIDs[ilist]) == fabs((*sis)->pdg_id()) && chargeconju) {
0166                   exclsis++;
0167                   checklistSis.push_back(ilist);
0168                 }
0169                 if (exclsisIDs[ilist] == (*sis)->pdg_id() && !chargeconju) {
0170                   exclsis++;
0171                   checklistSis.push_back(ilist);
0172                 }
0173               }
0174             }
0175           }
0176           //check sisters
0177           if (nsis != exclsis || nsis != nsisters)
0178             break;
0179           if ((*anc)->production_vertex()) {
0180             for (HepMC::GenVertex::particle_iterator granc =
0181                      (*anc)->production_vertex()->particles_begin(HepMC::parents);
0182                  granc != (*anc)->production_vertex()->particles_end(HepMC::parents);
0183                  ++granc) {
0184               if (fabs((*granc)->pdg_id()) != fabs(GrandMomID) && chargeconju)
0185                 continue;
0186               else if ((*granc)->pdg_id() != GrandMomID && !chargeconju)
0187                 continue;
0188               for (HepMC::GenVertex::particle_iterator aunt = (*granc)->end_vertex()->particles_begin(HepMC::children);
0189                    aunt != (*granc)->end_vertex()->particles_end(HepMC::children);
0190                    ++aunt) {
0191                 if ((*aunt)->pdg_id() == (*anc)->pdg_id())
0192                   continue;
0193                 if ((*aunt)->pdg_id() == 22 && !countQEDCorPhotons)
0194                   continue;
0195                 naunt++;
0196                 for (unsigned int ilist = 0; ilist < exclauntIDs.size(); ++ilist) {
0197                   if (AlreadyExcludedCheck(checklistAunt, ilist))
0198                     continue;
0199                   if (fabs(exclauntIDs[ilist]) == fabs((*aunt)->pdg_id()) && chargeconju) {
0200                     exclaunt++;
0201                     checklistAunt.push_back(ilist);
0202                   }
0203                   if (exclauntIDs[ilist] == (*aunt)->pdg_id() && !chargeconju) {
0204                     exclaunt++;
0205                     checklistAunt.push_back(ilist);
0206                   }
0207                 }
0208               }
0209             }
0210           }
0211           //check aunts
0212           if (naunt == exclaunt && naunt == naunts) {
0213             excludeTagParticle = true;
0214             break;
0215           }
0216         }
0217       }
0218     }
0219     if (excludeTagParticle)
0220       continue;
0221     accepted = true;
0222     break;
0223   }
0224 
0225   if (accepted) {
0226     return true;
0227   } else {
0228     return false;
0229   }
0230 }
[ 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!