RCNP Simulation

AtData/AtHit.h

@@ -88,7 +88,7 @@ class AtHit : public TObject {
    Double_t GetTimeStampCorrInter() const { return fTimeStampCorrInter; }
    const std::vector<AtHit::MCSimPoint> &GetMCSimPointArray() const { return fMCSimPointArray; }
 
-   static Bool_t SortHit(const AtHit &lhs, const AtHit &rhs) { return lhs.GetPadNum() < rhs.GetPadNum(); }
+   static Bool_t SortHit(const AtHit &lhs, const AtHit &rhs) { return lhs.GetPadNum() > rhs.GetPadNum(); }
    static Bool_t SortHit(const std::unique_ptr<AtHit> &lhs, const std::unique_ptr<AtHit> &rhs)
    {
       return SortHit(*lhs, *rhs);

AtData/AtTrack.h

@@ -80,7 +80,7 @@ class AtTrack : public TObject {
    const HitVector &GetHitArray() const { return fHitArray; }
 
    std::vector<AtHit> GetHitArrayObject() { return ContainerManip::GetObjectVector(fHitArray); }
-   // const std::vector<AtHit> &GetHitArrayConst() const { return fHitArray; }
+   //const std::vector<AtHit> &GetHitArrayConst() const { return fHitArray; }
    const AtPatterns::AtPattern *GetPattern() const { return fPattern.get(); }
 
    Double_t GetGeoTheta() const { return fGeoThetaAngle; }

AtReconstruction/AtFitter/AtGenfit.cxx

@@ -81,7 +81,7 @@ AtFITTER::AtGenfit::AtGenfit(Float_t magfield, Float_t minbrho, Float_t maxbrho,
    std::cout << " AtFITTER::AtGenfit::AtGenfit(): Checking materials that GENFIT will use "
              << "\n";
 
-   auto *gGeoMan = dynamic_cast<TGeoManager *>(gROOT->FindObject("FAIRGeom"));
+   /*auto *gGeoMan = dynamic_cast<TGeoManager *>(gROOT->FindObject("FAIRGeom"));
    TObjArray *volume_list = gGeoMan->GetListOfVolumes();
    if (!volume_list) {
       std::cout << cRED << " Warning! Null list of geometry volumes." << cNORMAL << "\n";
@@ -105,7 +105,7 @@ AtFITTER::AtGenfit::AtGenfit(Float_t magfield, Float_t minbrho, Float_t maxbrho,
       // Int_t mat_indx = mat->GetIndex();
       std::cout << cYELLOW << " - Material : " << mat->GetName() << cNORMAL << "\n";
    }
-
+   */
    // PDG definitions
    const Double_t kAu2Gev = 0.9314943228;
    const Double_t khSlash = 1.0545726663e-27;
@@ -609,7 +609,7 @@ std::vector<std::unique_ptr<AtFittedTrack>> AtFITTER::AtGenfit::ProcessTracks(st
                   // Backward extrapolation
                   try {
 
-                     for (auto iStep = 0; iStep < 200; ++iStep) {
+                     for (auto iStep = 0; iStep < 600; ++iStep) {
 
                         trackRep->extrapolateBy(fitState, stepXtr * iStep);
                         mom_ext_buff = fitState.getMom();

AtReconstruction/AtPatternRecognition/AtPRA.cxx

@@ -48,7 +48,28 @@ std::cout << " Processing track with " << track.GetHitArray().size() << " points
    RansacSmoothRadius.SetMinHitsPattern(0.1 * track.GetHitArray().size());
    RansacSmoothRadius.SetDistanceThreshold(6.0);
    RansacSmoothRadius.SetNumIterations(1000);
-   auto circularTracks = RansacSmoothRadius.Solve(ContainerManip::GetConstPointerVector(track.GetHitArray()))
+   track.SortHitArrayTime();
+   auto &hitArray = track.GetHitArray();
+   size_t numHits = hitArray.size();
+   size_t numHitsToUse;
+   if (numHits >= 50)
+   numHitsToUse = numHits * 10 / 10; // 50% of the hits
+
+   if (numHits < 50)
+   numHitsToUse = numHits; 
+
+   std::sort(hitArray.begin(), hitArray.end(), [](const std::unique_ptr<AtHit> &a, const std::unique_ptr<AtHit> &b) {
+        return a->GetTimeStamp() > b->GetTimeStamp();
+    });
+
+   std::vector<std::unique_ptr<AtHit>> first20PercentHits;
+   first20PercentHits.reserve(numHitsToUse);
+
+   for (size_t i = 0; i < numHitsToUse; ++i) {
+    first20PercentHits.push_back(std::make_unique<AtHit>(*hitArray[i]));
+   }
+
+   auto circularTracks = RansacSmoothRadius.Solve(ContainerManip::GetConstPointerVector(first20PercentHits)) // Only part of the spiral is used
                             .GetTrackCand(); // Only part of the spiral is used
                                              // This function also sets the coefficients
                                              // i.e. radius of curvature and center
@@ -60,7 +81,7 @@ std::cout << " Processing track with " << track.GetHitArray().size() << " points
       // auto circle = dynamic_cast<const AtPatterns::AtPatternCircle2D *>(circularTracks.at(0).GetPattern());
 
       auto circle = std::make_unique<AtPatterns::AtPatternCircle2D>();
-      circle->AtPattern::FitPattern(ContainerManip::GetConstPointerVector(track.GetHitArray()));
+      circle->AtPattern::FitPattern(ContainerManip::GetConstPointerVector(first20PercentHits));
 
       auto center = circle->GetCenter();
       auto radius = circle->GetRadius();

AtReconstruction/AtPatternRecognition/AtSampleConsensus.cxx

@@ -60,6 +60,8 @@ AtSampleConsensus::GeneratePatternFromHits(const std::vector<const AtHit *> &hit
 AtPatternEvent AtSampleConsensus::Solve(AtEvent *event)
 {
    auto hitVec = ContainerManip::GetConstPointerVector(event->GetHits());
+      std::cout<<"AtSampleConsensus::Solve"<<std::endl;
+
    return Solve(hitVec, event);
 }
 
@@ -71,6 +73,7 @@ AtPatternEvent AtSampleConsensus::Solve(const std::vector<AtHit> &hitArray, AtBa
 
 AtPatternEvent AtSampleConsensus::Solve(const std::vector<const AtHit *> &hitArray, AtBaseEvent *event)
 {
+
    // Return early if we were passed an event and it is marked bad
    if (event != nullptr && !event->IsGood())
       return {*event};

AtReconstruction/AtPatternRecognition/AtTrackFinderTC.cxx

@@ -4,7 +4,9 @@
 #include "AtHit.h"              // for AtHit
 #include "AtPatternEvent.h"     // for AtPatternEvent
 #include "AtTrack.h"            // for AtTrack
-#include "AtTrackTransformer.h" // for AtTrackTransformer
+#include "AtTrackTransformer.h"
+#include "AtFitter.h" // for AtTrackTransformer
+#include "AtTrackFinder.h"
 
 #include <Math/Point3D.h> // for PositionVector3D
 
@@ -25,6 +27,11 @@ constexpr auto cRED = "\033[1;31m";
 constexpr auto cYELLOW = "\033[1;33m";
 constexpr auto cNORMAL = "\033[0m";
 constexpr auto cGREEN = "\033[1;32m";
+std::vector<AtTrack *> candTrackPool;
+std::vector<AtTrack> mergedTrackPool;
+Bool_t fEnableSingleVertexTrack = kTRUE;
+Double_t fClusterSize = 20.0;
+
 
 AtPATTERN::AtTrackFinderTC::AtTrackFinderTC() : AtPATTERN::AtPRA() {}
 
@@ -87,12 +94,24 @@ std::unique_ptr<AtPatternEvent> AtPATTERN::AtTrackFinderTC::FindTracks(AtEvent &
    add_clusters(cloud_xyz, cl_group, opt_params.is_gnuplot());
 
    // Post processing
-   // process_pointcloud(cloud_xyz, 25, 0);
+   //process_pointcloud(cloud_xyz, 25, 0);
 
    // Adapt clusters to AtTrack
    return clustersToTrack(cloud_xyz, cl_group, event);
 }
 
+const double tolerance = 3.0;
+Bool_t SameTrack(AtTrack *trA, AtTrack *trB)
+{
+   auto theta0 = trA->GetGeoTheta();
+   auto theta = trB->GetGeoTheta();
+
+   if (std::abs((theta0 * TMath::RadToDeg()) - (theta * TMath::RadToDeg())) <= tolerance)
+      return true;
+   else
+      return false;
+}
+
 void AtPATTERN::AtTrackFinderTC::eventToClusters(AtEvent &event, PointCloud &cloud)
 {
    Int_t nHits = event.GetNumHits();
@@ -112,7 +131,6 @@ void AtPATTERN::AtTrackFinderTC::eventToClusters(AtEvent &event, PointCloud &clo
 std::unique_ptr<AtPatternEvent>
 AtPATTERN::AtTrackFinderTC::clustersToTrack(PointCloud &cloud, const std::vector<cluster_t> &clusters, AtEvent &event)
 {
-
    std::vector<AtTrack> tracks;
    // std::vector<Point> points = cloud;
    auto points = cloud;
@@ -160,11 +178,197 @@ AtPATTERN::AtTrackFinderTC::clustersToTrack(PointCloud &cloud, const std::vector
    for (const auto &point : points)
       retEvent->AddNoise(event.GetHit(point.GetID()));
 
+   for (auto &track : tracks) {
+      if (track.GetHitArray().size() > 0)
+         SetTrackInitialParameters(track);
+      //retEvent->AddTrack(std::move(track));
+   }
+
+   bool kMergeTracks = false;
+
+   std::vector<AtTrack> mergedTracks;
+   if(tracks.size() > 2){
+
+      //const double tolerance = 2.0;
+      std::vector<bool> processed(tracks.size(), false);
+
+      for(Int_t numtr = 0; numtr < tracks.size(); numtr++){
+         std::cout << "Theta " << tracks.at(numtr).GetGeoTheta() << std::endl;
+         if (processed[numtr]) continue;
+
+         auto track0 = tracks.at(numtr);  
+         auto hitArray0 = track0.GetHitArrayObject();
+         auto theta0 = track0.GetGeoTheta();
+
+         AtTrack newTrack;
+
+         for(auto &hit : hitArray0){
+            newTrack.AddHit(hit);
+
+         }
+
+         if (numtr + 1 < tracks.size()) {
+
+            for (size_t tr = numtr + 1; tr < tracks.size(); ++tr) {
+               if (processed[tr]) continue;
+
+               auto track = tracks.at(tr);
+               auto hitArray = track.GetHitArrayObject();
+               auto theta = track.GetGeoTheta();
+
+               if (std::abs((theta0 * TMath::RadToDeg()) - (theta * TMath::RadToDeg())) <= tolerance) {
+                  for(auto &hit : hitArray){
+                     newTrack.AddHit(hit);
+
+                  }
+                  processed[tr] = true;
+               }
+
+            }  
+         }
+
+         processed[numtr] = true;
+         fTrackTransformer->ClusterizeSmooth3D(newTrack, fClusterRadius, fClusterDistance);             
+         mergedTracks.push_back(newTrack);
+
+      }
+
+      kMergeTracks = true;  
+   }
+
+   if(kMergeTracks){
+    std::swap(tracks, mergedTracks);
+   }
+
    for (auto &track : tracks) {
       if (track.GetHitArray().size() > 0)
          SetTrackInitialParameters(track);
       retEvent->AddTrack(std::move(track));
    }
 
+
+
    return retEvent;
 }
+
+
+
+
+
+//First try of merging tracks
+ /*  bool kMergeTracks = false;
+
+   if(tracks.size() > 2){
+
+      const double tolerance = 2.0;
+      std::vector<bool> processed(tracks.size(), false);
+
+      for(Int_t numtr = 0; numtr < tracks.size(); numtr++){
+         std::cout << "Theta " << tracks.at(numtr).GetGeoTheta() << std::endl;
+         if (processed[numtr]) continue;
+
+         auto track0 = tracks.at(numtr);  
+         auto hitArray0 = track0.GetHitArrayObject();
+         auto theta0 = track0.GetGeoTheta();
+
+         AtTrack newTrack;
+
+         for(auto &hit : hitArray0){
+            newTrack.AddHit(hit);
+
+         }
+
+         if (numtr + 1 < tracks.size()) {
+
+            for (size_t tr = numtr + 1; tr < tracks.size(); ++tr) {
+               if (processed[tr]) continue;
+
+               auto track = tracks.at(tr);
+               auto hitArray = track.GetHitArrayObject();
+               auto theta = track.GetGeoTheta();
+
+               if (std::abs((theta0 * TMath::RadToDeg()) - (theta * TMath::RadToDeg())) <= tolerance) {
+                  for(auto &hit : hitArray){
+                     newTrack.AddHit(hit);
+
+                  }
+                  processed[tr] = true;
+               }
+
+            }  
+         }
+
+         processed[numtr] = true;
+         //fTrackTransformer->ClusterizeSmooth3D(newTrack, fClusterRadius, fClusterDistance);             
+         mergedTracks.push_back(newTrack);
+
+      }
+
+      kMergeTracks = true;  
+   }
+
+   if(kMergeTracks){
+    std::swap(tracks, mergedTracks);
+   }
+   */
+
+  /*
+   auto sp = std::unique_ptr<AtTrack[]>(new AtTrack[tracks.size()]);
+
+
+   for (auto iTrack = 0; iTrack < tracks.size(); ++iTrack) {
+      sp[iTrack] = tracks.at(iTrack);
+      candTrackPool.push_back(std::move(&sp[iTrack]));
+
+   }
+
+   std::vector<AtTrack *> candToMergePool;
+   AtTrackFinder merger;
+   for (auto itA = candTrackPool.begin(); itA != candTrackPool.end(); ++itA) {
+      if(candTrackPool.size() == 0) break;
+      AtTrack *trA = *(itA);
+      if(candTrackPool.size() == 1){
+         mergedTrackPool.push_back(*trA);
+         break;   
+      } 
+      candToMergePool.clear();
+
+      auto itB = std::copy_if(itA + 1, candTrackPool.end(), std::back_inserter(candToMergePool),
+                                 [&trA, this](AtTrack *track) { return SameTrack(trA, track); });
+
+                                std::cout << "candToMergePool size after copy_if: " << candToMergePool.size() << std::endl;  
+
+      if (candToMergePool.size() > 0) { // Merge if matches are found
+         candToMergePool.push_back(trA);
+         Bool_t merged = merger.MergeTracks(&candToMergePool, &mergedTrackPool, fEnableSingleVertexTrack, fClusterRadius,
+                                        fClusterSize);
+
+            std::cout << "Merged: " << merged << std::endl;
+            std::cout << "candTrackPool size before erase: " << candTrackPool.size() << std::endl;
+
+         itA = candTrackPool.erase(std::remove_if(itA, candTrackPool.end(),
+                  [&trA, this](AtTrack *track) { return SameTrack(trA, track); }),
+         candTrackPool.end());
+
+         } else {
+            mergedTrackPool.push_back(*trA);
+                    itA = candTrackPool.erase(itA); // Erase the track from the pool
+         }
+
+   }
+
+   std::swap(tracks, mergedTrackPool);
+
+   candToMergePool.clear();
+   candTrackPool.clear();
+   //mergedTrackPool.clear();
+
+ /*  for (auto &track : tracks) {
+      if (track.GetHitArray().size() > 0)
+         SetTrackInitialParameters(track);
+      retEvent->AddTrack(std::move(track));
+   }*/
+   //std::cout << "Number of tracks: " << tracks.size() << std::endl;
+
+
+

AtReconstruction/AtPatternRecognition/triplclust/src/directedgraph.cxx

@@ -9,6 +9,8 @@
 /** @file */
 
 #include "directedgraph.h"
+#include "orthogonallsq.h"
+#include "pointcloud.h"
 
 #include "util.h"
 #include <assert.h> /* assert */
@@ -41,17 +43,17 @@ Graph::Graph(PointCloud &cloud, std::vector<size_t> _indices)
    assert(std::set<size_t>(_indices.begin(), _indices.end()).size() == _indices.size()); // indices are unique
    assert((indices.back() < cloud.size()));                                              // checks range of indices
 
-   this->indices = indices;
-
+   this->indices = _indices;
+  
    // insert Nodes without Edges. Nodes have x,y,z Coordinates and the corresponding index
    for (unsigned int i = 0; i < indices.size(); i++) {
       graph.insert(std::pair<size_t, Node>(
          indices[i], Node(indices[i], cloud[indices[i]].x, cloud[indices[i]].y, cloud[indices[i]].z)));
    }
-
+    
    // Add first subtree: the first Point by time and the last Point by Index
    this->subtree_roots.push_back(&graph.at(indices[indices.size() - 1]));
-
+   
    // construct MST by adding Edges
    constructMST();
 }

AtReconstruction/AtPatternRecognition/triplclust/src/orthogonallsq.h

@@ -13,6 +13,8 @@
 
 #include "pointcloud.h"
 
+class PointCloud;
+
 //
 // orthogonal least squares fit with libeigen
 //   pc: points