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	onNewRun (const ldmx::RunHeader &rh) override
	onNewRun is the first function called for each processor after the conditions are fully configured and accessible.
void	onProcessEnd () override
	Callback for the EventProcessor to take any necessary action when the processing of events finishes, such as calculating job-summary quantities.
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	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.
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, float &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::array< float, 3 > momentum, std::array< float, 3 > 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.
std::vector< float >	trackProp (const ldmx::Tracks &tracks, ldmx::TrackStateType ts_type, const std::string &ts_title)
	Return a vector of parameters for a propagated recoil track.

Private Attributes
int	nevents_ {0}
float	processing_time_ {0.}
std::map< std::string, float >	profiling_map_
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< float > >	roc_range_values_
int	nEcalLayers_ {0}
int	nReadoutHits_ {0}
int	deepestLayerHit_ {0}
float	summedDet_ {0}
float	summedTightIso_ {0}
float	maxCellDep_ {0}
float	showerRMS_ {0}
float	xStd_ {0}
float	yStd_ {0}
float	avgLayerHit_ {0}
float	stdLayerHit_ {0}
float	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 as projected at ECAL.
float	epAngAtTarget_ {0}
	Angular separation between the projected photon and electron trajectories as at Target.
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 at Ecal.
float	epDotAtTarget_ {0}
	Dot product of the photon and electron momenta unit vectors at Target.
int	photonTerritoryHits_ {0}
	Number of hits in the photon territory.
float	bdtCutVal_ {0}
float	beamEnergyMeV_ {0}
bool	run_lin_reg_ {true}
float	linreg_radius_ {0}
std::string	bdtFileName_
std::string	rocFileName_
std::vector< float >	bdtFeatures_
std::string	featureListName_
std::string	sp_pass_name_
std::string	rec_pass_name_
std::string	rec_coll_name_
bool	recoil_from_tracking_
std::string	track_pass_name_
std::string	track_collection_
bool	inverse_skim_ {false}
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 35 of file EcalVetoProcessor.h.

Member Typedef Documentation

CellEnergyPair

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

Definition at line 37 of file EcalVetoProcessor.h.

XYCoords

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

Definition at line 39 of file EcalVetoProcessor.h.

Constructor & Destructor Documentation

EcalVetoProcessor()

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

inline

Definition at line 41 of file EcalVetoProcessor.h.

      : Producer(name, process) {}

~EcalVetoProcessor()

virtual ecal::EcalVetoProcessor::~EcalVetoProcessor ( )

inlinevirtual

Definition at line 44 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 45 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 159 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;
  epAngAtTarget_ = 0;
  epSep_ = 0;
  epDot_ = 0;
  epDotAtTarget_ = 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 103 of file EcalVetoProcessor.cxx.

                                                                       {
  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<float> values;
 
    // Read data, line by line
    while (std::getline(rocfile, line)) {
      std::stringstream ss(line);
      values.clear();
      while (std::getline(ss, value, ',')) {
        float 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");
  run_lin_reg_ = parameters.getParameter<bool>("run_lin_reg");
  linreg_radius_ = parameters.getParameter<double>("linreg_radius");
 
  // Set the collection name as defined in the configuration
  sp_pass_name_ = parameters.getParameter<std::string>("sp_pass_name", "");
  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");
  recoil_from_tracking_ = parameters.getParameter<bool>("recoil_from_tracking");
  track_collection_ = parameters.getParameter<std::string>("track_collection");
  track_pass_name_ =
      parameters.getParameter<std::string>("track_pass_name", "");
  inverse_skim_ = parameters.getParameter<bool>("inverse_skim");
}

References collectionName_.

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 1496 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 1483 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 1414 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 1423 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, float & showerRMS )

private

Definition at line 1367 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::array< float, 3 > momentum, std::array< float, 3 > position )

private

Definition at line 1466 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;
}

onNewRun()

void ecal::EcalVetoProcessor::onNewRun ( const ldmx::RunHeader & rh )

overridevirtual

onNewRun is the first function called for each processor after the conditions are fully configured and accessible.

This is where you could create single-processors, multi-event calculation objects.

Reimplemented from framework::EventProcessor.

Definition at line 31 of file EcalVetoProcessor.cxx.

                                                        {
  profiling_map_["setup"] = 0.;
  profiling_map_["recoil_electron"] = 0.;
  profiling_map_["trajectories"] = 0.;
  profiling_map_["roc_var"] = 0.;
  profiling_map_["fill_hitmaps"] = 0.;
  profiling_map_["containment_var"] = 0.;
  profiling_map_["mip_tracking_setup"] = 0.;
  profiling_map_["straight_tracks"] = 0.;
  profiling_map_["linreg_tracks"] = 0.;
  profiling_map_["set_variables"] = 0.;
  profiling_map_["bdt_variables"] = 0.;
}

