CMS 3D CMS Logo

Public Member Functions | Private Member Functions | Private Attributes

doubleEMEnrichingFilterAlgo Class Reference

#include <doubleEMEnrichingFilterAlgo.h>

List of all members.

Public Member Functions

 doubleEMEnrichingFilterAlgo (const edm::ParameterSet &)
bool filter (const edm::Event &iEvent, const edm::EventSetup &iSetup)
bool hasBCAncestors (reco::GenParticle gp)
 ~doubleEMEnrichingFilterAlgo ()

Private Member Functions

std::vector< reco::GenParticleapplyBFieldCurv (const std::vector< reco::GenParticle > &genPars, const edm::EventSetup &iSetup)
float deltaRxyAtEE (const reco::GenParticle &gp1, const reco::GenParticle &gp2)
int filterIsoGenPar (float etmin, float conesize, const reco::GenParticleCollection &gph, const reco::GenParticleCollection &gphCurved)
int filterPhotonElectronSeed (float clusterthreshold, float seedthreshold, float isoConeSize, float hOverEMax, float tkIsoMax, float caloIsoMax, bool requiretrackmatch, const std::vector< reco::GenParticle > &genPars, const std::vector< reco::GenParticle > &genParsCurved)
bool isBCBaryon (reco::GenParticle gp)
bool isBCHadron (reco::GenParticle gp)
bool isBCMeson (reco::GenParticle gp)

Private Attributes

float caloIsoMax_
float clusterThreshold_
float ECALBARRELMAXETA_
float ECALBARRELRADIUS_
float ECALENDCAPZ_
float eTThreshold_
float FILTER_CALOISOCUT_
float FILTER_ETA_MAX_
float FILTER_ETA_MIN_
float FILTER_TKISOCUT_
edm::InputTag genParSource_
float hOverEMax_
float isoConeSize_
float isoGenParConeSize_
float isoGenParETMin_
bool requireTrackMatch_
float seedThreshold_
std::vector< reco::GenParticlesel1seeds
std::vector< reco::GenParticlesel2seeds
std::vector< reco::GenParticleselBCtoEseeds
float tkIsoMax_

Detailed Description

Definition at line 26 of file doubleEMEnrichingFilterAlgo.h.


Constructor & Destructor Documentation

doubleEMEnrichingFilterAlgo::doubleEMEnrichingFilterAlgo ( const edm::ParameterSet iConfig)

Definition at line 26 of file doubleEMEnrichingFilterAlgo.cc.

References edm::ParameterSet::getParameter().

                                                                                       { 

  //set constants
  FILTER_TKISOCUT_=2;
  FILTER_CALOISOCUT_=2;
  FILTER_ETA_MIN_=0;
  FILTER_ETA_MAX_=2.5;
  ECALBARRELMAXETA_=1.479;
  ECALBARRELRADIUS_=129.0;
  ECALENDCAPZ_=304.5;

  // from bctoe
  eTThreshold_=(float)iConfig.getParameter<double>("eTThreshold");
  
  isoGenParETMin_=(float) iConfig.getParameter<double>("isoGenParETMin");
  isoGenParConeSize_=(float) iConfig.getParameter<double>("isoGenParConeSize");
  clusterThreshold_=(float) iConfig.getParameter<double>("clusterThreshold");
  seedThreshold_=(float) iConfig.getParameter<double>("seedThreshold");
  isoConeSize_=(float) iConfig.getParameter<double>("isoConeSize");
  hOverEMax_=(float) iConfig.getParameter<double>("hOverEMax");
  tkIsoMax_=(float) iConfig.getParameter<double>("tkIsoMax");
  caloIsoMax_=(float) iConfig.getParameter<double>("caloIsoMax");
  requireTrackMatch_=iConfig.getParameter<bool>("requireTrackMatch");
  genParSource_=iConfig.getParameter<edm::InputTag>("genParSource");

}
doubleEMEnrichingFilterAlgo::~doubleEMEnrichingFilterAlgo ( )

Definition at line 53 of file doubleEMEnrichingFilterAlgo.cc.

                                                          {
}

Member Function Documentation

std::vector< reco::GenParticle > doubleEMEnrichingFilterAlgo::applyBFieldCurv ( const std::vector< reco::GenParticle > &  genPars,
const edm::EventSetup iSetup 
) [private]

Definition at line 248 of file doubleEMEnrichingFilterAlgo.cc.

