From b36d9d27add80fc269d7a25bdf1b9f9265114fb3 Mon Sep 17 00:00:00 2001
From: adrien-matta <a.matta@surrey.ac.uk>
Date: Wed, 4 Nov 2020 15:39:24 +0100
Subject: [PATCH] * Working on FDC2
---
NPLib/TrackReconstruction/DCTest.cxx | 63 +++++++++++++++++++
.../NPDCReconstruction.cxx | 24 ++++---
2 files changed, 79 insertions(+), 8 deletions(-)
create mode 100644 NPLib/TrackReconstruction/DCTest.cxx
diff --git a/NPLib/TrackReconstruction/DCTest.cxx b/NPLib/TrackReconstruction/DCTest.cxx
new file mode 100644
index 000000000..89149ee6c
--- /dev/null
+++ b/NPLib/TrackReconstruction/DCTest.cxx
@@ -0,0 +1,63 @@
+#include"NPDCReconstruction.h"
+R__LOAD_LIBRARY(libNPTrackReconstruction.dylib)
+//
+void DrawTrack(double a, double b, int style);
+//
+void DrawCircle(double x, double z, double r);
+//
+void DCTest(){
+ // load the lib
+ auto c = new TCanvas();
+ TRandom3 rand;
+ // The ronconstruction utility object
+ NPL::DCReconstruction dcr;
+ vector<double> x,z,r;
+ // maximum drift length
+ double maxR = 10;
+ double stepZ = 50;
+ // infinite loop
+ while(true){
+ x.clear();z.clear();r.clear();
+ // create a random trajectory in 2D //
+ // a and b are the parameter of the 2D track
+ double a = rand.Uniform(-10,10);
+ double b = rand.Uniform(-10,10);
+ DrawTrack(a,b,2);
+ // Generate n points from the line
+ int n = (int) rand.Uniform(3,20);
+ for(unsigned int i = 0 ; i < n ; i++){
+ double Z = i*stepZ-100;
+ z.push_back(Z);
+ double R = rand.Uniform(-maxR,maxR);
+ r.push_back(R);
+ double X = (Z-b)/a+R;
+ x.push_back(X);
+ DrawCircle(X,Z,R);
+ }
+ double x0,x100,af,bf;
+ dcr.BuildTrack2D(x,z,r,x0,x100,af,bf);
+ DrawTrack(af,bf,1);
+ c->Update();
+ gPad->WaitPrimitive();
+ }
+}
+// draw a track based on z=ax+b
+void DrawTrack(double a, double b, int style){
+ double x[2]={-100,100};
+ double z[2]={a*-100+b,a*100+b};
+ auto l = new TGraph(2,x,z);
+ l->SetLineColor(kBlack);
+ if(style>1)
+ l->SetLineWidth(2);
+ l->SetLineStyle(style);
+ if(style==2)
+ l->Draw("al");
+ else
+ l->Draw("l");
+}
+// draw a circle
+void DrawCircle(double x, double z, double r){
+ auto e = new TEllipse(x,z,r,r);
+ e->SetLineColor(kBlack);
+ e->Draw("same");
+}
diff --git a/NPLib/TrackReconstruction/NPDCReconstruction.cxx b/NPLib/TrackReconstruction/NPDCReconstruction.cxx
index 7f8337a24..8221628fe 100644
--- a/NPLib/TrackReconstruction/NPDCReconstruction.cxx
+++ b/NPLib/TrackReconstruction/NPDCReconstruction.cxx
@@ -26,6 +26,11 @@ DCReconstruction::DCReconstruction(){
m_min=ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
m_func=ROOT::Math::Functor(this,&NPL::DCReconstruction::SumD,2);
m_min->SetFunction(m_func);
+ m_min->SetPrintLevel(0);
+ //m_min->SetMaxFunctionCalls(1000);
+ //m_min->SetMaxIterations(1000);
+ m_min->SetTolerance(1e64);
+ m_min->SetPrecision(1e-10);
}
////////////////////////////////////////////////////////////////////////////////
DCReconstruction::~DCReconstruction(){
@@ -42,17 +47,19 @@ double DCReconstruction::BuildTrack2D(const vector<double>& X,const vector<doubl
// Define the starting point of the fit: a straight line passing through the
// the first and last wire
// z = ax+b -> x=(z-b)/a
- ai = (Z[size-1]-Z[0])/(X[size-1]-R[size-1]-X[0]-R[0]);
- bi = Z[0]-ai*(X[0]+R[0]);
+// ai = (Z[size-1]-Z[0])/(X[size-1]-R[size-1]-X[0]-R[0]);
+// bi = Z[0]-ai*(X[0]+R[0]);
+ ai = (Z[size-1]-Z[0])/(X[size-1]-X[0]);
+ bi = Z[0]-ai*(X[0]);
+
parameter[0]=ai;
parameter[1]=bi;
+
+ m_min->Clear();
m_min->SetVariable(0,"a",parameter[0],1000);
m_min->SetVariable(1,"b",parameter[1],1000);
- m_min->SetTolerance(0.01);
-
// Perform minimisation
m_min->Minimize();
-
// access set of parameter that produce the minimum
const double *xs = m_min->X();
a=xs[0];
@@ -106,7 +113,7 @@ void DCReconstruction::ResolvePlane(const TVector3& L,const double& ThetaU ,cons
////////////////////////////////////////////////////////////////////////////////
double DCReconstruction::SumD(const double* parameter ){
- unsigned int size = fitX->size();
+ unsigned int sizeX = fitX->size();
// Compute the sum P of the distance between the circle and the track
P = 0;
a = parameter[0];
@@ -114,13 +121,14 @@ double DCReconstruction::SumD(const double* parameter ){
ab= a*b;
a2=a*a;
- for(unsigned int i = 0 ; i < size ; i++){
+ for(unsigned int i = 0 ; i < sizeX ; i++){
c = (*fitX)[i];
d = (*fitZ)[i];
r = (*fitR)[i];
x = (a*d-ab+c)/(1+a2);
z = a*x+b;
- P+= sqrt(abs( (x-c)*(x-c)+(z-d)*(z-d)-r*r));
+ //P+= sqrt(abs( (x-c)*(x-c)+(z-d)*(z-d)-r*r));
+ P+= abs( (x-c)*(x-c)+(z-d)*(z-d)-r*r);
}
// return normalized power
return P/size;
--
GitLab