Skip to content

Implementation of the fully hadronic modes using single-muon triggered data #4

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
amlyon wants to merge 11 commits into
base: main
Choose a base branch
from
Open

Implementation of the fully hadronic modes using single-muon triggered data #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
36d2a15
BTo4TrkBuilder added
May 12, 2025
6faf1fa
BsToPhiPhiTo4K_cff added
May 12, 2025
a8b64d8
DiTrackBuilder modified - updated loop strategy + sorting of candidates
May 12, 2025
83125e1
TrackMerger modified to be able to run on both dimuon and single-muon…
May 12, 2025
93e784d
tracks_cff modified to run on both dimuon and single-muon channels
May 12, 2025
826df64
track class name modified according to dimuon or single-muon channel
May 12, 2025
781253e
muons_cff modified to be able to run on both dimuon and single-muon c…
May 12, 2025
69519a9
nanoBPH_cff modified
May 12, 2025
ac6e193
handling of the muon producer modified to run consistently on mc
May 12, 2025
50865d4
fix
May 14, 2025
f97dfc1
run_bphNano_cfg modified
May 14, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
245 changes: 245 additions & 0 deletions BPHNano/plugins/BTo4TrkBuilder.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,245 @@
#include "FWCore/Framework/interface/global/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "FWCore/Utilities/interface/InputTag.h"

#include <vector>
#include <memory>
#include <map>
#include <string>
#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
#include "CommonTools/Utils/interface/StringCutObjectSelector.h"
#include "DataFormats/PatCandidates/interface/CompositeCandidate.h"
#include "DataFormats/Math/interface/deltaR.h"
#include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include "helper.h"
#include <limits>
#include <algorithm>
#include "KinVtxFitter.h"
#include "RecoVertex/VertexTools/interface/VertexDistance3D.h"

/*
* This class was originally designed to reconstruct the Bs -> phi(KK)phi/KK) decay.
* It can be generalised to other fully-hadronic B decays.
*/

class BTo4TrkBuilder : public edm::global::EDProducer<> {

public:
typedef std::vector<reco::TransientTrack> TransientTrackCollection;

explicit BTo4TrkBuilder(const edm::ParameterSet &cfg):
phi_candidates_{consumes<pat::CompositeCandidateCollection>( cfg.getParameter<edm::InputTag>("phis") )},
kaon_ttracks_{consumes<TransientTrackCollection>( cfg.getParameter<edm::InputTag>("phisTransientTracks") )},
pre_vtx_selection_ {cfg.getParameter<std::string>("pre_vtx_selection")},
post_vtx_selection_ {cfg.getParameter<std::string>("post_vtx_selection")},
beamspot_{consumes<reco::BeamSpot>(cfg.getParameter<edm::InputTag>("beamSpot"))}{
produces<pat::CompositeCandidateCollection>();
}

~BTo4TrkBuilder() override {}

void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;

static void fillDescriptions(edm::ConfigurationDescriptions &descriptions) {}

private:
const edm::EDGetTokenT<pat::CompositeCandidateCollection> phi_candidates_;
const edm::EDGetTokenT<TransientTrackCollection> kaon_ttracks_;
const StringCutObjectSelector<pat::CompositeCandidate> pre_vtx_selection_;
const StringCutObjectSelector<pat::CompositeCandidate> post_vtx_selection_;
const edm::EDGetTokenT<reco::BeamSpot> beamspot_;
};

