diff --git a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx
index a7bdf08f0fe3d65a9aa693764c02f2328e190e6f..aa1b5889cac81fedb75728bca85776845cf9f5b0 100644
--- a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx
+++ b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx
@@ -64,13 +64,14 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
static map<std::pair<unsigned int,double>, vector<double> > X ;
static map<std::pair<unsigned int,double>, vector<double> > Z ;
static map<std::pair<unsigned int,double>, vector<double> > R ;
+ static int det,layer,wire;
X.clear();Z.clear();R.clear();
unsigned int size = Detector.size();
for(unsigned int i = 0 ; i < size ; i++){
if( DriftLength[i] > DriftLowThreshold && DriftLength[i] < DriftUpThreshold){
- int det = Detector[i];
- int layer = Layer[i];
- int wire = Wire[i];
+ det = Detector[i];
+ layer = Layer[i];
+ wire = Wire[i];
SamuraiDCIndex idx(det,layer,wire);
std::pair<unsigned int, double> p(det,Wire_Angle[idx]);
X[p].push_back(Wire_X[idx]);
@@ -83,16 +84,15 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
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();
+ static map<std::pair<unsigned int,double>, double > D ;// the minimum distance
+ VX0.clear();VX100.clear(),D.clear();
for(auto it = X.begin();it!=X.end();++it){
- 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
-
+ D[it->first]=m_reconstruction.BuildTrack2D(X[it->first],Z[it->first],R[it->first],X0,X100,a,b);
+
+ // very large a means track perpendicular to the chamber, what happen when there is pile up
if(abs(a)>1000)
PileUp++;
- MultPlane[it->first] = X[it->first].size() ;
Mult+=X[it->first].size();
// Position at z=0
TVector3 P(X0,0,0);
@@ -107,15 +107,21 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
// Reconstruct the central position (z=0) for each detector
static map<unsigned int,vector<TVector3> > C ;
- C.clear();
+ static map<unsigned int,vector<double> > W ; // weight based on D
+ C.clear(),W.clear();
TVector3 P;
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
m_reconstruction.ResolvePlane(it1->second,it1->first.second,it2->second,it2->first.second,P);
- if(P.X()!=-10000)
+
+ if(P.X()!=-10000){
C[it1->first.first].push_back(P);
+ // Mean pos are weighted based on the the sum of distance from track
+ // to hit obtained during the minimisation
+ W[it1->first.first].push_back(1./(D[it1->first]+D[it2->first]));
+ }
}
}
}
@@ -136,37 +142,39 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
// Build the Reference position by averaging all possible pair
size = C[2].size();
- double PosX100,PosY100;
+ static double PosX100,PosY100,norm;
if(size){
PosX=0;
PosY=0;
PosX100=0;
PosY100=0;
+ norm=0;
for(unsigned int i = 0 ; i < size ; i++){
- PosX+= C[2][i].X();
- PosY+= C[2][i].Y();
- PosX100+= C100[2][i].X();
- PosY100+= C100[2][i].Y();
+ PosX+= C[2][i].X()*W[2][i];
+ PosY+= C[2][i].Y()*W[2][i];
+ PosX100+= C100[2][i].X()*W[2][i];
+ PosY100+= C100[2][i].Y()*W[2][i];
+ norm+=W[2][i];
// cout << C[2][i].X() << " (" << C[2][i].Y() << ") ";
}
// cout << endl;
MultMean=size;
// Mean position at Z=0
- PosX=PosX/size;
- PosY=PosY/size;
+ PosX=PosX/norm;
+ PosY=PosY/norm;
// Mean position at Z=100
- PosX100=PosX100/size;
- PosY100=PosY100/size;
+ PosX100=PosX100/norm;
+ PosY100=PosY100/norm;
devX=0;
devY=0;
for(unsigned int i = 0 ; i < size ; i++){
- devX+=(C[2][i].X()-PosX)*(C[2][i].X()-PosX);
- devY+=(C[2][i].Y()-PosY)*(C[2][i].Y()-PosY);
+ devX+=W[2][i]*(C[2][i].X()-PosX)*(C[2][i].X()-PosX);
+ devY+=W[2][i]*(C[2][i].Y()-PosY)*(C[2][i].Y()-PosY);
}
- devX=sqrt(devX/size);
- devY=sqrt(devY/size);
+ devX=sqrt(devX/((size-1)*norm));
+ devY=sqrt(devY/((size-1)*norm));
// Compute ThetaX, angle between the Direction vector projection in XZ with
// the Z axis
diff --git a/NPLib/TrackReconstruction/NPDCReconstruction.cxx b/NPLib/TrackReconstruction/NPDCReconstruction.cxx
index 29b633410fd572c2cd37c2a30a8fe3143b77a189..7f8337a24147ff9886d1c6dc6ee542cddaebd5cb 100644
--- a/NPLib/TrackReconstruction/NPDCReconstruction.cxx
+++ b/NPLib/TrackReconstruction/NPDCReconstruction.cxx
@@ -33,21 +33,22 @@ DCReconstruction::~DCReconstruction(){
}
////////////////////////////////////////////////////////////////////////////////
-void DCReconstruction::BuildTrack2D(const vector<double>& X,const vector<double>& Z,const vector<double>& R,double& X0,double& X100,double& a, double& b ){
+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
- unsigned int size = X.size();
+ 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};
+ ai = (Z[size-1]-Z[0])/(X[size-1]-R[size-1]-X[0]-R[0]);
+ bi = Z[0]-ai*(X[0]+R[0]);
+ parameter[0]=ai;
+ parameter[1]=bi;
m_min->SetVariable(0,"a",parameter[0],1000);
m_min->SetVariable(1,"b",parameter[1],1000);
- m_min->SetTolerance(0.1);
+ m_min->SetTolerance(0.01);
// Perform minimisation
m_min->Minimize();
@@ -58,29 +59,29 @@ void DCReconstruction::BuildTrack2D(const vector<double>& X,const vector<double>
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){
// direction of U and V wire
- TVector3 u = TVector3(0,1,0);
+ TVector3 u(0,1,0);
u.RotateZ(ThetaU);
- TVector3 v = TVector3(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
- long double av = v.Y()/v.X();
- long double bv = H.Y() - av*H.X();
+ av = v.Y()/v.X();
+ bv = H.Y() - av*H.X();
// du : y = au*x+bu
- long double au = u.Y()/u.X();
- long double 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:
- 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;
@@ -98,7 +99,7 @@ void DCReconstruction::ResolvePlane(const TVector3& L,const double& ThetaU ,cons
yM=-10000;
}
- PosXY=TVector3(xM,yM,0);
+ PosXY.SetXYZ(xM,yM,0);
}
@@ -107,12 +108,11 @@ void DCReconstruction::ResolvePlane(const TVector3& L,const double& ThetaU ,cons
double DCReconstruction::SumD(const double* parameter ){
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;
+ P = 0;
+ a = parameter[0];
+ b = parameter[1];
+ ab= a*b;
+ a2=a*a;
for(unsigned int i = 0 ; i < size ; i++){
c = (*fitX)[i];
@@ -120,9 +120,10 @@ double DCReconstruction::SumD(const double* parameter ){
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;
+ P+= sqrt(abs( (x-c)*(x-c)+(z-d)*(z-d)-r*r));
}
- return P;
+ // return normalized power
+ return P/size;
}
diff --git a/NPLib/TrackReconstruction/NPDCReconstruction.h b/NPLib/TrackReconstruction/NPDCReconstruction.h
index 830feb6517ca6b313980848652c7cb089efc66f5..5659aaf2cbd3c30329472d3f779c45ee7e37f76d 100644
--- a/NPLib/TrackReconstruction/NPDCReconstruction.h
+++ b/NPLib/TrackReconstruction/NPDCReconstruction.h
@@ -39,23 +39,23 @@
#include "Math/Minimizer.h"
#include "Math/Functor.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 );
+ // 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
+ double 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
ROOT::Math::Minimizer* m_min;
@@ -63,7 +63,18 @@ namespace NPL{
const std::vector<double>* fitX;
const std::vector<double>* fitZ;
const std::vector<double>* fitR;
- };
- }
+ // used by SumD
+ unsigned int size ;
+ double P,a,b,ab,a2,c,d,x,z,r;
+ // used by BuildTrack
+ double ai,bi;
+ double parameter[2];
+ // used by resolve plane
+ long double av,bv,au,bu;
+ double xM,yM;
+
+
+ };
+}
#endif