diff --git a/Projects/e793s/Analysis.cxx b/Projects/e793s/Analysis.cxx
index b2d69a906e58b281938f0d3fd3741316d2344913..630354c87dcbab42abbd472cb09467fb83a441fc 100755
--- a/Projects/e793s/Analysis.cxx
+++ b/Projects/e793s/Analysis.cxx
@@ -43,7 +43,6 @@ Analysis::~Analysis(){
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::Init() {
   ///////////////////////////////////////////////////////////////////////////////  
- 
 
   agata_zShift=51*mm;
   //BrhoRef=0.65;
@@ -60,6 +59,17 @@ void Analysis::Init() {
   // get reaction information
   reaction.ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
   OriginalBeamEnergy = reaction.GetBeamEnergy();
+  reaction.Print(); //TESTING PARSER
+
+  // get beam position from .reaction file
+  Beam = (NPL::Beam*) reaction.GetParticle1(); 
+  XBeam = Beam->GetMeanX();
+  YBeam = Beam->GetMeanY();
+
+  cout << " ---------------- Beam Position ---------------- " << endl;
+  cout << " \tX = " << XBeam << " mm\tY = " << YBeam  << " mm" << endl;
+  cout << " ----------------------------------------------- " << endl;
+
   // target thickness
   TargetThickness = m_DetectorManager->GetTargetThickness();
   string TargetMaterial = m_DetectorManager->GetTargetMaterial();
@@ -128,13 +138,8 @@ void Analysis::TreatEvent(){
     }
   }
 
-  double xbeam = 0.0; 
-  double ybeam = 0.0;
-
-  // double XTarget = CATS->GetPositionOnTarget().X();
-  // double YTarget = CATS->GetPositionOnTarget().Y();
-  TVector3 BeamDirection(xbeam,ybeam,1);
-  BeamImpact = TVector3(xbeam,ybeam,m_DetectorManager->GetTargetZ()); 
+  TVector3 BeamDirection(XBeam,YBeam,1);
+  BeamImpact = TVector3(XBeam,YBeam,m_DetectorManager->GetTargetZ()); 
 
   ParticleMult=M2->Si_E.size()+MG->DSSD_E.size();
   //ParticleMult=M2->Si_E.size();
@@ -167,7 +172,7 @@ void Analysis::TreatEvent(){
     Z.push_back(M2->GetPositionOfInteraction(countMust2).Z());
 
     ThetaM2Surface = HitDirection.Angle(- M2->GetTelescopeNormal(countMust2) );
-    ThetaNormalTarget = HitDirection.Angle( TVector3(xbeam,ybeam,1) ) ;
+    ThetaNormalTarget = HitDirection.Angle( TVector3(XBeam,YBeam,1) ) ;
 
    // cout<<"Must2 Znormal:"<<M2 -> GetTelescopeNormal(countMust2).Z()<<endl;
   //  cout<<"Must2 telescope:"<<M2->TelescopeNumber[countMust2]<<endl;
@@ -231,8 +236,6 @@ void Analysis::TreatEvent(){
     thetalab_tmp = 0;
     philab_tmp = 0;
     TVector3 HitDirection = MG->GetPositionOfInteraction(countMugast) - BeamImpact;
-//    cout << HitDirection.X() << "\t" << HitDirection.Y() << "\t" << HitDirection.Z() << "\t\t" << MG->GetPositionOfInteraction(countMugast).X() << "\t" << MG->GetPositionOfInteraction(countMugast).Y() << "\t" << MG->GetPositionOfInteraction(countMugast).Z() << "\t\t" << BeamImpact.X() << "\t" << BeamImpact.Y() << "\t" << BeamImpact.Z() << endl;
-    //HitDirection = [Xp, Yp, Zp] - [Xb, Yb, Zb] = [Xd, Td, Zd] IN MY CODE!
     thetalab_tmp = HitDirection.Angle(BeamDirection);
     philab_tmp = HitDirection.Phi();
 
@@ -242,7 +245,7 @@ void Analysis::TreatEvent(){
 
     //ThetaMGSurface = HitDirection.Angle( TVector3(0,0,1) ) ;
     ThetaMGSurface = HitDirection.Angle( MG -> GetTelescopeNormal(countMugast) );
-    ThetaNormalTarget = HitDirection.Angle( TVector3(xbeam,ybeam,-1) ) ;
+    ThetaNormalTarget = HitDirection.Angle( TVector3(XBeam,YBeam,-1) ) ;
 
     // Part 2 : Impact Energy
     Energy = elab_tmp = 0;
@@ -265,27 +268,6 @@ void Analysis::TreatEvent(){
      
     PhiLab.push_back(philab_tmp/deg);
 
-
-
-//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-/*
-cout << MG->TelescopeNumber[countMugast] << " "
-     << MG->GetPositionOfInteraction(countMugast).X() << " " 
-     << MG->GetPositionOfInteraction(countMugast).Y() << " " 
-     << MG->GetPositionOfInteraction(countMugast).Z() << " " 
-     << Energy <<  " " 
-     << endl; 
-*/
-
-///////////!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
-
-
-
-
-
     if(sizeMG==1){
       MG_T = MG->DSSD_T[0];
       MG_E = MG->DSSD_E[0];
@@ -602,7 +584,7 @@ cout << MG->TelescopeNumber[countMugast] << " "
     AddBack_EDC.clear();
     EAgata = -1000;
     ELab.clear();
-    RawEnergy.clear(); //CPx ADDITION
+    RawEnergy.clear(); 
     ThetaLab.clear();
     PhiLab.clear();
     ThetaCM.clear();
diff --git a/Projects/e793s/Analysis.h b/Projects/e793s/Analysis.h
index 3dfb5fcf308dd31110300d5987a2b1d199255437..3d3ed1bbc99fc02287a55c988593b85b56f06f0a 100755
--- a/Projects/e793s/Analysis.h
+++ b/Projects/e793s/Analysis.h
@@ -25,6 +25,7 @@
 #include"NPVAnalysis.h"
 #include"NPEnergyLoss.h"
 #include"NPReaction.h"
+#include"NPBeam.h"
 #include"RootOutput.h"
 #include"RootInput.h"
 #include "TMust2Physics.h"
@@ -64,7 +65,7 @@ class Analysis: public NPL::VAnalysis{
     std::vector<double> ELab;
     std::vector<double> Ex;
     std::vector<double> Ecm;
-    std::vector<double> RawEnergy; //CPx ADDITION
+    std::vector<double> RawEnergy;
     std::vector<double> ThetaLab;
     std::vector<double> PhiLab;
     std::vector<double> ThetaCM;
@@ -82,6 +83,10 @@ class Analysis: public NPL::VAnalysis{
     double OriginalBeamEnergy ; // AMEV
     double FinalBeamEnergy; 
 
+    // Beam Position
+    double XBeam;
+    double YBeam;
+
     // intermediate variable
     TVector3 BeamDirection;
     TVector3 BeamImpact;
@@ -211,6 +216,10 @@ class Analysis: public NPL::VAnalysis{
     TMugastPhysics* MG;
     //TCATSPhysics* CATS;
     TModularLeafPhysics* ML;
+
+    // Beam object
+    NPL::Beam* Beam;
+
     unsigned int GATCONF_MASTER;
 
     unsigned long long int count ;