Commit dbace9d3 authored by Johnston's avatar Johnston
Browse files

testing conformal method

parent d8a95009
......@@ -46,8 +46,8 @@ 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 );
std::vector<ROOT::Math::XYZTVector> UseConformalMap(
const::std::vector<ROOT::Math::XYZTVector>& hits) const;
TMultiGraph* FillHoughTransforms(const std::vector<ROOT::Math::XYZTVector>& hits);
......
......@@ -165,29 +165,19 @@ namespace genfind {
}
//_________________________________________________________________________
//
//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
//std::vector<double> HoughTransform::UseConformalMap(
std::vector<ROOT::Math::XYZTVector> HoughTransform::UseConformalMap(
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();
double y = ahit.Y();
double R = x*x + y*y;
res.push_back(ROOT::Math::XYZTVector{ 2.0*x/R, -2.0*y/R, ahit.Z(), ahit.T() });
});
return res;
}
//_________________________________________________________________________
......@@ -286,7 +276,8 @@ namespace genfind {
double rhoroot = f1->Eval(root);
//For better effect, make sure that there is a triple coincidence
fThresh = 3;
//fThresh = 3; // this works well for regular, not conformal hits
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];
......@@ -409,6 +400,7 @@ namespace genfind {
// find all points within a fairly relaxed threshold around max bin
// Accumulate list of all the hits that make up this group
//fThresh = 15; // this works for regular, not conformal hits
fThresh = 15;
std::vector< int > trackhitlist;
for(int i_func = 0; i_func < n_funcs; i_func++) {
......
void HT_test(
int i_event = 20,
int Ntracks = 3)
{
// test the efficiency of separating hits by track. when reconstructing, also fill a "truth" vector
// which holds which track each hit is supposed to go to. Compare against the result, check performance
//
gSystem->Load("libGenFind.so");
using namespace genfind;
using namespace ROOT::Math;
TFile * f = new TFile("simple_example_out.root","READ");
TTree * t = (TTree*)gROOT->FindObject("EVENT");
if(!t) { std::cout << " Tree not found " << std::endl; return;}
int nentries=t->GetEntries();
if(i_event+Ntracks>nentries){
std::cout << "Only " << nentries << " in file, this parameter set fails" << std::endl;
return;
}
std::vector<DD4hep::Simulation::Geant4Tracker::Hit*> * g4hits = nullptr;
std::vector<DD4hep::Simulation::Geant4Tracker::Hit*> * g4hits2 = nullptr;
t->SetBranchAddress("SiVertexBarrelHits", &g4hits);
t->SetBranchAddress("SiTrackerBarrelHits", &g4hits2);
HoughTransform* ht = new HoughTransform();
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", 100, -1010 ,1010, 100, -1010, 1010);
TH2F * huv = new TH2F("huv", "huv", 100, -.02 , .02, 100, -.02, .02);
//make composite tracks with real datapoints
std::vector<XYZTVector> hits ;
std::vector<int> tracktruth ;
double degree = TMath::Pi()/180.0;
for(int i_track = 0; i_track < Ntracks; i_track++){ // made tracks by combining events
t->GetEntry(i_event);
i_event++;
for( auto thit : (*g4hits) ) {
hits.push_back( {thit->position.X(), thit->position.Y(), thit->position.Z(), 0.0} );
tracktruth.push_back(i_track);
std::cout << thit->position.X() << " , " << thit->position.Y() << std::endl;
hxy->Fill(thit->position.X(), thit->position.Y());
hxz->Fill(thit->position.X(), thit->position.Z());
hyz->Fill(thit->position.Y(), thit->position.Z());
hxy2->Fill(thit->position.X(), thit->position.Y());
}
for( auto thit : (*g4hits2) ) {
hits.push_back( {thit->position.X(), thit->position.Y(), thit->position.Z(), 0.0} );
tracktruth.push_back(i_track);
std::cout << thit->position.X() << " , " << thit->position.Y() << std::endl;
hxy->Fill(thit->position.X(), thit->position.Y());
hxz->Fill(thit->position.X(), thit->position.Z());
hyz->Fill(thit->position.Y(), thit->position.Z());
hxy2->Fill(thit->position.X(), thit->position.Y());
}
}
std::cout << "Created 'event' by combining " << Ntracks << " real events from file\n";
//make track with artificial linear points.
/*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+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+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+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);*/
auto mg = ht->FillHoughTransforms(hits);
// use track truth and the tracknumber in the findtracks result to test performance.
// First, are all the hits properly separated,
// Second, are all the hits accounted for
std::vector<ROOT::Math::XYZTVector> conf_hits = ht->UseConformalMap(hits);
for (auto conf_hit : conf_hits ){
huv->Fill(conf_hit.X(),conf_hit.Y()) ;
}
std::vector<std::vector<std::tuple<int,double,double>>> FoundTracks = ht->FindTracks(conf_hits);
//std::vector<std::vector<std::tuple<int,double,double>>> FoundTracks = ht->FindTracks(hits);
std::cout << "FindTracks found " << FoundTracks.size() << " tracks " << endl;
std::vector<TH1F*> FoundTrackHists;
int tracknumber=0;
for( std::vector<std::tuple<int,double,double>> atrack : FoundTracks ) {
tracknumber++;
TH1F* singletrack = (TH1F*)hxy2->Clone();
singletrack->Reset();
singletrack->SetLineColor(2+FoundTrackHists.size());
int ref_truth_track = tracktruth[std::get<0>(atrack[0])];
int ref_diff = tracknumber - ref_truth_track;
//std::cout << "ref diff is " << ref_diff << std::endl;
int hit_diff=0;
int numhit=0;
for( std::tuple<int,double,double> thit : atrack ) {
hit_diff = hit_diff + tracknumber - tracktruth[std::get<0>(thit)]-ref_diff;
numhit++;
//std::cout << "track number " << tracknumber << " hit number " << std::get<0>(thit) << std::endl;
//std::cout << "track truth " << tracktruth[std::get<0>(thit)];
//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) );
singletrack->Fill( hits[std::get<0>(thit)].X(), hits[std::get<0>(thit)].Y() );
}
double grouped_test = hit_diff/numhit;
std::cout << "grouped test " << grouped_test << std::endl;
FoundTrackHists.push_back(singletrack);
}
// now, check for missed hits, hits that belonged to a track but were left out instead.
//
// doubled_hits = 0; // more complicated to check this, don't right now.
int missed_hits = 0;
int total_hits =0;
for( int hit_index=0; hit_index < tracktruth.size(); hit_index++ ) {
total_hits++;
//std::cout << " hit index " << hit_index << std::endl;
bool hit_listed = false;
for( std::vector<std::tuple<int,double,double>> atrack : FoundTracks ) {
for ( std::tuple<int,double,double> thit : atrack ) {
if (hit_index == std::get<0>(thit)) {
hit_listed = true;
}
}
}
if (!hit_listed) {
missed_hits++;
}
}
std::cout << "missed hits " << missed_hits << " total hits " << total_hits << std::endl;
TCanvas * c = new TCanvas();
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->cd(3);
huv->Draw();
c = new TCanvas();
c->Divide(2,2);
c->cd(1);
hxy2->Draw("box");
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");}
}
......@@ -27,11 +27,11 @@ void hough_transform2(
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", 100, -1010 ,1010, 100, -1010, 1010);
TH2F * hxy2 = new TH2F("hxy2", "hxy2", 1010, -1010 ,1010, 1010, -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);
TH2F * hpeaky2 = new TH2F("hpeakxy2", "hxy2", 50, -100 , 100, 50, -100, 100);
TH2F * hpeaky2 = new TH2F("hpeakxy2", "hxy2", 1000, -100 , 100, 1000, -100, 100);
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 );
......@@ -97,9 +97,9 @@ void hough_transform2(
//make track in a circle
int numhits = 7;
double radius = 3;
double centerx = 3;
int numhits = 200;
double radius = 100;
double centerx = 0;
double centery = 0;
for(int i_hit=0; i_hit < numhits; i_hit++){
......@@ -121,16 +121,18 @@ void hough_transform2(
for(const auto& fhit : hits){
//auto fhit = hits.at(0);
auto chits = ht->GetConformalCoordinates(fhit, hits);
//for(const auto& fhit : hits){
auto fhit = hits.at(50);
//auto chits = ht->GetConformalCoordinates(fhit, hits);
auto chits = ht->UseConformalMap(fhit, hits);
all_chits.insert(all_chits.end(), chits.begin(), chits.end());
//std::transform(hits.begin(), hits.end(), chits.begin(), master_hits.end(),
// [&](const auto& a, const auto& b){ return make_tuple(a,b,fhit);});
for(int ihit = 0; ihit < hits.size() ; ihit++){
master_hits.push_back(make_tuple(hits[ihit], chits[ihit], fhit));
}
}
//i}
/* std::cout << " Master " << master_hits.size() << std::endl;
if(master_hits.size()>2500){
......@@ -245,12 +247,12 @@ void hough_transform2(
//ht->fIntReduced->Draw();
c = new TCanvas();
huv->Draw("box");
huv->Draw("");
c = new TCanvas();
c->Divide(2,2);
c->cd(1);
hxy2->Draw("box");
hxy2->Draw("");
//for(auto ahist : FoundTrackHists){
// ahist->Draw("box,same");
//}
......
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