From 5b7aa9d1908b04725715ff57a83a7cdd32104fdc Mon Sep 17 00:00:00 2001
From: matta <matta@npt>
Date: Thu, 18 Feb 2010 06:06:11 +0000
Subject: [PATCH] *Change algorithm in TMust2Physics 	-Now test for
 compatibility in SiLi and CsI time 	-Allow time to be not coded, usefull
 for real data

---
 Inputs/EventGenerator/10He.reaction           |  2 +-
 Inputs/EventGenerator/60Fe.reaction           |  8 +-
 .../10He_Riken/include/ObjectManager.hh       | 60 ++++++------
 NPAnalysis/10He_Riken/macro/CrossSection.c    | 15 +--
 NPAnalysis/10He_Riken/src/Analysis.cc         | 98 +++++++++++++------
 NPAnalysis/e530/RunToTreat.txt                |  8 +-
 NPAnalysis/e530/pipo.txt                      |  2 +-
 NPAnalysis/e530/src/Analysis.cc               | 16 ++-
 NPLib/MUST2/TMust2Data.h                      |  2 +-
 NPLib/MUST2/TMust2Physics.cxx                 | 40 +++++---
 10 files changed, 156 insertions(+), 95 deletions(-)

diff --git a/Inputs/EventGenerator/10He.reaction b/Inputs/EventGenerator/10He.reaction
index 8567e32e3..b8a14bae8 100644
--- a/Inputs/EventGenerator/10He.reaction
+++ b/Inputs/EventGenerator/10He.reaction
@@ -7,7 +7,7 @@ TransfertToResonance
 	Target= 2H
 	Light= 3He
 	Heavy= 10He
-	ExcitationEnergy= 0.0
+	ExcitationEnergy= 5.0
 	BeamEnergy= 550
 	BeamEnergySpread= 0
 	SigmaThetaX= 0.6921330164
diff --git a/Inputs/EventGenerator/60Fe.reaction b/Inputs/EventGenerator/60Fe.reaction
index 467bd51d9..ea0874438 100644
--- a/Inputs/EventGenerator/60Fe.reaction
+++ b/Inputs/EventGenerator/60Fe.reaction
@@ -10,10 +10,10 @@ Transfert
 	ExcitationEnergy= 2.0
 	BeamEnergy= 1200
 	BeamEnergySpread= 0
-	SigmaX= 1
-	SigmaY= 1
-	SigmaThetaX= 2 
-	SigmaPhiY= 2
+	SigmaX= 3
+	SigmaY= 3
+	SigmaThetaX= 0.5 
+	SigmaPhiY= 0.5
 	CrossSectionPath= ni69_g7_01.n
 	ShootLight= 1
 	ShootHeavy= 0
diff --git a/NPAnalysis/10He_Riken/include/ObjectManager.hh b/NPAnalysis/10He_Riken/include/ObjectManager.hh
index c587cb506..47e09a3cb 100644
--- a/NPAnalysis/10He_Riken/include/ObjectManager.hh
+++ b/NPAnalysis/10He_Riken/include/ObjectManager.hh
@@ -75,7 +75,7 @@ namespace GRAPH
 	{
 		//	Declare your Spectra here:
 	
-			TH1F *myHist1D = new TH1F("Hist1D","Histogramm 1D ; x ; count", 1000 , -5 , 5 )					;
+			TH1F* myHist1D = new TH1F("Hist1D","Histogramm 1D ; x ; count", 100 , -40 , 40 )					;
 	
 			TH2F *myHist2D = new TH2F("Hist2D","Histogramm 2D ; x ; y ", 128 , 1 , 128 , 128 , 1 , 128 )	;
 
@@ -111,56 +111,56 @@ namespace ENERGYLOSS
 	//	3He Energy Loss
 			EnergyLoss He3TargetWind = EnergyLoss 	(	"He3_Mylar.G4table" 		,
 																								"G4Table",
-																									10000	 				);
+																									1000,
+																									3	 				);
 		
 			EnergyLoss He3TargetGaz = EnergyLoss 		(	"He3_D2.G4table" 	,
 																								"G4Table",
-																									10000	 				);
+																									1000,
+																									3	 				);
 			
 			EnergyLoss He3StripAl   = EnergyLoss 	(	"He3_Aluminium.G4table" 			,
 																							"G4Table",
-																							10000						);
-														
-			EnergyLoss He3StripSi   = EnergyLoss 	(	"3He_Si.txt" 			,
-																							"LISE",
-																							10000						,
-																							1					,
+																							1000,
 																							3						);
