diff --git a/Inputs/EventGenerator/10He.reaction b/Inputs/EventGenerator/10He.reaction
index 4e7f72f0d0faee9d2caef45ea37415efe0d8cd2a..4c7444abc5b71ac51da2a59bfcd031e1ca19e823 100644
--- a/Inputs/EventGenerator/10He.reaction
+++ b/Inputs/EventGenerator/10He.reaction
@@ -1,7 +1,18 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% Reaction file for 11Li(d,3He)10He reaction %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-Beam
+Isotropic
+ EnergyLow= 100
+ EnergyHigh= 100
+ HalfOpenAngleMin= 0
+ HalfOpenAngleMax= 90
+ x0= 0
+ y0= 0
+ z0= 0
+ Particle= 10He
+ ExcitationEnergy= 1.4
+
+%Beam
Particle= 11Li
Energy= 550
SigmaEnergy= 0
@@ -19,7 +30,7 @@ Beam
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-TwoBodyReaction
+%TwoBodyReaction
Beam= 11Li
Target= 2H
Light= 3He
@@ -32,9 +43,9 @@ TwoBodyReaction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%ParticleDecay 10He
+ParticleDecay 10He
Daughter= 9He n
- ExcitationEnergy= 0.5 0
+ ExcitationEnergy= 0 0
DifferentialCrossSection= flat.txt
shoot= 1 1
diff --git a/NPLib/InitialConditions/TInitialConditions.cxx b/NPLib/InitialConditions/TInitialConditions.cxx
index f19eb1df2ecd6f8b4d6d5b7bf68806128a2fee29..1c39613f876a5bf98bc444e1b65f16b0be49df86 100644
--- a/NPLib/InitialConditions/TInitialConditions.cxx
+++ b/NPLib/InitialConditions/TInitialConditions.cxx
@@ -53,7 +53,6 @@ void TInitialConditions::Clear(){
// emmitted particle
fIC_Particle_Name.clear();
- fIC_Process_Name.clear();
fIC_ThetaCM.clear();
fIC_Kinetic_Energy.clear();
fIC_Momentum_Direction_X.clear();
@@ -87,7 +86,6 @@ void TInitialConditions::Dump() const{
for(unsigned int i = 0 ; i < size; i ++){
cout << "\t ---- Particle " << i << " ---- " << endl;
cout << "\t Particle Name" << fIC_Particle_Name[i] << endl;
- // cout << "\t Process Name" << fIC_Process_Name[i] << endl;
cout << "\t Theta CM" << fIC_ThetaCM[i] << endl;
cout << "\t Energy" << fIC_Kinetic_Energy[i] << endl;
cout << "\t Momentum Direction: ( "
diff --git a/NPLib/InitialConditions/TInitialConditions.h b/NPLib/InitialConditions/TInitialConditions.h
index 4d8674cdaf846482ea36796089c474c8e872e8a7..b41dde71d863bd90f5fbd88ea2c1d2d5d19b681c 100644
--- a/NPLib/InitialConditions/TInitialConditions.h
+++ b/NPLib/InitialConditions/TInitialConditions.h
@@ -52,7 +52,6 @@ private:
// emmitted particle
vector<string> fIC_Particle_Name;
- vector<string> fIC_Process_Name;
vector<double> fIC_ThetaCM;
vector<double> fIC_Kinetic_Energy;
vector<double> fIC_Momentum_Direction_X;
@@ -82,7 +81,6 @@ public:
// emmitted particle
void SetParticleName (string Particle_Name) {fIC_Particle_Name.push_back(Particle_Name);}
- void SetProcessName (string Process_Name) {fIC_Process_Name.push_back(Process_Name);}
void SetThetaCM (double ThetaCM) {fIC_ThetaCM.push_back(ThetaCM);}
void SetKineticEnergy (double Kinetic_Energy) {fIC_Kinetic_Energy.push_back(Kinetic_Energy);}
void SetMomentumDirectionX (double Momentum_Direction_X) {fIC_Momentum_Direction_X.push_back(Momentum_Direction_X);}
@@ -103,7 +101,6 @@ public:
// emmitted particle
string GetParticleName (int i) const {return fIC_Particle_Name[i];}
- string GetProcessName (int i) const {return fIC_Process_Name[i];}
double GetThetaCM (int i) const {return fIC_ThetaCM[i];}
double GetKineticEnergy (int i) const {return fIC_Kinetic_Energy[i];}
double GetMomentumDirectionX (int i) const {return fIC_Momentum_Direction_X[i];}
diff --git a/NPSimulation/include/EventGeneratorIsotropic.hh b/NPSimulation/include/EventGeneratorIsotropic.hh
index f370a3ec4e4a8c936a4b5bd5bdf1d8a7c235299e..2da702dbf920b2490fa6e0e837f7e92102779117 100644
--- a/NPSimulation/include/EventGeneratorIsotropic.hh
+++ b/NPSimulation/include/EventGeneratorIsotropic.hh
@@ -23,37 +23,37 @@
*****************************************************************************/
// C++ header
#include <string>
+using namespace std;
// G4 headers
#include "G4Event.hh"
-#include "G4ParticleGun.hh"
-// NPTool headers
+// NPS headers
#include "VEventGenerator.hh"
-#include "TInitialConditions.h"
-
-using namespace std;
+#include "ParticleStack.hh"
-
-class EventGeneratorIsotropic : public VEventGenerator
-{
+class EventGeneratorIsotropic : public VEventGenerator{
public: // Constructor and destructor
- EventGeneratorIsotropic() ;
- virtual ~EventGeneratorIsotropic();
-
+ EventGeneratorIsotropic() ;
+ virtual ~EventGeneratorIsotropic();
+
public: // Inherit from VEventGenerator Class
- void ReadConfiguration(string) ;
- void GenerateEvent(G4Event*, G4ParticleGun*) ;
- void InitializeRootOutput() ;
-
+ void ReadConfiguration(string,int dump=0) ;
+ void GenerateEvent(G4Event*) ;
+ void InitializeRootOutput() ;
+
private: // Source parameter from input file
- G4double m_EnergyLow ; // Lower limit of energy range
- G4double m_EnergyHigh ; // Upper limit of energy range
- G4double m_HalfOpenAngleMin ; // Min Half open angle of the source
- G4double m_HalfOpenAngleMax ; // Max Half open angle of the source
- G4double m_x0 ; // Vertex Position X
- G4double m_y0 ; // Vertex Position Y
- G4double m_z0 ; // Vertex Position Z
- G4ParticleDefinition* m_particle ; // Kind of particle to shoot isotropically
+ G4double m_EnergyLow ; // Lower limit of energy range
+ G4double m_EnergyHigh ; // Upper limit of energy range
+ G4double m_HalfOpenAngleMin ; // Min Half open angle of the source
+ G4double m_HalfOpenAngleMax ; // Max Half open angle of the source
+ G4double m_x0 ; // Vertex Position X
+ G4double m_y0 ; // Vertex Position Y
+ G4double m_z0 ; // Vertex Position Z
+ G4ParticleDefinition* m_particle ; // Kind of particle to shoot isotropically
+ G4double m_ExcitationEnergy ; // Excitation energy of the emitted particle
+ string m_particleName ;
+ ParticleStack* m_ParticleStack ;
+
};
#endif
diff --git a/NPSimulation/src/EventGeneratorIsotropic.cc b/NPSimulation/src/EventGeneratorIsotropic.cc
index a1f1d714a77998df4795b034c7f5849422fbff56..e4a3d39f3f16f34ff024779430a954d760a95586 100644
--- a/NPSimulation/src/EventGeneratorIsotropic.cc
+++ b/NPSimulation/src/EventGeneratorIsotropic.cc
@@ -29,191 +29,203 @@
// for generating random numbers
#include "Randomize.hh"
-// NPTool headers
+// NPS headers
#include "EventGeneratorIsotropic.hh"
-#include "RootOutput.h"
-using namespace CLHEP;
+// NPl headers
+#include "RootOutput.h"
+#include "NPNucleus.h"
+using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-EventGeneratorIsotropic::EventGeneratorIsotropic()
-{
-
- m_EnergyLow = 0 ;
- m_EnergyHigh = 0 ;
- m_x0 = 0 ;
- m_y0 = 0 ;
- m_z0 = 0 ;
- m_particle = G4ParticleTable::GetParticleTable()->FindParticle("geantino");
+EventGeneratorIsotropic::EventGeneratorIsotropic(){
+ m_EnergyLow = 0 ;
+ m_EnergyHigh = 0 ;
+ m_x0 = 0 ;
+ m_y0 = 0 ;
+ m_z0 = 0 ;
+ m_particle = NULL;
+ m_ParticleStack = ParticleStack::getInstance();
+ m_ExcitationEnergy = 0 ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-EventGeneratorIsotropic::~EventGeneratorIsotropic()
-{
+EventGeneratorIsotropic::~EventGeneratorIsotropic(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorIsotropic::ReadConfiguration(string Path)
-{
- ////////General Reading needs////////
- string LineBuffer;
- string DataBuffer;
-
- bool ReadingStatus = false ;
- bool check_EnergyLow = false ;
- bool check_EnergyHigh = false ;
- bool check_HalfOpenAngleMin = false ;
- bool check_HalfOpenAngleMax = false ;
- bool check_x0 = false ;
- bool check_y0 = false ;
- bool check_z0 = false ;
- bool check_particle = false ;
-
- ////////Reaction Setting needs///////
- string particle ;
- //////////////////////////////////////////////////////////////////////////////////////////
- ifstream ReactionFile;
- ReactionFile.open(Path.c_str());
-
- if (ReactionFile.is_open()) {}
- else {
- return;
- }
-
- while (!ReactionFile.eof()) {
- //Pick-up next line
- getline(ReactionFile, LineBuffer);
-
- if (LineBuffer.compare(0, 9, "Isotropic") == 0) {
- G4cout << "Isotropic Source Found" << G4endl ;
- ReadingStatus = true;}
-
-
- while (ReadingStatus)
- {
- ReactionFile >> DataBuffer;
-
- //Search for comment Symbol %
- if (DataBuffer.compare(0, 1, "%") == 0) { ReactionFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );}
-
- else if (DataBuffer.compare(0, 10, "EnergyLow=") == 0) {
- check_EnergyLow = true ;
- ReactionFile >> DataBuffer;
- m_EnergyLow = atof(DataBuffer.c_str()) * MeV;
- G4cout << "Minimum energy " << m_EnergyLow / MeV << " MeV" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 11, "EnergyHigh=") == 0) {
- check_EnergyHigh = true ;
- ReactionFile >> DataBuffer;
- m_EnergyHigh = atof(DataBuffer.c_str()) * MeV;
- G4cout << "Maximum energy " << m_EnergyHigh / MeV << " MeV" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 17, "HalfOpenAngleMin=") == 0) {
- check_HalfOpenAngleMin = true ;
- ReactionFile >> DataBuffer;
- m_HalfOpenAngleMin = atof(DataBuffer.c_str()) * deg;
- G4cout << "HalfOpenAngleMin " << m_HalfOpenAngleMin / deg << " degree" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 17, "HalfOpenAngleMax=") == 0) {
- check_HalfOpenAngleMax = true ;
- ReactionFile >> DataBuffer;
- m_HalfOpenAngleMax = atof(DataBuffer.c_str()) * deg;
- G4cout << "HalfOpenAngleMax " << m_HalfOpenAngleMax / deg << " degree" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 3, "x0=") == 0) {
- check_x0 = true ;
- ReactionFile >> DataBuffer;
- m_x0 = atof(DataBuffer.c_str()) * mm;
- G4cout << "x0 " << m_x0 << " mm" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 3, "y0=") == 0) {
- check_y0 = true ;
- ReactionFile >> DataBuffer;
- m_y0 = atof(DataBuffer.c_str()) * mm;
- G4cout << "y0 " << m_y0 << " mm" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 3, "z0=") == 0) {
- check_z0 = true ;
- ReactionFile >> DataBuffer;
- m_y0 = atof(DataBuffer.c_str()) * mm;
- G4cout << "y0 " << m_y0 << " mm" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 9, "particle=") == 0) {
- check_particle = true ;
- ReactionFile >> DataBuffer;
- particle = DataBuffer;
- G4cout << "particle " << particle << G4endl;
- m_particle = G4ParticleTable::GetParticleTable()->FindParticle(particle);
- }
-
- // If no isotropic Token and no comment, toggle out
- 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 )
- ReadingStatus = false ;
-
- }
-
+void EventGeneratorIsotropic::ReadConfiguration(string Path,int dump){
+ dump= 0;
+ ////////General Reading needs////////
+ string LineBuffer;
+ string DataBuffer;
+
+ bool ReadingStatus = false ;
+ bool check_EnergyLow = false ;
+ bool check_EnergyHigh = false ;
+ bool check_HalfOpenAngleMin = false ;
+ bool check_HalfOpenAngleMax = false ;
+ bool check_x0 = false ;
+ bool check_y0 = false ;
+ bool check_z0 = false ;
+ bool check_particle = false ;
+ bool check_ExcitationEnergy = false ;
+
+ ////////Reaction Setting needs///////
+ string particle ;
+ //////////////////////////////////////////////////////////////////////////////////////////
+ ifstream ReactionFile;
+ ReactionFile.open(Path.c_str());
+
+ if (ReactionFile.is_open()) {}
+ else {
+ return;
+ }
+
+ while (!ReactionFile.eof()) {
+ //Pick-up next line
+ getline(ReactionFile, LineBuffer);
+
+ if (LineBuffer.compare(0, 9, "Isotropic") == 0) {
+ G4cout << "///////////////////////////////////////////////////" << G4endl ;
+ G4cout << "Isotropic Source Found" << G4endl ;
+ ReadingStatus = true;}
+
+
+ while (ReadingStatus)
+ {
+ ReactionFile >> DataBuffer;
+
+ //Search for comment Symbol %
+ if (DataBuffer.compare(0, 1, "%") == 0) { ReactionFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );}
+
+ else if (DataBuffer == "EnergyLow=") {
+ check_EnergyLow = true ;
+ ReactionFile >> DataBuffer;
+ m_EnergyLow = atof(DataBuffer.c_str()) * MeV;
+ G4cout << "Minimum energy " << m_EnergyLow / MeV << " MeV" << G4endl;
+ }
+
+ else if (DataBuffer == "EnergyHigh=") {
+ check_EnergyHigh = true ;
+ ReactionFile >> DataBuffer;
+ m_EnergyHigh = atof(DataBuffer.c_str()) * MeV;
+ G4cout << "Maximum energy " << m_EnergyHigh / MeV << " MeV" << G4endl;
+ }
+
+ else if (DataBuffer == "HalfOpenAngleMin=") {
+ check_HalfOpenAngleMin = true ;
+ ReactionFile >> DataBuffer;
+ m_HalfOpenAngleMin = atof(DataBuffer.c_str()) * deg;
+ G4cout << "HalfOpenAngleMin " << m_HalfOpenAngleMin / deg << " degree" << G4endl;
+ }
+
+ else if (DataBuffer == "HalfOpenAngleMax=") {
+ check_HalfOpenAngleMax = true ;
+ ReactionFile >> DataBuffer;
+ m_HalfOpenAngleMax = atof(DataBuffer.c_str()) * deg;
+ G4cout << "HalfOpenAngleMax " << m_HalfOpenAngleMax / deg << " degree" << G4endl;
}
-
- if( !check_EnergyLow || !check_EnergyHigh || !check_HalfOpenAngleMin || !check_HalfOpenAngleMax || !check_x0 || !check_y0 || !check_z0 || !check_particle )
- {cout << "WARNING : Token Sequence Incomplete, Isotropic definition could not be Fonctionnal" << endl ;}
-
- cout << "///////////////////////////////////////////////////" << endl << endl ;
-
+
+ else if (DataBuffer == "x0=") {
+ check_x0 = true ;
+ ReactionFile >> DataBuffer;
+ m_x0 = atof(DataBuffer.c_str()) * mm;
+ G4cout << "x0 " << m_x0 << " mm" << G4endl;
+ }
+
+ else if (DataBuffer == "y0=") {
+ check_y0 = true ;
+ ReactionFile >> DataBuffer;
+ m_y0 = atof(DataBuffer.c_str()) * mm;
+ G4cout << "y0 " << m_y0 << " mm" << G4endl;
+ }
+
+ else if (DataBuffer == "z0=" ) {
+ check_z0 = true ;
+ ReactionFile >> DataBuffer;
+ m_z0 = atof(DataBuffer.c_str()) * mm;
+ G4cout << "z0 " << m_z0 << " mm" << G4endl;
+ }
+
+ else if (DataBuffer=="Particle=") {
+ check_particle = true ;
+ ReactionFile >> m_particleName;
+ G4cout << "Particle : " << m_particleName << endl ;
+
+ }
+
+ else if (DataBuffer=="ExcitationEnergy=") {
+ check_ExcitationEnergy = true ;
+ ReactionFile >> DataBuffer;
+ m_ExcitationEnergy = atof(DataBuffer.c_str()) * MeV;
+ G4cout << "ExcitationEnergy : " << m_ExcitationEnergy << endl ;
+ }
+
+ // If no isotropic Token and no comment, toggle out
+ 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)
+ ReadingStatus = false ;
+
+ }
+
+ }
+
+ if( !check_EnergyLow || !check_EnergyHigh || !check_HalfOpenAngleMin || !check_HalfOpenAngleMax || !check_x0 || !check_y0 || !check_z0 || !check_particle )
+ {cout << "WARNING : Token Sequence Incomplete, Isotropic definition could not be Fonctionnal" << endl ;}
+
+
+
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorIsotropic::GenerateEvent(G4Event* anEvent, G4ParticleGun* particleGun)
-{
- // Clear TInitialConditions
-
- particleGun->SetParticleDefinition(m_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 = 3*pi/2; //RandFlat::shoot() * 2 * pi ;
- 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) ;
-
- // Set the gun to shoot
- particleGun->SetParticleMomentumDirection(G4ThreeVector(momentum_x, momentum_y, momentum_z)) ;
- particleGun->SetParticleEnergy(particle_energy) ;
- particleGun->SetParticlePosition(G4ThreeVector(m_x0, m_y0, m_z0)) ;
-
- //Shoot particle
- particleGun->GeneratePrimaryVertex(anEvent) ;
+void EventGeneratorIsotropic::GenerateEvent(G4Event* anEvent){
+
+ if(m_particle==NULL){
+ if(m_particleName!="gamma"){
+ NPL::Nucleus* N = new NPL::Nucleus(m_particleName);
+ m_particle = G4ParticleTable::GetParticleTable()->GetIon(N->GetZ(), N->GetA(),m_ExcitationEnergy);
+ delete N;
+ }
+ else{
+ m_particle = G4ParticleTable::GetParticleTable()->FindParticle("gamma");
+ }
+
+ }
+
+ // Clear TInitialConditions
+ 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) ;
+
+ Particle particle(m_particle, theta,particle_energy,G4ThreeVector(momentum_x, momentum_y, momentum_z),G4ThreeVector(m_x0, m_y0, m_z0));
+
+ m_ParticleStack->AddParticleToStack(particle);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorIsotropic::InitializeRootOutput()
-{
-
+void EventGeneratorIsotropic::InitializeRootOutput(){
+
}