Commit 67c59b15 authored by Johnston's avatar Johnston
Browse files

fix the merge conflicts

parents 4928755e d267ad98
......@@ -17,7 +17,6 @@ Proposal:
Simple transform (Hough or otherwise) to feature space given event hits,
binning, and maximum bin selection process to identify candidate hits.
Need:
Transforms suitable to geometry, B field. (Hough transform first)
......
......@@ -2,6 +2,8 @@
#define genfind_GenFindHits_HH
#include <memory>
#include <functional>
#include <algorithm>
#include "TObject.h"
#include "Math/Vector3D.h"
#include "Math/Vector4D.h"
......@@ -22,6 +24,8 @@ namespace genfind {
public:
Hit();
Hit(const ROOT::Math::XYZTVector& );
Hit(double x, double y, double z, double t=0);
Hit(const Hit&) = default;
Hit(Hit&&) = default;
Hit& operator=(const Hit&) = default;
......@@ -41,7 +45,6 @@ namespace genfind {
/** Conform space hits.
*/
class ConformalHit {
public:
ROOT::Math::XYZTVector fPosition = {0,0,0,0};
std::shared_ptr<genfind::Hit> fImageHit;
......@@ -63,7 +66,53 @@ namespace genfind {
ClassDef(ConformalHit,1)
};
/** @brief Compute the conformal transform of a hit.
*
* @param hit The hit that provides the position space coordinates.
* @param ref Hit providing the reference coordinates (x0,y0).
*
* \f$ u = \frac{x}{x^2+y^2} \f$
* \f$ v = \frac{y}{x^2+y^2} \f$
*
*/
ConformalHit compute_conformal_hit(std::shared_ptr<Hit> hit, std::shared_ptr<Hit> ref);
ConformalHit compute_conformal_hit(std::shared_ptr<Hit> hit);
template<class T>
ConformalHit compute_conformal_hit(const T& hit, const T& ref) {
ConformalHit chit;
chit.fImageHit = std::make_shared<genfind::Hit>(hit.X(), hit.Y(), 0.0) ;
chit.fImageRef = std::make_shared<genfind::Hit>(ref.X(), ref.Y(), 0.0) ;
double x = hit.X() - ref.X();
double y = hit.Y() - ref.Y();
double R = x*x + y*y;
if( !(hit == ref) ) {
chit.fPosition = { 2.0*x/R, 2.0*y/R, hit.Z(), hit.T() } ;
}
return chit;
}
template<class T>
std::vector<ConformalHit> compute_conformal_hits(const std::vector<T>& hits, const T& ref)
{
// Local class
struct Sum {
std::vector<ConformalHit> fChits;
T fRef;
Sum(const T& R) : fChits{}, fRef(R) { }
void operator()(const T& h){
fChits.push_back(compute_conformal_hit(h, fRef));
}
std::vector<ConformalHit> Get(){ return fChits; }
};
return std::for_each(hits.begin(),
hits.end(),
Sum(ref)
).Get();
}
}
......
......@@ -2,7 +2,9 @@
#define genfind_HoughTransform_HH
#include "TObject.h"
#include "TMath.h"
#include <array>
#include <functional>
#include <vector>
#include <utility>
#include "TH2F.h"
......@@ -11,9 +13,11 @@
#include "TSpectrum.h"
#include "TSpectrum2.h"
#include "TMultiGraph.h"
#include "GenFindHits.h"
namespace genfind {
class HoughTransform : public TObject {
public:
......@@ -39,7 +43,6 @@ namespace genfind {
public:
HoughTransform();
virtual ~HoughTransform();
std::vector<ROOT::Math::XYZTVector> GetConformalCoordinates(
......@@ -49,6 +52,61 @@ namespace genfind {
std::vector<ROOT::Math::XYZTVector> UseConformalMap(
const::std::vector<ROOT::Math::XYZTVector>& hits) const;
//====== delete the following from a merge if it conflicts with some of the other stuff.
//
/** @brief Method to get the result.
*
* @param hits Vector of hit type T which.
*
* Returns vector of vector of hits group into possilbe track candidates.
* Note: the reference hit for the conformal transform is at the origin.
*
*/
template<class T>
std::vector<std::vector<T>> operator()(const std::vector<T>& hits){
std::vector<std::vector<T>> res;
T ref = {0.0,0.0,0.0,0.0};
auto chits = genfind::compute_conformal_hits(hits, ref);
auto peaks = FindPhiPeaks(chits);
res = SelectPhiPeaks<T>(peaks, chits);
return res;
}
template<class T>
std::vector<std::vector<T>> SelectPhiPeaks(const std::vector<double>& peaks,
const std::vector<genfind::ConformalHit>& chits)
{
double degree = TMath::Pi()/180.0;
std::vector<std::vector<T>> res;
for(auto& p : peaks){
std::vector<T> p_hits;
for(const auto& ahit : chits){
double phi = TMath::ATan(ahit.Y()/ahit.X())/degree;
if( TMath::Abs(phi-p) < fPhiThresh ) {
p_hits.push_back(
{ahit.fImageHit->X(),
ahit.fImageHit->Y(),
ahit.fImageHit->Z(),
ahit.fImageHit->T()});
}
}
res.push_back(p_hits);
}
return res;
}
/** @brief Fill the Hough transform histograms for searching.
*
* @param chits Vector of ConformalHits used to find lines.
*
*/
void FillHistograms(const std::vector<genfind::ConformalHit>& chits);
std::vector<double> FindPhiPeaks(const std::vector<genfind::ConformalHit>& chits);
// next line, end deletes if the merge caused a conflict.
//>>>>>>> d267ad98e7036905dfceb221e8b3f316819b0990
TMultiGraph* FillHoughTransforms(const std::vector<ROOT::Math::XYZTVector>& hits);
std::vector<TF1*> FillFunctions(const std::vector<ROOT::Math::XYZTVector>& hits);
......
......@@ -6,6 +6,14 @@ namespace genfind {
{ }
//__________________________________________________________________________
Hit::Hit(const ROOT::Math::XYZTVector& p) : fPosition(p)
{ }
//__________________________________________________________________________
Hit::Hit(double x, double y, double z, double t) : fPosition({x,y,z,t})
{ }
//__________________________________________________________________________
Hit::~Hit()
{ }
//__________________________________________________________________________
......@@ -27,9 +35,18 @@ namespace genfind {
double y = hit->Y() - ref->Y();
double R = x*x + y*y;
if( !(hit == ref) ) {
chit.fPosition = { 2.0*x/R, -2.0*y/R, hit->Z(), hit->T() } ;
chit.fPosition = { 2.0*x/R, 2.0*y/R, hit->Z(), hit->T() } ;
}
return chit;
}
//__________________________________________________________________________
ConformalHit compute_conformal_hit(std::shared_ptr<Hit> hit)
{
// Use ref = {0,0};
return compute_conformal_hit(hit, std::make_shared<Hit>());
}
//__________________________________________________________________________
}
......@@ -71,93 +71,7 @@ namespace genfind {
double y3 = bhit.Y();
if( !(ahit == ref_point) ) {
//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;
}
res.push_back(ROOT::Math::XYZTVector{ 2.0*x/R, 2.0*y/R, ahit.Z(), ahit.T() });
}
}
};
......@@ -183,11 +97,73 @@ namespace genfind {
//_________________________________________________________________________
//
void HoughTransform::FillHistograms(const std::vector<genfind::ConformalHit>& chits)
{
double degree = TMath::Pi()/180.0;
double phi_min = 0.0;
double phi_max = 180.0;
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 : chits) {
temp0->Reset();
temp1->Reset();
temp2->Reset();
for(int i_theta = 0; i_theta<1000; i_theta++) {
double theta = 180.0*degree*double(i_theta)/double(500);//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);
}
delete temp0;
delete temp1;
delete temp2;
}
//__________________________________________________________________________
std::vector<double> HoughTransform::FindPhiPeaks(const std::vector<genfind::ConformalHit>& chits)
{
double degree = TMath::Pi()/180.0;
// First fill the phi histogram
fPhi->Reset();
for(const auto& ahit : chits){
double phi = TMath::ATan(ahit.Y()/ahit.X())/degree;
fPhi->Fill(phi);
}
double sigma = 2;
double thresh = 0.1;
int npeaks = fSpec1->Search(fPhi, sigma, "nodraw", thresh );
std::vector<double> peaks;
for(int i = 0; i<npeaks; i++) {
peaks.push_back( (fSpec1->GetPositionX())[i] );
}
return peaks;
}
//__________________________________________________________________________
TMultiGraph* HoughTransform::FillHoughTransforms(const std::vector<ROOT::Math::XYZTVector>& hits)
{
double degree = TMath::Pi()/180.0;
double phi_min = 0.0;
double phi_max = 180.0;
// This part is not really needed
TMultiGraph* mg = new TMultiGraph();
//std::vector<TF1*> funcs;
double avg = 0;
......@@ -316,7 +292,6 @@ namespace genfind {
}
int n_funcs = funcs.size();
//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 ){
......
#include "GenFind/GenFindHits.h"
void conformal_transform_test()
{
using namespace genfind;
std::vector<std::shared_ptr<Hit>> hits;
std::vector<ConformalHit> chits;
std::vector<ROOT::Math::XYZTVector> hits2;
std::vector<ConformalHit> chits2;
double degree = TMath::Pi()/180.0;
double r0 = 10.0;
double x0 = 0.0;
double y0 = r0;
ROOT::Math::XYZTVector ref_hit = {0.0,0.0,0.0,0.0};
for(int ith=0; ith < 10; ith++) {
double x = r0*TMath::Sin((5.0 + ith*3.0)*degree)- x0 ;
double y = TMath::Sqrt(r0*r0 - (x+x0)*(x+x0)) - y0 ;
auto ahit = std::make_shared<Hit>(x,y,0);
hits.push_back(ahit);
chits.push_back(compute_conformal_hit(ahit));
std::cout << " x " << x << ", y " << y << std::endl;
ROOT::Math::XYZTVector pos_hit = {x,y,0.0,0.0};
hits2.push_back(pos_hit);
//chits2.push_back(compute_conformal_hit(pos_hit, ref_hit));
}
chits2 = compute_conformal_hits(hits2,ref_hit);
TGraph* gr_hits = new TGraph();
int i_point = 0;
for(auto ahit : hits){
gr_hits->SetPoint(i_point, ahit->X(), ahit->Y());
i_point++;
}
TGraph* gr_chits = new TGraph();
i_point = 0;
for(auto ahit : chits){
gr_chits->SetPoint(i_point, ahit.X(), ahit.Y());
std::cout << " x " << ahit.X() << ", y " << ahit.Y() << std::endl;
i_point++;
}
TGraph* gr_hits2 = new TGraph();
i_point = 0;
for(auto ahit : hits2){
gr_hits2->SetPoint(i_point, ahit.X(), ahit.Y());
i_point++;
}
TGraph* gr_chits2 = new TGraph();
i_point = 0;
for(auto ahit : chits2){
gr_chits2->SetPoint(i_point, ahit.X(), ahit.Y());
std::cout << " x " << ahit.X() << ", y " << ahit.Y() << std::endl;
i_point++;
}
TCanvas* c = new TCanvas();
c->Divide(2,2);
c->cd(1);
gr_hits->SetMarkerStyle(20);
gr_hits->Draw("ap");
c->cd(2);
gr_chits->SetMarkerStyle(20);
gr_chits->Draw("ap");
c->cd(3);
gr_hits2->SetMarkerStyle(20);
gr_hits2->Draw("ap");
c->cd(4);
gr_chits2->SetMarkerStyle(20);
gr_chits2->Draw("ap");
}
......@@ -88,15 +88,6 @@ void fit_test1(
// smeared start state
genfit::MeasuredStateOnPlane stateSmeared(rep);
// approximate covariance
TMatrixDSym covM(6);
double res1 = 0.01;
for (int i = 0; i < 3; ++i)
covM(i,i) = res1*res1;
for (int i = 3; i < 6; ++i)
covM(i,i) = pow(res1/9.0/sqrt(3.0), 2);
const int detId(0); // detector ID
int planeId(0); // detector plane ID
int hitId(0); // hit ID
......@@ -117,6 +108,14 @@ void fit_test1(
}
auto first_hit_pos = (*((*g4hits)[0])).position;
// approximate covariance
TMatrixDSym covM(6);
double res1 = 0.01;
for (int i = 0; i < 3; ++i)
covM(i,i) = res1*res1;
for (int i = 3; i < 6; ++i)
covM(i,i) = pow(res1/9.0/sqrt(3.0), 2);
// start values for the fit, e.g. from pattern recognition
TVector3 pos(0.0,0.0,0.0);
......
......@@ -48,7 +48,7 @@ void hough_transform2(
t->GetEntry(i_event);
i_event++;
for( auto thit : (*g4hits) ) {
for( auto thit : (*g4hits2) ) {
hits.push_back( {thit->position.X(), thit->position.Y(), thit->position.Z(), 0.0} );
std::cout << thit->position.X() << " , " << thit->position.Y() << std::endl;
hxy->Fill(thit->position.X(), thit->position.Y());
......@@ -56,8 +56,7 @@ void hough_transform2(
hyz->Fill(thit->position.Y(), thit->position.Z());
hxy2->Fill(thit->position.X(), thit->position.Y());
}
for( auto thit : (*g4hits2) ) {
for( auto thit : (*g4hits) ) {
hits.push_back( {thit->position.X(), thit->position.Y(), thit->position.Z(), 0.0} );
std::cout << thit->position.X() << " , " << thit->position.Y() << std::endl;
hxy->Fill(thit->position.X(), thit->position.Y());
......@@ -65,6 +64,7 @@ void hough_transform2(
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";
*/
......@@ -119,7 +119,6 @@ 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(50);
......
void hough_transform3(
int i_event = 20,
int Ntracks = 1)
{
gSystem->Load("libGenFind.so");
using namespace genfind;
using namespace ROOT::Math;
double degree = TMath::Pi()/180.0;
TFile * f = new TFile("simple_example_out2.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;
}
// Setup branches
std::vector<DD4hep::Simulation::Geant4Tracker::Hit*> * g4hits = nullptr;
std::vector<DD4hep::Simulation::Geant4Tracker::Hit*> * g4hits2 = nullptr;
t->SetBranchAddress("SiVertexBarrelHits", &g4hits);
t->SetBranchAddress("SiTrackerBarrelHits", &g4hits2);
// Hough Transform class that has all methods. (Would be best to move some member functions to just functions)
HoughTransform* ht = new HoughTransform();
std::vector<XYZTVector> hits ;
// Get the track hits
for(int i_track = 0; i_track < Ntracks; i_track++){
t->GetEntry(i_event);
i_event++;
// SiTrackerBarrelHits
for( auto thit : (*g4hits2) ) {
hits.push_back( {thit->position.X(), thit->position.Y(), thit->position.Z(), 0.0} );
std::cout << thit->position.X() << " , " << thit->position.Y() << std::endl;
}
// SiVertexBarrelHits
for( auto thit : (*g4hits) ) {
hits.push_back( {thit->position.X(), thit->position.Y(), thit->position.Z(), 0.0} );
std::cout << thit->position.X() << " , " << thit->position.Y() << std::endl;
}
}
// ----------------------------------------------
// A point at the origin
XYZTVector zero_ref = {0,0,0,0};
// ----------------------------------------------
// This returns the hits in conformal space relative to zero_ref=(x0,y0)
// i.e.: u = (x-x0)/((x-x0)^2+(y-y0)^2), v = -(y-y0)/((x-x0)^2+(y-y0)^2)
auto chits = ht->GetConformalCoordinates(zero_ref, hits);