Adding new event generator reading GEF lmd files

Inputs/EventGenerator/GEFneutron0.source

+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%% 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.

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();

NPSimulation/EventGenerator/CMakeLists.txt

-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 )

NPSimulation/EventGenerator/EventGeneratorGEFReader.cc

+/*****************************************************************************
+ * 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;
+}

NPSimulation/EventGenerator/EventGeneratorGEFReader.hh

+#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