diff --git a/Inputs/DetectorConfiguration/Riken_65mm.detector b/Inputs/DetectorConfiguration/Riken_65mm.detector
index bee456b3d8bd5104bcb017232ee3a0cc5bc03615..af4c56547d29995c5fafbe03c332ddf04e3a9c64 100644
--- a/Inputs/DetectorConfiguration/Riken_65mm.detector
+++ b/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
diff --git a/Inputs/EventGenerator/10He.reaction b/Inputs/EventGenerator/10He.reaction
index 8567e32e3e2aad009e6f1f14e403149c23216e86..1b79c493cd70ded435c5fd5b141f9362efa3bb28 100644
--- a/Inputs/EventGenerator/10He.reaction
+++ b/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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Inputs/EventGenerator/34Si.reaction b/Inputs/EventGenerator/34Si.reaction
index edecfe2aa56337481e8871aedf2e9e364b0ce9cc..41a35ec8e0a24a92782a3bf5202d1700d4c628ac 100644
--- a/Inputs/EventGenerator/34Si.reaction
+++ b/Inputs/EventGenerator/34Si.reaction
@@ -1,21 +1,37 @@
+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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Inputs/EventGenerator/3He.source b/Inputs/EventGenerator/3He.source
index 8fa0b8698412b53f993a685d44e12862aabcb8ff..22957a57d0c16cf6c69813e071d726986534ad59 100644
--- a/Inputs/EventGenerator/3He.source
+++ b/Inputs/EventGenerator/3He.source
@@ -5,7 +5,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 0
- EnergyHigh= 21.6
+ EnergyHigh= 300
HalfOpenAngleMin= 0
HalfOpenAngleMax= 90
x0= 0
diff --git a/Inputs/EventGenerator/alpha.source b/Inputs/EventGenerator/alpha.source
index fd5f19e5259addc2b6bbccaeb2e2bcb97b08b3e1..66b3898dfd65a410b15cc81b6eaa35ba08d81aa8 100644
--- a/Inputs/EventGenerator/alpha.source
+++ b/Inputs/EventGenerator/alpha.source
@@ -5,7 +5,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 0
- EnergyHigh= 21.6
+ EnergyHigh= 300
HalfOpenAngleMin= 0
HalfOpenAngleMax= 90
x0= 0
diff --git a/NPAnalysis/10He_Riken/include/ObjectManager.hh b/NPAnalysis/10He_Riken/include/ObjectManager.hh
index c587cb506ac9b11d4a675b8660ffacc6dbc27679..fd9d4b26a378eb1bf7cdd13be6564d2fe29c458a 100644
--- a/NPAnalysis/10He_Riken/include/ObjectManager.hh
+++ b/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" ,
diff --git a/NPAnalysis/10He_Riken/macro/CrossSection.c b/NPAnalysis/10He_Riken/macro/CrossSection.c
index 5640efb446f5464a98882b5a40ea015ef2745171..56e071c820cdedaccfe71c6e3a1a7894070474cf 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 1c4954590a8caef6846113b1c18ffedd89d0b54b..0fe0d7c40e4efbc22e7d52b433013164740836da 100644
--- a/NPAnalysis/10He_Riken/src/Analysis.cc
+++ b/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() ;
diff --git a/NPSimulation/include/EventGeneratorTransfertToResonance.hh b/NPSimulation/include/EventGeneratorTransfertToResonance.hh
index 2a5f1c30f8c822ffae01140825a6fbd6c346b0fe..d42f08d0de2f4225662e8e7fab4a86486127f904 100644
--- a/NPSimulation/include/EventGeneratorTransfertToResonance.hh
+++ b/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 ,
diff --git a/NPSimulation/src/EventGeneratorTransfertToResonance.cc b/NPSimulation/src/EventGeneratorTransfertToResonance.cc
index dab91dec0c35eed7f27f43f8df34e9893adca48a..e09a894b23b660036f0a6b271e7fdf84f495c0d6 100644
--- a/NPSimulation/src/EventGeneratorTransfertToResonance.cc
+++ b/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 ;
diff --git a/NPSimulation/src/PrimaryGeneratorAction.cc b/NPSimulation/src/PrimaryGeneratorAction.cc
index a5d8a264e2eb54c33cb78404ac210f7ee0ab6ba4..ebae82b175efcd8037226d38ca9500495ddf3416 100644
--- a/NPSimulation/src/PrimaryGeneratorAction.cc
+++ b/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();
}