From 5c715898a5a239d42e676cd92c5be306eb61c68f Mon Sep 17 00:00:00 2001 From: tylerhoroho-UVA Date: Fri, 25 Oct 2024 11:33:57 -0500 Subject: [PATCH] First steps to build a visibles veto processor --- Hcal/include/Hcal/Event/VisiblesVetoResult.h | 79 ++++++++ Hcal/include/Hcal/VisiblesVetoProcessor.h | 64 ++++++ Hcal/src/Hcal/Event/VisiblesVetoResult.cxx | 55 +++++ Hcal/src/Hcal/VisiblesVetoProcessor.cxx | 202 +++++++++++++++++++ 4 files changed, 400 insertions(+) create mode 100644 Hcal/include/Hcal/Event/VisiblesVetoResult.h create mode 100644 Hcal/include/Hcal/VisiblesVetoProcessor.h create mode 100644 Hcal/src/Hcal/Event/VisiblesVetoResult.cxx create mode 100644 Hcal/src/Hcal/VisiblesVetoProcessor.cxx diff --git a/Hcal/include/Hcal/Event/VisiblesVetoResult.h b/Hcal/include/Hcal/Event/VisiblesVetoResult.h new file mode 100644 index 000000000..fcdcbc7a5 --- /dev/null +++ b/Hcal/include/Hcal/Event/VisiblesVetoResult.h @@ -0,0 +1,79 @@ +#ifndef EVENT_VISIBLESVETORESULT_H_ +#define EVENT_VISIBLESVETORESULT_H_ + +#include + +// ROOT // +#include + +namespace ldmx { + + class VisiblesVetoResult { + public: + /** Constructor */ + VisiblesVetoResult(); + + /** Destructor */ + virtual ~VisiblesVetoResult(); + + void setVariables(int nLayersHit, + double xStd, double yStd, double zStd, + double xMean, double yMean, double rMean, + int isoHits, double isoEnergy, + int nReadoutHits, double summedDet, + double rMeanFromPhotonProj); + + void Clear(); + + void Print() const; + + bool passesVeto() const { return passesVeto_; } + + double getDisc() const { return discValue_; } + + int getNLayersHit() const { return nLayersHit_; } + + double getXStd() const { return xStd_; } + + double getYStd() const { return yStd_; } + + double getZStd() const { return zStd_; } + + double getXMean() const { return xMean_; } + + double getYMean() const { return yMean_; } + + double getRMean() const { return rMean_; } + + int getIsoHits() const { return isoHits_; } + + double getIsoEnergy() const { return isoEnergy_; } + + int getNHits() const { return nReadoutHits_; } + + double getSummedDet() const { return summedDet_; } + + double getDistFromPhotonProj() const { return rMeanFromPhotonProj_; } + + private: + bool passesVeto_{false}; + + int nLayersHit_{0}; + double xStd_{0}; + double yStd_{0}; + double zStd_{0}; + double xMean_{0}; + double yMean_{0}; + double rMean_{0}; + int isoHits_{0}; + double isoEnergy_{0}; + int nReadoutHits_{0}; + double summedDet_{0}; + double rMeanFromPhotonProj_{0}; + + double discValue_{0}; + + }; +} // namespace ldmx + +#endif diff --git a/Hcal/include/Hcal/VisiblesVetoProcessor.h b/Hcal/include/Hcal/VisiblesVetoProcessor.h new file mode 100644 index 000000000..a9f96767d --- /dev/null +++ b/Hcal/include/Hcal/VisiblesVetoProcessor.h @@ -0,0 +1,64 @@ +#ifndef EVENTPROC_VISIBLESVETOPROCESSOR_H_ +#define EVENTPROC_VISIBLESVETOPROCESSOR_H_ + +// LDMX +#include "Hcal/Event/HcalHit.h" +#include "Hcal/Event/VisiblesVetoResult.h" +#include "Framework/Configure/Parameters.h" +#include "Framework/EventProcessor.h" +#include "Tools/ONNXRuntime.h" + +namespace hcal { + + class VisiblesVetoProcessor : public framework::Producer { + public: + EcalVetoProcessor(const std::string& name, framework::Process& process) + : Producer(name, process) {} + + virtual ~VisiblesVetoProcessor() {} + + void configure(framework::config::Parameters& parameters) override; + + void produce(framework::Event& event) override; + + private: + void clearProcessor(); + + void buildBDTFeatureVector(const ldmx::VisiblesVetoResult& result); + + int nLayersHit_{0}; + double xStd_{0}; + double yStd_{0}; + double zStd_{0}; + double xMean_{0}; + double yMean_{0}; + double rMean_{0}; + int isoHits_{0}; + double isoEnergy_{0}; + int nReadoutHits_{0}; + double summedDet_{0}; + double rMeanFromPhotonProj_{0}; + + double bdtCutVal_{0}; + + double beamEnergyMeV_{0}; + + bool verbose_{false}; + + std::string bdtFileName_; + std::string rocFileName_; + std::vector bdtFeatures_; + std::string featureListName_; + + std::string rec_pass_name_; + std::string rec_coll_name_; + + std::string collectionName_{"VisiblesVeto"}; + + std::unique_ptr rt_; + + }; + +} // namespace hcal + +#endif diff --git a/Hcal/src/Hcal/Event/VisiblesVetoResult.cxx b/Hcal/src/Hcal/Event/VisiblesVetoResult.cxx new file mode 100644 index 000000000..b8d44c9ce --- /dev/null +++ b/Hcal/src/Hcal/Event/VisiblesVetoResult.cxx @@ -0,0 +1,55 @@ +#include "Hcal/Event/VisiblesVetoResult.h" + +ClassImp(ldmx::VisiblesVetoResult); + +namespace ldmx { + VisiblesVetoResult::VisiblesVetoResult() {} + + VisiblesVetoResult::~VisiblesVetoResult() { Clear(); } + + void VisiblesVetoResult::Clear() { + passesVeto_ = false; + + nLayersHit_ = 0; + xStd_ = 0.; + yStd_ = 0.; + zStd_ = 0.; + xMean_ = 0.; + yMean_ = 0.; + rMean_ = 0.; + isoHits_ = 0; + isoEnergy_ = 0.; + nReadoutHits_ = 0; + summedDet_ = 0.; + rMeanFromPhotonProj_ = 0.; + + discValue_ = 0.; + } + + void VisiblesVetoResult::setVariables( + int nLayersHit, + double xStd, double yStd, double zStd, + double xMean, double yMean, double rMean, + int isoHits, double isoEnergy, + int nReadoutHits, double summedDet, + double rMeanFromPhotonProj) { + nLayersHit_ = nLayersHit; + xStd_ = xStd; + yStd_ = yStd; + zStd_ = zStd; + xMean_ = xMean; + yMean_ = yMean; + rMean_ = rMean; + isoHits_ = isoHits; + isoEnergy_ = isoEnergy; + nReadoutHits_ = nReadoutHits; + summedDet_ = summedDet; + rMeanFromPhotonProj_ = rMeanFromPhotonProj; + } + + void VisiblesVetoResult::Print() const { + std::cout << "[ VisiblesVetoResult ]:\n" + << "\t Passes veto : " << passesVeto_ << "\n" + << std::endl; + } +} // namespace ldmx diff --git a/Hcal/src/Hcal/VisiblesVetoProcessor.cxx b/Hcal/src/Hcal/VisiblesVetoProcessor.cxx new file mode 100644 index 000000000..d81a87427 --- /dev/null +++ b/Hcal/src/Hcal/VisiblesVetoProcessor.cxx @@ -0,0 +1,202 @@ +#include "Hcal/VisiblesVetoProcessor.h" + +// LDMX +#include "DetDescr/HcalID.h" +#include "DetDescr/SimSpecialID.h" +#include "SimCore/Event/SimParticle.h" +#include "SimCore/Event/SimTrackerHit.h" +#include "Hcal/Event/HcalHit.h" + +// C++ +#include +#include +#include +#include +#include + +namespace hcal { + void VisiblesVetoProcessor::buildBDTFeatureVector(const ldmx::VisiblesVetoResult &result) { + bdtFeatures_.push_back(result.getLayersHit()); + bdtFeatures_.push_back(result.getXStd()); + bdtFeatures_.push_back(result.getYStd()); + bdtFeatures_.push_back(result.getZStd()); + bdtFeatures_.push_back(result.getXMean()); + bdtFeatures_.push_back(result.getYMean()); + bdtFeatures_.push_back(result.getRMean()); + bdtFeatures_.push_back(result.getIsoHits()); + bdtFeatures_.push_back(result.getIsoEnergy()); + bdtFeatures_.push_back(result.getNReadoutHits()); + bdtFeatures_.push_back(result.getSummedDet()); + bdtFeatures_.push_back(result.getMeanDistFromPhotonProj()); + } + + void VisiblesVetoProcessor::configure(framework::config::Parameters ¶meters) { + verbose_ = parameters.getParameter("verbose"); + featureListName_ = parameters.getParameter("feature_list_name"); + // Load BDT ONNX file + rt_ = std::make_unique(parameters.getParameter("bdt_file")); + } + + void VisiblesVetoProcessor::clearProcessor() { + bdtFeatures_.clear(); + + nLayersHit_ = 0; + xStd_ = 0.; + yStd_ = 0.; + zStd_ = 0.; + xMean_ = 0.; + yMean_ = 0.; + rMean_ = 0.; + isoHits_ = 0; + isoEnergy_ = 0.; + nReadoutHits = 0; + summedDet_ = 0.; + rMeanFromPhotonProj_ = 0.; + } + + void VisiblesVetoProcessor::produce(framework::Event &event) { + + ldmx::VisiblesVetoResult result; + + clearProcessor(); + + auto particle_map{event.getMap("SimParticles")}; + + // Get target scoring plane hits for recoil electron + // Use this to calculate the projected photon line vector + // This currently uses truth-level information, but it should be replaced + // by reconstructed tracker information, when available + std::vector gamma_p(3); + std::vector gamma_x0(3); + if (event.exists("TargetScoringPlaneHits")) { + std::vector targetSPHits = + event.getCollection("TargetScoringPlaneHits"); + for (auto const &it : particle_map) { + for (auto const &sphit : targetSPHits) { + if (sphit.getPosition()[2] > 0) { + if (it.first == sphit.getTrackID()) { + if (it.second.getPdgID() == 11 && in_list(it.second.getParents(), 0)) { + gamma_x0 = sphit.getPosition(); + gamma_p[0] = -1*sphit.getMomentum()[0]; + gamma_p[1] = -1*sphit.getMomentum()[1]; + gamma_p[2] = 8000 - sphit.getMomentum()[2]; // hard-coded for 8 GeV + } + } + } + } + } + } + + // Get Hcal reconstructed hits and loop through them to build features + const std::vector hcalRecHits = + event.getCollection(rec_coll_name_, rec_pass_name_); + + double zMean = 0.; // need this when calculating zStd_ + std::vector layersHit; + for (const ldmx::HcalHit & hit : hcalRecHits) { + if (hit.getEnergy() > 0.) { + HcalID detID(hit.getID()); + if (detID.getSection() == 0) { // looking for hits in the back Hcal + nReadoutHits_ += 1; + double x = hit.getXPos(); + double y = hit.getYPos(); + double z = hit.getZPos(); + double r = sqrt(pow(x,2) + pow(y,2)); + + summedDet_ += hit.getEnergy(); + + xMean_ += x*hit.getEnergy(); + yMean_ += y*hit.getEnergy(); + zMean += z*hit.getEnergy(); + rMean_ += r*hit.getEnergy(); + + // check if this is a new layer + if (!(std::find(layersHit.begin(), layersHit.end(), detID.getLayerID()) != layersHit.end())) { + layersHit.push_back(detID.getLayerID()); + } + + double x_proj = gamma_x0[0] + (z - gamma_x0[2])*gamma_p[0]/gamma_p[2]; + double y_proj = gamma_x0[1] + (z - gamma_x0[2])*gamma_p[1]/gamma_p[2]; + + rMeanFromPhotonProj_ += hit.getEnergy()*sqrt(pow(x-x_proj,2) + pow(y-y_proj,2)); + + // Calculate isolated hits + double closestpoint = 9999.; + for (const ldmx::HcalHit &hit2 : hcalRecHits) { + if (hit2.getEnergy() > 0.) { + HcalID detID2(hit2.getID()); + if (detID2.getLayer() == detID.getLayer()) { + // Determine if a bar is vertical (along y-axis) or horizontal (along x-axis) + // Odd layers have horizontal strips + // Even layers have vertical strips + if (detID2.getLayer() % 2 == 0) { + if (abs(hit2.getYPos() - y) > 0) { + if (abs(hit2.getYPos() - y) < closestpoint) { + closestpoint = abs(hit2.getYPos() - y); + } + } + } + else { + if (abs(hit2.getXPos() - x) > 0) { + if (abs(hit2.getXPos() - x) < closestpoint) { + closestpoint = abs(hit2.getXPos() - x); + } + } + } + } + } + } + if (closestpoint > 50.) { + isoHits_ += 1; + isoEnergy_ += hit.getEnergy(); + } + } + } + } + + nLayersHit_ = layersHit.size(); + + if (summedDet_ > 0.) { + xMean_ /= summedDet_; + yMean_ /= summedDet_; + zMean /= summedDet_; + rMean_ /= summedDet_; + + rMeanFromPhotonProj_ /= summedDet_; + } + + for (const ldmx::HcalHit &hit : hcalRecHits) { + if (hit.getEnergy() > 0.) { + HcalID detID(hit.getID()); + if (detID.getSection() == 0) { + xStd_ += hit.getEnergy()*pow(hit.getXPos()-xMean_,2); + yStd_ += hit.getEnergy()*pow(hit.getYPos()-yMean_,2); + zStd_ += hit.getEnergy()*pow(hit.getZPos()-zMean ,2); + } + } + } + + if (summedDet_ > 0.) { + xStd_ = sqrt(xStd_/summedDet_); + yStd_ = sqrt(yStd_/summedDet_); + zStd_ = sqrt(zStd_/summedDet_); + } + + result.setVariables( + nLayersHit_, + xStd_, + yStd_, + zStd_, + xMean_, + yMean_, + rMean_, + isoHits_, + isoEnergy_, + nReadoutHits_, + summedDet_, + rMeanFromPhotonProj_); + + buildBDTFeatureVector(result); + } + +}