diff --git a/NanoAOD/plugins/MesonProducer.cc b/NanoAOD/plugins/MesonProducer.cc index dfad50c..3037153 100644 --- a/NanoAOD/plugins/MesonProducer.cc +++ b/NanoAOD/plugins/MesonProducer.cc @@ -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); @@ -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_; @@ -432,6 +438,8 @@ minDsPreselectMass_( iConfig.getParameter( "minDsPreselectMass" ) ), maxDsPreselectMass_( iConfig.getParameter( "maxDsPreselectMass" ) ), minDsMass_( iConfig.getParameter( "minDsMass" ) ), maxDsMass_( iConfig.getParameter( "maxDsMass" ) ), +minD0pi0PreselectMass_( iConfig.getParameter( "minD0pi0PreselectMass" ) ), +maxD0pi0PreselectMass_( iConfig.getParameter( "maxD0pi0PreselectMass" ) ), minD0PreselectMass_( iConfig.getParameter( "minD0PreselectMass" ) ), maxD0PreselectMass_( iConfig.getParameter( "maxD0PreselectMass" ) ), minD0Mass_( iConfig.getParameter( "minD0Mass" ) ), @@ -455,6 +463,7 @@ minVtxProb_( iConfig.getParameter( "minVtxProb" ) ) { produces("Ks"); produces("D0"); + produces("D0pi0"); produces("K0Star"); produces("Phi"); produces("Lambda"); @@ -890,6 +899,7 @@ void MesonProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { // Output collection auto kss = std::make_unique(); auto d0s = std::make_unique(); + auto d0pi0s = std::make_unique(); auto k0Stars = std::make_unique(); auto phis = std::make_unique(); auto lambdas = std::make_unique(); @@ -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 @@ -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"); @@ -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; @@ -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" ); @@ -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_); @@ -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.; @@ -1327,6 +1344,7 @@ 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" ); @@ -1334,34 +1352,157 @@ MesonProducer::getD0ToKPi(const edm::Event& iEvent, 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_ ); @@ -1372,6 +1513,7 @@ MesonProducer::getD0ToKPi(const edm::Event& iEvent, d0Cand.addUserFloat( "d0Star_photon2_phi", photon2_phi_ ); return d0Cand; + } @@ -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()=1 && abs(pdgId)<=5) || (abs(pdgId)>=11 && abs(pdgId)<=16) || pdgId==21 || pdgId==111 || abs(pdgId)==211 || abs(pdgId)==421 || abs(pdgId)==411 || (pdgId==22 && mass<1))?mass:0", float,precision="?((abs(pdgId)==6 || abs(pdgId)>1000000) && statusFlags().isLastCopy())?20:-1",doc="Mass stored for all particles with the exception of quarks (except top), leptons/neutrinos, photons with mass < 1 GeV, gluons, pi0(111), pi+(211), D0(421), and D+(411). For these particles, you can lookup the value from PDG."), + ) if triggerObjectTable.selections[1].id.value()==22: process.triggerObjectTable.selections[1].name=cms.string("Photon")