References alongMomentum, newFWLiteAna::build, reco::LeafCandidate::charge(), funct::cos(), reco::LeafCandidate::energy(), reco::LeafCandidate::et(), reco::LeafCandidate::eta(), edm::EventSetup::get(), AnalyticalPropagator::propagate(), LargeD0_PixelPairStep_cff::propagator, reco::LeafCandidate::px(), reco::LeafCandidate::py(), reco::LeafCandidate::pz(), reco::LeafCandidate::setP4(), funct::sin(), reco::LeafCandidate::status(), reco::LeafCandidate::vx(), reco::LeafCandidate::vy(), and reco::LeafCandidate::vz().

                                                                                                                                              {
  
  
  vector<reco::GenParticle> curvedPars;

  edm::ESHandle<MagneticField> magField;
  iSetup.get<IdealMagneticFieldRecord>().get(magField);

  Cylinder::CylinderPointer theBarrel=Cylinder::build(Surface::PositionType(0,0,0),Surface::RotationType(),ECALBARRELRADIUS_);
  Plane::PlanePointer endCapPlus=Plane::build(Surface::PositionType(0,0,ECALENDCAPZ_),Surface::RotationType());
  Plane::PlanePointer endCapMinus=Plane::build(Surface::PositionType(0,0,-1*ECALENDCAPZ_),Surface::RotationType());

  AnalyticalPropagator propagator(&(*magField),alongMomentum);  

  for (uint32_t ig=0;ig<genPars.size();ig++) {
    reco::GenParticle gp=genPars.at(ig);
    //don't bend trajectories of neutral particles, unstable particles, particles with < 1 GeV
    //particles with < ~0.9 GeV don't reach the barrel
    //so just put them as-is into the new vector
    if (gp.charge()==0 || gp.status()!=1 || gp.et()<1) {
      curvedPars.push_back(gp);
      continue;
    }
    GlobalPoint vertex(gp.vx(),gp.vy(),gp.vz());
    GlobalVector gvect(gp.px(),gp.py(),gp.pz());
    FreeTrajectoryState fts(vertex,gvect,gp.charge(),&(*magField));
    TrajectoryStateOnSurface impact;
    //choose to propagate to barrel, +Z endcap, or -Z endcap, according to eta
    if (fabs(gp.eta())<ECALBARRELMAXETA_) {
      impact=propagator.propagate(fts,*theBarrel);
    } else if (gp.eta()>0) {
      impact=propagator.propagate(fts,*endCapPlus);
    } else {
      impact=propagator.propagate(fts,*endCapMinus);
    }
    //in case the particle doesn't reach the barrel/endcap, just put it as-is into the new vector
    //it should reach though.
    if (!impact.isValid()) {
      curvedPars.push_back(gp);
      continue;
    }
    math::XYZTLorentzVector newp4;

    //the magnitude of p doesn't change, only the direction
    //NB I do get some small change in magnitude by the following,
    //I think it is a numerical precision issue
    float et=gp.et();
    float phinew=impact.globalPosition().phi();
    float pxnew=et*cos(phinew);
    float pynew=et*sin(phinew);
    newp4.SetPx(pxnew);
    newp4.SetPy(pynew);
    newp4.SetPz(gp.pz());
    newp4.SetE(gp.energy());
    reco::GenParticle gpnew=gp;
    gpnew.setP4(newp4);
    curvedPars.push_back(gpnew);
  }
  return curvedPars;


}
float doubleEMEnrichingFilterAlgo::deltaRxyAtEE ( const reco::GenParticle gp1,
const reco::GenParticle gp2 
) [private]

Definition at line 313 of file doubleEMEnrichingFilterAlgo.cc.

References funct::cos(), reco::LeafCandidate::phi(), reco::LeafCandidate::pz(), funct::sin(), mathSSE::sqrt(), funct::tan(), and reco::LeafCandidate::theta().

                                                                                                      {
  
  if (gp1.pz()*gp2.pz() < 0) return 9999.;
  
  float rxy1=ECALENDCAPZ_*tan(gp1.theta());
  float x1=cos(gp1.phi())*rxy1;
  float y1=sin(gp1.phi())*rxy1;

  float rxy2=ECALENDCAPZ_*tan(gp2.theta());
  float x2=cos(gp2.phi())*rxy2;
  float y2=sin(gp2.phi())*rxy2;
  
  float dxy=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
  return dxy;


}
bool doubleEMEnrichingFilterAlgo::filter ( const edm::Event iEvent,
const edm::EventSetup iSetup 
)

