BPHNano/plugins/BToV0LLBuilder.cc

0001 /////////////////////////////// BToV0LLBuilder ///////////////////////////////
0002 /// original authors: G Karathanasis (CERN),  G Melachroinos (NKUA)
0003 /// takes rebuilt V0 cands and a dilepton collection and produces B mothers
0004
0005 #include <algorithm>
0006 #include <limits>
0007 #include <map>
0008 #include <memory>
0009 #include <string>
0010 #include <vector>
0011
0012 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0013 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0014 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0015 #include "DataFormats/Candidate/interface/VertexCompositeCandidate.h"
0016 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0017 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidateFwd.h"
0018 #include "DataFormats/Math/interface/deltaR.h"
0019 #include "DataFormats/PatCandidates/interface/CompositeCandidate.h"
0020 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/global/EDProducer.h"
0023 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0026 #include "FWCore/Utilities/interface/InputTag.h"
0027 #include "KinVtxFitter.h"
0028 #include "MagneticField/Engine/interface/MagneticField.h"
0029 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0030 #include "RecoVertex/KinematicFit/interface/KinematicConstrainedVertexFitter.h"
0031 #include "RecoVertex/KinematicFit/interface/TwoTrackMassKinematicConstraint.h"
0032 #include "RecoVertex/KinematicFitPrimitives/interface/KinematicParticleFactoryFromTransientTrack.h"
0033 #include "RecoVertex/KinematicFitPrimitives/interface/MultiTrackKinematicConstraint.h"
0034 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0035 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0036 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0037 #include "helper.h"
0038
0039 class BToV0LLBuilder : public edm::global::EDProducer<> {
0040   // perhaps we need better structure here (begin run etc)
0041 public:
0042   typedef std::vector<reco::TransientTrack> TransientTrackCollection;
0043
0044   explicit BToV0LLBuilder(const edm::ParameterSet &cfg)
0045       : bFieldToken_{esConsumes<MagneticField, IdealMagneticFieldRecord>()},
0046         pre_vtx_selection_{cfg.getParameter<std::string>("preVtxSelection")},
0047         post_vtx_selection_{cfg.getParameter<std::string>("postVtxSelection")},
0048         dileptons_{consumes<pat::CompositeCandidateCollection>(cfg.getParameter<edm::InputTag>("dileptons"))},
0049         leptons_ttracks_{consumes<TransientTrackCollection>(cfg.getParameter<edm::InputTag>("leptonTransientTracks"))},
0050         v0s_{consumes<pat::CompositeCandidateCollection>(cfg.getParameter<edm::InputTag>("v0s"))},
0051         v0_ttracks_{consumes<TransientTrackCollection>(cfg.getParameter<edm::InputTag>("v0TransientTracks"))},
0052         pu_tracks_(consumes<pat::CompositeCandidateCollection>(cfg.getParameter<edm::InputTag>("PUtracks"))),
0053         beamspot_{consumes<reco::BeamSpot>(cfg.getParameter<edm::InputTag>("beamSpot"))},
0054         dilepton_constraint_{cfg.getParameter<bool>("dileptonMassContraint")} {
0055     produces<pat::CompositeCandidateCollection>();
0056   }
0057
0058   ~BToV0LLBuilder() override {}
0059
0060   void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override;
0061
0062 private:
0063   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> bFieldToken_;
0064
0065   // selection
0066   const StringCutObjectSelector<pat::CompositeCandidate> pre_vtx_selection_;
0067   const StringCutObjectSelector<pat::CompositeCandidate> post_vtx_selection_;
0068
0069   // input
0070   const edm::EDGetTokenT<pat::CompositeCandidateCollection> dileptons_;
0071   const edm::EDGetTokenT<TransientTrackCollection> leptons_ttracks_;
0072   const edm::EDGetTokenT<pat::CompositeCandidateCollection> v0s_;
0073   const edm::EDGetTokenT<TransientTrackCollection> v0_ttracks_;
0074   const edm::EDGetTokenT<pat::CompositeCandidateCollection> pu_tracks_;
0075   const edm::EDGetTokenT<reco::BeamSpot> beamspot_;
0076   const bool dilepton_constraint_;
0077 };
0078
0079 void BToV0LLBuilder::produce(edm::StreamID, edm::Event &evt, edm::EventSetup const &iSetup) const {
0080   // input
0081   edm::Handle<pat::CompositeCandidateCollection> dileptons;
0082   evt.getByToken(dileptons_, dileptons);
0083   edm::Handle<TransientTrackCollection> leptons_ttracks;
0084   evt.getByToken(leptons_ttracks_, leptons_ttracks);
0085
0086   edm::Handle<pat::CompositeCandidateCollection> v0s;
0087   evt.getByToken(v0s_, v0s);
0088   edm::Handle<TransientTrackCollection> v0_ttracks;
0089   evt.getByToken(v0_ttracks_, v0_ttracks);
0090
0091   edm::Handle<pat::CompositeCandidateCollection> pu_tracks;
0092   evt.getByToken(pu_tracks_, pu_tracks);
0093
0094   edm::Handle<reco::BeamSpot> beamspot;
0095   evt.getByToken(beamspot_, beamspot);
0096
0097   edm::ESHandle<MagneticField> fieldHandle;
0098   const auto &bField = iSetup.getData(bFieldToken_);
0099   AnalyticalImpactPointExtrapolator extrapolator(&bField);
0100
0101   // output
0102   std::unique_ptr<pat::CompositeCandidateCollection> ret_val(new pat::CompositeCandidateCollection());
0103
0104   // access V0
0105   for (size_t v0_idx = 0; v0_idx < v0s->size(); ++v0_idx) {
0106     edm::Ptr<pat::CompositeCandidate> v0_ptr(v0s, v0_idx);
0107
0108     // access ll
0109     for (size_t ll_idx = 0; ll_idx < dileptons->size(); ++ll_idx) {
0110       edm::Ptr<pat::CompositeCandidate> ll_ptr(dileptons, ll_idx);
0111       edm::Ptr<reco::Candidate> l1_ptr = ll_ptr->userCand("l1");
0112       edm::Ptr<reco::Candidate> l2_ptr = ll_ptr->userCand("l2");
0113       int l1_idx = ll_ptr->userInt("l1_idx");
0114       int l2_idx = ll_ptr->userInt("l2_idx");
0115
0116       pat::CompositeCandidate cand;
0117       cand.setP4(ll_ptr->p4() + v0_ptr->p4());
0118       cand.setCharge(ll_ptr->charge() + v0_ptr->charge());
0119
0120       cand.addUserInt("l1_idx", l1_idx);
0121       cand.addUserInt("l2_idx", l2_idx);
0122       cand.addUserInt("ll_idx", ll_idx);
0123       cand.addUserInt("v0_idx", v0_idx);
0124
0125       auto dr_info = bph::min_max_dr({l1_ptr, l2_ptr, v0_ptr});
0126       cand.addUserFloat("min_dr", dr_info.first);
0127       cand.addUserFloat("max_dr", dr_info.second);
0128
0129       // built B
0130       if (!pre_vtx_selection_(cand))
0131         continue;
0132
0133       KinVtxFitter fitter({leptons_ttracks->at(l1_idx), leptons_ttracks->at(l2_idx), v0_ttracks->at(v0_idx)},
0134                           {l1_ptr->mass(), l2_ptr->mass(), v0_ptr->mass()},
0135                           {bph::LEP_SIGMA, bph::LEP_SIGMA, v0_ptr->userFloat("massErr")});
0136
0137       if (!fitter.success())
0138         continue;
0139
0140       cand.setVertex(reco::Candidate::Point(fitter.fitted_vtx().x(), fitter.fitted_vtx().y(), fitter.fitted_vtx().z()));
0141
0142       cand.addUserFloat("sv_chi2", fitter.chi2());
0143       cand.addUserFloat("sv_ndof", fitter.dof());
0144       cand.addUserFloat("sv_prob", fitter.prob());
0145       cand.addUserFloat("fitted_mll", (fitter.daughter_p4(0) + fitter.daughter_p4(1)).mass());
0146       cand.addUserFloat("fitted_v0_mass", fitter.daughter_p4(2).mass());
0147
0148       auto fit_p4 = fitter.fitted_p4();
0149       cand.addUserFloat("fitted_pt", fit_p4.pt());
0150       cand.addUserFloat("fitted_eta", fit_p4.eta());
0151       cand.addUserFloat("fitted_phi", fit_p4.phi());
0152       cand.addUserFloat("fitted_mass", fitter.fitted_candidate().mass());
0153       cand.addUserFloat("fitted_massErr", sqrt(fitter.fitted_candidate().kinematicParametersError().matrix()(6, 6)));
0154       cand.addUserFloat("cos_theta_2D", bph::cos_theta_2D(fitter, *beamspot, cand.p4()));
0155       cand.addUserFloat("fitted_cos_theta_2D", bph::cos_theta_2D(fitter, *beamspot, fit_p4));
0156
0157       auto lxy = bph::l_xy(fitter, *beamspot);
0158       cand.addUserFloat("l_xy", lxy.value());
0159       cand.addUserFloat("l_xy_unc", lxy.error());
0160
0161       TrajectoryStateOnSurface tsos =
0162           extrapolator.extrapolate(v0_ttracks->at(v0_idx).impactPointState(), fitter.fitted_vtx());
0163       std::pair<bool, Measurement1D> cur2DIP =
0164           bph::signedTransverseImpactParameter(tsos, fitter.fitted_refvtx(), *beamspot);
0165       cand.addUserFloat("v0_svip2d", cur2DIP.second.value());
0166       cand.addUserFloat("v0_svip2d_err", cur2DIP.second.error());
0167
0168       if (!post_vtx_selection_(cand))
0169         continue;
0170
0171       cand.addUserFloat("vtx_x", cand.vx());
0172       cand.addUserFloat("vtx_y", cand.vy());
0173       cand.addUserFloat("vtx_z", cand.vz());
0174
0175       const auto &covMatrix = fitter.fitted_vtx_uncertainty();
0176       cand.addUserFloat("vtx_cxx", covMatrix.cxx());
0177       cand.addUserFloat("vtx_cyy", covMatrix.cyy());
0178       cand.addUserFloat("vtx_czz", covMatrix.czz());
0179       cand.addUserFloat("vtx_cyx", covMatrix.cyx());
0180       cand.addUserFloat("vtx_czx", covMatrix.czx());
0181       cand.addUserFloat("vtx_czy", covMatrix.czy());
0182
0183       // refitted daughters (leptons/tracks)
0184       std::vector<std::string> dnames{"l1", "l2", "v0"};
0185       for (size_t idaughter = 0; idaughter < dnames.size(); idaughter++) {
0186         cand.addUserFloat("fitted_" + dnames[idaughter] + "_pt", fitter.daughter_p4(idaughter).pt());
0187         cand.addUserFloat("fitted_" + dnames[idaughter] + "_eta", fitter.daughter_p4(idaughter).eta());
0188         cand.addUserFloat("fitted_" + dnames[idaughter] + "_phi", fitter.daughter_p4(idaughter).phi());
0189       }
0190
0191       // compute isolation
0192       std::vector<float> isos = bph::TrackerIsolation(pu_tracks, cand, dnames);
0193       for (size_t idaughter = 0; idaughter < dnames.size(); idaughter++) {
0194         cand.addUserFloat(dnames[idaughter] + "_iso04", isos[idaughter]);
0195       }
0196
0197       float cstr_pt = -99;
0198       float cstr_eta = -99;
0199       float cstr_phi = -99;
0200       float cstr_sv_prob = -99;
0201       float cstr_mass = -99;
0202       float cstr_massErr = -99;
0203       float cstr_vtx_x = -99;
0204       float cstr_vtx_y = -99;
0205       float cstr_vtx_z = -99;
0206       float cstr_vtx_cxx = -99;
0207       float cstr_vtx_cyy = -99;
0208       float cstr_vtx_czz = -99;
0209       float cstr_vtx_cyx = -99;
0210       float cstr_vtx_czx = -99;
0211       float cstr_vtx_czy = -99;
0212       float cstr_fitted_l1_pt = -99;
0213       float cstr_fitted_l1_eta = -99;
0214       float cstr_fitted_l1_phi = -99;
0215       float cstr_fitted_l2_pt = -99;
0216       float cstr_fitted_l2_eta = -99;
0217       float cstr_fitted_l2_phi = -99;
0218       float cstr_fitted_v0_pt = -99;
0219       float cstr_fitted_v0_eta = -99;
0220       float cstr_fitted_v0_phi = -99;
0221       float cstr_v0_mass = -99;
0222       float cstr_cos_theta_2D = -99;
0223
0224       const double dilepton_mass = ll_ptr->userFloat("fitted_mass");
0225       const double jpsi_bin[2] = {2.8, 3.35};
0226       const double psi2s_bin[2] = {3.45, 3.85};
0227       if (dilepton_constraint_ && ((dilepton_mass > jpsi_bin[0] && dilepton_mass < jpsi_bin[1]) ||
0228                                    (dilepton_mass > psi2s_bin[0] && dilepton_mass < psi2s_bin[1]))) {
0229         ParticleMass JPsi_mass = 3.0969;   // Jpsi mass 3.096900±0.000006
0230         ParticleMass Psi2S_mass = 3.6861;  // Psi2S mass 3.6861093±0.0000034
0231         ParticleMass mass_constraint = (dilepton_mass < jpsi_bin[1]) ? JPsi_mass : Psi2S_mass;
0232
0233         // Mass constraint is applied to the first two particles in the
0234         // "particles" vector Make sure that the first two particles are the
0235         // ones you want to constrain
0236         KinVtxFitter constraint_fitter(
0237             {leptons_ttracks->at(l1_idx), leptons_ttracks->at(l2_idx), v0_ttracks->at(v0_idx)},
0238             {l1_ptr->mass(), l2_ptr->mass(), v0_ptr->mass()},
0239             {bph::LEP_SIGMA, bph::LEP_SIGMA, v0_ptr->userFloat("massErr")},
0240             mass_constraint);
0241         if (constraint_fitter.success()) {
0242           auto constraint_p4 = constraint_fitter.fitted_p4();
0243           cstr_vtx_x = constraint_fitter.fitted_vtx().x();
0244           cstr_vtx_y = constraint_fitter.fitted_vtx().y();
0245           cstr_vtx_z = constraint_fitter.fitted_vtx().z();
0246
0247           const auto &constrained_covMatrix = constraint_fitter.fitted_vtx_uncertainty();
0248           cstr_vtx_cxx = constrained_covMatrix.cxx();
0249           cstr_vtx_cyy = constrained_covMatrix.cyy();
0250           cstr_vtx_czz = constrained_covMatrix.czz();
0251           cstr_vtx_cyx = constrained_covMatrix.cyx();
0252           cstr_vtx_czx = constrained_covMatrix.czx();
0253           cstr_vtx_czy = constrained_covMatrix.czy();
0254
0255           // refitted daughters (leptons/tracks)i
0256           cstr_fitted_l1_pt = constraint_fitter.daughter_p4(0).pt();
0257           cstr_fitted_l1_eta = constraint_fitter.daughter_p4(0).eta();
0258           cstr_fitted_l1_phi = constraint_fitter.daughter_p4(0).phi();
0259           cstr_fitted_l2_pt = constraint_fitter.daughter_p4(1).pt();
0260           cstr_fitted_l2_eta = constraint_fitter.daughter_p4(1).eta();
0261           cstr_fitted_l2_phi = constraint_fitter.daughter_p4(1).phi();
0262           cstr_fitted_v0_pt = constraint_fitter.daughter_p4(2).pt();
0263           cstr_fitted_v0_eta = constraint_fitter.daughter_p4(2).eta();
0264           cstr_fitted_v0_phi = constraint_fitter.daughter_p4(2).phi();
0265
0266           cstr_v0_mass = constraint_fitter.daughter_p4(2).mass();
0267
0268           cstr_sv_prob = constraint_fitter.prob();
0269           cstr_pt = constraint_p4.pt();
0270           cstr_eta = constraint_p4.eta();
0271           cstr_phi = constraint_p4.phi();
0272           cstr_mass = constraint_fitter.fitted_candidate().mass();
0273           cstr_massErr = sqrt(constraint_fitter.fitted_candidate().kinematicParametersError().matrix()(6, 6));
0274           cstr_cos_theta_2D = bph::cos_theta_2D(constraint_fitter, *beamspot, constraint_p4);
0275         }
0276       }
0277
0278       cand.addUserFloat("cstr_vtx_x", cstr_vtx_x);
0279       cand.addUserFloat("cstr_vtx_y", cstr_vtx_y);
0280       cand.addUserFloat("cstr_vtx_z", cstr_vtx_z);
0281
0282       cand.addUserFloat("cstr_vtx_cxx", cstr_vtx_cxx);
0283       cand.addUserFloat("cstr_vtx_cyy", cstr_vtx_cyy);
0284       cand.addUserFloat("cstr_vtx_czz", cstr_vtx_czz);
0285       cand.addUserFloat("cstr_vtx_cyx", cstr_vtx_cyx);
0286       cand.addUserFloat("cstr_vtx_czx", cstr_vtx_czx);
0287       cand.addUserFloat("cstr_vtx_czy", cstr_vtx_czy);
0288
0289       // refitted daughters (leptons/tracks)
0290       cand.addUserFloat("cstr_fitted_l1_pt", cstr_fitted_l1_pt);
0291       cand.addUserFloat("cstr_fitted_l1_eta", cstr_fitted_l1_eta);
0292       cand.addUserFloat("cstr_fitted_l1_phi", cstr_fitted_l1_phi);
0293       cand.addUserFloat("cstr_fitted_l2_pt", cstr_fitted_l2_pt);
0294       cand.addUserFloat("cstr_fitted_l2_eta", cstr_fitted_l2_eta);
0295       cand.addUserFloat("cstr_fitted_l2_phi", cstr_fitted_l2_phi);
0296       cand.addUserFloat("cstr_fitted_v0_pt", cstr_fitted_v0_pt);
0297       cand.addUserFloat("cstr_fitted_v0_eta", cstr_fitted_v0_eta);
0298       cand.addUserFloat("cstr_fitted_v0_phi", cstr_fitted_v0_phi);
0299       cand.addUserFloat("cstr_v0_mass", cstr_v0_mass);
0300
0301       cand.addUserFloat("cstr_sv_prob", cstr_sv_prob);
0302       cand.addUserFloat("cstr_pt", cstr_pt);
0303       cand.addUserFloat("cstr_eta", cstr_eta);
0304       cand.addUserFloat("cstr_phi", cstr_phi);
0305       cand.addUserFloat("cstr_mass", cstr_mass);
0306       cand.addUserFloat("cstr_massErr", cstr_massErr);
0307       cand.addUserFloat("cstr_cos_theta_2D", cstr_cos_theta_2D);
0308
0309       ret_val->push_back(cand);
0310     }
0311   }
0312   evt.put(std::move(ret_val));
0313 }
0314
0315 #include "FWCore/Framework/interface/MakerMacros.h"
0316 DEFINE_FWK_MODULE(BToV0LLBuilder);