void BTo4TrkBuilder::produce(edm::StreamID, edm::Event &evt, edm::EventSetup const &) const {

edm::Handle<pat::CompositeCandidateCollection> phi_candidates;
evt.getByToken(phi_candidates_, phi_candidates);

edm::Handle<TransientTrackCollection> kaon_ttracks;
evt.getByToken(kaon_ttracks_, kaon_ttracks);

edm::Handle<reco::BeamSpot> beamspot;
evt.getByToken(beamspot_, beamspot);

// output
std::unique_ptr<pat::CompositeCandidateCollection> ret_value(new pat::CompositeCandidateCollection());

// loop on first phi candidate (phi1)
for(size_t phi1_idx = 0; phi1_idx < phi_candidates->size(); ++phi1_idx) {
edm::Ptr<pat::CompositeCandidate> phi1_ptr(phi_candidates, phi1_idx);

// retrieve kaons
edm::Ptr<reco::Candidate> k1_ptr = phi1_ptr->userCand("trk1");
edm::Ptr<reco::Candidate> k2_ptr = phi1_ptr->userCand("trk2");
int k1_idx = phi1_ptr->userInt("trk1_idx");
int k2_idx = phi1_ptr->userInt("trk2_idx");

// loop on second phi candidate (phi2)
for(size_t phi2_idx = phi1_idx + 1; phi2_idx < phi_candidates->size(); ++phi2_idx) {
edm::Ptr<pat::CompositeCandidate> phi2_ptr(phi_candidates, phi2_idx);

// retrieve kaons
edm::Ptr<reco::Candidate> k3_ptr = phi2_ptr->userCand("trk1");
edm::Ptr<reco::Candidate> k4_ptr = phi2_ptr->userCand("trk2");
int k3_idx = phi2_ptr->userInt("trk1_idx");
int k4_idx = phi2_ptr->userInt("trk2_idx");

if(phi2_ptr->pt() > phi1_ptr->pt()){
std::cout << "not ordered in pt!" << std::endl;
}

if(k3_idx == k1_idx || k3_idx == k2_idx) {
continue;
}
if(k4_idx == k1_idx || k4_idx == k2_idx){
continue;
}

// build Bs candidate
pat::CompositeCandidate Bs_cand;
Bs_cand.setP4(phi1_ptr->p4() + phi2_ptr->p4());
Bs_cand.setCharge(phi1_ptr->charge() + phi2_ptr->charge());

// apply pre-fit selection on Bs candidate
if( !pre_vtx_selection_(Bs_cand) ) continue;

// save indices
Bs_cand.addUserInt("k1_idx", k1_idx);
Bs_cand.addUserInt("k2_idx", k2_idx);
Bs_cand.addUserInt("k3_idx", k3_idx);
Bs_cand.addUserInt("k4_idx", k4_idx);
Bs_cand.addUserInt("phi1_idx", phi1_idx);
Bs_cand.addUserInt("phi2_idx", phi2_idx);

// save candidates
Bs_cand.addUserCand("k1", k1_ptr);
Bs_cand.addUserCand("k2", k2_ptr);
Bs_cand.addUserCand("k3", k3_ptr);
Bs_cand.addUserCand("k4", k4_ptr);
Bs_cand.addUserCand("phi1", phi1_ptr);
Bs_cand.addUserCand("phi2", phi2_ptr);

// fit the four kaon tracks to a common vertex
KinVtxFitter fitter_Bs(
{kaon_ttracks->at(k1_idx), kaon_ttracks->at(k2_idx), kaon_ttracks->at(k3_idx), kaon_ttracks->at(k4_idx)},
{K_MASS, K_MASS, K_MASS, K_MASS}, // force kaon mass
{K_SIGMA, K_SIGMA, K_SIGMA, K_SIGMA}
);

if(!fitter_Bs.success()) continue;

Bs_cand.setVertex(
reco::Candidate::Point(
fitter_Bs.fitted_vtx().x(),
fitter_Bs.fitted_vtx().y(),
fitter_Bs.fitted_vtx().z()
)
);

Bs_cand.addUserFloat("Bs_vx", fitter_Bs.fitted_vtx().x());
Bs_cand.addUserFloat("Bs_vy", fitter_Bs.fitted_vtx().y());
Bs_cand.addUserFloat("Bs_vz", fitter_Bs.fitted_vtx().z());

const auto& covMatrix = fitter_Bs.fitted_vtx_uncertainty();
Bs_cand.addUserFloat("Bs_vtx_cxx", covMatrix.cxx());
Bs_cand.addUserFloat("Bs_vtx_cyy", covMatrix.cyy());
Bs_cand.addUserFloat("Bs_vtx_czz", covMatrix.czz());
Bs_cand.addUserFloat("Bs_vtx_cyx", covMatrix.cyx());
Bs_cand.addUserFloat("Bs_vtx_czx", covMatrix.czx());
Bs_cand.addUserFloat("Bs_vtx_czy", covMatrix.czy());

Bs_cand.addUserFloat("Bs_sv_chi2", fitter_Bs.chi2());
Bs_cand.addUserFloat("Bs_sv_ndof", fitter_Bs.dof());
Bs_cand.addUserFloat("Bs_sv_prob", fitter_Bs.prob());

auto fit_p4 = fitter_Bs.fitted_p4();
Bs_cand.addUserFloat("Bs_fitted_pt" , fit_p4.pt());
Bs_cand.addUserFloat("Bs_fitted_eta" , fit_p4.eta());
Bs_cand.addUserFloat("Bs_fitted_phi" , fit_p4.phi());
Bs_cand.addUserFloat("Bs_fitted_mass", fit_p4.mass());
Bs_cand.addUserFloat("Bs_fitted_mass_corr", fit_p4.mass() - (fitter_Bs.daughter_p4(0) + fitter_Bs.daughter_p4(1)).mass() + 1.0195 - (fitter_Bs.daughter_p4(2) + fitter_Bs.daughter_p4(3)).mass() + 1.0195);
Bs_cand.addUserFloat("Bs_fitted_massErr", sqrt(fitter_Bs.fitted_candidate().kinematicParametersError().matrix()(6,6)));
Bs_cand.addUserFloat("Bs_charge", Bs_cand.charge());
Bs_cand.addUserFloat("Bs_cos_theta_2D", cos_theta_2D(fitter_Bs, *beamspot, fit_p4));
auto Bs_lxy = l_xy(fitter_Bs, *beamspot);
Bs_cand.addUserFloat("Bs_lxy", Bs_lxy.value());
Bs_cand.addUserFloat("Bs_lxy_sig", Bs_lxy.value() / Bs_lxy.error());

auto k1_p4 = fitter_Bs.daughter_p4(0);
auto k2_p4 = fitter_Bs.daughter_p4(1);
auto k3_p4 = fitter_Bs.daughter_p4(2);
auto k4_p4 = fitter_Bs.daughter_p4(3);
auto phi1_p4 = k1_p4 + k2_p4;
auto phi2_p4 = k3_p4 + k4_p4;

// phi candidates
std::vector<std::string> phi_names{"phi1", "phi2"};
for (size_t i = 0; i < phi_names.size(); i++){
int idx1 = -99;
int idx2 = -99;
if(phi_names[i] == "phi1"){
idx1 = 0;
idx2 = 1;
}
else if(phi_names[i] == "phi2"){
idx1 = 2;
idx2 = 3;
}
Bs_cand.addUserFloat(phi_names[i] + "_fitted_mass", (fitter_Bs.daughter_p4(idx1) + fitter_Bs.daughter_p4(idx2)).mass());
Bs_cand.addUserFloat(phi_names[i] + "_fitted_pt", (fitter_Bs.daughter_p4(idx1) + fitter_Bs.daughter_p4(idx2)).pt());
Bs_cand.addUserFloat(phi_names[i] + "_fitted_eta", (fitter_Bs.daughter_p4(idx1) + fitter_Bs.daughter_p4(idx2)).eta());
Bs_cand.addUserFloat(phi_names[i] + "_fitted_phi", (fitter_Bs.daughter_p4(idx1) + fitter_Bs.daughter_p4(idx2)).phi());
}

// compute deltaR
Bs_cand.addUserFloat("deltaR_phi1phi2", reco::deltaR(phi1_p4, phi2_p4));
Bs_cand.addUserFloat("deltaR_k1k3", reco::deltaR(k1_p4, k3_p4));
Bs_cand.addUserFloat("deltaR_k1k4", reco::deltaR(k1_p4, k4_p4));
Bs_cand.addUserFloat("deltaR_k2k3", reco::deltaR(k2_p4, k3_p4));
Bs_cand.addUserFloat("deltaR_k2k4", reco::deltaR(k2_p4, k4_p4));

Bs_cand.addUserFloat("deltaR_k1k2", reco::deltaR(k1_p4, k2_p4));
Bs_cand.addUserFloat("deltaR_k3k4", reco::deltaR(k3_p4, k4_p4));

auto dr_info = min_max_dr({k1_ptr, k2_ptr, k3_ptr, k4_ptr});
Bs_cand.addUserFloat("deltaR_min", dr_info.first);
Bs_cand.addUserFloat("deltaR_max", dr_info.second);

// compute invariant masses
Bs_cand.addUserFloat("k1k3_mass", (k1_p4 + k3_p4).mass());
Bs_cand.addUserFloat("k1k4_mass", (k1_p4 + k4_p4).mass());
Bs_cand.addUserFloat("k2k3_mass", (k2_p4 + k3_p4).mass());
Bs_cand.addUserFloat("k2k4_mass", (k2_p4 + k4_p4).mass());

Bs_cand.addUserFloat("k1k3_pt", (k1_p4 + k3_p4).pt());
Bs_cand.addUserFloat("k1k4_pt", (k1_p4 + k4_p4).pt());
Bs_cand.addUserFloat("k2k3_pt", (k2_p4 + k3_p4).pt());
Bs_cand.addUserFloat("k2k4_pt", (k2_p4 + k4_p4).pt());

// kaons
std::vector<std::string> kaon_names{"k1", "k2", "k3", "k4"};

for (size_t i = 0; i < kaon_names.size(); i++){
Bs_cand.addUserFloat(kaon_names[i] + "_fitted_mass" ,fitter_Bs.daughter_p4(i).mass());
Bs_cand.addUserFloat(kaon_names[i] + "_fitted_pt" ,fitter_Bs.daughter_p4(i).pt());
Bs_cand.addUserFloat(kaon_names[i] + "_fitted_eta",fitter_Bs.daughter_p4(i).eta());
Bs_cand.addUserFloat(kaon_names[i] + "_fitted_phi",fitter_Bs.daughter_p4(i).phi());
}

// apply pots-fit selection on phi1 candidate
if(!post_vtx_selection_(Bs_cand)) continue;

// save candidate
ret_value->push_back(Bs_cand);
}
}

evt.put(std::move(ret_value));
}

