ecal::EcalVetoProcessor Class Reference

Determines if event is vetoable using ECAL hit information. More...

#include <EcalVetoProcessor.h>

Classes
struct	HitData

Public Types
typedef std::pair< ldmx::EcalID, float >	CellEnergyPair
typedef std::pair< float, float >	XYCoords

Public Member Functions
	EcalVetoProcessor (const std::string &name, framework::Process &process)
void	configure (framework::config::Parameters &parameters) override
	Configure the processor using the given user specified parameters.
void	produce (framework::Event &event) override
	Process the event and put new data products into it.
Public Member Functions inherited from framework::Producer
	Producer (const std::string &name, Process &process)
	Class constructor.
virtual void	beforeNewRun (ldmx::RunHeader &header)
	Handle allowing producers to modify run headers before the run begins.
Public Member Functions inherited from framework::EventProcessor
	EventProcessor (const std::string &name, Process &process)
	Class constructor.
virtual	~EventProcessor ()
	Class destructor.
virtual void	onNewRun (const ldmx::RunHeader &runHeader)
	Callback for the EventProcessor to take any necessary action when the run being processed changes.
virtual void	onFileOpen (EventFile &eventFile)
	Callback for the EventProcessor to take any necessary action when a new event input ROOT file is opened.
virtual void	onFileClose (EventFile &eventFile)
	Callback for the EventProcessor to take any necessary action when a event input ROOT file is closed.
virtual void	onProcessStart ()
	Callback for the EventProcessor to take any necessary action when the processing of events starts, such as creating histograms.
virtual void	onProcessEnd ()
	Callback for the EventProcessor to take any necessary action when the processing of events finishes, such as calculating job-summary quantities.
template<class T >
const T &	getCondition (const std::string &condition_name)
	Access a conditions object for the current event.
TDirectory *	getHistoDirectory ()
	Access/create a directory in the histogram file for this event processor to create histograms and analysis tuples.
void	setStorageHint (framework::StorageControl::Hint hint)
	Mark the current event as having the given storage control hint from this module.
void	setStorageHint (framework::StorageControl::Hint hint, const std::string &purposeString)
	Mark the current event as having the given storage control hint from this module and the given purpose string.
int	getLogFrequency () const
	Get the current logging frequency from the process.
int	getRunNumber () const
	Get the run number from the process.
std::string	getName () const
	Get the processor name.
void	createHistograms (const std::vector< framework::config::Parameters > &histos)
	Internal function which is used to create histograms passed from the python configuration @parma histos vector of Parameters that configure histograms to create.

Private Member Functions
void	clearProcessor ()
ldmx::EcalID	GetShowerCentroidIDAndRMS (const std::vector< ldmx::EcalHit > &ecalRecHits, double &showerRMS)
void	fillHitMap (const std::vector< ldmx::EcalHit > &ecalRecHits, std::map< ldmx::EcalID, float > &cellMap_)
	Function to load up empty vector of hit maps.
void	fillIsolatedHitMap (const std::vector< ldmx::EcalHit > &ecalRecHits, ldmx::EcalID globalCentroid, std::map< ldmx::EcalID, float > &cellMap, std::map< ldmx::EcalID, float > &cellMapIso, bool doTight=false)
std::vector< XYCoords >	getTrajectory (std::vector< double > momentum, std::vector< float > position)
void	buildBDTFeatureVector (const ldmx::EcalVetoResult &result)
float	distTwoLines (TVector3 v1, TVector3 v2, TVector3 w1, TVector3 w2)
	Returns the distance between the lines v and w, with v defined to pass through the points (v1,v2) (and similarly for w).
float	distPtToLine (TVector3 h1, TVector3 p1, TVector3 p2)
	Return the minimum distance between the point h1 and the line passing through points p1 and p2.

Private Attributes
std::map< ldmx::EcalID, float >	cellMap_
std::map< ldmx::EcalID, float >	cellMapTightIso_
std::vector< float >	ecalLayerEdepRaw_
std::vector< float >	ecalLayerEdepReadout_
std::vector< float >	ecalLayerTime_
std::vector< std::vector< double > >	roc_range_values_
int	nEcalLayers_ {0}
int	nReadoutHits_ {0}
int	deepestLayerHit_ {0}
double	summedDet_ {0}
double	summedTightIso_ {0}
double	maxCellDep_ {0}
double	showerRMS_ {0}
double	xStd_ {0}
double	yStd_ {0}
double	avgLayerHit_ {0}
double	stdLayerHit_ {0}
double	ecalBackEnergy_ {0}
int	nStraightTracks_ {0}
	Number of "straight" tracks found in the event.
int	nLinregTracks_ {0}
	Number of "linreg" tracks found in the event.
int	firstNearPhLayer_ {0}
	First ECal layer in which a hit is found near the photon.
int	nNearPhHits_ {0}
	Number of hits near the photon trajectory.
float	epAng_ {0}
	Angular separation between the projected photon and electron trajectories (currently unused)
float	epSep_ {0}
	Distance between the projected photon and electron trajectories at the ECal face.
float	epDot_ {0}
	Dot product of the photon and electron momenta unit vectors.
int	photonTerritoryHits_ {0}
	Number of hits in the photon territory.
double	bdtCutVal_ {0}
double	beamEnergyMeV_ {0}
double	linreg_radius_ {0}
bool	verbose_ {false}
std::string	bdtFileName_
std::string	rocFileName_
std::vector< float >	bdtFeatures_
std::string	featureListName_
std::string	rec_pass_name_
std::string	rec_coll_name_
std::string	collectionName_ {"EcalVeto"}
	Name of the collection which will containt the results.
std::unique_ptr< ldmx::Ort::ONNXRuntime >	rt_
const ldmx::EcalGeometry *	geometry_
	handle to current geometry (to share with member functions)

Additional Inherited Members
Static Public Member Functions inherited from framework::EventProcessor
static void	declare (const std::string &classname, int classtype, EventProcessorMaker *)
	Internal function which is part of the PluginFactory machinery.
Static Public Attributes inherited from framework::Producer
static const int	CLASSTYPE {1}
	Constant used to track EventProcessor types by the PluginFactory.
Protected Member Functions inherited from framework::EventProcessor
void	abortEvent ()
	Abort the event immediately.
Protected Attributes inherited from framework::EventProcessor
HistogramHelper	histograms_
	Interface class for making and filling histograms.
NtupleManager &	ntuple_ {NtupleManager::getInstance()}
	Manager for any ntuples.
logging::logger	theLog_
	The logger for this EventProcessor.

Detailed Description

Determines if event is vetoable using ECAL hit information.

Definition at line 32 of file EcalVetoProcessor.h.

Member Typedef Documentation

CellEnergyPair

typedef std::pair<ldmx::EcalID, float> ecal::EcalVetoProcessor::CellEnergyPair

Definition at line 34 of file EcalVetoProcessor.h.

XYCoords

typedef std::pair<float, float> ecal::EcalVetoProcessor::XYCoords

Definition at line 36 of file EcalVetoProcessor.h.

Constructor & Destructor Documentation

EcalVetoProcessor()

ecal::EcalVetoProcessor::EcalVetoProcessor ( const std::string &  name, framework::Process &  process  )

inline

Definition at line 38 of file EcalVetoProcessor.h.

      : Producer(name, process) {}

~EcalVetoProcessor()

virtual ecal::EcalVetoProcessor::~EcalVetoProcessor ( )

inlinevirtual

Definition at line 41 of file EcalVetoProcessor.h.

{}

Member Function Documentation

buildBDTFeatureVector()

void ecal::EcalVetoProcessor::buildBDTFeatureVector ( const ldmx::EcalVetoResult &  result )

private

Electron RoC variables

Photon RoC variables

Outside RoC variables

