diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index e081c2958d2390dc24d811c44a48cd79fc78305f..52037c6d0e4a14d28b8359d4b1b982be482688e5 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -35,6 +35,7 @@
#include "NPBeam.h"
#include "NPInputParser.h"
#include "NPParticle.h"
+
using namespace NPL;
// ROOT header
@@ -86,9 +87,9 @@ namespace NPL {
private:
Beam fParticle1; // Beam
- Particle fParticle2; // Target
- Particle fParticle3; // Light ejectile
- Particle fParticle4; // Heavy ejectile
+ NPL::Particle fParticle2; // Target
+ NPL::Particle fParticle3; // Light ejectile
+ NPL::Particle fParticle4; // Heavy ejectile
double fQValue; // Q-value in MeV
double fEcm; // Ecm in MeV
double fBeamEnergy; // Beam energy in MeV
@@ -156,14 +157,14 @@ namespace NPL {
double GetExcitation4() const { return fExcitation4; }
double GetQValue() const { return fQValue; }
double GetEcm() const { return fEcm; }
- Particle* GetParticle1() { return &fParticle1; }
- Particle* GetParticle2() { return &fParticle2; }
- Particle* GetParticle3() { return &fParticle3; }
- Particle* GetParticle4() { return &fParticle4; }
- Particle* GetNucleus1() { return GetParticle1(); }
- Particle* GetNucleus2() { return GetParticle2(); }
- Particle* GetNucleus3() { return GetParticle3(); }
- Particle* GetNucleus4() { return GetParticle4(); }
+ NPL::Particle* GetParticle1() { return &fParticle1; }
+ NPL::Particle* GetParticle2() { return &fParticle2; }
+ NPL::Particle* GetParticle3() { return &fParticle3; }
+ NPL::Particle* GetParticle4() { return &fParticle4; }
+ NPL::Particle* GetNucleus1() { return GetParticle1(); }
+ NPL::Particle* GetNucleus2() { return GetParticle2(); }
+ NPL::Particle* GetNucleus3() { return GetParticle3(); }
+ NPL::Particle* GetNucleus4() { return GetParticle4(); }
TH1D* GetCrossSectionHist() const { return fCrossSectionHist; }
int GetVerboseLevel() const { return fVerboseLevel; }
diff --git a/NPSimulation/Core/Particle.cc b/NPSimulation/Core/Particle.cc
index 99bd7ae4771025178bdaa89fd79481c4345776a8..1915b4f32336d908986f398f3debbc5056260875 100644
--- a/NPSimulation/Core/Particle.cc
+++ b/NPSimulation/Core/Particle.cc
@@ -23,6 +23,7 @@
#include "Particle.hh"
+using namespace NPS;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Particle::Particle(){
m_ParticleDefinition = NULL;
diff --git a/NPSimulation/Core/Particle.hh b/NPSimulation/Core/Particle.hh
index 764f1ffd8f5a0fc0ecf8aefde46083c208e60e5d..f7e46bbc2bc2369bcaa246bfbf5ed939f853ab43 100644
--- a/NPSimulation/Core/Particle.hh
+++ b/NPSimulation/Core/Particle.hh
@@ -28,37 +28,39 @@
#include"G4ThreeVector.hh"
#include"CLHEP/Units/SystemOfUnits.h"
using namespace CLHEP;
-class Particle{
-
-public: // Constructor and Destructor
- // Empty constructor (return geantino at zero degree)
- Particle();
- ~Particle();
- // Constructor to be used
- Particle(G4ParticleDefinition* particle,double ThetaCM,double T,G4ThreeVector Direction, G4ThreeVector Position,bool ShootStatus=true);
-
-private: // Private Member
- G4ParticleDefinition* m_ParticleDefinition;
- double m_ThetaCM;
- double m_T ;
- G4ThreeVector m_Direction;
- G4ThreeVector m_Position;
- bool m_ShootStatus;
-
-public: // Setter and Getter
- // Getter
- G4ParticleDefinition* GetParticleDefinition();
- double GetParticleThetaCM();
- double GetParticleKineticEnergy();
- G4ThreeVector GetParticleMomentumDirection();
- G4ThreeVector GetParticlePosition();
- bool GetShootStatus();
- // Setter
- void SetParticleDefinition(G4ParticleDefinition*);
- void SetParticleThetaCM(double);
- void SetParticleKineticEnergy(double);
- void SetParticlePosition(G4ThreeVector);
- void SetParticleMomentumDirection(G4ThreeVector);
- void SetShootStatus(bool);
-};
+namespace NPS{
+ class Particle{
+
+ public: // Constructor and Destructor
+ // Empty constructor (return geantino at zero degree)
+ Particle();
+ ~Particle();
+ // Constructor to be used
+ Particle(G4ParticleDefinition* particle,double ThetaCM,double T,G4ThreeVector Direction, G4ThreeVector Position,bool ShootStatus=true);
+
+ private: // Private Member
+ G4ParticleDefinition* m_ParticleDefinition;
+ double m_ThetaCM;
+ double m_T ;
+ G4ThreeVector m_Direction;
+ G4ThreeVector m_Position;
+ bool m_ShootStatus;
+
+ public: // Setter and Getter
+ // Getter
+ G4ParticleDefinition* GetParticleDefinition();
+ double GetParticleThetaCM();
+ double GetParticleKineticEnergy();
+ G4ThreeVector GetParticleMomentumDirection();
+ G4ThreeVector GetParticlePosition();
+ bool GetShootStatus();
+ // Setter
+ void SetParticleDefinition(G4ParticleDefinition*);
+ void SetParticleThetaCM(double);
+ void SetParticleKineticEnergy(double);
+ void SetParticlePosition(G4ThreeVector);
+ void SetParticleMomentumDirection(G4ThreeVector);
+ void SetShootStatus(bool);
+ };
+}
#endif
diff --git a/NPSimulation/Core/ParticleStack.cc b/NPSimulation/Core/ParticleStack.cc
index 9c30eba23af31d0ba6b9cf464f1d1ae55b4d26b8..d71dc8767705f86b1a87dd5666ede2287623b858 100644
--- a/NPSimulation/Core/ParticleStack.cc
+++ b/NPSimulation/Core/ParticleStack.cc
@@ -76,17 +76,17 @@ void ParticleStack::AttachInitialConditions(){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-vector<Particle> ParticleStack::GetParticleStack(){
+vector<NPS::Particle> ParticleStack::GetParticleStack(){
return m_ParticleStack;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void ParticleStack::SetParticleStack(vector<Particle> particle_stack){
+void ParticleStack::SetParticleStack(vector<NPS::Particle> particle_stack){
m_ParticleStack = particle_stack;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void ParticleStack::AddParticleToStack(Particle& particle){
+void ParticleStack::AddParticleToStack(NPS::Particle& particle){
// If the particle is the first one to be added, then the IC are cleared
if(m_First)
@@ -103,7 +103,7 @@ void ParticleStack::AddParticleToStack(Particle& particle){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void ParticleStack::AddBeamParticleToStack(Particle& particle){
+void ParticleStack::AddBeamParticleToStack(NPS::Particle& particle){
// If the particle is the first one to be added, then the IC are cleared
if(m_First)
m_InitialConditions->Clear();
@@ -141,148 +141,148 @@ void ParticleStack::AddBeamParticleToStack(Particle& particle){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-Particle ParticleStack::SearchAndRemoveParticle(string name){
-
- unsigned int size = m_ParticleStack.size();
-
- for(unsigned int i = 0 ; i < size ; i++){
- string ParticleName = m_ParticleStack[i].GetParticleDefinition()->GetParticleName();
- if(ParticleName.compare(0, name.length(), name) == 0){
- Particle my_Particule = m_ParticleStack[i];
- m_ParticleStack.erase(m_ParticleStack.begin()+i);
- return my_Particule;
- }
+NPS::Particle ParticleStack::SearchAndRemoveParticle(string name){
+
+ unsigned int size = m_ParticleStack.size();
+
+ for(unsigned int i = 0 ; i < size ; i++){
+ string ParticleName = m_ParticleStack[i].GetParticleDefinition()->GetParticleName();
+ if(ParticleName.compare(0, name.length(), name) == 0){
+ NPS::Particle my_Particule = m_ParticleStack[i];
+ m_ParticleStack.erase(m_ParticleStack.begin()+i);
+ return my_Particule;
}
-
- return Particle();
+ }
+
+ return NPS::Particle();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
string ParticleStack::ChangeNameToG4Standard(string OriginalName){
- string NumberOfMass ;
- string Nucleid;
-
- for (unsigned int i = 0; i < OriginalName.length(); i++) {
- ostringstream character;
- character << OriginalName[i];
- if (character.str()=="0") NumberOfMass+="0";
- else if (character.str()=="1") NumberOfMass+="1";
- else if (character.str()=="2") NumberOfMass+="2";
- else if (character.str()=="3") NumberOfMass+="3";
- else if (character.str()=="4") NumberOfMass+="4";
- else if (character.str()=="5") NumberOfMass+="5";
- else if (character.str()=="6") NumberOfMass+="6";
- else if (character.str()=="7") NumberOfMass+="7";
- else if (character.str()=="8") NumberOfMass+="8";
- else if (character.str()=="9") NumberOfMass+="9";
-
- else if (character.str()=="A") Nucleid+="A";
- else if (character.str()=="B") Nucleid+="B";
- else if (character.str()=="C") Nucleid+="C";
- else if (character.str()=="D") Nucleid+="D";
- else if (character.str()=="E") Nucleid+="E";
- else if (character.str()=="F") Nucleid+="F";
- else if (character.str()=="G") Nucleid+="G";
- else if (character.str()=="H") Nucleid+="H";
- else if (character.str()=="I") Nucleid+="I";
- else if (character.str()=="J") Nucleid+="J";
- else if (character.str()=="K") Nucleid+="K";
- else if (character.str()=="L") Nucleid+="L";
- else if (character.str()=="M") Nucleid+="M";
- else if (character.str()=="N") Nucleid+="N";
- else if (character.str()=="O") Nucleid+="O";
- else if (character.str()=="P") Nucleid+="P";
- else if (character.str()=="Q") Nucleid+="Q";
- else if (character.str()=="R") Nucleid+="R";
- else if (character.str()=="S") Nucleid+="S";
- else if (character.str()=="T") Nucleid+="T";
- else if (character.str()=="U") Nucleid+="U";
- else if (character.str()=="V") Nucleid+="V";
- else if (character.str()=="W") Nucleid+="W";
- else if (character.str()=="X") Nucleid+="X";
- else if (character.str()=="Y") Nucleid+="Y";
- else if (character.str()=="Z") Nucleid+="Z";
-
- else if (character.str()=="a") Nucleid+="a";
- else if (character.str()=="b") Nucleid+="b";
- else if (character.str()=="c") Nucleid+="c";
- else if (character.str()=="d") Nucleid+="d";
- else if (character.str()=="e") Nucleid+="e";
- else if (character.str()=="f") Nucleid+="f";
- else if (character.str()=="g") Nucleid+="g";
- else if (character.str()=="h") Nucleid+="h";
- else if (character.str()=="i") Nucleid+="i";
- else if (character.str()=="j") Nucleid+="j";
- else if (character.str()=="k") Nucleid+="k";
- else if (character.str()=="l") Nucleid+="l";
- else if (character.str()=="m") Nucleid+="m";
- else if (character.str()=="n") Nucleid+="n";
- else if (character.str()=="o") Nucleid+="o";
- else if (character.str()=="p") Nucleid+="p";
- else if (character.str()=="q") Nucleid+="q";
- else if (character.str()=="r") Nucleid+="r";
- else if (character.str()=="s") Nucleid+="s";
- else if (character.str()=="t") Nucleid+="t";
- else if (character.str()=="u") Nucleid+="u";
- else if (character.str()=="v") Nucleid+="v";
- else if (character.str()=="w") Nucleid+="w";
- else if (character.str()=="x") Nucleid+="x";
- else if (character.str()=="y") Nucleid+="y";
- else if (character.str()=="z") Nucleid+="z";
- }
-
- // Special case for light particles
- string FinalName=Nucleid+NumberOfMass;
- if (FinalName=="H1") FinalName="proton";
- // else if (FinalName=="H2") FinalName="deuteron";
- else if (FinalName=="H3") FinalName="triton";
- else if (FinalName=="He4") FinalName="alpha";
- else if (FinalName=="p") FinalName="proton";
- else if (FinalName=="d") FinalName="deuteron";
- else if (FinalName=="t") FinalName="triton";
- else if (FinalName=="a") FinalName="alpha";
- else if (FinalName=="proton") FinalName="proton";
- else if (FinalName=="deuteron") FinalName="deuteron";
- else if (FinalName=="triton") FinalName="triton";
- else if (FinalName=="alpha") FinalName="alpha";
- else if (FinalName=="mu+") FinalName="mu+";
- else if (FinalName=="mu-") FinalName="mu-";
- else if (FinalName=="n") FinalName="neutron";
- else if (FinalName=="neutron") FinalName="neutron";
- return(FinalName);
+ string NumberOfMass ;
+ string Nucleid;
+
+ for (unsigned int i = 0; i < OriginalName.length(); i++) {
+ ostringstream character;
+ character << OriginalName[i];
+ if (character.str()=="0") NumberOfMass+="0";
+ else if (character.str()=="1") NumberOfMass+="1";
+ else if (character.str()=="2") NumberOfMass+="2";
+ else if (character.str()=="3") NumberOfMass+="3";
+ else if (character.str()=="4") NumberOfMass+="4";
+ else if (character.str()=="5") NumberOfMass+="5";
+ else if (character.str()=="6") NumberOfMass+="6";
+ else if (character.str()=="7") NumberOfMass+="7";
+ else if (character.str()=="8") NumberOfMass+="8";
+ else if (character.str()=="9") NumberOfMass+="9";
+
+ else if (character.str()=="A") Nucleid+="A";
+ else if (character.str()=="B") Nucleid+="B";
+ else if (character.str()=="C") Nucleid+="C";
+ else if (character.str()=="D") Nucleid+="D";
+ else if (character.str()=="E") Nucleid+="E";
+ else if (character.str()=="F") Nucleid+="F";
+ else if (character.str()=="G") Nucleid+="G";
+ else if (character.str()=="H") Nucleid+="H";
+ else if (character.str()=="I") Nucleid+="I";
+ else if (character.str()=="J") Nucleid+="J";
+ else if (character.str()=="K") Nucleid+="K";
+ else if (character.str()=="L") Nucleid+="L";
+ else if (character.str()=="M") Nucleid+="M";
+ else if (character.str()=="N") Nucleid+="N";
+ else if (character.str()=="O") Nucleid+="O";
+ else if (character.str()=="P") Nucleid+="P";
+ else if (character.str()=="Q") Nucleid+="Q";
+ else if (character.str()=="R") Nucleid+="R";
+ else if (character.str()=="S") Nucleid+="S";
+ else if (character.str()=="T") Nucleid+="T";
+ else if (character.str()=="U") Nucleid+="U";
+ else if (character.str()=="V") Nucleid+="V";
+ else if (character.str()=="W") Nucleid+="W";
+ else if (character.str()=="X") Nucleid+="X";
+ else if (character.str()=="Y") Nucleid+="Y";
+ else if (character.str()=="Z") Nucleid+="Z";
+
+ else if (character.str()=="a") Nucleid+="a";
+ else if (character.str()=="b") Nucleid+="b";
+ else if (character.str()=="c") Nucleid+="c";
+ else if (character.str()=="d") Nucleid+="d";
+ else if (character.str()=="e") Nucleid+="e";
+ else if (character.str()=="f") Nucleid+="f";
+ else if (character.str()=="g") Nucleid+="g";
+ else if (character.str()=="h") Nucleid+="h";
+ else if (character.str()=="i") Nucleid+="i";
+ else if (character.str()=="j") Nucleid+="j";
+ else if (character.str()=="k") Nucleid+="k";
+ else if (character.str()=="l") Nucleid+="l";
+ else if (character.str()=="m") Nucleid+="m";
+ else if (character.str()=="n") Nucleid+="n";
+ else if (character.str()=="o") Nucleid+="o";
+ else if (character.str()=="p") Nucleid+="p";
+ else if (character.str()=="q") Nucleid+="q";
+ else if (character.str()=="r") Nucleid+="r";
+ else if (character.str()=="s") Nucleid+="s";
+ else if (character.str()=="t") Nucleid+="t";
+ else if (character.str()=="u") Nucleid+="u";
+ else if (character.str()=="v") Nucleid+="v";
+ else if (character.str()=="w") Nucleid+="w";
+ else if (character.str()=="x") Nucleid+="x";
+ else if (character.str()=="y") Nucleid+="y";
+ else if (character.str()=="z") Nucleid+="z";
+ }
+
+ // Special case for light particles
+ string FinalName=Nucleid+NumberOfMass;
+ if (FinalName=="H1") FinalName="proton";
+ // else if (FinalName=="H2") FinalName="deuteron";
+ else if (FinalName=="H3") FinalName="triton";
+ else if (FinalName=="He4") FinalName="alpha";
+ else if (FinalName=="p") FinalName="proton";
+ else if (FinalName=="d") FinalName="deuteron";
+ else if (FinalName=="t") FinalName="triton";
+ else if (FinalName=="a") FinalName="alpha";
+ else if (FinalName=="proton") FinalName="proton";
+ else if (FinalName=="deuteron") FinalName="deuteron";
+ else if (FinalName=="triton") FinalName="triton";
+ else if (FinalName=="alpha") FinalName="alpha";
+ else if (FinalName=="mu+") FinalName="mu+";
+ else if (FinalName=="mu-") FinalName="mu-";
+ else if (FinalName=="n") FinalName="neutron";
+ else if (FinalName=="neutron") FinalName="neutron";
+ return(FinalName);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void ParticleStack::ShootAllParticle(G4Event* anEvent){
- unsigned int size = m_ParticleStack.size();
-
- for(unsigned int i = 0 ; i < size ; i++){
-
- if( m_ParticleStack[i].GetShootStatus()){
- // Write the DEDX table for charged particle and
- // all material used in the simulation
- if( anEvent->GetEventID()==0
- && G4RunManager::GetRunManager()->GetCurrentRun()->GetRunID()==0
- && m_ParticleStack[i].GetParticleDefinition()->GetPDGCharge()!=0){
- MaterialManager::getInstance()
- ->WriteDEDXTable(m_ParticleStack[i].GetParticleDefinition(),
- 1e-3*eV,m_ParticleStack[i].GetParticleKineticEnergy());
- }
-
- m_particleGun->SetParticleDefinition(m_ParticleStack[i].GetParticleDefinition());
- m_particleGun->SetParticleEnergy(m_ParticleStack[i].GetParticleKineticEnergy());
- m_particleGun->SetParticleMomentumDirection(m_ParticleStack[i].GetParticleMomentumDirection());
- m_particleGun->SetParticlePosition(m_ParticleStack[i].GetParticlePosition());
- m_particleGun->GeneratePrimaryVertex(anEvent);
- //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
- m_InitialConditions-> SetParticleName ( m_ParticleStack[i].GetParticleDefinition()->GetParticleName()) ;
- m_InitialConditions-> SetThetaCM ( m_ParticleStack[i].GetParticleThetaCM()/deg);
- m_InitialConditions-> SetKineticEnergy ( m_ParticleStack[i].GetParticleKineticEnergy());
- m_InitialConditions-> SetMomentumDirectionX ( m_ParticleStack[i].GetParticleMomentumDirection().x());
- m_InitialConditions-> SetMomentumDirectionY ( m_ParticleStack[i].GetParticleMomentumDirection().y());
- m_InitialConditions-> SetMomentumDirectionZ ( m_ParticleStack[i].GetParticleMomentumDirection().z());
- }
+ unsigned int size = m_ParticleStack.size();
+
+ for(unsigned int i = 0 ; i < size ; i++){
+
+ if( m_ParticleStack[i].GetShootStatus()){
+ // Write the DEDX table for charged particle and
+ // all material used in the simulation
+ if( anEvent->GetEventID()==0
+ && G4RunManager::GetRunManager()->GetCurrentRun()->GetRunID()==0
+ && m_ParticleStack[i].GetParticleDefinition()->GetPDGCharge()!=0){
+ MaterialManager::getInstance()
+ ->WriteDEDXTable(m_ParticleStack[i].GetParticleDefinition(),
+ 1e-3*eV,m_ParticleStack[i].GetParticleKineticEnergy());
+ }
+
+ m_particleGun->SetParticleDefinition(m_ParticleStack[i].GetParticleDefinition());
+ m_particleGun->SetParticleEnergy(m_ParticleStack[i].GetParticleKineticEnergy());
+ m_particleGun->SetParticleMomentumDirection(m_ParticleStack[i].GetParticleMomentumDirection());
+ m_particleGun->SetParticlePosition(m_ParticleStack[i].GetParticlePosition());
+ m_particleGun->GeneratePrimaryVertex(anEvent);
+ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+ m_InitialConditions-> SetParticleName ( m_ParticleStack[i].GetParticleDefinition()->GetParticleName()) ;
+ m_InitialConditions-> SetThetaCM ( m_ParticleStack[i].GetParticleThetaCM()/deg);
+ m_InitialConditions-> SetKineticEnergy ( m_ParticleStack[i].GetParticleKineticEnergy());
+ m_InitialConditions-> SetMomentumDirectionX ( m_ParticleStack[i].GetParticleMomentumDirection().x());
+ m_InitialConditions-> SetMomentumDirectionY ( m_ParticleStack[i].GetParticleMomentumDirection().y());
+ m_InitialConditions-> SetMomentumDirectionZ ( m_ParticleStack[i].GetParticleMomentumDirection().z());
}
- m_ParticleStack.clear();
- m_First=true;
+ }
+ m_ParticleStack.clear();
+ m_First=true;
}
diff --git a/NPSimulation/Core/ParticleStack.hh b/NPSimulation/Core/ParticleStack.hh
index d11d4a9bfa1248eecb8c6d3fa7657a4b0dc7b3b4..33e2cea64a26bcab8ae05c5ed2f84162163f17ab 100644
--- a/NPSimulation/Core/ParticleStack.hh
+++ b/NPSimulation/Core/ParticleStack.hh
@@ -79,22 +79,22 @@ class ParticleStack{
bool m_EventZero;
private: // Private Member
- vector<Particle> m_ParticleStack;
+ vector<NPS::Particle> m_ParticleStack;
public: // Getter and Setter
- vector<Particle> GetParticleStack();
- void SetParticleStack(vector<Particle>);
+ vector<NPS::Particle> GetParticleStack();
+ void SetParticleStack(vector<NPS::Particle>);
G4ParticleGun* GetParticleGun() {return m_particleGun;};
public: // Particle management and shooting method
// EventGenerator use this method to add particle in the stack
- void AddParticleToStack(Particle&);
+ void AddParticleToStack(NPS::Particle&);
// Add the Particle to the stack and fill Initial Condition Beam field
- void AddBeamParticleToStack(Particle&);
+ void AddBeamParticleToStack(NPS::Particle&);
// Search for a specific particle in the stack
- Particle SearchAndRemoveParticle(string);
+ NPS::Particle SearchAndRemoveParticle(string);
// Transform the particle name to G4 standard: i.e: 10He -> He10
string ChangeNameToG4Standard(string);
diff --git a/NPSimulation/EventGenerator/EventGeneratorBeam.cc b/NPSimulation/EventGenerator/EventGeneratorBeam.cc
index 4e33ce3cf296bcfc2ae0d443f2409589355d3ae9..a8684d28b867724853c9f0616783b3dd56515e42 100644
--- a/NPSimulation/EventGenerator/EventGeneratorBeam.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorBeam.cc
@@ -107,7 +107,7 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent){
///////////////////////////////////////////////////////
///// Add the Beam particle to the particle Stack /////
///////////////////////////////////////////////////////
- Particle BeamParticle( m_particle,
+ NPS::Particle BeamParticle( m_particle,
InitialBeamEnergy,
InitialBeamEnergy,
BeamDir.unit(),
diff --git a/NPSimulation/EventGenerator/EventGeneratorCosmic.cc b/NPSimulation/EventGenerator/EventGeneratorCosmic.cc
index c38c08ff0eebef286b7ba491581084755ba17335..8a55fdbb083200349f4db4b7a5c40bb2a65127c8 100644
--- a/NPSimulation/EventGenerator/EventGeneratorCosmic.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorCosmic.cc
@@ -226,7 +226,7 @@ void EventGeneratorCosmic::GenerateEvent(G4Event*){
//DistribMomZMomX->Fill(momentum_x,momentum_z);
//DistribCosmicAngle->Fill(CosmicAngle);
- Particle particle(par.m_particle, 0,particle_energy,G4ThreeVector(momentum_x, momentum_y, momentum_z),G4ThreeVector(x0, par.m_y0, z0));
+ NPS::Particle particle(par.m_particle, 0,particle_energy,G4ThreeVector(momentum_x, momentum_y, momentum_z),G4ThreeVector(x0, par.m_y0, z0));
m_ParticleStack->AddParticleToStack(particle);
diff --git a/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc b/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc
index 21b514ed86e842f9c2e4b2edecd49f8dcb9498a1..70b703588901d5fe9975d8acc9e64c280ea962ea 100644
--- a/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc
@@ -48,7 +48,9 @@ EventGeneratorGEFReader::EventGeneratorGEFReader(){
m_ParticleStack = ParticleStack::getInstance();
event_ID=0;
+ m_isTwoBody = false;
HasInputDataFile = false;
+ m_FissionConditions = new TFissionConditions();
}
@@ -71,8 +73,19 @@ EventGeneratorGEFReader::SourceParameters::SourceParameters(){
m_direction = 'z' ;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void EventGeneratorGEFReader::AttachFissionConditions(){
+ if(RootOutput::getInstance()->GetTree()->FindBranch("FissionConditions"))
+ RootOutput::getInstance()->GetTree()->SetBranchAddress("FissionConditions", &m_FissionConditions);
+}
+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorGEFReader::ReadConfiguration(NPL::InputParser parser){
+
+ AttachFissionConditions();
+ if(!RootOutput::getInstance()->GetTree()->FindBranch("FissionConditions")){
+ RootOutput::getInstance()->GetTree()->Branch("FissionConditions", "TFissionConditions", &m_FissionConditions);
+ }
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("GEFReader");
m_Parameters.reserve(blocks.size());
if(NPOptionManager::getInstance()->GetVerboseLevel())
@@ -85,45 +98,46 @@ void EventGeneratorGEFReader::ReadConfiguration(NPL::InputParser parser){
const vector<SourceParameters>::reverse_iterator it = m_Parameters.rbegin();
if(blocks[i]->HasToken("GEFversion"))
- {
- it->m_GEFversion=blocks[i]->GetDouble("GEFversion","");
- if(it->m_GEFversion>=2023.32 && it->m_GEFversion<=2024.33) version_shift=0;
- else if (it->m_GEFversion>=2023.11 && it->m_GEFversion<=2023.12) version_shift=1;
- else {
- cout << "!!!!!!!!!!!!!!!!! ERROR: GEF version not verified !!!!!!!!!!!!!!!!!!" << endl;
- cout << "Check the column positions of your GEF file and" << endl;
- cout << "add the version to the correct shift in EventGeneratorGEFReader.cc" << endl;
- exit(1);
- }
- }
+ {
+ it->m_GEFversion=blocks[i]->GetDouble("GEFversion","");
+ if(it->m_GEFversion>=2023.32 && it->m_GEFversion<=2024.33) version_shift=0;
+ else if (it->m_GEFversion>=2023.11 && it->m_GEFversion<=2023.12) version_shift=1;
+ else {
+ cout << "!!!!!!!!!!!!!!!!! ERROR: GEF version not verified !!!!!!!!!!!!!!!!!!" << endl;
+ cout << "Check the column positions of your GEF file and" << endl;
+ cout << "add the version to the correct shift in EventGeneratorGEFReader.cc" << endl;
+ exit(1);
+ }
+ }
if(blocks[i]->HasToken("InputDataFile")){
- vector<string> file = blocks[i]->GetVectorString("InputDataFile");
- fInputDataFile.open(file.at(0).c_str());
- HasInputDataFile = true;
-
- string line;
- string comment="*";
- double data=0;
-
- while(comment=="*" || data==0)
- {
- getline(fInputDataFile,line);
- istringstream iss(line);
- //cout << "line : " << line << endl;
-
- iss >> data;
- comment=to_string(data);
- if(comment!="*" && data>0)
- {// Storing the first line to be read by GenerateEvent
- LastLine.push_back(data);
- while(iss >> data) LastLine.push_back(data);
- break;
- }
- }
+ vector<string> file = blocks[i]->GetVectorString("InputDataFile");
+ fInputDataFile.open(file.at(0).c_str());
+ HasInputDataFile = true;
+
+ string line;
+ string comment="*";
+ double data=0;
+
+ while(comment=="*" || data==0)
+ {
+ getline(fInputDataFile,line);
+ istringstream iss(line);
+ //cout << "line : " << line << endl;
+
+ iss >> data;
+ comment=to_string(data);
+ if(comment!="*" && data>0)
+ {// Storing the first line to be read by GenerateEvent
+ LastLine.push_back(data);
+ while(iss >> data) LastLine.push_back(data);
+ break;
+ }
+ }
}
-
- if(blocks[i]->HasToken("Direction"))
+
+ if(blocks[i]->HasToken("Direction")){
it->m_direction=blocks[i]->GetString("Direction");
+ }
it->m_x0 = blocks[i]->GetDouble("x0","mm");
it->m_y0 = blocks[i]->GetDouble("y0","mm");
it->m_z0 = blocks[i]->GetDouble("z0","mm");
@@ -134,28 +148,33 @@ void EventGeneratorGEFReader::ReadConfiguration(NPL::InputParser parser){
else if(particleName[j]=="neutron") {it->m_particleName.push_back("neutron") ;}
else it->m_particleName.push_back(particleName[j]);
}
- if(blocks[i]->HasToken("KineticEnergy_FS"))
- it->m_Boost=blocks[i]->GetDouble("KineticEnergy_FS","MeV");
+ if(blocks[i]->HasToken("TwoBodyReaction")){
+ string my_reaction = blocks[i]->GetString("TwoBodyReaction");
+ m_TwoBodyReaction = new NPL::Reaction(my_reaction);
+ m_isTwoBody = true;
+ }
+
+ //if(blocks[i]->HasToken("KineticEnergy_FS"))
+ //it->m_Boost=blocks[i]->GetDouble("KineticEnergy_FS","MeV");
if(blocks[i]->HasToken("FissioningSystem"))
- {
- it->m_FissioningSystem=blocks[i]->GetString("FissioningSystem");
- FissioningSystem=new NPL::Particle(it->m_FissioningSystem);
- FissioningSystem->SetKineticEnergy(it->m_Boost,0,0);
- }
+ {
+ it->m_FissioningSystemName=blocks[i]->GetString("FissioningSystem");
+ m_FissioningSystem=new NPL::Particle(it->m_FissioningSystemName);
+ }
if(blocks[i]->HasToken("SigmaX"))
it->m_SigmaX=blocks[i]->GetDouble("SigmaX","mm");
if(blocks[i]->HasToken("SigmaY"))
it->m_SigmaY=blocks[i]->GetDouble("SigmaY","mm");
if(blocks[i]->HasToken("SigmaZ"))
it->m_SigmaZ=blocks[i]->GetDouble("SigmaZ","mm");
-
+
}
else{
cout << "ERROR: check your input file formatting \033[0m" << endl;
exit(1);
}
}
-
+
for(auto& par : m_Parameters) {
for(auto partName: par.m_particleName){
AllowedParticles.push_back(partName);
@@ -168,11 +187,25 @@ void EventGeneratorGEFReader::ReadConfiguration(NPL::InputParser parser){
}
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void EventGeneratorGEFReader::GetBoostFromTwoBodyReaction(double Ex){
+ double Theta3, E3;
+ double Theta4, E4;
+ double ThetaCM = RandFlat::shoot() * 2 * pi;
+ m_TwoBodyReaction->SetExcitation4(Ex);
+ m_TwoBodyReaction->SetThetaCM(ThetaCM);
+ m_TwoBodyReaction->KineRelativistic(Theta3,E3,Theta4,E4);
+
+ double Phi4 = RandFlat::shoot() * 2 * pi;
+ m_FissioningSystem->SetKineticEnergy(E4,Theta4,Phi4);
+
+}
+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorGEFReader::GenerateEvent(G4Event*){
for(auto& par : m_Parameters) {
int fileParticleMultiplicity=-1;
-
+
vector<double> ELab;
vector<double> CosThetaLab;
vector<double> PhiLab;
@@ -187,7 +220,6 @@ void EventGeneratorGEFReader::GenerateEvent(G4Event*){
vector<TVector3> FragmentBoostforN;
double CN_ExcitationEnergy;
- TVector3 FissioningSystemBoost=FissioningSystem->GetEnergyImpulsion().BoostVector();
TVector3 LightFragmentDirection(0,0,1.);
TVector3 HeavyFragmentDirection(0,0,1.);
TLorentzVector ImpulsionFrag[2];
@@ -196,269 +228,292 @@ void EventGeneratorGEFReader::GenerateEvent(G4Event*){
NPL::Particle *Proton=new NPL::Particle("1H");
NPL::Particle *Gamma=new NPL::Particle("gamma");
NPL::Particle *Fragment;
-
+ TVector3 FissioningSystemBoost;
int it=0;
if(HasInputDataFile && !fInputDataFile.eof())
+ {
+ fileParticleMultiplicity=0;
+ // for(int GEFline=0;GEFline<4;GEFline++)
+ while(it==0 || LastLine.size()==1 || (LastLine.size()>0 && LastLine.at(1)<50))
{
- fileParticleMultiplicity=0;
- // for(int GEFline=0;GEFline<4;GEFline++)
- while(it==0 || LastLine.size()==1 || (LastLine.size()>0 && LastLine.at(1)<50))
- {
- it++;
- string line;
- getline(fInputDataFile,line);
- //cout << line << endl;
-
- istringstream iss(line);
- double data;
- vector<double> DataLine;
-
- // cout << "it=" << it << " line=" << line << endl;
-
- DataLine = LastLine;
- LastLine.clear();
-
- while(iss >> data) LastLine.push_back(data);
-
- //cout << "line starts with " << DataLine.at(0) << " length : " << DataLine.size() << endl;
- //for (int d=0; d<DataLine.size(); d++){
- // cout << d << " " << DataLine.at(d) << endl;
- //}
-
- if(DataLine.size()>=26 && DataLine.at(1)>50)
- {
- // Fragment emission
- for(int frag=0;frag<2;frag++)
- {
- int Zf = DataLine.at(3 + frag);
- int Af = DataLine.at(7 + frag); // Post-scission
- double ELabf = DataLine.at(13 + 3*frag);
- double cos_thetaf = DataLine.at(14 + 3*frag);
- double Phif = DataLine.at(15 + 3*frag) *M_PI/180;
-
- Fragment=new NPL::Particle(Zf,Af);
- Fragment->EnergyToEnergyImpulsion(ELabf,acos(cos_thetaf),Phif);
-
- ImpulsionFrag[frag] = Fragment->GetEnergyImpulsion();
- ImpulsionFrag[frag].Boost(FissioningSystemBoost);
-
- if(frag==0) LightFragmentDirection.SetMagThetaPhi(1.,acos(cos_thetaf),Phif);
- if(frag==1) HeavyFragmentDirection.SetMagThetaPhi(1.,acos(cos_thetaf),Phif);
- TVector3 Boost=Fragment->GetEnergyImpulsion().BoostVector();
-
- ELabf = ImpulsionFrag[frag].E() - Fragment->Mass();
- cos_thetaf = ImpulsionFrag[frag].CosTheta();
-
- ZFrag .push_back( Zf );
- AFrag .push_back( Af );
- ELabFrag .push_back( ELabf );
- cos_thetaFrag .push_back( cos_thetaf);
- PhiFrag .push_back( Phif);
- NneutronsEmission.push_back( DataLine.at(21 + frag));
- FragmentBoostforN.push_back(Boost);
- //delete Fragment;
- }
- CN_ExcitationEnergy = DataLine.at(25);
-
- if(DataLine.size()>26)
- {
- // Particle emission by the CN before fission
- string ParticleList = to_string(DataLine.at(26));
- int idx=0;
- for(char type: ParticleList)
- if((type=='1' || type=='3' || type=='5')
- && (AllowedParticles.size()==0 || std::find(AllowedParticles.begin(), AllowedParticles.end(), "neutron") != AllowedParticles.end()))
- {// Pre-fission neutron treatment:
-
- double ELabn = DataLine.at(27 + idx*3);
- double cos_thetan = DataLine.at(28 + idx*3);
- double Phin = RandFlat::shoot() * 2 * pi;
- Neutron->SetKineticEnergy(ELabn,acos(cos_thetan),Phin);
-
- TLorentzVector ImpulsionNeutron = Neutron->GetEnergyImpulsion();
- ImpulsionNeutron.Boost(FissioningSystemBoost);
- ELabn = ImpulsionNeutron.E() - Neutron->Mass();
- cos_thetan = ImpulsionNeutron.CosTheta();
-
- ELab .push_back(ELabn);
- CosThetaLab.push_back(cos_thetan);
- PhiLab .push_back(Phin);
- vParticle .push_back(G4ParticleTable::GetParticleTable()->FindParticle("neutron"));
-
- fileParticleMultiplicity++;
- idx++;
- }
- else if((type=='2' || type=='4')
- && (AllowedParticles.size()==0 || std::find(AllowedParticles.begin(), AllowedParticles.end(), "proton") != AllowedParticles.end()))
- {// Pre-fission proton treatment:
-
-
- double ELabp = DataLine.at(27 + idx*3);
- double cos_thetap = DataLine.at(28 + idx*3);
- double Phip = RandFlat::shoot() * 2 * pi;
- Proton->SetKineticEnergy(ELabp,acos(cos_thetap),Phip);
-
- TLorentzVector ImpulsionProton = Proton->GetEnergyImpulsion();
- ImpulsionProton.Boost(FissioningSystemBoost);
- ELabp = ImpulsionProton.E() - Proton->Mass();
- cos_thetap = ImpulsionProton.CosTheta();
-
- ELab .push_back(ELabp);
- CosThetaLab.push_back(cos_thetap);
- PhiLab .push_back(RandFlat::shoot() * 2 * pi);
- vParticle .push_back(G4ParticleTable::GetParticleTable()->FindParticle("proton"));
-
- fileParticleMultiplicity++;
- idx++;
- }
- else if(type=='0') break;
- }
- }
- else if(DataLine.size()>0 && DataLine.at(0)==0
- && (AllowedParticles.size()==0 || std::find(AllowedParticles.begin(), AllowedParticles.end(), "neutron") != AllowedParticles.end()))
- for(int it=0;it<(DataLine.size()-1)/3;it++)
- {// Prompt fission neutron treatment (just get the angles)
-
- double ELabn = DataLine.at(1+3*it);
- double cos_thetan = DataLine.at(2+3*it);
- double Phin = DataLine.at(3+3*it) *M_PI/180.;
-
- //std::cout << " Elab 0 " << ELabn << std::endl;
-
- ELab .push_back(ELabn);
- CosThetaLab.push_back(cos_thetan);
- PhiLab .push_back(Phin);
- vParticle .push_back(G4ParticleTable::GetParticleTable()->FindParticle("neutron"));
-
- fileParticleMultiplicity++;
- }
-
- else if(DataLine.size()>0 && DataLine.at(0)==1
- && (AllowedParticles.size()==0 || std::find(AllowedParticles.begin(), AllowedParticles.end(), "neutron") != AllowedParticles.end()))
- for(int it=0;it<(DataLine.size()-1)/2;it++)
- {// LIGHT fragment prompt fission neutron treatment (apply boost by linking neutron to fragment emitting it)
-
- double ELabn_for_ID = DataLine.at(2+2*it);
-
- //std::cout << " Elab for ID " << ELabn_for_ID << std::endl;
-
- // Find the index of ELabn_for_ID in Elab
- auto iterator = std::find(ELab.begin(), ELab.end(), ELabn_for_ID);
-
- int index;
- if (iterator != ELab.end()) {
- // Calculate the index
- index = std::distance(ELab.begin(), iterator);
- //std::cout << "Index of " << ELabn_for_ID << " is " << index << std::endl;
- } else {
- //std::cout << " " << ELabn_for_ID << " not found in ELab." << std::endl;
- }
-
- double ELabn = ELab.at(index);
- double cos_thetan = CosThetaLab.at(index);
- double Phin = PhiLab.at(index);
-
- TVector3 NeutronLabAngle(0,0,1.);
- NeutronLabAngle.SetMagThetaPhi(1.,acos(cos_thetan),Phin);
- NeutronLabAngle.RotateUz(LightFragmentDirection);
- Neutron->SetKineticEnergy(ELabn,NeutronLabAngle.Theta(),NeutronLabAngle.Phi());
-
- TLorentzVector ImpulsionNeutron = Neutron->GetEnergyImpulsion();
- ImpulsionNeutron.Boost(FragmentBoostforN.at(0));
- ELabn = ImpulsionNeutron.E()-Neutron->Mass();
- cos_thetan = ImpulsionNeutron.CosTheta();
- Phin = NeutronLabAngle.Phi();
-
- ELab.at(index) = ELabn;
- CosThetaLab.at(index) = cos_thetan;
- PhiLab.at(index) = Phin;
- }
-
- else if(DataLine.size()>0 && DataLine.at(0)==2
- && (AllowedParticles.size()==0 || std::find(AllowedParticles.begin(), AllowedParticles.end(), "neutron") != AllowedParticles.end()))
- for(int it=0;it<(DataLine.size()-1)/2;it++)
- {// HEAVY fragment prompt fission neutron treatment (apply boost by linking neutron to fragment emitting it)
-
- double ELabn_for_ID = DataLine.at(2+2*it);
-
- // Find the index of ELabn_for_ID in Elab
- auto iterator = std::find(ELab.begin(), ELab.end(), ELabn_for_ID);
-
- int index;
- if (iterator != ELab.end()) {
- // Calculate the index
- index = std::distance(ELab.begin(), iterator);
- //std::cout << "Index of " << ELabn_for_ID << " is " << index << std::endl;
- } else {
- //std::cout << " " << ELabn_for_ID << " not found in ELab." << std::endl;
- }
-
- double ELabn = ELab.at(index);
- double cos_thetan = CosThetaLab.at(index);
- double Phin = PhiLab.at(index);
-
- TVector3 NeutronLabAngle(0,0,1.);
- NeutronLabAngle.SetMagThetaPhi(1.,acos(cos_thetan),Phin);
- NeutronLabAngle.RotateUz(HeavyFragmentDirection);
- Neutron->SetKineticEnergy(ELabn,NeutronLabAngle.Theta(),NeutronLabAngle.Phi());
-
- TLorentzVector ImpulsionNeutron = Neutron->GetEnergyImpulsion();
- ImpulsionNeutron.Boost(FragmentBoostforN.at(1));
- ELabn = ImpulsionNeutron.E()-Neutron->Mass();
- cos_thetan = ImpulsionNeutron.CosTheta();
- Phin = NeutronLabAngle.Phi();
-
- ELab.at(index) = ELabn;
- CosThetaLab.at(index) = cos_thetan;
- PhiLab.at(index) = Phin;
- }
-
- else if(DataLine.size()>0 && (DataLine.at(0)>=3 && DataLine.at(0)<=8)
- && (AllowedParticles.size()==0 || std::find(AllowedParticles.begin(), AllowedParticles.end(), "gamma") != AllowedParticles.end()))
- for(int it=0;it<(DataLine.size()-1);it++)
- {// Prompt fission gamma treatment
-
- if((DataLine.at(0)==5 || DataLine.at(0)==8) && (to_string(DataLine.at(it+1))=="GS" || to_string(DataLine.at(it+1)).at(0)=='M'))
- {
- cout << "DataLine.at(it+1)=" << DataLine.at(it+1) << endl;
- break;
- }
-
- double ELabg = DataLine.at(it+1);
- double cos_thetag = RandFlat::shoot(-1.,1.);
- double Phig = RandFlat::shoot() * 2 * pi;
-
- Gamma->SetKineticEnergy(ELabg,acos(cos_thetag),Phig);
- TLorentzVector ImpulsionGamma = Gamma->GetEnergyImpulsion();
- if(DataLine.at(0)<=5)
- ImpulsionGamma.Boost(ImpulsionFrag[0].BoostVector());
- else if(DataLine.at(0)<=8)
- ImpulsionGamma.Boost(ImpulsionFrag[1].BoostVector());
-
- ImpulsionGamma.Boost(FissioningSystemBoost);
- ELabg = ImpulsionGamma.E();
- cos_thetag = ImpulsionGamma.CosTheta();
-
- ELab .push_back( ELabg);
- CosThetaLab.push_back( cos_thetag); // Need to add the fragment boost to the gamma emission.
- PhiLab .push_back( Phig);
- vParticle .push_back(G4ParticleTable::GetParticleTable()->FindParticle("gamma"));
-
- fileParticleMultiplicity++;
- }
- }
- if(fInputDataFile.eof())
- {
- cout << "Problem with Input data file ? no more data to put in" << endl;
- return;
- }
+ it++;
+ string line;
+ getline(fInputDataFile,line);
+ //cout << line << endl;
+
+ istringstream iss(line);
+ double data;
+ vector<double> DataLine;
+
+ // cout << "it=" << it << " line=" << line << endl;
+
+ DataLine = LastLine;
+ LastLine.clear();
+
+ while(iss >> data) LastLine.push_back(data);
+
+ //cout << "line starts with " << DataLine.at(0) << " length : " << DataLine.size() << endl;
+ //for (int d=0; d<DataLine.size(); d++){
+ // cout << d << " " << DataLine.at(d) << endl;
+ //}
+
+ if(DataLine.size()>=26 && DataLine.at(1)>50)
+ {
+ double Ex = DataLine.at(25); // Excitation energy at fission
+ if(m_isTwoBody){
+ GetBoostFromTwoBodyReaction(Ex);
+ }
+ FissioningSystemBoost = m_FissioningSystem->GetEnergyImpulsion().BoostVector();
+
+ int Z_CN = DataLine.at(1);
+ int A_CN = DataLine.at(2);
+
+ m_FissionConditions->SetZ_CN(Z_CN);
+ m_FissionConditions->SetA_CN(A_CN);
+ m_FissionConditions->SetEx_CN(Ex);
+ m_FissionConditions->SetELab_CN(m_FissioningSystem->GetEnergy());
+ m_FissionConditions->SetThetaLab_CN(m_FissioningSystem->GetEnergyImpulsion().Theta());
+
+ // Fragment emission
+ for(int frag=0;frag<2;frag++)
+ {
+ int Zf = DataLine.at(3 + frag);
+ int Af = DataLine.at(7 + frag); // Post-scission
+ double ELabf = DataLine.at(13 + 3*frag);
+ double cos_thetaf = DataLine.at(14 + 3*frag);
+ double Phif = DataLine.at(15 + 3*frag) *M_PI/180;
+ Fragment=new NPL::Particle(Zf,Af);
+ Fragment->SetKineticEnergy(ELabf);
+ Fragment->EnergyToEnergyImpulsion(ELabf,acos(cos_thetaf),Phif);
+
+ ImpulsionFrag[frag] = Fragment->GetEnergyImpulsion();
+ ImpulsionFrag[frag].Boost(FissioningSystemBoost);
+
+ if(frag==0) LightFragmentDirection.SetMagThetaPhi(1.,acos(cos_thetaf),Phif);
+ if(frag==1) HeavyFragmentDirection.SetMagThetaPhi(1.,acos(cos_thetaf),Phif);
+ TVector3 Boost=Fragment->GetEnergyImpulsion().BoostVector();
+
+ ELabf = ImpulsionFrag[frag].E() - Fragment->Mass();
+ cos_thetaf = ImpulsionFrag[frag].CosTheta();
+
+ ZFrag .push_back( Zf );
+ AFrag .push_back( Af );
+ ELabFrag .push_back( ELabf );
+ cos_thetaFrag .push_back( cos_thetaf);
+ PhiFrag .push_back( Phif);
+ NneutronsEmission.push_back( DataLine.at(21 + frag));
+ FragmentBoostforN.push_back(Boost);
+ //delete Fragment;
+
+ // Filling fission conditions
+ m_FissionConditions->SetFragmentZ(Zf);
+ m_FissionConditions->SetFragmentA(Af);
+ m_FissionConditions->SetFragmentKineticEnergy(ELabf);
+ m_FissionConditions->SetFragmentBrho(Fragment->GetBrho());
+ m_FissionConditions->SetFragmentTheta(acos(cos_thetaf));
+ m_FissionConditions->SetFragmentPhi(Phif);
+ }
+ CN_ExcitationEnergy = DataLine.at(25);
+
+ if(DataLine.size()>26)
+ {
+ // Particle emission by the CN before fission
+ string ParticleList = to_string(DataLine.at(26));
+ int idx=0;
+ for(char type: ParticleList)
+ if((type=='1' || type=='3' || type=='5')
+ && (AllowedParticles.size()==0 || std::find(AllowedParticles.begin(), AllowedParticles.end(), "neutron") != AllowedParticles.end()))
+ {// Pre-fission neutron treatment:
+
+ double ELabn = DataLine.at(27 + idx*3);
+ double cos_thetan = DataLine.at(28 + idx*3);
+ double Phin = RandFlat::shoot() * 2 * pi;
+ Neutron->SetKineticEnergy(ELabn,acos(cos_thetan),Phin);
+
+ TLorentzVector ImpulsionNeutron = Neutron->GetEnergyImpulsion();
+ ImpulsionNeutron.Boost(FissioningSystemBoost);
+ ELabn = ImpulsionNeutron.E() - Neutron->Mass();
+ cos_thetan = ImpulsionNeutron.CosTheta();
+
+ ELab .push_back(ELabn);
+ CosThetaLab.push_back(cos_thetan);
+ PhiLab .push_back(Phin);
+ vParticle .push_back(G4ParticleTable::GetParticleTable()->FindParticle("neutron"));
+
+ fileParticleMultiplicity++;
+ idx++;
+ }
+ else if((type=='2' || type=='4')
+ && (AllowedParticles.size()==0 || std::find(AllowedParticles.begin(), AllowedParticles.end(), "proton") != AllowedParticles.end()))
+ {// Pre-fission proton treatment:
+
+
+ double ELabp = DataLine.at(27 + idx*3);
+ double cos_thetap = DataLine.at(28 + idx*3);
+ double Phip = RandFlat::shoot() * 2 * pi;
+ Proton->SetKineticEnergy(ELabp,acos(cos_thetap),Phip);
+
+ TLorentzVector ImpulsionProton = Proton->GetEnergyImpulsion();
+ ImpulsionProton.Boost(FissioningSystemBoost);
+ ELabp = ImpulsionProton.E() - Proton->Mass();
+ cos_thetap = ImpulsionProton.CosTheta();
+
+ ELab .push_back(ELabp);
+ CosThetaLab.push_back(cos_thetap);
+ PhiLab .push_back(RandFlat::shoot() * 2 * pi);
+ vParticle .push_back(G4ParticleTable::GetParticleTable()->FindParticle("proton"));
+
+ fileParticleMultiplicity++;
+ idx++;
+ }
+ else if(type=='0') break;
+ }
+ }
+ else if(DataLine.size()>0 && DataLine.at(0)==0
+ && (AllowedParticles.size()==0 || std::find(AllowedParticles.begin(), AllowedParticles.end(), "neutron") != AllowedParticles.end()))
+ for(int it=0;it<(DataLine.size()-1)/3;it++)
+ {// Prompt fission neutron treatment (just get the angles)
+
+ double ELabn = DataLine.at(1+3*it);
+ double cos_thetan = DataLine.at(2+3*it);
+ double Phin = DataLine.at(3+3*it) *M_PI/180.;
+
+ //std::cout << " Elab 0 " << ELabn << std::endl;
+
+ ELab .push_back(ELabn);
+ CosThetaLab.push_back(cos_thetan);
+ PhiLab .push_back(Phin);
+ vParticle .push_back(G4ParticleTable::GetParticleTable()->FindParticle("neutron"));
+
+ fileParticleMultiplicity++;
+ }
+
+ else if(DataLine.size()>0 && DataLine.at(0)==1
+ && (AllowedParticles.size()==0 || std::find(AllowedParticles.begin(), AllowedParticles.end(), "neutron") != AllowedParticles.end()))
+ for(int it=0;it<(DataLine.size()-1)/2;it++)
+ {// LIGHT fragment prompt fission neutron treatment (apply boost by linking neutron to fragment emitting it)
+
+ double ELabn_for_ID = DataLine.at(2+2*it);
+
+ //std::cout << " Elab for ID " << ELabn_for_ID << std::endl;
+
+ // Find the index of ELabn_for_ID in Elab
+ auto iterator = std::find(ELab.begin(), ELab.end(), ELabn_for_ID);
+
+ int index;
+ if (iterator != ELab.end()) {
+ // Calculate the index
+ index = std::distance(ELab.begin(), iterator);
+ //std::cout << "Index of " << ELabn_for_ID << " is " << index << std::endl;
+ } else {
+ //std::cout << " " << ELabn_for_ID << " not found in ELab." << std::endl;
+ }
+
+ double ELabn = ELab.at(index);
+ double cos_thetan = CosThetaLab.at(index);
+ double Phin = PhiLab.at(index);
+
+ TVector3 NeutronLabAngle(0,0,1.);
+ NeutronLabAngle.SetMagThetaPhi(1.,acos(cos_thetan),Phin);
+ NeutronLabAngle.RotateUz(LightFragmentDirection);
+ Neutron->SetKineticEnergy(ELabn,NeutronLabAngle.Theta(),NeutronLabAngle.Phi());
+
+ TLorentzVector ImpulsionNeutron = Neutron->GetEnergyImpulsion();
+ ImpulsionNeutron.Boost(FragmentBoostforN.at(0));
+ ELabn = ImpulsionNeutron.E()-Neutron->Mass();
+ cos_thetan = ImpulsionNeutron.CosTheta();
+ Phin = NeutronLabAngle.Phi();
+
+ ELab.at(index) = ELabn;
+ CosThetaLab.at(index) = cos_thetan;
+ PhiLab.at(index) = Phin;
+ }
+
+ else if(DataLine.size()>0 && DataLine.at(0)==2
+ && (AllowedParticles.size()==0 || std::find(AllowedParticles.begin(), AllowedParticles.end(), "neutron") != AllowedParticles.end()))
+ for(int it=0;it<(DataLine.size()-1)/2;it++)
+ {// HEAVY fragment prompt fission neutron treatment (apply boost by linking neutron to fragment emitting it)
+
+ double ELabn_for_ID = DataLine.at(2+2*it);
+
+ // Find the index of ELabn_for_ID in Elab
+ auto iterator = std::find(ELab.begin(), ELab.end(), ELabn_for_ID);
+
+ int index;
+ if (iterator != ELab.end()) {
+ // Calculate the index
+ index = std::distance(ELab.begin(), iterator);
+ //std::cout << "Index of " << ELabn_for_ID << " is " << index << std::endl;
+ } else {
+ //std::cout << " " << ELabn_for_ID << " not found in ELab." << std::endl;
+ }
+
+ double ELabn = ELab.at(index);
+ double cos_thetan = CosThetaLab.at(index);
+ double Phin = PhiLab.at(index);
+
+ TVector3 NeutronLabAngle(0,0,1.);
+ NeutronLabAngle.SetMagThetaPhi(1.,acos(cos_thetan),Phin);
+ NeutronLabAngle.RotateUz(HeavyFragmentDirection);
+ Neutron->SetKineticEnergy(ELabn,NeutronLabAngle.Theta(),NeutronLabAngle.Phi());
+
+ TLorentzVector ImpulsionNeutron = Neutron->GetEnergyImpulsion();
+ ImpulsionNeutron.Boost(FragmentBoostforN.at(1));
+ ELabn = ImpulsionNeutron.E()-Neutron->Mass();
+ cos_thetan = ImpulsionNeutron.CosTheta();
+ Phin = NeutronLabAngle.Phi();
+
+ ELab.at(index) = ELabn;
+ CosThetaLab.at(index) = cos_thetan;
+ PhiLab.at(index) = Phin;
+ }
+
+ else if(DataLine.size()>0 && (DataLine.at(0)>=3 && DataLine.at(0)<=8)
+ && (AllowedParticles.size()==0 || std::find(AllowedParticles.begin(), AllowedParticles.end(), "gamma") != AllowedParticles.end()))
+ for(int it=0;it<(DataLine.size()-1);it++)
+ {// Prompt fission gamma treatment
+
+ if((DataLine.at(0)==5 || DataLine.at(0)==8) && (to_string(DataLine.at(it+1))=="GS" || to_string(DataLine.at(it+1)).at(0)=='M'))
+ {
+ cout << "DataLine.at(it+1)=" << DataLine.at(it+1) << endl;
+ break;
+ }
+
+ double ELabg = DataLine.at(it+1);
+ double cos_thetag = RandFlat::shoot(-1.,1.);
+ double Phig = RandFlat::shoot() * 2 * pi;
+
+ Gamma->SetKineticEnergy(ELabg,acos(cos_thetag),Phig);
+ TLorentzVector ImpulsionGamma = Gamma->GetEnergyImpulsion();
+ if(DataLine.at(0)<=5)
+ ImpulsionGamma.Boost(ImpulsionFrag[0].BoostVector());
+ else if(DataLine.at(0)<=8)
+ ImpulsionGamma.Boost(ImpulsionFrag[1].BoostVector());
+
+ ImpulsionGamma.Boost(FissioningSystemBoost);
+ ELabg = ImpulsionGamma.E();
+ cos_thetag = ImpulsionGamma.CosTheta();
+
+ ELab .push_back( ELabg);
+ CosThetaLab.push_back( cos_thetag); // Need to add the fragment boost to the gamma emission.
+ PhiLab .push_back( Phig);
+ vParticle .push_back(G4ParticleTable::GetParticleTable()->FindParticle("gamma"));
+
+ fileParticleMultiplicity++;
+ }
}
- else
+ if(fInputDataFile.eof())
{
- cout << "Problem with Input data file ? no data to put in" << endl;
- return;
+ cout << "Problem with Input data file ? no more data to put in" << endl;
+ return;
}
+ }
+ else
+ {
+ cout << "Problem with Input data file ? no data to put in" << endl;
+ return;
+ }
G4double x0, y0, z0;
TVector3 Momentum;
@@ -468,31 +523,31 @@ void EventGeneratorGEFReader::GenerateEvent(G4Event*){
z0 = RandGauss::shoot(par.m_z0,par.m_SigmaZ);
for(int i_m=0;i_m<fileParticleMultiplicity;i_m++)
- {
- G4double theta = acos(CosThetaLab.at(i_m)) ;
- G4double phi = PhiLab.at(i_m) ;
- G4double particle_energy = ELab.at(i_m) ;
- G4ParticleDefinition* GeneratedParticle = vParticle.at(i_m) ;
-
- Momentum = ShootParticle(theta,phi,par.m_direction);
-
- Particle particle(GeneratedParticle, theta,particle_energy,G4ThreeVector(Momentum.x(), Momentum.y(), Momentum.z()),G4ThreeVector(x0, y0, z0));
- if(particle_energy<=20.0) // NeutronHP crashes for neutrons of higher energies. Need to use a different physics list.
- m_ParticleStack->AddParticleToStack(particle);
- }
+ {
+ G4double theta = acos(CosThetaLab.at(i_m)) ;
+ G4double phi = PhiLab.at(i_m) ;
+ G4double particle_energy = ELab.at(i_m) ;
+ G4ParticleDefinition* GeneratedParticle = vParticle.at(i_m) ;
+
+ Momentum = ShootParticle(theta,phi,par.m_direction);
+
+ NPS::Particle particle(GeneratedParticle, theta,particle_energy,G4ThreeVector(Momentum.x(), Momentum.y(), Momentum.z()),G4ThreeVector(x0, y0, z0));
+ if(particle_energy<=20.0) // NeutronHP crashes for neutrons of higher energies. Need to use a different physics list.
+ m_ParticleStack->AddParticleToStack(particle);
+ }
if(AllowedParticles.size()==0 || std::find(AllowedParticles.begin(), AllowedParticles.end(), "fragments") != AllowedParticles.end())
for(int i_m=0;i_m<2;i_m++)
- {
- G4double theta = acos(cos_thetaFrag.at(i_m)) ;
- G4double phi = PhiFrag.at(i_m) ;
- G4double particle_energy = ELabFrag.at(i_m) ;
- G4ParticleDefinition* GeneratedParticle = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(ZFrag.at(i_m),AFrag.at(i_m),0.);
-
- Momentum = ShootParticle(theta,phi,par.m_direction);
- //cout << "particle_energy=" << particle_energy << endl;
- Particle particle(GeneratedParticle, theta,particle_energy,G4ThreeVector(Momentum.x(), Momentum.y(), Momentum.z()),G4ThreeVector(x0, y0, z0));
- m_ParticleStack->AddParticleToStack(particle);
- }
+ {
+ G4double theta = acos(cos_thetaFrag.at(i_m)) ;
+ G4double phi = PhiFrag.at(i_m) ;
+ G4double particle_energy = ELabFrag.at(i_m) ;
+ G4ParticleDefinition* GeneratedParticle = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(ZFrag.at(i_m),AFrag.at(i_m),0.);
+
+ Momentum = ShootParticle(theta,phi,par.m_direction);
+ //cout << "particle_energy=" << particle_energy << endl;
+ NPS::Particle particle(GeneratedParticle, theta,particle_energy,G4ThreeVector(Momentum.x(), Momentum.y(), Momentum.z()),G4ThreeVector(x0, y0, z0));
+ m_ParticleStack->AddParticleToStack(particle);
+ }
}
}
@@ -502,31 +557,31 @@ void EventGeneratorGEFReader::InitializeRootOutput(){
}
TVector3 EventGeneratorGEFReader::ShootParticle(double theta,double phi,TString direction){
-
+
G4double momentum_x, momentum_y, momentum_z;
if(direction == 'z')
- {
- // Direction of particle, energy and laboratory angle
- momentum_x = sin(theta) * cos(phi) ;
- momentum_y = sin(theta) * sin(phi) ;
- momentum_z = cos(theta) ;
+ {
+ // Direction of particle, energy and laboratory angle
+ momentum_x = sin(theta) * cos(phi) ;
+ momentum_y = sin(theta) * sin(phi) ;
+ momentum_z = cos(theta) ;
- }
+ }
else if(direction == 'y')
- {
- // Direction of particle, energy and laboratory angle
- momentum_z = sin(theta) * cos(phi) ;
- momentum_x = sin(theta) * sin(phi) ;
- momentum_y = cos(theta) ;
- }
+ {
+ // Direction of particle, energy and laboratory angle
+ momentum_z = sin(theta) * cos(phi) ;
+ momentum_x = sin(theta) * sin(phi) ;
+ momentum_y = cos(theta) ;
+ }
else // = 'x'
- {
- // Direction of particle, energy and laboratory angle
- momentum_y = sin(theta) * cos(phi) ;
- momentum_z = sin(theta) * sin(phi) ;
- momentum_x = cos(theta) ;
- }
-
+ {
+ // Direction of particle, energy and laboratory angle
+ momentum_y = sin(theta) * cos(phi) ;
+ momentum_z = sin(theta) * sin(phi) ;
+ momentum_x = cos(theta) ;
+ }
+
TVector3 Momentum(momentum_x,momentum_y,momentum_z);
return Momentum;
diff --git a/NPSimulation/EventGenerator/EventGeneratorGEFReader.hh b/NPSimulation/EventGenerator/EventGeneratorGEFReader.hh
index 3c29ba31fdabb60f675080e1f41572e35a18470d..c4c9319830a6d0dd858eb02943541d5d3888165b 100644
--- a/NPSimulation/EventGenerator/EventGeneratorGEFReader.hh
+++ b/NPSimulation/EventGenerator/EventGeneratorGEFReader.hh
@@ -35,8 +35,13 @@ using namespace CLHEP;
// NPS headers
#include "VEventGenerator.hh"
#include "ParticleStack.hh"
+#include "Particle.hh"
+
+// NPL headers
+#include "NPReaction.h"
#include "NPInputParser.h"
#include "NPParticle.h"
+#include "TFissionConditions.h"
// ROOT headers
#include "TString.h"
@@ -44,43 +49,50 @@ using namespace CLHEP;
#include "TH1.h"
class EventGeneratorGEFReader : public NPS::VEventGenerator{
-public: // Constructor and destructor
+ public: // Constructor and destructor
EventGeneratorGEFReader() ;
virtual ~EventGeneratorGEFReader();
-public: // Inherit from VEventGenerator Class
+ public: // Inherit from VEventGenerator Class
void ReadConfiguration(NPL::InputParser);
void GenerateEvent(G4Event*) ;
void InitializeRootOutput() ;
TVector3 ShootParticle(double,double, TString);
+ void GetBoostFromTwoBodyReaction(double Ex);
-private: // Source parameter from input file
+ private: // Source parameter from input file
G4int event_ID;
- struct SourceParameters {
- SourceParameters() ;
- G4double m_x0 ; // Vertex Position X
- G4double m_y0 ; // Vertex Position Y
- G4double m_z0 ; // Vertex Position Z
- G4double m_SigmaX ;
- G4double m_SigmaY ;
- G4double m_SigmaZ ;
- G4double m_Boost ;
- TString m_direction ;
- vector<string> m_particleName ;
- string m_FissioningSystem ;
- double m_GEFversion ;
- };
- vector<SourceParameters> m_Parameters ;
+ struct SourceParameters {
+ SourceParameters() ;
+ G4double m_x0 ; // Vertex Position X
+ G4double m_y0 ; // Vertex Position Y
+ G4double m_z0 ; // Vertex Position Z
+ G4double m_SigmaX ;
+ G4double m_SigmaY ;
+ G4double m_SigmaZ ;
+ G4double m_Boost ;
+ TString m_direction ;
+ vector<string> m_particleName ;
+ string m_FissioningSystemName ;
+ double m_GEFversion ;
+ };
+ vector<SourceParameters> m_Parameters ;
ParticleStack* m_ParticleStack ;
ifstream fInputDataFile;
bool HasInputDataFile;
- NPL::Particle *FissioningSystem;
+ NPL::Particle* m_FissioningSystem;
+ NPL::Reaction* m_TwoBodyReaction;
+ bool m_isTwoBody;
int version_shift;
vector<double> LastLine;
vector<string> AllowedParticles;
+ private:
+ TFissionConditions* m_FissionConditions;
+ public:
+ void AttachFissionConditions();
};
#endif
diff --git a/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc b/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc
index 24342b9650b7769df9ea75e202281b8908d8cf9b..fb258d50fd0223a516f150dbf5503961a177b6a3 100644
--- a/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc
@@ -217,7 +217,7 @@ void EventGeneratorIsotropic::GenerateEvent(G4Event*){
}
- Particle particle(par.m_particle, theta,particle_energy,G4ThreeVector(momentum_x, momentum_y, momentum_z),G4ThreeVector(x0, y0, z0));
+ NPS::Particle particle(par.m_particle, theta,particle_energy,G4ThreeVector(momentum_x, momentum_y, momentum_z),G4ThreeVector(x0, y0, z0));
m_ParticleStack->AddParticleToStack(particle);
}
}
diff --git a/NPSimulation/EventGenerator/EventGeneratorMultipleParticle.cc b/NPSimulation/EventGenerator/EventGeneratorMultipleParticle.cc
index b906069387357d26827a8cf2492d1a584975560d..c3692131ae720d84d48b2c16f75abf48b1a2733c 100644
--- a/NPSimulation/EventGenerator/EventGeneratorMultipleParticle.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorMultipleParticle.cc
@@ -205,7 +205,7 @@ void EventGeneratorMultipleParticle::GenerateEvent(G4Event* evt){
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));
+ NPS::Particle particle(m_particle, theta,particle_energy,G4ThreeVector(momentum_x, momentum_y, momentum_z),G4ThreeVector(x0, y0, m_z0));
m_ParticleStack->AddParticleToStack(particle);
}