onProcessEnd()

void ecal::EcalVetoProcessor::onProcessEnd ( )

overridevirtual

Callback for the EventProcessor to take any necessary action when the processing of events finishes, such as calculating job-summary quantities.

Reimplemented from framework::EventProcessor.

Definition at line 1319 of file EcalVetoProcessor.cxx.

                                     {
  ldmx_log(info) << "AVG Time/Event: " << std::fixed << std::setprecision(2)
                 << processing_time_ / nevents_ << " ms";
 
  ldmx_log(info) << "Breakdown::";
 
  ldmx_log(info) << "setup                  Avg Time/Event = " << std::fixed
                 << std::setprecision(3) << profiling_map_["setup"] / nevents_
                 << " ms";
 
  ldmx_log(info) << "recoil_electron        Avg Time/Event = " << std::fixed
                 << std::setprecision(3)
                 << profiling_map_["recoil_electron"] / nevents_ << " ms";
  ldmx_log(info) << "trajectories           Avg Time/Event = " << std::fixed
                 << std::setprecision(3)
                 << profiling_map_["trajectories"] / nevents_ << " ms";
 
  ldmx_log(info) << "roc_var                Avg Time/Event = " << std::fixed
                 << std::setprecision(3) << profiling_map_["roc_var"] / nevents_
                 << " ms";
 
  ldmx_log(info) << "fill_hitmaps           Avg Time/Event = " << std::fixed
                 << std::setprecision(3)
                 << profiling_map_["fill_hitmaps"] / nevents_ << " ms";
  ldmx_log(info) << "containment_var        Avg Time/Event = " << std::fixed
                 << std::setprecision(3)
                 << profiling_map_["containment_var"] / nevents_ << " ms";
  ldmx_log(info) << "mip_tracking_setup     Avg Time/Event = " << std::fixed
                 << std::setprecision(3)
                 << profiling_map_["mip_tracking_setup"] / nevents_ << " ms";
 
  ldmx_log(info) << "straight_tracks        Avg Time/Event = " << std::fixed
                 << std::setprecision(3)
                 << profiling_map_["straight_tracks"] / nevents_ << " ms";
 
  ldmx_log(info) << "linreg_tracks          Avg Time/Event = " << std::fixed
                 << std::setprecision(3)
                 << profiling_map_["linreg_tracks"] / nevents_ << " ms";
 
  ldmx_log(info) << "set_variables           Avg Time/Event = " << std::fixed
                 << std::setprecision(3)
                 << profiling_map_["set_variables"] / nevents_ << " ms";
 
  ldmx_log(info) << "bdt_variables           Avg Time/Event = " << std::fixed
                 << std::setprecision(3)
                 << profiling_map_["bdt_variables"] / nevents_ << " ms";
}

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 192 of file EcalVetoProcessor.cxx.

                                                     {
  auto start = std::chrono::high_resolution_clock::now();
  nevents_++;
 
  // 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::array<float, 3> recoilP = {0., 0., 0.};
  std::array<float, 3> recoilPos = {-9999., -9999., -9999.};
  std::array<float, 3> recoilPAtTarget = {0., 0., 0.};
  std::array<float, 3> recoilPosAtTarget = {-9999., -9999., -9999.};
 
  auto setup = std::chrono::high_resolution_clock::now();
  profiling_map_["setup"] +=
      std::chrono::duration<float, std::milli>(setup - start).count();
 
  if (!recoil_from_tracking_ &&
      event.exists("EcalScoringPlaneHits", sp_pass_name_)) {
    ldmx_log(trace) << "   Loop through all of the sim particles and find the "
                       "recoil electron";
    //
    // 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", sp_pass_name_)};
    float pmax = 0;
    for (ldmx::SimTrackerHit &spHit : ecalSpHits) {
      ldmx::SimSpecialID hit_id(spHit.getID());
      auto ecal_sp_momentum = spHit.getMomentum();
      auto ecal_sp_position = spHit.getPosition();
      if (hit_id.plane() != 31 || ecal_sp_momentum[2] <= 0) continue;
 
      if (spHit.getTrackID() == recoilTrackID) {
        // A*A is faster than pow(A,2)
        if (sqrt((ecal_sp_momentum[0] * ecal_sp_momentum[0]) +
                 (ecal_sp_momentum[1] * ecal_sp_momentum[1]) +
                 (ecal_sp_momentum[2] * ecal_sp_momentum[2])) > pmax) {
          recoilP = {static_cast<float>(ecal_sp_momentum[0]),
                     static_cast<float>(ecal_sp_momentum[1]),
                     static_cast<float>(ecal_sp_momentum[2])};
          recoilPos = {(ecal_sp_position[0]), (ecal_sp_position[1]),
                       (ecal_sp_position[2])};
          pmax = sqrt(recoilP[0] * recoilP[0] + recoilP[1] * recoilP[1] +
                      recoilP[2] * recoilP[2]);
        }
      }
    }
 
    // Find target SP hit for recoil electron
    if (event.exists("TargetScoringPlaneHits", sp_pass_name_)) {
      std::vector<ldmx::SimTrackerHit> targetSpHits =
          event.getCollection<ldmx::SimTrackerHit>("TargetScoringPlaneHits",
                                                   sp_pass_name_);
      pmax = 0;
      for (ldmx::SimTrackerHit &spHit : targetSpHits) {
        ldmx::SimSpecialID hit_id(spHit.getID());
        auto target_sp_momentum = spHit.getMomentum();
        auto target_sp_position = spHit.getPosition();
        if (hit_id.plane() != 1 || target_sp_momentum[2] <= 0) continue;
 
        if (spHit.getTrackID() == recoilTrackID) {
          if (sqrt((target_sp_momentum[0] * target_sp_momentum[0]) +
                   (target_sp_momentum[1] * target_sp_momentum[1]) +
                   (target_sp_momentum[2] * target_sp_momentum[2])) > pmax) {
            recoilPAtTarget = {static_cast<float>(target_sp_momentum[0]),
                               static_cast<float>(target_sp_momentum[1]),
                               static_cast<float>(target_sp_momentum[2])};
            recoilPosAtTarget = {target_sp_position[0], target_sp_position[1],
                                 target_sp_position[2]};
            // (A*A) is faster than pow(A,2)
            pmax = sqrt((recoilPAtTarget[0] * recoilPAtTarget[0]) +
                        (recoilPAtTarget[1] * recoilPAtTarget[1]) +
                        (recoilPAtTarget[2] * recoilPAtTarget[2]));
          }
        }
      }  // end loop on target SP hits
    }    // end condition on target SP
  }      // end condition on ecal SP
 
  // Get recoilPos using recoil tracking
  bool fiducial_in_tracker{false};
  if (recoil_from_tracking_) {
    ldmx_log(trace) << "  Get recoil tracks collection";
 
    // Get the recoil track collection
    auto recoil_tracks{
        event.getCollection<ldmx::Track>(track_collection_, track_pass_name_)};
 
    ldmx_log(trace) << "  Propagate the recoil ele to the ECAL";
    ldmx::TrackStateType ts_type = ldmx::TrackStateType::AtECAL;
    auto recoil_track_states_ecal = trackProp(recoil_tracks, ts_type, "ecal");
    ldmx_log(trace) << "  Propagate the recoil ele to the Target";
    ldmx::TrackStateType ts_type_target = ldmx::TrackStateType::AtTarget;
    auto recoil_track_states_target =
        trackProp(recoil_tracks, ts_type_target, "target");
 
    ldmx_log(trace) << "  Set recoilPos and recoilP";
    // Redefining recoilPos now to come from the track state
    // track_state_loc0 is recoilPos[0] and track_state_loc1 is recoilPos[1]
    if (!recoil_track_states_ecal.empty()) {
      fiducial_in_tracker = true;
      recoilPos = {recoil_track_states_ecal[0], recoil_track_states_ecal[1],
                   recoil_track_states_ecal[2]};
      recoilP = {(recoil_track_states_ecal[3]), (recoil_track_states_ecal[4]),
                 (recoil_track_states_ecal[5])};
    } else {
      ldmx_log(trace) << "  No recoil track at ECAL";
      fiducial_in_tracker = false;
    }
    ldmx_log(debug) << "  Set recoilP = (" << recoilP[0] << ", " << recoilP[1]
                    << ", " << recoilP[2] << ") and recoilPos = ("
                    << recoilPos[0] << ", " << recoilPos[1] << ", "
                    << recoilPos[2] << ")";
 
    // Repeat the above but now for the taget states
    if (!recoil_track_states_target.empty()) {
      recoilPosAtTarget = {(recoil_track_states_target[0]),
                           (recoil_track_states_target[1]),
                           (recoil_track_states_target[2])};
      recoilPAtTarget = {recoil_track_states_target[3],
                         recoil_track_states_target[4],
                         recoil_track_states_target[5]};
    }
    ldmx_log(debug) << "  Set recoilPAtTarget = (" << recoilPAtTarget[0] << ", "
                    << recoilPAtTarget[1] << ", " << recoilPAtTarget[2]
                    << ") and recoilPosAtTarget = (" << recoilPosAtTarget[0]
                    << ", " << recoilPosAtTarget[1] << ", "
                    << recoilPosAtTarget[2] << ")";
  }  // condition to do recoil information from tracking
 
  ldmx_log(trace) << "   Get projected trajectories for electron and photon";
 
  auto recoil_electron = std::chrono::high_resolution_clock::now();
  profiling_map_["recoil_electron"] +=
      std::chrono::duration<float, std::milli>(recoil_electron - setup).count();
 
  // Get projected trajectories for electron and photon
  std::vector<XYCoords> ele_trajectory, photon_trajectory,
      ele_trajectory_at_target;
  // Require that z-momentum is positive (which will also exclude the default
  // initializaton) Require that the positions are not the default initializaton
  if ((recoilP[2] > 0.) && (recoilPAtTarget[2] > 0.) &&
      (recoilPos[0] != -9999.) && (recoilPosAtTarget[0] != -9999.)) {
    ele_trajectory = getTrajectory(recoilP, recoilPos);
    // Get the photon projection. This does not require that the photon exists
    // tho
    std::array<float, 3> photon_proj_momentum = {
        -recoilPAtTarget[0], -recoilPAtTarget[1],
        beamEnergyMeV_ - recoilPAtTarget[2]};
    photon_trajectory = getTrajectory(photon_proj_momentum, recoilPosAtTarget);
  } else {
    ldmx_log(trace) << "Ele / photon trajectory cannot be determined, pZ = "
                    << recoilP[2] << " pZAtTarget = " << recoilPAtTarget[2]
                    << " X = " << recoilPos[0]
                    << " XAtTarget = " << recoilPosAtTarget[0];
  }
 
  float recoilPMag = (recoilP[2] > 0.) ? sqrt((recoilP[0] * recoilP[0]) +
                                              (recoilP[1] * recoilP[1]) +
                                              (recoilP[2] * recoilP[2]))
                                       : -1.0;
  float recoilTheta =
      recoilPMag > 0 ? acos(recoilP[2] / recoilPMag) * 180.0 / M_PI : -1.0;
 
  ldmx_log(trace) << "   Build Radii of containment (ROC)";
 
  auto trajectories = std::chrono::high_resolution_clock::now();
  profiling_map_["trajectories"] +=
      std::chrono::duration<float, std::milli>(trajectories - recoil_electron)
          .count();
 
  // Use the appropriate containment radii for the recoil electron
  std::vector<float> roc_values_bin0(roc_range_values_[0].begin() + 4,
                                     roc_range_values_[0].end());
  std::vector<float> ele_radii = roc_values_bin0;
  float 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<float> 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<float> photon_radii = roc_values_bin0;
 
  auto roc_var = std::chrono::high_resolution_clock::now();
  profiling_map_["roc_var"] +=
      std::chrono::duration<float, std::milli>(roc_var - trajectories).count();
 
  // 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);
 
  auto fill_hitmaps = std::chrono::high_resolution_clock::now();
  profiling_map_["fill_hitmaps"] +=
      std::chrono::duration<float, std::milli>(fill_hitmaps - roc_var).count();
 
  // 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));
 
  auto containment_var = std::chrono::high_resolution_clock::now();
  profiling_map_["containment_var"] +=
      std::chrono::duration<float, std::milli>(containment_var - fill_hitmaps)
          .count();
 
  // 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;
 
  ldmx_log(trace)
      << "   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) {
      ldmx_log(fatal)
          << "ECal hit has negative or zero energy, this should never happen, "
             "check the threshold settings in HgcrocEmulator";
      continue;
    }
    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]);
    }
  }
 
  ldmx_log(trace) << "   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[2] > 0.) {
    ldmx_log(trace) << "   Recoil electron pX = " << recoilP[0]
                    << " pY = " << recoilP[1] << " pZ = " << recoilP[2];
    ldmx_log(trace) << "   Recoil electron PosX = " << recoilPos[0]
                    << " PosY = " << recoilPos[1] << " PosZ = " << recoilPos[2];
    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_ecal_cell{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_ecal_cell = true;
    }
  }
 
  ldmx_log(info) << "   Is this event is fiducial in ECAL? "
                 << inside_ecal_cell;
 
  // ------------------------------------------------------
  // 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 e_traj_target_start;
  TVector3 e_traj_target_end;
  TVector3 p_traj_start;
  TVector3 p_traj_end;
  if (!ele_trajectory.empty() && !photon_trajectory.empty()) {
    // 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();
 
    // Calculate the angle between the projected electron at Ecal and the photon
    // (at target)
    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));
    // Calculate the electron trajectory with positions and momentum as measured
    // at the target
    if (!ele_trajectory_at_target.empty()) {
      e_traj_target_start.SetXYZ(ele_trajectory_at_target[0].first,
                                 ele_trajectory_at_target[0].second,
                                 geometry_->getZPosition(0));
      e_traj_target_end.SetXYZ(ele_trajectory_at_target[(0)].first,
                               ele_trajectory_at_target[(0)].second,
                               geometry_->getZPosition((nEcalLayers_ - 1)));
      // Now calculate the ep angle at the target
      TVector3 evec_target = e_traj_target_end - e_traj_target_start;
      TVector3 e_norm_target = evec_target.Unit();
      epDotAtTarget_ = e_norm_target.Dot(p_norm);
      epAngAtTarget_ = acos(epDotAtTarget_) * 180.0 / M_PI;
    }
    ldmx_log(trace) << "   Electron trajectory calculated";
  } 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).
    ldmx_log(trace) << "   Electron trajectory is missing";
    e_traj_start = TVector3(999, 999, geometry_->getZPosition(0));
    e_traj_end =
        TVector3(999, 999, geometry_->getZPosition((nEcalLayers_ - 1)));
    p_traj_start = TVector3(1000, 1000, geometry_->getZPosition(0));
    p_traj_end =
        TVector3(1000, 1000, geometry_->getZPosition((nEcalLayers_ - 1)));
    /*ensures event will not be vetoed by angle/separation cut */
    epAng_ = 999.;
    epAngAtTarget_ = 999.;
    epSep_ = 999.;
    epDot_ = 999.;
    epDotAtTarget_ = 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 no photon trajectory, leave this at the default (ECal back)
  if (!photon_trajectory.empty()) {
    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.empty()) {
    for (auto &hitData : trackingHitList) {
      TVector3 hitPos = hitData.pos;
      TVector3 hitPrime = hitPos - origin;
      if (hitPrime.Dot(gToe) <= 0) {
        photonTerritoryHits_++;
      }
    }
  } else {
    photonTerritoryHits_ = nReadoutHits_;
  }
 
  auto mip_tracking_setup = std::chrono::high_resolution_clock::now();
  profiling_map_["mip_tracking_setup"] +=
      std::chrono::duration<float, std::milli>(mip_tracking_setup -
                                               containment_var)
          .count();
 
  // ------------------------------------------------------
  // 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
 
  ldmx_log(trace) << "====== Tracking hit list (original) length "
                  << trackingHitList.size() << " ======";
  for (int i = 0; i < trackingHitList.size(); i++) {
    ldmx_log(trace) << "   [" << trackingHitList[i].pos.X() << ", "
                    << trackingHitList[i].pos.Y() << ", "
                    << trackingHitList[i].layer << "]";
  }
  ldmx_log(trace) << "====== 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;
    ldmx_log(trace) << "====== After rejection for MIP tracking ======";
    ldmx_log(trace) << "current hit: [" << trackingHitList[iHit].pos.X() << ", "
                    << trackingHitList[iHit].pos.Y() << ", "
                    << trackingHitList[iHit].layer << "]";
 
    for (int k = 0; k < trackLen; k++) {
      ldmx_log(trace) << "   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
      ldmx_log(trace) << "====== Tracking hit list (after erase) length "
                      << trackingHitList.size() << " ======";
      for (int i = 0; i < trackingHitList.size(); i++) {
        ldmx_log(trace) << "   [" << trackingHitList[i].pos.X() << ", "
                        << trackingHitList[i].pos.Y() << ", "
                        << trackingHitList[i].layer << "] ";
      }
      ldmx_log(trace) << "====== 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();
  for (int iTrack = 0; iTrack < track_list.size(); iTrack++) {
    ldmx_log(trace) << "Track " << iTrack << ":";
    for (int iHit = 0; iHit < track_list[iTrack].size(); iHit++) {
      ldmx_log(trace) << "  Hit " << iHit << ": ["
                      << track_list[iTrack][iHit].pos.X() << ", "
                      << track_list[iTrack][iHit].pos.Y() << ", "
                      << track_list[iTrack][iHit].layer << "]";
    }
  }
 
  // 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
    ldmx_log(trace) << "  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()
        ldmx_log(trace) << "     ** Compatible track found at index "
                        << track_j;
        ldmx_log(trace) << "     Tail xylayer: " << head_hitdata.pos.X() << ","
                        << head_hitdata.pos.Y() << "," << head_hitdata.layer;
        ldmx_log(trace) << "     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 << "]";
    }
  }
 
  auto straight_tracks = std::chrono::high_resolution_clock::now();
  profiling_map_["straight_tracks"] += std::chrono::duration<float, std::milli>(
                                           straight_tracks - mip_tracking_setup)
                                           .count();
 
  // ------------------------------------------------------
  // Linreg tracking:
  ldmx_log(info) << "Finding linreg tracks from a total of "
                 << trackingHitList.size() << " hits using a radius of "
                 << linreg_radius_ << " mm";
 
  int max_lin_reg_hit{0};
  if (run_lin_reg_) max_lin_reg_hit = trackingHitList.size();
  for (int iHit = 0; iHit < max_lin_reg_hit; 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};
 
    ldmx_log(trace) << "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
        float 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++) {
      if (!trackingHitList.empty()) {
        trackingHitList.erase(trackingHitList.begin() + bestHitNums[lHit]);
      }
    }
    iHit--;
  }  // end loop on all hits
 
  auto linreg_tracks = std::chrono::high_resolution_clock::now();
  profiling_map_["linreg_tracks"] +=
      std::chrono::duration<float, std::milli>(linreg_tracks - straight_tracks)
          .count();
 
  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_, epAngAtTarget_, epSep_, epDot_,
      epDotAtTarget_, 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);
 
  auto set_variables = std::chrono::high_resolution_clock::now();
  profiling_map_["set_variables"] +=
      std::chrono::duration<float, std::milli>(set_variables - linreg_tracks)
          .count();
 
  buildBDTFeatureVector(result);
  ldmx::Ort::FloatArrays inputs({bdtFeatures_});
  float pred = rt_->run({featureListName_}, inputs, {"probabilities"})[0].at(1);
  ldmx_log(info) << " BDT was ran, score is " << pred;
  // 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(info) << " The pred > bdtCutVal = " << (pred > bdtCutVal_)
                 << " and MIP tracking passed = " << passesTrackingVeto;
 
  // Persist in the event if the recoil ele is fiducial
  result.setFiducial(inside_ecal_cell);
  result.setTrackingFiducial(fiducial_in_tracker);
 
  // If the event passes the veto, keep it. Otherwise,
  // drop the event.
  if (!inverse_skim_) {
    if (result.passesVeto()) {
      setStorageHint(framework::hint_shouldKeep);
    } else {
      setStorageHint(framework::hint_shouldDrop);
    }
  } else {
    // Invert the skimming logic
    if (result.passesVeto()) {
      setStorageHint(framework::hint_shouldDrop);
    } else {
      setStorageHint(framework::hint_shouldKeep);
    }
  }
 
  event.add(collectionName_, result);
 
  auto bdt_variables = std::chrono::high_resolution_clock::now();
  profiling_map_["bdt_variables"] +=
      std::chrono::duration<float, std::milli>(bdt_variables - set_variables)
          .count();
 
  auto end = std::chrono::high_resolution_clock::now();
  auto time_diff = end - start;
  processing_time_ +=
      std::chrono::duration<float, std::milli>(time_diff).count();
}

