diff --git a/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx b/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx
index 75ae729ff4f7032199e117626b571cb056109888..9017c919d68a4ae77948733498a935f7f4f2f195 100644
--- a/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx
+++ b/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.cxx
@@ -154,7 +154,7 @@ void TSamuraiBDCPhysics::BuildPhysicalEvent(){
             //cout << "BDC" << endl;
             m_reconstruction.ResolvePlane(it1->second,it1->first,it2->second,it2->first,P);
             // cout << "done " << D[it1->first] << " " << D[it2->first] << endl;
-            if(P.X()!=-10000 && D[it1->first]<PowerThreshold&& D[it2->first]<PowerThreshold){
+            if(P.X()!=-10000 /*&& D[it1->first]<PowerThreshold&& D[it2->first]<PowerThreshold*/){
               C.push_back(P);
               // Mean pos are weighted based on the the sum of distance from track
               // to hit obtained during the minimisation
@@ -172,7 +172,7 @@ void TSamuraiBDCPhysics::BuildPhysicalEvent(){
           if(it1!=it2 ){// different plane, same detector
             m_reconstruction.ResolvePlane(it1->second,it1->first,it2->second,it2->first,P);
 
-            if(P.X()!=-10000&& D[it1->first]<PowerThreshold && D[it2->first]<PowerThreshold)
+            if(P.X()!=-10000/*&& D[it1->first]<PowerThreshold && D[it2->first]<PowerThreshold*/)
               C100.push_back(P);
           }
         }
@@ -272,7 +272,10 @@ void TSamuraiBDCPhysics::PreTreat(){
       if(etime && time && etime-time>ToTThreshold_L && etime-time<ToTThreshold_H){
         channel="SamuraiBDC"+NPL::itoa(det)+"/L" + NPL::itoa(layer);
         SamuraiDCIndex idx(det,layer,wire);
-        m_DCHit[det].push_back(DCHit(det,layer,wire,time,etime-time,2.5-Cal->ApplySigmoid(channel,etime)));
+        if(!m_invertD[det])
+          m_DCHit[det].push_back(DCHit(det,layer,wire,time,etime-time,Cal->ApplySigmoid(channel,etime)));
+        else
+          m_DCHit[det].push_back(DCHit(det,layer,wire,time,etime-time,2.5-Cal->ApplySigmoid(channel,etime)));
       }
     }
   }
@@ -303,7 +306,7 @@ void TSamuraiBDCPhysics::ReadConfiguration(NPL::InputParser parser){
   if(NPOptionManager::getInstance()->GetVerboseLevel())
     cout << "//// " << blocks.size() << " detector(s) found " << endl; 
 
-  vector<string> token= {"XML","Offset","InvertX","InvertY"};
+  vector<string> token= {"XML","Offset","InvertX","InvertY","InvertD"};
 
   for(unsigned int i = 0 ; i < blocks.size() ; i++){
     if(blocks[i]->HasTokenList(token)){
@@ -316,9 +319,11 @@ void TSamuraiBDCPhysics::ReadConfiguration(NPL::InputParser parser){
       TVector3 offset = blocks[i]->GetTVector3("Offset","mm"); 
       bool invertX = blocks[i]->GetInt("InvertX"); 
       bool invertY = blocks[i]->GetInt("InvertY"); 
+      bool invertD = blocks[i]->GetInt("InvertD"); 
       m_offset[det] = offset;
       m_invertX[det] = invertX;
       m_invertY[det] = invertY;
+      m_invertD[det] = invertD;
     }
   }
 
diff --git a/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.h b/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.h
index bfa9c5b9d0135ce3e7f024acce646239e87d49e9..37219d89cf01762b1bbeacb8831c184f9b7e1a1e 100644
--- a/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.h
+++ b/NPLib/Detectors/Samurai/TSamuraiBDCPhysics.h
@@ -199,6 +199,7 @@ class TSamuraiBDCPhysics : public TObject, public NPL::VDetector{
     std::map<unsigned int, TVector3> m_offset;//!
     std::map<unsigned int, bool> m_invertX;//!
     std::map<unsigned int, bool> m_invertY;//!
+    std::map<unsigned int, bool> m_invertD;//!
 
   private: // Spectra Class
     // TSamuraiBDCSpectra* m_Spectra; // !
diff --git a/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.cxx b/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.cxx
index cb344120d7e46db5ff1d3151f77af8fdd2c9964a..b1538925e436aa0884dbf3863d915501a989166b 100644
--- a/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.cxx
+++ b/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.cxx
@@ -46,11 +46,14 @@ ClassImp(TSamuraiFDC0Physics)
     m_PreTreatedData    = new TSamuraiFDC0Data ;
     m_EventPhysics      = this ;
     //m_Spectra           = NULL;
+//    ToTThreshold_L = 0;
+//    ToTThreshold_H = 180;
     ToTThreshold_L = 0;
-    ToTThreshold_H = 180;
+    ToTThreshold_H = 10000;
+ 
     DriftLowThreshold=0.1 ;
     DriftUpThreshold=2.4;
-    PowerThreshold=5;
+    PowerThreshold=15;
   }
 
 ///////////////////////////////////////////////////////////////////////////
@@ -107,9 +110,9 @@ void TSamuraiFDC0Physics::BuildPhysicalEvent(){
 #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);
+   { std::ofstream f("distance.txt", std::ios::app);
    f<< D[it->first] << endl;
    f.close();
    }
@@ -141,7 +144,7 @@ void TSamuraiFDC0Physics::BuildPhysicalEvent(){
         //cout << "FDC0" << endl;
         m_reconstruction.ResolvePlane(it1->second,it1->first,it2->second,it2->first,P);
         // cout << "done " << D[it1->first] << " " << D[it2->first] << endl;
-        if(P.X()!=-10000 && D[it1->first]<PowerThreshold&& D[it2->first]<PowerThreshold){
+        if(P.X()!=-10000/* && D[it1->first]<PowerThreshold&& D[it2->first]<PowerThreshold*/){
           C.push_back(P);
           // Mean pos are weighted based on the the sum of distance from track
           // to hit obtained during the minimisation
@@ -158,7 +161,7 @@ void TSamuraiFDC0Physics::BuildPhysicalEvent(){
       if(it1!=it2 ){// different plane, same detector
         m_reconstruction.ResolvePlane(it1->second,it1->first,it2->second,it2->first,P);
 
-        if(P.X()!=-10000&& D[it1->first]<PowerThreshold && D[it2->first]<PowerThreshold)
+        if(P.X()!=-10000/*&& D[it1->first]<PowerThreshold && D[it2->first]<PowerThreshold*/)
           C100.push_back(P);
       }
     }
@@ -180,10 +183,12 @@ void TSamuraiFDC0Physics::BuildPhysicalEvent(){
       PosY100+= C100[i].Y()*W[i]; 
       norm+=W[i];
     } 
+
     MultMean=size;
     // Mean position at Z=0
     PosX=PosX/norm; 
     PosY=PosY/norm; 
+
     // Mean position at Z=100
     PosX100=PosX100/norm; 
     PosY100=PosY100/norm; 
@@ -196,15 +201,28 @@ void TSamuraiFDC0Physics::BuildPhysicalEvent(){
     }
     devX=sqrt(devX/((size-1)*norm));
     devY=sqrt(devY/((size-1)*norm));
+   
+    if(m_invertX){
+      PosX*=-1;
+      PosX100*=-1;
+    }
+
+    if(m_invertY){
+      PosY*=-1;
+      PosY100*=-1;
+    }
     // Compute ThetaX, angle between the Direction vector projection in XZ with
     // the Z axis
-    //ThetaX=atan((PosX100-PosX)/100.);
-    ThetaX = (PosX100-PosX)/100.;
+    ThetaX=atan((PosX100-PosX)/100.);
+    //ThetaX = (PosX100-PosX)/100.;
     // Compute PhiY, angle between the Direction vector projection in YZ with
     // the Z axis
-    //PhiY=atan((PosY100-PosY)/100.);
-    PhiY=(PosY100-PosY)/100.;
+    PhiY=atan((PosY100-PosY)/100.);
+    //PhiY=(PosY100-PosY)/100.;
     Dir=TVector3(PosX100-PosX,PosY100-PosY,100).Unit();
+    PosX+=m_offset.X();
+    PosY+=m_offset.Y();
+
   }
 
   return;
@@ -253,7 +271,10 @@ void TSamuraiFDC0Physics::PreTreat(){
         Time.push_back(time);
         ToT.push_back(etime-time);
         channel="SamuraiFDC0/L" + NPL::itoa(layer);
-        DriftLength.push_back(2.5-Cal->ApplySigmoid(channel,etime));
+        if(!m_invertD)
+          DriftLength.push_back(Cal->ApplySigmoid(channel,etime));
+        else
+          DriftLength.push_back(2.5-Cal->ApplySigmoid(channel,etime));
       }
     }
 
@@ -306,6 +327,7 @@ void TSamuraiFDC0Physics::Clear(){
   PileUp=0;
   Mult=0;
   PosX=PosY=-10000;
+  ThetaX=PhiY=-10000;
   devX=devY=-10000;
   DriftLength.clear();
   Detector.clear();
@@ -327,7 +349,7 @@ void TSamuraiFDC0Physics::ReadConfiguration(NPL::InputParser parser){
   if(NPOptionManager::getInstance()->GetVerboseLevel())
     cout << "//// " << blocks.size() << " detector(s) found " << endl; 
 
-  vector<string> token= {"XML"};
+  vector<string> token= {"XML","Offset","InvertX","InvertY","InvertD"};
 
   for(unsigned int i = 0 ; i < blocks.size() ; i++){
     cout << endl << "////  Samurai FDC0 (" << i+1 << ")" << endl;
@@ -335,6 +357,10 @@ void TSamuraiFDC0Physics::ReadConfiguration(NPL::InputParser parser){
     NPL::XmlParser xml;
     xml.LoadFile(xmlpath);
     AddDC("SAMURAIFDC0",xml);
+    m_offset = blocks[i]->GetTVector3("Offset","mm"); 
+    m_invertX = blocks[i]->GetInt("InvertX"); 
+    m_invertY = blocks[i]->GetInt("InvertY"); 
+    m_invertD = blocks[i]->GetInt("InvertD"); 
   }
 
 #if __cplusplus > 199711L && NPMULTITHREADING 
diff --git a/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.h b/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.h
index 161c587ce0dd1cd94cb0fcff536262eac8ad779d..ae43f580d56268cc9f617f23c2c923ecbfb27913 100644
--- a/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.h
+++ b/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.h
@@ -181,6 +181,11 @@ class TSamuraiFDC0Physics : public TObject, public NPL::VDetector{
     TSamuraiFDC0Data*         m_PreTreatedData;//!
     TSamuraiFDC0Physics*      m_EventPhysics;//!
 
+  private: // offset and inversion 
+    TVector3 m_offset;//!
+    bool m_invertX;//!
+    bool m_invertY;//!
+    bool m_invertD;//!
 
   private: // Spectra Class
    // TSamuraiFDC0Spectra* m_Spectra; // !
diff --git a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx
index 6385567618500aabb2204c0b1dcabe94a4514ab3..7f2b5a7fdc99ad3c85a2d9b684461f5010c01ca5 100644
--- a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx
+++ b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx
@@ -44,11 +44,15 @@ ClassImp(TSamuraiFDC2Physics)
     m_EventData         = new TSamuraiFDC2Data ;
     m_EventPhysics      = this ;
     //m_Spectra           = NULL;
-    ToTThreshold_L = 180;
-    ToTThreshold_H = 1000;
-    DriftLowThreshold=0.4;
-    DriftUpThreshold=9.3;
-    PowerThreshold=14;
+ //   ToTThreshold_L = 180;
+ //   ToTThreshold_H = 1000;
+    ToTThreshold_L = 0;
+    ToTThreshold_H = 10000; 
+    DriftLowThreshold=0.1;
+    DriftUpThreshold=9.9;
+    //PowerThreshold=15;
+
+    PowerThreshold=1e6;
   }
 
 ///////////////////////////////////////////////////////////////////////////
@@ -110,14 +114,14 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
   D[it->first]=m_reconstruction.GetResults(uid++,X0,X100,a,b); 
 #endif
 
-  /* // for Debug, write a file of 
+/*   // for Debug, write a file of 
    { std::ofstream f("distance.txt", std::ios::app);
    f<< D[it->first] << endl;
    f.close();
    }
-   */
-    // very large a means track perpendicular to the chamber, what happen when there is pile up
-    if(abs(a)>5000)
+*/    
+   // very large a means track perpendicular to the chamber, what happen when there is pile up
+   if(abs(a)>5000)
       PileUp++;
 
     Mult+=X[it->first].size();
@@ -141,8 +145,7 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
     for(auto it2 = it1;it2!=VX0.end();++it2){
       if(it1!=it2){// different plane, same detector
         m_reconstruction.ResolvePlane(it1->second,it1->first,it2->second,it2->first,P);
-        
-        if(P.X()!=-10000 && D[it1->first]<PowerThreshold && D[it2->first]<PowerThreshold){
+        if(P.X()!=-10000 /*&& D[it1->first]<PowerThreshold && D[it2->first]<PowerThreshold*/){
           C.push_back(P);
           // Mean pos are weighted based on the the sum of distance from track
           // to hit obtained during the minimisation
@@ -151,7 +154,7 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
       }
     }
   }
-
+  
   // Reconstruct the position at z=100 for each detector
   static vector<TVector3> C100 ;  
   C100.clear();
@@ -160,7 +163,7 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
       if(it1!=it2){// different plane
         m_reconstruction.ResolvePlane(it1->second,it1->first,it2->second,it2->first,P);
 
-        if(P.X()!=-10000&& D[it1->first]<PowerThreshold && D[it2->first]<PowerThreshold)
+        if(P.X()!=-10000/*&& D[it1->first]<PowerThreshold && D[it2->first]<PowerThreshold*/)
           C100.push_back(P);
       }
     }
@@ -202,6 +205,16 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
 
     devX=sqrt(devX/((size-1)*norm));
     devY=sqrt(devY/((size-1)*norm));
+   
+    if(m_invertX){
+      PosX*=-1;
+      PosX100*=-1;
+    }
+
+    if(m_invertY){
+      PosY*=-1;
+      PosY100*=-1;
+    }
 
     // Compute ThetaX, angle between the Direction vector projection in XZ with
     // the Z axis
@@ -212,8 +225,13 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
     //PhiY=atan((PosY100-PosY)/100.);
     PhiY=(PosY100-PosY)/100.;
     Dir=TVector3(PosX100-PosX,PosY100-PosY,100).Unit();
+    PosX+=m_offset.X();
+    PosY+=m_offset.Y();
   }
-//cout <<size << " " << PosX << endl;
+/*  if(PosX==-10000)
+    cout << " bad " <<  Detector.size()<< " " << size << endl;
+  else
+    cout << " okay" <<  Detector.size()<< " " << size << endl;*/
   return;
 }
 ///////////////////////////////////////////////////////////////////////////
@@ -258,7 +276,10 @@ void TSamuraiFDC2Physics::PreTreat(){
           channel="SamuraiFDC2/L" + NPL::itoa(layer);
           // rescalling is needed because calib are bad.
           // to be fixed
-          DriftLength.push_back(10-Cal->ApplySigmoid(channel,etime));
+          if(!m_invertD)
+            DriftLength.push_back(Cal->ApplySigmoid(channel,etime));
+          else
+            DriftLength.push_back(10-Cal->ApplySigmoid(channel,etime));
         }
       }
     }
@@ -329,6 +350,7 @@ void TSamuraiFDC2Physics::Clear(){
   PileUp=0;
   Mult=0;
   PosX=PosY=-10000;
+  ThetaX=PhiY=-10000;
   devX=devY=-10000;
   DriftLength.clear();
   Detector.clear();
@@ -350,7 +372,7 @@ void TSamuraiFDC2Physics::ReadConfiguration(NPL::InputParser parser){
   if(NPOptionManager::getInstance()->GetVerboseLevel())
     cout << "//// " << blocks.size() << " detector(s) found " << endl; 
 
-  vector<string> token= {"XML"};
+  vector<string> token= {"XML","Offset","InvertX","InvertY","InvertD"};
 
   for(unsigned int i = 0 ; i < blocks.size() ; i++){
     cout << endl << "////  Samurai FDC2 (" << i+1 << ")" << endl;
@@ -358,6 +380,10 @@ void TSamuraiFDC2Physics::ReadConfiguration(NPL::InputParser parser){
     NPL::XmlParser xml;
     xml.LoadFile(xmlpath);
     AddDC("SAMURAIFDC2",xml);
+    m_offset = blocks[i]->GetTVector3("Offset","mm"); 
+    m_invertX = blocks[i]->GetInt("InvertX"); 
+    m_invertY = blocks[i]->GetInt("InvertY"); 
+    m_invertD = blocks[i]->GetInt("InvertD"); 
   }
 
 
diff --git a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h
index b601cd71b050f3329290ba98caf0bab6497985da..99c03a43b478000246707d5702446cfa61cadbc0 100644
--- a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h
+++ b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h
@@ -176,6 +176,11 @@ class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{
     TSamuraiFDC2Data*         m_EventData;//!
     TSamuraiFDC2Physics*      m_EventPhysics;//!
 
+  private: // offset and inversion 
+    TVector3 m_offset;//!
+    bool m_invertX;//!
+    bool m_invertY;//!
+    bool m_invertD;//!
 
   private: // Spectra Class
    // TSamuraiFDC2Spectra* m_Spectra; // !
diff --git a/NPLib/Physics/NPParticle.h b/NPLib/Physics/NPParticle.h
index 9a470fb1aa8d115810483fd524a08cda24ceced2..fba2e1deef973628296c0780cb1c01d75a85ac08 100644
--- a/NPLib/Physics/NPParticle.h
+++ b/NPLib/Physics/NPParticle.h
@@ -126,8 +126,8 @@ namespace NPL {
     void				SetSpinParity(const char* spinparity)	{fSpinParity = spinparity;}
     void				SetSpin(double spin) {fSpin = spin;}
     void				SetParity(const char* parity)	{fParity = parity;}
-    void          SetLifeTime(double LifeTime) {fLifeTime=LifeTime;}
-    void          SetLifeTimeError(double LifeTimeErr) {fLifeTimeErr=LifeTimeErr;}
+    void        SetLifeTime(double LifeTime) {fLifeTime=LifeTime;}
+    void        SetLifeTimeError(double LifeTimeErr) {fLifeTimeErr=LifeTimeErr;}
     void				SetKineticEnergy(double energy)	{fKineticEnergy = energy; EnergyToBrho(); EnergyToTof(); EnergyToBeta(); BetaToGamma();BetaToVelocity();}
     void				SetBrho(double brho) {fBrho = brho; BrhoToEnergy(); BrhoToTof(); EnergyToBeta(); BetaToGamma();BetaToVelocity();}
     void				SetTimeOfFlight(double tof) {fTimeOfFlight = tof; TofToEnergy(); TofToBrho(); EnergyToBeta(); BetaToGamma();BetaToVelocity();}
diff --git a/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx b/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx
index 8e6b3c62ce4f11ed07c8c9e3fb95de563391c698..73c443aeca852ef06f3b37fc1e67a563fe7a9eb4 100644
--- a/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx
+++ b/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx
@@ -236,13 +236,9 @@ 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
-      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++;
-      }
+        ai = ((*fitZ[id])[sizeX[id]-1]-(*fitZ[id])[0])/((*fitX[id])[sizeX[id]-1]-(*fitX[id])[0]+(*fitR[id])[sizeX[id]-1]-(*fitR[id])[0]);
+        bi = (*fitZ[id])[0]-ai*((*fitX[id])[0]+(*fitR[id])[0]);
+
       if(isinf(ai)){ // then there is no two point in different layer
         m_a[uid]=-10000;
         m_b[uid]=-10000;
@@ -250,6 +246,7 @@ void DCReconstructionMT::StartThread(unsigned int id){
         m_X100[uid]=-10000;
         m_minimum[uid] = 10000;
       }
+
       else{
         mini->Clear(); 
         mini->SetVariable(0,"a",ai,1);
diff --git a/Projects/S034/Analysis.cxx b/Projects/S034/Analysis.cxx
index 7c82b7244b88e9f2ce8d90a3b20a8cd0fc4dc03c..9094323d62dc7f379e8f5af977d281d87fcb00f5 100644
--- a/Projects/S034/Analysis.cxx
+++ b/Projects/S034/Analysis.cxx
@@ -54,13 +54,19 @@ void Analysis::TreatEvent(){
   //cout << Trigger << " " ; 
   Trigger=Trigger&0x00ff;
 //cout << Trigger << endl;
-  // Compute Broh 
-  if(FDC2->PosX>-10000 && FDC0->PosX>-10000){ // if both are correctly build
-  double broh_param[6]={FDC0->PosX, FDC0->PosY, tan(FDC0->ThetaX), tan(FDC0->PhiY), FDC2->PosX, FDC2->ThetaX};
+  // Compute Brho 
+  if(FDC2->PosX>-10000 && FDC0->PosX>-10000 ){ // if both are correctly build
+   // Compute ThetaX and PhiY using Minos vertex and FDC0 XY
+   double FDC0_ThetaX = FDC0->ThetaX;
+   double FDC0_PhiY   = FDC0->PhiY;
 
-  Broh=r_fit(broh_param);
+   if(Minos->Z_Vertex>0){
+    FDC0_ThetaX = atan((FDC0->PosX-Minos->X_Vertex)/(1283.7-Minos->Z_Vertex));
+    FDC0_PhiY   = atan((FDC0->PosY-Minos->Y_Vertex)/(1283.7-Minos->Z_Vertex));
+   } 
+    double brho_param[6]={FDC0->PosX/*+1.77*/, FDC0->PosY, tan(FDC0_ThetaX), tan(FDC0_PhiY), FDC2->PosX/*-252.55*/, FDC2->ThetaX};
+    Brho=r_fit(brho_param);
   }
-
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -69,13 +75,13 @@ void Analysis::End(){
 
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::Clear(){
-  Broh=-1000;
+  Brho=-1000;
   Beta_f=-1000;
 }
 
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::InitOutputBranch() {
-  RootOutput::getInstance()->GetTree()->Branch("Broh",&Broh,"Broh/D");
+  RootOutput::getInstance()->GetTree()->Branch("Brho",&Brho,"Brho/D");
   RootOutput::getInstance()->GetTree()->Branch("Beta_f",&Beta_f,"Beta_f/D");
   RootOutput::getInstance()->GetTree()->Branch("Trigger",&Trigger,"Trigger/I");
 } 
diff --git a/Projects/S034/Analysis.h b/Projects/S034/Analysis.h
index e1e1136b0229229870079d613d69e271de76e536..6c893a7b90ce7a5c4b6ee1eb726e5077e4dddfee 100644
--- a/Projects/S034/Analysis.h
+++ b/Projects/S034/Analysis.h
@@ -50,7 +50,7 @@ class Analysis: public NPL::VAnalysis{
     TSamuraiHodoscopePhysics* Hodo;
  
   private: // output variable
-    double Broh;
+    double Brho;
     double Beta_f;
     int    Trigger;
   public:
diff --git a/Projects/S034/macro/IDz.cxx b/Projects/S034/macro/IDz.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..124bc11b3a807b32e4e773035c2999a23969c74e
--- /dev/null
+++ b/Projects/S034/macro/IDz.cxx
@@ -0,0 +1,35 @@
+void IDz(){
+  auto fz = new TFile("root/zaihong/run0582_RIG20210424_6He.root");
+  auto tz = (TTree*) fz->FindObjectAny("rig");
+ 
+  tz->Draw("rig/1000.>>h(1000,0,8)","","");  
+  tz->Draw("rig/1000.>>hn(1000,0,8)","FragID==-1","same");  
+  tz->Draw("rig/1000.>>hd(1000,0,8)","FragID==12","same");  
+  tz->Draw("rig/1000.>>ht(1000,0,8)","FragID==13","same");  
+  tz->Draw("rig/1000.>>h4He(1000,0,8)","FragID==24","same");  
+  tz->Draw("rig/1000.>>h6He(1000,0,8)","FragID==26","same");  
+  auto hn = (TH1*) gDirectory->FindObjectAny("hn");
+  auto hd = (TH1*) gDirectory->FindObjectAny("hd");
+  auto ht = (TH1*) gDirectory->FindObjectAny("ht");
+  auto h4He = (TH1*) gDirectory->FindObjectAny("h4He");
+  auto h6He = (TH1*) gDirectory->FindObjectAny("h6He");
+
+  hn->SetFillStyle(1001);
+  hd ->SetFillStyle(1001);
+  ht  ->SetFillStyle(1001);
+  h4He ->SetFillStyle(1001);
+  h6He ->SetFillStyle(1001);
+
+  hn->SetFillColor(kBlack);
+  hd ->SetFillColor(kAzure+7);
+  ht  ->SetFillColor(kOrange+7);
+  h4He ->SetFillColor(kGreen-3);
+  h6He ->SetFillColor(kMagenta-3);
+  hn->SetLineColor(kBlack);
+  hd ->SetLineColor(kAzure+7);
+  ht  ->SetLineColor(kOrange+7);
+  h4He ->SetLineColor(kGreen-3);
+  h6He ->SetLineColor(kMagenta-3);
+  gPad->Update();
+
+}
diff --git a/Projects/S034/macro/compare.cxx b/Projects/S034/macro/compare.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..07ee2a00dbaa72f28a7390f3ffca8c0f74eb124b
--- /dev/null
+++ b/Projects/S034/macro/compare.cxx
@@ -0,0 +1,101 @@
+void compare(){
+
+ auto fz = new TFile("root/zaihong/run0582_RIG20210424_6He.root");
+ auto tz = (TTree*) fz->FindObjectAny("rig");
+ auto fl = new TFile("root/analysis/test582.root");
+ auto tl = (TTree*) fl->FindObjectAny("PhysicsTree");
+ tl->AddFriend(tz);
+
+ auto cfdc0= new TCanvas();
+ cfdc0->Divide(2,4);
+ cfdc0->cd(1);
+ tz->Draw("FDC0_X>>h1(1000,-80,80)","FDC0_X>-80 && FDC0_X<80");
+ cfdc0->cd(2);
+ tz->Draw("FDC0_Y>>h2(1000,-80,80)","FDC0_X>-80 && FDC0_Y<80");
+ cfdc0->cd(3);
+ tl->Draw("SamuraiFDC0.PosX>>h3(1000,-80,80)","SamuraiFDC0.PosX>-80 && SamuraiFDC0.PosX<80");
+ cfdc0->cd(4);
+ tl->Draw("SamuraiFDC0.PosY>>h4(1000,-80,80)","SamuraiFDC0.PosY>-80 && SamuraiFDC0.PosY<80");
+ cfdc0->cd(5);
+ tl->Draw("SamuraiFDC0.PosX:FDC0_X>>h5(1000,-80,80,1000,-80,80)","SamuraiFDC0.PosX>-80 && SamuraiFDC0.PosX<80","colz");
+ auto l = new TLine(-80,-80,80,80);
+ l->SetLineColor(kOrange+7);
+ l->SetLineWidth(2);
+ l->Draw();
+ cfdc0->cd(6);
+ tl->Draw("SamuraiFDC0.PosY:FDC0_Y>>h6(1000,-80,80,1000,-80,80)","SamuraiFDC0.PosY>-80 && SamuraiFDC0.PosY<80","colz");
+ l->Draw();
+
+ cfdc0->cd(7);
+ tl->Draw("SamuraiFDC0.PosX-FDC0_X>>h7(1000,-5,5)","SamuraiFDC0.PosX>-80 && SamuraiFDC0.PosX<80","");
+ cfdc0->cd(8);
+ tl->Draw("SamuraiFDC0.PosY-FDC0_Y>>h8(1000,-5,5)","SamuraiFDC0.PosY>-80 && SamuraiFDC0.PosY<80","");
+ 
+ auto cfdc2= new TCanvas();
+ cfdc2->Divide(2,4);
+ cfdc2->cd(1);
+ tz->Draw("FDC2_X>>g1(1000,-2000,2000)","FDC2_X>-2000 && FDC2_X<2000");
+ cfdc2->cd(2);
+ tz->Draw("FDC2_Y>>g2(1000,-2000,2000)","FDC2_X>-2000 && FDC2_Y<2000");
+ cfdc2->cd(3);
+ tl->Draw("SamuraiFDC2.PosX>>g3(1000,-2000,2000)","SamuraiFDC2.PosX>-2000 && SamuraiFDC2.PosX<2000&&SamuraiFDC2.devX<10");
+ cfdc2->cd(4);
+ tl->Draw("SamuraiFDC2.PosY>>g4(1000,-2000,2000)","SamuraiFDC2.PosY>-2000 && SamuraiFDC2.PosY<2000&&SamuraiFDC2.devY<10");
+ cfdc2->cd(5);
+ tl->Draw("SamuraiFDC2.PosX:FDC2_X>>g5(1000,-2000,2000,1000,-2000,2000)","SamuraiFDC2.PosX>-2000 && SamuraiFDC2.PosX<2000&&SamuraiFDC2.devX<10","colz");
+ l = new TLine(-2000,-2000,2000,2000);
+ l->SetLineColor(kOrange+7);
+ l->SetLineWidth(2);
+ l->Draw();
+
+ cfdc2->cd(6);
+ tl->Draw("SamuraiFDC2.PosY:FDC2_Y>>g6(1000,-2000,2000,1000,-2000,2000)","SamuraiFDC2.PosY>-2000 && SamuraiFDC2.PosY<2000&&SamuraiFDC2.devX<10","colz");
+ l->Draw();
+
+ cfdc2->cd(7);
+ tl->Draw("SamuraiFDC2.PosX-FDC2_X>>g7(1000,-5,5)","SamuraiFDC2.PosX>-2000 && SamuraiFDC2.PosX<2000&&SamuraiFDC2.devX<10","");
+ cfdc2->cd(8);
+ tl->Draw("SamuraiFDC2.PosY-FDC2_Y>>g8(1000,-5,5)","SamuraiFDC2.PosY>-2000 && SamuraiFDC2.PosY<2000&&SamuraiFDC2.devY<10","");
+ 
+ auto cfdc2b= new TCanvas();
+ cfdc2b->Divide(2,4);
+ cfdc2b->cd(1);
+ tz->Draw("FDC2_ThetaX>>i1(1000,-2,2)","FDC2_ThetaX>-1000");
+ cfdc2b->cd(2);
+ tz->Draw("FDC2_ThetaY>>i2(1000,-2,2)","FDC2_ThetaY>-1000");
+ cfdc2b->cd(3);
+ tl->Draw("SamuraiFDC2.ThetaX>>i3(1000,-2,2)","SamuraiFDC2.ThetaX>-10000");
+ cfdc2b->cd(4);
+ tl->Draw("SamuraiFDC2.PhiY>>i4(1000,-2,2)","SamuraiFDC2.PhiY>-1000");
+ cfdc2b->cd(5);
+ tl->Draw("SamuraiFDC2.ThetaX:FDC2_ThetaX>>i5(1000,-2,2,1000,-2,2)","FDC2_ThetaX>-1000&&SamuraiFDC2.ThetaX>-10000","colz");
+ l = new TLine(-2,-2,2,2);
+ l->SetLineColor(kOrange+7);
+ l->SetLineWidth(2);
+ l->Draw();
+
+ cfdc2b->cd(6);
+ tl->Draw("SamuraiFDC2.PhiY:FDC2_ThetaY>>i6(1000,-2,2,1000,-2,2)","FDC2_ThetaY>-1000&&SamuraiFDC2.PhiY>-1000","colz");
+ l->Draw();
+
+ cfdc2b->cd(7);
+ tl->Draw("SamuraiFDC2.ThetaX-FDC2_ThetaX>>i7(1000,-0.1,0.1)","FDC2_ThetaX>-1000&&SamuraiFDC2.ThetaX>-10000","");
+ cfdc2b->cd(8);
+ tl->Draw("SamuraiFDC2.PhiY-FDC2_ThetaY>>i8(1000,-0.1,0.1)","FDC2_ThetaY>-1000&&SamuraiFDC2.PhiY>-1000","");
+
+ auto cbroh= new TCanvas();
+ cbroh->Divide(2,2);
+ cbroh->cd(1);
+ tl->Draw("Brho>>b1(1000,0,8)","Brho>0 && SamuraiFDC2.devX<10 && SamuraiFDC2.devX<10");
+ cbroh->cd(2);
+ tz->Draw("rig>>b2(1000,0,8000)","rig>0");
+ cbroh->cd(3);
+ tl->Draw("Brho:rig/1000.>>b3(1000,0,8,1000,0,8)","rig>0 && Brho>0","colz");
+ l = new TLine(0,0,8,8);
+ l->SetLineColor(kOrange+7);
+ l->SetLineWidth(2);
+ l->Draw();
+ cbroh->cd(4);
+ tz->Draw("rig/1000.:FragID>>b4(1000,-2000,2000,31,-1,30)","","colz");
+
+}
diff --git a/Projects/S034/macro/dist.cxx b/Projects/S034/macro/dist.cxx
index dadc60f8eb1a0aa2ae8748a6b0d128db588f4dc0..fd8eff3a0d8ad8c7079ab60fe4ba95c6105ea49d 100644
--- a/Projects/S034/macro/dist.cxx
+++ b/Projects/S034/macro/dist.cxx
@@ -2,7 +2,7 @@ void dist(){
 
   ifstream f("distance.txt");
   double d;
-  auto h = new TH1D("h","h",1000,0,10);
+  auto h = new TH1D("h","h",1000,0,100);
 
   while(f>>d){
     h->Fill(d);
diff --git a/Projects/S034/macro/rigz.cxx b/Projects/S034/macro/rigz.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..029b37dc0b9a28ec2c1b8ae42aa055bd347b588a
--- /dev/null
+++ b/Projects/S034/macro/rigz.cxx
@@ -0,0 +1,1048 @@
+double r_fit(double *x) ;
+
+void rigz(){
+
+ auto fz = new TFile("root/zaihong/run0582_RIG20210424_6He.root");
+ auto tz = (TTree*) fz->FindObjectAny("rig");
+ auto fl = new TFile("root/analysis/test582.root");
+ auto tl = (TTree*) fl->FindObjectAny("PhysicsTree");
+ 
+ double FDC0_X,FDC0_Y,FDC2_X,FDC2_ThetaX;
+ tz->SetBranchAddress("FDC0_X",&FDC0_X);
+ tz->SetBranchAddress("FDC0_Y",&FDC0_Y);
+ tz->SetBranchAddress("FDC2_X",&FDC2_X);
+ tz->SetBranchAddress("FDC2_ThetaX",&FDC2_ThetaX);
+ auto h = new TH1D("brho","brho",1000,0,8);
+ unsigned int entries = tz->GetEntries();
+
+ for(unsigned int i = 0 ; i < 100000 /*entries*/ ; i++){
+  tz->GetEntry(i);
+  double brho_param[6]={FDC0_X/*+1.77*/, FDC0_Y, 0, 0, FDC2_X/*-252.55*/, FDC2_ThetaX};
+  double Brho=r_fit(brho_param);
+  if(Brho>0)
+    h->Fill(Brho);
+  }
+  h->Scale(1./h->Integral());
+  h->Draw();
+  cout << tl->Draw("Brho>>g","Brho>0","same") << endl;
+  auto g = (TH1*) gDirectory->FindObjectAny("g");
+  g->Scale(1./g->Integral());
+}
+// EOF for ../Geant4/samurai/simtree_energy.C
+// -*- mode: c++ -*-
+// 
+// File ../Geant4/samurai/simtree_momentum.C generated by TMultiDimFit::MakeRealCode
+// on Tue Jun  9 23:41:17 2020
+// ROOT version 5.32/00
+//
+// This file contains the function 
+//
+//    double  p_fit(double *x); 
+//
+// For evaluating the parameterization obtained
+// from TMultiDimFit and the point x
+// 
+// See TMultiDimFit class documentation for more information 
+// 
+//
+// Static data variables
+//
+static int    p_gNVariables    = 6;
+static int    p_gNCoefficients = 57;
+static double p_gDMean         = 524.342;
+// Assignment to mean vector.
+static double p_gXMean[] = {
+  0.0025732, -0.284362, -0.566756, -0.946845, 0.361982, -0.52295 };
+
+// Assignment to minimum vector.
+static double p_gXMin[] = {
+  -331.189, -196.066, -257.522, -304.613, -333.804, -292.129 };
+
+// Assignment to maximum vector.
+static double p_gXMax[] = {
+  329.322, 351.904, 931.885, 11176.9, 1147.99, 933.394 };
+
+// Assignment to coefficients vector.
+static double p_gCoefficient[] = {
+  -4150.11,
+  -7349.66,
+  -1376.07,
+  726.65,
+  -604.755,
+  157.175,
+  -7402.61,
+  137.906,
+  1266,
+  607.703,
+  9432.11,
+  2336.63,
+  -12820,
+  -981.145,
+  -2219,
+  -1146.1,
+  1069.02,
+  -9607.95,
+  -135.98,
+  236.78,
+  -5477.97,
+  2213.85,
+  6282.03,
+  89.2823,
+  -10.974,
+  -104.594,
+  -62.7333,
+  -9.89892,
+  -21.8523,
+  -18.8223,
+  16.9863,
+  3546.43,
+  -461.655,
+  -17.9853,
+  -2026.01,
+  -134.173,
+  5565.29,
+  18164.2,
+  5239.56,
+  -4192.18,
+  22.9582,
+  -15479.1,
+  -114.231,
+  2974.51,
+  -10402.5,
+  2516.63,
+  25.4881,
+  303.465,
+  17.2719,
+  1202.72,
+  50.8832,
+  -1.17495,
+  1284.76,
+  17.0545,
+  -0.894684,
+  -136.872,
+  -134.371
+ };
+
+// Assignment to error coefficients vector.
+static double p_gCoefficientRMS[] = {
+  5.0702e-12,
+  9.69107e-12,
+  3.89769e-11,
+  1.13183e-11,
+  1.53099e-11,
+  5.35431e-12,
+  8.94257e-12,
+  1.12034e-11,
+  2.16374e-11,
+  6.32208e-12,
+  6.87455e-11,
+  2.70028e-11,
+  1.70926e-11,
+  2.9263e-11,
+  1.4199e-11,
+  5.24502e-12,
+  6.53991e-12,
+  4.1161e-11,
+  1.19525e-11,
+  1.99626e-11,
+  3.47403e-11,
+  9.32059e-11,
+  1.19932e-10,
+  6.39059e-12,
+  2.49985e-11,
+  4.91273e-11,
+  1.92969e-11,
+  2.95574e-10,
+  1.39773e-11,
+  6.12734e-11,
+  3.16965e-11,
+  1.26652e-10,
+  1.15348e-11,
+  1.5818e-11,
+  1.00252e-11,
+  1.09154e-10,
+  2.29228e-10,
+  1.31396e-10,
+  5.16128e-11,
+  2.71397e-11,
+  2.1153e-10,
+  7.86728e-11,
+  9.25091e-12,
+  1.78175e-10,
+  6.64139e-11,
+  1.25002e-11,
+  8.61292e-11,
+  3.38293e-11,
+  1.16568e-11,
+  1.64392e-10,
+  1.42658e-11,
+  4.97717e-11,
+  1.20839e-11,
+  4.82003e-11,
+  1.9685e-11,
+  4.28747e-11,
+  1.11506e-11
+ };
+
+// Assignment to powers vector.
+// The powers are stored row-wise, that is
+//  p_ij = p_gPower[i * NVariables + j];
+static int    p_gPower[] = {
+  1,  1,  1,  1,  1,  1,
+  1,  1,  1,  1,  1,  2,
+  2,  1,  1,  1,  1,  1,
+  1,  1,  1,  1,  2,  1,
+  1,  2,  1,  1,  1,  1,
+  1,  1,  1,  2,  1,  1,
+  1,  1,  2,  1,  1,  1,
+  1,  1,  1,  1,  1,  3,
+  1,  1,  1,  1,  2,  2,
+  1,  3,  1,  1,  1,  1,
+  2,  1,  2,  1,  1,  1,
+  1,  2,  2,  1,  1,  1,
+  1,  1,  2,  1,  1,  2,
+  1,  2,  1,  1,  1,  2,
+  1,  1,  3,  1,  1,  1,
+  3,  1,  1,  1,  1,  1,
+  1,  1,  1,  1,  3,  1,
+  2,  1,  1,  2,  1,  1,
+  1,  1,  1,  2,  2,  1,
+  1,  1,  2,  1,  2,  1,
+  1,  2,  1,  1,  2,  1,
+  2,  1,  1,  1,  2,  1,
+  2,  2,  1,  1,  1,  1,
+  1,  1,  1,  3,  1,  1,
+  1,  1,  1,  1,  2,  3,
+  2,  1,  1,  3,  1,  1,
+  1,  2,  1,  3,  1,  1,
+  2,  2,  1,  1,  2,  1,
+  1,  3,  1,  1,  2,  1,
+  1,  2,  2,  1,  2,  1,
+  1,  1,  3,  1,  2,  1,
+  2,  2,  1,  2,  1,  1,
+  1,  1,  2,  1,  3,  1,
+  3,  2,  1,  1,  1,  1,
+  3,  1,  1,  1,  1,  2,
+  2,  1,  3,  1,  1,  1,
+  2,  2,  1,  1,  1,  2,
+  2,  1,  2,  1,  1,  2,
+  1,  2,  2,  1,  1,  2,
+  1,  1,  3,  1,  1,  2,
+  2,  2,  2,  1,  1,  1,
+  2,  1,  1,  2,  1,  2,
+  3,  1,  2,  1,  1,  1,
+  2,  1,  1,  1,  2,  2,
+  1,  2,  1,  1,  2,  2,
+  1,  1,  1,  1,  3,  2,
+  2,  1,  1,  1,  1,  3,
+  1,  2,  1,  1,  1,  3,
+  3,  1,  1,  1,  2,  1,
+  2,  1,  2,  1,  2,  1,
+  1,  1,  1,  3,  2,  1,
+  2,  1,  1,  1,  3,  1,
+  1,  3,  1,  1,  1,  2,
+  2,  3,  1,  1,  1,  1,
+  1,  2,  1,  1,  3,  1,
+  1,  2,  3,  1,  1,  1,
+  1,  3,  2,  1,  1,  1
+};
+
+// 
+// The function   double p_fit(double *x)
+// 
+double p_fit(double *x) {
+  double returnValue = p_gDMean;
+  int    i = 0, j = 0, k = 0;
+  for (i = 0; i < p_gNCoefficients ; i++) {
+    // Evaluate the ith term in the expansion
+    double term = p_gCoefficient[i];
+    for (j = 0; j < p_gNVariables; j++) {
+      // Evaluate the polynomial in the jth variable.
+      int power = p_gPower[p_gNVariables * i + j]; 
+      double p1 = 1, p2 = 0, p3 = 0, r = 0;
+      double v =  1 + 2. / (p_gXMax[j] - p_gXMin[j]) * (x[j] - p_gXMax[j]);
+      // what is the power to use!
+      switch(power) {
+      case 1: r = 1; break; 
+      case 2: r = v; break; 
+      default: 
+        p2 = v; 
+        for (k = 3; k <= power; k++) { 
+          p3 = p2 * v;
+          p3 = 2 * v * p2 - p1; 
+          p1 = p2; p2 = p3; 
+        }
+        r = p3;
+      }
+      // multiply this term by the poly in the jth var
+      term *= r; 
+    }
+    // Add this term to the final result
+    returnValue += term;
+  }
+  return returnValue;
+}
+
+// EOF for ../Geant4/samurai/simtree_momentum.C
+// -*- mode: c++ -*-
+// 
+// File ../Geant4/samurai/simtree_rigidity.C generated by TMultiDimFit::MakeRealCode
+// on Tue Jun  9 23:41:17 2020
+// ROOT version 5.32/00
+//
+// This file contains the function 
+//
+//    double  r_fit(double *x); 
+//
+// For evaluating the parameterization obtained
+// from TMultiDimFit and the point x
+// 
+// See TMultiDimFit class documentation for more information 
+// 
+//
+// Static data variables
+//
+static int    r_gNVariables    = 6;
+static int    r_gNCoefficients = 57;
+static double r_gDMean         = 5.24705;
+// Assignment to mean vector.
+static double r_gXMean[] = {
+  0.0025732, -0.284362, -0.566756, -0.946845, 0.361982, -0.52295 };
+
+// Assignment to minimum vector.
+static double r_gXMin[] = {
+  -331.189, -196.066, -257.522, -304.613, -333.804, -292.129 };
+
+// Assignment to maximum vector.
+static double r_gXMax[] = {
+  329.322, 351.904, 931.885, 11176.9, 1147.99, 933.394 };
+
+// Assignment to coefficients vector.
+static double r_gCoefficient[] = {
+  -41.5298,
+  -73.5475,
+  -13.7703,
+  7.27153,
+  -6.05174,
+  1.57283,
+  -74.0774,
+  1.38001,
+  12.6688,
+  6.08124,
+  94.3864,
+  23.3825,
+  -128.289,
+  -9.81824,
+  -22.2054,
+  -11.4689,
+  10.6976,
+  -96.146,
+  -1.36074,
+  2.36944,
+  -54.8176,
+  22.1538,
+  62.8638,
+  0.893442,
+  -0.109816,
+  -1.04666,
+  -0.627767,
+  -0.0990577,
+  -0.218674,
+  -0.188353,
+  0.16998,
+  35.4888,
+  -4.61974,
+  -0.179978,
+  -20.2741,
+  -1.34266,
+  55.6915,
+  181.768,
+  52.4319,
+  -41.9508,
+  0.229741,
+  -154.898,
+  -1.1431,
+  29.7657,
+  -104.097,
+  25.1837,
+  0.255057,
+  3.03675,
+  0.172838,
+  12.0355,
+  0.509184,
+  -0.0117577,
+  12.8565,
+  0.170664,
+  -0.00895304,
+  -1.36967,
+  -1.34464
+ };
+
+// Assignment to error coefficients vector.
+static double r_gCoefficientRMS[] = {
+  5.0702e-12,
+  9.69107e-12,
+  3.89769e-11,
+  1.13183e-11,
+  1.53099e-11,
+  5.35431e-12,
+  8.94257e-12,
+  1.12034e-11,
+  2.16374e-11,
+  6.32208e-12,
+  6.87455e-11,
+  2.70028e-11,
+  1.70926e-11,
+  2.9263e-11,
+  1.4199e-11,
+  5.24502e-12,
+  6.53991e-12,
+  4.1161e-11,
+  1.19525e-11,
+  1.99626e-11,
+  3.47403e-11,
+  9.32059e-11,
+  1.19932e-10,
+  6.39059e-12,
+  2.49985e-11,
+  4.91273e-11,
+  1.92969e-11,
+  2.95574e-10,
+  1.39773e-11,
+  6.12734e-11,
+  3.16965e-11,
+  1.26652e-10,
+  1.15348e-11,
+  1.5818e-11,
+  1.00252e-11,
+  1.09154e-10,
+  2.29228e-10,
+  1.31396e-10,
+  5.16128e-11,
+  2.71397e-11,
+  2.1153e-10,
+  7.86728e-11,
+  9.25091e-12,
+  1.78175e-10,
+  6.64139e-11,
+  1.25002e-11,
+  8.61292e-11,
+  3.38293e-11,
+  1.16568e-11,
+  1.64392e-10,
+  1.42658e-11,
+  4.97717e-11,
+  1.20839e-11,
+  4.82003e-11,
+  1.9685e-11,
+  4.28747e-11,
+  1.11506e-11
+ };
+
+// Assignment to powers vector.
+// The powers are stored row-wise, that is
+//  p_ij = r_gPower[i * NVariables + j];
+static int    r_gPower[] = {
+  1,  1,  1,  1,  1,  1,
+  1,  1,  1,  1,  1,  2,
+  2,  1,  1,  1,  1,  1,
+  1,  1,  1,  1,  2,  1,
+  1,  2,  1,  1,  1,  1,
+  1,  1,  1,  2,  1,  1,
+  1,  1,  2,  1,  1,  1,
+  1,  1,  1,  1,  1,  3,
+  1,  1,  1,  1,  2,  2,
+  1,  3,  1,  1,  1,  1,
+  2,  1,  2,  1,  1,  1,
+  1,  2,  2,  1,  1,  1,
+  1,  1,  2,  1,  1,  2,
+  1,  2,  1,  1,  1,  2,
+  1,  1,  3,  1,  1,  1,
+  3,  1,  1,  1,  1,  1,
+  1,  1,  1,  1,  3,  1,
+  2,  1,  1,  2,  1,  1,
+  1,  1,  1,  2,  2,  1,
+  1,  1,  2,  1,  2,  1,
+  1,  2,  1,  1,  2,  1,
+  2,  1,  1,  1,  2,  1,
+  2,  2,  1,  1,  1,  1,
+  1,  1,  1,  3,  1,  1,
+  1,  1,  1,  1,  2,  3,
+  2,  1,  1,  3,  1,  1,
+  1,  2,  1,  3,  1,  1,
+  2,  2,  1,  1,  2,  1,
+  1,  3,  1,  1,  2,  1,
+  1,  2,  2,  1,  2,  1,
+  1,  1,  3,  1,  2,  1,
+  2,  2,  1,  2,  1,  1,
+  1,  1,  2,  1,  3,  1,
+  3,  2,  1,  1,  1,  1,
+  3,  1,  1,  1,  1,  2,
+  2,  1,  3,  1,  1,  1,
+  2,  2,  1,  1,  1,  2,
+  2,  1,  2,  1,  1,  2,
+  1,  2,  2,  1,  1,  2,
+  1,  1,  3,  1,  1,  2,
+  2,  2,  2,  1,  1,  1,
+  2,  1,  1,  2,  1,  2,
+  3,  1,  2,  1,  1,  1,
+  2,  1,  1,  1,  2,  2,
+  1,  2,  1,  1,  2,  2,
+  1,  1,  1,  1,  3,  2,
+  2,  1,  1,  1,  1,  3,
+  1,  2,  1,  1,  1,  3,
+  3,  1,  1,  1,  2,  1,
+  2,  1,  2,  1,  2,  1,
+  1,  1,  1,  3,  2,  1,
+  2,  1,  1,  1,  3,  1,
+  1,  3,  1,  1,  1,  2,
+  2,  3,  1,  1,  1,  1,
+  1,  2,  1,  1,  3,  1,
+  1,  2,  3,  1,  1,  1,
+  1,  3,  2,  1,  1,  1
+};
+
+// 
+// The function   double r_fit(double *x)
+// 
+double r_fit(double *x) {
+  double returnValue = r_gDMean;
+  int    i = 0, j = 0, k = 0;
+  for (i = 0; i < r_gNCoefficients ; i++) {
+    // Evaluate the ith term in the expansion
+    double term = r_gCoefficient[i];
+    for (j = 0; j < r_gNVariables; j++) {
+      // Evaluate the polynomial in the jth variable.
+      int power = r_gPower[r_gNVariables * i + j]; 
+      double p1 = 1, p2 = 0, p3 = 0, r = 0;
+      double v =  1 + 2. / (r_gXMax[j] - r_gXMin[j]) * (x[j] - r_gXMax[j]);
+      // what is the power to use!
+      switch(power) {
+      case 1: r = 1; break; 
+      case 2: r = v; break; 
+      default: 
+        p2 = v; 
+        for (k = 3; k <= power; k++) { 
+          p3 = p2 * v;
+          p3 = 2 * v * p2 - p1; 
+          p1 = p2; p2 = p3; 
+        }
+        r = p3;
+      }
+      // multiply this term by the poly in the jth var
+      term *= r; 
+    }
+    // Add this term to the final result
+    returnValue += term;
+  }
+  return returnValue;
+}
+
+// EOF for ../Geant4/samurai/simtree_rigidity.C
+// -*- mode: c++ -*-
+// 
+// File ../Geant4/samurai/simtree_length.C generated by TMultiDimFit::MakeRealCode
+// on Tue Jun  9 23:41:17 2020
+// ROOT version 5.32/00
+//
+// This file contains the function 
+//
+//    double  l_fit(double *x); 
+//
+// For evaluating the parameterization obtained
+// from TMultiDimFit and the point x
+// 
+// See TMultiDimFit class documentation for more information 
+// 
+//
+// Static data variables
+//
+static int    l_gNVariables    = 6;
+static int    l_gNCoefficients = 56;
+static double l_gDMean         = 9302.09;
+// Assignment to mean vector.
+static double l_gXMean[] = {
+  0.0025732, -0.284362, -0.566756, -0.946845, 0.361982, -0.52295 };
+
+// Assignment to minimum vector.
+static double l_gXMin[] = {
+  -331.189, -196.066, -257.522, -304.613, -333.804, -292.129 };
+
+// Assignment to maximum vector.
+static double l_gXMax[] = {
+  329.322, 351.904, 931.885, 11176.9, 1147.99, 933.394 };
+
+// Assignment to coefficients vector.
+static double l_gCoefficient[] = {
+  -5063.74,
+  -8210.4,
+  8531.82,
+  -734.412,
+  -360.971,
+  -8616.61,
+  331.099,
+  -854.689,
+  13989.9,
+  -14949,
+  -2334.33,
+  23725.3,
+  -232.89,
+  863.204,
+  -5135.67,
+  -671.923,
+  193.744,
+  -7892.44,
+  5470.48,
+  2106.04,
+  220.14,
+  392.627,
+  -750.825,
+  -13.1656,
+  -15.4854,
+  -15.0127,
+  -692.948,
+  -46.1235,
+  -10.895,
+  -2288.18,
+  13.9869,
+  -7.15341,
+  0.581274,
+  42.6349,
+  -8.22167,
+  -116.588,
+  8175.36,
+  1080.47,
+  129.214,
+  26910.7,
+  7002.79,
+  -4598.69,
+  -14.7257,
+  2944.04,
+  -467.226,
+  -15008,
+  2105.01,
+  6.06541,
+  361.147,
+  -1073.29,
+  3745.34,
+  -284.423,
+  1666.15,
+  -499.767,
+  8786.59,
+  119.522
+ };
+
+// Assignment to error coefficients vector.
+static double l_gCoefficientRMS[] = {
+  1092.36,
+  1874.66,
+  2316.91,
+  992.765,
+  293.062,
+  1327.01,
+  164.035,
+  1469.78,
+  1822.87,
+  2291.78,
+  587.934,
+  2842.24,
+  222.135,
+  300.266,
+  2772.37,
+  341.844,
+  284.914,
+  849.513,
+  1225.2,
+  5295.19,
+  203.834,
+  313.389,
+  883.734,
+  23.072,
+  38.4748,
+  15.8304,
+  233.869,
+  387.214,
+  119.356,
+  5419.4,
+  456.337,
+  6.95297,
+  5.94522,
+  60.7779,
+  168.339,
+  150.975,
+  1134.29,
+  329.863,
+  198.355,
+  3547.52,
+  1330.81,
+  1077.78,
+  656.951,
+  1360.04,
+  385.463,
+  1653.15,
+  541.211,
+  56.6718,
+  276.195,
+  707.914,
+  681.038,
+  733.694,
+  700.282,
+  173.163,
+  2230.31,
+  136.351
+ };
+
+// Assignment to powers vector.
+// The powers are stored row-wise, that is
+//  p_ij = l_gPower[i * NVariables + j];
+static int    l_gPower[] = {
+  1,  1,  1,  1,  1,  1,
+  1,  1,  1,  1,  1,  2,
+  2,  1,  1,  1,  1,  1,
+  1,  1,  1,  1,  2,  1,
+  1,  1,  1,  2,  1,  1,
+  1,  1,  2,  1,  1,  1,
+  1,  1,  1,  1,  1,  3,
+  1,  1,  1,  1,  2,  2,
+  2,  1,  2,  1,  1,  1,
+  1,  1,  2,  1,  1,  2,
+  1,  1,  3,  1,  1,  1,
+  2,  1,  1,  1,  1,  2,
+  3,  1,  1,  1,  1,  1,
+  1,  1,  1,  1,  3,  1,
+  2,  1,  1,  2,  1,  1,
+  1,  1,  1,  2,  2,  1,
+  1,  1,  2,  1,  2,  1,
+  1,  2,  1,  1,  2,  1,
+  2,  1,  1,  1,  2,  1,
+  2,  2,  1,  1,  1,  1,
+  1,  1,  1,  3,  1,  1,
+  1,  1,  1,  1,  2,  3,
+  2,  1,  1,  3,  1,  1,
+  3,  1,  1,  1,  2,  1,
+  2,  2,  1,  1,  2,  1,
+  1,  3,  1,  1,  2,  1,
+  1,  3,  1,  2,  1,  1,
+  1,  2,  2,  1,  2,  1,
+  1,  1,  3,  1,  2,  1,
+  2,  2,  1,  2,  1,  1,
+  1,  1,  1,  3,  2,  1,
+  2,  1,  1,  1,  3,  1,
+  1,  2,  1,  1,  3,  1,
+  3,  2,  1,  1,  1,  1,
+  1,  2,  3,  1,  1,  1,
+  2,  1,  3,  1,  1,  1,
+  2,  2,  1,  1,  1,  2,
+  1,  3,  1,  1,  1,  2,
+  1,  3,  2,  1,  1,  1,
+  2,  1,  2,  1,  1,  2,
+  1,  2,  2,  1,  1,  2,
+  1,  1,  3,  1,  1,  2,
+  2,  2,  2,  1,  1,  1,
+  1,  2,  1,  2,  1,  2,
+  3,  1,  2,  1,  1,  1,
+  1,  2,  1,  1,  2,  2,
+  1,  1,  1,  1,  3,  2,
+  2,  3,  1,  1,  1,  1,
+  2,  1,  1,  1,  1,  3,
+  1,  2,  1,  1,  1,  1,
+  1,  2,  2,  1,  1,  1,
+  1,  2,  1,  3,  1,  1,
+  2,  1,  2,  1,  2,  1,
+  1,  1,  2,  1,  3,  1,
+  2,  1,  1,  1,  2,  2,
+  1,  2,  1,  1,  1,  3
+};
+
+// 
+// The function   double l_fit(double *x)
+// 
+double l_fit(double *x) {
+  double returnValue = l_gDMean;
+  int    i = 0, j = 0, k = 0;
+  for (i = 0; i < l_gNCoefficients ; i++) {
+    // Evaluate the ith term in the expansion
+    double term = l_gCoefficient[i];
+    for (j = 0; j < l_gNVariables; j++) {
+      // Evaluate the polynomial in the jth variable.
+      int power = l_gPower[l_gNVariables * i + j]; 
+      double p1 = 1, p2 = 0, p3 = 0, r = 0;
+      double v =  1 + 2. / (l_gXMax[j] - l_gXMin[j]) * (x[j] - l_gXMax[j]);
+      // what is the power to use!
+      switch(power) {
+      case 1: r = 1; break; 
+      case 2: r = v; break; 
+      default: 
+        p2 = v; 
+        for (k = 3; k <= power; k++) { 
+          p3 = p2 * v;
+          p3 = 2 * v * p2 - p1; 
+          p1 = p2; p2 = p3; 
+        }
+        r = p3;
+      }
+      // multiply this term by the poly in the jth var
+      term *= r; 
+    }
+    // Add this term to the final result
+    returnValue += term;
+  }
+  return returnValue;
+}
+
+
+// EOF for ../Geant4/samurai/simtree_length.C
+// -*- mode: c++ -*-
+// 
+// File ../Geant4/samurai/simtree_tof.C generated by TMultiDimFit::MakeRealCode
+// on Tue Jun  9 23:41:17 2020
+// ROOT version 5.32/00
+//
+// This file contains the function 
+//
+//    double  t_fit(double *x); 
+//
+// For evaluating the parameterization obtained
+// from TMultiDimFit and the point x
+// 
+// See TMultiDimFit class documentation for more information 
+// 
+//
+// Static data variables
+//
+static int    t_gNVariables    = 6;
+static int    t_gNCoefficients = 56;
+static double t_gDMean         = 63.5106;
+// Assignment to mean vector.
+static double t_gXMean[] = {
+  0.0025732, -0.284362, -0.566756, -0.946845, 0.361982, -0.52295 };
+
+// Assignment to minimum vector.
+static double t_gXMin[] = {
+  -331.189, -196.066, -257.522, -304.613, -333.804, -292.129 };
+
+// Assignment to maximum vector.
+static double t_gXMax[] = {
+  329.322, 351.904, 931.885, 11176.9, 1147.99, 933.394 };
+
+// Assignment to coefficients vector.
+static double t_gCoefficient[] = {
+  296.575,
+  509.773,
+  -1298.15,
+  -177.102,
+  38.3448,
+  495.622,
+  -7.92753,
+  -88.3852,
+  -44.7681,
+  -204.021,
+  916.644,
+  158.841,
+  -954.907,
+  54.8928,
+  -71.5493,
+  -1160.79,
+  -151.56,
+  -6.48142,
+  -20.8229,
+  -95.551,
+  21.8388,
+  4.53363,
+  -368.211,
+  -11.1978,
+  -0.122623,
+  0.899953,
+  0.664196,
+  -4.26981,
+  -2.76381,
+  -353.964,
+  -41.6085,
+  -0.104552,
+  28.9895,
+  2.23953,
+  8.78264,
+  123.105,
+  3.83715,
+  -338.603,
+  -92.1278,
+  -1248.05,
+  -432.176,
+  299.722,
+  10.2068,
+  17.7278,
+  -20.2405,
+  -135.529,
+  644.109,
+  -168.691,
+  -0.0126702,
+  -23.4315,
+  -633.876,
+  -171.644,
+  -48.5535,
+  7.73714,
+  -3.93291,
+  0.0801397
+ };
+
+// Assignment to error coefficients vector.
+static double t_gCoefficientRMS[] = {
+  90.0633,
+  153.552,
+  286.58,
+  57.4084,
+  10.6353,
+  98.8219,
+  12.2663,
+  110.128,
+  18.1831,
+  52.623,
+  177.01,
+  44.7459,
+  220.758,
+  29.2743,
+  22.5181,
+  280.17,
+  36.8213,
+  66.015,
+  23.1365,
+  90.5589,
+  10.4606,
+  23.5416,
+  99.1486,
+  50.5072,
+  1.90998,
+  3.21059,
+  1.31771,
+  31.4861,
+  9.39573,
+  67.4273,
+  33.8813,
+  0.497038,
+  12.9849,
+  4.9474,
+  13.3374,
+  41.0747,
+  11.7911,
+  86.8346,
+  27.0137,
+  280.477,
+  98.0117,
+  82.1288,
+  53.8616,
+  98.7791,
+  31.4828,
+  165.676,
+  125.204,
+  40.4042,
+  4.27406,
+  11.177,
+  143.484,
+  55.4336,
+  52.6846,
+  15.5019,
+  22.798,
+  0.581552
+ };
+
+// Assignment to powers vector.
+// The powers are stored row-wise, that is
+//  p_ij = t_gPower[i * NVariables + j];
+static int    t_gPower[] = {
+  1,  1,  1,  1,  1,  1,
+  1,  1,  1,  1,  1,  2,
+  2,  1,  1,  1,  1,  1,
+  1,  1,  1,  1,  2,  1,
+  1,  1,  1,  2,  1,  1,
+  1,  1,  2,  1,  1,  1,
+  1,  1,  1,  1,  1,  3,
+  1,  1,  1,  1,  2,  2,
+  1,  3,  1,  1,  1,  1,
+  1,  2,  2,  1,  1,  1,
+  1,  1,  2,  1,  1,  2,
+  1,  1,  3,  1,  1,  1,
+  2,  1,  1,  1,  1,  2,
+  3,  1,  1,  1,  1,  1,
+  1,  1,  1,  1,  3,  1,
+  2,  1,  1,  2,  1,  1,
+  1,  1,  1,  2,  2,  1,
+  1,  2,  1,  2,  1,  1,
+  1,  1,  2,  1,  2,  1,
+  2,  1,  1,  1,  2,  1,
+  1,  1,  1,  3,  1,  1,
+  1,  1,  1,  1,  2,  3,
+  2,  1,  1,  3,  1,  1,
+  1,  2,  1,  3,  1,  1,
+  3,  1,  1,  1,  2,  1,
+  2,  2,  1,  1,  2,  1,
+  1,  3,  1,  1,  2,  1,
+  1,  2,  2,  1,  2,  1,
+  1,  1,  3,  1,  2,  1,
+  1,  2,  1,  2,  2,  1,
+  1,  1,  1,  3,  2,  1,
+  1,  2,  1,  1,  3,  1,
+  1,  1,  2,  1,  3,  1,
+  3,  2,  1,  1,  1,  1,
+  1,  2,  3,  1,  1,  1,
+  3,  1,  1,  1,  1,  2,
+  2,  1,  3,  1,  1,  1,
+  2,  2,  1,  1,  1,  2,
+  1,  3,  1,  1,  1,  2,
+  2,  1,  2,  1,  1,  2,
+  1,  2,  2,  1,  1,  2,
+  1,  1,  3,  1,  1,  2,
+  2,  2,  2,  1,  1,  1,
+  1,  2,  1,  2,  1,  2,
+  3,  1,  2,  1,  1,  1,
+  2,  1,  1,  1,  2,  2,
+  1,  2,  1,  1,  2,  2,
+  1,  1,  1,  1,  3,  2,
+  2,  3,  1,  1,  1,  1,
+  1,  2,  1,  1,  1,  3,
+  2,  1,  2,  1,  1,  1,
+  2,  2,  1,  1,  1,  1,
+  2,  1,  2,  1,  2,  1,
+  1,  3,  2,  1,  1,  1,
+  2,  1,  1,  1,  1,  3,
+  2,  1,  1,  1,  3,  1
+};
+
+// 
+// The function   double t_fit(double *x)
+// 
+double t_fit(double *x) {
+  double returnValue = t_gDMean;
+  int    i = 0, j = 0, k = 0;
+  for (i = 0; i < t_gNCoefficients ; i++) {
+    // Evaluate the ith term in the expansion
+    double term = t_gCoefficient[i];
+    for (j = 0; j < t_gNVariables; j++) {
+      // Evaluate the polynomial in the jth variable.
+      int power = t_gPower[t_gNVariables * i + j]; 
+      double p1 = 1, p2 = 0, p3 = 0, r = 0;
+      double v =  1 + 2. / (t_gXMax[j] - t_gXMin[j]) * (x[j] - t_gXMax[j]);
+      // what is the power to use!
+      switch(power) {
+      case 1: r = 1; break; 
+      case 2: r = v; break; 
+      default: 
+        p2 = v; 
+        for (k = 3; k <= power; k++) { 
+          p3 = p2 * v;
+          p3 = 2 * v * p2 - p1; 
+          p1 = p2; p2 = p3; 
+        }
+        r = p3;
+      }
+      // multiply this term by the poly in the jth var
+      term *= r; 
+    }
+    // Add this term to the final result
+    returnValue += term;
+  }
+  return returnValue;
+}
+
+// EOF for ../Geant4/samurai/simtree_tof.C
+
+
diff --git a/Projects/S034/s034.detector b/Projects/S034/s034.detector
index 19ce7fa1c552d815a5d5a88dd5f83a591b0623b6..6d99cb8d8220a1dd9233417b483fb2f34c2c22e0 100644
--- a/Projects/S034/s034.detector
+++ b/Projects/S034/s034.detector
@@ -4,6 +4,7 @@ SAMURAIBDC 1
   Offset= 0 0 0 mm
   InvertX= 0 
   InvertY= 0
+  InvertD= 1
 
 %%%%%%%%%%%%
 SAMURAIBDC 2
@@ -11,6 +12,7 @@ SAMURAIBDC 2
   Offset= 0 0 0 mm
   InvertX= 0 
   InvertY= 0
+  IverteD= 1
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Minos
@@ -27,12 +29,20 @@ Minos
  XML= db/MINOS.xml
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-SAMURAIFDC2
-  XML= db/SAMURAIFDC2.xml
+SAMURAIFDC0
+ XML= db/SAMURAIFDC0_20200109.xml
+ Offset= 3.005 1.864 0 mm
+ InvertX= 1 
+ InvertY= 0
+ InvertD= 1
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-SAMURAIFDC0
-  XML= db/SAMURAIFDC0_20200109.xml
+SAMURAIFDC2
+ XML= db/SAMURAIFDC2.xml
+ Offset= 0 0 0 mm
+ InvertX= 0 
+ InvertY= 1
+ InvertD= 1
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 SAMURAIHOD