-														
 
-	
-//		//	3He Energy Loss
 //			EnergyLoss He3TargetWind = EnergyLoss 	(	"3He_Mylar.txt" 		,
-//														10000	 				,
-//														1					,
-//														3						);
+//																								"LISE",
+//																							10000						,
+//																							1					,
+//																							3	);					
 //		
-//			EnergyLoss He3TargetGaz = EnergyLoss 	(	"3He_D2gaz_1b_26K.txt" 	,
-//														10000	 				,
-//														1						,
-//														3						);
+//			EnergyLoss He3TargetGaz = EnergyLoss 		(	"3He_D2gaz_1b_26K.txt" 	,
+//																								"LISE",
+//																							10000						,
+//																							1					,
+//																							3					);	
 //			
-//			EnergyLoss He3StripAl   = EnergyLoss 	(	"3He_Al.txt" 			,
-//														10000						,
-//														1						,
-//														3						);
-//														
-//			EnergyLoss He3StripSi   = EnergyLoss 	(	"3He_Si.txt" 			,
-//														10000						,
-//														1					,
-//														3						);
+//			EnergyLoss He3StripAl   = EnergyLoss 	(		"3He_Al.txt"			,
+//																							"LISE",
+//																							10000						,
+//																							1					,
+//																							3			);			
+
+
+														
+			EnergyLoss He3StripSi   = EnergyLoss 	(	"3He_Si.txt" 			,
+																							"LISE",
+																							100						,
+																							3					,
+																							1						);
 														
 														
 		//	proton Energy Loss
 			EnergyLoss protonTargetWind = EnergyLoss 	(	"proton_Mylar.txt"	 		,
 															"LISE",
-															1000		 				,
+															100		 				,
 															1						,
 															1							);
 		
 			EnergyLoss protonTargetGaz = EnergyLoss 	(	"proton_D2gaz_1b_26K.txt" 	,
 																"LISE",
-															1000		 				,
+															100		 				,
 															1						,
 															1							);
 			
diff --git a/NPAnalysis/10He_Riken/macro/CrossSection.c b/NPAnalysis/10He_Riken/macro/CrossSection.c
index 5640efb44..56e071c82 100644
--- a/NPAnalysis/10He_Riken/macro/CrossSection.c
+++ b/NPAnalysis/10He_Riken/macro/CrossSection.c
@@ -11,16 +11,8 @@
 	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();
-	
+TFile *file0 = TFile::Open("./Result/thetaCM.root");
+
 	cCM = new TCanvas("cCm" , "Cross Section (CM)" , 100 , 100 , 900, 900) ;
 	hThetaCM->Draw();
 	
@@ -79,8 +71,7 @@ TFile *file0 = TFile::Open("./Result/myResult.root");
 	
 	
 	//Normalisation:
-	//Int_t Maximum_Bin 		= hCrossSection->GetMaximumBin()			;
-	Int_t Maximum_Bin 		= 3 										;
+	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)	;
diff --git a/NPAnalysis/10He_Riken/src/Analysis.cc b/NPAnalysis/10He_Riken/src/Analysis.cc
index 8243ec186..893c375b0 100644
--- a/NPAnalysis/10He_Riken/src/Analysis.cc
+++ b/NPAnalysis/10He_Riken/src/Analysis.cc
@@ -17,7 +17,6 @@ int main(int argc,char** argv)
 	string calibrationfileName 	= argv[2]	;
 	string runToReadfileName 		= argv[3]	;
 	
-	/////////////////////////////////////////////////////////////////////////////////////////////////////
 	//	First of All instantiate RootInput and Output
 	//	Detector will be attached later
 	RootInput:: getInstance(runToReadfileName)	;
@@ -30,40 +29,40 @@ int main(int argc,char** argv)
 	NPL::Reaction*  Li10Reaction = new Reaction								;
 	Li10Reaction	->	ReadConfigurationFile("9Li-dp-10Li.reaction")		;
 
