Implement dEdx mass estimator to do PID with recoil tracker (#1510)

* Create event class TrackDeDxMassEstimate
* Mass estimator with SimHits
* Create DQM module for dEdx mass estimation

---------

Co-authored-by: Tamas Vami <Tamas.Almos.Vami@cern.ch>
Co-authored-by: github-actions[bot] <github-actions[bot]@users.noreply.github.com>

.github/validation_samples/ecal_pn/config.py

@@ -185,6 +185,7 @@
 seed_recoil_dqm.title = ""
 seed_recoil_dqm.buildHistograms()
 
+
 recoil_dqm = tkdqm.TrackingRecoDQM("RecoilDQM")
 recoil_dqm.track_collection = tracking_recoil.out_trk_collection
 recoil_dqm.truth_collection = "RecoilTruthTracks"

DQM/include/DQM/TrkDeDxMassEstFeatures.h

@@ -0,0 +1,46 @@
+#ifndef DQM_TRKDEDXMASSESTFEATURES_H
+#define DQM_TRKDEDXMASSESTFEATURES_H
+
+// LDMX Framework
+#include "Framework/Configure/Parameters.h"  // Needed to import parameters from configuration file
+#include "Framework/EventProcessor.h"  //Needed to declare processor
+
+namespace dqm {
+
+/**
+ * @class TrkDeDxMassEstFeatures
+ * @brief Generate histograms to check tracker dE/dx mass estimate features
+ */
+class TrkDeDxMassEstFeatures : public framework::Analyzer {
+ public:
+  /**
+   * Constructor
+   *
+   * Blank Analyzer constructor
+   */
+  TrkDeDxMassEstFeatures(const std::string& name, framework::Process& process)
+      : framework::Analyzer(name, process) {}
+
+  /**
+   * Input python configuration parameters
+   */
+  virtual void configure(framework::config::Parameters& ps);
+
+  /**
+   * Fills histograms
+   */
+  virtual void analyze(const framework::Event& event);
+
+  /// Method executed before processing of events begins.
+  void onProcessStart() override;
+
+ private:
+  /// Collection Name for mass estimate object
+  std::string mass_estimate_name_;
+
+  /// Pass Name for mass estimate object
+  std::string mass_estimate_pass_;
+};
+}  // namespace dqm
+
+#endif /* DQM_TRKDEDXMASSESTFEATURES_H */

DQM/python/dqm.py

@@ -552,6 +552,30 @@ def __init__(self,name='RecoilTracker') :
             self.build1DHistogram("tpy%s" % t, "Recoil e^{-} Truth p_{y} (MeV)", 100, -10, 10)
             self.build1DHistogram("tpz%s" % t, "Recoil e^{-} Truth p_{z} (MeV)", 260, -100, 2500)
 
+class TrkDeDxMassEstFeatures(ldmxcfg.Analyzer) :
+    """Configured TrkDeDxMassEstFeatures python object
+    
+    Contains an instance of TrkDeDxMassEstFeatures that
+    has already been configured.
+    
+    Builds the necessary histograms as well.
+    
+    Examples
+    --------
+        from LDMX.DQM import dqm
+        p.sequence.append( dqm.TrkDeDxMassEstFeatures() )
+    """
+
+    def __init__(self,name='TrkDeDxMassEstFeatures') :
+        super().__init__(name, "dqm::TrkDeDxMassEstFeatures",'DQM')
+
+        self.mass_estimate_name = "TrackDeDxMassEstimate"
+        self.mass_estimate_pass = ""
+
+        self.build1DHistogram("mass_estimate", "Mass Estimate [MeV]", 100, 0, 2000)
+        self.build1DHistogram("track_type", "Track Type", 3, 0, 3)
+
+
 class TrigScintSimDQM(ldmxcfg.Analyzer) :
     """Configured TrigScintSimDQM python object
     
@@ -714,6 +738,9 @@ def __init__(self, name='SampleValidation') :
         RecoilTrackerDQM()
         ]
 
+dEdx_dqm = [
+        TrkDeDxMassEstFeatures()
+        ]
 
 trigScint_dqm = [
     TrigScintSimDQM('TrigScintSimPad1','TriggerPad1SimHits','pad1'),

DQM/src/DQM/TrkDeDxMassEstFeatures.cxx

@@ -0,0 +1,40 @@
+
+#include "DQM/TrkDeDxMassEstFeatures.h"
+
+#include "Recon/Event/TrackDeDxMassEstimate.h"
+
+namespace dqm {
+
+void TrkDeDxMassEstFeatures::configure(framework::config::Parameters &ps) {
+  mass_estimate_name_ = ps.getParameter<std::string>("mass_estimate_name");
+  mass_estimate_pass_ = ps.getParameter<std::string>("mass_estimate_pass");
+
+  return;
+}
+
+void TrkDeDxMassEstFeatures::analyze(const framework::Event &event) {
+  auto mass_estimates_{event.getCollection<ldmx::TrackDeDxMassEstimate>(
+      mass_estimate_name_, mass_estimate_pass_)};
+
+  for (const auto &mass_est_ : mass_estimates_) {
+    histograms_.fill("mass_estimate", mass_est_.getMass());
+    histograms_.fill("track_type", mass_est_.getTrackType());
+  }
+
+  return;
+}
+
+void TrkDeDxMassEstFeatures::onProcessStart() {
+  std::vector<std::string> labels = {"Other",   // 0
+                                     "Tagger",  // 1
+                                     "Recoil",  // 2
+                                     ""};
+  TH1 *hist = histograms_.get("track_type");
+  for (int ilabel{1}; ilabel < labels.size(); ++ilabel) {
+    hist->GetXaxis()->SetBinLabel(ilabel, labels[ilabel - 1].c_str());
+  }
+}
+
+}  // namespace dqm
+
+DECLARE_ANALYZER_NS(dqm, TrkDeDxMassEstFeatures);

Recon/CMakeLists.txt

@@ -24,17 +24,19 @@ if(BUILD_EVENT_ONLY)
                          class "TriggerResult" )
   register_event_object( module_path "Recon/Event" namespace "ldmx" 
                          class "HgcrocDigiCollection" )
-                       register_event_object( module_path "Recon/Event" 
+  register_event_object( module_path "Recon/Event" 
                          namespace "ldmx" class "HgcrocTrigDigi" 
                          type "collection" )
-                       register_event_object( module_path "Recon/Event" 
+  register_event_object( module_path "Recon/Event" 
                          namespace "ldmx" class "CaloTrigPrim" 
                          type "collection" )
-		       register_event_object(module_path "Recon/Event" namespace "ldmx" 
+  register_event_object(module_path "Recon/Event" namespace "ldmx" 
                         class "PFCandidate" type "collection" )
-		      register_event_object( module_path "Recon/Event" namespace "ldmx"
+  register_event_object( module_path "Recon/Event" namespace "ldmx"
   			class "BeamElectronTruth" type "collection" )
-   setup_library(module Recon name Event
+  register_event_object( module_path "Recon/Event" namespace "ldmx" 
+                         class "TrackDeDxMassEstimate" type "collection" )
+ setup_library(module Recon name Event
                  dependencies ROOT::Core
                  register_target)
 

Recon/include/Recon/Event/TrackDeDxMassEstimate.h

@@ -0,0 +1,107 @@
+/**
+ * @file TrackDeDxMassEstimate.h
+ * @brief Class that represents the estimated mass of a particle
+ * using tracker dE/dx information
+ * @author Danyi Zhang, Tamas Almos Vami (UCSB)
+ */
+
+#ifndef RECON_TRACKDEDXMASSESTIMATE_H_
+#define RECON_TRACKDEDXMASSESTIMATE_H_
+
+#include <iostream>
+
+// ROOT
+#include "TObject.h"  //For ClassDef
+
+// LDMX
+// #include "Event/Track.h"
+// #include "Event/Measurement.h"
+
+namespace ldmx {
+/**
+ * @class TrackDeDxMassEstimate
+ * @brief Represents the estimated mass of a particle
+ * using tracker dE/dx information
+ * @note This class represents the estimated mass information
+ * from a tracker including mass, track index, and the track type
+ */
+
+class TrackDeDxMassEstimate {
+ public:
+  /**
+   * Class constructor.
+   */
+  TrackDeDxMassEstimate();
+
+  /**
+   * Class destructor.
+   */
+  virtual ~TrackDeDxMassEstimate();
+
+  /**
+   * Clear the data in the object.
+   */
+  void Clear();
+
+  /**
+   * Print out the object.
+   */
+  void Print() const;
+
+  /**
+   * Set the estimated mass of the particle/track.
+   * @param mass The estimated mass of the particle/track.
+   */
+  void setMass(double mass) { mass_ = mass; }
+
+  /**
+   * Set the index of the track.
+   * @param track_index The index of the track.
+   */
+  void setTrackIndex(int track_index) { track_index_ = track_index; }
+
+  /**
+   * Set the type of the track.
+   * @param track_type The type of the track.
+   * 1: tagger track, 2: recoil track
+   * Possibly consider truth with 0
+   * and ECAL with 3
+   */
+  void setTrackType(int track_type) { track_type_ = track_type; }
+
+  /**
+   * Get the estimated mass of the particle/track.
+   * @return The estimated mass of the particle/track.
+   */
+  double getMass() const { return mass_; }
+
+  /**
+   * Get the index of the track.
+   * @return The index of the track.
+   */
+  int getTrackIndex() const { return track_index_; }
+
+  /**
+   * Get the type of the track.
+   * @return The type of the track.
+   */
+  int getTrackType() const { return track_type_; }
+
+ private:
+  /* The estimated mass of the particle/track */
+  double mass_{0.};
+
+  /* The index of the track */
+  int track_index_{-1};
+
+  /* The type of the track */
+  int track_type_{-1};
+
+  /**
+   * The ROOT class definition.
+   */
+  ClassDef(TrackDeDxMassEstimate, 1);
+};
+}  // namespace ldmx
+
+#endif  // RECON_TRACKDEDXMASSESTIMATE_H_

Recon/include/Recon/TrackDeDxMassEstimator.h

@@ -0,0 +1,50 @@
+/**
+ * @file TrackDeDxMassEstimator.h
+ * @brief Class that estimates the mass of a particle using tracker dE/dx
+ * information
+ * @author Danyi Zhang, Tamas Almos Vami (UCSB)
+ */
+
+#ifndef RECON_TRACKDEDXMASSESTIMATOR_H_
+#define RECON_TRACKDEDXMASSESTIMATOR_H_
+
+// LDMX Framework
+#include "Framework/Configure/Parameters.h"
+#include "Framework/EventProcessor.h"
+#include "Recon/Event/TrackDeDxMassEstimate.h"
+#include "Tracking/Event/Measurement.h"
+#include "Tracking/Event/Track.h"
+
+namespace recon {
+
+/**
+ * @class TrackDeDxMassEstimator
+ * @brief
+ */
+class TrackDeDxMassEstimator : public framework::Producer {
+ public:
+  TrackDeDxMassEstimator(const std::string& name, framework::Process& process)
+      : framework::Producer(name, process) {}
+
+  virtual void configure(framework::config::Parameters& ps) override;
+
+  virtual void produce(framework::Event& event) override;
+
+ private:
+  // specific verbosity of this producer
+  int verbose_{0};
+
+  float fit_res_C_{0.};
+  float fit_res_K_{-9999.};
+
+  // name of input track collection
+  std::string track_collection_;
+
+  // name of input simhit collection
+  std::string simhit_collection_;
+
+};  // TrackDeDxMassEstimator
+
+}  // namespace recon
+
+#endif  // RECON_TRACKDEDXMASSESTIMATOR_H_

Recon/python/trackDeDxMassEstimator.py

@@ -0,0 +1,40 @@
+"""Configuration for TrackDeDxMassEstimator
+
+The tracker dE/dx vs momentum 2d distribution profile 
+histogram can be fitted using an approximated 
+Bethe-Bloch parametrization at low relativistic regime:
+    dE/dx = K * m^2 / p^2 + C,
+where m is the particle mass, p is the particle momentum,
+K and C are fit parameters.
+Using the fitted result of the parameters, the mass of 
+the particle can be calculated from its dE/dx and momentum.
+
+Attributes:
+------------- 
+track_collection : string
+    Name of the track collection used as input
+fit_res_C : float
+    The fitted result of the constant term C (unit: MeV/mm)
+fit_res_K : float
+    The fitted result of the factor K of the quadratic term (dimentionless)
+
+Examples
+--------
+    from LDMX.Recon.trackDeDxMassEstimator import reoilTrackMassEstimator
+    p.sequence.append( reoilTrackMassEstimator )
+"""
+
+from LDMX.Framework import ldmxcfg
+
+class trackDeDxMassEstimator(ldmxcfg.Producer) :
+    """Configuration for the mass estimator from tracker dEdx"""
+
+    def __init__(self, name="TrackDeDxMassEstimator") :
+        super().__init__(name,'recon::TrackDeDxMassEstimator','Recon')
+
+        self.track_collection = "RecoilTruthTracks"
+        self.fit_res_C = 3.094
+        self.fit_res_K = 1.862
+
+recoilTrackMassEstimator = trackDeDxMassEstimator("RecoilTrackMassEstimator")
+recoilTrackMassEstimator.track_collection = "RecoilTruthTracks"

Recon/src/Recon/Event/TrackDeDxMassEstimate.cxx

@@ -0,0 +1,23 @@
+#include "Recon/Event/TrackDeDxMassEstimate.h"
+
+ClassImp(ldmx::TrackDeDxMassEstimate);
+
+namespace ldmx {
+TrackDeDxMassEstimate::TrackDeDxMassEstimate() {}
+
+TrackDeDxMassEstimate::~TrackDeDxMassEstimate() { Clear(); }
+
+void TrackDeDxMassEstimate::Clear() {
+  mass_ = 0.;
+  track_index_ = -1;
+  track_type_ = -1;
+}
+
+void TrackDeDxMassEstimate::Print() const {
+  std::cout << "TrackDeDxMassEstimate { "
+            << "Mass: " << mass_ << ", "
+            << "Track Index: " << track_index_ << ", "
+            << "Track Type: " << track_type_ << " }" << std::endl;
+}
+
+}  // namespace ldmx