BCtoE filter modified to store particle seeds

Definition at line 57 of file doubleEMEnrichingFilterAlgo.cc.

References abs, reco::LeafCandidate::et(), reco::LeafCandidate::eta(), edm::Event::getByLabel(), reco::LeafCandidate::pdgId(), query::result, and reco::LeafCandidate::status().

                                                                                             {

  
  //----------

  Handle<reco::GenParticleCollection> genParsHandle;
  iEvent.getByLabel(genParSource_,genParsHandle);
  reco::GenParticleCollection genPars=*genParsHandle;

  bool result = false;
  int BCtoEgood=0;
  int PhotElecgood=0;
  int IsoGenPargood=0;

  // cleaning
  sel1seeds.clear();
  sel2seeds.clear();
  selBCtoEseeds.clear();

  //bending of traj. of charged particles under influence of B-field
  std::vector<reco::GenParticle> genParsCurved=applyBFieldCurv(genPars,iSetup);

  PhotElecgood=filterPhotonElectronSeed(clusterThreshold_,
                                        seedThreshold_,
                                        isoConeSize_,
                                        hOverEMax_,
                                        tkIsoMax_,
                                        caloIsoMax_,
                                        requireTrackMatch_,
                                        genPars,
                                        genParsCurved);
  
  IsoGenPargood=filterIsoGenPar(isoGenParETMin_,isoGenParConeSize_,genPars,genParsCurved);


  for (uint32_t ig=0;ig<genParsCurved.size();ig++) {
    reco::GenParticle gp=genParsCurved.at(ig);
    if (gp.status()==1 && abs(gp.pdgId())==11 && gp.et()>eTThreshold_ && fabs(gp.eta())<FILTER_ETA_MAX_) {
      if (hasBCAncestors(gp)) {
        BCtoEgood++;
        selBCtoEseeds.push_back(gp);
      }
    }
  }

  // we want 2 different em candidates
  if ( PhotElecgood> 1) { 
    result = true;
  }
  else if ( IsoGenPargood> 1) {
    result = true;  
  }
  else if ( BCtoEgood > 1) {
    result = true;
  }
  else if ( PhotElecgood==1 && IsoGenPargood==1) {
    if ( (sel1seeds.at(0).eta()!=sel2seeds.at(0).eta()) && (sel1seeds.at(0).phi()!=sel2seeds.at(0).phi()) && (sel1seeds.at(0).et()!=sel2seeds.at(0).et()) ) result = true;
  }
  else if ( PhotElecgood==1 && BCtoEgood==1) {
    if ( (sel1seeds.at(0).eta()!=selBCtoEseeds.at(0).eta()) && (sel1seeds.at(0).phi()!=selBCtoEseeds.at(0).phi()) && (sel1seeds.at(0).et()!=selBCtoEseeds.at(0).et()) ) result = true;
  }
  else if ( BCtoEgood==1 && IsoGenPargood==1) {
    if ( (selBCtoEseeds.at(0).eta()!=sel2seeds.at(0).eta()) && (selBCtoEseeds.at(0).phi()!=sel2seeds.at(0).phi()) && (selBCtoEseeds.at(0).et()!=sel2seeds.at(0).et()) ) result = true;
  }
  
  return result;
        
}
int doubleEMEnrichingFilterAlgo::filterIsoGenPar ( float  etmin,
float  conesize,
const reco::GenParticleCollection gph,
const reco::GenParticleCollection gphCurved 
) [private]

Definition at line 338 of file doubleEMEnrichingFilterAlgo.cc.

