* Improving EventGeneratorTransfertToResonance

	- Now the weight of each event is stored instead of using a rejection methods when making the phase space decay

* Adding a pure phase space Event Generator
	- Using an input file, one can simulate any many body phase space (up to 18 body)
	- Same weight methods is used as in TransfertToResonance

Inputs/DetectorConfiguration/Riken_65mm.detector

@@ -22,7 +22,7 @@ CryoTarget
 	TEMPERATURE= 26
 	PRESSURE= 1
 	MATERIAL= D2	
-	WINDOWSTHICKNESS= 15
+	WINDOWSTHICKNESS= 10
 	WINDOWSMATERIAL= Mylar
 	X= 0
 	Y= 0

Inputs/EventGenerator/10He.reaction

@@ -7,19 +7,19 @@ TransfertToResonance
 	Target= 2H
 	Light= 3He
 	Heavy= 10He
-	ExcitationEnergy= 0.0
+	ExcitationEnergy= 1.0
 	BeamEnergy= 550
 	BeamEnergySpread= 0
 	SigmaThetaX= 0.6921330164
 	SigmaPhiY= 0.963142053
 	SigmaX= 6.232
 	SigmaY= 9.069
-	ResonanceWidth= 0
+	ResonanceWidth= 1
 	ResonanceDecayZ= 2
 	ResonanceDecayA= 8
 	CrossSectionPath= 11Li(d,3He)10He.txt	
 	ShootLight= 1
-	ShootHeavy= 1
+	ShootHeavy= 0
 	ShootDecayProduct= 0
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 

Inputs/EventGenerator/34Si.reaction

+Transfert
+       Beam= 19C
+       Target= 208Pb
+       Light= 208Pb
+       Heavy= 19C
+       ExcitationEnergy= 5.9
+       BeamEnergy= 722
+       BeamEnergySpread= 0
+       SigmaThetaX= 0.6921330164
+       SigmaPhiY= 0.963142053
+       SigmaX= 6.232
+       SigmaY= 9.069
+       CrossSectionPath= flat.txt
+       ShootLight= 0
+       ShootHeavy= 1
+
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%% Reaction file for 11Li(d,3He)10He reaction %%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%Beam energy given in MeV ; Excitation in MeV 748  612
-Transfert
-	Beam= 34Si
-	Target= 2H
-	Light= 1H
-	Heavy= 35Si
-	ExcitationEnergy= 2.0
-	BeamEnergy= 748
-	BeamEnergySpread= 0
-	SigmaX= 0.5
-	SigmaY= 0.5
-	SigmaThetaX= 0
-	SigmaPhiY= 0
-	BeamSpreadX= 0
-	CrossSectionPath= ni69_g7_01.n
-	ShootLight= 1
-	ShootHeavy= 0
+%Transfert
+%	Beam= 34Si
+%	Target= 2H
+%	Light= 1H
+%	Heavy= 35Si
+%	ExcitationEnergy= 2.0
+%	BeamEnergy= 748
+%	BeamEnergySpread= 0
+%	SigmaX= 0.5
+%	SigmaY= 0.5
+%	SigmaThetaX= 0
+%	SigmaPhiY= 0
+%	CrossSectionPath= flat.txt
+%	ShootLight= 1
+%	ShootHeavy= 0
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Inputs/EventGenerator/3He.source

@@ -5,7 +5,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Isotropic
 	EnergyLow=  0	
-	EnergyHigh= 21.6
+	EnergyHigh= 300
 	HalfOpenAngleMin= 0
 	HalfOpenAngleMax= 90
 	x0= 0	

Inputs/EventGenerator/alpha.source

@@ -5,7 +5,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Isotropic
 	EnergyLow= 0
-	EnergyHigh= 21.6
+	EnergyHigh= 300
 	HalfOpenAngleMin= 0
 	HalfOpenAngleMax= 90
 	x0= 0	

NPAnalysis/10He_Riken/include/ObjectManager.hh

@@ -120,6 +120,26 @@ namespace ENERGYLOSS
 			EnergyLoss He3StripAl   = EnergyLoss 	(	"He3_Aluminium.G4table" 			,
 																							"G4Table",
 																							10000						);
