diff --git a/NPLib/TrackReconstruction/NPDCReconstruction.cxx b/NPLib/TrackReconstruction/NPDCReconstruction.cxx
index fd24aee035b727ed639ed0f31a0cc12eadef849b..0f3b4811f40f64a4708768504155aab17256ec48 100644
--- a/NPLib/TrackReconstruction/NPDCReconstruction.cxx
+++ b/NPLib/TrackReconstruction/NPDCReconstruction.cxx
@@ -16,7 +16,7 @@
* This class have all the method needed to analyse Drift Chambers *
*****************************************************************************/
-#include"NPDCReconstruction.h"
+#include "NPDCReconstruction.h"
#include "Math/Factory.h"
#include "TError.h"
#include "TGraph.h"
@@ -24,149 +24,146 @@ using namespace std;
using namespace NPL;
////////////////////////////////////////////////////////////////////////////////
-DCReconstruction::DCReconstruction(){
- m_min=ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
- m_func=ROOT::Math::Functor(this,&NPL::DCReconstruction::SumD,2);
+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);
- // this avoid error
+ // this avoid error
gErrorIgnoreLevel = kError;
- //m_min->SetMaxFunctionCalls(1000);
- //m_min->SetMaxIterations(1000);
- //m_min->SetTolerance(1);
- //m_min->SetPrecision(1e-10);
+ // m_min->SetMaxFunctionCalls(1000);
+ // m_min->SetMaxIterations(1000);
+ // m_min->SetTolerance(1);
+ // m_min->SetPrecision(1e-10);
}
////////////////////////////////////////////////////////////////////////////////
-DCReconstruction::~DCReconstruction(){
- delete m_min;
-}
+DCReconstruction::~DCReconstruction() { delete m_min; }
////////////////////////////////////////////////////////////////////////////////
-double DCReconstruction::BuildTrack2D(const vector<double>& X,const vector<double>& Z,const vector<double>& R,double& X0,double& X100,double& a, double& b ){
- fitX=&X;
- fitZ=&Z;
- fitR=&R;
+double DCReconstruction::BuildTrack2D(const vector<double>& X, const vector<double>& Z, const vector<double>& R,
+ double& X0, double& X100, double& a, double& b) {
+ fitX = &X;
+ fitZ = &Z;
+ fitR = &R;
// assume all X,Z,R of same size
sizeX = X.size();
- // Define the starting point of the fit: a straight line passing through the
+ // 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
unsigned int i = 1;
- ai=1/0.;
- while(isinf(ai)&&i!=sizeX){
- ai = (Z[sizeX-i]-Z[0])/(X[sizeX-i]-X[0]);
- bi = Z[0]-ai*(X[0]);
+ ai = 1 / 0.;
+ while (std::isinf(ai) && i != sizeX) {
+ ai = (Z[sizeX - i] - Z[0]) / (X[sizeX - i] - X[0]);
+ bi = Z[0] - ai * (X[0]);
i++;
}
- if(isinf(ai)){ // then there is no two point in different layer
- a=-10000;
- b=-10000;
- X0=-10000;
- X100=-10000;
+ if (std::isinf(ai)) { // then there is no two point in different layer
+ a = -10000;
+ b = -10000;
+ X0 = -10000;
+ X100 = -10000;
return 10000;
}
- m_min->Clear();
- m_min->SetVariable(0,"a",ai,1);
- m_min->SetVariable(1,"b",bi,1);
+ m_min->Clear();
+ m_min->SetVariable(0, "a", ai, 1);
+ m_min->SetVariable(1, "b", bi, 1);
// Perform minimisation
- m_min->Minimize();
- //std::cout << "EDM:" << m_min->Edm() << std::endl;
- // access set of parameter that produce the minimum
- const double *xs = m_min->X();
- a=xs[0];
- b=xs[1];
- X0=-b/a;
- X100=(100-b)/a;
- return m_min->MinValue() ;
+ m_min->Minimize();
+ // std::cout << "EDM:" << m_min->Edm() << std::endl;
+ // access set of parameter that produce the minimum
+ const double* xs = m_min->X();
+ a = xs[0];
+ b = xs[1];
+ X0 = -b / a;
+ X100 = (100 - b) / a;
+ return m_min->MinValue();
}
////////////////////////////////////////////////////////////////////////////////
-void DCReconstruction::ResolvePlane(const TVector3& L,const double& ThetaU ,const TVector3& H, const double& ThetaV, TVector3& PosXY){
+void DCReconstruction::ResolvePlane(const TVector3& L, const double& ThetaU, const TVector3& H, const double& ThetaV,
+ TVector3& PosXY) {
// direction of U and V wire
- TVector3 u(0,1,0);
+ TVector3 u(0, 1, 0);
u.RotateZ(ThetaU);
- TVector3 v(0,1,0);
+ TVector3 v(0, 1, 0);
v.RotateZ(ThetaV);
-
// Compute the coeff of the two line of vecotr u (v) going through H (L)
// dv : y = av*x+bv
- av = v.Y()/v.X();
- bv = H.Y() - av*H.X();
+ av = v.Y() / v.X();
+ bv = H.Y() - av * H.X();
// du : y = au*x+bu
- au = u.Y()/u.X();
- bu = L.Y() - au*L.X();
+ au = u.Y() / u.X();
+ bu = L.Y() - au * L.X();
// We look for M(xM, yM) that intersect du and dv:
- if(!isinf(au) && !isinf(av)){ // au and av are not inf, i.e. not vertical line
- xM = (bv-bu)/(au-av);
- yM = au*xM+bu;
+ if (!std::isinf(au) && !std::isinf(av)) { // au and av are not inf, i.e. not vertical line
+ xM = (bv - bu) / (au - av);
+ yM = au * xM + bu;
}
- else if(isinf(av)){// av is inf, so v is along Y axis, H is direct measure of X
+ else if (std::isinf(av)) { // av is inf, so v is along Y axis, H is direct measure of X
xM = H.X();
- yM = au*xM+bu;
+ yM = au * xM + bu;
}
- else if (isinf(au)){//au is inf, so u is along Y axis, L is direct measure of X
+ else if (std::isinf(au)) { // au is inf, so u is along Y axis, L is direct measure of X
xM = L.X();
- yM = av*xM+bv;
+ yM = av * xM + bv;
}
- else{ // all is lost
- xM=-10000;
- yM=-10000;
+ else { // all is lost
+ xM = -10000;
+ yM = -10000;
}
- PosXY.SetXYZ(xM,yM,0);
+ PosXY.SetXYZ(xM, yM, 0);
}
-
////////////////////////////////////////////////////////////////////////////////
-double DCReconstruction::SumD(const double* parameter ){
+double DCReconstruction::SumD(const double* parameter) {
// Compute the sum P of the distance between the circle and the track
P = 0;
a = parameter[0];
b = parameter[1];
- ab= a*b;
- a2=a*a;
- for(unsigned int i = 0 ; i < sizeX ; i++){
+ ab = a * b;
+ a2 = a * a;
+ 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= (x-c)*(x-c)+(z-d)*(z-d)-r*r;
+ x = (a * d - ab + c) / (1 + a2);
+ z = a * x + b;
+ p = (x - c) * (x - c) + (z - d) * (z - d) - r * r;
// numerical trick to have a smooth derivative instead of using abs
- P+= sqrt(p*p+0.1);
+ P += sqrt(p * p + 0.1);
}
// return normalized power
- return P/sizeX;
+ return P / sizeX;
}
-
////////////////////////////////////////////////////////////////////////////////
-TGraph* DCReconstruction::Scan(double a, double b, int tovary, double minV, double maxV){
- vector<double> x,y;
- unsigned int sizeT=1000;
- double step = (maxV-minV)/sizeT;
- double p[2]={a,b};
- for(unsigned int i = 0 ; i < sizeT ; i++){
- if(!tovary){
- p[0]=minV+step*i;
+TGraph* DCReconstruction::Scan(double a, double b, int tovary, double minV, double maxV) {
+ vector<double> x, y;
+ unsigned int sizeT = 1000;
+ double step = (maxV - minV) / sizeT;
+ double p[2] = {a, b};
+ for (unsigned int i = 0; i < sizeT; i++) {
+ if (!tovary) {
+ p[0] = minV + step * i;
x.push_back(p[0]);
y.push_back(SumD(p));
}
- else{
- p[1]=minV+step*i;
+ else {
+ p[1] = minV + step * i;
x.push_back(p[1]);
y.push_back(SumD(p));
}
}
- TGraph* g = new TGraph(x.size(),&x[0],&y[0]);
+ TGraph* g = new TGraph(x.size(), &x[0], &y[0]);
return g;
}
diff --git a/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx b/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx
index 390a6272281c105481431a3cdfa4a2620d90b9f0..cca5634728e8015af0a8f9b6ecc6d1cf98310e10 100644
--- a/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx
+++ b/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx
@@ -1,4 +1,4 @@
-#if __cplusplus > 199711L // require c++11
+#if __cplusplus > 199711L // require c++11
/*****************************************************************************
* Copyright (C) 2009-2020 this file is part of the NPTool Project *
* *
@@ -17,183 +17,180 @@
* This class have all the method needed to analyse Drift Chambers *
*****************************************************************************/
-#include"NPDCReconstructionMT.h"
+#include "NPDCReconstructionMT.h"
// ROOT
#include "Math/Factory.h"
-#include "Math/Minimizer.h"
#include "Math/Functor.h"
+#include "Math/Minimizer.h"
#include "TError.h"
#include "TGraph.h"
-#include "TVector3.h"
#include "TROOT.h"
+#include "TVector3.h"
using namespace std;
using namespace NPL;
////////////////////////////////////////////////////////////////////////////////
-DCReconstructionMT::DCReconstructionMT(unsigned int number_thread){
+DCReconstructionMT::DCReconstructionMT(unsigned int number_thread) {
ROOT::EnableThreadSafety();
- m_nbr_thread= number_thread;
+ m_nbr_thread = number_thread;
// force loading of the minimizer plugin ahead
- ROOT::Math::Minimizer* mini=ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
+ ROOT::Math::Minimizer* mini = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
delete mini;
}
////////////////////////////////////////////////////////////////////////////////
-DCReconstructionMT::~DCReconstructionMT(){
- StopThread();
-}
+DCReconstructionMT::~DCReconstructionMT() { StopThread(); }
////////////////////////////////////////////////////////////////////////////////
-double DCReconstructionMT::GetResults(unsigned int ID,double& X0,double& X100,double& a, double& b){
- if(m_X0.find(ID)!=m_X0.end()){
- X0=m_X0[ID];
- X100=m_X100[ID];
- a=m_a[ID];
- b=m_b[ID];
+double DCReconstructionMT::GetResults(unsigned int ID, double& X0, double& X100, double& a, double& b) {
+ if (m_X0.find(ID) != m_X0.end()) {
+ X0 = m_X0[ID];
+ X100 = m_X100[ID];
+ a = m_a[ID];
+ b = m_b[ID];
return m_minimum[ID];
}
return -1;
}
////////////////////////////////////////////////////////////////////////////////
-void DCReconstructionMT::AddPlan(unsigned int ID, const vector<double>& X, const vector<double>& Z, const vector<double>& R){
+void DCReconstructionMT::AddPlan(unsigned int ID, const vector<double>& X, const vector<double>& Z,
+ const vector<double>& R) {
// Select a free thread
unsigned sizeR = m_Ready.size();
unsigned int free_thread;
- bool found_thread=false;
- while(!found_thread){
- for(unsigned int i = 0 ; i < sizeR ; i++){
- if(!m_Ready[i]){
- free_thread=i;
- found_thread=true;
+ bool found_thread = false;
+ while (!found_thread) {
+ for (unsigned int i = 0; i < sizeR; i++) {
+ if (!m_Ready[i]) {
+ free_thread = i;
+ found_thread = true;
break;
}
}
}
- fitX[free_thread]=&X;
- fitZ[free_thread]=&Z;
- fitR[free_thread]=&R;
- m_uid[free_thread]=ID;
+ fitX[free_thread] = &X;
+ fitZ[free_thread] = &Z;
+ fitR[free_thread] = &R;
+ m_uid[free_thread] = ID;
// assume all X,Z,R of same size
sizeX[free_thread] = X.size();
- m_Ready[free_thread]=true;
+ m_Ready[free_thread] = true;
return;
}
////////////////////////////////////////////////////////////////////////////////
-void DCReconstructionMT::BuildTrack2D(){
- while(!IsDone())
+void DCReconstructionMT::BuildTrack2D() {
+ while (!IsDone())
std::this_thread::yield();
return;
}
////////////////////////////////////////////////////////////////////////////////
-bool DCReconstructionMT::IsDone(){
- std::vector<bool>::iterator begin = m_Ready.begin();
+bool DCReconstructionMT::IsDone() {
+ std::vector<bool>::iterator begin = m_Ready.begin();
std::vector<bool>::iterator end = m_Ready.end();
- for(auto it = begin ; it!=end ; it++){
- if((*it))
+ for (auto it = begin; it != end; it++) {
+ if ((*it))
return false;
}
return true;
}
////////////////////////////////////////////////////////////////////////////////
-void DCReconstructionMT::ResolvePlane(const TVector3& L,const double& ThetaU ,const TVector3& H, const double& ThetaV, TVector3& PosXY){
+void DCReconstructionMT::ResolvePlane(const TVector3& L, const double& ThetaU, const TVector3& H, const double& ThetaV,
+ TVector3& PosXY) {
// direction of U and V wire
- TVector3 u(0,1,0);
+ TVector3 u(0, 1, 0);
u.RotateZ(ThetaU);
- TVector3 v(0,1,0);
+ TVector3 v(0, 1, 0);
v.RotateZ(ThetaV);
-
// Compute the coeff of the two line of vecotr u (v) going through H (L)
// dv : y = av*x+bv
- av = v.Y()/v.X();
- bv = H.Y() - av*H.X();
+ av = v.Y() / v.X();
+ bv = H.Y() - av * H.X();
// du : y = au*x+bu
- au = u.Y()/u.X();
- bu = L.Y() - au*L.X();
+ au = u.Y() / u.X();
+ bu = L.Y() - au * L.X();
// We look for M(xM, yM) that intersect du and dv:
- if(!isinf(au) && !isinf(av)){ // au and av are not inf, i.e. not vertical line
- xM = (bv-bu)/(au-av);
- yM = au*xM+bu;
+ if (!std::isinf(au) && !std::isinf(av)) { // au and av are not inf, i.e. not vertical line
+ xM = (bv - bu) / (au - av);
+ yM = au * xM + bu;
}
- else if(isinf(av)){// av is inf, so v is along Y axis, H is direct measure of X
+ else if (std::isinf(av)) { // av is inf, so v is along Y axis, H is direct measure of X
xM = H.X();
- yM = au*xM+bu;
+ yM = au * xM + bu;
}
- else if (isinf(au)){//au is inf, so u is along Y axis, L is direct measure of X
+ else if (std::isinf(au)) { // au is inf, so u is along Y axis, L is direct measure of X
xM = L.X();
- yM = av*xM+bv;
+ yM = av * xM + bv;
}
- else{ // all is lost
- xM=-10000;
- yM=-10000;
+ else { // all is lost
+ xM = -10000;
+ yM = -10000;
}
- PosXY.SetXYZ(xM,yM,0);
+ PosXY.SetXYZ(xM, yM, 0);
}
-
////////////////////////////////////////////////////////////////////////////////
-double DCReconstructionMT::SumD(const double* parameter ){
+double DCReconstructionMT::SumD(const double* parameter) {
// Compute the sum P of the distance between the circle and the track
double P = 0;
double a = parameter[0];
double b = parameter[1];
- double ab= a*b;
- double a2=a*a;
+ double ab = a * b;
+ double a2 = a * a;
unsigned int id = parameter[2];
- unsigned int size = sizeX[id];
- const std::vector<double>* X=fitX[id];
- const std::vector<double>* Z=fitZ[id];
- const std::vector<double>* R=fitR[id];
+ unsigned int size = sizeX[id];
+ const std::vector<double>* X = fitX[id];
+ const std::vector<double>* Z = fitZ[id];
+ const std::vector<double>* R = fitR[id];
- double c,d,r,x,z,p;
- for(unsigned int i = 0 ; i < size ; i++){
+ double c, d, r, x, z, p;
+ for (unsigned int i = 0; i < size; i++) {
c = (*X)[i];
d = (*Z)[i];
r = (*R)[i];
- x = (a*d-ab+c)/(1+a2);
- z = a*x+b;
- p= (x-c)*(x-c)+(z-d)*(z-d)-r*r;
+ x = (a * d - ab + c) / (1 + a2);
+ z = a * x + b;
+ p = (x - c) * (x - c) + (z - d) * (z - d) - r * r;
// numerical trick to have a smooth derivative instead of using abs
- P+= sqrt(p*p+0.1);
+ P += sqrt(p * p + 0.1);
}
// return normalized power
- return P/size;
+ return P / size;
}
-
////////////////////////////////////////////////////////////////////////////////
-TGraph* DCReconstructionMT::Scan(double a, double b, int tovary, double minV, double maxV){
- vector<double> x,y;
- unsigned int sizeT=1000;
- double step = (maxV-minV)/sizeT;
- double p[2]={a,b};
- for(unsigned int i = 0 ; i < sizeT ; i++){
- if(!tovary){
- p[0]=minV+step*i;
+TGraph* DCReconstructionMT::Scan(double a, double b, int tovary, double minV, double maxV) {
+ vector<double> x, y;
+ unsigned int sizeT = 1000;
+ double step = (maxV - minV) / sizeT;
+ double p[2] = {a, b};
+ for (unsigned int i = 0; i < sizeT; i++) {
+ if (!tovary) {
+ p[0] = minV + step * i;
x.push_back(p[0]);
y.push_back(SumD(p));
}
- else{
- p[1]=minV+step*i;
+ else {
+ p[1] = minV + step * i;
x.push_back(p[1]);
y.push_back(SumD(p));
}
}
- TGraph* g = new TGraph(x.size(),&x[0],&y[0]);
+ TGraph* g = new TGraph(x.size(), &x[0], &y[0]);
return g;
}
////////////////////////////////////////////////////////////////////////////////
-void NPL::DCReconstructionMT::InitThreadPool(){
+void NPL::DCReconstructionMT::InitThreadPool() {
// this avoid error printout during fitting
gErrorIgnoreLevel = kError;
@@ -201,28 +198,28 @@ void NPL::DCReconstructionMT::InitThreadPool(){
StopThread();
m_ThreadPool.clear();
m_Ready.clear();
- m_Ready.resize(m_nbr_thread,false);
- for (unsigned int i=0; i < m_nbr_thread; i++) {
+ m_Ready.resize(m_nbr_thread, false);
+ for (unsigned int i = 0; i < m_nbr_thread; i++) {
// Register minimiser for futur deletion
- m_ThreadPool.push_back( std::thread(&NPL::DCReconstructionMT::StartThread,this,i) );
+ m_ThreadPool.push_back(std::thread(&NPL::DCReconstructionMT::StartThread, this, i));
}
m_stop = false;
- for(auto& th: m_ThreadPool){
+ for (auto& th : m_ThreadPool) {
th.detach();
}
}
////////////////////////////////////////////////////////////////////////////////
-void DCReconstructionMT::StartThread(unsigned int id){
+void DCReconstructionMT::StartThread(unsigned int id) {
// usefull variable
- double ai,bi;
+ double ai, bi;
unsigned int uid;
const double* xs;
- // create the functor
- // each threads needs its own or the minisation is not thread safe
- ROOT::Math::Functor* func= new ROOT::Math::Functor(this,&NPL::DCReconstructionMT::SumD,3);
- //Create the minimiser (deleted by the thread)
+ // create the functor
+ // each threads needs its own or the minisation is not thread safe
+ ROOT::Math::Functor* func = new ROOT::Math::Functor(this, &NPL::DCReconstructionMT::SumD, 3);
+ // Create the minimiser (deleted by the thread)
m_mtx.lock();
- ROOT::Math::Minimizer* mini=ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
+ ROOT::Math::Minimizer* mini = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
m_mtx.unlock();
mini->SetFunction(*func);
@@ -230,63 +227,64 @@ void DCReconstructionMT::StartThread(unsigned int id){
// Let the main thread start
std::this_thread::sleep_for(std::chrono::milliseconds(250));
- while(true){
+ while (true) {
// Do the job if possible
- if(m_Ready[id]){
- // Define the starting point of the fit: a straight line passing through the
+ if (m_Ready[id]) {
+ // 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 = ((*fitZ[id])[sizeX[id]-1]-(*fitZ[id])[0])/((*fitX[id])[sizeX[id]-1]-(*fitX[id])[0]+(*fitR[id])[sizeX[id]-1]-(*fitR[id])[0]);
- bi = (*fitZ[id])[0]-ai*((*fitX[id])[0]+(*fitR[id])[0]);
+ ai = ((*fitZ[id])[sizeX[id] - 1] - (*fitZ[id])[0]) /
+ ((*fitX[id])[sizeX[id] - 1] - (*fitX[id])[0] + (*fitR[id])[sizeX[id] - 1] - (*fitR[id])[0]);
+ bi = (*fitZ[id])[0] - ai * ((*fitX[id])[0] + (*fitR[id])[0]);
- if(isinf(ai)){ // then there is no two point in different layer
- m_a[uid]=-10000;
- m_b[uid]=-10000;
- m_X0[uid]=-10000;
- m_X100[uid]=-10000;
+ if (std::isinf(ai)) { // then there is no two point in different layer
+ m_a[uid] = -10000;
+ m_b[uid] = -10000;
+ m_X0[uid] = -10000;
+ m_X100[uid] = -10000;
m_minimum[uid] = 10000;
}
- else{
- mini->Clear();
- mini->SetVariable(0,"a",ai,1);
- mini->SetVariable(1,"b",bi,1);
- mini->SetFixedVariable(2,"id",id);
+ else {
+ mini->Clear();
+ mini->SetVariable(0, "a", ai, 1);
+ mini->SetVariable(1, "b", bi, 1);
+ mini->SetFixedVariable(2, "id", id);
// Perform minimisation
- mini->Minimize();
+ mini->Minimize();
// access set of parameter that produce the minimum
xs = mini->X();
- uid = m_uid[id];
- m_a[uid]=xs[0];
- m_b[uid]=xs[1];
- m_X0[uid]=-m_b[uid]/m_a[uid];
- m_X100[uid]=(100-m_b[uid])/m_a[uid];
+ uid = m_uid[id];
+ m_a[uid] = xs[0];
+ m_b[uid] = xs[1];
+ m_X0[uid] = -m_b[uid] / m_a[uid];
+ m_X100[uid] = (100 - m_b[uid]) / m_a[uid];
m_minimum[uid] = mini->MinValue();
}
// notify main thread job is done
- m_mtx.lock();// make sure no other thread is reading/writing to the map
+ m_mtx.lock(); // make sure no other thread is reading/writing to the map
m_Ready[id].flip();
m_mtx.unlock();
// Let other thread move up in the queu
std::this_thread::yield();
}
- else{
+ else {
std::this_thread::yield();
}
// Return if stopped
- if(m_stop){
+ if (m_stop) {
delete mini;
delete func;
return;
}
- }
+ }
}
////////////////////////////////////////////////////////////////////////////////
-void DCReconstructionMT::StopThread(){
+void DCReconstructionMT::StopThread() {
// make sure the last thread are schedule before stopping;
std::this_thread::yield();
- m_stop=true;
+ m_stop = true;
std::this_thread::yield();
}
#endif