* 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

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();
 }
 ///////////////////////////////////////////////////////////////////////////
 

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]);

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();