File indexing completed on 2024-04-06 12:13:41
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/Framework/interface/Frameworkfwd.h"
0015 #include "FWCore/Framework/interface/global/EDFilter.h"
0016 #include "FWCore/Framework/interface/MakerMacros.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/Utilities/interface/EDGetToken.h"
0019 #include "FWCore/Utilities/interface/ESGetToken.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 PythiaFilterGammaJetWithOutBg : public edm::global::EDFilter<> {
0031 public:
0032 explicit PythiaFilterGammaJetWithOutBg(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 const edm::ESGetToken<ParticleDataTable, PDTRecord> particleDataTableToken_;
0039
0040 const double etaMax;
0041 const double ptSeed;
0042 const double ptMin;
0043 const double ptMax;
0044 const double dphiMin;
0045 const double detaMax;
0046 const double etaPhotonCut2;
0047
0048 const double cone;
0049 const double ebEtaMax;
0050 const double deltaEB;
0051 const double deltaEE;
0052 };
0053
0054 namespace {
0055
0056 double deltaR2(double eta0, double phi0, double eta, double phi) {
0057 double dphi = phi - phi0;
0058 if (dphi > M_PI)
0059 dphi -= 2 * M_PI;
0060 else if (dphi <= -M_PI)
0061 dphi += 2 * M_PI;
0062 return dphi * dphi + (eta - eta0) * (eta - eta0);
0063 }
0064
0065 double deltaPhi(double phi0, double phi) {
0066 double dphi = phi - phi0;
0067 if (dphi > M_PI)
0068 dphi -= 2 * M_PI;
0069 else if (dphi <= -M_PI)
0070 dphi += 2 * M_PI;
0071 return dphi;
0072 }
0073
0074 class ParticlePtGreater {
0075 public:
0076 int operator()(const HepMC::GenParticle* p1, const HepMC::GenParticle* p2) const {
0077 return p1->momentum().perp() > p2->momentum().perp();
0078 }
0079 };
0080 }
0081
0082 PythiaFilterGammaJetWithOutBg::PythiaFilterGammaJetWithOutBg(const edm::ParameterSet& iConfig)
0083 : token_(consumes<edm::HepMCProduct>(
0084 edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0085 particleDataTableToken_(esConsumes<ParticleDataTable, PDTRecord>()),
0086 etaMax(iConfig.getUntrackedParameter<double>("MaxPhotonEta", 2.8)),
0087 ptSeed(iConfig.getUntrackedParameter<double>("PhotonSeedPt", 5.)),
0088 ptMin(iConfig.getUntrackedParameter<double>("MinPhotonPt")),
0089 ptMax(iConfig.getUntrackedParameter<double>("MaxPhotonPt")),
0090 dphiMin(iConfig.getUntrackedParameter<double>("MinDeltaPhi", -1) / 180 * M_PI),
0091 detaMax(iConfig.getUntrackedParameter<double>("MaxDeltaEta", 10.)),
0092 etaPhotonCut2(iConfig.getUntrackedParameter<double>("MinPhotonEtaForwardJet", 1.3)),
0093 cone(0.5),
0094 ebEtaMax(1.479),
0095 deltaEB(0.01745 / 2 * 5),
0096 deltaEE(2.93 / 317 / 2 * 5) {}
0097
0098 bool PythiaFilterGammaJetWithOutBg::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0099 bool accepted = false;
0100 edm::Handle<edm::HepMCProduct> evt;
0101 iEvent.getByToken(token_, evt);
0102
0103 std::list<const HepMC::GenParticle*> seeds;
0104 const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0105
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
0115 for (std::list<const HepMC::GenParticle*>::const_iterator is = seeds.begin(); is != seeds.end(); is++) {
0116 double etaPhoton = (*is)->momentum().eta();
0117 double phiPhoton = (*is)->momentum().phi();
0118
0119 HepMC::GenEvent::particle_const_iterator ppp = myGenEvent->particles_begin();
0120 for (int i = 0; i < 6; ++i)
0121 ppp++;
0122 HepMC::GenParticle* particle7 = (*ppp);
0123 ppp++;
0124 HepMC::GenParticle* particle8 = (*ppp);
0125
0126 double dphi7 = std::abs(deltaPhi(phiPhoton, particle7->momentum().phi()));
0127 double dphi8 = std::abs(deltaPhi(phiPhoton, particle8->momentum().phi()));
0128
0129 double etaJet = dphi7 > dphi8 ? particle7->momentum().eta() : particle8->momentum().eta();
0130
0131 double eta1 = etaJet - detaMax;
0132 double eta2 = etaJet + detaMax;
0133 if (eta1 > etaPhotonCut2)
0134 eta1 = etaPhotonCut2;
0135 if (eta2 < -etaPhotonCut2)
0136 eta2 = -etaPhotonCut2;
0137
0138 if (etaPhoton < eta1 || etaPhoton > eta2) {
0139 continue;
0140 }
0141 bool inEB(false);
0142 double tgx(0);
0143 double tgy(0);
0144 if (std::abs(etaPhoton) < ebEtaMax)
0145 inEB = true;
0146 else {
0147 tgx = (*is)->momentum().px() / (*is)->momentum().pz();
0148 tgy = (*is)->momentum().py() / (*is)->momentum().pz();
0149 }
0150
0151 double etPhoton = 0;
0152 double ptMaxHadron = 0;
0153
0154 for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0155 ++p) {
0156 if ((*p)->status() != 1)
0157 continue;
0158 int apid = std::abs((*p)->pdg_id());
0159 if (apid > 11 && apid < 21)
0160 continue;
0161 double eta = (*p)->momentum().eta();
0162 double phi = (*p)->momentum().phi();
0163 if (deltaR2(etaPhoton, phiPhoton, eta, phi) > cone * cone)
0164 continue;
0165 double pt = (*p)->momentum().perp();
0166
0167
0168 if (inEB) {
0169 if (std::abs(eta - etaPhoton) > deltaEB || std::abs(deltaPhi(phi, phiPhoton)) > deltaEB)
0170 continue;
0171 } else if (std::abs((*p)->momentum().px() / (*p)->momentum().pz() - tgx) > deltaEE ||
0172 std::abs((*p)->momentum().py() / (*p)->momentum().pz() - tgy) > deltaEE)
0173 continue;
0174
0175 etPhoton += pt;
0176 if (apid > 100 && apid != 310 && pt > ptMaxHadron)
0177 ptMaxHadron = pt;
0178 }
0179
0180 if (etPhoton < ptMin || etPhoton > ptMax) {
0181 continue;
0182 }
0183
0184 accepted = true;
0185 break;
0186
0187 }
0188 return accepted;
0189 }
0190
0191 DEFINE_FWK_MODULE(PythiaFilterGammaJetWithOutBg);