Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:24:09 

0001 #ifndef PhysicsTools_PatUtils_interface_MuonVPlusJetsIDSelectionFunctor_h
0002 #define PhysicsTools_PatUtils_interface_MuonVPlusJetsIDSelectionFunctor_h
0003 
0004 #ifndef __GCCXML__
0005 #include "FWCore/Framework/interface/ConsumesCollector.h"
0006 #endif
0007 #include "DataFormats/PatCandidates/interface/Muon.h"
0008 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0009 #include "DataFormats/VertexReco/interface/Vertex.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 
0012 #include "PhysicsTools/SelectorUtils/interface/Selector.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 
0015 #include <iostream>
0016 
0017 class MuonVPlusJetsIDSelectionFunctor : public Selector<pat::Muon> {
0018 public:  // interface
0019   bool verbose_;
0020 
0021   enum Version_t { SUMMER08, FIRSTDATA, SPRING10, FALL10, N_VERSIONS, KITQCD };
0022 
0023   MuonVPlusJetsIDSelectionFunctor() {}
0024 
0025 #ifndef __GCCXML__
0026   MuonVPlusJetsIDSelectionFunctor(edm::ParameterSet const& parameters, edm::ConsumesCollector& iC)
0027       : MuonVPlusJetsIDSelectionFunctor(parameters) {
0028     beamLineSrcToken_ = iC.consumes<reco::BeamSpot>(beamLineSrc_);
0029     pvSrcToken_ = iC.consumes<std::vector<reco::Vertex> >(pvSrc_);
0030   }
0031 #endif
0032 
0033   MuonVPlusJetsIDSelectionFunctor(edm::ParameterSet const& parameters) {
0034     verbose_ = false;
0035 
0036     std::string versionStr = parameters.getParameter<std::string>("version");
0037 
0038     Version_t version = N_VERSIONS;
0039 
0040     if (versionStr == "SUMMER08") {
0041       version = SUMMER08;
0042     } else if (versionStr == "FIRSTDATA") {
0043       version = FIRSTDATA;
0044     } else if (versionStr == "SPRING10") {
0045       version = SPRING10;
0046     } else if (versionStr == "FALL10") {
0047       version = FALL10;
0048       if (verbose_)
0049         std::cout << "\nMUON SELECTION - you are using FALL10 Selection" << std::endl;
0050     } else if (versionStr == "KITQCD") {
0051       version = KITQCD;
0052       if (verbose_)
0053         std::cout << "\nMUON SELECTION - you are using KITQCD Selection" << std::endl;
0054     } else {
0055       throw cms::Exception("InvalidInput")
0056           << "Expect version to be one of SUMMER08, FIRSTDATA, SPRING10, FALL10" << std::endl;
0057     }
0058 
0059     initialize(version,
0060                parameters.getParameter<double>("Chi2"),
0061                parameters.getParameter<double>("D0"),
0062                parameters.getParameter<double>("ED0"),
0063                parameters.getParameter<double>("SD0"),
0064                parameters.getParameter<int>("NHits"),
0065                parameters.getParameter<int>("NValMuHits"),
0066                parameters.getParameter<double>("ECalVeto"),
0067                parameters.getParameter<double>("HCalVeto"),
0068                parameters.getParameter<double>("RelIso"),
0069                parameters.getParameter<double>("LepZ"),
0070                parameters.getParameter<int>("nPixelHits"),
0071                parameters.getParameter<int>("nMatchedStations"));
0072     if (parameters.exists("cutsToIgnore"))
0073       setIgnoredCuts(parameters.getParameter<std::vector<std::string> >("cutsToIgnore"));
0074 
0075     retInternal_ = getBitTemplate();
0076 
0077     recalcDBFromBSp_ = parameters.getParameter<bool>("RecalcFromBeamSpot");
0078     pvSrc_ = parameters.getParameter<edm::InputTag>("pvSrc");
0079   }
0080 
0081   MuonVPlusJetsIDSelectionFunctor(Version_t version,
0082                                   double chi2 = 10.0,
0083                                   double d0 = 0.2,
0084                                   double ed0 = 999.0,
0085                                   double sd0 = 999.0,
0086                                   int nhits = 11,
0087                                   int nValidMuonHits = 0,
0088                                   double ecalveto = 4.0,
0089                                   double hcalveto = 6.0,
0090                                   double reliso = 0.05,
0091                                   double maxLepZ = 1.0,
0092                                   int minPixelHits = 1,
0093                                   int minNMatches = 1)
0094       : recalcDBFromBSp_(false) {
0095     initialize(version,
0096                chi2,
0097                d0,
0098                ed0,
0099                sd0,
0100                nhits,
0101                nValidMuonHits,
0102                ecalveto,
0103                hcalveto,
0104                reliso,
0105                maxLepZ,
0106                minPixelHits,
0107                minNMatches);
0108 
0109     retInternal_ = getBitTemplate();
0110   }
0111 
0112   void initialize(Version_t version,
0113                   double chi2 = 10.0,
0114                   double d0 = 999.0,
0115                   double ed0 = 999.0,
0116                   double sd0 = 3.0,
0117                   int nhits = 11,
0118                   int nValidMuonHits = 0,
0119                   double ecalveto = 4.0,
0120                   double hcalveto = 6.0,
0121                   double reliso = 0.05,
0122                   double maxLepZ = 1.0,
0123                   int minPixelHits = 1,
0124                   int minNMatches = 1) {
0125     version_ = version;
0126 
0127     push_back("Chi2", chi2);
0128     push_back("D0", d0);
0129     push_back("ED0", ed0);
0130     push_back("SD0", sd0);
0131     push_back("NHits", nhits);
0132     push_back("NValMuHits", nValidMuonHits);
0133     push_back("ECalVeto", ecalveto);
0134     push_back("HCalVeto", hcalveto);
0135     push_back("RelIso", reliso);
0136     push_back("LepZ", maxLepZ);
0137     push_back("nPixelHits", minPixelHits);
0138     push_back("nMatchedStations", minNMatches);
0139 
0140     set("Chi2");
0141     set("D0");
0142     set("ED0");
0143     set("SD0");
0144     set("NHits");
0145     set("NValMuHits");
0146     set("ECalVeto");
0147     set("HCalVeto");
0148     set("RelIso");
0149     set("LepZ");
0150     set("nPixelHits");
0151     set("nMatchedStations");
0152 
0153     indexChi2_ = index_type(&bits_, "Chi2");
0154     indexD0_ = index_type(&bits_, "D0");
0155     indexED0_ = index_type(&bits_, "ED0");
0156     indexSD0_ = index_type(&bits_, "SD0");
0157     indexNHits_ = index_type(&bits_, "NHits");
0158     indexNValMuHits_ = index_type(&bits_, "NValMuHits");
0159     indexECalVeto_ = index_type(&bits_, "ECalVeto");
0160     indexHCalVeto_ = index_type(&bits_, "HCalVeto");
0161     indexRelIso_ = index_type(&bits_, "RelIso");
0162     indexLepZ_ = index_type(&bits_, "LepZ");
0163     indexPixHits_ = index_type(&bits_, "nPixelHits");
0164     indexStations_ = index_type(&bits_, "nMatchedStations");
0165 
0166     if (version == FALL10) {
0167       set("ED0", false);
0168       set("SD0", false);
0169       set("ECalVeto", false);
0170       set("HCalVeto", false);
0171     } else if (version == SPRING10) {
0172       set("ED0", false);
0173       set("SD0", false);
0174       set("ECalVeto", false);
0175       set("HCalVeto", false);
0176       set("LepZ", false);
0177       set("nPixelHits", false);
0178       set("nMatchedStations", false);
0179     } else if (version_ == FIRSTDATA) {
0180       set("D0", false);
0181       set("ED0", false);
0182       set("NValMuHits", false);
0183       set("LepZ", false);
0184       set("nPixelHits", false);
0185       set("nMatchedStations", false);
0186     } else if (version == SUMMER08) {
0187       set("SD0", false);
0188       set("NValMuHits", false);
0189       set("LepZ", false);
0190       set("nPixelHits", false);
0191       set("nMatchedStations", false);
0192     }
0193   }
0194 
0195   // Allow for multiple definitions of the cuts.
0196   bool operator()(const pat::Muon& muon, edm::EventBase const& event, pat::strbitset& ret) override {
0197     if (version_ == FALL10)
0198       return fall10Cuts(muon, event, ret);
0199     else if (version_ == SPRING10)
0200       return spring10Cuts(muon, event, ret);
0201     else if (version_ == SUMMER08)
0202       return summer08Cuts(muon, ret);
0203     else if (version_ == FIRSTDATA)
0204       return firstDataCuts(muon, ret);
0205     else if (version_ == KITQCD) {
0206       if (verbose_)
0207         std::cout << "Calling KIT selection method" << std::endl;
0208       return kitQCDCuts(muon, event, ret);
0209     } else {
0210       return false;
0211     }
0212   }
0213 
0214   // For compatibility with older versions.
0215   bool operator()(const pat::Muon& muon, pat::strbitset& ret) override {
0216     if (version_ == SPRING10 || version_ == FALL10)
0217       throw cms::Exception("LogicError") << "MuonVPlusJetsSelectionFunctor SPRING10 and FALL10 versions needs the "
0218                                             "event! Call operator()(muon,event,ret)"
0219                                          << std::endl;
0220 
0221     else if (version_ == SUMMER08)
0222       return summer08Cuts(muon, ret);
0223     else if (version_ == FIRSTDATA)
0224       return firstDataCuts(muon, ret);
0225     else {
0226       return false;
0227     }
0228   }
0229 
0230   using Selector<pat::Muon>::operator();
0231 
0232   // cuts based on craft 08 analysis.
0233   bool summer08Cuts(const pat::Muon& muon, pat::strbitset& ret) {
0234     ret.set(false);
0235 
0236     double norm_chi2 = muon.normChi2();
0237     double corr_d0 = muon.dB();
0238     int nhits = static_cast<int>(muon.numberOfValidHits());
0239 
0240     double ecalVeto = muon.isolationR03().emVetoEt;
0241     double hcalVeto = muon.isolationR03().hadVetoEt;
0242 
0243     double hcalIso = muon.hcalIso();
0244     double ecalIso = muon.ecalIso();
0245     double trkIso = muon.trackIso();
0246     double pt = muon.pt();
0247 
0248     double relIso = (ecalIso + hcalIso + trkIso) / pt;
0249 
0250     if (norm_chi2 < cut(indexChi2_, double()) || ignoreCut(indexChi2_))
0251       passCut(ret, indexChi2_);
0252     if (fabs(corr_d0) < cut(indexD0_, double()) || ignoreCut(indexD0_))
0253       passCut(ret, indexD0_);
0254     if (nhits >= cut(indexNHits_, int()) || ignoreCut(indexNHits_))
0255       passCut(ret, indexNHits_);
0256     if (hcalVeto < cut(indexHCalVeto_, double()) || ignoreCut(indexHCalVeto_))
0257       passCut(ret, indexHCalVeto_);
0258     if (ecalVeto < cut(indexECalVeto_, double()) || ignoreCut(indexECalVeto_))
0259       passCut(ret, indexECalVeto_);
0260     if (relIso < cut(indexRelIso_, double()) || ignoreCut(indexRelIso_))
0261       passCut(ret, indexRelIso_);
0262 
0263     setIgnored(ret);
0264 
0265     return (bool)ret;
0266   }
0267 
0268   // cuts based on craft 08 analysis.
0269   bool firstDataCuts(const pat::Muon& muon, pat::strbitset& ret) {
0270     ret.set(false);
0271 
0272     double norm_chi2 = muon.normChi2();
0273     double corr_d0 = muon.dB();
0274     double corr_ed0 = muon.edB();
0275     double corr_sd0 = (corr_ed0 > 0.000000001) ? corr_d0 / corr_ed0 : 999.0;
0276     int nhits = static_cast<int>(muon.numberOfValidHits());
0277 
0278     double ecalVeto = muon.isolationR03().emVetoEt;
0279     double hcalVeto = muon.isolationR03().hadVetoEt;
0280 
0281     double hcalIso = muon.hcalIso();
0282     double ecalIso = muon.ecalIso();
0283     double trkIso = muon.trackIso();
0284     double pt = muon.pt();
0285 
0286     double relIso = (ecalIso + hcalIso + trkIso) / pt;
0287 
0288     if (norm_chi2 < cut(indexChi2_, double()) || ignoreCut(indexChi2_))
0289       passCut(ret, indexChi2_);
0290     if (fabs(corr_d0) < cut(indexD0_, double()) || ignoreCut(indexD0_))
0291       passCut(ret, indexD0_);
0292     if (fabs(corr_ed0) < cut(indexED0_, double()) || ignoreCut(indexED0_))
0293       passCut(ret, indexED0_);
0294     if (fabs(corr_sd0) < cut(indexSD0_, double()) || ignoreCut(indexSD0_))
0295       passCut(ret, indexSD0_);
0296     if (nhits >= cut(indexNHits_, int()) || ignoreCut(indexNHits_))
0297       passCut(ret, indexNHits_);
0298     if (hcalVeto < cut(indexHCalVeto_, double()) || ignoreCut(indexHCalVeto_))
0299       passCut(ret, indexHCalVeto_);
0300     if (ecalVeto < cut(indexECalVeto_, double()) || ignoreCut(indexECalVeto_))
0301       passCut(ret, indexECalVeto_);
0302     if (relIso < cut(indexRelIso_, double()) || ignoreCut(indexRelIso_))
0303       passCut(ret, indexRelIso_);
0304 
0305     setIgnored(ret);
0306 
0307     return (bool)ret;
0308   }
0309 
0310   // cuts based on top group L+J synchronization exercise
0311   bool spring10Cuts(const pat::Muon& muon, edm::EventBase const& event, pat::strbitset& ret) {
0312     ret.set(false);
0313 
0314     double norm_chi2 = muon.normChi2();
0315     double corr_d0 = muon.dB();
0316     double corr_ed0 = muon.edB();
0317     double corr_sd0 = (corr_ed0 > 0.000000001) ? corr_d0 / corr_ed0 : 999.0;
0318 
0319     //If required, recalculate the impact parameter using the beam spot
0320     if (recalcDBFromBSp_) {
0321       //Get the beam spot
0322       reco::TrackBase::Point beamPoint(0, 0, 0);
0323       reco::BeamSpot beamSpot;
0324       edm::Handle<reco::BeamSpot> beamSpotHandle;
0325       event.getByLabel(beamLineSrc_, beamSpotHandle);
0326 
0327       if (beamSpotHandle.isValid()) {
0328         beamSpot = *beamSpotHandle;
0329       } else {
0330         edm::LogError("DataNotAvailable")
0331             << "No beam spot available from EventSetup, not adding high level selection \n";
0332       }
0333       beamPoint = reco::TrackBase::Point(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
0334 
0335       //Use the beamspot to correct the impact parameter and uncertainty
0336       reco::TrackRef innerTrack = muon.innerTrack();
0337       if (innerTrack.isNonnull() && innerTrack.isAvailable()) {
0338         corr_d0 = -1.0 * innerTrack->dxy(beamPoint);
0339         corr_ed0 =
0340             sqrt(innerTrack->d0Error() * innerTrack->d0Error() + 0.5 * beamSpot.BeamWidthX() * beamSpot.BeamWidthX() +
0341                  0.5 * beamSpot.BeamWidthY() * beamSpot.BeamWidthY());
0342         corr_sd0 = (corr_ed0 > 0.000000001) ? corr_d0 / corr_ed0 : 999.0;
0343 
0344       } else {
0345         corr_d0 = 999.;
0346         corr_ed0 = 999.;
0347       }
0348     }
0349 
0350     int nhits = static_cast<int>(muon.numberOfValidHits());
0351     int nValidMuonHits = static_cast<int>(muon.globalTrack()->hitPattern().numberOfValidMuonHits());
0352 
0353     double ecalVeto = muon.isolationR03().emVetoEt;
0354     double hcalVeto = muon.isolationR03().hadVetoEt;
0355 
0356     double hcalIso = muon.hcalIso();
0357     double ecalIso = muon.ecalIso();
0358     double trkIso = muon.trackIso();
0359     double pt = muon.pt();
0360 
0361     double relIso = (ecalIso + hcalIso + trkIso) / pt;
0362 
0363     if (norm_chi2 < cut(indexChi2_, double()) || ignoreCut(indexChi2_))
0364       passCut(ret, indexChi2_);
0365     if (fabs(corr_d0) < cut(indexD0_, double()) || ignoreCut(indexD0_))
0366       passCut(ret, indexD0_);
0367     if (fabs(corr_ed0) < cut(indexED0_, double()) || ignoreCut(indexED0_))
0368       passCut(ret, indexED0_);
0369     if (fabs(corr_sd0) < cut(indexSD0_, double()) || ignoreCut(indexSD0_))
0370       passCut(ret, indexSD0_);
0371     if (nhits >= cut(indexNHits_, int()) || ignoreCut(indexNHits_))
0372       passCut(ret, indexNHits_);
0373     if (nValidMuonHits > cut(indexNValMuHits_, int()) || ignoreCut(indexNValMuHits_))
0374       passCut(ret, indexNValMuHits_);
0375     if (hcalVeto < cut(indexHCalVeto_, double()) || ignoreCut(indexHCalVeto_))
0376       passCut(ret, indexHCalVeto_);
0377     if (ecalVeto < cut(indexECalVeto_, double()) || ignoreCut(indexECalVeto_))
0378       passCut(ret, indexECalVeto_);
0379     if (relIso < cut(indexRelIso_, double()) || ignoreCut(indexRelIso_))
0380       passCut(ret, indexRelIso_);
0381 
0382     setIgnored(ret);
0383 
0384     return (bool)ret;
0385   }
0386 
0387   // cuts based on top group L+J synchronization exercise
0388   bool fall10Cuts(const pat::Muon& muon, edm::EventBase const& event, pat::strbitset& ret) {
0389     ret.set(false);
0390 
0391     double norm_chi2 = muon.normChi2();
0392     double corr_d0 = muon.dB();
0393     double corr_ed0 = muon.edB();
0394     double corr_sd0 = (corr_ed0 > 0.000000001) ? corr_d0 / corr_ed0 : 999.0;
0395 
0396     // Get the PV for the muon z requirement
0397     edm::Handle<std::vector<reco::Vertex> > pvtxHandle_;
0398     event.getByLabel(pvSrc_, pvtxHandle_);
0399 
0400     double zvtx = -999;
0401     if (!pvtxHandle_->empty()) {
0402       zvtx = pvtxHandle_->at(0).z();
0403     } else {
0404       throw cms::Exception("InvalidInput")
0405           << " There needs to be at least one primary vertex in the event." << std::endl;
0406     }
0407 
0408     //If required, recalculate the impact parameter using the beam spot
0409     if (recalcDBFromBSp_) {
0410       //Get the beam spot
0411       reco::TrackBase::Point beamPoint(0, 0, 0);
0412       reco::BeamSpot beamSpot;
0413       edm::Handle<reco::BeamSpot> beamSpotHandle;
0414       event.getByLabel(beamLineSrc_, beamSpotHandle);
0415 
0416       if (beamSpotHandle.isValid()) {
0417         beamSpot = *beamSpotHandle;
0418       } else {
0419         edm::LogError("DataNotAvailable")
0420             << "No beam spot available from EventSetup, not adding high level selection \n";
0421       }
0422       beamPoint = reco::TrackBase::Point(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
0423 
0424       //Use the beamspot to correct the impact parameter and uncertainty
0425       reco::TrackRef innerTrack = muon.innerTrack();
0426       if (innerTrack.isNonnull() && innerTrack.isAvailable()) {
0427         corr_d0 = -1.0 * innerTrack->dxy(beamPoint);
0428         corr_ed0 =
0429             sqrt(innerTrack->d0Error() * innerTrack->d0Error() + 0.5 * beamSpot.BeamWidthX() * beamSpot.BeamWidthX() +
0430                  0.5 * beamSpot.BeamWidthY() * beamSpot.BeamWidthY());
0431         corr_sd0 = (corr_ed0 > 0.000000001) ? corr_d0 / corr_ed0 : 999.0;
0432 
0433       } else {
0434         corr_d0 = 999.;
0435         corr_ed0 = 999.;
0436       }
0437     }
0438 
0439     int nhits = static_cast<int>(muon.numberOfValidHits());
0440     int nValidMuonHits = static_cast<int>(muon.globalTrack()->hitPattern().numberOfValidMuonHits());
0441 
0442     double ecalVeto = muon.isolationR03().emVetoEt;
0443     double hcalVeto = muon.isolationR03().hadVetoEt;
0444 
0445     double hcalIso = muon.hcalIso();
0446     double ecalIso = muon.ecalIso();
0447     double trkIso = muon.trackIso();
0448     double pt = muon.pt();
0449 
0450     double relIso = (ecalIso + hcalIso + trkIso) / pt;
0451 
0452     double z_mu = muon.vertex().z();
0453 
0454     int nPixelHits = muon.innerTrack()->hitPattern().pixelLayersWithMeasurement();
0455 
0456     int nMatchedStations = muon.numberOfMatches();
0457 
0458     if (norm_chi2 < cut(indexChi2_, double()) || ignoreCut(indexChi2_))
0459       passCut(ret, indexChi2_);
0460     if (fabs(corr_d0) < cut(indexD0_, double()) || ignoreCut(indexD0_))
0461       passCut(ret, indexD0_);
0462     if (fabs(corr_ed0) < cut(indexED0_, double()) || ignoreCut(indexED0_))
0463       passCut(ret, indexED0_);
0464     if (fabs(corr_sd0) < cut(indexSD0_, double()) || ignoreCut(indexSD0_))
0465       passCut(ret, indexSD0_);
0466     if (nhits >= cut(indexNHits_, int()) || ignoreCut(indexNHits_))
0467       passCut(ret, indexNHits_);
0468     if (nValidMuonHits > cut(indexNValMuHits_, int()) || ignoreCut(indexNValMuHits_))
0469       passCut(ret, indexNValMuHits_);
0470     if (hcalVeto < cut(indexHCalVeto_, double()) || ignoreCut(indexHCalVeto_))
0471       passCut(ret, indexHCalVeto_);
0472     if (ecalVeto < cut(indexECalVeto_, double()) || ignoreCut(indexECalVeto_))
0473       passCut(ret, indexECalVeto_);
0474     if (relIso < cut(indexRelIso_, double()) || ignoreCut(indexRelIso_))
0475       passCut(ret, indexRelIso_);
0476     if (fabs(z_mu - zvtx) < cut(indexLepZ_, double()) || ignoreCut(indexLepZ_))
0477       passCut(ret, indexLepZ_);
0478     if (nPixelHits > cut(indexPixHits_, int()) || ignoreCut(indexPixHits_))
0479       passCut(ret, indexPixHits_);
0480     if (nMatchedStations > cut(indexStations_, int()) || ignoreCut(indexStations_))
0481       passCut(ret, indexStations_);
0482 
0483     setIgnored(ret);
0484 
0485     return (bool)ret;
0486   }
0487 
0488   // cuts based on top group L+J synchronization exercise
0489   // this is a copy of fall 10 cuts
0490   // with a hack to include a double-sided reliso cut
0491 
0492   bool kitQCDCuts(const pat::Muon& muon, edm::EventBase const& event, pat::strbitset& ret) {
0493     ret.set(false);
0494 
0495     double norm_chi2 = muon.normChi2();
0496     double corr_d0 = muon.dB();
0497     double corr_ed0 = muon.edB();
0498     double corr_sd0 = (corr_ed0 > 0.000000001) ? corr_d0 / corr_ed0 : 999.0;
0499 
0500     // Get the PV for the muon z requirement
0501     edm::Handle<std::vector<reco::Vertex> > pvtxHandle_;
0502     event.getByLabel(pvSrc_, pvtxHandle_);
0503 
0504     double zvtx = -999;
0505     if (!pvtxHandle_->empty()) {
0506       zvtx = pvtxHandle_->at(0).z();
0507     } else {
0508       throw cms::Exception("InvalidInput")
0509           << " There needs to be at least one primary vertex in the event." << std::endl;
0510     }
0511 
0512     //If required, recalculate the impact parameter using the beam spot
0513     if (recalcDBFromBSp_) {
0514       //Get the beam spot
0515       reco::TrackBase::Point beamPoint(0, 0, 0);
0516       reco::BeamSpot beamSpot;
0517       edm::Handle<reco::BeamSpot> beamSpotHandle;
0518       event.getByLabel(beamLineSrc_, beamSpotHandle);
0519 
0520       if (beamSpotHandle.isValid()) {
0521         beamSpot = *beamSpotHandle;
0522       } else {
0523         edm::LogError("DataNotAvailable")
0524             << "No beam spot available from EventSetup, not adding high level selection \n";
0525       }
0526       beamPoint = reco::TrackBase::Point(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
0527 
0528       //Use the beamspot to correct the impact parameter and uncertainty
0529       reco::TrackRef innerTrack = muon.innerTrack();
0530       if (innerTrack.isNonnull() && innerTrack.isAvailable()) {
0531         corr_d0 = -1.0 * innerTrack->dxy(beamPoint);
0532         corr_ed0 =
0533             sqrt(innerTrack->d0Error() * innerTrack->d0Error() + 0.5 * beamSpot.BeamWidthX() * beamSpot.BeamWidthX() +
0534                  0.5 * beamSpot.BeamWidthY() * beamSpot.BeamWidthY());
0535         corr_sd0 = (corr_ed0 > 0.000000001) ? corr_d0 / corr_ed0 : 999.0;
0536 
0537       } else {
0538         corr_d0 = 999.;
0539         corr_ed0 = 999.;
0540       }
0541     }
0542 
0543     int nhits = static_cast<int>(muon.numberOfValidHits());
0544     int nValidMuonHits = static_cast<int>(muon.globalTrack()->hitPattern().numberOfValidMuonHits());
0545 
0546     double ecalVeto = muon.isolationR03().emVetoEt;
0547     double hcalVeto = muon.isolationR03().hadVetoEt;
0548 
0549     double hcalIso = muon.hcalIso();
0550     double ecalIso = muon.ecalIso();
0551     double trkIso = muon.trackIso();
0552     double pt = muon.pt();
0553 
0554     double relIso = (ecalIso + hcalIso + trkIso) / pt;
0555 
0556     double z_mu = muon.vertex().z();
0557 
0558     int nPixelHits = muon.innerTrack()->hitPattern().pixelLayersWithMeasurement();
0559 
0560     int nMatchedStations = muon.numberOfMatches();
0561 
0562     if (norm_chi2 < cut(indexChi2_, double()) || ignoreCut(indexChi2_))
0563       passCut(ret, indexChi2_);
0564     if (fabs(corr_d0) < cut(indexD0_, double()) || ignoreCut(indexD0_))
0565       passCut(ret, indexD0_);
0566     if (fabs(corr_ed0) < cut(indexED0_, double()) || ignoreCut(indexED0_))
0567       passCut(ret, indexED0_);
0568     if (fabs(corr_sd0) < cut(indexSD0_, double()) || ignoreCut(indexSD0_))
0569       passCut(ret, indexSD0_);
0570     if (nhits >= cut(indexNHits_, int()) || ignoreCut(indexNHits_))
0571       passCut(ret, indexNHits_);
0572     if (nValidMuonHits > cut(indexNValMuHits_, int()) || ignoreCut(indexNValMuHits_))
0573       passCut(ret, indexNValMuHits_);
0574     if (hcalVeto < cut(indexHCalVeto_, double()) || ignoreCut(indexHCalVeto_))
0575       passCut(ret, indexHCalVeto_);
0576     if (ecalVeto < cut(indexECalVeto_, double()) || ignoreCut(indexECalVeto_))
0577       passCut(ret, indexECalVeto_);
0578     if (fabs(z_mu - zvtx) < cut(indexLepZ_, double()) || ignoreCut(indexLepZ_))
0579       passCut(ret, indexLepZ_);
0580     if (nPixelHits > cut(indexPixHits_, int()) || ignoreCut(indexPixHits_))
0581       passCut(ret, indexPixHits_);
0582     if (nMatchedStations > cut(indexStations_, int()) || ignoreCut(indexStations_))
0583       passCut(ret, indexStations_);
0584 
0585     ////////////////////////////////////////////////////////////////
0586     //
0587     // JMS Dec 13 2010
0588     // HACK
0589     // Need double-sided relIso cut to implement data-driven QCD
0590     //
0591     //
0592     ///////////////////////////////////////////////////////////////
0593     if (((relIso > 0.2) && (relIso < 0.75)) || ignoreCut(indexRelIso_))
0594       passCut(ret, indexRelIso_);
0595 
0596     setIgnored(ret);
0597 
0598     return (bool)ret;
0599   }
0600 
0601 private:  // member variables
0602   Version_t version_;
0603   bool recalcDBFromBSp_;
0604   edm::InputTag beamLineSrc_;
0605 #ifndef __GCCXML__
0606   edm::EDGetTokenT<reco::BeamSpot> beamLineSrcToken_;
0607 #endif
0608   edm::InputTag pvSrc_;
0609 #ifndef __GCCXML__
0610   edm::EDGetTokenT<std::vector<reco::Vertex> > pvSrcToken_;
0611 #endif
0612 
0613   index_type indexChi2_;
0614   index_type indexD0_;
0615   index_type indexED0_;
0616   index_type indexSD0_;
0617   index_type indexNHits_;
0618   index_type indexNValMuHits_;
0619   index_type indexECalVeto_;
0620   index_type indexHCalVeto_;
0621   index_type indexRelIso_;
0622   index_type indexLepZ_;
0623   index_type indexPixHits_;
0624   index_type indexStations_;
0625 };
0626 
0627 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!