Adding tracker param init from cluster.

- All charges (-1,0,1) are added to the candidates
- The initial vector points to the EM cluster, set with momentum equal
to the energy where we assume a massless particle
- Added new seeding algorithm which uses Ecal clusters and Vertex
tracker hits.

JugReco/CMakeLists.txt

@@ -34,6 +34,7 @@ gaudi_add_module(JugRecoPlugins
   src/components/TestACTSLogger.cpp 
   src/components/TrackParamTruthInit.cpp 
   src/components/TrackParamClusterInit.cpp 
+  src/components/TrackParamVertexClusterInit.cpp 
   src/components/SimpleClustering.cpp 
   src/components/ParticlesFromTrackFit.cpp 
   src/components/CalorimeterIslandCluster.cpp 

JugReco/src/components/ParticlesFromTrackFit.cpp

@@ -15,13 +15,19 @@
 #include "JugBase/DataHandle.h"
 #include "JugBase/IGeoSvc.h"
 
+#include "Acts/EventData/MultiTrajectory.hpp"
+#include "Acts/EventData/MultiTrajectoryHelpers.hpp"
+ 
 // Event Model related classes
-//#include "GaudiExamples/MyTrack.h"
 #include "eicd/ParticleCollection.h"
 #include "eicd/TrackerHitCollection.h"
+#include "eicd/TrackParametersCollection.h"
 #include "JugReco/SourceLinks.h"
 #include "JugReco/Track.hpp"
 
+#include "Acts/Utilities/Helpers.hpp"
+
+
 namespace Jug {
   namespace Reco {
   
@@ -31,8 +37,9 @@ namespace Jug {
    class ParticlesFromTrackFit : public GaudiAlgorithm {
    public:
     //DataHandle<eic::RawTrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
-    DataHandle<TrajectoryContainer>      m_inputTrajectories{"inputTrajectories", Gaudi::DataHandle::Reader, this};
+    DataHandle<TrajectoryContainer>     m_inputTrajectories{"inputTrajectories", Gaudi::DataHandle::Reader, this};
     DataHandle<eic::ParticleCollection> m_outputParticles{"outputParticles", Gaudi::DataHandle::Writer, this};
+    DataHandle<eic::TrackParametersCollection> m_outputTrackParameters{"outputTrackParameters", Gaudi::DataHandle::Writer, this};
 
    public:
     //  ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
@@ -40,6 +47,7 @@ namespace Jug {
         : GaudiAlgorithm(name, svcLoc) {
           declareProperty("inputTrajectories", m_inputTrajectories,"");
           declareProperty("outputParticles", m_outputParticles, "");
+          declareProperty("outputTrackParameters", m_outputTrackParameters, "ACTS Track Parameters");
         }
 
     StatusCode initialize() override {
@@ -53,12 +61,68 @@ namespace Jug {
       const TrajectoryContainer* trajectories = m_inputTrajectories.get();
       // create output collections
       auto rec_parts = m_outputParticles.createAndPut();
+      auto track_pars = m_outputTrackParameters.createAndPut();
 
       for(const auto& traj : *trajectories) {
         //traj.trajectory().first
+        const auto& [trackTips, mj] = traj.trajectory();
+        if (trackTips.empty()) {
+          debug() << "Empty multiTrajectory." << endmsg;
+          continue;
+        }
+
+            // Get the entry index for the single trajectory
+    auto& trackTip = trackTips.front();
+
+    // Collect the trajectory summary info
+    auto trajState       = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
+    int m_nMeasurements = trajState.nMeasurements;
+    int m_nStates = trajState.nStates;
+
+            // Get the fitted track parameter
+    bool  m_hasFittedParams = false;
+    if (traj.hasTrackParameters(trackTip)) {
+      m_hasFittedParams = true;
+      const auto& boundParam = traj.trackParameters(trackTip);
+      const auto& parameter = boundParam.parameters();
+      const auto& covariance = *boundParam.covariance();
+      debug() << "loc 0 = " << parameter[Acts::eBoundLoc0]   << endmsg;
+      debug() << "loc 1 = " << parameter[Acts::eBoundLoc1]   << endmsg;
+      debug() << "phi   = " << parameter[Acts::eBoundPhi]    << endmsg;
+      debug() << "theta = " << parameter[Acts::eBoundTheta]  << endmsg;
+      debug() << "q/p   = " << parameter[Acts::eBoundQOverP] << endmsg;
+      debug() << "p     = " << 1.0/parameter[Acts::eBoundQOverP] << endmsg;
+
+      debug() << "err phi = " << sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi))      << endmsg;
+      debug() << "err th  = " << sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta))  << endmsg;
+      debug() << "err q/p = " << sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP))<< endmsg;
+
+      debug() << " chi2 = " << trajState.chi2Sum << endmsg;
+
+      eic::TrackParameters pars({
+          parameter[Acts::eBoundLoc0], parameter[Acts::eBoundLoc1], parameter[Acts::eBoundPhi],
+          parameter[Acts::eBoundTheta], parameter[Acts::eBoundQOverP],parameter[Acts::eBoundTime],
+          sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)), sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)),
+          sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)), sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)),
+          sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)),sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime))});
+      track_pars->push_back(pars);
+
+      //m_ePHI_fit = parameter[Acts::eBoundPhi];
+      //m_eTHETA_fit = parameter[Acts::eBoundTheta];
+      //m_eQOP_fit = parameter[Acts::eBoundQOverP];
+      //m_eT_fit = parameter[Acts::eBoundTime];
+      //m_err_eLOC0_fit = sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0));
+      //m_err_eLOC1_fit = sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1));
+      //m_err_ePHI_fit = sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi));
+      //m_err_eTHETA_fit = sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta));
+      //m_err_eQOP_fit = sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP));
+      //m_err_eT_fit = sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime));
+    }
+
         auto tsize = traj.trajectory().first.size();
         debug() << "# fitted parameters : " << tsize << endmsg;
         if(tsize == 0 ) continue;
