diff --git a/Inputs/EventGenerator/np.source b/Inputs/EventGenerator/np.source
new file mode 100644
index 0000000000000000000000000000000000000000..0fc3b7342f98e48010358a8560e1fd6ef7932c3d
--- /dev/null
+++ b/Inputs/EventGenerator/np.source
@@ -0,0 +1,17 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%% An Isotropic Source to be used as EventGenerator %%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Energy are given in MeV , Position in mm %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Isotropic
+ EnergyLow= 10
+ EnergyHigh= 250
+ HalfOpenAngleMin= 0
+ HalfOpenAngleMax= 10
+ x0= 0
+ y0= 0
+ z0= 0
+ Multiplicity= 1 1
+ Particle= neutron proton
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Supported particle type: proton, neutron, deuton, triton, He3 , alpha
diff --git a/NPSimulation/Core/EventGeneratorIsotropic.cc b/NPSimulation/Core/EventGeneratorIsotropic.cc
index 2f83bb9b000dc61d519c1e60dc75a5c9d3ae7a5f..0523ff8f064682c371304f184053dca169fbd450 100644
--- a/NPSimulation/Core/EventGeneratorIsotropic.cc
+++ b/NPSimulation/Core/EventGeneratorIsotropic.cc
@@ -47,7 +47,7 @@ EventGeneratorIsotropic::EventGeneratorIsotropic(){
m_z0 = 0 ;
m_SigmaX = 0;
m_SigmaY = 0;
- m_Multiplicty = 1;
+ m_Multiplicty.clear();
m_particle = NULL;
m_ParticleStack = ParticleStack::getInstance();
m_ExcitationEnergy = 0 ;
@@ -66,6 +66,7 @@ void EventGeneratorIsotropic::ReadConfiguration(string Path,int){
////////General Reading needs////////
string LineBuffer;
string DataBuffer;
+ istringstream LineStream;
bool ReadingStatus = false ;
bool check_EnergyLow = false ;
@@ -170,25 +171,45 @@ void EventGeneratorIsotropic::ReadConfiguration(string Path,int){
else if (DataBuffer == "Multiplicity=" ) {
check_Multiplicity = true;
- ReactionFile >> DataBuffer;
- m_Multiplicty = atof(DataBuffer.c_str()) * mm;
- G4cout << "Multiplicity= " << m_Multiplicty << " " << G4endl;
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+
+ getline(ReactionFile, LineBuffer);
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ G4cout << " Multiplicity: " ;
+ while(LineStream >> DataBuffer){
+ m_Multiplicty.push_back(atof(DataBuffer.c_str()));
+ G4cout << DataBuffer << " " ;
+ }
+ G4cout << G4endl;
+ //ReactionFile >> DataBuffer;
+ //m_Multiplicty = atof(DataBuffer.c_str());
+ //G4cout << "Multiplicity= " << m_Multiplicty << " " << G4endl;
}
else if (DataBuffer=="Particle=" || DataBuffer=="particle=") {
check_particle = true ;
- ReactionFile >> m_particleName;
- G4cout << "Particle : " << m_particleName << G4endl ;
-
- // Case of light particle
- if(m_particleName=="proton"){ m_particleName="1H" ; check_ExcitationEnergy = true ;}
- else if(m_particleName=="deuton"){ m_particleName="2H" ; check_ExcitationEnergy = true ;}
- else if(m_particleName=="triton"){ m_particleName="3H" ; check_ExcitationEnergy = true ;}
- else if(m_particleName=="alpha") { m_particleName="4He" ; check_ExcitationEnergy = true ;}
- else if(m_particleName=="gamma") { check_ExcitationEnergy = true ;}
- else if(m_particleName=="neutron") { check_ExcitationEnergy = true ;}
- else { check_ExcitationEnergy = true ;}
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ getline(ReactionFile, LineBuffer);
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ G4cout << " Particle: " ;
+ while(LineStream >> DataBuffer){
+ //m_particleName.push_back(DataBuffer);
+ G4cout << DataBuffer << " " ;
+ // Case of light particle
+ if(DataBuffer=="proton"){ m_particleName.push_back("1H") ; check_ExcitationEnergy = true ;}
+ else if(DataBuffer=="deuton"){ m_particleName.push_back("2H") ; check_ExcitationEnergy = true ;}
+ else if(DataBuffer=="triton"){ m_particleName.push_back("3H") ; check_ExcitationEnergy = true ;}
+ else if(DataBuffer=="alpha") { m_particleName.push_back("4He") ; check_ExcitationEnergy = true ;}
+ else if(DataBuffer=="gamma") { m_particleName.push_back("gamma") ; check_ExcitationEnergy = true ;}
+ else if(DataBuffer=="neutron") {m_particleName.push_back("neutron") ; check_ExcitationEnergy = true ;}
+ else { check_ExcitationEnergy = true ;}
+ }
+ G4cout << G4endl;
}
else if (DataBuffer=="ExcitationEnergy=") {
@@ -202,6 +223,7 @@ void EventGeneratorIsotropic::ReadConfiguration(string Path,int){
else
{ReadingStatus = false; G4cout << "WARNING : Wrong Token Sequence: Getting out " << G4endl ;}
+
///////////////////////////////////////////////////
// If all Token found toggle out
if( check_EnergyLow && check_EnergyHigh && check_HalfOpenAngleMin && check_HalfOpenAngleMax && check_x0 && check_y0 && check_z0 && check_particle && check_ExcitationEnergy)
@@ -210,12 +232,15 @@ void EventGeneratorIsotropic::ReadConfiguration(string Path,int){
}
}
+ if(check_Multiplicity==false) {m_Multiplicty.push_back(1);check_Multiplicity=true;}
+
+ if(m_Multiplicty.size()!=m_particleName.size()) {
+ G4cout << "WARNING : Size of m_particleName and m_Multiplicty are different" << G4endl ;
+ G4cout << "Assuming m_Multiplicty=1 for all particles" << G4endl ;
+ }
if( !check_EnergyLow || !check_EnergyHigh || !check_HalfOpenAngleMin || !check_HalfOpenAngleMax || !check_x0 || !check_y0 || !check_z0 || !check_particle)
{cout << "ERROR : Token Sequence Incomplete, Isotropic definition can not be Fonctionnal" << G4endl ; exit(1);}
-
-
-
}
@@ -223,38 +248,42 @@ void EventGeneratorIsotropic::ReadConfiguration(string Path,int){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorIsotropic::GenerateEvent(G4Event*){
- for(int i=0; i<m_Multiplicty; i++){
- if(m_particle==NULL){
- if(m_particleName=="gamma" || m_particleName=="neutron" || m_particleName=="opticalphoton"){
- m_particle = G4ParticleTable::GetParticleTable()->FindParticle(m_particleName.c_str());
- }
- else{
- NPL::Nucleus* N = new NPL::Nucleus(m_particleName);
- m_particle = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(N->GetZ(), N->GetA(),m_ExcitationEnergy);
- delete N;
+ for(unsigned int p=0; p<m_particleName.size(); p++){
+ for(int i=0; i<m_Multiplicty[p]; i++){
+ m_particle=NULL;
+ if(m_particle==NULL){
+
+ if(m_particleName[p]=="gamma" || m_particleName[p]=="neutron" || m_particleName[p]=="opticalphoton"){
+ m_particle = G4ParticleTable::GetParticleTable()->FindParticle(m_particleName[p].c_str());
+ }
+ else{
+ NPL::Nucleus* N = new NPL::Nucleus(m_particleName[p]);
+ m_particle = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(N->GetZ(), N->GetA(),m_ExcitationEnergy);
+ delete N;
+ }
+
}
+ G4double cos_theta_min = cos(m_HalfOpenAngleMin);
+ G4double cos_theta_max = cos(m_HalfOpenAngleMax);
+ G4double cos_theta = cos_theta_min + (cos_theta_max - cos_theta_min) * RandFlat::shoot();
+ G4double theta = acos(cos_theta) ;
+ G4double phi = RandFlat::shoot() * 2 * pi ;
+ G4double particle_energy = m_EnergyLow + RandFlat::shoot() * (m_EnergyHigh - m_EnergyLow) ;
+
+ // Direction of particle, energy and laboratory angle
+ G4double momentum_x = sin(theta) * cos(phi) ;
+ G4double momentum_y = sin(theta) * sin(phi) ;
+ G4double momentum_z = cos(theta) ;
+
+ G4double x0 = RandGauss::shoot(m_x0,m_SigmaX);
+ G4double y0 = RandGauss::shoot(m_y0,m_SigmaY);
+
+ Particle particle(m_particle, theta,particle_energy,G4ThreeVector(momentum_x, momentum_y, momentum_z),G4ThreeVector(x0, y0, m_z0));
+
+
+ m_ParticleStack->AddParticleToStack(particle);
}
-
- G4double cos_theta_min = cos(m_HalfOpenAngleMin);
- G4double cos_theta_max = cos(m_HalfOpenAngleMax);
- G4double cos_theta = cos_theta_min + (cos_theta_max - cos_theta_min) * RandFlat::shoot();
- G4double theta = acos(cos_theta) ;
- G4double phi = RandFlat::shoot() * 2 * pi ;
- G4double particle_energy = m_EnergyLow + RandFlat::shoot() * (m_EnergyHigh - m_EnergyLow) ;
-
- // Direction of particle, energy and laboratory angle
- G4double momentum_x = sin(theta) * cos(phi) ;
- G4double momentum_y = sin(theta) * sin(phi) ;
- G4double momentum_z = cos(theta) ;
-
- G4double x0 = RandGauss::shoot(m_x0,m_SigmaX);
- G4double y0 = RandGauss::shoot(m_y0,m_SigmaY);
-
- Particle particle(m_particle, theta,particle_energy,G4ThreeVector(momentum_x, momentum_y, momentum_z),G4ThreeVector(x0, y0, m_z0));
-
-
- m_ParticleStack->AddParticleToStack(particle);
}
}
diff --git a/NPSimulation/Core/EventGeneratorIsotropic.hh b/NPSimulation/Core/EventGeneratorIsotropic.hh
index 5c5e9514ecbcb59d509a69b1906f4463e1073673..19870fb82172d651856d5c6efe2027bf3739df59 100644
--- a/NPSimulation/Core/EventGeneratorIsotropic.hh
+++ b/NPSimulation/Core/EventGeneratorIsotropic.hh
@@ -23,6 +23,7 @@
*****************************************************************************/
// C++ header
#include <string>
+#include <cmath>
using namespace std;
using namespace CLHEP;
@@ -56,9 +57,9 @@ private: // Source parameter from input file
G4double m_SigmaY ;
G4ParticleDefinition* m_particle ; // Kind of particle to shoot isotropically
G4double m_ExcitationEnergy ; // Excitation energy of the emitted particle
- string m_particleName ;
+ vector<string> m_particleName ;
ParticleStack* m_ParticleStack ;
- G4int m_Multiplicty;//Adding by Pierre
+ vector<G4int> m_Multiplicty;
};
#endif