#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(BTo4TrkBuilder);
95 changes: 61 additions & 34 deletions BPHNano/plugins/DiTrackBuilder.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,6 @@




class DiTrackBuilder : public edm::global::EDProducer<> {


Expand Down Expand Up @@ -81,9 +80,10 @@ void DiTrackBuilder::produce(edm::StreamID, edm::Event &evt, edm::EventSetup con


// output
std::unique_ptr<pat::CompositeCandidateCollection> kstar_out(new pat::CompositeCandidateCollection());

std::unique_ptr<pat::CompositeCandidateCollection> cand_out(new pat::CompositeCandidateCollection());

// needed to sort in pt
std::vector<pat::CompositeCandidate> vector_candidates;

// main loop
for (size_t trk1_idx = 0; trk1_idx < pfcands->size(); ++trk1_idx ) {
Expand All @@ -98,74 +98,101 @@ void DiTrackBuilder::produce(edm::StreamID, edm::Event &evt, edm::EventSetup con
if (!trk2_selection_(*trk2_ptr)) continue;

bool UsedAgain = false;

std::vector< std::pair<double, double> > list_masses;
if (trk1_mass_ == trk2_mass_) {
list_masses.push_back(std::pair<double, double>(trk1_mass_, trk2_mass_));
}
else{
list_masses.push_back(std::pair<double, double>(trk1_mass_, trk2_mass_));
list_masses.push_back(std::pair<double, double>(trk2_mass_, trk1_mass_));
}

// Loop in all possible hypothesis
for ( std::pair<double, double> masses : { std::pair<double, double>(trk1_mass_, trk2_mass_), std::pair<double, double>(trk2_mass_, trk1_mass_) } ) {
// create a K* candidate; add first quantities that can be used for pre fit selection
pat::CompositeCandidate kstar_cand;
for ( std::pair<double, double> masses : list_masses ) {
// create a candidate; add first quantities that can be used for pre fit selection
pat::CompositeCandidate cand;

auto trk1_p4 = trk1_ptr->polarP4();
auto trk2_p4 = trk2_ptr->polarP4();
trk1_p4.SetM(masses.first);
trk2_p4.SetM(masses.second);

//adding stuff for pre fit selection
kstar_cand.setP4(trk1_p4 + trk2_p4);
kstar_cand.setCharge(trk1_ptr->charge() + trk2_ptr->charge());
kstar_cand.addUserFloat("trk_deltaR", reco::deltaR(*trk1_ptr, *trk2_ptr));
cand.setP4(trk1_p4 + trk2_p4);
cand.setCharge(trk1_ptr->charge() + trk2_ptr->charge());
cand.addUserFloat("trk_deltaR", reco::deltaR(*trk1_ptr, *trk2_ptr));

// save indices
kstar_cand.addUserInt("trk1_idx", trk1_idx );
kstar_cand.addUserInt("trk2_idx", trk2_idx );
kstar_cand.addUserFloat("trk1_mass", masses.first);
kstar_cand.addUserFloat("trk2_mass", masses.second);
cand.addUserInt("trk1_idx", trk1_idx );
cand.addUserInt("trk2_idx", trk2_idx );
cand.addUserFloat("trk1_mass", masses.first);
cand.addUserFloat("trk2_mass", masses.second);

// save cands
kstar_cand.addUserCand("trk1", trk1_ptr );
kstar_cand.addUserCand("trk2", trk2_ptr );
cand.addUserCand("trk1", trk1_ptr );
cand.addUserCand("trk2", trk2_ptr );

// selection before fit
if ( !pre_vtx_selection_(kstar_cand) ) continue;
if ( !pre_vtx_selection_(cand) ) continue;
//std::cout<<"trk1 "<<trk1_idx<<" trk2 "<<trk2_idx<<" dr "<< reco::deltaR(*trk1_ptr, *trk2_ptr)<<" pt1 "<<trk1_p4.pt()<<" pt2 "<<trk2_p4.pt()<<" mass "<<(trk1_p4+trk2_p4).mass()<<std::endl;

KinVtxFitter fitter(
{ttracks->at(trk1_idx), ttracks->at(trk2_idx)},
{ttracks->at(trk1_idx), ttracks->at(trk2_idx)},
{ masses.first, masses.second },
{K_SIGMA, K_SIGMA} //K and PI sigma equal...
);

if ( !fitter.success() ) continue;
kstar_cand.setVertex(
cand.setVertex(
reco::Candidate::Point(
fitter.fitted_vtx().x(),
fitter.fitted_vtx().y(),
fitter.fitted_vtx().z()
)
);
// save quantities after fit
kstar_cand.addUserInt("sv_ok", fitter.success() ? 1 : 0);
kstar_cand.addUserFloat("sv_chi2", fitter.chi2());
kstar_cand.addUserFloat("sv_ndof", fitter.dof());
kstar_cand.addUserFloat("sv_prob", fitter.prob());
kstar_cand.addUserFloat("fitted_mass", fitter.fitted_candidate().mass() );
kstar_cand.addUserFloat("fitted_pt",
cand.addUserInt("sv_ok", fitter.success() ? 1 : 0);
cand.addUserFloat("sv_chi2", fitter.chi2());
cand.addUserFloat("sv_ndof", fitter.dof());
cand.addUserFloat("sv_prob", fitter.prob());
cand.addUserFloat("fitted_mass", fitter.fitted_candidate().mass() );
cand.addUserFloat("fitted_pt",
fitter.fitted_candidate().globalMomentum().perp() );

kstar_cand.addUserFloat("fitted_eta",
cand.addUserFloat("fitted_eta",
fitter.fitted_candidate().globalMomentum().eta() );

kstar_cand.addUserFloat("fitted_phi",
cand.addUserFloat("fitted_phi",
fitter.fitted_candidate().globalMomentum().phi() );

kstar_cand.addUserInt("second_mass_hypothesis", UsedAgain );
kstar_cand.addUserFloat("vtx_x", kstar_cand.vx());
kstar_cand.addUserFloat("vtx_y", kstar_cand.vy());
kstar_cand.addUserFloat("vtx_z", kstar_cand.vz());
cand.addUserInt("second_mass_hypothesis", UsedAgain );
cand.addUserFloat("vtx_x", cand.vx());
cand.addUserFloat("vtx_y", cand.vy());
cand.addUserFloat("vtx_z", cand.vz());
cand.addUserFloat("deltaR_postfit", reco::deltaR(fitter.daughter_p4(0), fitter.daughter_p4(1)));
cand.addUserFloat("fitted_k1_pt", fitter.daughter_p4(0).pt());
cand.addUserFloat("fitted_k2_pt", fitter.daughter_p4(1).pt());

// after fit selection
if ( !post_vtx_selection_(kstar_cand) ) continue;
kstar_out->emplace_back(kstar_cand);
if ( !post_vtx_selection_(cand) ) continue;
vector_candidates.emplace_back(cand);
UsedAgain = true;
if (masses.first == masses.second) break;

} // end for ( auto & masses:
} // end for(size_t trk2_idx = trk1_idx + 1
} //for(size_t trk1_idx = 0

evt.put(std::move(kstar_out));
// sort candidate collection in pt
std::sort(vector_candidates.begin(), vector_candidates.end(),
[] (auto & cand1, auto & cand2) ->
bool {return (cand1.pt() > cand2.pt());}
);

for (auto & cand: vector_candidates){
cand_out->emplace_back(cand);
}
evt.put(std::move(cand_out));
}

#include "FWCore/Framework/interface/MakerMacros.h"
Expand Down
Loading