-	//	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)		;
-	
-	//	Ask the detector manager to load the parameter added by the detector in the calibrationfileName
-	myCalibration->LoadParameterFromFile() ;
-	/////////////////////////////////////////////////////////////////////////////////////////////////////
-	
-
 
+	//	Instantiate the Calibration Manger using a file
+	CalibrationManager* myCalibration = CalibrationManager::getInstance(calibrationfileName) ;
+	
 	//	Attach more branch to the output
 	
 	double ELab[2], ExcitationEnergy[2]	;
 	double ThetaLab[2]	, ThetaCM[2]		;
 	double X[2] , Y[2] 									;
+	double EventWeight=0;
 
 	// Exitation Energy
 	RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy",ExcitationEnergy,"ExcitationEnergy[2]/D") ;
 
 	//E lab et Theta lab
-	RootOutput::getInstance()->GetTree()->Branch("ThetaCM",ThetaCM,"ThetaCM[2]/D") 			;
-	RootOutput::getInstance()->GetTree()->Branch("ELab",ELab,"ELab[2]/D") 							;
-	RootOutput::getInstance()->GetTree()->Branch("ThetaLab",ThetaLab,"ThetaLab[2]/D") 	;
-	RootOutput::getInstance()->GetTree()->Branch("X",X,"X/D[2]") 												;
-	RootOutput::getInstance()->GetTree()->Branch("Y",Y,"Y/D[2]")									 			;
+	RootOutput::getInstance()->GetTree()->Branch("ThetaCM",ThetaCM,"ThetaCM[2]/D") 	;
+	RootOutput::getInstance()->GetTree()->Branch("ELab",ELab,"ELab[2]/D") 					;
+	RootOutput::getInstance()->GetTree()->Branch("ThetaLab",ThetaLab,"ThetaLab[2]/D") 					;
+	RootOutput::getInstance()->GetTree()->Branch("X",X,"X/D[2]") 										;
+	RootOutput::getInstance()->GetTree()->Branch("Y",Y,"Y/D[2]")									 	;
 	
+	// For phasespace event
+	RootOutput::getInstance()->GetTree()->Branch("EventWeight",&EventWeight,"EventWeight/D")									 	;
 	
 	//	Get the formed Chained Tree and Treat it
 	TChain* Chain = RootInput:: getInstance() -> GetChain()	;
 
 	// Open the ThinSi Branch
+	Chain->SetBranchStatus("EventWeight",true);
+	Chain->SetBranchAddress("EventWeight"	,&EventWeight		);
+	
 	Chain->SetBranchStatus("InitialConditions",true)	;
 	Chain->SetBranchStatus("fIC_*",true)	; 
  		
@@ -79,8 +78,38 @@ int main(int argc,char** argv)
 	TMust2Physics* M2 		= (TMust2Physics*) 			myDetector -> m_Detector["MUST2"] 	;
 	TPlasticPhysics* Pl 	= (TPlasticPhysics*) 		myDetector -> m_Detector["Plastic"] ;
 	TSSSDPhysics* ThinSi 	= (TSSSDPhysics*) 			myDetector -> m_Detector["SSSD"] 		;
