GenFilters/plugins/PhotonGenFilter.cc

0001 #include "GeneratorInterface/GenFilters/plugins/PhotonGenFilter.h"
0002 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0003 #include <iostream>
0004
0005 using namespace edm;
0006 using namespace std;
0007
0008 PhotonGenFilter::PhotonGenFilter(const edm::ParameterSet &iConfig)
0009     : token_(consumes<edm::HepMCProduct>(
0010           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))) {
0011   // Constructor implementation
0012   ptMin = iConfig.getUntrackedParameter<double>("MinPt", 20.);
0013   etaMin = iConfig.getUntrackedParameter<double>("MinEta", -2.4);
0014   etaMax = iConfig.getUntrackedParameter<double>("MaxEta", 2.4);
0015   drMin = iConfig.getUntrackedParameter<double>("drMin", 0.1);
0016   ptThreshold = iConfig.getUntrackedParameter<double>("ptThreshold", 2.);
0017 }
0018
0019 PhotonGenFilter::~PhotonGenFilter() {
0020   // Destructor implementation
0021 }
0022
0023 bool PhotonGenFilter::filter(edm::StreamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const {
0024   using namespace edm;
0025   Handle<HepMCProduct> evt;
0026   iEvent.getByToken(token_, evt);
0027   const HepMC::GenEvent *myGenEvent = evt->GetEvent();
0028   bool accepted_event = false;
0029
0030   for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0031        ++p) {  // loop through all particles
0032     if ((*p)->pdg_id() == 22 && !hasAncestor(*p, [](int x) {
0033           return x > 30 && x != 2212;
0034         })) {  // loop through all photons that don't come from hadrons
0035       if ((*p)->momentum().perp() > ptMin && (*p)->status() == 1 && (*p)->momentum().eta() > etaMin &&
0036           (*p)->momentum().eta() < etaMax) {  // check if photon passes pt and eta cuts
0037         bool good_photon = true;
0038         double phi = (*p)->momentum().phi();
0039         double eta = (*p)->momentum().eta();
0040         double pt = (*p)->momentum().perp();
0041         double frixione_isolation_coefficient = pt / (1 - cos(drMin));
0042         vector<double> particles_pt, particles_deltar;
0043         for (HepMC::GenEvent::particle_const_iterator q = myGenEvent->particles_begin();
0044              q != myGenEvent->particles_end();
0045              ++q) {        // loop through all particles to compute frixione isolation
0046           if (&p != &q) {  // don't compare the photon to itself
0047             if ((*q)->momentum().perp() > ptThreshold && (*q)->pdg_id() != 22 && (*q)->pdg_id() != 12 &&
0048                 (*q)->pdg_id() != 14 && (*q)->pdg_id() != 16 &&
0049                 (*q)->status() == 1) {  // check if particle passes pt and status cuts and is not a neutrino or photon
0050               double phi2 = (*q)->momentum().phi();
0051               double deltaphi = fabs(phi - phi2);
0052               if (deltaphi > M_PI)
0053                 deltaphi = 2. * M_PI - deltaphi;
0054               double eta2 = (*q)->momentum().eta();
0055               double deltaeta = fabs(eta - eta2);
0056               double deltaR = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
0057               if (deltaR < drMin)  // check if particle is within drMin of photon, for isolation
0058               {
0059                 particles_pt.push_back((*q)->momentum().perp());
0060                 particles_deltar.push_back(deltaR);
0061               }
0062             }
0063           }
0064         }
0065         //order particles_pt and particles_deltar according to increasing particles_deltar
0066         for (unsigned int i = 0; i < particles_deltar.size(); i++) {
0067           for (unsigned int j = i + 1; j < particles_deltar.size(); j++) {
0068             if (particles_deltar[i] > particles_deltar[j]) {
0069               double temp = particles_deltar[i];
0070               particles_deltar[i] = particles_deltar[j];
0071               particles_deltar[j] = temp;
0072               temp = particles_pt[i];
0073               particles_pt[i] = particles_pt[j];
0074               particles_pt[j] = temp;
0075             }
0076           }
0077         }
0078         //calculate frixione isolation
0079         double total_pt = 0;
0080         for (unsigned int i = 0; i < particles_deltar.size(); i++) {
0081           total_pt += particles_pt[i];
0082           if (total_pt > frixione_isolation_coefficient * (1 - cos(particles_deltar[i]))) {
0083             good_photon =
0084                 false;  // if for some delta R the isolation condition is not satisfied, the photon is not good
0085             break;
0086           }
0087         }
0088         if (good_photon) {
0089           if (hasAncestor(*p, [](int x) { return x == 11 || x == 13 || x == 15; }) ||
0090               hasAncestor(
0091                   *p,
0092                   [](int x) { return x == 24 || x == 5; },
0093                   true,
0094                   false)) {  // check if photon is from decay and defines the event as acceptable
0095             accepted_event = true;
0096           }
0097           if (!hasAncestor(*p, [](int x) { return x == 11 || x == 13 || x == 15; }) &&
0098               hasAncestor(
0099                   *p,
0100                   [](int x) { return x == 24 || x == 5; },
0101                   false,
0102                   true)) {  // check if photon comes from the hard process and discards it, in case it is
0103             return false;
0104           }
0105         }
0106       }
0107     }
0108   }
0109
0110   // Implementation for event filtering
0111   return accepted_event;  // Accept if it has found at least one good photon from decay and none from hard process
0112 }
0113
0114 void PhotonGenFilter::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0115   edm::ParameterSetDescription desc;
0116   desc.addUntracked<double>("MaxEta", 2.4);
0117   desc.addUntracked<double>("MinEta", -2.4);
0118   desc.addUntracked<double>("MinPt", 20.);
0119   desc.addUntracked<double>("drMin", 0.1);
0120   desc.addUntracked<double>("ptThreshold", 2.);
0121
0122   descriptions.add("PhotonGenFilter", desc);
0123 }
0124
0125 bool PhotonGenFilter::hasAncestor(
0126     HepMC::GenParticle *particle,
0127     function<bool(int)> check,
0128     bool isWorBFromDecayCheck,
0129     bool isWorBPromptCheck) const {  // function to check if a particle has a certain ancestor
0130   if (!particle)
0131     return false;
0132
0133   HepMC::GenVertex *prodVertex = particle->production_vertex();
0134
0135   // If the particle doesn't have a production vertex, it has no parents.
0136   if (!prodVertex)
0137     return false;
0138
0139   // Loop over all parents (incoming particles) of the vertex
0140   for (auto parent = prodVertex->particles_begin(HepMC::parents); parent != prodVertex->particles_end(HepMC::parents);
0141        ++parent) {
0142     int pdgId = abs((*parent)->pdg_id());
0143
0144     // Check if the PDG ID respects a check
0145     if (check(pdgId)) {
0146       if (isWorBFromDecayCheck) {
0147         return hasAncestor(
0148             *parent,
0149             [](int x) { return x == 6; },
0150             false,
0151             false);  // if the photon has a W or b quark ancestor, check if it comes from a top quark
0152       } else if (isWorBPromptCheck) {
0153         return !hasAncestor(
0154             *parent,
0155             [](int x) { return x == 6; },
0156             false,
0157             false);  // if the photon has a W or b quark ancestor, check that it doesn't come from a top quark decay (therefore it comes form the hard process)
0158       } else
0159         return true;
0160     }
0161
0162     return hasAncestor(
0163         *parent, check, isWorBFromDecayCheck, isWorBPromptCheck);  // Recursively check the ancestry of the parent
0164   }
0165
0166   return false;
0167 }
0168
0169 DEFINE_FWK_MODULE(PhotonGenFilter);