CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/PhysicsTools/SelectorUtils/interface/MuonVPlusJetsIDSelectionFunctor.h

Go to the documentation of this file.
00001 #ifndef PhysicsTools_PatUtils_interface_MuonVPlusJetsIDSelectionFunctor_h
00002 #define PhysicsTools_PatUtils_interface_MuonVPlusJetsIDSelectionFunctor_h
00003 
00004 #include "DataFormats/PatCandidates/interface/Muon.h"
00005 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007 
00008 #include "PhysicsTools/SelectorUtils/interface/Selector.h"
00009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00010 
00011 #include <iostream>
00012 
00013 class MuonVPlusJetsIDSelectionFunctor : public Selector<pat::Muon> {
00014 
00015  public: // interface
00016 
00017   bool verbose_;
00018   
00019   enum Version_t { SUMMER08, FIRSTDATA, SPRING10, FALL10, N_VERSIONS, KITQCD };
00020 
00021   MuonVPlusJetsIDSelectionFunctor() {}
00022 
00023   MuonVPlusJetsIDSelectionFunctor( edm::ParameterSet const & parameters ) {
00024 
00025     verbose_ = false;
00026     
00027     std::string versionStr = parameters.getParameter<std::string>("version");
00028 
00029     Version_t version = N_VERSIONS;
00030 
00031     if ( versionStr == "SUMMER08" ) {
00032       version = SUMMER08;
00033     }
00034     else if ( versionStr == "FIRSTDATA" ) {
00035       version = FIRSTDATA;
00036     }
00037     else if ( versionStr == "SPRING10" ) {
00038       version = SPRING10;
00039     }
00040     else if ( versionStr == "FALL10" ) {
00041       version = FALL10;
00042       if (verbose_) std::cout << "\nMUON SELECTION - you are using FALL10 Selection" << std::endl;
00043     }
00044     else if (versionStr == "KITQCD") {
00045       version = KITQCD;
00046       if (verbose_) std::cout << "\nMUON SELECTION - you are using KITQCD Selection" << std::endl;
00047     }
00048     else {
00049       throw cms::Exception("InvalidInput") << "Expect version to be one of SUMMER08, FIRSTDATA, SPRING10, FALL10" << std::endl;
00050     }
00051 
00052     initialize( version, 
00053                 parameters.getParameter<double>("Chi2"),
00054                 parameters.getParameter<double>("D0")  ,
00055                 parameters.getParameter<double>("ED0")  ,
00056                 parameters.getParameter<double>("SD0")  ,
00057                 parameters.getParameter<int>   ("NHits")   ,
00058                 parameters.getParameter<int>   ("NValMuHits"),
00059                 parameters.getParameter<double>("ECalVeto")   ,
00060                 parameters.getParameter<double>("HCalVeto")   ,
00061                 parameters.getParameter<double>("RelIso"),
00062                 parameters.getParameter<double>("LepZ"), 
00063                 parameters.getParameter<int>("nPixelHits"),
00064                 parameters.getParameter<int>("nMatchedStations")
00065                 );
00066     if ( parameters.exists("cutsToIgnore") )
00067       setIgnoredCuts( parameters.getParameter<std::vector<std::string> >("cutsToIgnore") );
00068         
00069     retInternal_ = getBitTemplate();
00070 
00071     recalcDBFromBSp_ = parameters.getParameter<bool>("RecalcFromBeamSpot");
00072     beamLineSrc_ = parameters.getParameter<edm::InputTag>("beamLineSrc");
00073     pvSrc_ = parameters.getParameter<edm::InputTag>("pvSrc");
00074   }
00075 
00076 
00077   MuonVPlusJetsIDSelectionFunctor( Version_t version,
00078                                    double chi2 = 10.0,
00079                                    double d0 = 0.2,
00080                                    double ed0 = 999.0,
00081                                    double sd0 = 999.0,
00082                                    int nhits = 11,
00083                                    int nValidMuonHits = 0,
00084                                    double ecalveto = 4.0,
00085                                    double hcalveto = 6.0,
00086                                    double reliso = 0.05,
00087                                    double maxLepZ = 1.0,
00088                                    int minPixelHits = 1,
00089                                    int minNMatches = 1
00090                                    ) : recalcDBFromBSp_(false) {
00091     initialize( version, chi2, d0, ed0, sd0, nhits, nValidMuonHits, ecalveto, hcalveto, reliso,
00092                 maxLepZ, minPixelHits, minNMatches );
00093 
00094     retInternal_ = getBitTemplate();
00095   }
00096   
00097   
00098 
00099 
00100   void initialize( Version_t version,
00101                    double chi2 = 10.0,
00102                    double d0 = 999.0,
00103                    double ed0 = 999.0,
00104                    double sd0 = 3.0,
00105                    int nhits = 11,
00106                    int nValidMuonHits = 0,
00107                    double ecalveto = 4.0,
00108                    double hcalveto = 6.0,
00109                    double reliso = 0.05,
00110                    double maxLepZ = 1.0,
00111                    int minPixelHits = 1,
00112                    int minNMatches = 1 )
00113   {
00114     version_ = version; 
00115 
00116     push_back("Chi2",      chi2   );
00117     push_back("D0",        d0     );
00118     push_back("ED0",       ed0    );
00119     push_back("SD0",       sd0    );
00120     push_back("NHits",     nhits  );
00121     push_back("NValMuHits",     nValidMuonHits  );
00122     push_back("ECalVeto",  ecalveto);
00123     push_back("HCalVeto",  hcalveto);
00124     push_back("RelIso",    reliso );
00125     push_back("LepZ",      maxLepZ);
00126     push_back("nPixelHits",minPixelHits);
00127     push_back("nMatchedStations", minNMatches);
00128 
00129     set("Chi2");
00130     set("D0");
00131     set("ED0");
00132     set("SD0");
00133     set("NHits");
00134     set("NValMuHits");
00135     set("ECalVeto");
00136     set("HCalVeto");
00137     set("RelIso");   
00138     set("LepZ");
00139     set("nPixelHits");
00140     set("nMatchedStations");  
00141 
00142     indexChi2_          = index_type(&bits_, "Chi2"         );
00143     indexD0_            = index_type(&bits_, "D0"           );
00144     indexED0_           = index_type(&bits_, "ED0"          );
00145     indexSD0_           = index_type(&bits_, "SD0"          );
00146     indexNHits_         = index_type(&bits_, "NHits"        );
00147     indexNValMuHits_    = index_type(&bits_, "NValMuHits"   );
00148     indexECalVeto_      = index_type(&bits_, "ECalVeto"     );
00149     indexHCalVeto_      = index_type(&bits_, "HCalVeto"     );
00150     indexRelIso_        = index_type(&bits_, "RelIso"       );
00151     indexLepZ_          = index_type( &bits_, "LepZ");
00152     indexPixHits_       = index_type( &bits_, "nPixelHits");
00153     indexStations_      = index_type( &bits_, "nMatchedStations");
00154 
00155     if ( version == FALL10) {
00156       set("ED0", false );
00157       set("SD0", false);
00158       set("ECalVeto",false);
00159       set("HCalVeto",false);
00160     } else if ( version == SPRING10) {
00161       set("ED0", false );
00162       set("SD0", false);
00163       set("ECalVeto",false);
00164       set("HCalVeto",false);
00165       set("LepZ", false);
00166       set("nPixelHits", false);
00167       set("nMatchedStations", false);  
00168     } else if ( version_ == FIRSTDATA ) {
00169       set("D0", false );
00170       set("ED0", false );
00171       set("NValMuHits",false);
00172       set("LepZ", false);
00173       set("nPixelHits", false);
00174       set("nMatchedStations", false);  
00175     } else if (version == SUMMER08 ) {
00176       set("SD0", false);
00177       set("NValMuHits",false);      
00178       set("LepZ", false);
00179       set("nPixelHits", false);
00180       set("nMatchedStations", false);  
00181 
00182     }
00183 
00184   }
00185 
00186   // Allow for multiple definitions of the cuts. 
00187   bool operator()( const pat::Muon & muon, edm::EventBase const & event, pat::strbitset & ret ) 
00188   { 
00189 
00190     if (version_ == FALL10 ) return fall10Cuts(muon, event, ret);
00191     else if (version_ == SPRING10 ) return spring10Cuts(muon, event, ret);
00192     else if ( version_ == SUMMER08 ) return summer08Cuts( muon, ret );
00193     else if ( version_ == FIRSTDATA ) return firstDataCuts( muon, ret );
00194     else if ( version_ == KITQCD ) {
00195       if (verbose_) std::cout << "Calling KIT selection method" << std::endl;
00196       return kitQCDCuts (muon, event, ret);
00197     }      
00198     else {
00199       return false;
00200     }
00201   }
00202 
00203   // For compatibility with older versions.
00204   bool operator()( const pat::Muon & muon, pat::strbitset & ret ) 
00205   { 
00206 
00207     if (version_ == SPRING10 || version_ == FALL10 ) throw cms::Exception("LogicError") 
00208       << "MuonVPlusJetsSelectionFunctor SPRING10 and FALL10 versions needs the event! Call operator()(muon,event,ret)"
00209       <<std::endl;
00210 
00211     else if ( version_ == SUMMER08 ) return summer08Cuts( muon, ret );
00212     else if ( version_ == FIRSTDATA ) return firstDataCuts( muon, ret );
00213     else {
00214       return false;
00215     }
00216   }
00217 
00218 
00219   using Selector<pat::Muon>::operator();
00220 
00221   // cuts based on craft 08 analysis. 
00222   bool summer08Cuts( const pat::Muon & muon, pat::strbitset & ret)
00223   {
00224 
00225     ret.set(false);
00226 
00227     double norm_chi2 = muon.normChi2();
00228     double corr_d0 = muon.dB();
00229     int nhits = static_cast<int>( muon.numberOfValidHits() );
00230     
00231     double ecalVeto = muon.isolationR03().emVetoEt;
00232     double hcalVeto = muon.isolationR03().hadVetoEt;
00233         
00234     double hcalIso = muon.hcalIso();
00235     double ecalIso = muon.ecalIso();
00236     double trkIso  = muon.trackIso();
00237     double pt      = muon.pt() ;
00238 
00239     double relIso = (ecalIso + hcalIso + trkIso) / pt;
00240 
00241     if ( norm_chi2     <  cut(indexChi2_,   double()) || ignoreCut(indexChi2_)    ) passCut(ret, indexChi2_   );
00242     if ( fabs(corr_d0) <  cut(indexD0_,     double()) || ignoreCut(indexD0_)      ) passCut(ret, indexD0_     );
00243     if ( nhits         >= cut(indexNHits_,  int()   ) || ignoreCut(indexNHits_)   ) passCut(ret, indexNHits_  );
00244     if ( hcalVeto      <  cut(indexHCalVeto_,double())|| ignoreCut(indexHCalVeto_)) passCut(ret, indexHCalVeto_);
00245     if ( ecalVeto      <  cut(indexECalVeto_,double())|| ignoreCut(indexECalVeto_)) passCut(ret, indexECalVeto_);
00246     if ( relIso        <  cut(indexRelIso_, double()) || ignoreCut(indexRelIso_)  ) passCut(ret, indexRelIso_ );
00247 
00248     setIgnored(ret);
00249 
00250     return (bool)ret;
00251   }
00252 
00253 
00254 
00255   // cuts based on craft 08 analysis. 
00256   bool firstDataCuts( const pat::Muon & muon, pat::strbitset & ret)
00257   {
00258 
00259     ret.set(false);
00260 
00261     double norm_chi2 = muon.normChi2();
00262     double corr_d0 = muon.dB();
00263     double corr_ed0 = muon.edB();
00264     double corr_sd0 = ( corr_ed0 > 0.000000001 ) ? corr_d0 / corr_ed0 : 999.0;
00265     int nhits = static_cast<int>( muon.numberOfValidHits() );
00266     
00267     double ecalVeto = muon.isolationR03().emVetoEt;
00268     double hcalVeto = muon.isolationR03().hadVetoEt;
00269         
00270     double hcalIso = muon.hcalIso();
00271     double ecalIso = muon.ecalIso();
00272     double trkIso  = muon.trackIso();
00273     double pt      = muon.pt() ;
00274 
00275     double relIso = (ecalIso + hcalIso + trkIso) / pt;
00276 
00277     if ( norm_chi2     <  cut(indexChi2_,   double()) || ignoreCut(indexChi2_)    ) passCut(ret, indexChi2_   );
00278     if ( fabs(corr_d0) <  cut(indexD0_,     double()) || ignoreCut(indexD0_)      ) passCut(ret, indexD0_     );
00279     if ( fabs(corr_ed0)<  cut(indexED0_,    double()) || ignoreCut(indexED0_)     ) passCut(ret, indexED0_    );
00280     if ( fabs(corr_sd0)<  cut(indexSD0_,    double()) || ignoreCut(indexSD0_)     ) passCut(ret, indexSD0_    );
00281     if ( nhits         >= cut(indexNHits_,  int()   ) || ignoreCut(indexNHits_)   ) passCut(ret, indexNHits_  );
00282     if ( hcalVeto      <  cut(indexHCalVeto_,double())|| ignoreCut(indexHCalVeto_)) passCut(ret, indexHCalVeto_);
00283     if ( ecalVeto      <  cut(indexECalVeto_,double())|| ignoreCut(indexECalVeto_)) passCut(ret, indexECalVeto_);
00284     if ( relIso        <  cut(indexRelIso_, double()) || ignoreCut(indexRelIso_)  ) passCut(ret, indexRelIso_ );
00285 
00286     setIgnored(ret);
00287     
00288     return (bool)ret;
00289   }
00290 
00291   // cuts based on top group L+J synchronization exercise
00292   bool spring10Cuts( const pat::Muon & muon, edm::EventBase const & event, pat::strbitset & ret)
00293   {
00294 
00295     ret.set(false);
00296 
00297     double norm_chi2 = muon.normChi2();
00298     double corr_d0 = muon.dB();
00299     double corr_ed0 = muon.edB();
00300     double corr_sd0 = ( corr_ed0 > 0.000000001 ) ? corr_d0 / corr_ed0 : 999.0;
00301 
00302     //If required, recalculate the impact parameter using the beam spot
00303     if (recalcDBFromBSp_) {
00304 
00305       //Get the beam spot
00306       reco::TrackBase::Point beamPoint(0,0,0);
00307       reco::BeamSpot beamSpot;
00308       edm::Handle<reco::BeamSpot> beamSpotHandle;
00309       event.getByLabel(beamLineSrc_, beamSpotHandle);
00310       
00311       if( beamSpotHandle.isValid() ){
00312         beamSpot = *beamSpotHandle;
00313       } else{
00314         edm::LogError("DataNotAvailable")
00315           << "No beam spot available from EventSetup, not adding high level selection \n";
00316       }
00317       beamPoint = reco::TrackBase::Point ( beamSpot.x0(), beamSpot.y0(), beamSpot.z0() );
00318       
00319       //Use the beamspot to correct the impact parameter and uncertainty
00320       reco::TrackRef innerTrack = muon.innerTrack();
00321       if ( innerTrack.isNonnull() && innerTrack.isAvailable() ) {
00322         corr_d0 = -1.0 * innerTrack->dxy( beamPoint );
00323         corr_ed0 = sqrt( innerTrack->d0Error() * innerTrack->d0Error() 
00324                          + 0.5* beamSpot.BeamWidthX()*beamSpot.BeamWidthX() 
00325                          + 0.5* beamSpot.BeamWidthY()*beamSpot.BeamWidthY() );
00326         corr_sd0 = ( corr_ed0 > 0.000000001 ) ? corr_d0 / corr_ed0 : 999.0;
00327 
00328       } else {
00329         corr_d0 =  999.;
00330         corr_ed0 = 999.;
00331       }
00332     }
00333 
00334     int nhits = static_cast<int>( muon.numberOfValidHits() );
00335     int nValidMuonHits = static_cast<int> (muon.globalTrack()->hitPattern().numberOfValidMuonHits());
00336     
00337     double ecalVeto = muon.isolationR03().emVetoEt;
00338     double hcalVeto = muon.isolationR03().hadVetoEt;
00339         
00340     double hcalIso = muon.hcalIso();
00341     double ecalIso = muon.ecalIso();
00342     double trkIso  = muon.trackIso();
00343     double pt      = muon.pt() ;
00344 
00345     double relIso = (ecalIso + hcalIso + trkIso) / pt;
00346 
00347     if ( norm_chi2     <  cut(indexChi2_,   double()) || ignoreCut(indexChi2_)    ) passCut(ret, indexChi2_   );
00348     if ( fabs(corr_d0) <  cut(indexD0_,     double()) || ignoreCut(indexD0_)      ) passCut(ret, indexD0_     );
00349     if ( fabs(corr_ed0)<  cut(indexED0_,    double()) || ignoreCut(indexED0_)     ) passCut(ret, indexED0_    );
00350     if ( fabs(corr_sd0)<  cut(indexSD0_,    double()) || ignoreCut(indexSD0_)     ) passCut(ret, indexSD0_    );
00351     if ( nhits         >= cut(indexNHits_,  int()   ) || ignoreCut(indexNHits_)   ) passCut(ret, indexNHits_  );
00352     if ( nValidMuonHits> cut(indexNValMuHits_,int()) || ignoreCut(indexNValMuHits_)) passCut(ret, indexNValMuHits_  );
00353     if ( hcalVeto      <  cut(indexHCalVeto_,double())|| ignoreCut(indexHCalVeto_)) passCut(ret, indexHCalVeto_);
00354     if ( ecalVeto      <  cut(indexECalVeto_,double())|| ignoreCut(indexECalVeto_)) passCut(ret, indexECalVeto_);
00355     if ( relIso        <  cut(indexRelIso_, double()) || ignoreCut(indexRelIso_)  ) passCut(ret, indexRelIso_ );
00356 
00357     setIgnored(ret);
00358     
00359     return (bool)ret;
00360   }
00361 
00362 
00363 
00364 
00365   // cuts based on top group L+J synchronization exercise
00366   bool fall10Cuts( const pat::Muon & muon, edm::EventBase const & event, pat::strbitset & ret)
00367   {
00368 
00369     ret.set(false);
00370 
00371     double norm_chi2 = muon.normChi2();
00372     double corr_d0 = muon.dB();
00373     double corr_ed0 = muon.edB();
00374     double corr_sd0 = ( corr_ed0 > 0.000000001 ) ? corr_d0 / corr_ed0 : 999.0;
00375 
00376     // Get the PV for the muon z requirement
00377     edm::Handle<std::vector<reco::Vertex> > pvtxHandle_;
00378     event.getByLabel( pvSrc_, pvtxHandle_ );
00379 
00380     double zvtx = -999;
00381     if ( pvtxHandle_->size() > 0 ) {
00382       zvtx = pvtxHandle_->at(0).z();
00383     } else {
00384       throw cms::Exception("InvalidInput") << " There needs to be at least one primary vertex in the event." << std::endl;
00385     }
00386 
00387     //If required, recalculate the impact parameter using the beam spot
00388     if (recalcDBFromBSp_) {
00389 
00390       //Get the beam spot
00391       reco::TrackBase::Point beamPoint(0,0,0);
00392       reco::BeamSpot beamSpot;
00393       edm::Handle<reco::BeamSpot> beamSpotHandle;
00394       event.getByLabel(beamLineSrc_, beamSpotHandle);
00395       
00396       if( beamSpotHandle.isValid() ){
00397         beamSpot = *beamSpotHandle;
00398       } else{
00399         edm::LogError("DataNotAvailable")
00400           << "No beam spot available from EventSetup, not adding high level selection \n";
00401       }
00402       beamPoint = reco::TrackBase::Point ( beamSpot.x0(), beamSpot.y0(), beamSpot.z0() );
00403       
00404       //Use the beamspot to correct the impact parameter and uncertainty
00405       reco::TrackRef innerTrack = muon.innerTrack();
00406       if ( innerTrack.isNonnull() && innerTrack.isAvailable() ) {
00407         corr_d0 = -1.0 * innerTrack->dxy( beamPoint );
00408         corr_ed0 = sqrt( innerTrack->d0Error() * innerTrack->d0Error() 
00409                          + 0.5* beamSpot.BeamWidthX()*beamSpot.BeamWidthX() 
00410                          + 0.5* beamSpot.BeamWidthY()*beamSpot.BeamWidthY() );
00411         corr_sd0 = ( corr_ed0 > 0.000000001 ) ? corr_d0 / corr_ed0 : 999.0;
00412 
00413       } else {
00414         corr_d0 =  999.;
00415         corr_ed0 = 999.;
00416       }
00417     }
00418 
00419     int nhits = static_cast<int>( muon.numberOfValidHits() );
00420     int nValidMuonHits = static_cast<int> (muon.globalTrack()->hitPattern().numberOfValidMuonHits());
00421     
00422     double ecalVeto = muon.isolationR03().emVetoEt;
00423     double hcalVeto = muon.isolationR03().hadVetoEt;
00424         
00425     double hcalIso = muon.hcalIso();
00426     double ecalIso = muon.ecalIso();
00427     double trkIso  = muon.trackIso();
00428     double pt      = muon.pt() ;
00429 
00430     double relIso = (ecalIso + hcalIso + trkIso) / pt;
00431 
00432 
00433     double z_mu = muon.vertex().z();
00434 
00435     int nPixelHits = muon.innerTrack()->hitPattern().pixelLayersWithMeasurement();
00436 
00437     int nMatchedStations = muon.numberOfMatches();
00438 
00439     if ( norm_chi2     <  cut(indexChi2_,   double()) || ignoreCut(indexChi2_)    ) passCut(ret, indexChi2_   );
00440     if ( fabs(corr_d0) <  cut(indexD0_,     double()) || ignoreCut(indexD0_)      ) passCut(ret, indexD0_     );
00441     if ( fabs(corr_ed0)<  cut(indexED0_,    double()) || ignoreCut(indexED0_)     ) passCut(ret, indexED0_    );
00442     if ( fabs(corr_sd0)<  cut(indexSD0_,    double()) || ignoreCut(indexSD0_)     ) passCut(ret, indexSD0_    );
00443     if ( nhits         >= cut(indexNHits_,  int()   ) || ignoreCut(indexNHits_)   ) passCut(ret, indexNHits_  );
00444     if ( nValidMuonHits> cut(indexNValMuHits_,int()) || ignoreCut(indexNValMuHits_)) passCut(ret, indexNValMuHits_  );
00445     if ( hcalVeto      <  cut(indexHCalVeto_,double())|| ignoreCut(indexHCalVeto_)) passCut(ret, indexHCalVeto_);
00446     if ( ecalVeto      <  cut(indexECalVeto_,double())|| ignoreCut(indexECalVeto_)) passCut(ret, indexECalVeto_);
00447     if ( relIso        <  cut(indexRelIso_, double()) || ignoreCut(indexRelIso_)  ) passCut(ret, indexRelIso_ );
00448     if ( fabs(z_mu-zvtx)<  cut(indexLepZ_, double()) || ignoreCut(indexLepZ_)  ) passCut(ret, indexLepZ_ );
00449     if ( nPixelHits    >  cut(indexPixHits_,int())    || ignoreCut(indexPixHits_))  passCut(ret, indexPixHits_);
00450     if ( nMatchedStations> cut(indexStations_,int())    || ignoreCut(indexStations_))  passCut(ret, indexStations_);
00451 
00452     setIgnored(ret);
00453     
00454     return (bool)ret;
00455   }
00456 
00457 
00458   // cuts based on top group L+J synchronization exercise
00459   // this is a copy of fall 10 cuts
00460   // with a hack to include a double-sided reliso cut
00461   
00462   bool kitQCDCuts( const pat::Muon & muon, edm::EventBase const & event, pat::strbitset & ret)
00463   {
00464 
00465     ret.set(false);
00466 
00467     double norm_chi2 = muon.normChi2();
00468     double corr_d0 = muon.dB();
00469     double corr_ed0 = muon.edB();
00470     double corr_sd0 = ( corr_ed0 > 0.000000001 ) ? corr_d0 / corr_ed0 : 999.0;
00471 
00472     // Get the PV for the muon z requirement
00473     edm::Handle<std::vector<reco::Vertex> > pvtxHandle_;
00474     event.getByLabel( pvSrc_, pvtxHandle_ );
00475 
00476     double zvtx = -999;
00477     if ( pvtxHandle_->size() > 0 ) {
00478       zvtx = pvtxHandle_->at(0).z();
00479     } else {
00480       throw cms::Exception("InvalidInput") << " There needs to be at least one primary vertex in the event." << std::endl;
00481     }
00482 
00483     //If required, recalculate the impact parameter using the beam spot
00484     if (recalcDBFromBSp_) {
00485 
00486       //Get the beam spot
00487       reco::TrackBase::Point beamPoint(0,0,0);
00488       reco::BeamSpot beamSpot;
00489       edm::Handle<reco::BeamSpot> beamSpotHandle;
00490       event.getByLabel(beamLineSrc_, beamSpotHandle);
00491       
00492       if( beamSpotHandle.isValid() ){
00493         beamSpot = *beamSpotHandle;
00494       } else{
00495         edm::LogError("DataNotAvailable")
00496           << "No beam spot available from EventSetup, not adding high level selection \n";
00497       }
00498       beamPoint = reco::TrackBase::Point ( beamSpot.x0(), beamSpot.y0(), beamSpot.z0() );
00499       
00500       //Use the beamspot to correct the impact parameter and uncertainty
00501       reco::TrackRef innerTrack = muon.innerTrack();
00502       if ( innerTrack.isNonnull() && innerTrack.isAvailable() ) {
00503         corr_d0 = -1.0 * innerTrack->dxy( beamPoint );
00504         corr_ed0 = sqrt( innerTrack->d0Error() * innerTrack->d0Error() 
00505                          + 0.5* beamSpot.BeamWidthX()*beamSpot.BeamWidthX() 
00506                          + 0.5* beamSpot.BeamWidthY()*beamSpot.BeamWidthY() );
00507         corr_sd0 = ( corr_ed0 > 0.000000001 ) ? corr_d0 / corr_ed0 : 999.0;
00508 
00509       } else {
00510         corr_d0 =  999.;
00511         corr_ed0 = 999.;
00512       }
00513     }
00514 
00515     int nhits = static_cast<int>( muon.numberOfValidHits() );
00516     int nValidMuonHits = static_cast<int> (muon.globalTrack()->hitPattern().numberOfValidMuonHits());
00517     
00518     double ecalVeto = muon.isolationR03().emVetoEt;
00519     double hcalVeto = muon.isolationR03().hadVetoEt;
00520         
00521     double hcalIso = muon.hcalIso();
00522     double ecalIso = muon.ecalIso();
00523     double trkIso  = muon.trackIso();
00524     double pt      = muon.pt() ;
00525 
00526     double relIso = (ecalIso + hcalIso + trkIso) / pt;
00527 
00528 
00529     double z_mu = muon.vertex().z();
00530 
00531     int nPixelHits = muon.innerTrack()->hitPattern().pixelLayersWithMeasurement();
00532 
00533     int nMatchedStations = muon.numberOfMatches();
00534 
00535     if ( norm_chi2     <  cut(indexChi2_,   double()) || ignoreCut(indexChi2_)    ) passCut(ret, indexChi2_   );
00536     if ( fabs(corr_d0) <  cut(indexD0_,     double()) || ignoreCut(indexD0_)      ) passCut(ret, indexD0_     );
00537     if ( fabs(corr_ed0)<  cut(indexED0_,    double()) || ignoreCut(indexED0_)     ) passCut(ret, indexED0_    );
00538     if ( fabs(corr_sd0)<  cut(indexSD0_,    double()) || ignoreCut(indexSD0_)     ) passCut(ret, indexSD0_    );
00539     if ( nhits         >= cut(indexNHits_,  int()   ) || ignoreCut(indexNHits_)   ) passCut(ret, indexNHits_  );
00540     if ( nValidMuonHits> cut(indexNValMuHits_,int()) || ignoreCut(indexNValMuHits_)) passCut(ret, indexNValMuHits_  );
00541     if ( hcalVeto      <  cut(indexHCalVeto_,double())|| ignoreCut(indexHCalVeto_)) passCut(ret, indexHCalVeto_);
00542     if ( ecalVeto      <  cut(indexECalVeto_,double())|| ignoreCut(indexECalVeto_)) passCut(ret, indexECalVeto_);
00543     if ( fabs(z_mu-zvtx)<  cut(indexLepZ_, double()) || ignoreCut(indexLepZ_)  ) passCut(ret, indexLepZ_ );
00544     if ( nPixelHits    >  cut(indexPixHits_,int())    || ignoreCut(indexPixHits_))  passCut(ret, indexPixHits_);
00545     if ( nMatchedStations> cut(indexStations_,int())    || ignoreCut(indexStations_))  passCut(ret, indexStations_);
00546 
00547 
00549     //
00550     // JMS Dec 13 2010
00551     // HACK 
00552     // Need double-sided relIso cut to implement data-driven QCD
00553     //
00554     //
00556     if ( ((relIso > 0.2) && (relIso < 0.75))
00557          || ignoreCut(indexRelIso_)  ) passCut(ret, indexRelIso_ );
00558 
00559 
00560 
00561     
00562     setIgnored(ret);
00563     
00564     return (bool)ret;
00565   }
00566 
00567 
00568   
00569   
00570  private: // member variables
00571   
00572   Version_t version_;
00573   bool recalcDBFromBSp_;
00574   edm::InputTag beamLineSrc_;
00575   edm::InputTag pvSrc_;
00576 
00577   index_type indexChi2_;          
00578   index_type indexD0_;            
00579   index_type indexED0_;           
00580   index_type indexSD0_;           
00581   index_type indexNHits_;         
00582   index_type indexNValMuHits_;    
00583   index_type indexECalVeto_;      
00584   index_type indexHCalVeto_;      
00585   index_type indexRelIso_;        
00586   index_type indexLepZ_;
00587   index_type indexPixHits_;
00588   index_type indexStations_;
00589 
00590 
00591 };
00592 
00593 #endif