From 5717ffbc7446b705ad8bdce4e968272043f314ad Mon Sep 17 00:00:00 2001 From: Rocco Ardino Date: Wed, 11 Sep 2024 00:51:43 +0200 Subject: [PATCH] Add H->Q+gamma and H->QQ analyses --- .../Phase2/plugins/ScPhase2PuppiH2PhiDemo.cc | 300 +++++++++++++++ .../Phase2/plugins/ScPhase2PuppiH2RhoDemo.cc | 300 +++++++++++++++ .../plugins/ScPhase2PuppiHJPsiGammaDemo.cc | 333 ++++++++++++++++ .../plugins/ScPhase2PuppiHPhiGammaDemo.cc | 335 +++++++++++++++++ .../plugins/ScPhase2PuppiHPhiJPsiDemo.cc | 300 +++++++++++++++ .../plugins/ScPhase2PuppiHRhoGammaDemo.cc | 335 +++++++++++++++++ .../Phase2/test/runScoutingPhase2HRare_cfg.py | 355 ++++++++++++++++++ 7 files changed, 2258 insertions(+) create mode 100644 L1TriggerScouting/Phase2/plugins/ScPhase2PuppiH2PhiDemo.cc create mode 100644 L1TriggerScouting/Phase2/plugins/ScPhase2PuppiH2RhoDemo.cc create mode 100644 L1TriggerScouting/Phase2/plugins/ScPhase2PuppiHJPsiGammaDemo.cc create mode 100644 L1TriggerScouting/Phase2/plugins/ScPhase2PuppiHPhiGammaDemo.cc create mode 100644 L1TriggerScouting/Phase2/plugins/ScPhase2PuppiHPhiJPsiDemo.cc create mode 100644 L1TriggerScouting/Phase2/plugins/ScPhase2PuppiHRhoGammaDemo.cc create mode 100644 L1TriggerScouting/Phase2/test/runScoutingPhase2HRare_cfg.py diff --git a/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiH2PhiDemo.cc b/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiH2PhiDemo.cc new file mode 100644 index 0000000000000..cb7351e0232ec --- /dev/null +++ b/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiH2PhiDemo.cc @@ -0,0 +1,300 @@ +#include +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/StreamID.h" + +#include "DataFormats/L1Scouting/interface/OrbitCollection.h" +#include "DataFormats/L1Scouting/interface/OrbitFlatTable.h" +#include "DataFormats/L1TParticleFlow/interface/L1ScoutingPuppi.h" +#include "DataFormats/L1TParticleFlow/interface/L1ScoutingTkEm.h" +#include "L1TriggerScouting/Utilities/interface/BxOffsetsFiller.h" + +#include +#include +#include +#include +#include +#include +#include + +class ScPhase2PuppiH2PhiDemo : public edm::stream::EDProducer<> { +public: + explicit ScPhase2PuppiH2PhiDemo(const edm::ParameterSet &); + ~ScPhase2PuppiH2PhiDemo() override; + static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); + +private: + void beginStream(edm::StreamID) override; + void produce(edm::Event &, const edm::EventSetup &) override; + void endStream() override; + template + void runObj(const OrbitCollection &src, + edm::Event &out, + unsigned long &nTry, + unsigned long &nPass, + const std::string &bxLabel); + + bool doStruct_; + edm::EDGetTokenT> structToken_; + + struct Cuts { + float minptD = 5; + float minptQ = 10; + float maxdeltarD2 = 0.40 * 0.40; + float minmassH = 100; + float maxmassH = 150; + float minmassQ = 0.95; + float maxmassQ = 1.25; + float mindr2 = 0.05 * 0.05; + float maxdr2 = 0.25 * 0.25; + float maxiso = 0.25; + } cuts; + + template + bool isolationQ( + unsigned int pidex1, unsigned int pidex2, const T *cands, unsigned int size) const; + + std::tuple deltar(float eta1, float eta2, float phi1, float phi2) const; + + template + static float pairmass(const std::array &t, const T *cands, const std::array &massD); + + template + static float quadrupletmass(const std::array &t, const T *cands, const std::array &massD); + + unsigned long countStruct_; + unsigned long passStruct_; +}; + +ScPhase2PuppiH2PhiDemo::ScPhase2PuppiH2PhiDemo(const edm::ParameterSet &iConfig) + : doStruct_(iConfig.getParameter("runStruct")) { + if (doStruct_) { + structToken_ = consumes>(iConfig.getParameter("src")); + produces>("selectedBx"); + produces("h2phi"); + } +} + +ScPhase2PuppiH2PhiDemo::~ScPhase2PuppiH2PhiDemo(){}; + +void ScPhase2PuppiH2PhiDemo::beginStream(edm::StreamID) { + countStruct_ = 0; + passStruct_ = 0; +} + +void ScPhase2PuppiH2PhiDemo::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) { + if (doStruct_) { + edm::Handle> src; + iEvent.getByToken(structToken_, src); + + runObj(*src, iEvent, countStruct_, passStruct_, ""); + } +} + +void ScPhase2PuppiH2PhiDemo::endStream() { + if (doStruct_) + std::cout << "H2Phi Struct analysis: " << countStruct_ << " -> " << passStruct_ << std::endl; +} + +template +void ScPhase2PuppiH2PhiDemo::runObj(const OrbitCollection &src, + edm::Event &iEvent, + unsigned long &nTry, + unsigned long &nPass, + const std::string &label) { + l1ScoutingRun3::BxOffsetsFillter bxOffsetsFiller; + bxOffsetsFiller.start(); + auto ret = std::make_unique>(); + std::vector masses; + std::vector i0s, i1s, i2s, i3s; + ROOT::RVec ix; // + std::array bestPair1, bestPair2; + bool bestPair1Found, bestPair2Found; + float bestPair1Score, bestPair2Score; + for (unsigned int bx = 1; bx <= OrbitCollection::NBX; ++bx) { + nTry++; + auto range = src.bxIterator(bx); + const T *cands = &range.front(); + auto size = range.size(); + + ix.clear(); + for (unsigned int i = 0; i < size; ++i) { //make list of all hadrons + if ((std::abs(cands[i].pdgId()) == 211 or std::abs(cands[i].pdgId()) == 11)) { + if (cands[i].pt() >= cuts.minptD) + ix.push_back(i); + } + } + unsigned int ndaus = ix.size(); + if (ndaus < 4) + continue; + + // Q1 candidate from closest OS pair with mass compatible with mQ + bestPair1Found = false; + bestPair1Score = 999; + for (unsigned int i1 = 0; i1 < ndaus; ++i1) { + if (cands[ix[i1]].pt() < cuts.minptD) + continue; // D1 pt cut + for (unsigned int i2 = 0; i2 < ndaus; ++i2) { + if (i2 == i1 || cands[ix[i2]].pt() < cuts.minptD) + continue; // D2 pt cut + + if (!(cands[ix[i1]].charge()*cands[ix[i2]].charge() < 0)) + continue; + + auto mass2 = pairmass({{ix[i1], ix[i2]}}, cands, {{0.4937, 0.4937}}); + if (mass2 >= cuts.minmassQ and mass2 <= cuts.maxmassQ) + continue; + + auto [drcond, drQ] = deltar(cands[ix[i1]].eta(), cands[ix[i2]].eta(), cands[ix[i1]].phi(), cands[ix[i2]].phi()); + if (!drcond) + continue; // angular sep of top 2 tracks + + std::array pair{{ix[i1], ix[i2]}}; // pair of indices + if (drQ < bestPair1Score) { + std::copy_n(pair.begin(), 2, bestPair1.begin()); + bestPair1Score = drQ; + if (bestPair1Score*bestPair1Score < cuts.maxdeltarD2) + bestPair1Found = true; + } + } + } + if (!bestPair1Found) + continue; // pair was found + auto ptQ = (cands[bestPair1[0]].p4() + cands[bestPair1[1]].p4()).pt(); + if (ptQ < cuts.minptQ) + continue; // Q pt + if (!isolationQ(bestPair1[0], bestPair1[1], cands, size)) + continue; // Q isolation + + // Q2 candidate from closest OS pair with mass compatible with mQ + bestPair2Found = false; + bestPair2Score = 999; + for (unsigned int i3 = 0; i3 < ndaus; ++i3) { + if (cands[ix[i3]].pt() < cuts.minptD) + continue; // D1 pt cut + if (ix[i3]==bestPair1[0] or ix[i3]==bestPair1[1]) + continue; // don't reuse candidates from previous pair + for (unsigned int i4 = 0; i4 < ndaus; ++i4) { + if (i4 == i3 || cands[ix[i4]].pt() < cuts.minptD) + continue; // D2 pt cut + if (ix[i4]==bestPair1[0] or ix[i4]==bestPair1[1]) + continue; // don't reuse candidates from previous pair + if (!(cands[ix[i3]].charge()*cands[ix[i4]].charge() < 0)) + continue; // OS pair + auto mass2 = pairmass({{ix[i3], ix[i4]}}, cands, {{0.4937, 0.4937}}); // (cands[ix[i3]].p4() + cands[ix[i4]].p4()).mass(); + if (mass2 >= cuts.minmassQ and mass2 <= cuts.maxmassQ) + continue; // Q mass + auto [drcond, drQ] = deltar(cands[ix[i3]].eta(), cands[ix[i4]].eta(), cands[ix[i3]].phi(), cands[ix[i4]].phi()); + if (!drcond) + continue; // angular sep of top 2 tracks + + std::array pair{{ix[i3], ix[i4]}}; // pair of indices + if (drQ < bestPair2Score) { + std::copy_n(pair.begin(), 2, bestPair2.begin()); + bestPair2Score = drQ; + if (bestPair2Score*bestPair2Score < cuts.maxdeltarD2) + bestPair2Found = true; + } + } + } + if (!bestPair2Found) + continue; // pair was found + ptQ = (cands[bestPair2[0]].p4() + cands[bestPair2[1]].p4()).pt(); + if (ptQ < cuts.minptQ) + continue; // Q pt + if (!isolationQ(bestPair2[0], bestPair2[1], cands, size)) + continue; // Q isolation + + std::array bestQuadruplet{{bestPair1[0], bestPair1[1], bestPair2[0], bestPair2[1]}}; + // H mass + auto mass = quadrupletmass(bestQuadruplet, cands, {{0.4937, 0.4937, 0.4937, 0.4937}}); + if (!(mass >= cuts.minmassH and mass <= cuts.maxmassH)) + continue; + + ret->emplace_back(bx); + nPass++; + masses.push_back(mass); + i0s.push_back(bestQuadruplet[0]); + i1s.push_back(bestQuadruplet[1]); + i2s.push_back(bestQuadruplet[2]); + i3s.push_back(bestQuadruplet[3]); + bxOffsetsFiller.addBx(bx, 1); + } // loop on BXs + + iEvent.put(std::move(ret), "selectedBx" + label); + // now we make the table + auto bxOffsets = bxOffsetsFiller.done(); + auto tab = std::make_unique(bxOffsets, "H2Phi" + label, true); + tab->addColumn("mass", masses, "4 kaons invariant mass"); + tab->addColumn("i0", i0s, "1st kaon (phi1)"); + tab->addColumn("i1", i1s, "2nd kaon (phi1)"); + tab->addColumn("i2", i2s, "1st kaon (phi2)"); + tab->addColumn("i3", i3s, "2nd kaon (phi2)"); + iEvent.put(std::move(tab), "h2phi" + label); +} + +//TEST functions +template +bool ScPhase2PuppiH2PhiDemo::isolationQ( + unsigned int pidex1, unsigned int pidex2, const T *cands, unsigned int size) const { + bool passed = false; + float psum = 0; + float eta = cands[pidex1].eta(); //center cone around leading track + float phi = cands[pidex1].phi(); + for (unsigned int j = 0u; j < size; ++j) { //loop over other particles + if (pidex1 == j or pidex2 == j) + continue; + float deta = eta - cands[j].eta(), dphi = ROOT::VecOps::DeltaPhi(phi, cands[j].phi()); + float dr2 = deta * deta + dphi * dphi; + if (dr2 >= cuts.mindr2 && dr2 <= cuts.maxdr2) + psum += cands[j].pt(); + } + if (psum <= cuts.maxiso * (cands[pidex1].pt() + cands[pidex2].pt())) + passed = true; + return passed; +} + +std::tuple ScPhase2PuppiH2PhiDemo::deltar(float eta1, float eta2, float phi1, float phi2) const { + bool passed = true; + float deta = eta1 - eta2; + float dphi = ROOT::VecOps::DeltaPhi(phi1, phi2); + float dr2 = deta * deta + dphi * dphi; + if (dr2 > cuts.maxdeltarD2) { + passed = false; + return std::tuple(passed, dr2); + } + return std::tuple(passed, dr2); +} + +template +float ScPhase2PuppiH2PhiDemo::pairmass(const std::array &t, + const T *cands, + const std::array &massD) { + ROOT::Math::PtEtaPhiMVector p1(cands[t[0]].pt(), cands[t[0]].eta(), cands[t[0]].phi(), massD[0]); + ROOT::Math::PtEtaPhiMVector p2(cands[t[1]].pt(), cands[t[1]].eta(), cands[t[1]].phi(), massD[1]); + float mass = (p1 + p2).M(); + return mass; +} + +template +float ScPhase2PuppiH2PhiDemo::quadrupletmass(const std::array &t, + const T *cands, + const std::array &massD) { + ROOT::Math::PtEtaPhiMVector p1(cands[t[0]].pt(), cands[t[0]].eta(), cands[t[0]].phi(), massD[0]); + ROOT::Math::PtEtaPhiMVector p2(cands[t[1]].pt(), cands[t[1]].eta(), cands[t[1]].phi(), massD[1]); + ROOT::Math::PtEtaPhiMVector p3(cands[t[2]].pt(), cands[t[2]].eta(), cands[t[2]].phi(), massD[2]); + ROOT::Math::PtEtaPhiMVector p4(cands[t[3]].pt(), cands[t[3]].eta(), cands[t[3]].phi(), massD[3]); + float mass = (p1 + p2 + p3 + p4).M(); + return mass; +} + +void ScPhase2PuppiH2PhiDemo::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + +DEFINE_FWK_MODULE(ScPhase2PuppiH2PhiDemo); diff --git a/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiH2RhoDemo.cc b/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiH2RhoDemo.cc new file mode 100644 index 0000000000000..54d8f235f75ad --- /dev/null +++ b/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiH2RhoDemo.cc @@ -0,0 +1,300 @@ +#include +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/StreamID.h" + +#include "DataFormats/L1Scouting/interface/OrbitCollection.h" +#include "DataFormats/L1Scouting/interface/OrbitFlatTable.h" +#include "DataFormats/L1TParticleFlow/interface/L1ScoutingPuppi.h" +#include "DataFormats/L1TParticleFlow/interface/L1ScoutingTkEm.h" +#include "L1TriggerScouting/Utilities/interface/BxOffsetsFiller.h" + +#include +#include +#include +#include +#include +#include +#include + +class ScPhase2PuppiH2RhoDemo : public edm::stream::EDProducer<> { +public: + explicit ScPhase2PuppiH2RhoDemo(const edm::ParameterSet &); + ~ScPhase2PuppiH2RhoDemo() override; + static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); + +private: + void beginStream(edm::StreamID) override; + void produce(edm::Event &, const edm::EventSetup &) override; + void endStream() override; + template + void runObj(const OrbitCollection &src, + edm::Event &out, + unsigned long &nTry, + unsigned long &nPass, + const std::string &bxLabel); + + bool doStruct_; + edm::EDGetTokenT> structToken_; + + struct Cuts { + float minptD = 5; + float minptQ = 10; + float maxdeltarD2 = 0.40 * 0.40; + float minmassH = 100; + float maxmassH = 150; + float minmassQ = 0.40; + float maxmassQ = 1.30; + float mindr2 = 0.05 * 0.05; + float maxdr2 = 0.25 * 0.25; + float maxiso = 0.25; + } cuts; + + template + bool isolationQ( + unsigned int pidex1, unsigned int pidex2, const T *cands, unsigned int size) const; + + std::tuple deltar(float eta1, float eta2, float phi1, float phi2) const; + + template + static float pairmass(const std::array &t, const T *cands, const std::array &massD); + + template + static float quadrupletmass(const std::array &t, const T *cands, const std::array &massD); + + unsigned long countStruct_; + unsigned long passStruct_; +}; + +ScPhase2PuppiH2RhoDemo::ScPhase2PuppiH2RhoDemo(const edm::ParameterSet &iConfig) + : doStruct_(iConfig.getParameter("runStruct")) { + if (doStruct_) { + structToken_ = consumes>(iConfig.getParameter("src")); + produces>("selectedBx"); + produces("h2rho"); + } +} + +ScPhase2PuppiH2RhoDemo::~ScPhase2PuppiH2RhoDemo(){}; + +void ScPhase2PuppiH2RhoDemo::beginStream(edm::StreamID) { + countStruct_ = 0; + passStruct_ = 0; +} + +void ScPhase2PuppiH2RhoDemo::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) { + if (doStruct_) { + edm::Handle> src; + iEvent.getByToken(structToken_, src); + + runObj(*src, iEvent, countStruct_, passStruct_, ""); + } +} + +void ScPhase2PuppiH2RhoDemo::endStream() { + if (doStruct_) + std::cout << "H2Rho Struct analysis: " << countStruct_ << " -> " << passStruct_ << std::endl; +} + +template +void ScPhase2PuppiH2RhoDemo::runObj(const OrbitCollection &src, + edm::Event &iEvent, + unsigned long &nTry, + unsigned long &nPass, + const std::string &label) { + l1ScoutingRun3::BxOffsetsFillter bxOffsetsFiller; + bxOffsetsFiller.start(); + auto ret = std::make_unique>(); + std::vector masses; + std::vector i0s, i1s, i2s, i3s; + ROOT::RVec ix; // + std::array bestPair1, bestPair2; + bool bestPair1Found, bestPair2Found; + float bestPair1Score, bestPair2Score; + for (unsigned int bx = 1; bx <= OrbitCollection::NBX; ++bx) { + nTry++; + auto range = src.bxIterator(bx); + const T *cands = &range.front(); + auto size = range.size(); + + ix.clear(); + for (unsigned int i = 0; i < size; ++i) { //make list of all hadrons + if ((std::abs(cands[i].pdgId()) == 211 or std::abs(cands[i].pdgId()) == 11)) { + if (cands[i].pt() >= cuts.minptD) + ix.push_back(i); + } + } + unsigned int ndaus = ix.size(); + if (ndaus < 4) + continue; + + // Q1 candidate from closest OS pair with mass compatible with mQ + bestPair1Found = false; + bestPair1Score = 999; + for (unsigned int i1 = 0; i1 < ndaus; ++i1) { + if (cands[ix[i1]].pt() < cuts.minptD) + continue; // D1 pt cut + for (unsigned int i2 = 0; i2 < ndaus; ++i2) { + if (i2 == i1 || cands[ix[i2]].pt() < cuts.minptD) + continue; // D2 pt cut + + if (!(cands[ix[i1]].charge()*cands[ix[i2]].charge() < 0)) + continue; + + auto mass2 = pairmass({{ix[i1], ix[i2]}}, cands, {{0.1396, 0.1396}}); + if (mass2 >= cuts.minmassQ and mass2 <= cuts.maxmassQ) + continue; + + auto [drcond, drQ] = deltar(cands[ix[i1]].eta(), cands[ix[i2]].eta(), cands[ix[i1]].phi(), cands[ix[i2]].phi()); + if (!drcond) + continue; // angular sep of top 2 tracks + + std::array pair{{ix[i1], ix[i2]}}; // pair of indices + if (drQ < bestPair1Score) { + std::copy_n(pair.begin(), 2, bestPair1.begin()); + bestPair1Score = drQ; + if (bestPair1Score*bestPair1Score < cuts.maxdeltarD2) + bestPair1Found = true; + } + } + } + if (!bestPair1Found) + continue; // pair was found + auto ptQ = (cands[bestPair1[0]].p4() + cands[bestPair1[1]].p4()).pt(); + if (ptQ < cuts.minptQ) + continue; // Q pt + if (!isolationQ(bestPair1[0], bestPair1[1], cands, size)) + continue; // Q isolation + + // Q2 candidate from closest OS pair with mass compatible with mQ + bestPair2Found = false; + bestPair2Score = 999; + for (unsigned int i3 = 0; i3 < ndaus; ++i3) { + if (cands[ix[i3]].pt() < cuts.minptD) + continue; // D1 pt cut + if (ix[i3]==bestPair1[0] or ix[i3]==bestPair1[1]) + continue; // don't reuse candidates from previous pair + for (unsigned int i4 = 0; i4 < ndaus; ++i4) { + if (i4 == i3 || cands[ix[i4]].pt() < cuts.minptD) + continue; // D2 pt cut + if (ix[i4]==bestPair1[0] or ix[i4]==bestPair1[1]) + continue; // don't reuse candidates from previous pair + if (!(cands[ix[i3]].charge()*cands[ix[i4]].charge() < 0)) + continue; // OS pair + auto mass2 = pairmass({{ix[i3], ix[i4]}}, cands, {{0.1396, 0.1396}}); + if (mass2 >= cuts.minmassQ and mass2 <= cuts.maxmassQ) + continue; // Q mass + auto [drcond, drQ] = deltar(cands[ix[i3]].eta(), cands[ix[i4]].eta(), cands[ix[i3]].phi(), cands[ix[i4]].phi()); + if (!drcond) + continue; // angular sep of top 2 tracks + + std::array pair{{ix[i3], ix[i4]}}; // pair of indices + if (drQ < bestPair2Score) { + std::copy_n(pair.begin(), 2, bestPair2.begin()); + bestPair2Score = drQ; + if (bestPair2Score*bestPair2Score < cuts.maxdeltarD2) + bestPair2Found = true; + } + } + } + if (!bestPair2Found) + continue; // pair was found + ptQ = (cands[bestPair2[0]].p4() + cands[bestPair2[1]].p4()).pt(); + if (ptQ < cuts.minptQ) + continue; // Q pt + if (!isolationQ(bestPair2[0], bestPair2[1], cands, size)) + continue; // Q isolation + + std::array bestQuadruplet{{bestPair1[0], bestPair1[1], bestPair2[0], bestPair2[1]}}; + // H mass + auto mass = quadrupletmass(bestQuadruplet, cands, {{0.1396, 0.1396, 0.1396, 0.1396}}); + if (!(mass >= cuts.minmassH and mass <= cuts.maxmassH)) + continue; + + ret->emplace_back(bx); + nPass++; + masses.push_back(mass); + i0s.push_back(bestQuadruplet[0]); + i1s.push_back(bestQuadruplet[1]); + i2s.push_back(bestQuadruplet[2]); + i3s.push_back(bestQuadruplet[3]); + bxOffsetsFiller.addBx(bx, 1); + } // loop on BXs + + iEvent.put(std::move(ret), "selectedBx" + label); + // now we make the table + auto bxOffsets = bxOffsetsFiller.done(); + auto tab = std::make_unique(bxOffsets, "H2Rho" + label, true); + tab->addColumn("mass", masses, "4 pions invariant mass"); + tab->addColumn("i0", i0s, "1st pion (rho1)"); + tab->addColumn("i1", i1s, "2nd pion (rho1)"); + tab->addColumn("i2", i2s, "1st pion (rho2)"); + tab->addColumn("i3", i3s, "2nd pion (rho2)"); + iEvent.put(std::move(tab), "h2rho" + label); +} + +//TEST functions +template +bool ScPhase2PuppiH2RhoDemo::isolationQ( + unsigned int pidex1, unsigned int pidex2, const T *cands, unsigned int size) const { + bool passed = false; + float psum = 0; + float eta = cands[pidex1].eta(); //center cone around leading track + float phi = cands[pidex1].phi(); + for (unsigned int j = 0u; j < size; ++j) { //loop over other particles + if (pidex1 == j or pidex2 == j) + continue; + float deta = eta - cands[j].eta(), dphi = ROOT::VecOps::DeltaPhi(phi, cands[j].phi()); + float dr2 = deta * deta + dphi * dphi; + if (dr2 >= cuts.mindr2 && dr2 <= cuts.maxdr2) + psum += cands[j].pt(); + } + if (psum <= cuts.maxiso * (cands[pidex1].pt() + cands[pidex2].pt())) + passed = true; + return passed; +} + +std::tuple ScPhase2PuppiH2RhoDemo::deltar(float eta1, float eta2, float phi1, float phi2) const { + bool passed = true; + float deta = eta1 - eta2; + float dphi = ROOT::VecOps::DeltaPhi(phi1, phi2); + float dr2 = deta * deta + dphi * dphi; + if (dr2 > cuts.maxdeltarD2) { + passed = false; + return std::tuple(passed, dr2); + } + return std::tuple(passed, dr2); +} + +template +float ScPhase2PuppiH2RhoDemo::pairmass(const std::array &t, + const T *cands, + const std::array &massD) { + ROOT::Math::PtEtaPhiMVector p1(cands[t[0]].pt(), cands[t[0]].eta(), cands[t[0]].phi(), massD[0]); + ROOT::Math::PtEtaPhiMVector p2(cands[t[1]].pt(), cands[t[1]].eta(), cands[t[1]].phi(), massD[1]); + float mass = (p1 + p2).M(); + return mass; +} + +template +float ScPhase2PuppiH2RhoDemo::quadrupletmass(const std::array &t, + const T *cands, + const std::array &massD) { + ROOT::Math::PtEtaPhiMVector p1(cands[t[0]].pt(), cands[t[0]].eta(), cands[t[0]].phi(), massD[0]); + ROOT::Math::PtEtaPhiMVector p2(cands[t[1]].pt(), cands[t[1]].eta(), cands[t[1]].phi(), massD[1]); + ROOT::Math::PtEtaPhiMVector p3(cands[t[2]].pt(), cands[t[2]].eta(), cands[t[2]].phi(), massD[2]); + ROOT::Math::PtEtaPhiMVector p4(cands[t[3]].pt(), cands[t[3]].eta(), cands[t[3]].phi(), massD[3]); + float mass = (p1 + p2 + p3 + p4).M(); + return mass; +} + +void ScPhase2PuppiH2RhoDemo::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + +DEFINE_FWK_MODULE(ScPhase2PuppiH2RhoDemo); diff --git a/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiHJPsiGammaDemo.cc b/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiHJPsiGammaDemo.cc new file mode 100644 index 0000000000000..a6c214494b1a6 --- /dev/null +++ b/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiHJPsiGammaDemo.cc @@ -0,0 +1,333 @@ +#include +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/StreamID.h" + +#include "DataFormats/L1Scouting/interface/OrbitCollection.h" +#include "DataFormats/L1Scouting/interface/OrbitFlatTable.h" +#include "DataFormats/L1TParticleFlow/interface/L1ScoutingPuppi.h" +#include "DataFormats/L1TParticleFlow/interface/L1ScoutingTkEm.h" +#include "L1TriggerScouting/Utilities/interface/BxOffsetsFiller.h" + +#include +#include +#include +#include +#include +#include +#include + +class ScPhase2PuppiHJPsiGammaDemo : public edm::stream::EDProducer<> { +public: + explicit ScPhase2PuppiHJPsiGammaDemo(const edm::ParameterSet &); + ~ScPhase2PuppiHJPsiGammaDemo() override; + static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); + +private: + void beginStream(edm::StreamID) override; + void produce(edm::Event &, const edm::EventSetup &) override; + void endStream() override; + template + void runObj(const OrbitCollection &src, + const OrbitCollection &src2, + edm::Event &out, + unsigned long &nTry, + unsigned long &nPass, + const std::string &bxLabel); + + bool doStruct_; + edm::EDGetTokenT> structToken_; + edm::EDGetTokenT> struct2Token_; + + struct Cuts { + float minpt1 = 5; + float minpt2 = 5; + float minpt3 = 10; + float minptQ = 10; + float maxdeltar2 = 0.40 * 0.40; + float minmass = 100; + float maxmass = 150; + float minmass2 = 2.50; + float maxmass2 = 3.50; + float mindr2 = 0.05 * 0.05; + float maxdr2 = 0.25 * 0.25; + float maxiso = 0.25; + float mindr2tkem = 0.05 * 0.05; + float maxdr2tkem = 0.25 * 0.25; + float maxisotkem = 0.25; + } cuts; + + template + bool isolationQ( + unsigned int pidex1, unsigned int pidex2, const T *cands, unsigned int size) const; + + template + bool isolationTkEm(float pt, float eta, float phi, const T *cands, unsigned int size) const; + + std::tuple deltar(float eta1, float eta2, float phi1, float phi2) const; + + template + static float pairmass(const std::array &t, const T *cands, const std::array &massD); + + template + float tripletmass(const std::array &t, const T *cands, const U *cands2, const std::array &masses); + + unsigned long countStruct_; + unsigned long passStruct_; +}; + +ScPhase2PuppiHJPsiGammaDemo::ScPhase2PuppiHJPsiGammaDemo(const edm::ParameterSet &iConfig) + : doStruct_(iConfig.getParameter("runStruct")) { + if (doStruct_) { + structToken_ = consumes>(iConfig.getParameter("src")); + struct2Token_ = consumes>(iConfig.getParameter("src2")); + produces>("selectedBx"); + produces("hjpsigamma"); + } +} + +ScPhase2PuppiHJPsiGammaDemo::~ScPhase2PuppiHJPsiGammaDemo(){}; + +void ScPhase2PuppiHJPsiGammaDemo::beginStream(edm::StreamID) { + countStruct_ = 0; + passStruct_ = 0; +} + +void ScPhase2PuppiHJPsiGammaDemo::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) { + if (doStruct_) { + edm::Handle> src; + iEvent.getByToken(structToken_, src); + + edm::Handle> src2; + iEvent.getByToken(struct2Token_, src2); + + runObj(*src, *src2, iEvent, countStruct_, passStruct_, ""); + } +} + +void ScPhase2PuppiHJPsiGammaDemo::endStream() { + if (doStruct_) + std::cout << "HJPsiGamma Struct analysis: " << countStruct_ << " -> " << passStruct_ << std::endl; +} + +template +void ScPhase2PuppiHJPsiGammaDemo::runObj(const OrbitCollection &src, + const OrbitCollection &src2, + edm::Event &iEvent, + unsigned long &nTry, + unsigned long &nPass, + const std::string &label) { + l1ScoutingRun3::BxOffsetsFillter bxOffsetsFiller; + bxOffsetsFiller.start(); + auto ret = std::make_unique>(); + std::vector masses; + std::vector i0s, i1s, i2s; // i2s is the photon + ROOT::RVec ix; // + ROOT::RVec ig; // photons + std::array bestPair; + std::array bestTriplet; + bool bestPairFound; + float bestPairScore; + float bestPhoScore; + for (unsigned int bx = 1; bx <= OrbitCollection::NBX; ++bx) { + nTry++; + auto range = src.bxIterator(bx); + const T *cands = &range.front(); + auto size = range.size(); + + auto range2 = src2.bxIterator(bx); + const U *cands2 = &range2.front(); + auto size2 = range2.size(); + + ix.clear(); + int highcut = 0; + for (unsigned int i = 0; i < size; ++i) { // make list of puppi with a certain threshold + if (cands[i].pt() >= cuts.minpt1) { + ix.push_back(i); + if (cands[i].pt() >= cuts.minpt2) + highcut++; + } + } + unsigned int ndaus = ix.size(); + if (highcut < 1 || ndaus < 2) + continue; + + ig.clear(); + for (unsigned int i = 0; i < size2; ++i) { // make list of all photons + if (cands2[i].pt() >= cuts.minpt3) { + ig.push_back(i); + } + } + unsigned int ngammas = ig.size(); + if (ngammas < 1) + continue; + + // Q candidate from closest OS pair with mass compatible with mQ + bestPairFound = false; + bestPairScore = 999; + for (unsigned int i1 = 0; i1 < ndaus; ++i1) { + if (cands[ix[i1]].pt() < cuts.minpt2) + continue; // high pt cut + for (unsigned int i2 = 0; i2 < ndaus; ++i2) { + if (i2 == i1 || cands[ix[i2]].pt() < cuts.minpt1) + continue; + + if (!(cands[ix[i1]].charge()*cands[ix[i2]].charge() < 0)) + continue; + + auto mass2 = pairmass({{ix[i1], ix[i2]}}, cands, {{0.1057, 0.1057}}); + if (mass2 >= cuts.minmass2 and mass2 <= cuts.maxmass2) + continue; + + auto [drcond, drQ] = deltar(cands[ix[i1]].eta(), cands[ix[i2]].eta(), cands[ix[i1]].phi(), cands[ix[i2]].phi()); + if (!drcond) + continue; //angular sep of top 2 tracks + + std::array pair{{ix[i1], ix[i2]}}; // pair of indices + if (drQ < bestPairScore) { + std::copy_n(pair.begin(), 2, bestPair.begin()); + bestPairScore = drQ; + if (bestPairScore*bestPairScore < cuts.maxdeltar2) + bestPairFound = true; + } + } + } + if (!bestPairFound) + continue; + + // Q pt + auto ptQ = (cands[bestPair[0]].p4() + cands[bestPair[1]].p4()).pt(); + if (ptQ < cuts.minptQ) + continue; + + // Q isolation + if (!isolationQ(bestPair[0], bestPair[1], cands, size)) + continue; + + // photon + bestPhoScore = 0; + for (unsigned int i3 = 0; i3 < ngammas; ++i3) { + if (cands2[ig[i3]].pt() < cuts.minpt3) + continue; // photon pt cut + + std::array tr{{bestPair[0], bestPair[1], ig[i3]}}; // triplet of indices + + if (cands2[ig[i3]].pt() > bestPhoScore) { + std::copy_n(tr.begin(), 3, bestTriplet.begin()); + bestPhoScore = cands2[ig[i3]].pt(); + } + } + if (bestPhoScore < 0) + continue; + + // photon isolation + bool isop = isolationTkEm( + cands2[bestTriplet[2]].pt(), cands2[bestTriplet[2]].eta(), cands2[bestTriplet[2]].phi(), cands, size); + if (!isop) + continue; + + // H mass + auto mass = tripletmass(bestTriplet, cands, cands2, {{0.1057, 0.1057, 0.}}); + if (!(mass >= cuts.minmass and mass <= cuts.maxmass)) + continue; + + ret->emplace_back(bx); + nPass++; + masses.push_back(mass); + i0s.push_back(bestTriplet[0]); + i1s.push_back(bestTriplet[1]); + i2s.push_back(bestTriplet[2]); + bxOffsetsFiller.addBx(bx, 1); + } // loop on BXs + + iEvent.put(std::move(ret), "selectedBx" + label); + // now we make the table + auto bxOffsets = bxOffsetsFiller.done(); + auto tab = std::make_unique(bxOffsets, "HJPsiGamma" + label, true); + tab->addColumn("mass", masses, "2 muons plus photon invariant mass"); + tab->addColumn("i0", i0s, "leading muon"); + tab->addColumn("i1", i1s, "subleading muon"); + tab->addColumn("i2", i2s, "photon"); + iEvent.put(std::move(tab), "hjpsigamma" + label); +} + +//TEST functions +template +bool ScPhase2PuppiHJPsiGammaDemo::isolationQ( + unsigned int pidex1, unsigned int pidex2, const T *cands, unsigned int size) const { + bool passed = false; + float psum = 0; + float eta = cands[pidex1].eta(); //center cone around leading track + float phi = cands[pidex1].phi(); + for (unsigned int j = 0u; j < size; ++j) { //loop over other particles + if (pidex1 == j or pidex2 == j) + continue; + float deta = eta - cands[j].eta(), dphi = ROOT::VecOps::DeltaPhi(phi, cands[j].phi()); + float dr2 = deta * deta + dphi * dphi; + if (dr2 >= cuts.mindr2 && dr2 <= cuts.maxdr2) + psum += cands[j].pt(); + } + if (psum <= cuts.maxiso * (cands[pidex1].pt() + cands[pidex2].pt())) + passed = true; + return passed; +} + +template +bool ScPhase2PuppiHJPsiGammaDemo::isolationTkEm(float pt, float eta, float phi, const T *cands, unsigned int size) const { + bool passed = false; + float psum = 0; + for (unsigned int j = 0u; j < size; ++j) { //loop over other particles + float deta = eta - cands[j].eta(), dphi = ROOT::VecOps::DeltaPhi(phi, cands[j].phi()); + float dr2 = deta * deta + dphi * dphi; + if (dr2 >= cuts.mindr2tkem && dr2 <= cuts.maxdr2tkem) + psum += cands[j].pt(); + } + if (psum <= cuts.maxisotkem * pt) + passed = true; + return passed; +} + +std::tuple ScPhase2PuppiHJPsiGammaDemo::deltar(float eta1, float eta2, float phi1, float phi2) const { + bool passed = true; + float deta = eta1 - eta2; + float dphi = ROOT::VecOps::DeltaPhi(phi1, phi2); + float dr2 = deta * deta + dphi * dphi; + if (dr2 > cuts.maxdeltar2) { + passed = false; + return std::tuple(passed, dr2); + } + return std::tuple(passed, dr2); +} + +template +float ScPhase2PuppiHJPsiGammaDemo::pairmass(const std::array &t, + const T *cands, + const std::array &massD) { + ROOT::Math::PtEtaPhiMVector p1(cands[t[0]].pt(), cands[t[0]].eta(), cands[t[0]].phi(), massD[0]); + ROOT::Math::PtEtaPhiMVector p2(cands[t[1]].pt(), cands[t[1]].eta(), cands[t[1]].phi(), massD[1]); + float mass = (p1 + p2).M(); + return mass; +} + +template +float ScPhase2PuppiHJPsiGammaDemo::tripletmass(const std::array &t, + const T *cands, + const U *cands2, + const std::array &masses) { + ROOT::Math::PtEtaPhiMVector p1(cands[t[0]].pt(), cands[t[0]].eta(), cands[t[0]].phi(), masses[0]); + ROOT::Math::PtEtaPhiMVector p2(cands[t[1]].pt(), cands[t[1]].eta(), cands[t[1]].phi(), masses[1]); + ROOT::Math::PtEtaPhiMVector p3(cands2[t[2]].pt(), cands2[t[2]].eta(), cands2[t[2]].phi(), masses[2]); + float mass = (p1 + p2 + p3).M(); + return mass; +} + +void ScPhase2PuppiHJPsiGammaDemo::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + +DEFINE_FWK_MODULE(ScPhase2PuppiHJPsiGammaDemo); diff --git a/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiHPhiGammaDemo.cc b/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiHPhiGammaDemo.cc new file mode 100644 index 0000000000000..24f7eb65384af --- /dev/null +++ b/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiHPhiGammaDemo.cc @@ -0,0 +1,335 @@ +#include +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/StreamID.h" + +#include "DataFormats/L1Scouting/interface/OrbitCollection.h" +#include "DataFormats/L1Scouting/interface/OrbitFlatTable.h" +#include "DataFormats/L1TParticleFlow/interface/L1ScoutingPuppi.h" +#include "DataFormats/L1TParticleFlow/interface/L1ScoutingTkEm.h" +#include "L1TriggerScouting/Utilities/interface/BxOffsetsFiller.h" + +#include +#include +#include +#include +#include +#include +#include + +class ScPhase2PuppiHPhiGammaDemo : public edm::stream::EDProducer<> { +public: + explicit ScPhase2PuppiHPhiGammaDemo(const edm::ParameterSet &); + ~ScPhase2PuppiHPhiGammaDemo() override; + static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); + +private: + void beginStream(edm::StreamID) override; + void produce(edm::Event &, const edm::EventSetup &) override; + void endStream() override; + template + void runObj(const OrbitCollection &src, + const OrbitCollection &src2, + edm::Event &out, + unsigned long &nTry, + unsigned long &nPass, + const std::string &bxLabel); + + bool doStruct_; + edm::EDGetTokenT> structToken_; + edm::EDGetTokenT> struct2Token_; + + struct Cuts { + float minpt1 = 5; + float minpt2 = 5; + float minpt3 = 10; + float minptQ = 10; + float maxdeltar2 = 0.40 * 0.40; + float minmass = 100; + float maxmass = 150; + float minmass2 = 0.95; + float maxmass2 = 1.25; + float mindr2 = 0.05 * 0.05; + float maxdr2 = 0.25 * 0.25; + float maxiso = 0.25; + float mindr2tkem = 0.05 * 0.05; + float maxdr2tkem = 0.25 * 0.25; + float maxisotkem = 0.25; + } cuts; + + template + bool isolationQ( + unsigned int pidex1, unsigned int pidex2, const T *cands, unsigned int size) const; + + template + bool isolationTkEm(float pt, float eta, float phi, const T *cands, unsigned int size) const; + + std::tuple deltar(float eta1, float eta2, float phi1, float phi2) const; + + template + static float pairmass(const std::array &t, const T *cands, const std::array &massD); + + template + float tripletmass(const std::array &t, const T *cands, const U *cands2, const std::array &masses); + + unsigned long countStruct_; + unsigned long passStruct_; +}; + +ScPhase2PuppiHPhiGammaDemo::ScPhase2PuppiHPhiGammaDemo(const edm::ParameterSet &iConfig) + : doStruct_(iConfig.getParameter("runStruct")) { + if (doStruct_) { + structToken_ = consumes>(iConfig.getParameter("src")); + struct2Token_ = consumes>(iConfig.getParameter("src2")); + produces>("selectedBx"); + produces("hphigamma"); + } +} + +ScPhase2PuppiHPhiGammaDemo::~ScPhase2PuppiHPhiGammaDemo(){}; + +void ScPhase2PuppiHPhiGammaDemo::beginStream(edm::StreamID) { + countStruct_ = 0; + passStruct_ = 0; +} + +void ScPhase2PuppiHPhiGammaDemo::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) { + if (doStruct_) { + edm::Handle> src; + iEvent.getByToken(structToken_, src); + + edm::Handle> src2; + iEvent.getByToken(struct2Token_, src2); + + runObj(*src, *src2, iEvent, countStruct_, passStruct_, ""); + } +} + +void ScPhase2PuppiHPhiGammaDemo::endStream() { + if (doStruct_) + std::cout << "HPhiGamma Struct analysis: " << countStruct_ << " -> " << passStruct_ << std::endl; +} + +template +void ScPhase2PuppiHPhiGammaDemo::runObj(const OrbitCollection &src, + const OrbitCollection &src2, + edm::Event &iEvent, + unsigned long &nTry, + unsigned long &nPass, + const std::string &label) { + l1ScoutingRun3::BxOffsetsFillter bxOffsetsFiller; + bxOffsetsFiller.start(); + auto ret = std::make_unique>(); + std::vector masses; + std::vector i0s, i1s, i2s; // i2s is the photon + ROOT::RVec ix; // + ROOT::RVec ig; // photons + std::array bestPair; + std::array bestTriplet; + bool bestPairFound; + float bestPairScore; + float bestPhoScore; + for (unsigned int bx = 1; bx <= OrbitCollection::NBX; ++bx) { + nTry++; + auto range = src.bxIterator(bx); + const T *cands = &range.front(); + auto size = range.size(); + + auto range2 = src2.bxIterator(bx); + const U *cands2 = &range2.front(); + auto size2 = range2.size(); + + ix.clear(); + int highcut = 0; + for (unsigned int i = 0; i < size; ++i) { //make list of all hadrons + if ((std::abs(cands[i].pdgId()) == 211 or std::abs(cands[i].pdgId()) == 11)) { + if (cands[i].pt() >= cuts.minpt1) { + ix.push_back(i); + if (cands[i].pt() >= cuts.minpt2) + highcut++; + } + } + } + unsigned int ndaus = ix.size(); + if (highcut < 1 || ndaus < 2) + continue; + + ig.clear(); + for (unsigned int i = 0; i < size2; ++i) { // make list of all photons + if (cands2[i].pt() >= cuts.minpt3) { + ig.push_back(i); + } + } + unsigned int ngammas = ig.size(); + if (ngammas < 1) + continue; + + // Q candidate from closest OS pair with mass compatible with mQ + bestPairFound = false; + bestPairScore = 999; + for (unsigned int i1 = 0; i1 < ndaus; ++i1) { + if (cands[ix[i1]].pt() < cuts.minpt2) + continue; // high pt cut + for (unsigned int i2 = 0; i2 < ndaus; ++i2) { + if (i2 == i1 || cands[ix[i2]].pt() < cuts.minpt1) + continue; + + if (!(cands[ix[i1]].charge()*cands[ix[i2]].charge() < 0)) + continue; + + auto mass2 = pairmass({{ix[i1], ix[i2]}}, cands, {{0.4937, 0.4937}}); + if (mass2 >= cuts.minmass2 and mass2 <= cuts.maxmass2) + continue; + + auto [drcond, drQ] = deltar(cands[ix[i1]].eta(), cands[ix[i2]].eta(), cands[ix[i1]].phi(), cands[ix[i2]].phi()); + if (!drcond) + continue; //angular sep of top 2 tracks + + std::array pair{{ix[i1], ix[i2]}}; // pair of indices + if (drQ < bestPairScore) { + std::copy_n(pair.begin(), 2, bestPair.begin()); + bestPairScore = drQ; + if (bestPairScore*bestPairScore < cuts.maxdeltar2) + bestPairFound = true; + } + } + } + if (!bestPairFound) + continue; + + // Q pt + auto ptQ = (cands[bestPair[0]].p4() + cands[bestPair[1]].p4()).pt(); + if (ptQ < cuts.minptQ) + continue; + + // Q isolation + if (!isolationQ(bestPair[0], bestPair[1], cands, size)) + continue; + + // photon + bestPhoScore = 0; + for (unsigned int i3 = 0; i3 < ngammas; ++i3) { + if (cands2[ig[i3]].pt() < cuts.minpt3) + continue; // photon pt cut + + std::array tr{{bestPair[0], bestPair[1], ig[i3]}}; // triplet of indices + + if (cands2[ig[i3]].pt() > bestPhoScore) { + std::copy_n(tr.begin(), 3, bestTriplet.begin()); + bestPhoScore = cands2[ig[i3]].pt(); + } + } + if (bestPhoScore < 0) + continue; + + // photon isolation + bool isop = isolationTkEm( + cands2[bestTriplet[2]].pt(), cands2[bestTriplet[2]].eta(), cands2[bestTriplet[2]].phi(), cands, size); + if (!isop) + continue; + + // H mass + auto mass = tripletmass(bestTriplet, cands, cands2, {{0.4937, 0.4937, 0.}}); + if (!(mass >= cuts.minmass and mass <= cuts.maxmass)) + continue; + + ret->emplace_back(bx); + nPass++; + masses.push_back(mass); + i0s.push_back(bestTriplet[0]); + i1s.push_back(bestTriplet[1]); + i2s.push_back(bestTriplet[2]); + bxOffsetsFiller.addBx(bx, 1); + } // loop on BXs + + iEvent.put(std::move(ret), "selectedBx" + label); + // now we make the table + auto bxOffsets = bxOffsetsFiller.done(); + auto tab = std::make_unique(bxOffsets, "HPhiGamma" + label, true); + tab->addColumn("mass", masses, "2 kaons plus photon invariant mass"); + tab->addColumn("i0", i0s, "leading kaon"); + tab->addColumn("i1", i1s, "subleading kaon"); + tab->addColumn("i2", i2s, "photon"); + iEvent.put(std::move(tab), "hphigamma" + label); +} + +//TEST functions +template +bool ScPhase2PuppiHPhiGammaDemo::isolationQ( + unsigned int pidex1, unsigned int pidex2, const T *cands, unsigned int size) const { + bool passed = false; + float psum = 0; + float eta = cands[pidex1].eta(); //center cone around leading track + float phi = cands[pidex1].phi(); + for (unsigned int j = 0u; j < size; ++j) { //loop over other particles + if (pidex1 == j or pidex2 == j) + continue; + float deta = eta - cands[j].eta(), dphi = ROOT::VecOps::DeltaPhi(phi, cands[j].phi()); + float dr2 = deta * deta + dphi * dphi; + if (dr2 >= cuts.mindr2 && dr2 <= cuts.maxdr2) + psum += cands[j].pt(); + } + if (psum <= cuts.maxiso * (cands[pidex1].pt() + cands[pidex2].pt())) + passed = true; + return passed; +} + +template +bool ScPhase2PuppiHPhiGammaDemo::isolationTkEm(float pt, float eta, float phi, const T *cands, unsigned int size) const { + bool passed = false; + float psum = 0; + for (unsigned int j = 0u; j < size; ++j) { //loop over other particles + float deta = eta - cands[j].eta(), dphi = ROOT::VecOps::DeltaPhi(phi, cands[j].phi()); + float dr2 = deta * deta + dphi * dphi; + if (dr2 >= cuts.mindr2tkem && dr2 <= cuts.maxdr2tkem) + psum += cands[j].pt(); + } + if (psum <= cuts.maxisotkem * pt) + passed = true; + return passed; +} + +std::tuple ScPhase2PuppiHPhiGammaDemo::deltar(float eta1, float eta2, float phi1, float phi2) const { + bool passed = true; + float deta = eta1 - eta2; + float dphi = ROOT::VecOps::DeltaPhi(phi1, phi2); + float dr2 = deta * deta + dphi * dphi; + if (dr2 > cuts.maxdeltar2) { + passed = false; + return std::tuple(passed, dr2); + } + return std::tuple(passed, dr2); +} + +template +float ScPhase2PuppiHPhiGammaDemo::pairmass(const std::array &t, + const T *cands, + const std::array &massD) { + ROOT::Math::PtEtaPhiMVector p1(cands[t[0]].pt(), cands[t[0]].eta(), cands[t[0]].phi(), massD[0]); + ROOT::Math::PtEtaPhiMVector p2(cands[t[1]].pt(), cands[t[1]].eta(), cands[t[1]].phi(), massD[1]); + float mass = (p1 + p2).M(); + return mass; +} + +template +float ScPhase2PuppiHPhiGammaDemo::tripletmass(const std::array &t, + const T *cands, + const U *cands2, + const std::array &masses) { + ROOT::Math::PtEtaPhiMVector p1(cands[t[0]].pt(), cands[t[0]].eta(), cands[t[0]].phi(), masses[0]); + ROOT::Math::PtEtaPhiMVector p2(cands[t[1]].pt(), cands[t[1]].eta(), cands[t[1]].phi(), masses[1]); + ROOT::Math::PtEtaPhiMVector p3(cands2[t[2]].pt(), cands2[t[2]].eta(), cands2[t[2]].phi(), masses[2]); + float mass = (p1 + p2 + p3).M(); + return mass; +} + +void ScPhase2PuppiHPhiGammaDemo::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + +DEFINE_FWK_MODULE(ScPhase2PuppiHPhiGammaDemo); diff --git a/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiHPhiJPsiDemo.cc b/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiHPhiJPsiDemo.cc new file mode 100644 index 0000000000000..bc7fb3e95ebba --- /dev/null +++ b/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiHPhiJPsiDemo.cc @@ -0,0 +1,300 @@ +#include +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/StreamID.h" + +#include "DataFormats/L1Scouting/interface/OrbitCollection.h" +#include "DataFormats/L1Scouting/interface/OrbitFlatTable.h" +#include "DataFormats/L1TParticleFlow/interface/L1ScoutingPuppi.h" +#include "DataFormats/L1TParticleFlow/interface/L1ScoutingTkEm.h" +#include "L1TriggerScouting/Utilities/interface/BxOffsetsFiller.h" + +#include +#include +#include +#include +#include +#include +#include + +class ScPhase2PuppiHPhiJPsiDemo : public edm::stream::EDProducer<> { +public: + explicit ScPhase2PuppiHPhiJPsiDemo(const edm::ParameterSet &); + ~ScPhase2PuppiHPhiJPsiDemo() override; + static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); + +private: + void beginStream(edm::StreamID) override; + void produce(edm::Event &, const edm::EventSetup &) override; + void endStream() override; + template + void runObj(const OrbitCollection &src, + edm::Event &out, + unsigned long &nTry, + unsigned long &nPass, + const std::string &bxLabel); + + bool doStruct_; + edm::EDGetTokenT> structToken_; + + struct Cuts { + float minptD = 5; + float minptQ = 10; + float maxdeltarD2 = 0.40 * 0.40; + float minmassH = 100; + float maxmassH = 150; + float minmassQ1 = 0.95; + float maxmassQ1 = 1.25; + float minmassQ2 = 2.50; + float maxmassQ2 = 3.50; + float mindr2 = 0.05 * 0.05; + float maxdr2 = 0.25 * 0.25; + float maxiso = 0.25; + } cuts; + + template + bool isolationQ( + unsigned int pidex1, unsigned int pidex2, const T *cands, unsigned int size) const; + + std::tuple deltar(float eta1, float eta2, float phi1, float phi2) const; + + template + static float pairmass(const std::array &t, const T *cands, const std::array &massD); + + template + static float quadrupletmass(const std::array &t, const T *cands, const std::array &massD); + + unsigned long countStruct_; + unsigned long passStruct_; +}; + +ScPhase2PuppiHPhiJPsiDemo::ScPhase2PuppiHPhiJPsiDemo(const edm::ParameterSet &iConfig) + : doStruct_(iConfig.getParameter("runStruct")) { + if (doStruct_) { + structToken_ = consumes>(iConfig.getParameter("src")); + produces>("selectedBx"); + produces("hphijpsi"); + } +} + +ScPhase2PuppiHPhiJPsiDemo::~ScPhase2PuppiHPhiJPsiDemo(){}; + +void ScPhase2PuppiHPhiJPsiDemo::beginStream(edm::StreamID) { + countStruct_ = 0; + passStruct_ = 0; +} + +void ScPhase2PuppiHPhiJPsiDemo::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) { + if (doStruct_) { + edm::Handle> src; + iEvent.getByToken(structToken_, src); + + runObj(*src, iEvent, countStruct_, passStruct_, ""); + } +} + +void ScPhase2PuppiHPhiJPsiDemo::endStream() { + if (doStruct_) + std::cout << "HPhiJPsi Struct analysis: " << countStruct_ << " -> " << passStruct_ << std::endl; +} + +template +void ScPhase2PuppiHPhiJPsiDemo::runObj(const OrbitCollection &src, + edm::Event &iEvent, + unsigned long &nTry, + unsigned long &nPass, + const std::string &label) { + l1ScoutingRun3::BxOffsetsFillter bxOffsetsFiller; + bxOffsetsFiller.start(); + auto ret = std::make_unique>(); + std::vector masses; + std::vector i0s, i1s, i2s, i3s; + ROOT::RVec ix; // + std::array bestPair1, bestPair2; + bool bestPair1Found, bestPair2Found; + float bestPair1Score, bestPair2Score; + for (unsigned int bx = 1; bx <= OrbitCollection::NBX; ++bx) { + nTry++; + auto range = src.bxIterator(bx); + const T *cands = &range.front(); + auto size = range.size(); + + ix.clear(); + for (unsigned int i = 0; i < size; ++i) { //make list of all hadrons + if (cands[i].pt() >= cuts.minptD) + ix.push_back(i); + } + unsigned int ndaus = ix.size(); + if (ndaus < 4) + continue; + + // Q1 candidate from closest OS pair with mass compatible with mQ + bestPair1Found = false; + bestPair1Score = 999; + for (unsigned int i1 = 0; i1 < ndaus; ++i1) { + if (cands[ix[i1]].pt() < cuts.minptD) + continue; // D1 pt cut + for (unsigned int i2 = 0; i2 < ndaus; ++i2) { + if (i2 == i1 || cands[ix[i2]].pt() < cuts.minptD) + continue; // D2 pt cut + + if (!(cands[ix[i1]].charge()*cands[ix[i2]].charge() < 0)) + continue; + + auto mass2 = pairmass({{ix[i1], ix[i2]}}, cands, {{0.4937, 0.4937}}); + if (mass2 >= cuts.minmassQ1 and mass2 <= cuts.maxmassQ1) + continue; + + auto [drcond, drQ] = deltar(cands[ix[i1]].eta(), cands[ix[i2]].eta(), cands[ix[i1]].phi(), cands[ix[i2]].phi()); + if (!drcond) + continue; // angular sep of top 2 tracks + + std::array pair{{ix[i1], ix[i2]}}; // pair of indices + if (drQ < bestPair1Score) { + std::copy_n(pair.begin(), 2, bestPair1.begin()); + bestPair1Score = drQ; + if (bestPair1Score*bestPair1Score < cuts.maxdeltarD2) + bestPair1Found = true; + } + } + } + if (!bestPair1Found) + continue; // pair was found + auto ptQ = (cands[bestPair1[0]].p4() + cands[bestPair1[1]].p4()).pt(); + if (ptQ < cuts.minptQ) + continue; // Q pt + if (!isolationQ(bestPair1[0], bestPair1[1], cands, size)) + continue; // Q isolation + + // Q2 candidate from closest OS pair with mass compatible with mQ + bestPair2Found = false; + bestPair2Score = 999; + for (unsigned int i3 = 0; i3 < ndaus; ++i3) { + if (cands[ix[i3]].pt() < cuts.minptD) + continue; // D1 pt cut + if (ix[i3]==bestPair1[0] or ix[i3]==bestPair1[1]) + continue; // don't reuse candidates from previous pair + for (unsigned int i4 = 0; i4 < ndaus; ++i4) { + if (i4 == i3 || cands[ix[i4]].pt() < cuts.minptD) + continue; // D2 pt cut + if (ix[i4]==bestPair1[0] or ix[i4]==bestPair1[1]) + continue; // don't reuse candidates from previous pair + if (!(cands[ix[i3]].charge()*cands[ix[i4]].charge() < 0)) + continue; // OS pair + auto mass2 = pairmass({{ix[i3], ix[i4]}}, cands, {{0.1057, 0.1057}}); // (cands[ix[i3]].p4() + cands[ix[i4]].p4()).mass(); + if (mass2 >= cuts.minmassQ2 and mass2 <= cuts.maxmassQ2) + continue; // Q mass + auto [drcond, drQ] = deltar(cands[ix[i3]].eta(), cands[ix[i4]].eta(), cands[ix[i3]].phi(), cands[ix[i4]].phi()); + if (!drcond) + continue; // angular sep of top 2 tracks + + std::array pair{{ix[i3], ix[i4]}}; // pair of indices + if (drQ < bestPair2Score) { + std::copy_n(pair.begin(), 2, bestPair2.begin()); + bestPair2Score = drQ; + if (bestPair2Score*bestPair2Score < cuts.maxdeltarD2) + bestPair2Found = true; + } + } + } + if (!bestPair2Found) + continue; // pair was found + ptQ = (cands[bestPair2[0]].p4() + cands[bestPair2[1]].p4()).pt(); + if (ptQ < cuts.minptQ) + continue; // Q pt + if (!isolationQ(bestPair2[0], bestPair2[1], cands, size)) + continue; // Q isolation + + std::array bestQuadruplet{{bestPair1[0], bestPair1[1], bestPair2[0], bestPair2[1]}}; + // H mass + auto mass = quadrupletmass(bestQuadruplet, cands, {{0.4937, 0.4937, 0.1057, 0.1057}}); + if (!(mass >= cuts.minmassH and mass <= cuts.maxmassH)) + continue; + + ret->emplace_back(bx); + nPass++; + masses.push_back(mass); + i0s.push_back(bestQuadruplet[0]); + i1s.push_back(bestQuadruplet[1]); + i2s.push_back(bestQuadruplet[2]); + i3s.push_back(bestQuadruplet[3]); + bxOffsetsFiller.addBx(bx, 1); + } // loop on BXs + + iEvent.put(std::move(ret), "selectedBx" + label); + // now we make the table + auto bxOffsets = bxOffsetsFiller.done(); + auto tab = std::make_unique(bxOffsets, "HPhiJPsi" + label, true); + tab->addColumn("mass", masses, "4 kaons invariant mass"); + tab->addColumn("i0", i0s, "1st kaon (phi)"); + tab->addColumn("i1", i1s, "2nd kaon (phi)"); + tab->addColumn("i2", i2s, "1st muon (jpsi)"); + tab->addColumn("i3", i3s, "2nd muon (jpsi)"); + iEvent.put(std::move(tab), "hphijpsi" + label); +} + +//TEST functions +template +bool ScPhase2PuppiHPhiJPsiDemo::isolationQ( + unsigned int pidex1, unsigned int pidex2, const T *cands, unsigned int size) const { + bool passed = false; + float psum = 0; + float eta = cands[pidex1].eta(); //center cone around leading track + float phi = cands[pidex1].phi(); + for (unsigned int j = 0u; j < size; ++j) { //loop over other particles + if (pidex1 == j or pidex2 == j) + continue; + float deta = eta - cands[j].eta(), dphi = ROOT::VecOps::DeltaPhi(phi, cands[j].phi()); + float dr2 = deta * deta + dphi * dphi; + if (dr2 >= cuts.mindr2 && dr2 <= cuts.maxdr2) + psum += cands[j].pt(); + } + if (psum <= cuts.maxiso * (cands[pidex1].pt() + cands[pidex2].pt())) + passed = true; + return passed; +} + +std::tuple ScPhase2PuppiHPhiJPsiDemo::deltar(float eta1, float eta2, float phi1, float phi2) const { + bool passed = true; + float deta = eta1 - eta2; + float dphi = ROOT::VecOps::DeltaPhi(phi1, phi2); + float dr2 = deta * deta + dphi * dphi; + if (dr2 > cuts.maxdeltarD2) { + passed = false; + return std::tuple(passed, dr2); + } + return std::tuple(passed, dr2); +} + +template +float ScPhase2PuppiHPhiJPsiDemo::pairmass(const std::array &t, + const T *cands, + const std::array &massD) { + ROOT::Math::PtEtaPhiMVector p1(cands[t[0]].pt(), cands[t[0]].eta(), cands[t[0]].phi(), massD[0]); + ROOT::Math::PtEtaPhiMVector p2(cands[t[1]].pt(), cands[t[1]].eta(), cands[t[1]].phi(), massD[1]); + float mass = (p1 + p2).M(); + return mass; +} + +template +float ScPhase2PuppiHPhiJPsiDemo::quadrupletmass(const std::array &t, + const T *cands, + const std::array &massD) { + ROOT::Math::PtEtaPhiMVector p1(cands[t[0]].pt(), cands[t[0]].eta(), cands[t[0]].phi(), massD[0]); + ROOT::Math::PtEtaPhiMVector p2(cands[t[1]].pt(), cands[t[1]].eta(), cands[t[1]].phi(), massD[1]); + ROOT::Math::PtEtaPhiMVector p3(cands[t[2]].pt(), cands[t[2]].eta(), cands[t[2]].phi(), massD[2]); + ROOT::Math::PtEtaPhiMVector p4(cands[t[3]].pt(), cands[t[3]].eta(), cands[t[3]].phi(), massD[3]); + float mass = (p1 + p2 + p3 + p4).M(); + return mass; +} + +void ScPhase2PuppiHPhiJPsiDemo::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + +DEFINE_FWK_MODULE(ScPhase2PuppiHPhiJPsiDemo); diff --git a/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiHRhoGammaDemo.cc b/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiHRhoGammaDemo.cc new file mode 100644 index 0000000000000..7fb9e6efedcf1 --- /dev/null +++ b/L1TriggerScouting/Phase2/plugins/ScPhase2PuppiHRhoGammaDemo.cc @@ -0,0 +1,335 @@ +#include +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/StreamID.h" + +#include "DataFormats/L1Scouting/interface/OrbitCollection.h" +#include "DataFormats/L1Scouting/interface/OrbitFlatTable.h" +#include "DataFormats/L1TParticleFlow/interface/L1ScoutingPuppi.h" +#include "DataFormats/L1TParticleFlow/interface/L1ScoutingTkEm.h" +#include "L1TriggerScouting/Utilities/interface/BxOffsetsFiller.h" + +#include +#include +#include +#include +#include +#include +#include + +class ScPhase2PuppiHRhoGammaDemo : public edm::stream::EDProducer<> { +public: + explicit ScPhase2PuppiHRhoGammaDemo(const edm::ParameterSet &); + ~ScPhase2PuppiHRhoGammaDemo() override; + static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); + +private: + void beginStream(edm::StreamID) override; + void produce(edm::Event &, const edm::EventSetup &) override; + void endStream() override; + template + void runObj(const OrbitCollection &src, + const OrbitCollection &src2, + edm::Event &out, + unsigned long &nTry, + unsigned long &nPass, + const std::string &bxLabel); + + bool doStruct_; + edm::EDGetTokenT> structToken_; + edm::EDGetTokenT> struct2Token_; + + struct Cuts { + float minpt1 = 5; + float minpt2 = 5; + float minpt3 = 10; + float minptQ = 10; + float maxdeltar2 = 0.40 * 0.40; + float minmass = 100; + float maxmass = 150; + float minmass2 = 0.40; + float maxmass2 = 1.30; + float mindr2 = 0.05 * 0.05; + float maxdr2 = 0.25 * 0.25; + float maxiso = 0.25; + float mindr2tkem = 0.05 * 0.05; + float maxdr2tkem = 0.25 * 0.25; + float maxisotkem = 0.25; + } cuts; + + template + bool isolationQ( + unsigned int pidex1, unsigned int pidex2, const T *cands, unsigned int size) const; + + template + bool isolationTkEm(float pt, float eta, float phi, const T *cands, unsigned int size) const; + + std::tuple deltar(float eta1, float eta2, float phi1, float phi2) const; + + template + static float pairmass(const std::array &t, const T *cands, const std::array &massD); + + template + float tripletmass(const std::array &t, const T *cands, const U *cands2, const std::array &masses); + + unsigned long countStruct_; + unsigned long passStruct_; +}; + +ScPhase2PuppiHRhoGammaDemo::ScPhase2PuppiHRhoGammaDemo(const edm::ParameterSet &iConfig) + : doStruct_(iConfig.getParameter("runStruct")) { + if (doStruct_) { + structToken_ = consumes>(iConfig.getParameter("src")); + struct2Token_ = consumes>(iConfig.getParameter("src2")); + produces>("selectedBx"); + produces("hrhogamma"); + } +} + +ScPhase2PuppiHRhoGammaDemo::~ScPhase2PuppiHRhoGammaDemo(){}; + +void ScPhase2PuppiHRhoGammaDemo::beginStream(edm::StreamID) { + countStruct_ = 0; + passStruct_ = 0; +} + +void ScPhase2PuppiHRhoGammaDemo::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) { + if (doStruct_) { + edm::Handle> src; + iEvent.getByToken(structToken_, src); + + edm::Handle> src2; + iEvent.getByToken(struct2Token_, src2); + + runObj(*src, *src2, iEvent, countStruct_, passStruct_, ""); + } +} + +void ScPhase2PuppiHRhoGammaDemo::endStream() { + if (doStruct_) + std::cout << "HRhoGamma Struct analysis: " << countStruct_ << " -> " << passStruct_ << std::endl; +} + +template +void ScPhase2PuppiHRhoGammaDemo::runObj(const OrbitCollection &src, + const OrbitCollection &src2, + edm::Event &iEvent, + unsigned long &nTry, + unsigned long &nPass, + const std::string &label) { + l1ScoutingRun3::BxOffsetsFillter bxOffsetsFiller; + bxOffsetsFiller.start(); + auto ret = std::make_unique>(); + std::vector masses; + std::vector i0s, i1s, i2s; // i2s is the photon + ROOT::RVec ix; // + ROOT::RVec ig; // photons + std::array bestPair; + std::array bestTriplet; + bool bestPairFound; + float bestPairScore; + float bestPhoScore; + for (unsigned int bx = 1; bx <= OrbitCollection::NBX; ++bx) { + nTry++; + auto range = src.bxIterator(bx); + const T *cands = &range.front(); + auto size = range.size(); + + auto range2 = src2.bxIterator(bx); + const U *cands2 = &range2.front(); + auto size2 = range2.size(); + + ix.clear(); + int highcut = 0; + for (unsigned int i = 0; i < size; ++i) { //make list of all hadrons + if ((std::abs(cands[i].pdgId()) == 211 or std::abs(cands[i].pdgId()) == 11)) { + if (cands[i].pt() >= cuts.minpt1) { + ix.push_back(i); + if (cands[i].pt() >= cuts.minpt2) + highcut++; + } + } + } + unsigned int ndaus = ix.size(); + if (highcut < 1 || ndaus < 2) + continue; + + ig.clear(); + for (unsigned int i = 0; i < size2; ++i) { // make list of all photons + if (cands2[i].pt() >= cuts.minpt3) { + ig.push_back(i); + } + } + unsigned int ngammas = ig.size(); + if (ngammas < 1) + continue; + + // Q candidate from closest OS pair with mass compatible with mQ + bestPairFound = false; + bestPairScore = 999; + for (unsigned int i1 = 0; i1 < ndaus; ++i1) { + if (cands[ix[i1]].pt() < cuts.minpt2) + continue; // high pt cut + for (unsigned int i2 = 0; i2 < ndaus; ++i2) { + if (i2 == i1 || cands[ix[i2]].pt() < cuts.minpt1) + continue; + + if (!(cands[ix[i1]].charge()*cands[ix[i2]].charge() < 0)) + continue; + + auto mass2 = pairmass({{ix[i1], ix[i2]}}, cands, {{0.1396, 0.1396}}); + if (mass2 >= cuts.minmass2 and mass2 <= cuts.maxmass2) + continue; + + auto [drcond, drQ] = deltar(cands[ix[i1]].eta(), cands[ix[i2]].eta(), cands[ix[i1]].phi(), cands[ix[i2]].phi()); + if (!drcond) + continue; //angular sep of top 2 tracks + + std::array pair{{ix[i1], ix[i2]}}; // pair of indices + if (drQ < bestPairScore) { + std::copy_n(pair.begin(), 2, bestPair.begin()); + bestPairScore = drQ; + if (bestPairScore*bestPairScore < cuts.maxdeltar2) + bestPairFound = true; + } + } + } + if (!bestPairFound) + continue; + + // Q pt + auto ptQ = (cands[bestPair[0]].p4() + cands[bestPair[1]].p4()).pt(); + if (ptQ < cuts.minptQ) + continue; + + // Q isolation + if (!isolationQ(bestPair[0], bestPair[1], cands, size)) + continue; + + // photon + bestPhoScore = 0; + for (unsigned int i3 = 0; i3 < ngammas; ++i3) { + if (cands2[ig[i3]].pt() < cuts.minpt3) + continue; // photon pt cut + + std::array tr{{bestPair[0], bestPair[1], ig[i3]}}; // triplet of indices + + if (cands2[ig[i3]].pt() > bestPhoScore) { + std::copy_n(tr.begin(), 3, bestTriplet.begin()); + bestPhoScore = cands2[ig[i3]].pt(); + } + } + if (bestPhoScore < 0) + continue; + + // photon isolation + bool isop = isolationTkEm( + cands2[bestTriplet[2]].pt(), cands2[bestTriplet[2]].eta(), cands2[bestTriplet[2]].phi(), cands, size); + if (!isop) + continue; + + // H mass + auto mass = tripletmass(bestTriplet, cands, cands2, {{0.1396, 0.1396, 0.}}); + if (!(mass >= cuts.minmass and mass <= cuts.maxmass)) + continue; + + ret->emplace_back(bx); + nPass++; + masses.push_back(mass); + i0s.push_back(bestTriplet[0]); + i1s.push_back(bestTriplet[1]); + i2s.push_back(bestTriplet[2]); + bxOffsetsFiller.addBx(bx, 1); + } // loop on BXs + + iEvent.put(std::move(ret), "selectedBx" + label); + // now we make the table + auto bxOffsets = bxOffsetsFiller.done(); + auto tab = std::make_unique(bxOffsets, "HRhoGamma" + label, true); + tab->addColumn("mass", masses, "2 pions plus photon invariant mass"); + tab->addColumn("i0", i0s, "leading pion"); + tab->addColumn("i1", i1s, "subleading pion"); + tab->addColumn("i2", i2s, "photon"); + iEvent.put(std::move(tab), "hrhogamma" + label); +} + +//TEST functions +template +bool ScPhase2PuppiHRhoGammaDemo::isolationQ( + unsigned int pidex1, unsigned int pidex2, const T *cands, unsigned int size) const { + bool passed = false; + float psum = 0; + float eta = cands[pidex1].eta(); //center cone around leading track + float phi = cands[pidex1].phi(); + for (unsigned int j = 0u; j < size; ++j) { //loop over other particles + if (pidex1 == j or pidex2 == j) + continue; + float deta = eta - cands[j].eta(), dphi = ROOT::VecOps::DeltaPhi(phi, cands[j].phi()); + float dr2 = deta * deta + dphi * dphi; + if (dr2 >= cuts.mindr2 && dr2 <= cuts.maxdr2) + psum += cands[j].pt(); + } + if (psum <= cuts.maxiso * (cands[pidex1].pt() + cands[pidex2].pt())) + passed = true; + return passed; +} + +template +bool ScPhase2PuppiHRhoGammaDemo::isolationTkEm(float pt, float eta, float phi, const T *cands, unsigned int size) const { + bool passed = false; + float psum = 0; + for (unsigned int j = 0u; j < size; ++j) { //loop over other particles + float deta = eta - cands[j].eta(), dphi = ROOT::VecOps::DeltaPhi(phi, cands[j].phi()); + float dr2 = deta * deta + dphi * dphi; + if (dr2 >= cuts.mindr2tkem && dr2 <= cuts.maxdr2tkem) + psum += cands[j].pt(); + } + if (psum <= cuts.maxisotkem * pt) + passed = true; + return passed; +} + +std::tuple ScPhase2PuppiHRhoGammaDemo::deltar(float eta1, float eta2, float phi1, float phi2) const { + bool passed = true; + float deta = eta1 - eta2; + float dphi = ROOT::VecOps::DeltaPhi(phi1, phi2); + float dr2 = deta * deta + dphi * dphi; + if (dr2 > cuts.maxdeltar2) { + passed = false; + return std::tuple(passed, dr2); + } + return std::tuple(passed, dr2); +} + +template +float ScPhase2PuppiHRhoGammaDemo::pairmass(const std::array &t, + const T *cands, + const std::array &massD) { + ROOT::Math::PtEtaPhiMVector p1(cands[t[0]].pt(), cands[t[0]].eta(), cands[t[0]].phi(), massD[0]); + ROOT::Math::PtEtaPhiMVector p2(cands[t[1]].pt(), cands[t[1]].eta(), cands[t[1]].phi(), massD[1]); + float mass = (p1 + p2).M(); + return mass; +} + +template +float ScPhase2PuppiHRhoGammaDemo::tripletmass(const std::array &t, + const T *cands, + const U *cands2, + const std::array &masses) { + ROOT::Math::PtEtaPhiMVector p1(cands[t[0]].pt(), cands[t[0]].eta(), cands[t[0]].phi(), masses[0]); + ROOT::Math::PtEtaPhiMVector p2(cands[t[1]].pt(), cands[t[1]].eta(), cands[t[1]].phi(), masses[1]); + ROOT::Math::PtEtaPhiMVector p3(cands2[t[2]].pt(), cands2[t[2]].eta(), cands2[t[2]].phi(), masses[2]); + float mass = (p1 + p2 + p3).M(); + return mass; +} + +void ScPhase2PuppiHRhoGammaDemo::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + +DEFINE_FWK_MODULE(ScPhase2PuppiHRhoGammaDemo); diff --git a/L1TriggerScouting/Phase2/test/runScoutingPhase2HRare_cfg.py b/L1TriggerScouting/Phase2/test/runScoutingPhase2HRare_cfg.py new file mode 100644 index 0000000000000..4d596c2edacc0 --- /dev/null +++ b/L1TriggerScouting/Phase2/test/runScoutingPhase2HRare_cfg.py @@ -0,0 +1,355 @@ +from __future__ import print_function +import FWCore.ParameterSet.Config as cms +import FWCore.ParameterSet.VarParsing as VarParsing +import os + +options = VarParsing.VarParsing ('analysis') +options.register ('runNumber', + 37, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.int, # string, int, or float + "Run Number") + +options.register ('lumiNumber', + 1, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.int, # string, int, or float + "Run Number") + +options.register ('daqSourceMode', + 'ScoutingPhase2', # default value + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, # string, int, or float + "DAQ source data mode") + +options.register ('broker', + 'none', # default value + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, # string, int, or float + "Broker: 'none' or 'hostname:port'") + +options.register ('buBaseDir', + '/dev/shm/ramdisk', # default value + VarParsing.VarParsing.multiplicity.list, + VarParsing.VarParsing.varType.string, # string, int, or float + "BU base directory") + +options.register ('buNumStreams', + "2", # default value = 2: puppi, tkem + VarParsing.VarParsing.multiplicity.list, + VarParsing.VarParsing.varType.int, # string, int, or float + "Number of input streams (i.e. files) used simultaneously for each BU directory") + +options.register ('timeslices', + 1, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.int, # string, int, or float + "Number of timeslices") + +options.register ('tmuxPeriod', + 1, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.int, # string, int, or float + "Time multiplex period") + +options.register ('puppiStreamIDs', + [], # default value + VarParsing.VarParsing.multiplicity.list, + VarParsing.VarParsing.varType.int, # string, int, or float + "Stream IDs for the Puppi inputs") + +options.register ('tkEmStreamIDs', + [], # default value + VarParsing.VarParsing.multiplicity.list, + VarParsing.VarParsing.varType.int, # string, int, or float + "Stream IDs for the TkEm inputs") + +options.register ('fuBaseDir', + '/dev/shm/data', # default value + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, # string, int, or float + "BU base directory") + +options.register ('fffBaseDir', + '/dev/shm', # default value + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, # string, int, or float + "FFF base directory") + +options.register ('numThreads', + 1, # default value + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.int, # string, int, or float + "Number of CMSSW threads") + +options.register ('numFwkStreams', + 1, # default value + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.int, # string, int, or float + "Number of CMSSW streams") + +options.register ('run', + 'both', # default value + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, # string, int, or float + "'inclusive', 'selected', 'both' (default).") + +options.register ('analyses', + [], # default value + VarParsing.VarParsing.multiplicity.list, + VarParsing.VarParsing.varType.string, # string, int, or float + "analyses: any list of 'w3pi', 'wdsg'.") + +options.register ('outMode', + 'none', # default value + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, # string, int, or float + "output (none, nanoSelected, nanoInclusive, nanoBoth)") + +options.register ('outFile', + "NanoOutput.root", + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, + "Sub lumisection number to process") + + +options.parseArguments() +analyses = options.analyses if options.analyses else ["hrhog", "hphig", "hjpsig", "h2rho", "h2phi"] +print(f"Analyses set to {analyses}") + +process = cms.Process("SCPU") +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(options.maxEvents) +) + +process.options = cms.untracked.PSet( + numberOfThreads = cms.untracked.uint32(options.numThreads), + numberOfStreams = cms.untracked.uint32(options.numFwkStreams), + numberOfConcurrentLuminosityBlocks = cms.untracked.uint32(1), + wantSummary = cms.untracked.bool(True) +) +process.load("FWCore.MessageService.MessageLogger_cfi") +process.MessageLogger.cerr.FwkReport.reportEvery = 100 + +if len(options.buNumStreams) != len(options.buBaseDir): + raise RuntimeError("Mismatch between buNumStreams (%d) and buBaseDirs (%d)" % (len(options.buNumStreams), len(options.buBaseDir))) + +if options.puppiStreamIDs == [] and options.tkEmStreamIDs == []: + nStreamsTot = sum(options.buNumStreams) + puppiStreamIDs = list(range(nStreamsTot//2)) # take first half + tkEmStreamIDs = list(range(nStreamsTot//2, nStreamsTot)) # take second half +else: + puppiStreamIDs = options.puppiStreamIDs + tkEmStreamIDs = options.tkEmStreamIDs + +process.EvFDaqDirector = cms.Service("EvFDaqDirector", + useFileBroker = cms.untracked.bool(options.broker != "none"), + fileBrokerHostFromCfg = cms.untracked.bool(False), + fileBrokerHost = cms.untracked.string(options.broker.split(":")[0] if options.broker != "none" else "htcp40.cern.ch"), + fileBrokerPort = cms.untracked.string(options.broker.split(":")[1] if options.broker != "none" else "8080"), + runNumber = cms.untracked.uint32(options.runNumber), + baseDir = cms.untracked.string(options.fuBaseDir), + buBaseDir = cms.untracked.string(options.buBaseDir[0]), + buBaseDirsAll = cms.untracked.vstring(*options.buBaseDir), + buBaseDirsNumStreams = cms.untracked.vint32(*options.buNumStreams), + directorIsBU = cms.untracked.bool(False), +) +process.FastMonitoringService = cms.Service("FastMonitoringService") + +process.load( "HLTrigger.Timer.FastTimerService_cfi" ) +process.FastTimerService.writeJSONSummary = cms.untracked.bool(True) +process.FastTimerService.jsonFileName = cms.untracked.string(f'resources.{os.uname()[1]}.json') +process.MessageLogger.cerr.FastReport = cms.untracked.PSet( limit = cms.untracked.int32( 10000000 ) ) + +fuDir = options.fuBaseDir+("/run%06d" % options.runNumber) +buDirs = [b+("/run%06d" % options.runNumber) for b in options.buBaseDir] +for d in [fuDir, options.fuBaseDir] + buDirs + options.buBaseDir: + if not os.path.isdir(d): + os.makedirs(d) + +process.source = cms.Source("DAQSource", + testing = cms.untracked.bool(True), + dataMode = cms.untracked.string(options.daqSourceMode), + verifyChecksum = cms.untracked.bool(True), + useL1EventID = cms.untracked.bool(False), + eventChunkBlock = cms.untracked.uint32(2 * 1024), + eventChunkSize = cms.untracked.uint32(2 * 1024), + maxChunkSize = cms.untracked.uint32(4 * 1024), + numBuffers = cms.untracked.uint32(4), + maxBufferedFiles = cms.untracked.uint32(4), + fileListMode = cms.untracked.bool(options.broker == "none"), + fileNames = cms.untracked.vstring( + buDirs[0] + "/" + "run%06d_ls%04d_index%06d_stream00.raw" % (options.runNumber, options.lumiNumber, 1), + ) +) +os.system("touch " + buDirs[0] + "/" + "fu.lock") + +## test pluging +process.scPhase2PuppiRawToDigiStruct = cms.EDProducer('ScPhase2PuppiRawToDigi', + src = cms.InputTag('rawDataCollector'), + fedIDs = cms.vuint32(*puppiStreamIDs), +) + +process.scPhase2TkEmRawToDigiStruct = cms.EDProducer('ScPhase2TkEmRawToDigi', + src = cms.InputTag('rawDataCollector'), + fedIDs = cms.vuint32(*tkEmStreamIDs), +) + +process.goodOrbitsByNBX = cms.EDFilter("GoodOrbitNBxSelector", + unpackers = cms.VInputTag( + cms.InputTag("scPhase2PuppiRawToDigiStruct"), + cms.InputTag("scPhase2TkEmRawToDigiStruct"), + ), + nbxMin = cms.uint32(3564 * options.timeslices // options.tmuxPeriod) +) + +process.hrhogStruct = cms.EDProducer("ScPhase2PuppiHRhoGammaDemo", + src = cms.InputTag("scPhase2PuppiRawToDigiStruct"), + src2 = cms.InputTag("scPhase2TkEmRawToDigiStruct"), + runStruct = cms.bool(True) +) + +process.hphigStruct = cms.EDProducer("ScPhase2PuppiHPhiGammaDemo", + src = cms.InputTag("scPhase2PuppiRawToDigiStruct"), + src2 = cms.InputTag("scPhase2TkEmRawToDigiStruct"), + runStruct = cms.bool(True) +) + +process.hjpsigStruct = cms.EDProducer("ScPhase2PuppiHJPsiGammaDemo", + src = cms.InputTag("scPhase2PuppiRawToDigiStruct"), + src2 = cms.InputTag("scPhase2TkEmRawToDigiStruct"), + runStruct = cms.bool(True) +) + +process.h2rhoStruct = cms.EDProducer("ScPhase2PuppiH2RhoDemo", + src = cms.InputTag("scPhase2PuppiRawToDigiStruct"), + runStruct = cms.bool(True) +) + +process.h2phiStruct = cms.EDProducer("ScPhase2PuppiH2PhiDemo", + src = cms.InputTag("scPhase2PuppiRawToDigiStruct"), + runStruct = cms.bool(True) +) + +process.hphijpsiStruct = cms.EDProducer("ScPhase2PuppiHPhiJPsiDemo", + src = cms.InputTag("scPhase2PuppiRawToDigiStruct"), + runStruct = cms.bool(True) +) + +process.scPhase2SelectedBXs = cms.EDFilter("FinalBxSelector", + analysisLabels = cms.VInputTag([cms.InputTag(f"{a}Struct", "selectedBx") for a in analyses]), +) + +process.scPhase2PuppiMasked = cms.EDProducer("MaskOrbitBxScoutingPuppi", + dataTag = cms.InputTag("scPhase2PuppiRawToDigiStruct"), + selectBxs = cms.InputTag("scPhase2SelectedBXs","SelBx"), +) + +process.scPhase2TkEmMasked = cms.EDProducer("MaskOrbitBxScoutingTkEm", + dataTag = cms.InputTag("scPhase2TkEmRawToDigiStruct"), + selectBxs = cms.InputTag("scPhase2SelectedBXs","SelBx"), +) + +process.scPhase2TkEleMasked = cms.EDProducer("MaskOrbitBxScoutingTkEle", + dataTag = cms.InputTag("scPhase2TkEmRawToDigiStruct"), + selectBxs = cms.InputTag("scPhase2SelectedBXs","SelBx"), +) + +process.scPhase2PuppiStructToTable = cms.EDProducer("ScPuppiToOrbitFlatTable", + src = cms.InputTag("scPhase2PuppiRawToDigiStruct"), + name = cms.string("L1Puppi"), + doc = cms.string("L1Puppi candidates from Correlator Layer 2"), +) + +process.scPhase2PuppiMaskedStructToTable = process.scPhase2PuppiStructToTable.clone( + src = "scPhase2PuppiMasked" +) + +process.scPhase2TkEmStructToTable = cms.EDProducer("ScTkEmToOrbitFlatTable", + src = cms.InputTag("scPhase2TkEmRawToDigiStruct"), + name = cms.string("L1TkEm"), + doc = cms.string("L1TkEm candidates"), +) + +process.scPhase2TkEmMaskedStructToTable = process.scPhase2TkEmStructToTable.clone( + src = "scPhase2TkEmMasked" +) + +process.scPhase2TkEleStructToTable = cms.EDProducer("ScTkEleToOrbitFlatTable", + src = cms.InputTag("scPhase2TkEmRawToDigiStruct"), + name = cms.string("L1TkEle"), + doc = cms.string("L1TkEle candidates"), +) + +process.scPhase2TkEleMaskedStructToTable = process.scPhase2TkEleStructToTable.clone( + src = "scPhase2TkEleMasked" +) + +process.s_unpackers = cms.Sequence( + process.scPhase2PuppiRawToDigiStruct + + process.scPhase2TkEmRawToDigiStruct + + process.goodOrbitsByNBX +) + +process.p_inclusive = cms.Path( + process.s_unpackers + + process.scPhase2PuppiStructToTable + + process.scPhase2TkEmStructToTable + + process.scPhase2TkEleStructToTable +) + +analysisModules = [getattr(process,f"{a}Struct") for a in analyses] +process.s_analyses = cms.Sequence(sum(analysisModules[1:], analysisModules[0])) + +process.p_selected = cms.Path( + process.s_unpackers + + process.s_analyses + + process.scPhase2SelectedBXs + + process.scPhase2PuppiMasked + + process.scPhase2TkEmMasked + + process.scPhase2TkEleMasked + + process.scPhase2PuppiMaskedStructToTable + + process.scPhase2TkEmMaskedStructToTable + + process.scPhase2TkEleMaskedStructToTable +) + +process.scPhase2NanoAll = cms.OutputModule("OrbitNanoAODOutputModule", + fileName = cms.untracked.string(options.outFile.replace(".root","")+".inclusive.root"), + SelectEvents = cms.untracked.PSet(SelectEvents = cms.vstring('p_inclusive')), + outputCommands = cms.untracked.vstring("drop *", + "keep l1ScoutingRun3OrbitFlatTable_scPhase2PuppiStructToTable_*_*", + "keep l1ScoutingRun3OrbitFlatTable_scPhase2TkEmStructToTable_*_*", + "keep l1ScoutingRun3OrbitFlatTable_scPhase2TkEleStructToTable_*_*"), + compressionLevel = cms.untracked.int32(4), + compressionAlgorithm = cms.untracked.string("LZ4"), +) + +process.scPhase2PuppiNanoSelected = cms.OutputModule("OrbitNanoAODOutputModule", + fileName = cms.untracked.string(options.outFile.replace(".root","")+".selected.root"), + SelectEvents = cms.untracked.PSet(SelectEvents = cms.vstring('p_selected')), + selectedBx = cms.InputTag("scPhase2SelectedBXs","SelBx"), + outputCommands = cms.untracked.vstring("drop *", + "keep l1ScoutingRun3OrbitFlatTable_scPhase2PuppiMaskedStructToTable_*_*", + "keep l1ScoutingRun3OrbitFlatTable_scPhase2TkEmMaskedStructToTable_*_*", + "keep l1ScoutingRun3OrbitFlatTable_scPhase2TkEleMaskedStructToTable_*_*", + "keep *_hrhogStruct_*_*", + "keep *_hphigStruct_*_*", + "keep *_hjpsigStruct_*_*", + "keep *_h2rhoStruct_*_*", + "keep *_h2phiStruct_*_*", + "keep *_scPhase2SelectedBXs_*_*" + ), + compressionLevel = cms.untracked.int32(4), + compressionAlgorithm = cms.untracked.string("LZ4"), +) +process.o_nanoInclusive = cms.EndPath(process.scPhase2NanoAll) +process.o_nanoSelected = cms.EndPath(process.scPhase2PuppiNanoSelected) +process.o_nanoBoth = cms.EndPath(process.scPhase2NanoAll + process.scPhase2PuppiNanoSelected) + +sched = [ process.p_inclusive, process.p_selected ] +if options.run == "inclusive": sched = [ process.p_inclusive ] +if options.run == "selected": sched = [ process.p_selected ] + +if options.outMode != "none": + sched.append(getattr(process, "o_"+options.outMode)) +process.schedule = cms.Schedule(*sched)