Definition at line 30 of file EcalVetoProcessor.cxx.

                                      {
  // Base variables
  bdtFeatures_.push_back(result.getNReadoutHits());
  bdtFeatures_.push_back(result.getSummedDet());
  bdtFeatures_.push_back(result.getSummedTightIso());
  bdtFeatures_.push_back(result.getMaxCellDep());
  bdtFeatures_.push_back(result.getShowerRMS());
  bdtFeatures_.push_back(result.getXStd());
  bdtFeatures_.push_back(result.getYStd());
  bdtFeatures_.push_back(result.getAvgLayerHit());
  bdtFeatures_.push_back(result.getStdLayerHit());
  bdtFeatures_.push_back(result.getDeepestLayerHit());
  bdtFeatures_.push_back(result.getEcalBackEnergy());
  // MIP tracking
  bdtFeatures_.push_back(result.getNStraightTracks());
  // bdtFeatures_.push_back(result.getNLinregTracks());
  bdtFeatures_.push_back(result.getFirstNearPhLayer());
  bdtFeatures_.push_back(result.getNNearPhHits());
  bdtFeatures_.push_back(result.getPhotonTerritoryHits());
  bdtFeatures_.push_back(result.getEPSep());
  bdtFeatures_.push_back(result.getEPDot());
  // Longitudinal segment variables
  bdtFeatures_.push_back(result.getEnergySeg()[0]);
  bdtFeatures_.push_back(result.getXMeanSeg()[0]);
  bdtFeatures_.push_back(result.getYMeanSeg()[0]);
  bdtFeatures_.push_back(result.getLayerMeanSeg()[0]);
  bdtFeatures_.push_back(result.getEnergySeg()[1]);
  bdtFeatures_.push_back(result.getYMeanSeg()[2]);
  bdtFeatures_.push_back(result.getEleContEnergy()[0][0]);
  bdtFeatures_.push_back(result.getEleContEnergy()[1][0]);
  bdtFeatures_.push_back(result.getEleContYMean()[0][0]);
  bdtFeatures_.push_back(result.getEleContEnergy()[0][1]);
  bdtFeatures_.push_back(result.getEleContEnergy()[1][1]);
  bdtFeatures_.push_back(result.getEleContYMean()[0][1]);
  bdtFeatures_.push_back(result.getPhContNHits()[0][0]);
  bdtFeatures_.push_back(result.getPhContYMean()[0][0]);
  bdtFeatures_.push_back(result.getPhContNHits()[0][1]);
  bdtFeatures_.push_back(result.getOutContEnergy()[0][0]);
  bdtFeatures_.push_back(result.getOutContEnergy()[1][0]);
  bdtFeatures_.push_back(result.getOutContEnergy()[2][0]);
  bdtFeatures_.push_back(result.getOutContNHits()[0][0]);
  bdtFeatures_.push_back(result.getOutContXMean()[0][0]);
  bdtFeatures_.push_back(result.getOutContYMean()[0][0]);
  bdtFeatures_.push_back(result.getOutContYMean()[1][0]);
  bdtFeatures_.push_back(result.getOutContYStd()[0][0]);
  bdtFeatures_.push_back(result.getOutContEnergy()[0][1]);
  bdtFeatures_.push_back(result.getOutContEnergy()[1][1]);
  bdtFeatures_.push_back(result.getOutContEnergy()[2][1]);
  bdtFeatures_.push_back(result.getOutContLayerMean()[0][1]);
  bdtFeatures_.push_back(result.getOutContLayerStd()[0][1]);
  bdtFeatures_.push_back(result.getOutContEnergy()[0][2]);
  bdtFeatures_.push_back(result.getOutContLayerMean()[0][2]);
}

References ldmx::EcalVetoResult::getNStraightTracks().

Referenced by produce().

clearProcessor()

void ecal::EcalVetoProcessor::clearProcessor ( )

private

Definition at line 138 of file EcalVetoProcessor.cxx.

                                       {
  cellMap_.clear();
  cellMapTightIso_.clear();
  bdtFeatures_.clear();
 
  nReadoutHits_ = 0;
  summedDet_ = 0;
  summedTightIso_ = 0;
  maxCellDep_ = 0;
  showerRMS_ = 0;
  xStd_ = 0;
  yStd_ = 0;
  avgLayerHit_ = 0;
  stdLayerHit_ = 0;
  deepestLayerHit_ = 0;
  ecalBackEnergy_ = 0;
  // MIP tracking
  nStraightTracks_ = 0;
  nLinregTracks_ = 0;
  firstNearPhLayer_ = 0;
  nNearPhHits_ = 0;
  epAng_ = 0;
  epSep_ = 0;
  epDot_ = 0;
  photonTerritoryHits_ = 0;
 
  std::fill(ecalLayerEdepRaw_.begin(), ecalLayerEdepRaw_.end(), 0);
  std::fill(ecalLayerEdepReadout_.begin(), ecalLayerEdepReadout_.end(), 0);
  std::fill(ecalLayerTime_.begin(), ecalLayerTime_.end(), 0);
}

configure()

void ecal::EcalVetoProcessor::configure ( framework::config::Parameters &  parameters )

overridevirtual

Configure the processor using the given user specified parameters.

Parameters

parameters

Set of parameters used to configure this processor.

Reimplemented from framework::EventProcessor.

Definition at line 88 of file EcalVetoProcessor.cxx.

                                                                       {
  verbose_ = parameters.getParameter<bool>("verbose");
  featureListName_ = parameters.getParameter<std::string>("feature_list_name");
  // Load BDT ONNX file
  rt_ = std::make_unique<ldmx::Ort::ONNXRuntime>(
      parameters.getParameter<std::string>("bdt_file"));
 
  // Read in arrays holding 68% containment radius per layer
  // for different bins in momentum/angle
  rocFileName_ = parameters.getParameter<std::string>("roc_file");
  if (!std::ifstream(rocFileName_).good()) {
    EXCEPTION_RAISE(
        "EcalVetoProcessor",
        "The specified RoC file '" + rocFileName_ + "' does not exist!");
  } else {
    std::ifstream rocfile(rocFileName_);
    std::string line, value;
 
    // Extract the first line in the file
    std::getline(rocfile, line);
    std::vector<double> values;
 
    // Read data, line by line
    while (std::getline(rocfile, line)) {
      std::stringstream ss(line);
      values.clear();
      while (std::getline(ss, value, ',')) {
        double f_value = (value != "") ? std::stof(value) : -1.0;
        values.push_back(f_value);
      }
      roc_range_values_.push_back(values);
    }
  }
 
  nEcalLayers_ = parameters.getParameter<int>("num_ecal_layers");
 
  bdtCutVal_ = parameters.getParameter<double>("disc_cut");
  ecalLayerEdepRaw_.resize(nEcalLayers_, 0);
  ecalLayerEdepReadout_.resize(nEcalLayers_, 0);
  ecalLayerTime_.resize(nEcalLayers_, 0);
 
  beamEnergyMeV_ = parameters.getParameter<double>("beam_energy");
  linreg_radius_ = parameters.getParameter<double>("linreg_radius");
 
  // Set the collection name as defined in the configuration
  collectionName_ = parameters.getParameter<std::string>("collection_name");
  rec_pass_name_ = parameters.getParameter<std::string>("rec_pass_name");
  rec_coll_name_ = parameters.getParameter<std::string>("rec_coll_name");
}

References collectionName_, and framework::config::Parameters::getParameter().

distPtToLine()

float ecal::EcalVetoProcessor::distPtToLine ( TVector3  h1, TVector3  p1, TVector3  p2  )

private

Return the minimum distance between the point h1 and the line passing through points p1 and p2.

Parameters

[in]	h1	Point to find the distance to
[in]	p1	An arbitrary point on the line
[in]	p2	A second, distinct point on the line

Returns

Minimum distance between h1 and the line

Definition at line 1273 of file EcalVetoProcessor.cxx.

                                                                           {
  return ((h1 - p1).Cross(h1 - p2)).Mag() / (p1 - p2).Mag();
}

Referenced by produce().

distTwoLines()

float ecal::EcalVetoProcessor::distTwoLines ( TVector3  v1, TVector3  v2, TVector3  w1, TVector3  w2  )

private

Returns the distance between the lines v and w, with v defined to pass through the points (v1,v2) (and similarly for w).

Parameters

[in]	v1	An arbitrary point on line v
[in]	v2	A second, distinct point on line v
[in]	w1	An arbitrary point on line w
[in]	w2	A second, distinct point on line w

Returns

Closest distance of approach of lines u and v

Definition at line 1260 of file EcalVetoProcessor.cxx.

                                                   {
  TVector3 e1 = v1 - v2;
  TVector3 e2 = w1 - w2;
  TVector3 crs = e1.Cross(e2);
  if (crs.Mag() == 0) {
    return 100.0;  // arbitrary large number; edge case that shouldn't cause
                   // problems.
  } else {
    return std::abs(crs.Dot(v1 - w1) / crs.Mag());
  }
}

Referenced by produce().

fillHitMap()

void ecal::EcalVetoProcessor::fillHitMap ( const std::vector< ldmx::EcalHit > &  ecalRecHits, std::map< ldmx::EcalID, float > &  cellMap_  )

private

Function to load up empty vector of hit maps.

Definition at line 1191 of file EcalVetoProcessor.cxx.

                                        {
  for (const ldmx::EcalHit &hit : ecalRecHits) {
    ldmx::EcalID id(hit.getID());
    cellMap.emplace(id, hit.getEnergy());
  }
}

Referenced by produce().

fillIsolatedHitMap()

void ecal::EcalVetoProcessor::fillIsolatedHitMap ( const std::vector< ldmx::EcalHit > &  ecalRecHits, ldmx::EcalID  globalCentroid, std::map< ldmx::EcalID, float > &  cellMap, std::map< ldmx::EcalID, float > &  cellMapIso, bool  doTight = false  )

private

Definition at line 1200 of file EcalVetoProcessor.cxx.

                                                         {
  for (const ldmx::EcalHit &hit : ecalRecHits) {
    auto isolatedHit = std::make_pair(true, ldmx::EcalID());
    ldmx::EcalID id(hit.getID());
    if (doTight) {
      // Disregard hits that are on the centroid.
      if (id == globalCentroid) continue;
 
      // Skip hits that are on centroid inner ring
      if (geometry_->isNN(globalCentroid, id)) {
        continue;
      }
    }
 
    // Skip hits that have a readout neighbor
    // Get neighboring cell id's and try to look them up in the full cell map
    // (constant speed algo.)
    //  these ideas are only cell/module (must ignore layer)
    std::vector<ldmx::EcalID> cellNbrIds = geometry_->getNN(id);
 
    for (int k = 0; k < cellNbrIds.size(); k++) {
      // update neighbor ID to the current layer
      cellNbrIds[k] = ldmx::EcalID(id.layer(), cellNbrIds[k].module(),
                                   cellNbrIds[k].cell());
      // look in cell hit map to see if it is there
      if (cellMap.find(cellNbrIds[k]) != cellMap.end()) {
        isolatedHit = std::make_pair(false, cellNbrIds[k]);
        break;
      }
    }
    if (!isolatedHit.first) {
      continue;
    }
    // Insert isolated hit
    cellMapIso.emplace(id, hit.getEnergy());
  }
}