References buildBDTFeatureVector(), collectionName_, distPtToLine(), distTwoLines(), epAng_, epAngAtTarget_, epDot_, epDotAtTarget_, epSep_, framework::Event::exists(), fillHitMap(), firstNearPhLayer_, geometry_, ldmx::EcalGeometry::getCellMaxR(), framework::EventProcessor::getCondition(), 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(), ldmx::EcalVetoResult::setVariables(), and trackProp().

trackProp()

std::vector< float > ecal::EcalVetoProcessor::trackProp ( const ldmx::Tracks & tracks, ldmx::TrackStateType ts_type, const std::string & ts_title )

private

Return a vector of parameters for a propagated recoil track.

Parameters

[in]	tracks	The track collection
[in]	ts_type	The track state type, i.e. tracks state at the ECAL face
[in]	ts_title	The track state title, most likely "ecal"

Returns

Vector of parameters for a propagated recoil track

Definition at line 1500 of file EcalVetoProcessor.cxx.

                                                                           {
  // Vector to hold the new track state variables
  std::vector<float> new_track_states;
 
  // Return if no tracks
  if (tracks.empty()) return new_track_states;
 
  // Otherwise loop on the tracks
  for (auto &track : tracks) {
    // Get track state for ts_type
    auto trk_ts = track.getTrackState(ts_type);
    // Continue if there's no value
    if (!trk_ts.has_value()) continue;
    ldmx::Track::TrackState &ecal_track_state = trk_ts.value();
 
    // Check that the track state is filled
    if (ecal_track_state.params.size() < 5) continue;
 
    float track_state_loc0 = static_cast<float>(ecal_track_state.params[0]);
    float track_state_loc1 = static_cast<float>(ecal_track_state.params[1]);
    // param 2 = phi (azimuthal), param 3 = theta (polar)
    // param 4 = QoP
    // ACTS (local)  to  LDMX (global) coordinates: (y,z,x)->  (x,y,z)
    // convert qop [1/GeV] to p [MeV]
    float p_track_state = (-1 / ecal_track_state.params[4]) * 1000;
    // p * sin(theta) * sin(phi)
    float recoil_mom_x = p_track_state * sin(ecal_track_state.params[3]) *
                         sin(ecal_track_state.params[2]);
    // p * cos(theta)
    float recoil_mom_y = p_track_state * cos(ecal_track_state.params[3]);
    // p * sin(theta) * cos(phi)
    float recoil_mom_z = p_track_state * sin(ecal_track_state.params[3]) *
                         cos(ecal_track_state.params[2]);
 
    // Store the new track state variables
    new_track_states.push_back(track_state_loc0);
    new_track_states.push_back(track_state_loc1);
    // z-position at the ECAL (4) or Target (1)
    if (ts_type == 4) {
      // this should match `ECAL_SCORING_PLANE` in CKFProcessor
      new_track_states.push_back(240.5);
    } else if (ts_type == 1) {
      // This should match `target_surface` in CKFProcessor
      new_track_states.push_back(0.0);
    }
 
    new_track_states.push_back(recoil_mom_x);
    new_track_states.push_back(recoil_mom_y);
    new_track_states.push_back(recoil_mom_z);
 
    // Break after getting the first valid track state
    // TODO: interface this with CLUE to make sure the propageted track
    //       has an associated cluster in the ECAL
    break;
  }
 
  return new_track_states;
}

