MuonVPlusJetsIDSelectionFunctor.h

Go to the documentation of this file.

#ifndef PhysicsTools_PatUtils_interface_MuonVPlusJetsIDSelectionFunctor_h
#define PhysicsTools_PatUtils_interface_MuonVPlusJetsIDSelectionFunctor_h

#ifndef __GCCXML__
#include "FWCore/Framework/interface/ConsumesCollector.h"
#endif
#include "DataFormats/PatCandidates/interface/Muon.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "PhysicsTools/SelectorUtils/interface/Selector.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include <iostream>

class MuonVPlusJetsIDSelectionFunctor : public Selector<pat::Muon> {

 public: // interface

  bool verbose_;

  enum Version_t { SUMMER08, FIRSTDATA, SPRING10, FALL10, N_VERSIONS, KITQCD };

  MuonVPlusJetsIDSelectionFunctor() {}

#ifndef __GCCXML__
  MuonVPlusJetsIDSelectionFunctor( edm::ParameterSet const & parameters, edm::ConsumesCollector& iC ) :
    MuonVPlusJetsIDSelectionFunctor(parameters)
  {
    beamLineSrcToken_ = iC.consumes<reco::BeamSpot>(beamLineSrc_);
    pvSrcToken_ = iC.consumes<std::vector<reco::Vertex> >(pvSrc_);
  }
#endif

  MuonVPlusJetsIDSelectionFunctor( edm::ParameterSet const & parameters ) {

    verbose_ = false;

    std::string versionStr = parameters.getParameter<std::string>("version");

    Version_t version = N_VERSIONS;

    if ( versionStr == "SUMMER08" ) {
      version = SUMMER08;
    }
    else if ( versionStr == "FIRSTDATA" ) {
      version = FIRSTDATA;
    }
    else if ( versionStr == "SPRING10" ) {
      version = SPRING10;
    }
    else if ( versionStr == "FALL10" ) {
      version = FALL10;
      if (verbose_) std::cout << "\nMUON SELECTION - you are using FALL10 Selection" << std::endl;
    }
    else if (versionStr == "KITQCD") {
      version = KITQCD;
      if (verbose_) std::cout << "\nMUON SELECTION - you are using KITQCD Selection" << std::endl;
    }
    else {
      throw cms::Exception("InvalidInput") << "Expect version to be one of SUMMER08, FIRSTDATA, SPRING10, FALL10" << std::endl;
    }

    initialize( version,
        parameters.getParameter<double>("Chi2"),
        parameters.getParameter<double>("D0")  ,
        parameters.getParameter<double>("ED0")  ,
        parameters.getParameter<double>("SD0")  ,
        parameters.getParameter<int>   ("NHits")   ,
        parameters.getParameter<int>   ("NValMuHits"),
        parameters.getParameter<double>("ECalVeto")   ,
        parameters.getParameter<double>("HCalVeto")   ,
        parameters.getParameter<double>("RelIso"),
        parameters.getParameter<double>("LepZ"),
        parameters.getParameter<int>("nPixelHits"),
        parameters.getParameter<int>("nMatchedStations")
        );
    if ( parameters.exists("cutsToIgnore") )
      setIgnoredCuts( parameters.getParameter<std::vector<std::string> >("cutsToIgnore") );

    retInternal_ = getBitTemplate();

    recalcDBFromBSp_ = parameters.getParameter<bool>("RecalcFromBeamSpot");
    pvSrc_ = parameters.getParameter<edm::InputTag>("pvSrc");
  }


  MuonVPlusJetsIDSelectionFunctor( Version_t version,
                   double chi2 = 10.0,
                   double d0 = 0.2,
                   double ed0 = 999.0,
                   double sd0 = 999.0,
                   int nhits = 11,
                   int nValidMuonHits = 0,
                   double ecalveto = 4.0,
                   double hcalveto = 6.0,
                   double reliso = 0.05,
                   double maxLepZ = 1.0,
                   int minPixelHits = 1,
                   int minNMatches = 1
                   ) : recalcDBFromBSp_(false) {
    initialize( version, chi2, d0, ed0, sd0, nhits, nValidMuonHits, ecalveto, hcalveto, reliso,
        maxLepZ, minPixelHits, minNMatches );

    retInternal_ = getBitTemplate();
  }




