Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-04-11 22:25:11 

0001 // -*- C++ -*-
0002 //
0003 // Package:    MuMu/DisappearingMuonsSkimming
0004 // Class:      DisappearingMuonsSkimming
0005 //
0006 /**\class DisappearingMuonsSkimming DisappearingMuonsSkimming_AOD.cc MuMu/DisappearingMuonsSkimming/plugins/DisappearingMuonsSkimming_AOD.cc
0007 
0008  Description: [one line class summary]
0009 
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Michael Revering
0015 //         Created:  Mon, 18 Jun 2018 21:22:23 GMT
0016 //
0017 //
0018 // system include files
0019 #include <memory>
0020 #include "Math/VectorUtil.h"
0021 
0022 // user include files
0023 #include "Configuration/Skimming/interface/DisappearingMuonsSkimming.h"
0024 #include "DataFormats/Math/interface/deltaR.h"
0025 #include "FWCore/Common/interface/TriggerNames.h"
0026 #include "FWCore/Framework/interface/MakerMacros.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/Frameworkfwd.h"
0029 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0030 #include "FWCore/Utilities/interface/InputTag.h"
0031 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0032 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0033 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0034 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0035 
0036 DisappearingMuonsSkimming::DisappearingMuonsSkimming(const edm::ParameterSet& iConfig)
0037     : recoMuonToken_(consumes<reco::MuonCollection>(iConfig.getParameter<edm::InputTag>("recoMuons"))),
0038       standaloneMuonToken_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("StandaloneTracks"))),
0039       trackCollectionToken_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"))),
0040       primaryVerticesToken_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("primaryVertices"))),
0041       reducedEndcapRecHitCollectionToken_(
0042           consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("EERecHits"))),
0043       reducedBarrelRecHitCollectionToken_(
0044           consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("EBRecHits"))),
0045       trigResultsToken_(consumes<edm::TriggerResults>(iConfig.getParameter<edm::InputTag>("TriggerResultsTag"))),
0046       transientTrackToken_(
0047           esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"))),
0048       geometryToken_(esConsumes<CaloGeometry, CaloGeometryRecord>(edm::ESInputTag{})),
0049       muonPathsToPass_(iConfig.getParameter<std::vector<std::string>>("muonPathsToPass")),
0050       minMuPt_(iConfig.getParameter<double>("minMuPt")),
0051       maxMuEta_(iConfig.getParameter<double>("maxMuEta")),
0052       minTrackEta_(iConfig.getParameter<double>("minTrackEta")),
0053       maxTrackEta_(iConfig.getParameter<double>("maxTrackEta")),
0054       minTrackPt_(iConfig.getParameter<double>("minTrackPt")),
0055       maxTransDCA_(iConfig.getParameter<double>("maxTransDCA")),
0056       maxLongDCA_(iConfig.getParameter<double>("maxLongDCA")),
0057       maxVtxChi_(iConfig.getParameter<double>("maxVtxChi")),
0058       minInvMass_(iConfig.getParameter<double>("minInvMass")),
0059       maxInvMass_(iConfig.getParameter<double>("maxInvMass")),
0060       trackIsoConesize_(iConfig.getParameter<double>("trackIsoConesize")),
0061       trackIsoInnerCone_(iConfig.getParameter<double>("trackIsoInnerCone")),
0062       ecalIsoConesize_(iConfig.getParameter<double>("ecalIsoConesize")),
0063       minEcalHitE_(iConfig.getParameter<double>("minEcalHitE")),
0064       maxTrackIso_(iConfig.getParameter<double>("maxTrackIso")),
0065       maxEcalIso_(iConfig.getParameter<double>("maxEcalIso")),
0066       minSigInvMass_(iConfig.getParameter<double>("minSigInvMass")),
0067       maxSigInvMass_(iConfig.getParameter<double>("maxSigInvMass")),
0068       minStandaloneDr_(iConfig.getParameter<double>("minStandaloneDr")),
0069       maxStandaloneDE_(iConfig.getParameter<double>("maxStandaloneDE")),
0070       keepOffPeak_(iConfig.getParameter<bool>("keepOffPeak")),
0071       keepSameSign_(iConfig.getParameter<bool>("keepSameSign")),
0072       keepTotalRegion_(iConfig.getParameter<bool>("keepTotalRegion")),
0073       keepPartialRegion_(iConfig.getParameter<bool>("keepPartialRegion")) {}
0074 
0075 //
0076 // member functions
0077 //
0078 
0079 // ------------ method called for each event  ------------
0080 bool DisappearingMuonsSkimming::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0081   using namespace edm;
0082   using namespace std;
0083   using namespace reco;
0084 
0085   bool totalRegion = false;
0086   bool sameSign = false;
0087   bool offPeak = false;
0088   bool partialRegion = false;
0089 
0090   const auto& staTracks = iEvent.get(standaloneMuonToken_);
0091   const auto& vertices = iEvent.get(primaryVerticesToken_);
0092   const auto& recoMuons = iEvent.get(recoMuonToken_);
0093   const auto& thePATTracks = iEvent.get(trackCollectionToken_);
0094   const auto& triggerResults = iEvent.get(trigResultsToken_);
0095 
0096   // this wraps tracks with additional methods that are used in vertex-calculation
0097   const TransientTrackBuilder* transientTrackBuilder = &iSetup.getData(transientTrackToken_);
0098 
0099   if (!passTriggers(iEvent, triggerResults, muonPathsToPass_))
0100     return false;
0101 
0102   int nMuonTrackCand = 0;
0103   float MuonTrackMass = 0.;
0104 
0105   //Looping over the reconstructed Muons
0106   for (const auto& iMuon : recoMuons) {
0107     if (!(iMuon.isPFMuon() && iMuon.isGlobalMuon()))
0108       continue;
0109     if (!(iMuon.passed(reco::Muon::CutBasedIdTight)))
0110       continue;
0111     if (!(iMuon.passed(reco::Muon::PFIsoTight)))
0112       continue;
0113     if (iMuon.pt() < minMuPt_ || std::abs(iMuon.eta()) > maxMuEta_)
0114       continue;
0115 
0116     //Looping over tracks for any good muon
0117     int indx(0);
0118     for (const auto& iTrack : thePATTracks) {
0119       reco::TrackRef trackRef = reco::TrackRef(&thePATTracks, indx);
0120       indx++;
0121 
0122       if (!iTrack.quality(reco::Track::qualityByName("highPurity")))
0123         continue;
0124       if (std::abs(iTrack.eta()) > maxTrackEta_ || std::abs(iTrack.eta()) < minTrackEta_)
0125         continue;
0126       if (iTrack.pt() < minTrackPt_)
0127         continue;
0128       //Check if the track belongs to a primary vertex for isolation calculation
0129 
0130       unsigned int vtxIndex;
0131       unsigned int tkIndex;
0132       bool foundtrack = this->findTrackInVertices(trackRef, vertices, vtxIndex, tkIndex);
0133 
0134       if (!foundtrack)
0135         continue;
0136 
0137       reco::VertexRef vtx(&vertices, vtxIndex);
0138       GlobalPoint tkVtx = GlobalPoint(vtx->x(), vtx->y(), vtx->z());
0139 
0140       reco::TransientTrack tk = transientTrackBuilder->build(iTrack);
0141       TrajectoryStateClosestToPoint traj = tk.trajectoryStateClosestToPoint(tkVtx);
0142       double transDCA = traj.perigeeParameters().transverseImpactParameter();
0143       double longDCA = traj.perigeeParameters().longitudinalImpactParameter();
0144       if (std::abs(longDCA) > maxLongDCA_)
0145         continue;
0146       if (std::abs(transDCA) > maxTransDCA_)
0147         continue;
0148       // make a pair of TransientTracks to feed the vertexer
0149       std::vector<reco::TransientTrack> tracksToVertex;
0150       tracksToVertex.push_back(tk);
0151       tracksToVertex.push_back(transientTrackBuilder->build(iMuon.globalTrack()));
0152       // try to fit these two tracks to a common vertex
0153       KalmanVertexFitter vertexFitter;
0154       CachingVertex<5> fittedVertex = vertexFitter.vertex(tracksToVertex);
0155       double vtxChi = 0;
0156       // some poor fits will simply fail to find a common vertex
0157       if (fittedVertex.isValid() && fittedVertex.totalChiSquared() >= 0. && fittedVertex.degreesOfFreedom() > 0) {
0158         // others we can exclude by their poor fit
0159         vtxChi = fittedVertex.totalChiSquared() / fittedVertex.degreesOfFreedom();
0160 
0161         if (vtxChi < maxVtxChi_) {
0162           // important! evaluate momentum vectors AT THE VERTEX
0163           TrajectoryStateClosestToPoint one_TSCP =
0164               tracksToVertex[0].trajectoryStateClosestToPoint(fittedVertex.position());
0165           TrajectoryStateClosestToPoint two_TSCP =
0166               tracksToVertex[1].trajectoryStateClosestToPoint(fittedVertex.position());
0167           GlobalVector one_momentum = one_TSCP.momentum();
0168           GlobalVector two_momentum = two_TSCP.momentum();
0169 
0170           double total_energy = sqrt(one_momentum.mag2() + 0.106 * 0.106) + sqrt(two_momentum.mag2() + 0.106 * 0.106);
0171           double total_px = one_momentum.x() + two_momentum.x();
0172           double total_py = one_momentum.y() + two_momentum.y();
0173           double total_pz = one_momentum.z() + two_momentum.z();
0174           MuonTrackMass = sqrt(pow(total_energy, 2) - pow(total_px, 2) - pow(total_py, 2) - pow(total_pz, 2));
0175         } else {
0176           continue;
0177         }
0178       } else {
0179         continue;
0180       }
0181       if (MuonTrackMass < minInvMass_ || MuonTrackMass > maxInvMass_)
0182         continue;
0183 
0184       double trackIso = getTrackIsolation(trackRef, vertices);
0185       //Track iso returns -1 when it fails to find the vertex containing the track (already checked in track selection, but might as well make sure)
0186       if (trackIso < 0)
0187         continue;
0188       double ecalIso = getECALIsolation(iEvent, iSetup, transientTrackBuilder->build(iTrack));
0189       if (trackIso > maxTrackIso_ || ecalIso > maxEcalIso_)
0190         continue;
0191 
0192       //A good tag/probe pair has been selected, now check for control or signal regions
0193       if (iMuon.charge() == iTrack.charge()) {
0194         sameSign = true;
0195       }
0196 
0197       //If not same sign CR, need to check standalone muons for signal regions
0198       double staMinDr2 = 1000;
0199       double staMinDEoverE = -10;
0200       if (!staTracks.empty()) {
0201         for (const auto& staTrack : staTracks) {
0202           reco::TransientTrack track = transientTrackBuilder->build(staTrack);
0203           double dR2 = deltaR2(track.impactPointTSCP().momentum().eta(),
0204                                track.impactPointTSCP().momentum().phi(),
0205                                iTrack.eta(),
0206                                iTrack.phi());
0207           double staDE = (std::sqrt(track.impactPointTSCP().momentum().mag2()) - iTrack.p()) / iTrack.p();
0208           if (dR2 < staMinDr2) {
0209             staMinDr2 = dR2;
0210           }
0211           //Pick the largest standalone muon within the cone
0212           if (dR2 < minStandaloneDr_ * minStandaloneDr_ && staDE > staMinDEoverE) {
0213             staMinDEoverE = staDE;
0214           }
0215         }
0216       }
0217       if (staMinDr2 > minStandaloneDr_ * minStandaloneDr_) {
0218         if (MuonTrackMass < minSigInvMass_ || MuonTrackMass > maxSigInvMass_) {
0219           offPeak = true;
0220         } else {
0221           totalRegion = true;
0222         }
0223       } else {
0224         if (staMinDEoverE < maxStandaloneDE_) {
0225           if (MuonTrackMass < minSigInvMass_ || MuonTrackMass > maxSigInvMass_) {
0226             offPeak = true;
0227           } else {
0228             partialRegion = true;
0229           }
0230         }
0231       }
0232       nMuonTrackCand++;
0233     }
0234   }
0235 
0236   if (nMuonTrackCand == 0)
0237     return false;
0238   bool passes = false;
0239   //Pass all same sign CR events
0240   if (sameSign && keepSameSign_) {
0241     passes = true;
0242   }
0243   //Pass all total disappearance events
0244   else if (totalRegion && keepTotalRegion_) {
0245     passes = true;
0246   }
0247   //Pass all partial disappearkance off-peak events
0248   else if (offPeak && keepOffPeak_) {
0249     passes = true;
0250   }
0251   //Pass partial region events that pass minimum standalone muon DE/E.
0252   else if (partialRegion && keepPartialRegion_) {
0253     passes = true;
0254   }
0255   return passes;
0256 }
0257 
0258 bool DisappearingMuonsSkimming::passTriggers(const edm::Event& iEvent,
0259                                              const edm::TriggerResults& triggerResults,
0260                                              const std::vector<std::string>& m_muonPathsToPass) {
0261   bool passTriggers = false;
0262   const edm::TriggerNames& trigNames = iEvent.triggerNames(triggerResults);
0263   for (size_t i = 0; i < trigNames.size(); ++i) {
0264     const std::string& name = trigNames.triggerName(i);
0265     for (auto& pathName : m_muonPathsToPass) {
0266       if ((name.find(pathName) != std::string::npos)) {
0267         if (triggerResults.accept(i)) {
0268           passTriggers = true;
0269           break;
0270         }
0271       }
0272     }
0273   }
0274   return passTriggers;
0275 }
0276 
0277 bool DisappearingMuonsSkimming::findTrackInVertices(const reco::TrackRef& tkToMatch,
0278                                                     const reco::VertexCollection& vertices,
0279                                                     unsigned int& vtxIndex,
0280                                                     unsigned int& trackIndex) {
0281   // initialize indices
0282   vtxIndex = -1;
0283   trackIndex = -1;
0284 
0285   bool foundtrack{false};
0286   unsigned int idx = 0;
0287   for (const auto& vtx : vertices) {
0288     if (!vtx.isValid()) {
0289       idx++;
0290       continue;
0291     }
0292 
0293     const auto& thePVtracks{vtx.tracks()};
0294     std::vector<unsigned int> thePVkeys;
0295     thePVkeys.reserve(thePVtracks.size());
0296     for (const auto& tv : thePVtracks) {
0297       thePVkeys.push_back(tv.key());
0298     }
0299 
0300     auto result = std::find_if(
0301         thePVkeys.begin(), thePVkeys.end(), [tkToMatch](const auto& tkRef) { return tkRef == tkToMatch.key(); });
0302     if (result != thePVkeys.end()) {
0303       foundtrack = true;
0304       trackIndex = std::distance(thePVkeys.begin(), result);
0305       vtxIndex = idx;
0306     }
0307     if (foundtrack)
0308       break;
0309 
0310     idx++;
0311   }
0312   return foundtrack;
0313 }
0314 
0315 double DisappearingMuonsSkimming::getTrackIsolation(const reco::TrackRef& tkToMatch,
0316                                                     const reco::VertexCollection& vertices) {
0317   unsigned int vtxindex;
0318   unsigned int trackIndex;
0319   bool foundtrack = this->findTrackInVertices(tkToMatch, vertices, vtxindex, trackIndex);
0320 
0321   LogDebug("DisappearingMuonsSkimming") << "getTrackIsolation vtx Index: " << vtxindex
0322                                         << " track Index: " << trackIndex;
0323 
0324   if (!foundtrack) {
0325     return -1.;
0326   }
0327 
0328   reco::VertexRef primaryVtx(&vertices, vtxindex);
0329 
0330   double Isolation = 0;
0331   for (unsigned int i = 0; i < primaryVtx->tracksSize(); i++) {
0332     if (i == trackIndex)
0333       continue;
0334     reco::TrackBaseRef secondarytrack = primaryVtx->trackRefAt(i);
0335     if (deltaR2(tkToMatch.get()->eta(), tkToMatch.get()->phi(), secondarytrack->eta(), secondarytrack->phi()) >
0336             trackIsoConesize_ * trackIsoConesize_ ||
0337         deltaR2(tkToMatch.get()->eta(), tkToMatch.get()->phi(), secondarytrack->eta(), secondarytrack->phi()) <
0338             trackIsoInnerCone_ * trackIsoInnerCone_)
0339       continue;
0340     Isolation += secondarytrack->pt();
0341   }
0342 
0343   return Isolation / tkToMatch.get()->pt();
0344 }
0345 
0346 double DisappearingMuonsSkimming::getECALIsolation(const edm::Event& iEvent,
0347                                                    const edm::EventSetup& iSetup,
0348                                                    const reco::TransientTrack& track) {
0349   edm::Handle<EcalRecHitCollection> rechitsEE;
0350   iEvent.getByToken(reducedEndcapRecHitCollectionToken_, rechitsEE);
0351 
0352   edm::Handle<EcalRecHitCollection> rechitsEB;
0353   iEvent.getByToken(reducedBarrelRecHitCollectionToken_, rechitsEB);
0354 
0355   const CaloGeometry& caloGeom = iSetup.getData(geometryToken_);
0356   TrajectoryStateClosestToPoint t0 = track.impactPointTSCP();
0357   double eDR = 0;
0358 
0359   for (EcalRecHitCollection::const_iterator hit = rechitsEE->begin(); hit != rechitsEE->end(); hit++) {
0360     const DetId id = (*hit).detid();
0361     const GlobalPoint hitPos = caloGeom.getSubdetectorGeometry(id)->getGeometry(id)->getPosition();
0362     //Check if hit and track trajectory ar in the same endcap (transient track projects both ways)
0363     if ((hitPos.eta() * t0.momentum().eta()) < 0) {
0364       continue;
0365     }
0366     TrajectoryStateClosestToPoint traj = track.trajectoryStateClosestToPoint(hitPos);
0367     math::XYZVector idPositionRoot(hitPos.x(), hitPos.y(), hitPos.z());
0368     math::XYZVector trajRoot(traj.position().x(), traj.position().y(), traj.position().z());
0369     if (deltaR2(idPositionRoot, trajRoot) < ecalIsoConesize_ * ecalIsoConesize_ && (*hit).energy() > minEcalHitE_) {
0370       eDR += (*hit).energy();
0371     }
0372   }
0373   for (EcalRecHitCollection::const_iterator hit = rechitsEB->begin(); hit != rechitsEB->end(); hit++) {
0374     const DetId id = (*hit).detid();
0375     const GlobalPoint hitPos = caloGeom.getSubdetectorGeometry(id)->getGeometry(id)->getPosition();
0376     if ((hitPos.eta() * t0.momentum().eta()) < 0) {
0377       continue;
0378     }
0379     TrajectoryStateClosestToPoint traj = track.trajectoryStateClosestToPoint(hitPos);
0380     math::XYZVector idPositionRoot(hitPos.x(), hitPos.y(), hitPos.z());
0381     math::XYZVector trajRoot(traj.position().x(), traj.position().y(), traj.position().z());
0382     if (deltaR2(idPositionRoot, trajRoot) < ecalIsoConesize_ * ecalIsoConesize_ && (*hit).energy() > minEcalHitE_) {
0383       eDR += (*hit).energy();
0384     }
0385   }
0386 
0387   return eDR;
0388 }
0389 
0390 void DisappearingMuonsSkimming::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0391   edm::ParameterSetDescription desc;
0392 
0393   desc.add<edm::InputTag>("recoMuons", edm::InputTag("muons"));
0394   desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0395   desc.add<edm::InputTag>("StandaloneTracks", edm::InputTag("standAloneMuons"));
0396   desc.add<edm::InputTag>("primaryVertices", edm::InputTag("offlinePrimaryVertices"));
0397   desc.add<edm::InputTag>("EERecHits", edm::InputTag("reducedEcalRecHitsEE"));
0398   desc.add<edm::InputTag>("EBRecHits", edm::InputTag("reducedEcalRecHitsEB"));
0399   desc.add<edm::InputTag>("TriggerResultsTag", edm::InputTag("TriggerResults", "", "HLT"));
0400   desc.add<std::vector<std::string>>("muonPathsToPass",
0401                                      {
0402                                          "HLT_IsoMu24_v",
0403                                          "HLT_IsoMu27_v",
0404                                      });
0405   desc.add<double>("minMuPt", 26);
0406   desc.add<double>("maxMuEta", 2.4);
0407   desc.add<double>("minTrackEta", 0);
0408   desc.add<double>("maxTrackEta", 2.4);
0409   desc.add<double>("minTrackPt", 20);
0410   desc.add<double>("maxTransDCA", 0.005);
0411   desc.add<double>("maxLongDCA", 0.05);
0412   desc.add<double>("maxVtxChi", 3.0);
0413   desc.add<double>("minInvMass", 50);
0414   desc.add<double>("maxInvMass", 150);
0415   desc.add<double>("trackIsoConesize", 0.3);
0416   desc.add<double>("trackIsoInnerCone", 0.01);
0417   desc.add<double>("ecalIsoConesize", 0.4);
0418   desc.add<double>("minEcalHitE", 0.3);
0419   desc.add<double>("maxTrackIso", 0.05);
0420   desc.add<double>("maxEcalIso", 10);
0421   desc.add<double>("minSigInvMass", 76);
0422   desc.add<double>("maxSigInvMass", 106);
0423   desc.add<double>("minStandaloneDr", 1.0);
0424   desc.add<double>("maxStandaloneDE", -0.5);
0425   desc.add<bool>("keepOffPeak", true);
0426   desc.add<bool>("keepSameSign", true);
0427   desc.add<bool>("keepTotalRegion", true);
0428   desc.add<bool>("keepPartialRegion", true);
0429   descriptions.addWithDefaultLabel(desc);
0430 }
0431 
0432 //define this as a plug-in
0433 DEFINE_FWK_MODULE(DisappearingMuonsSkimming);
[ 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!