GetShowerCentroidIDAndRMS()

ldmx::EcalID ecal::EcalVetoProcessor::GetShowerCentroidIDAndRMS ( const std::vector< ldmx::EcalHit > &  ecalRecHits, double &  showerRMS  )

private

Definition at line 1144 of file EcalVetoProcessor.cxx.

                                                                  {
  auto wgtCentroidCoords = std::make_pair<float, float>(0., 0.);
  float sumEdep = 0;
  ldmx::EcalID returnCellId;
 
  // Calculate Energy Weighted Centroid
  for (const ldmx::EcalHit &hit : ecalRecHits) {
    ldmx::EcalID id(hit.getID());
    CellEnergyPair cell_energy_pair = std::make_pair(id, hit.getEnergy());
    auto [x, y, z] = geometry_->getPosition(id);
    XYCoords centroidCoords = std::make_pair(x, y);
    wgtCentroidCoords.first = wgtCentroidCoords.first +
                              centroidCoords.first * cell_energy_pair.second;
    wgtCentroidCoords.second = wgtCentroidCoords.second +
                               centroidCoords.second * cell_energy_pair.second;
    sumEdep += cell_energy_pair.second;
  }
  wgtCentroidCoords.first = (sumEdep > 1E-6) ? wgtCentroidCoords.first / sumEdep
                                             : wgtCentroidCoords.first;
  wgtCentroidCoords.second = (sumEdep > 1E-6)
                                 ? wgtCentroidCoords.second / sumEdep
                                 : wgtCentroidCoords.second;
  // Find Nearest Cell to Centroid
  float maxDist = 1e6;
  for (const ldmx::EcalHit &hit : ecalRecHits) {
    auto [x, y, z] = geometry_->getPosition(hit.getID());
    XYCoords centroidCoords = std::make_pair(x, y);
 
    float deltaR =
        pow(pow((centroidCoords.first - wgtCentroidCoords.first), 2) +
                pow((centroidCoords.second - wgtCentroidCoords.second), 2),
            .5);
    showerRMS += deltaR * hit.getEnergy();
    if (deltaR < maxDist) {
      maxDist = deltaR;
      returnCellId = ldmx::EcalID(hit.getID());
    }
  }
  if (sumEdep > 0) showerRMS = showerRMS / sumEdep;
  // flatten
  return ldmx::EcalID(0, returnCellId.module(), returnCellId.cell());
}

getTrajectory()

std::vector< std::pair< float, float > > ecal::EcalVetoProcessor::getTrajectory ( std::vector< double >  momentum, std::vector< float >  position  )

private

Definition at line 1243 of file EcalVetoProcessor.cxx.

                                                           {
  std::vector<XYCoords> positions;
  for (int iLayer = 0; iLayer < nEcalLayers_; iLayer++) {
    float posX =
        position[0] + (momentum[0] / momentum[2]) *
                          (geometry_->getZPosition(iLayer) - position[2]);
    float posY =
        position[1] + (momentum[1] / momentum[2]) *
                          (geometry_->getZPosition(iLayer) - position[2]);
    positions.push_back(std::make_pair(posX, posY));
  }
  return positions;
}

produce()

void ecal::EcalVetoProcessor::produce ( framework::Event &  event )

overridevirtual

Process the event and put new data products into it.

Parameters

event

The Event to process.

Implements framework::Producer.

