Commit 237dee56 authored by Johnston's avatar Johnston
Browse files

works for lines pretty good

parent d958b1b1
......@@ -50,16 +50,26 @@ namespace genfind {
std::vector<TF1*> FillFunctions(const std::vector<ROOT::Math::XYZTVector>& hits);
int FindIntersections( const std::vector<ROOT::Math::XYZTVector>& hits);
std::vector<std::tuple<int,int,double,double>> FindIntersections( const std::vector<ROOT::Math::XYZTVector>& hits);
int GetReferenceHit( const std::vector<ROOT::Math::XYZTVector>& hits);
int GetReferenceHit( const std::vector<ROOT::Math::XYZTVector>& hits, std::vector<std::tuple<int,int,double,double>> Intersections );
std::vector< int > FindTrack( const std::vector<ROOT::Math::XYZTVector>& hits, int ref );
std::vector< int > FindTrack( const std::vector<ROOT::Math::XYZTVector>& hits, std::vector<std::tuple<int,int,double,double>> Intersections, int ref );
int ReduceIntersections( std::vector< int > trackhitlist );
std::vector<std::tuple<int,int,double,double>> ReduceIntersections(std::vector<std::tuple<int,int,double,double>> Intersections, std::vector< int > trackhitlist );
int FindPeaks(const std::vector<ROOT::Math::XYZTVector>& hits);
std::vector<std::vector<std::tuple<int,double,double>>> FindTracks( const std::vector<ROOT::Math::XYZTVector>& hits );
//return type for find tracks
// innermost, hit element of a tracklist, the hit index, plus the hit:
// std::tuple<int,std::vector<XYZTVector>> typename HTL
// Next, one full tracklist. Vector of HTLs
// std::vector<std::tuple<int,std::vector<XYZTVector>>> typename TL
// next, full findtracks output. Vector of tracklists.
// std::vector<std::vector<std::tuple<int,std::vector<XYZTVector>>>> typename FT
//
// crude idea:
//struct PreTrack {
// double phi_peak;
......
......@@ -80,8 +80,9 @@ namespace genfind {
mg->Add(new TGraph(f),"l");
}
std::cout << "evaluate fFuncs[0](7) " << fFuncs[0]->Eval(7) << std::endl;
return mg;
//fFuncs as member changes in ways I don't understand.
//Not constant, can't be trusted in other class methods.
}
//_______________________________________________________________________________
......@@ -99,43 +100,46 @@ namespace genfind {
f->SetParameter(0, ahit.X());
f->SetParameter(1, ahit.Y());
funcs.push_back(f);
std::cout << "evaluated " << f->Eval(7) << std::endl;
}
std::cout << "evaluate funcs[0](7) " << funcs[0]->Eval(7) << std::endl;
return funcs;
}
//_______________________________________________________________________________
int HoughTransform::FindIntersections( const std::vector<ROOT::Math::XYZTVector>& hits ) {
// for some reason, the returned vector from this is not constant.
// Changes in ways I don't understand.
}
//______________________________________________________________________________
//
// this function finds the crossings for every pair of hit-functions.
// These are only kept in the list if a third hit-function passes very near.
std::vector<std::tuple<int,int,double,double>> HoughTransform::FindIntersections( const std::vector<ROOT::Math::XYZTVector>& hits ) {
//create a list of rho/theta functions from the hit parameters
double degree = TMath::Pi()/180.0;
double phi_min = 0.0;
double phi_max = 180.0;
std::vector<TF1*> funcs;
std::vector<std::tuple<int,int,double,double>> Intersections;
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->SetParameter(0, ahit.X());
f->SetParameter(1, ahit.Y());
funcs.push_back(f);
std::cout << "evaluated " << f->Eval(7) << std::endl;
}
int n_funcs = funcs.size();
//Find the crossings for every pair, by solving roots of difference functions.
int n_roots = 0;
int n_funcs = funcs.size();
for(int i_func = 0; i_func < n_funcs; i_func++) {
for(int j_func = i_func+1; j_func < n_funcs; j_func++) {
auto f1 = funcs[i_func];
if(j_func==i_func+1) {std::cout << "i_func(7) is " << f1->Eval(7) << std::endl; }
auto f2 = funcs[j_func];
TF1 * fdiff = new TF1("fdiff",[&](double*x, double *p){
double val1 = f1->Eval(x[0]);
double val2 = f2->Eval(x[0]);
double dif = val1-val2;
//std::cout << " v1 - v2 = " << val1 << " - " << val2 << " = " << dif << std::endl;
return dif;
}, phi_min, phi_max, 0);
ROOT::Math::WrappedTF1 wf1(*fdiff);
......@@ -144,13 +148,11 @@ namespace genfind {
brf.Solve();
double root = brf.Root();
fHTRoots.push_back( std::make_tuple(i_func, j_func, root) );
//avg += root;
n_roots++;
double rhoroot = f1->Eval(root);
//For better effect, make sure that there is a triple coincidence
//only fill fIntersect in the case of at least 3 nearly 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;}
......@@ -158,151 +160,131 @@ namespace genfind {
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;
Intersections.push_back( std::make_tuple(i_func, j_func, root, rhoroot) );
break;
}
}
}
}
return n_roots;
return Intersections;
}
//_______________________________________________________________________________
//
// This function selects one reference hit, using a histogram
// This hit is in the bin with the most intersections
// It isn't actually a needed function if FindTrack just uses the intersection method itself...
int HoughTransform::GetReferenceHit( const std::vector<ROOT::Math::XYZTVector>& hits, std::vector<std::tuple<int,int,double,double>> Intersections ) {
int HoughTransform::GetReferenceHit( const std::vector<ROOT::Math::XYZTVector>& hits ) {
//recreate the damn funcs vector since it's not stable otherwise...
double degree = TMath::Pi()/180.0;
double phi_min = 0.0;
double phi_max = 180.0;
std::vector<TF1*> funcs;
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->SetParameter(0, ahit.X());
f->SetParameter(1, ahit.Y());
funcs.push_back(f);
std::cout << "evaluated " << f->Eval(7) << std::endl;
}
std::cout << "evaluate funcs[0](7) " << funcs[0]->Eval(7) << std::endl;
// This function selects one reference hit,
// This hit is in the bin with the most intersections
// To allow multiple tracks to be found, make this callable with
// fIntersect as a variable. Or maybe the fHTIntersections tuple var.
int n_funcs = funcs.size();
std::cout << "evaluate funcs[0](7) " << funcs[0]->Eval(7) << std::endl;
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;
std::cout << "evaluate funcs[0](7) again " << funcs[0]->Eval(7) << std::endl;
//fIntersect->GetBinXYZ(max, binx, biny, binz); some bizarre reason this fucks up fFuncs
std::cout << "evaluate funcs[0](7) " << funcs[0]->Eval(7) << std::endl;
double xmaxbin = 7; //fIntersect->GetXaxis()->GetBinCenter(binx);
double ymaxbin = 4; //fIntersect->GetYaxis()->GetBinCenter(biny);
//find the location of the most intersections, using hist bins.
TH2F* hist_intersect = new TH2F{"hist_intersect","hist_intersect",90,0, 180, 75,-150,150 };
for(auto t : Intersections ){
hist_intersect->Fill(std::get<2>(t),std::get<3>(t));
}
double max = hist_intersect->GetMaximumBin();
double value = hist_intersect->GetBinContent(max);
Int_t binx,biny,binz;
hist_intersect->GetBinXYZ(max, binx, biny, binz);
double xmaxbin = hist_intersect->GetXaxis()->GetBinCenter(binx);
double ymaxbin = hist_intersect->GetYaxis()->GetBinCenter(biny);
std::cout << "binx " << xmaxbin << "biny " << ymaxbin << std::endl;
std::cout << "evaluate funcs[0](7) again " << funcs[0]->Eval(7) << 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; // instead, probably look for the minimum difference rather than set threshold.
std::vector< int > hitlist ;
std::cout << "evaluate funcs[0](7) " << funcs[0]->Eval(7) << std::endl;
delete hist_intersect;
//find one hit with function that goes closest to the maximum bin (
int min = abs(funcs[0]->Eval(xmaxbin)-ymaxbin);
int index = -1;
for(int i = 0; i < n_funcs; i++) {
auto f1 = funcs[i];
std::cout << "evaluate funcs[0](7) " << funcs[0]->Eval(7) << std::endl;
std::cout << "i_func(7) is " << f1->Eval(7) << std::endl;
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);
if (abs(dif) < min) {
min = abs(dif);
index = i;
}
}
std::cout << "hitlist size " << hitlist.size() << std::endl;
return hitlist[0];
std::cout << "reference " << index << std::endl;
return index;
}
//_____________________________________________________________________________
//
//This finds a list of candidate hits constituting a single track.
//It requires a list of hits and intersections, created by HoughTransform::FindIntersections
//It also requires a reference hit. Suitable hit found by HoughTransform::GetReferenceHit
//It doesn't actually need to require a reference hit, just does right now.
//The track is the longest candidate track in the list of hit intersections given
std::vector< int > HoughTransform::FindTrack( const std::vector<ROOT::Math::XYZTVector>& hits, std::vector<std::tuple<int,int,double,double>> Intersections, int ref ) {
std::vector< int > HoughTransform::FindTrack( const std::vector<ROOT::Math::XYZTVector>& hits, int ref ) {
//recreate the damn funcs vector since it's not stable otherwise...
double degree = TMath::Pi()/180.0;
double phi_min = 0.0;
double phi_max = 180.0;
std::vector<TF1*> funcs;
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->SetParameter(0, ahit.X());
f->SetParameter(1, ahit.Y());
funcs.push_back(f);
std::cout << "evaluated " << f->Eval(7) << std::endl;
}
int n_funcs = funcs.size();
//use fHTIntersections to find candidate tracks.
std::vector<std::tuple<int,int,double,double>> fHTIntReduced;
std::vector<std::tuple<int,int,double,double>> fHTCandidate;
std::cout << "reference " << ref << std::endl;
for(auto t:fHTIntersections) {
int index1 = std::get<0>(t);
int index2 = std::get<1>(t);
//std::cout << "index 1 " << index1 << " and 2 " << index2 << std::endl;
if( (index1 != ref) && (index2 != ref) ){
// neither index matches the candidate hit, enter this into the fHTIntReduced;
// I don't want to fill it at this point, because I haven't selected the densest point for this hit
// Need to remove from this listing ALL points corresponding to the max grouping for this hit
// THEN the remainder is ready for continued searching
//std::cout << "filling reduced with " << std::get<2>(t) << " " << std::get<3>(t) << std::endl;
fHTIntReduced.push_back(t);
//fIntReduced->Fill(std::get<2>(t),std::get<3>(t));
} else {
// one of the indexes matches the candidate hit-- fill fHTCandidate
fHTCandidate.push_back(t);
fIntReduced->Fill(std::get<2>(t),std::get<3>(t));
std::cout << "filling candidate with " << std::get<2>(t) << " " << std::get<3>(t) << std::endl;
}
//check that the ref value isn't the -1 error value...
if(ref<0){
std::cout << "The reference is wrong " << std::endl;
//also, this function doesn't actually use reference. rethink.
//It uses the maximum histogram bin position instead
}
int reducedmaxbin = fIntReduced->GetMaximumBin();
Int_t binxR;
Int_t binyR;
Int_t binzR;
fIntReduced->GetBinXYZ(reducedmaxbin, binxR, binyR, binzR);
double xmaxbinR = fIntersect->GetXaxis()->GetBinCenter(binxR);
double ymaxbinR = fIntersect->GetYaxis()->GetBinCenter(binyR);
std::cout << "evaluate funcs[0](7) " << funcs[0]->Eval(7) << std::endl;
std::cout << "binxR " << xmaxbinR << "binyR " << ymaxbinR << std::endl;
// use the fIntReduced, which now holds candidate points,
//find the location of the most intersections, using hist bins.
TH2F* hist_intersect = new TH2F{"hist_intersect","hist_intersect",90,0, 180, 75,-150,150 };
for(auto t : Intersections ){
hist_intersect->Fill(std::get<2>(t),std::get<3>(t));
}
double max = hist_intersect->GetMaximumBin();
double value = hist_intersect->GetBinContent(max);
Int_t binx,biny,binz;
hist_intersect->GetBinXYZ(max, binx, biny, binz);
double xmaxbin = hist_intersect->GetXaxis()->GetBinCenter(binx);
double ymaxbin = hist_intersect->GetYaxis()->GetBinCenter(biny);
std::cout << "binx " << xmaxbin << "biny " << ymaxbin << std::endl;
delete hist_intersect;
// find all points within a fairly relaxed threshold around max bin
// Accumulate list of all the hits that make up this group
// Then make the fHTIntReduced vector/tuple of the remaining hits.
fThresh = 15;
std::vector< int > trackhitlist;
for(int i_func = 0; i_func < n_funcs; i_func++) {
auto f1 = funcs[i_func];
double val = f1->Eval(7); //xmaxbinR);
std::cout << "val is " << val << std::endl;
std::cout << "evaluate funcs[0](7) " << funcs[0]->Eval(7) << std::endl;
// first time, the above line outputs sensible. Every subsequent loop,
// outputs fucking insanity.
double dif = ymaxbinR-val;
double val = f1->Eval(xmaxbin);
double dif = ymaxbin-val;
std::cout << "index: " << i_func << " dif is " << dif << std::endl;
if (abs(dif) < fThresh) {
// include this hit in the pretrack for most likely track
//
trackhitlist.push_back(i_func);
std::cout << "filling trackhitlist with " << i_func << std::endl;
}
......@@ -313,8 +295,13 @@ namespace genfind {
//__________________________________________________________________________________
//
int HoughTransform::ReduceIntersections( std::vector< int > trackhitlist ) {
for(auto t:fHTIntersections) {
//Taking an initial set of function intersections, and a candidate track list,
//This function checks each intersection and culls it if either hit is on the track list.
std::vector<std::tuple<int,int,double,double>> HoughTransform::ReduceIntersections( std::vector<std::tuple<int,int,double,double>> Intersections, std::vector< int > trackhitlist ) {
std::vector<std::tuple<int,int,double,double>> ReducedInt;
for(auto t:Intersections) {
int index1 = std::get<0>(t);
int index2 = std::get<1>(t);
//std::cout << "index " << index1 << " and " << index2 << std::endl;
......@@ -322,18 +309,51 @@ namespace genfind {
for( int i = 0; i<trackhitlist.size(); i++ ) {
if ( (index1 == trackhitlist[i]) || (index2 == trackhitlist[i]) ){
addflag = 0;
std::cout << trackhitlist[i] << std::endl;
//std::cout << trackhitlist[i] << std::endl;
}
}
if ( addflag ==1 ) {
fIntLeft->Fill(std::get<2>(t),std::get<3>(t));
//std::cout << "filling tuple " << std::endl;
ReducedInt.push_back(t);
std::cout << "filling reduced tuple " << std::get<0>(t) << " " << std::get<1>(t) << std::endl;
}
}
return 0;
}
return ReducedInt;
}
//__________________________________________________________________________
//
//should return vector of vector of tuple of int and vector of xyztvectors
//should continue calling until all tracks are found.
std::vector<std::vector<std::tuple<int,double,double>>> HoughTransform::FindTracks( const std::vector<ROOT::Math::XYZTVector>& hits ){
int num_tracks = 0;
int max_tracks = 3;
std::vector<std::vector<std::tuple<int,double,double>>> FoundTracks;
std::vector<std::tuple<int,int,double,double>> fi = this->FindIntersections(hits);
while (fi.size()>0&&num_tracks<max_tracks){
int gr = this->GetReferenceHit(hits,fi);
std::cout<< "intersections have tracks # " << fi.size() << std::endl;
std::vector< int > trackhitlist = this->FindTrack( hits, fi, gr );
std::vector<std::tuple<int,double,double>> TL;
for(int i=0; i<trackhitlist.size(); i++){
int tracknumber = trackhitlist[i];
std::tuple<int,double,double> HTL = std::make_tuple(tracknumber,hits[tracknumber].X(),hits[tracknumber].Y());
//push back std::tuple<int, HIT> ? or some other more clever solution here
TL.push_back(HTL);
}
FoundTracks.push_back(TL);
num_tracks++;
fi = this->ReduceIntersections( fi, trackhitlist );
std::cout<< " intersections size " << fi.size() << std::endl;
}
//return num_tracks;
return FoundTracks;
}
//__________________________________________________________________________
//
......
......@@ -126,19 +126,34 @@ void hough_transform2(
auto mg = ht->FillHoughTransforms(hits);
std::vector<TF1*> ff = ht->FillFunctions(hits);
int fi = ht->FindIntersections(hits);
/*int fi = ht->FindIntersections(hits);
int gr = ht->GetReferenceHit(hits);
std::vector< int > trackhitlist = ht->FindTrack( hits,gr );
int ri = ht->ReduceIntersections( trackhitlist );
*/
//int numtracks = ht->FindTracks(hits);
std::vector<std::vector<std::tuple<int,double,double>>> FoundTracks = ht->FindTracks(hits);
std::cout << "FindTracks found " << FoundTracks.size() << " tracks " << endl;
std::cout << " n peaks : " << ht->FindPeaks(hits) << std::endl;
std::vector<TH1F*> FoundTrackHists;
int tracknumber=0;
for( std::vector<std::tuple<int,double,double>> atrack : FoundTracks ) {
tracknumber++;
TH1F* singletrack = (TH1F*)htrack1->Clone();
singletrack->Reset();
singletrack->SetLineColor(2+FoundTrackHists.size());
for( std::tuple<int,double,double> thit : atrack ) {
std::cout << "track number " << tracknumber << " hit number " << std::get<0>(thit) << std::endl;
std::cout << " x value " << std::get<1>(thit) << " y value " << std::get<2>(thit) << std::endl;
singletrack->Fill( std::get<1>(thit), std::get<2>(thit) );
}
FoundTrackHists.push_back(singletrack);
}
//std::cout << " n peaks : " << ht->FindPeaks(all_chits) << std::endl;
auto pretrack_hits = ht->GetPreTrackHits(master_hits,4.0);
std::vector<TH1F*> hists;
//auto pretrack_hits = ht->GetPreTrackHits(master_hits,4.0);
/*std::vector<TH1F*> hists;
std::vector<TH1F*> alltrack;
std::vector<TH1F*> vertextrack;
for( auto atrack : pretrack_hits ) {
......@@ -161,7 +176,8 @@ void hough_transform2(
hists.push_back(fPhi2);
alltrack.push_back(singletrack);
vertextrack.push_back(singlevertextrack);
}
}*/
// auto trackone = pretrack_hits.at(0);
// for( auto thit : trackone ) {
// htrack1->Fill( std::get<0>(thit).X(), std::get<0>(thit).Y());
......@@ -178,18 +194,30 @@ void hough_transform2(
}
TCanvas * c = new TCanvas();
ht->fIntersect->Draw();
//ht->fIntersect->Draw();
c = new TCanvas();
ht->fIntReduced->Draw();
c->Divide(2,2);
c->cd(1);
mg->Draw("a");
std::cout << "multigraph list is " << mg->GetListOfGraphs()->GetSize() << endl;
std::cout << mg->GetListOfGraphs() << endl;
std::cout << mg->GetListOfFunctions() << endl;
c->cd(2);
mg->GetYaxis()->SetRangeUser(-50,50);
mg->Draw("a");
//c = new TCanvas();
//ht->fIntReduced->Draw();
//c = new TCanvas();
//ht->fIntLeft->Draw();
c = new TCanvas();
ht->fIntLeft->Draw();
c = new TCanvas();
c->Divide(2,2);
c->cd(1);
hxy2->Draw("box");
//for(auto ahist : alltrack){
//for(auto ahist : FoundTrackHists){
// ahist->Draw("box,same");
//}
//htrack1->SetLineColor(2);
......@@ -198,9 +226,11 @@ void hough_transform2(
//htrack2->Draw("box,same");
//hpeaky->Draw("box,same");
c->cd(2);
huv->Draw("box");
//huv->Draw("box");
if(FoundTrackHists.size()>0) {FoundTrackHists[0]->Draw("box");}
c->cd(3);
hrphi->Draw("");
if(FoundTrackHists.size()>1) {FoundTrackHists[1]->Draw("box");}
//hrphi->Draw("");
/*hxy2->Draw("box");
for(auto ahist : vertextrack){
ahist->Draw("box,same");
......@@ -208,13 +238,14 @@ void hough_transform2(
//hpeaky2->SetLineColor(2);
//hpeaky2->Draw("box,same");
c->cd(4);
hphi->Draw();
if(FoundTrackHists.size()>2) {FoundTrackHists[2]->Draw("box");}
//hphi->Draw();
//hrphi->Draw("lego2");
c = new TCanvas();
/*c = new TCanvas();
c->Divide(2,2);
int counter = 1;
for( auto ahist : alltrack ){
for( auto ahist : FoundTrackHists ){
c->cd(counter);
ahist->Draw("box");
counter++;
......@@ -229,8 +260,8 @@ void hough_transform2(
hxz->Draw("box");
c->cd(3);
hyz->Draw("box");
*/
/*
c = new TCanvas();
c->Divide(2,2);
c->cd(1);
......@@ -241,17 +272,9 @@ void hough_transform2(
//ht->fRhoTheta0->Draw("colz");
c->cd(2);
htheta->Draw();
*/
c->cd(3);
mg->Draw("a");
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 = new TCanvas();
c->Divide(2,2);
c->cd(1);
ht->fRhoTheta0->SetLineColor(2);
......@@ -260,5 +283,5 @@ void hough_transform2(
ht->fRhoTheta1->Draw("colz");
c->cd(3);
ht->fRhoTheta2->Draw("colz");
*/
}
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