diff --git a/Inputs/DetectorConfiguration/e530.detector b/Inputs/DetectorConfiguration/e530.detector
index 93a4408ed63fcd9853486411b498e3f6db6c706d..8bf4603f47974ff7e7fc522b03b5a745455c2968 100644
--- a/Inputs/DetectorConfiguration/e530.detector
+++ b/Inputs/DetectorConfiguration/e530.detector
@@ -25,7 +25,7 @@ Target
ANGLE= 0
X= 0
Y= 0
- Z= -40
+ Z= 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
AnnularS1
Z= -158.7
diff --git a/NPAnalysis/must2/Makefile b/NPAnalysis/must2/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..215794b2868b46e37601e95c6dedb6b1f84d7c11
--- /dev/null
+++ b/NPAnalysis/must2/Makefile
@@ -0,0 +1,10 @@
+
+Analyse:
+ make -C ./src
+
+clean:
+ make clean -C ./src
+
+distclean:
+ make clean -C ./src
+ rm Analysis
diff --git a/NPAnalysis/must2/RunToTreat.txt b/NPAnalysis/must2/RunToTreat.txt
new file mode 100644
index 0000000000000000000000000000000000000000..698f27aa2f5d83de671b3370b1ef6ccf2a561d21
--- /dev/null
+++ b/NPAnalysis/must2/RunToTreat.txt
@@ -0,0 +1,4 @@
+TTreeName
+ SimulatedTree
+RootFileName
+ ../../Outputs/Simulation/myResult.root
diff --git a/NPAnalysis/must2/configs/ConfigMust2.dat b/NPAnalysis/must2/configs/ConfigMust2.dat
new file mode 100644
index 0000000000000000000000000000000000000000..2e760ba65c8e895d2747401508988924f637c1d1
--- /dev/null
+++ b/NPAnalysis/must2/configs/ConfigMust2.dat
@@ -0,0 +1,6 @@
+ConfigMust2
+ MAX_STRIP_MULTIPLICITY 1
+ SI_X_E_RAW_THRESHOLD 0
+ SILI_E_RAW_THRESHOLD 0
+ SILI_E_THRESHOLD 0
+
diff --git a/NPAnalysis/must2/include/ObjectManager.hh b/NPAnalysis/must2/include/ObjectManager.hh
new file mode 100644
index 0000000000000000000000000000000000000000..101c92289bbefcae8adb905885c04b987f3481b0
--- /dev/null
+++ b/NPAnalysis/must2/include/ObjectManager.hh
@@ -0,0 +1,217 @@
+// 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 ---------------------------------------
+
+// NPA
+#include "DetectorManager.h"
+#include "NPOptionManager.h"
+#include "TMust2Physics.h"
+
+// STL C++
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <cmath>
+#include <cstdlib>
+#include <time.h>
+
+// ROOT
+#include <TROOT.h>
+#include <TCutG.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TLeaf.h>
+#include <TVector3.h>
+#include <TRandom.h>
+
+// NPL
+#include "TPlasticData.h"
+#include "NPReaction.h"
+#include "RootInput.h"
+#include "RootOutput.h"
+#include "TInitialConditions.h"
+
+// Use CLHEP System of unit and Physical Constant
+#include "CLHEP/Units/GlobalSystemOfUnits.h"
+#include "CLHEP/Units/PhysicalConstants.h"
+
+
+// ----------------------------------------------------------------------------------------------
+double ThetaCalculation (TVector3 A , TVector3 B) ;
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// ----------------------------------- DOUBLE, INT, BOOL AND MORE -------------------------------
+namespace VARIABLE
+ {
+ // Declare your Variable here:
+
+ double X1,Y1,Z1 ;
+ int N1,N2 = 0 ;
+ bool check= false ;
+
+ // A Usefull Simple Random Generator
+ TRandom Rand;
+ }
+
+using namespace VARIABLE ;
+// ----------------------------------------------------------------------------------------------
+
+
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// -----------------------------------GRAPH------------------------------------------------------
+#include <TObject.h>
+#include <TH1.h>
+#include <TH1F.h>
+#include <TH2.h>
+#include <TH2F.h>
+#include <TGraph2D.h>
+
+namespace GRAPH
+ {
+ // Declare your Spectra here:
+
+ TH1F *myHist1D = new TH1F("Hist1D","Histogramm 1D ; x ; count", 1000 , -5 , 5 ) ;
+
+ TH2F *myHist2D = new TH2F("Hist2D","Histogramm 2D ; x ; y ", 128 , 1 , 128 , 128 , 1 , 128 ) ;
+
+ }
+
+using namespace GRAPH ;
+// --------------------------------------------------------------------------------------------
+
+
+
+///////////////////////////////////////////////////////////////////////////////////////////////
+// -----------------------------------CUT------------------------------------------------------
+
+#include <TCutG.h>
+namespace CUT
+ {
+ //open the ROOT file for cuts
+ // TFile *FCuts = new TFile("cut_protons.root");
+ // TCutG *cut_protons = (TCutG*) FCuts->Get("protons");
+ }
+
+using namespace CUT ;
+// --------------------------------------------------------------------------------------------
+
+
+
+////////////////////////////////////////////////////////////////////////////////////////////////
+// -----------------------------------ENERGY LOSS----------------------------------------------
+#include "NPEnergyLoss.h"
+using namespace NPL ;
+namespace ENERGYLOSS
+ {
+ // Beam Energy Loss
+
+ EnergyLoss Fe60TargetCD2 = EnergyLoss ("Fe60[0.0]_CD2.G4table",
+ "G4Table" ,
+ 1000 );
+
+ /*
+ EnergyLoss Fe60TargetCD2 = EnergyLoss ("60Fe_CD2_new.txt",
+ "LISE" ,
+ 1000 ,
+ 3 ,
+ 59.934 );
+
+ EnergyLoss Fe60TargetCarbone= EnergyLoss ("60Fe_carbone.txt",
+ "LISE" ,
+ 1000 ,
+ 3 ,
+ 59.934 );
+ */
+ /*
+ EnergyLoss S34TargetCD2 = EnergyLoss ("S34[0.0]_CD2.G4table",
+ "G4Table" ,
+ 1000 );
+ */
+ /*
+ EnergyLoss Si34TargetCD2 = EnergyLoss ("34Si_CD2.txt",
+ "LISE" ,
+ 1000 ,
+ 1 ,
+ 34 );
+
+ EnergyLoss Si34TargetCarbone= EnergyLoss ("34Si_carbone.txt",
+ "LISE" ,
+ 1000 ,
+ 1 ,
+ 34 );
+
+ EnergyLoss S36TargetCD2 = EnergyLoss ("36S_CD2.txt",
+ "LISE" ,
+ 1000 ,
+ 1 ,
+ 36 );
+
+ EnergyLoss S36TargetCarbone= EnergyLoss ("36S_carbone.txt",
+ "LISE" ,
+ 1000 ,
+ 1 ,
+ 36 );
+
+ */
+ // proton Energy Loss
+
+ EnergyLoss protonTargetCD2 = EnergyLoss ("proton_CD2.G4table",
+ "G4Table" ,
+ 1000 );
+
+ EnergyLoss protonTargetWind = EnergyLoss ("proton_Mylar.txt",
+ "LISE" ,
+ 1000 ,
+ 3,
+ 1.008 );
+
+ EnergyLoss protonTargetC = EnergyLoss ("proton_carbone.txt" ,
+ "LISE" ,
+ 1000 ,
+ 3 ,
+ 1.008 );
+
+ /*
+ EnergyLoss protonTargetCD2 = EnergyLoss ("proton_cd2.txt" ,
+ "LISE" ,
+ 1000 ,
+ 3 ,
+ 1.008 );
+
+ */
+ EnergyLoss protonStripAl = EnergyLoss ("proton_Al.txt" ,
+ "LISE" ,
+ 1000 ,
+ 3 ,
+ 1.008 );
+
+ EnergyLoss protonPadSi = EnergyLoss ("proton_Si.txt",
+ "LISE" ,
+ 1000 ,
+ 3 ,
+ 1.008 );
+
+
+ }
+
+using namespace ENERGYLOSS ;
+// ----------------------------------------------------------------------------------------------
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// -----------------------------------Random Engine----------------------------------------------
+#include "TRandom3.h"
+namespace RANDOMENGINE
+ {
+
+ TRandom3 RandomEngine = TRandom3();
+
+ }
+
+using namespace RANDOMENGINE ;
+// ----------------------------------------------------------------------------------------------
+/////////////////////////////////////////////////////////////////////////////////////////////////
+
diff --git a/NPAnalysis/must2/macro/CrossSection.c b/NPAnalysis/must2/macro/CrossSection.c
new file mode 100644
index 0000000000000000000000000000000000000000..5640efb446f5464a98882b5a40ea015ef2745171
--- /dev/null
+++ b/NPAnalysis/must2/macro/CrossSection.c
@@ -0,0 +1,121 @@
+{
+ gROOT->Reset();
+ gStyle->SetOptStat(1);
+ gStyle->SetPalette(51,0);
+ ///////////////////////
+ ///////////////////////
+
+ double Pi = 3.141592654 ;
+
+///Usefull myAngleInRad * RadToDeg = myAngle In deg :p
+ double DegToRad = Pi/180. ; // 2Pi/360 = Pi/180
+ double RadToDeg = 180./Pi ; // 360/2Pi = 180/Pi
+
+TFile *file0 = TFile::Open("./Result/myResult.root");
+
+ cEA = new TCanvas("cEA","Kinematic Line" ,100,100,900,900);
+ hEA->Draw("COLZ");
+ cEx = new TCanvas("cEx","Excitation Energy" ,100,100,600,600);
+ hEx->Draw();
+
+ cEHexa = new TCanvas("cEHexa","Hexaneutron bound Energy" ,100,100,600,600);
+ hEHexa->Draw();
+
+ cCM = new TCanvas("cCm" , "Cross Section (CM)" , 100 , 100 , 900, 900) ;
+ hThetaCM->Draw();
+
+ vector<double> relative_error;
+ double NumberOfHitX=0;
+ ifstream efficFile;
+ efficFile.open("./17cm.efficiency");
+ string EffString, ThetaString;
+ double xxx = hThetaLabCM->GetNbinsY();
+
+ for(Int_t i = 1 ; i < hThetaLabCM->GetNbinsX() ; i++)
+ {
+ // Calculate Relative error
+ NumberOfHitX=hThetaLabCM->Integral(i, i, 1, xxx );
+ if(NumberOfHitX!=0)relative_error.push_back( 1. / ( sqrt(NumberOfHitX) ) );
+ else {relative_error.push_back(0) ;}
+ }
+
+ for(Int_t j = 1 ; j < hThetaLabCM->GetNbinsY() ; j++)
+ {
+ //Efficacity correction
+ efficFile >> ThetaString >> EffString;
+ double Eff = atoi(EffString.c_str());
+ for(Int_t i = 1 ; i < hThetaLabCM->GetNbinsX() ; i++)
+ {
+ if (Eff!=0)
+ {
+ double pipo = (hThetaLabCM->GetCellContent(i,j))/(Eff) ;
+ hThetaLabCM->SetCellContent( i, j, pipo );
+ }
+ else hThetaLabCM->SetCellContent( i, j, 0 );
+ }
+ }
+
+ TH1D* hCrossSection = new TH1D();
+ hCrossSection = hThetaLabCM->ProjectionX("Cross Section", 0, -1, "") ;
+
+ for(Int_t i = 1 ; i < hCrossSection->GetNbinsX() ; i++)
+ {
+ // Calculate Error
+ hCrossSection->SetBinContent(i,
+ NumberOfHitX=hCrossSection->GetBinContent(i)/(sin(hCrossSection->GetBinCenter(i)*DegToRad) ));
+ double error=NumberOfHitX*relative_error[i];
+ hCrossSection->SetBinError(i,error);
+ }
+
+ hCrossSection->Draw();
+
+ string Path="../Data/CrossSection/11Li(d,3He)10He.txt";
+ ifstream CrossSectionFile;
+ CrossSectionFile.open(Path.c_str());
+ if( CrossSectionFile.is_open() )
+ cout << " file " << Path << " loading " << endl;
+ else{
+ cout << " Error, no file " << Path << " found" << endl; return;}
+
+
+ //Normalisation:
+ //Int_t Maximum_Bin = hCrossSection->GetMaximumBin() ;
+ Int_t Maximum_Bin = 3 ;
+ Double_t Maximum_theta = hCrossSection->GetBinCenter(Maximum_Bin) ;
+ Double_t Bin_Width = hCrossSection->GetBinWidth(Maximum_Bin) ;
+ Double_t Maximum = hCrossSection->GetBinContent(Maximum_Bin) ;
+
+ Double_t Normalisation_Factor = 0 ;
+ Double_t ThetaVal = 0 ;
+ Double_t CSVal = 0 ;
+ TMarker marker ;
+ string theta,CS ;
+
+ while( !CrossSectionFile.eof() && ThetaVal<60)
+ {
+ CrossSectionFile >> theta >> CS ;
+
+ ThetaVal=atof(theta.c_str()) ;
+
+ if(ThetaVal>(Maximum_theta-Bin_Width/1000) && ThetaVal<(Maximum_theta+Bin_Width/1000) )
+ {
+ CSVal=atof(CS.c_str()) ;
+ Normalisation_Factor = Maximum/CSVal;
+ }
+ }
+ CrossSectionFile.close();
+ CrossSectionFile.open(Path.c_str());
+
+ //Reading Cross Section
+
+ ThetaVal=0;
+ while( !CrossSectionFile.eof() && ThetaVal<60)
+ {
+
+ CrossSectionFile >> theta >> CS ;
+ CSVal=atof(CS.c_str()) ;
+ ThetaVal=atof(theta.c_str()) ;
+ marker->DrawMarker(ThetaVal,CSVal*Normalisation_Factor ) ;
+ }
+
+}
diff --git a/NPAnalysis/must2/macro/TimeOfFlight.c b/NPAnalysis/must2/macro/TimeOfFlight.c
new file mode 100644
index 0000000000000000000000000000000000000000..916932fc04c540832f14c0c4cf14b4889d8972be
--- /dev/null
+++ b/NPAnalysis/must2/macro/TimeOfFlight.c
@@ -0,0 +1,14 @@
+{
+
+ gROOT->Reset();
+ gStyle->SetOptStat(1);
+ gStyle->SetPalette(51,0);
+ //gStyle->SetPalette(1);
+ ///////////////////////
+ ///////////////////////
+TFile *file0 = TFile::Open("./Result/myResult.root");
+
+ TCanvas* cTOF = new TCanvas("cTOF","Time of Flight" ,10,10,800,800);
+ hTOF1234->Draw("COLZ");
+
+}
diff --git a/NPAnalysis/must2/macro/affich.c b/NPAnalysis/must2/macro/affich.c
new file mode 100644
index 0000000000000000000000000000000000000000..2c8537ea62be3b91ec797d490b630dd8f4ecc6eb
--- /dev/null
+++ b/NPAnalysis/must2/macro/affich.c
@@ -0,0 +1,291 @@
+{
+
+ gROOT->Reset();
+ gStyle->SetOptStat(1);
+ //gStyle->SetPalette(51,0);
+ gStyle->SetPalette(1);
+ ///////////////////////
+ ///////////////////////
+TFile *file0 = TFile::Open("./Result/myResult.root");
+
+/* TCanvas* cTOF = new TCanvas("cTOF","Time of Flight" ,10,10,800,800);
+ cTOF->Divide(3,2);
+
+ cTOF->cd(1);
+ hTOF1->Draw("COLZ");
+
+ cTOF->cd(2);
+ hTOF2->Draw("COLZ");
+
+ cTOF->cd(3);
+ hTOF3->Draw("COLZ");
+
+ cTOF->cd(4);
+ hTOF4->Draw("COLZ");
+
+ cTOF->cd(5);
+ hTOF5->Draw("COLZ"); */
+
+ double Pi = 3.141592654 ;
+ double DegToRad = Pi/180. ; // 2Pi/360 = Pi/180
+
+/* vector<double> relative_error;
+ double NumberOfHitX=0;
+ cTheta2D = new TCanvas("cTheta2d","Cross Section" ,100,100,900,900);
+ hThetaLabCM->Draw("COLZ");*/
+ cTheta = new TCanvas("cTheta","Cross Section" ,100,100,900,900);
+ hTheta->Draw();
+/* ifstream efficFile;
+
+ efficFile.open("./efficiency/20cm2.efficiency");
+ string EffString, ThetaString;
+ double xxx = hThetaLabCM->GetNbinsY();
+
+ for(Int_t i = 1 ; i < hThetaLabCM->GetNbinsX() ; i++)
+ {
+ // Calculate Relative error
+ NumberOfHitX=hThetaLabCM->Integral(i, i, 1, xxx );
+ if(NumberOfHitX!=0)relative_error.push_back( ( sqrt(NumberOfHitX) ) / (NumberOfHitX) );
+ else {relative_error.push_back(0) ;}
+ }
+
+ for(Int_t j = 1 ; j < hThetaLabCM->GetNbinsY() ; j++)
+ {
+ //Efficacity correction
+ efficFile >> ThetaString >> EffString;
+ double Eff = atoi(EffString.c_str());
+ for(Int_t i = 1 ; i < hThetaLabCM->GetNbinsX() ; i++)
+ {
+ if (Eff!=0)
+ {
+ double pipo = (hThetaLabCM->GetCellContent(i,j))/(Eff/100) ;
+ hThetaLabCM->SetCellContent( i, j, pipo );
+ }
+ else hThetaLabCM->SetCellContent( i, j, 0 );
+ }
+ }
+
+
+ TH1D* hCrossSection = new TH1D();
+ hCrossSection = hThetaLabCM->ProjectionX("Cross Section", 0, -1, "") ;
+
+
+ for(Int_t i = 1 ; i < hCrossSection->GetNbinsX() ; i++)
+ {
+ // Calculate Error
+ NumberOfHitX=hCrossSection->GetBinContent(i);
+ double error=NumberOfHitX*relative_error[i];
+ hCrossSection->SetBinError(i,error);
+ }
+
+ hCrossSection->Draw();
+
+/* ofstream efficFile;
+ efficFile.open("20cm.efficiency");
+
+ for(Int_t i = 1 ; i < hTheta->GetNbinsX() ; i++)
+ {
+
+ Double_t Flux = 100000*sin(hTheta->GetBinCenter(i)*DegToRad)*2*DegToRad/2;
+
+ if(Flux!=0)
+ {
+ Double_t Efficiency = ( ( Flux - hTheta->GetBinContent(i) ) / Flux ) * 100 ;
+
+ hTheta->SetBinContent(i, 100-Efficiency ) ;
+
+ efficFile << hTheta->GetBinCenter(i) << " " << 100-Efficiency << endl;
+ }
+
+ else hTheta->SetBinContent(i, 0 );
+ }
+
+
+ //hTheta->Draw();*/
+
+ //efficiency correction and error bar calculation
+/* ifstream efficFile;
+ efficFile.open("./efficiency/20cmCM.efficiency");
+ string EffString, ThetaString;
+ double Eff, relative_error, error, toto ;
+
+ for(Int_t i = 1 ; i < hThetaCM->GetNbinsX() ; i++)
+ {
+
+ efficFile >> ThetaString >> EffString;
+ Eff = atoi(EffString.c_str());
+ if(hThetaCM->GetBinContent(i)!=0)relative_error = sqrt(hThetaCM->GetBinContent(i))/(hThetaCM->GetBinContent(i));
+
+ toto=hThetaCM->GetBinContent(i)*(100-Eff) ;
+ hThetaCM->SetBinContent(i, toto) ;
+
+ error=hThetaCM->GetBinContent(i)*relative_error;
+ hThetaCM->SetBinError(i,error);
+ }*/
+
+
+
+
+ //hThetaCM->Draw();
+
+/* string Path="/home/Adrien/Desktop/geant/Simulation/Data/CrossSection/11Li(d,3He)10He.txt";
+ ifstream CrossSectionFile;
+ CrossSectionFile.open(Path.c_str());
+ if( CrossSectionFile.is_open() )
+ cout << " file " << Path << " loading " << endl;
+ else{
+ cout << " Error, no file " << Path << " found" << endl; return;}
+
+
+ //Normalisation:
+ Int_t Maximum_Bin = hCrossSection->GetMaximumBin() ;
+ Double_t Maximum_theta = hCrossSection->GetBinCenter(Maximum_Bin) ;
+ Double_t Bin_Width = hCrossSection->GetBinWidth(Maximum_Bin) ;
+ Double_t Maximum = hCrossSection->GetBinContent(Maximum_Bin) ;
+
+ Double_t Normalisation_Factor=0;
+
+
+
+ Double_t ThetaVal=0 ;
+ Double_t CSVal=0 ;
+ TMarker marker ;
+ string theta,CS ;
+
+ while( !CrossSectionFile.eof() && ThetaVal<60)
+ {
+ CrossSectionFile >> theta >> CS ;
+
+ ThetaVal=atof(theta.c_str()) ;
+
+ if(ThetaVal>(Maximum_theta-Bin_Width/4) && ThetaVal<(Maximum_theta+Bin_Width/4) )
+ {
+ CSVal=atof(CS.c_str()) ;
+ Normalisation_Factor = Maximum/CSVal;
+ }
+ }
+ CrossSectionFile.close();
+ CrossSectionFile.open(Path.c_str());
+
+ //Reading Cross Section
+
+ ThetaVal=0;
+ while( !CrossSectionFile.eof() && ThetaVal<60)
+ {
+
+ CrossSectionFile >> theta >> CS ;
+ CSVal=atof(CS.c_str()) ;
+ ThetaVal=atof(theta.c_str()) ;
+ marker->DrawMarker(ThetaVal,CSVal*Normalisation_Factor) ;
+ }
+
+
+
+/*
+ //Normalisation:
+ Double_t Maximum_Bin = hTheta->GetMaximumBin() ;
+ Double_t Maximum = hTheta->GetBinContent(Maximum_Bin) ;
+
+ Double_t Normalisation_Factor = Maximum/3.;
+
+
+
+ //Reading Cross Section
+ Double_t ThetaVal=0 ;
+ Double_t CSVal=0 ;
+ TMarker marker ;
+ string theta,CS ;
+
+ while( !CrossSectionFile.eof() && ThetaVal<60)
+ {
+
+ CrossSectionFile >> theta >> CS ;
+ CSVal=atof(CS.c_str()) ;
+ ThetaVal=atof(theta.c_str()) ;
+
+ ThetaVal=0.000479597*pow(ThetaVal,3)-0.0670771*pow(ThetaVal,2)+3.08149*ThetaVal-0.339958 ;
+
+ marker->DrawMarker(ThetaVal,CSVal*Normalisation_Factor) ;
+ }
+// cTheta->SetLogy(1);
+
+
+
+
+ /*
+ cKine = new TCanvas("cKine","Kinematics lines" ,10,10,800,800);
+ hKine->Draw("COLZ");*/
+
+/* cKineInit = new TCanvas("cKineInit","Initial Kinematics lines" ,100,100,600,600);
+ hKineInit->Draw("COLZ");
+
+ cEDE = new TCanvas("EDE","EDE add Strip" ,100,100,600,600);
+ hEDEAddStrip->Draw("COLZ");
+/* cG = new TCanvas("cG","Strip position",500,100,800,600);
+ cG->Divide(3,2);
+
+ cG->cd(1);
+ Agraph2D->Draw("P0") ;
+ cG->cd(2);
+ Bgraph2D->Draw("P0") ;
+ cG->cd(3);
+ Cgraph2D->Draw("P0") ;
+ cG->cd(4);
+ Dgraph2D->Draw("P0") ;
+ cG->cd(5);
+ Egraph2D->Draw("P0") ;
+
+ cG2 = new TCanvas("cG2","all Strip position",500,100,800,600);
+ TOTgraph2D->Draw("P0") ;
+ */
+ cH = new TCanvas("cH","Hit density",500,100,1000,800);
+ cH->Divide(3,2);
+
+ cH->cd(1);
+ hHIT4->Draw("COLZ");
+ cH->cd(2);
+ hHIT2->Draw("COLZ");
+ cH->cd(3);
+ hHIT5->Draw("COLZ");
+ cH->cd(4);
+ hHIT1->Draw("COLZ");
+ cH->cd(5);
+ hHIT3->Draw("COLZ");
+ cH->cd(6);
+ hThetaHeavy->Draw();
+
+/* cH->cd(5);
+ hXZ->Draw("COLZ");
+ cH->cd(6);
+ hXY->Draw("COLZ");
+
+ cEDE = new TCanvas("cEDE","EDE indentification Spectra",500,100,1000,800);
+ cEDE->Divide(3,2);
+
+ cEDE->cd(1);
+ hEDE1->Draw("COLZ");
+ cEDE->cd(2);
+ hEDE2->Draw("COLZ");
+ cEDE->cd(3);
+ hEDE3->Draw("COLZ");
+ cEDE->cd(4);
+ hEDE4->Draw("COLZ");
+ cEDE->cd(5);
+ hEDE5->Draw("COLZ");
+ cEDE->cd(6);
+ TH2F* hEDET= new TH2F("hEDET","",4000,-1,600, 400, -1, 25) ;
+ hEDET->Add(hEDE1);
+ hEDET->Add(hEDE2);
+ hEDET->Add(hEDE3);
+ hEDET->Add(hEDE4);
+ hEDET->Draw("CLOZ");*/
+
+ cEx = new TCanvas("cEx","Excitation Energy" ,100,100,600,600);
+ hEx->Draw();
+ /*cE = new TCanvas("cE","Light Energy" ,500,100,1000,800);
+ cE->Divide(2);
+ cE->cd(1);
+ hE1234->Draw();
+ cE->cd(2);
+ hE5->Draw();*/
+}
diff --git a/NPAnalysis/must2/macro/efficiency.c b/NPAnalysis/must2/macro/efficiency.c
new file mode 100644
index 0000000000000000000000000000000000000000..9a04af185f75853d0bce5bae2009cd8cd01cae1a
--- /dev/null
+++ b/NPAnalysis/must2/macro/efficiency.c
@@ -0,0 +1,58 @@
+{
+ gROOT->Reset();
+ gStyle->SetOptStat(1);
+ gStyle->SetPalette(51,0);
+ ///////////////////////
+ ///////////////////////
+TFile *file0 = TFile::Open("./Result/myResult.root");
+
+ double Pi = 3.141592654 ;
+ double DegToRad = Pi/180. ; // 2Pi/360 = Pi/180
+
+
+ ofstream efficFile;
+ efficFile.open("12cm.efficiency");
+
+ for(Int_t i = 1 ; i < hTheta->GetNbinsX() ; i++)
+ {
+ Double_t Flux = 1000000./2. * sin(hTheta->GetBinCenter(i)*DegToRad) * (2*DegToRad) ;
+
+ if(Flux!=0)
+ {
+ Double_t Efficiency = 100*( (hTheta->GetBinContent(i)) / Flux );
+ hTheta->SetBinContent(i, Efficiency ) ;
+ }
+
+ else hTheta->SetBinContent(i, 0 );
+
+ efficFile << hTheta->GetBinCenter(i) << " " << Efficiency << endl;
+ }
+
+ cEff = new TCanvas("cEff","Efficiency" ,100,100,600,600);
+ hTheta->Draw();
+
+ cHit = new TCanvas("cHit","Hit" ,100,100,600,600);
+ hXY->Draw("COLZ");
+
+ cHit2 = new TCanvas("cHit2","Hit" ,100,100,600,600);
+ cHit2->Divide(3,2);
+
+ cHit2->cd(1);
+ hXY1->Draw("COLZ");
+
+ cHit2->cd(2);
+ hXY2->Draw("COLZ");
+
+ cHit2->cd(3);
+ hXY3->Draw("COLZ");
+
+ cHit2->cd(4);
+ hXY4->Draw("COLZ");
+
+ cHit2->cd(5);
+ hXY5->Draw("COLZ");
+
+ cHit2->cd(6);
+ hXY6->Draw("COLZ");
+
+}
diff --git a/NPAnalysis/must2/macro/hit.c b/NPAnalysis/must2/macro/hit.c
new file mode 100644
index 0000000000000000000000000000000000000000..50db19da5e7c54b1cba977ffe34232fb328dcdd0
--- /dev/null
+++ b/NPAnalysis/must2/macro/hit.c
@@ -0,0 +1,44 @@
+{
+ gROOT->Reset();
+ gStyle->SetOptStat(1);
+ gStyle->SetPalette(51,0);
+ ///////////////////////
+ ///////////////////////
+TFile *file0 = TFile::Open("./Result/myResult.root");
+
+ cEff = new TCanvas("cEff","Theta Distribution" ,100,100,600,600);
+ hTheta->Draw();
+
+
+ cCM = new TCanvas("cCm" , "ThetaCM" , 100 , 100 , 600, 600) ;
+ hThetaCM->Draw();
+
+ cHit = new TCanvas("cHit","Hit" ,100,100,600,600);
+ hXY->Draw("COLZ");
+
+
+ cHit2 = new TCanvas("cHit2","Hit" ,100,100,600,600);
+ cHit2->Divide(3,2);
+
+ cHit2->cd(1);
+ hXY1->Draw("COLZ");
+
+ cHit2->cd(2);
+ hXY2->Draw("COLZ");
+
+ cHit2->cd(3);
+ hXY3->Draw("COLZ");
+
+ cHit2->cd(4);
+ hXY4->Draw("COLZ");
+
+ cHit2->cd(5);
+ hXY5->Draw("COLZ");
+
+ cHit2->cd(6);
+ hXY6->Draw("COLZ");
+
+ cEx = new TCanvas("cEx","Excitation Energy" ,100,100,300,300);
+ hEx->Draw();
+
+}
diff --git a/NPAnalysis/must2/src/Analysis.cc b/NPAnalysis/must2/src/Analysis.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1cefa90e30428e969a44f82fb98b459d7b169f4b
--- /dev/null
+++ b/NPAnalysis/must2/src/Analysis.cc
@@ -0,0 +1,205 @@
+#include "ObjectManager.hh"
+using namespace std;
+
+#define MM_DEADLAYER 0.4 // Must2 aluminium dead layer in um
+
+double ThetaCalculation (TVector3 A , TVector3 B);
+double PhiCalculation (TVector3 A , TVector3 B);
+
+int main(int argc,char** argv)
+{
+ // command line parsing
+ NPOptionManager* myOptionManager = NPOptionManager::getInstance(argc,argv);
+ string reactionfileName = myOptionManager->GetReactionFilePath();
+ string detectorfileName = myOptionManager->GetDetectorFilePath();
+ string calibrationfileName = myOptionManager->GetCalibrationFilePath();
+ string runToReadfileName = myOptionManager->GetRunToReadFilePath();
+ string OutputfileName = myOptionManager->GetOutputFilePath();
+
+ /////////////////////////////////////////////////////////////////////////////////////////////////////
+ // First of All instantiate RootInput and Output
+ // Detector will be attached later
+ RootInput:: getInstance(runToReadfileName);
+ RootOutput::getInstance("Analysis/Must2_AnalysedData", "AnalysedTree");
+
+ // Instantiate some Reaction
+ cout << endl << "/////////// Event generator ///////////" << endl;
+ NPL::Reaction* myReaction = new Reaction();
+ myReaction->ReadConfigurationFile(reactionfileName);
+
+ // Instantiate the Calibration Manger using a file (WARNING:prior to the detector instantiation)
+ CalibrationManager* myCalibration = CalibrationManager::getInstance(calibrationfileName) ;
+
+ // Instantiate the detector using a file
+ NPA::DetectorManager* myDetector = new DetectorManager();
+ myDetector->ReadConfigurationFile(detectorfileName);
+
+ // Calculate beam energy at target middle
+ // Target informations
+ cout << endl;
+ cout << "/////////// Target information ///////////" << endl;
+ cout << "Thickness (um): " << myDetector->GetTargetThickness() << endl;
+ // Get nominal beam energy
+ Double_t BeamEnergyNominal = myReaction->GetBeamEnergy() * MeV;
+ cout << "Nominal beam energy (MeV): " << BeamEnergyNominal << endl;
+ // Slow beam at target middle
+ Double_t BeamEnergy = Fe60TargetCD2.Slow(BeamEnergyNominal, myDetector->GetTargetThickness()/2 * micrometer, 0);
+ cout << "Middle-target beam energy (MeV): " << BeamEnergy << endl;
+ // Set energy beam at target middle
+ myReaction->SetBeamEnergy(BeamEnergy);
+ myReaction->Print();
+
+ // Ask the detector manager to load the parameter added by the detector in the calibrationfileName
+ myCalibration->LoadParameterFromFile() ;
+
+ // Attach more branch to the output
+ const int mult = 2;
+ double ELab[mult], ThetaLab[mult], ThetaCM[mult], PhiLab[mult], ExcitationEnergy[mult];
+
+ // some initializations
+ for (int i = 0; i < mult; i++) {
+ ELab[i] = -10000; ThetaLab[i] = -10000; ThetaCM[i] = -10000; ExcitationEnergy[i] = -10000; PhiLab[i] = -10000;
+ }
+
+ RootOutput::getInstance()->GetTree()->Branch("ELab", &ELab, "ELab[2]/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaLab", &ThetaLab, "ThetaLab[2]/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaCM", &ThetaCM, "ThetaCM[2]/D");
+ RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy", &ExcitationEnergy, "ExcitationEnergy[2]/D");
+
+ // Get the formed Chained Tree and Treat it
+ TChain* Chain = RootInput::getInstance()->GetChain();
+
+ // Connect TInitialConditions branch
+ TInitialConditions *initCond = 0;
+ Chain->SetBranchAddress("InitialConditions", &initCond);
+ Chain->SetBranchStatus("InitialConditions", 1);
+
+ // Get TMust2Physics pointer
+ TMust2Physics *M2 = (TMust2Physics*) myDetector -> m_Detector["MUST2"] ;
+
+ // define user spectra
+ TH2F* DE_E_protons = new TH2F("DE_E_protons", "DE_E et cut protons", 1000, 0, 25000, 1000, 0, 25000);
+ TH2F* DE_E = new TH2F("DE_E", "DE_E", 1000, 0, 25000, 1000, 0, 25000);
+ TH2F* E_Theta = new TH2F("E_Theta", "E_Theta", 60, 100, 180, 500, 0, 20);
+ TH2F* Position_M2 = new TH2F("Must2Positions", "Must2Positions", 2000, -200, 200, 2000, -200, 200);
+ TH2F* Position_M2_tot = new TH2F("Must2Positions_tot", "Must2Positions_tot", 130, -1, 129, 130, -1, 130);
+
+ cout << " ///////// Starting Analysis ///////// "<< endl << endl ;
+ int N = Chain -> GetEntries();
+ cout << " Number of Event to be treated : " << N << endl ;
+ clock_t begin=clock();
+ clock_t end=begin;
+
+ for (int i = 0; i < N; i++) {
+ if (i%10000 == 0 && i!=0) {
+ cout.precision(5);
+ end = clock();
+ double TimeElapsed = (end-begin) / CLOCKS_PER_SEC;
+ double percent = (double)i/N ;
+ double TimeToWait = (TimeElapsed/percent) - TimeElapsed;
+ cout << "\r Progression:" << percent*100 << " % \t | \t Remaining time : ~" << TimeToWait <<"s"<< flush;
+ }
+ else if (i==N-1) cout << "\r Progression:" << " 100% " <<endl;
+
+ // Get raw data
+ Chain->GetEntry(i);
+
+ // Clear Previous Event
+ myDetector->ClearEventPhysics();
+
+ // Build the new event
+ myDetector->BuildPhysicalEvent();
+
+ // Get target interaction position from initial conditions
+ double XTarget = initCond->GetICPositionX(0);
+ double YTarget = initCond->GetICPositionY(0);
+ double ZTarget = initCond->GetICPositionZ(0);
+
+ // Calculate beam direction
+ double BeamTheta = initCond->GetICIncidentAngleTheta(0)*deg; double BeamPhi = initCond->GetICIncidentAnglePhi(0)*deg;
+ TVector3 BeamDirection = TVector3(cos(BeamPhi)*sin(BeamTheta) , sin(BeamPhi)*sin(BeamTheta) , cos(BeamTheta));
+
+ for (int hit = 0; hit < M2->Si_E.size(); hit++) {
+ ExcitationEnergy[hit] = -10000;
+ ELab[hit] = -10000 ; ThetaLab[hit] = -10000;
+
+ // Get hit direction
+ TVector3 HitDirection = M2->GetPositionOfInteraction(hit) - TVector3(XTarget, YTarget, myDetector->GetTargetZ());
+
+ // Angle between beam and particle
+ ThetaLab[hit] = ThetaCalculation(HitDirection , BeamDirection) * rad;
+ PhiLab[hit] = (M2->GetPositionOfInteraction(hit)).Phi();
+
+ // Angle between particule and z axis (target Normal)
+ double ThetaN = ThetaCalculation(HitDirection , TVector3(0,0,1));
+ // Angle between particule and Must2 Si surface
+ double ThetaMM2Surface = ThetaCalculation(HitDirection, M2->GetTelescopeNormal(hit));
+
+ if ((M2->Si_E[hit] > 0)) {
+ Position_M2->Fill((M2->GetPositionOfInteraction(hit)).X(), (M2->GetPositionOfInteraction(hit)).Y());
+ }
+
+ if ((M2->Si_E[hit] > 0) && (M2->SiLi_E[hit]) < 0) {
+ ELab[hit] = M2->Si_E[hit] * MeV;
+ ELab[hit] = protonStripAl.EvaluateInitialEnergy(ELab[hit], MM_DEADLAYER*micrometer, ThetaMM2Surface) /MeV;
+ ELab[hit] = protonTargetCD2.EvaluateInitialEnergy(ELab[hit], myDetector -> GetTargetThickness() /2 *micrometer , ThetaLab[hit]) /MeV;
+
+ ExcitationEnergy[hit] = myReaction->ReconstructRelativistic(ELab[hit]/MeV, ThetaLab[hit]) /MeV;
+
+ ThetaLab[hit] = ThetaLab[hit] /deg;
+ PhiLab[hit] = PhiLab[hit] /deg;
+ }
+
+ else if ((M2->Si_E[hit] > 0) && (M2->SiLi_E[hit]) > 0) {
+ Position_M2_tot->Fill(M2->Si_X[hit], M2->Si_Y[hit]);
+ DE_E ->Fill(M2->SiLi_E[hit], M2->Si_E[hit]);
+
+ ELab[hit] = M2->SiLi_E[hit]*MeV;
+ ELab[hit] = protonStripAl.EvaluateInitialEnergy(ELab[hit]/MeV, MM_DEADLAYER*micrometer, ThetaMM2Surface) /MeV;
+ ELab[hit] += M2->Si_E[hit]*MeV;
+ ELab[hit] = protonStripAl.EvaluateInitialEnergy(ELab[hit]/MeV, MM_DEADLAYER*micrometer, ThetaMM2Surface) /MeV;
+ ELab[hit] = protonTargetCD2.EvaluateInitialEnergy(ELab[hit]/MeV,myDetector -> GetTargetThickness() /2 *micrometer,ThetaLab[hit]) /MeV;
+
+ ExcitationEnergy[hit] = myReaction -> ReconstructRelativistic( ELab[hit]/MeV, ThetaLab[hit]) /MeV ;
+ ThetaLab[hit] = ThetaLab[hit] /deg;
+ PhiLab[hit] = PhiLab[hit] /deg;
+ }
+
+ E_Theta->Fill(ThetaLab[hit],ELab[hit]/MeV);
+
+ } // end of for loop over hit
+
+ RootOutput::getInstance()->GetTree()->Fill();
+ } // end of for loop over events
+
+
+ cout << " A total of " << N << " event has been analysed " << endl ;
+ cout << endl << " ///////////////////////////////////// "<< endl<< endl ;
+ RootOutput::getInstance()->Destroy();
+
+ return 0;
+}
+
+
+
+double ThetaCalculation (TVector3 A , TVector3 B)
+{
+ double Theta = acos( (A.Dot(B)) / (A.Mag()*B.Mag()) ) ;
+
+ return Theta*rad;
+}
+
+
+
+double PhiCalculation (TVector3 A, TVector3 B)
+{
+ double Theta = acos((A.Dot(B)) / (A.Mag()*B.Mag()));
+ double R = A.Mag();
+
+ double Phi = 0;
+ double test = A.X() / R * sin(Theta);
+
+ Phi = acos(test);
+
+ return Phi*rad;
+}
diff --git a/NPAnalysis/must2/src/GNUmakefile b/NPAnalysis/must2/src/GNUmakefile
new file mode 100644
index 0000000000000000000000000000000000000000..28c404622cec4d40b2e93aa0009fc4fc4fa97d26
--- /dev/null
+++ b/NPAnalysis/must2/src/GNUmakefile
@@ -0,0 +1,42 @@
+###Make file made by Adrien MATTA/ Institut de Physique Nucleaire d'Orsay IPNO###
+# Made to compile the ROOT Analyser for MUST2 experiment
+
+CPP=g++
+EXEC=Analysis
+
+# local includes
+NPAINCLUDES = ../include
+
+# ROOT includes
+CXXFLAGS += `root-config --cflags`
+
+# CLHEP includes
+CXXFLAGS += -I$(CLHEP_INCLUDE_DIR)
+CXXFLAGS += -I$(NPAINCLUDES)
+CXXFLAGS += -I$(NPLIB)/include
+
+LDFLAGS = `root-config --libs` -lMathMore
+LDFLAGS+= `$(NPLIB)/liblist`
+LDFLAGS+= -L$(CLHEP_LIB_DIR) -l$(CLHEP_LIB)
+
+SRC= $(wildcard *.cc)
+INC= $(wildcard $(NPAINCLUDES)/*.hh)
+OBJ=$(SRC:.cc=.o)
+
+#all:$(EXEC)
+# @$(CPP) -o $@ -c $< $(CXXFLAGS)
+
+Analysis:$(OBJ) $(INC)
+ @$(CPP) $(CXXFLAGS) -o $@ $^ $(LDFLAGS)
+ mv Analysis ../Analysis
+
+%.o: %.cc
+ @$(CPP) $(CXXFLAGS) -o $@ -c $<
+
+.PHONY: clean mrproper
+
+clean:
+ rm -rf *.o
+
+mrproper: clean
+ rm -rf $(EXEC)