Use recoil tracking's extrapolated position in EcalVeto (#1495)

* Use recoil tracking to determine projections in EcalVeto
* The feature is off by default
* Use tracking in the ecal pn validation config

.github/validation_samples/ecal_pn/config.py

@@ -194,6 +194,10 @@
 recoil_dqm.truth_hit_collection="RecoilSimHits"
 recoil_dqm.buildHistograms()
 
+# Load ecal veto and use tracking in it
+ecalVeto = ecal_vetos.EcalVetoProcessor()
+ecalVeto.recoil_from_tracking = True
+
 p.sequence.extend([
         digi_tagger,
         digi_recoil,
@@ -204,7 +208,7 @@
         tracking_recoil,
         ecal_digi.EcalDigiProducer(),
         ecal_digi.EcalRecProducer(), 
-        ecal_vetos.EcalVetoProcessor(),
+        ecalVeto,
         hcal_digi.HcalDigiProducer(),
         hcal_digi.HcalRecProducer(),
         *ts_digis,

Ecal/include/Ecal/EcalVetoProcessor.h

@@ -19,6 +19,9 @@
 // ROOT (MIP tracking)
 #include "TVector3.h"
 
+// For recoil tracking
+#include "Tracking/Event/Track.h"
+
 // C++
 #include <map>
 #include <memory>
@@ -102,6 +105,17 @@ class EcalVetoProcessor : public framework::Producer {
    */
   float distPtToLine(TVector3 h1, TVector3 p1, TVector3 p2);
 
+  /**
+   * Return a vector of parameters for a propagated recoil track
+   * @param[in] tracks The track collection
+   * @param[in] ts_type The track state type, i.e. tracks state at the ECAL face
+   * @param[in] ts_title The track state title, most likely "ecal"
+   * @returns Vector of parameters for a propagated recoil track
+   */
+  std::vector<float> trackProp(const ldmx::Tracks& tracks,
+                               ldmx::TrackStateType ts_type,
+                               const std::string& ts_title);
+
  private:
   std::map<ldmx::EcalID, float> cellMap_;
   std::map<ldmx::EcalID, float> cellMapTightIso_;
@@ -159,6 +173,8 @@ class EcalVetoProcessor : public framework::Producer {
 
   std::string rec_pass_name_;
   std::string rec_coll_name_;
+  bool recoil_from_tracking_;
+  std::string track_collection_;
 
   /** Name of the collection which will containt the results. */
   std::string collectionName_{"EcalVeto"};

Ecal/python/vetos.py

@@ -26,6 +26,8 @@ def __init__(self,name = 'ecalVeto') :
         self.collection_name = "EcalVeto"
         self.rec_coll_name = 'EcalRecHits'
         self.rec_pass_name = ''
+        self.recoil_from_tracking = False # Will be True soon
+        self.track_collection = 'RecoilTracks'
 
 
 class DNNEcalVetoProcessor(ldmxcfg.Producer) :

Ecal/src/Ecal/EcalVetoProcessor.cxx

@@ -133,6 +133,8 @@ void EcalVetoProcessor::configure(framework::config::Parameters &parameters) {
   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");
 }
 
 void EcalVetoProcessor::clearProcessor() {
@@ -242,13 +244,28 @@ void EcalVetoProcessor::produce(framework::Event &event) {
       }
     }
   }
+  // Get recoilPos using recoil tracking
+  if (recoil_from_tracking_) {
+    std::vector<float> recoil_track_states;
+    if (verbose_) {
+      ldmx_log(debug) << "   Propagate recoil tracks to ECAL face";
+    }
+    // Get the recoil track collection
+    auto recoil_tracks{event.getCollection<ldmx::Track>(track_collection_)};
+
+    ldmx::TrackStateType ts_type = ldmx::TrackStateType::AtECAL;
+    recoil_track_states = trackProp(recoil_tracks, ts_type, "ecal");
+    // Redefining recoilPos now to come from the track state
+    // track_state_loc0 is recoilPos[0] and track_state_loc1 is recoilPos[1]
+    recoilPos = recoil_track_states;
+  }
 
   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) {
+  if (!recoilP.empty() && !recoilPos.empty()) {
     ele_trajectory = getTrajectory(recoilP, recoilPos);
     std::vector<double> pvec = recoilPAtTarget.size()
                                    ? recoilPAtTarget
@@ -657,7 +674,7 @@ void EcalVetoProcessor::produce(framework::Event &event) {
   const float dz_from_face{7.932};
   float drifted_recoil_x{-9999.};
   float drifted_recoil_y{-9999.};
-  if (recoilP.size() > 0) {
+  if (!recoilP.empty()) {
     drifted_recoil_x =
         (dz_from_face * (recoilP[0] / recoilP[2])) + recoilPos[0];
     drifted_recoil_y =
@@ -697,7 +714,7 @@ void EcalVetoProcessor::produce(framework::Event &event) {
   TVector3 e_traj_end;
   TVector3 p_traj_start;
   TVector3 p_traj_end;
-  if (ele_trajectory.size() > 0 && photon_trajectory.size() > 0) {
+  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,
@@ -741,8 +758,8 @@ void EcalVetoProcessor::produce(framework::Event &event) {
   // segmipBDT
   firstNearPhLayer_ = nEcalLayers_ - 1;
 
-  if (photon_trajectory.size() !=
-      0) {  // If no photon trajectory, leave this at the default (ECal back)
+  // 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 =
@@ -760,7 +777,7 @@ void EcalVetoProcessor::produce(framework::Event &event) {
   // 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) {
+  if (!ele_trajectory.empty()) {
     for (auto &hitData : trackingHitList) {
       TVector3 hitPos = hitData.pos;
       TVector3 hitPrime = hitPos - origin;
@@ -1118,6 +1135,7 @@ void EcalVetoProcessor::produce(framework::Event &event) {
   buildBDTFeatureVector(result);
   ldmx::Ort::FloatArrays inputs({bdtFeatures_});
   float pred = rt_->run({featureListName_}, inputs, {"probabilities"})[0].at(1);
+  ldmx_log(debug) << "  BDT was ran, score is " << pred;
   // Other considerations were (nLinregTracks_ == 0)  && (firstNearPhLayer_ >=
   // 6)
   // && (epAng_ > 3.0 && epAng_ < 900 || epSep_ > 10.0 && epSep_ < 900)
@@ -1274,6 +1292,44 @@ float EcalVetoProcessor::distPtToLine(TVector3 h1, TVector3 p1, TVector3 p2) {
   return ((h1 - p1).Cross(h1 - p2)).Mag() / (p1 - p2).Mag();
 }
 
+std::vector<float> EcalVetoProcessor::trackProp(const ldmx::Tracks &tracks,
+                                                ldmx::TrackStateType ts_type,
+                                                const std::string &ts_title) {
+  // 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]);
+
+    // 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 scoring plane
+    new_track_states.push_back(239.999);
+
+    // 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;
+}
+
 }  // namespace ecal
 
 DECLARE_PRODUCER_NS(ecal, EcalVetoProcessor);

Ecal/src/Ecal/Event/EcalVetoResult.cxx

@@ -170,6 +170,8 @@ void EcalVetoResult::setVariables(
     recoilPx_ = recoilP[0];
     recoilPy_ = recoilP[1];
     recoilPz_ = recoilP[2];
+  }
+  if (!recoilPos.empty()) {
     recoilX_ = recoilPos[0];
     recoilY_ = recoilPos[1];
   }