Referenced by produce().

Member Data Documentation

avgLayerHit_

float ecal::EcalVetoProcessor::avgLayerHit_ {0}

private

Definition at line 156 of file EcalVetoProcessor.h.

{0};

bdtCutVal_

float ecal::EcalVetoProcessor::bdtCutVal_ {0}

private

Definition at line 184 of file EcalVetoProcessor.h.

{0};

bdtFeatures_

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

private

Definition at line 192 of file EcalVetoProcessor.h.

bdtFileName_

std::string ecal::EcalVetoProcessor::bdtFileName_

private

Definition at line 190 of file EcalVetoProcessor.h.

beamEnergyMeV_

float ecal::EcalVetoProcessor::beamEnergyMeV_ {0}

private

Definition at line 186 of file EcalVetoProcessor.h.

{0};

cellMap_

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

private

Definition at line 137 of file EcalVetoProcessor.h.

cellMapTightIso_

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

private

Definition at line 138 of file EcalVetoProcessor.h.

collectionName_

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

private

Name of the collection which will containt the results.

Definition at line 205 of file EcalVetoProcessor.h.

{"EcalVeto"};

Referenced by configure(), and produce().

deepestLayerHit_

int ecal::EcalVetoProcessor::deepestLayerHit_ {0}

private

Definition at line 148 of file EcalVetoProcessor.h.

