OfflineValidation/plugins/SagittaBiasNtuplizer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    Alignment/OfflineValidation
0004 // Class:      SagittaBiasNtuplizer
0005 //
0006 /*
0007  *\class SagittaBiasNtuplizer SagittaBiasNtuplizer.cc Alignment/OfflineValidation/plugins/SagittaBiasNtuplizer.cc
0008
0009  Description: This module is meant to create a simple ntuple to be able to validate the tracker alignment for the presence of Sagitta Bias.
0010
0011  Implementation: the implementation is straightforward.
0012
0013 */
0014 //
0015 // Original Author:  Marco Musich, C. Alexe
0016 //         Created:  Fri, 05 Jan 2023 11:41:00 GMT
0017 //
0018 //
0019
0020 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0021 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0022 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0023 #include "DataFormats/Math/interface/deltaR.h"
0024 #include "DataFormats/MuonReco/interface/Muon.h"
0025 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0026 #include "DataFormats/TrackReco/interface/Track.h"
0027 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0028 #include "DataFormats/VertexReco/interface/Vertex.h"
0029 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/Framework/interface/EventSetup.h"
0032 #include "FWCore/Framework/interface/Frameworkfwd.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034 #include "FWCore/Framework/interface/MakerMacros.h"
0035 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0036 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0037 #include "FWCore/ServiceRegistry/interface/Service.h"
0038
0039 #include "TH1F.h"
0040 #include "TH2I.h"
0041 #include "TTree.h"
0042 #include "TLorentzVector.h"
0043
0044 class SagittaBiasNtuplizer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0045 public:
0046   explicit SagittaBiasNtuplizer(const edm::ParameterSet&);
0047   ~SagittaBiasNtuplizer() override = default;
0048
0049   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0050
0051 private:
0052   void analyze(const edm::Event& event, const edm::EventSetup& setup) override;
0053   double openingAngle(const reco::Track* track1, const reco::Track* track2);
0054   std::pair<unsigned int, reco::Vertex> findClosestVertex(const reco::Track* track1,
0055                                                           const reco::VertexCollection& vertices);
0056   const bool useReco_;
0057   const bool doGen_;
0058
0059   const double muonEtaCut_;
0060   const double muonPtCut_;
0061   const double muondxySigCut_;
0062   const double minMassWindowCut_;
0063   const double maxMassWindowCut_;
0064   const double d0CompatibilityCut_;
0065   const double z0CompatibilityCut_;
0066
0067   std::vector<double> pTthresholds_;
0068
0069   // either on or the other!
0070   //used to select what muon tracks to read from configuration file
0071   edm::EDGetTokenT<reco::MuonCollection> muonsToken_;
0072   //used to select what tracks to read from configuration file
0073   edm::EDGetTokenT<reco::TrackCollection> alcaRecoToken_;
0074
0075   //used to select what tracks to read from configuration file
0076   edm::EDGetTokenT<reco::TrackCollection> tracksToken_;
0077
0078   // for associated genParticles
0079   edm::EDGetTokenT<edm::View<reco::Candidate>> genParticlesToken_;
0080   //edm::EDGetTokenT<std::vector<reco::GenParticle>> genParticlesToken_;
0081
0082   const edm::EDGetTokenT<reco::VertexCollection> vtxToken_;
0083   const edm::EDGetTokenT<reco::BeamSpot> bsToken_;
0084
0085   static double constexpr muMass = 0.1056583745;  //!< Muon mass [GeV]
0086
0087   TTree* tree_;
0088   float mass_;
0089   float posTrackDz_;
0090   float negTrackDz_;
0091   float posTrackD0_;
0092   float negTrackD0_;
0093   float posTrackEta_;
0094   float negTrackEta_;
0095   float posTrackPhi_;
0096   float negTrackPhi_;
0097   float posTrackPt_;
0098   float negTrackPt_;
0099
0100   // for the gen-level info
0101   float genPosMuonDz_;
0102   float genNegMuonDz_;
0103   float genPosMuonD0_;
0104   float genNegMuonD0_;
0105   float genPosMuonEta_;
0106   float genNegMuonEta_;
0107   float genPosMuonPhi_;
0108   float genNegMuonPhi_;
0109   float genPosMuonPt_;
0110   float genNegMuonPt_;
0111
0112   // control histograms
0113   TH2I* h_VertexMatrix;
0114   TH1F* h_cutFlow;
0115   TH1F* h_DeltaD0;
0116   TH1F* h_DeltaDz;
0117   TH1F* h_CosOpeningAngle;
0118 };
0119
0120 //_________________________________________________________________________________
0121 SagittaBiasNtuplizer::SagittaBiasNtuplizer(const edm::ParameterSet& iConfig)
0122     : useReco_(iConfig.getParameter<bool>("useReco")),
0123       doGen_(iConfig.getParameter<bool>("doGen")),
0124       muonEtaCut_(iConfig.getParameter<double>("muonEtaCut")),
0125       muonPtCut_(iConfig.getParameter<double>("muonPtCut")),
0126       muondxySigCut_(iConfig.getParameter<double>("muondxySigCut")),
0127       minMassWindowCut_(iConfig.getParameter<double>("minMassWindowCut")),
0128       maxMassWindowCut_(iConfig.getParameter<double>("maxMassWindowCut")),
0129       d0CompatibilityCut_(iConfig.getParameter<double>("d0CompatibilityCut")),
0130       z0CompatibilityCut_(iConfig.getParameter<double>("z0CompatibilityCut")),
0131       pTthresholds_(iConfig.getParameter<std::vector<double>>("pTThresholds")),
0132       vtxToken_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))),
0133       bsToken_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpot"))),
0134       mass_{0},
0135       posTrackDz_{0},
0136       negTrackDz_{0},
0137       posTrackD0_{0},
0138       negTrackD0_{0},
0139       posTrackEta_{0},
0140       negTrackEta_{0},
0141       posTrackPhi_{0},
0142       negTrackPhi_{0},
0143       posTrackPt_{0},
0144       negTrackPt_{0},
0145       genPosMuonEta_{-99.},
0146       genNegMuonEta_{-99.},
0147       genPosMuonPhi_{-99.},
0148       genNegMuonPhi_{-99.},
0149       genPosMuonPt_{-99.},
0150       genNegMuonPt_{-99.} {
0151   if (useReco_) {
0152     tracksToken_ = mayConsume<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"));
0153     muonsToken_ = mayConsume<reco::MuonCollection>(iConfig.getParameter<edm::InputTag>("muons"));
0154   } else {
0155     alcaRecoToken_ = mayConsume<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("muonTracks"));
0156   }
0157
0158   if (doGen_) {
0159     genParticlesToken_ = consumes<edm::View<reco::Candidate>>(iConfig.getParameter<edm::InputTag>("genParticles"));
0160   }
0161
0162   usesResource(TFileService::kSharedResource);
0163
0164   // sort the vector of thresholds
0165   std::sort(pTthresholds_.begin(), pTthresholds_.end(), [](const double& lhs, const double& rhs) { return lhs > rhs; });
0166
0167   edm::LogInfo("SagittaBiasNtuplizer") << __FUNCTION__;
0168   for (const auto& thr : pTthresholds_) {
0169     edm::LogInfo("SagittaBiasNtuplizer") << " Threshold: " << thr << " ";
0170   }
0171
0172   edm::Service<TFileService> fs;
0173   h_cutFlow = fs->make<TH1F>("cutFlow", "cutFlow;cut;remaining events", 9, -0.5, 8.5);
0174   std::string labels[9] = {"all events",
0175                            "common vertex",
0176                            "d0 cut",
0177                            "p_{T} cut",
0178                            "#eta cut",
0179                            "mass window",
0180                            "#delta d_{0}",
0181                            "#delta d_{z}",
0182                            "opening angle"};
0183   unsigned int count{0};
0184   for (const auto& label : labels) {
0185     count++;
0186     h_cutFlow->GetXaxis()->SetBinLabel(count, label.c_str());
0187   }
0188
0189   h_VertexMatrix = fs->make<TH2I>("VertexMatrix",
0190                                   ";index of closest vertex to #mu_{0};index of closest vertex to #mu_{1}",
0191                                   100,
0192                                   0,
0193                                   100,
0194                                   100,
0195                                   0,
0196                                   100);
0197   h_DeltaD0 = fs->make<TH1F>("DeltaD0", "#Deltad_{0};#Deltad_{0} [cm];events", 100, -0.5, 0.5);
0198   h_DeltaDz = fs->make<TH1F>("DeltaDz", "#Deltad_{z};#Deltad_{z} [cm];events", 100, -1, 1);
0199   h_CosOpeningAngle = fs->make<TH1F>("OpeningAngle", "cos(#gamma(#mu^{+},#mu^{-}));events", 100, -1., 1.);
0200
0201   tree_ = fs->make<TTree>("tree", "My TTree");
0202   tree_->Branch("mass", &mass_, "mass/F");
0203   tree_->Branch("posTrackDz", &posTrackDz_, "posTrackDz/F");
0204   tree_->Branch("negTrackDz", &negTrackDz_, "negTrackDz/F");
0205   tree_->Branch("posTrackD0", &posTrackD0_, "posTrackD0/F");
0206   tree_->Branch("negTrackD0", &negTrackD0_, "negTrackD0/F");
0207   tree_->Branch("posTrackEta", &posTrackEta_, "posTrackEta/F");
0208   tree_->Branch("negTrackEta", &negTrackEta_, "negTrackEta/F");
0209   tree_->Branch("posTrackPhi", &posTrackPhi_, "posTrackPhi/F");
0210   tree_->Branch("negTrackPhi", &negTrackPhi_, "negTrackPhi/F");
0211   tree_->Branch("posTrackPt", &posTrackPt_, "posTrackPt/F");
0212   tree_->Branch("negTrackPt", &negTrackPt_, "negTrackPt/F");
0213
0214   if (doGen_) {
0215     tree_->Branch("genPosMuonDz", &genPosMuonDz_, "genPosMuonDz/F");
0216     tree_->Branch("genNegMuonDz", &genNegMuonDz_, "genNegMuonDz/F");
0217     tree_->Branch("genPosMuonD0", &genPosMuonD0_, "genPosMuonD0/F");
0218     tree_->Branch("genNegMuonD0", &genNegMuonD0_, "genNegMuonD0/F");
0219     tree_->Branch("genPosMuonEta", &genPosMuonEta_, "genPosMuonEta/F");
0220     tree_->Branch("genNegMuonEta", &genNegMuonEta_, "genNegMuonEta/F");
0221     tree_->Branch("genPosMuonPhi", &genPosMuonPhi_, "genPosMuonPhi/F");
0222     tree_->Branch("genNegMuonPhi", &genNegMuonPhi_, "genNegMuonPhi/F");
0223     tree_->Branch("genPosMuonPt", &genPosMuonPt_, "genPosMuonPt/F");
0224     tree_->Branch("genNegMuonPt", &genNegMuonPt_, "genNegMuonPt/F");
0225   }
0226 }
0227
0228 //_________________________________________________________________________________
0229 void SagittaBiasNtuplizer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0230   edm::ParameterSetDescription desc;
0231   desc.ifValue(
0232           edm::ParameterDescription<bool>("useReco", true, true),
0233           true >> edm::ParameterDescription<edm::InputTag>("muons", edm::InputTag("muons"), true) or
0234               false >> edm::ParameterDescription<edm::InputTag>("muonTracks", edm::InputTag("ALCARECOTkAlZMuMu"), true))
0235       ->setComment("If useReco is true need to specify the muon tracks, otherwise take the ALCARECO Inner tracks");
0236   desc.add<bool>("doGen", false);
0237   desc.add<edm::InputTag>("genParticles", edm::InputTag("genParticles"));
0238   desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0239   desc.add<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
0240   desc.add<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
0241   desc.add<double>("muonEtaCut", 2.5)->setComment("muon system acceptance");
0242   desc.add<double>("muonPtCut", 12.)->setComment("in GeV");
0243   desc.add<double>("muondxySigCut", 4.)->setComment("significance of the d0 compatibility with closest vertex");
0244   desc.add<double>("minMassWindowCut", 70.)->setComment("in GeV");
0245   desc.add<double>("maxMassWindowCut", 110.)->setComment("in GeV");
0246   desc.add<double>("d0CompatibilityCut", 0.01)->setComment("d0 compatibility between the two muons");
0247   desc.add<double>("z0CompatibilityCut", 0.06)->setComment("z0 compatibility between the two muons");
0248   desc.add<std::vector<double>>("pTThresholds", {30., 10.});
0249   descriptions.addWithDefaultLabel(desc);
0250 }
0251
0252 //_________________________________________________________________________________
0253 void SagittaBiasNtuplizer::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0254   h_cutFlow->Fill(0);
0255
0256   // the di-muon tracks
0257   std::vector<const reco::Track*> myTracks;
0258
0259   // if we have to start from scratch from RECO data-tier
0260   if (useReco_) {
0261     // select the good muons
0262     std::vector<const reco::Muon*> myGoodMuonVector;
0263     for (const auto& muon : event.get(muonsToken_)) {
0264       const reco::TrackRef t = muon.innerTrack();
0265       if (!t.isNull()) {
0266         if (t->quality(reco::TrackBase::highPurity)) {
0267           if (t->chi2() / t->ndof() <= 2.5 && t->numberOfValidHits() >= 5 &&
0268               t->hitPattern().numberOfValidPixelHits() >= 2 && t->quality(reco::TrackBase::highPurity))
0269             myGoodMuonVector.emplace_back(&muon);
0270         }
0271       }
0272     }
0273
0274     LogDebug("SagittaBiasNtuplizer") << "myGoodMuonVector size: " << myGoodMuonVector.size() << std::endl;
0275     std::sort(myGoodMuonVector.begin(), myGoodMuonVector.end(), [](const reco::Muon*& lhs, const reco::Muon*& rhs) {
0276       return lhs->pt() > rhs->pt();
0277     });
0278
0279     // just check the ordering
0280     for (const auto& muon : myGoodMuonVector) {
0281       LogDebug("SagittaBiasNtuplizer") << "pT: " << muon->pt() << " ";
0282     }
0283     LogDebug("SagittaBiasNtuplizer") << std::endl;
0284
0285     // reject if there's no Z
0286     if (myGoodMuonVector.size() < 2)
0287       return;
0288
0289     if ((myGoodMuonVector[0]->pt()) < pTthresholds_[0] || (myGoodMuonVector[1]->pt() < pTthresholds_[1]))
0290       return;
0291
0292     if (myGoodMuonVector[0]->charge() * myGoodMuonVector[1]->charge() > 0)
0293       return;
0294
0295     //const auto& m1 = myGoodMuonVector[1]->p4();
0296     //const auto& m0 = myGoodMuonVector[0]->p4();
0297     //const auto& mother = m0 + m1;
0298
0299     // just copy the top two muons
0300     std::vector<const reco::Muon*> theZMuonVector;
0301     theZMuonVector.reserve(2);
0302     theZMuonVector.emplace_back(myGoodMuonVector[1]);
0303     theZMuonVector.emplace_back(myGoodMuonVector[0]);
0304
0305     // do the matching of Z muons with inner tracks
0306     unsigned int i = 0;
0307     for (const auto& muon : theZMuonVector) {
0308       i++;
0309       float minD = 1e6;
0310       const reco::Track* theMatch = nullptr;
0311       for (const auto& track : event.get(tracksToken_)) {
0312         float D = ::deltaR2(muon->eta(), muon->phi(), track.eta(), track.phi());
0313         if (D < minD) {
0314           minD = D;
0315           theMatch = &track;
0316         }
0317       }
0318       LogDebug("SagittaBiasNtuplizer") << "pushing new track: " << i << std::endl;
0319       myTracks.emplace_back(theMatch);
0320     }
0321   } else {
0322     // we start directly with the pre-selected ALCARECO tracks
0323     for (const auto& muon : event.get(alcaRecoToken_)) {
0324       myTracks.emplace_back(&muon);
0325     }
0326   }
0327
0328   const reco::VertexCollection& vertices = event.get(vtxToken_);
0329   const reco::BeamSpot& beamSpot = event.get(bsToken_);
0330   math::XYZPoint bs(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
0331
0332   //edm::LogPrint("SagittaBiasNtuplizer") << " track size:" << myTracks.size() << " vertex size:" << vertices.size() << std::endl;
0333
0334   if ((myTracks.size() != 2)) {
0335     LogTrace("SagittaBiasNtuplizer") << "Found " << myTracks.size() << " muons in the event. Skipping";
0336     return;
0337   }
0338
0339   bool passSameVertex{true};
0340   bool passD0sigCut{true};
0341   bool passPtCut{true};
0342   bool passEtaCut{true};
0343   bool passMassWindow{true};
0344   bool passDeltaD0{true};
0345   bool passDeltaDz{true};
0346   bool passOpeningAngle{true};
0347   double d0[2] = {0., 0.};
0348   double dz[2] = {0., 0.};
0349   unsigned int vtxIndex[2] = {999, 999};
0350
0351   unsigned int i = 0;
0352   for (const auto& track : myTracks) {
0353     if (track->pt() < muonPtCut_) {
0354       passPtCut = false;
0355       continue;
0356     }
0357
0358     if (std::abs(track->eta()) > muonEtaCut_) {
0359       passEtaCut = false;
0360       continue;
0361     }
0362
0363     const auto& closestVertex = this->findClosestVertex(track, vertices);
0364     vtxIndex[i] = closestVertex.first;
0365     d0[i] = track->dxy(closestVertex.second.position());
0366     dz[i] = track->dz(closestVertex.second.position());
0367
0368     if (d0[i] / track->dxyError() > muondxySigCut_) {
0369       passD0sigCut = false;
0370       continue;
0371     }
0372
0373     if (track->charge() > 0) {
0374       posTrackDz_ = dz[i];
0375       posTrackD0_ = d0[i];
0376       posTrackEta_ = track->eta();
0377       posTrackPhi_ = track->phi();
0378       posTrackPt_ = track->pt();
0379     } else {
0380       negTrackDz_ = dz[i];
0381       negTrackD0_ = d0[i];
0382       negTrackEta_ = track->eta();
0383       negTrackPhi_ = track->phi();
0384       negTrackPt_ = track->pt();
0385     }
0386     i++;
0387   }
0388
0389   h_VertexMatrix->Fill(vtxIndex[0], vtxIndex[1]);
0390   // check if the two muons have the same vertex
0391   passSameVertex = (vtxIndex[0] == vtxIndex[1]);
0392   if (!passSameVertex)
0393     return;
0394   h_cutFlow->Fill(1);
0395
0396   // checks if both muons pass the IP signficance cut (w.r.t BeamSpot)
0397   if (!passD0sigCut)
0398     return;
0399   h_cutFlow->Fill(2);
0400
0401   // checks if both muons pass the pT cut
0402   if (!passPtCut)
0403     return;
0404   h_cutFlow->Fill(3);
0405
0406   // check if both muons pass the eta cut
0407   if (!passEtaCut)
0408     return;
0409   h_cutFlow->Fill(4);
0410
0411   // compute the invariant mass of the system
0412   TLorentzVector posTrack, negTrack, mother;
0413   posTrack.SetPtEtaPhiM(posTrackPt_, posTrackEta_, posTrackPhi_, muMass);  // assume muon mass for tracks
0414   negTrack.SetPtEtaPhiM(negTrackPt_, negTrackEta_, negTrackPhi_, muMass);
0415   mother = posTrack + negTrack;
0416   mass_ = mother.M();
0417
0418   // checks if invariant mass of the system lies in the fiducial window
0419   passMassWindow = (mass_ > minMassWindowCut_ && mass_ < maxMassWindowCut_);
0420   if (!passMassWindow)
0421     return;
0422   h_cutFlow->Fill(5);
0423
0424   // checks if the di-muon system passes the d0 compatibility cut
0425   passDeltaD0 = (std::abs(d0[0] - d0[1]) < d0CompatibilityCut_);
0426   h_DeltaD0->Fill(d0[0] - d0[1]);
0427   h_DeltaDz->Fill(dz[0] - dz[1]);
0428   if (!passDeltaD0)
0429     return;
0430   h_cutFlow->Fill(6);
0431
0432   // checks if the di-muon system passes the z0 compatibility cut
0433   passDeltaDz = (std::abs(dz[0] - dz[1]) < z0CompatibilityCut_);
0434   if (!passDeltaDz)
0435     return;
0436   h_cutFlow->Fill(7);
0437
0438   // checks if the di-muon system passes the opening angle cut
0439   double openingAngle = this->openingAngle(myTracks[0], myTracks[1]);
0440   h_CosOpeningAngle->Fill(openingAngle);
0441   passOpeningAngle = true;  //(openingAngle > M_PI/4.);
0442
0443   if (!passOpeningAngle)
0444     return;
0445   h_cutFlow->Fill(8);
0446
0447   if (doGen_) {
0448     const edm::View<reco::Candidate>* genPartCollection = &event.get(genParticlesToken_);
0449
0450     // loop on the reconstructed tracks
0451     for (const auto& track : myTracks) {
0452       float drmin = 0.01;
0453       // loop on the gen particles
0454       for (auto g = genPartCollection->begin(); g != genPartCollection->end(); ++g) {
0455         if (g->status() != 1)
0456           continue;
0457
0458         if (std::abs(g->pdgId()) != 13)
0459           continue;
0460
0461         if (g->charge() != track->charge())
0462           continue;
0463
0464         float dR = reco::deltaR2(*g, *track);
0465
0466         auto const& vtx = g->vertex();
0467         auto const& myBeamSpot = beamSpot.position(vtx.z());
0468         const auto& theptinv2 = 1 / g->pt() * g->pt();
0469
0470         if (dR < drmin) {
0471           drmin = dR;
0472
0473           if (g->charge() > 0) {
0474             genPosMuonPt_ = g->pt();
0475             genPosMuonEta_ = g->eta();
0476             genPosMuonPhi_ = g->phi();
0477             //d0
0478             genPosMuonD0_ = -(-(vtx.x() - myBeamSpot.x()) * g->py() + (vtx.y() - myBeamSpot.y()) * g->px()) / g->pt();
0479             //dz
0480             genPosMuonDz_ =
0481                 (vtx.z() - myBeamSpot.z()) -
0482                 ((vtx.x() - myBeamSpot.x()) * g->px() + (vtx.y() - myBeamSpot.y()) * g->py()) * g->pz() * theptinv2;
0483           } else {
0484             genNegMuonPt_ = g->pt();
0485             genNegMuonEta_ = g->eta();
0486             genNegMuonPhi_ = g->phi();
0487             //d0
0488             genNegMuonD0_ = (-(vtx.x() - myBeamSpot.x()) * g->py() + (vtx.y() - myBeamSpot.y()) * g->px()) / g->pt();
0489             //dz
0490             genNegMuonDz_ =
0491                 (vtx.z() - myBeamSpot.z()) -
0492                 ((vtx.x() - myBeamSpot.x()) * g->px() + (vtx.y() - myBeamSpot.y()) * g->py()) * g->pz() * theptinv2;
0493           }
0494         } else {
0495           continue;
0496         }
0497       }
0498     }
0499   }
0500
0501   // if all cuts are passed fill the TTree
0502   tree_->Fill();
0503 }
0504
0505 //_________________________________________________________________________________
0506 double SagittaBiasNtuplizer::openingAngle(const reco::Track* trk1, const reco::Track* trk2) {
0507   math::XYZVector vec1(trk1->px(), trk1->py(), trk1->pz());
0508   math::XYZVector vec2(trk2->px(), trk2->py(), trk2->pz());
0509   return vec1.Dot(vec2) / (vec1.R() * vec2.R());
0510 }
0511
0512 /*
0513 double SagittaBiasNtuplizer::openingAngle(const reco::Track& track1, const reco::Track& track2) {
0514   double dPhi = std::abs(track1.phi() - track2.phi());
0515   double dEta = std::abs(track1.eta() - track2.eta());
0516   double deltaR2 = reco::deltaR2(track1, track2);
0517   double openingAngle = 2 * std::atan2(std::sqrt(std::sin(dEta / 2) * std::sin(dEta / 2) + std::sinh(dPhi / 2) * std::sinh(dPhi / 2)), std::sqrt(1 - deltaR2));
0518   return openingAngle;
0519 }
0520 */
0521
0522 //_________________________________________________________________________________
0523 std::pair<unsigned int, reco::Vertex> SagittaBiasNtuplizer::findClosestVertex(const reco::Track* track,
0524                                                                               const reco::VertexCollection& vertices) {
0525   // Initialize variables for minimum distance and closest vertex
0526   double minDistance = std::numeric_limits<double>::max();
0527   reco::Vertex closestVertex;
0528
0529   unsigned int index{0}, theIndex{999};
0530   // Loop over the primary vertices and calculate the distance to the track
0531   for (const auto& vertex : vertices) {
0532     const math::XYZPoint& vertexPosition = vertex.position();
0533
0534     // Calculate the distance to the track
0535     const auto& trackMomentum = track->momentum();
0536     const auto& vertexToPoint = vertexPosition - track->referencePoint();
0537     double distance = vertexToPoint.Cross(trackMomentum).R() / trackMomentum.R();
0538
0539     // Check if this is the closest vertex so far
0540     if (distance < minDistance) {
0541       minDistance = distance;
0542       closestVertex = vertex;
0543       theIndex = index;
0544     }
0545     index++;
0546   }
0547   return std::make_pair(theIndex, closestVertex);
0548 }
0549
0550 DEFINE_FWK_MODULE(SagittaBiasNtuplizer);