+
+//			EnergyLoss He3TargetWind = EnergyLoss 	(	"3He_Mylar.txt" 		,
+//																								"LISE",
+//																							10000						,
+//																							1					,
+//																							3	);					
+//		
+//			EnergyLoss He3TargetGaz = EnergyLoss 		(	"3He_D2gaz_1b_26K.txt" 	,
+//																								"LISE",
+//																							10000						,
+//																							1					,
+//																							3					);	
+//			
+//			EnergyLoss He3StripAl   = EnergyLoss 	(		"3He_Al.txt"			,
+//																							"LISE",
+//																							10000						,
+//																							1					,
+//																							3			);			
+
+
 														
 			EnergyLoss He3StripSi   = EnergyLoss 	(	"3He_Si.txt" 			,
 																							"LISE",
@@ -127,29 +147,6 @@ namespace ENERGYLOSS
 																							1					,
 																							3						);
 														
-
-	
-//		//	3He Energy Loss
-//			EnergyLoss He3TargetWind = EnergyLoss 	(	"3He_Mylar.txt" 		,
-//														10000	 				,
-//														1					,
-//														3						);
-//		
-//			EnergyLoss He3TargetGaz = EnergyLoss 	(	"3He_D2gaz_1b_26K.txt" 	,
-//														10000	 				,
-//														1						,
-//														3						);
-//			
-//			EnergyLoss He3StripAl   = EnergyLoss 	(	"3He_Al.txt" 			,
-//														10000						,
-//														1						,
-//														3						);
-//														
-//			EnergyLoss He3StripSi   = EnergyLoss 	(	"3He_Si.txt" 			,
-//														10000						,
-//														1					,
-//														3						);
-														
 														
 		//	proton Energy Loss
 			EnergyLoss protonTargetWind = EnergyLoss 	(	"proton_Mylar.txt"	 		,

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)	;

NPAnalysis/10He_Riken/src/Analysis.cc

@@ -117,7 +117,7 @@ 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 ++)
@@ -152,7 +152,7 @@ int main(int argc,char** argv)
 											}
 
 										ELab[hit]= He3TargetWind.EvaluateInitialEnergy( 	ELab[hit] 				, // Energy of the detected particle
-																																			15*micrometer			, // Target Thickness at 0 degree
+																																			10*micrometer			, // Target Thickness at 0 degree
 																																			ThetaN						);
 																			
 										ELab[hit]= He3TargetGaz.EvaluateInitialEnergy(		ELab[hit] 				, // Energy of the detected particle
@@ -161,7 +161,6 @@ 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 ) 		;
 									X[hit] = HitDirection . X();
 									Y[hit] = HitDirection . Y();	
 									ThetaLab[hit] = ThetaLab[hit] / deg ;
@@ -194,7 +193,7 @@ int main(int argc,char** argv)
 										}
 											
 									ELab[hit]= He3TargetWind.EvaluateInitialEnergy( 	ELab[hit] 					, // Energy of the detected particle
-																																		15*micrometer				, // Target Thickness at 0 degree
+																																		10*micrometer				, // Target Thickness at 0 degree
 																																		ThetaN							);
 									
 									ELab[hit]= He3TargetGaz.EvaluateInitialEnergy(		ELab[hit] 					, // Energy of the detected particle
@@ -203,7 +202,6 @@ 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 ) 		;
 									X[hit] = HitDirection . X();
 									Y[hit] = HitDirection . Y();	
 									ThetaLab[hit] = ThetaLab[hit] / deg ;
@@ -216,7 +214,6 @@ int main(int argc,char** argv)
 					/*else if(M2 -> GetPositionOfInteraction(hit).Z()<0)
 						{}	*/
 
-
 				}			
 			
 			RootOutput::getInstance()->GetTree()->Fill()	;

NPSimulation/include/EventGeneratorTransfertToResonance.hh

@@ -99,13 +99,17 @@ class EventGeneratorTransfertToResonance : public VEventGenerator
 		Reaction*	m_Reaction				;
 
 	private: // Resonance decay channel