{0};

ecalBackEnergy_

float ecal::EcalVetoProcessor::ecalBackEnergy_ {0}

private

Definition at line 158 of file EcalVetoProcessor.h.

{0};

ecalLayerEdepRaw_

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

private

Definition at line 140 of file EcalVetoProcessor.h.

ecalLayerEdepReadout_

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

private

Definition at line 141 of file EcalVetoProcessor.h.

ecalLayerTime_

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

private

Definition at line 142 of file EcalVetoProcessor.h.

epAng_

float ecal::EcalVetoProcessor::epAng_ {0}

private

Angular separation between the projected photon and electron trajectories as projected at ECAL.

Definition at line 170 of file EcalVetoProcessor.h.

{0};

Referenced by produce().

epAngAtTarget_

float ecal::EcalVetoProcessor::epAngAtTarget_ {0}

private

Angular separation between the projected photon and electron trajectories as at Target.

Definition at line 173 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 at Ecal.

Definition at line 178 of file EcalVetoProcessor.h.

{0};

Referenced by produce().

epDotAtTarget_

float ecal::EcalVetoProcessor::epDotAtTarget_ {0}

private

Dot product of the photon and electron momenta unit vectors at Target.

Definition at line 180 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 176 of file EcalVetoProcessor.h.

{0};

