diff --git a/NPAnalysis/MUGAST/Analysis.cxx b/NPAnalysis/MUGAST/Analysis.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..4849a27d93ec319a547ed1f02744becf92f1f8ea
--- /dev/null
+++ b/NPAnalysis/MUGAST/Analysis.cxx
@@ -0,0 +1,283 @@
+#include "Analysis.h"
+
+int main(int argc, char** argv){
+ // command line parsing
+ NPOptionManager* myOptionManager = NPOptionManager::getInstance(argc,argv);
+
+ // Instantiate RootInput
+ string runToReadfileName = myOptionManager->GetRunToReadFile();
+ RootInput:: getInstance("RunToTreat.txt");
+ TChain* Chain = RootInput:: getInstance()->GetChain();
+ // if input files are not given, use those from TAsciiFile
+ if (myOptionManager->IsDefault("DetectorConfiguration")) {
+ string name = RootInput::getInstance()->DumpAsciiFile("DetectorConfiguration");
+ myOptionManager->SetDetectorFile(name);
+ }
+
+ if (myOptionManager->IsDefault("EventGenerator")) {
+ string name = RootInput::getInstance()->DumpAsciiFile("EventGenerator");
+ myOptionManager->SetReactionFile(name);
+ }
+
+ // get input files from NPOptionManager
+ string detectorfileName = myOptionManager->GetDetectorFile();
+ string OutputfileName = myOptionManager->GetOutputFile();
+
+ // Instantiate RootOutput
+ RootOutput::getInstance("Analysis/"+OutputfileName, "ResultTree");
+ // RootOutput::getInstance()->GetFile()->SetCompressionLevel(0);
+ // Instantiate the detector using a file
+ NPA::DetectorManager* myDetector = new DetectorManager();
+ myDetector->ReadConfigurationFile(detectorfileName);
+ // Attach new branch
+ InitOutputBranch();
+ InitInputBranch();
+
+ // Instantiate the Reaction
+ NPL::Reaction* myReaction = new Reaction ;
+ myReaction -> ReadConfigurationFile(myOptionManager->GetReactionFile());
+
+ ////////////////////////////////////////////////////////
+ // Get pointer to the different detector
+
+ TMust2Physics* M2 = (TMust2Physics*) myDetector -> GetDetector("MUST2");
+ GaspardTracker* GD = (GaspardTracker*) myDetector -> GetDetector("GASPARD");
+
+ // Beam Energy
+ double BeamEnergy = 4.0* 74; // AMEV
+
+ // intermediate variable
+ TRandom3 Rand = TRandom3();
+ int DetectorNumber = 0 ;
+ double ThetaNormalTarget = 0 ;
+ double ThetaM2Surface = 0;
+ double X_M2 = 0 ;
+ double Y_M2 = 0 ;
+ double Z_M2 = 0 ;
+ double Si_E_M2 = 0 ;
+ double CsI_E_M2 = 0 ;
+ double Energy = 0;
+ double E_M2 = 0;
+ double Si_X_M2 = 0;
+ double Si_Y_M2 = 0;
+ double ZTarget = 0;
+ double TargetThickness = 9*micrometer;
+
+ double ThetaGDSurface = 0;
+ double X_GD = 0 ;
+ double Y_GD = 0 ;
+ double Z_GD = 0 ;
+ double Si_E_GD = 0 ;
+ double E_GD = 0;
+ double Si_X_GD = 0;
+ double Si_Y_GD = 0;
+
+
+ // Get number of events to treat
+ cout << endl << "///////// Starting Analysis ///////// "<< endl;
+ int nentries = Chain->GetEntries();
+ cout << " Number of Event to be treated : " << nentries << endl;
+ clock_t begin = clock();
+ clock_t end = begin;
+ cout.precision(5);
+ //////////////////////////////////////////////////////////////////////////////
+ // main loop on entries //
+ for (int i = 0 ; i < nentries; i++) {
+ if (i%10000 == 0 && i!=0) {
+ end = clock();
+ long double TimeElapsed = (long double) (end-begin) / CLOCKS_PER_SEC;
+ double percent = (double)i/nentries;
+ double TimeToWait = (TimeElapsed/percent) - TimeElapsed;
+ cout << " "<< flush;
+ cout << "\r Progression:"
+ << percent*100 << " % \t | \t Remaining time : ~"
+ << TimeToWait <<"s | Analysis Rate : "
+ << (double) i/TimeElapsed << flush;
+ }
+ else if (i == nentries-1) cout << "\r Progression:" << " 100% " <<endl;
+ // Get the raw Data
+ Chain -> GetEntry(i);
+ // Clear previous Event
+ myDetector->ClearEventPhysics();
+ // Build the current event
+ myDetector->BuildPhysicalEvent();
+ // Reinitiate calculated variable
+ ReInitValue();
+
+ // Get the Init information on beam position and energy
+ // and apply by hand the experimental resolution
+ // This is because the beam diagnosis are not simulated
+ // PPAC position resolution on target is assumed to be 1mm
+ double XTarget = 0;
+ double YTarget = 0;
+ TVector3 BeamDirection = TVector3(0,0,1);
+
+ // Beam energy is measured using F3 and F2 plastic TOF
+ BeamEnergy = BeamCD2.Slow(BeamEnergy,TargetThickness/2.,0);
+ myReaction->SetBeamEnergy(BeamEnergy);
+
+ //////////////////////////// LOOP on MUST2 //////////////////
+ for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
+ /************************************************/
+ //Part 0 : Get the usefull Data
+ // MUST2
+ int X = M2->Si_X[countMust2];
+ int Y = M2->Si_Y[countMust2];
+ int TelescopeNumber = M2->TelescopeNumber[countMust2];
+ Si_X_M2 = X ;
+ Si_Y_M2 = Y ;
+
+ /************************************************/
+ // Part 1 : Impact Angle
+ ThetaM2Surface = 0;
+ ThetaNormalTarget = 0;
+ if(XTarget>-1000 && YTarget>-1000){
+ TVector3 BeamImpact(XTarget,YTarget,0);
+ TVector3 HitDirection = M2 -> GetPositionOfInteraction(countMust2) - BeamImpact ;
+ ThetaLab = HitDirection.Angle( BeamDirection );
+
+ ThetaM2Surface = HitDirection.Angle(- M2 -> GetTelescopeNormal(countMust2) );
+ ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
+ X_M2 = M2 -> GetPositionOfInteraction(countMust2).X() ;
+ Y_M2 = M2 -> GetPositionOfInteraction(countMust2).Y() ;
+ Z_M2 = M2 -> GetPositionOfInteraction(countMust2).Z() ;
+ }
+
+ else{
+ BeamDirection = TVector3(-1000,-1000,-1000);
+ ThetaM2Surface = -1000 ;
+ ThetaNormalTarget = -1000 ;
+ }
+
+ /************************************************/
+
+ /************************************************/
+ // Part 2 : Impact Energy
+ Energy = ELab = 0;
+ Si_E_M2 = M2->Si_E[countMust2];
+ CsI_E_M2= M2->CsI_E[countMust2];
+
+ // if CsI
+ if(CsI_E_M2>0 ){
+ // The energy in CsI is calculate form dE/dx Table because
+ // 20um resolution is poor
+ Energy =
+ LightSi.EvaluateEnergyFromDeltaE(Si_E_M2,300*micrometer,
+ ThetaM2Surface, 0.01*MeV,
+ 450.*MeV,0.001*MeV ,1000);
+ E_M2=CsI_E_M2;
+ }
+
+ else
+ Energy = Si_E_M2;
+
+ E_M2 += Si_E_M2;
+
+ // Evaluate energy using the thickness
+ ELab = LightAl.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaM2Surface);
+ // Target Correction
+ ELab = LightCD2.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
+ /************************************************/
+
+ /************************************************/
+ // Part 3 : Excitation Energy Calculation
+ Ex = myReaction -> ReconstructRelativistic( ELab , ThetaLab );
+ ThetaLab=ThetaLab/deg;
+ /************************************************/
+
+ /************************************************/
+ // Part 4 : Theta CM Calculation
+ ThetaCM = myReaction -> EnergyLabToThetaCM( ELab , 0)/deg;
+ /************************************************/
+ }//end loop MUST2
+
+ ////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////
+ //////////////////////////// LOOP on GASPARD //////////////////
+ if(GD->GetEnergyDeposit()>0){
+ /************************************************/
+ // Part 1 : Impact Angle
+ ThetaGDSurface = 0;
+ ThetaNormalTarget = 0;
+ if(XTarget>-1000 && YTarget>-1000){
+ TVector3 BeamImpact(XTarget,YTarget,0);
+ TVector3 HitDirection = GD -> GetPositionOfInteraction() - BeamImpact ;
+ ThetaLab = HitDirection.Angle( BeamDirection );
+
+ ThetaGDSurface = HitDirection.Angle( TVector3(0,0,1) ) ;
+ ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
+ X_GD = GD -> GetPositionOfInteraction().X() ;
+ Y_GD = GD -> GetPositionOfInteraction().Y() ;
+ Z_GD = GD -> GetPositionOfInteraction().Z() ;
+ }
+
+ else{
+ BeamDirection = TVector3(-1000,-1000,-1000);
+ ThetaGDSurface = -1000 ;
+ ThetaNormalTarget = -1000 ;
+ }
+
+ /************************************************/
+
+ /************************************************/
+ // Part 2 : Impact Energy
+ Energy = ELab = 0;
+ Energy = GD->GetEnergyDeposit();
+
+ ELab = LightAl.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaGDSurface);
+ // Target Correction
+ ELab = LightCD2.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
+ /************************************************/
+
+ /************************************************/
+ // Part 3 : Excitation Energy Calculation
+ Ex = myReaction -> ReconstructRelativistic( ELab , ThetaLab );
+ if(Ex!=Ex){
+ cout << ELab << " " << ThetaLab << " " << Ex<< endl;
+ }
+ ThetaLab=ThetaLab/deg;
+ /************************************************/
+
+ /************************************************/
+ // Part 4 : Theta CM Calculation
+ ThetaCM = myReaction -> EnergyLabToThetaCM( ELab , 0)/deg;
+ /************************************************/
+ }//end loop GASPARD
+
+ if(ELab>0)
+ RootOutput::getInstance()->GetTree()->Fill();
+ }// loop over events
+
+ cout << "A total of " << nentries << " event has been annalysed " << endl ;
+
+ RootOutput::getInstance()->Destroy();
+ RootInput::getInstance()->Destroy();
+ NPOptionManager::getInstance()->Destroy();
+ /////////////////////////////////////////////////////////////////////////////
+ return 0 ;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void InitOutputBranch() {
+ RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void InitInputBranch(){
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&Init );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
+}
+////////////////////////////////////////////////////////////////////////////////
+void ReInitValue(){
+ Ex = -1000 ;
+ ELab = -1000;
+ ThetaLab = -1000;
+ ThetaCM = -1000;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
diff --git a/NPAnalysis/MUGAST/Analysis.h b/NPAnalysis/MUGAST/Analysis.h
new file mode 100755
index 0000000000000000000000000000000000000000..4e36754d2768993edc66b6412c673760b983ee2e
--- /dev/null
+++ b/NPAnalysis/MUGAST/Analysis.h
@@ -0,0 +1,67 @@
+// You can use this file to declare your spectra, file, energy loss , ... and whatever you want.
+// This way you can remove all unnecessary declaration in the main programm.
+// In order to help debugging and organizing we use Name Space.
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// -------------------------------------- VARIOUS INCLUDE ---------------------------------------
+
+// NPL
+#include "DetectorManager.h"
+#include "NPOptionManager.h"
+#include "NPReaction.h"
+#include "RootInput.h"
+#include "RootOutput.h"
+#include "TMust2Physics.h"
+#include "GaspardTracker.h"
+
+#include "TInitialConditions.h"
+#include "NPEnergyLoss.h"
+using namespace NPL ;
+// STL C++
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <cmath>
+#include <cstdlib>
+using namespace std;
+// ROOT
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TVector3.h>
+#include <TRandom3.h>
+#include <TMath.h>
+#include <TObject.h>
+
+// ----------------------------------------------------------------------------------------------
+void InitOutputBranch() ;
+void InitInputBranch() ;
+void ReInitValue() ;
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// ----------------------------------- DOUBLE, INT, BOOL AND MORE -------------------------------
+namespace VARIABLE{
+ double Ex;
+ double ELab;
+ double ThetaLab;
+ double ThetaCM;
+ TInitialConditions* Init = new TInitialConditions();
+}
+
+using namespace VARIABLE ;
+// ----------------------------------------------------------------------------------------------
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// -----------------------------------ENERGY LOSS----------------------------------------------
+namespace ENERGYLOSS{
+ // Energy loss table: the G4Table are generated by the simulation
+ EnergyLoss LightCD2 = EnergyLoss("proton_CD2.G4table","G4Table",100 );
+ EnergyLoss LightAl = EnergyLoss("proton_Al.G4table","G4Table",10);
+ EnergyLoss LightSi = EnergyLoss("proton_Si.G4table","G4Table",10);
+ EnergyLoss BeamCD2 = EnergyLoss("Kr74[0.0]_CD2.G4table","G4Table",100);
+}
+
+using namespace ENERGYLOSS ;
+// ----------------------------------------------------------------------------------------------
+/////////////////////////////////////////////////////////////////////////////////////////////////
+
+
diff --git a/NPAnalysis/MUGAST/Makefile b/NPAnalysis/MUGAST/Makefile
new file mode 100755
index 0000000000000000000000000000000000000000..abdac126fca99a6b559dbcf1b3f71fbd0baa1c94
--- /dev/null
+++ b/NPAnalysis/MUGAST/Makefile
@@ -0,0 +1,31 @@
+# include same architecture file than for NPLib
+# so that consistency is ensured
+include $(NPLIB)/Makefile.arch
+
+# additional libraries
+LIBRARY = `$(NPLIB)/liblist`
+LIBRARY += -L$(CLHEP_LIB_DIR) -l$(CLHEP_LIB)
+
+PROGRAMS = NPAnalysis
+
+all: $(PROGRAMS)
+
+NPAnalysis: Analysis.o
+ $(LD) $(LDFLAGS) $^ $(LIBS) $(LIBRARY) $(OutPutOpt) $@
+ @echo "$@ done"
+
+
+# rule for creating .o from .cxx
+.SUFFIXES: .$(SrcSuf)
+.$(SrcSuf).$(ObjSuf):
+ $(CXX) $(CXXFLAGS) $(INCLUDE) -c $<
+
+# some cleaning
+clean:
+ rm -rf *.o
+
+distclean:
+ make clean; rm $(PROGRAMS)
+
+# dependences
+Analysis.o: Analysis.cxx Analysis.h
diff --git a/NPAnalysis/MUGAST/NPAnalysis b/NPAnalysis/MUGAST/NPAnalysis
new file mode 100755
index 0000000000000000000000000000000000000000..0d3882cc25d0c134a2037c81bb354721d0ab6dd2
Binary files /dev/null and b/NPAnalysis/MUGAST/NPAnalysis differ
diff --git a/NPAnalysis/MUGAST/RunToTreat.txt b/NPAnalysis/MUGAST/RunToTreat.txt
new file mode 100755
index 0000000000000000000000000000000000000000..237ac3af20e33924cd9305e8124184068af7fe29
--- /dev/null
+++ b/NPAnalysis/MUGAST/RunToTreat.txt
@@ -0,0 +1,5 @@
+TTreeName
+ SimulatedTree
+RootFileName
+ ../../Outputs/Simulation/myResult.root
+
diff --git a/NPAnalysis/MUGAST/configs/ConfigMust2.dat b/NPAnalysis/MUGAST/configs/ConfigMust2.dat
new file mode 100755
index 0000000000000000000000000000000000000000..2c36d1d62aba8284349f5afb674c5f63153299da
--- /dev/null
+++ b/NPAnalysis/MUGAST/configs/ConfigMust2.dat
@@ -0,0 +1,26 @@
+ConfigMust2
+ MAX_STRIP_MULTIPLICITY 1
+ STRIP_ENERGY_MATCHING_NUMBER_OF_SIGMA 5
+ STRIP_ENERGY_MATCHING_SIGMA 0.02
+ DISABLE_CHANNEL MM1STRY12
+ DISABLE_CHANNEL MM2STRY12
+ DISABLE_CHANNEL MM3STRY12
+ DISABLE_CHANNEL MM4STRY12
+ DISABLE_CHANNEL MM1STRX12
+ DISABLE_CHANNEL MM2STRX12
+ DISABLE_CHANNEL MM3STRX12
+ DISABLE_CHANNEL MM4STRX12
+ DISABLE_CHANNEL MM1STRY124
+ DISABLE_CHANNEL MM4STRX3
+ DISABLE_CHANNEL MM4STRX4
+ DISABLE_CHANNEL MM4STRX5
+ DISABLE_CHANNEL MM4STRX6
+ DISABLE_CHANNEL MM4STRX8
+ DISABLE_CHANNEL MM4STRX31
+ DISABLE_ALL MM5
+ DISABLE_ALL MM6
+ DISABLE_ALL MM7
+ DISABLE_ALL MM8
+ SI_X_E_RAW_THRESHOLD 0
+ CSI_E_RAW_THRESHOLD 0
+ CSI_SIZE 256
diff --git a/NPAnalysis/MUGAST/configs/ConfigSSSD.dat b/NPAnalysis/MUGAST/configs/ConfigSSSD.dat
new file mode 100755
index 0000000000000000000000000000000000000000..a56fa0b905431080044528396495a410ba90dc7b
--- /dev/null
+++ b/NPAnalysis/MUGAST/configs/ConfigSSSD.dat
@@ -0,0 +1,4 @@
+ConfigSSSD
+ PEDESTAL_THRESHOLD 0
+ DISABLE_CHANNEL SI2STR1
+ DISABLE_CHANNEL SI4STR13