References abs, reco::LeafCandidate::charge(), deltaR(), reco::LeafCandidate::energy(), reco::LeafCandidate::et(), reco::LeafCandidate::eta(), reco::LeafCandidate::pdgId(), createTree::pp, and reco::LeafCandidate::status().

                                             {

  int passed = 0;
  for (uint32_t ip=0;ip<gph.size();ip++) {

    reco::GenParticle gp=gph.at(ip);
    reco::GenParticle gpCurved=gphCurved.at(ip);    
    int gpabsid=abs(gp.pdgId());
    //find potential candidates
    if (gp.et()<=etmin || gp.status()!=1) continue;
    if (gpabsid!=11 && gpabsid != 211 && gpabsid!= 321) continue;
    if (fabs(gp.eta()) < FILTER_ETA_MIN_) continue;
    if (fabs(gp.eta()) > FILTER_ETA_MAX_) continue;
    
    //check isolation
    float tkiso=0;
    float caloiso=0;
    for (uint32_t jp=0;jp<gph.size();jp++) {
      if (jp==ip) continue;
      reco::GenParticle pp=gph.at(jp);
      reco::GenParticle ppCurved=gphCurved.at(jp);
      if (pp.status() != 1) continue;
      float dr=deltaR(gp,pp);
      float drCurved=deltaR(gpCurved,ppCurved);
      if (abs(pp.charge())==1 && pp.et()>2 && dr<conesize) {
        tkiso+=pp.et();
      }
      if (pp.et()>2 && abs(pp.pdgId())!=22 && drCurved<conesize) {
        caloiso+=pp.energy();
      }
    }
    if (tkiso<FILTER_TKISOCUT_ && caloiso<FILTER_CALOISOCUT_) { 
      sel2seeds.push_back(gpCurved); 
      passed++;
    }
  }
  return passed;
}
int doubleEMEnrichingFilterAlgo::filterPhotonElectronSeed ( float  clusterthreshold,
float  seedthreshold,
float  isoConeSize,
float  hOverEMax,
float  tkIsoMax,
float  caloIsoMax,
bool  requiretrackmatch,
const std::vector< reco::GenParticle > &  genPars,
const std::vector< reco::GenParticle > &  genParsCurved 
) [private]

Definition at line 135 of file doubleEMEnrichingFilterAlgo.cc.

References abs, reco::LeafCandidate::charge(), Geom::deltaPhi(), deltaR(), reco::LeafCandidate::et(), reco::LeafCandidate::eta(), reco::LeafCandidate::pdgId(), reco::LeafCandidate::phi(), and reco::LeafCandidate::status().

                                                                                                      {

  float conesizeendcap=15;

  int retval=0;
  
  vector<reco::GenParticle> seeds;
  //find electron and photon seeds - must have E>seedthreshold GeV
  for (uint32_t is=0;is<genParsCurved.size();is++) {
    reco::GenParticle gp=genParsCurved.at(is);
    if (gp.status()!=1 || fabs(gp.eta()) > FILTER_ETA_MAX_ || fabs(gp.eta()) < FILTER_ETA_MIN_) continue;
    int absid=abs(gp.pdgId());
    if (absid!=11 && absid!=22) continue;
    if (gp.et()>seedthreshold) seeds.push_back(gp);
  }
  
  bool matchtrack=false;
  
  //for every seed, try to cluster stable particles about it in cone
  for (uint32_t is=0;is<seeds.size();is++) {
    float eTInCone=0;//eT associated to the electron cluster
    float tkIsoET=0;//tracker isolation energy
    float caloIsoET=0;//calorimeter isolation energy
    float hadET=0;//isolation energy from heavy hadrons that goes in the same area as the "electron" - so contributes to H/E
    bool isBarrel=fabs(seeds.at(is).eta())<ECALBARRELMAXETA_;
    for (uint32_t ig=0;ig<genParsCurved.size();ig++) {
      reco::GenParticle gp=genParsCurved.at(ig);
      reco::GenParticle gpUnCurv=genPars.at(ig);//for tk isolation, p at vertex
      if (gp.status()!=1) continue;
      int gpabsid=abs(gp.pdgId());
      if (gp.et()<1) continue;//ignore very soft particles
      //BARREL
      if (isBarrel) {
        float dr=deltaR(seeds.at(is),gp);
        float dphi=deltaPhi(seeds.at(is).phi(),gp.phi());
        float deta=fabs(seeds.at(is).eta()-gp.eta());
        if (deta<0.03 && dphi<0.2) {
          if (gpabsid==22 || gpabsid==11 || gpabsid==211 || gpabsid==321) {
            //contributes to electron
            eTInCone+=gp.et();
            //check for a matched track with at least 5 GeV
            if ((gpabsid==11 || gpabsid==211 || gpabsid==321) && gp.et()>5) matchtrack=true;
          } else {
            //contributes to H/E
            hadET+=gp.et();
          }
        } else {
          float drUnCurv=deltaR(seeds.at(is),gpUnCurv);
          if ((gp.charge()==0 && dr<isoConeSize && gpabsid!=22) ||
              (gp.charge()!=0 && drUnCurv<isoConeSize)) {      
            //contributes to calo isolation energy
            caloIsoET+=gp.et();
          }
          if (gp.charge()!=0 && drUnCurv<isoConeSize) {
            //contributes to track isolation energy
            tkIsoET+=gp.et();
          }
        }
        //ENDCAP
      } else {
        float drxy=deltaRxyAtEE(seeds.at(is),gp);
        float dr=deltaR(seeds.at(is),gp);//the isolation is done in dR
        if (drxy<conesizeendcap) {
          if (gpabsid==22 || gpabsid==11 || gpabsid==211 || gpabsid==321) {
            //contributes to electron
            eTInCone+=gp.et();
            //check for a matched track with at least 5 GeV
            if ((gpabsid==11 || gpabsid==211 || gpabsid==321) && gp.et()>5) matchtrack=true;
          } else {
            //contributes to H/E
            hadET+=gp.et();
          }
        } else {
          float drUnCurv=deltaR(seeds.at(is),gpUnCurv);
          if ((gp.charge()==0 && dr<isoConeSize && gpabsid!=22) ||
              (gp.charge()!=0 && drUnCurv<isoConeSize)) {
            //contributes to calo isolation energy
            caloIsoET+=gp.et();
          }
          if (gp.charge()!=0 && drUnCurv<isoConeSize) {
            //contributes to track isolation energy
            tkIsoET+=gp.et();
          }
        }
      }
    }

    if (eTInCone>clusterthreshold && (!requiretrackmatch || matchtrack)) {
      //       cout <<"isoET: "<<isoET<<endl;
      if (hadET/eTInCone<hOverEMax && tkIsoET<tkIsoMax && caloIsoET<caloIsoMax) { 
        retval=retval+1;
        sel1seeds.push_back(seeds[is]);
        // break;
      }
    }
  }
  
  return retval;
}
bool doubleEMEnrichingFilterAlgo::hasBCAncestors ( reco::GenParticle  gp)