  void initialize( Version_t version,
           double chi2 = 10.0,
           double d0 = 999.0,
           double ed0 = 999.0,
           double sd0 = 3.0,
           int nhits = 11,
           int nValidMuonHits = 0,
           double ecalveto = 4.0,
           double hcalveto = 6.0,
           double reliso = 0.05,
           double maxLepZ = 1.0,
           int minPixelHits = 1,
           int minNMatches = 1 )
  {
    version_ = version;

    push_back("Chi2",      chi2   );
    push_back("D0",        d0     );
    push_back("ED0",       ed0    );
    push_back("SD0",       sd0    );
    push_back("NHits",     nhits  );
    push_back("NValMuHits",     nValidMuonHits  );
    push_back("ECalVeto",  ecalveto);
    push_back("HCalVeto",  hcalveto);
    push_back("RelIso",    reliso );
    push_back("LepZ",      maxLepZ);
    push_back("nPixelHits",minPixelHits);
    push_back("nMatchedStations", minNMatches);

    set("Chi2");
    set("D0");
    set("ED0");
    set("SD0");
    set("NHits");
    set("NValMuHits");
    set("ECalVeto");
    set("HCalVeto");
    set("RelIso");
    set("LepZ");
    set("nPixelHits");
    set("nMatchedStations");

    indexChi2_          = index_type(&bits_, "Chi2"         );
    indexD0_            = index_type(&bits_, "D0"           );
    indexED0_           = index_type(&bits_, "ED0"          );
    indexSD0_           = index_type(&bits_, "SD0"          );
    indexNHits_         = index_type(&bits_, "NHits"        );
    indexNValMuHits_    = index_type(&bits_, "NValMuHits"   );
    indexECalVeto_      = index_type(&bits_, "ECalVeto"     );
    indexHCalVeto_      = index_type(&bits_, "HCalVeto"     );
    indexRelIso_        = index_type(&bits_, "RelIso"       );
    indexLepZ_          = index_type( &bits_, "LepZ");
    indexPixHits_       = index_type( &bits_, "nPixelHits");
    indexStations_      = index_type( &bits_, "nMatchedStations");

    if ( version == FALL10) {
      set("ED0", false );
      set("SD0", false);
      set("ECalVeto",false);
      set("HCalVeto",false);
    } else if ( version == SPRING10) {
      set("ED0", false );
      set("SD0", false);
      set("ECalVeto",false);
      set("HCalVeto",false);
      set("LepZ", false);
      set("nPixelHits", false);
      set("nMatchedStations", false);
    } else if ( version_ == FIRSTDATA ) {
      set("D0", false );
      set("ED0", false );
      set("NValMuHits",false);
      set("LepZ", false);
      set("nPixelHits", false);
      set("nMatchedStations", false);
    } else if (version == SUMMER08 ) {
      set("SD0", false);
      set("NValMuHits",false);
      set("LepZ", false);
      set("nPixelHits", false);
      set("nMatchedStations", false);

    }

  }

  // Allow for multiple definitions of the cuts.
  bool operator()( const pat::Muon & muon, edm::EventBase const & event, pat::strbitset & ret ) override
  {

    if (version_ == FALL10 ) return fall10Cuts(muon, event, ret);
    else if (version_ == SPRING10 ) return spring10Cuts(muon, event, ret);
    else if ( version_ == SUMMER08 ) return summer08Cuts( muon, ret );
    else if ( version_ == FIRSTDATA ) return firstDataCuts( muon, ret );
    else if ( version_ == KITQCD ) {
      if (verbose_) std::cout << "Calling KIT selection method" << std::endl;
      return kitQCDCuts (muon, event, ret);
    }
    else {
      return false;
    }
  }

