* Now putting the DC reconstruction into the TrackReconstruction

Library

NPLib/Detectors/Samurai/CMakeLists.txt

@@ -4,6 +4,6 @@ add_custom_command(OUTPUT TSamuraiFDC2PhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR
 
 add_library(NPSamurai SHARED TSamuraiFDC2Data.cxx TSamuraiFDC2DataDict.cxx TSamuraiFDC2Physics.cxx TSamuraiFDC2PhysicsDict.cxx)
 
-target_link_libraries(NPSamurai ${ROOT_LIBRARIES} NPCore) 
+target_link_libraries(NPSamurai ${ROOT_LIBRARIES} NPCore NPTrackReconstruction) 
 install(FILES TSamuraiFDC2Data.h TSamuraiFDC2Physics.h  DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
 

NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx

@@ -35,11 +35,6 @@
 #include "NPDetectorFactory.h"
 #include "NPSystemOfUnits.h"
 //   ROOT
-#include "TChain.h"
-#include "TVector2.h"
-#include "Math/Minimizer.h"
-#include "Math/Functor.h"
-#include "Math/Factory.h"
 using namespace NPUNITS;
 ///////////////////////////////////////////////////////////////////////////
 
@@ -83,13 +78,13 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
   }
 
   // Reconstruct the vector for each of the plane of each of the detector
-  static double X0,X100;
+  static double X0,X100,a,b; // store the BuildTrack2D results
   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>, int > MultPlane ;  
   VX0.clear();VX100.clear();
   for(auto it = X.begin();it!=X.end();++it){
-    double a = Track2D(X[it->first],Z[it->first],R[it->first],X0,X100); 
+    m_reconstruction.BuildTrack2D(X[it->first],Z[it->first],R[it->first],X0,X100,a,b); 
     // very small a means track perpendicular to the chamber, what happen when there is pile up
 
     if(abs(a)>1000)
@@ -115,7 +110,7 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
   for(auto it1 = VX0.begin();it1!=VX0.end();++it1){
     for(auto it2 = it1;it2!=VX0.end();++it2){
       if(it1!=it2 && it1->first.first==it2->first.first){// different plane, same detector
-        ResolvePlane(it1->second,it1->first.second,it2->second,it2->first.second,P);
+        m_reconstruction.ResolvePlane(it1->second,it1->first.second,it2->second,it2->first.second,P);
         if(P.X()!=-10000)
           C[it1->first.first].push_back(P);
       }
@@ -127,7 +122,7 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
   for(auto it1 = VX100.begin();it1!=VX100.end();++it1){
     for(auto it2 = it1;it2!=VX100.end();++it2){
       if(it1!=it2 && it1->first.first==it2->first.first){// different plane, same detector
-        ResolvePlane(it1->second,it1->first.second,it2->second,it2->first.second,P);
+        m_reconstruction.ResolvePlane(it1->second,it1->first.second,it2->second,it2->first.second,P);
 
         if(P.X()!=-10000)
           C100[it1->first.first].push_back(P);
@@ -168,119 +163,6 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
 
   return;
 }
-////////////////////////////////////////////////////////////////////////////////
-void TSamuraiFDC2Physics::ResolvePlane(const TVector3& L,const double& ThetaU ,const TVector3& H, const double& ThetaV, TVector3& PosXY){
-  // direction of U and V wire 
-  TVector3 u = TVector3(0,1,0);
-  u.RotateZ(ThetaU); 
-
-  TVector3 v = TVector3(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
-  double av = v.Y()/v.X();
-  double bv = H.Y() - av*H.X();   
-  
-  // du : y = au*x+bv
-  double au = u.Y()/u.X();
-  double bu = L.Y() - au*L.X();   
-
-  // We look for M(xM, yM) that intersect du and dv:
-  double xM,yM;
-  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=TVector3(xM,yM,0);
-};
-///////////////////////////////////////////////////////////////////////////
-double TSamuraiFDC2Physics::Track2D(const vector<double>& X,const vector<double>& Z,const vector<double>& R,double& X0,double& X100 ){
-  fitX=X;
-  fitZ=Z;
-  fitR=R;
-  // assume all X,Z,R of same size
-  unsigned int size = X.size();
-/*
-    ofstream out("data.txt");
-      for(unsigned int i = 0 ; i < size ; i++){
-      out << fitX[i] << " " << fitZ[i] << " "<< fitR[i] << endl;
-      }
-      out.close();
-  */      
-  // 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
-  double a = (fitZ[size-1]-fitZ[0])/(fitX[size-1]-fitR[size-1]-fitX[0]-fitR[0]);
-  double b = fitZ[0]-a*(fitX[0]+fitR[0]);
-  double parameter[2]={a,b};
-/*    out.open("line.txt");
-      out << a << " " << b <<endl;
-      out.close();
-      */
-  //cout << a << " - " << b << " " << SumD(parameter) <<endl;
-
-  static ROOT::Math::Minimizer* min = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
-  static ROOT::Math::Functor f(this,&TSamuraiFDC2Physics::SumD,2); 
-  min->SetFunction(f);
-  min->SetVariable(0,"a",parameter[0],1000);
-  min->SetVariable(1,"b",parameter[1],1000);
-  min->SetTolerance(0.1);
-  // cout <<"start " << endl;
-  min->Minimize(); 
-  //cout << "start" << endl;
-  const double *xs = min->X();
-  X0=-xs[1]/xs[0];
-  X100=(100-xs[1])/xs[0];
-/*  
-     out.open("fit.txt");
-     out << xs[0] << " " << xs[1] <<endl;
-     out.close();
-     if(X0==0){
-     exit(0);
-     }
-  */   
-  return xs[0];
-  //  cout << " done " << SumD(xs) <<endl;  
-  //cout << xs[0] << " / " << xs[1] << " " << SumD(xs) <<endl;
-}
-////////////////////////////////////////////////////////////////////////////////
-double TSamuraiFDC2Physics::SumD(const double* parameter){
-  //cout << " "<<parameter[0] << " h " << parameter[1] ;
-  unsigned int size = fitX.size();
-  // 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 c,d,x,z;
-
-  for(unsigned int i = 0 ; i < size ; i++){
-    c = fitX[i];
-    d = fitZ[i];
-    x = (a*d-ab+c)/(1+a2);
-    z = a*x+b;
-    //P+= sqrt(abs( (x-c)*(x-c)+(z-d)*(z-d)-fitR[i]*fitR[i]))/fitR[i];
-    P+= abs( (x-c)*(x-c)+(z-d)*(z-d)-fitR[i]*fitR[i])/fitR[i];
-    //  cout << sqrt(abs((x-c)*(x-c)+(z-d)*(z-d)-fitR[i]*fitR[i] ))/fitR[i]<< " ";
-  }
-  //cout << endl;
-  return P;
-}
 ///////////////////////////////////////////////////////////////////////////
 void TSamuraiFDC2Physics::PreTreat(){
   ClearPreTreatedData();

NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h

@@ -33,12 +33,9 @@
 #include "NPVDetector.h"
 #include "NPInputParser.h"
 #include "NPXmlParser.h"
+#include "NPDCReconstruction.h"
 // ROOT 
-#include "TVector2.h" 
 #include "TVector3.h" 
-#include "TObject.h"
-#include "TCanvas.h"
-#include "TRandom3.h"
 // Forward declaration
 //class TSamuraiFDC2Spectra;
 
@@ -106,6 +103,7 @@ class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{
     int Mult;
     int MultMean;
     int PileUp;
+
   public:
     // Projected position at given Z plan
     TVector3 ProjectedPosition(double Z);
@@ -121,13 +119,8 @@ class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{
     void RemoveNoise();
     // Construct the 2D track and ref position at Z=0 and Z=100 based on X,Z and Radius provided
 
-    double Track2D(const vector<double>& X,const vector<double>& Z,const vector<double>& R,double& X0,double& X100 );
-    // Compute X and Y of interaction point based on drift vector of two different wire plane
-    void ResolvePlane(const TVector3& PosU,const double& ThetaU ,const TVector3& PosV, const double& ThetaV, TVector3& PosXY);
-    double SumD(const double* parameter );
-    vector<double> fitX;
-    vector<double> fitZ;
-    vector<double> fitR;
+    // Object use to perform the DC reconstruction
+    NPL::DCReconstruction m_reconstruction;//!
 
   public: //   Innherited from VDetector Class
 

NPLib/TrackReconstruction/CMakeLists.txt

@@ -4,7 +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)
+add_library(NPTrackReconstruction SHARED NPRansac.cxx NPCluster.cxx NPTrack.cxx Tracking.cxx NPRansacDict.cxx NPClusterDict.cxx TrackingDict.cxx NPDCReconstruction.cxx)
+
 target_link_libraries(NPTrackReconstruction ${ROOT_LIBRARIES} NPCore)
 
-install(FILES NPRansac.h NPCluster.h NPTrack.h Tracking.h NPTrackingUtility.h   DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+install(FILES NPRansac.h NPCluster.h NPTrack.h Tracking.h NPTrackingUtility.h NPDCReconstruction.h  DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})

NPLib/TrackReconstruction/NPDCReconstruction.cxx

+#include"NPDCReconstruction.h"
+#include "Math/Minimizer.h"
+#include "Math/Functor.h"
+#include "Math/Factory.h"
+
+using namespace std;
+using namespace NPL;
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+void 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
+  unsigned int size = X.size();
+  // 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
+  double ai = (Z[size-1]-Z[0])/(X[size-1]-R[size-1]-X[0]-R[0]);
+  double bi = Z[0]-ai*(X[0]+R[0]);
+  double parameter[2]={ai,bi};
+
+  static ROOT::Math::Minimizer* min = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
+  static ROOT::Math::Functor f(this,&NPL::DCReconstruction::SumD,2); 
+  min->SetFunction(f);
+  min->SetVariable(0,"a",parameter[0],1000);
+  min->SetVariable(1,"b",parameter[1],1000);
+  min->SetTolerance(0.1);
+  // Perform minimisation
+  min->Minimize(); 
+  // access set of parameter that produce the minimum
+  const double *xs = min->X();
+  a=xs[0];
+  b=xs[1];
+  X0=-b/a;
+  X100=(100-b)/a;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+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 = TVector3(0,1,0);
+  u.RotateZ(ThetaU);
+
+  TVector3 v = TVector3(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
+  double av = v.Y()/v.X();
+  double bv = H.Y() - av*H.X();
+
+  // du : y = au*x+bv
+  double au = u.Y()/u.X();
+  double bu = L.Y() - au*L.X();
+
+  // We look for M(xM, yM) that intersect du and dv:
+  double xM,yM;
+  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=TVector3(xM,yM,0);
+
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+double DCReconstruction::SumD(const double* parameter ){
+  //cout << " "<<parameter[0] << " h " << parameter[1] ;
+  unsigned int size = fitX->size();
+  // 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 c,d,x,z,r;
+
+  for(unsigned int i = 0 ; i < size ; i++){
+    c = (*fitX)[i];
+    d = (*fitZ)[i];
+    r = (*fitR)[i];
+    x = (a*d-ab+c)/(1+a2);
+    z = a*x+b;
+    P+= abs( (x-c)*(x-c)+(z-d)*(z-d)-r*r)/r;
+  }
+  return P;
+
+}
+
+

NPLib/TrackReconstruction/NPDCReconstruction.h

+#ifndef NPDCRECONSTRUCTION_H
+#define NPDCRECONSTRUCTION_H
+/*****************************************************************************
+ * 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 :  Pierre MORFOUACE 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<vector>
+#include"TVector3.h"
+namespace NPL{
+  
+  class DCReconstruction{
+    public:
+      DCReconstruction(){};
+      ~DCReconstruction(){};
+    
+    public:
+    // 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
+    void BuildTrack2D(const std::vector<double>& X,const std::vector<double>& Z,const std::vector<double>& R,double& X0,double& X100,double& a, double& b );
+    
+    // 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 );
+
+    private: // private member used by SumD
+      const std::vector<double>* fitX;
+      const std::vector<double>* fitZ;
+      const std::vector<double>* fitR;
+    };
+  }
+
+#endif