No commit message

No commit message

NPAnalysis/10He_Riken/Makefile

+
+Analyse:
+	make -C ./src
+	
+clean:
+	make clean -C ./src

NPAnalysis/10He_Riken/RunToTreat.txt

+TTreeName 
+	SimulatedTree
+RootFileName 
+	../../Outputs/Simulation/mySimul.root

NPAnalysis/10He_Riken/include/DetectorManager.hh

+#ifndef DetectorManager_h
+#define DetectorManager_h
+
+//	NPL
+#include "VDetector.h"
+
+//	STL
+#include <string>
+#include <map>
+
+using namespace std ;
+using namespace NPA ;
+
+// This class manage a map of virtual detector
+
+class DetectorManager
+	{
+		public:
+			DetectorManager()		;
+			~DetectorManager()		;
+			
+		public:
+			//	Read stream at Path and pick-up Token declaration of Detector
+			void	ReadConfigurationFile(string Path)	;
+			void	BuildPhysicalEvent()				;
+			void	BuildSimplePhysicalEvent()			;
+			void	InitializeRootInput()				;
+			void	InitializeRootOutput()				;
+			void	AddDetector(string,VDetector*)		;
+			void	ClearEventPhysics()					;
+			void	ClearEventData()					;
+
+		public:	//	The map containning all detectors
+			//	Using a Map one can access to any detector using its name
+			map<string,VDetector*>	m_Detector	;
+				
+	};
+
+#endif 

NPAnalysis/10He_Riken/include/ObjectManager.hh

+// 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.hh"
+#include "Must2Array.h"
+
+// STL C++
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <cmath>
+#include <cstdlib>
+
+// ROOT
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TLeaf.h>
+#include <TVector3.h>
+#include <TRandom.h>
+
+// NPL
+#include "TMust2Data.h"
+#include "TMust2Physics.h"
+#include "NPReaction.h"
+#include "RootInput.h"
+#include "RootOutput.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
+	{
+		//	Declare your Cut here:
+
+	}
+
+using namespace CUT ;
+// --------------------------------------------------------------------------------------------
+
+
+
+////////////////////////////////////////////////////////////////////////////////////////////////
+// -----------------------------------ENERGY LOSS----------------------------------------------
+#include "NPEnergyLoss.h"
+using namespace NPL ;
+namespace ENERGYLOSS
+	{
+	
+		//	Declare your Energy loss here	:
+	/*		EnergyLoss ProtonTarget = EnergyLoss 	(	"CD2.txt" 	,
+														100 		,
+														1			);
+	*/
+	}
+	
+using namespace ENERGYLOSS ;
+// ----------------------------------------------------------------------------------------------
+/////////////////////////////////////////////////////////////////////////////////////////////////
+
+

NPAnalysis/10He_Riken/macro/CrossSection.c

+{
+	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	)	;
+	}
+	
+}

NPAnalysis/10He_Riken/macro/TimeOfFlight.c

+{
+
+	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"); 
+	
+}

NPAnalysis/10He_Riken/macro/affich.c

+{
+
+	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();*/
+}

NPAnalysis/10He_Riken/macro/efficiency.c

+{
+	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");
+	
+}

NPAnalysis/10He_Riken/macro/hit.c

+{
+	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();
+	
+}

NPAnalysis/10He_Riken/src/Analysis.cc

