Commit 4adf68cf authored by Adrien Matta's avatar Adrien Matta
Browse files

* adding Multithreading to DC reconstruction

parent 5450012f
Pipeline #92302 failed with stages
in 6 minutes and 1 second
......@@ -33,8 +33,6 @@
#include <iostream>
#if __cplusplus > 199711L
#include <thread>
#include <mutex>
#include <condition_variable>
#endif
......
......@@ -50,6 +50,14 @@ ClassImp(TSamuraiFDC0Physics)
DriftLowThreshold=0.1 ;
DriftUpThreshold=2.4;
PowerThreshold=5;
#if __cplusplus > 199711L && NPMULTITHREADING
// one thread for each plan X,Y = 2
// ! more than that this will not help !
m_reconstruction.SetNumberOfThread(2);
m_reconstruction.InitThreadPool();
#endif
}
///////////////////////////////////////////////////////////////////////////
......@@ -88,9 +96,26 @@ void TSamuraiFDC0Physics::BuildPhysicalEvent(){
static map<std::pair<unsigned int,double>, TVector3 > VX0 ;
static map<std::pair<unsigned int,double>, TVector3 > VX100 ;
static map<std::pair<unsigned int,double>, double > D ;// the minimum distance
static unsigned int uid=0;
VX0.clear();VX100.clear(),D.clear();
for(auto it = X.begin();it!=X.end();++it){
#if __cplusplus > 199711L && NPMULTITHREADING
m_reconstruction.AddPlan(uid++,X[it->first],Z[it->first],R[it->first]);
#else
D[it->first]=m_reconstruction.BuildTrack2D(X[it->first],Z[it->first],R[it->first],X0,X100,a,b);
#endif
}
#if __cplusplus > 199711L && NPMULTITHREADING
// do all plan at once in parallele, return when all plan are done
m_reconstruction.BuildTrack2D();
uid=0;
#endif
for(auto it = X.begin();it!=X.end();++it){
#if __cplusplus > 199711L && NPMULTITHREADING
D[it->first]=m_reconstruction.GetResults(uid++,X0,X100,a,b);
#endif
// for Debug, write a file of
/* { std::ofstream f("distance.txt", std::ios::app);
......
......@@ -34,7 +34,13 @@
#include "NPVDetector.h"
#include "NPInputParser.h"
#include "NPXmlParser.h"
#if __cplusplus > 199711L && NPMULTITHREADING
#include "NPDCReconstructionMT.h"
#else
#include "NPDCReconstruction.h"
#endif
// ROOT
#include "TVector3.h"
......@@ -97,7 +103,11 @@ class TSamuraiFDC0Physics : public TObject, public NPL::VDetector{
// Construct the 2D track and ref position at Z=0 and Z=100 based on X,Z and Radius provided
// Object use to perform the DC reconstruction
#if __cplusplus > 199711L && NPMULTITHREADING
NPL::DCReconstructionMT m_reconstruction;//!
#else
NPL::DCReconstruction m_reconstruction;//!
#endif
public: // Innherited from VDetector Class
......
......@@ -50,6 +50,13 @@ ClassImp(TSamuraiFDC2Physics)
DriftLowThreshold=0.4 ;
DriftUpThreshold=9.3;
PowerThreshold=5;
#if __cplusplus > 199711L && NPMULTITHREADING
// one thread for each plan X,U,V = 3
// ! more than this will not help !
m_reconstruction.SetNumberOfThread(3);
m_reconstruction.InitThreadPool();
#endif
}
///////////////////////////////////////////////////////////////////////////
......@@ -87,9 +94,26 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
static map<std::pair<unsigned int,double>, TVector3 > VX0 ;
static map<std::pair<unsigned int,double>, TVector3 > VX100 ;
static map<std::pair<unsigned int,double>, double > D ;// the minimum distance
static unsigned int uid=0;
VX0.clear();VX100.clear(),D.clear();
for(auto it = X.begin();it!=X.end();++it){
#if __cplusplus > 199711L && NPMULTITHREADING
m_reconstruction.AddPlan(uid++,X[it->first],Z[it->first],R[it->first]);
#else
D[it->first]=m_reconstruction.BuildTrack2D(X[it->first],Z[it->first],R[it->first],X0,X100,a,b);
#endif
}
#if __cplusplus > 199711L && NPMULTITHREADING
// do all plan at once in parallele, return when all plan are done
m_reconstruction.BuildTrack2D();
uid=0;
#endif
for(auto it = X.begin();it!=X.end();++it){
#if __cplusplus > 199711L && NPMULTITHREADING
D[it->first]=m_reconstruction.GetResults(uid++,X0,X100,a,b);
#endif
/* // for Debug, write a file of
{ std::ofstream f("distance.txt", std::ios::app);
......@@ -110,7 +134,6 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
TVector3 P100= TVector3(X100,0,0);
P100.RotateZ(it->first.second);
VX100[it->first]=P100;
}
// Reconstruct the central position (z=0) for each detector
......
......@@ -34,7 +34,13 @@
#include "NPVDetector.h"
#include "NPInputParser.h"
#include "NPXmlParser.h"
#if __cplusplus > 199711L && NPMULTITHREADING
#include "NPDCReconstructionMT.h"
#else
#include "NPDCReconstruction.h"
#endif
// ROOT
#include "TVector3.h"
// Forward declaration
......@@ -96,7 +102,13 @@ class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{
// Construct the 2D track and ref position at Z=0 and Z=100 based on X,Z and Radius provided
// Object use to perform the DC reconstruction
#if __cplusplus > 199711L && NPMULTITHREADING
NPL::DCReconstructionMT m_reconstruction;//!
#else
NPL::DCReconstruction m_reconstruction;//!
#endif
public: // Innherited from VDetector Class
......
......@@ -4,8 +4,8 @@ add_custom_command(OUTPUT NPClusterDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/
add_custom_command(OUTPUT TrackingDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh Tracking.h TrackingDict.cxx Tracking.rootmap libNPTrackReconstruction.so NPTrackReconstructionLinkDef.h DEPENDS Tracking.h)
add_library(NPTrackReconstruction SHARED NPRansac.cxx NPCluster.cxx NPTrack.cxx Tracking.cxx NPRansacDict.cxx NPClusterDict.cxx TrackingDict.cxx NPDCReconstruction.cxx)
add_library(NPTrackReconstruction SHARED NPRansac.cxx NPCluster.cxx NPTrack.cxx Tracking.cxx NPRansacDict.cxx NPClusterDict.cxx TrackingDict.cxx NPDCReconstruction.cxx NPDCReconstructionMT.cxx)
target_link_libraries(NPTrackReconstruction ${ROOT_LIBRARIES} NPCore)
target_link_libraries(NPTrackReconstruction ${ROOT_LIBRARIES} -lMinuit2 NPCore)
install(FILES NPRansac.h NPCluster.h NPTrack.h Tracking.h NPTrackingUtility.h NPDCReconstruction.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
install(FILES NPRansac.h NPCluster.h NPTrack.h Tracking.h NPTrackingUtility.h NPDCReconstruction.h NPDCReconstructionMT.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
#if __cplusplus > 199711L // require c++11
/*****************************************************************************
* Copyright (C) 2009-2020 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 : Adrien Matta contact address: matta@lpcaen.in2p3.fr *
* *
* Creation Date : October 2020 *
* Last update : October 2020 *
*---------------------------------------------------------------------------*
* Decription: *
* This class have all the method needed to analyse Drift Chambers *
*****************************************************************************/
#include"NPDCReconstructionMT.h"
#include "Math/Factory.h"
#include "TError.h"
#include "TGraph.h"
using namespace std;
using namespace NPL;
////////////////////////////////////////////////////////////////////////////////
DCReconstructionMT::DCReconstructionMT(unsigned int number_thread){
// this avoid error printout
gErrorIgnoreLevel = kError;
m_nbr_thread= number_thread;
}
////////////////////////////////////////////////////////////////////////////////
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];
return m_minimum[ID];
}
return -1;
}
////////////////////////////////////////////////////////////////////////////////
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;
break;
}
}
}
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;
return;
}
////////////////////////////////////////////////////////////////////////////////
void DCReconstructionMT::BuildTrack2D(){
while(!IsDone())
std::this_thread::yield();
return;
}
////////////////////////////////////////////////////////////////////////////////
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))
return false;
}
return true;
}
////////////////////////////////////////////////////////////////////////////////
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);
u.RotateZ(ThetaU);
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();
// du : y = au*x+bu
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;
}
else if(isinf(av)){// av is inf, so v is along Y axis, H is direct measure of X
xM = H.X();
yM = au*xM+bu;
}
else if (isinf(au)){//au is inf, so u is along Y axis, L is direct measure of X
xM = L.X();
yM = av*xM+bv;
}
else{ // all is lost
xM=-10000;
yM=-10000;
}
PosXY.SetXYZ(xM,yM,0);
}
////////////////////////////////////////////////////////////////////////////////
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;
unsigned int id = parameter[2];
double c,d,r,x,z,p;
for(unsigned int i = 0 ; i < sizeX[id] ; i++){
c = (*fitX[id])[i];
d = (*fitZ[id])[i];
r = (*fitR[id])[i];
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);
}
// return normalized power
return P/sizeX[id];
}
////////////////////////////////////////////////////////////////////////////////
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;
x.push_back(p[1]);
y.push_back(SumD(p));
}
}
TGraph* g = new TGraph(x.size(),&x[0],&y[0]);
return g;
}
////////////////////////////////////////////////////////////////////////////////
void NPL::DCReconstructionMT::InitThreadPool(){
StopThread();
m_ThreadPool.clear();
m_Ready.clear();
for (unsigned int i=0; i < m_nbr_thread; i++) {
//Create the minimiser (deleted by the thread)
ROOT::Math::Minimizer* mini=ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
// Register minimiser for futur deletion
m_ThreadPool.push_back( std::thread(&NPL::DCReconstructionMT::StartThread,this,mini,i) );
m_Ready.push_back(false);
}
m_stop = false;
for(auto& th: m_ThreadPool){
th.detach();
}
}
////////////////////////////////////////////////////////////////////////////////
void DCReconstructionMT::StartThread(ROOT::Math::Minimizer* mini,unsigned int id){
// Let the main thread start
std::this_thread::sleep_for(std::chrono::milliseconds(100));
// 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);
mini->SetFunction(*func);
mini->SetPrintLevel(0);
double ai,bi;
unsigned int uid;
const double* xs;
while(true){
// Do the job if possible
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]);
bi = (*fitZ[id])[0]-ai*((*fitX[id])[0]);
mini->Clear();
mini->SetVariable(0,"a",ai,1);
mini->SetVariable(1,"b",bi,1);
mini->SetFixedVariable(2,"id",id);
// Perform minimisation
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];
m_minimum[uid] = mini->MinValue();
// notify main thread job is done
m_Ready[id].flip();
// Let other thread move up in the queu
std::this_thread::yield();
}
else{
std::this_thread::yield();
}
// Return if stopped
if(m_stop){
delete mini;
delete func;
return;
}
}
}
////////////////////////////////////////////////////////////////////////////////
void DCReconstructionMT::StopThread(){
// make sure the last thread are schedule before stopping;
std::this_thread::yield();
m_stop=true;
std::this_thread::yield();
}
#endif
#ifndef NPDCRECONSTRUCTIONMT_H
#define NPDCRECONSTRUCTIONMT_H
#if __cplusplus > 199711L // require c++11
/*****************************************************************************
* Copyright (C) 2009-2020 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 : Adrien Matta contact address: matta@lpcaen.in2p3.fr *
* *
* Creation Date : October 2020 *
* Last update : October 2020 *
*---------------------------------------------------------------------------*
* Decription: *
* This class have all the method needed to analyse Drift Chambers *
* Notation: *
* - The Wire are along the Y axis, there fore providing a measure of the *
* X positions *
* - First we find a 2D track in the X-Z plane based on the measurement of*
* several drift radius by minimisation: *
* - Let be Di the distance between a trial track and one of the drift *
* circle of coordinate Xi Zi and radius Ri. *
* - We compute for the track P= sum_i(Di^2/Ri) *
* - We minimise P to find the correct track *
* - Warning : this algo assume only one track *
* - Once the track found, we return the X position at Z=0 and Z=100 to *
* the user. *
* - User can rotate found position depending on the real position of the *
* Wires. *
* - Once the position found for all wire plan, one can look for the *
* intersaction in any given plane. This is done using ResolvePlane. *
* - Resolving plane for two Z plane will provide a point and a direction *
*****************************************************************************/
#include <vector>
#include <map>
#include "TVector3.h"
#include "Math/Minimizer.h"
#include "Math/Functor.h"
#include <thread>
class TGraph;
namespace NPL{
class DCReconstructionMT{
public:
DCReconstructionMT(unsigned int number_thread=1);
~DCReconstructionMT();
public:
// Set number of thread
// require to stop and reinit the thread
void SetNumberOfThread(unsigned int number_thread){m_nbr_thread=number_thread;};
void AddPlan(unsigned int ID,const std::vector<double>& X, const std::vector<double>& Z, const std::vector<double>& R);
double GetResults(unsigned int ID,double& X0,double& X100,double& a, double& b);
// Build a track in 2D based on drift circle of Radius R and position X,Z
// return X0(X100) the X position at Z=0 (Z=100)
// return a and b the coeff of the 2D line
// when all thread are done
void BuildTrack2D();
// Compute X and Y crossing coordinate of 2 plane of Wire
void ResolvePlane(const TVector3& L,const double& ThetaU ,const TVector3& H, const double& ThetaV, TVector3& PosXY);
// Function used by the minimizer in BuildTrack2D
double SumD(const double* parameter );
// For debugging/optimisation
// Scan Sumd versus parameter a or b (tovary =0 for a, 1 for b)
// return a TGraph for display
TGraph* Scan(double a, double b, int tovary, double minV, double maxV);
private: // private member used by SumD
// data to minize index by thread ID
std::map<unsigned int,unsigned int> sizeX;
std::map<unsigned int,unsigned int> m_uid; // match thread id and user id
std::map<unsigned int,const std::vector<double>*> fitX;
std::map<unsigned int,const std::vector<double>*> fitZ;
std::map<unsigned int,const std::vector<double>*> fitR;
// Computed value indexed by user ID
std::map<unsigned int,double> m_minimum;
std::map<unsigned int,double> m_X0;
std::map<unsigned int,double> m_X100;
std::map<unsigned int,double> m_a;
std::map<unsigned int,double> m_b;
private: // Thread Pool defined if C++11 is available
unsigned int m_nbr_thread;
std::vector<std::thread> m_ThreadPool;
std::vector<bool> m_Ready;
bool m_stop;
public: // Init the Thread Pool
void StopThread();
void StartThread(ROOT::Math::Minimizer*,unsigned int);
void InitThreadPool();
bool IsDone();
// used by SumD
// used by resolve plane
long double av,bv,au,bu;
double xM,yM;
};
}
#endif//c++11
#endif//ndef
......@@ -1560,7 +1560,7 @@
<geo>16</geo>
<ch>7</ch>
</SAMURAIFDC0>
<!--<SAMURAIFDC0>
<SAMURAIFDC0>
<ID>105</ID>
<NAME>FDC0_3_8</NAME>
<FPL>13</FPL>
......@@ -1604,7 +1604,7 @@
<det>28</det>
<geo>16</geo>
<ch>10</ch>
</SAMURAIFDC0>-->
</SAMURAIFDC0>
<SAMURAIFDC0>
<ID>108</ID>
<NAME>FDC0_3_11</NAME>
......
Markdown is supported
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