Skip to content
GitLab
Commit 4fbde59f authored by Johnston's avatar Johnston
Browse files

using intersections in rho/theta space

parent 6caebbf0
Hide whitespace changes
Inline Side-by-side
src/GenFind/include/HoughTransform.h
View file @ 4fbde59f
......@@ -21,6 +21,7 @@ namespace genfind {
TH2F* fRhoTheta0 = nullptr;
TH2F* fRhoTheta1 = nullptr;
TH2F* fRhoTheta2 = nullptr;
TH2F* fIntersect = nullptr;
TSpectrum* fSpec1 = nullptr;
TSpectrum2* fSpec = nullptr;
......@@ -29,6 +30,7 @@ namespace genfind {
int fMaxTracks = 20;
std::vector<std::tuple<int,int,double>> fHTRoots;
std::vector<std::tuple<int,int,double,double>> fHTIntersections;
std::vector<double> fPhiPeaks;
public:
......@@ -42,6 +44,8 @@ namespace genfind {
TMultiGraph* FillHoughTransforms(const std::vector<ROOT::Math::XYZTVector>& hits);
int FindIntersections( TMultiGraph* mg );
int FindPeaks(const std::vector<ROOT::Math::XYZTVector>& hits);
// crude idea:
......
src/GenFind/src/HoughTransform.cxx
View file @ 4fbde59f
#include "HoughTransform.h"
#include "TMath.h"
#include "Math/WrappedTF1.h"
#include "Math/BrentRootFinder.h"
#include <algorithm>
#include "TRandom3.h"
......@@ -7,6 +9,8 @@
#include "Math/WrappedTF1.h"
#include "Math/BrentRootFinder.h"
#include "TGraph.h"
#include "TList.h"
#include "TObject.h"
namespace genfind {
......@@ -15,9 +19,10 @@ namespace genfind {
bool status = TH1::AddDirectoryStatus();
TH1::AddDirectory(false);
fPhi = new TH1F{"hPhi","hPhi",90,-90, 90};
fRhoTheta0 = new TH2F{"hRhoTheta0","hRhoTheta0",120,0, 180, 100,-0.1,0.1};
fRhoTheta0 = new TH2F{"hRhoTheta0","hRhoTheta0",100,-1010, 1010, 100,-1010,1010};
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};
fIntersect = new TH2F{"hIntersect","hIntersect",90,0, 180, 75,-150,150 };
fSpec = new TSpectrum2(2*fMaxTracks);
fSpec1 = new TSpectrum(2*fMaxTracks);
TH1::AddDirectory(status);
......@@ -30,6 +35,7 @@ namespace genfind {
delete fRhoTheta0;
delete fRhoTheta1;
delete fRhoTheta2;
delete fIntersect;
if(fSpec) delete fSpec; fSpec = nullptr;
if(fSpec1) delete fSpec1; fSpec1 = nullptr;
}
......@@ -60,6 +66,8 @@ namespace genfind {
TMultiGraph* mg = new TMultiGraph();
std::vector<TF1*> funcs;
double avg = 0;
for(const auto& ahit : hits) {
TF1 * f = new TF1("f",[&](double*x, double *p){ return p[0]*TMath::Cos(x[0]*degree) + p[1]*TMath::Sin(x[0]*degree); }, phi_min, phi_max, 2);
f->SetLineColor(1);
......@@ -68,8 +76,10 @@ namespace genfind {
f->SetParameter(1, ahit.Y());
funcs.push_back(f);
mg->Add(new TGraph(f),"l");
}
// eventually move this to another function call?
int n_roots = 0;
int n_funcs = funcs.size();
for(int i_func = 0; i_func < n_funcs; i_func++) {
......@@ -91,11 +101,81 @@ namespace genfind {
fHTRoots.push_back( std::make_tuple(i_func, j_func, root) );
avg += root;
n_roots++;
double rhoroot = f1->Eval(root);
//fIntersect->Fill(root,rhoroot);
//fHTIntersections.push_back( std::make_tuple(i_func, j_func, root, rhoroot) );
//std::cout << " root is at " << root << std::endl;
//
//For better effect, make sure that there is a triple coincidence
//only fill fIntersect in the case of at least 3 linear points
fThresh = 3;
for(int k_func = 0; k_func < n_funcs; k_func++) {
if (k_func==i_func || k_func==j_func){continue;}
auto f3 = funcs[k_func];
double checkval = f3->Eval(root);
if (abs(checkval-rhoroot)<fThresh) {
fIntersect->Fill(root,rhoroot);
std::cout << "filled Intersect! " << std::endl;
fHTIntersections.push_back( std::make_tuple(i_func, j_func, root, rhoroot) );
std::cout << "root (theta) is at " << root << " rhoval " << rhoroot << std::endl;
break;
}
}
}
}
TRandom3 rand;
// eventually move this out to another function call
double max = fIntersect->GetMaximumBin();
double value = fIntersect->GetBinContent(max);
std::cout << "maximum bin value is " << value << std::endl;
Int_t binx;
Int_t biny;
Int_t binz;
fIntersect->GetBinXYZ(max, binx, biny, binz);
double xmaxbin = fIntersect->GetXaxis()->GetBinCenter(binx);
double ymaxbin = fIntersect->GetYaxis()->GetBinCenter(biny);
std::cout << "binx " << xmaxbin << "biny " << ymaxbin << std::endl;
//find which hits have functions that go close enough to the maximum bin (generalize to all peaks later)
// helps because the indexing of the functions in funcs is the same as the hit number
// have not solved the case for if the root isn't at the center of the bin and the slope in the bin is high
// something like this should be like "pretrackhits" so that the pretrack hits can be gotten as vectors
fThresh = 1;
std::vector< int > hitlist ;
for(int i_func = 0; i_func < n_funcs; i_func++) {
auto f1 = funcs[i_func];
double val = f1->Eval(xmaxbin);
double dif = ymaxbin-val;
std::cout << "dif is " << dif << std::endl;
if (abs(dif) < fThresh) {
// include this hit in the pretrack for most likely track
//
hitlist.push_back(i_func);
}
}
std::cout << "hitlist size " << hitlist.size() << std::endl;
int hitindex = 0;
for(const auto& ahit : hits) {
//std::cout << "hit index " << hitindex << std::endl;
//fill fRhoTheta0 as placeholder to see "likeliest track"
for(int m = 0; m < hitlist.size(); m++ ){
//std::cout << "hitindex " << hitindex << "hitlist " << hitlist[m] << std::endl;
int filled = 0;
if (hitlist[m]==hitindex) {
fRhoTheta0->Fill(ahit.X(),ahit.Y());
std::cout << "hitindex " << hitindex << "hitlist " << hitlist[m] << std::endl;
filled = 1;
}
if (filled == 1) {break;}
}
hitindex++;
}
/*TRandom3 rand;
TH1F* temp0 = (TH1F*)fRhoTheta0->Clone();
TH1F* temp1 = (TH1F*)fRhoTheta1->Clone();
TH1F* temp2 = (TH1F*)fRhoTheta2->Clone();
......@@ -122,8 +202,10 @@ namespace genfind {
fRhoTheta0->Add(temp0);
fRhoTheta1->Add(temp1);
fRhoTheta2->Add(temp2);
}
}*/
// This makes sense if you are forcing all tracks through 0,0, probably reasonable given vertexing
//std::cout << " average : " << avg/double(n_roots) << std::endl;
//double z_max = fRhoTheta1->GetBinContent(fRhoTheta1->GetMaximumBin());
//double z_mean = fRhoTheta1->GetMean(3);
......@@ -140,7 +222,61 @@ namespace genfind {
return mg;
}
//__________________________________________________________________________
int HoughTransform::FindIntersections( TMultiGraph* mg)
{ // Using fIntercept or the tuple version,
// find peaks corresponding to lines in original image
//
// The next step is to find all the points that correspond to
// rho/theta functions going through the peak points
// Group these points into track seeds.
double max = fIntersect->GetMaximumBin();
double value = fIntersect->GetBinContent(max);
std::cout << "maximum bin value is " << value << std::endl;
Int_t binx;
Int_t biny;
Int_t binz;
fIntersect->GetBinXYZ(max, binx, biny, binz);
double x = fIntersect->GetXaxis()->GetBinCenter(binx);
double y = fIntersect->GetYaxis()->GetBinCenter(biny);
std::cout << "binx " << x << "biny " << y << std::endl;
//use FHTIntersections to find the identity of the hits corresponding to the functions that cross in that bin
//
return max;
/*int size = mg->GetListOfGraphs()->GetSize();// lots of nonsense here accessing multigraphs
TList *funclist = mg->GetListOfFunctions();
TList *graphlist = mg->GetListOfGraphs();
TObjOptLink *lnk = (TObjOptLink*)graphlist->FirstLink();
TObject *obj = lnk->GetObject();
const char * name = obj->GetName();
const char * classname = obj->ClassName();
const char * title = obj->GetTitle();
std::cout << "object name is " << name << std::endl;
std::cout << "class name is " << classname << std::endl;
std::cout << "object title is " << title << std::endl;
lnk = (TObjOptLink*)lnk->Next();
name = obj->GetName();
classname = obj->ClassName();
title = obj->GetTitle();
std::cout << "object name is " << name << std::endl;
std::cout << "class name is " << classname << std::endl;
std::cout << "object title is " << title << std::endl;
return size; */
}
//__________________________________________________________________________
int HoughTransform::FindPeaks(const std::vector<ROOT::Math::XYZTVector>& hits)
{
double degree = TMath::Pi()/180.0;
......
src/GenFind/tests/hough_transform2.cxx
View file @ 4fbde59f
......@@ -24,10 +24,10 @@ void hough_transform2(
HoughTransform* ht = new HoughTransform();
TH2F * hxy = new TH2F("hxy", "hxy", 100, -1000, 1000, 100, -1000, 1000);
TH2F * hxy = new TH2F("hxy", "hxy", 100, -200, 200, 100, -200, 200);
TH2F * hxz = new TH2F("hxz", "hxz", 100, -1000, 1000, 100, -1000, 1000);
TH2F * hyz = new TH2F("hyz", "hyz", 100, -1000, 1000, 100, -1000, 1000);
TH2F * hxy2 = new TH2F("hxy2", "hxy2", 50, -100 , 100, 50, -100, 100);
TH2F * hxy2 = new TH2F("hxy2", "hxy2", 100, -1010 ,1010, 100, -1010, 1010);
TH2F * htrack2 = new TH2F("htrack2", "ht2", 50, -100, 100, 50, -1000, 1000);
TH2F * htrack1 = new TH2F("htrack1", "ht1", 100, -1000, 1000, 100, -1000, 1000);
TH2F * hpeaky = new TH2F("hpeakxy", "hxy", 100, -1000, 1000, 100, -1000, 1000);
......@@ -41,7 +41,7 @@ void hough_transform2(
double degree = TMath::Pi()/180.0;
/*for(int i_track = 0; i_track < Ntracks; i_track++){ // made tracks by combining events
for(int i_track = 0; i_track < Ntracks; i_track++){ // made tracks by combining events
t->GetEntry(i_event);
i_event++;
......@@ -65,21 +65,32 @@ void hough_transform2(
}
}
std::cout << "Created 'event' by combining " << Ntracks << " real events from file\n";
*/
for(int i_hit = 0; i_hit < 100; i_hit++){ // make test hits in straight line
/*for(int i_hit = 0; i_hit < 100; i_hit++){ // make test hits in straight line
double x_incr = 1;
double y_incr = 1;
double z_incr = 1;
hits.push_back( {x_incr*i_hit , y_incr*i_hit, z_incr*i_hit , 0.0} );
hits.push_back( {x_incr*i_hit+10 , y_incr*i_hit, z_incr*i_hit , 0.0} );
std::cout << x_incr*i_hit << " , " << y_incr*i_hit << std::endl;
hxy->Fill(x_incr*i_hit, y_incr*i_hit);
hxz->Fill(x_incr*i_hit, z_incr*i_hit);
hxy->Fill(x_incr*i_hit+10, y_incr*i_hit);
hxz->Fill(x_incr*i_hit+10, z_incr*i_hit);
hyz->Fill(y_incr*i_hit, z_incr*i_hit);
hxy2->Fill(x_incr*i_hit, y_incr*i_hit);
hxy2->Fill(x_incr*i_hit+10, y_incr*i_hit);
}
double x_extra = 0.0;
double y_extra = 50.0;
double z_extra = 0.0;
hits.push_back( {x_extra, y_extra, z_extra, 0.0} );
std::cout << x_extra << " , " << y_extra << std::endl;
hxy->Fill(x_extra, y_extra);
hxz->Fill(x_extra, z_extra);
hyz->Fill(y_extra, z_extra);
hxy2->Fill(x_extra, y_extra);*/
std::vector<std::tuple<XYZTVector,XYZTVector,XYZTVector>> master_hits ;
std::vector<XYZTVector> all_chits ;
......@@ -103,7 +114,7 @@ void hough_transform2(
//for(auto ahit : all_chits) { // find lines in real space first, then conformal space
for(auto ahit : hits) {
//std::cout << ahit.X() << " , " << ahit.Y() << std::endl;
//std::cout << ahit.X() << " , " << ahit.Y() << std::en1l;
huv->Fill( ahit.X(), ahit.Y() );
hphi->Fill(TMath::ATan(ahit.Y()/ahit.X()));
for(int i_theta = 0; i_theta<180; i_theta++) {
......@@ -162,13 +173,16 @@ void hough_transform2(
}
TCanvas * c = new TCanvas();
ht->fIntersect->Draw();
c = new TCanvas();
c->Divide(2,2);
c->cd(1);
hxy->Draw("box");
for(auto ahist : alltrack){
ahist->Draw("box,same");
}
hxy2->Draw("box");
//for(auto ahist : alltrack){
// ahist->Draw("box,same");
//}
//htrack1->SetLineColor(2);
//htrack1->Draw("box,same");
//htrack2->SetLineColor(3);
......@@ -221,16 +235,22 @@ void hough_transform2(
c->cd(3);
mg->Draw("a");
mg->GetYaxis()->SetRangeUser(-0.10,0.10);
std::cout << "multigraph list is " << mg->GetListOfGraphs()->GetSize() << endl;
std::cout << mg->GetListOfGraphs() << endl;
std::cout << mg->GetListOfFunctions() << endl;
c->cd(4);
mg->GetYaxis()->SetRangeUser(-50,50);
mg->Draw("a");
c = new TCanvas();
c->Divide(2,2);
c->cd(1);
ht->fRhoTheta0->SetLineColor(2);
ht->fRhoTheta0->Draw("colz");
c->cd(2);
ht->fRhoTheta1->Draw("colz");
c->cd(3);
ht->fRhoTheta2->Draw("colz");
std::cout << "find intersections " << ht->FindIntersections(mg) << endl;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment