diff --git a/NPLib/Detectors/Sofia/TSofFissionFragment.cxx b/NPLib/Detectors/Sofia/TSofFissionFragment.cxx
index 53be42bac83cd1d2b1727beddc71a5ac6b8f5c96..d8e0195971e97303991c0ee9717372cdb76f7b81 100644
--- a/NPLib/Detectors/Sofia/TSofFissionFragment.cxx
+++ b/NPLib/Detectors/Sofia/TSofFissionFragment.cxx
@@ -45,6 +45,7 @@ TSofFissionFragment::~TSofFissionFragment() {
 void TSofFissionFragment::Clear() {
   fFF_Z.clear();
   fFF_iZ.clear();
+  fFF_Plastic.clear();
   fFF_AoQ.clear();
   fFF_A.clear();
   fFF_Beta.clear();
@@ -53,6 +54,7 @@ void TSofFissionFragment::Clear() {
   fFF_Brho.clear();
   fFF_DT.clear();
   fFF_ThetaIn.clear();
+  fFF_ThetaOut.clear();
   fFF_TofPosX.clear();
   fFF_TofPosY.clear();
   fFF_PosX1.clear();
@@ -63,7 +65,7 @@ void TSofFissionFragment::Clear() {
   fFF_PosY3.clear();
 
   fFF_Zsum = -1;
-  fFF_IntZsum = -1;
+  fFF_iZsum = -1;
 }
 
 
diff --git a/NPLib/Detectors/Sofia/TSofFissionFragment.h b/NPLib/Detectors/Sofia/TSofFissionFragment.h
index c265a6ccaf4da076973a8a0fb059595cdb9ced02..48bd0f0808f261448b2a9e80a1a0769e29c8167b 100644
--- a/NPLib/Detectors/Sofia/TSofFissionFragment.h
+++ b/NPLib/Detectors/Sofia/TSofFissionFragment.h
@@ -36,6 +36,7 @@ class TSofFissionFragment : public TObject {
   private:
     vector<double> fFF_Z; 
     vector<int> fFF_iZ;
+    vector<int> fFF_Plastic;
     vector<double> fFF_AoQ;
     vector<double> fFF_A;
     vector<double> fFF_Beta;
@@ -44,6 +45,7 @@ class TSofFissionFragment : public TObject {
     vector<double> fFF_Brho;
     vector<double> fFF_DT;
     vector<double> fFF_ThetaIn;
+    vector<double> fFF_ThetaOut;
     vector<double> fFF_TofPosX;
     vector<double> fFF_TofPosY;
     vector<double> fFF_PosX1;
@@ -53,7 +55,7 @@ class TSofFissionFragment : public TObject {
     vector<double> fFF_PosY2;
     vector<double> fFF_PosY3;
     double fFF_Zsum;
-    int fFF_IntZsum;
+    int fFF_iZsum;
 
   //////////////////////////////////////////////////////////////
   // Constructor and destructor
@@ -78,7 +80,8 @@ class TSofFissionFragment : public TObject {
   public:
     //////////////////////    SETTERS    ////////////////////////
     inline void SetZsum(double val){fFF_Zsum = val;};//!
-    inline void SetIntZsum(int val){fFF_IntZsum = val;};//!
+    inline void SetiZsum(int val){fFF_iZsum = val;};//!
+    inline void SetPlastic(int val){fFF_Plastic.push_back(val);};//!
     inline void SetZ(double val){fFF_Z.push_back(val);};//!
     inline void SetiZ(int val){fFF_iZ.push_back(val);};//!
     inline void SetAoQ(double val){fFF_AoQ.push_back(val);};//!
@@ -89,6 +92,7 @@ class TSofFissionFragment : public TObject {
     inline void SetBrho(double val){fFF_Brho.push_back(val);};//!
     inline void SetDT(double val){fFF_DT.push_back(val);};//!
     inline void SetThetaIn(double val){fFF_ThetaIn.push_back(val);};//!
+    inline void SetThetaOut(double val){fFF_ThetaOut.push_back(val);};//!
     inline void SetTofPosX(double val){fFF_TofPosX.push_back(val);};//!
     inline void SetTofPosY(double val){fFF_TofPosY.push_back(val);};//!
     inline void SetPosX1(double val){fFF_PosX1.push_back(val);};//!
@@ -106,7 +110,8 @@ class TSofFissionFragment : public TObject {
     int GetMultMwpc2() {return fFF_PosY2.size();}//!
     int GetMultMwpc3() {return fFF_PosY3.size();}//!
     inline double GetZsum() const {return fFF_Zsum;}//! 
-    inline int GetIntZsum() const {return fFF_IntZsum;}//! 
+    inline int GetiZsum() const {return fFF_iZsum;}//! 
+    inline int GetPlastic(int i) const {return fFF_Plastic[i];}//! 
     inline double GetZ(int i) const {return fFF_Z[i];}//! 
     inline int GetiZ(int i) const {return fFF_iZ[i];}//! 
     inline double GetAoQ(int i) const {return fFF_AoQ[i];}//! 
@@ -117,6 +122,7 @@ class TSofFissionFragment : public TObject {
     inline double GetBrho(int i) const {return fFF_Brho[i];}//! 
     inline double GetDT(int i) const {return fFF_DT[i];}//! 
     inline double GetThetaIn(int i) const {return fFF_ThetaIn[i];}//! 
+    inline double GetThetaOut(int i) const {return fFF_ThetaOut[i];}//! 
     inline double GetTofPosX(int i) const {return fFF_TofPosX[i];}//! 
     inline double GetTofPosY(int i) const {return fFF_TofPosY[i];}//! 
     inline double GetPosX1(int i) const {return fFF_PosX1[i];}//! 
diff --git a/NPLib/Detectors/Sofia/TSofMwpcPhysics.cxx b/NPLib/Detectors/Sofia/TSofMwpcPhysics.cxx
index a7351c1f4706fe5cd7254ed080b3a246d21400ea..2bf422432ae545051d531c7c1b1765b8d05d6eff 100644
--- a/NPLib/Detectors/Sofia/TSofMwpcPhysics.cxx
+++ b/NPLib/Detectors/Sofia/TSofMwpcPhysics.cxx
@@ -209,40 +209,52 @@ void TSofMwpcPhysics::BuildPhysicalEvent() {
       }
 
 
-      if(fClusterX.size()==fClusterY.size()){
-        int size = fClusterX.size();
-        vector<double> Xpos;
-        vector<double> Ypos;
-        for(unsigned int i=0; i<size; i++){
-          // *** strip X *** //
-          vector<pair<int,int>> hitX;
-          hitX = fClusterX[i];
-          double x = -1000;
-          int qleft = hitX[0].second;
-          int qmax = hitX[1].second;
-          int qright = hitX[2].second;
-          int padmax = hitX[2].first;
-          if(padmax>0 && padmax+1<288 && qmax>0 && qleft>0 && qright>0){
-            x = GetPositionX(det_num, qmax, padmax, qleft, qright);
-            Xpos.push_back(x);
-          }
-
-          // *** strip Y *** //
-          vector<pair<int,int>> hitY;
-          hitY = fClusterY[i];
-          double y = -1000;
-          int qdown = hitY[0].second;
-          qmax = hitY[1].second;
-          int qup = hitY[2].second;
-          padmax = hitY[2].first;
-          if(padmax>0 && padmax+1<120 && qmax>0 && qdown>0 && qup>0){
-            y = GetPositionY(det_num, qmax, padmax, qdown, qup);
-            Ypos.push_back(y);
-          }
+      // *** strip X *** //
+      int sizeX = fClusterX.size();
+      vector<double> Xpos;
+      for(unsigned int i=0; i<sizeX; i++){
+        vector<pair<int,int>> hitX;
+        hitX = fClusterX[i];
+        double x = -1000;
+        int qleft = hitX[0].second;
+        int qmax = hitX[1].second;
+        int qright = hitX[2].second;
+        int padmax = hitX[2].first;
+        if(padmax>0 && padmax+1<288 && qmax>0 && qleft>0 && qright>0){
+          x = GetPositionX(det_num, qmax, padmax, qleft, qright);
+          Xpos.push_back(x);
+        }
+      }
+
+      // *** strip Y *** //
+      int sizeY = fClusterY.size();
+      vector<double> Ypos;
+      for(unsigned int i=0; i<sizeY; i++){
+        vector<pair<int,int>> hitY;
+        hitY = fClusterY[i];
+        double y = -1000;
+        int qdown = hitY[0].second;
+        int qmax = hitY[1].second;
+        int qup = hitY[2].second;
+        int padmax = hitY[2].first;
+        if(padmax>0 && padmax+1<120 && qmax>0 && qdown>0 && qup>0){
+          y = GetPositionY(det_num, qmax, padmax, qdown, qup);
+          Ypos.push_back(y);
+        }
+      }
+
+      for(unsigned int i=0; i<Xpos.size(); i++){
+        if(Xpos.size()==Ypos.size()){
+          DetectorNbr.push_back(det_num);
+          PositionX1.push_back(Xpos[i]);
+          PositionX2.push_back(-1000);
+          PositionY.push_back(Ypos[i]);
+        }
+        else{
           DetectorNbr.push_back(det_num);
-          PositionX1.push_back(x);
+          PositionX1.push_back(Xpos[i]);
           PositionX2.push_back(-1000);
-          PositionY.push_back(y);
+          PositionY.push_back(-1000);
         }
       }
     }
diff --git a/NPLib/Detectors/Sofia/TSofTwimPhysics.cxx b/NPLib/Detectors/Sofia/TSofTwimPhysics.cxx
index 8904b0635ef71c468ed1f29411e63e92fc09c7d4..98aa802c84c78536b891dbac138a4f7cc1e0a0be 100644
--- a/NPLib/Detectors/Sofia/TSofTwimPhysics.cxx
+++ b/NPLib/Detectors/Sofia/TSofTwimPhysics.cxx
@@ -93,10 +93,15 @@ void TSofTwimPhysics::BuildPhysicalEvent() {
   double Esec2=0;
   double Esec3=0;
   double Esec4=0;
+  double Theta1=-10;
+  double Theta2=-10;
+  double Theta3=-10;
+  double Theta4=-10;
   double DTsec1=-1e6;
   double DTsec2=-1e6;
   double DTsec3=-1e6;
   double DTsec4=-1e6;
+  double AnodeDriftTime[4][16];
 
   unsigned int mysizeE = m_PreTreatedData->GetMultiplicity();
   for (UShort_t e = 0; e < mysizeE ; e++) {
@@ -112,6 +117,8 @@ void TSofTwimPhysics::BuildPhysicalEvent() {
     AnodeEnergy.push_back(Energy);
     AnodeDT.push_back(DT);
 
+    AnodeDriftTime[SectionNbr-1][AnodeNumber-1] = DT;
+
     if(SectionNbr==1){
       anode_energy_sec1.push_back(Energy);
       if(AnodeNumber==8)
@@ -166,24 +173,36 @@ void TSofTwimPhysics::BuildPhysicalEvent() {
     EnergySection.push_back(Esec1);
     DriftTime.push_back(DTsec1);
     SectionNbr.push_back(1);
+
+    Theta1 = asin((AnodeDriftTime[0][11]-AnodeDriftTime[0][3])/(7*31.));
+    Theta.push_back(Theta1);
   }
   if(Esec2>0){
     Esec2 = Cal->ApplyCalibration("SofTwim/SEC2_ALIGN",Esec2);
     EnergySection.push_back(Esec2);
     DriftTime.push_back(DTsec2);
     SectionNbr.push_back(2);
+
+    Theta2 = asin((AnodeDriftTime[1][11]-AnodeDriftTime[1][3])/(7*31.));
+    Theta.push_back(Theta2);
   }
   if(Esec3>0){
     Esec3 = Cal->ApplyCalibration("SofTwim/SEC3_ALIGN",Esec3);
     EnergySection.push_back(Esec3);
     DriftTime.push_back(DTsec3);
     SectionNbr.push_back(3);
+
+    Theta3 = asin((AnodeDriftTime[2][11]-AnodeDriftTime[2][3])/(7*31.));
+    Theta.push_back(Theta3);
   }
   if(Esec4>0){
     Esec4 = Cal->ApplyCalibration("SofTwim/SEC4_ALIGN",Esec4);
     EnergySection.push_back(Esec4);
     DriftTime.push_back(DTsec4);
     SectionNbr.push_back(4);
+
+    Theta4 = asin((AnodeDriftTime[3][11]-AnodeDriftTime[3][3])/(7*31.));
+    Theta.push_back(Theta4);
   }
 
   m_Beta = -1;
@@ -310,6 +329,7 @@ void TSofTwimPhysics::Clear() {
   AnodeSecNbr.clear();
   AnodeEnergy.clear();
   AnodeDT.clear();
+  Theta.clear();
 
   mult1 = 0;
   mult2 = 0;
diff --git a/NPLib/Detectors/Sofia/TSofTwimPhysics.h b/NPLib/Detectors/Sofia/TSofTwimPhysics.h
index 865d2f761de938098ac2790560c9d8ee82f8d4c8..ac10c26c971c73e9ef29cf39fc15f9339c4708e6 100644
--- a/NPLib/Detectors/Sofia/TSofTwimPhysics.h
+++ b/NPLib/Detectors/Sofia/TSofTwimPhysics.h
@@ -63,6 +63,7 @@ class TSofTwimPhysics : public TObject, public NPL::VDetector {
     vector<int>      SectionNbr;
     vector<double>   EnergySection;
     vector<double>   DriftTime;
+    vector<double>   Theta;
     vector<int>      AnodeNbr;
     vector<int>      AnodeSecNbr;
     vector<double>   AnodeEnergy;
diff --git a/Projects/s455/Analysis.cxx b/Projects/s455/Analysis.cxx
index 07fd371ecde3f332e45b8c2b0390babc1798f944..42e6afa0c08b41aef1439abef6f2201658cb6d23 100644
--- a/Projects/s455/Analysis.cxx
+++ b/Projects/s455/Analysis.cxx
@@ -27,6 +27,33 @@ using namespace std;
 #include"NPDetectorManager.h"
 #include"NPPhysicalConstants.h"
 #include"NPGlobalSystemOfUnits.h"
+
+
+////////////////////////////////////////////////////////////////////////////////
+struct TofPair
+{
+  double x = -1000;
+  double y = -1000;
+  double tof = -1000;
+  double velocity = -1;
+  double beta = -1;
+  double theta_in = -10;
+  double theta_out = -10;
+  int plastic = -1;
+  //
+  int isLorR  = -1;
+  // *** isLorR = 1 -> Left
+  // *** isLorR = 2 -> Right 
+  //
+  int isUorD  = -1;
+  // *** isUorD = 1 -> Up
+  // *** isUorD = 2 -> Down
+  int section = -1;
+  double Esec = -1;
+  double DT = -100;
+};
+
+
 ////////////////////////////////////////////////////////////////////////////////
 Analysis::Analysis(){
 }
@@ -66,6 +93,42 @@ void Analysis::TreatEvent(){
     SofTofW->BuildPhysicalEvent();
 
     FissionFragmentAnalysis();
+    //BeamFragmentAnalysis();
+  }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::BeamFragmentAnalysis(){
+  unsigned int softofw_size = SofTofW->PlasticNbr.size();
+
+  double L_CC = 8.45;
+  TofPair TofHit;
+  if(softofw_size==1){
+    TofHit.plastic  = SofTofW->PlasticNbr[0];
+    TofHit.x        = SofTofW->CalPosX[0];
+    TofHit.y        = SofTofW->CalPosY[0];
+    TofHit.tof      = SofTofW->CalTof[0];
+    TofHit.velocity = L_CC/TofHit.tof;
+    TofHit.beta     = TofHit.velocity * m/ns / NPUNITS::c_light;
+
+    double Brho = 9.62543 + 0.0076642*TofHit.x;
+    double Lfactor = 9.17/L_CC;
+    double Beta = TofHit.beta*Lfactor;
+    double Gamma1 = 1. / sqrt(1 - Beta * Beta);
+
+    double AoQ = Brho / (3.10761 * Beta * Gamma1);
+
+    SofFF->SetTOF(TofHit.tof);
+    SofFF->SetTofPosX(TofHit.x);
+    SofFF->SetTofPosY(TofHit.y);
+    SofFF->SetPlastic(TofHit.plastic);
+
+    SofFF->SetBeta(Beta);
+    SofFF->SetGamma(Gamma1);
+    SofFF->SetAoQ(AoQ);
+    SofFF->SetBrho(Brho);
+
+
   }
 }
 
@@ -74,14 +137,6 @@ void Analysis::FissionFragmentAnalysis(){
   unsigned int softofw_size = SofTofW->PlasticNbr.size();
   unsigned int softwim_size = SofTwim->SectionNbr.size();
 
-  double TOF_CC[2];
-  double Plastic[2];
-  double Plastic_left = -1;
-  double Plastic_right = -1;
-  double TOF_left = -1;
-  double TOF_right = -1;
-  double TOF_up = -1;
-  double TOF_down = -1;
   double E1 = -1;
   double E2 = -1;
   double E3 = -1;
@@ -90,15 +145,11 @@ void Analysis::FissionFragmentAnalysis(){
   double DT2 = -1000;
   double DT3 = -1000;
   double DT4 = -1000;
+  double Theta1 = -1000;
+  double Theta2 = -1000;
+  double Theta3 = -1000;
+  double Theta4 = -1000;
   double L_CC = 8.45;
-  double Beta_left = -1;
-  double Beta_right = -1;
-  double Beta_up = -1;
-  double Beta_down = -1;
-  double Beta1 = -1;
-  double Beta2 = -1;
-  double Beta3 = -1;
-  double Beta4 = -1;
   double Beta_norm = 0.745;
   double Gamma1 = -1;
   double Gamma2 = -1;
@@ -116,25 +167,34 @@ void Analysis::FissionFragmentAnalysis(){
   int iZ1 = -1;
   int iZ2 = -1;
   int iZsum = -1;
-  double ThetaIn1 = -1000;
-  double ThetaIn2 = -1000;
-
-  for(int i = 0; i<2; i++){
-    TOF_CC[i] = -1;
-    Plastic[i] = -1;
-  }
 
-  vector<double> PosY;
-  vector<double> PosX;
+  TofPair TofHit[2];
   if(softofw_size==2){
     for(unsigned int i=0; i<softofw_size; i++){
-      TOF_CC[i] = SofTofW->CalTof[i];
-      Plastic[i] = SofTofW->PlasticNbr[i];
-      PosX.push_back(SofTofW->CalPosX[i]);
-      PosY.push_back(SofTofW->CalPosY[i]);
+      TofHit[i].plastic  = SofTofW->PlasticNbr[i];
+      TofHit[i].x        = SofTofW->CalPosX[i];
+      TofHit[i].y        = SofTofW->CalPosY[i];
+      TofHit[i].tof      = SofTofW->CalTof[i];
+      TofHit[i].velocity = L_CC/TofHit[i].tof;
+      TofHit[i].beta     = TofHit[i].velocity * m/ns / NPUNITS::c_light;
+    }
 
-      SofFF->SetTofPosX(SofTofW->CalPosX[i]);
-      SofFF->SetTofPosY(SofTofW->CalPosY[i]);
+    if(TofHit[0].x>TofHit[1].x){
+      TofHit[0].isLorR = 1;
+      TofHit[1].isLorR = 2;
+    }
+    else if(TofHit[0].x<TofHit[1].x){
+      TofHit[0].isLorR = 2;
+      TofHit[1].isLorR = 1;
+    }
+
+    if(TofHit[0].y>TofHit[1].y){
+      TofHit[0].isUorD = 1;
+      TofHit[1].isUorD = 2;
+    }
+    else if(TofHit[0].y<TofHit[1].y){
+      TofHit[0].isUorD = 2;
+      TofHit[1].isUorD = 1;
     }
   }
 
@@ -154,38 +214,31 @@ void Analysis::FissionFragmentAnalysis(){
       SofFF->SetPosY2(SofMwpc->PositionY[i]);
     }
     if(SofMwpc->DetectorNbr[i]==4){
-      X3.push_back(SofMwpc->PositionX1[i]);
-      Y3.push_back(SofMwpc->PositionY[i]);
-
-      SofFF->SetPosX3(SofMwpc->PositionX1[i]);
-      SofFF->SetPosY3(SofMwpc->PositionY[i]);
+      if(SofMwpc->PositionX1[i]!=-1000)
+        X3.push_back(SofMwpc->PositionX1[i]);
+      
+      if(SofMwpc->PositionY[i]!=-1000)
+        Y3.push_back(SofMwpc->PositionY[i]);
+
+      //SofFF->SetPosX3(SofMwpc->PositionX1[i]);
+      //SofFF->SetPosY3(SofMwpc->PositionY[i]);
     }
   }
+  
 
   vector<double> good_posx;
   vector<double> good_posy;
-  for(unsigned int i=0; i<PosX.size(); i++){
-    double tofx = PosX[i];
-    double tofy = PosY[i];
+  for(unsigned int i=0; i<2; i++){
+    double tofx = TofHit[i].x;
+    //double tofy = TofHit[i].y;
     for(unsigned int k=0; k<X3.size(); k++){
       double posx = X3[k];
       if(abs(posx-tofx) < 100){
         good_posx.push_back(posx);
-        good_posy.push_back(tofy);   
       }
     }
   }
 
-
-  /*if(good_posx.size()==2 && good_posy.size()==2){
-    SofFF->SetTofPosX(good_posx[0]);
-    SofFF->SetTofPosX(good_posx[1]);
-
-    SofFF->SetTofPosY(good_posy[0]);
-    SofFF->SetTofPosY(good_posy[1]);
-    }*/
-
-
   int mult1 = SofTwim->mult1;
   int mult2 = SofTwim->mult2;
   int mult3 = SofTwim->mult3;
@@ -202,18 +255,22 @@ void Analysis::FissionFragmentAnalysis(){
         if(sec==1){
           E1 = SofTwim->EnergySection[i];
           DT1 = SofTwim->DriftTime[i];
+          Theta1 = SofTwim->Theta[i];
         }
         if(sec==2){
           E2 = SofTwim->EnergySection[i];
           DT2 = SofTwim->DriftTime[i];
+          Theta2 = SofTwim->Theta[i];
         }
         if(sec==3){
           E3 = SofTwim->EnergySection[i]; 
           DT3 = SofTwim->DriftTime[i];
+          Theta3 = SofTwim->Theta[i];
         }
         if(sec==4){
           E4 = SofTwim->EnergySection[i];     
           DT4 = SofTwim->DriftTime[i];
+          Theta4 = SofTwim->Theta[i];
         }
       }
 
@@ -237,122 +294,239 @@ void Analysis::FissionFragmentAnalysis(){
       if(E4>0)
         E4 = E4/16;
 
-      if(Plastic[0]<Plastic[1]){
-        Plastic_left = Plastic[0];
-        Plastic_right = Plastic[1];
-        TOF_left = TOF_CC[0];
-        TOF_right = TOF_CC[1];
-      }
-      if(Plastic[0]>Plastic[1]){
-        Plastic_left = Plastic[1];
-        Plastic_right = Plastic[0];
-        TOF_left = TOF_CC[1];
-        TOF_right = TOF_CC[0];
-      }
 
-      if(PosY.size()==2){
-        if(PosY[0]>PosY[1]){
-          TOF_up = TOF_CC[0];
-          TOF_down = TOF_CC[1];
+      // *** case 1 *** //
+      if(E1!=-1 && E2!=-1){
+        if(TofHit[0].isUorD==1 && TofHit[1].isUorD==2){
+          TofHit[0].Esec = E2;
+          TofHit[1].Esec = E1;
+
+          TofHit[0].theta_in = Theta2;
+          TofHit[1].theta_in = Theta1;
+ 
+          TofHit[0].DT = DT2;
+          TofHit[1].DT = DT1;
+
+          TofHit[0].section = 2;
+          TofHit[1].section = 1;
         }
-        if(PosY[0]<PosY[1]){
-          TOF_up = TOF_CC[1];
-          TOF_down = TOF_CC[0];
+        if(TofHit[0].isUorD==2 && TofHit[1].isUorD==1){
+          TofHit[0].Esec = E1;
+          TofHit[1].Esec = E2;
+
+          TofHit[0].theta_in = Theta1;
+          TofHit[1].theta_in = Theta2;
+ 
+          TofHit[0].DT = DT1;
+          TofHit[1].DT = DT2;
+
+          TofHit[0].section = 1;
+          TofHit[1].section = 2;
         }
       }
 
-      double velocity_left = L_CC/TOF_left;
-      double velocity_right = L_CC/TOF_right;
-      Beta_left = velocity_left * m/ns / NPUNITS::c_light;
-      Beta_right = velocity_right * m/ns / NPUNITS::c_light;
-
-      double velocity_down = L_CC/TOF_down;
-      double velocity_up = L_CC/TOF_up;
-      Beta_down = velocity_down * m/ns / NPUNITS::c_light;
-      Beta_up = velocity_up * m/ns / NPUNITS::c_light;
-
-      if(E1!=-1 && E2!=-1){
-        Beta1 = Beta_down;
-        Beta2 = Beta_up;
-      }
+      // *** case 2 *** //
       if(E1!=-1 && E3!=-1){
-        Beta1 = Beta_down;
-        Beta3 = Beta_up;
+        if(TofHit[0].isUorD==1 && TofHit[1].isUorD==2){
+          TofHit[0].Esec = E3;
+          TofHit[1].Esec = E1;
+
+          TofHit[0].theta_in = Theta3;
+          TofHit[1].theta_in = Theta1;
+ 
+          TofHit[0].DT = DT3;
+          TofHit[1].DT = DT1;
+
+          TofHit[0].section = 3;
+          TofHit[1].section = 1;
+        }
+        if(TofHit[0].isUorD==2 && TofHit[1].isUorD==1){
+          TofHit[0].Esec = E1;
+          TofHit[1].Esec = E3;
+
+          TofHit[0].DT = DT1;
+          TofHit[1].DT = DT3;
+
+          TofHit[0].theta_in = Theta1;
+          TofHit[1].theta_in = Theta3;
+ 
+          TofHit[0].section = 1;
+          TofHit[1].section = 3; 
+        }
       }
+
+      // *** case 3 *** //
       if(E1!=-1 && E4!=-1){
-        Beta1 = Beta_left;
-        Beta4 = Beta_right;
+        if(TofHit[0].isLorR==1 && TofHit[1].isLorR==2){
+          TofHit[0].Esec = E1;
+          TofHit[1].Esec = E4;
+
+          TofHit[0].theta_in = Theta1;
+          TofHit[1].theta_in = Theta4;
+          
+          TofHit[0].DT = DT1;
+          TofHit[1].DT = DT4;
+
+          TofHit[0].section = 1;
+          TofHit[1].section = 4;
+        }
+        if(TofHit[0].isLorR==2 && TofHit[1].isLorR==1){
+          TofHit[0].Esec = E4;
+          TofHit[1].Esec = E1;
+
+          TofHit[0].theta_in = Theta4;
+          TofHit[1].theta_in = Theta1;
+           
+          TofHit[0].DT = DT4;
+          TofHit[1].DT = DT1;
+
+          TofHit[0].section = 4;
+          TofHit[1].section = 1;
+        }
       }
+
+      // *** case 4 *** //
       if(E2!=-1 && E3!=-1){
-        Beta2 = Beta_left;
-        Beta3 = Beta_right;
+        if(TofHit[0].isLorR==1 && TofHit[1].isLorR==2){
+          TofHit[0].Esec = E2;
+          TofHit[1].Esec = E3;
+
+          TofHit[0].theta_in = Theta2;
+          TofHit[1].theta_in = Theta3;
+           
+          TofHit[0].DT = DT2;
+          TofHit[1].DT = DT3;
+
+          TofHit[0].section = 2;
+          TofHit[1].section = 3;
+        }
+        if(TofHit[0].isLorR==2 && TofHit[1].isLorR==1){
+          TofHit[0].Esec = E3;
+          TofHit[1].Esec = E2;
+
+          TofHit[0].theta_in = Theta3;
+          TofHit[1].theta_in = Theta2;
+ 
+          TofHit[0].DT = DT3;
+          TofHit[1].DT = DT2;
+
+          TofHit[0].section = 3;
+          TofHit[1].section = 2;
+        }
       }
+
+      // *** case 5 *** //
       if(E2!=-1 && E4!=-1){
-        Beta2 = Beta_up;
-        Beta4 = Beta_down;
+        if(TofHit[0].isUorD==1 && TofHit[1].isUorD==2){
+          TofHit[0].Esec = E2;
+          TofHit[1].Esec = E4;
+
+          TofHit[0].theta_in = Theta2;
+          TofHit[1].theta_in = Theta4;
+ 
+          TofHit[0].DT = DT2;
+          TofHit[1].DT = DT4;
+
+          TofHit[0].section = 2;
+          TofHit[1].section = 4;
+        }
+        if(TofHit[0].isUorD==2 && TofHit[1].isUorD==1){
+          TofHit[0].Esec = E4;
+          TofHit[1].Esec = E2;
+
+          TofHit[0].theta_in = Theta4;
+          TofHit[1].theta_in = Theta2;
+ 
+          TofHit[0].DT = DT4;
+          TofHit[1].DT = DT2;
+
+          TofHit[0].section = 4;
+          TofHit[1].section = 2;
+        }
       }
+
+      // *** case 6 *** //
       if(E3!=-1 && E4!=-1){
-        Beta3 = Beta_up;
-        Beta4 = Beta_down;
+        if(TofHit[0].isUorD==1 && TofHit[1].isUorD==2){
+          TofHit[0].Esec = E3;
+          TofHit[1].Esec = E4;
+
+          TofHit[0].theta_in = Theta3;
+          TofHit[1].theta_in = Theta4;
+ 
+          TofHit[0].DT = DT3;
+          TofHit[1].DT = DT4;
+
+          TofHit[0].section = 3;
+          TofHit[1].section = 4;
+        }
+        if(TofHit[0].isUorD==2 && TofHit[1].isUorD==1){
+          TofHit[0].Esec = E4;
+          TofHit[1].Esec = E3;
+
+          TofHit[0].theta_in = Theta4;
+          TofHit[1].theta_in = Theta3;
+ 
+          TofHit[0].DT = DT4;
+          TofHit[1].DT = DT3;
+
+          TofHit[0].section = 4;
+          TofHit[1].section = 3;
+        }
       }
 
 
-      if(Beta1!=-1 && E1>0){ 
-        E1 = E1 / fcorr_z_beta[0]->Eval(Beta1) * fcorr_z_beta[0]->Eval(Beta_norm);
-        E1 = E1 / fcorr_z_dt[0]->Eval(DT1) * fcorr_z_dt[0]->Eval(55);
-      }
-      if(Beta2!=-1 && E2>0){ 
-        E2 = E2 / fcorr_z_beta[1]->Eval(Beta2) * fcorr_z_beta[1]->Eval(Beta_norm);
-        E2 = E2 / fcorr_z_dt[1]->Eval(DT2) * fcorr_z_dt[1]->Eval(55);
-      }
-      if(Beta3!=-1 && E3>0){ 
-        E3 = E3 / fcorr_z_beta[2]->Eval(Beta3) * fcorr_z_beta[2]->Eval(Beta_norm);
-        E3 = E3 / fcorr_z_dt[2]->Eval(DT3) * fcorr_z_dt[2]->Eval(-55);
-      }
-      if(Beta4!=-1 && E4>0){ 
-        E4 = E4 / fcorr_z_beta[3]->Eval(Beta4) * fcorr_z_beta[3]->Eval(Beta_norm);
-        E4 = E4 / fcorr_z_dt[3]->Eval(DT4) * fcorr_z_dt[3]->Eval(-55);
-      }
+      // *** spline correction *** //
+      for(int i=0; i<2; i++){
+        int section = TofHit[i].section;
 
+        double drift_time;
+        if(section==1) drift_time = DT1;
+        if(section==2) drift_time = DT2;
+        if(section==3) drift_time = DT3;
+        if(section==4) drift_time = DT4;
 
-      // Z calibration //
-      if(E1>0 && E2>0 && E3==-1 && E4==-1){
-        Z1 = E1;
-        Z2 = E2;
-        Beta_Z1 = Beta1;
-        Beta_Z2 = Beta2;
-      }
-      if(E1>0 && E2==-1 && E3>0 && E4==-1){
-        Z1 = E1;
-        Z2 = E3;
-        Beta_Z1 = Beta1;
-        Beta_Z2 = Beta3;
-      }
-      if(E1>0 && E2==-1 && E3==-1 && E4>0){
-        Z1 = E1;
-        Z2 = E4;
-        Beta_Z1 = Beta1;
-        Beta_Z2 = Beta4;
-      }
-      if(E1==-1 && E2>0 && E3>0 && E4==-1){
-        Z1 = E2;
-        Z2 = E3;
-        Beta_Z1 = Beta2;
-        Beta_Z2 = Beta3;
-      }
-      if(E1==-1 && E2>0 && E3==-1 && E4>0){
-        Z1 = E2;
-        Z2 = E4;
-        Beta_Z1 = Beta2;
-        Beta_Z2 = Beta4;
+        double DT_eval;
+        if(section<3) DT_eval=55;
+        if(section>2) DT_eval=-55;
+        if(section>0){
+          TofHit[i].Esec = TofHit[i].Esec / fcorr_z_beta[section-1]->Eval(TofHit[i].beta) * fcorr_z_beta[section-1]->Eval(Beta_norm);
+
+          TofHit[i].Esec = TofHit[i].Esec / fcorr_z_dt[section-1]->Eval(drift_time) * fcorr_z_dt[section-1]->Eval(DT_eval);
+        }
       }
-      if(E1==-1 && E2==-1 && E3>0 && E4>0){
-        Z1 = E3;
-        Z2 = E4;
-        Beta_Z1 = Beta3;
-        Beta_Z2 = Beta4;
+
+
+      // *** Calculation Theta_out *** //
+      double XA;
+      double ZA = 2328.;
+      double XC;
+      double ZC = 4434.;
+      double XMW3 = -1436;
+      double ZMW3 = 8380;
+      double X3lab;
+      double Z3lab;
+      double Theta0 = 20.*deg;
+      TVector3 vOut;
+      TVector3 vZ = TVector3(0,0,1);
+      for(int i=0; i<2; i++){
+        XA = TofHit[i].DT;
+        XC = XA+(ZC-ZA)*tan(TofHit[i].theta_in);
+
+        X3lab = TofHit[i].x*cos(Theta0) + XMW3;
+        Z3lab = ZMW3 + TofHit[i].x*cos(Theta0);
+
+        vOut = TVector3(X3lab-XC,0,Z3lab-ZC);
+
+        double angle = vZ.Angle(vOut);
+
+        TofHit[i].theta_out = angle;
       }
 
+      Z1 = TofHit[0].Esec;
+      Z2 = TofHit[1].Esec;
+
+
       if(Z1>0 && Z2>0){
         Z1 = fZff_p0 + fZff_p1*Z1 + fZff_p2*Z1*Z1;
         Z2 = fZff_p0 + fZff_p1*Z2 + fZff_p2*Z2*Z2;
@@ -367,6 +541,20 @@ void Analysis::FissionFragmentAnalysis(){
         iZsum = iZ1 + iZ2;
       }
 
+      //double Bfactor = 2185./2413.;
+      //Brho1 = Bfactor*(9.62543 + 0.0076642*TofHit[0].x);
+      //Brho2 = Bfactor*(9.62543 + 0.0076642*TofHit[1].x);
+      
+      double MagB = 2185*2.2/3584;
+      double Leff = 2.067;
+      double Tilt = 14*deg;
+      double rho1 = Leff/abs(2*sin(0.5*(TofHit[0].theta_out-TofHit[0].theta_in))*cos(Tilt-0.5*(TofHit[0].theta_out-TofHit[0].theta_in)));
+      double rho2 = Leff/abs(2*sin(0.5*(TofHit[1].theta_out-TofHit[1].theta_in))*cos(Tilt-0.5*(TofHit[1].theta_out-TofHit[1].theta_in)));
+      double Brho1 = MagB*rho1;
+      double Brho2 = MagB*rho2;
+      double Lfactor = 1.;//9.5/L_CC;
+      Beta_Z1 = TofHit[0].beta*Lfactor;
+      Beta_Z2 = TofHit[1].beta*Lfactor;
       Gamma1 = 1. / sqrt(1 - Beta_Z1 * Beta_Z1);
       Gamma2 = 1. / sqrt(1 - Beta_Z2 * Beta_Z2);
 
@@ -377,8 +565,24 @@ void Analysis::FissionFragmentAnalysis(){
       A2 = AoQ2 * iZ2;
 
       // *** Filling the Fission Fragment Tree *** //
-      SofFF->SetTOF(TOF_left);
-      SofFF->SetTOF(TOF_right);
+      SofFF->SetTOF(TofHit[0].tof);
+      SofFF->SetTOF(TofHit[1].tof);
+      SofFF->SetTofPosX(TofHit[0].x);
+      SofFF->SetTofPosX(TofHit[1].x);
+      SofFF->SetTofPosY(TofHit[0].y);
+      SofFF->SetTofPosY(TofHit[1].y);
+      SofFF->SetPlastic(TofHit[0].plastic);
+      SofFF->SetPlastic(TofHit[1].plastic);
+
+      for(int i=0; i<good_posx.size(); i++){
+        SofFF->SetPosX3(good_posx[i]);
+      }
+      SofFF->SetThetaIn(TofHit[0].theta_in);
+      SofFF->SetThetaIn(TofHit[1].theta_in);
+      SofFF->SetThetaOut(TofHit[0].theta_out);
+      SofFF->SetThetaOut(TofHit[1].theta_out);
+
+
       SofFF->SetBeta(Beta_Z1);
       SofFF->SetBeta(Beta_Z2);
       SofFF->SetGamma(Gamma1);
@@ -388,7 +592,7 @@ void Analysis::FissionFragmentAnalysis(){
       SofFF->SetZ(Z1);
       SofFF->SetZ(Z2);
       SofFF->SetAoQ(AoQ1);
-      SofFF->SetAoQ(AoQ1);
+      SofFF->SetAoQ(AoQ2);
       SofFF->SetA(A1);
       SofFF->SetA(A2);
       SofFF->SetBrho(Brho1);
@@ -400,7 +604,7 @@ void Analysis::FissionFragmentAnalysis(){
       SofFF->SetDT(DT4);
 
       SofFF->SetZsum(Zsum);
-      SofFF->SetIntZsum(iZsum);
+      SofFF->SetiZsum(iZsum);
     }
   }
 }
@@ -452,6 +656,7 @@ void Analysis::BeamAnalysis(){
       double Y_p1 = 12.362;
       Zbeam = Zbeam/(Y_p0 + Y_p1*YCC)*Y_p0;
 
+
       // Z calibration //
       Zbeam = fZbeam_p0 + fZbeam_p1*Zbeam + fZbeam_p2*Zbeam*Zbeam;
       Zbeam = sqrt(Zbeam);
@@ -524,11 +729,15 @@ void Analysis::InitParameter(){
   fDCC   = -10000;
   fK_LS2 = -30e-8;
 
-  fRunID = 12;
+  fBrho0 = 12.3255;
+  fRunID = 6;
 
   // Beam parameter //
   fZBeta_p0 = 1;
   fZBeta_p1 = 0;
+  fZbeam_p0 = 1651.57;
+  fZbeam_p1 = 0.0876127;
+  fZbeam_p2 = 4.02563e-6;
 
   // FF parameter //
   fZff_p0 = 2.80063;
@@ -646,8 +855,11 @@ void Analysis::InitParameter(){
     fZbeam_p0 = 186.892;
     fZbeam_p1 = 0.20739;
     fZbeam_p2 = 1.61797e-6;
-
   }
+  if(fRunID==15){
+    fBrho0 = 12.3352;
+  }
+
 }
 
 ////////////////////////////////////////////////////////////////////////////////
diff --git a/Projects/s455/Analysis.h b/Projects/s455/Analysis.h
index db6b3ce1afe7f6f6a28e31bb7bbaa4e5dffc95da..dbd4a351ce058f06b718563e9fc48eb0fcf339e2 100644
--- a/Projects/s455/Analysis.h
+++ b/Projects/s455/Analysis.h
@@ -50,6 +50,7 @@ class Analysis: public NPL::VAnalysis{
     void ReInitValue();
     void BeamAnalysis();
     void FissionFragmentAnalysis();
+    void BeamFragmentAnalysis();
 
     static NPL::VAnalysis* Construct();
 
diff --git a/Projects/s455/RunToTreat.txt b/Projects/s455/RunToTreat.txt
index 9c7e8b6b60cedc31fc59528a06318122575a34aa..c0a52b7fec0f9963b5c0b05cf1fd5094f6383d57 100644
--- a/Projects/s455/RunToTreat.txt
+++ b/Projects/s455/RunToTreat.txt
@@ -2,6 +2,9 @@ TTreeName
   RawTree
 RootFileName
   
+  %/media/sofia/s455/raw/run_raw_0367.root
+  %/media/sofia/s455/raw/run_raw_0368.root
+  %/media/sofia/s455/raw/run_raw_0369.root
   %/media/sofia/s455/raw/run_raw_0419.root
   %/media/sofia/s455/raw/run_raw_0422.root
   
@@ -46,10 +49,10 @@ RootFileName
   %/media/sofia/s455/raw/run_raw_0399.root
 
   %% 189Pb
-  %/media/pierre/proton/raw/run_raw_0400.root
+  /media/sofia/s455/raw/run_raw_0400.root
   
   %% 182Hg
-  /media/pierre/proton/raw/run_raw_0401.root
+  %/media/sofia/s455/raw/run_raw_0401.root
   
   %% 187Pb
   %/media/sofia/s455/raw/run_raw_0402.root
diff --git a/Projects/s455/calibration/SofSci/SofSci_physics.cal b/Projects/s455/calibration/SofSci/SofSci_physics.cal
index 8778ee911a3131c0216c0cc52ad439b9029ffaff..b5d5734487c506da72a4aeca02535cd84b096df9 100644
--- a/Projects/s455/calibration/SofSci/SofSci_physics.cal
+++ b/Projects/s455/calibration/SofSci/SofSci_physics.cal
@@ -1,4 +1,4 @@
-SofSci_TOF2INV_V -8.0515 0.00737389
+SofSci_TOF2INV_V -8.04267 0.00737389
 SofSci_LENGTH_S2 135.614
-SofSci_DET1_POSPAR 45.0 80.0
+SofSci_DET1_POSPAR 90 54.6
 SofSci_DET2_POSPAR 943.287 86.652
diff --git a/Projects/s455/calibration/SofTrim/SofTrim_SectionAlign.cal b/Projects/s455/calibration/SofTrim/SofTrim_SectionAlign.cal
index 37d3efaaa6803eb42745b0b186822b5f6513519a..8a3e3df454ce31ee8ad1bf110603e457707558ff 100644
--- a/Projects/s455/calibration/SofTrim/SofTrim_SectionAlign.cal
+++ b/Projects/s455/calibration/SofTrim/SofTrim_SectionAlign.cal
@@ -1,3 +1,3 @@
-SofTrim_SEC1_ALIGN -457.863 1.03809
+SofTrim_SEC1_ALIGN 310.216 1.00792
 SofTrim_SEC2_ALIGN 0 1
-SofTrim_SEC3_ALIGN -2599.29 1.18757
+SofTrim_SEC3_ALIGN -1382.53 1.13744
diff --git a/Projects/s455/calibration/SofTrim/spline/spline_section_angle_26062021.root b/Projects/s455/calibration/SofTrim/spline/spline_section_angle_26062021.root
new file mode 100644
index 0000000000000000000000000000000000000000..d246cb81becf321f2c03df5e9da6fdb51e9eefe2
Binary files /dev/null and b/Projects/s455/calibration/SofTrim/spline/spline_section_angle_26062021.root differ
diff --git a/Projects/s455/calibration/SofTrim/spline/spline_section_dt_26062021.root b/Projects/s455/calibration/SofTrim/spline/spline_section_dt_26062021.root
new file mode 100644
index 0000000000000000000000000000000000000000..18385be0794ccf796409baf33ecd092ca8e3a903
Binary files /dev/null and b/Projects/s455/calibration/SofTrim/spline/spline_section_dt_26062021.root differ