Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/main'
Browse files Browse the repository at this point in the history
  • Loading branch information
@mariadalfonso
mariadalfonso committed Aug 4, 2023
2 parents 13147d2 + e1e303a commit d06fd15
Show file tree
Hide file tree
Showing 3 changed files with 226 additions and 58 deletions.
253 changes: 198 additions & 55 deletions NanoAOD/plugins/MesonProducer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -266,6 +266,10 @@ class MesonProducer : public edm::EDProducer {
const pat::PackedCandidate& kaonCand,
const pat::PackedCandidate& pion);
pat::CompositeCandidate
getD0ToKPiPi0(const edm::Event& iEvent,
const pat::PackedCandidate& kaonCand,
const pat::PackedCandidate& pion);
pat::CompositeCandidate
getPhiToKK(const edm::Event& iEvent,
const pat::PackedCandidate& pfCand1,
const pat::PackedCandidate& pfCand2);
Expand Down Expand Up @@ -377,6 +381,8 @@ class MesonProducer : public edm::EDProducer {
double maxDsPreselectMass_;
double minDsMass_;
double maxDsMass_;
double minD0pi0PreselectMass_;
double maxD0pi0PreselectMass_;
double minD0PreselectMass_;
double maxD0PreselectMass_;
double minD0Mass_;
Expand Down Expand Up @@ -432,6 +438,8 @@ minDsPreselectMass_( iConfig.getParameter<double>( "minDsPreselectMass" ) ),
maxDsPreselectMass_( iConfig.getParameter<double>( "maxDsPreselectMass" ) ),
minDsMass_( iConfig.getParameter<double>( "minDsMass" ) ),
maxDsMass_( iConfig.getParameter<double>( "maxDsMass" ) ),
minD0pi0PreselectMass_( iConfig.getParameter<double>( "minD0pi0PreselectMass" ) ),
maxD0pi0PreselectMass_( iConfig.getParameter<double>( "maxD0pi0PreselectMass" ) ),
minD0PreselectMass_( iConfig.getParameter<double>( "minD0PreselectMass" ) ),
maxD0PreselectMass_( iConfig.getParameter<double>( "maxD0PreselectMass" ) ),
minD0Mass_( iConfig.getParameter<double>( "minD0Mass" ) ),
Expand All @@ -455,6 +463,7 @@ minVtxProb_( iConfig.getParameter<double>( "minVtxProb" ) )
{
produces<pat::CompositeCandidateCollection>("Ks");
produces<pat::CompositeCandidateCollection>("D0");
produces<pat::CompositeCandidateCollection>("D0pi0");
produces<pat::CompositeCandidateCollection>("K0Star");
produces<pat::CompositeCandidateCollection>("Phi");
produces<pat::CompositeCandidateCollection>("Lambda");
Expand Down Expand Up @@ -890,6 +899,7 @@ void MesonProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
// Output collection
auto kss = std::make_unique<pat::CompositeCandidateCollection>();
auto d0s = std::make_unique<pat::CompositeCandidateCollection>();
auto d0pi0s = std::make_unique<pat::CompositeCandidateCollection>();
auto k0Stars = std::make_unique<pat::CompositeCandidateCollection>();
auto phis = std::make_unique<pat::CompositeCandidateCollection>();
auto lambdas = std::make_unique<pat::CompositeCandidateCollection>();
Expand Down Expand Up @@ -970,6 +980,22 @@ void MesonProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
}
}

//D0->K-pi+pi0: Br=14% <-- this is the target
// D0ToKPi (1864) D0*ToD0Gammas (2007)
if ( pfCand1.charge() < 0 ) {
auto d0pi0Cand1 = getD0ToKPiPi0(iEvent, pfCand1, pfCand2);
if (d0pi0Cand1.numberOfDaughters() > 0 and d0pi0Cand1.pt() > 5){
d0pi0Cand1.addUserFloat( "doca", tt_doca);
d0pi0s->push_back(d0pi0Cand1);
}
} else if ( pfCand2.charge() < 0 ) {
auto d0pi0Cand2 = getD0ToKPiPi0(iEvent, pfCand2, pfCand1);
if (d0pi0Cand2.numberOfDaughters() > 0 and d0pi0Cand2.pt() > 5){
d0pi0Cand2.addUserFloat( "doca", tt_doca);
d0pi0s->push_back(d0pi0Cand2);
}
}