Definition at line 382 of file doubleEMEnrichingFilterAlgo.cc.

References reco::CompositeRefCandidateT< D >::motherRef(), and reco::CompositeRefCandidateT< D >::numberOfMothers().

                                                                   {

  //stopping condition: this particle is a b or c hadron
  if (isBCHadron(gp)) return true;
  //stopping condition: this particle has no mothers
  if (gp.numberOfMothers()==0) return false;
  //otherwise continue
  bool retval=false;
  for (uint32_t im=0;im<gp.numberOfMothers();im++) {
    retval=retval || hasBCAncestors(*gp.motherRef(im));
  }
  return retval;
  
}
bool doubleEMEnrichingFilterAlgo::isBCBaryon ( reco::GenParticle  gp) [private]

Definition at line 413 of file doubleEMEnrichingFilterAlgo.cc.

References abs, and reco::LeafCandidate::pdgId().

                                                               {
  
  uint32_t pdgid=abs(gp.pdgId());
  uint32_t thousands=pdgid%10000;
  if (thousands>=4000 && thousands <6000) {
    return true;
  } else {
    return false;
  }

}
bool doubleEMEnrichingFilterAlgo::isBCHadron ( reco::GenParticle  gp) [private]

Definition at line 397 of file doubleEMEnrichingFilterAlgo.cc.

                                                               {
  return isBCMeson(gp) || isBCBaryon(gp);
}
bool doubleEMEnrichingFilterAlgo::isBCMeson ( reco::GenParticle  gp) [private]

Definition at line 401 of file doubleEMEnrichingFilterAlgo.cc.

References abs, and reco::LeafCandidate::pdgId().

                                                              {
  
  uint32_t pdgid=abs(gp.pdgId());
  uint32_t hundreds=pdgid%1000;
  if (hundreds>=400 && hundreds<600) {
    return true;
  } else {
    return false;
  }

}

Member Data Documentation

Definition at line 77 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 72 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 65 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 66 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 67 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 78 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 62 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 64 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 63 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 61 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 80 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 75 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 74 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 71 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 70 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 79 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 73 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 83 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 84 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 85 of file doubleEMEnrichingFilterAlgo.h.

Definition at line 76 of file doubleEMEnrichingFilterAlgo.h.