  // For compatibility with older versions.
  bool operator()( const pat::Muon & muon, pat::strbitset & ret ) override
  {

    if (version_ == SPRING10 || version_ == FALL10 ) throw cms::Exception("LogicError")
      << "MuonVPlusJetsSelectionFunctor SPRING10 and FALL10 versions needs the event! Call operator()(muon,event,ret)"
      <<std::endl;

    else if ( version_ == SUMMER08 ) return summer08Cuts( muon, ret );
    else if ( version_ == FIRSTDATA ) return firstDataCuts( muon, ret );
    else {
      return false;
    }
  }


  using Selector<pat::Muon>::operator();

  // cuts based on craft 08 analysis.
  bool summer08Cuts( const pat::Muon & muon, pat::strbitset & ret)
  {

    ret.set(false);

    double norm_chi2 = muon.normChi2();
    double corr_d0 = muon.dB();
    int nhits = static_cast<int>( muon.numberOfValidHits() );

    double ecalVeto = muon.isolationR03().emVetoEt;
    double hcalVeto = muon.isolationR03().hadVetoEt;

    double hcalIso = muon.hcalIso();
    double ecalIso = muon.ecalIso();
    double trkIso  = muon.trackIso();
    double pt      = muon.pt() ;

    double relIso = (ecalIso + hcalIso + trkIso) / pt;

    if ( norm_chi2     <  cut(indexChi2_,   double()) || ignoreCut(indexChi2_)    ) passCut(ret, indexChi2_   );
    if ( fabs(corr_d0) <  cut(indexD0_,     double()) || ignoreCut(indexD0_)      ) passCut(ret, indexD0_     );
    if ( nhits         >= cut(indexNHits_,  int()   ) || ignoreCut(indexNHits_)   ) passCut(ret, indexNHits_  );
    if ( hcalVeto      <  cut(indexHCalVeto_,double())|| ignoreCut(indexHCalVeto_)) passCut(ret, indexHCalVeto_);
    if ( ecalVeto      <  cut(indexECalVeto_,double())|| ignoreCut(indexECalVeto_)) passCut(ret, indexECalVeto_);
    if ( relIso        <  cut(indexRelIso_, double()) || ignoreCut(indexRelIso_)  ) passCut(ret, indexRelIso_ );

    setIgnored(ret);

    return (bool)ret;
  }



  // cuts based on craft 08 analysis.
  bool firstDataCuts( const pat::Muon & muon, pat::strbitset & ret)
  {

    ret.set(false);

    double norm_chi2 = muon.normChi2();
    double corr_d0 = muon.dB();
    double corr_ed0 = muon.edB();
    double corr_sd0 = ( corr_ed0 > 0.000000001 ) ? corr_d0 / corr_ed0 : 999.0;
    int nhits = static_cast<int>( muon.numberOfValidHits() );

    double ecalVeto = muon.isolationR03().emVetoEt;
    double hcalVeto = muon.isolationR03().hadVetoEt;

    double hcalIso = muon.hcalIso();
    double ecalIso = muon.ecalIso();
    double trkIso  = muon.trackIso();
    double pt      = muon.pt() ;

    double relIso = (ecalIso + hcalIso + trkIso) / pt;

    if ( norm_chi2     <  cut(indexChi2_,   double()) || ignoreCut(indexChi2_)    ) passCut(ret, indexChi2_   );
    if ( fabs(corr_d0) <  cut(indexD0_,     double()) || ignoreCut(indexD0_)      ) passCut(ret, indexD0_     );
    if ( fabs(corr_ed0)<  cut(indexED0_,    double()) || ignoreCut(indexED0_)     ) passCut(ret, indexED0_    );
    if ( fabs(corr_sd0)<  cut(indexSD0_,    double()) || ignoreCut(indexSD0_)     ) passCut(ret, indexSD0_    );
    if ( nhits         >= cut(indexNHits_,  int()   ) || ignoreCut(indexNHits_)   ) passCut(ret, indexNHits_  );
    if ( hcalVeto      <  cut(indexHCalVeto_,double())|| ignoreCut(indexHCalVeto_)) passCut(ret, indexHCalVeto_);
    if ( ecalVeto      <  cut(indexECalVeto_,double())|| ignoreCut(indexECalVeto_)) passCut(ret, indexECalVeto_);
    if ( relIso        <  cut(indexRelIso_, double()) || ignoreCut(indexRelIso_)  ) passCut(ret, indexRelIso_ );

    setIgnored(ret);

    return (bool)ret;
  }

