Commit d958b1b1 authored by Johnston's avatar Johnston
Browse files

separated functions

parent 46ffbc01
......@@ -32,6 +32,7 @@ namespace genfind {
int fThresh = 3;
int fMaxTracks = 20;
std::vector<TF1*> fFuncs;
std::vector<std::tuple<int,int,double>> fHTRoots;
std::vector<std::tuple<int,int,double,double>> fHTIntersections;
std::vector<double> fPhiPeaks;
......@@ -47,7 +48,15 @@ namespace genfind {
TMultiGraph* FillHoughTransforms(const std::vector<ROOT::Math::XYZTVector>& hits);
int FindIntersections( TMultiGraph* mg );
std::vector<TF1*> FillFunctions(const std::vector<ROOT::Math::XYZTVector>& hits);
int FindIntersections( const std::vector<ROOT::Math::XYZTVector>& hits);
int GetReferenceHit( const std::vector<ROOT::Math::XYZTVector>& hits);
std::vector< int > FindTrack( const std::vector<ROOT::Math::XYZTVector>& hits, int ref );
int ReduceIntersections( std::vector< int > trackhitlist );
int FindPeaks(const std::vector<ROOT::Math::XYZTVector>& hits);
......
......@@ -66,7 +66,7 @@ namespace genfind {
double phi_min = 0.0;
double phi_max = 180.0;
TMultiGraph* mg = new TMultiGraph();
std::vector<TF1*> funcs;
//std::vector<TF1*> funcs;
double avg = 0;
......@@ -76,17 +76,60 @@ namespace genfind {
f->SetLineWidth(1);
f->SetParameter(0, ahit.X());
f->SetParameter(1, ahit.Y());
funcs.push_back(f);
fFuncs.push_back(f);
mg->Add(new TGraph(f),"l");
}
std::cout << "evaluate fFuncs[0](7) " << fFuncs[0]->Eval(7) << std::endl;
return mg;
}
//_______________________________________________________________________________
std::vector<TF1*> HoughTransform::FillFunctions(const std::vector<ROOT::Math::XYZTVector>& hits)
{
double degree = TMath::Pi()/180.0;
double phi_min = 0.0;
double phi_max = 180.0;
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->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 ) {
// 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.
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;
}
// 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++) {
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]);
......@@ -101,16 +144,13 @@ namespace genfind {
brf.Solve();
double root = brf.Root();
fHTRoots.push_back( std::make_tuple(i_func, j_func, root) );
avg += 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
//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;}
......@@ -118,16 +158,40 @@ namespace genfind {
double checkval = f3->Eval(root);
if (abs(checkval-rhoroot)<fThresh) {
fIntersect->Fill(root,rhoroot);
std::cout << "filled Intersect! " << std::endl;
//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;
//std::cout << "root (theta) is at " << root << " rhoval " << rhoroot << std::endl;
break;
}
}
}
}
return n_roots;
}
//_______________________________________________________________________________
//
// eventually move this out to another function call
int HoughTransform::GetReferenceHit( const std::vector<ROOT::Math::XYZTVector>& hits ) {
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;
......@@ -135,63 +199,69 @@ namespace genfind {
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 << "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);
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;
fThresh = 1; // instead, probably look for the minimum difference rather than set threshold.
std::vector< int > hitlist ;
for(int i_func = 0; i_func < n_funcs; i_func++) {
auto f1 = funcs[i_func];
std::cout << "evaluate funcs[0](7) " << funcs[0]->Eval(7) << std::endl;
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_func);
hitlist.push_back(i);
}
}
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()); // I think I should just do hit[proper index]
std::cout << "hitindex " << hitindex << " hitlist " << hitlist[m] << std::endl;
filled = 1;
}
if (filled == 1) {break;}
}
hitindex++;
}
return hitlist[0];
}
//_____________________________________________________________________________
//
std::vector< int > HoughTransform::FindTrack( const std::vector<ROOT::Math::XYZTVector>& hits, int ref ) {
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.
//TH2F* fIntReduced = new TH2F{"hIntReduced","hIntReduced",90,0, 180, 75,-150,150 };
std::vector<std::tuple<int,int,double,double>> fHTIntReduced;
std::vector<std::tuple<int,int,double,double>> fHTCandidate;
//for(int i_func = 0; i_func < n_funcs; i_func++) {
//first element of hitlist should be in the maximum bin...
//this is the index I'm interested in at the moment
std::cout << "hitlist first element " << hitlist[0] << std::endl;
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 != hitlist[0]) && (index2 != hitlist[0]) ){
//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
......@@ -204,31 +274,30 @@ namespace genfind {
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;
}
}
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 << "binxR " << xmaxbin << "binyR " << ymaxbin << std::endl;
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 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(xmaxbinR);
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;
std::cout << "index: " << i_func << " dif is " << dif << std::endl;
if (abs(dif) < fThresh) {
......@@ -238,12 +307,17 @@ namespace genfind {
std::cout << "filling trackhitlist with " << i_func << std::endl;
}
}
return trackhitlist;
}
//__________________________________________________________________________________
//
int HoughTransform::ReduceIntersections( std::vector< int > trackhitlist ) {
for(auto t:fHTIntersections) {
int index1 = std::get<0>(t);
int index2 = std::get<1>(t);
std::cout << "index " << index1 << " and " << index2 << std::endl;
//std::cout << "index " << index1 << " and " << index2 << std::endl;
int addflag = 1;
for( int i = 0; i<trackhitlist.size(); i++ ) {
if ( (index1 == trackhitlist[i]) || (index2 == trackhitlist[i]) ){
......@@ -253,116 +327,16 @@ namespace genfind {
}
if ( addflag ==1 ) {
fIntLeft->Fill(std::get<2>(t),std::get<3>(t));
std::cout << "filling tuple " << std::endl;
//std::cout << "filling tuple " << std::endl;
}
}
//find densest group in fHTCandidate
/*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);
}*/
// 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);
//std::cout << z_max << std::endl;
//std::cout << z_mean << std::endl;
//double aZ_max = fRhoTheta1->GetMaximum(z_max);
//fThresh = 0.9*z_max;
//int bin = fRhoTheta1->FindFirstBinAbove(fThresh,2);
//while( bin > 0 ) {
// std::cout << fRhoTheta1->GetBinContent(bin) << std::endl;
// //aZ_max = fRhoTheta1->GetMaximum(aZ_max-0.0001);
// bin = fRhoTheta1->FindFirstBinAbove(fThresh);
//}
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; */
return 0;
}
//__________________________________________________________________________
//
int HoughTransform::FindPeaks(const std::vector<ROOT::Math::XYZTVector>& hits)
{
double degree = TMath::Pi()/180.0;
......
......@@ -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,9 +65,9 @@ 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;
......@@ -89,7 +89,7 @@ void hough_transform2(
hxy->Fill(x_extra, y_extra);
hxz->Fill(x_extra, z_extra);
hyz->Fill(y_extra, z_extra);
hxy2->Fill(x_extra, y_extra);
hxy2->Fill(x_extra, y_extra);*/
std::vector<std::tuple<XYZTVector,XYZTVector,XYZTVector>> master_hits ;
std::vector<XYZTVector> all_chits ;
......@@ -125,6 +125,11 @@ void hough_transform2(
//auto mg = ht->FillHoughTransforms(all_chits);
auto mg = ht->FillHoughTransforms(hits);
std::vector<TF1*> ff = ht->FillFunctions(hits);
int fi = ht->FindIntersections(hits);
int gr = ht->GetReferenceHit(hits);
std::vector< int > trackhitlist = ht->FindTrack( hits,gr );
int ri = ht->ReduceIntersections( trackhitlist );
std::cout << " n peaks : " << ht->FindPeaks(hits) << std::endl;
......@@ -256,5 +261,4 @@ void hough_transform2(
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