Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:29:58

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/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 }  // namespace
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),       // delta_eta, delta_phi
0096       deltaEE(2.93 / 317 / 2 * 5) {}  // delta_x/z, delta_y/z
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   //  std::cout<<" Number of seeds "<<seeds.size()<<" ProccessID "<<myGenEvent->signal_process_id()<<std::endl;
0116   //  std::cout<<" ParticleId 7= "<<myGenEvent->particle(7)->pdg_id()
0117   //  <<" pT "<<myGenEvent->particle(7)->momentum().perp()
0118   //  <<" Eta "<<myGenEvent->particle(7)->momentum().eta()
0119   //  <<" Phi "<<myGenEvent->particle(7)->momentum().phi()<<std::endl;
0120   //  std::cout<<" ParticleId 8= "<<myGenEvent->particle(8)->pdg_id()<<" pT "<<myGenEvent->particle(8)->momentum().perp()<<" Eta "<<myGenEvent->particle(8)->momentum().eta()<<" Phi "<<myGenEvent->particle(8)->momentum().phi()<<std::endl;
0121 
0122   for (std::list<const HepMC::GenParticle*>::const_iterator is = seeds.begin(); is != seeds.end(); is++) {
0123     double etaPhoton = (*is)->momentum().eta();
0124     double phiPhoton = (*is)->momentum().phi();
0125 
0126     /*
0127     double dphi7=std::abs(deltaPhi(phiPhoton,
0128                    myGenEvent->particle(7)->momentum().phi()));
0129     double dphi=std::abs(deltaPhi(phiPhoton,
0130                   myGenEvent->particle(8)->momentum().phi()));
0131     */
0132 
0133     HepMC::GenEvent::particle_const_iterator ppp = myGenEvent->particles_begin();
0134     for (int i = 0; i < 6; ++i)
0135       ppp++;
0136     HepMC::GenParticle* particle7 = (*ppp);
0137     ppp++;
0138     HepMC::GenParticle* particle8 = (*ppp);
0139 
0140     double dphi7 = std::abs(deltaPhi(phiPhoton, particle7->momentum().phi()));
0141     double dphi = std::abs(deltaPhi(phiPhoton, particle8->momentum().phi()));
0142 
0143     int jetline = 8;
0144     if (dphi7 > dphi) {
0145       dphi = dphi7;
0146       jetline = 7;
0147     }
0148 
0149     //    std::cout<<" Dphi "<<dphi<<" "<<dphiMin<<std::endl;
0150     //    if(dphi<dphiMin) {
0151     //      std::cout<<"Reject dphi"<<std::endl;
0152     //      continue;
0153     //    }
0154 
0155     //    double etaJet= myGenEvent->particle(jetline)->momentum().eta();
0156     double etaJet = 0.0;
0157     if (jetline == 8)
0158       etaJet = particle8->momentum().eta();
0159     else
0160       etaJet = particle7->momentum().eta();
0161 
0162     double eta1 = etaJet - detaMax;
0163     double eta2 = etaJet + detaMax;
0164     if (eta1 > etaPhotonCut2)
0165       eta1 = etaPhotonCut2;
0166     if (eta2 < -etaPhotonCut2)
0167       eta2 = -etaPhotonCut2;
0168 
0169     //    std::cout<<" Etaphoton "<<etaPhoton<<" "<<eta1<<" "<<eta2<<std::endl;
0170 
0171     if (etaPhoton < eta1 || etaPhoton > eta2) {
0172       //       std::cout<<"Reject eta"<<std::endl;
0173       continue;
0174     }
0175     bool inEB(false);
0176     double tgx(0);
0177     double tgy(0);
0178     if (std::abs(etaPhoton) < ebEtaMax)
0179       inEB = true;
0180     else {
0181       tgx = (*is)->momentum().px() / (*is)->momentum().pz();
0182       tgy = (*is)->momentum().py() / (*is)->momentum().pz();
0183     }
0184 
0185     double etPhoton = 0;
0186     double etPhotonCharged = 0;
0187     double etCone = 0;
0188     double etConeCharged = 0;
0189     double ptMaxHadron = 0;
0190 
0191     for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0192          ++p) {
0193       if ((*p)->status() != 1)
0194         continue;
0195       int pid = (*p)->pdg_id();
0196       int apid = std::abs(pid);
0197       if (apid > 11 && apid < 21)
0198         continue;  //get rid of muons and neutrinos
0199       double eta = (*p)->momentum().eta();
0200       double phi = (*p)->momentum().phi();
0201       if (deltaR2(etaPhoton, phiPhoton, eta, phi) > cone * cone)
0202         continue;
0203       double pt = (*p)->momentum().perp();
0204 
0205       //       double charge=(*p)->particledata().charge();
0206       //      int charge3=(*p)->particleID().threeCharge();
0207 
0208       //***
0209       ParticleDataTable const& pdt = iSetup.getData(particleDataTableToken_);
0210 
0211       int charge3 = ((pdt.particle((*p)->pdg_id()))->ID().threeCharge());
0212       //***
0213 
0214       etCone += pt;
0215       if (charge3 && pt < 2)
0216         etConeCharged += pt;
0217 
0218       //select particles matching a crystal array centered on photon
0219       if (inEB) {
0220         if (std::abs(eta - etaPhoton) > deltaEB || std::abs(deltaPhi(phi, phiPhoton)) > deltaEB)
0221           continue;
0222       } else if (std::abs((*p)->momentum().px() / (*p)->momentum().pz() - tgx) > deltaEE ||
0223                  std::abs((*p)->momentum().py() / (*p)->momentum().pz() - tgy) > deltaEE)
0224         continue;
0225 
0226       etPhoton += pt;
0227       if (charge3 && pt < 2)
0228         etPhotonCharged += pt;
0229       if (apid > 100 && apid != 310 && pt > ptMaxHadron)
0230         ptMaxHadron = pt;
0231     }
0232     //    std::cout<<" etPhoton "<<etPhoton<<" "<<ptMin<<" "<<ptMax<<std::endl;
0233     if (etPhoton < ptMin || etPhoton > ptMax) {
0234       //     std::cout<<" Reject etPhoton "<<std::endl;
0235       continue;
0236     }
0237     //isolation cuts
0238 
0239     //    double isocut1 = 5+etPhoton/20-etPhoton*etPhoton/1e4;
0240     //    double isocut2 = 3+etPhoton/20-etPhoton*etPhoton*etPhoton/1e6;
0241     //    double isocut3 = 4.5+etPhoton/40;
0242     //    if (etPhoton>165.)
0243     //    {
0244     //      isocut1 = 5.+165./20.-165.*165./1e4;
0245     //      isocut2 = 3.+165./20.-165.*165.*165./1e6;
0246     //      isocut3 = 4.5+165./40.;
0247     //    }
0248     //    std::cout<<" etCone "<<etCone<<" "<<etPhoton<<" "<<etCone-etPhoton<<" "<<isocut1<<std::endl;
0249     //    std::cout<<"Second cut on iso "<<etCone-etPhoton-(etConeCharged-etPhotonCharged)<<" cut value "<<isocut2<<" etPhoton "<<etPhoton<<std::endl;
0250     //    std::cout<<" PtHadron "<<ptMaxHadron<<" "<<4.5+etPhoton/40<<std::endl;
0251     //    if(etCone-etPhoton> 5+etPhoton/20-etPhoton*etPhoton/1e4) continue;
0252     //    std::cout<<" Accept 1"<<std::endl;
0253     //    if(etCone-etPhoton-(etConeCharged-etPhotonCharged) >
0254     //       isocut2) continue;
0255     //     std::cout<<" Accept 2"<<std::endl;
0256     //    if(ptMaxHadron > isocut3) continue;
0257     //    std::cout<<"Accept event "<<std::endl;
0258 
0259     accepted = true;
0260     break;
0261 
0262   }  //loop over seeds
0263   return accepted;
0264 }
0265 
0266 DEFINE_FWK_MODULE(PythiaFilterGammaJetWithOutBg);