Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-09-08 03:21:53

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  */
0005 
0006 #include <fmt/printf.h>
0007 
0008 #include "DQMOffline/Alignment/interface/DiMuonVertexMonitor.h"
0009 #include "DQMServices/Core/interface/DQMStore.h"
0010 #include "DataFormats/Math/interface/deltaR.h"
0011 #include "DataFormats/TrackReco/interface/Track.h"
0012 #include "DataFormats/TrackReco/interface/TrackBase.h"
0013 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FWCore/Utilities/interface/InputTag.h"
0016 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0017 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0018 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0019 #include "RecoVertex/VertexTools/interface/VertexDistanceXY.h"
0020 #include "TrackingTools/IPTools/interface/IPTools.h"
0021 
0022 #include "TLorentzVector.h"
0023 
0024 namespace {
0025   constexpr float cmToum = 10e4;
0026   constexpr float mumass2 = 0.105658367 * 0.105658367;  //mu mass squared (GeV^2/c^4)
0027 }  // namespace
0028 
0029 DiMuonVertexMonitor::DiMuonVertexMonitor(const edm::ParameterSet& iConfig)
0030     : ttbESToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0031       tracksToken_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("muonTracks"))),
0032       vertexToken_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))),
0033       motherName_(iConfig.getParameter<std::string>("decayMotherName")),
0034       MEFolderName_(iConfig.getParameter<std::string>("FolderName")),
0035       useClosestVertex_(iConfig.getParameter<bool>("useClosestVertex")),
0036       maxSVdist_(iConfig.getParameter<double>("maxSVdist")) {
0037   if (motherName_.find('Z') != std::string::npos) {
0038     massLimits_ = std::make_pair(50., 120);
0039   } else if (motherName_.find("J/#psi") != std::string::npos) {
0040     massLimits_ = std::make_pair(2.7, 3.4);
0041   } else if (motherName_.find("#Upsilon") != std::string::npos) {
0042     massLimits_ = std::make_pair(8.9, 9.9);
0043   } else {
0044     edm::LogError("DiMuonVertexMonitor") << " unrecognized decay mother particle: " << motherName_
0045                                          << " setting the default for the Z->mm (50.,120.)" << std::endl;
0046     massLimits_ = std::make_pair(50., 120);
0047   }
0048 }
0049 
0050 void DiMuonVertexMonitor::bookHistograms(DQMStore::IBooker& iBooker, edm::Run const&, edm::EventSetup const&) {
0051   iBooker.setCurrentFolder(MEFolderName_ + "/DiMuonVertexMonitor");
0052 
0053   // clang-format off
0054   std::string ts = fmt::sprintf(";%s vertex probability;N(#mu#mu pairs)", motherName_);
0055   std::string ps = "N(#mu#mu pairs)";
0056   hSVProb_ = iBooker.book1D("VtxProb", ts, 100, 0., 1.);
0057 
0058   ts = fmt::sprintf("#chi^{2} of the %s vertex; #chi^{2} of the %s vertex; %s", motherName_, motherName_, ps);
0059   hSVChi2_ = iBooker.book1D("VtxChi2", ts, 200, 0., 200.);
0060 
0061   ts = fmt::sprintf("#chi^{2}/ndf of the %s vertex; #chi^{2}/ndf of %s vertex; %s", motherName_, motherName_, ps);
0062   hSVNormChi2_ = iBooker.book1D("VtxNormChi2", ts, 100, 0., 20.);
0063 
0064   std::string histTit = motherName_ + " #rightarrow #mu^{+}#mu^{-}";
0065   ts = fmt::sprintf("%s;PV- %sV xy distance [#mum];%s", histTit, motherName_, ps);
0066   hSVDist_ = iBooker.book1D("VtxDist", ts, 100, 0., 300.);
0067 
0068   ts = fmt::sprintf("%s;PV-%sV xy distance error [#mum];%s", histTit, motherName_, ps);
0069   hSVDistErr_ = iBooker.book1D("VtxDistErr", ts, 100, 0., 1000.);
0070 
0071   ts = fmt::sprintf("%s;PV-%sV xy distance signficance;%s", histTit, motherName_, ps);
0072   hSVDistSig_ = iBooker.book1D("VtxDistSig", ts, 100, 0., 5.);
0073 
0074   ts = fmt::sprintf("compatibility of %s vertex; compatibility of %s vertex; %s", motherName_, motherName_, ps);
0075   hSVCompatibility_ = iBooker.book1D("VtxCompatibility", ts, 100, 0., 100.);
0076 
0077   ts = fmt::sprintf("%s;PV-%sV 3D distance [#mum];%s", histTit, motherName_, ps);
0078   hSVDist3D_ = iBooker.book1D("VtxDist3D", ts, 100, 0., 300.);
0079 
0080   ts = fmt::sprintf("%s;PV-%sV 3D distance error [#mum];%s", histTit, motherName_, ps);
0081   hSVDist3DErr_ = iBooker.book1D("VtxDist3DErr", ts, 100, 0., 1000.);
0082 
0083   ts = fmt::sprintf("%s;PV-%sV 3D distance signficance;%s", histTit, motherName_, ps);
0084   hSVDist3DSig_ = iBooker.book1D("VtxDist3DSig", ts, 100, 0., 5.);
0085 
0086   ts = fmt::sprintf("3D compatibility of %s vertex;3D compatibility of %s vertex; %s", motherName_, motherName_, ps);
0087   hSVCompatibility3D_ = iBooker.book1D("VtxCompatibility3D", ts, 100, 0., 100.);
0088 
0089   hInvMass_ = iBooker.book1D("InvMass", fmt::sprintf("%s;M(#mu,#mu) [GeV];%s", histTit, ps), 70., massLimits_.first, massLimits_.second);
0090   hCosPhi_ = iBooker.book1D("CosPhi", fmt::sprintf("%s;cos(#phi_{xy});%s", histTit, ps), 50, -1., 1.);
0091   hCosPhi3D_ = iBooker.book1D("CosPhi3D", fmt::sprintf("%s;cos(#phi_{3D});%s", histTit, ps), 50, -1., 1.);
0092   hCosPhiInv_ = iBooker.book1D("CosPhiInv", fmt::sprintf("%s;inverted cos(#phi_{xy});%s", histTit, ps), 50, -1., 1.);
0093   hCosPhiInv3D_ = iBooker.book1D("CosPhiInv3D", fmt::sprintf("%s;inverted cos(#phi_{3D});%s", histTit, ps), 50, -1., 1.);
0094 
0095   hdxy_ = iBooker.book1D("dxy", fmt::sprintf("%s;muon track d_{xy}(PV) [#mum];muon tracks", histTit), 150, -300, 300);
0096   hdz_ = iBooker.book1D("dz", fmt::sprintf("%s;muon track d_{z}(PV) [#mum];muon tracks", histTit), 150, -300, 300);
0097   hdxyErr_ = iBooker.book1D("dxyErr", fmt::sprintf("%s;muon track err_{dxy} [#mum];muon tracks", histTit), 250, 0., 500.);
0098   hdzErr_ = iBooker.book1D("dzErr", fmt::sprintf("%s;muon track err_{dz} [#mum];muon tracks", histTit), 250, 0., 500.);
0099   hIP2d_ = iBooker.book1D("IP2d", fmt::sprintf("%s;muon track IP_{2D} [#mum];muon tracks", histTit), 150, -300, 300);
0100   hIP3d_ = iBooker.book1D("IP3d", fmt::sprintf("%s;muon track IP_{3D} [#mum];muon tracks", histTit), 150, -300, 300);
0101   hIP2dsig_ = iBooker.book1D("IP2Dsig", fmt::sprintf("%s;muon track IP_{2D} significance;muon tracks", histTit), 100, 0., 5.);
0102   hIP3dsig_ = iBooker.book1D("IP3Dsig", fmt::sprintf("%s;muon track IP_{3D} significance;muon tracks", histTit), 100, 0., 5.);
0103   // clang-format on
0104 }
0105 
0106 void DiMuonVertexMonitor::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0107   std::vector<const reco::Track*> myTracks;
0108   const auto trackHandle = iEvent.getHandle(tracksToken_);
0109   if (!trackHandle.isValid()) {
0110     edm::LogError("DiMuonVertexMonitor") << "invalid track collection encountered!";
0111     return;
0112   }
0113 
0114   for (const auto& muonTrk : *trackHandle) {
0115     myTracks.emplace_back(&muonTrk);
0116   }
0117 
0118   const TransientTrackBuilder* theB = &iSetup.getData(ttbESToken_);
0119   TransientVertex mumuTransientVtx;
0120   std::vector<reco::TransientTrack> tks;
0121 
0122   if (myTracks.size() != 2) {
0123     edm::LogWarning("DiMuonVertexMonitor") << "There are not enough tracks to monitor!";
0124     return;
0125   }
0126 
0127   const auto& t1 = myTracks[1]->momentum();
0128   const auto& t0 = myTracks[0]->momentum();
0129   const auto& ditrack = t1 + t0;
0130 
0131   const auto& tplus = myTracks[0]->charge() > 0 ? myTracks[0] : myTracks[1];
0132   const auto& tminus = myTracks[0]->charge() < 0 ? myTracks[0] : myTracks[1];
0133 
0134   TLorentzVector p4_tplus(tplus->px(), tplus->py(), tplus->pz(), sqrt((tplus->p() * tplus->p()) + mumass2));
0135   TLorentzVector p4_tminus(tminus->px(), tminus->py(), tminus->pz(), sqrt((tminus->p() * tminus->p()) + mumass2));
0136 
0137   const auto& Zp4 = p4_tplus + p4_tminus;
0138   float track_invMass = Zp4.M();
0139   hInvMass_->Fill(track_invMass);
0140 
0141   // creat the pair of TLorentVectors used to make the plos
0142   std::pair<TLorentzVector, TLorentzVector> tktk_p4 = std::make_pair(p4_tplus, p4_tminus);
0143 
0144   math::XYZPoint ZpT(ditrack.x(), ditrack.y(), 0);
0145   math::XYZPoint Zp(ditrack.x(), ditrack.y(), ditrack.z());
0146 
0147   for (const auto& track : myTracks) {
0148     reco::TransientTrack trajectory = theB->build(track);
0149     tks.push_back(trajectory);
0150   }
0151 
0152   KalmanVertexFitter kalman(true);
0153   mumuTransientVtx = kalman.vertex(tks);
0154 
0155   double SVProb = TMath::Prob(mumuTransientVtx.totalChiSquared(), (int)mumuTransientVtx.degreesOfFreedom());
0156   hSVProb_->Fill(SVProb);
0157   hSVChi2_->Fill(mumuTransientVtx.totalChiSquared());
0158   hSVNormChi2_->Fill(mumuTransientVtx.totalChiSquared() / (int)mumuTransientVtx.degreesOfFreedom());
0159 
0160   if (!mumuTransientVtx.isValid())
0161     return;
0162 
0163   const reco::Vertex* theClosestVertex;
0164   // get collection of reconstructed vertices from event
0165   edm::Handle<reco::VertexCollection> vertexHandle = iEvent.getHandle(vertexToken_);
0166   if (vertexHandle.isValid()) {
0167     const reco::VertexCollection* vertices = vertexHandle.product();
0168     theClosestVertex = this->findClosestVertex(mumuTransientVtx, vertices);
0169   } else {
0170     edm::LogWarning("DiMuonVertexMonitor") << "invalid vertex collection encountered Skipping event!";
0171     return;
0172   }
0173 
0174   reco::Vertex theMainVtx;
0175   if (!useClosestVertex_ || theClosestVertex == nullptr) {
0176     theMainVtx = *theClosestVertex;
0177   } else {
0178     theMainVtx = vertexHandle.product()->front();
0179   }
0180 
0181   const math::XYZPoint theMainVtxPos(theMainVtx.position().x(), theMainVtx.position().y(), theMainVtx.position().z());
0182   const math::XYZPoint myVertex(
0183       mumuTransientVtx.position().x(), mumuTransientVtx.position().y(), mumuTransientVtx.position().z());
0184   const math::XYZPoint deltaVtx(
0185       theMainVtxPos.x() - myVertex.x(), theMainVtxPos.y() - myVertex.y(), theMainVtxPos.z() - myVertex.z());
0186 
0187   if (theMainVtx.isValid()) {
0188     // fill the impact parameter plots
0189     for (const auto& track : myTracks) {
0190       hdxy_->Fill(track->dxy(theMainVtxPos) * cmToum);
0191       hdz_->Fill(track->dz(theMainVtxPos) * cmToum);
0192       hdxyErr_->Fill(track->dxyError() * cmToum);
0193       hdzErr_->Fill(track->dzError() * cmToum);
0194 
0195       const auto& ttrk = theB->build(track);
0196       Global3DVector dir(track->px(), track->py(), track->pz());
0197       const auto& ip2d = IPTools::signedTransverseImpactParameter(ttrk, dir, theMainVtx);
0198       const auto& ip3d = IPTools::signedImpactParameter3D(ttrk, dir, theMainVtx);
0199 
0200       hIP2d_->Fill(ip2d.second.value() * cmToum);
0201       hIP3d_->Fill(ip3d.second.value() * cmToum);
0202       hIP2dsig_->Fill(ip2d.second.significance());
0203       hIP3dsig_->Fill(ip3d.second.significance());
0204     }
0205 
0206     // Z Vertex distance in the xy plane
0207     VertexDistanceXY vertTool;
0208     double distance = vertTool.distance(mumuTransientVtx, theMainVtx).value();
0209     double dist_err = vertTool.distance(mumuTransientVtx, theMainVtx).error();
0210     float compatibility = 0.;
0211 
0212     try {
0213       compatibility = vertTool.compatibility(mumuTransientVtx, theMainVtx);
0214     } catch (cms::Exception& er) {
0215       LogTrace("DiMuonVertexMonitor") << "caught std::exception " << er.what() << std::endl;
0216     }
0217 
0218     hSVCompatibility_->Fill(compatibility);
0219     hSVDist_->Fill(distance * cmToum);
0220     hSVDistErr_->Fill(dist_err * cmToum);
0221     hSVDistSig_->Fill(distance / dist_err);
0222 
0223     // Z Vertex distance in 3D
0224     VertexDistance3D vertTool3D;
0225     double distance3D = vertTool3D.distance(mumuTransientVtx, theMainVtx).value();
0226     double dist3D_err = vertTool3D.distance(mumuTransientVtx, theMainVtx).error();
0227     float compatibility3D = 0.;
0228 
0229     try {
0230       compatibility3D = vertTool3D.compatibility(mumuTransientVtx, theMainVtx);
0231     } catch (cms::Exception& er) {
0232       LogTrace("DiMuonVertexMonitor") << "caught std::exception " << er.what() << std::endl;
0233     }
0234 
0235     hSVCompatibility3D_->Fill(compatibility3D);
0236     hSVDist3D_->Fill(distance3D * cmToum);
0237     hSVDist3DErr_->Fill(dist3D_err * cmToum);
0238     hSVDist3DSig_->Fill(distance3D / dist3D_err);
0239 
0240     // cut on the PV - SV distance
0241     if (distance * cmToum < maxSVdist_) {
0242       double cosphi = (ZpT.x() * deltaVtx.x() + ZpT.y() * deltaVtx.y()) /
0243                       (sqrt(ZpT.x() * ZpT.x() + ZpT.y() * ZpT.y()) *
0244                        sqrt(deltaVtx.x() * deltaVtx.x() + deltaVtx.y() * deltaVtx.y()));
0245 
0246       double cosphi3D = (Zp.x() * deltaVtx.x() + Zp.y() * deltaVtx.y() + Zp.z() * deltaVtx.z()) /
0247                         (sqrt(Zp.x() * Zp.x() + Zp.y() * Zp.y() + Zp.z() * Zp.z()) *
0248                          sqrt(deltaVtx.x() * deltaVtx.x() + deltaVtx.y() * deltaVtx.y() + deltaVtx.z() * deltaVtx.z()));
0249 
0250       hCosPhi_->Fill(cosphi);
0251       hCosPhi3D_->Fill(cosphi3D);
0252       // inverted
0253       hCosPhiInv_->Fill(-cosphi);
0254       hCosPhiInv3D_->Fill(-cosphi3D);
0255     }
0256   } else {
0257     edm::LogWarning("DiMuonVertexMonitor") << "hardest primary vertex in the event is not valid!";
0258   }
0259 }
0260 
0261 // compute the closest vertex to di-lepton ------------------------------------
0262 const reco::Vertex* DiMuonVertexMonitor::findClosestVertex(const TransientVertex aTransVtx,
0263                                                            const reco::VertexCollection* vertices) const {
0264   reco::Vertex* defaultVtx = nullptr;
0265 
0266   if (!aTransVtx.isValid())
0267     return defaultVtx;
0268 
0269   // find the closest vertex to the secondary vertex in 3D
0270   VertexDistance3D vertTool3D;
0271   float minD = 9999.;
0272   int closestVtxIndex = 0;
0273   int counter = 0;
0274   for (const auto& vtx : *vertices) {
0275     double dist3D = vertTool3D.distance(aTransVtx, vtx).value();
0276     if (dist3D < minD) {
0277       minD = dist3D;
0278       closestVtxIndex = counter;
0279     }
0280     counter++;
0281   }
0282 
0283   if ((*vertices).at(closestVtxIndex).isValid()) {
0284     return &(vertices->at(closestVtxIndex));
0285   } else {
0286     return defaultVtx;
0287   }
0288 }
0289 
0290 void DiMuonVertexMonitor::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0291   edm::ParameterSetDescription desc;
0292   desc.add<edm::InputTag>("muonTracks", edm::InputTag("ALCARECOTkAlDiMuon"));
0293   desc.add<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
0294   desc.add<std::string>("FolderName", "DiMuonVertexMonitor");
0295   desc.add<std::string>("decayMotherName", "Z");
0296   desc.add<bool>("useClosestVertex", true);
0297   desc.add<double>("maxSVdist", 50.);
0298   descriptions.addWithDefaultLabel(desc);
0299 }
0300 
0301 DEFINE_FWK_MODULE(DiMuonVertexMonitor);