Alignment/src/DiMuonMassBiasMonitor.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  */
0005
0006 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0007 #include "DQMOffline/Alignment/interface/DiMuonMassBiasMonitor.h"
0008 #include "DQMServices/Core/interface/DQMStore.h"
0009 #include "DataFormats/Histograms/interface/DQMToken.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; /* unused for now */
0026   constexpr float mumass2 = 0.105658367 * 0.105658367;  //mu mass squared (GeV^2/c^4)
0027 }  // namespace
0028
0029 DiMuonMassBiasMonitor::DiMuonMassBiasMonitor(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       beamSpotToken_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpot"))),
0034       MEFolderName_(iConfig.getParameter<std::string>("FolderName")),
0035       decayMotherName_(iConfig.getParameter<std::string>("decayMotherName")),
0036       distanceScaleFactor_(iConfig.getParameter<double>("distanceScaleFactor")),
0037       DiMuMassConfiguration_(iConfig.getParameter<edm::ParameterSet>("DiMuMassConfig")) {}
0038
0039 void DiMuonMassBiasMonitor::bookHistograms(DQMStore::IBooker& iBooker, edm::Run const&, edm::EventSetup const&) {
0040   iBooker.setCurrentFolder(MEFolderName_ + "/DiMuonMassBiasMonitor/MassBias");
0041   ZMassPlots.bookFromPSet(iBooker, DiMuMassConfiguration_);
0042
0043   iBooker.setCurrentFolder(MEFolderName_ + "/DiMuonMassBiasMonitor");
0044
0045   // retrieve the mass bins
0046   const auto& mass_bins = DiMuMassConfiguration_.getParameter<int32_t>("NyBins");
0047   const auto& mass_min = DiMuMassConfiguration_.getParameter<double>("ymin");
0048   const auto& mass_max = DiMuMassConfiguration_.getParameter<double>("ymax");
0049
0050   bookDecayHists(
0051       iBooker, histosZmm, decayMotherName_, "#mu^{+}#mu^{-}", mass_bins, mass_min, mass_max, distanceScaleFactor_);
0052
0053   iBooker.setCurrentFolder(MEFolderName_ + "/DiMuonMassBiasMonitor/components");
0054   bookDecayComponentHistograms(iBooker, histosZmm);
0055 }
0056
0057 void DiMuonMassBiasMonitor::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0058   std::vector<const reco::Track*> myTracks;
0059   const auto trackHandle = iEvent.getHandle(tracksToken_);
0060   if (!trackHandle.isValid()) {
0061     edm::LogError("DiMuonMassBiasMonitor") << "invalid track collection encountered!";
0062     return;
0063   }
0064
0065   for (const auto& muonTrk : *trackHandle) {
0066     myTracks.emplace_back(&muonTrk);
0067   }
0068
0069   if (myTracks.size() != 2) {
0070     edm::LogWarning("DiMuonMassBiasMonitor") << "There are not enough tracks to monitor!";
0071     return;
0072   }
0073
0074   const auto& t1 = myTracks[1]->momentum();
0075   const auto& t0 = myTracks[0]->momentum();
0076   const auto& ditrack = t1 + t0;
0077
0078   const auto& tplus = myTracks[0]->charge() > 0 ? myTracks[0] : myTracks[1];
0079   const auto& tminus = myTracks[0]->charge() < 0 ? myTracks[0] : myTracks[1];
0080
0081   TLorentzVector p4_tplus(tplus->px(), tplus->py(), tplus->pz(), sqrt((tplus->p() * tplus->p()) + mumass2));
0082   TLorentzVector p4_tminus(tminus->px(), tminus->py(), tminus->pz(), sqrt((tminus->p() * tminus->p()) + mumass2));
0083
0084   const auto& Zp4 = p4_tplus + p4_tminus;
0085   float track_invMass = Zp4.M();
0086
0087   // creat the pair of TLorentVectors used to make the plos
0088   std::pair<TLorentzVector, TLorentzVector> tktk_p4 = std::make_pair(p4_tplus, p4_tminus);
0089
0090   // fill the z->mm mass plots
0091   ZMassPlots.fillPlots(track_invMass, tktk_p4);
0092
0093   math::XYZPoint ZpT(ditrack.x(), ditrack.y(), 0);
0094   math::XYZPoint Zp(ditrack.x(), ditrack.y(), ditrack.z());
0095
0096   // get collection of reconstructed vertices from event
0097   const auto& vertexHandle = iEvent.getHandle(vertexToken_);
0098   if (!vertexHandle.isValid()) {
0099     edm::LogError("DiMuonMassBiasMonitor") << "invalid vertex collection encountered!";
0100     return;
0101   }
0102
0103   // get the vertices from the event
0104   const auto& vertices = vertexHandle.product();
0105
0106   // fill the decay vertex plots
0107   auto decayVertex = fillDecayHistograms(histosZmm, myTracks, vertices, iSetup);
0108
0109   // get the beamspot from the event
0110   const reco::BeamSpot* bs;
0111   const auto& beamSpotHandle = iEvent.getHandle(beamSpotToken_);
0112   if (!beamSpotHandle.isValid()) {
0113     bs = nullptr;
0114   } else {
0115     bs = beamSpotHandle.product();
0116   }
0117
0118   // fill the components plots
0119   fillComponentHistograms(histosZmm.decayComponents[0], tplus, bs, decayVertex);
0120   fillComponentHistograms(histosZmm.decayComponents[1], tminus, bs, decayVertex);
0121 }
0122
0123 void DiMuonMassBiasMonitor::bookDecayComponentHistograms(DQMStore::IBooker& ibook, DecayHists& histos) const {
0124   bookComponentHists(ibook, histos, "mu_plus", 0.1);
0125   bookComponentHists(ibook, histos, "mu_minus", 0.1);
0126 }
0127
0128 void DiMuonMassBiasMonitor::bookComponentHists(DQMStore::IBooker& ibook,
0129                                                DecayHists& histos,
0130                                                TString const& componentName,
0131                                                float distanceScaleFactor) const {
0132   ComponentHists comp;
0133
0134   comp.h_pt = ibook.book1D(componentName + "_pt", "track momentum ;p_{T} [GeV]", 100, 0., 100.);
0135   comp.h_eta = ibook.book1D(componentName + "_eta", "track rapidity;#eta", 100, -2.5, 2.5);
0136   comp.h_phi = ibook.book1D(componentName + "_phi", "track azimuth;#phi", 100, -M_PI, M_PI);
0137   comp.h_dxy = ibook.book1D(componentName + "_dxyBS",
0138                             "TIP w.r.t BS;d_{xy}(BS) [cm]",
0139                             100,
0140                             -0.5 * distanceScaleFactor,
0141                             0.5 * distanceScaleFactor);
0142   comp.h_exy =
0143       ibook.book1D(componentName + "_exy", "Error on TIP ;#sigma(d_{xy}(BS)) [cm]", 100, 0, 0.05 * distanceScaleFactor);
0144   comp.h_dz = ibook.book1D(componentName + "_dzPV",
0145                            "LIP w.r.t PV;d_{z}(PV) [cm]",
0146                            100,
0147                            -0.5 * distanceScaleFactor,
0148                            0.5 * distanceScaleFactor);
0149   comp.h_ez =
0150       ibook.book1D(componentName + "_ez", "Error on LIP;#sigma(d_{z}(PV)) [cm]", 100, 0, 0.5 * distanceScaleFactor);
0151   comp.h_chi2 = ibook.book1D(componentName + "_chi2", ";#chi^{2}", 100, 0, 20);
0152
0153   histos.decayComponents.push_back(comp);
0154 }
0155
0156 void DiMuonMassBiasMonitor::bookDecayHists(DQMStore::IBooker& ibook,
0157                                            DecayHists& decayHists,
0158                                            std::string const& name,
0159                                            std::string const& products,
0160                                            int nMassBins,
0161                                            float massMin,
0162                                            float massMax,
0163                                            float distanceScaleFactor) const {
0164   std::string histTitle = name + " #rightarrow " + products + ";";
0165
0166   decayHists.h_mass = ibook.book1D("h_mass", histTitle + "M(" + products + ") [GeV]", nMassBins, massMin, massMax);
0167   decayHists.h_pt = ibook.book1D("h_pt", histTitle + "p_{T} [GeV]", 100, 0.00, 200.0);
0168   decayHists.h_eta = ibook.book1D("h_eta", histTitle + "#eta", 100, -5., 5.);
0169   decayHists.h_phi = ibook.book1D("h_phi", histTitle + "#varphi [rad]", 100, -M_PI, M_PI);
0170   decayHists.h_displ2D =
0171       ibook.book1D("h_displ2D", histTitle + "vertex 2D displacement [cm]", 100, 0.00, 0.1 * distanceScaleFactor);
0172   decayHists.h_sign2D =
0173       ibook.book1D("h_sign2D", histTitle + "vertex 2D displ. significance", 100, 0.00, 100.0 * distanceScaleFactor);
0174   decayHists.h_ct = ibook.book1D("h_ct", histTitle + "c#tau [cm]", 100, 0.00, 0.4 * distanceScaleFactor);
0175   decayHists.h_pointing = ibook.book1D("h_pointing", histTitle + "cos( 2D pointing angle )", 100, -1, 1);
0176   decayHists.h_vertNormChi2 = ibook.book1D("h_vertNormChi2", histTitle + "vertex #chi^{2}/ndof", 100, 0.00, 10);
0177   decayHists.h_vertProb = ibook.book1D("h_vertProb", histTitle + "vertex prob.", 100, 0.00, 1.0);
0178 }
0179
0180 reco::Vertex const* DiMuonMassBiasMonitor::fillDecayHistograms(DecayHists const& histos,
0181                                                                std::vector<const reco::Track*> const& tracks,
0182                                                                const reco::VertexCollection* const& pvs,
0183                                                                const edm::EventSetup& iSetup) const {
0184   if (tracks.size() != 2) {
0185     edm::LogWarning("DiMuonVertexMonitor") << "There are not enough tracks to construct a vertex!";
0186     return nullptr;
0187   }
0188
0189   const TransientTrackBuilder* theB = &iSetup.getData(ttbESToken_);
0190   TransientVertex mumuTransientVtx;
0191   std::vector<reco::TransientTrack> tks;
0192
0193   for (const auto& track : tracks) {
0194     reco::TransientTrack trajectory = theB->build(track);
0195     tks.push_back(trajectory);
0196   }
0197
0198   KalmanVertexFitter kalman(true);
0199   mumuTransientVtx = kalman.vertex(tks);
0200
0201   auto svtx = reco::Vertex(mumuTransientVtx);
0202   if (not svtx.isValid()) {
0203     return nullptr;
0204   }
0205
0206   const auto& mom_t1 = tracks[1]->momentum();
0207   const auto& mom_t0 = tracks[0]->momentum();
0208   const auto& momentum = mom_t1 + mom_t0;
0209
0210   TLorentzVector p4_t0(mom_t0.x(), mom_t0.y(), mom_t0.z(), sqrt(mom_t0.mag2() + mumass2));
0211   TLorentzVector p4_t1(mom_t1.x(), mom_t1.y(), mom_t1.z(), sqrt(mom_t1.mag2() + mumass2));
0212
0213   const auto& p4 = p4_t0 + p4_t1;
0214   float mass = p4.M();
0215
0216   auto pvtx = std::min_element(pvs->begin(), pvs->end(), [svtx](reco::Vertex const& pv1, reco::Vertex const& pv2) {
0217     return abs(pv1.z() - svtx.z()) < abs(pv2.z() - svtx.z());
0218   });
0219
0220   if (pvtx == pvs->end()) {
0221     return nullptr;
0222   }
0223
0224   VertexDistanceXY vdistXY;
0225   Measurement1D distXY = vdistXY.distance(svtx, *pvtx);
0226
0227   auto pvtPos = pvtx->position();
0228   const auto& svtPos = svtx.position();
0229
0230   math::XYZVector displVect2D(svtPos.x() - pvtPos.x(), svtPos.y() - pvtPos.y(), 0);
0231   auto cosAlpha = displVect2D.Dot(momentum) / (displVect2D.Rho() * momentum.rho());
0232
0233   auto ct = distXY.value() * cosAlpha * mass / momentum.rho();
0234
0235   histos.h_pointing->Fill(cosAlpha);
0236
0237   histos.h_mass->Fill(mass);
0238
0239   histos.h_pt->Fill(momentum.rho());
0240   histos.h_eta->Fill(momentum.eta());
0241   histos.h_phi->Fill(momentum.phi());
0242
0243   histos.h_ct->Fill(ct);
0244
0245   histos.h_displ2D->Fill(distXY.value());
0246   histos.h_sign2D->Fill(distXY.significance());
0247
0248   if (svtx.chi2() >= 0) {
0249     histos.h_vertNormChi2->Fill(svtx.chi2() / svtx.ndof());
0250     histos.h_vertProb->Fill(ChiSquaredProbability(svtx.chi2(), svtx.ndof()));
0251   }
0252
0253   return &*pvtx;
0254 }
0255
0256 void DiMuonMassBiasMonitor::fillComponentHistograms(ComponentHists const& histos,
0257                                                     const reco::Track* const& component,
0258                                                     reco::BeamSpot const* bs,
0259                                                     reco::Vertex const* pv) const {
0260   histos.h_pt->Fill(component->pt());
0261   histos.h_eta->Fill(component->eta());
0262   histos.h_phi->Fill(component->phi());
0263
0264   math::XYZPoint zero(0, 0, 0);
0265   math::Error<3>::type zeroCov;  // needed for dxyError
0266   if (bs) {
0267     histos.h_dxy->Fill(component->dxy(*bs));
0268     histos.h_exy->Fill(component->dxyError(*bs));
0269   } else {
0270     histos.h_dxy->Fill(component->dxy(zero));
0271     histos.h_exy->Fill(component->dxyError(zero, zeroCov));
0272   }
0273   if (pv) {
0274     histos.h_dz->Fill(component->dz(pv->position()));
0275   } else {
0276     histos.h_dz->Fill(component->dz(zero));
0277   }
0278   histos.h_ez->Fill(component->dzError());
0279
0280   histos.h_chi2->Fill(component->chi2() / component->ndof());
0281 }
0282
0283 void DiMuonMassBiasMonitor::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0284   edm::ParameterSetDescription desc;
0285   desc.add<edm::InputTag>("muonTracks", edm::InputTag("ALCARECOTkAlDiMuon"));
0286   desc.add<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
0287   desc.add<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
0288   desc.add<std::string>("FolderName", "DiMuonMassBiasMonitor");
0289   desc.add<std::string>("decayMotherName", "Z");
0290   desc.add<double>("distanceScaleFactor", 0.1);
0291
0292   {
0293     edm::ParameterSetDescription psDiMuMass;
0294     psDiMuMass.add<std::string>("name", "DiMuMass");
0295     psDiMuMass.add<std::string>("title", "M(#mu#mu)");
0296     psDiMuMass.add<std::string>("yUnits", "[GeV]");
0297     psDiMuMass.add<int>("NxBins", 24);
0298     psDiMuMass.add<int>("NyBins", 50);
0299     psDiMuMass.add<double>("ymin", 65.);
0300     psDiMuMass.add<double>("ymax", 115.);
0301     psDiMuMass.add<double>("maxDeltaEta", 3.7);
0302     desc.add<edm::ParameterSetDescription>("DiMuMassConfig", psDiMuMass);
0303   }
0304
0305   descriptions.addWithDefaultLabel(desc);
0306 }
0307
0308 DEFINE_FWK_MODULE(DiMuonMassBiasMonitor);