diff --git a/Inputs/EventGenerator/GEFneutron0.source b/Inputs/EventGenerator/GEFneutron0.source
new file mode 100755
index 0000000000000000000000000000000000000000..23ff732723c7621acb8a865c6d44f53f205ffc61
--- /dev/null
+++ b/Inputs/EventGenerator/GEFneutron0.source
@@ -0,0 +1,19 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%% An Isotropic Source to be used as EventGenerator %%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Energy are given in MeV , Position in mm %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+GEFReader
+ GEFversion= 2023.33
+ EnergyLow= 0.01
+ EnergyHigh= 20
+ InputDataFile= /home/benoit/Physics/nptool/nptoolv3/Inputs/InputData/test_ng.txt
+ x0= 0
+ y0= 0
+ z0= 0 mm
+ SigmaX= 12 mm
+ SigmaY= 12 mm
+ Particle= fragments
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Supported particle type: proton, neutron, gamma, fragments.
+% If you want to see all particles of the file, put nothing in argument.
diff --git a/NPSimulation/Core/PrimaryGeneratorAction.cc b/NPSimulation/Core/PrimaryGeneratorAction.cc
index 260e99523735b6fd08660fee68d261dd72d5b9d4..c061ec19943462ef39d7df6ff7f6b0df0d180045 100644
--- a/NPSimulation/Core/PrimaryGeneratorAction.cc
+++ b/NPSimulation/Core/PrimaryGeneratorAction.cc
@@ -41,6 +41,7 @@
// Event Generator Class
#include "EventGeneratorIsotropic.hh"
+#include "EventGeneratorGEFReader.hh"
#include "EventGeneratorCosmic.hh"
#include "EventGeneratorMultipleParticle.hh"
#include "EventGeneratorBeam.hh"
@@ -90,6 +91,14 @@ void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path){
m_EventGenerator.push_back(myEventGenerator);
}
blocks.clear();
+ blocks = parser.GetAllBlocksWithToken("GEFReader");
+ if (blocks.size()>0) {
+ NPS::VEventGenerator* myEventGenerator = new EventGeneratorGEFReader();
+ myEventGenerator->ReadConfiguration(parser);
+ myEventGenerator->InitializeRootOutput();
+ m_EventGenerator.push_back(myEventGenerator);
+ }
+ blocks.clear();
blocks = parser.GetAllBlocksWithToken("MultipleParticle");
if (blocks.size()>0) {
NPS::VEventGenerator* myEventGenerator = new EventGeneratorMultipleParticle();
diff --git a/NPSimulation/EventGenerator/CMakeLists.txt b/NPSimulation/EventGenerator/CMakeLists.txt
index 3bcc1361a6f2285bca1cbc0c19b6d88e650c4752..5d5f47b34a146c2b066964e8ffd826aba667b033 100644
--- a/NPSimulation/EventGenerator/CMakeLists.txt
+++ b/NPSimulation/EventGenerator/CMakeLists.txt
@@ -1,2 +1,2 @@
-add_library(NPSEventGenerator OBJECT EventGeneratorBeam.cc EventGeneratorMultipleParticle.cc EventGeneratorCosmic.cc EventGeneratorIsotropic.cc VEventGenerator.cc)
+add_library(NPSEventGenerator OBJECT EventGeneratorBeam.cc EventGeneratorMultipleParticle.cc EventGeneratorCosmic.cc EventGeneratorIsotropic.cc EventGeneratorGEFReader.cc VEventGenerator.cc)
#target_link_libraries(NPSEventGenerator ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPInitialConditions -lNPInteractionCoordinates )
diff --git a/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc b/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc
new file mode 100644
index 0000000000000000000000000000000000000000..0c179868757432e4ef1346c285510fcc3eb5df9c
--- /dev/null
+++ b/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc
@@ -0,0 +1,356 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2016 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Benoit MAUSS contact address: benoit.mauss@cea.fr *
+ * *
+ * Creation Date : October 2024 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This event Generator is used to simulate fission events based on *
+ * data files generated by the GEF model code (General Description of *
+ * Fission Observables) *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+// C++
+#include<limits>
+#include<iostream>
+#include<sstream>
+
+// G4 headers
+#include "G4ParticleTable.hh"
+#include "G4IonTable.hh"
+// G4 headers including CLHEP headers
+// for generating random numbers
+#include "Randomize.hh"
+
+// NPS headers
+#include "EventGeneratorGEFReader.hh"
+
+// NPL headers
+#include "RootOutput.h"
+#include "NPNucleus.h"
+#include "NPOptionManager.h"
+#include "NPFunction.h"
+using namespace CLHEP;
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+EventGeneratorGEFReader::EventGeneratorGEFReader(){
+ m_ParticleStack = ParticleStack::getInstance();
+ event_ID=0;
+
+ HasInputDataFile = false;
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+EventGeneratorGEFReader::~EventGeneratorGEFReader(){
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+EventGeneratorGEFReader::SourceParameters::SourceParameters(){
+
+ m_x0 = 0 ;
+ m_y0 = 0 ;
+ m_z0 = 0 ;
+ m_SigmaX = 0 ;
+ m_SigmaY = 0 ;
+ m_SigmaZ = 0 ;
+ m_direction = 'z' ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void EventGeneratorGEFReader::ReadConfiguration(NPL::InputParser parser){
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("GEFReader");
+ m_Parameters.reserve(blocks.size());
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "\033[1;35m//// GEFReader source found " << endl;
+
+ vector<string> token = {"x0","y0","z0","Particle"};
+ for(unsigned int i = 0 ; i < blocks.size() ; i++){
+ if(blocks[i]->HasTokenList(token)){
+ m_Parameters.push_back(SourceParameters());
+ const vector<SourceParameters>::reverse_iterator it = m_Parameters.rbegin();
+
+/* if(blocks[i]->HasToken("GEFversion"))
+ {
+ it->m_GEFversion=blocks[i]->GetDouble("GEFversion");
+ if(version>=2023.32 && version<=2023.33) version_shift=0;
+ }*/
+ if(blocks[i]->HasToken("InputDataFile")){
+ vector<string> file = blocks[i]->GetVectorString("InputDataFile");
+ cout << "Input file to read: " << file[0] << endl;
+ 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);
+
+ iss >> data;
+ comment=to_string(data);
+ cout << "comment=" << comment << " data=" << data << endl;
+ if(comment!="*" && data>0)
+ {// Storing the first line to be read by GenerateEvent
+ LastLine.push_back(data);
+ while(iss >> data) LastLine.push_back(data);
+ }
+ }
+ }
+
+ 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");
+ vector<string> particleName = blocks[i]->GetVectorString("Particle");
+ for(unsigned int j = 0 ; j < particleName.size() ; j++){
+ if(particleName[j]=="proton"){ it->m_particleName.push_back("proton") ;}
+ else if(particleName[j]=="gamma") { it->m_particleName.push_back("gamma") ;}
+ else if(particleName[j]=="neutron") {it->m_particleName.push_back("neutron") ;}
+ else it->m_particleName.push_back(particleName[j]);
+ }
+
+ 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);
+ }
+ cout << "///////// Warning: Only ";
+ for(auto particle: AllowedParticles) cout << particle << ", " ;
+ cout << "will be simulated" << endl;
+
+ if(AllowedParticles.size()==0) cout << "//////// All particles of the file will be simulated" << endl;
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void EventGeneratorGEFReader::GenerateEvent(G4Event*){
+
+ for(auto& par : m_Parameters) {
+ int fileParticleMultiplicity=-1;
+
+ vector<double> ELab;
+ vector<double> ThetaLab;
+ vector<double> PhiLab;
+ vector<G4ParticleDefinition*> vParticle;
+
+ vector<int> ZFrag;
+ vector<int> AFrag;
+ vector<double> ELabFrag;
+ vector<double> cos_thetaFrag;
+ vector<double> PhiFrag;
+ vector<int> NneutronsEmission; // may be used to correct the energy of the fragments
+ double CN_ExcitationEnergy;
+ int it=0;
+
+ if(HasInputDataFile && !fInputDataFile.eof())
+ {
+ fileParticleMultiplicity=0;
+ // for(int GEFline=0;GEFline<4;GEFline++)
+ while(LastLine.size()==1 || LastLine.at(1)<50 || it==0)
+ {
+ 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);
+
+ if(DataLine.size()>=26 && DataLine.at(1)==92)
+ {
+ // Fragment emission
+ for(int frag=0;frag<2;frag++)
+ {
+ ZFrag .push_back( DataLine.at(3 + frag));
+ AFrag .push_back( DataLine.at(7 + frag));
+ ELabFrag .push_back( DataLine.at(13 + 3*frag));
+ cos_thetaFrag .push_back( DataLine.at(14 + 3*frag));
+ PhiFrag .push_back( DataLine.at(15 + 3*frag));
+ NneutronsEmission.push_back( DataLine.at(21+frag));
+ }
+ CN_ExcitationEnergy = DataLine.at(25);
+
+ if(DataLine.size()>26)
+ {
+ // Particle emission
+ 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:
+
+ ELab .push_back(DataLine.at(27 + idx*3));
+ ThetaLab .push_back(DataLine.at(28 + idx*3));
+ PhiLab .push_back(RandFlat::shoot() * 2 * pi);
+ 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:
+
+ ELab .push_back(DataLine.at(27 + idx*3));
+ ThetaLab .push_back(DataLine.at(28 + idx*3));
+ 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)==1 || DataLine.at(0)==2)
+ && (AllowedParticles.size()==0 || std::find(AllowedParticles.begin(), AllowedParticles.end(), "neutron") != AllowedParticles.end()))
+ for(int it=0;it<(DataLine.size()-1)/4;it++)
+ {// Promt fission neutron treatment
+
+ ELab .push_back(DataLine.at(3+4*it));
+ ThetaLab .push_back(DataLine.at(4+4*it));
+ PhiLab .push_back(RandFlat::shoot() * 2 * pi);
+ vParticle.push_back(G4ParticleTable::GetParticleTable()->FindParticle("neutron"));
+
+ fileParticleMultiplicity++;
+ }
+ 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
+
+ G4double cos_theta = RandFlat::shoot(-1.,1.);
+
+ if((DataLine.at(0)==5 || DataLine.at(0)==8) && to_string(DataLine.at(it+1))=="GS")
+ break;
+ ELab .push_back(DataLine.at(it+1));
+ ThetaLab .push_back( cos_theta ); // Need to add the fragment boost to the gamma emission.
+ PhiLab .push_back(RandFlat::shoot() * 2 * pi);
+ 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;
+ }
+ }
+ else
+ {
+ cout << "Problem with Input data file ? no data to put in" << endl;
+ return;
+ }
+
+ G4double x0, y0, z0;
+ TVector3 Momentum;
+
+ x0 = RandGauss::shoot(par.m_x0,par.m_SigmaX);
+ y0 = RandGauss::shoot(par.m_y0,par.m_SigmaY);
+ z0 = RandGauss::shoot(par.m_z0,par.m_SigmaZ);
+
+ for(int i_m=0;i_m<fileParticleMultiplicity;i_m++)
+ {
+ G4double theta = acos(ThetaLab.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);
+ }
+ 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);
+ }
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+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) ;
+
+ }
+ 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) ;
+ }
+ 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) ;
+ }
+
+ TVector3 Momentum(momentum_x,momentum_y,momentum_z);
+
+ return Momentum;
+}
diff --git a/NPSimulation/EventGenerator/EventGeneratorGEFReader.hh b/NPSimulation/EventGenerator/EventGeneratorGEFReader.hh
new file mode 100644
index 0000000000000000000000000000000000000000..82325809dd307b9d3e03f9efeadd8a023769fdfe
--- /dev/null
+++ b/NPSimulation/EventGenerator/EventGeneratorGEFReader.hh
@@ -0,0 +1,84 @@
+#ifndef EventGeneratorGEFReader_h
+#define EventGeneratorGEFReader_h
+/*****************************************************************************
+ * Copyright (C) 2009-2016 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: matta@lpccaen.in2p3.fr *
+ * *
+ * Creation Date : January 2009 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This event Generator is used to simulated Isotropic ion Source *
+ * Very usefull to figure out Geometric Efficacity of experimental Set-Up *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+// C++ header
+#include <string>
+#include <cmath>
+#include <fstream>
+using namespace std;
+
+using namespace CLHEP;
+
+// G4 headers
+#include "G4Event.hh"
+
+// NPS headers
+#include "VEventGenerator.hh"
+#include "ParticleStack.hh"
+#include "NPInputParser.h"
+
+// ROOT headers
+#include "TString.h"
+#include "TF1.h"
+#include "TH1.h"
+
+class EventGeneratorGEFReader : public NPS::VEventGenerator{
+public: // Constructor and destructor
+ EventGeneratorGEFReader() ;
+ virtual ~EventGeneratorGEFReader();
+
+public: // Inherit from VEventGenerator Class
+ void ReadConfiguration(NPL::InputParser);
+ void GenerateEvent(G4Event*) ;
+ void InitializeRootOutput() ;
+ TVector3 ShootParticle(double,double, TString);
+
+private: // Source parameter from input file
+ G4int event_ID;
+ struct SourceParameters {
+ SourceParameters() ;
+ G4double m_EnergyLow ; // Lower limit of energy range
+ G4double m_EnergyHigh ; // Upper limit of energy range
+ 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 ;
+ TString m_direction ;
+ vector<string> m_particleName ;
+ double m_GEFversion ;
+ };
+ vector<SourceParameters> m_Parameters ;
+ ParticleStack* m_ParticleStack ;
+
+ ifstream fInputDataFile;
+ bool HasInputDataFile;
+
+ int version_shift;
+ vector<double> LastLine;
+
+ vector<string> AllowedParticles;
+
+};
+#endif