Definition at line 169 of file EcalVetoProcessor.cxx.

                                                     {
  // Get the Ecal Geometry
  geometry_ = &getCondition<ldmx::EcalGeometry>(
      ldmx::EcalGeometry::CONDITIONS_OBJECT_NAME);
 
  ldmx::EcalVetoResult result;
 
  clearProcessor();
 
  // Get the collection of Ecal scoring plane hits. If it doesn't exist,
  // don't bother adding any truth tracking information.
 
  std::vector<double> recoilP;
  std::vector<float> recoilPos;
  std::vector<double> recoilPAtTarget;
  std::vector<float> recoilPosAtTarget;
 
  if (verbose_) {
    ldmx_log(debug) << "   Loop through all of the sim particles and find the "
                       "recoil electron";
  }
 
  if (event.exists("EcalScoringPlaneHits")) {
    //
    // Loop through all of the sim particles and find the recoil electron.
    //
 
    // Get the collection of simulated particles from the event
    auto particleMap{event.getMap<int, ldmx::SimParticle>("SimParticles")};
 
    // Search for the recoil electron
    auto [recoilTrackID, recoilElectron] = Analysis::getRecoil(particleMap);
 
    // Find ECAL SP hit for recoil electron
    auto ecalSpHits{
        event.getCollection<ldmx::SimTrackerHit>("EcalScoringPlaneHits")};
    float pmax = 0;
    for (ldmx::SimTrackerHit &spHit : ecalSpHits) {
      ldmx::SimSpecialID hit_id(spHit.getID());
      if (hit_id.plane() != 31 || spHit.getMomentum()[2] <= 0) continue;
 
      if (spHit.getTrackID() == recoilTrackID) {
        if (sqrt(pow(spHit.getMomentum()[0], 2) +
                 pow(spHit.getMomentum()[1], 2) +
                 pow(spHit.getMomentum()[2], 2)) > pmax) {
          recoilP = spHit.getMomentum();
          recoilPos = spHit.getPosition();
          pmax = sqrt(pow(recoilP[0], 2) + pow(recoilP[1], 2) +
                      pow(recoilP[2], 2));
        }
      }
    }
 
    // Find target SP hit for recoil electron
    if (event.exists("TargetScoringPlaneHits")) {
      std::vector<ldmx::SimTrackerHit> targetSpHits =
          event.getCollection<ldmx::SimTrackerHit>("TargetScoringPlaneHits");
      pmax = 0;
      for (ldmx::SimTrackerHit &spHit : targetSpHits) {
        ldmx::SimSpecialID hit_id(spHit.getID());
        if (hit_id.plane() != 1 || spHit.getMomentum()[2] <= 0) continue;
 
        if (spHit.getTrackID() == recoilTrackID) {
          if (sqrt(pow(spHit.getMomentum()[0], 2) +
                   pow(spHit.getMomentum()[1], 2) +
                   pow(spHit.getMomentum()[2], 2)) > pmax) {
            recoilPAtTarget = spHit.getMomentum();
            recoilPosAtTarget = spHit.getPosition();
            pmax =
                sqrt(pow(recoilPAtTarget[0], 2) + pow(recoilPAtTarget[1], 2) +
                     pow(recoilPAtTarget[2], 2));
          }
        }
      }
    }
  }
 
  if (verbose_) {
    ldmx_log(debug) << "   Get projected trajectories for electron and photon";
  }
  // Get projected trajectories for electron and photon
  std::vector<XYCoords> ele_trajectory, photon_trajectory;
  if (recoilP.size() > 0) {
    ele_trajectory = getTrajectory(recoilP, recoilPos);
    std::vector<double> pvec = recoilPAtTarget.size()
                                   ? recoilPAtTarget
                                   : std::vector<double>{0.0, 0.0, 0.0};
    std::vector<float> posvec = recoilPosAtTarget.size()
                                    ? recoilPosAtTarget
                                    : std::vector<float>{0.0, 0.0, 0.0};
    photon_trajectory =
        getTrajectory({-pvec[0], -pvec[1], beamEnergyMeV_ - pvec[2]}, posvec);
  }
 
  float recoilPMag =
      recoilP.size()
          ? sqrt(pow(recoilP[0], 2) + pow(recoilP[1], 2) + pow(recoilP[2], 2))
          : -1.0;
  float recoilTheta =
      recoilPMag > 0 ? acos(recoilP[2] / recoilPMag) * 180.0 / M_PI : -1.0;
 
  if (verbose_) {
    ldmx_log(debug) << "   Build Radii of containment (ROC)";
  }
  // Use the appropriate containment radii for the recoil electron
  std::vector<double> roc_values_bin0(roc_range_values_[0].begin() + 4,
                                      roc_range_values_[0].end());
  std::vector<double> ele_radii = roc_values_bin0;
  double theta_min, theta_max, p_min, p_max;
  bool inrange;
 
  for (int i = 0; i < roc_range_values_.size(); i++) {
    theta_min = roc_range_values_[i][0];
    theta_max = roc_range_values_[i][1];
    p_min = roc_range_values_[i][2];
    p_max = roc_range_values_[i][3];
    inrange = true;
 
    if (theta_min != -1.0) {
      inrange = inrange && (recoilTheta >= theta_min);
    }
    if (theta_max != -1.0) {
      inrange = inrange && (recoilTheta < theta_max);
    }
    if (p_min != -1.0) {
      inrange = inrange && (recoilPMag >= p_min);
    }
    if (p_max != -1.0) {
      inrange = inrange && (recoilPMag < p_max);
    }
    if (inrange) {
      std::vector<double> roc_values_bini(roc_range_values_[i].begin() + 4,
                                          roc_range_values_[i].end());
      ele_radii = roc_values_bini;
    }
  }
  // Use default RoC bin for photon
  std::vector<double> photon_radii = roc_values_bin0;
 
  // Get the collection of digitized Ecal hits from the event.
  const std::vector<ldmx::EcalHit> ecalRecHits =
      event.getCollection<ldmx::EcalHit>(rec_coll_name_, rec_pass_name_);
 
  ldmx::EcalID globalCentroid =
      GetShowerCentroidIDAndRMS(ecalRecHits, showerRMS_);
  /* ~~ Fill the hit map ~~ O(n)  */
  fillHitMap(ecalRecHits, cellMap_);
  bool doTight = true;
  /* ~~ Fill the isolated hit maps ~~ O(n)  */
  fillIsolatedHitMap(ecalRecHits, globalCentroid, cellMap_, cellMapTightIso_,
                     doTight);
 
  // Loop over the hits from the event to calculate the rest of the important
  // quantities
 
  float wavgLayerHit = 0;
  float xMean = 0;
  float yMean = 0;
 
  // Containment variables
  unsigned int nregions = 5;
  std::vector<float> electronContainmentEnergy(nregions, 0.0);
  std::vector<float> photonContainmentEnergy(nregions, 0.0);
  std::vector<float> outsideContainmentEnergy(nregions, 0.0);
  std::vector<int> outsideContainmentNHits(nregions, 0);
  std::vector<float> outsideContainmentXmean(nregions, 0.0);
  std::vector<float> outsideContainmentYmean(nregions, 0.0);
  std::vector<float> outsideContainmentXstd(nregions, 0.0);
  std::vector<float> outsideContainmentYstd(nregions, 0.0);
  // Longitudinal segmentation
  std::vector<int> segLayers = {0, 6, 17, 34};
  unsigned int nsegments = segLayers.size() - 1;
  std::vector<float> energySeg(nsegments, 0.0);
  std::vector<float> xMeanSeg(nsegments, 0.0);
  std::vector<float> xStdSeg(nsegments, 0.0);
  std::vector<float> yMeanSeg(nsegments, 0.0);
  std::vector<float> yStdSeg(nsegments, 0.0);
  std::vector<float> layerMeanSeg(nsegments, 0.0);
  std::vector<float> layerStdSeg(nsegments, 0.0);
  std::vector<std::vector<float>> eContEnergy(
      nregions, std::vector<float>(nsegments, 0.0));
  std::vector<std::vector<float>> eContXMean(
      nregions, std::vector<float>(nsegments, 0.0));
  std::vector<std::vector<float>> eContYMean(
      nregions, std::vector<float>(nsegments, 0.0));
  std::vector<std::vector<float>> gContEnergy(
      nregions, std::vector<float>(nsegments, 0.0));
  std::vector<std::vector<int>> gContNHits(nregions,
                                           std::vector<int>(nsegments, 0));
  std::vector<std::vector<float>> gContXMean(
      nregions, std::vector<float>(nsegments, 0.0));
  std::vector<std::vector<float>> gContYMean(
      nregions, std::vector<float>(nsegments, 0.0));
  std::vector<std::vector<float>> oContEnergy(
      nregions, std::vector<float>(nsegments, 0.0));
  std::vector<std::vector<int>> oContNHits(nregions,
                                           std::vector<int>(nsegments, 0));
  std::vector<std::vector<float>> oContXMean(
      nregions, std::vector<float>(nsegments, 0.0));
  std::vector<std::vector<float>> oContYMean(
      nregions, std::vector<float>(nsegments, 0.0));
  std::vector<std::vector<float>> oContXStd(nregions,
                                            std::vector<float>(nsegments, 0.0));
  std::vector<std::vector<float>> oContYStd(nregions,
                                            std::vector<float>(nsegments, 0.0));
  std::vector<std::vector<float>> oContLayerMean(
      nregions, std::vector<float>(nsegments, 0.0));
  std::vector<std::vector<float>> oContLayerStd(
      nregions, std::vector<float>(nsegments, 0.0));
 
  // MIP tracking:  vector of hits to be used in the MIP tracking algorithm. All
  // hits inside the electron ROC (or all hits in the ECal if the event is
  // missing an electron) will be included.
  std::vector<HitData> trackingHitList;
 
  if (verbose_) {
    ldmx_log(debug)
        << "   Loop over the hits from the event to calculate the BDT features";
  }
 
  for (const ldmx::EcalHit &hit : ecalRecHits) {
    // Layer-wise quantities
    ldmx::EcalID id(hit.getID());
    ecalLayerEdepRaw_[id.layer()] =
        ecalLayerEdepRaw_[id.layer()] + hit.getEnergy();
    if (id.layer() >= 20) ecalBackEnergy_ += hit.getEnergy();
    if (maxCellDep_ < hit.getEnergy()) maxCellDep_ = hit.getEnergy();
    if (hit.getEnergy() > 0) {
      nReadoutHits_++;
      ecalLayerEdepReadout_[id.layer()] += hit.getEnergy();
      ecalLayerTime_[id.layer()] += (hit.getEnergy()) * hit.getTime();
      auto [x, y, z] = geometry_->getPosition(id);
      xMean += x * hit.getEnergy();
      yMean += y * hit.getEnergy();
      avgLayerHit_ += id.layer();
      wavgLayerHit += id.layer() * hit.getEnergy();
      if (deepestLayerHit_ < id.layer()) {
        deepestLayerHit_ = id.layer();
      }
      XYCoords xy_pair = std::make_pair(x, y);
      float distance_ele_trajectory =
          ele_trajectory.size()
              ? sqrt(
                    pow((xy_pair.first - ele_trajectory[id.layer()].first), 2) +
                    pow((xy_pair.second - ele_trajectory[id.layer()].second),
                        2))
              : -1.0;
      float distance_photon_trajectory =
          photon_trajectory.size()
              ? sqrt(
                    pow((xy_pair.first - photon_trajectory[id.layer()].first),
                        2) +
                    pow((xy_pair.second - photon_trajectory[id.layer()].second),
                        2))
              : -1.0;
 
      // Decide which longitudinal segment the hit is in and add to sums
      for (unsigned int iseg = 0; iseg < nsegments; iseg++) {
        if (id.layer() >= segLayers[iseg] &&
            id.layer() <= segLayers[iseg + 1] - 1) {
          energySeg[iseg] += hit.getEnergy();
          xMeanSeg[iseg] += xy_pair.first * hit.getEnergy();
          yMeanSeg[iseg] += xy_pair.second * hit.getEnergy();
          layerMeanSeg[iseg] += id.layer() * hit.getEnergy();
 
          // Decide which containment region the hit is in and add to sums
          for (unsigned int ireg = 0; ireg < nregions; ireg++) {
            if (distance_ele_trajectory >= ireg * ele_radii[id.layer()] &&
                distance_ele_trajectory < (ireg + 1) * ele_radii[id.layer()]) {
              eContEnergy[ireg][iseg] += hit.getEnergy();
              eContXMean[ireg][iseg] += xy_pair.first * hit.getEnergy();
              eContYMean[ireg][iseg] += xy_pair.second * hit.getEnergy();
            }
            if (distance_photon_trajectory >= ireg * photon_radii[id.layer()] &&
                distance_photon_trajectory <
                    (ireg + 1) * photon_radii[id.layer()]) {
              gContEnergy[ireg][iseg] += hit.getEnergy();
              gContNHits[ireg][iseg] += 1;
              gContXMean[ireg][iseg] += xy_pair.first * hit.getEnergy();
              gContYMean[ireg][iseg] += xy_pair.second * hit.getEnergy();
            }
            if (distance_ele_trajectory > (ireg + 1) * ele_radii[id.layer()] &&
                distance_photon_trajectory >
                    (ireg + 1) * photon_radii[id.layer()]) {
              oContEnergy[ireg][iseg] += hit.getEnergy();
              oContNHits[ireg][iseg] += 1;
              oContXMean[ireg][iseg] += xy_pair.first * hit.getEnergy();
              oContYMean[ireg][iseg] += xy_pair.second * hit.getEnergy();
              oContLayerMean[ireg][iseg] += id.layer() * hit.getEnergy();
            }
          }
        }
      }
 
      // Decide which containment region the hit is in and add to sums
      for (unsigned int ireg = 0; ireg < nregions; ireg++) {
        if (distance_ele_trajectory >= ireg * ele_radii[id.layer()] &&
            distance_ele_trajectory < (ireg + 1) * ele_radii[id.layer()])
          electronContainmentEnergy[ireg] += hit.getEnergy();
        if (distance_photon_trajectory >= ireg * photon_radii[id.layer()] &&
            distance_photon_trajectory < (ireg + 1) * photon_radii[id.layer()])
          photonContainmentEnergy[ireg] += hit.getEnergy();
        if (distance_ele_trajectory > (ireg + 1) * ele_radii[id.layer()] &&
            distance_photon_trajectory >
                (ireg + 1) * photon_radii[id.layer()]) {
          outsideContainmentEnergy[ireg] += hit.getEnergy();
          outsideContainmentNHits[ireg] += 1;
          outsideContainmentXmean[ireg] += xy_pair.first * hit.getEnergy();
          outsideContainmentYmean[ireg] += xy_pair.second * hit.getEnergy();
        }
      }
 
      // MIP tracking:  Decide whether hit should be added to trackingHitList
      // If inside e- RoC or if etraj is missing, use the hit for tracking:
      if (distance_ele_trajectory >= ele_radii[id.layer()] ||
          distance_ele_trajectory == -1.0) {
        HitData hd;
        hd.pos = TVector3(xy_pair.first, xy_pair.second,
                          geometry_->getZPosition(id.layer()));
        hd.layer = id.layer();
        trackingHitList.push_back(hd);
      }
    }
  }  // end loop over rechits
 
  for (const auto &[id, energy] : cellMapTightIso_) {
    if (energy > 0) summedTightIso_ += energy;
  }
 
  for (int iLayer = 0; iLayer < ecalLayerEdepReadout_.size(); iLayer++) {
    ecalLayerTime_[iLayer] =
        ecalLayerTime_[iLayer] / ecalLayerEdepReadout_[iLayer];
    summedDet_ += ecalLayerEdepReadout_[iLayer];
  }
 
  if (nReadoutHits_ > 0) {
    avgLayerHit_ /= nReadoutHits_;
    wavgLayerHit /= summedDet_;
    xMean /= summedDet_;
    yMean /= summedDet_;
  } else {
    wavgLayerHit = 0;
    avgLayerHit_ = 0;
    xMean = 0;
    yMean = 0;
  }
 
  // If necessary, quotient out the total energy from the means
  for (unsigned int iseg = 0; iseg < nsegments; iseg++) {
    if (energySeg[iseg] > 0) {
      xMeanSeg[iseg] /= energySeg[iseg];
      yMeanSeg[iseg] /= energySeg[iseg];
      layerMeanSeg[iseg] /= energySeg[iseg];
    }
    for (unsigned int ireg = 0; ireg < nregions; ireg++) {
      if (eContEnergy[ireg][iseg] > 0) {
        eContXMean[ireg][iseg] /= eContEnergy[ireg][iseg];
        eContYMean[ireg][iseg] /= eContEnergy[ireg][iseg];
      }
      if (gContEnergy[ireg][iseg] > 0) {
        gContXMean[ireg][iseg] /= gContEnergy[ireg][iseg];
        gContYMean[ireg][iseg] /= gContEnergy[ireg][iseg];
      }
      if (oContEnergy[ireg][iseg] > 0) {
        oContXMean[ireg][iseg] /= oContEnergy[ireg][iseg];
        oContYMean[ireg][iseg] /= oContEnergy[ireg][iseg];
        oContLayerMean[ireg][iseg] /= oContEnergy[ireg][iseg];
      }
    }
  }
 
  for (unsigned int ireg = 0; ireg < nregions; ireg++) {
    if (outsideContainmentEnergy[ireg] > 0) {
      outsideContainmentXmean[ireg] /= outsideContainmentEnergy[ireg];
      outsideContainmentYmean[ireg] /= outsideContainmentEnergy[ireg];
    }
  }
 
  // Loop over hits a second time to find the standard deviations.
  for (const ldmx::EcalHit &hit : ecalRecHits) {
    ldmx::EcalID id(hit.getID());
    auto [x, y, z] = geometry_->getPosition(id);
    if (hit.getEnergy() > 0) {
      xStd_ += pow((x - xMean), 2) * hit.getEnergy();
      yStd_ += pow((y - yMean), 2) * hit.getEnergy();
      stdLayerHit_ += pow((id.layer() - wavgLayerHit), 2) * hit.getEnergy();
    }
    XYCoords xy_pair = std::make_pair(x, y);
    float distance_ele_trajectory =
        ele_trajectory.size()
            ? sqrt(pow((xy_pair.first - ele_trajectory[id.layer()].first), 2) +
                   pow((xy_pair.second - ele_trajectory[id.layer()].second), 2))
            : -1.0;
    float distance_photon_trajectory =
        photon_trajectory.size()
            ? sqrt(pow((xy_pair.first - photon_trajectory[id.layer()].first),
                       2) +
                   pow((xy_pair.second - photon_trajectory[id.layer()].second),
                       2))
            : -1.0;
 
    for (unsigned int iseg = 0; iseg < nsegments; iseg++) {
      if (id.layer() >= segLayers[iseg] &&
          id.layer() <= segLayers[iseg + 1] - 1) {
        xStdSeg[iseg] +=
            pow((xy_pair.first - xMeanSeg[iseg]), 2) * hit.getEnergy();
        yStdSeg[iseg] +=
            pow((xy_pair.second - yMeanSeg[iseg]), 2) * hit.getEnergy();
        layerStdSeg[iseg] +=
            pow((id.layer() - layerMeanSeg[iseg]), 2) * hit.getEnergy();
 
        for (unsigned int ireg = 0; ireg < nregions; ireg++) {
          if (distance_ele_trajectory > (ireg + 1) * ele_radii[id.layer()] &&
              distance_photon_trajectory >
                  (ireg + 1) * photon_radii[id.layer()]) {
            oContXStd[ireg][iseg] +=
                pow((xy_pair.first - oContXMean[ireg][iseg]), 2) *
                hit.getEnergy();
            oContYStd[ireg][iseg] +=
                pow((xy_pair.second - oContYMean[ireg][iseg]), 2) *
                hit.getEnergy();
            oContLayerStd[ireg][iseg] +=
                pow((id.layer() - oContLayerMean[ireg][iseg]), 2) *
                hit.getEnergy();
          }
        }
      }
    }
 
    for (unsigned int ireg = 0; ireg < nregions; ireg++) {
      if (distance_ele_trajectory > (ireg + 1) * ele_radii[id.layer()] &&
          distance_photon_trajectory > (ireg + 1) * photon_radii[id.layer()]) {
        outsideContainmentXstd[ireg] +=
            pow((xy_pair.first - outsideContainmentXmean[ireg]), 2) *
            hit.getEnergy();
        outsideContainmentYstd[ireg] +=
            pow((xy_pair.second - outsideContainmentYmean[ireg]), 2) *
            hit.getEnergy();
      }
    }
  }  // end loop over rechits (2nd time)
 
  if (nReadoutHits_ > 0) {
    xStd_ = sqrt(xStd_ / summedDet_);
    yStd_ = sqrt(yStd_ / summedDet_);
    stdLayerHit_ = sqrt(stdLayerHit_ / summedDet_);
  } else {
    xStd_ = 0;
    yStd_ = 0;
    stdLayerHit_ = 0;
  }
 
  // Quotient out the total energies from the standard deviations if possible
  // and take root
  for (unsigned int iseg = 0; iseg < nsegments; iseg++) {
    if (energySeg[iseg] > 0) {
      xStdSeg[iseg] = sqrt(xStdSeg[iseg] / energySeg[iseg]);
      yStdSeg[iseg] = sqrt(yStdSeg[iseg] / energySeg[iseg]);
      layerStdSeg[iseg] = sqrt(layerStdSeg[iseg] / energySeg[iseg]);
    }
    for (unsigned int ireg = 0; ireg < nregions; ireg++) {
      if (oContEnergy[ireg][iseg] > 0) {
        oContXStd[ireg][iseg] =
            sqrt(oContXStd[ireg][iseg] / oContEnergy[ireg][iseg]);
        oContYStd[ireg][iseg] =
            sqrt(oContYStd[ireg][iseg] / oContEnergy[ireg][iseg]);
        oContLayerStd[ireg][iseg] =
            sqrt(oContLayerStd[ireg][iseg] / oContEnergy[ireg][iseg]);
      }
    }
  }
 
  for (unsigned int ireg = 0; ireg < nregions; ireg++) {
    if (outsideContainmentEnergy[ireg] > 0) {
      outsideContainmentXstd[ireg] =
          sqrt(outsideContainmentXstd[ireg] / outsideContainmentEnergy[ireg]);
      outsideContainmentYstd[ireg] =
          sqrt(outsideContainmentYstd[ireg] / outsideContainmentEnergy[ireg]);
    }
  }
 
  if (verbose_) {
    ldmx_log(debug) << "   Find out if the recoil electron is fiducial";
  }
 
  // Find the location of the recoil electron
  // Ecal face is not where the first layer starts,
  // defined in DetDescr/python/EcalGeometry.py
  const float dz_from_face{7.932};
  float drifted_recoil_x{-9999.};
  float drifted_recoil_y{-9999.};
  if (recoilP.size() > 0) {
    drifted_recoil_x =
        (dz_from_face * (recoilP[0] / recoilP[2])) + recoilPos[0];
    drifted_recoil_y =
        (dz_from_face * (recoilP[1] / recoilP[2])) + recoilPos[1];
  }
  const int recoil_layer_index = 0;
 
  // Check if it's fiducial
  bool inside{false};
  // At module level
  const auto ecalID = geometry_->getID(drifted_recoil_x, drifted_recoil_y,
                                       recoil_layer_index, true);
  if (!ecalID.null()) {
    // If fiducial at module level, check at cell level
    const auto cellID =
        geometry_->getID(drifted_recoil_x, drifted_recoil_y, recoil_layer_index,
                         ecalID.getModuleID(), true);
    if (!cellID.null()) {
      inside = true;
    }
  }
 
  if (!inside) {
    ldmx_log(info) << "This event is non-fiducial in ECAL";
  }
 
  // ------------------------------------------------------
  // MIP tracking starts here
 
  /* Goal:  Calculate
   *  nStraightTracks (self-explanatory),
   *  nLinregTracks (tracks found by linreg algorithm),
   */
 
  // Find epAng and epSep, and prepare EP trajectory vectors:
  TVector3 e_traj_start;
  TVector3 e_traj_end;
  TVector3 p_traj_start;
  TVector3 p_traj_end;
  if (ele_trajectory.size() > 0 && photon_trajectory.size() > 0) {
    // Create TVector3s marking the start and endpoints of each projected
    // trajectory
    e_traj_start.SetXYZ(ele_trajectory[0].first, ele_trajectory[0].second,
                        geometry_->getZPosition(0));
    e_traj_end.SetXYZ(ele_trajectory[(nEcalLayers_ - 1)].first,
                      ele_trajectory[(nEcalLayers_ - 1)].second,
                      geometry_->getZPosition((nEcalLayers_ - 1)));
    p_traj_start.SetXYZ(photon_trajectory[0].first, photon_trajectory[0].second,
                        geometry_->getZPosition(0));
    p_traj_end.SetXYZ(photon_trajectory[(nEcalLayers_ - 1)].first,
                      photon_trajectory[(nEcalLayers_ - 1)].second,
                      geometry_->getZPosition((nEcalLayers_ - 1)));
 
    TVector3 evec = e_traj_end - e_traj_start;
    TVector3 e_norm = evec.Unit();
    TVector3 pvec = p_traj_end - p_traj_start;
    TVector3 p_norm = pvec.Unit();
    epDot_ = e_norm.Dot(p_norm);
    epAng_ = acos(epDot_) * 180.0 / M_PI;
    epSep_ = sqrt(pow(e_traj_start.X() - p_traj_start.X(), 2) +
                  pow(e_traj_start.Y() - p_traj_start.Y(), 2));
  } else {
    // Electron trajectory is missing, so all hits in the Ecal are fair game.
    // Pick e/ptraj so that they won't restrict the tracking algorithm (place
    // them far outside the ECal).
    e_traj_start = TVector3(999, 999, geometry_->getZPosition(0));  // 0);
    e_traj_end = TVector3(
        999, 999, geometry_->getZPosition((nEcalLayers_ - 1)));       // 999);
    p_traj_start = TVector3(1000, 1000, geometry_->getZPosition(0));  // 0);
    p_traj_end = TVector3(
        1000, 1000, geometry_->getZPosition((nEcalLayers_ - 1)));  // 1000);
    /*ensures event will not be vetoed by angle/separation cut */
    epAng_ = 999.;
    epSep_ = 999.;
    epDot_ = 999.;
  }
 
  // Near photon step:  Find the first layer of the ECal where a hit near the
  // projected photon trajectory is found Currently unusued pending further
  // study; performance has dropped between v9 and v12. Currently used in
  // segmipBDT
  firstNearPhLayer_ = nEcalLayers_ - 1;
 
  if (photon_trajectory.size() !=
      0) {  // If no photon trajectory, leave this at the default (ECal back)
    for (std::vector<HitData>::iterator it = trackingHitList.begin();
         it != trackingHitList.end(); ++it) {
      float ehDist =
          sqrt(pow((*it).pos.X() - photon_trajectory[(*it).layer].first, 2) +
               pow((*it).pos.Y() - photon_trajectory[(*it).layer].second, 2));
      if (ehDist < 8.7) {
        nNearPhHits_++;
        if ((*it).layer < firstNearPhLayer_) {
          firstNearPhLayer_ = (*it).layer;
        }
      }
    }
  }
 
  // Territories limited to trackingHitList
  TVector3 gToe = (e_traj_start - p_traj_start).Unit();
  TVector3 origin = p_traj_start + 0.5 * 8.7 * gToe;
  if (ele_trajectory.size() > 0) {
    for (auto &hitData : trackingHitList) {
      TVector3 hitPos = hitData.pos;
      TVector3 hitPrime = hitPos - origin;
      if (hitPrime.Dot(gToe) <= 0) {
        photonTerritoryHits_++;
      }
    }
  } else {
    photonTerritoryHits_ = nReadoutHits_;
  }
 
  // ------------------------------------------------------
  // Find straight MIP tracks:
 
  std::sort(trackingHitList.begin(), trackingHitList.end(),
            [](HitData ha, HitData hb) { return ha.layer > hb.layer; });
  // For merging tracks:  Need to keep track of existing tracks
  // Candidate tracks to merge in will always be in front of the current track
  // (lower z), so only store the last hit 3-layer vector:  each track = vector
  // of 3-tuples (xy+layer).
  std::vector<std::vector<HitData>> track_list;
 
  // print trackingHitList
  if (verbose_) {
    ldmx_log(debug) << "====== Tracking hit list (original) length "
                    << trackingHitList.size() << " ======";
    for (int i = 0; i < trackingHitList.size(); i++) {
      std::cout << "[" << trackingHitList[i].pos.X() << ", "
                << trackingHitList[i].pos.Y() << ", "
                << trackingHitList[i].layer << "], ";
    }
    std::cout << std::endl;
    ldmx_log(debug) << "====== END OF Tracking hit list ======";
  }
 
  // in v14 minR is 4.17 mm
  // while maxR is 4.81 mm
  float cellWidth = 2 * geometry_->getCellMaxR();
  for (int iHit = 0; iHit < trackingHitList.size(); iHit++) {
    // list of hit numbers in track (34 = maximum theoretical length)
    int track[34];
    int currenthit{iHit};
    int trackLen{1};
 
    track[0] = iHit;
 
    // Search for hits to add to the track:
    // repeatedly find hits in the front two layers with same x & y positions
    // but since v14 the odd layers are offset, so we allow half a cellWidth
    // deviation and then add to track until no more hits are found
    int jHit = iHit;
    while (jHit < trackingHitList.size()) {
      if ((trackingHitList[jHit].layer ==
               trackingHitList[currenthit].layer - 1 ||
           trackingHitList[jHit].layer ==
               trackingHitList[currenthit].layer - 2) &&
          abs(trackingHitList[jHit].pos.X() -
              trackingHitList[currenthit].pos.X()) <= 0.5 * cellWidth &&
          abs(trackingHitList[jHit].pos.Y() -
              trackingHitList[currenthit].pos.Y()) <= 0.5 * cellWidth) {
        track[trackLen] = jHit;
        trackLen++;
        currenthit = jHit;
      }
      jHit++;
    }
 
    // Confirm that the track is valid:
    if (trackLen < 2) continue;  // Track must contain at least 2 hits
    float closest_e = distTwoLines(trackingHitList[track[0]].pos,
                                   trackingHitList[track[trackLen - 1]].pos,
                                   e_traj_start, e_traj_end);
    float closest_p = distTwoLines(trackingHitList[track[0]].pos,
                                   trackingHitList[track[trackLen - 1]].pos,
                                   p_traj_start, p_traj_end);
    // Make sure that the track is near the photon trajectory and away from the
    // electron trajectory Details of these constraints may be revised
    if (closest_p > cellWidth and closest_e < 2 * cellWidth) continue;
    if (trackLen < 4 and closest_e > closest_p) continue;
    if (verbose_) {
      ldmx_log(debug) << "====== After rejection for MIP tracking ======";
      ldmx_log(debug) << "current hit: [" << trackingHitList[iHit].pos.X()
                      << ", " << trackingHitList[iHit].pos.Y() << ", "
                      << trackingHitList[iHit].layer << "]";
 
      for (int k = 0; k < trackLen; k++) {
        ldmx_log(debug) << "track[" << k << "] position = ["
                        << trackingHitList[track[k]].pos.X() << ", "
                        << trackingHitList[track[k]].pos.Y() << ", "
                        << trackingHitList[track[k]].layer << "]";
      }
    }
 
    // if track found, increment nStraightTracks and remove all hits in track
    // from future consideration
    if (trackLen >= 2) {
      std::vector<HitData> temp_track_list;
      int n_remove = 0;
      for (int kHit = 0; kHit < trackLen; kHit++) {
        temp_track_list.push_back(trackingHitList[track[kHit] - n_remove]);
        trackingHitList.erase(trackingHitList.begin() + track[kHit] - n_remove);
        n_remove++;
      }
      // print trackingHitList
      if (verbose_) {
        ldmx_log(debug) << "====== Tracking hit list (after erase) length "
                        << trackingHitList.size() << " ======";
        for (int i = 0; i < trackingHitList.size(); i++) {
          std::cout << "[" << trackingHitList[i].pos.X() << ", "
                    << trackingHitList[i].pos.Y() << ", "
                    << trackingHitList[i].layer << "] ";
        }
        std::cout << std::endl;
        ldmx_log(debug) << "====== END OF Tracking hit list ======";
      }
 
      track_list.push_back(temp_track_list);
      // The *current* hit will have been removed, so iHit is currently pointing
      // to the next hit. Decrement iHit so no hits will get skipped by iHit++
      iHit--;
    }
  }
 
  ldmx_log(debug) << "Straight tracks found (before merge): "
                  << track_list.size();
  if (verbose_) {
    for (int iTrack = 0; iTrack < track_list.size(); iTrack++) {
      ldmx_log(debug) << "Track " << iTrack << ":";
      for (int iHit = 0; iHit < track_list[iTrack].size(); iHit++) {
        std::cout << "  Hit " << iHit << ": ["
                  << track_list[iTrack][iHit].pos.X() << ", "
                  << track_list[iTrack][iHit].pos.Y() << ", "
                  << track_list[iTrack][iHit].layer << "]" << std::endl;
      }
      std::cout << std::endl;
    }
  }
 
  // Optional addition:  Merge nearby straight tracks.  Not necessary for veto.
  // Criteria:  consider tail of track.  Merge if head of next track is 1/2
  // layers behind, within 1 cell of xy position.
  ldmx_log(debug) << "Beginning track merging using " << track_list.size()
                  << " tracks";
 
  for (int track_i = 0; track_i < track_list.size(); track_i++) {
    // for each track, check the remainder of the track list for compatible
    // tracks
    std::vector<HitData> base_track = track_list[track_i];
    HitData tail_hitdata = base_track.back();  // xylayer of last hit in track
    if (verbose_) ldmx_log(debug) << "  Considering track " << track_i;
    for (int track_j = track_i + 1; track_j < track_list.size(); track_j++) {
      std::vector<HitData> checking_track = track_list[track_j];
      HitData head_hitdata = checking_track.front();
      // if 1-2 layers behind, and xy within one cell...
      if ((head_hitdata.layer == tail_hitdata.layer + 1 ||
           head_hitdata.layer == tail_hitdata.layer + 2) &&
          pow(pow(head_hitdata.pos.X() - tail_hitdata.pos.X(), 2) +
                  pow(head_hitdata.pos.Y() - tail_hitdata.pos.Y(), 2),
              0.5) <= cellWidth) {
        // ...then append the second track to the first one and delete it
        // NOTE:  TO ADD:  (trackingHitList[iHit].pos -
        // trackingHitList[jHit].pos).Mag()
        if (verbose_) {
          ldmx_log(debug) << "     ** Compatible track found at index "
                          << track_j;
          ldmx_log(debug) << "     Tail xylayer: " << head_hitdata.pos.X()
                          << "," << head_hitdata.pos.Y() << ","
                          << head_hitdata.layer;
          ldmx_log(debug) << "     Head xylayer: " << tail_hitdata.pos.X()
                          << "," << tail_hitdata.pos.Y() << ","
                          << tail_hitdata.layer;
        }
        for (int hit_k = 0; hit_k < checking_track.size(); hit_k++) {
          base_track.push_back(track_list[track_j][hit_k]);
        }
        track_list[track_i] = base_track;
        track_list.erase(track_list.begin() + track_j);
        break;
      }
    }
  }
  nStraightTracks_ = track_list.size();
  // print the track list
  ldmx_log(debug) << "Straight tracks found (after merge): "
                  << nStraightTracks_;
  for (int track_i = 0; track_i < track_list.size(); track_i++) {
    ldmx_log(debug) << "Track " << track_i << ":";
    for (int hit_i = 0; hit_i < track_list[track_i].size(); hit_i++) {
      ldmx_log(debug) << "  Hit " << hit_i << ": ["
                      << track_list[track_i][hit_i].pos.X() << ", "
                      << track_list[track_i][hit_i].pos.Y() << ", "
                      << track_list[track_i][hit_i].layer << "]";
    }
  }
 
  // ------------------------------------------------------
  // Linreg tracking:
  ldmx_log(debug) << "Finding linreg tracks from a total of "
                  << trackingHitList.size() << " hits using a radius of "
                  << linreg_radius_ << " mm";
 
  for (int iHit = 0; iHit < trackingHitList.size(); iHit++) {
    // Hits being considered at a given time
    std::vector<int> hitsInRegion;
    TMatrixD Vm(3, 3);
    TMatrixD hdt(3, 3);
    TVector3 slopeVec;
    TVector3 hmean;
    TVector3 hpoint;
    float r_corr_best{0.0};
    // Temp array having 3 potential hits
    int hitNums[3];
    // From the above which are passing the correlation reqs
    int bestHitNums[3];
 
    hitsInRegion.push_back(iHit);
    // Find all hits within 2 cells of the primary hit:
    for (int jHit = 0; jHit < trackingHitList.size(); jHit++) {
      // Dont try to put hits on the same layer to the lin-reg track
      if (trackingHitList[iHit].pos(2) == trackingHitList[jHit].pos(2)) {
        continue;
      }
      float dstToHit =
          (trackingHitList[iHit].pos - trackingHitList[jHit].pos).Mag();
      // This distance optimized to give the best significance
      // it used to be 2*cellWidth, i.e. 4.81 mm
      // note, the layers in the back have a separation of 22.3
      if (dstToHit <= 2 * linreg_radius_) {
        hitsInRegion.push_back(jHit);
      }
    }
    // Found a track that passed the lin-reg reqs
    bool bestLinRegFound{false};
    if (verbose_) {
      ldmx_log(debug) << "There are " << hitsInRegion.size()
                      << " hits within a radius of " << linreg_radius_ << " mm";
    }
    // Look at combinations of hits within the region (do not consider the same
    // combination twice):
    hitNums[0] = iHit;
    for (int jHitInReg = 1; jHitInReg < hitsInRegion.size() - 1; jHitInReg++) {
      // We require (exactly) 3 hits for the lin-reg track building
      if (hitsInRegion.size() < 3) break;
      hitNums[1] = hitsInRegion[jHitInReg];
      for (int kHitReg = jHitInReg + 1; kHitReg < hitsInRegion.size();
           kHitReg++) {
        hitNums[2] = hitsInRegion[kHitReg];
        for (int hInd = 0; hInd < 3; hInd++) {
          // hmean = geometric mean, subtract off from hits to improve SVD
          // performance
          hmean(hInd) = (trackingHitList[hitNums[0]].pos(hInd) +
                         trackingHitList[hitNums[1]].pos(hInd) +
                         trackingHitList[hitNums[2]].pos(hInd)) /
                        3.0;
        }
        for (int hInd = 0; hInd < 3; hInd++) {
          for (int lInd = 0; lInd < 3; lInd++) {
            hdt(hInd, lInd) =
                trackingHitList[hitNums[hInd]].pos(lInd) - hmean(lInd);
          }
        }
 
        // Perform "linreg" on selected points
        // Calculate the determinant of the matrix
        double determinant =
            hdt(0, 0) * (hdt(1, 1) * hdt(2, 2) - hdt(1, 2) * hdt(2, 1)) -
            hdt(0, 1) * (hdt(1, 0) * hdt(2, 2) - hdt(1, 2) * hdt(2, 0)) +
            hdt(0, 2) * (hdt(1, 0) * hdt(2, 1) - hdt(1, 1) * hdt(2, 0));
        // Exit early if the matrix is singular (i.e. det = 0)
        if (determinant == 0) continue;
        // Perform matrix decomposition with SVD
        TDecompSVD svdObj(hdt);
        bool decomposed = svdObj.Decompose();
        if (!decomposed) continue;
 
        // First col of V matrix is the slope of the best-fit line
        Vm = svdObj.GetV();
        for (int hInd = 0; hInd < 3; hInd++) {
          slopeVec(hInd) = Vm[0][hInd];
        }
        // hmean, hpoint are points on the best-fit line
        hpoint = slopeVec + hmean;
        // linreg complete:  Now have best-fit line for 3 hits under
        // consideration Check whether the track is valid:  r^2 must be high,
        // and the track must plausibly originate from the photon
        float closest_e = distTwoLines(hmean, hpoint, e_traj_start, e_traj_end);
        float closest_p = distTwoLines(hmean, hpoint, p_traj_start, p_traj_end);
        // Projected track must be close to the photon; details may change after
        // future study.
        if (closest_p > cellWidth or closest_e < 1.5 * cellWidth) continue;
        // find r^2
        // ~variance
        float vrnc = (trackingHitList[hitNums[0]].pos - hmean).Mag() +
                     (trackingHitList[hitNums[1]].pos - hmean).Mag() +
                     (trackingHitList[hitNums[2]].pos - hmean).Mag();
        // sum of |errors|
        float sumerr =
            distPtToLine(trackingHitList[hitNums[0]].pos, hmean, hpoint) +
            distPtToLine(trackingHitList[hitNums[1]].pos, hmean, hpoint) +
            distPtToLine(trackingHitList[hitNums[2]].pos, hmean, hpoint);
        float r_corr = 1 - sumerr / vrnc;
        // Check whether r^2 exceeds a low minimum r_corr:  "Fake" tracks are
        // still much more common in background, so making the algorithm
        // oversensitive doesn't lower performance significantly
        if (r_corr > r_corr_best and r_corr > .6) {
          r_corr_best = r_corr;
          // Only looking for 3-hit tracks currently
          bestLinRegFound = true;
          for (int k = 0; k < 3; k++) {
            bestHitNums[k] = hitNums[k];
          }
        }
      }  // end loop on hits in the region
    }    // end 2nd loop on hits in the region
 
    // Continue early if not hits on track
    if (!bestLinRegFound) continue;
    // Otherwise increase the number of lin-reg tracks
    nLinregTracks_++;
    ldmx_log(debug) << " Lin-reg track " << nLinregTracks_;
    for (int finalHitIndx = 0; finalHitIndx < 3; finalHitIndx++) {
      ldmx_log(debug) << "   Hit " << finalHitIndx << " ["
                      << trackingHitList[bestHitNums[finalHitIndx]].pos(0)
                      << ", "
                      << trackingHitList[bestHitNums[finalHitIndx]].pos(1)
                      << ", "
                      << trackingHitList[bestHitNums[finalHitIndx]].pos(2)
                      << "] ";
    }
 
    // Exclude all hits in a found track from further consideration:
    for (int lHit = 0; lHit < 3; lHit++) {
      trackingHitList.erase(trackingHitList.begin() + bestHitNums[lHit]);
    }
    iHit--;
  }  // end loop on all hits
 
  ldmx_log(info) << " MIP tracking completed; found " << nStraightTracks_
                 << " straight tracks and " << nLinregTracks_
                 << " lin-reg tracks";
 
  result.setVariables(
      nReadoutHits_, deepestLayerHit_, summedDet_, summedTightIso_, maxCellDep_,
      showerRMS_, xStd_, yStd_, avgLayerHit_, stdLayerHit_, ecalBackEnergy_,
      nStraightTracks_, nLinregTracks_, firstNearPhLayer_, nNearPhHits_,
      photonTerritoryHits_, epAng_, epSep_, epDot_, electronContainmentEnergy,
      photonContainmentEnergy, outsideContainmentEnergy,
      outsideContainmentNHits, outsideContainmentXstd, outsideContainmentYstd,
      energySeg, xMeanSeg, yMeanSeg, xStdSeg, yStdSeg, layerMeanSeg,
      layerStdSeg, eContEnergy, eContXMean, eContYMean, gContEnergy, gContNHits,
      gContXMean, gContYMean, oContEnergy, oContNHits, oContXMean, oContYMean,
      oContXStd, oContYStd, oContLayerMean, oContLayerStd,
      ecalLayerEdepReadout_, recoilP, recoilPos);
 
  buildBDTFeatureVector(result);
  ldmx::Ort::FloatArrays inputs({bdtFeatures_});
  float pred = rt_->run({featureListName_}, inputs, {"probabilities"})[0].at(1);
  // Other considerations were (nLinregTracks_ == 0)  && (firstNearPhLayer_ >=
  // 6)
  // && (epAng_ > 3.0 && epAng_ < 900 || epSep_ > 10.0 && epSep_ < 900)
  bool passesTrackingVeto = (nStraightTracks_ < 3);
  result.setVetoResult(pred > bdtCutVal_ && passesTrackingVeto);
  result.setDiscValue(pred);
  ldmx_log(debug) << " The pred > bdtCutVal = " << (pred > bdtCutVal_);
 
  // Persist in the event if the recoil ele is fiducial
  result.setFiducial(inside);
 
  // If the event passes the veto, keep it. Otherwise,
  // drop the event.
  if (result.passesVeto()) {
    setStorageHint(framework::hint_shouldKeep);
  } else {
    setStorageHint(framework::hint_shouldDrop);
  }
 
  event.add(collectionName_, result);
}

