diff --git a/Projects/e870/Analysis.cxx b/Projects/e870/Analysis.cxx
index af55fe7ba7ff505fa1267fe2f1d7178ab39cf748..928df1b414da66b2a7ad63bc850b23242be8c63d 100644
--- a/Projects/e870/Analysis.cxx
+++ b/Projects/e870/Analysis.cxx
@@ -21,7 +21,7 @@
  *****************************************************************************/
 #include <iostream>
 using namespace std;
-
+#include "NPVAnalysis.h"
 #include "Analysis.h"
 #include "NPAnalysisFactory.h"
 #include "NPDetectorManager.h"
@@ -91,9 +91,6 @@ void Analysis::Init() {
   CsI_E_M2 = 0;
   Energy = 0;
   ThetaGDSurface = 0;
-  X = 0;
-  Y = 0;
-  Z = 0;
   dE = 0;
   BeamDirection = TVector3(0, 0, 1);
   nbTrack = 0;
@@ -128,19 +125,23 @@ void Analysis::TreatEvent() {
     int TelescopeNumber = M2->TelescopeNumber[countMust2];
 
     /************************************************/
+
     // Part 1 : Impact Angle
+    //THIS IS THE NEW PART
+
     ThetaM2Surface = 0;
     ThetaNormalTarget = 0;
     TVector3 HitDirection = M2->GetPositionOfInteraction(countMust2) - BeamImpact;
+    TVector3 Coords = M2->GetPositionOfInteraction(countMust2);
     double Theta = HitDirection.Angle(BeamDirection);
-
-    X = M2->GetPositionOfInteraction(countMust2).X();
-    Y = M2->GetPositionOfInteraction(countMust2).Y();
-    Z = M2->GetPositionOfInteraction(countMust2).Z();
+    X.push_back(Coords.x());
+    Y.push_back(Coords.y());
+    Z.push_back(Coords.z());
 
     ThetaM2Surface = HitDirection.Angle(-M2->GetTelescopeNormal(countMust2));
     ThetaNormalTarget = HitDirection.Angle(TVector3(0, 0, 1));
 
+
     /************************************************/
 
     /************************************************/
@@ -167,9 +168,11 @@ void Analysis::TreatEvent() {
     Energy = LightTarget.EvaluateInitialEnergy(Energy, TargetThickness * 0.5, Theta);
 
     // What is written in the tree
+ 
     ThetaLab.push_back(Theta / deg);
     Ex.push_back(reaction.ReconstructRelativistic(Energy, Theta));
     ELab.push_back(Energy);
+    
     /************************************************/
 
   } // end loop MUST2
@@ -185,9 +188,12 @@ void Analysis::InitOutputBranch() {
   RootOutput::getInstance()->GetTree()->Branch("ThetaLab", &ThetaLab);
   RootOutput::getInstance()->GetTree()->Branch("ThetaCM", &ThetaCM, "ThetaCM/D");
   RootOutput::getInstance()->GetTree()->Branch("Run", &Run, "Run/I");
-  RootOutput::getInstance()->GetTree()->Branch("X", &X, "X/D");
-  RootOutput::getInstance()->GetTree()->Branch("Y", &Y, "Y/D");
-  RootOutput::getInstance()->GetTree()->Branch("Z", &Z, "Z/D");
+  // RootOutput::getInstance()->GetTree()->Branch("X", &X, "X/D"); Like before
+  // RootOutput::getInstance()->GetTree()->Branch("Y", &Y, "Y/D"); 
+  // RootOutput::getInstance()->GetTree()->Branch("Z", &Z, "Z/D");
+  RootOutput::getInstance()->GetTree()->Branch("X", &X);
+  RootOutput::getInstance()->GetTree()->Branch("Y", &Y);
+  RootOutput::getInstance()->GetTree()->Branch("Z", &Z);
   RootOutput::getInstance()->GetTree()->Branch("dE", &dE, "dE/D");
   if (!simulation) {
   }
@@ -222,9 +228,12 @@ void Analysis::ReInitValue() {
   EDC = -1000;
   BeamEnergy = -1000;
   ThetaCM = -1000;
-  X = -1000;
-  Y = -1000;
-  Z = -1000;
+  // X = -1000;
+  // Y = -1000;
+  // Z = -1000;
+  X.clear();
+  Y.clear();
+  Z.clear();
   dE = -1000;
   ELab.clear();
   ThetaLab.clear();
diff --git a/Projects/e870/Analysis.h b/Projects/e870/Analysis.h
index 753d9fcdbfb81f25e6efff8ba3727f67561884b3..6ba81025d7d41cb800feaf2604192a9a06b1152b 100644
--- a/Projects/e870/Analysis.h
+++ b/Projects/e870/Analysis.h
@@ -93,9 +93,9 @@ class Analysis: public NPL::VAnalysis{
   double BeamEnergy;
   
   double ThetaGDSurface ;
-  double X ;
-  double Y ;
-  double Z ;
+  std::vector<double> X ;
+  std::vector<double> Y ;
+  std::vector<double> Z ;
   // Vamos Branches
   unsigned long long int LTS;
   
diff --git a/Projects/e870/DetectorConfiguration/MUST2_E805.detector b/Projects/e870/DetectorConfiguration/MUST2_E805.detector
index 4c78f7563ad7c4d4783c632bee6c0a61c96cc24f..bc14fd00cb6f8c2b1fd2902b90aac644f6608fd6 100644
--- a/Projects/e870/DetectorConfiguration/MUST2_E805.detector
+++ b/Projects/e870/DetectorConfiguration/MUST2_E805.detector
@@ -58,3 +58,18 @@ M2Telescope
  VIS= all  
  CsIOffset= 1
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+CATSDetector
+ CATSNumber = 1
+ X1_Y1= 43.71  -41.76  -1587.1 mm
+ X28_Y1= -27.41  -41.76  -1587.1 mm
+ X1_Y28= 43.71  29.36  -1587.1 mm
+ X28_Y28= -27.41 29.36  -1587.1 mm
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+CATSDetector
+ CATSNumber = 2
+ X1_Y1= 33.94    -37.06   -1090.1 mm
+ X28_Y1= -37.18  -37.06   -1090.1 mm
+ X1_Y28= 33.94    34.06  -1090.1 mm 
+ X28_Y28= -37.18  34.06   -1090.1 mm
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
diff --git a/Projects/e870/Makefile b/Projects/e870/Makefile
index 719c7c7a4e2b71e65030d87930757eb33ae68a19..223d5e7e565eed065829403cf49c6917558f3848 100644
--- a/Projects/e870/Makefile
+++ b/Projects/e870/Makefile
@@ -1,44 +1,35 @@
 # CMAKE generated file: DO NOT EDIT!
-# Generated by "Unix Makefiles" Generator, CMake Version 3.26
+# Generated by "Unix Makefiles" Generator, CMake Version 3.16
 
 # Default target executed when no arguments are given to make.
 default_target: all
+
 .PHONY : default_target
 
 # Allow only one "make -f Makefile2" at a time, but pass parallelism.
 .NOTPARALLEL:
 
+
 #=============================================================================
 # Special targets provided by cmake.
 
 # Disable implicit rules so canonical targets will work.
 .SUFFIXES:
 
-# Disable VCS-based implicit rules.
-% : %,v
-
-# Disable VCS-based implicit rules.
-% : RCS/%
-
-# Disable VCS-based implicit rules.
-% : RCS/%,v
-
-# Disable VCS-based implicit rules.
-% : SCCS/s.%
 
-# Disable VCS-based implicit rules.
-% : s.%
+# Remove some rules from gmake that .SUFFIXES does not remove.
+SUFFIXES =
 
 .SUFFIXES: .hpux_make_needs_suffix_list
 
-# Command-line flag to silence nested $(MAKE).
-$(VERBOSE)MAKESILENT = -s
 
-#Suppress display of executed commands.
+# Suppress display of executed commands.
 $(VERBOSE).SILENT:
 
+
 # A target that is always out of date.
 cmake_force:
+
 .PHONY : cmake_force
 
 #=============================================================================
@@ -48,67 +39,70 @@ cmake_force:
 SHELL = /bin/sh
 
 # The CMake executable.
-CMAKE_COMMAND = /usr/local/Cellar/cmake/3.26.4/bin/cmake
+CMAKE_COMMAND = /usr/bin/cmake
 
 # The command to remove a file.
-RM = /usr/local/Cellar/cmake/3.26.4/bin/cmake -E rm -f
+RM = /usr/bin/cmake -E remove -f
 
 # Escaping for special characters.
 EQUALS = =
 
 # The top-level source directory on which CMake was run.
-CMAKE_SOURCE_DIR = /Users/valerian/Software/nptool_gitlab/nptool/Projects/e870
+CMAKE_SOURCE_DIR = /home/tania/Software/nptool/Projects/e870
 
 # The top-level build directory on which CMake was run.
-CMAKE_BINARY_DIR = /Users/valerian/Software/nptool_gitlab/nptool/Projects/e870
+CMAKE_BINARY_DIR = /home/tania/Software/nptool/Projects/e870
 
 #=============================================================================
 # Targets provided globally by CMake.
 
-# Special rule for the target edit_cache
-edit_cache:
-	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake cache editor..."
-	/usr/local/Cellar/cmake/3.26.4/bin/ccmake -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR)
-.PHONY : edit_cache
-
-# Special rule for the target edit_cache
-edit_cache/fast: edit_cache
-.PHONY : edit_cache/fast
-
 # Special rule for the target rebuild_cache
 rebuild_cache:
 	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake to regenerate build system..."
-	/usr/local/Cellar/cmake/3.26.4/bin/cmake --regenerate-during-build -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR)
+	/usr/bin/cmake -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR)
 .PHONY : rebuild_cache
 
 # Special rule for the target rebuild_cache
 rebuild_cache/fast: rebuild_cache
+
 .PHONY : rebuild_cache/fast
 
+# Special rule for the target edit_cache
+edit_cache:
+	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake cache editor..."
+	/usr/bin/ccmake -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR)
+.PHONY : edit_cache
+
+# Special rule for the target edit_cache
+edit_cache/fast: edit_cache
+
+.PHONY : edit_cache/fast
+
 # The main all target
 all: cmake_check_build_system
-	$(CMAKE_COMMAND) -E cmake_progress_start /Users/valerian/Software/nptool_gitlab/nptool/Projects/e870/CMakeFiles /Users/valerian/Software/nptool_gitlab/nptool/Projects/e870//CMakeFiles/progress.marks
-	$(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 all
-	$(CMAKE_COMMAND) -E cmake_progress_start /Users/valerian/Software/nptool_gitlab/nptool/Projects/e870/CMakeFiles 0
+	$(CMAKE_COMMAND) -E cmake_progress_start /home/tania/Software/nptool/Projects/e870/CMakeFiles /home/tania/Software/nptool/Projects/e870/CMakeFiles/progress.marks
+	$(MAKE) -f CMakeFiles/Makefile2 all
+	$(CMAKE_COMMAND) -E cmake_progress_start /home/tania/Software/nptool/Projects/e870/CMakeFiles 0
 .PHONY : all
 
 # The main clean target
 clean:
-	$(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 clean
+	$(MAKE) -f CMakeFiles/Makefile2 clean
 .PHONY : clean
 
 # The main clean target
 clean/fast: clean
+
 .PHONY : clean/fast
 
 # Prepare targets for installation.
 preinstall: all
-	$(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 preinstall
+	$(MAKE) -f CMakeFiles/Makefile2 preinstall
 .PHONY : preinstall
 
 # Prepare targets for installation.
 preinstall/fast:
-	$(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 preinstall
+	$(MAKE) -f CMakeFiles/Makefile2 preinstall
 .PHONY : preinstall/fast
 
 # clear depends
@@ -121,36 +115,39 @@ depend:
 
 # Build rule for target.
 NPAnalysis: cmake_check_build_system
-	$(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 NPAnalysis
+	$(MAKE) -f CMakeFiles/Makefile2 NPAnalysis
 .PHONY : NPAnalysis
 
 # fast build rule for target.
 NPAnalysis/fast:
-	$(MAKE) $(MAKESILENT) -f CMakeFiles/NPAnalysis.dir/build.make CMakeFiles/NPAnalysis.dir/build
+	$(MAKE) -f CMakeFiles/NPAnalysis.dir/build.make CMakeFiles/NPAnalysis.dir/build
 .PHONY : NPAnalysis/fast
 
 Analysis.o: Analysis.cxx.o
+
 .PHONY : Analysis.o
 
 # target to build an object file
 Analysis.cxx.o:
-	$(MAKE) $(MAKESILENT) -f CMakeFiles/NPAnalysis.dir/build.make CMakeFiles/NPAnalysis.dir/Analysis.cxx.o
+	$(MAKE) -f CMakeFiles/NPAnalysis.dir/build.make CMakeFiles/NPAnalysis.dir/Analysis.cxx.o
 .PHONY : Analysis.cxx.o
 
 Analysis.i: Analysis.cxx.i
+
 .PHONY : Analysis.i
 
 # target to preprocess a source file
 Analysis.cxx.i:
-	$(MAKE) $(MAKESILENT) -f CMakeFiles/NPAnalysis.dir/build.make CMakeFiles/NPAnalysis.dir/Analysis.cxx.i
+	$(MAKE) -f CMakeFiles/NPAnalysis.dir/build.make CMakeFiles/NPAnalysis.dir/Analysis.cxx.i
 .PHONY : Analysis.cxx.i
 
 Analysis.s: Analysis.cxx.s
+
 .PHONY : Analysis.s
 
 # target to generate assembly for a file
 Analysis.cxx.s:
-	$(MAKE) $(MAKESILENT) -f CMakeFiles/NPAnalysis.dir/build.make CMakeFiles/NPAnalysis.dir/Analysis.cxx.s
+	$(MAKE) -f CMakeFiles/NPAnalysis.dir/build.make CMakeFiles/NPAnalysis.dir/Analysis.cxx.s
 .PHONY : Analysis.cxx.s
 
 # Help Target
@@ -159,8 +156,8 @@ help:
 	@echo "... all (the default if no target is provided)"
 	@echo "... clean"
 	@echo "... depend"
-	@echo "... edit_cache"
 	@echo "... rebuild_cache"
+	@echo "... edit_cache"
 	@echo "... NPAnalysis"
 	@echo "... Analysis.o"
 	@echo "... Analysis.i"
diff --git a/Projects/e870/macros/DrawExInvariant.cxx b/Projects/e870/macros/DrawExInvariant.cxx
index 0d191812285a13367b9dcfdda01464145e69b49d..0d2aa1f1147f99ba71d06917362d50913d612e4e 100644
--- a/Projects/e870/macros/DrawExInvariant.cxx
+++ b/Projects/e870/macros/DrawExInvariant.cxx
@@ -1,61 +1,66 @@
-void DrawDalitz() {
+void DrawExInvariant() {
   ////////////////////////////////////////////////////////////////////////////////
-  TLorentzVector LV_alpha, LV_proton, LV_fragment;
+  TLorentzVector LV_alpha, LV_7Li;
+  // TLorentzVector LV_alpha, LV_proton, LV_fragment;
   double mu = 931.5;                   // MeV/c2
-  double m_proton = mu + 7.289;        // Mev/c2
   double m_alpha = 4 * mu + 2.425;     // Mev/c2
-  double m_fragment = 44 * mu - 3.755; // Mev/c2
-  double m_5Li = 5 * mu + 11.68;       // MeV/c2
+  double m_7Li = 7 * mu + 14.91;       // MeV/c2
+  // double m_proton = mu + 7.289;        // Mev/c2
+  // double m_fragment = 44 * mu - 3.755; // Mev/c2
+  // double m_5Li = 5 * mu + 11.68;       // MeV/c2
 
-  TChain* t = new TChain("PhysicsTree"); // SimulatedTree
-  t->Add("NPRootA/NPr0371_*a.root");     // Name of simulated file
+  TChain* t = new TChain("PhysicsTree");
+  t->Add("../../Outputs/Analysis/testnew.root");     // Name of simulated file
 
   // YOUR ALPHA AND Li CUT
-  TFile* fcuta = new TFile("CUT_e805/CUT_alpha_must.root", "read");
-  TCutG* cuta = (TCutG*)gDirectory->FindObjectAny("CUT_alpha_must");
-  TFile* fcutp = new TFile("CUT_e805/CUT_proton_must.root", "read");
-  TCutG* cutp = (TCutG*)gDirectory->FindObjectAny("CUT_proton_must");
+  TFile* fcuta = new TFile("cuts/CUT_alpha.root", "read");
+  TCutG* cuta = (TCutG*)gDirectory->FindObjectAny("CUT_alpha");
+  TFile* fcutli = new TFile("cuts/CUT_7Li.root", "read");
+  TCutG* cutli = (TCutG*)gDirectory->FindObjectAny("CUT_7Li");
 
   ////////////////////////////////////////////////////////////////////////////////
   TMust2Physics* M2 = new TMust2Physics();
   std::vector<double>* ELab = 0;
   std::vector<double>* ThetaLab = 0;
-  std::vector<double>* M2_X = 0;
-  std::vector<double>* M2_Y = 0;
-  std::vector<double>* M2_Z = 0;
+  std::vector<double>* X = 0;
+  std::vector<double>* Y = 0;
+  std::vector<double>* Z = 0;
+  std::vector<double>* Ex = 0;
   TBranch* b_ELab;
   TBranch* b_ThetaLab;
-  TBranch* b_M2_X;
-  TBranch* b_M2_Y;
-  TBranch* b_M2_Z;
+  TBranch* b_X;
+  TBranch* b_Y;
+  TBranch* b_Z;
+  TBranch* b_Ex;
 
-  float CATS1_X, CATS1_Y, CATS2_X, CATS2_Y, Xf, Yf;
+  // float CATS1_X, CATS1_Y, CATS2_X, CATS2_Y, Xf, Yf;
 
   t->SetBranchStatus("*", false);
   t->SetBranchStatus("MUST2", true);
   t->SetBranchAddress("MUST2", &M2);
-  t->SetBranchStatus("M2_ELab", true);
-  t->SetBranchAddress("M2_ELab", &ELab, &b_ELab);
-  t->SetBranchStatus("M2_ThetaLab", true);
-  t->SetBranchAddress("M2_ThetaLab", &ThetaLab, &b_ThetaLab);
-  t->SetBranchStatus("CATS1_X", "true");
-  t->SetBranchAddress("CATS1_X", &CATS1_X);
-  t->SetBranchStatus("CATS1_Y", "true");
-  t->SetBranchAddress("CATS1_Y", &CATS1_Y);
-  t->SetBranchStatus("CATS2_X", "true");
-  t->SetBranchAddress("CATS2_X", &CATS2_X);
-  t->SetBranchStatus("CATS2_Y", "true");
-  t->SetBranchAddress("CATS2_Y", &CATS2_Y);
-  t->SetBranchStatus("M2_X", "true");
-  t->SetBranchAddress("M2_X", &M2_X, &b_M2_X);
-  t->SetBranchStatus("M2_Y", "true");
-  t->SetBranchAddress("M2_Y", &M2_Y, &b_M2_Y);
-  t->SetBranchStatus("M2_Z", "true");
-  t->SetBranchAddress("M2_Z", &M2_Z, &b_M2_Z);
+  t->SetBranchStatus("ELab", true);
+  t->SetBranchAddress("ELab", &ELab, &b_ELab);
+  t->SetBranchStatus("ThetaLab", true);
+  t->SetBranchAddress("ThetaLab", &ThetaLab, &b_ThetaLab);
+  // t->SetBranchStatus("CATS1_X", "true");
+  // t->SetBranchAddress("CATS1_X", &CATS1_X);
+  // t->SetBranchStatus("CATS1_Y", "true");
+  // t->SetBranchAddress("CATS1_Y", &CATS1_Y);
+  // t->SetBranchStatus("CATS2_X", "true");
+  // t->SetBranchAddress("CATS2_X", &CATS2_X);
+  // t->SetBranchStatus("CATS2_Y", "true");
+  // t->SetBranchAddress("CATS2_Y", &CATS2_Y);
+  t->SetBranchStatus("X", "true");
+  t->SetBranchAddress("X", &X, &b_X);
+  t->SetBranchStatus("Y", "true");
+  t->SetBranchAddress("Y", &Y, &b_Y);
+  t->SetBranchStatus("Z", "true");
+  t->SetBranchAddress("Z", &Z, &b_Z);
+  t->SetBranchStatus("Ex", "true");
+  t->SetBranchAddress("Ex", &Ex, &b_Ex);
 
   TH1F* hExTot = new TH1F("hExTot", "hExTot", 1000, -100, 100);
-  TH1F* hEx5Li = new TH1F("hEx5Li", "hEx5Li", 1000, -100, 100);
-
+  TH1F* hEx = new TH1F("hEx", "hEx", 1000, -100, 100);
   TH2F* hE1E2 = new TH2F("hE1E2", "hE1E2", 1000, 0, 500, 1000, 0, 500);
 
   ////////////////////////////////////////////////////////////////////////////////
@@ -68,61 +73,57 @@ void DrawDalitz() {
     double theta_alpha;
     double phi_alpha;
 
-    double e_proton;
-    double theta_proton;
-    double phi_proton;
-
-    double e_fragment;
-    double theta_fragment;
-    double phi_fragment;
+    double e_li;
+    double theta_li;
+    double phi_li;
 
     double xtarget = 0;
     double ytarget = 0;
-    double CATSX_diff = 0;
-    double CATSY_diff = 0;
-
-    if (CATS1_X > -1500 && CATS2_X > -1500 && Xf > -1500) {
-      CATSX_diff = -(CATS2_X - CATS1_X);
-      xtarget = -Xf;
-    }
-    if (CATS1_Y > -1500 && CATS2_X > -1500 && Yf > -1500) {
-      ytarget = Yf;
-      CATSY_diff = CATS2_X - CATS1_X;
-    }
+    // double CATSX_diff = 0;
+    // double CATSY_diff = 0;
+
+    // if (CATS1_X > -1500 && CATS2_X > -1500 && Xf > -1500) {
+    //   CATSX_diff = -(CATS2_X - CATS1_X);
+    //   xtarget = -Xf;
+    // }
+    // if (CATS1_Y > -1500 && CATS2_X > -1500 && Yf > -1500) {
+    //   ytarget = Yf;
+    //   CATSY_diff = CATS2_X - CATS1_X;
+    // }
 
     TVector3 BeamImpact(xtarget, ytarget, 0);
-    TVector3 BeamDirection(CATSX_diff, CATSY_diff, 497.1);
+    TVector3 BeamDirection(0,0,1);  
 
     if (M2->Si_E.size() == 2) {
-      if ((cuta->IsInside(M2->CsI_E[0], M2->Si_E[0]) && cutp->IsInside(M2->CsI_E[1], M2->Si_E[1]))) {
+      if ((cuta->IsInside(M2->CsI_E[0], M2->Si_E[0]) && cutli->IsInside(M2->CsI_E[1], M2->Si_E[1]))) {
         // Particle1 = alpha
-        TVector3 PositionOfInteraction_alpha(M2_X->at(0), M2_Y->at(0), M2_Z->at(0));
+        TVector3 PositionOfInteraction_alpha(X->at(0), Y->at(0), Z->at(0));
         TVector3 HitDirection_alpha = PositionOfInteraction_alpha - BeamImpact;
         e_alpha = ELab->at(0);
         theta_alpha = HitDirection_alpha.Angle(BeamDirection);
         phi_alpha = HitDirection_alpha.Phi();
 
-        // Particle2 = proton
-        TVector3 PositionOfInteraction_proton(M2_X->at(1), M2_Y->at(1), M2_Z->at(1));
-        TVector3 HitDirection_proton = PositionOfInteraction_proton - BeamImpact;
-        e_proton = ELab->at(1);
-        theta_proton = HitDirection_proton.Angle(BeamDirection);
-        phi_proton = HitDirection_proton.Phi();
+        // Particle2 = 7li
+        TVector3 PositionOfInteraction_li(X->at(1), Y->at(1), Z->at(1));
+        TVector3 HitDirection_li = PositionOfInteraction_li - BeamImpact;
+        e_li = ELab->at(1);
+        theta_li = HitDirection_li.Angle(BeamDirection);
+        phi_li = HitDirection_li.Phi();
       }
-      else if (cuta->IsInside(M2->CsI_E[1], M2->Si_E[1]) && cutp->IsInside(M2->CsI_E[0], M2->Si_E[0])) {
+      else if (cuta->IsInside(M2->CsI_E[1], M2->Si_E[1]) && cutli->IsInside(M2->CsI_E[0], M2->Si_E[0])) {
         // Particle2 = alpha
-        TVector3 PositionOfInteraction_alpha(M2_X->at(1), M2_Y->at(1), M2_Z->at(1));
+        TVector3 PositionOfInteraction_alpha(X->at(1), Y->at(1), Z->at(1));
         TVector3 HitDirection_alpha = PositionOfInteraction_alpha - BeamImpact;
         e_alpha = ELab->at(1);
         theta_alpha = HitDirection_alpha.Angle(BeamDirection);
         phi_alpha = HitDirection_alpha.Phi();
 
-        // Particle1 = proton
-        TVector3 PositionOfInteraction_proton(M2_X->at(0), M2_Y->at(0), M2_Z->at(0));
-        TVector3 HitDirection_proton = PositionOfInteraction_proton - BeamImpact;
-        e_proton = ELab->at(0);
-        theta_proton = HitDirection_proton.Angle(BeamDirection);
-        phi_proton = HitDirection_proton.Phi();
+        // Particle1 = 7li
+        TVector3 PositionOfInteraction_li(X->at(0), Y->at(0), Z->at(0));
+        TVector3 HitDirection_li = PositionOfInteraction_li - BeamImpact;
+        e_li = ELab->at(0);
+        theta_li = HitDirection_li.Angle(BeamDirection);
+        phi_li = HitDirection_li.Phi();
       }
 
       ////////////////////////////////////////////////////////////////////////////////
@@ -138,29 +139,48 @@ void DrawDalitz() {
       LV_alpha.SetPxPyPzE(p_alpha.x(), p_alpha.y(), p_alpha.z(), etot_alpha);
 
       ////////////////////////////////////////////////////////////////////////////////
-      double gamma_proton = e_proton / m_proton + 1;
-      double beta_proton = sqrt(1 - 1 / (pow(gamma_proton, 2.)));
-      double p_proton_mag = gamma_proton * m_proton * beta_proton;
-      double etot_proton = sqrt(pow(p_proton_mag, 2) + pow(m_proton, 2));
-
-      TVector3 p_proton(sin(theta_proton) * cos(phi_proton), // px
-                        sin(theta_proton) * sin(phi_proton), // py
-                        cos(theta_proton));                  // pz
-      p_proton.SetMag(p_proton_mag);
-      LV_proton.SetPxPyPzE(p_proton.x(), p_proton.y(), p_proton.z(), etot_proton);
-      TLorentzVector LV_total = LV_alpha + LV_proton;
-
-      double ExTot = LV_total.Mag() - m_alpha - m_proton;
+      double gamma_li = e_li / m_7Li + 1;
+      double beta_li = sqrt(1 - 1 / (pow(gamma_li, 2.)));
+      double p_li_mag = gamma_li * m_7Li * beta_li;
+      double etot_li = sqrt(pow(p_li_mag, 2) + pow(m_7Li, 2));
+
+      TVector3 p_li(sin(theta_li) * cos(phi_li), // px
+                        sin(theta_li) * sin(phi_li), // py
+                        cos(theta_li));                  // pz
+      p_li.SetMag(p_li_mag);
+      LV_7Li.SetPxPyPzE(p_li.x(), p_li.y(), p_li.z(), etot_li);
+      TLorentzVector LV_total = LV_alpha + LV_7Li;
+
+      double ExTot = LV_total.Mag() - m_alpha - m_7Li;
       hExTot->Fill(ExTot);
-      if (Ex.size() >= 2) {
-        hEx->Fill(Ex.at(0));
-        hEx->Fill(Ex.at(1));
+      if (Ex->size() >= 2) {  
+        hEx->Fill(Ex->at(0));
+        hEx->Fill(Ex->at(1));
       }
-      hE1E2->Fill(e_alpha, e_proton);
+      hE1E2->Fill(e_alpha, e_li);
     }
   }
   TCanvas* c1 = new TCanvas();
+  gStyle->SetOptStat(0);
   hExTot->Draw();
+  hEx->SetLineColor(kRed);
   hEx->Draw("same");
+
+  TF1* fitExTot = new TF1("fitExTot", "gaus", hExTot->GetXaxis()->GetXmin(), hExTot->GetXaxis()->GetXmax());
+  hExTot->Fit(fitExTot, "Q");
+  double meanExTot = fitExTot->GetParameter(1);
+  std::cout << "Center of hExTot: " << meanExTot << " MeV" << std::endl; 
+
+  double fitRangeMin = 0.0; 
+  double fitRangeMax = 15.0; 
+  TF1* fitEx = new TF1("fitEx", "gaus", fitRangeMin, fitRangeMax);
+  hEx->Fit(fitEx);
+  double meanEx = fitEx->GetParameter(1);
+  std::cout << "Center of hEx: " << meanEx << " MeV" << std::endl;
+
+  TLegend* legend = new TLegend(0.75, 0.75, 0.85, 0.85);
+  legend->AddEntry(hExTot, "hExTot", "l"); // "l" for line
+  legend->AddEntry(hEx, "hex", "l");
+  legend->Draw();
 }
-~
+
diff --git a/Projects/e870/reaction/10Bepalpha.reaction b/Projects/e870/reaction/10Bepalpha.reaction
index a4b5f079592b21fa8e3deaba9be601711f87b4f3..3945abf5819ef29662dfec91a4a517eb9ac3a662 100644
--- a/Projects/e870/reaction/10Bepalpha.reaction
+++ b/Projects/e870/reaction/10Bepalpha.reaction
@@ -13,9 +13,9 @@ Beam
   SigmaX= 6.216 mm
   SigmaY= 6.109 mm
   % Playground:
-  %SigmaEnergy= ?? MeV
-  %SigmaX= ?? mm
-  %SigmaY= ?? mm
+  %SigmaEnergy= 23.2 MeV
+  %SigmaX= 7 mm
+  %SigmaY= 6.8 mm
 
   MeanThetaX= 0 deg
   MeanPhiY= 0 deg