From 0ad4ac94889b065d8f45d124f95469ce5c268a81 Mon Sep 17 00:00:00 2001
From: Adrien Matta <matta@lpccaen.in2p3.fr>
Date: Fri, 21 May 2021 09:05:12 +0200
Subject: [PATCH] * Improving DC reconstruction algorithm - In many
cases the first and last wire had the same X causing the initial
value to be inf and the fit to return nan
---
.../Detectors/Samurai/TSamuraiBDCPhysics.cxx | 19 ++++----
.../NPDCReconstruction.cxx | 42 ++++++++++-------
.../NPDCReconstructionMT.cxx | 47 ++++++++++++-------
3 files changed, 66 insertions(+), 42 deletions(-)
diff --git a/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx b/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx
index 7c94e762b..4d22d3d42 100644
--- a/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx
+++ b/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx
@@ -48,8 +48,8 @@ ClassImp(TSamuraiBDCPhysics)
//m_Spectra = NULL;
ToTThreshold_L = 0;
ToTThreshold_H = 1000;
- DriftLowThreshold=0.1 ;
- DriftUpThreshold=2.4;
+ DriftLowThreshold=0 ;
+ DriftUpThreshold=2.5;
PowerThreshold=5;
}
@@ -75,6 +75,7 @@ void TSamuraiBDCPhysics::BuildPhysicalEvent(){
static unsigned int uid; uid=0;
static vector<TVector3> C ;
static vector<double > W ; // weight based on D
+ static double PosX100,PosY100,norm;
unsigned int count = 0 ;
for(auto it = m_DCHit.begin(); it!=m_DCHit.end(); it++){
// Each entry in the map is a detector
@@ -134,7 +135,7 @@ void TSamuraiBDCPhysics::BuildPhysicalEvent(){
// very large "a" means track perpendicular to the chamber, what happen when there is pile up
if(abs(a)>5000)
PileUp[count]++;
-
+//cout << a << " " << b << endl;
// Position at z=0
TVector3 P(X0,0,0);
P.RotateZ(it->first);
@@ -177,9 +178,8 @@ void TSamuraiBDCPhysics::BuildPhysicalEvent(){
}
// Build the Reference position by averaging all possible pair
size = C.size();
- static double PosX100,PosY100,norm;
if(size){
- norm=0;
+ norm=0;PosX100=0;PosY100=0;
for(unsigned int i = 0 ; i < size ; i++){
PosX[count]+= C[i].X()*W[i];
PosY[count]+= C[i].Y()*W[i];
@@ -199,19 +199,20 @@ void TSamuraiBDCPhysics::BuildPhysicalEvent(){
devX[count]+=W[i]*(C[i].X()-PosX[count])*(C[i].X()-PosX[count]);
devY[count]+=W[i]*(C[i].Y()-PosY[count])*(C[i].Y()-PosY[count]);
}
+
devX[count]=sqrt(devX[count]/((size-1)*norm));
devY[count]=sqrt(devY[count]/((size-1)*norm));
+
// Compute ThetaX, angle between the Direction vector projection in XZ with
// the Z axis
- //ThetaX=atan((PosX100-PosX)/100.);
- ThetaX[count] = (PosX100-PosX[count])/100.;
+ ThetaX[count]=(PosX100-PosX[count])/100.;
// Compute PhiY, angle between the Direction vector projection in YZ with
// the Z axis
- //PhiY=atan((PosY100-PosY)/100.);
PhiY[count]=(PosY100-PosY[count])/100.;
Dir[count]=TVector3(PosX100-PosX[count],PosY100-PosY[count],100).Unit();
}
}
+
if(PosX[count]==0){
PosX[count]=-10000;
PosY[count]=-10000;
@@ -262,7 +263,6 @@ void TSamuraiBDCPhysics::PreTreat(){
m_DCHit[det].push_back(DCHit(det,layer,wire,time,etime-time,2.5-Cal->ApplySigmoid(channel,etime)));
}
}
-
}
return;
}
@@ -278,6 +278,7 @@ void TSamuraiBDCPhysics::Clear(){
devY.clear();
Dir.clear();
PileUp.clear();
+ Detector.clear();
}
///////////////////////////////////////////////////////////////////////////
diff --git a/NPLib/TrackReconstruction/NPDCReconstruction.cxx b/NPLib/TrackReconstruction/NPDCReconstruction.cxx
index c1978c899..fd24aee03 100644
--- a/NPLib/TrackReconstruction/NPDCReconstruction.cxx
+++ b/NPLib/TrackReconstruction/NPDCReconstruction.cxx
@@ -31,7 +31,7 @@ DCReconstruction::DCReconstruction(){
m_min->SetPrintLevel(0);
// this avoid error
- gErrorIgnoreLevel = kError;
+ gErrorIgnoreLevel = kError;
//m_min->SetMaxFunctionCalls(1000);
//m_min->SetMaxIterations(1000);
//m_min->SetTolerance(1);
@@ -52,10 +52,20 @@ double DCReconstruction::BuildTrack2D(const vector<double>& X,const vector<doubl
// 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 = (Z[size-1]-Z[0])/(X[size-1]-R[size-1]-X[0]-R[0]);
- // bi = Z[0]-ai*(X[0]+R[0]);
- ai = (Z[sizeX-1]-Z[0])/(X[sizeX-1]-X[0]);
- bi = Z[0]-ai*(X[0]);
+ 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]);
+ i++;
+ }
+ if(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);
@@ -144,17 +154,17 @@ TGraph* DCReconstruction::Scan(double a, double b, int tovary, double minV, doub
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));
- }
+ 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]);
diff --git a/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx b/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx
index 3f4d5ae8f..8e6b3c62c 100644
--- a/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx
+++ b/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx
@@ -236,24 +236,37 @@ void DCReconstructionMT::StartThread(unsigned int 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();
+ unsigned int i = 1;
+ ai=1/0.;
+ while(isinf(ai)&&i!=sizeX[id]){
+ ai = ((*fitZ[id])[sizeX[id]-i]-(*fitZ[id])[0])/((*fitX[id])[sizeX[id]-i]-(*fitX[id])[0]);
+ bi = (*fitZ[id])[0]-ai*((*fitX[id])[0]);
+ i++;
+ }
+ 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;
+ m_minimum[uid] = 10000;
+ }
+ else{
+ 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();
+ // 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_mtx.lock();// make sure no other thread is reading/writing to the map
m_Ready[id].flip();
--
GitLab