GenFilters/plugins/PythiaFilterGammaJet.cc

0001 /** \class PythiaFilterGammaJet
0002  *
0003  *  PythiaFilterGammaJet filter implements generator-level preselections
0004  *  for photon+jet like events to be used in jet energy calibration.
0005  *  Ported from fortran code written by V.Konoplianikov.
0006  *
0007  * \author A.Ulyanov, ITEP
0008  *
0009  ************************************************************/
0010
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "FWCore/Utilities/interface/ESGetToken.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 #include "FWCore/Framework/interface/Frameworkfwd.h"
0016 #include "FWCore/Framework/interface/global/EDFilter.h"
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/Utilities/interface/EDGetToken.h"
0020 #include "FWCore/Utilities/interface/InputTag.h"
0021 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0022 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0023 #include "SimGeneral/HepPDTRecord/interface/PDTRecord.h"
0024
0025 #include <cmath>
0026 #include <cstdlib>
0027 #include <list>
0028 #include <string>
0029
0030 class PythiaFilterGammaJet : public edm::global::EDFilter<> {
0031 public:
0032   explicit PythiaFilterGammaJet(const edm::ParameterSet&);
0033
0034   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0035
0036 private:
0037   const edm::EDGetTokenT<edm::HepMCProduct> token_;
0038   edm::ESGetToken<ParticleDataTable, PDTRecord> particleDataTableToken_;
0039   const double etaMax;
0040   const double ptSeed;
0041   const double ptMin;
0042   const double ptMax;
0043   const double dphiMin;
0044   const double detaMax;
0045   const double etaPhotonCut2;
0046
0047   const double cone;
0048   const double ebEtaMax;
0049   const double deltaEB;
0050   const double deltaEE;
0051 };
0052
0053 namespace {
0054
0055   double deltaR2(double eta0, double phi0, double eta, double phi) {
0056     double dphi = phi - phi0;
0057     if (dphi > M_PI)
0058       dphi -= 2 * M_PI;
0059     else if (dphi <= -M_PI)
0060       dphi += 2 * M_PI;
0061     return dphi * dphi + (eta - eta0) * (eta - eta0);
0062   }
0063
0064   double deltaPhi(double phi0, double phi) {
0065     double dphi = phi - phi0;
0066     if (dphi > M_PI)
0067       dphi -= 2 * M_PI;
0068     else if (dphi <= -M_PI)
0069       dphi += 2 * M_PI;
0070     return dphi;
0071   }
0072
0073   class ParticlePtGreater {
0074   public:
0075     int operator()(const HepMC::GenParticle* p1, const HepMC::GenParticle* p2) const {
0076       return p1->momentum().perp() > p2->momentum().perp();
0077     }
0078   };
0079 }  // namespace
0080
0081 PythiaFilterGammaJet::PythiaFilterGammaJet(const edm::ParameterSet& iConfig)
0082     : token_(consumes<edm::HepMCProduct>(
0083           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0084       particleDataTableToken_(esConsumes<ParticleDataTable, PDTRecord>()),
0085       etaMax(iConfig.getUntrackedParameter<double>("MaxPhotonEta", 2.8)),
0086       ptSeed(iConfig.getUntrackedParameter<double>("PhotonSeedPt", 5.)),
0087       ptMin(iConfig.getUntrackedParameter<double>("MinPhotonPt")),
0088       ptMax(iConfig.getUntrackedParameter<double>("MaxPhotonPt")),
0089       dphiMin(iConfig.getUntrackedParameter<double>("MinDeltaPhi", -1) / 180 * M_PI),
0090       detaMax(iConfig.getUntrackedParameter<double>("MaxDeltaEta", 10.)),
0091       etaPhotonCut2(iConfig.getUntrackedParameter<double>("MinPhotonEtaForwardJet", 1.3)),
0092       cone(0.5),
0093       ebEtaMax(1.479),
0094       deltaEB(0.01745 / 2 * 5),       // delta_eta, delta_phi
0095       deltaEE(2.93 / 317 / 2 * 5) {}  // delta_x/z, delta_y/z
0096
0097 bool PythiaFilterGammaJet::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0098   bool accepted = false;
0099   edm::Handle<edm::HepMCProduct> evt;
0100   iEvent.getByToken(token_, evt);
0101
0102   std::list<const HepMC::GenParticle*> seeds;
0103   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0104
0105   if (myGenEvent->signal_process_id() == 14 || myGenEvent->signal_process_id() == 29) {
0106     for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0107          ++p) {
0108       if ((*p)->pdg_id() == 22 && (*p)->status() == 1 && (*p)->momentum().perp() > ptSeed &&
0109           std::abs((*p)->momentum().eta()) < etaMax)
0110         seeds.push_back(*p);
0111     }
0112
0113     seeds.sort(ParticlePtGreater());
0114     for (std::list<const HepMC::GenParticle*>::const_iterator is = seeds.begin(); is != seeds.end(); is++) {
0115       double etaPhoton = (*is)->momentum().eta();
0116       double phiPhoton = (*is)->momentum().phi();
0117
0118       HepMC::GenEvent::particle_const_iterator ppp = myGenEvent->particles_begin();
0119       for (int i = 0; i < 6; ++i)
0120         ppp++;
0121       HepMC::GenParticle* particle7 = (*ppp);
0122       ppp++;
0123       HepMC::GenParticle* particle8 = (*ppp);
0124
0125       double dphi7 = std::abs(deltaPhi(phiPhoton, particle7->momentum().phi()));
0126       double dphi = std::abs(deltaPhi(phiPhoton, particle8->momentum().phi()));
0127       int jetline = 8;
0128       if (dphi7 > dphi) {
0129         dphi = dphi7;
0130         jetline = 7;
0131       }
0132       if (dphi < dphiMin)
0133         continue;
0134       //double etaJet= myGenEvent->particle(jetline)->momentum().eta();
0135       double etaJet = 0.0;
0136       if (jetline == 8)
0137         etaJet = particle8->momentum().eta();
0138       else
0139         etaJet = particle7->momentum().eta();
0140
0141       double eta1 = etaJet - detaMax;
0142       double eta2 = etaJet + detaMax;
0143       if (eta1 > etaPhotonCut2)
0144         eta1 = etaPhotonCut2;
0145       if (eta2 < -etaPhotonCut2)
0146         eta2 = -etaPhotonCut2;
0147       if (etaPhoton < eta1 || etaPhoton > eta2)
0148         continue;
0149
0150       bool inEB(false);
0151       double tgx(0);
0152       double tgy(0);
0153       if (std::abs(etaPhoton) < ebEtaMax)
0154         inEB = true;
0155       else {
0156         tgx = (*is)->momentum().px() / (*is)->momentum().pz();
0157         tgy = (*is)->momentum().py() / (*is)->momentum().pz();
0158       }
0159
0160       double etPhoton = 0;
0161       double etPhotonCharged = 0;
0162       double etCone = 0;
0163       double etConeCharged = 0;
0164       double ptMaxHadron = 0;
0165
0166       for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0167            ++p) {
0168         if ((*p)->status() != 1)
0169           continue;
0170         int pid = (*p)->pdg_id();
0171         int apid = std::abs(pid);
0172         if (apid > 11 && apid < 20)
0173           continue;  //get rid of muons and neutrinos
0174         double eta = (*p)->momentum().eta();
0175         double phi = (*p)->momentum().phi();
0176         if (deltaR2(etaPhoton, phiPhoton, eta, phi) > cone * cone)
0177           continue;
0178         double pt = (*p)->momentum().perp();
0179
0180         ParticleDataTable const& pdt = iSetup.getData(particleDataTableToken_);
0181
0182         //       double charge=(*p)->particledata().charge();
0183         //int charge3=(*p)->particleID().threeCharge();
0184
0185         int charge3 = ((pdt.particle((*p)->pdg_id()))->ID().threeCharge());
0186         etCone += pt;
0187         if (charge3 && pt < 2)
0188           etConeCharged += pt;
0189
0190         //select particles matching a crystal array centered on photon
0191         if (inEB) {
0192           if (std::abs(eta - etaPhoton) > deltaEB || std::abs(deltaPhi(phi, phiPhoton)) > deltaEB)
0193             continue;
0194         } else if (std::abs((*p)->momentum().px() / (*p)->momentum().pz() - tgx) > deltaEE ||
0195                    std::abs((*p)->momentum().py() / (*p)->momentum().pz() - tgy) > deltaEE)
0196           continue;
0197
0198         etPhoton += pt;
0199         if (charge3 && pt < 2)
0200           etPhotonCharged += pt;
0201         if (apid > 100 && apid != 310 && pt > ptMaxHadron)
0202           ptMaxHadron = pt;
0203       }
0204
0205       if (etPhoton < ptMin || etPhoton > ptMax)
0206         continue;
0207
0208       //isolation cuts
0209       if (etCone - etPhoton > 5 + etPhoton / 20 - etPhoton * etPhoton / 1e4)
0210         continue;
0211       if (etCone - etPhoton - (etConeCharged - etPhotonCharged) >
0212           3 + etPhoton / 20 - etPhoton * etPhoton * etPhoton / 1e6)
0213         continue;
0214       if (ptMaxHadron > 4.5 + etPhoton / 40)
0215         continue;
0216
0217       accepted = true;
0218       break;
0219
0220     }  //loop over seeds
0221
0222   } else {
0223     // end of if(gammajetevent)
0224     return true;
0225     // accept all non-gammajet events
0226   }
0227   return accepted;
0228 }
0229
0230 DEFINE_FWK_MODULE(PythiaFilterGammaJet);