+#include "ObjectManager.hh"
+
+using namespace std;
+
+int main(int argc,char** argv)
+{	
+	
+	if(argc!=4) 
+		{
+			cout << 
+			"you need to specify both a Reaction file and a Detector file such as : Analysis 	myReaction.reaction myDetector.detector runToRead.run" 
+			<< endl;
+			return 0;
+		}
+	
+	string reactionfileName 	= argv[1]	;
+	string detectorfileName 	= argv[2]	;
+	string runToReadfileName 	= argv[3]	;
+	
+	//	First of All instantiate RootInput and Output
+	//	Detector will be attached later
+	RootInput:: getInstance(runToReadfileName)	;
+	RootOutput::getInstance("Analysis/10HeRiken_AnalyzedData", "AnalyzedTree")					;
+	
+	//	Instantiate a Reaction
+	NPL::Reaction* myReaction = new Reaction					;
+	myReaction	->	ReadConfigurationFile(reactionfileName)		;
+
+		//	Instantiate the detector using a file 
+	DetectorManager* myDetector = new DetectorManager 			;
+	myDetector	->	ReadConfigurationFile(detectorfileName)		;
+	
+	//	Attach more branch to the output
+	double Ex = 0 ; double EE = 0 ; double TT = 0 ; double X = 0 ; double Y = 0 ; int det ;
+	RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy",&Ex,"Ex/D") ;
+	RootOutput::getInstance()->GetTree()->Branch("E",&EE,"EE/D") ;
+	RootOutput::getInstance()->GetTree()->Branch("A",&TT,"TT/D") ;
+	RootOutput::getInstance()->GetTree()->Branch("X",&X,"X/D") ;
+	RootOutput::getInstance()->GetTree()->Branch("Y",&Y,"Y/D") ;
+	
+	// Get Must2 Pointer:
+	MUST2Array* M2 = (MUST2Array*) myDetector -> m_Detector["MUST2"] ;
+	//	Get the formed Chained Tree and Treat it
+	TChain* Chain = RootInput:: getInstance() -> GetChain()	;
+	int i;	
+	for ( i = 0 ; i < Chain -> GetEntries() ; i ++ )
+		{
+			if( i%10000 == 0 && i!=0) cout << i << " Event annalysed " << endl ;						
+			Chain -> GetEntry(i);
+			
+			myDetector -> ClearEventPhysics()				;
+			myDetector -> BuildPhysicalEvent()				;
+			
+			
+			double E = M2 -> GetEnergyDeposit();
+			TVector3 A = M2 -> GetPositionOfInteraction();
+			double Theta = ThetaCalculation ( A , TVector3(0,0,1) ) ;
+			if(E>-1000)		Ex = myReaction -> ReconstructRelativistic( E , Theta ) ;
+			else Ex = -100 ;
+			EE = E ; TT = Theta/deg ;
+			if(E>-1000){
+			X = A . X();
+			Y = A . Y();}
+			else{X = -1000 ; Y = -1000;}
+			
+			RootOutput::getInstance()->GetTree()->Fill()	;
+		}
+	cout << "A total of " << i << " event has been annalysed " << endl ;
+	
+	RootOutput::getInstance()->Destroy();
+	return 0	;
+}
+
+
+double ThetaCalculation (TVector3 A , TVector3 B)
+	{
+		double Theta = acos( (A.Dot(B)) / (A.Mag()*B.Mag()) ) ;
+		return Theta ;
+	}
+

NPAnalysis/10He_Riken/src/DetectorManager.cc