-	cout <<  " ///////// Starting Analysis ///////// "<< endl << endl ;
-	
+
+
+RootOutput::getInstance()->GetList()->Add(myHist1D);
+
+   TCutG *cut3He_MUST2 = new TCutG("Cut3HeMUST2",11);
+   cut3He_MUST2->SetPoint(0,-1.49426,24.2781);
+   cut3He_MUST2->SetPoint(1,15.3161,14.3151);
+   cut3He_MUST2->SetPoint(2,47.069,7.6732);
+   cut3He_MUST2->SetPoint(3,110.575,4.35222);
+   cut3He_MUST2->SetPoint(4,308.563,2.25477);
+   cut3He_MUST2->SetPoint(5,310.431,1.20604);
+   cut3He_MUST2->SetPoint(6,232.917,1.11864);
+   cut3He_MUST2->SetPoint(7,89.0948,2.42955);
+   cut3He_MUST2->SetPoint(8,32.1264,6.71186);
+   cut3He_MUST2->SetPoint(9,1.30747,11.9555);
+   cut3He_MUST2->SetPoint(10,-1.49426,24.2781);
+   
+   TCutG *cut3He_M2_SSSD = new TCutG("Cut3HeM2SSSD",11);
+   cut3He_M2_SSSD->SetPoint(0,7.44252,1.45432);
+   cut3He_M2_SSSD->SetPoint(1,40.6322,0.684984);
+   cut3He_M2_SSSD->SetPoint(2,95.9483,0.406912);
+   cut3He_M2_SSSD->SetPoint(3,218.649,0.216896);
+   cut3He_M2_SSSD->SetPoint(4,362.471,0.101033);
+   cut3He_M2_SSSD->SetPoint(5,355.431,0.0315148);
+   cut3He_M2_SSSD->SetPoint(6,100.977,0.0315148);
+   cut3He_M2_SSSD->SetPoint(7,25.546,0.314221);
+   cut3He_M2_SSSD->SetPoint(8,3.41954,0.995498);
+   cut3He_M2_SSSD->SetPoint(9,3.41954,1.40797);
+   cut3He_M2_SSSD->SetPoint(10,7.44252,1.45432);
+
+
+	cout <<  " ///////// Starting Analysis ///////// "<< endl << endl ;	
 	int i ,N=Chain -> GetEntries();
 	
 	cout << " Number of Event to be treated : " << N << endl ;
@@ -91,8 +120,8 @@ int main(int argc,char** argv)
 			// Clear local branch
 			for(int hh = 0 ; hh <2 ; hh++)
 				{
-					ELab[hh] = -1 	; ExcitationEnergy[hh] = -1 ; ThetaLab[hh] = -1 ;	
-					X[hh]	 = -1000	; Y[hh]					= -1000 ; ThetaCM[hh]  = -1 ;
+					ELab[hh] = -10000 	; ExcitationEnergy[hh] = -10000 ; ThetaLab[hh] = -10000 ;	
+					X[hh]	 = -10000	; Y[hh]					= -10000 ; ThetaCM[hh]  = -10000 ;
 				}
 
 			// Minimum code
@@ -124,12 +153,11 @@ int main(int argc,char** argv)
 			//	YTarget = RandomEngine.Gaus(Init->GetICPositionY(0),1);
 			BeamTheta = Init->GetICIncidentAngleTheta(0)*deg ; BeamPhi = Init->GetICIncidentAnglePhi(0)*deg ; 
 			TVector3 BeamDirection = TVector3(cos(BeamPhi)*sin(BeamTheta) , sin(BeamPhi)*sin(BeamTheta) , cos(BeamTheta)) ;
