diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index 66eea8347260e59cb66e2af6b34b07d1435fb0c1..14f8868f68b5eb2a55057dd379ae6e246a58fa29 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -141,10 +141,6 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
   m_ReactionConditions->Clear();
   m_shoot=true;
   }
-  else if (((in-m_StepSize) <= 1E-20) && m_shoot) { // last step
-    return true;
-  }
-
 
   // If the condition is met, the event is generated
   if (ratio < rand) {
@@ -153,13 +149,17 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
     if(m_ReactionType=="QFSReaction"){
         if ( m_shoot && m_QFS.IsAllowed() ) {
             return true;
-        } 
+        }
+        else
+          m_shoot=false;
     }
    
     else if(m_ReactionType=="TwoBodyReaction"){
         if ( m_shoot && m_Reaction.IsAllowed(PrimaryTrack->GetKineticEnergy()) ) {
             return true;
         } 
+        else
+          m_shoot=false;
     }
   }
 
@@ -170,8 +170,6 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
     m_PreviousEnergy = PrimaryTrack->GetKineticEnergy();
     m_PreviousDirection = PrimaryTrack->GetMomentumDirection();
   }
-  else
-    return true;
 
   return false;
 }
@@ -202,6 +200,7 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
     double rand   = G4RandFlat::shoot();
     double length = (1-rand) * (m_PreviousLength);
     double reac_energy = energy + (1-rand) * (m_PreviousEnergy - energy);
+
     G4ThreeVector ldir = m_PreviousDirection;
     ldir *= length;
     localPosition = localPosition - ldir;
diff --git a/Projects/MUGAST/16Odp17O_gs.reaction b/Projects/MUGAST/16Odp17O_gs.reaction
new file mode 100644
index 0000000000000000000000000000000000000000..2fdc85b1b2c1449e24eb1e0963e110e06817bf10
--- /dev/null
+++ b/Projects/MUGAST/16Odp17O_gs.reaction
@@ -0,0 +1,29 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%% Reaction file: 56Ni(d,p)57Ni at VAMOS-AGATA %%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+  Particle= 16O
+  Energy= 96 MeV
+  SigmaEnergy= 0 MeV
+  SigmaX= 0 mm
+  SigmaY= 0 mm
+  SigmaThetaX= 0.0 deg
+  SigmaPhiY= 0.0 deg
+  MeanThetaX= 0 mm
+  MeanPhiY= 0 mm
+  MeanX= 0 mm
+  MeanY= 0 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+TwoBodyReaction
+ Beam= 16O
+ Target= 2H
+ Light= 1H
+ Heavy= 17O
+ ExcitationEnergyLight= 0.0 MeV
+ ExcitationEnergyHeavy= 0.0 MeV
+ CrossSectionPath= 16Odp17O_s12_6AMeV.dat CS
+ ShootLight= 1
+ ShootHeavy= 1
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Projects/MUGAST/Analysis.cxx b/Projects/MUGAST/Analysis.cxx
index 4c86b827c1cfcf406fa6e03d761020e8e296ec9c..d01b32e73c3d8549381af650b288b0fd7680c33b 100644
--- a/Projects/MUGAST/Analysis.cxx
+++ b/Projects/MUGAST/Analysis.cxx
@@ -76,9 +76,10 @@ void Analysis::Init() {
   LightTarget = NPL::EnergyLoss(light+"_"+TargetMaterial+".G4table","G4Table",100 );
   LightAl = NPL::EnergyLoss(light+"_Al.G4table","G4Table",100);
   LightSi = NPL::EnergyLoss(light+"_Si.G4table","G4Table",100);
-  BeamCD2 = NPL::EnergyLoss(beam+"_"+TargetMaterial+".G4table","G4Table",100);
+  BeamCD2 = NPL::EnergyLoss(beam+"_"+TargetMaterial+".G4table","G4Table",00);
 
   FinalBeamEnergy = BeamCD2.Slow(OriginalBeamEnergy,TargetThickness*0.5,0); 
+  //FinalBeamEnergy = OriginalBeamEnergy; 
   cout << "Original beam energy: " << OriginalBeamEnergy  << " MeV      Mid-target beam energy: " << FinalBeamEnergy << "MeV " << endl; 
   reaction.SetBeamEnergy(FinalBeamEnergy);
 
diff --git a/Projects/MUGAST/ShowResults.C b/Projects/MUGAST/ShowResults.C
index a2545748330069d630b22ab6eb54f93527a0bb32..619b3a2381dabc254431c77964c09518ef4d3424 100644
--- a/Projects/MUGAST/ShowResults.C
+++ b/Projects/MUGAST/ShowResults.C
@@ -72,10 +72,10 @@ void ShowResults(){
   reaction->GetKinematicLine3()->Draw("c");
 
   c1->cd(2);
-  TH1F* hEx = new TH1F("hEx", "hEx",1000, -1, 5);
+  TH1F* hEx = new TH1F("hEx", "hEx",240, -1, 5);
   t->Draw("Ex>>hEx","ThetaLab>100 && ThetaLab<156","col");
   hEx->GetXaxis()->SetTitle("Ex (MeV)");
-  hEx->GetYaxis()->SetTitle("counts/100 keV");
+  hEx->GetYaxis()->SetTitle("counts/25 keV");
  
   c1->cd(3);
   TH1F *hCM = new TH1F("hCM", "hCM", 36, 0, 180); 
@@ -87,7 +87,6 @@ void ShowResults(){
     }
   }
 
-/*  
   TCanvas *c2 = new TCanvas("c2", "Control", 1000, 1000);
   c2->Divide(2,2);
   c2->cd(1);
@@ -100,7 +99,6 @@ void ShowResults(){
   t->Draw("ELab-OriginalELab>>hcE","","col");
   TLine* lE = new TLine(0,0,60,60);
   //lE->Draw();
-*/
 
 }
 
diff --git a/Projects/MUGAST/e793.detector b/Projects/MUGAST/e793.detector
index 1bb1133dc9fbab2031012fdd6b18b505317468f7..e25775f0f4671a9507cda0cc9db788cd76a0ed82 100644
--- a/Projects/MUGAST/e793.detector
+++ b/Projects/MUGAST/e793.detector
@@ -1,13 +1,14 @@
 %%%%%%%%%%Detector%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
 Target
- THICKNESS= 0.01 micrometer
+ THICKNESS= 9.43 micrometer
  ANGLE= 0 deg
  RADIUS= 10 mm
  MATERIAL= CD2
- X= 0
- Y= 0
- Z= 0
+ X= 0 mm
+ Y= 0 mm
+ Z= 0 mm
+ NbSlices= 100
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 CATSDetector
diff --git a/Projects/MUGAST/run.mac b/Projects/MUGAST/run.mac
new file mode 100644
index 0000000000000000000000000000000000000000..02b80e5efd512013cc3fcf9bb00160a53811fec6
--- /dev/null
+++ b/Projects/MUGAST/run.mac
@@ -0,0 +1,3 @@
+/run/beamOn 10000
+/gen/open 16Odp17O_gs.reaction
+/run/beamOn 10000