Referenced by produce().

featureListName_

std::string ecal::EcalVetoProcessor::featureListName_

private

Definition at line 193 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 165 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 210 of file EcalVetoProcessor.h.

Referenced by produce().

inverse_skim_

bool ecal::EcalVetoProcessor::inverse_skim_ {false}

private

Definition at line 202 of file EcalVetoProcessor.h.

{false};

linreg_radius_

float ecal::EcalVetoProcessor::linreg_radius_ {0}

private

Definition at line 188 of file EcalVetoProcessor.h.

{0};

maxCellDep_

float ecal::EcalVetoProcessor::maxCellDep_ {0}

private

Definition at line 152 of file EcalVetoProcessor.h.

{0};

nEcalLayers_

int ecal::EcalVetoProcessor::nEcalLayers_ {0}

private

Definition at line 146 of file EcalVetoProcessor.h.

{0};

nevents_

int ecal::EcalVetoProcessor::nevents_ {0}

private

Definition at line 133 of file EcalVetoProcessor.h.

{0};

nLinregTracks_

int ecal::EcalVetoProcessor::nLinregTracks_ {0}

private

Number of "linreg" tracks found in the event.

Definition at line 163 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 167 of file EcalVetoProcessor.h.

{0};

Referenced by produce().

nReadoutHits_

