Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-12 23:01:04 

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 // IP tools
0054 #include "TrackingTools/IPTools/interface/IPTools.h"
0055 
0056 // ROOT
0057 #include "TH1F.h"
0058 #include "TH2F.h"
0059 
0060 //#define LogDebug(X) std::cout << X
0061 
0062 //
0063 // class declaration
0064 //
0065 
0066 class DiMuonVertexValidation : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0067 public:
0068   explicit DiMuonVertexValidation(const edm::ParameterSet&);
0069   ~DiMuonVertexValidation() override;
0070 
0071   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0072 
0073 private:
0074   void beginJob() override;
0075   void analyze(const edm::Event&, const edm::EventSetup&) override;
0076   const reco::Vertex* findClosestVertex(const TransientVertex aTransVtx, const reco::VertexCollection* vertices) const;
0077   void endJob() override;
0078 
0079   // ----------member data ---------------------------
0080   DiLeptonHelp::Counts myCounts;
0081 
0082   const std::string motherName_;
0083   const bool useReco_;
0084   const bool useClosestVertex_;
0085   std::pair<float, float> massLimits_; /* for the mass plot x-range */
0086   std::vector<double> pTthresholds_;
0087   const float maxSVdist_;
0088 
0089   // plot configurations
0090   const edm::ParameterSet CosPhiConfiguration_;
0091   const edm::ParameterSet CosPhi3DConfiguration_;
0092   const edm::ParameterSet VtxProbConfiguration_;
0093   const edm::ParameterSet VtxDistConfiguration_;
0094   const edm::ParameterSet VtxDist3DConfiguration_;
0095   const edm::ParameterSet VtxDistSigConfiguration_;
0096   const edm::ParameterSet VtxDist3DSigConfiguration_;
0097   const edm::ParameterSet DiMuMassConfiguration_;
0098 
0099   // control plots
0100   TH1F* hSVProb_;
0101   TH1F* hSVChi2_;
0102   TH1F* hSVNormChi2_;
0103 
0104   TH1F* hSVDist_;
0105   TH1F* hSVDistErr_;
0106   TH1F* hSVDistSig_;
0107   TH1F* hSVCompatibility_;
0108 
0109   TH1F* hSVDist3D_;
0110   TH1F* hSVDist3DErr_;
0111   TH1F* hSVDist3DSig_;
0112   TH1F* hSVCompatibility3D_;
0113 
0114   TH1F* hCosPhi_;
0115   TH1F* hCosPhi3D_;
0116   TH1F* hCosPhiInv_;
0117   TH1F* hCosPhiInv3D_;
0118   TH1F* hCosPhiUnbalance_;
0119   TH1F* hCosPhi3DUnbalance_;
0120 
0121   TH1F* hInvMass_;
0122   TH1F* hTrackInvMass_;
0123 
0124   TH1F* hCutFlow_;
0125 
0126   // impact parameters information
0127   TH1F* hdxy_;
0128   TH1F* hdz_;
0129   TH1F* hdxyErr_;
0130   TH1F* hdzErr_;
0131   TH1F* hIP2d_;
0132   TH1F* hIP3d_;
0133   TH1F* hIP2dsig_;
0134   TH1F* hIP3dsig_;
0135 
0136   // 2D maps
0137 
0138   DiLeptonHelp::PlotsVsKinematics CosPhiPlots = DiLeptonHelp::PlotsVsKinematics(DiLeptonHelp::MM);
0139   DiLeptonHelp::PlotsVsKinematics CosPhi3DPlots = DiLeptonHelp::PlotsVsKinematics(DiLeptonHelp::MM);
0140   DiLeptonHelp::PlotsVsKinematics VtxProbPlots = DiLeptonHelp::PlotsVsKinematics(DiLeptonHelp::MM);
0141   DiLeptonHelp::PlotsVsKinematics VtxDistPlots = DiLeptonHelp::PlotsVsKinematics(DiLeptonHelp::MM);
0142   DiLeptonHelp::PlotsVsKinematics VtxDist3DPlots = DiLeptonHelp::PlotsVsKinematics(DiLeptonHelp::MM);
0143   DiLeptonHelp::PlotsVsKinematics VtxDistSigPlots = DiLeptonHelp::PlotsVsKinematics(DiLeptonHelp::MM);
0144   DiLeptonHelp::PlotsVsKinematics VtxDist3DSigPlots = DiLeptonHelp::PlotsVsKinematics(DiLeptonHelp::MM);
0145   DiLeptonHelp::PlotsVsKinematics ZMassPlots = DiLeptonHelp::PlotsVsKinematics(DiLeptonHelp::MM);
0146 
0147   // plots vs region
0148   DiLeptonHelp::PlotsVsDiLeptonRegion CosPhi3DInEtaBins = DiLeptonHelp::PlotsVsDiLeptonRegion(1.5);
0149   DiLeptonHelp::PlotsVsDiLeptonRegion InvMassInEtaBins = DiLeptonHelp::PlotsVsDiLeptonRegion(1.5);
0150 
0151   const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> ttbESToken_;
0152 
0153   //used to select what tracks to read from configuration file
0154   edm::EDGetTokenT<reco::TrackCollection> tracksToken_;
0155   //used to select what vertices to read from configuration file
0156   const edm::EDGetTokenT<reco::VertexCollection> vertexToken_;
0157 
0158   // either on or the other!
0159   edm::EDGetTokenT<reco::MuonCollection> muonsToken_;      // used to select tracks to read from configuration file
0160   edm::EDGetTokenT<reco::TrackCollection> alcaRecoToken_;  //used to select muon tracks to read from configuration file
0161 };
0162 
0163 //
0164 // constants, enums and typedefs
0165 //
0166 
0167 static constexpr float cmToum = 10e4;
0168 static constexpr float mumass2 = 0.105658367 * 0.105658367;  //mu mass squared (GeV^2/c^4)
0169 
0170 //
0171 // static data member definitions
0172 //
0173 
0174 //
0175 // constructors and destructor
0176 //
0177 DiMuonVertexValidation::DiMuonVertexValidation(const edm::ParameterSet& iConfig)
0178     : motherName_(iConfig.getParameter<std::string>("decayMotherName")),
0179       useReco_(iConfig.getParameter<bool>("useReco")),
0180       useClosestVertex_(iConfig.getParameter<bool>("useClosestVertex")),
0181       pTthresholds_(iConfig.getParameter<std::vector<double>>("pTThresholds")),
0182       maxSVdist_(iConfig.getParameter<double>("maxSVdist")),
0183       CosPhiConfiguration_(iConfig.getParameter<edm::ParameterSet>("CosPhiConfig")),
0184       CosPhi3DConfiguration_(iConfig.getParameter<edm::ParameterSet>("CosPhi3DConfig")),
0185       VtxProbConfiguration_(iConfig.getParameter<edm::ParameterSet>("VtxProbConfig")),
0186       VtxDistConfiguration_(iConfig.getParameter<edm::ParameterSet>("VtxDistConfig")),
0187       VtxDist3DConfiguration_(iConfig.getParameter<edm::ParameterSet>("VtxDist3DConfig")),
0188       VtxDistSigConfiguration_(iConfig.getParameter<edm::ParameterSet>("VtxDistSigConfig")),
0189       VtxDist3DSigConfiguration_(iConfig.getParameter<edm::ParameterSet>("VtxDist3DSigConfig")),
0190       DiMuMassConfiguration_(iConfig.getParameter<edm::ParameterSet>("DiMuMassConfig")),
0191       ttbESToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0192       vertexToken_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))) {
0193   if (useReco_) {
0194     tracksToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"));
0195     muonsToken_ = consumes<reco::MuonCollection>(iConfig.getParameter<edm::InputTag>("muons"));
0196   } else {
0197     alcaRecoToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("muonTracks"));
0198   }
0199 
0200   usesResource(TFileService::kSharedResource);
0201 
0202   // sort the vector of thresholds
0203   std::sort(pTthresholds_.begin(), pTthresholds_.end(), [](const double& lhs, const double& rhs) { return lhs > rhs; });
0204 
0205   edm::LogInfo("DiMuonVertexValidation") << __FUNCTION__;
0206   for (const auto& thr : pTthresholds_) {
0207     edm::LogInfo("DiMuonVertexValidation") << " Threshold: " << thr << " ";
0208   }
0209   edm::LogInfo("DiMuonVertexValidation") << "Max SV distance: " << maxSVdist_ << " ";
0210 
0211   // set the limits for the mass plots
0212   if (motherName_.find('Z') != std::string::npos) {
0213     massLimits_ = std::make_pair(50., 120);
0214   } else if (motherName_.find("J/#psi") != std::string::npos) {
0215     massLimits_ = std::make_pair(2.7, 3.4);
0216   } else if (motherName_.find("#Upsilon") != std::string::npos) {
0217     massLimits_ = std::make_pair(8.9, 9.9);
0218   } else {
0219     edm::LogError("DiMuonVertexValidation") << " unrecognized decay mother particle: " << motherName_
0220                                             << " setting the default for the Z->mm (50.,120.)" << std::endl;
0221     massLimits_ = std::make_pair(50., 120);
0222   }
0223 }
0224 
0225 DiMuonVertexValidation::~DiMuonVertexValidation() = default;
0226 
0227 //
0228 // member functions
0229 //
0230 
0231 // ------------ method called for each event  ------------
0232 void DiMuonVertexValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0233   using namespace edm;
0234 
0235   myCounts.eventsTotal++;
0236 
0237   // the di-muon tracks
0238   std::vector<const reco::Track*> myTracks;
0239 
0240   // if we have to start from scratch from RECO data-tier
0241   if (useReco_) {
0242     // select the good muons
0243     std::vector<const reco::Muon*> myGoodMuonVector;
0244     for (const auto& muon : iEvent.get(muonsToken_)) {
0245       const reco::TrackRef t = muon.innerTrack();
0246       if (!t.isNull()) {
0247         if (t->quality(reco::TrackBase::highPurity)) {
0248           if (t->chi2() / t->ndof() <= 2.5 && t->numberOfValidHits() >= 5 &&
0249               t->hitPattern().numberOfValidPixelHits() >= 2 && t->quality(reco::TrackBase::highPurity))
0250             myGoodMuonVector.emplace_back(&muon);
0251         }
0252       }
0253     }
0254 
0255     LogDebug("DiMuonVertexValidation") << "myGoodMuonVector size: " << myGoodMuonVector.size() << std::endl;
0256     std::sort(myGoodMuonVector.begin(), myGoodMuonVector.end(), [](const reco::Muon*& lhs, const reco::Muon*& rhs) {
0257       return lhs->pt() > rhs->pt();
0258     });
0259 
0260     // just check the ordering
0261     for (const auto& muon : myGoodMuonVector) {
0262       LogDebug("DiMuonVertexValidation") << "pT: " << muon->pt() << " ";
0263     }
0264     LogDebug("DiMuonVertexValidation") << std::endl;
0265 
0266     // reject if there's no Z
0267     if (myGoodMuonVector.size() < 2)
0268       return;
0269 
0270     myCounts.eventsAfterMult++;
0271 
0272     if ((myGoodMuonVector[0]->pt()) < pTthresholds_[0] || (myGoodMuonVector[1]->pt() < pTthresholds_[1]))
0273       return;
0274 
0275     myCounts.eventsAfterPt++;
0276     myCounts.eventsAfterEta++;
0277 
0278     if (myGoodMuonVector[0]->charge() * myGoodMuonVector[1]->charge() > 0)
0279       return;
0280 
0281     const auto& m1 = myGoodMuonVector[1]->p4();
0282     const auto& m0 = myGoodMuonVector[0]->p4();
0283     const auto& mother = m0 + m1;
0284 
0285     float invMass = mother.M();
0286     hInvMass_->Fill(invMass);
0287 
0288     // just copy the top two muons
0289     std::vector<const reco::Muon*> theZMuonVector;
0290     theZMuonVector.reserve(2);
0291     theZMuonVector.emplace_back(myGoodMuonVector[1]);
0292     theZMuonVector.emplace_back(myGoodMuonVector[0]);
0293 
0294     // do the matching of Z muons with inner tracks
0295     unsigned int i = 0;
0296     for (const auto& muon : theZMuonVector) {
0297       i++;
0298       float minD = 1000.;
0299       const reco::Track* theMatch = nullptr;
0300       for (const auto& track : iEvent.get(tracksToken_)) {
0301         float D = ::deltaR(muon->eta(), muon->phi(), track.eta(), track.phi());
0302         if (D < minD) {
0303           minD = D;
0304           theMatch = &track;
0305         }
0306       }
0307       LogDebug("DiMuonVertexValidation") << "pushing new track: " << i << std::endl;
0308       myTracks.emplace_back(theMatch);
0309     }
0310   } else {
0311     // we start directly with the pre-selected ALCARECO tracks
0312     for (const auto& muon : iEvent.get(alcaRecoToken_)) {
0313       myTracks.emplace_back(&muon);
0314     }
0315   }
0316 
0317 #ifdef EDM_ML_DEBUG
0318   for (const auto& track : myTracks) {
0319     edm::LogVerbatim("DiMuonVertexValidation") << __PRETTY_FUNCTION__ << " pT: " << track->pt() << " GeV"
0320                                                << " , pT error: " << track->ptError() << " GeV"
0321                                                << " , eta: " << track->eta() << " , phi: " << track->phi() << std::endl;
0322   }
0323 #endif
0324 
0325   LogDebug("DiMuonVertexValidation") << "selected tracks: " << myTracks.size() << std::endl;
0326 
0327   const TransientTrackBuilder* theB = &iSetup.getData(ttbESToken_);
0328   TransientVertex aTransientVertex;
0329   std::vector<reco::TransientTrack> tks;
0330 
0331   if (myTracks.size() != 2)
0332     return;
0333 
0334   const auto& t1 = myTracks[1]->momentum();
0335   const auto& t0 = myTracks[0]->momentum();
0336   const auto& ditrack = t1 + t0;
0337 
0338   const auto& tplus = myTracks[0]->charge() > 0 ? myTracks[0] : myTracks[1];
0339   const auto& tminus = myTracks[0]->charge() < 0 ? myTracks[0] : myTracks[1];
0340 
0341   TLorentzVector p4_tplus(tplus->px(), tplus->py(), tplus->pz(), sqrt((tplus->p() * tplus->p()) + mumass2));
0342   TLorentzVector p4_tminus(tminus->px(), tminus->py(), tminus->pz(), sqrt((tminus->p() * tminus->p()) + mumass2));
0343 
0344 #ifdef EDM_ML_DEBUG
0345   // Define a lambda function to convert TLorentzVector to a string
0346   auto tLorentzVectorToString = [](const TLorentzVector& vector) {
0347     return std::to_string(vector.Px()) + " " + std::to_string(vector.Py()) + " " + std::to_string(vector.Pz()) + " " +
0348            std::to_string(vector.E());
0349   };
0350 
0351   edm::LogVerbatim("DiMuonVertexValidation") << "mu+" << tLorentzVectorToString(p4_tplus) << std::endl;
0352   edm::LogVerbatim("DiMuonVertexValidation") << "mu-" << tLorentzVectorToString(p4_tminus) << std::endl;
0353 #endif
0354 
0355   const auto& Zp4 = p4_tplus + p4_tminus;
0356   float track_invMass = Zp4.M();
0357   hTrackInvMass_->Fill(track_invMass);
0358 
0359   // creat the pair of TLorentVectors used to make the plos
0360   std::pair<TLorentzVector, TLorentzVector> tktk_p4 = std::make_pair(p4_tplus, p4_tminus);
0361 
0362   // fill the z->mm mass plots
0363   ZMassPlots.fillPlots(track_invMass, tktk_p4);
0364   InvMassInEtaBins.fillTH1Plots(track_invMass, tktk_p4);
0365 
0366   math::XYZPoint ZpT(ditrack.x(), ditrack.y(), 0);
0367   math::XYZPoint Zp(ditrack.x(), ditrack.y(), ditrack.z());
0368 
0369   for (const auto& track : myTracks) {
0370     reco::TransientTrack trajectory = theB->build(track);
0371     tks.push_back(trajectory);
0372   }
0373 
0374   KalmanVertexFitter kalman(true);
0375   aTransientVertex = kalman.vertex(tks);
0376 
0377   double SVProb = TMath::Prob(aTransientVertex.totalChiSquared(), (int)aTransientVertex.degreesOfFreedom());
0378 
0379   LogDebug("DiMuonVertexValidation") << " vertex prob: " << SVProb << std::endl;
0380 
0381   hSVProb_->Fill(SVProb);
0382   hSVChi2_->Fill(aTransientVertex.totalChiSquared());
0383   hSVNormChi2_->Fill(aTransientVertex.totalChiSquared() / (int)aTransientVertex.degreesOfFreedom());
0384 
0385   LogDebug("DiMuonVertexValidation") << " vertex norm chi2: "
0386                                      << (aTransientVertex.totalChiSquared() / (int)aTransientVertex.degreesOfFreedom())
0387                                      << std::endl;
0388 
0389   if (!aTransientVertex.isValid())
0390     return;
0391 
0392   myCounts.eventsAfterVtx++;
0393 
0394   // fill the VtxProb plots
0395   VtxProbPlots.fillPlots(SVProb, tktk_p4);
0396 
0397   math::XYZPoint mainVtxPos(0, 0, 0);
0398   const reco::Vertex* theClosestVertex = nullptr;
0399   // get collection of reconstructed vertices from event
0400   edm::Handle<reco::VertexCollection> vertexHandle = iEvent.getHandle(vertexToken_);
0401   if (vertexHandle.isValid()) {
0402     const reco::VertexCollection* vertices = vertexHandle.product();
0403     theClosestVertex = this->findClosestVertex(aTransientVertex, vertices);
0404   } else {
0405     edm::LogWarning("DiMuonVertexMonitor") << "invalid vertex collection encountered Skipping event!";
0406     return;
0407   }
0408 
0409   reco::Vertex theMainVertex;
0410   if (!useClosestVertex_ || theClosestVertex == nullptr) {
0411     // if the closest vertex is not available, or explicitly not chosen
0412     theMainVertex = vertexHandle.product()->front();
0413   } else {
0414     theMainVertex = *theClosestVertex;
0415   }
0416 
0417   mainVtxPos.SetXYZ(theMainVertex.position().x(), theMainVertex.position().y(), theMainVertex.position().z());
0418   const math::XYZPoint myVertex(
0419       aTransientVertex.position().x(), aTransientVertex.position().y(), aTransientVertex.position().z());
0420   const math::XYZPoint deltaVtx(
0421       mainVtxPos.x() - myVertex.x(), mainVtxPos.y() - myVertex.y(), mainVtxPos.z() - myVertex.z());
0422 
0423 #ifdef EDM_ML_DEBUG
0424   edm::LogVerbatim("DiMuonVertexValidation")
0425       << "mm vertex position:" << aTransientVertex.position().x() << "," << aTransientVertex.position().y() << ","
0426       << aTransientVertex.position().z();
0427 
0428   edm::LogVerbatim("DiMuonVertexValidation") << "main vertex position:" << theMainVertex.position().x() << ","
0429                                              << theMainVertex.position().y() << "," << theMainVertex.position().z();
0430 #endif
0431 
0432   if (theMainVertex.isValid()) {
0433     // fill the impact parameter plots
0434     for (const auto& track : myTracks) {
0435       hdxy_->Fill(track->dxy(mainVtxPos) * cmToum);
0436       hdz_->Fill(track->dz(mainVtxPos) * cmToum);
0437       hdxyErr_->Fill(track->dxyError() * cmToum);
0438       hdzErr_->Fill(track->dzError() * cmToum);
0439 
0440       const auto& ttrk = theB->build(track);
0441       Global3DVector dir(track->px(), track->py(), track->pz());
0442       const auto& ip2d = IPTools::signedTransverseImpactParameter(ttrk, dir, theMainVertex);
0443       const auto& ip3d = IPTools::signedImpactParameter3D(ttrk, dir, theMainVertex);
0444 
0445       hIP2d_->Fill(ip2d.second.value() * cmToum);
0446       hIP3d_->Fill(ip3d.second.value() * cmToum);
0447       hIP2dsig_->Fill(ip2d.second.significance());
0448       hIP3dsig_->Fill(ip3d.second.significance());
0449     }
0450 
0451     LogDebug("DiMuonVertexValidation") << " after filling the IP histograms " << std::endl;
0452 
0453     // Z Vertex distance in the xy plane
0454     VertexDistanceXY vertTool;
0455     double distance = vertTool.distance(aTransientVertex, theMainVertex).value();
0456     double dist_err = vertTool.distance(aTransientVertex, theMainVertex).error();
0457     float compatibility = 0.;
0458 
0459     try {
0460       compatibility = vertTool.compatibility(aTransientVertex, theMainVertex);
0461     } catch (cms::Exception& er) {
0462       LogTrace("DiMuonVertexValidation") << "caught std::exception " << er.what() << std::endl;
0463     }
0464 
0465     hSVCompatibility_->Fill(compatibility);
0466     hSVDist_->Fill(distance * cmToum);
0467     hSVDistErr_->Fill(dist_err * cmToum);
0468     hSVDistSig_->Fill(distance / dist_err);
0469 
0470     // fill the VtxDist plots
0471     VtxDistPlots.fillPlots(distance * cmToum, tktk_p4);
0472 
0473     // fill the VtxDisSig plots
0474     VtxDistSigPlots.fillPlots(distance / dist_err, tktk_p4);
0475 
0476     // Z Vertex distance in 3D
0477 
0478     VertexDistance3D vertTool3D;
0479     double distance3D = vertTool3D.distance(aTransientVertex, theMainVertex).value();
0480     double dist3D_err = vertTool3D.distance(aTransientVertex, theMainVertex).error();
0481     float compatibility3D = 0.;
0482 
0483     try {
0484       compatibility3D = vertTool3D.compatibility(aTransientVertex, theMainVertex);
0485     } catch (cms::Exception& er) {
0486       LogTrace("DiMuonVertexMonitor") << "caught std::exception " << er.what() << std::endl;
0487     }
0488 
0489     hSVCompatibility3D_->Fill(compatibility3D);
0490     hSVDist3D_->Fill(distance3D * cmToum);
0491     hSVDist3DErr_->Fill(dist3D_err * cmToum);
0492     hSVDist3DSig_->Fill(distance3D / dist3D_err);
0493 
0494     // fill the VtxDist3D plots
0495     VtxDist3DPlots.fillPlots(distance3D * cmToum, tktk_p4);
0496 
0497     // fill the VtxDisSig plots
0498     VtxDist3DSigPlots.fillPlots(distance3D / dist3D_err, tktk_p4);
0499 
0500     LogDebug("DiMuonVertexValidation") << "distance: " << distance << "+/-" << dist_err << std::endl;
0501     // cut on the PV - SV distance
0502     if (distance * cmToum < maxSVdist_) {
0503       myCounts.eventsAfterDist++;
0504 
0505       double cosphi = (ZpT.x() * deltaVtx.x() + ZpT.y() * deltaVtx.y()) /
0506                       (sqrt(ZpT.x() * ZpT.x() + ZpT.y() * ZpT.y()) *
0507                        sqrt(deltaVtx.x() * deltaVtx.x() + deltaVtx.y() * deltaVtx.y()));
0508 
0509       double cosphi3D = (Zp.x() * deltaVtx.x() + Zp.y() * deltaVtx.y() + Zp.z() * deltaVtx.z()) /
0510                         (sqrt(Zp.x() * Zp.x() + Zp.y() * Zp.y() + Zp.z() * Zp.z()) *
0511                          sqrt(deltaVtx.x() * deltaVtx.x() + deltaVtx.y() * deltaVtx.y() + deltaVtx.z() * deltaVtx.z()));
0512 
0513       LogDebug("DiMuonVertexValidation") << "cos(phi) = " << cosphi << std::endl;
0514 
0515       hCosPhi_->Fill(cosphi);
0516       hCosPhi3D_->Fill(cosphi3D);
0517 
0518 #ifdef EDM_ML_DEBUG
0519       edm::LogVerbatim("DiMuonVertexValidation")
0520           << "distance " << distance * cmToum << " cosphi3D:" << cosphi3D << std::endl;
0521 #endif
0522 
0523       // unbalance
0524       hCosPhiUnbalance_->Fill(cosphi, 1.);
0525       hCosPhiUnbalance_->Fill(-cosphi, -1.);
0526       hCosPhi3DUnbalance_->Fill(cosphi3D, 1.);
0527       hCosPhi3DUnbalance_->Fill(-cosphi3D, -1.);
0528 
0529       // fill the cosphi plots
0530       CosPhiPlots.fillPlots(cosphi, tktk_p4);
0531 
0532       // fill the cosphi3D plots
0533       CosPhi3DPlots.fillPlots(cosphi3D, tktk_p4);
0534 
0535       // fill the cosphi3D plots in eta bins
0536       CosPhi3DInEtaBins.fillTH1Plots(cosphi3D, tktk_p4);
0537     }
0538   }
0539 }
0540 
0541 // ------------ method called once each job just before starting event loop  ------------
0542 void DiMuonVertexValidation::beginJob() {
0543   edm::Service<TFileService> fs;
0544 
0545   // clang-format off
0546   TH1F::SetDefaultSumw2(kTRUE);
0547   hSVProb_ = fs->make<TH1F>("VtxProb", ";#mu^{+}#mu^{-} vertex probability;N(#mu#mu pairs)", 100, 0., 1.);
0548 
0549   auto extractRangeValues = [](const edm::ParameterSet& PSetConfiguration_) -> std::pair<double, double> {
0550     double min = PSetConfiguration_.getParameter<double>("ymin");
0551     double max = PSetConfiguration_.getParameter<double>("ymax");
0552     return {min, max};
0553   };
0554 
0555   std::string ts = fmt::sprintf(";%s vertex probability;N(#mu#mu pairs)", motherName_);
0556   std::string ps = "N(#mu#mu pairs)";
0557   ts = fmt::sprintf("#chi^{2} of the %s vertex; #chi^{2} of the %s vertex; %s", motherName_, motherName_, ps);
0558   hSVChi2_ = fs->make<TH1F>("VtxChi2", ts.c_str(), 200, 0., 200.);
0559 
0560   ts = fmt::sprintf("#chi^{2}/ndf of the %s vertex; #chi^{2}/ndf of %s vertex; %s", motherName_, motherName_, ps);
0561   hSVNormChi2_ = fs->make<TH1F>("VtxNormChi2", ts.c_str(), 100, 0., 20.);
0562 
0563   // take the range from the 2D histograms
0564   const auto& svDistRng = extractRangeValues(VtxDistConfiguration_);
0565   hSVDist_ = fs->make<TH1F>("VtxDist", ";PV-#mu^{+}#mu^{-} vertex xy distance [#mum];N(#mu#mu pairs)", 100, svDistRng.first, svDistRng.second);
0566 
0567   std::string histTit = motherName_ + " #rightarrow #mu^{+}#mu^{-}";
0568   ts = fmt::sprintf("%s;PV-%sV xy distance error [#mum];%s", histTit, motherName_, ps);
0569   hSVDistErr_ = fs->make<TH1F>("VtxDistErr", ts.c_str(), 100, 0., 1000.);
0570 
0571   // take the range from the 2D histograms
0572   const auto& svDistSigRng = extractRangeValues(VtxDistSigConfiguration_);
0573   hSVDistSig_ = fs->make<TH1F>("VtxDistSig", ";PV-#mu^{+}#mu^{-} vertex xy distance signficance;N(#mu#mu pairs)", 100, svDistSigRng.first, svDistSigRng.second);
0574 
0575   // take the range from the 2D histograms
0576   const auto& svDist3DRng = extractRangeValues(VtxDist3DConfiguration_);
0577   hSVDist3D_ = fs->make<TH1F>("VtxDist3D", ";PV-#mu^{+}#mu^{-} vertex 3D distance [#mum];N(#mu#mu pairs)", 100, svDist3DRng.first, svDist3DRng.second);
0578 
0579   ts = fmt::sprintf("%s;PV-%sV 3D distance error [#mum];%s", histTit, motherName_, ps);
0580   hSVDist3DErr_ = fs->make<TH1F>("VtxDist3DErr", ts.c_str(), 100, 0., 1000.);
0581 
0582   // take the range from the 2D histograms
0583   const auto& svDist3DSigRng = extractRangeValues(VtxDist3DSigConfiguration_);
0584   hSVDist3DSig_ = fs->make<TH1F>("VtxDist3DSig", ";PV-#mu^{+}#mu^{-} vertex 3D distance signficance;N(#mu#mu pairs)", 100, svDist3DSigRng.first, svDist3DSigRng.second);
0585 
0586   ts = fmt::sprintf("compatibility of %s vertex; compatibility of %s vertex; %s", motherName_, motherName_, ps);
0587   hSVCompatibility_ = fs->make<TH1F>("VtxCompatibility", ts.c_str(), 100, 0., 100.);
0588 
0589   ts = fmt::sprintf("3D compatibility of %s vertex;3D compatibility of %s vertex; %s", motherName_, motherName_, ps);
0590   hSVCompatibility3D_ = fs->make<TH1F>("VtxCompatibility3D", ts.c_str(), 100, 0., 100.);
0591 
0592   // take the range from the 2D histograms
0593   const auto& massRng = extractRangeValues(DiMuMassConfiguration_);
0594   hInvMass_ = fs->make<TH1F>("InvMass", ";M(#mu#mu) [GeV];N(#mu#mu pairs)", 70., massRng.first, massRng.second);
0595   hTrackInvMass_ = fs->make<TH1F>("TkTkInvMass", ";M(tk,tk) [GeV];N(tk tk pairs)", 70., massRng.first, massRng.second);
0596 
0597   hCosPhi_ = fs->make<TH1F>("CosPhi", ";cos(#phi_{xy});N(#mu#mu pairs)", 50, -1., 1.);
0598   hCosPhi3D_ = fs->make<TH1F>("CosPhi3D", ";cos(#phi_{3D});N(#mu#mu pairs)", 50, -1., 1.);
0599 
0600   hCosPhiInv_ = fs->make<TH1F>("CosPhiInv", ";inverted cos(#phi_{xy});N(#mu#mu pairs)", 50, -1., 1.);
0601   hCosPhiInv3D_ = fs->make<TH1F>("CosPhiInv3D", ";inverted cos(#phi_{3D});N(#mu#mu pairs)", 50, -1., 1.);
0602 
0603   hCosPhiUnbalance_ = fs->make<TH1F>("CosPhiUnbalance", fmt::sprintf("%s;cos(#phi_{xy}) unbalance;#Delta%s", histTit, ps).c_str(), 50, -1.,1.);
0604   hCosPhi3DUnbalance_ = fs->make<TH1F>("CosPhi3DUnbalance", fmt::sprintf("%s;cos(#phi_{3D}) unbalance;#Delta%s", histTit, ps).c_str(), 50, -1., 1.);
0605 
0606   hdxy_ = fs->make<TH1F>("dxy", fmt::sprintf("%s;muon track d_{xy}(PV) [#mum];muon tracks", histTit).c_str(), 150, -300, 300);
0607   hdz_ = fs->make<TH1F>("dz", fmt::sprintf("%s;muon track d_{z}(PV) [#mum];muon tracks", histTit).c_str(), 150, -300, 300);
0608   hdxyErr_ = fs->make<TH1F>("dxyErr", fmt::sprintf("%s;muon track err_{dxy} [#mum];muon tracks", histTit).c_str(), 250, 0., 500.);
0609   hdzErr_ = fs->make<TH1F>("dzErr", fmt::sprintf("%s;muon track err_{dz} [#mum];muon tracks", histTit).c_str(), 250, 0., 500.);
0610   hIP2d_ = fs->make<TH1F>("IP2d", fmt::sprintf("%s;muon track IP_{2D} [#mum];muon tracks", histTit).c_str(), 150, -300, 300);
0611   hIP3d_ = fs->make<TH1F>("IP3d", fmt::sprintf("%s;muon track IP_{3D} [#mum];muon tracks", histTit).c_str(), 150, -300, 300);
0612   hIP2dsig_ = fs->make<TH1F>("IP2Dsig", fmt::sprintf("%s;muon track IP_{2D} significance;muon tracks", histTit).c_str(), 100, 0., 5.);
0613   hIP3dsig_ = fs->make<TH1F>("IP3Dsig", fmt::sprintf("%s;muon track IP_{3D} significance;muon tracks", histTit).c_str(), 100, 0., 5.);
0614   // clang-format on
0615 
0616   // 2D Maps
0617 
0618   TFileDirectory dirCosPhi = fs->mkdir("CosPhiPlots");
0619   CosPhiPlots.bookFromPSet(dirCosPhi, CosPhiConfiguration_);
0620 
0621   TFileDirectory dirCosPhi3D = fs->mkdir("CosPhi3DPlots");
0622   CosPhi3DPlots.bookFromPSet(dirCosPhi3D, CosPhi3DConfiguration_);
0623 
0624   TFileDirectory dirVtxProb = fs->mkdir("VtxProbPlots");
0625   VtxProbPlots.bookFromPSet(dirVtxProb, VtxProbConfiguration_);
0626 
0627   TFileDirectory dirVtxDist = fs->mkdir("VtxDistPlots");
0628   VtxDistPlots.bookFromPSet(dirVtxDist, VtxDistConfiguration_);
0629 
0630   TFileDirectory dirVtxDist3D = fs->mkdir("VtxDist3DPlots");
0631   VtxDist3DPlots.bookFromPSet(dirVtxDist3D, VtxDist3DConfiguration_);
0632 
0633   TFileDirectory dirVtxDistSig = fs->mkdir("VtxDistSigPlots");
0634   VtxDistSigPlots.bookFromPSet(dirVtxDistSig, VtxDistSigConfiguration_);
0635 
0636   TFileDirectory dirVtxDist3DSig = fs->mkdir("VtxDist3DSigPlots");
0637   VtxDist3DSigPlots.bookFromPSet(dirVtxDist3DSig, VtxDist3DSigConfiguration_);
0638 
0639   TFileDirectory dirInvariantMass = fs->mkdir("InvariantMassPlots");
0640   ZMassPlots.bookFromPSet(dirInvariantMass, DiMuMassConfiguration_);
0641 
0642   // CosPhi3D in eta bins
0643   TFileDirectory dirCosphi3DEta = fs->mkdir("CosPhi3DInEtaBins");
0644   CosPhi3DInEtaBins.bookSet(dirCosphi3DEta, hCosPhi3D_);
0645 
0646   // Z-> mm mass in eta bins
0647   TFileDirectory dirResMassEta = fs->mkdir("TkTkMassInEtaBins");
0648   InvMassInEtaBins.bookSet(dirResMassEta, hTrackInvMass_);
0649 
0650   // cut flow
0651 
0652   hCutFlow_ = fs->make<TH1F>("hCutFlow", "cut flow;cut step;events left", 6, -0.5, 5.5);
0653   std::string names[6] = {"Total", "Mult.", ">pT", "<eta", "hasVtx", "VtxDist"};
0654   for (unsigned int i = 0; i < 6; i++) {
0655     hCutFlow_->GetXaxis()->SetBinLabel(i + 1, names[i].c_str());
0656   }
0657 
0658   myCounts.zeroAll();
0659 }
0660 
0661 // ------------ method called once each job just after ending the event loop  ------------
0662 void DiMuonVertexValidation::endJob() {
0663   myCounts.printCounts();
0664 
0665   hCutFlow_->SetBinContent(1, myCounts.eventsTotal);
0666   hCutFlow_->SetBinContent(2, myCounts.eventsAfterMult);
0667   hCutFlow_->SetBinContent(3, myCounts.eventsAfterPt);
0668   hCutFlow_->SetBinContent(4, myCounts.eventsAfterEta);
0669   hCutFlow_->SetBinContent(5, myCounts.eventsAfterVtx);
0670   hCutFlow_->SetBinContent(6, myCounts.eventsAfterDist);
0671 
0672   TH1F::SetDefaultSumw2(kTRUE);
0673   const unsigned int nBinsX = hCosPhi_->GetNbinsX();
0674   for (unsigned int i = 1; i <= nBinsX; i++) {
0675     //float binContent = hCosPhi_->GetBinContent(i);
0676     float invertedBinContent = hCosPhi_->GetBinContent(nBinsX + 1 - i);
0677     float invertedBinError = hCosPhi_->GetBinError(nBinsX + 1 - i);
0678     hCosPhiInv_->SetBinContent(i, invertedBinContent);
0679     hCosPhiInv_->SetBinError(i, invertedBinError);
0680 
0681     //float binContent3D = hCosPhi3D_->GetBinContent(i);
0682     float invertedBinContent3D = hCosPhi3D_->GetBinContent(nBinsX + 1 - i);
0683     float invertedBinError3D = hCosPhi3D_->GetBinError(nBinsX + 1 - i);
0684     hCosPhiInv3D_->SetBinContent(i, invertedBinContent3D);
0685     hCosPhiInv3D_->SetBinError(i, invertedBinError3D);
0686   }
0687 }
0688 
0689 // compute the closest vertex to di-lepton ------------------------------------
0690 const reco::Vertex* DiMuonVertexValidation::findClosestVertex(const TransientVertex aTransVtx,
0691                                                               const reco::VertexCollection* vertices) const {
0692   reco::Vertex* defaultVtx = nullptr;
0693 
0694   if (!aTransVtx.isValid())
0695     return defaultVtx;
0696 
0697   // find the closest vertex to the secondary vertex in 3D
0698   VertexDistance3D vertTool3D;
0699   float minD = 9999.;
0700   int closestVtxIndex = 0;
0701   int counter = 0;
0702   for (const auto& vtx : *vertices) {
0703     double dist3D = vertTool3D.distance(aTransVtx, vtx).value();
0704     if (dist3D < minD) {
0705       minD = dist3D;
0706       closestVtxIndex = counter;
0707     }
0708     counter++;
0709   }
0710 
0711   if ((*vertices).at(closestVtxIndex).isValid()) {
0712     return &(vertices->at(closestVtxIndex));
0713   } else {
0714     return defaultVtx;
0715   }
0716 }
0717 
0718 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0719 void DiMuonVertexValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0720   edm::ParameterSetDescription desc;
0721   desc.ifValue(edm::ParameterDescription<bool>("useReco", true, true),
0722                true >> edm::ParameterDescription<edm::InputTag>("muons", edm::InputTag("muons"), true) or
0723                    false >> edm::ParameterDescription<edm::InputTag>(
0724                                 "muonTracks", edm::InputTag("ALCARECOTkAlDiMuon"), true))
0725       ->setComment("If useReco is true need to specify the muon tracks, otherwise take the ALCARECO Inner tracks");
0726   //desc.add<bool>("useReco",true);
0727   //desc.add<edm::InputTag>("muons", edm::InputTag("muons"));
0728   //desc.add<edm::InputTag>("muonTracks", edm::InputTag("ALCARECOTkAlDiMuon"));
0729   desc.add<std::string>("decayMotherName", "Z");
0730   desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0731   desc.add<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
0732   desc.add<std::vector<double>>("pTThresholds", {30., 10.});
0733   desc.add<bool>("useClosestVertex", true);
0734   desc.add<double>("maxSVdist", 50.);
0735 
0736   {
0737     edm::ParameterSetDescription psDiMuMass;
0738     psDiMuMass.add<std::string>("name", "DiMuMass");
0739     psDiMuMass.add<std::string>("title", "M(#mu#mu)");
0740     psDiMuMass.add<std::string>("yUnits", "[GeV]");
0741     psDiMuMass.add<int>("NxBins", 24);
0742     psDiMuMass.add<int>("NyBins", 50);
0743     psDiMuMass.add<double>("ymin", 70.);
0744     psDiMuMass.add<double>("ymax", 120.);
0745     desc.add<edm::ParameterSetDescription>("DiMuMassConfig", psDiMuMass);
0746   }
0747   {
0748     edm::ParameterSetDescription psCosPhi;
0749     psCosPhi.add<std::string>("name", "CosPhi");
0750     psCosPhi.add<std::string>("title", "cos(#phi_{xy})");
0751     psCosPhi.add<std::string>("yUnits", "");
0752     psCosPhi.add<int>("NxBins", 50);
0753     psCosPhi.add<int>("NyBins", 50);
0754     psCosPhi.add<double>("ymin", -1.);
0755     psCosPhi.add<double>("ymax", 1.);
0756     desc.add<edm::ParameterSetDescription>("CosPhiConfig", psCosPhi);
0757   }
0758   {
0759     edm::ParameterSetDescription psCosPhi3D;
0760     psCosPhi3D.add<std::string>("name", "CosPhi3D");
0761     psCosPhi3D.add<std::string>("title", "cos(#phi_{3D})");
0762     psCosPhi3D.add<std::string>("yUnits", "");
0763     psCosPhi3D.add<int>("NxBins", 50);
0764     psCosPhi3D.add<int>("NyBins", 50);
0765     psCosPhi3D.add<double>("ymin", -1.);
0766     psCosPhi3D.add<double>("ymax", 1.);
0767     desc.add<edm::ParameterSetDescription>("CosPhi3DConfig", psCosPhi3D);
0768   }
0769   {
0770     edm::ParameterSetDescription psVtxProb;
0771     psVtxProb.add<std::string>("name", "VtxProb");
0772     psVtxProb.add<std::string>("title", "Prob(#chi^{2}_{SV})");
0773     psVtxProb.add<std::string>("yUnits", "");
0774     psVtxProb.add<int>("NxBins", 50);
0775     psVtxProb.add<int>("NyBins", 50);
0776     psVtxProb.add<double>("ymin", 0);
0777     psVtxProb.add<double>("ymax", 1.);
0778     desc.add<edm::ParameterSetDescription>("VtxProbConfig", psVtxProb);
0779   }
0780   {
0781     edm::ParameterSetDescription psVtxDist;
0782     psVtxDist.add<std::string>("name", "VtxDist");
0783     psVtxDist.add<std::string>("title", "d_{xy}(PV,SV)");
0784     psVtxDist.add<std::string>("yUnits", "[#mum]");
0785     psVtxDist.add<int>("NxBins", 50);
0786     psVtxDist.add<int>("NyBins", 100);
0787     psVtxDist.add<double>("ymin", 0);
0788     psVtxDist.add<double>("ymax", 300.);
0789     desc.add<edm::ParameterSetDescription>("VtxDistConfig", psVtxDist);
0790   }
0791   {
0792     edm::ParameterSetDescription psVtxDist3D;
0793     psVtxDist3D.add<std::string>("name", "VtxDist3D");
0794     psVtxDist3D.add<std::string>("title", "d_{3D}(PV,SV)");
0795     psVtxDist3D.add<std::string>("yUnits", "[#mum]");
0796     psVtxDist3D.add<int>("NxBins", 50);
0797     psVtxDist3D.add<int>("NyBins", 250);
0798     psVtxDist3D.add<double>("ymin", 0);
0799     psVtxDist3D.add<double>("ymax", 500.);
0800     desc.add<edm::ParameterSetDescription>("VtxDist3DConfig", psVtxDist3D);
0801   }
0802   {
0803     edm::ParameterSetDescription psVtxDistSig;
0804     psVtxDistSig.add<std::string>("name", "VtxDistSig");
0805     psVtxDistSig.add<std::string>("title", "d_{xy}(PV,SV)/#sigma_{dxy}(PV,SV)");
0806     psVtxDistSig.add<std::string>("yUnits", "");
0807     psVtxDistSig.add<int>("NxBins", 50);
0808     psVtxDistSig.add<int>("NyBins", 100);
0809     psVtxDistSig.add<double>("ymin", 0);
0810     psVtxDistSig.add<double>("ymax", 5.);
0811     desc.add<edm::ParameterSetDescription>("VtxDistSigConfig", psVtxDistSig);
0812   }
0813   {
0814     edm::ParameterSetDescription psVtxDist3DSig;
0815     psVtxDist3DSig.add<std::string>("name", "VtxDist3DSig");
0816     psVtxDist3DSig.add<std::string>("title", "d_{3D}(PV,SV)/#sigma_{d3D}(PV,SV)");
0817     psVtxDist3DSig.add<std::string>("yUnits", "");
0818     psVtxDist3DSig.add<int>("NxBins", 50);
0819     psVtxDist3DSig.add<int>("NyBins", 100);
0820     psVtxDist3DSig.add<double>("ymin", 0);
0821     psVtxDist3DSig.add<double>("ymax", 5.);
0822     desc.add<edm::ParameterSetDescription>("VtxDist3DSigConfig", psVtxDist3DSig);
0823   }
0824 
0825   descriptions.addWithDefaultLabel(desc);
0826 }
0827 
0828 //define this as a plug-in
0829 DEFINE_FWK_MODULE(DiMuonVertexValidation);
[ 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!