GenFilters/plugins/EMEnrichingFilterAlgo.cc

0001 #include "GeneratorInterface/GenFilters/plugins/EMEnrichingFilterAlgo.h"
0002 #include "DataFormats/Common/interface/Handle.h"
0003 #include "DataFormats/Math/interface/deltaR.h"
0004 #include "DataFormats/Math/interface/deltaPhi.h"
0005 #include "FWCore/Framework/interface/ConsumesCollector.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/Framework/interface/EventSetup.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Utilities/interface/InputTag.h"
0010
0011 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0012 #include "DataFormats/GeometrySurface/interface/Plane.h"
0013 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0014 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0015
0016 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0017 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0018 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0019
0020 #include <cmath>
0021 #include <cstdint>
0022 #include <cstdlib>
0023
0024 using namespace edm;
0025 using namespace std;
0026
0027 EMEnrichingFilterAlgo::EMEnrichingFilterAlgo(const edm::ParameterSet &iConfig, edm::ConsumesCollector &&iConsumes)
0028     //set constants
0029     : FILTER_TKISOCUT_(4),
0030       FILTER_CALOISOCUT_(7),
0031       FILTER_ETA_MIN_(0),
0032       FILTER_ETA_MAX_(2.4),
0033       ECALBARRELMAXETA_(1.479),
0034       ECALBARRELRADIUS_(129.0),
0035       ECALENDCAPZ_(304.5),
0036       isoGenParETMin_((float)iConfig.getParameter<double>("isoGenParETMin")),
0037       isoGenParConeSize_((float)iConfig.getParameter<double>("isoGenParConeSize")),
0038       clusterThreshold_((float)iConfig.getParameter<double>("clusterThreshold")),
0039       isoConeSize_((float)iConfig.getParameter<double>("isoConeSize")),
0040       hOverEMax_((float)iConfig.getParameter<double>("hOverEMax")),
0041       tkIsoMax_((float)iConfig.getParameter<double>("tkIsoMax")),
0042       caloIsoMax_((float)iConfig.getParameter<double>("caloIsoMax")),
0043       requireTrackMatch_(iConfig.getParameter<bool>("requireTrackMatch")),
0044       genParSourceToken_(
0045           iConsumes.consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("genParSource"))),
0046       magneticFieldToken_(iConsumes.esConsumes<MagneticField, IdealMagneticFieldRecord>()) {}
0047
0048 bool EMEnrichingFilterAlgo::filter(const edm::Event &iEvent, const edm::EventSetup &iSetup) const {
0049   Handle<reco::GenParticleCollection> genParsHandle;
0050   iEvent.getByToken(genParSourceToken_, genParsHandle);
0051   reco::GenParticleCollection genPars = *genParsHandle;
0052
0053   //bending of traj. of charged particles under influence of B-field
0054   std::vector<reco::GenParticle> genParsCurved = applyBFieldCurv(genPars, iSetup);
0055
0056   bool result1 = filterPhotonElectronSeed(
0057       clusterThreshold_, isoConeSize_, hOverEMax_, tkIsoMax_, caloIsoMax_, requireTrackMatch_, genPars, genParsCurved);
0058
0059   bool result2 = filterIsoGenPar(isoGenParETMin_, isoGenParConeSize_, genPars, genParsCurved);
0060
0061   bool result = result1 || result2;
0062
0063   return result;
0064 }
0065
0066 //filter that uses clustering around photon and electron seeds
0067 //only electrons, photons, charged pions, and charged kaons are clustered
0068 //additional requirements:
0069
0070 //seed threshold, total threshold, and cone size/shape are specified separately for the barrel and the endcap
0071 bool EMEnrichingFilterAlgo::filterPhotonElectronSeed(float clusterthreshold,
0072                                                      float isoConeSize,
0073                                                      float hOverEMax,
0074                                                      float tkIsoMax,
0075                                                      float caloIsoMax,
0076                                                      bool requiretrackmatch,
0077                                                      const std::vector<reco::GenParticle> &genPars,
0078                                                      const std::vector<reco::GenParticle> &genParsCurved) const {
0079   float seedthreshold = 5;
0080   float conesizeendcap = 15;
0081
0082   bool retval = false;
0083
0084   vector<reco::GenParticle> seeds;
0085   //find electron and photon seeds - must have E>seedthreshold GeV
0086   for (uint32_t is = 0; is < genParsCurved.size(); is++) {
0087     const reco::GenParticle &gp = genParsCurved.at(is);
0088     if (gp.status() != 1 || fabs(gp.eta()) > FILTER_ETA_MAX_ || fabs(gp.eta()) < FILTER_ETA_MIN_)
0089       continue;
0090     int absid = abs(gp.pdgId());
0091     if (absid != 11 && absid != 22)
0092       continue;
0093     if (gp.et() > seedthreshold)
0094       seeds.push_back(gp);
0095   }
0096
0097   bool matchtrack = false;
0098
0099   //for every seed, try to cluster stable particles about it in cone
0100   for (uint32_t is = 0; is < seeds.size(); is++) {
0101     float eTInCone = 0;   //eT associated to the electron cluster
0102     float tkIsoET = 0;    //tracker isolation energy
0103     float caloIsoET = 0;  //calorimeter isolation energy
0104     float hadET =
0105         0;  //isolation energy from heavy hadrons that goes in the same area as the "electron" - so contributes to H/E
0106     bool isBarrel = fabs(seeds.at(is).eta()) < ECALBARRELMAXETA_;
0107     for (uint32_t ig = 0; ig < genParsCurved.size(); ig++) {
0108       const reco::GenParticle &gp = genParsCurved.at(ig);
0109       const reco::GenParticle &gpUnCurv = genPars.at(ig);  //for tk isolation, p at vertex
0110       if (gp.status() != 1)
0111         continue;
0112       int gpabsid = abs(gp.pdgId());
0113       if (gp.et() < 1)
0114         continue;  //ignore very soft particles
0115       //BARREL
0116       if (isBarrel) {
0117         float dr = deltaR(seeds.at(is), gp);
0118         float dphi = deltaPhi(seeds.at(is).phi(), gp.phi());
0119         float deta = fabs(seeds.at(is).eta() - gp.eta());
0120         if (deta < 0.03 && dphi < 0.2) {
0121           if (gpabsid == 22 || gpabsid == 11 || gpabsid == 211 || gpabsid == 321) {
0122             //contributes to electron
0123             eTInCone += gp.et();
0124             //check for a matched track with at least 5 GeV
0125             if ((gpabsid == 11 || gpabsid == 211 || gpabsid == 321) && gp.et() > 5)
0126               matchtrack = true;
0127           } else {
0128             //contributes to H/E
0129             hadET += gp.et();
0130           }
0131         } else {
0132           float drUnCurv = deltaR(seeds.at(is), gpUnCurv);
0133           if ((gp.charge() == 0 && dr < isoConeSize && gpabsid != 22) || (gp.charge() != 0 && drUnCurv < isoConeSize)) {
0134             //contributes to calo isolation energy
0135             caloIsoET += gp.et();
0136           }
0137           if (gp.charge() != 0 && drUnCurv < isoConeSize) {
0138             //contributes to track isolation energy
0139             tkIsoET += gp.et();
0140           }
0141         }
0142         //ENDCAP
0143       } else {
0144         float drxy = deltaRxyAtEE(seeds.at(is), gp);
0145         float dr = deltaR(seeds.at(is), gp);  //the isolation is done in dR
0146         if (drxy < conesizeendcap) {
0147           if (gpabsid == 22 || gpabsid == 11 || gpabsid == 211 || gpabsid == 321) {
0148             //contributes to electron
0149             eTInCone += gp.et();
0150             //check for a matched track with at least 5 GeV
0151             if ((gpabsid == 11 || gpabsid == 211 || gpabsid == 321) && gp.et() > 5)
0152               matchtrack = true;
0153           } else {
0154             //contributes to H/E
0155             hadET += gp.et();
0156           }
0157         } else {
0158           float drUnCurv = deltaR(seeds.at(is), gpUnCurv);
0159           if ((gp.charge() == 0 && dr < isoConeSize && gpabsid != 22) || (gp.charge() != 0 && drUnCurv < isoConeSize)) {
0160             //contributes to calo isolation energy
0161             caloIsoET += gp.et();
0162           }
0163           if (gp.charge() != 0 && drUnCurv < isoConeSize) {
0164             //contributes to track isolation energy
0165             tkIsoET += gp.et();
0166           }
0167         }
0168       }
0169     }
0170
0171     if (eTInCone > clusterthreshold && (!requiretrackmatch || matchtrack)) {
0172       //       cout <<"isoET: "<<isoET<<endl;
0173       if (hadET / eTInCone < hOverEMax && tkIsoET < tkIsoMax && caloIsoET < caloIsoMax)
0174         retval = true;
0175       break;
0176     }
0177   }
0178
0179   return retval;
0180 }
0181
0182 //make new genparticles vector taking into account the bending of charged particles in the b field
0183 //only stable-final-state (status==1) particles, with ET>=1 GeV, have their trajectories bent
0184 std::vector<reco::GenParticle> EMEnrichingFilterAlgo::applyBFieldCurv(const std::vector<reco::GenParticle> &genPars,
0185                                                                       const edm::EventSetup &iSetup) const {
0186   vector<reco::GenParticle> curvedPars;
0187
0188   MagneticField const *magField = &iSetup.getData(magneticFieldToken_);
0189
0190   Cylinder::CylinderPointer theBarrel =
0191       Cylinder::build(Surface::PositionType(0, 0, 0), Surface::RotationType(), ECALBARRELRADIUS_);
0192   Plane::PlanePointer endCapPlus = Plane::build(Surface::PositionType(0, 0, ECALENDCAPZ_), Surface::RotationType());
0193   Plane::PlanePointer endCapMinus =
0194       Plane::build(Surface::PositionType(0, 0, -1 * ECALENDCAPZ_), Surface::RotationType());
0195
0196   AnalyticalPropagator propagator(magField, alongMomentum);
0197
0198   for (uint32_t ig = 0; ig < genPars.size(); ig++) {
0199     const reco::GenParticle &gp = genPars.at(ig);
0200     //don't bend trajectories of neutral particles, unstable particles, particles with < 1 GeV
0201     //particles with < ~0.9 GeV don't reach the barrel
0202     //so just put them as-is into the new vector
0203     if (gp.charge() == 0 || gp.status() != 1 || gp.et() < 1) {
0204       curvedPars.push_back(gp);
0205       continue;
0206     }
0207     GlobalPoint vertex(gp.vx(), gp.vy(), gp.vz());
0208     GlobalVector gvect(gp.px(), gp.py(), gp.pz());
0209     FreeTrajectoryState fts(vertex, gvect, gp.charge(), magField);
0210     TrajectoryStateOnSurface impact;
0211     //choose to propagate to barrel, +Z endcap, or -Z endcap, according to eta
0212     if (fabs(gp.eta()) < ECALBARRELMAXETA_) {
0213       impact = propagator.propagate(fts, *theBarrel);
0214     } else if (gp.eta() > 0) {
0215       impact = propagator.propagate(fts, *endCapPlus);
0216     } else {
0217       impact = propagator.propagate(fts, *endCapMinus);
0218     }
0219     //in case the particle doesn't reach the barrel/endcap, just put it as-is into the new vector
0220     //it should reach though.
0221     if (!impact.isValid()) {
0222       curvedPars.push_back(gp);
0223       continue;
0224     }
0225     math::XYZTLorentzVector newp4;
0226
0227     //the magnitude of p doesn't change, only the direction
0228     //NB I do get some small change in magnitude by the following,
0229     //I think it is a numerical precision issue
0230     float et = gp.et();
0231     float phinew = impact.globalPosition().phi();
0232     float pxnew = et * cos(phinew);
0233     float pynew = et * sin(phinew);
0234     newp4.SetPx(pxnew);
0235     newp4.SetPy(pynew);
0236     newp4.SetPz(gp.pz());
0237     newp4.SetE(gp.energy());
0238     reco::GenParticle gpnew = gp;
0239     gpnew.setP4(newp4);
0240     curvedPars.push_back(gpnew);
0241   }
0242   return curvedPars;
0243 }
0244
0245 //calculate the difference in the xy-plane positions of gp1 and gp1 at the ECAL endcap
0246 //if they go in different z directions returns 9999.
0247 float EMEnrichingFilterAlgo::deltaRxyAtEE(const reco::GenParticle &gp1, const reco::GenParticle &gp2) const {
0248   if (gp1.pz() * gp2.pz() < 0)
0249     return 9999.;
0250
0251   float rxy1 = ECALENDCAPZ_ * tan(gp1.theta());
0252   float x1 = cos(gp1.phi()) * rxy1;
0253   float y1 = sin(gp1.phi()) * rxy1;
0254
0255   float rxy2 = ECALENDCAPZ_ * tan(gp2.theta());
0256   float x2 = cos(gp2.phi()) * rxy2;
0257   float y2 = sin(gp2.phi()) * rxy2;
0258
0259   float dxy = sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1));
0260   return dxy;
0261 }
0262
0263 //filter looking for isolated charged pions, charged kaons, and electrons.
0264 //isolation done in cone of given size, looking at charged particles and neutral hadrons
0265 //photons aren't counted in the isolation requirements
0266
0267 //need to have both the the curved and uncurved genpar collections
0268 //because tracker iso has to be treated differently than calo iso
0269 bool EMEnrichingFilterAlgo::filterIsoGenPar(float etmin,
0270                                             float conesize,
0271                                             const reco::GenParticleCollection &gph,
0272                                             const reco::GenParticleCollection &gphCurved) const {
0273   for (uint32_t ip = 0; ip < gph.size(); ip++) {
0274     const reco::GenParticle &gp = gph.at(ip);
0275     const reco::GenParticle &gpCurved = gphCurved.at(ip);
0276     int gpabsid = abs(gp.pdgId());
0277     //find potential candidates
0278     if (gp.et() <= etmin || gp.status() != 1)
0279       continue;
0280     if (gpabsid != 11 && gpabsid != 211 && gpabsid != 321)
0281       continue;
0282     if (fabs(gp.eta()) < FILTER_ETA_MIN_)
0283       continue;
0284     if (fabs(gp.eta()) > FILTER_ETA_MAX_)
0285       continue;
0286
0287     //check isolation
0288     float tkiso = 0;
0289     float caloiso = 0;
0290     for (uint32_t jp = 0; jp < gph.size(); jp++) {
0291       if (jp == ip)
0292         continue;
0293       const reco::GenParticle &pp = gph.at(jp);
0294       const reco::GenParticle &ppCurved = gphCurved.at(jp);
0295       if (pp.status() != 1)
0296         continue;
0297       float dr = deltaR(gp, pp);
0298       float drCurved = deltaR(gpCurved, ppCurved);
0299       if (abs(pp.charge()) == 1 && pp.et() > 2 && dr < conesize) {
0300         tkiso += pp.et();
0301       }
0302       if (pp.et() > 2 && abs(pp.pdgId()) != 22 && drCurved < conesize) {
0303         caloiso += pp.energy();
0304       }
0305     }
0306     if (tkiso < FILTER_TKISOCUT_ && caloiso < FILTER_CALOISOCUT_)
0307       return true;
0308   }
0309   return false;
0310 }