added originIDs to cluster collection

Ecal/include/Ecal/EcalClusterProducer.h

@@ -70,6 +70,9 @@ class EcalClusterProducer : public framework::Producer {
 
   std::string recHitCollName_;
   std::string recHitPassName_;
+  std::string simHitCollName_;
+  std::string simHitPassName_;
+
   std::string algoCollName_;
   std::string clusterCollName_;
 

Ecal/include/Ecal/Event/EcalCluster.h

@@ -10,6 +10,8 @@
 // ldmx-sw
 #include "Ecal/Event/EcalHit.h"
 #include "Recon/Event/CaloCluster.h"
+#include "SimCore/Event/SimCalorimeterHit.h"
+
 
 namespace ldmx {
 
@@ -40,6 +42,8 @@ class EcalCluster : public ldmx::CaloCluster {
 
   void addFirstLayerHits(const std::vector<ldmx::EcalHit> hitsVec);
 
+  void findHitOrigins(const std::vector<ldmx::SimCalorimeterHit>& ecalSimHits);
+
   bool operator<(const EcalCluster& rhs) const {
     return this->getEnergy() < rhs.getEnergy();
   }
@@ -62,6 +66,8 @@ class EcalCluster : public ldmx::CaloCluster {
   // Could add further ECal-specific info here...
 
   std::vector<unsigned int> firstLayerHitIDs_;
+  std::vector<unsigned int> hitOriginIDs_;
+
 
   double firstLayerCentroidX_{0};
   double firstLayerCentroidY_{0};

Ecal/python/ecalClusters.py

@@ -17,6 +17,9 @@ def __init__(self,name='ecalClusters') :
         self.recHitCollName = 'EcalRecHits'
         self.recHitPassName = ''
 
+        self.simHitCollName = 'EcalSimHits'
+        self.simHitPassName = ''
+
         # Name of the cluster collection to make
         self.clusterCollName = "ecalClusters"
 

Ecal/src/Ecal/EcalClusterProducer.cxx

@@ -27,6 +27,10 @@ void EcalClusterProducer::configure(framework::config::Parameters& parameters) {
 
   recHitCollName_ = parameters.getParameter<std::string>("recHitCollName");
   recHitPassName_ = parameters.getParameter<std::string>("recHitPassName");
+
+  simHitCollName_ = parameters.getParameter<std::string>("simHitCollName");
+  simHitPassName_ = parameters.getParameter<std::string>("simHitPassName");
+
   algoCollName_ = parameters.getParameter<std::string>("algoCollName");
   algoName_ = parameters.getParameter<std::string>("algoName");
   clusterCollName_ = parameters.getParameter<std::string>("clusterCollName");
@@ -40,6 +44,10 @@ void EcalClusterProducer::configure(framework::config::Parameters& parameters) {
 void EcalClusterProducer::produce(framework::Event& event) {
   std::vector<ldmx::EcalHit> ecalHits =
       event.getCollection<ldmx::EcalHit>(recHitCollName_, recHitPassName_);
+
+  std::vector<ldmx::SimCalorimeterHit> ecalSimHits =
+      event.getCollection<ldmx::SimCalorimeterHit>(simHitCollName_, simHitPassName_);
+
 
   // Don't do anything if there are no ECal digis!
   if (!(ecalHits.size() > 0)) {
@@ -73,6 +81,7 @@ void EcalClusterProducer::produce(framework::Event& event) {
       cluster.setNHits(wcVec[aWC].getHits().size());
       cluster.addHits(wcVec[aWC].getHits());
       cluster.addFirstLayerHits(fWcVec[aWC].getHits());
+      cluster.findHitOrigins(ecalSimHits);
 
       histograms_.fill("nHits", wcVec[aWC].getHits().size());
       histograms_.fill("cluster_energy", wcVec[aWC].centroid().E());

Ecal/src/Ecal/EcalVetoProcessor.cxx

@@ -1231,6 +1231,7 @@ std::vector<std::pair<float, float>> EcalVetoProcessor::getTrajectory(
     std::vector<double> momentum, std::vector<float> position) {
   std::vector<XYCoords> positions;
   for (int iLayer = 0; iLayer < nEcalLayers_; iLayer++) {
+    //std::cout << iLayer << " " << geometry_->getZPosition(iLayer) << std::endl;
     float posX =
         position[0] + (momentum[0] / momentum[2]) *
                           (geometry_->getZPosition(iLayer) - position[2]);

Ecal/src/Ecal/Event/EcalCluster.cxx

@@ -30,4 +30,41 @@ ClassImp(ldmx::EcalCluster)
     }
     firstLayerHitIDs_ = vecIDs;
   }
+
+  void EcalCluster::findHitOrigins(const std::vector<ldmx::SimCalorimeterHit>& ecalSimHits) {
+
+    std::vector<unsigned int> vecIDs;
+    for (const auto& id : this->getHitIDs()) {
+      int tag = 0;
+      auto it = std::find_if(ecalSimHits.begin(), ecalSimHits.end(),
+        [&](const auto& simHit) { return simHit.getID() == id; });
+      if (it != ecalSimHits.end()) {
+        int ancestor = 0;
+        int prevAncestor = 0;
+        bool tagged = false;
+        tag = 0;
+        for (int i = 0; i < it->getNumberOfContribs(); i++) {
+          // for each contrib in this simhit
+          const auto& c = it->getContrib(i);
+          // get origin electron ID
+          ancestor = c.originID;
+          if (!tagged && i != 0 && prevAncestor != ancestor) {
+            // if origin electron ID does not match previous origin electron ID
+            // this hit has contributions from several electrons, ie mixed case
+            tag = 0;
+            tagged = true;
+          }
+          prevAncestor = ancestor;
+        }
+        if (!tagged) {
+          // if not tagged, hit was from a single electron
+          tag = prevAncestor;
+        }
+      }
+      else {tag = -1;}
+      vecIDs.push_back(tag);
+    }
+    hitOriginIDs_ = vecIDs;           
+  }
+
 }  // namespace ldmx

SimCore/include/SimCore/Event/SimCalorimeterHit.h

@@ -241,6 +241,7 @@ class SimCalorimeterHit {
    * @param pdgCode The PDG code of the actual track.
    * @param edep The energy deposition of the hit [MeV].
    * @param time The time of the hit [ns].
+   * @param originID the Geant4 track ID for the particle's electron parent.
    */
   void addContrib(int incidentID, int trackID, int pdgCode, float edep,
                   float time, int originID = -1);

SimCore/src/SimCore/Event/SimCalorimeterHit.cxx

@@ -16,6 +16,7 @@ ClassImp(ldmx::SimCalorimeterHit)
     pdgCodeContribs_.clear();
     edepContribs_.clear();
     timeContribs_.clear();
+    originContribs_.clear();
 
     nContribs_ = 0;
     id_ = 0;

Tracking/include/Tracking/Event/Track.h

@@ -205,6 +205,25 @@ class Track {
 
   void addTrackState(const ldmx::Track::TrackState& ts) {
     trackStates_.push_back(ts);
+
+    if (ts.ts_type == 1) {
+      trackStateTargetLoc_[0] = ts.refX;
+      trackStateTargetLoc_[1] = ts.refY;
+      trackStateTargetLoc_[2] = ts.refZ;
+
+      trackStateTargetParams_ = ts.params;
+      trackStateTargetCov_ = ts.cov;
+    }
+
+    if (ts.ts_type == 4) {
+      trackStateECalLoc_[0] = ts.refX;
+      trackStateECalLoc_[1] = ts.refY;
+      trackStateECalLoc_[2] = ts.refZ;
+
+      trackStateECalParams_ = ts.params;
+      trackStateECalCov_ = ts.cov;
+
+    }
   };
 
   std::vector<TrackState> getTrackStates() const { return trackStates_; }
@@ -268,6 +287,15 @@ class Track {
   // pdgID
   int pdgID_{0};
 
+    // TrackStates leaves
+  std::vector<double> trackStateTargetLoc_{-999., -999., -999.};
+  std::vector<double> trackStateTargetParams_{-999., -999., -999.,-999., -999., -999.};
+  std::vector<double> trackStateTargetCov_;
+
+  std::vector<double> trackStateECalLoc_{-999., -999., -999.};
+  std::vector<double> trackStateECalParams_{-999., -999., -999.,-999., -999., -999.};
+  std::vector<double> trackStateECalCov_;
+
   // Track States
   std::vector<TrackState> trackStates_;