int ecal::EcalVetoProcessor::nReadoutHits_ {0}

private

Definition at line 147 of file EcalVetoProcessor.h.

{0};

nStraightTracks_

int ecal::EcalVetoProcessor::nStraightTracks_ {0}

private

Number of "straight" tracks found in the event.

Definition at line 161 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 182 of file EcalVetoProcessor.h.

{0};

Referenced by produce().

processing_time_

float ecal::EcalVetoProcessor::processing_time_ {0.}

private

Definition at line 134 of file EcalVetoProcessor.h.

{0.};

profiling_map_

std::map<std::string, float> ecal::EcalVetoProcessor::profiling_map_

private

Definition at line 136 of file EcalVetoProcessor.h.

rec_coll_name_

std::string ecal::EcalVetoProcessor::rec_coll_name_

private

Definition at line 198 of file EcalVetoProcessor.h.

rec_pass_name_

std::string ecal::EcalVetoProcessor::rec_pass_name_

private

Definition at line 197 of file EcalVetoProcessor.h.

recoil_from_tracking_

bool ecal::EcalVetoProcessor::recoil_from_tracking_

private

Definition at line 199 of file EcalVetoProcessor.h.

roc_range_values_

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

private

Definition at line 144 of file EcalVetoProcessor.h.

rocFileName_

