diff --git a/Projects/e870/Analysis.cxx b/Projects/e870/Analysis.cxx
index a3345f779b1543846b1114d8a7ffbaf9cc231b40..30b63ee0c7c0ecdc5a5854a8383e78d7d15a47f9 100644
--- a/Projects/e870/Analysis.cxx
+++ b/Projects/e870/Analysis.cxx
@@ -21,12 +21,12 @@
  *****************************************************************************/
 #include <iostream>
 using namespace std;
-#include "NPVAnalysis.h"
 #include "Analysis.h"
 #include "NPAnalysisFactory.h"
 #include "NPDetectorManager.h"
 #include "NPFunction.h"
 #include "NPOptionManager.h"
+#include "NPVAnalysis.h"
 ////////////////////////////////////////////////////////////////////////////////
 Analysis::Analysis() {}
 ////////////////////////////////////////////////////////////////////////////////
@@ -60,6 +60,7 @@ void Analysis::Init() {
 
   // get reaction information
   reaction.ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
+  reaction2.ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
   OriginalBeamEnergy = reaction.GetBeamEnergy();
   std::cout << OriginalBeamEnergy << std::endl;
   // target thickness
@@ -83,7 +84,6 @@ void Analysis::Init() {
   reaction.SetBeamEnergy(FinalBeamEnergy);
 
   // initialize various parameters
-  Rand = TRandom3();
   DetectorNumber = 0;
   ThetaNormalTarget = 0;
   ThetaM2Surface = 0;
@@ -94,6 +94,7 @@ void Analysis::Init() {
   dE = 0;
   BeamDirection = TVector3(0, 0, 1);
   nbTrack = 0;
+  Rand = new TRandom3();
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -104,15 +105,21 @@ void Analysis::TreatEvent() {
   TVector3 BeamDirection;
 
   // For simulation
-  XTarget = 0;
-  YTarget = 0;
   BeamDirection = TVector3(0, 0, 1);
-
+  XTarget = ReactionConditions->GetVertexPositionX();
+  YTarget = ReactionConditions->GetVertexPositionY();
+  XTarget = Rand->Gaus(XTarget * mm, 1 * mm);
+  YTarget = Rand->Gaus(YTarget * mm, 1 * mm);
   BeamImpact = TVector3(XTarget, YTarget, 0);
+  BeamImpact2 = TVector3(0, 0, 0);
+
+  BeamEnergy = ReactionConditions->GetBeamEnergy();
+  BeamEnergy = Rand->Gaus(BeamEnergy, 0.1 * BeamEnergy / 2.35);
+
   // determine beam energy for a randomized interaction point in target
   // 1% FWHM randominastion (E/100)/2.35
   // reaction.SetBeamEnergy(Rand.Gaus(BeamEnergy, BeamEnergy * 1. / 100. / 2.35));
-  // reaction.SetBeamEnergy(BeamEnergy);
+  reaction.SetBeamEnergy(BeamEnergy);
 
   ////////////////////////////////////////////////////////////////////////////
   //////////////////////////////// LOOP on MUST2  ////////////////////////////
@@ -127,14 +134,18 @@ void Analysis::TreatEvent() {
     /************************************************/
 
     // Part 1 : Impact Angle
-    //THIS IS THE NEW PART
+    // 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);
     double Phi = HitDirection.Phi();
+
+    TVector3 HitDirection2 = M2->GetPositionOfInteraction(countMust2) - BeamImpact2;
+    double ThetaNoCATS =HitDirection2.Angle(BeamDirection); 
+
+    TVector3 Coords = M2->GetPositionOfInteraction(countMust2);
     X.push_back(Coords.x());
     Y.push_back(Coords.y());
     Z.push_back(Coords.z());
@@ -142,7 +153,6 @@ void Analysis::TreatEvent() {
     ThetaM2Surface = HitDirection.Angle(-M2->GetTelescopeNormal(countMust2));
     ThetaNormalTarget = HitDirection.Angle(TVector3(0, 0, 1));
 
-
     /************************************************/
 
     /************************************************/
@@ -169,12 +179,13 @@ void Analysis::TreatEvent() {
     Energy = LightTarget.EvaluateInitialEnergy(Energy, TargetThickness * 0.5, Theta);
 
     // What is written in the tree
- 
+
     ThetaLab.push_back(Theta / deg);
     PhiLab.push_back(Phi / deg);
     Ex.push_back(reaction.ReconstructRelativistic(Energy, Theta));
+    ExNoCATS.push_back(reaction2.ReconstructRelativistic(Energy, Theta));
     ELab.push_back(Energy);
-    
+
     /************************************************/
 
   } // end loop MUST2
@@ -184,27 +195,17 @@ void Analysis::TreatEvent() {
 void Analysis::End() {}
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::InitOutputBranch() {
-  // RootOutput::getInstance()->GetTree()->Branch("Exogam", &ExogamData);
   RootOutput::getInstance()->GetTree()->Branch("Ex", &Ex);
+  RootOutput::getInstance()->GetTree()->Branch("ExNoCATS", &ExNoCATS);
   RootOutput::getInstance()->GetTree()->Branch("ELab", &ELab);
   RootOutput::getInstance()->GetTree()->Branch("ThetaLab", &ThetaLab);
   RootOutput::getInstance()->GetTree()->Branch("PhiLab", &PhiLab);
   RootOutput::getInstance()->GetTree()->Branch("ThetaCM", &ThetaCM, "ThetaCM/D");
   RootOutput::getInstance()->GetTree()->Branch("Run", &Run, "Run/I");
-  // 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) {
-  }
-  else {
-    RootOutput::getInstance()->GetTree()->Branch("OriginalELab", &OriginalELab, "OriginalELab/D");
-    RootOutput::getInstance()->GetTree()->Branch("OriginalThetaLab", &OriginalThetaLab, "OriginalThetaLab/D");
-    RootOutput::getInstance()->GetTree()->Branch("BeamEnergy", &BeamEnergy, "BeamEnergy/D");
-  }
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -220,10 +221,7 @@ void Analysis::InitInputBranch() {
     RootInput::getInstance()->GetChain()->SetBranchStatus("fRC_*", true);
     RootInput::getInstance()->GetChain()->SetBranchAddress("ReactionConditions", &ReactionConditions);
   }
-  // RootInput::getInstance()->GetChain()->SetBranchStatus("Exogam", true);
   RootInput::getInstance()->GetChain()->SetBranchStatus("fExo*", true);
-  // RootInput::getInstance()->GetChain()->SetBranchAddress("Exogam", &ExogamData);
-  // RootInput::getInstance()->GetChain()->SetBranchStatus("ZDD", true);
 }
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::ReInitValue() {
@@ -231,9 +229,6 @@ void Analysis::ReInitValue() {
   EDC = -1000;
   BeamEnergy = -1000;
   ThetaCM = -1000;
-  // X = -1000;
-  // Y = -1000;
-  // Z = -1000;
   X.clear();
   Y.clear();
   Z.clear();
@@ -242,6 +237,7 @@ void Analysis::ReInitValue() {
   ThetaLab.clear();
   PhiLab.clear();
   Ex.clear();
+  ExNoCATS.clear();
 }
 
 ////////////////////////////////////////////////////////////////////////////////
diff --git a/Projects/e870/Analysis.h b/Projects/e870/Analysis.h
index afe2d339f85511b3826f370f188c4b3248655e80..fcf71b91ab7a3e667b769772378c57750ceac521 100644
--- a/Projects/e870/Analysis.h
+++ b/Projects/e870/Analysis.h
@@ -1,4 +1,4 @@
-#ifndef Analysis_h 
+#ifndef Analysis_h
 #define Analysis_h
 /*****************************************************************************
  * Copyright (C) 2009-2014    this file is part of the NPTool Project        *
@@ -21,38 +21,39 @@
  *                                                                           *
  *                                                                           *
  *****************************************************************************/
-#include"NPVAnalysis.h"
-#include"NPEnergyLoss.h"
-#include"NPReaction.h"
-#include"RootOutput.h"
-#include"RootInput.h"
+#include "NPEnergyLoss.h"
+#include "NPReaction.h"
+#include "NPVAnalysis.h"
+#include "RootInput.h"
+#include "RootOutput.h"
 #include "TInitialConditions.h"
-#include "TReactionConditions.h"
-#include "TMust2Physics.h"
 #include "TMugastPhysics.h"
+#include "TMust2Physics.h"
+#include "TReactionConditions.h"
 // #include "TCATSPhysics.h"
 #include "TExogamData.h"
+#include <TMath.h>
 #include <TRandom3.h>
 #include <TVector3.h>
-#include <TMath.h>
 
-class Analysis: public NPL::VAnalysis{
-  public:
-    Analysis();
-    ~Analysis();
+class Analysis : public NPL::VAnalysis {
+ public:
+  Analysis();
+  ~Analysis();
 
-  public: 
-    void Init();
-    void TreatEvent();
-    void End();
+ public:
+  void Init();
+  void TreatEvent();
+  void End();
 
   void InitOutputBranch();
   void InitInputBranch();
   void ReInitValue();
   static NPL::VAnalysis* Construct();
- 
-  private:
-  std::vector<double>Ex;
+
+ private:
+  std::vector<double> Ex;
+  std::vector<double> ExNoCATS;
   double ExNoBeam;
   double ExNoProton;
   double EDC;
@@ -64,7 +65,8 @@ class Analysis: public NPL::VAnalysis{
   double OriginalThetaLab;
 
   NPL::Reaction reaction;
-    //	Energy loss table: the G4Table are generated by the simulation
+  NPL::Reaction reaction2;
+  //	Energy loss table: the G4Table are generated by the simulation
   NPL::EnergyLoss LightTarget;
   NPL::EnergyLoss LightMyl;
   NPL::EnergyLoss LightAl;
@@ -73,49 +75,51 @@ class Analysis: public NPL::VAnalysis{
   NPL::EnergyLoss* BeamWindow;
   NPL::EnergyLoss* LightWindow;
 
-  double TargetThickness ;
+  TRandom3* Rand;
+
+  double TargetThickness;
   double WindowsThickness;
   // Beam Energy
-  double OriginalBeamEnergy ; // AMEV
-  double FinalBeamEnergy; 
-  
+  double OriginalBeamEnergy; // AMEV
+  double FinalBeamEnergy;
+
   // intermediate variable
   TVector3 BeamDirection;
   TVector3 BeamImpact;
-  TRandom3 Rand ;
+  TVector3 BeamImpact2;
   int Run;
-  int DetectorNumber  ;
+  int DetectorNumber;
   double ThetaNormalTarget;
-  double ThetaM2Surface ;
-  double ThetaMGSurface ;
-  double Si_E_M2 ;
-  double CsI_E_M2  ;
-  double Energy ;
+  double ThetaM2Surface;
+  double ThetaMGSurface;
+  double Si_E_M2;
+  double CsI_E_M2;
+  double Energy;
   double BeamEnergy;
-  
-  double ThetaGDSurface ;
-  std::vector<double> X ;
-  std::vector<double> Y ;
-  std::vector<double> Z ;
+
+  double ThetaGDSurface;
+  std::vector<double> X;
+  std::vector<double> Y;
+  std::vector<double> Z;
   // Vamos Branches
   unsigned long long int LTS;
-  
+
   // Agata branches
-   double agata_zShift;
-   unsigned long long int TStrack;
-   int nbHits;
-   int nbTrack;
-   float *trackE= new float(100);
-   float *trackX1= new float(100);
-   float *trackY1= new float(100);
-   float *trackZ1= new float(100);
-   float *trackT= new float(100);
-   int *trackCrystalID = new int(100);
-   int nbCores;
-   int *coreId= new int(100);
-   ULong64_t *coreTS= new ULong64_t(100);
-   float *coreE0= new float(100);
-   //
+  double agata_zShift;
+  unsigned long long int TStrack;
+  int nbHits;
+  int nbTrack;
+  float* trackE = new float(100);
+  float* trackX1 = new float(100);
+  float* trackY1 = new float(100);
+  float* trackZ1 = new float(100);
+  float* trackT = new float(100);
+  int* trackCrystalID = new int(100);
+  int nbCores;
+  int* coreId = new int(100);
+  ULong64_t* coreTS = new ULong64_t(100);
+  float* coreE0 = new float(100);
+  //
   double dE;
   double dTheta;
   // Branches and detectors
@@ -126,6 +130,5 @@ class Analysis: public NPL::VAnalysis{
   TInitialConditions* Initial;
   TReactionConditions* ReactionConditions;
   TExogamData* ExogamData;
-
 };
 #endif
diff --git a/Projects/e870/macros/DrawExInvariant.cxx b/Projects/e870/macros/DrawExInvariant.cxx
index c624a907a909411fbe7cbb646db2af37e78d5106..52ea1c5fc68cddbe0ab8814178ab4c8c3cfd9652 100644
--- a/Projects/e870/macros/DrawExInvariant.cxx
+++ b/Projects/e870/macros/DrawExInvariant.cxx
@@ -23,6 +23,7 @@ void DrawExInvariant() {
   std::vector<double>* ThetaLab = 0;
   std::vector<double>* PhiLab = 0;
   std::vector<double>* Ex = 0;
+  std::vector<double>* ExNoCATS = 0;
   TBranch* b_ELab;
   TBranch* b_ThetaLab;
   TBranch* b_PhiLab;
@@ -30,6 +31,7 @@ void DrawExInvariant() {
   TBranch* b_Y;
   TBranch* b_Z;
   TBranch* b_Ex;
+  TBranch* b_ExNoCATS;
 
   t->SetBranchStatus("*", false);
   t->SetBranchStatus("MUST2", true);
@@ -42,10 +44,13 @@ void DrawExInvariant() {
   t->SetBranchAddress("PhiLab", &PhiLab, &b_PhiLab);
   t->SetBranchStatus("Ex", "true");
   t->SetBranchAddress("Ex", &Ex, &b_Ex);
+  t->SetBranchStatus("ExNoCATS", "true");
+  t->SetBranchAddress("ExNoCATS", &ExNoCATS, &b_ExNoCATS);
 
-  TH1F* hExInv = new TH1F("hExInv", "hExInv", 1000, -20, 20);
-  TH1F* hExMM = new TH1F("hExMM", "hExMM", 1000, -20, 20);
-  TH2F* hE1E2 = new TH2F("hE1E2", "hE1E2", 1000, 0, 500, 1000, 0, 500);
+  TH1F* hExInv = new TH1F("hExInv", "hExInv", 200, -5, 5);
+  TH1F* hExMM = new TH1F("hExMM", "hExMM", 200, -5, 5);
+  TH1F* hExMMNoCATS = new TH1F("hExMMNoCATS", "hExMMNoCATS", 200, -5, 5);
+  TH2F* hE1E2 = new TH2F("hE1E2", "hE1E2", 200, 0, 500, 1000, 0, 500);
   TH2F* hELabThetaLab = new TH2F("hELabThetLab", "hELabThetaLab", 1000, 0, 50, 1000, 0, 500);
 
   ////////////////////////////////////////////////////////////////////////////////
@@ -64,6 +69,8 @@ void DrawExInvariant() {
 
     if (M2->Si_E.size() == 2) {
       if ((cuta->IsInside(M2->CsI_E[0], M2->Si_E[0])) && (cutli->IsInside(M2->CsI_E[1], M2->Si_E[1]))) {
+        hExMM->Fill(Ex->at(0));
+        hExMMNoCATS->Fill(ExNoCATS->at(0));
         // Particle1 = alpha
         e_alpha = ELab->at(0);
         theta_alpha = ThetaLab->at(0) * M_PI / 180.;
@@ -83,10 +90,10 @@ void DrawExInvariant() {
         else
           phi_li = phi_li + 180;
         phi_li = phi_li * M_PI / 180.;
-        hExMM->Fill(Ex->at(0));
       }
       else if ((cuta->IsInside(M2->CsI_E[1], M2->Si_E[1])) && (cutli->IsInside(M2->CsI_E[0], M2->Si_E[0]))) {
         hExMM->Fill(Ex->at(1));
+        hExMMNoCATS->Fill(ExNoCATS->at(1));
         // Particle2 = alpha
         e_alpha = ELab->at(1);
         theta_alpha = ThetaLab->at(1) * M_PI / 180.;
@@ -112,8 +119,6 @@ void DrawExInvariant() {
         ////////////////////////////////////////////////////////////////////////////////
         double gamma_alpha = e_alpha / m_alpha + 1;
         double beta_alpha = sqrt(1 - 1 / (pow(gamma_alpha, 2.)));
-        cout << i << endl;
-        cout << "alpha: " << e_alpha << " " << beta_alpha << " " << theta_alpha << endl;
         double p_alpha_mag = gamma_alpha * m_alpha * beta_alpha;
         double etot_alpha = sqrt(pow(p_alpha_mag, 2) + pow(m_alpha, 2));
 
@@ -131,7 +136,6 @@ void DrawExInvariant() {
         ////////////////////////////////////////////////////////////////////////////////
         double gamma_li = e_li / m_7Li + 1;
         double beta_li = sqrt(1 - 1 / (pow(gamma_li, 2.)));
-        cout << "li: " << e_li << " " << beta_li << " " << theta_li << endl;
         double p_li_mag = gamma_li * m_7Li * beta_li;
         double etot_li = sqrt(pow(p_li_mag, 2) + pow(m_7Li, 2));
 
@@ -151,7 +155,6 @@ void DrawExInvariant() {
         TLorentzVector LV_total = LV_alpha + LV_7Li - LV_p;
 
         double ExTot = LV_total.Mag() - m_10Be;
-        cout << ExTot << endl;
         hExInv->Fill(ExTot);
         hE1E2->Fill(e_alpha, e_li);
         hELabThetaLab->Fill(theta_alpha * 180 / M_PI, e_alpha);
@@ -164,18 +167,28 @@ void DrawExInvariant() {
   hExInv->Draw();
   hExMM->SetLineColor(kRed);
   hExMM->Draw("same");
+  hExMMNoCATS->SetLineColor(kGreen);
+  hExMMNoCATS->Draw("same");
 
   TF1* fitExTot = new TF1("fitExTot", "gaus", hExInv->GetXaxis()->GetXmin(), hExInv->GetXaxis()->GetXmax());
+  fitExTot->SetNpx(10000);
   hExInv->Fit(fitExTot, "Q");
-  double meanExTot = fitExTot->GetParameter(1);
-  std::cout << "Center of hExInv: " << meanExTot << " MeV" << std::endl;
+  double sigmaExTot = fitExTot->GetParameter(2);
+  std::cout << "FWHM of hExInv: " << sigmaExTot*2.35 << " MeV" << std::endl;
 
   double fitRangeMin = 0.0;
   double fitRangeMax = 15.0;
   TF1* fitEx = new TF1("fitEx", "gaus", fitRangeMin, fitRangeMax);
-  hExMM->Fit(fitEx, "", "", -10, 10);
-  double meanEx = fitEx->GetParameter(1);
-  std::cout << "Center of hEx: " << meanEx << " MeV" << std::endl;
+  fitEx->SetNpx(10000);
+  hExMM->Fit(fitEx, "Q", "", -10, 10);
+  double sigmaEx = fitEx->GetParameter(2);
+  std::cout << "FWHM of hEx: " << sigmaEx*2.35 << " MeV" << std::endl;
+
+  TF1* fitExNoCATS = new TF1("fitExNoCATS", "gaus", fitRangeMin, fitRangeMax);
+  fitExNoCATS->SetNpx(10000);
+  hExMMNoCATS->Fit(fitExNoCATS, "Q", "", -10, 10);
+  double sigmaExNoCATS = fitExNoCATS->GetParameter(2);
+  std::cout << "FWHM of hExNoCATS: " << sigmaExNoCATS*2.35 << " MeV" << std::endl;
 
   TLegend* legend = new TLegend(0.75, 0.75, 0.85, 0.85);
   legend->AddEntry(hExInv, "hExInv", "l"); // "l" for line