-		double			m_ResonanceWidth	;
-		double			m_ResonanceMean		;
-	   	int            	m_ResonanceDecayZ   ;
-	   	int            	m_ResonanceDecayA   ;
+		double			m_ResonanceWidth		;
+		double			m_ResonanceMean			;
+	   	int       m_ResonanceDecayZ   ;
+	   	int       m_ResonanceDecayA   ;
+	   	 	
+	  //	When the Phase Space Generator is called, the weight of the current configuration is return and stored in this variable
+	  //	Spectrum then need to be weighted by this paramater to be realistic
+	  //	NB: This procedure avoid long calculation time of the rejection methods previously instantiate and therefore allow simulation of manybody phase space decay 	 	
+	  double 			m_EventWeight				;
 
 	   //Other
-	   void     Print() const        ;
 	   void     InitializeRootOutput()  ;
 	   void     ResonanceDecay(	G4double EnergyHeavy    ,
 					         	G4double ThetaHeavy     ,

NPSimulation/src/EventGeneratorTransfertToResonance.cc

@@ -46,32 +46,34 @@
 //Root Headers
 #include "TGenPhaseSpace.h"
 
+//CLHEP
+#include "CLHEP/Random/RandBreitWigner.h"
+
 using namespace std;
 using namespace CLHEP;
 
-
-
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 EventGeneratorTransfertToResonance::EventGeneratorTransfertToResonance()
 {
    //------------- Default Constructor -------------
 	m_InitConditions	= new TInitialConditions()	;
-	m_Reaction = new Reaction() ;
-	m_Target = new Target();
-	m_SigmaX       		=  0 ;
-	m_SigmaY       		=  0 ;
-	m_SigmaThetaX     =  0 ;
-	m_SigmaPhiY 		=  0 ;
-	m_ResonanceDecayZ 		=  0 ;
-	m_ResonanceDecayA 		=  0 ;
+	m_Reaction 				= new Reaction() 						;
+	m_Target 					= new Target()							;
+	m_SigmaX       		=  0 												;
+	m_SigmaY       		=  0 												;
+	m_SigmaThetaX     =  0 												;
+	m_SigmaPhiY 			=  0 												;
+	m_ResonanceDecayZ = 0 												;
+	m_ResonanceDecayA = 0 												;
+	m_EventWeight     =	0 												;
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 EventGeneratorTransfertToResonance::~EventGeneratorTransfertToResonance()
 {
   //------------- Default Destructor ------------
-	delete m_InitConditions;
-	delete m_Reaction ;
+	delete m_InitConditions	;
+	delete m_Reaction 			;
 }
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 void	EventGeneratorTransfertToResonance::SetTarget(Target* Target) 
@@ -119,6 +121,8 @@ EventGeneratorTransfertToResonance::EventGeneratorTransfertToResonance(	  string
 		           	ShootHeavy        	,
 		           	ShootDecayProduct   ,
 		            Path				);        
+		            
+    m_EventWeight     =	0 												;
 
 }
 
@@ -129,13 +133,7 @@ void EventGeneratorTransfertToResonance::InitializeRootOutput()
    RootOutput *pAnalysis = RootOutput::getInstance();
    TTree *pTree = pAnalysis->GetTree();
    pTree->Branch("InitialConditions", "TInitialConditions", &m_InitConditions);
-}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorTransfertToResonance::Print() const
-{
-
+   pTree->Branch("EventWeight",&m_EventWeight,"EventWeigt/D");
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -188,7 +186,7 @@ void EventGeneratorTransfertToResonance::ReadConfiguration(string Path)
 
       
 
-      if (LineBuffer.compare(0, 9, "Transfert") == 0) { ReadingStatus = true ;}
+      if (LineBuffer.compare(0, 20, "TransfertToResonance") == 0) { ReadingStatus = true ;}
 
 
 while(ReadingStatus){
@@ -345,24 +343,24 @@ while(ReadingStatus){
 
 	}
    
-   SetEverything(	Beam            	,
-			         Target            	,
-			         Light          	,
-			         Heavy          	,
-			         BeamEnergy       	,
-			         ExcitationEnergy  	,
-			         BeamEnergySpread  	,
-			         SigmaX         	,
-			         SigmaY         	,
-			         SigmaThetaX ,
-			         SigmaPhiY 	,
-			         ResonanceWidth,
-			         ResonanceDecayZ    ,
-			         ResonanceDecayA    ,
-			         ShootLight        	,
-			         ShootHeavy        	,
-			         ShootDecayProduct  ,
-			         CrossSectionPath	);
+   SetEverything(	Beam            		,
+			         		Target            	,
+			         		Light          			,
+			         		Heavy          			,
+			         		BeamEnergy       		,
+			         		ExcitationEnergy  	,
+			         		BeamEnergySpread  	,
+			         		SigmaX         			,
+			        		SigmaY         			,
+			         		SigmaThetaX				 	,
+			         		SigmaPhiY 					,
+			         		ResonanceWidth			,
+			         		ResonanceDecayZ    	,
+			         		ResonanceDecayA    	,
+			         		ShootLight        	,
+			         		ShootHeavy        	,
+			        		ShootDecayProduct  	,
+			         		CrossSectionPath		);
 
    		ReactionFile.close();
 }
@@ -407,21 +405,19 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part
    G4int HeavyZ = m_Reaction->GetNucleus4()->GetZ() ;
    G4int HeavyA = m_Reaction->GetNucleus4()->GetA() ;
 
-   G4ParticleDefinition* HeavyName
-   = G4ParticleTable::GetParticleTable()->GetIon(HeavyZ, HeavyA, m_Reaction->GetExcitation()*MeV);
+
+	 // Shoot the Resonance energy following the mean and width value
+   	double EXX = -10 ;
+		EXX = RandBreitWigner::shoot(m_ResonanceMean,m_ResonanceWidth) ;
+
+  	m_Reaction->SetExcitation( EXX );
 
    // Beam
    G4int BeamZ = m_Reaction->GetNucleus1()->GetZ();
    G4int BeamA = m_Reaction->GetNucleus1()->GetA();
    G4ParticleDefinition* BeamName = G4ParticleTable::GetParticleTable()->GetIon(BeamZ, BeamA, 0);
 
-   // Shoot the Resonance energy following the mean and width value
-   double EXX = -10 ;
-   
-   while(EXX<0)
-   EXX = RandGauss::shoot(m_ResonanceMean,m_ResonanceWidth) ;
 
-	m_Reaction->SetExcitation( EXX  );
 
 	 ///////////////////////////////////////////////////////////////////////
    ///// Calculate the incident beam direction as well as the vertex /////
@@ -482,11 +478,6 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part
    ///// Angles for emitted particles following Cross Section //////
    ///// Angles are in the beam frame                         //////
    /////////////////////////////////////////////////////////////////
-   // Beam incident energy
-   G4double NominalBeamEnergy = m_BeamEnergy;
-   G4double IncidentBeamEnergy = RandGauss::shoot(NominalBeamEnergy, m_BeamEnergySpread);
-   m_Reaction->SetBeamEnergy(IncidentBeamEnergy);
-   m_InitConditions->SetICIncidentEnergy(IncidentBeamEnergy / MeV);
    // Angles
    RandGeneral CrossSectionShoot(m_Reaction->GetCrossSection(), m_Reaction->GetCrossSectionSize());
    G4double ThetaCM = CrossSectionShoot.shoot() * (180*deg);
@@ -545,12 +536,6 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part
    }
    
    // Case of recoil particle
