modified: tests/simple_example_out.root

src/GenFind/include/HoughTransform.h

@@ -44,6 +44,19 @@ namespace GenFind {
 
       int FindPeaks(const std::vector<ROOT::Math::XYZTVector>& hits);
 
+      // crude idea:
+      //struct PreTrack {
+      //  double phi_peak;
+      //  ROOT::Math::XYZTVector ref_hit; // image or conformal space?
+      //  std::vector<ROOT::Math::XYZTVector> hits; // image space
+      //  std::vector<ROOT::Math::XYZTVector> chits; // conformal space
+      //}
+
+      //PreTrack GetPreTrack(
+      //    const std::vector<std::tuple<ROOT::Math::XYZTVector,ROOT::Math::XYZTVector>>& hits,
+      //    ROOT::Math::XYZTVector ref_hit,
+      //    double dphi = 5.0 /*degrees*/);
+
       std::vector<std::vector<std::tuple<ROOT::Math::XYZTVector,ROOT::Math::XYZTVector,ROOT::Math::XYZTVector>>> GetPreTrackHits(
           const std::vector<std::tuple<ROOT::Math::XYZTVector,ROOT::Math::XYZTVector,ROOT::Math::XYZTVector>>& hits,
           double dphi = 5.0 /*degrees*/);

src/GenFind/src/HoughTransform.cxx

@@ -14,10 +14,10 @@ namespace GenFind {
   {
     bool status = TH1::AddDirectoryStatus();
     TH1::AddDirectory(false);
-    fPhi       = new TH1F{"hPhi","hPhi",45,-90, 90};
-    fRhoTheta0 = new TH2F{"hRhoTheta0","hRhoTheta0",120,0, 360, 100,-0.02,0.02};
-    fRhoTheta1 = new TH2F{"hRhoTheta1","hRhoTheta1",90, 0, 360, 100,-0.02,0.02};
-    fRhoTheta2 = new TH2F{"hRhoTheta2","hRhoTheta2",40, 0, 360, 100,-0.02,0.02};
+    fPhi       = new TH1F{"hPhi","hPhi",90,-90, 90};
+    fRhoTheta0 = new TH2F{"hRhoTheta0","hRhoTheta0",120,0, 180, 100,-0.1,0.1};
+    fRhoTheta1 = new TH2F{"hRhoTheta1","hRhoTheta1",90, 0, 180, 100,-0.1,0.1};
+    fRhoTheta2 = new TH2F{"hRhoTheta2","hRhoTheta2",40, 0, 180, 10,-0.1,0.1};
     fSpec  = new TSpectrum2(2*fMaxTracks);
     fSpec1 = new TSpectrum(2*fMaxTracks);
     TH1::AddDirectory(status);
@@ -94,12 +94,35 @@ namespace GenFind {
         //std::cout << " root is at " << root << std::endl;
       }
     }
-      //  for(int i_theta = 0; i_theta<500; i_theta++) {
-      //    double theta = rand.Uniform(0.0,360.0)*degree;
-      //    fRhoTheta0->Fill(theta/degree, ahit.X()*TMath::Cos(theta) + ahit.Y()*TMath::Sin(theta) );
-      //    fRhoTheta1->Fill(theta/degree, ahit.X()*TMath::Cos(theta) + ahit.Y()*TMath::Sin(theta) );
-      //    fRhoTheta2->Fill(theta/degree, ahit.X()*TMath::Cos(theta) + ahit.Y()*TMath::Sin(theta) );
-      //  }
+
+    TRandom3 rand;
+    TH1F* temp0 = (TH1F*)fRhoTheta0->Clone();
+    TH1F* temp1 = (TH1F*)fRhoTheta1->Clone(); 
+    TH1F* temp2 = (TH1F*)fRhoTheta2->Clone(); 
+    fRhoTheta0->Reset();
+    fRhoTheta1->Reset();
+    fRhoTheta2->Reset();
+    for(const auto& ahit : hits) {
+      temp0->Reset();
+      temp1->Reset();
+      temp2->Reset();
+      for(int i_theta = 0; i_theta<500; i_theta++) {
+
+        double theta = rand.Uniform(0.0,180.0)*degree;
+
+        int bin = temp0->FindBin(theta/degree, ahit.X()*TMath::Cos(theta) + ahit.Y()*TMath::Sin(theta) );
+        temp0->SetBinContent(bin, 1);
+
+        bin = temp1->FindBin(theta/degree, ahit.X()*TMath::Cos(theta) + ahit.Y()*TMath::Sin(theta) );
+        temp1->SetBinContent(bin, 1);
+
+        bin = temp2->FindBin(theta/degree, ahit.X()*TMath::Cos(theta) + ahit.Y()*TMath::Sin(theta) );
+        temp2->SetBinContent(bin, 1);
+      }
+      fRhoTheta0->Add(temp0);
+      fRhoTheta1->Add(temp1);
+      fRhoTheta2->Add(temp2);
+    }
 
     //std::cout << " average : " << avg/double(n_roots) << std::endl;
     //double z_max  = fRhoTheta1->GetBinContent(fRhoTheta1->GetMaximumBin());

src/GenFind/tests/hough_transform2.cxx

@@ -25,7 +25,7 @@ void hough_transform2(
   TH2F * huv   = new TH2F("huv",   "huv",   100,   -0.12, 0.12  , 100,   -0.12, 0.12  );
   TH2F * hrphi = new TH2F("hrphi", "hrphi",  90,    0, 180, 100, -0.2, 0.2);
   TH1F * hphi  = new TH1F("hphi",  "hphi",   50,   -4, 4 );
-  TH1F * htheta  = new TH1F("htheta",  "htheta", 90,   0, 180 );
+  TH1F * htheta  = new TH1F("htheta",  "htheta", 180,   0, 180 );
 
   std::vector<XYZTVector> hits ;
 
@@ -67,16 +67,16 @@ void hough_transform2(
   std::cout << " Master " << master_hits.size() << std::endl;
 
 
-  //for(auto ahit : all_chits) {
-  //  std::cout << ahit.X() << " , " <<  ahit.Y() << std::endl;
-  //  huv->Fill( ahit.X(), ahit.Y() );
-  //  hphi->Fill(TMath::ATan(ahit.Y()/ahit.X()));
-  //  for(int i_theta = 0; i_theta<180; i_theta++) {
-  //    double theta = 2.0*double(i_theta)*CLHEP::degree;
-  //    hrphi->Fill(theta/CLHEP::degree, ahit.X()*TMath::Cos(theta) + ahit.Y()*TMath::Sin(theta) );
-  //  }
-  //}
-  //auto mg = ht->FillHoughTransforms(all_chits);
+  for(auto ahit : all_chits) {
+    std::cout << ahit.X() << " , " <<  ahit.Y() << std::endl;
+    huv->Fill( ahit.X(), ahit.Y() );
+    hphi->Fill(TMath::ATan(ahit.Y()/ahit.X()));
+    for(int i_theta = 0; i_theta<180; i_theta++) {
+      double theta = 2.0*double(i_theta)*CLHEP::degree;
+      hrphi->Fill(theta/CLHEP::degree, ahit.X()*TMath::Cos(theta) + ahit.Y()*TMath::Sin(theta) );
+    }
+  }
+  auto mg = ht->FillHoughTransforms(all_chits);
 
   std::cout << " n peaks : " << ht->FindPeaks(all_chits) << std::endl;
 
@@ -124,15 +124,22 @@ void hough_transform2(
   for(auto ahist : hists){
   ahist->Draw("same");
   }
-  //mg->Draw("a");
-  //mg->GetYaxis()->SetRangeUser(-0.10,0.10);
-  //mg->Draw("a");
   //ht->fRhoTheta0->Draw("colz");
   c->cd(2);
   htheta->Draw();
-  //ht->fRhoTheta1->Draw("colz");
-  //c->cd(3);
-  //ht->fRhoTheta2->Draw("colz");
 
+  c->cd(3);
+  mg->Draw("a");
+  mg->GetYaxis()->SetRangeUser(-0.10,0.10);
+  mg->Draw("a");
+
+  c = new TCanvas();
+  c->Divide(2,2);
+  c->cd(1);
+  ht->fRhoTheta0->Draw("colz");
+  c->cd(2);
+  ht->fRhoTheta1->Draw("colz");
+  c->cd(3);
+  ht->fRhoTheta2->Draw("colz");
 
 }