References buildBDTFeatureVector(), collectionName_, distPtToLine(), distTwoLines(), epAng_, epDot_, epSep_, framework::Event::exists(), fillHitMap(), firstNearPhLayer_, geometry_, ldmx::EcalGeometry::getCellMaxR(), ldmx::EcalGeometry::getID(), ldmx::EcalGeometry::getPosition(), ldmx::EcalGeometry::getZPosition(), framework::hint_shouldDrop, framework::hint_shouldKeep, nLinregTracks_, nNearPhHits_, nStraightTracks_, ldmx::EcalVetoResult::passesVeto(), photonTerritoryHits_, ldmx::SimSpecialID::plane(), framework::EventProcessor::setStorageHint(), and ldmx::EcalVetoResult::setVariables().

Member Data Documentation

avgLayerHit_

double ecal::EcalVetoProcessor::avgLayerHit_ {0}

private

Definition at line 125 of file EcalVetoProcessor.h.

{0};

bdtCutVal_

double ecal::EcalVetoProcessor::bdtCutVal_ {0}

private

Definition at line 148 of file EcalVetoProcessor.h.

{0};

bdtFeatures_

std::vector<float> ecal::EcalVetoProcessor::bdtFeatures_

private

Definition at line 157 of file EcalVetoProcessor.h.