+#include "DetectorManager.hh"
+
+//	STL
+#include <iostream>
+#include <fstream>
+#include <cstdlib>
+
+//	Detector	
+#include "Must2Array.h"
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+//	Default Constructor and Destructor
+DetectorManager::DetectorManager()	
+	{}
+	
+/////////////////////////////////////////////////////////////////////////////////////////////////	
+DetectorManager::~DetectorManager()
+	{}
+
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+//	Read stream at ConfigFile and pick-up Token declaration of Detector
+void	DetectorManager::ReadConfigurationFile(string Path)	
+	{
+		////////General Reading needs////////
+	   string LineBuffer;
+	   string DataBuffer;
+
+	   /////////Boolean////////////////////
+	   bool MUST2           = false  ;
+	   bool AddThinSi       = false  ;
+	   bool GeneralTarget   = false  ;
+	   bool GPDTracker      = false  ;	// Gaspard Tracker
+	   //////////////////////////////////////////////////////////////////////////////////////////
+	   // added by Nicolas [07/05/09]
+	   string GlobalPath = getenv("NPTOOL");
+	   Path = GlobalPath + "/Inputs/DetectorConfiguration/" + Path;
+	   ifstream ConfigFile;
+	   ConfigFile.open(Path.c_str());
+
+	   if (ConfigFile.is_open())
+	   	{
+	   		cout << "/////////////////////////////" << endl;
+	      	cout << " Configuration file " << Path << " loading " << endl;
+	   	}	
+	   	
+	   else 
+	   	{
+	      cout << " Error, no configuration file" << Path << " found" << endl;
+	      return;
+	   	}
+
+
+	   while (!ConfigFile.eof()) {
+	      //Pick-up next line
+	      getline(ConfigFile, LineBuffer);
+	      //Search for comment Symbol: %
+	      if (LineBuffer.compare(0, 1, "%") == 0) {   /*Do  Nothing*/ ;}
+
+	    /*  ////////////////////////////////////////////
+	      //////////// Search for Gaspard ////////////
+	      ////////////////////////////////////////////
+	      else if (LineBuffer.compare(0, 14, "GaspardTracker") == 0 && GPDTracker == false) {
+	         GPDTracker = true ;
+	         cout << "//////// Gaspard Tracker ////////" << endl   ;
+
+	         // Instantiate the new array as a VDetector Object
+	         VDetector* myDetector = new GaspardTracker()                  ;
+
+	         // Read Position of Telescope
+	         ConfigFile.close()                                 ;
+	         myDetector->ReadConfiguration(Path)                   ;
+	         ConfigFile.open(Path.c_str())                      ;
+
+	         // Add array to the VDetector Vector
+	         AddDetector(myDetector)                            ;
+	      }*/
+
+	      ////////////////////////////////////////////
+	      //////// Search for MUST2 Array  ////////
+	      ////////////////////////////////////////////
+	      else if (LineBuffer.compare(0, 10, "MUST2Array") == 0 && MUST2 == false) {
+	         MUST2 = true ;
+	         cout << "//////// MUST2 Array ////////" << endl   << endl   ;
+
+	         // Instantiate the new array as a VDetector Object
+	         VDetector* myDetector = new MUST2Array()                 ;
+
+	         // Read Position of Telescope
+	         ConfigFile.close()                                    ;
+	         myDetector->ReadConfiguration(Path)                      ;
+	         ConfigFile.open(Path.c_str())                         ;
+
+	         // Add array to the VDetector Vector
+	         AddDetector("MUST2" , myDetector)                               ;
+	      }
+
+	/*      ////////////////////////////////////////////
+	      ////////// Search for Add.ThinSi ///////////
+	      ////////////////////////////////////////////
+	      else if (LineBuffer.compare(0, 9, "AddThinSi") == 0 && AddThinSi == false) {
+	         AddThinSi = true ;
+	         cout << "//////// Thin Si ////////" << endl << endl   ;
+
+	         // Instantiate the new array as a VDetector Object
+	         VDetector* myDetector = new ThinSi()                  ;
+
+	         // Read Position of Telescope
+	         ConfigFile.close()                                 ;
+	         myDetector->ReadConfiguration(Path)                   ;
+	         ConfigFile.open(Path.c_str())                      ;
+
+	         // Add array to the VDetector Vector
+	         AddDetector(myDetector)                            ;
+	      }
+
+	      ////////////////////////////////////////////
+	      //////////// Search for Target /////////////
+	      ////////////////////////////////////////////
+
+	      else if (LineBuffer.compare(0, 13, "GeneralTarget") == 0 && GeneralTarget == false) {
+	         GeneralTarget = true ;
+	         cout << "////////// Target ///////////" << endl   << endl   ;
+
+	         // Instantiate the new array as a VDetector Objects
+	         VDetector* myDetector = new Target()                     ;
+
+	         // Read Position and target specification
+	         ConfigFile.close()                                    ;
+	         myDetector->ReadConfiguration(Path)                      ;
+	         ConfigFile.open(Path.c_str())                         ;
+
+	         m_TargetThickness = ((Target*)myDetector)->GetTargetThickness()   ;
+	         m_TargetRadius    = ((Target*)myDetector)->GetTargetRadius()   ;
+	         m_TargetX        = ((Target*)myDetector)->GetTargetX()         ;
+	         m_TargetY        = ((Target*)myDetector)->GetTargetY()         ;
+	         m_TargetZ        = ((Target*)myDetector)->GetTargetZ()         ;
+
+	         // Add target to the VDetector Vector
+	         AddDetector(myDetector)                               ;
+	      }*/
+
+	      //Nothing understandable
+	      //else ;
+	   }
+
+		ConfigFile.close();
+		
+		InitializeRootInput();
+		InitializeRootOutput();
+		
+		return   ;
+	}
+	
+/////////////////////////////////////////////////////////////////////////////////////////////////	
+
+void	DetectorManager::BuildPhysicalEvent()						
+	{	
+		map<string,VDetector*>::iterator it	;
+		
+		for( it = m_Detector.begin() ; it != m_Detector.end() ; ++it) 
+			{
+				it->second->BuildPhysicalEvent()	;	
+			}           
+	}
+/////////////////////////////////////////////////////////////////////////////////////////////////
+	
+void	DetectorManager::BuildSimplePhysicalEvent()					
+	{
+		map<string,VDetector*>::iterator it	;
+		
+		for( it = m_Detector.begin() ; it != m_Detector.end() ; ++it) 
+			{
+				it->second->BuildSimplePhysicalEvent()	;	
+			}    
+	}
+/////////////////////////////////////////////////////////////////////////////////////////////////
+	
+void	DetectorManager::InitializeRootInput()						
+	{
+		map<string,VDetector*>::iterator it	;
+		
+		for( it = m_Detector.begin() ; it != m_Detector.end() ; ++it) 
+			{
+				it->second->InitializeRootInput()	;	
+			}    
+	}
+/////////////////////////////////////////////////////////////////////////////////////////////////	
+
+void	DetectorManager::InitializeRootOutput()						
+	{
+		map<string,VDetector*>::iterator it	;
+		
+		for( it = m_Detector.begin() ; it != m_Detector.end() ; ++it) 
+			{
+				it->second->InitializeRootOutput()	;	
+			}   
+		
+	}
+/////////////////////////////////////////////////////////////////////////////////////////////////
+	
+void	DetectorManager::AddDetector(string DetectorName , VDetector* newDetector)	
+	{
+		m_Detector["MUST2"] = newDetector ;
+	}
+/////////////////////////////////////////////////////////////////////////////////////////////////
+void	DetectorManager::ClearEventPhysics()	
+	{
+		map<string,VDetector*>::iterator it	;
+		
+		for( it = m_Detector.begin() ; it != m_Detector.end() ; ++it) 
+			{
+				it->second->ClearEventPhysics()	;	
+			}   
+		
+	}
+/////////////////////////////////////////////////////////////////////////////////////////////////
+void	DetectorManager::ClearEventData()
+	{
+		map<string,VDetector*>::iterator it	;
+		
+		for( it = m_Detector.begin() ; it != m_Detector.end() ; ++it) 
+			{
+				it->second->ClearEventData()	;	
+			}   
+		
+	}
+			