  // cuts based on top group L+J synchronization exercise
  bool spring10Cuts( const pat::Muon & muon, edm::EventBase const & event, pat::strbitset & ret)
  {

    ret.set(false);

    double norm_chi2 = muon.normChi2();
    double corr_d0 = muon.dB();
    double corr_ed0 = muon.edB();
    double corr_sd0 = ( corr_ed0 > 0.000000001 ) ? corr_d0 / corr_ed0 : 999.0;

    //If required, recalculate the impact parameter using the beam spot
    if (recalcDBFromBSp_) {

      //Get the beam spot
      reco::TrackBase::Point beamPoint(0,0,0);
      reco::BeamSpot beamSpot;
      edm::Handle<reco::BeamSpot> beamSpotHandle;
      event.getByLabel(beamLineSrc_, beamSpotHandle);

      if( beamSpotHandle.isValid() ){
    beamSpot = *beamSpotHandle;
      } else{
    edm::LogError("DataNotAvailable")
      << "No beam spot available from EventSetup, not adding high level selection \n";
      }
      beamPoint = reco::TrackBase::Point ( beamSpot.x0(), beamSpot.y0(), beamSpot.z0() );

      //Use the beamspot to correct the impact parameter and uncertainty
      reco::TrackRef innerTrack = muon.innerTrack();
      if ( innerTrack.isNonnull() && innerTrack.isAvailable() ) {
    corr_d0 = -1.0 * innerTrack->dxy( beamPoint );
    corr_ed0 = sqrt( innerTrack->d0Error() * innerTrack->d0Error()
             + 0.5* beamSpot.BeamWidthX()*beamSpot.BeamWidthX()
             + 0.5* beamSpot.BeamWidthY()*beamSpot.BeamWidthY() );
    corr_sd0 = ( corr_ed0 > 0.000000001 ) ? corr_d0 / corr_ed0 : 999.0;

      } else {
    corr_d0 =  999.;
    corr_ed0 = 999.;
      }
    }

    int nhits = static_cast<int>( muon.numberOfValidHits() );
    int nValidMuonHits = static_cast<int> (muon.globalTrack()->hitPattern().numberOfValidMuonHits());

    double ecalVeto = muon.isolationR03().emVetoEt;
    double hcalVeto = muon.isolationR03().hadVetoEt;

    double hcalIso = muon.hcalIso();
    double ecalIso = muon.ecalIso();
    double trkIso  = muon.trackIso();
    double pt      = muon.pt() ;

    double relIso = (ecalIso + hcalIso + trkIso) / pt;

    if ( norm_chi2     <  cut(indexChi2_,   double()) || ignoreCut(indexChi2_)    ) passCut(ret, indexChi2_   );
    if ( fabs(corr_d0) <  cut(indexD0_,     double()) || ignoreCut(indexD0_)      ) passCut(ret, indexD0_     );
    if ( fabs(corr_ed0)<  cut(indexED0_,    double()) || ignoreCut(indexED0_)     ) passCut(ret, indexED0_    );
    if ( fabs(corr_sd0)<  cut(indexSD0_,    double()) || ignoreCut(indexSD0_)     ) passCut(ret, indexSD0_    );
    if ( nhits         >= cut(indexNHits_,  int()   ) || ignoreCut(indexNHits_)   ) passCut(ret, indexNHits_  );
    if ( nValidMuonHits> cut(indexNValMuHits_,int()) || ignoreCut(indexNValMuHits_)) passCut(ret, indexNValMuHits_  );
    if ( hcalVeto      <  cut(indexHCalVeto_,double())|| ignoreCut(indexHCalVeto_)) passCut(ret, indexHCalVeto_);
    if ( ecalVeto      <  cut(indexECalVeto_,double())|| ignoreCut(indexECalVeto_)) passCut(ret, indexECalVeto_);
    if ( relIso        <  cut(indexRelIso_, double()) || ignoreCut(indexRelIso_)  ) passCut(ret, indexRelIso_ );

    setIgnored(ret);

    return (bool)ret;
  }




