Added complete finder using the operator() to HoughTransform

This method uses only the phi distribution to find peaks.
TODO: reevaluate how the hits are passed around maybe it is better to
pass around the index?

	modified:   include/HoughTransform.h
	modified:   src/HoughTransform.cxx
	modified:   tests/ht_test.cxx

src/GenFind/include/HoughTransform.h

@@ -2,6 +2,7 @@
 #define genfind_HoughTransform_HH
 
 #include "TObject.h"
+#include "TMath.h"
 #include <array>
 #include <vector>
 #include <utility>
@@ -30,6 +31,8 @@ namespace genfind {
       int fThresh    = 3;
       int fMaxTracks = 20;
 
+      double fPhiThresh = 5.0; /// degrees
+
       std::vector<std::tuple<int,int,double>> fHTRoots;
       std::vector<double>    fPhiPeaks;
 
@@ -54,7 +57,31 @@ namespace genfind {
         std::vector<std::vector<T>> res; 
         T ref = {0.0,0.0,0.0,0.0};
         auto chits = genfind::compute_conformal_hits(hits, ref);
-        res.push_back(hits);
+        auto peaks = FindPhiPeaks(chits);
+        res = SelectPhiPeaks<T>(peaks, chits);
+        return res;
+      }
+
+      template<class T>
+      std::vector<std::vector<T>> SelectPhiPeaks(const std::vector<double>& peaks,
+                                                 const std::vector<genfind::ConformalHit>& chits)
+      {
+        double degree = TMath::Pi()/180.0;
+        std::vector<std::vector<T>> res;
+        for(auto& p : peaks){
+          std::vector<T> p_hits;
+          for(const auto& ahit : chits){
+            double phi  = TMath::ATan(ahit.Y()/ahit.X())/degree;
+            if( TMath::Abs(phi-p) < fPhiThresh ) {
+              p_hits.push_back(
+                {ahit.fImageHit->X(),
+                  ahit.fImageHit->Y(), 
+                  ahit.fImageHit->Z(),
+                  ahit.fImageHit->T()});
+            }
+          }
+          res.push_back(p_hits); 
+        }
         return res;
       }
 
@@ -65,6 +92,7 @@ namespace genfind {
        */
       void FillHistograms(const std::vector<genfind::ConformalHit>& chits);
       
+      std::vector<double> FindPhiPeaks(const std::vector<genfind::ConformalHit>& chits);
 
       TMultiGraph* FillHoughTransforms(const std::vector<ROOT::Math::XYZTVector>& hits);
 

src/GenFind/src/HoughTransform.cxx

@@ -92,6 +92,26 @@ namespace genfind {
   }
   //__________________________________________________________________________
 
+  std::vector<double> HoughTransform::FindPhiPeaks(const std::vector<genfind::ConformalHit>& chits)
+  {
+    double degree = TMath::Pi()/180.0;
+    // First fill the phi histogram
+    fPhi->Reset();
+    for(const auto& ahit : chits){
+      double phi  = TMath::ATan(ahit.Y()/ahit.X())/degree;
+      fPhi->Fill(phi);
+    }
+    double sigma  = 2;
+    double thresh = 0.1;
+    int    npeaks = fSpec1->Search(fPhi, sigma, "nodraw", thresh );
+    std::vector<double> peaks;
+    for(int i = 0; i<npeaks; i++) {
+      peaks.push_back( (fSpec1->GetPositionX())[i] );
+    }
+    return peaks;
+  }
+  //__________________________________________________________________________
+
   TMultiGraph* HoughTransform::FillHoughTransforms(const std::vector<ROOT::Math::XYZTVector>& hits)
   {
     double degree = TMath::Pi()/180.0;

src/GenFind/tests/ht_test.cxx

@@ -31,12 +31,18 @@ void ht_test()
 
   HoughTransform ht;
   auto thits = ht(hits2);
+
   std::cout << "ht test\n";
   for(auto ahit : thits[0]){
     std::cout << " x " << ahit.X() << ", y " << ahit.Y() << std::endl;
   }
   std::cout << "ht test\n";
 
+  auto peaks = ht.FindPhiPeaks(chits2);
+  for(auto ahit : peaks){
+    std::cout << " phi " << ahit << std::endl;
+  }
+
   TGraph* gr_hits = new TGraph();
   int i_point = 0;
   for(auto ahit : hits){
@@ -54,7 +60,7 @@ void ht_test()
 
   TGraph* gr_hits2 = new TGraph();
   i_point = 0;
-  for(auto ahit : hits2){
+  for(auto ahit : thits[0]){
      gr_hits2->SetPoint(i_point, ahit.X(), ahit.Y());
      i_point++;
   }
@@ -68,23 +74,28 @@ void ht_test()
   }
   
   TCanvas* c = new TCanvas();
+  TMultiGraph* mg = new TMultiGraph();
   c->Divide(2,2);
 
   c->cd(1);
   gr_hits->SetMarkerStyle(20);
-  gr_hits->Draw("ap");
+  gr_hits2->SetMarkerStyle(21);
+  gr_hits2->SetMarkerColor(2);
+  mg->Add(gr_hits,"p");
+  mg->Add(gr_hits2,"p");
+  mg->Draw("a");
 
   c->cd(2);
   gr_chits->SetMarkerStyle(20);
   gr_chits->Draw("ap");
 
   c->cd(3);
-  gr_hits2->SetMarkerStyle(20);
   gr_hits2->Draw("ap");
 
   c->cd(4);
-  gr_chits2->SetMarkerStyle(20);
-  gr_chits2->Draw("ap");
+  ht.fPhi->DrawCopy();
+  //gr_chits2->SetMarkerStyle(20);
+  //gr_chits2->Draw("ap");
 
 
 }