first pass at nontrivial pflow

LDMX-Software · Jul 21, 2023 · 16fa557 · 16fa557
1 parent 9e392d7
commit 16fa557
Show file tree

Hide file tree

Showing 5 changed files with 371 additions and 152 deletions.
diff --git a/Recon/include/Recon/ParticleFlow.h b/Recon/include/Recon/ParticleFlow.h
@@ -13,6 +13,13 @@
 #include "Framework/EventProcessor.h"  //Needed to declare processor
 #include "TGraph.h"
 
+#include "SimCore/Event/SimParticle.h"
+#include "SimCore/Event/SimTrackerHit.h"
+#include "Recon/Event/CaloCluster.h"
+#include "Ecal/Event/EcalCluster.h"
+#include "Hcal/Event/HcalCluster.h"
+#include "Recon/Event/PFCandidate.h"
+
 namespace recon {
 
 /**
@@ -36,6 +43,10 @@ class ParticleFlow : public framework::Producer {
 
   virtual void onProcessEnd();
 
+  void FillCandTrack(ldmx::PFCandidate &cand, const ldmx::SimTrackerHit &tk);
+  void FillCandEMCalo(ldmx::PFCandidate &cand, const ldmx::CaloCluster &em);
+  void FillCandHadCalo(ldmx::PFCandidate &cand, const ldmx::CaloCluster &had);
+
  private:
   // specific verbosity of this producer
   int verbose_{0};
@@ -49,6 +60,8 @@ class ParticleFlow : public framework::Producer {
   std::string inputTrackCollName_;
   // name of collection for PF outputs
   std::string outputCollName_;
+  // configuration
+  bool singleParticle_;
 };
 }  // namespace recon
 

diff --git a/Recon/python/pfReco.py b/Recon/python/pfReco.py
@@ -49,6 +49,7 @@ def __init__(self, name='PFlow') :
         self.inputHcalCollName  = 'PFHcalClusters'
         self.inputTrackCollName = 'PFTracks'
         self.outputCollName     = 'PFCandidates'
+        self.singleParticle     = False
 
 class pfTruthProducer(ldmxcfg.Producer) :
     """Configuration for track selector for particle reco"""

diff --git a/Recon/src/Recon/DBScanClusterBuilder.cxx b/Recon/src/Recon/DBScanClusterBuilder.cxx
@@ -82,9 +82,9 @@ namespace recon {
   float e(0),x(0),y(0),z(0),xx(0),yy(0),zz(0),n(0);
   float w = 1; // weight
   float sumw = 0;
-  std::vector<float> xvals{};
-  std::vector<float> yvals{};
-  std::vector<float> zvals{};
+  // std::vector<float> xvals{};
+  // std::vector<float> yvals{};
+  // std::vector<float> zvals{};
   std::vector<float> raw_xvals{};
   std::vector<float> raw_yvals{};
   std::vector<float> raw_zvals{};
@@ -102,9 +102,9 @@ namespace recon {
     zz += w * h->getZPos() * h->getZPos();
     n += 1;
     sumw += w;
-    xvals.push_back(x);
-    yvals.push_back(y);
-    zvals.push_back(z);
+    // xvals.push_back(x); // unused
+    // yvals.push_back(y);
+    // zvals.push_back(z);
     raw_xvals.push_back(h->getXPos());
     raw_yvals.push_back(h->getYPos());
     raw_zvals.push_back(h->getZPos());
@@ -128,20 +128,28 @@ namespace recon {
   cl->setHitValsZ(raw_zvals);
   cl->setHitValsE(raw_evals);
 
-  if (xvals.size()>2){
-    for(int i=0; i<xvals.size();i++){ // mean subtract
-      xvals[i] = xvals[i] - x;
-      yvals[i] = yvals[i] - y;
-      zvals[i] = zvals[i] - z;
+  if (raw_xvals.size()>2){
+    // skip fits for 'vertical' clusters
+    std::vector<float> sortedZ = raw_zvals;
+    std::sort(sortedZ.begin(), sortedZ.end());
+    if (sortedZ[sortedZ.size()-1] - sortedZ[0] > 1e3){
+
+      for(int i=0; i<raw_xvals.size();i++){ // mean subtract
+	raw_xvals[i] = raw_xvals[i] - x;
+	raw_yvals[i] = raw_yvals[i] - y;
+	raw_zvals[i] = raw_zvals[i] - z;
+      }
+
+      TGraph gxz(raw_zvals.size(), raw_zvals.data(), raw_xvals.data());
+      auto r_xz = gxz.Fit("pol1","SQ"); // p0 + x*p1
+      cl->setDXDZ( r_xz->Value(1) );
+      cl->setEDXDZ( r_xz->ParError(1) );
+
+      TGraph gyz(raw_zvals.size(), raw_zvals.data(), raw_yvals.data());
+      auto r_yz = gyz.Fit("pol1","SQ"); // p0 + x*p1
+      cl->setDYDZ( r_yz->Value(1) );
+      cl->setEDYDZ( r_yz->ParError(1) );
     }
-    TGraph gxz(zvals.size(), zvals.data(), xvals.data());
-    auto r_xz = gxz.Fit("pol1","SQ"); // p0 + x*p1
-    cl->setDXDZ( r_xz->Value(1) );
-    cl->setEDXDZ( r_xz->ParError(1) );
-    TGraph gyz(zvals.size(), zvals.data(), yvals.data());
-    auto r_yz = gyz.Fit("pol1","SQ"); // p0 + x*p1
-    cl->setDYDZ( r_yz->Value(1) );
-    cl->setEDYDZ( r_yz->ParError(1) );
   }
   return;
 }

diff --git a/Recon/src/Recon/PFTrackProducer.cxx b/Recon/src/Recon/PFTrackProducer.cxx
@@ -10,6 +10,11 @@ void PFTrackProducer::configure(framework::config::Parameters& ps) {
   outputTrackCollName_ = ps.getParameter<std::string>("outputTrackCollName"); 
 }
 
+double getP(const ldmx::SimTrackerHit &tk){
+  std::vector<double> pxyz = tk.getMomentum();
+  return sqrt(pow(pxyz[0],2)+pow(pxyz[1],2)+pow(pxyz[2],2));
+}
+
 void PFTrackProducer::produce(framework::Event& event) {
 
   if (!event.exists(inputTrackCollName_)) return;
@@ -26,7 +31,7 @@ void PFTrackProducer::produce(framework::Event& event) {
   }
   std::sort(pfTracks.begin(), pfTracks.end(),
 	    [](ldmx::SimTrackerHit a, ldmx::SimTrackerHit b) {
-	      return a.getEnergy() > b.getEnergy();
+	      return getP(a) > getP(b);
 	    });
   event.add(outputTrackCollName_, pfTracks);
 }