diff --git a/Inputs/DetectorConfiguration/MUGAST_Manu.detector b/Inputs/DetectorConfiguration/MUGAST_Manu.detector
index c7f70ebf7c9aad5c9ea812c73af99100c8d4b7d3..9135644bd2b4a39353040e5aff69c3f756869b35 100644
--- a/Inputs/DetectorConfiguration/MUGAST_Manu.detector
+++ b/Inputs/DetectorConfiguration/MUGAST_Manu.detector
@@ -9,6 +9,7 @@ Target
 	X= 0
 	Y= 0
 	Z= 0
+	NBLAYERS= 100
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Inputs/EventGenerator/30Pdp.reaction b/Inputs/EventGenerator/30Pdp.reaction
index 07e4253bdba34fc320b84569e476ebd2712a1ce2..fcbada0df7b5eaaead3f74c5297a857169139dc3 100644
--- a/Inputs/EventGenerator/30Pdp.reaction
+++ b/Inputs/EventGenerator/30Pdp.reaction
@@ -3,7 +3,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Beam
    Particle= 30P
-   Energy= 180
+   Energy= 360
    SigmaEnergy= 0
    SigmaX= 0
    SigmaY= 0
diff --git a/Projects/MUGAST/Analysis.cxx b/Projects/MUGAST/Analysis.cxx
index 42349c3a5482f16f229cb128a53ef77ff8a564f4..88c65f9902eb41c4a044de6207ecd03d0667af9b 100644
--- a/Projects/MUGAST/Analysis.cxx
+++ b/Projects/MUGAST/Analysis.cxx
@@ -21,10 +21,11 @@
  *****************************************************************************/
 #include<iostream>
 using namespace std;
-#include"Analysis.h"
-#include"NPAnalysisFactory.h"
-#include"NPDetectorManager.h"
-#include"NPOptionManager.h"
+
+#include "Analysis.h"
+#include "NPAnalysisFactory.h"
+#include "NPDetectorManager.h"
+#include "NPOptionManager.h"
 ////////////////////////////////////////////////////////////////////////////////
 Analysis::Analysis(){
 }
@@ -33,48 +34,65 @@ Analysis::~Analysis(){
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-void Analysis::Init(){
+void Analysis::Init() {
+  // initialize input and output branches
   InitOutputBranch();
   InitInputBranch();
   
-  M2  = (TMust2Physics*) m_DetectorManager -> GetDetector("MUST2Array");
-  GD = (GaspardTracker*)  m_DetectorManager -> GetDetector("GaspardTracker");
+  // get MUST2 and Gaspard objects
+  M2 = (TMust2Physics*)  m_DetectorManager -> GetDetector("MUST2Array");
+  GD = (GaspardTracker*) m_DetectorManager -> GetDetector("GaspardTracker");
+
+  // energy losses
   LightCD2 = EnergyLoss("proton_CD2.G4table","G4Table",100 );
   LightAl = EnergyLoss("proton_Al.G4table","G4Table",100);
   LightSi = EnergyLoss("proton_Si.G4table","G4Table",100);
-  BeamCD2 = EnergyLoss("Na24[0.0]_CD2.G4table","G4Table",100);
+//  BeamCD2 = EnergyLoss("Na24[0.0]_CD2.G4table","G4Table",100);
+  BeamCD2 = EnergyLoss("P30_CD2.G4table","G4Table",100);
+
+  // get reaction information
   myReaction = new NPL::Reaction();
   myReaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
-   TargetThickness = 18*micrometer;
   OriginalBeamEnergy = myReaction->GetBeamEnergy();
-   Rand = TRandom3();
-   DetectorNumber = 0 ;
-   ThetaNormalTarget = 0 ;
-   ThetaM2Surface = 0;
-   Si_E_M2 = 0 ;
-   CsI_E_M2 = 0 ;
-   Energy = 0;
-   E_M2 = 0;
-  
-   ThetaGDSurface = 0;
-   X_GD = 0 ;
-   Y_GD = 0 ;
-   Z_GD = 0 ;
-   Si_E_GD = 0 ;
-   E_GD = 0;
-   Si_X_GD = 0;
-   Si_Y_GD = 0;
+
+  // target thickness
+  TargetThickness = m_DetectorManager->GetTargetThickness()*micrometer;
+
+  // initialize various parameters
+  Rand = TRandom3();
+  DetectorNumber = 0;
+  ThetaNormalTarget = 0;
+  ThetaM2Surface = 0;
+  Si_E_M2 = 0;
+  CsI_E_M2 = 0;
+  Energy = 0;
+  E_M2 = 0;
+
+  ThetaGDSurface = 0;
+  X_GD = 0;
+  Y_GD = 0;
+  Z_GD = 0;
+  Si_E_GD = 0;
+  E_GD = 0;
+  Si_X_GD = 0;
+  Si_Y_GD = 0;
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-void Analysis::TreatEvent(){
+void Analysis::TreatEvent() {
   // Reinitiate calculated variable
   ReInitValue();
+
+  // assume beam is centered and parallel to z axis
   double XTarget = 0;
   double YTarget = 0;
   TVector3 BeamDirection = TVector3(0,0,1);
-  double BeamEnergy = BeamCD2.Slow(OriginalBeamEnergy,Rand.Uniform(0,TargetThickness),0);
+
+  // determine beam energy for a randomized interaction point in target
+//  double BeamEnergy = BeamCD2.Slow(OriginalBeamEnergy, Rand.Uniform(0,TargetThickness), 0);
+  double BeamEnergy = BeamCD2.Slow(OriginalBeamEnergy, TargetThickness/2., 0);
   myReaction->SetBeamEnergy(BeamEnergy);
+
   //////////////////////////// LOOP on MUST2 //////////////////
   for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
     /************************************************/
diff --git a/Projects/MUGAST/ShowResults.C b/Projects/MUGAST/ShowResults.C
new file mode 100644
index 0000000000000000000000000000000000000000..aa9a4893ad1b55d42640f9e22432baaf5f7b2e5d
--- /dev/null
+++ b/Projects/MUGAST/ShowResults.C
@@ -0,0 +1,22 @@
+void ShowResults()
+{
+   // get tree   
+   TFile *f = new TFile("../../Outputs/Analysis/PhysicsTree.root");
+   TTree *t = (TTree*) f->Get("PhysicsTree");
+
+   // draw 
+   TCanvas *c1 = new TCanvas("c1", "kinematic information", 600, 600);
+   c1->Draw();
+   TH2F *hk = new TH2F("hk", "hk", 90, 90, 180, 200, 0, 12);
+   hk->GetXaxis()->SetTitle("#Theta_{lab} (deg)");
+   hk->GetYaxis()->SetTitle("E_{p} (MeV)");
+   t->Draw("ELab:ThetaLab>>hk");
+   // inset
+   TPad *pad = new TPad("pad1", "excitation energy", 0.45, 0.45, 0.87, 0.87);
+   pad->Draw();
+   pad->cd();
+   TH1F *hx = new TH1F("hx", "hx", 150, 5, 8);
+   hx->SetXTitle("E_{X} (MeV)");
+   hx->SetYTitle("counts / (20 keV)");
+   t->Draw("Ex>>hx");
+}