+
         traj.trajectory().second.visitBackwards(tsize-1, [&](auto&& trackstate) {
           //debug() << trackstate.hasPredicted() << endmsg;
           //debug() << trackstate.predicted() << endmsg;

JugReco/src/components/TrackFindingAlgorithm.cpp

@@ -72,7 +72,7 @@ namespace Jug::Reco {
     m_fieldctx              = BFieldVariant(m_BField);
 
     // chi2 and #sourclinks per surface cutoffs
-    m_sourcelinkSelectorCfg = { {Acts::GeometryIdentifier(), {3, 3}},
+    m_sourcelinkSelectorCfg = { {Acts::GeometryIdentifier(), {15, 10}},
     };
 
     findTracks = TrackFindingAlgorithm::makeTrackFinderFunction(m_geoSvc->trackingGeometry(), m_BField);
@@ -113,8 +113,8 @@ namespace Jug::Reco {
       //                                                 &(*pSurface));
 
       debug() << "Invoke track finding seeded by truth particle " << iseed << endmsg;
-
       auto result = findTracks(*src_links, initialParams, ckfOptions);
+      debug() << "finding done." << endmsg;
       if (result.ok()) {
         // Get the track finding output object
         const auto& trackFindingOutput = result.value();
@@ -174,6 +174,7 @@ namespace Jug::Reco {
           using SourceLinkSelector = Acts::CKFSourceLinkSelector;
           using CKF                = Acts::CombinatorialKalmanFilter<Propagator, Updater, Smoother, SourceLinkSelector>;
 
+          std::cout << " finding ...\n";
           // construct all components for the track finder
           MagneticField field(std::move(inputField));
           Stepper       stepper(std::move(field));
@@ -181,6 +182,7 @@ namespace Jug::Reco {
           navigator.resolvePassive   = false;
           navigator.resolveMaterial  = true;
           navigator.resolveSensitive = true;
+          std::cout << " propagator\n" ;
           Propagator propagator(std::move(stepper), std::move(navigator));
           CKF        trackFinder(std::move(propagator));
 

JugReco/src/components/TrackParamClusterInit.cpp

@@ -80,7 +80,7 @@ namespace Jug::Reco {
         cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = 1.0 * mm*1.0 * mm;
         cov(Acts::eBoundPhi, Acts::eBoundPhi)     = M_PI / 180.0;
         cov(Acts::eBoundTheta, Acts::eBoundTheta) = M_PI / 180.0;
-        cov(Acts::eBoundQOverP, Acts::eBoundQOverP)     = 1.0 / (0.05*0.05*p*p);
+        cov(Acts::eBoundQOverP, Acts::eBoundQOverP)     = 1.0 / (p*p);
         cov(Acts::eBoundTime, Acts::eBoundTime)         = Acts::UnitConstants::ns;
 
         // add all charges to the track candidate...
@@ -90,9 +90,9 @@ namespace Jug::Reco {
         init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
                                       Acts::Vector3D(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 1,
                                       std::make_optional(cov));
-        init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
-                                      Acts::Vector3D(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 0,
-                                      std::make_optional(cov));
+        //init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
+        //                              Acts::Vector3D(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 0,
+        //                              std::make_optional(cov));
         debug() << "Invoke track finding seeded by truth particle with p = " << p/GeV  << " GeV" << endmsg;
       }
       return StatusCode::SUCCESS;

JugReco/src/components/TrackParamTruthInit.cpp

@@ -73,7 +73,7 @@ namespace Jug::Reco {
         using Acts::UnitConstants::GeV;
         using Acts::UnitConstants::mm;
 
-        double p = std::hypot( part.psx() * MeV, part.psy() * MeV, part.psz() * MeV);
+        double p = std::hypot( part.psx() * GeV, part.psy() * GeV, part.psz() * GeV);
 
         // build some track cov matrix
         Acts::BoundSymMatrix cov        = Acts::BoundSymMatrix::Zero();
@@ -85,7 +85,7 @@ namespace Jug::Reco {
         cov(Acts::eBoundTime, Acts::eBoundTime)         = Acts::UnitConstants::ns;
 
         init_trk_params->emplace_back(Acts::Vector4D(part.vsx() * mm, part.vsy() * mm, part.vsz() * mm, part.time() * Acts::UnitConstants::ns),
-                                      Acts::Vector3D(part.psx() * MeV, part.psy() * MeV, part.psz() * MeV),
+                                      Acts::Vector3D(part.psx() * GeV, part.psy() * GeV, part.psz() * GeV),
                                       p+200*MeV,
                                       ((part.pdgID() > 0) ? -1 : 1),
                                       std::make_optional(std::move(cov))

JugReco/src/components/TrackParamVertexClusterInit.cpp

+#include <cmath>
+// Gaudi
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "GaudiAlg/Transformer.h"
+#include "GaudiAlg/GaudiTool.h"
+#include "GaudiKernel/RndmGenerators.h"
+#include "GaudiKernel/Property.h"
+
+#include "JugBase/DataHandle.h"
+#include "JugBase/IGeoSvc.h"
+#include "JugReco/Track.hpp"
+#include "Acts/Utilities/Units.hpp"
+
+#include "eicd/TrackerHitCollection.h"
+#include "eicd/ClusterCollection.h"
+
+
+  ///// (Reconstructed) track parameters e.g. close to the vertex.
+  //using TrackParameters = Acts::CurvilinearTrackParameters;
+
+  ///// Container of reconstructed track states for multiple tracks.
+  //using TrackParametersContainer = std::vector<TrackParameters>;
+
+  ///// MultiTrajectory definition
+  //using Trajectory = Acts::MultiTrajectory<SourceLink>;
+
+  ///// Container for the truth fitting/finding track(s)
+  //using TrajectoryContainer = std::vector<SimMultiTrajectory>;
+
+namespace Jug::Reco {
+
+  /** Initial Track parameters from clusters and vertex hits.
+   *
+   *
+   * The momentum of the initial track is estimated from the cluster  energy and 
+   * the direction is set using the vertex hits.
+   *
+   */
+  class TrackParamVertexClusterInit : public GaudiAlgorithm {
+  public:
+    using Clusters =  eic::ClusterCollection;
+    using VertexHits =  eic::TrackerHitCollection;
+
+    DataHandle<VertexHits>               m_inputVertexHits{"inputVertexHits", Gaudi::DataHandle::Reader, this};
+    DataHandle<Clusters>                 m_inputClusters{"inputClusters", Gaudi::DataHandle::Reader, this};
+    DataHandle<TrackParametersContainer> m_outputInitialTrackParameters{"outputInitialTrackParameters",
+                                                                        Gaudi::DataHandle::Writer, this};
+
+  public:
+    TrackParamVertexClusterInit(const std::string& name, ISvcLocator* svcLoc)
+        : GaudiAlgorithm(name, svcLoc) {
+      declareProperty("inputVertexHits", m_inputVertexHits, "Vertex tracker hits");
+      declareProperty("inputClusters", m_inputClusters, "Input clusters");
+      declareProperty("outputInitialTrackParameters", m_outputInitialTrackParameters, "");
+    }
+
+    StatusCode initialize() override {
+      if (GaudiAlgorithm::initialize().isFailure())
+        return StatusCode::FAILURE;
+      IRndmGenSvc* randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
+      if(!randSvc)
+        return StatusCode::FAILURE;
+      return StatusCode::SUCCESS;
+    }
+
+    StatusCode execute() override {
+      // input collection
+      const Clusters*   clusters = m_inputClusters.get();
+      const VertexHits* vtx_hits = m_inputVertexHits.get();
+      // Create output collections
+      auto init_trk_params = m_outputInitialTrackParameters.createAndPut();
+
+      for(const auto& c : *clusters) {
+
+        using Acts::UnitConstants::GeV;
+        using Acts::UnitConstants::MeV;
+        using Acts::UnitConstants::mm;
+        using Acts::UnitConstants::ns;
+
+        double p = c.energy()*GeV;
+        if( p < 1.0) {
+          debug() << " skipping cluster with energy " << p/GeV << " GeV" << endmsg;
+          continue;
+        }
+
+        double len =  std::hypot( c.x() , c.y() , c.z() );
+
+        // build some track cov matrix
+        Acts::BoundSymMatrix cov        = Acts::BoundSymMatrix::Zero();
+        cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = 1.0 * mm*1.0 * mm;
+        cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = 1.0 * mm*1.0 * mm;
+        cov(Acts::eBoundPhi, Acts::eBoundPhi)     = M_PI / 180.0;
+        cov(Acts::eBoundTheta, Acts::eBoundTheta) = M_PI / 180.0;
+        cov(Acts::eBoundQOverP, Acts::eBoundQOverP)     = 1.0 / (p*p);
+        cov(Acts::eBoundTime, Acts::eBoundTime)         = Acts::UnitConstants::ns;
+
+        // add all charges to the track candidate...
+        init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
+                                      Acts::Vector3D(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, -1,
+                                      std::make_optional(cov));
+        init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
+                                      Acts::Vector3D(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 1,
+                                      std::make_optional(cov));
+        //init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
+        //                              Acts::Vector3D(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 0,
+        //                              std::make_optional(cov));
+        debug() << "Invoke track finding seeded by truth particle with p = " << p/GeV  << " GeV" << endmsg;
+      }
+      return StatusCode::SUCCESS;
+    }
+  };
+  DECLARE_COMPONENT(TrackParamVertexClusterInit)
+
+} // namespace Jug::reco
+

JugReco/src/components/TrackerSourceLinker.cpp

@@ -100,7 +100,7 @@ namespace Jug::Reco {
 
         auto vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.cellID());
         auto vol_id = vol_ctx->identifier;
-        debug() << " hit          : \n" <<  ahit << endmsg;
+        //debug() << " hit          : \n" <<  ahit << endmsg;
         //debug() << "cell ID : " << ahit.cellID() << endmsg;
         //debug() << " position : (" <<  ahit.position(0) << ", " <<  ahit.position(1) << ", "<<  ahit.position(2) << ") " << endmsg;
         debug() << " vol_id       : " <<  vol_id << endmsg;