Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-02-08 02:59:45

0001 // -*- C++ -*-
0002 //
0003 // Package:    Alignment/OfflineValidation
0004 // Class:      DiMuonVertexValidation
0005 //
0006 /**\class DiMuonVertexValidation DiMuonVertexValidation.cc Alignment/OfflineValidation/plugins/DiMuonVertexValidation.cc
0007 
0008  Description: Class to perform validation Tracker Alignment Validations by means of a PV and the SV constructed with a di-muon pair
0009 
0010 */
0011 //
0012 // Original Author:  Marco Musich
0013 //         Created:  Wed, 21 Apr 2021 09:06:25 GMT
0014 //
0015 //
0016 
0017 // system include files
0018 #include <memory>
0019 #include <utility>
0020 
0021 // user include files
0022 #include "FWCore/Framework/interface/Frameworkfwd.h"
0023 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/MakerMacros.h"
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027 #include "FWCore/Utilities/interface/InputTag.h"
0028 
0029 // muons
0030 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0031 #include "DataFormats/MuonReco/interface/Muon.h"
0032 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0033 
0034 // utils
0035 #include "Alignment/OfflineValidation/interface/DiLeptonVertexHelpers.h"
0036 #include "DataFormats/Math/interface/deltaR.h"
0037 #include "TLorentzVector.h"
0038 
0039 // tracks
0040 #include "DataFormats/TrackReco/interface/Track.h"
0041 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0042 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0043 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0044 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0045 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0046 
0047 // vertices
0048 #include "RecoVertex/VertexTools/interface/VertexDistanceXY.h"
0049 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0050 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0051 #include "FWCore/ServiceRegistry/interface/Service.h"
0052 
0053 // ROOT
0054 #include "TH1F.h"
0055 #include "TH2F.h"
0056 
0057 //#define LogDebug(X) std::cout << X <<
0058 
0059 //
0060 // class declaration
0061 //
0062 
0063 class DiMuonVertexValidation : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0064 public:
0065   explicit DiMuonVertexValidation(const edm::ParameterSet&);
0066   ~DiMuonVertexValidation() override;
0067 
0068   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0069 
0070 private:
0071   void beginJob() override;
0072   void analyze(const edm::Event&, const edm::EventSetup&) override;
0073   const reco::Vertex* findClosestVertex(const TransientVertex aTransVtx, const reco::VertexCollection* vertices) const;
0074   void endJob() override;
0075 
0076   // ----------member data ---------------------------
0077   DiLeptonHelp::Counts myCounts;
0078 
0079   const bool useReco_;
0080   const bool useClosestVertex_;
0081   std::vector<double> pTthresholds_;
0082   const float maxSVdist_;
0083 
0084   // plot configurations
0085   const edm::ParameterSet CosPhiConfiguration_;
0086   const edm::ParameterSet CosPhi3DConfiguration_;
0087   const edm::ParameterSet VtxProbConfiguration_;
0088   const edm::ParameterSet VtxDistConfiguration_;
0089   const edm::ParameterSet VtxDist3DConfiguration_;
0090   const edm::ParameterSet VtxDistSigConfiguration_;
0091   const edm::ParameterSet VtxDist3DSigConfiguration_;
0092   const edm::ParameterSet DiMuMassConfiguration_;
0093 
0094   // control plots
0095 
0096   TH1F* hSVProb_;
0097   TH1F* hSVDist_;
0098   TH1F* hSVDistSig_;
0099   TH1F* hSVDist3D_;
0100   TH1F* hSVDist3DSig_;
0101 
0102   TH1F* hCosPhi_;
0103   TH1F* hCosPhi3D_;
0104   TH1F* hCosPhiInv_;
0105   TH1F* hCosPhiInv3D_;
0106 
0107   TH1F* hInvMass_;
0108   TH1F* hTrackInvMass_;
0109 
0110   TH1F* hCutFlow_;
0111 
0112   // 2D maps
0113 
0114   DiLeptonHelp::PlotsVsKinematics CosPhiPlots = DiLeptonHelp::PlotsVsKinematics(DiLeptonHelp::MM);
0115   DiLeptonHelp::PlotsVsKinematics CosPhi3DPlots = DiLeptonHelp::PlotsVsKinematics(DiLeptonHelp::MM);
0116   DiLeptonHelp::PlotsVsKinematics VtxProbPlots = DiLeptonHelp::PlotsVsKinematics(DiLeptonHelp::MM);
0117   DiLeptonHelp::PlotsVsKinematics VtxDistPlots = DiLeptonHelp::PlotsVsKinematics(DiLeptonHelp::MM);
0118   DiLeptonHelp::PlotsVsKinematics VtxDist3DPlots = DiLeptonHelp::PlotsVsKinematics(DiLeptonHelp::MM);
0119   DiLeptonHelp::PlotsVsKinematics VtxDistSigPlots = DiLeptonHelp::PlotsVsKinematics(DiLeptonHelp::MM);
0120   DiLeptonHelp::PlotsVsKinematics VtxDist3DSigPlots = DiLeptonHelp::PlotsVsKinematics(DiLeptonHelp::MM);
0121   DiLeptonHelp::PlotsVsKinematics ZMassPlots = DiLeptonHelp::PlotsVsKinematics(DiLeptonHelp::MM);
0122 
0123   // plots vs region
0124   DiLeptonHelp::PlotsVsDiLeptonRegion CosPhi3DInEtaBins = DiLeptonHelp::PlotsVsDiLeptonRegion(1.5);
0125   DiLeptonHelp::PlotsVsDiLeptonRegion InvMassInEtaBins = DiLeptonHelp::PlotsVsDiLeptonRegion(1.5);
0126 
0127   const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> ttbESToken_;
0128 
0129   //used to select what tracks to read from configuration file
0130   edm::EDGetTokenT<reco::TrackCollection> tracksToken_;
0131   //used to select what vertices to read from configuration file
0132   const edm::EDGetTokenT<reco::VertexCollection> vertexToken_;
0133 
0134   // either on or the other!
0135   edm::EDGetTokenT<reco::MuonCollection> muonsToken_;      // used to select tracks to read from configuration file
0136   edm::EDGetTokenT<reco::TrackCollection> alcaRecoToken_;  //used to select muon tracks to read from configuration file
0137 };
0138 
0139 //
0140 // constants, enums and typedefs
0141 //
0142 
0143 static constexpr float cmToum = 10e4;
0144 static constexpr float mumass2 = 0.105658367 * 0.105658367;  //mu mass squared (GeV^2/c^4)
0145 
0146 //
0147 // static data member definitions
0148 //
0149 
0150 //
0151 // constructors and destructor
0152 //
0153 DiMuonVertexValidation::DiMuonVertexValidation(const edm::ParameterSet& iConfig)
0154     : useReco_(iConfig.getParameter<bool>("useReco")),
0155       useClosestVertex_(iConfig.getParameter<bool>("useClosestVertex")),
0156       pTthresholds_(iConfig.getParameter<std::vector<double>>("pTThresholds")),
0157       maxSVdist_(iConfig.getParameter<double>("maxSVdist")),
0158       CosPhiConfiguration_(iConfig.getParameter<edm::ParameterSet>("CosPhiConfig")),
0159       CosPhi3DConfiguration_(iConfig.getParameter<edm::ParameterSet>("CosPhi3DConfig")),
0160       VtxProbConfiguration_(iConfig.getParameter<edm::ParameterSet>("VtxProbConfig")),
0161       VtxDistConfiguration_(iConfig.getParameter<edm::ParameterSet>("VtxDistConfig")),
0162       VtxDist3DConfiguration_(iConfig.getParameter<edm::ParameterSet>("VtxDist3DConfig")),
0163       VtxDistSigConfiguration_(iConfig.getParameter<edm::ParameterSet>("VtxDistSigConfig")),
0164       VtxDist3DSigConfiguration_(iConfig.getParameter<edm::ParameterSet>("VtxDist3DSigConfig")),
0165       DiMuMassConfiguration_(iConfig.getParameter<edm::ParameterSet>("DiMuMassConfig")),
0166       ttbESToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0167       vertexToken_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))) {
0168   if (useReco_) {
0169     tracksToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"));
0170     muonsToken_ = consumes<reco::MuonCollection>(iConfig.getParameter<edm::InputTag>("muons"));
0171   } else {
0172     alcaRecoToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("muonTracks"));
0173   }
0174 
0175   usesResource(TFileService::kSharedResource);
0176 
0177   // sort the vector of thresholds
0178   std::sort(pTthresholds_.begin(), pTthresholds_.end(), [](const double& lhs, const double& rhs) { return lhs > rhs; });
0179 
0180   edm::LogInfo("DiMuonVertexValidation") << __FUNCTION__;
0181   for (const auto& thr : pTthresholds_) {
0182     edm::LogInfo("DiMuonVertexValidation") << " Threshold: " << thr << " ";
0183   }
0184   edm::LogInfo("DiMuonVertexValidation") << "Max SV distance: " << maxSVdist_ << " ";
0185 }
0186 
0187 DiMuonVertexValidation::~DiMuonVertexValidation() = default;
0188 
0189 //
0190 // member functions
0191 //
0192 
0193 // ------------ method called for each event  ------------
0194 void DiMuonVertexValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0195   using namespace edm;
0196 
0197   myCounts.eventsTotal++;
0198 
0199   // the di-muon tracks
0200   std::vector<const reco::Track*> myTracks;
0201 
0202   // if we have to start from scratch from RECO data-tier
0203   if (useReco_) {
0204     // select the good muons
0205     std::vector<const reco::Muon*> myGoodMuonVector;
0206     for (const auto& muon : iEvent.get(muonsToken_)) {
0207       const reco::TrackRef t = muon.innerTrack();
0208       if (!t.isNull()) {
0209         if (t->quality(reco::TrackBase::highPurity)) {
0210           if (t->chi2() / t->ndof() <= 2.5 && t->numberOfValidHits() >= 5 &&
0211               t->hitPattern().numberOfValidPixelHits() >= 2 && t->quality(reco::TrackBase::highPurity))
0212             myGoodMuonVector.emplace_back(&muon);
0213         }
0214       }
0215     }
0216 
0217     LogDebug("DiMuonVertexValidation") << "myGoodMuonVector size: " << myGoodMuonVector.size() << std::endl;
0218     std::sort(myGoodMuonVector.begin(), myGoodMuonVector.end(), [](const reco::Muon*& lhs, const reco::Muon*& rhs) {
0219       return lhs->pt() > rhs->pt();
0220     });
0221 
0222     // just check the ordering
0223     for (const auto& muon : myGoodMuonVector) {
0224       LogDebug("DiMuonVertexValidation") << "pT: " << muon->pt() << " ";
0225     }
0226     LogDebug("DiMuonVertexValidation") << std::endl;
0227 
0228     // reject if there's no Z
0229     if (myGoodMuonVector.size() < 2)
0230       return;
0231 
0232     myCounts.eventsAfterMult++;
0233 
0234     if ((myGoodMuonVector[0]->pt()) < pTthresholds_[0] || (myGoodMuonVector[1]->pt() < pTthresholds_[1]))
0235       return;
0236 
0237     myCounts.eventsAfterPt++;
0238     myCounts.eventsAfterEta++;
0239 
0240     if (myGoodMuonVector[0]->charge() * myGoodMuonVector[1]->charge() > 0)
0241       return;
0242 
0243     const auto& m1 = myGoodMuonVector[1]->p4();
0244     const auto& m0 = myGoodMuonVector[0]->p4();
0245     const auto& mother = m0 + m1;
0246 
0247     float invMass = mother.M();
0248     hInvMass_->Fill(invMass);
0249 
0250     // just copy the top two muons
0251     std::vector<const reco::Muon*> theZMuonVector;
0252     theZMuonVector.reserve(2);
0253     theZMuonVector.emplace_back(myGoodMuonVector[1]);
0254     theZMuonVector.emplace_back(myGoodMuonVector[0]);
0255 
0256     // do the matching of Z muons with inner tracks
0257     unsigned int i = 0;
0258     for (const auto& muon : theZMuonVector) {
0259       i++;
0260       float minD = 1000.;
0261       const reco::Track* theMatch = nullptr;
0262       for (const auto& track : iEvent.get(tracksToken_)) {
0263         float D = ::deltaR(muon->eta(), muon->phi(), track.eta(), track.phi());
0264         if (D < minD) {
0265           minD = D;
0266           theMatch = &track;
0267         }
0268       }
0269       LogDebug("DiMuonVertexValidation") << "pushing new track: " << i << std::endl;
0270       myTracks.emplace_back(theMatch);
0271     }
0272   } else {
0273     // we start directly with the pre-selected ALCARECO tracks
0274     for (const auto& muon : iEvent.get(alcaRecoToken_)) {
0275       myTracks.emplace_back(&muon);
0276     }
0277   }
0278 
0279 #ifdef EDM_ML_DEBUG
0280   for (const auto& track : myTracks) {
0281     edm::LogVerbatim("DiMuonVertexValidation") << __PRETTY_FUNCTION__ << " pT: " << track->pt() << " GeV"
0282                                                << " , pT error: " << track->ptError() << " GeV"
0283                                                << " , eta: " << track->eta() << " , phi: " << track->phi() << std::endl;
0284   }
0285 #endif
0286 
0287   LogDebug("DiMuonVertexValidation") << "selected tracks: " << myTracks.size() << std::endl;
0288 
0289   const TransientTrackBuilder* theB = &iSetup.getData(ttbESToken_);
0290   TransientVertex aTransientVertex;
0291   std::vector<reco::TransientTrack> tks;
0292 
0293   if (myTracks.size() != 2)
0294     return;
0295 
0296   const auto& t1 = myTracks[1]->momentum();
0297   const auto& t0 = myTracks[0]->momentum();
0298   const auto& ditrack = t1 + t0;
0299 
0300   const auto& tplus = myTracks[0]->charge() > 0 ? myTracks[0] : myTracks[1];
0301   const auto& tminus = myTracks[0]->charge() < 0 ? myTracks[0] : myTracks[1];
0302 
0303   TLorentzVector p4_tplus(tplus->px(), tplus->py(), tplus->pz(), sqrt((tplus->p() * tplus->p()) + mumass2));
0304   TLorentzVector p4_tminus(tminus->px(), tminus->py(), tminus->pz(), sqrt((tminus->p() * tminus->p()) + mumass2));
0305 
0306 #ifdef EDM_ML_DEBUG
0307   // Define a lambda function to convert TLorentzVector to a string
0308   auto tLorentzVectorToString = [](const TLorentzVector& vector) {
0309     return std::to_string(vector.Px()) + " " + std::to_string(vector.Py()) + " " + std::to_string(vector.Pz()) + " " +
0310            std::to_string(vector.E());
0311   };
0312 
0313   edm::LogVerbatim("DiMuonVertexValidation") << "mu+" << tLorentzVectorToString(p4_tplus) << std::endl;
0314   edm::LogVerbatim("DiMuonVertexValidation") << "mu-" << tLorentzVectorToString(p4_tminus) << std::endl;
0315 #endif
0316 
0317   const auto& Zp4 = p4_tplus + p4_tminus;
0318   float track_invMass = Zp4.M();
0319   hTrackInvMass_->Fill(track_invMass);
0320 
0321   // creat the pair of TLorentVectors used to make the plos
0322   std::pair<TLorentzVector, TLorentzVector> tktk_p4 = std::make_pair(p4_tplus, p4_tminus);
0323 
0324   // fill the z->mm mass plots
0325   ZMassPlots.fillPlots(track_invMass, tktk_p4);
0326   InvMassInEtaBins.fillTH1Plots(track_invMass, tktk_p4);
0327 
0328   math::XYZPoint ZpT(ditrack.x(), ditrack.y(), 0);
0329   math::XYZPoint Zp(ditrack.x(), ditrack.y(), ditrack.z());
0330 
0331   for (const auto& track : myTracks) {
0332     reco::TransientTrack trajectory = theB->build(track);
0333     tks.push_back(trajectory);
0334   }
0335 
0336   KalmanVertexFitter kalman(true);
0337   aTransientVertex = kalman.vertex(tks);
0338 
0339   double SVProb = TMath::Prob(aTransientVertex.totalChiSquared(), (int)aTransientVertex.degreesOfFreedom());
0340 
0341   LogDebug("DiMuonVertexValidation") << " vertex prob: " << SVProb << std::endl;
0342 
0343   hSVProb_->Fill(SVProb);
0344 
0345   if (!aTransientVertex.isValid())
0346     return;
0347 
0348   myCounts.eventsAfterVtx++;
0349 
0350   // fill the VtxProb plots
0351   VtxProbPlots.fillPlots(SVProb, tktk_p4);
0352 
0353   math::XYZPoint mainVtxPos(0, 0, 0);
0354   const reco::Vertex* theClosestVertex = nullptr;
0355   // get collection of reconstructed vertices from event
0356   edm::Handle<reco::VertexCollection> vertexHandle = iEvent.getHandle(vertexToken_);
0357   if (vertexHandle.isValid()) {
0358     const reco::VertexCollection* vertices = vertexHandle.product();
0359     theClosestVertex = this->findClosestVertex(aTransientVertex, vertices);
0360   } else {
0361     edm::LogWarning("DiMuonVertexMonitor") << "invalid vertex collection encountered Skipping event!";
0362     return;
0363   }
0364 
0365   reco::Vertex theMainVertex;
0366   if (!useClosestVertex_ || theClosestVertex == nullptr) {
0367     // if the closest vertex is not available, or explicitly not chosen
0368     theMainVertex = vertexHandle.product()->front();
0369   } else {
0370     theMainVertex = *theClosestVertex;
0371   }
0372 
0373   mainVtxPos.SetXYZ(theMainVertex.position().x(), theMainVertex.position().y(), theMainVertex.position().z());
0374   const math::XYZPoint myVertex(
0375       aTransientVertex.position().x(), aTransientVertex.position().y(), aTransientVertex.position().z());
0376   const math::XYZPoint deltaVtx(
0377       mainVtxPos.x() - myVertex.x(), mainVtxPos.y() - myVertex.y(), mainVtxPos.z() - myVertex.z());
0378 
0379 #ifdef EDM_ML_DEBUG
0380   edm::LogVerbatim("DiMuonVertexValidation")
0381       << "mm vertex position:" << aTransientVertex.position().x() << "," << aTransientVertex.position().y() << ","
0382       << aTransientVertex.position().z();
0383 
0384   edm::LogVerbatim("DiMuonVertexValidation") << "main vertex position:" << theMainVertex.position().x() << ","
0385                                              << theMainVertex.position().y() << "," << theMainVertex.position().z();
0386 #endif
0387 
0388   if (theMainVertex.isValid()) {
0389     // Z Vertex distance in the xy plane
0390 
0391     VertexDistanceXY vertTool;
0392     double distance = vertTool.distance(aTransientVertex, theMainVertex).value();
0393     double dist_err = vertTool.distance(aTransientVertex, theMainVertex).error();
0394 
0395     hSVDist_->Fill(distance * cmToum);
0396     hSVDistSig_->Fill(distance / dist_err);
0397 
0398     // fill the VtxDist plots
0399     VtxDistPlots.fillPlots(distance * cmToum, tktk_p4);
0400 
0401     // fill the VtxDisSig plots
0402     VtxDistSigPlots.fillPlots(distance / dist_err, tktk_p4);
0403 
0404     // Z Vertex distance in 3D
0405 
0406     VertexDistance3D vertTool3D;
0407     double distance3D = vertTool3D.distance(aTransientVertex, theMainVertex).value();
0408     double dist3D_err = vertTool3D.distance(aTransientVertex, theMainVertex).error();
0409 
0410     hSVDist3D_->Fill(distance3D * cmToum);
0411     hSVDist3DSig_->Fill(distance3D / dist3D_err);
0412 
0413     // fill the VtxDist3D plots
0414     VtxDist3DPlots.fillPlots(distance3D * cmToum, tktk_p4);
0415 
0416     // fill the VtxDisSig plots
0417     VtxDist3DSigPlots.fillPlots(distance3D / dist3D_err, tktk_p4);
0418 
0419     LogDebug("DiMuonVertexValidation") << "distance: " << distance << "+/-" << dist_err << std::endl;
0420     // cut on the PV - SV distance
0421     if (distance * cmToum < maxSVdist_) {
0422       myCounts.eventsAfterDist++;
0423 
0424       double cosphi = (ZpT.x() * deltaVtx.x() + ZpT.y() * deltaVtx.y()) /
0425                       (sqrt(ZpT.x() * ZpT.x() + ZpT.y() * ZpT.y()) *
0426                        sqrt(deltaVtx.x() * deltaVtx.x() + deltaVtx.y() * deltaVtx.y()));
0427 
0428       double cosphi3D = (Zp.x() * deltaVtx.x() + Zp.y() * deltaVtx.y() + Zp.z() * deltaVtx.z()) /
0429                         (sqrt(Zp.x() * Zp.x() + Zp.y() * Zp.y() + Zp.z() * Zp.z()) *
0430                          sqrt(deltaVtx.x() * deltaVtx.x() + deltaVtx.y() * deltaVtx.y() + deltaVtx.z() * deltaVtx.z()));
0431 
0432       LogDebug("DiMuonVertexValidation") << "cos(phi) = " << cosphi << std::endl;
0433 
0434       hCosPhi_->Fill(cosphi);
0435       hCosPhi3D_->Fill(cosphi3D);
0436 
0437 #ifdef EDM_ML_DEBUG
0438       edm::LogVerbatim("DiMuonVertexValidation")
0439           << "distance " << distance * cmToum << " cosphi3D:" << cosphi3D << std::endl;
0440 #endif
0441 
0442       // fill the cosphi plots
0443       CosPhiPlots.fillPlots(cosphi, tktk_p4);
0444 
0445       // fill the cosphi3D plots
0446       CosPhi3DPlots.fillPlots(cosphi3D, tktk_p4);
0447 
0448       // fill the cosphi3D plots in eta bins
0449       CosPhi3DInEtaBins.fillTH1Plots(cosphi3D, tktk_p4);
0450     }
0451   }
0452 }
0453 
0454 // ------------ method called once each job just before starting event loop  ------------
0455 void DiMuonVertexValidation::beginJob() {
0456   edm::Service<TFileService> fs;
0457 
0458   // clang-format off
0459   TH1F::SetDefaultSumw2(kTRUE);
0460   hSVProb_ = fs->make<TH1F>("VtxProb", ";#mu^{+}#mu^{-} vertex probability;N(#mu#mu pairs)", 100, 0., 1.);
0461 
0462   auto extractRangeValues = [](const edm::ParameterSet& PSetConfiguration_) -> std::pair<double, double> {
0463     double min = PSetConfiguration_.getParameter<double>("ymin");
0464     double max = PSetConfiguration_.getParameter<double>("ymax");
0465     return {min, max};
0466   };
0467 
0468   // take the range from the 2D histograms
0469   const auto& svDistRng = extractRangeValues(VtxDistConfiguration_);
0470   hSVDist_ = fs->make<TH1F>("VtxDist", ";PV-#mu^{+}#mu^{-} vertex xy distance [#mum];N(#mu#mu pairs)", 100, svDistRng.first, svDistRng.second);
0471 
0472   // take the range from the 2D histograms
0473   const auto& svDistSigRng = extractRangeValues(VtxDistSigConfiguration_);
0474   hSVDistSig_ = fs->make<TH1F>("VtxDistSig", ";PV-#mu^{+}#mu^{-} vertex xy distance signficance;N(#mu#mu pairs)", 100, svDistSigRng.first, svDistRng.second);
0475 
0476   // take the range from the 2D histograms
0477   const auto& svDist3DRng = extractRangeValues(VtxDist3DConfiguration_);
0478   hSVDist3D_ = fs->make<TH1F>("VtxDist3D", ";PV-#mu^{+}#mu^{-} vertex 3D distance [#mum];N(#mu#mu pairs)", 100, svDist3DRng.first, svDist3DRng.second);
0479 
0480   // take the range from the 2D histograms
0481   const auto& svDist3DSigRng = extractRangeValues(VtxDist3DSigConfiguration_);
0482   hSVDist3DSig_ = fs->make<TH1F>("VtxDist3DSig", ";PV-#mu^{+}#mu^{-} vertex 3D distance signficance;N(#mu#mu pairs)", 100, svDist3DSigRng.first, svDist3DSigRng.second);
0483 
0484   // take the range from the 2D histograms
0485   const auto& massRng = extractRangeValues(DiMuMassConfiguration_);
0486   hInvMass_ = fs->make<TH1F>("InvMass", ";M(#mu#mu) [GeV];N(#mu#mu pairs)", 70., massRng.first, massRng.second);
0487   hTrackInvMass_ = fs->make<TH1F>("TkTkInvMass", ";M(tk,tk) [GeV];N(tk tk pairs)", 70., massRng.first, massRng.second);
0488 
0489   hCosPhi_ = fs->make<TH1F>("CosPhi", ";cos(#phi_{xy});N(#mu#mu pairs)", 50, -1., 1.);
0490   hCosPhi3D_ = fs->make<TH1F>("CosPhi3D", ";cos(#phi_{3D});N(#mu#mu pairs)", 50, -1., 1.);
0491 
0492   hCosPhiInv_ = fs->make<TH1F>("CosPhiInv", ";inverted cos(#phi_{xy});N(#mu#mu pairs)", 50, -1., 1.);
0493   hCosPhiInv3D_ = fs->make<TH1F>("CosPhiInv3D", ";inverted cos(#phi_{3D});N(#mu#mu pairs)", 50, -1., 1.);
0494   // clang-format on
0495 
0496   // 2D Maps
0497 
0498   TFileDirectory dirCosPhi = fs->mkdir("CosPhiPlots");
0499   CosPhiPlots.bookFromPSet(dirCosPhi, CosPhiConfiguration_);
0500 
0501   TFileDirectory dirCosPhi3D = fs->mkdir("CosPhi3DPlots");
0502   CosPhi3DPlots.bookFromPSet(dirCosPhi3D, CosPhi3DConfiguration_);
0503 
0504   TFileDirectory dirVtxProb = fs->mkdir("VtxProbPlots");
0505   VtxProbPlots.bookFromPSet(dirVtxProb, VtxProbConfiguration_);
0506 
0507   TFileDirectory dirVtxDist = fs->mkdir("VtxDistPlots");
0508   VtxDistPlots.bookFromPSet(dirVtxDist, VtxDistConfiguration_);
0509 
0510   TFileDirectory dirVtxDist3D = fs->mkdir("VtxDist3DPlots");
0511   VtxDist3DPlots.bookFromPSet(dirVtxDist3D, VtxDist3DConfiguration_);
0512 
0513   TFileDirectory dirVtxDistSig = fs->mkdir("VtxDistSigPlots");
0514   VtxDistSigPlots.bookFromPSet(dirVtxDistSig, VtxDistSigConfiguration_);
0515 
0516   TFileDirectory dirVtxDist3DSig = fs->mkdir("VtxDist3DSigPlots");
0517   VtxDist3DSigPlots.bookFromPSet(dirVtxDist3DSig, VtxDist3DSigConfiguration_);
0518 
0519   TFileDirectory dirInvariantMass = fs->mkdir("InvariantMassPlots");
0520   ZMassPlots.bookFromPSet(dirInvariantMass, DiMuMassConfiguration_);
0521 
0522   // CosPhi3D in eta bins
0523   TFileDirectory dirCosphi3DEta = fs->mkdir("CosPhi3DInEtaBins");
0524   CosPhi3DInEtaBins.bookSet(dirCosphi3DEta, hCosPhi3D_);
0525 
0526   // Z-> mm mass in eta bins
0527   TFileDirectory dirResMassEta = fs->mkdir("TkTkMassInEtaBins");
0528   InvMassInEtaBins.bookSet(dirResMassEta, hTrackInvMass_);
0529 
0530   // cut flow
0531 
0532   hCutFlow_ = fs->make<TH1F>("hCutFlow", "cut flow;cut step;events left", 6, -0.5, 5.5);
0533   std::string names[6] = {"Total", "Mult.", ">pT", "<eta", "hasVtx", "VtxDist"};
0534   for (unsigned int i = 0; i < 6; i++) {
0535     hCutFlow_->GetXaxis()->SetBinLabel(i + 1, names[i].c_str());
0536   }
0537 
0538   myCounts.zeroAll();
0539 }
0540 
0541 // ------------ method called once each job just after ending the event loop  ------------
0542 void DiMuonVertexValidation::endJob() {
0543   myCounts.printCounts();
0544 
0545   hCutFlow_->SetBinContent(1, myCounts.eventsTotal);
0546   hCutFlow_->SetBinContent(2, myCounts.eventsAfterMult);
0547   hCutFlow_->SetBinContent(3, myCounts.eventsAfterPt);
0548   hCutFlow_->SetBinContent(4, myCounts.eventsAfterEta);
0549   hCutFlow_->SetBinContent(5, myCounts.eventsAfterVtx);
0550   hCutFlow_->SetBinContent(6, myCounts.eventsAfterDist);
0551 
0552   TH1F::SetDefaultSumw2(kTRUE);
0553   const unsigned int nBinsX = hCosPhi_->GetNbinsX();
0554   for (unsigned int i = 1; i <= nBinsX; i++) {
0555     //float binContent = hCosPhi_->GetBinContent(i);
0556     float invertedBinContent = hCosPhi_->GetBinContent(nBinsX + 1 - i);
0557     float invertedBinError = hCosPhi_->GetBinError(nBinsX + 1 - i);
0558     hCosPhiInv_->SetBinContent(i, invertedBinContent);
0559     hCosPhiInv_->SetBinError(i, invertedBinError);
0560 
0561     //float binContent3D = hCosPhi3D_->GetBinContent(i);
0562     float invertedBinContent3D = hCosPhi3D_->GetBinContent(nBinsX + 1 - i);
0563     float invertedBinError3D = hCosPhi3D_->GetBinError(nBinsX + 1 - i);
0564     hCosPhiInv3D_->SetBinContent(i, invertedBinContent3D);
0565     hCosPhiInv3D_->SetBinError(i, invertedBinError3D);
0566   }
0567 }
0568 
0569 // compute the closest vertex to di-lepton ------------------------------------
0570 const reco::Vertex* DiMuonVertexValidation::findClosestVertex(const TransientVertex aTransVtx,
0571                                                               const reco::VertexCollection* vertices) const {
0572   reco::Vertex* defaultVtx = nullptr;
0573 
0574   if (!aTransVtx.isValid())
0575     return defaultVtx;
0576 
0577   // find the closest vertex to the secondary vertex in 3D
0578   VertexDistance3D vertTool3D;
0579   float minD = 9999.;
0580   int closestVtxIndex = 0;
0581   int counter = 0;
0582   for (const auto& vtx : *vertices) {
0583     double dist3D = vertTool3D.distance(aTransVtx, vtx).value();
0584     if (dist3D < minD) {
0585       minD = dist3D;
0586       closestVtxIndex = counter;
0587     }
0588     counter++;
0589   }
0590 
0591   if ((*vertices).at(closestVtxIndex).isValid()) {
0592     return &(vertices->at(closestVtxIndex));
0593   } else {
0594     return defaultVtx;
0595   }
0596 }
0597 
0598 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0599 void DiMuonVertexValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0600   edm::ParameterSetDescription desc;
0601   desc.ifValue(edm::ParameterDescription<bool>("useReco", true, true),
0602                true >> edm::ParameterDescription<edm::InputTag>("muons", edm::InputTag("muons"), true) or
0603                    false >> edm::ParameterDescription<edm::InputTag>(
0604                                 "muonTracks", edm::InputTag("ALCARECOTkAlDiMuon"), true))
0605       ->setComment("If useReco is true need to specify the muon tracks, otherwise take the ALCARECO Inner tracks");
0606   //desc.add<bool>("useReco",true);
0607   //desc.add<edm::InputTag>("muons", edm::InputTag("muons"));
0608   //desc.add<edm::InputTag>("muonTracks", edm::InputTag("ALCARECOTkAlDiMuon"));
0609   desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0610   desc.add<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
0611   desc.add<std::vector<double>>("pTThresholds", {30., 10.});
0612   desc.add<bool>("useClosestVertex", true);
0613   desc.add<double>("maxSVdist", 50.);
0614 
0615   {
0616     edm::ParameterSetDescription psDiMuMass;
0617     psDiMuMass.add<std::string>("name", "DiMuMass");
0618     psDiMuMass.add<std::string>("title", "M(#mu#mu)");
0619     psDiMuMass.add<std::string>("yUnits", "[GeV]");
0620     psDiMuMass.add<int>("NxBins", 24);
0621     psDiMuMass.add<int>("NyBins", 50);
0622     psDiMuMass.add<double>("ymin", 70.);
0623     psDiMuMass.add<double>("ymax", 120.);
0624     desc.add<edm::ParameterSetDescription>("DiMuMassConfig", psDiMuMass);
0625   }
0626   {
0627     edm::ParameterSetDescription psCosPhi;
0628     psCosPhi.add<std::string>("name", "CosPhi");
0629     psCosPhi.add<std::string>("title", "cos(#phi_{xy})");
0630     psCosPhi.add<std::string>("yUnits", "");
0631     psCosPhi.add<int>("NxBins", 50);
0632     psCosPhi.add<int>("NyBins", 50);
0633     psCosPhi.add<double>("ymin", -1.);
0634     psCosPhi.add<double>("ymax", 1.);
0635     desc.add<edm::ParameterSetDescription>("CosPhiConfig", psCosPhi);
0636   }
0637   {
0638     edm::ParameterSetDescription psCosPhi3D;
0639     psCosPhi3D.add<std::string>("name", "CosPhi3D");
0640     psCosPhi3D.add<std::string>("title", "cos(#phi_{3D})");
0641     psCosPhi3D.add<std::string>("yUnits", "");
0642     psCosPhi3D.add<int>("NxBins", 50);
0643     psCosPhi3D.add<int>("NyBins", 50);
0644     psCosPhi3D.add<double>("ymin", -1.);
0645     psCosPhi3D.add<double>("ymax", 1.);
0646     desc.add<edm::ParameterSetDescription>("CosPhi3DConfig", psCosPhi3D);
0647   }
0648   {
0649     edm::ParameterSetDescription psVtxProb;
0650     psVtxProb.add<std::string>("name", "VtxProb");
0651     psVtxProb.add<std::string>("title", "Prob(#chi^{2}_{SV})");
0652     psVtxProb.add<std::string>("yUnits", "");
0653     psVtxProb.add<int>("NxBins", 50);
0654     psVtxProb.add<int>("NyBins", 50);
0655     psVtxProb.add<double>("ymin", 0);
0656     psVtxProb.add<double>("ymax", 1.);
0657     desc.add<edm::ParameterSetDescription>("VtxProbConfig", psVtxProb);
0658   }
0659   {
0660     edm::ParameterSetDescription psVtxDist;
0661     psVtxDist.add<std::string>("name", "VtxDist");
0662     psVtxDist.add<std::string>("title", "d_{xy}(PV,SV)");
0663     psVtxDist.add<std::string>("yUnits", "[#mum]");
0664     psVtxDist.add<int>("NxBins", 50);
0665     psVtxDist.add<int>("NyBins", 100);
0666     psVtxDist.add<double>("ymin", 0);
0667     psVtxDist.add<double>("ymax", 300.);
0668     desc.add<edm::ParameterSetDescription>("VtxDistConfig", psVtxDist);
0669   }
0670   {
0671     edm::ParameterSetDescription psVtxDist3D;
0672     psVtxDist3D.add<std::string>("name", "VtxDist3D");
0673     psVtxDist3D.add<std::string>("title", "d_{3D}(PV,SV)");
0674     psVtxDist3D.add<std::string>("yUnits", "[#mum]");
0675     psVtxDist3D.add<int>("NxBins", 50);
0676     psVtxDist3D.add<int>("NyBins", 250);
0677     psVtxDist3D.add<double>("ymin", 0);
0678     psVtxDist3D.add<double>("ymax", 500.);
0679     desc.add<edm::ParameterSetDescription>("VtxDist3DConfig", psVtxDist3D);
0680   }
0681   {
0682     edm::ParameterSetDescription psVtxDistSig;
0683     psVtxDistSig.add<std::string>("name", "VtxDistSig");
0684     psVtxDistSig.add<std::string>("title", "d_{xy}(PV,SV)/#sigma_{dxy}(PV,SV)");
0685     psVtxDistSig.add<std::string>("yUnits", "");
0686     psVtxDistSig.add<int>("NxBins", 50);
0687     psVtxDistSig.add<int>("NyBins", 100);
0688     psVtxDistSig.add<double>("ymin", 0);
0689     psVtxDistSig.add<double>("ymax", 5.);
0690     desc.add<edm::ParameterSetDescription>("VtxDistSigConfig", psVtxDistSig);
0691   }
0692   {
0693     edm::ParameterSetDescription psVtxDist3DSig;
0694     psVtxDist3DSig.add<std::string>("name", "VtxDist3DSig");
0695     psVtxDist3DSig.add<std::string>("title", "d_{3D}(PV,SV)/#sigma_{d3D}(PV,SV)");
0696     psVtxDist3DSig.add<std::string>("yUnits", "");
0697     psVtxDist3DSig.add<int>("NxBins", 50);
0698     psVtxDist3DSig.add<int>("NyBins", 100);
0699     psVtxDist3DSig.add<double>("ymin", 0);
0700     psVtxDist3DSig.add<double>("ymax", 5.);
0701     desc.add<edm::ParameterSetDescription>("VtxDist3DSigConfig", psVtxDist3DSig);
0702   }
0703 
0704   descriptions.addWithDefaultLabel(desc);
0705 }
0706 
0707 //define this as a plug-in
0708 DEFINE_FWK_MODULE(DiMuonVertexValidation);