bdtFileName_

std::string ecal::EcalVetoProcessor::bdtFileName_

private

Definition at line 155 of file EcalVetoProcessor.h.

beamEnergyMeV_

double ecal::EcalVetoProcessor::beamEnergyMeV_ {0}

private

Definition at line 150 of file EcalVetoProcessor.h.

{0};

cellMap_

std::map<ldmx::EcalID, float> ecal::EcalVetoProcessor::cellMap_

private

Definition at line 106 of file EcalVetoProcessor.h.

cellMapTightIso_

std::map<ldmx::EcalID, float> ecal::EcalVetoProcessor::cellMapTightIso_

private

Definition at line 107 of file EcalVetoProcessor.h.

collectionName_

std::string ecal::EcalVetoProcessor::collectionName_ {"EcalVeto"}

private

Name of the collection which will containt the results.

Definition at line 164 of file EcalVetoProcessor.h.

{"EcalVeto"};

Referenced by configure(), and produce().

deepestLayerHit_

int ecal::EcalVetoProcessor::deepestLayerHit_ {0}

private

Definition at line 117 of file EcalVetoProcessor.h.

{0};

ecalBackEnergy_

double ecal::EcalVetoProcessor::ecalBackEnergy_ {0}

private

Definition at line 127 of file EcalVetoProcessor.h.