-   /*   // Particle type
-      particleGun->SetParticleDefinition(HeavyName);
-      // Particle energy
-      particleGun->SetParticleEnergy(EnergyHeavy);
-      // Particle vertex position
-      particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0));
       // Particle direction
       // Kinematical angles in the beam frame are transformed 
       // to the "world" frame*/
@@ -560,14 +545,15 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part
       G4double phi_world   = momentum_kine_world.phi();
       if (phi_world < 1e-6) phi_world += 2*pi;
       
+      if(m_ShootHeavy || m_ShootDecayProduct)
       EventGeneratorTransfertToResonance::ResonanceDecay(	EnergyHeavy    	,
-      														theta_world     ,
-      														phi_world      	,
-      														x0            	,
-      														y0            	,
-      														z0             	,
-      														anEvent        	,
-      														particleGun		);
+      																										theta_world     ,
+												      														phi_world      	,
+												      														x0            	,
+												      														y0            	,
+												      														z0             	,
+												      														anEvent        	,
+												      														particleGun		  );
       
    
 }
@@ -598,7 +584,7 @@ void EventGeneratorTransfertToResonance::ResonanceDecay(  G4double EnergyHeavy
    else {
 		//Obtain Mass of daughter Nuclei
 		G4ParticleDefinition* parent
-		= G4ParticleTable::GetParticleTable()->GetIon(parentZ, parentA, m_Reaction->GetExcitation())     ;
+		= G4ParticleTable::GetParticleTable()->GetIon(parentZ, parentA, 0 )     ;
 
 		G4ParticleDefinition* daughter
 		= G4ParticleTable::GetParticleTable()->GetIon(m_ResonanceDecayZ, m_ResonanceDecayA, 0.) ;
@@ -609,7 +595,7 @@ void EventGeneratorTransfertToResonance::ResonanceDecay(  G4double EnergyHeavy
 		G4ParticleDefinition* proton
 		=  G4ParticleTable::GetParticleTable()->FindParticle("proton");
 
-		G4double M  = parent   -> GetPDGMass()     ;
+		G4double M  = parent   -> GetPDGMass() + m_Reaction->GetExcitation()     ;
 		G4double md = daughter -> GetPDGMass()     ;
 		G4double mn = neutron  -> GetPDGMass()     ;
 		G4double mp = proton   -> GetPDGMass()     ;
@@ -638,15 +624,9 @@ void EventGeneratorTransfertToResonance::ResonanceDecay(  G4double EnergyHeavy
 		TGenPhaseSpace TPhaseSpace ;
 
 		if( !TPhaseSpace.SetDecay(Initial, NumberOfDecayProducts+1, masses) ) cout << "Warning: Phase Space Decay forbiden by kinematic, or more than 18 particles "<<endl;
-		double MaxWt=TPhaseSpace.GetWtMax() ;
-		double Weight = 0 	;
-		double Rand   = 1	; 
-
-		while( Rand > Weight )
-			{  
-				Weight = TPhaseSpace.Generate() 		;
-				Rand = CLHEP::RandFlat::shoot()*MaxWt	; 
-			}
+
+		//	Generate an event and store is weight in the Output Tree
+		m_EventWeight = TPhaseSpace.Generate() 		;
 		
 		
 		TLorentzVector* daugterLV ;

NPSimulation/src/PrimaryGeneratorAction.cc

@@ -35,6 +35,7 @@
 #include "EventGeneratorTransfertToResonance.hh"
 #include "EventGeneratorIsotropic.hh"
 #include "EventGeneratorBeam.hh"
+#include "EventGeneratorPhaseSpace.hh"
 
 #include <cstdlib>
 
@@ -83,10 +84,11 @@ void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path)
    ifstream EventGeneratorFile;
    EventGeneratorFile.open(Path.c_str());
 
-   bool check_TransfertToResonance = false   ;
-   bool check_Isotropic        = false ;
-   bool check_Transfert        = false ;
-   bool check_Beam                = false ;
+   bool check_TransfertToResonance 	= false   ;
+   bool check_PhaseSpace				 		= false   ;
+   bool check_Isotropic        			= false 	;
+   bool check_Transfert        			= false 	;
+   bool check_Beam                	= false 	;
 
    if (EventGeneratorFile.is_open()) cout << " Event Generator file " << Path << " loading " << endl  ;
    else                      {
@@ -151,6 +153,18 @@ void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path)
          myEventGenerator->SetTarget(m_detector->GetTarget());
          m_EventGenerator = myEventGenerator                            ;
       }
+      
+      //Search for Transfert To Resonance
+      else if (LineBuffer.compare(0, 10, "PhaseSpace") == 0 && !check_PhaseSpace) {
+         check_PhaseSpace = true                             								 	;
+         VEventGenerator* myEventGenerator = new EventGeneratorPhaseSpace()  	;
+         EventGeneratorFile.close()                                  					;
+         myEventGenerator->ReadConfiguration(Path)                      			;	
+         EventGeneratorFile.open(Path.c_str())                          			;
+         myEventGenerator->InitializeRootOutput()           									;
+         myEventGenerator->SetTarget(m_detector->GetTarget())									;
+         m_EventGenerator = myEventGenerator                            			;
+      }
    }
    EventGeneratorFile.close();
 }