Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-05-09 22:40:13 

0001 ////////////////////////////// BToTrkTrkLLBuilder //////////////////////////////
0002 /// original authors: G Karathanasis (CERN),  G Melachroinos (NKUA)
0003 // takes the ditrack collection and a dilepton collection and produces B moth
0004 // - ers using a four-track vertex
0005 
0006 #include <algorithm>
0007 #include <limits>
0008 #include <map>
0009 #include <memory>
0010 #include <string>
0011 #include <vector>
0012 
0013 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0014 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0015 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0016 #include "DataFormats/Candidate/interface/Candidate.h"
0017 #include "DataFormats/Math/interface/deltaR.h"
0018 #include "DataFormats/PatCandidates/interface/CompositeCandidate.h"
0019 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0020 #include "FWCore/Framework/interface/Event.h"
0021 #include "FWCore/Framework/interface/global/EDProducer.h"
0022 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0025 #include "FWCore/Utilities/interface/InputTag.h"
0026 #include "KinVtxFitter.h"
0027 #include "MagneticField/Engine/interface/MagneticField.h"
0028 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0029 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0030 #include "helper.h"
0031 
0032 class BToTrkTrkLLBuilder : public edm::global::EDProducer<> {
0033   // perhaps we need better structure here (begin run etc)
0034 public:
0035   typedef std::vector<reco::TransientTrack> TransientTrackCollection;
0036 
0037   explicit BToTrkTrkLLBuilder(const edm::ParameterSet &cfg)
0038       : bFieldToken_{esConsumes<MagneticField, IdealMagneticFieldRecord>()},
0039         // selections
0040         pre_vtx_selection_{cfg.getParameter<std::string>("preVtxSelection")},
0041         post_vtx_selection_{cfg.getParameter<std::string>("postVtxSelection")},
0042         // inputs
0043         dileptons_{consumes<pat::CompositeCandidateCollection>(cfg.getParameter<edm::InputTag>("dileptons"))},
0044         ditracks_{consumes<pat::CompositeCandidateCollection>(cfg.getParameter<edm::InputTag>("ditracks"))},
0045         leptons_ttracks_{consumes<TransientTrackCollection>(cfg.getParameter<edm::InputTag>("leptonTransientTracks"))},
0046         ditracks_ttracks_{consumes<TransientTrackCollection>(cfg.getParameter<edm::InputTag>("transientTracks"))},
0047         pu_tracks_(consumes<pat::CompositeCandidateCollection>(cfg.getParameter<edm::InputTag>("PUtracks"))),
0048         beamspot_{consumes<reco::BeamSpot>(cfg.getParameter<edm::InputTag>("beamSpot"))},
0049         dilepton_constraint_{cfg.getParameter<bool>("dileptonMassContraint")} {
0050     // output
0051     produces<pat::CompositeCandidateCollection>("SelectedBToTrkTrkMuMu");
0052     produces<pat::CompositeCandidateCollection>("SelectedTrkTrk");
0053   }
0054 
0055   ~BToTrkTrkLLBuilder() override {}
0056 
0057   void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override;
0058 
0059 private:
0060   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> bFieldToken_;
0061 
0062   // selections
0063   const StringCutObjectSelector<pat::CompositeCandidate> pre_vtx_selection_;
0064   const StringCutObjectSelector<pat::CompositeCandidate> post_vtx_selection_;
0065 
0066   // inputs
0067   const edm::EDGetTokenT<pat::CompositeCandidateCollection> dileptons_;
0068   const edm::EDGetTokenT<pat::CompositeCandidateCollection> ditracks_;
0069   const edm::EDGetTokenT<TransientTrackCollection> leptons_ttracks_;
0070   const edm::EDGetTokenT<TransientTrackCollection> ditracks_ttracks_;
0071   const edm::EDGetTokenT<pat::CompositeCandidateCollection> pu_tracks_;
0072   const edm::EDGetTokenT<reco::BeamSpot> beamspot_;
0073   const bool dilepton_constraint_;
0074 };
0075 
0076 void BToTrkTrkLLBuilder::produce(edm::StreamID, edm::Event &evt, edm::EventSetup const &iSetup) const {
0077   // input
0078   edm::Handle<pat::CompositeCandidateCollection> dileptons;
0079   evt.getByToken(dileptons_, dileptons);
0080   edm::Handle<TransientTrackCollection> leptons_ttracks;
0081   evt.getByToken(leptons_ttracks_, leptons_ttracks);
0082 
0083   edm::Handle<pat::CompositeCandidateCollection> ditracks;
0084   evt.getByToken(ditracks_, ditracks);
0085   edm::Handle<TransientTrackCollection> ditracks_ttracks;
0086   evt.getByToken(ditracks_ttracks_, ditracks_ttracks);
0087 
0088   edm::Handle<pat::CompositeCandidateCollection> pu_tracks;
0089   evt.getByToken(pu_tracks_, pu_tracks);
0090 
0091   edm::Handle<reco::BeamSpot> beamspot;
0092   evt.getByToken(beamspot_, beamspot);
0093 
0094   edm::ESHandle<MagneticField> fieldHandle;
0095   const auto &bField = iSetup.getData(bFieldToken_);
0096   AnalyticalImpactPointExtrapolator extrapolator(&bField);
0097 
0098   // output
0099   std::unique_ptr<pat::CompositeCandidateCollection> ret_val(new pat::CompositeCandidateCollection());
0100   std::unique_ptr<pat::CompositeCandidateCollection> ditrack_out(new pat::CompositeCandidateCollection());
0101 
0102   for (size_t ditracks_idx = 0; ditracks_idx < ditracks->size(); ++ditracks_idx) {
0103     // both k*,phi and lep pair already passed cuts; no need for more
0104     // preselection
0105     edm::Ptr<pat::CompositeCandidate> ditracks_ptr(ditracks, ditracks_idx);
0106     edm::Ptr<reco::Candidate> trk1_ptr = ditracks_ptr->userCand("trk1");
0107     edm::Ptr<reco::Candidate> trk2_ptr = ditracks_ptr->userCand("trk2");
0108     int trk1_idx = ditracks_ptr->userInt("trk1_idx");
0109     int trk2_idx = ditracks_ptr->userInt("trk2_idx");
0110 
0111     for (size_t ll_idx = 0; ll_idx < dileptons->size(); ++ll_idx) {
0112       edm::Ptr<pat::CompositeCandidate> ll_ptr(dileptons, ll_idx);
0113       edm::Ptr<reco::Candidate> l1_ptr = ll_ptr->userCand("l1");
0114       edm::Ptr<reco::Candidate> l2_ptr = ll_ptr->userCand("l2");
0115       int l1_idx = ll_ptr->userInt("l1_idx");
0116       int l2_idx = ll_ptr->userInt("l2_idx");
0117 
0118       // B0 candidate
0119       pat::CompositeCandidate cand;
0120       cand.setP4(ll_ptr->p4() + ditracks_ptr->p4());
0121       cand.setCharge(l1_ptr->charge() + l2_ptr->charge() + trk1_ptr->charge() + trk2_ptr->charge());
0122 
0123       // save daughters - unfitted
0124       cand.addUserCand("l1", l1_ptr);
0125       cand.addUserCand("l2", l2_ptr);
0126       cand.addUserCand("trk1", trk1_ptr);
0127       cand.addUserCand("trk2", trk2_ptr);
0128       cand.addUserCand("ditrack", ditracks_ptr);
0129       cand.addUserCand("dilepton", ll_ptr);
0130 
0131       // save indices
0132       cand.addUserInt("l1_idx", l1_idx);
0133       cand.addUserInt("l2_idx", l2_idx);
0134       cand.addUserInt("trk1_idx", trk1_idx);
0135       cand.addUserInt("trk2_idx", trk2_idx);
0136       cand.addUserInt("ditrack_idx", ditracks_idx);
0137 
0138       auto lep1_p4 = l1_ptr->polarP4();
0139       auto lep2_p4 = l2_ptr->polarP4();
0140       lep1_p4.SetM(l1_ptr->mass());
0141       lep2_p4.SetM(l2_ptr->mass());
0142 
0143       auto trk1_p4 = trk1_ptr->polarP4();
0144       auto trk2_p4 = trk2_ptr->polarP4();
0145 
0146       trk1_p4.SetM(bph::K_MASS);
0147       trk2_p4.SetM(bph::K_MASS);
0148       cand.addUserFloat("unfitted_B_mass_KK", (trk1_p4 + trk2_p4 + lep1_p4 + lep2_p4).M());
0149       trk1_p4.SetM(bph::K_MASS);
0150       trk2_p4.SetM(bph::PI_MASS);
0151       cand.addUserFloat("unfitted_B_mass_Kpi", (trk1_p4 + trk2_p4 + lep1_p4 + lep2_p4).M());
0152       trk2_p4.SetM(bph::K_MASS);
0153       trk1_p4.SetM(bph::PI_MASS);
0154       cand.addUserFloat("unfitted_B_mass_piK", (trk1_p4 + trk2_p4 + lep1_p4 + lep2_p4).M());
0155 
0156       auto dr_info = bph::min_max_dr({l1_ptr, l2_ptr, trk1_ptr, trk2_ptr});
0157       cand.addUserFloat("min_dr", dr_info.first);
0158       cand.addUserFloat("max_dr", dr_info.second);
0159 
0160       // check if pass pre vertex cut
0161       if (!pre_vtx_selection_(cand))
0162         continue;
0163 
0164       KinVtxFitter fitter({leptons_ttracks->at(l1_idx),
0165                            leptons_ttracks->at(l2_idx),
0166                            ditracks_ttracks->at(trk1_idx),
0167                            ditracks_ttracks->at(trk2_idx)},
0168                           {l1_ptr->mass(), l2_ptr->mass(), bph::K_MASS, bph::K_MASS},
0169                           {bph::LEP_SIGMA, bph::LEP_SIGMA, bph::K_SIGMA, bph::K_SIGMA});
0170       if (!fitter.success())
0171         continue;
0172       KinVtxFitter fitter_Kpi({leptons_ttracks->at(l1_idx),
0173                                leptons_ttracks->at(l2_idx),
0174                                ditracks_ttracks->at(trk1_idx),
0175                                ditracks_ttracks->at(trk2_idx)},
0176                               {l1_ptr->mass(), l2_ptr->mass(), bph::K_MASS, bph::PI_MASS},
0177                               {bph::LEP_SIGMA, bph::LEP_SIGMA, bph::K_SIGMA, bph::K_SIGMA});
0178       if (!fitter_Kpi.success())
0179         continue;
0180       KinVtxFitter fitter_piK({leptons_ttracks->at(l1_idx),
0181                                leptons_ttracks->at(l2_idx),
0182                                ditracks_ttracks->at(trk1_idx),
0183                                ditracks_ttracks->at(trk2_idx)},
0184                               {l1_ptr->mass(), l2_ptr->mass(), bph::PI_MASS, bph::K_MASS},
0185                               {bph::LEP_SIGMA, bph::LEP_SIGMA, bph::K_SIGMA, bph::K_SIGMA});
0186       if (!fitter_piK.success())
0187         continue;
0188 
0189       // B0 position
0190       cand.setVertex(reco::Candidate::Point(fitter.fitted_vtx().x(), fitter.fitted_vtx().y(), fitter.fitted_vtx().z()));
0191 
0192       // vertex vars
0193       cand.addUserFloat("sv_chi2", fitter.chi2());
0194       cand.addUserFloat("sv_ndof", fitter.dof());
0195       cand.addUserFloat("sv_prob", fitter.prob());
0196 
0197       // refitted kinematic vars
0198       cand.addUserFloat("fitted_ditrack_mass_KK", (fitter.daughter_p4(2) + fitter.daughter_p4(3)).mass());
0199       cand.addUserFloat("fitted_ditrack_mass_Kpi", (fitter_Kpi.daughter_p4(2) + fitter_Kpi.daughter_p4(3)).mass());
0200       cand.addUserFloat("fitted_ditrack_mass_piK", (fitter_piK.daughter_p4(2) + fitter_piK.daughter_p4(3)).mass());
0201       cand.addUserFloat("fitted_massErr_KK", sqrt(fitter.fitted_candidate().kinematicParametersError().matrix()(6, 6)));
0202       cand.addUserFloat("fitted_massErr_Kpi",
0203                         sqrt(fitter_Kpi.fitted_candidate().kinematicParametersError().matrix()(6, 6)));
0204       cand.addUserFloat("fitted_massErr_piK",
0205                         sqrt(fitter_piK.fitted_candidate().kinematicParametersError().matrix()(6, 6)));
0206 
0207       cand.addUserFloat("fitted_mll", (fitter.daughter_p4(0) + fitter.daughter_p4(1)).mass());
0208 
0209       auto fit_p4 = fitter.fitted_p4();
0210       cand.addUserFloat("fitted_pt", fit_p4.pt());
0211       cand.addUserFloat("fitted_eta", fit_p4.eta());
0212       cand.addUserFloat("fitted_phi", fit_p4.phi());
0213 
0214       cand.addUserFloat("fitted_mass_KK", fit_p4.mass());
0215       cand.addUserFloat("fitted_mass_Kpi", fitter_Kpi.fitted_p4().mass());
0216       cand.addUserFloat("fitted_mass_piK", fitter_piK.fitted_p4().mass());
0217 
0218       // other vars
0219       cand.addUserFloat("cos_theta_2D", bph::cos_theta_2D(fitter, *beamspot, cand.p4()));
0220 
0221       cand.addUserFloat("fitted_cos_theta_2D", bph::cos_theta_2D(fitter, *beamspot, fit_p4));
0222 
0223       auto lxy = bph::l_xy(fitter, *beamspot);
0224       cand.addUserFloat("l_xy", lxy.value());
0225       cand.addUserFloat("l_xy_unc", lxy.error());
0226       // track impact parameter from dilepton SV
0227 
0228       TrajectoryStateOnSurface tsos1 =
0229           extrapolator.extrapolate(ditracks_ttracks->at(trk1_idx).impactPointState(), fitter.fitted_vtx());
0230       std::pair<bool, Measurement1D> cur2DIP1 =
0231           bph::signedTransverseImpactParameter(tsos1, fitter.fitted_refvtx(), *beamspot);
0232       cand.addUserFloat("trk1_svip2d", cur2DIP1.second.value());
0233       cand.addUserFloat("trk1_svip2d_err", cur2DIP1.second.error());
0234 
0235       TrajectoryStateOnSurface tsos2 =
0236           extrapolator.extrapolate(ditracks_ttracks->at(trk2_idx).impactPointState(), fitter.fitted_vtx());
0237       std::pair<bool, Measurement1D> cur2DIP2 =
0238           bph::signedTransverseImpactParameter(tsos2, fitter.fitted_refvtx(), *beamspot);
0239       cand.addUserFloat("trk2_svip2d", cur2DIP2.second.value());
0240       cand.addUserFloat("trk2_svip2d_err", cur2DIP2.second.error());
0241 
0242       // post fit selection
0243       if (!post_vtx_selection_(cand))
0244         continue;
0245 
0246       cand.addUserFloat("vtx_x", cand.vx());
0247       cand.addUserFloat("vtx_y", cand.vy());
0248       cand.addUserFloat("vtx_z", cand.vz());
0249 
0250       const auto &covMatrix = fitter.fitted_vtx_uncertainty();
0251       cand.addUserFloat("vtx_cxx", covMatrix.cxx());
0252       cand.addUserFloat("vtx_cyy", covMatrix.cyy());
0253       cand.addUserFloat("vtx_czz", covMatrix.czz());
0254       cand.addUserFloat("vtx_cyx", covMatrix.cyx());
0255       cand.addUserFloat("vtx_czx", covMatrix.czx());
0256       cand.addUserFloat("vtx_czy", covMatrix.czy());
0257 
0258       // refitted daughters (leptons/tracks)
0259       std::vector<std::string> dnames{"l1", "l2", "trk1", "trk2"};
0260 
0261       for (size_t idaughter = 0; idaughter < dnames.size(); idaughter++) {
0262         cand.addUserFloat("fitted_" + dnames[idaughter] + "_pt", fitter.daughter_p4(idaughter).pt());
0263         cand.addUserFloat("fitted_" + dnames[idaughter] + "_eta", fitter.daughter_p4(idaughter).eta());
0264         cand.addUserFloat("fitted_" + dnames[idaughter] + "_phi", fitter.daughter_p4(idaughter).phi());
0265       }
0266 
0267       // compute isolation
0268       std::vector<float> isos = bph::TrackerIsolation(pu_tracks, cand, dnames);
0269       for (size_t idaughter = 0; idaughter < dnames.size(); idaughter++) {
0270         cand.addUserFloat(dnames[idaughter] + "_iso04", isos[idaughter]);
0271       }
0272 
0273       float constraint_sv_prob = -9;
0274       float constraint_pt = -9;
0275       float constraint_eta = -9;
0276       float constraint_phi = -9;
0277       float constraint_mass_KK = -9;
0278       float constraint_mass_Kpi = -9;
0279       float constraint_mass_piK = -9;
0280       float constraint_massErr_KK = -9;
0281       float constraint_massErr_Kpi = -9;
0282       float constraint_massErr_piK = -9;
0283       float constraint_mll = -9;
0284 
0285       const double dilepton_mass = ll_ptr->userFloat("fitted_mass");
0286       const double jpsi_bin[2] = {2.8, 3.35};
0287       const double psi2s_bin[2] = {3.45, 3.85};
0288 
0289       if (dilepton_constraint_ && ((dilepton_mass > jpsi_bin[0] && dilepton_mass < jpsi_bin[1]) ||
0290                                    (dilepton_mass > psi2s_bin[0] && dilepton_mass < psi2s_bin[1]))) {
0291         ParticleMass JPsi_mass = 3.0969;   // Jpsi mass 3.096900±0.000006
0292         ParticleMass Psi2S_mass = 3.6861;  // Psi2S mass 3.6861093±0.0000034
0293         ParticleMass mass_constraint = (dilepton_mass < jpsi_bin[1]) ? JPsi_mass : Psi2S_mass;
0294 
0295         // Mass constraint is applied to the first two particles in the
0296         // "particles" vector Make sure that the first two particles are the
0297         // ones you want to constrain
0298 
0299         KinVtxFitter constraint_fitter_KK({leptons_ttracks->at(l1_idx),
0300                                            leptons_ttracks->at(l2_idx),
0301                                            ditracks_ttracks->at(trk1_idx),
0302                                            ditracks_ttracks->at(trk2_idx)},
0303                                           {l1_ptr->mass(), l2_ptr->mass(), bph::K_MASS, bph::K_MASS},
0304                                           {bph::LEP_SIGMA, bph::LEP_SIGMA, bph::K_SIGMA, bph::K_SIGMA},
0305                                           mass_constraint);
0306         if (!constraint_fitter_KK.success())
0307           continue;
0308         KinVtxFitter constraint_fitter_Kpi({leptons_ttracks->at(l1_idx),
0309                                             leptons_ttracks->at(l2_idx),
0310                                             ditracks_ttracks->at(trk1_idx),
0311                                             ditracks_ttracks->at(trk2_idx)},
0312                                            {l1_ptr->mass(), l2_ptr->mass(), bph::K_MASS, bph::PI_MASS},
0313                                            {bph::LEP_SIGMA, bph::LEP_SIGMA, bph::K_SIGMA, bph::K_SIGMA},
0314                                            mass_constraint);
0315         if (!constraint_fitter_Kpi.success())
0316           continue;
0317         KinVtxFitter constraint_fitter_piK({leptons_ttracks->at(l1_idx),
0318                                             leptons_ttracks->at(l2_idx),
0319                                             ditracks_ttracks->at(trk1_idx),
0320                                             ditracks_ttracks->at(trk2_idx)},
0321                                            {l1_ptr->mass(), l2_ptr->mass(), bph::PI_MASS, bph::K_MASS},
0322                                            {bph::LEP_SIGMA, bph::LEP_SIGMA, bph::K_SIGMA, bph::K_SIGMA},
0323                                            mass_constraint);
0324         if (!constraint_fitter_piK.success())
0325           continue;
0326 
0327         if (constraint_fitter_KK.success()) {
0328           auto constraint_p4 = constraint_fitter_KK.fitted_p4();
0329           constraint_sv_prob = constraint_fitter_KK.prob();
0330           constraint_pt = constraint_p4.pt();
0331           constraint_eta = constraint_p4.eta();
0332           constraint_phi = constraint_p4.phi();
0333           constraint_mass_KK = constraint_fitter_KK.fitted_candidate().mass();
0334           constraint_massErr_KK =
0335               sqrt(constraint_fitter_KK.fitted_candidate().kinematicParametersError().matrix()(6, 6));
0336           constraint_mass_Kpi = constraint_fitter_Kpi.fitted_candidate().mass();
0337           constraint_massErr_Kpi =
0338               sqrt(constraint_fitter_Kpi.fitted_candidate().kinematicParametersError().matrix()(6, 6));
0339           constraint_mass_piK = constraint_fitter_piK.fitted_candidate().mass();
0340           constraint_massErr_piK =
0341               sqrt(constraint_fitter_piK.fitted_candidate().kinematicParametersError().matrix()(6, 6));
0342           constraint_mll = (constraint_fitter_KK.daughter_p4(0) + constraint_fitter_KK.daughter_p4(1)).mass();
0343         }
0344       }
0345       cand.addUserFloat("constraint_sv_prob", constraint_sv_prob);
0346       cand.addUserFloat("constraint_pt", constraint_pt);
0347       cand.addUserFloat("constraint_eta", constraint_eta);
0348       cand.addUserFloat("constraint_phi", constraint_phi);
0349       cand.addUserFloat("constraint_mass_KK", constraint_mass_KK);
0350       cand.addUserFloat("constraint_mass_Kpi", constraint_mass_Kpi);
0351       cand.addUserFloat("constraint_mass_piK", constraint_mass_piK);
0352       cand.addUserFloat("constraint_massErr_KK", constraint_massErr_KK);
0353       cand.addUserFloat("constraint_massErr_Kpi", constraint_massErr_Kpi);
0354       cand.addUserFloat("constraint_massErr_piK", constraint_massErr_piK);
0355       cand.addUserFloat("constraint_mll", constraint_mll);
0356 
0357       ret_val->emplace_back(cand);
0358       ditrack_out->emplace_back(*ditracks_ptr);
0359 
0360     }  // for(size_t ll_idx = 0; ll_idx < dileptons->size(); ++ll_idx) {
0361 
0362   }  // for(size_t k_idx = 0; k_idx < ditracks->size(); ++k_idx)
0363 
0364   evt.put(std::move(ret_val), "SelectedBToTrkTrkMuMu");
0365   evt.put(std::move(ditrack_out), "SelectedTrkTrk");
0366 }
0367 
0368 #include "FWCore/Framework/interface/MakerMacros.h"
0369 DEFINE_FWK_MODULE(BToTrkTrkLLBuilder);
[ 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!