-			////
+			//// 
 			
 			// Must 2 And ThinSi //
 			for(int hit = 0; hit < M2 -> Si_E.size() ; hit ++)
 				{
-					ELab[hit] = -1 ; ThetaLab[hit] = -1;
 					//	Get Hit Direction
 					TVector3 HitDirection  = M2 -> GetPositionOfInteraction(hit) - TVector3(XTarget,YTarget,0);
 					// Angle between beam and particle
@@ -142,21 +170,24 @@ int main(int argc,char** argv)
 
 					if(M2 -> GetPositionOfInteraction(hit).Z() > 0)
 						{
-							if( M2 -> CsI_E[hit] == 0 && M2 -> Si_E[hit] > 0)
+							if( M2 -> CsI_E[hit] == 0 && M2 -> Si_E[hit] > 0 )
 								{  
 										ELab[hit] = M2 -> Si_E[hit]  ; 
 										
 										ELab[hit]= He3StripAl.EvaluateInitialEnergy(			ELab[hit] 				, // Energy of the detected particle
-																																			2*0.4*micrometer	, // Target Thickness at 0 degree
+																																			0.4*micrometer	, // Target Thickness at 0 degree
 																																			ThetaMM2Surface		);														
 										
 										if(ThinSi -> Energy.size() > 0)
 											{
+												ELab[hit]= He3StripAl.EvaluateInitialEnergy(	ELab[hit] 				, // Energy of the detected particle
+																																			0.4*micrometer		, // Target Thickness at 0 degree
+																																			ThetaMM2Surface		);
 												ELab[hit] += ThinSi-> Energy[hit];
 												ELab[hit]= He3StripAl.EvaluateInitialEnergy(	ELab[hit] 				, // Energy of the detected particle
 																																			0.4*micrometer		, // Target Thickness at 0 degree
 																																			ThetaMM2Surface		);
-											}
+											} 
 
 										ELab[hit]= He3TargetWind.EvaluateInitialEnergy( 	ELab[hit] 				, // Energy of the detected particle
 																																			15*micrometer			, // Target Thickness at 0 degree
@@ -168,7 +199,11 @@ int main(int argc,char** argv)
 																		 				
 									ThetaCM[hit] = He10Reaction -> EnergyLabToThetaCM( ELab[hit] ) /deg 	;
 									ExcitationEnergy[hit] = He10Reaction -> ReconstructRelativistic( ELab[hit] , ThetaLab[hit] ) 		;	
-//								  	ExcitationEnergy[hit] = He10Reaction -> ReconstructRelativistic( ELab[hit] , TVRAI ) 		;
+									
+									if(ThinSi -> Energy.size() > 0 )
+								  	if(cut3He_M2_SSSD->IsInside(ThinSi -> Energy[hit], M2 -> Si_E[hit]) )
+											myHist1D->Fill(ExcitationEnergy[hit],EventWeight);
+									
 									X[hit] = HitDirection . X();
 									Y[hit] = HitDirection . Y();	
 									ThetaLab[hit] = ThetaLab[hit] / deg ;
@@ -194,6 +229,9 @@ int main(int argc,char** argv)
 								
 									if(ThinSi -> Energy.size() > 0)
 										{
+											ELab[hit]= He3StripAl.EvaluateInitialEnergy(	ELab[hit] 					, // Energy of the detected particle
+																																		0.4*micrometer			, // Target Thickness at 0 degree
+																																		ThetaMM2Surface			);
 											ELab[hit] += ThinSi-> Energy[hit];
 											ELab[hit]= He3StripAl.EvaluateInitialEnergy(	ELab[hit] 					, // Energy of the detected particle
 																																		0.4*micrometer			, // Target Thickness at 0 degree
@@ -207,10 +245,13 @@ int main(int argc,char** argv)
 									ELab[hit]= He3TargetGaz.EvaluateInitialEnergy(		ELab[hit] 					, // Energy of the detected particle
 																																		1.5*mm							, // Target Thickness at 0 degree
 																																		ThetaN							);
-																				
+																																		
 									ThetaCM[hit]= He10Reaction -> EnergyLabToThetaCM( ELab[hit] ) /deg 	;	
 									ExcitationEnergy[hit] = He10Reaction -> ReconstructRelativistic( ELab[hit], ThetaLab[hit]) ;	
-//									ExcitationEnergy[hit] = He10Reaction -> ReconstructRelativistic( ELab[hit] , TVRAI ) 		;
+										
+									if( cut3He_MUST2->IsInside(M2 -> Si_E[hit], M2 -> CsI_E[hit]) )
+										myHist1D->Fill(ExcitationEnergy[hit],EventWeight);
+										
 									X[hit] = HitDirection . X();
 									Y[hit] = HitDirection . Y();	
 									ThetaLab[hit] = ThetaLab[hit] / deg ;
@@ -223,14 +264,13 @@ int main(int argc,char** argv)
 					/*else if(M2 -> GetPositionOfInteraction(hit).Z()<0)
 						{}	*/
 
-
 				}			
-			
 			RootOutput::getInstance()->GetTree()->Fill()	;
 		}
 
 	cout << " A total of " << i << " event has been annalysed " << endl ;
 	cout << endl << "/////////////////////////////////"<< endl<< endl ;
+	myHist1D->Write();
 	RootOutput::getInstance()->Destroy();
 	return 0	;
 }
diff --git a/NPAnalysis/e530/RunToTreat.txt b/NPAnalysis/e530/RunToTreat.txt
index 1a9cc7e28..c56d1ce44 100644
--- a/NPAnalysis/e530/RunToTreat.txt
+++ b/NPAnalysis/e530/RunToTreat.txt
@@ -1,5 +1,7 @@
 TTreeName 
-	AutoTree
+	SimulatedTree
 RootFileName 
-	../../../../root/run_0003.root
-
+	%../../Outputs/Simulation/3He_source.root
+	%../../Outputs/Simulation/alpha_source.root
+	%../../Outputs/Simulation/REALIST_DATA.root
+	../../Outputs/Simulation/mySimul.root
diff --git a/NPAnalysis/e530/pipo.txt b/NPAnalysis/e530/pipo.txt
index 853f2240a..a6ffecd1b 100644
--- a/NPAnalysis/e530/pipo.txt
+++ b/NPAnalysis/e530/pipo.txt
@@ -1,2 +1,2 @@
 CalibrationFilePath
-	/home/Adrien/Desktop/NPTool/NPTool.dev.CalibMust2/NPAnalysis/e530/calibPipo.txt
+	/home/Adrien/Desktop/NPTool/NPTool.dev.CalibMust2/NPAnalysis/e530/calibrrrPipo.txt
diff --git a/NPAnalysis/e530/src/Analysis.cc b/NPAnalysis/e530/src/Analysis.cc
index c0b02c5d5..2561d14e4 100644
--- a/NPAnalysis/e530/src/Analysis.cc
+++ b/NPAnalysis/e530/src/Analysis.cc
@@ -41,7 +41,10 @@ int main(int argc,char** argv)
 	/////////////////////////////////////////////////////////////////////////////////////////////////////
 
   //	Attach more branch to the output
-	
+  double ELab=0,ThetaLab=0,ExcitationEnergy=0;
+	RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D") 							;
+	RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D") 	;
+	RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy",&ExcitationEnergy,"ExcitationEnergy/D") ;
   //	Get the formed Chained Tree and Treat it
   TChain* Chain = RootInput:: getInstance() -> GetChain()	;
 
@@ -80,6 +83,17 @@ int main(int argc,char** argv)
       ////
 		
 			
+			// Must 2
+			for(int hit = 0; hit < M2 -> Si_E.size() ; hit ++)
+				{
+					ELab = -1 ; ThetaLab = -1;
+					//	Get Hit Direction
+					TVector3 HitDirection  = M2 -> GetPositionOfInteraction(hit) - TVector3(0,0,-40);
+					// Angle between beam and particle
+					ThetaLab  = ThetaCalculation ( HitDirection , TVector3(0,0,1)   ) ;	
+					ELab = M2 -> Si_E[hit] + M2 -> SiLi_E[hit]	;
+			  }
+			
       RootOutput::getInstance()->GetTree()->Fill()	;
     }
 
diff --git a/NPLib/MUST2/TMust2Data.h b/NPLib/MUST2/TMust2Data.h
index 654e999cd..511cfc0ca 100644
--- a/NPLib/MUST2/TMust2Data.h
+++ b/NPLib/MUST2/TMust2Data.h
@@ -151,7 +151,7 @@ class TMust2Data : public TObject {
 
 		// CsI 
 		//(E)
-		UShort_t	GetMMCsIEMult()					{return fMM_CsIE_DetectorNbr.size();}
+		UShort_t	GetMMCsIEMult()									{return fMM_CsIE_DetectorNbr.size();}
 		UShort_t	GetMMCsIEDetectorNbr(Int_t i)   {return fMM_CsIE_DetectorNbr.at(i);}
 		UShort_t	GetMMCsIECristalNbr(Int_t i)    {return fMM_CsIE_CristalNbr.at(i);}
 		Double_t	GetMMCsIEEnergy(Int_t i)        {return fMM_CsIE_Energy.at(i);}
diff --git a/NPLib/MUST2/TMust2Physics.cxx b/NPLib/MUST2/TMust2Physics.cxx
index 07243af1b..1cf3262a3 100644
--- a/NPLib/MUST2/TMust2Physics.cxx
+++ b/NPLib/MUST2/TMust2Physics.cxx
@@ -16,8 +16,6 @@
  *                                                                           *
  *---------------------------------------------------------------------------*
  * Comment:                                                                  *
- *  Only multiplicity one and multiplicity 2 are down.                       *
- *  Improvment needed                                                        *
  *                                                                           *
  *****************************************************************************/
 #include "TMust2Physics.h"
@@ -95,14 +93,24 @@ void TMust2Physics::BuildPhysicalEvent()
 							{
 								if(EventData->GetMMSiLiEDetectorNbr(j)==N)
 									{
-										//	if SiLi energy is above threshold check the compatibility
+										// SiLi energy is above threshold check the compatibility
 										if( fSiLi_E(EventData , j)>SiLi_E_Threshold )
 											{
+												// pad vs strip number match
 												if( Match_Si_SiLi( X, Y , EventData->GetMMSiLiEPadNbr(j) ) )
 												{
 													SiLi_N.push_back(EventData->GetMMSiLiEPadNbr(j))	;
 													SiLi_E.push_back(fSiLi_E(EventData , j))					;
-													SiLi_T.push_back(fSiLi_T(EventData , j))		;
+													
+													// Look for associate energy
+													// Note: in case of use of SiLi "Orsay" time is not coded.
+													for(int k =0 ; k  < EventData->GetMMSiLiTMult() ; k ++)
+														{
+															// Same Pad, same Detector
+															if( EventData->GetMMSiLiEPadNbr(j)==EventData->GetMMSiLiEPadNbr(k) && EventData->GetMMSiLiEDetectorNbr(j)==EventData->GetMMSiLiTDetectorNbr(k) )
+																{SiLi_T.push_back(fSiLi_T(EventData , k))		; break ;}
+														}
+													
 													check_SILI = true ;
 													
 												}
@@ -114,14 +122,20 @@ void TMust2Physics::BuildPhysicalEvent()
 							{
 								if(EventData->GetMMCsIEDetectorNbr(j)==N)
 									{
-										//	ifCsI energy is above threshold check the compatibility
+										//	CsI energy is above threshold check the compatibility
 										if( fCsI_T(EventData , j)>CsI_E_Threshold )
 											{
-												if( Match_Si_CsI( X, Y , EventData->GetMMCsIECristalNbr(j) ) )
+												if(Match_Si_CsI( X, Y , EventData->GetMMCsIECristalNbr(j) ) )
 													{
 														CsI_N.push_back(EventData->GetMMCsIECristalNbr(j))	;
-														CsI_E.push_back(fCsI_E(EventData , j))			;
-														CsI_T.push_back(fCsI_T(EventData , j))				;
+														CsI_E.push_back(fCsI_E(EventData , j))							;
+														
+														for(int k =0 ; k  < EventData->GetMMCsITMult() ; k ++)
+															{
+																if( EventData->GetMMCsIECristalNbr(j)==EventData->GetMMCsITCristalNbr(k) && EventData->GetMMCsIEDetectorNbr(j)==EventData->GetMMCsITDetectorNbr(k) )
+																	{CsI_T.push_back(fCsI_T(EventData , k))	; break ;}
+															}
+														
 														check_CSI = true ;
 													}
 											}
@@ -132,15 +146,15 @@ void TMust2Physics::BuildPhysicalEvent()
 						if(!check_SILI)
 							{
 								SiLi_N.push_back(0)	;
-								SiLi_E.push_back(0)	;
-								SiLi_T.push_back(0)	;
+								SiLi_E.push_back(-10000)	;
+								SiLi_T.push_back(-10000)	;
 							}
 
 						if(!check_CSI) 
 							{
 								CsI_N.push_back(0)	;
-								CsI_E.push_back(0)	;
-								CsI_T.push_back(0)	;
+								CsI_E.push_back(-10000)	;
+								CsI_T.push_back(-10000)	;
 							}
 					
 					}
@@ -196,7 +210,7 @@ vector < TVector2 > TMust2Physics :: Match_X_Y()
 										//	if same detector check energy
 										if ( EventData->GetMMStripXEDetectorNbr(i) == EventData->GetMMStripYEDetectorNbr(j) )
 											{
-													//	Look if energy match
+													//	Look if energy match (within 10%)
 													if( ( fSi_X_E(EventData , i) - fSi_Y_E(EventData , j) ) / fSi_X_E(EventData , i) < 0.1	)
 														ArrayOfGoodCouple . push_back ( TVector2(i,j) ) ;	
 											}
-- 
GitLab