{0};

ecalLayerEdepRaw_

std::vector<float> ecal::EcalVetoProcessor::ecalLayerEdepRaw_

private

Definition at line 109 of file EcalVetoProcessor.h.

ecalLayerEdepReadout_

std::vector<float> ecal::EcalVetoProcessor::ecalLayerEdepReadout_

private

Definition at line 110 of file EcalVetoProcessor.h.

ecalLayerTime_

std::vector<float> ecal::EcalVetoProcessor::ecalLayerTime_

private

Definition at line 111 of file EcalVetoProcessor.h.

epAng_

float ecal::EcalVetoProcessor::epAng_ {0}

private

Angular separation between the projected photon and electron trajectories (currently unused)

Definition at line 139 of file EcalVetoProcessor.h.

{0};

Referenced by produce().

epDot_

float ecal::EcalVetoProcessor::epDot_ {0}

private

Dot product of the photon and electron momenta unit vectors.

Definition at line 144 of file EcalVetoProcessor.h.

{0};

Referenced by produce().

epSep_

float ecal::EcalVetoProcessor::epSep_ {0}

private

Distance between the projected photon and electron trajectories at the ECal face.

Definition at line 142 of file EcalVetoProcessor.h.

{0};

Referenced by produce().

featureListName_

std::string ecal::EcalVetoProcessor::featureListName_