  // cuts based on top group L+J synchronization exercise
  bool fall10Cuts( const pat::Muon & muon, edm::EventBase const & event, pat::strbitset & ret)
  {

    ret.set(false);

    double norm_chi2 = muon.normChi2();
    double corr_d0 = muon.dB();
    double corr_ed0 = muon.edB();
    double corr_sd0 = ( corr_ed0 > 0.000000001 ) ? corr_d0 / corr_ed0 : 999.0;

    // Get the PV for the muon z requirement
    edm::Handle<std::vector<reco::Vertex> > pvtxHandle_;
    event.getByLabel( pvSrc_, pvtxHandle_ );

    double zvtx = -999;
    if ( !pvtxHandle_->empty() ) {
      zvtx = pvtxHandle_->at(0).z();
    } else {
      throw cms::Exception("InvalidInput") << " There needs to be at least one primary vertex in the event." << std::endl;
    }

    //If required, recalculate the impact parameter using the beam spot
    if (recalcDBFromBSp_) {

      //Get the beam spot
      reco::TrackBase::Point beamPoint(0,0,0);
      reco::BeamSpot beamSpot;
      edm::Handle<reco::BeamSpot> beamSpotHandle;
      event.getByLabel(beamLineSrc_, beamSpotHandle);

      if( beamSpotHandle.isValid() ){
    beamSpot = *beamSpotHandle;
      } else{
    edm::LogError("DataNotAvailable")
      << "No beam spot available from EventSetup, not adding high level selection \n";
      }
      beamPoint = reco::TrackBase::Point ( beamSpot.x0(), beamSpot.y0(), beamSpot.z0() );

      //Use the beamspot to correct the impact parameter and uncertainty
      reco::TrackRef innerTrack = muon.innerTrack();
      if ( innerTrack.isNonnull() && innerTrack.isAvailable() ) {
    corr_d0 = -1.0 * innerTrack->dxy( beamPoint );
    corr_ed0 = sqrt( innerTrack->d0Error() * innerTrack->d0Error()
             + 0.5* beamSpot.BeamWidthX()*beamSpot.BeamWidthX()
             + 0.5* beamSpot.BeamWidthY()*beamSpot.BeamWidthY() );
    corr_sd0 = ( corr_ed0 > 0.000000001 ) ? corr_d0 / corr_ed0 : 999.0;

      } else {
    corr_d0 =  999.;
    corr_ed0 = 999.;
      }
    }

    int nhits = static_cast<int>( muon.numberOfValidHits() );
    int nValidMuonHits = static_cast<int> (muon.globalTrack()->hitPattern().numberOfValidMuonHits());

    double ecalVeto = muon.isolationR03().emVetoEt;
    double hcalVeto = muon.isolationR03().hadVetoEt;

    double hcalIso = muon.hcalIso();
    double ecalIso = muon.ecalIso();
    double trkIso  = muon.trackIso();
    double pt      = muon.pt() ;

    double relIso = (ecalIso + hcalIso + trkIso) / pt;


    double z_mu = muon.vertex().z();

    int nPixelHits = muon.innerTrack()->hitPattern().pixelLayersWithMeasurement();

    int nMatchedStations = muon.numberOfMatches();

    if ( norm_chi2     <  cut(indexChi2_,   double()) || ignoreCut(indexChi2_)    ) passCut(ret, indexChi2_   );
    if ( fabs(corr_d0) <  cut(indexD0_,     double()) || ignoreCut(indexD0_)      ) passCut(ret, indexD0_     );
    if ( fabs(corr_ed0)<  cut(indexED0_,    double()) || ignoreCut(indexED0_)     ) passCut(ret, indexED0_    );
    if ( fabs(corr_sd0)<  cut(indexSD0_,    double()) || ignoreCut(indexSD0_)     ) passCut(ret, indexSD0_    );
    if ( nhits         >= cut(indexNHits_,  int()   ) || ignoreCut(indexNHits_)   ) passCut(ret, indexNHits_  );
    if ( nValidMuonHits> cut(indexNValMuHits_,int()) || ignoreCut(indexNValMuHits_)) passCut(ret, indexNValMuHits_  );
    if ( hcalVeto      <  cut(indexHCalVeto_,double())|| ignoreCut(indexHCalVeto_)) passCut(ret, indexHCalVeto_);
    if ( ecalVeto      <  cut(indexECalVeto_,double())|| ignoreCut(indexECalVeto_)) passCut(ret, indexECalVeto_);
    if ( relIso        <  cut(indexRelIso_, double()) || ignoreCut(indexRelIso_)  ) passCut(ret, indexRelIso_ );
    if ( fabs(z_mu-zvtx)<  cut(indexLepZ_, double()) || ignoreCut(indexLepZ_)  ) passCut(ret, indexLepZ_ );
    if ( nPixelHits    >  cut(indexPixHits_,int())    || ignoreCut(indexPixHits_))  passCut(ret, indexPixHits_);
    if ( nMatchedStations> cut(indexStations_,int())    || ignoreCut(indexStations_))  passCut(ret, indexStations_);

    setIgnored(ret);

    return (bool)ret;
  }


