Commit d8a95009 authored by Johnston's avatar Johnston
Browse files

work on the conformal map method

 Changes to be committed:
	modified:   ../include/HoughTransform.h
	modified:   HoughTransform.cxx
	modified:   ../tests/hough_transform2.cxx
parent 36c32766
......@@ -46,6 +46,9 @@ namespace genfind {
const ROOT::Math::XYZTVector ref_point,
const std::vector<ROOT::Math::XYZTVector>& hits) const;
std::vector<double> UseConformalMap( double a, double b);
//const::std::vector<ROOT::Math::XYZTVector>& hits );
TMultiGraph* FillHoughTransforms(const std::vector<ROOT::Math::XYZTVector>& hits);
std::vector<TF1*> FillFunctions(const std::vector<ROOT::Math::XYZTVector>& hits);
......
......@@ -42,23 +42,156 @@ namespace genfind {
if(fSpec1) delete fSpec1; fSpec1 = nullptr;
}
//__________________________________________________________________________
//
// the conformal map, 1/z, for complex value z, transforms a circle into a line
// However, it only does this for a circle positioned on (radius,0).
// The method finds the center and radius for a hit and a reference hit,
// assuming the track segment comes from the vertex
//
// the conformal map is only saved if there exists a fourth hit, bhit,
// which also lies on the defined circle.
//
// If the points are such that the radius is too large, the points are checked for linearity.
std::vector<ROOT::Math::XYZTVector> HoughTransform::GetConformalCoordinates(
const ROOT::Math::XYZTVector ref_point,
const std::vector<ROOT::Math::XYZTVector>& hits) const
{
std::vector<ROOT::Math::XYZTVector> res;
std::for_each(hits.begin(), hits.end(), [&](const auto& ahit) {
double x = ahit.X() - ref_point.X();
double y = ahit.Y() - ref_point.Y();
double R = x*x + y*y;
//std::for_each(hits.begin(), hits.end(), [&](const auto& ahit) {
for(const auto& ahit : hits) {
for(const auto & bhit : hits) {
// we are assuming an x0,y0 that is the vertex.
double x1 = ahit.X();
double x2 = ref_point.X();
double x3 = bhit.X();
double y1 = ahit.Y();
double y2 = ref_point.Y();
double y3 = bhit.Y();
if( !(ahit == ref_point) ) {
res.push_back(ROOT::Math::XYZTVector{ 2.0*x/R, -2.0*y/R, ahit.Z(), ahit.T() });
//find the center and radius of the circle described by the three points, (p,q), and rad
double mag1 = x1*x1+y1*y1;
double mag2 = x2*x2+y2*y2;
double mag3 = x3*x3+y3*y3;
//double q = (mag3-mag2)*(x2-x1)-(mag2-mag1)*(x3-x2)/( 2*( (y1-y2)*(x3-x2)-(y2-y3)*(x2-x1) ) );
//double p = (mag3-mag2+2*q*(y2-y3))/( 2*(x3-x2) );
//double rad = std::sqrt( (x3-p)*(x3-p) + (y3-q)*(y3-q) );
//the following assume the track goes through the vertex
double divide = 2*(x2*y1-x1*y2);
double p=0;
double q=0;
double rad=0;
double check = 100;
bool linear = false;
if ( divide > .001 ){
p = -( mag1*y2-mag2*y1 ) / divide;
q = ( mag1*x2-mag2*x1 ) / divide;
rad = std::sqrt(p*p+q*q);
std::cout << "p is " << p << " q is " << q << " r is " << rad << std::endl;
// confirm that there exists an additional hit, bhit, which follows this equation
check = (x3-p)*(x3-p)+(y3-p)*(y3-p);
if ( abs(check - rad*rad) > .5*rad ) {
//this combination of hits is not worth adding to the return set, continue and check another bhit.
continue;
}
}else{
// vertex, ahit, and refhit did not give a proper equation.
// see if instead of a circle, these points are a line
double m;
double b;
if ( (x2-x1) > .001 ) {
m = (y2-y1)/(x2-x1);
b = y1=x1*m;
// ok, have an equation for a line, now check if either x3 or vertex is on it.
check = abs(y3-(m*x3+b));
if ( check < .05*fThresh ) {
linear = true;
}else {
if ( abs(b) < .05*fThresh ) {
linear = true;
}
}
} else {
//this might be a vertical line.
if ( (x2-x1) < .05*fThresh ){
// check if either x3 or vertex is vertical with x2/x1
if ( (x3-(x2+x1)/2 ) < .05*fThresh ) {
linear = true;
if ( (x2+x1)/2 < .05*fThresh ) {
linear = true;
}
}
}
}
}
// if linear = true, don't use the conformal map 1/z, just include the points as is.
if ( linear == false ) {
// this conformal map 1/z in complex plane,
// only works to transform a circle with a center on (r,0)
// where r is the radius of that circle. Transform point 1 to that position
double x_trans = x1-p+rad;
double y_trans = y1-q;
double R = x_trans*x_trans + y_trans*y_trans;
res.push_back(ROOT::Math::XYZTVector{ x_trans/R, y_trans/R, ahit.Z(), ahit.T() });
//std::cout << "ref hit at (" << x2 << "," << y2 << ")" << std::endl;
//std::cout << "hit is at (" << x1 << "," << y1 << ")" << std::endl;
std::cout << "center is at (" << p << "," << q << ")" << std::endl;
std::cout << "radius is at " << rad << std::endl;
std::cout << std::endl;
} else if ( linear == true ) {
// linear case is sensible enough, just include untransformed points.
res.push_back(ROOT::Math::XYZTVector{ x1, y1, ahit.Z(), ahit.T() });
std::cout << " 3 points are linear, not transforming these coordinates " << std::endl;
}
}
});
}
};
return res;
}
//__________________________________________________________________________
//_________________________________________________________________________
//
//GetConformalCoordinates does not actually transform a circle into a line at all
//Trying another equation to test.
//
//std::vector<ROOT::Math::XYZTVector> HoughTransform::UseConformalMap(
std::vector<double> HoughTransform::UseConformalMap(
// const::std::vector<ROOT::Math::XYZTVector>& hits
double a, double b
)
{
std::vector<double> res;
//std::vector<ROOT::Math::XYZTVector> res;
//std::for_each(hits.begin(), hits.end(), [&](const auto& ahit) {
double x = a/(a*a+b*b);
double y = -b/(a*a+b*b);
std::cout << "transformed x is " << x << " transformed y is " << y << std::endl;
res.push_back(x);
res.push_back(y);
//});
//test some basic values
return res;
}
//_________________________________________________________________________
//
TMultiGraph* HoughTransform::FillHoughTransforms(const std::vector<ROOT::Math::XYZTVector>& hits)
{
......
......@@ -97,17 +97,18 @@ void hough_transform2(
//make track in a circle
int numhits = 20;
double radius = 1000;
double centerx = -1000;
int numhits = 7;
double radius = 3;
double centerx = 3;
double centery = 0;
for(int i_hit=0; i_hit < numhits; i_hit++){
double theta = 30./numhits*degree*i_hit;
double theta = 360./numhits*degree*i_hit+5;
double xval = TMath::Cos(theta)*radius+centerx;
double yval = TMath::Sin(theta)*radius+centery;
hits.push_back( {xval , yval , 0 , 0} );
std::cout << " x val " << xval << " y val " << yval << std::endl;
std::cout << " x val " << xval << " y val " << yval << std::endl;
hxy->Fill(xval,yval);
hxz->Fill(xval,0);
hyz->Fill(yval,0);
......@@ -118,6 +119,8 @@ void hough_transform2(
std::vector<std::tuple<XYZTVector,XYZTVector,XYZTVector>> master_hits ;
std::vector<XYZTVector> all_chits ;
for(const auto& fhit : hits){
//auto fhit = hits.at(0);
auto chits = ht->GetConformalCoordinates(fhit, hits);
......@@ -138,7 +141,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::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++) {
......@@ -156,7 +159,12 @@ void hough_transform2(
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;
......@@ -174,6 +182,8 @@ void hough_transform2(
}
FoundTrackHists.push_back(singletrack);
}
*/
//std::cout << " n peaks : " << ht->FindPeaks(all_chits) << std::endl;
//auto pretrack_hits = ht->GetPreTrackHits(master_hits,4.0);
......@@ -249,12 +259,16 @@ void hough_transform2(
//htrack2->SetLineColor(3);
//htrack2->Draw("box,same");
//hpeaky->Draw("box,same");
/*
c->cd(2);
if(FoundTrackHists.size()>0) {FoundTrackHists[0]->Draw("box");}
c->cd(3);
if(FoundTrackHists.size()>1) {FoundTrackHists[1]->Draw("box");}
c->cd(4);
if(FoundTrackHists.size()>2) {FoundTrackHists[2]->Draw("box");}
*/
/*c = new TCanvas();
c->Divide(2,2);
......
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