private

Definition at line 158 of file EcalVetoProcessor.h.

firstNearPhLayer_

int ecal::EcalVetoProcessor::firstNearPhLayer_ {0}

private

First ECal layer in which a hit is found near the photon.

Definition at line 134 of file EcalVetoProcessor.h.

{0};

Referenced by produce().

geometry_

const ldmx::EcalGeometry* ecal::EcalVetoProcessor::geometry_

private

handle to current geometry (to share with member functions)

Definition at line 169 of file EcalVetoProcessor.h.

Referenced by produce().

linreg_radius_

double ecal::EcalVetoProcessor::linreg_radius_ {0}

private

Definition at line 151 of file EcalVetoProcessor.h.

{0};

maxCellDep_

double ecal::EcalVetoProcessor::maxCellDep_ {0}

private

Definition at line 121 of file EcalVetoProcessor.h.

{0};

nEcalLayers_

int ecal::EcalVetoProcessor::nEcalLayers_ {0}

private

Definition at line 115 of file EcalVetoProcessor.h.

{0};

nLinregTracks_

int ecal::EcalVetoProcessor::nLinregTracks_ {0}

private

Number of "linreg" tracks found in the event.

Definition at line 132 of file EcalVetoProcessor.h.

{0};

Referenced by produce().

nNearPhHits_

int ecal::EcalVetoProcessor::nNearPhHits_ {0}

private

Number of hits near the photon trajectory.

Definition at line 136 of file EcalVetoProcessor.h.

{0};

Referenced by produce().

nReadoutHits_

int ecal::EcalVetoProcessor::nReadoutHits_ {0}

private

Definition at line 116 of file EcalVetoProcessor.h.

{0};

nStraightTracks_

int ecal::EcalVetoProcessor::nStraightTracks_ {0}

private

Number of "straight" tracks found in the event.

Definition at line 130 of file EcalVetoProcessor.h.

{0};

Referenced by produce().

photonTerritoryHits_

int ecal::EcalVetoProcessor::photonTerritoryHits_ {0}

private

Number of hits in the photon territory.

Definition at line 146 of file EcalVetoProcessor.h.

{0};

Referenced by produce().

rec_coll_name_

std::string ecal::EcalVetoProcessor::rec_coll_name_

private

Definition at line 161 of file EcalVetoProcessor.h.

rec_pass_name_

std::string ecal::EcalVetoProcessor::rec_pass_name_

private

Definition at line 160 of file EcalVetoProcessor.h.

roc_range_values_

std::vector<std::vector<double> > ecal::EcalVetoProcessor::roc_range_values_

private

Definition at line 113 of file EcalVetoProcessor.h.

rocFileName_

std::string ecal::EcalVetoProcessor::rocFileName_

private

Definition at line 156 of file EcalVetoProcessor.h.

rt_

std::unique_ptr<ldmx::Ort::ONNXRuntime> ecal::EcalVetoProcessor::rt_

private

Definition at line 166 of file EcalVetoProcessor.h.

showerRMS_

double ecal::EcalVetoProcessor::showerRMS_ {0}

private

Definition at line 122 of file EcalVetoProcessor.h.

{0};

stdLayerHit_

double ecal::EcalVetoProcessor::stdLayerHit_ {0}

private

Definition at line 126 of file EcalVetoProcessor.h.

{0};

summedDet_

double ecal::EcalVetoProcessor::summedDet_ {0}

private

Definition at line 119 of file EcalVetoProcessor.h.

{0};

summedTightIso_

double ecal::EcalVetoProcessor::summedTightIso_ {0}

private

Definition at line 120 of file EcalVetoProcessor.h.

{0};

verbose_

bool ecal::EcalVetoProcessor::verbose_ {false}

private

Definition at line 153 of file EcalVetoProcessor.h.

{false};

xStd_

double ecal::EcalVetoProcessor::xStd_ {0}

private

Definition at line 123 of file EcalVetoProcessor.h.

{0};

yStd_

double ecal::EcalVetoProcessor::yStd_ {0}

private

Definition at line 124 of file EcalVetoProcessor.h.

{0};

---

The documentation for this class was generated from the following files:

• Ecal/include/Ecal/EcalVetoProcessor.h
• Ecal/src/Ecal/EcalVetoProcessor.cxx