NPAnalysis/10He_Riken/src/GNUmakefile

+###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+= -L$(NPLIB)/lib -lVDetector -lIORoot -lReaction -lEnergyLoss \
+			-lMust2Data -lMust2Physics \
+			-lAnnularS1Data -lGaspardData \
+			-lInitialConditions -lInteractionCoordinates
+LDFLAGS+= -L$(CLHEP_LIB_DIR) -l$(CLHEP_LIB) 
+
+SRC= $(wildcard *.cc)
+OBJ=$(SRC:.cc=.o)
+
+#all:$(EXEC)
+#	@$(CPP) -o $@ -c $< $(CXXFLAGS)
+
+Analysis:$(OBJ)
+	@$(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)

NPSimulation/src/ThinSi.cc

@@ -44,15 +44,15 @@ namespace THINSI
    const G4double ResoEnergy  = 0.042 ;// = 100keV of Resolution   //   Unit is MeV/2.35
 
    // Geometry
-   const G4double DetectorSize   = 70*mm           ;
+   const G4double DetectorSize   	= 70*mm           ;
    const G4double SiliconThickness  = 20*micrometer      ;
-   const G4double SiliconSize    = 50*mm           ;
+   const G4double SiliconSize    	= 50*mm           ;
    const G4double AluThickness      = 0.4*micrometer  ;
-   const G4int  NumberOfStrip    = 32           ;
+   const G4int  NumberOfStrip    	= 32           ;
 
-   const G4double AluStripFront_PosZ = -0.5*SiliconThickness - 0.5*AluThickness ;
-   const G4double Si_PosZ        = 0                                 ;
-   const G4double AluStripBack_PosZ  = 0.5*SiliconThickness + 0.5*AluThickness  ;
+   const G4double AluStripFront_PosZ 	= -0.5*SiliconThickness - 0.5*AluThickness ;
+   const G4double Si_PosZ        		= 0                                 ;
+   const G4double AluStripBack_PosZ  	= 0.5*SiliconThickness + 0.5*AluThickness  ;
 
    const G4double Si_PosX_Shift = 4*mm ;
    const G4double Si_PosY_Shift = 2*mm ;