  // cuts based on top group L+J synchronization exercise
  // this is a copy of fall 10 cuts
  // with a hack to include a double-sided reliso cut

  bool kitQCDCuts( const pat::Muon & muon, edm::EventBase const & event, pat::strbitset & ret)
  {

    ret.set(false);

    double norm_chi2 = muon.normChi2();
    double corr_d0 = muon.dB();
    double corr_ed0 = muon.edB();
    double corr_sd0 = ( corr_ed0 > 0.000000001 ) ? corr_d0 / corr_ed0 : 999.0;

    // Get the PV for the muon z requirement
    edm::Handle<std::vector<reco::Vertex> > pvtxHandle_;
    event.getByLabel( pvSrc_, pvtxHandle_ );

    double zvtx = -999;
    if ( !pvtxHandle_->empty() ) {
      zvtx = pvtxHandle_->at(0).z();
    } else {
      throw cms::Exception("InvalidInput") << " There needs to be at least one primary vertex in the event." << std::endl;
    }

    //If required, recalculate the impact parameter using the beam spot
    if (recalcDBFromBSp_) {

      //Get the beam spot
      reco::TrackBase::Point beamPoint(0,0,0);
      reco::BeamSpot beamSpot;
      edm::Handle<reco::BeamSpot> beamSpotHandle;
      event.getByLabel(beamLineSrc_, beamSpotHandle);

      if( beamSpotHandle.isValid() ){
    beamSpot = *beamSpotHandle;
      } else{
    edm::LogError("DataNotAvailable")
      << "No beam spot available from EventSetup, not adding high level selection \n";
      }
      beamPoint = reco::TrackBase::Point ( beamSpot.x0(), beamSpot.y0(), beamSpot.z0() );

      //Use the beamspot to correct the impact parameter and uncertainty
      reco::TrackRef innerTrack = muon.innerTrack();
      if ( innerTrack.isNonnull() && innerTrack.isAvailable() ) {
    corr_d0 = -1.0 * innerTrack->dxy( beamPoint );
    corr_ed0 = sqrt( innerTrack->d0Error() * innerTrack->d0Error()
             + 0.5* beamSpot.BeamWidthX()*beamSpot.BeamWidthX()
             + 0.5* beamSpot.BeamWidthY()*beamSpot.BeamWidthY() );
    corr_sd0 = ( corr_ed0 > 0.000000001 ) ? corr_d0 / corr_ed0 : 999.0;

      } else {
    corr_d0 =  999.;
    corr_ed0 = 999.;
      }
    }

    int nhits = static_cast<int>( muon.numberOfValidHits() );
    int nValidMuonHits = static_cast<int> (muon.globalTrack()->hitPattern().numberOfValidMuonHits());

    double ecalVeto = muon.isolationR03().emVetoEt;
    double hcalVeto = muon.isolationR03().hadVetoEt;

    double hcalIso = muon.hcalIso();
    double ecalIso = muon.ecalIso();
    double trkIso  = muon.trackIso();
    double pt      = muon.pt() ;

    double relIso = (ecalIso + hcalIso + trkIso) / pt;


    double z_mu = muon.vertex().z();

    int nPixelHits = muon.innerTrack()->hitPattern().pixelLayersWithMeasurement();

    int nMatchedStations = muon.numberOfMatches();

    if ( norm_chi2     <  cut(indexChi2_,   double()) || ignoreCut(indexChi2_)    ) passCut(ret, indexChi2_   );
    if ( fabs(corr_d0) <  cut(indexD0_,     double()) || ignoreCut(indexD0_)      ) passCut(ret, indexD0_     );
    if ( fabs(corr_ed0)<  cut(indexED0_,    double()) || ignoreCut(indexED0_)     ) passCut(ret, indexED0_    );
    if ( fabs(corr_sd0)<  cut(indexSD0_,    double()) || ignoreCut(indexSD0_)     ) passCut(ret, indexSD0_    );
    if ( nhits         >= cut(indexNHits_,  int()   ) || ignoreCut(indexNHits_)   ) passCut(ret, indexNHits_  );
    if ( nValidMuonHits> cut(indexNValMuHits_,int()) || ignoreCut(indexNValMuHits_)) passCut(ret, indexNValMuHits_  );
    if ( hcalVeto      <  cut(indexHCalVeto_,double())|| ignoreCut(indexHCalVeto_)) passCut(ret, indexHCalVeto_);
    if ( ecalVeto      <  cut(indexECalVeto_,double())|| ignoreCut(indexECalVeto_)) passCut(ret, indexECalVeto_);
    if ( fabs(z_mu-zvtx)<  cut(indexLepZ_, double()) || ignoreCut(indexLepZ_)  ) passCut(ret, indexLepZ_ );
    if ( nPixelHits    >  cut(indexPixHits_,int())    || ignoreCut(indexPixHits_))  passCut(ret, indexPixHits_);
    if ( nMatchedStations> cut(indexStations_,int())    || ignoreCut(indexStations_))  passCut(ret, indexStations_);


    //
    // JMS Dec 13 2010
    // HACK
    // Need double-sided relIso cut to implement data-driven QCD
    //
    //
    if ( ((relIso > 0.2) && (relIso < 0.75))
         || ignoreCut(indexRelIso_)  ) passCut(ret, indexRelIso_ );




    setIgnored(ret);

    return (bool)ret;
  }




 private: // member variables

  Version_t version_;
  bool recalcDBFromBSp_;
  edm::InputTag beamLineSrc_;
#ifndef __GCCXML__
  edm::EDGetTokenT<reco::BeamSpot> beamLineSrcToken_;
#endif
  edm::InputTag pvSrc_;
#ifndef __GCCXML__
  edm::EDGetTokenT<std::vector<reco::Vertex> > pvSrcToken_;
#endif

  index_type indexChi2_;
  index_type indexD0_;
  index_type indexED0_;
  index_type indexSD0_;
  index_type indexNHits_;
  index_type indexNValMuHits_;
  index_type indexECalVeto_;
  index_type indexHCalVeto_;
  index_type indexRelIso_;
  index_type indexLepZ_;
  index_type indexPixHits_;
  index_type indexStations_;


};

#endif