std::string ecal::EcalVetoProcessor::rocFileName_

private

Definition at line 191 of file EcalVetoProcessor.h.

rt_

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

private

Definition at line 207 of file EcalVetoProcessor.h.

run_lin_reg_

bool ecal::EcalVetoProcessor::run_lin_reg_ {true}

private

Definition at line 187 of file EcalVetoProcessor.h.

{true};

showerRMS_

float ecal::EcalVetoProcessor::showerRMS_ {0}

private

Definition at line 153 of file EcalVetoProcessor.h.

{0};

sp_pass_name_

std::string ecal::EcalVetoProcessor::sp_pass_name_

private

Definition at line 196 of file EcalVetoProcessor.h.

stdLayerHit_

float ecal::EcalVetoProcessor::stdLayerHit_ {0}

private

Definition at line 157 of file EcalVetoProcessor.h.

{0};

summedDet_

float ecal::EcalVetoProcessor::summedDet_ {0}

private

Definition at line 150 of file EcalVetoProcessor.h.

{0};

summedTightIso_

float ecal::EcalVetoProcessor::summedTightIso_ {0}

private

Definition at line 151 of file EcalVetoProcessor.h.

{0};

track_collection_

std::string ecal::EcalVetoProcessor::track_collection_

private

Definition at line 201 of file EcalVetoProcessor.h.

track_pass_name_

std::string ecal::EcalVetoProcessor::track_pass_name_

private

Definition at line 200 of file EcalVetoProcessor.h.

xStd_

float ecal::EcalVetoProcessor::xStd_ {0}

private

Definition at line 154 of file EcalVetoProcessor.h.

{0};

yStd_

float ecal::EcalVetoProcessor::yStd_ {0}

private

Definition at line 155 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