// Look for V0s built from displaced tracks

if ( not displacedTrack(pfCand1) or
Expand Down Expand Up @@ -1000,6 +1026,7 @@ void MesonProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {

iEvent.put(std::move(kss), "Ks");
iEvent.put(std::move(d0s), "D0");
iEvent.put(std::move(d0pi0s), "D0pi0");
iEvent.put(std::move(k0Stars), "K0Star");
iEvent.put(std::move(phis), "Phi");
iEvent.put(std::move(lambdas), "Lambda");
Expand Down Expand Up @@ -1128,10 +1155,6 @@ MesonProducer::getOmegasToPiPiPi0(const edm::Event& iEvent,
) return pat::CompositeCandidate();

// look for a Photon
float massOmegaFullCand = 0.f;
float ptOmegaFullCand = 0.f;
float etaOmegaFullCand = -99.f;
float phiOmegaFullCand = -99.f;
int nPhotons = 0;
float photon_pt_ = -1.f;
float photon_eta_ = 0.f;
Expand All @@ -1158,6 +1181,7 @@ MesonProducer::getOmegasToPiPiPi0(const edm::Event& iEvent,
photon.setMass(0.);
photon.setVertex(math::XYZPoint(vtx_point.x(), vtx_point.y(), vtx_point.z()));

// Omega: aim to get two photons from Pi0
nPhotons++;
if (nPhotons==1) omegaFullCand.addDaughter( photon , "photon" );
if (nPhotons==2) omegaFullCand.addDaughter( photon , "photon2" );
Expand All @@ -1166,41 +1190,38 @@ MesonProducer::getOmegasToPiPiPi0(const edm::Event& iEvent,

addP4.set(omegaFullCand);

LorentzVector tot_4;
tot_4.SetPxPyPzE(0,0,0,0);

unsigned int ndau = omegaFullCand.numberOfDaughters();

if ( omegaFullCand.mass() < minRhosPreselectMass_ or omegaFullCand.mass() > maxRhosPreselectMass_ )
return pat::CompositeCandidate();

tot_4 = pfCand1.p4()+pfCand2.p4();
TLorentzVector tot_4;
tot_4.SetPtEtaPhiM(omegasVtxFit.p3().perp(),omegasVtxFit.p3().eta(),omegasVtxFit.p3().phi(),omegasVtxFit.mass());

massOmegaFullCand = omegaFullCand.mass();
ptOmegaFullCand = omegaFullCand.pt();
etaOmegaFullCand = omegaFullCand.eta();
phiOmegaFullCand = omegaFullCand.phi();
unsigned int ndau = omegaFullCand.numberOfDaughters();

if (ndau>2) {
const reco::Candidate* ph1 = omegaFullCand.daughter("photon");
photon_pt_ = ph1->pt();
photon_eta_ = ph1->eta();
photon_phi_ = ph1->phi();
tot_4 = tot_4 + ph1->p4();
TLorentzVector ph1_;
ph1_.SetPtEtaPhiM(photon_pt_, photon_eta_, photon_phi_, 0.);
tot_4 = tot_4 + ph1_;
}
if (ndau>3) {
const reco::Candidate* ph2 = omegaFullCand.daughter("photon2");
photon2_pt_ = ph2->pt();
photon2_eta_ = ph2->eta();
photon2_phi_ = ph2->phi();
tot_4 = tot_4 + ph2->p4();
TLorentzVector ph2_;
ph2_.SetPtEtaPhiM(photon2_pt_, photon2_eta_, photon2_phi_, 0.);
tot_4 = tot_4 + ph2_;
}

omegasCand.addUserInt( "Nphotons", nPhotons );
omegasCand.addUserFloat( "Nbody_mass", massOmegaFullCand );
omegasCand.addUserFloat( "Nbody_pt", ptOmegaFullCand );
omegasCand.addUserFloat( "Nbody_eta", etaOmegaFullCand );
omegasCand.addUserFloat( "Nbody_phi", phiOmegaFullCand );
omegasCand.addUserFloat( "Nbody_mass", tot_4.M() );
omegasCand.addUserFloat( "Nbody_pt", tot_4.Pt() );
omegasCand.addUserFloat( "Nbody_eta", tot_4.Eta() );
omegasCand.addUserFloat( "Nbody_phi", tot_4.Phi() );

omegasCand.addUserFloat( "photon_pt", photon_pt_ );
omegasCand.addUserFloat( "photon_eta", photon_eta_);
Expand Down Expand Up @@ -1288,19 +1309,15 @@ MesonProducer::getD0ToKPi(const edm::Event& iEvent,

auto d0VtxFit = fillInfo(d0Cand, iEvent, kaon, pion, "kaon", "pion");

if ( not d0VtxFit.valid() or
if ( not d0VtxFit.valid() or
d0VtxFit.vtxProb() < minVtxProb_ or
d0VtxFit.mass() < minD0Mass_ or
d0VtxFit.mass() < minD0Mass_ or
d0VtxFit.mass() > maxD0Mass_ ) return pat::CompositeCandidate();
// d0VtxFit.lxy > maxLxy_ or
// d0VtxFit.sigLxy < minSigLxy_ ) return pat::CompositeCandidate();

// look for a Photon
// const pat::PackedCandidate* d0Star_photon(nullptr);
float massd0Star = -1.;
float ptd0Star = 0.;
float etad0Star = -99.;
float phid0Star = -99.;
int nPhotons = 0;
float photon_pt_ = 0.;
float photon_eta_ = -99.;
Expand All @@ -1327,41 +1344,165 @@ MesonProducer::getD0ToKPi(const edm::Event& iEvent,
photon.setMass(0.);
photon.setVertex(math::XYZPoint(vtx_point.x(), vtx_point.y(), vtx_point.z()));

// D0->k/pi; there the two photons are aiming to get one or two photons from Pi0 from the D0*->D0
nPhotons++;
if (nPhotons==1) d0StarCand.addDaughter( photon , "photon" );
if (nPhotons==2) d0StarCand.addDaughter( photon , "photon2" );
}

addP4.set( d0StarCand);

/// fixme
if ( d0StarCand.mass() > minD0StarMass_ and
d0StarCand.mass() < maxD0StarMass_ ) {
massd0Star = d0StarCand.mass();
ptd0Star = d0StarCand.pt();
etad0Star = d0StarCand.eta();
phid0Star = d0StarCand.phi();

unsigned int ndau = d0StarCand.numberOfDaughters();

if (ndau>2) {
const reco::Candidate* ph1 = d0StarCand.daughter("photon");
photon_pt_ = ph1->pt();
photon_eta_ = ph1->eta();
photon_phi_ = ph1->phi();
}
if (ndau>3) {
const reco::Candidate* ph2 = d0StarCand.daughter("photon2");
photon2_pt_ = ph2->pt();
photon2_eta_ = ph2->eta();
photon2_phi_ = ph2->phi();
}
// keep the lower bound to catch the cases with 0 photons
if ( d0StarCand.mass() < minD0Mass_ or d0StarCand.mass() > maxD0StarMass_ )
return pat::CompositeCandidate();

TLorentzVector tot_4;
tot_4.SetPtEtaPhiM(d0VtxFit.p3().perp(),d0VtxFit.p3().eta(),d0VtxFit.p3().phi(),d0VtxFit.mass());

unsigned int ndau = d0StarCand.numberOfDaughters();

if (ndau>2) {
const reco::Candidate* ph1 = d0StarCand.daughter("photon");
photon_pt_ = ph1->pt();
photon_eta_ = ph1->eta();
photon_phi_ = ph1->phi();
TLorentzVector ph1_;
ph1_.SetPtEtaPhiM(photon_pt_, photon_eta_, photon_phi_, 0.);
tot_4 = tot_4 + ph1_;
}
if (ndau>3) {
const reco::Candidate* ph2 = d0StarCand.daughter("photon2");
photon2_pt_ = ph2->pt();
photon2_eta_ = ph2->eta();
photon2_phi_ = ph2->phi();
TLorentzVector ph2_;
ph2_.SetPtEtaPhiM(photon2_pt_, photon2_eta_, photon2_phi_, 0.);
tot_4 = tot_4 + ph2_;
}

d0Cand.addUserFloat( "d0Star_Nbody_mass", tot_4.M() );
d0Cand.addUserFloat( "d0Star_Nbody_pt", tot_4.Pt() );
d0Cand.addUserFloat( "d0Star_Nbody_eta", tot_4.Eta() );
d0Cand.addUserFloat( "d0Star_Nbody_phi", tot_4.Phi() );

d0Cand.addUserInt( "d0Star_Nphotons", nPhotons );
d0Cand.addUserFloat( "d0Star_photon_pt", photon_pt_ );
d0Cand.addUserFloat( "d0Star_photon_eta", photon_eta_ );
d0Cand.addUserFloat( "d0Star_photon_phi", photon_phi_ );
d0Cand.addUserFloat( "d0Star_photon2_pt", photon2_pt_ );
d0Cand.addUserFloat( "d0Star_photon2_eta", photon2_eta_ );
d0Cand.addUserFloat( "d0Star_photon2_phi", photon2_phi_ );

return d0Cand;
}

pat::CompositeCandidate
MesonProducer::getD0ToKPiPi0(const edm::Event& iEvent,
const pat::PackedCandidate& ikaon,
const pat::PackedCandidate& ipion)
{
float signalCone_ = 0.05; // to collect the pi0
pat::CompositeCandidate d0Cand;
pat::PackedCandidate kaon(ikaon);
kaon.setMass(kaon_mass_);
pat::PackedCandidate pion(ipion);
pion.setMass(pion_mass_);
d0Cand.addDaughter( kaon , "kaon" );
d0Cand.addDaughter( pion , "pion" );
AddFourMomenta addP4;
addP4.set( d0Cand );

if ( d0Cand.mass() < minD0pi0PreselectMass_ or d0Cand.mass() > maxD0pi0PreselectMass_ )
return pat::CompositeCandidate();

d0Cand.addUserFloat( "kaon_pt", kaon.pt() );
d0Cand.addUserFloat( "kaon_eta", kaon.eta() );
d0Cand.addUserFloat( "kaon_phi", kaon.phi() );
d0Cand.addUserFloat( "kaon_charge", kaon.charge() );
d0Cand.addUserFloat( "pion_pt", pion.pt() );
d0Cand.addUserFloat( "pion_eta", pion.eta() );
d0Cand.addUserFloat( "pion_phi", pion.phi() );
d0Cand.addUserFloat( "pion_charge", pion.charge() );
d0Cand.addUserFloat( "kaon_sip", trackImpactParameterSignificance(kaon) );
d0Cand.addUserFloat( "pion_sip", trackImpactParameterSignificance(pion) );
// d0Cand.addUserInt( "kaon_mu_index", match_to_muon(kaon, *muonHandle_));
// d0Cand.addUserInt( "pion_mu_index", match_to_muon(pion, *muonHandle_));

auto d0VtxFit = fillInfo(d0Cand, iEvent, kaon, pion, "kaon", "pion");

if ( not d0VtxFit.valid() or
d0VtxFit.vtxProb() < minVtxProb_ or
d0VtxFit.mass() < minD0pi0PreselectMass_ or
d0VtxFit.mass() > maxD0pi0PreselectMass_ ) return pat::CompositeCandidate();
// d0VtxFit.lxy > maxLxy_ or
// d0VtxFit.sigLxy < minSigLxy_ ) return pat::CompositeCandidate();

// look for a Photon
// const pat::PackedCandidate* d0Star_photon(nullptr);
int nPhotons = 0;
float photon_pt_ = 0.;
float photon_eta_ = -99.;
float photon_phi_ = -99.;
float photon2_pt_ = 0.;
float photon2_eta_ = -99.;
float photon2_phi_ = -99.;

const auto & vtx_point = d0VtxFit.refitVertex->vertexState().position();
pat::CompositeCandidate d0StarCand;
d0StarCand.addDaughter( kaon , "kaon" );
d0StarCand.addDaughter( pion , "pion" );

for (unsigned int k=0; k < pfCandHandle_->size(); ++k){
auto iphoton(pfCandHandle_->at(k));

if (iphoton.charge() != 0 ) continue;
if (iphoton.mass() > 1 ) continue; // some uninitialized mass
if (abs(iphoton.pdgId()) !=22) continue; // otherwise some "KL" enters
if (iphoton.pt() < 5. ) continue;
if (deltaR((kaon.p4()+pion.p4()),iphoton.p4()) > signalCone_) continue; // photon should be collimated

pat::PackedCandidate photon(iphoton);
photon.setMass(0.);
photon.setVertex(math::XYZPoint(vtx_point.x(), vtx_point.y(), vtx_point.z()));

nPhotons++;
if (nPhotons==1) d0StarCand.addDaughter( photon , "photon" );
if (nPhotons==2) d0StarCand.addDaughter( photon , "photon2" );
}

d0Cand.addUserFloat( "d0Star_Nbody_mass", massd0Star );
d0Cand.addUserFloat( "d0Star_Nbody_pt", ptd0Star );
d0Cand.addUserFloat( "d0Star_Nbody_eta", etad0Star );
d0Cand.addUserFloat( "d0Star_Nbody_phi", phid0Star );
addP4.set( d0StarCand);

TLorentzVector tot_4;
tot_4.SetPtEtaPhiM(d0VtxFit.p3().perp(),d0VtxFit.p3().eta(),d0VtxFit.p3().phi(),d0VtxFit.mass());

if ( d0StarCand.mass() < minD0StarMass_ and d0StarCand.mass() > maxD0StarMass_ )
return pat::CompositeCandidate();

unsigned int ndau = d0StarCand.numberOfDaughters();

if (ndau>2) {
const reco::Candidate* ph1 = d0StarCand.daughter("photon");
photon_pt_ = ph1->pt();
photon_eta_ = ph1->eta();
photon_phi_ = ph1->phi();
TLorentzVector ph1_;
ph1_.SetPtEtaPhiM(photon_pt_, photon_eta_, photon_phi_, 0.);
tot_4 = tot_4 + ph1_;
}
if (ndau>3) {
const reco::Candidate* ph2 = d0StarCand.daughter("photon2");
photon2_pt_ = ph2->pt();
photon2_eta_ = ph2->eta();
photon2_phi_ = ph2->phi();
TLorentzVector ph2_;
ph2_.SetPtEtaPhiM(photon2_pt_, photon2_eta_, photon2_phi_, 0.);
tot_4 = tot_4 + ph2_;
}

d0Cand.addUserFloat( "d0Star_Nbody_mass", tot_4.M() );
d0Cand.addUserFloat( "d0Star_Nbody_pt", tot_4.Pt() );
d0Cand.addUserFloat( "d0Star_Nbody_eta", tot_4.Eta() );
d0Cand.addUserFloat( "d0Star_Nbody_phi", tot_4.Phi() );

d0Cand.addUserInt( "d0Star_Nphotons", nPhotons );
d0Cand.addUserFloat( "d0Star_photon_pt", photon_pt_ );
Expand All @@ -1372,6 +1513,7 @@ MesonProducer::getD0ToKPi(const edm::Event& iEvent,
d0Cand.addUserFloat( "d0Star_photon2_phi", photon2_phi_ );

return d0Cand;

}


Expand Down Expand Up @@ -1693,10 +1835,11 @@ MesonProducer::computeCandIsolation(const pat::PackedCandidate& pfCand1, const p
if (particleType==0) {
if (pfCandIso.charge() == 0 ) continue; // for chargedIsolation
if (!pfCandIso.hasTrackDetails()) continue;
} else if (particleType==22) {
if (pfCandIso.charge() != 0 ) continue; // for neutralIsolation MARIA
} else if (particleType==22 and pfCandIso.pdgId()==22) {
if (pfCandIso.charge() != 0 ) continue; // only photons for neutralIsolation
} else {
if (pfCandIso.charge() != 0 ) continue; // for neutralIsolation MARIA
if (pfCandIso.pdgId() == 22) continue; // only neutral hadrons
if (pfCandIso.charge() != 0 ) continue; // for neutralIsolation
}
if (pfCandIso.pt()<minPt) continue;
if (pfCandIso.vertexRef().key()!=primaryVertexIndex) continue;
Expand Down
Loading

0 comments on commit d06fd15

Please sign in to comment.