diff --git a/NPLib/TrackReconstruction/NPLinearRansac3D.cxx b/NPLib/TrackReconstruction/NPLinearRansac3D.cxx
index 8ef0839a43c82fd5b263a9b10b5351f813b596d1..6b2c17eaa6a38ceed4c138368a892f124bf0c1cb 100644
--- a/NPLib/TrackReconstruction/NPLinearRansac3D.cxx
+++ b/NPLib/TrackReconstruction/NPLinearRansac3D.cxx
@@ -17,8 +17,14 @@
*****************************************************************************/
#include "NPLinearRansac3D.h"
+#include<iostream>
using namespace std;
using namespace NPL;
+
+
+ClassImp(LinearCluster3D);
+
+ClassImp(LinearRansac3D);
////////////////////////////////////////////////////////////////////////////////
void LinearCluster3D::LinearFit() {
@@ -36,6 +42,7 @@ void LinearCluster3D::LinearFit() {
double meanZ = 0;
double w;
+
for (unsigned int i = 0; i < sizeI; i++) {
TVector3 Pi((*m_X)[m_index[i]], (*m_Y)[m_index[i]], (*m_Z)[m_index[i]]);
for (unsigned int j = i + 1; j < sizeI; j++) {
@@ -69,11 +76,200 @@ void LinearCluster3D::LinearFit() {
m_P0 = TVector3(meanX, meanY, meanZ);
m_Dir = (TVector3(meanDX, meanDY, meanDZ)).Unit();
+ m_isFitted = true;
// Cluster dir is arbitrary. we choose to always have positive Z dir
if (m_Dir.Z() < 0)
m_Dir = -m_Dir;
};
+void LinearCluster3D::FastLinearFit() {
+
+double xx=0, yy=0, zz=0, xx2=0, yy2=0, zz2=0, xy=0, xz=0, yz=0, qq=0;
+unsigned int sizeI = m_index.size();
+ for(unsigned int i = 0; i < sizeI; i++){
+ xx+=(*m_X)[m_index[i]];
+ yy+=(*m_Y)[m_index[i]];
+ zz+=(*m_Z)[m_index[i]];
+ xx2+=(*m_X)[m_index[i]]*(*m_X)[m_index[i]];
+ yy2+=(*m_Y)[m_index[i]]*(*m_Y)[m_index[i]];
+ zz2+=(*m_Z)[m_index[i]]*(*m_Z)[m_index[i]];
+ xy+=(*m_X)[m_index[i]]*(*m_Y)[m_index[i]];
+ xz+=(*m_X)[m_index[i]]*(*m_Z)[m_index[i]];
+ yz+=(*m_Y)[m_index[i]]*(*m_Z)[m_index[i]];
+ qq+=1;
+ }
+
+ double Xm = xx/qq, Ym = yy/qq, Zm = zz/qq;
+ double Sxx = -Xm*Xm+xx2/qq, Syy = -Ym*Ym+yy2/qq, Szz = -Zm*Zm+zz2/qq;
+ double Sxy = -Xm*Ym+xy/qq, Sxz = -Xm*Zm+xz/qq, Syz = -Ym*Zm+yz/qq;
+ double theta = 0.5*atan(2*Sxy/(Sxx-Syy));
+
+ double ct = cos(theta), st = sin(theta), ct2 = ct*ct, st2 = st*st, cst = cos(theta)*sin(theta);
+ double K11 = (Syy+Szz)*ct2+(Sxx+Szz)*st2-2*Sxy*cst;
+ double K22 = (Syy+Szz)*st2+(Sxx+Szz)*ct2+2*Sxy*cst;
+ //double K12 = -Sxy*(ct2-st2)+(Sxx-Syy)*cst; //=0 by construction
+ double K10 = Sxz*ct+Syz*st;
+ double K01 = -Sxz*st+Syz*ct;
+ double K00 = Sxx+Syy;
+
+ double C2 = -K00-K11-K22;
+ double C1 = K00*K11+K00*K22+K11*K22-K01*K01-K10*K10;
+ double C0 = K01*K01*K11+K10*K10*K22-K00*K11*K22;
+
+ double p = C1-C2*C2/3, q = 2*C2*C2*C2/27-C1*C2/3+C0;
+ double R = q*q/4+p*p*p/27;
+ double dm2 = 0;
+ if(R>0)
+ {
+ dm2 = -C2/3+pow(-q/2+sqrt(R), 1./3.)+pow(-q/2-sqrt(R), 1./3.);
+ }
+ else
+ {
+ double rho = sqrt(-p*p*p/27), phi = acos(-q/(2*rho)); double R = q*q/4+p*p*p/27;
+
+ double sqrt3rho = pow(rho, 1./3.);
+ double d1 = -C2/3+2*sqrt3rho*cos(phi/3);
+ double d2 = -C2/3+2*sqrt3rho*cos((phi+2*M_PI)/3);
+ double d3 = -C2/3+2*sqrt3rho*cos((phi+4*M_PI)/3);
+ dm2 = std::min(std::min(d1, d2), d3);
+ }
+ double a = -K10*ct/(K11-dm2)+K01*st/(K22-dm2);
+ double b = -K10*st/(K11-dm2)-K01*ct/(K22-dm2);
+
+ double denom = 1./(1+a*a+b*b);
+ double u = ((1+b*b)*Xm-a*b*Ym+a*Zm)*denom;
+ double v = (-a*b*Xm+(1+a*a)*Ym+b*Zm)*denom;
+ double w = (a*Xm+b*Ym+(a*a+b*b)*Zm)*denom;
+
+ m_P0 = TVector3(Xm, Ym, Zm);
+ m_Dir = (TVector3(Xm-u, Ym-v, Zm-w)).Unit();
+}
+
+void LinearCluster3D::FastLinearFitShort() {
+
+double xx=0, yy=0, zz=0, xx2=0, yy2=0, zz2=0, xy=0, xz=0, yz=0, qq=0;
+unsigned int sizeI = m_index_short.size();
+ for(unsigned int i = 0; i < sizeI; i++){
+ xx+=(*m_X)[m_index_short[i]];
+ yy+=(*m_Y)[m_index_short[i]];
+ zz+=(*m_Z)[m_index_short[i]];
+ xx2+=(*m_X)[m_index_short[i]]*(*m_X)[m_index_short[i]];
+ yy2+=(*m_Y)[m_index_short[i]]*(*m_Y)[m_index_short[i]];
+ zz2+=(*m_Z)[m_index_short[i]]*(*m_Z)[m_index_short[i]];
+ xy+=(*m_X)[m_index_short[i]]*(*m_Y)[m_index_short[i]];
+ xz+=(*m_X)[m_index_short[i]]*(*m_Z)[m_index_short[i]];
+ yz+=(*m_Y)[m_index_short[i]]*(*m_Z)[m_index_short[i]];
+ qq+=1;
+ }
+
+ double Xm = xx/qq, Ym = yy/qq, Zm = zz/qq;
+ double Sxx = -Xm*Xm+xx2/qq, Syy = -Ym*Ym+yy2/qq, Szz = -Zm*Zm+zz2/qq;
+ double Sxy = -Xm*Ym+xy/qq, Sxz = -Xm*Zm+xz/qq, Syz = -Ym*Zm+yz/qq;
+ double theta = 0.5*atan(2*Sxy/(Sxx-Syy));
+
+ double ct = cos(theta), st = sin(theta), ct2 = ct*ct, st2 = st*st, cst = cos(theta)*sin(theta);
+ double K11 = (Syy+Szz)*ct2+(Sxx+Szz)*st2-2*Sxy*cst;
+ double K22 = (Syy+Szz)*st2+(Sxx+Szz)*ct2+2*Sxy*cst;
+ //double K12 = -Sxy*(ct2-st2)+(Sxx-Syy)*cst; //=0 by construction
+ double K10 = Sxz*ct+Syz*st;
+ double K01 = -Sxz*st+Syz*ct;
+ double K00 = Sxx+Syy;
+
+ double C2 = -K00-K11-K22;
+ double C1 = K00*K11+K00*K22+K11*K22-K01*K01-K10*K10;
+ double C0 = K01*K01*K11+K10*K10*K22-K00*K11*K22;
+
+ double p = C1-C2*C2/3, q = 2*C2*C2*C2/27-C1*C2/3+C0;
+ double R = q*q/4+p*p*p/27;
+ double dm2 = 0;
+ if(R>0)
+ {
+ dm2 = -C2/3+pow(-q/2+sqrt(R), 1./3.)+pow(-q/2-sqrt(R), 1./3.);
+ }
+ else
+ {
+ double rho = sqrt(-p*p*p/27), phi = acos(-q/(2*rho)); double R = q*q/4+p*p*p/27;
+
+ double sqrt3rho = pow(rho, 1./3.);
+ double d1 = -C2/3+2*sqrt3rho*cos(phi/3);
+ double d2 = -C2/3+2*sqrt3rho*cos((phi+2*M_PI)/3);
+ double d3 = -C2/3+2*sqrt3rho*cos((phi+4*M_PI)/3);
+ dm2 = std::min(std::min(d1, d2), d3);
+ }
+ double a = -K10*ct/(K11-dm2)+K01*st/(K22-dm2);
+ double b = -K10*st/(K11-dm2)-K01*ct/(K22-dm2);
+
+ double denom = 1./(1+a*a+b*b);
+ double u = ((1+b*b)*Xm-a*b*Ym+a*Zm)*denom;
+ double v = (-a*b*Xm+(1+a*a)*Ym+b*Zm)*denom;
+ double w = (a*Xm+b*Ym+(a*a+b*b)*Zm)*denom;
+
+ m_P0 = TVector3(Xm, Ym, Zm);
+ m_Dir = (TVector3(Xm-u, Ym-v, Zm-w)).Unit();
+}
+
+void LinearCluster3D::CheckDirection(const TVector3 & vertex){
+ if(m_FurthestPoint == TVector3(0, 0, 0)){
+ double max_distance = 0;
+ for(int i = 0; i < m_index.size(); i++)
+ {
+ TVector3 aPoint((*m_X)[m_index[i]], (*m_Y)[m_index[i]], (*m_Z)[m_index[i]]);
+ double distance = (aPoint-vertex).Mag();
+ if(distance > max_distance){
+ max_distance = distance;
+ m_FurthestPoint = aPoint;
+ }
+ }
+ }
+ if((m_FurthestPoint-vertex).X()*m_Dir.X() < 0) this->Inverse();
+}
+
+void LinearCluster3D::CalculateTrackLength(const TVector3 & vertex, const double & percentRange){
+ TH1D* Hrange = new TH1D("Range", "range", 256, 0, 512);
+ TH1D* Hrangecumule = new TH1D("RangeCumulated", "rangecumule", 256, 0, 512);
+
+ double total_charge = 0;
+ double cumulated_charge = 0;
+ double max_distance = 0;
+
+
+ for(int i = 0; i < m_index.size(); i++)
+ {
+ TVector3 aPoint((*m_X)[m_index[i]], (*m_Y)[m_index[i]], (*m_Z)[m_index[i]]);
+ double distance = (aPoint-vertex).Mag();
+ total_charge+=(*m_Q)[m_index[i]];
+ Hrange->Fill(distance, (*m_Q)[m_index[i]]);
+ if(distance > max_distance){
+ max_distance = distance;
+ m_FurthestPoint = aPoint;
+ }
+ if(distance < m_ShortTrackMaxLength){
+ m_index_short.push_back(m_index[i]);
+ }
+ }
+ m_Charge = total_charge;
+ //Inverse the track if it is not going in the good direction
+ if((m_FurthestPoint-vertex).X()*m_Dir.X() < 0) this->Inverse();
+
+ for(int bin = 1; bin < Hrange->GetNbinsX(); bin++)
+ {
+ cumulated_charge+=Hrange->GetBinContent(bin);
+ Hrangecumule->SetBinContent(bin, cumulated_charge/total_charge);
+ }
+ Int_t nbin = 1;
+ while(Hrangecumule->GetBinContent(nbin) < percentRange && nbin < Hrangecumule->GetNbinsX()-1)
+ {
+ nbin++;
+ }
+ Double_t delta, deltabinap, deltabinav;
+ delta = Hrangecumule->GetBinContent(nbin)-Hrangecumule->GetBinContent(nbin-1);
+ deltabinap = Hrangecumule->GetBinContent(nbin)-percentRange;
+ deltabinav = percentRange-Hrangecumule->GetBinContent(nbin-1);
+ m_TrackLength = deltabinap*Hrangecumule->GetBinCenter(nbin-1)/delta+deltabinav*Hrangecumule->GetBinCenter(nbin)/delta;
+ delete Hrange;
+ delete Hrangecumule;
+}
+
////////////////////////////////////////////////////////////////////////////////
vector<LinearCluster3D> LinearRansac3D::TreatEvent(vector<double>& X, vector<double>& Y, vector<double>& Z) {
cluster_id.clear();
@@ -92,7 +288,7 @@ vector<LinearCluster3D> LinearRansac3D::TreatEvent(vector<double>& X, vector<dou
m_iteration_d = 0;
m_iteration = 0;
- while (m_iteration++ < m_max_iteration && m_assigned.size() < sizeX) {
+ while (m_iteration++ < m_max_iteration/* && m_assigned.size() < sizeX*/) {
LinearCluster3D cluster(&X, &Y, &Z);
// take 2 distant point randomly that has not been match before
d = 0;
@@ -123,7 +319,7 @@ vector<LinearCluster3D> LinearRansac3D::TreatEvent(vector<double>& X, vector<dou
}
} // while
-
+
// loop over the cluster starting with the biggest
unsigned int current_cluster = 0;
unsigned int index;
diff --git a/NPLib/TrackReconstruction/NPLinearRansac3D.h b/NPLib/TrackReconstruction/NPLinearRansac3D.h
index b824e6e21a93cb6d3f27a874165dc8b06217a402..cdbf0b08b9e1952256069dfc4ca0283ddcfd0cab 100644
--- a/NPLib/TrackReconstruction/NPLinearRansac3D.h
+++ b/NPLib/TrackReconstruction/NPLinearRansac3D.h
@@ -24,15 +24,25 @@
//root
#include "TVector3.h"
#include "TRandom2.h"
+#include "TH1.h"
// nptool
#include "NPTrackingUtility.h"
namespace NPL{
/////////////////////////////////
class LinearCluster3D{
public:
- LinearCluster3D(){};
+ LinearCluster3D(){
+ m_isFitted = false;
+ m_FurthestPoint = TVector3(0, 0, 0);
+ };
+
LinearCluster3D(std::vector<double>* X, std::vector<double>* Y, std::vector<double>* Z){
- m_X=X; m_Y=Y; m_Z=Z;
+ m_X=X; m_Y=Y; m_Z=Z; m_isFitted = false;
+ m_FurthestPoint = TVector3(0, 0, 0);
+ }
+ LinearCluster3D(std::vector<double>* X, std::vector<double>* Y, std::vector<double>* Z, std::vector<double>* Q){
+ m_X=X; m_Y=Y; m_Z=Z; m_Q=Q; m_isFitted = false;
+ m_FurthestPoint = TVector3(0, 0, 0);
}
~LinearCluster3D(){};
@@ -40,22 +50,48 @@ namespace NPL{
void AddIndex(const unsigned int& i) {m_index.push_back(i);};
unsigned int size() const {return m_index.size();};
unsigned int GetIndex (const unsigned int& i) const {return m_index[i];};
+ void SetShortTrackMaxLength(const double & l) {m_ShortTrackMaxLength = l;}
+ //Inverse the direction
+ void Inverse() {m_Dir = -m_Dir;}; //!
// Fit the cluster
- void LinearFit();
+ void LinearFit(); //!
+ void FastLinearFit(); //!
+ void FastLinearFitShort(); //!
+
+ //Calculate the track length from the interaction vertex, and for a given percent of charge deposited along the way
+ //It also inverses the track if needed, and regroups the points below a certain distance from the vertex to refit it without the deviation
+ void CalculateTrackLength(const TVector3 & vertex, const double & percentRange = 0.95); //!
+ void CheckDirection(const TVector3 & vertex);
private:
std::vector<double>* m_X;//!
std::vector<double>* m_Y;//!
std::vector<double>* m_Z;//!
+ std::vector<double>* m_Q;//!
+
std::vector<unsigned int> m_index;//!
- TVector3 m_P0;// a point on the line
- TVector3 m_Dir;// direction of the line
+ std::vector<unsigned int> m_index_short;//!
+ TVector3 m_P0; // a point on the line
+ TVector3 m_Dir; // direction of the line
+
+ double m_TrackLength; //!
+ double m_Charge; //!
+ // check if the track has been fitted with the "LinearFit" function
+ bool m_isFitted; //!
+ //Furthest point compared to vertex
+ TVector3 m_FurthestPoint; //!
+ double m_ShortTrackMaxLength; //!
+
public:
std::vector<unsigned int> GetIndex() const {return m_index;};
TVector3 GetP0()const {return m_P0;};
TVector3 GetDir()const {return m_Dir;};
+ double GetTrackLength()const {return m_TrackLength;}
+ bool IsFitted() const {return m_isFitted;}
+ TVector3 GetFurthestPoint()const {return m_FurthestPoint;}
+ double GetCharge()const {return m_Charge;}
public: //operator
// reversed logic so the cluster are ordered from largest to smallest
@@ -69,6 +105,8 @@ namespace NPL{
friend bool operator<(const LinearCluster3D& A, const LinearCluster3D& B){
return (A.size()>B.size());
}
+ ClassDef(LinearCluster3D, 0);
+
};
/////////////////////////////////
class LinearRansac3D{
@@ -109,6 +147,8 @@ namespace NPL{
public:
std::vector<LinearCluster3D> TreatEvent(std::vector<double>& X, std::vector<double>&Y, std::vector<double>&Z);
+
+ ClassDef(LinearRansac3D, 0);
};
diff --git a/NPLib/TrackReconstruction/NPRansacACTAR.cxx b/NPLib/TrackReconstruction/NPRansacACTAR.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..d8a9d035f5be393b7d49c2a8413584c793c1b62d
--- /dev/null
+++ b/NPLib/TrackReconstruction/NPRansacACTAR.cxx
@@ -0,0 +1,184 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2016 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * *
+ * Original Author : Damien THISSE contact address: damien.thisse@cea.fr *
+ * *
+ * Creation Date : October 2024 *
+ * Last update : October 2024 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class deal with finding all the track event by event *
+ *****************************************************************************/
+
+#include "NPRansacACTAR.h"
+#include "NPLinearRansac3D.h"
+
+using namespace NPL;
+using namespace std;
+
+/////////////////////////////////////////////////
+RansacACTAR::RansacACTAR()
+{
+ fRANSACMaxIteration = 1000;
+ fRANSACThreshold = 16;//100;
+ fRANSACPointThreshold = 0.05;//0.07;
+ fRANSACChargeThreshold = 0.05;//0.07;
+ fRANSACDistance = 7;//7;
+ fMAXBarycenterDistance = 10;
+ fAngleMax = 2;//degree
+
+ fNumberOfPadsX = 128;
+ fNumberOfPadsY = 128;
+
+ Rand = new TRandom3();
+}
+
+vector<LinearCluster3D> RansacACTAR::SimpleRansac(vector<double> & vX, vector<double> & vY, vector<double> & vZ, vector<double> & vQ){
+ fTotalCharge = 0;
+ fOriginalCloudSize = vX.size();
+ vector<int> index;
+ vector<LinearCluster3D> all_clusters;
+
+ for(int QQ = 0; QQ < vQ.size(); QQ++){
+ fTotalCharge+= vQ.at(QQ);
+ index.push_back(QQ);
+ }
+ double RemainingCharge = fTotalCharge;
+ int loop_counter = 0;
+ // cout << "Total charge is " << fTotalCharge << " ---> Threshold to end RANSAC is thus " << fTotalCharge*fRANSACChargeThreshold << endl;
+ // cout << "Total voxels is " << index.size() << " ---> Threshold to end RANSAC is thus " << fOriginalCloudSize*fRANSACPointThreshold << endl;
+ do{
+ loop_counter++;
+ //cout << "Loop " << loop_counter << ": remaining charge is " << RemainingCharge << "; number of points remaining : " << index.size() << endl;
+ LinearCluster3D foundCluster(&vX, &vY, &vZ, &vQ);
+ int highest_inliers = 0;
+ vector<int> index_to_remove;
+
+ for(int i = 0; i < fRANSACMaxIteration; i++){
+ LinearCluster3D aCluster(&vX, &vY, &vZ, &vQ);
+ vector<int> pos_index;
+ UInt_t rand1 = Rand->Integer(index.size());
+ UInt_t rand2;
+ do{
+ rand2 = Rand->Integer(index.size());
+ } while(rand2 == rand1);
+
+ TVector3 pt1(vX[index[rand1]], vY[index[rand1]], vZ[index[rand1]]);
+ TVector3 pt2(vX[index[rand2]], vY[index[rand2]], vZ[index[rand2]]);
+
+ for(int j = 0; j < index.size(); j++){
+ TVector3 pt(vX[index[j]], vY[index[j]], vZ[index[j]]);
+ if(MinimumDistancePointLine(pt1, pt2, pt) < fRANSACDistance){
+ aCluster.AddIndex(index[j]);
+ pos_index.push_back(j);
+
+ }
+ }
+ if(aCluster.size() > highest_inliers){
+ foundCluster = aCluster;
+ highest_inliers = aCluster.size();
+ index_to_remove = pos_index;
+ }
+ }
+
+ if(highest_inliers > fRANSACThreshold){
+ all_clusters.push_back(foundCluster);
+ for(int i = index_to_remove.size()-1; i >-1; i--){
+ RemainingCharge-= vQ.at(index[index_to_remove[i]]);
+ index.erase(index.begin()+index_to_remove.at(i));
+ }
+ // RemainingCharge = 0;
+ // for(int QQ = 0; QQ < index.size(); QQ++){
+ // RemainingCharge+= vQ.at(index.at(QQ));
+ // }
+ }
+ else break;
+
+ } while(RemainingCharge > fTotalCharge*fRANSACChargeThreshold && index.size() > fOriginalCloudSize*fRANSACPointThreshold);
+
+ return all_clusters;
+}
+
+void RansacACTAR::ReadParameterValue(string filename)
+{
+ bool ReadingStatus = false;
+
+ ifstream RansacParamFile;
+ RansacParamFile.open(filename.c_str());
+
+ if(!RansacParamFile.is_open()){
+ cout << "No Paramter File found for Ransac parameters" << endl;
+ cout << "Using the default parameter:" << endl;
+ cout << "RANSAC MaxIteration= " << fRANSACMaxIteration << endl;
+ cout << "RANSAC Threshold=" << fRANSACThreshold << endl;
+ cout << "RANSAC ChargeThreshold= " << fRANSACChargeThreshold << endl;
+ cout << "RANSAC Distance= " << fRANSACDistance << endl;
+ cout << "RANSAC PointThreshold= " << fRANSACPointThreshold << endl;
+ cout << "MAX Barycenter Distance= " << fMAXBarycenterDistance << endl;
+ cout << "Angle Max to merge tracks= " << fAngleMax << endl;
+ }
+ else{
+ string LineBuffer, whatToDo, DataBuffer;
+ while(!RansacParamFile.eof()){
+ getline(RansacParamFile,LineBuffer);
+ string name = "ConfigRansac";
+ if(LineBuffer.compare(0, name.length(), name) == 0){
+ cout << endl;
+ cout << "**** ConfigRansac found" << endl;
+ ReadingStatus = true;
+ }
+ while(ReadingStatus){
+ whatToDo="";
+ RansacParamFile >> whatToDo;
+ // Search for comment symbol (%)
+ if (whatToDo.compare(0, 1, "%") == 0) {
+ RansacParamFile.ignore(numeric_limits<streamsize>::max(), '\n' );
+ }
+ else if (whatToDo.compare(0,19,"RANSACMaxIteration=") == 0) {
+ RansacParamFile >> DataBuffer;
+ fRANSACMaxIteration = atoi(DataBuffer.c_str());
+ cout << "/// RANSAC MaxIteration= " << " " << fRANSACMaxIteration << " ///" << endl;
+ }
+ else if (whatToDo.compare(0,15,"RANSACDistance=") == 0) {
+ RansacParamFile >> DataBuffer;
+ fRANSACDistance = atoi(DataBuffer.c_str());
+ cout << "/// RANSAC Distance= " << " " << fRANSACDistance << " ///" << endl;
+ }
+ else if (whatToDo.compare(0,22,"RANSACChargeThreshold=") == 0) {
+ RansacParamFile >> DataBuffer;
+ fRANSACChargeThreshold = atof(DataBuffer.c_str());
+ cout << "/// RANSAC ChargeThreshold= " << " " << fRANSACChargeThreshold << " ///" << endl;
+ }
+ else if (whatToDo.compare(0,21,"RANSACPointThreshold=") == 0) {
+ RansacParamFile >> DataBuffer;
+ fRANSACPointThreshold = atof(DataBuffer.c_str());
+ cout << "/// RANSAC PointThreshold= " << " " << fRANSACPointThreshold << " ///" << endl;
+ }
+ else if (whatToDo.compare(0,16,"RANSACThreshold=") == 0) {
+ RansacParamFile >> DataBuffer;
+ fRANSACThreshold = atoi(DataBuffer.c_str());
+ cout << "/// RANSAC Threshold= " << " " << fRANSACThreshold << " ///" << endl;
+ }
+ else if (whatToDo.compare(0,22,"MAXBarycenterDistance=") == 0) {
+ RansacParamFile >> DataBuffer;
+ fMAXBarycenterDistance = atoi(DataBuffer.c_str());
+ cout << "/// Max Barycenter Distance= " << " " << fMAXBarycenterDistance << " ///" << endl;
+ }
+ else if (whatToDo.compare(0,16,"AngleMaxToMerge=") == 0) {
+ RansacParamFile >> DataBuffer;
+ fAngleMax = atoi(DataBuffer.c_str());
+ cout << "/// Angle Max to merge tracks= " << " " << fAngleMax << " ///" << endl;
+ }
+ else{
+ ReadingStatus = false;
+ }
+ }
+ }
+ }
+}
\ No newline at end of file
diff --git a/NPLib/TrackReconstruction/NPRansacACTAR.h b/NPLib/TrackReconstruction/NPRansacACTAR.h
new file mode 100644
index 0000000000000000000000000000000000000000..e7def0bccadce2b27884b245d7cac10d10c3b639
--- /dev/null
+++ b/NPLib/TrackReconstruction/NPRansacACTAR.h
@@ -0,0 +1,79 @@
+#ifndef __RANSACACTAR__
+#define __RANSACACTAR__
+/*****************************************************************************
+ * Copyright (C) 2009-2016 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * *
+ * Original Author : Damien THISSE contact address: damien.thisse@cea.fr *
+ * *
+ * Creation Date : October 2024 *
+ * Last update : October 2024 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class deal with finding all the track event by event *
+ *****************************************************************************/
+
+//C++ Header
+#include <stdio.h>
+#include <vector>
+#include <sstream>
+#include <iostream>
+#include <fstream>
+#include <cmath>
+#include <stdlib.h>
+#include <limits>
+#include <string>
+#include <algorithm>
+
+//NPL
+#include "NPTrack.h"
+#include "NPLinearRansac3D.h"
+
+using namespace NPL;
+
+// ROOT Headers
+#include <TH2F.h>
+#include <TMath.h>
+#include <TVector3.h>
+#include <TRandom3.h>
+#include <TVector.h>
+
+using namespace std;
+
+namespace NPL{
+
+ class RansacACTAR
+ {
+ public:
+ RansacACTAR();
+ ~RansacACTAR() {};
+
+ public:
+ void ReadParameterValue(string filename);
+ vector<NPL::LinearCluster3D> SimpleRansac(vector<double> & vX, vector<double> & vY, vector<double> & vZ, vector<double> & vQ);
+
+ private:
+ TRandom3* Rand;
+
+ private:
+ float fRANSACThreshold;
+ float fRANSACPointThreshold;
+ float fRANSACChargeThreshold;
+ float fRANSACDistance;
+ float fRANSACMaxIteration;
+ float fMAXBarycenterDistance;
+ float fAngleMax;
+ int fNumberOfTracksMax;
+ int fOriginalCloudSize;
+ double fTotalCharge;
+ int fNumberOfPadsX;
+ int fNumberOfPadsY;
+
+ };
+}
+#endif