Skip to content
Snippets Groups Projects
Commit ae7a0a1b authored by Morfouace's avatar Morfouace
Browse files

* Update of EventGeneratorIsotropic to be able to simulate different 

* particle for the same event. Example of token is in
* Inputs/EventGenerator/np.source No need to modify token of
* proton.source or deutron.source... It reads m_particleName as a
* vector<string> now and loop over the size of m_particleName
parent 5e684ffa
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Inputs/EventGenerator/np.source 0 → 100644
+ 17
− 0
View file @ ae7a0a1b
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% 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
NPSimulation/Core/EventGeneratorIsotropic.cc
+ 76
− 47
View file @ ae7a0a1b
......@@ -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);
}
}
......
NPSimulation/Core/EventGeneratorIsotropic.hh
+ 3
− 2
View file @ ae7a0a1b
......@@ -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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment