diff --git a/NPLib/Detectors/Sofia/TSofFissionFragment.cxx b/NPLib/Detectors/Sofia/TSofFissionFragment.cxx
index f5c75a644712852b85358909e5d558e39223020c..5f107d4e3bfa6a19a42338a9c9291979276d228f 100644
--- a/NPLib/Detectors/Sofia/TSofFissionFragment.cxx
+++ b/NPLib/Detectors/Sofia/TSofFissionFragment.cxx
@@ -51,8 +51,10 @@ void TSofFissionFragment::Clear() {
   fFF_TOF.clear();
   fFF_Gamma.clear();
   fFF_Brho.clear();
+  fFF_DT.clear();
   
   fFF_Zsum = -1;
+  fFF_IntZsum = -1;
 }
 
 
diff --git a/NPLib/Detectors/Sofia/TSofFissionFragment.h b/NPLib/Detectors/Sofia/TSofFissionFragment.h
index 069a7f7a58155340825e42c021204467a910f9e3..0fbed2f6ec5a5247586a02ace578809a64dcdc1c 100644
--- a/NPLib/Detectors/Sofia/TSofFissionFragment.h
+++ b/NPLib/Detectors/Sofia/TSofFissionFragment.h
@@ -42,7 +42,9 @@ class TSofFissionFragment : public TObject {
     vector<double> fFF_TOF;
     vector<double> fFF_Gamma;
     vector<double> fFF_Brho;
+    vector<double> fFF_DT;
     double fFF_Zsum;
+    int fFF_IntZsum;
 
   //////////////////////////////////////////////////////////////
   // Constructor and destructor
@@ -67,6 +69,7 @@ class TSofFissionFragment : public TObject {
   public:
     //////////////////////    SETTERS    ////////////////////////
     inline void SetZsum(double val){fFF_Zsum = val;};//!
+    inline void SetIntZsum(int val){fFF_IntZsum = val;};//!
     inline void SetZ(double val){fFF_Z.push_back(val);};//!
     inline void SetQmax(double val){fFF_Qmax.push_back(val);};//!
     inline void SetAoQ(double val){fFF_AoQ.push_back(val);};//!
@@ -75,10 +78,12 @@ class TSofFissionFragment : public TObject {
     inline void SetTOF(double val){fFF_TOF.push_back(val);};//!
     inline void SetGamma(double val){fFF_Gamma.push_back(val);};//!
     inline void SetBrho(double val){fFF_Brho.push_back(val);};//!
+    inline void SetDT(double val){fFF_DT.push_back(val);};//!
 
     //////////////////////    GETTERS    ////////////////////////
     int GetMult() {return fFF_Z.size();}//!
     inline double GetZsum() const {return fFF_Zsum;}//! 
+    inline int GetIntZsum() const {return fFF_IntZsum;}//! 
     inline double GetZ(int i) const {return fFF_Z[i];}//! 
     inline double GetQmax(int i) const {return fFF_Qmax[i];}//! 
     inline double GetAoQ(int i) const {return fFF_AoQ[i];}//! 
@@ -87,6 +92,7 @@ class TSofFissionFragment : public TObject {
     inline double GetTOF(int i) const {return fFF_TOF[i];}//! 
     inline double GetGamma(int i) const {return fFF_Gamma[i];}//! 
     inline double GetBrho(int i) const {return fFF_Brho[i];}//! 
+    inline double GetDT(int i) const {return fFF_DT[i];}//! 
 
   //////////////////////////////////////////////////////////////
   // Required for ROOT dictionnary
diff --git a/NPLib/Detectors/Sofia/TSofSciPhysics.cxx b/NPLib/Detectors/Sofia/TSofSciPhysics.cxx
index 7535e20e96fae4fb8a74b22f38bf5ec00056d00a..94236bb7d94833c3156a839d4e5c056e737d57a2 100644
--- a/NPLib/Detectors/Sofia/TSofSciPhysics.cxx
+++ b/NPLib/Detectors/Sofia/TSofSciPhysics.cxx
@@ -134,6 +134,9 @@ void TSofSciPhysics::BuildPhysicalEvent() {
     }
   }
 
+  if(S2_pmtR.size()==0 || S2_pmtL.size()==0 || CC_pmtR.size()==0 || CC_pmtL.size()==0)
+    return;
+
   if(m_NumberOfDetectors==2){
     double CC_rawpos;
     double S2_rawpos;
diff --git a/NPLib/Detectors/Sofia/TSofTwimPhysics.cxx b/NPLib/Detectors/Sofia/TSofTwimPhysics.cxx
index 0988a97ead1e28c1c1e394e23bb295f5fa8ee888..8904b0635ef71c468ed1f29411e63e92fc09c7d4 100644
--- a/NPLib/Detectors/Sofia/TSofTwimPhysics.cxx
+++ b/NPLib/Detectors/Sofia/TSofTwimPhysics.cxx
@@ -82,24 +82,21 @@ void TSofTwimPhysics::BuildSimplePhysicalEvent() {
 ///////////////////////////////////////////////////////////////////////////
 void TSofTwimPhysics::BuildPhysicalEvent() {
 
-  //if(m_Beta<0)
-  //  return;
-
   // apply thresholds and calibration
   PreTreat();
-  /*if(m_PreTreatedData->GetMultiplicity() != 32){
-  //m_Beta = -1;
-  return;
-  }*/
 
   vector<double> anode_energy_sec1;
   vector<double> anode_energy_sec2;
   vector<double> anode_energy_sec3;
   vector<double> anode_energy_sec4;
-  vector<double> anode_dt_sec1;
-  vector<double> anode_dt_sec2;
-  vector<double> anode_dt_sec3;
-  vector<double> anode_dt_sec4;
+  double Esec1=0;
+  double Esec2=0;
+  double Esec3=0;
+  double Esec4=0;
+  double DTsec1=-1e6;
+  double DTsec2=-1e6;
+  double DTsec3=-1e6;
+  double DTsec4=-1e6;
 
   unsigned int mysizeE = m_PreTreatedData->GetMultiplicity();
   for (UShort_t e = 0; e < mysizeE ; e++) {
@@ -107,6 +104,8 @@ void TSofTwimPhysics::BuildPhysicalEvent() {
     int AnodeNumber = m_PreTreatedData->GetAnodeNbr(e);
     double Energy   = m_PreTreatedData->GetEnergy(e);
     double DT       = m_PreTreatedData->GetDriftTime(e);
+    int PU          = m_PreTreatedData->GetPileUp(e);
+    int OV          = m_PreTreatedData->GetOverflow(e);
 
     AnodeSecNbr.push_back(SectionNbr);
     AnodeNbr.push_back(AnodeNumber);
@@ -115,80 +114,73 @@ void TSofTwimPhysics::BuildPhysicalEvent() {
 
     if(SectionNbr==1){
       anode_energy_sec1.push_back(Energy);
-      anode_dt_sec1.push_back(DT);
+      if(AnodeNumber==8)
+        DTsec1 = DT;
     }
     if(SectionNbr==2){
       anode_energy_sec2.push_back(Energy);
-      anode_dt_sec2.push_back(DT);
+      if(AnodeNumber==8)
+        DTsec2 = DT;
     }
     if(SectionNbr==3){
       anode_energy_sec3.push_back(Energy);
-      anode_dt_sec3.push_back(DT);
+      if(AnodeNumber==8)
+        DTsec3 = DT;
     }
     if(SectionNbr==4){
       anode_energy_sec4.push_back(Energy);
-      anode_dt_sec4.push_back(DT);
+      if(AnodeNumber==8)
+        DTsec4 = DT;
     }
   }
 
-  static CalibrationManager* Cal = CalibrationManager::getInstance();
-  double Esec1=0;
-  double Esec2=0;
-  double Esec3=0;
-  double Esec4=0;
-  double DTsec1=0;
-  double DTsec2=0;
-  double DTsec3=0;
-  double DTsec4=0;
-  for(int i=0; i<anode_energy_sec1.size(); i++){
-    Esec1 += anode_energy_sec1[i];
-    DTsec1 += anode_dt_sec1[i];
+  mult1 = anode_energy_sec1.size();
+  mult2 = anode_energy_sec2.size();
+  mult3 = anode_energy_sec3.size();
+  mult4 = anode_energy_sec4.size();
+
+  if(mult1<17){
+    for(int i=0; i<mult1; i++){
+      Esec1 += anode_energy_sec1[i];
+    }
   }
-  for(int i=0; i<anode_energy_sec2.size(); i++){
-    Esec2 += anode_energy_sec2[i];
-    DTsec2 += anode_dt_sec2[i];
+  if(mult2<17){
+    for(int i=0; i<mult2; i++){
+      Esec2 += anode_energy_sec2[i];
+    }
   }
-  for(int i=0; i<anode_energy_sec3.size(); i++){
-    Esec3 += anode_energy_sec3[i];
-    DTsec3 += anode_dt_sec3[i];
+  if(mult3<17){
+    for(int i=0; i<mult3; i++){
+      Esec3 += anode_energy_sec3[i];
+    }
   }
-  for(int i=0; i<anode_energy_sec4.size(); i++){
-    Esec4 += anode_energy_sec4[i];
-    DTsec4 += anode_dt_sec4[i];
+  if(mult4<17){
+    for(int i=0; i<mult4; i++){
+      Esec4 += anode_energy_sec4[i];
+    }
   }
 
-  if(Esec1>0 && anode_energy_sec1.size()>2){
-    Esec1 = Esec1 / anode_energy_sec1.size();
+  static CalibrationManager* Cal = CalibrationManager::getInstance();
+  if(Esec1>0){
     Esec1 = Cal->ApplyCalibration("SofTwim/SEC1_ALIGN",Esec1);
-
-    DTsec1 = DTsec1 / anode_dt_sec1.size();
     EnergySection.push_back(Esec1);
     DriftTime.push_back(DTsec1);
     SectionNbr.push_back(1);
   }
-  if(Esec2>0 && anode_energy_sec2.size()>2){
-    Esec2 = Esec2 / anode_energy_sec2.size();
+  if(Esec2>0){
     Esec2 = Cal->ApplyCalibration("SofTwim/SEC2_ALIGN",Esec2);
-
-    DTsec2 = DTsec2 / anode_dt_sec2.size();
     EnergySection.push_back(Esec2);
     DriftTime.push_back(DTsec2);
     SectionNbr.push_back(2);
   }
-  if(Esec3>0 && anode_energy_sec3.size()>2){
-    Esec3 = Esec3 / anode_energy_sec3.size();
+  if(Esec3>0){
     Esec3 = Cal->ApplyCalibration("SofTwim/SEC3_ALIGN",Esec3);
-
-    DTsec3 = DTsec3 / anode_dt_sec3.size();
     EnergySection.push_back(Esec3);
     DriftTime.push_back(DTsec3);
     SectionNbr.push_back(3);
   }
-  if(Esec4>0 && anode_energy_sec4.size()>2){
-    Esec4 = Esec4 / anode_energy_sec4.size();
+  if(Esec4>0){
     Esec4 = Cal->ApplyCalibration("SofTwim/SEC4_ALIGN",Esec4);
-
-    DTsec4 = DTsec4 / anode_dt_sec4.size();
     EnergySection.push_back(Esec4);
     DriftTime.push_back(DTsec4);
     SectionNbr.push_back(4);
@@ -318,6 +310,11 @@ void TSofTwimPhysics::Clear() {
   AnodeSecNbr.clear();
   AnodeEnergy.clear();
   AnodeDT.clear();
+
+  mult1 = 0;
+  mult2 = 0;
+  mult3 = 0;
+  mult4 = 0;
 }
 
 
diff --git a/NPLib/Detectors/Sofia/TSofTwimPhysics.h b/NPLib/Detectors/Sofia/TSofTwimPhysics.h
index bbdedc69d2f0580c38164d6702046c71656d16fa..865d2f761de938098ac2790560c9d8ee82f8d4c8 100644
--- a/NPLib/Detectors/Sofia/TSofTwimPhysics.h
+++ b/NPLib/Detectors/Sofia/TSofTwimPhysics.h
@@ -67,6 +67,10 @@ class TSofTwimPhysics : public TObject, public NPL::VDetector {
     vector<int>      AnodeSecNbr;
     vector<double>   AnodeEnergy;
     vector<double>   AnodeDT;
+    int mult1;
+    int mult2;
+    int mult3;
+    int mult4;
 
   /// A usefull method to bundle all operation to add a detector
   void AddDetector(TVector3 POS); 
diff --git a/Projects/s455/Analysis.cxx b/Projects/s455/Analysis.cxx
index b9613b15f9b28943dfd86fc7c843afdb10ca37ac..9b0c02b6bebcc1893b33752eba8951ebd3df73b2 100644
--- a/Projects/s455/Analysis.cxx
+++ b/Projects/s455/Analysis.cxx
@@ -85,23 +85,23 @@ void Analysis::FissionFragmentAnalysis(){
   double TOF_right = -1;
   double TOF_up = -1;
   double TOF_down = -1;
-  double E_left = -1;
-  double E_right = -1;
   double E1 = -1;
   double E2 = -1;
   double E3 = -1;
   double E4 = -1;
-  double DT1 = -1;
-  double DT2 = -1;
-  double DT3 = -1;
-  double DT4 = -1;
-  double E_up = -1;
-  double E_down = -1;
+  double DT1 = -1000;
+  double DT2 = -1000;
+  double DT3 = -1000;
+  double DT4 = -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;
 
   for(int i = 0; i<2; i++){
@@ -118,382 +118,438 @@ void Analysis::FissionFragmentAnalysis(){
     }
   }
 
+  int mult1 = SofTwim->mult1;
+  int mult2 = SofTwim->mult2;
+  int mult3 = SofTwim->mult3;
+  int mult4 = SofTwim->mult4;
+
+  int multL = mult1 + mult2;
+  int multR = mult3 + mult4;
+
   if(softwim_size>1){
-    for(unsigned int i=0; i< softwim_size; i++){
-      int sec = SofTwim->SectionNbr[i];
+    //if( (multL>15 && multR>15) || (multR==32) || (multL==32)){
+    if( (mult1>1 && mult1<17) || (mult2>1 && mult2<17) || (mult3>1 && mult3<17) || (mult4>1 && mult4<17)){
+      for(unsigned int i=0; i< softwim_size; i++){
+        int sec = SofTwim->SectionNbr[i];
+
+        if(sec==1){
+          E1 = SofTwim->EnergySection[i];
+          DT1 = SofTwim->DriftTime[i];
+        }
+        if(sec==2){
+          E2 = SofTwim->EnergySection[i];
+          DT2 = SofTwim->DriftTime[i];
+        }
+        if(sec==3){
+          E3 = SofTwim->EnergySection[i]; 
+          DT3 = SofTwim->DriftTime[i];
+        }
+        if(sec==4){
+          E4 = SofTwim->EnergySection[i];     
+          DT4 = SofTwim->DriftTime[i];
+        }
+      }
 
-      if(sec==1){
-        E1 = SofTwim->EnergySection[i];
-        DT1 = SofTwim->DriftTime[i];
+      if(softwim_size>2){
+        if(E1>0 && E2>0){
+          E1 = E1+E2;
+          E2 = -1;
+        }
+        if(E3>0 && E4>0){
+          E3 = E3+E4;
+          E4 = -1;
+        }
       }
-      else if(sec==2){
-        E2 = SofTwim->EnergySection[i];
-        DT2 = SofTwim->DriftTime[i];
+
+      if(E1>0)
+        E1 = E1/16;
+      if(E2>0)
+        E2 = E2/16;
+      if(E3>0)
+        E3 = E3/16;
+      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];
       }
-      else if(sec==3){
-        E3 = SofTwim->EnergySection[i]; 
-        DT3 = SofTwim->DriftTime[i];
+      if(Plastic[0]>Plastic[1]){
+        Plastic_left = Plastic[1];
+        Plastic_right = Plastic[0];
+        TOF_left = TOF_CC[1];
+        TOF_right = TOF_CC[0];
       }
-      else if(sec==4){
-        E4 = SofTwim->EnergySection[i];     
-        DT4 = SofTwim->DriftTime[i];
+
+      if(PosY[0]>PosY[1]){
+        TOF_up = TOF_CC[0];
+        TOF_down = TOF_CC[1];
+      }
+      if(PosY[0]<PosY[1]){
+        TOF_up = TOF_CC[1];
+        TOF_down = TOF_CC[0];
       }
-    }
-  }
 
-  if(E1>0 && E2>0 && E3>0 && E4>0){
-    E1 = (E1+E2)/2;
-    E2 = -1;
-    E3 = (E3+E4)/2;
-    E4 = -1;
-  }
-  if(E1>0 && E2>0 && (E3>0||E4>0)){
-    E1 = (E1+E2)/2;
-    E2 = -1;
-  }
-  if(E3>0 && E4>0 && (E1>0||E2>0)){
-    E3 = (E3+E4)/2;
-    E4 = -1;
-  }
+      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;
 
-  if(Plastic[0]<Plastic[1]){
-    Plastic_left = Plastic[0];
-    Plastic_right = Plastic[1];
-    TOF_left = TOF_CC[0];
-    TOF_right = TOF_CC[1];
-  }
-  else if(Plastic[0]>Plastic[1]){
-    Plastic_left = Plastic[1];
-    Plastic_right = Plastic[0];
-    TOF_left = TOF_CC[1];
-    TOF_right = TOF_CC[0];
-  }
+      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(PosY[0]>PosY[1]){
-    TOF_up = TOF_CC[0];
-    TOF_down = TOF_CC[1];
-  }
-  else if(PosY[0]<PosY[1]){
-    TOF_up = TOF_CC[1];
-    TOF_down = TOF_CC[0];
-  }
+      if(E1!=-1 && E2!=-1){
+        Beta1 = Beta_down;
+        Beta2 = Beta_up;
+      }
+      if(E1!=-1 && E3!=-1){
+        Beta1 = Beta_down;
+        Beta3 = Beta_up;
+      }
+      if(E1!=-1 && E4!=-1){
+        Beta1 = Beta_left;
+        Beta4 = Beta_right;
+      }
+      if(E2!=-1 && E3!=-1){
+        Beta2 = Beta_left;
+        Beta3 = Beta_right;
+      }
+      if(E2!=-1 && E4!=-1){
+        Beta2 = Beta_up;
+        Beta4 = Beta_down;
+      }
+      if(E3!=-1 && E4!=-1){
+        Beta3 = Beta_up;
+        Beta4 = Beta_down;
+      }
 
-  if(TOF_left != -1 && TOF_right != -1){
-    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;
+      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);
+      }
 
-    if(E1 != -1 && E2==-1){
-      E1 = E1 / fcorr_z_beta[0]->Eval(Beta_left) * fcorr_z_beta[0]->Eval(Beta_norm);
-      E1 = E1 / fcorr_z_dt[0]->Eval(DT1) * fcorr_z_dt[0]->Eval(0);
-      E_left = E1;
-    }
-    if(E2 != -1 && E1==-1){
-      E2 = E2 / fcorr_z_beta[1]->Eval(Beta_left) * fcorr_z_beta[1]->Eval(Beta_norm);
-      E2 = E2 / fcorr_z_dt[1]->Eval(DT2) * fcorr_z_dt[1]->Eval(0);
-      E_left = E2;
-    }
-    if(E3 != -1 && E4==-1){
-      E3 = E3 / fcorr_z_beta[2]->Eval(Beta_right) * fcorr_z_beta[2]->Eval(Beta_norm);
-      E3 = E3 / fcorr_z_dt[2]->Eval(DT3) * fcorr_z_dt[2]->Eval(0);
-      E_right = E3;
-    }
-    if(E4 != -1 && E3==-1){
-      E4 = E4 / fcorr_z_beta[3]->Eval(Beta_right) * fcorr_z_beta[3]->Eval(Beta_norm);
-      E4 = E4 / fcorr_z_dt[3]->Eval(DT4) * fcorr_z_dt[3]->Eval(0);
-      E_right = E4;
-    }
-  }
-  
-  if(TOF_up != -1 && TOF_down != -1){
-    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>0 && E2>0 && E3==-1 && E4==-1){
-      E1 = E1 / fcorr_z_beta[0]->Eval(Beta_down) * fcorr_z_beta[0]->Eval(Beta_norm);
-      E2 = E2 / fcorr_z_beta[1]->Eval(Beta_up) * fcorr_z_beta[1]->Eval(Beta_norm);
-      
-      E_up = E2;
-      E_down = E1;
-    }
-    else if(E1==-1 && E2==-1 && E3>0 && E4>0){
-      E3 = E3 / fcorr_z_beta[2]->Eval(Beta_up) * fcorr_z_beta[2]->Eval(Beta_norm);
-      E4 = E4 / fcorr_z_beta[3]->Eval(Beta_down) * fcorr_z_beta[3]->Eval(Beta_norm);
-      
-      E_up = E3;
-      E_down = E4;
-    }
-  }
 
+      // Z calibration //
+      double p0 = 2.80063;
+      double p1 = 6.91985e-2;
+      double p2 = 1.01598e-7;
+      double Z1=-1;
+      double Z2=-1;
+      double Zsum=-1;
+      int iZsum=-1;
+
+      if(E1>0 && E2>0 && E3==-1 && E4==-1){
+        Z1 = E1;
+        Z2 = E2;
+      }
+      if(E1>0 && E2==-1 && E3>0 && E4==-1){
+        Z1 = E1;
+        Z2 = E3;
+      }
+      if(E1>0 && E2==-1 && E3==-1 && E4>0){
+        Z1 = E1;
+        Z2 = E4;
+      }
+      if(E1==-1 && E2>0 && E3>0 && E4==-1){
+        Z1 = E2;
+        Z2 = E3;
+      }
+      if(E1==-1 && E2>0 && E3==-1 && E4>0){
+        Z1 = E2;
+        Z2 = E4;
+      }
+      if(E1==-1 && E2==-1 && E3>0 && E4>0){
+        Z1 = E3;
+        Z2 = E4;
+      }
+
+      if(Z1>0 && Z2>0){
+        Z1 = p0 + p1*Z1 + p2*Z1*Z1;
+        Z2 = p0 + p1*Z2 + p2*Z2*Z2;
+        
+        Z1 = sqrt(Z1);
+        Z2 = sqrt(Z2);
 
-  // Z calibration //
-  double p0 = -0.787072;
-  double p1 = 0.273064;
-  double Z1=-1;
-  double Z2=-1;
-  double Zsum=-1;
-  if(E_left!=-1 && E_right!=-1){
-    Z1 = p0 + p1*sqrt(E_left);
-    Z2 = p0 + p1*sqrt(E_right);
-    Zsum = Z1+Z2;
-  }
-  if(E_up!=-1 && E_down!=-1){
-    double Z1 = p0 + p1*sqrt(E_down);
-    double Z2 = p0 + p1*sqrt(E_up);
-    Zsum = Z1+Z2;
-  }
+        Zsum = Z1+Z2;
 
-  if(E_left != -1 && E_right != -1){
-    SofFF->SetTOF(TOF_left);
-    SofFF->SetTOF(TOF_right);
-    SofFF->SetBeta(Beta_left);
-    SofFF->SetBeta(Beta_right);
-    SofFF->SetZ(p0+p1*sqrt(E1));
-    SofFF->SetZ(p0+p1*sqrt(E2));
-    SofFF->SetZ(p0+p1*sqrt(E3));
-    SofFF->SetZ(p0+p1*sqrt(E4));
-    SofFF->SetZsum(Zsum);
+        int iZ1 = (int) round(Z1);
+        int iZ2 = (int) round(Z2);
+        iZsum = iZ1 + iZ2;
+      }
+
+      SofFF->SetTOF(TOF_left);
+      SofFF->SetTOF(TOF_right);
+      SofFF->SetBeta(Beta1);
+      SofFF->SetBeta(Beta2);
+      SofFF->SetBeta(Beta3);
+      SofFF->SetBeta(Beta4);
+
+      /*SofFF->SetZ(E1);
+      SofFF->SetZ(E2);
+      SofFF->SetZ(E3);
+      SofFF->SetZ(E4);*/
+
+      /*SofFF->SetZ( sqrt(p0+p1*E1+p2*E1*E1) );
+      SofFF->SetZ( sqrt(p0+p1*E2+p2*E2*E2) );
+      SofFF->SetZ( sqrt(p0+p1*E3+p2*E3*E3) );
+      SofFF->SetZ( sqrt(p0+p1*E4+p2*E4*E4) );*/
+      SofFF->SetZ(Z1);
+      SofFF->SetZ(Z2);
+
+      SofFF->SetDT(DT1);
+      SofFF->SetDT(DT2);
+      SofFF->SetDT(DT3);
+      SofFF->SetDT(DT4);
+
+      SofFF->SetZsum(Zsum);
+      SofFF->SetIntZsum(iZsum);
+    }
   }
-  if(E_up!=-1 && E_down!=-1){
-    SofFF->SetTOF(TOF_down);
-    SofFF->SetTOF(TOF_up);
-    SofFF->SetBeta(Beta_down);
-    SofFF->SetBeta(Beta_up);
-    SofFF->SetZ(p0+p1*sqrt(E1));
-    SofFF->SetZ(p0+p1*sqrt(E2));
-    SofFF->SetZ(p0+p1*sqrt(E3));
-    SofFF->SetZ(p0+p1*sqrt(E4));
-    SofFF->SetZsum(Zsum);
   }
-}
 
-////////////////////////////////////////////////////////////////////////////////
-void Analysis::BeamAnalysis(){
-  unsigned int sofsci_size = SofSci->DetectorNbr.size();
-  double Y_p0 = 23943.8;
-  double Y_p1 = 12.362;
-
-  double Z_p0 = -9.31873;
-  double Z_p1 = 0.564459;
-  
-  if(sofsci_size==2){
-    double beta = SofSci->Beta[0];
-    //cout << "Set beta to " << beta << endl;
-    SofTrim->SetBeta(beta);
-    SofTrim->BuildSimplePhysicalEvent();
-    double Zbeam,Qmax,Theta;
-    if(SofTrim->EnergySection.size()>0){
-      Zbeam = SofTrim->GetMaxEnergySection();
-      Qmax = DetermineQmax();
-      Theta = SofTrim->Theta[0];
-
-      double TofFromS2    = SofSci->CalTof[0];
-      double velocity_mns = SofSci->VelocityMNs[0];
-      double Beta         = SofSci->Beta[0];
-      double XS2          = SofSci->PosMm[0];
-      //double XCC          = SofSci->PosMm[1];
-      double XCC=0;
-      double YCC=0;
-      for(unsigned int i=0; i<SofMwpc->DetectorNbr.size(); i++){
-        if(SofMwpc->DetectorNbr[i]==1){
-          XCC = SofMwpc->PositionX1[i];
-          YCC = SofMwpc->PositionY[i];
+  ////////////////////////////////////////////////////////////////////////////////
+  void Analysis::BeamAnalysis(){
+    unsigned int sofsci_size = SofSci->DetectorNbr.size();
+    double Y_p0 = 23943.8;
+    double Y_p1 = 12.362;
+
+    double alpha = 2.09862;
+    double beta = 0.972689;
+    double Z_p0 = alpha+beta*-9.31873;
+    double Z_p1 = beta*0.564459;
+
+    if(sofsci_size==2){
+      double beta = SofSci->Beta[0];
+      //cout << "Set beta to " << beta << endl;
+      SofTrim->SetBeta(beta);
+      SofTrim->BuildSimplePhysicalEvent();
+      double Zbeam,Qmax,Theta;
+      if(SofTrim->EnergySection.size()>0){
+        Zbeam = SofTrim->GetMaxEnergySection();
+        Qmax = DetermineQmax();
+        Theta = SofTrim->Theta[0];
+
+        double TofFromS2    = SofSci->CalTof[0];
+        double velocity_mns = SofSci->VelocityMNs[0];
+        double Beta         = SofSci->Beta[0];
+        double XS2          = SofSci->PosMm[0];
+        //double XCC          = SofSci->PosMm[1];
+        double XCC=0;
+        double YCC=0;
+        for(unsigned int i=0; i<SofMwpc->DetectorNbr.size(); i++){
+          if(SofMwpc->DetectorNbr[i]==1){
+            XCC = SofMwpc->PositionX1[i];
+            YCC = SofMwpc->PositionY[i];
+          }
         }
+
+        double LS2;
+        LS2 = fLS2_0;//*(1 + fK_LS2*Theta);
+        velocity_mns = LS2/TofFromS2;
+        Beta = velocity_mns * m/ns / NPUNITS::c_light;
+        double Gamma        = 1./(TMath::Sqrt(1 - TMath::Power(Beta,2)));
+        double Brho = fBrho0 * (1 - XS2/fDS2 - XCC/fDCC);
+        double AoQ  = Brho / (3.10716*Gamma*Beta);
+        double A    = AoQ * Qmax;
+
+        // Y dependence correction //
+        Zbeam = Zbeam/(Y_p0 + Y_p1*YCC)*Y_p0;
+
+        // Z calibration //
+        Zbeam = Z_p0 + Z_p1*sqrt(Zbeam);
+        //Zbeam = 2.09862 + 0.972689*Zbeam;
+
+        // Filling Beam tree
+        SofBeamID->SetZbeam(Zbeam);
+        SofBeamID->SetQmax(rand.Gaus(Qmax,0.15));
+        SofBeamID->SetAoQ(AoQ);
+        SofBeamID->SetAbeam(A);
+        SofBeamID->SetBeta(Beta);
+        SofBeamID->SetGamma(Gamma);
+        SofBeamID->SetBrho(Brho);
+        SofBeamID->SetXS2(XS2);
+        SofBeamID->SetXCC(XCC);
+        SofBeamID->SetYCC(YCC);
       }
-  
-      double LS2;
-      LS2 = fLS2_0;//*(1 + fK_LS2*Theta);
-      velocity_mns = LS2/TofFromS2;
-      Beta = velocity_mns * m/ns / NPUNITS::c_light;
-      double Gamma        = 1./(TMath::Sqrt(1 - TMath::Power(Beta,2)));
-      double Brho = fBrho0 * (1 - XS2/fDS2 - XCC/fDCC);
-      double AoQ  = Brho / (3.10716*Gamma*Beta);
-      double A    = AoQ * Qmax;
-
-      // Y dependence correction //
-      Zbeam = Zbeam/(Y_p0 + Y_p1*YCC)*Y_p0;
-      
-      // Z calibration //
-      Zbeam = Z_p0 + Z_p1*sqrt(Zbeam);
-
-      // Filling Beam tree
-      SofBeamID->SetZbeam(Zbeam);
-      SofBeamID->SetQmax(rand.Gaus(Qmax,0.15));
-      SofBeamID->SetAoQ(AoQ);
-      SofBeamID->SetAbeam(A);
-      SofBeamID->SetBeta(Beta);
-      SofBeamID->SetGamma(Gamma);
-      SofBeamID->SetBrho(Brho);
-      SofBeamID->SetXS2(XS2);
-      SofBeamID->SetXCC(XCC);
-      SofBeamID->SetYCC(YCC);
     }
   }
-}
 
-////////////////////////////////////////////////////////////////////////////////
-void Analysis::LoadCut(){
-  TString input_path = "./calibration/SofTrim/cut/";
-
-  TString rootfile;
-  TString cutfile;
-  TFile* file;
-  for(int i=0; i<3; i++){
-    // Q=77
-    rootfile = Form("cutsec%iQ77.root",i+1);
-    cutfile = input_path + rootfile;
-    file = new TFile(cutfile);
-    cutQ77[i] = (TCutG*) file->Get(Form("cutsec%iQ77",i+1));
-    // Q=78
-    rootfile = Form("cutsec%iQ78.root",i+1);
-    cutfile = input_path + rootfile;
-    file = new TFile(cutfile);
-    cutQ78[i] = (TCutG*) file->Get(Form("cutsec%iQ78",i+1));
-    // Q=79
-    rootfile = Form("cutsec%iQ79.root",i+1);
-    cutfile = input_path + rootfile;
-    file = new TFile(cutfile);
-    cutQ79[i] = (TCutG*) file->Get(Form("cutsec%iQ79",i+1));
-    // Q=80
-    rootfile = Form("cutsec%iQ80.root",i+1);
-    cutfile = input_path + rootfile;
-    file = new TFile(cutfile);
-    cutQ80[i] = (TCutG*) file->Get(Form("cutsec%iQ80",i+1));
+  ////////////////////////////////////////////////////////////////////////////////
+  void Analysis::LoadCut(){
+    TString input_path = "./calibration/SofTrim/cut/";
+
+    TString rootfile;
+    TString cutfile;
+    TFile* file;
+    for(int i=0; i<3; i++){
+      // Q=77
+      rootfile = Form("cutsec%iQ77.root",i+1);
+      cutfile = input_path + rootfile;
+      file = new TFile(cutfile);
+      cutQ77[i] = (TCutG*) file->Get(Form("cutsec%iQ77",i+1));
+      // Q=78
+      rootfile = Form("cutsec%iQ78.root",i+1);
+      cutfile = input_path + rootfile;
+      file = new TFile(cutfile);
+      cutQ78[i] = (TCutG*) file->Get(Form("cutsec%iQ78",i+1));
+      // Q=79
+      rootfile = Form("cutsec%iQ79.root",i+1);
+      cutfile = input_path + rootfile;
+      file = new TFile(cutfile);
+      cutQ79[i] = (TCutG*) file->Get(Form("cutsec%iQ79",i+1));
+      // Q=80
+      rootfile = Form("cutsec%iQ80.root",i+1);
+      cutfile = input_path + rootfile;
+      file = new TFile(cutfile);
+      cutQ80[i] = (TCutG*) file->Get(Form("cutsec%iQ80",i+1));
+    }
   }
-}
 
 
-////////////////////////////////////////////////////////////////////////////////
-void Analysis::LoadSpline(){
-  TString input_path = "./calibration/SofTwim/spline/";
-
-  TString rootfile = input_path + "spline_beta.root";
-  TFile* ifile = new TFile(rootfile,"read");
-
-  TString splinename;
-  if(ifile->IsOpen()){
-    cout << "Loading Beta spline for fission fragment analysis..." << endl;
-    for(int i=0; i<4; i++){
-      splinename = Form("spline_beta_sec%i",i+1);
-      fcorr_z_beta[i] = (TSpline3*) ifile->FindObjectAny(splinename);
+  ////////////////////////////////////////////////////////////////////////////////
+  void Analysis::LoadSpline(){
+    TString input_path = "./calibration/SofTwim/spline/";
+
+    TString rootfile = input_path + "spline_beta.root";
+    TFile* ifile = new TFile(rootfile,"read");
+
+    TString splinename;
+    if(ifile->IsOpen()){
+      cout << "Loading Beta spline for fission fragment analysis..." << endl;
+      for(int i=0; i<4; i++){
+        splinename = Form("spline_beta_sec%i",i+1);
+        fcorr_z_beta[i] = (TSpline3*) ifile->FindObjectAny(splinename);
+      }
+      ifile->Close();
     }
-    ifile->Close();
-  }
-  else
-    cout << "File " << rootfile << " not found!" << endl;
-
-  //*** ***//
-  rootfile = input_path + "spline_dt.root";
-  ifile = new TFile(rootfile,"read");
-
-  if(ifile->IsOpen()){
-    cout << "Loading DT spline for fission fragment analysis..." << endl;
-    for(int i=0; i<4; i++){
-      splinename = Form("spline_dt_sec%i",i+1);
-      fcorr_z_dt[i] = (TSpline3*) ifile->FindObjectAny(splinename);
+    else
+      cout << "File " << rootfile << " not found!" << endl;
+
+    //*** ***//
+    rootfile = input_path + "spline_dt.root";
+    ifile = new TFile(rootfile,"read");
+
+    if(ifile->IsOpen()){
+      cout << "Loading DT spline for fission fragment analysis..." << endl;
+      for(int i=0; i<4; i++){
+        splinename = Form("spline_dt_sec%i",i+1);
+        fcorr_z_dt[i] = (TSpline3*) ifile->FindObjectAny(splinename);
+      }
+      ifile->Close();
     }
-    ifile->Close();
-  }
-  else
-    cout << "File " << rootfile << " not found!" << endl;
+    else
+      cout << "File " << rootfile << " not found!" << endl;
 
-}
+  }
 
-////////////////////////////////////////////////////////////////////////////////
-int Analysis::DetermineQmax(){
-  int Qmax;
-  int Qsec[3];
-
-  unsigned int size = SofTrim->EnergySection.size();
-  for(unsigned int i=0; i<size; i++){
-    int SecNbr   = SofTrim->SectionNbr[i];
-    double Theta = SofTrim->Theta[i];
-    double Esec  = SofTrim->EnergySection[i];
-
-
-    if(cutQ77[SecNbr-1]->IsInside(Theta,Esec))
-      Qsec[SecNbr-1] = 77;
-    else if(cutQ78[SecNbr-1]->IsInside(Theta,Esec))
-      Qsec[SecNbr-1] = 78;
-    else if(cutQ79[SecNbr-1]->IsInside(Theta,Esec))
-      Qsec[SecNbr-1] = 79;
-    else if(cutQ80[SecNbr-1]->IsInside(Theta,Esec))
-      Qsec[SecNbr-1] = 80;
-    //else if(cutQ81[SecNbr-1]->IsInside(Theta,Esec))
+  ////////////////////////////////////////////////////////////////////////////////
+  int Analysis::DetermineQmax(){
+    int Qmax;
+    int Qsec[3];
+
+    unsigned int size = SofTrim->EnergySection.size();
+    for(unsigned int i=0; i<size; i++){
+      int SecNbr   = SofTrim->SectionNbr[i];
+      double Theta = SofTrim->Theta[i];
+      double Esec  = SofTrim->EnergySection[i];
+
+
+      if(cutQ77[SecNbr-1]->IsInside(Theta,Esec))
+        Qsec[SecNbr-1] = 77;
+      else if(cutQ78[SecNbr-1]->IsInside(Theta,Esec))
+        Qsec[SecNbr-1] = 78;
+      else if(cutQ79[SecNbr-1]->IsInside(Theta,Esec))
+        Qsec[SecNbr-1] = 79;
+      else if(cutQ80[SecNbr-1]->IsInside(Theta,Esec))
+        Qsec[SecNbr-1] = 80;
+      //else if(cutQ81[SecNbr-1]->IsInside(Theta,Esec))
       //Qsec[SecNbr-1] = 81;
-  }
+    }
 
-  Qmax = max(Qsec[0],Qsec[1]);
-  Qmax = max(Qsec[2],Qmax);
+    Qmax = max(Qsec[0],Qsec[1]);
+    Qmax = max(Qsec[2],Qmax);
 
-  return Qmax;
-}
+    return Qmax;
+  }
 
-////////////////////////////////////////////////////////////////////////////////
-void Analysis::End(){
-}
+  ////////////////////////////////////////////////////////////////////////////////
+  void Analysis::End(){
+  }
 
-////////////////////////////////////////////////////////////////////////////////
-void Analysis::InitParameter(){
-  fLS2_0 = 135.614;
-  fDS2   = 8000;
-  fDCC   = -10000;
-  fK_LS2 = -30e-8;
-  
-
-  //fBrho0 = 14.1008; // 238U run 367
-  //fBrho0 = 12.8719; // 238U run 368
-  //fBrho0 = 12.3255; // 238U run 369
-  //fBrho0 = 12.3352; // 238U run 419
-  //fBrho0 = 10.9476; // 189Pb
-  //fBrho0 = 10.9558; // 184Hg
-  //fBrho0 = 10.8183; // 182Hg
-  //fBrho0 = 10.6814; // 180Hg
-  //fBrho0 = 10.8138; // 187Pb
-  //fBrho0 = 11.3418; // 216Th
-  //fBrho0 = 11.0864; // 204Fr
-  fBrho0 = 11.2712; // 207Fr
-  //fBrho0 = 10.6814; // 175Pt
-  //fBrho0 = 11.5067; // 221Pa
-  //fBrho0 = 11.0955; // 199At run 423 & 424
-  //fBrho0 = 10.9970; // 199At run 425 & 426
-  //fBrho0 = 10.8697; //197At
-}
+  ////////////////////////////////////////////////////////////////////////////////
+  void Analysis::InitParameter(){
+    fLS2_0 = 135.614;
+    fDS2   = 8000;
+    fDCC   = -10000;
+    fK_LS2 = -30e-8;
+
+
+    //fBrho0 = 14.1008; // 238U run 367
+    //fBrho0 = 12.8719; // 238U run 368
+    //fBrho0 = 12.3255; // 238U run 369
+    //fBrho0 = 12.3352; // 238U run 422
+    //fBrho0 = 10.9476; // 189Pb
+    //fBrho0 = 10.9558; // 184Hg
+    //fBrho0 = 10.8183; // 182Hg
+    //fBrho0 = 10.6814; // 180Hg
+    //fBrho0 = 10.8138; // 187Pb
+    //fBrho0 = 11.3418; // 216Th
+    fBrho0 = 11.0864; // 204Fr
+    //fBrho0 = 11.2712; // 207Fr
+    //fBrho0 = 10.6814; // 175Pt
+    //fBrho0 = 11.5067; // 221Pa
+    //fBrho0 = 11.0955; // 199At run 423 & 424
+    //fBrho0 = 10.9970; // 199At run 425 & 426
+    //fBrho0 = 10.8697; //197At
+  }
 
-////////////////////////////////////////////////////////////////////////////////
-void Analysis::ReInitValue(){
-  SofBeamID->Clear();
-  SofFF->Clear();
-}
-////////////////////////////////////////////////////////////////////////////////
-void Analysis::InitOutputBranch(){
-  //RootOutput::getInstance()->GetTree()->Branch("Zbeam",&Zbeam,"Zbeam/D");
-  RootOutput::getInstance()->GetTree()->Branch("SofBeamID","TSofBeamID",&SofBeamID);
-  RootOutput::getInstance()->GetTree()->Branch("SofFissionFragment","TSofFissionFragment",&SofFF);
+  ////////////////////////////////////////////////////////////////////////////////
+  void Analysis::ReInitValue(){
+    SofBeamID->Clear();
+    SofFF->Clear();
+  }
+  ////////////////////////////////////////////////////////////////////////////////
+  void Analysis::InitOutputBranch(){
+    //RootOutput::getInstance()->GetTree()->Branch("Zbeam",&Zbeam,"Zbeam/D");
+    RootOutput::getInstance()->GetTree()->Branch("SofBeamID","TSofBeamID",&SofBeamID);
+    RootOutput::getInstance()->GetTree()->Branch("SofFissionFragment","TSofFissionFragment",&SofFF);
 
-}
-////////////////////////////////////////////////////////////////////////////////
-//            Construct Method to be pass to the DetectorFactory              //
-////////////////////////////////////////////////////////////////////////////////
-NPL::VAnalysis* Analysis::Construct(){
-  return (NPL::VAnalysis*) new Analysis();
-}
+  }
+  ////////////////////////////////////////////////////////////////////////////////
+  //            Construct Method to be pass to the DetectorFactory              //
+  ////////////////////////////////////////////////////////////////////////////////
+  NPL::VAnalysis* Analysis::Construct(){
+    return (NPL::VAnalysis*) new Analysis();
+  }
 
-////////////////////////////////////////////////////////////////////////////////
-//            Registering the construct method to the factory                 //
-////////////////////////////////////////////////////////////////////////////////
-extern "C"{
-  class proxy{
-    public:
-      proxy(){
-        NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
-      }
-  };
+  ////////////////////////////////////////////////////////////////////////////////
+  //            Registering the construct method to the factory                 //
+  ////////////////////////////////////////////////////////////////////////////////
+  extern "C"{
+    class proxy{
+      public:
+        proxy(){
+          NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
+        }
+    };
 
-  proxy p;
-}
+    proxy p;
+  }