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;

  m_isTwoBody = false;
  HasInputDataFile = false;
  m_FissionConditions = new TFissionConditions();
}



//....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_Boost        =  0  ;
  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())
    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     (2023.11 <= it->m_GEFversion && it->m_GEFversion <= 2023.22) {
          version_shift=0;
          version_ff_shift=0;
        }
        else if(2023.31 <= it->m_GEFversion && it->m_GEFversion <= 2023.33) {
          version_shift=1;
          version_ff_shift=0;
        }
        else if(2024.11 <= it->m_GEFversion && it->m_GEFversion <= 2024.11){
          version_shift=1;
          version_ff_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;
        int    count=0;

        while(count==0)
        {
          getline(fInputDataFile,line);
          //cout << "line : " << line << endl;
          if (line.find("*")==0) continue;
          
          istringstream iss(line);
          count = std::distance(std::istream_iterator<std::string>(iss), std::istream_iterator<std::string>());
          //cout << "line : " << line << endl;
          
          if(count>0) 
          {// Storing the first line to be read by GenerateEvent
            iss.clear();              
            iss.seekg(0, iss.beg);    
            double data;
            while(iss >> data) {
              //cout << "data = " << data << endl;
              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("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_FissioningSystemName=blocks[i]->GetString("FissioningSystem");
        m_FissioningSystem=new NPL::Particle(it->m_FissioningSystemName);
        m_FissioningSystem->SetKineticEnergy(0,0,0);
      }
      if(blocks[i]->HasToken("BeamProfile"))
        it->m_BeamProfile=blocks[i]->GetString("BeamProfile");
      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","um");  // attention in um

    }
    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::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;
  double Phi3 = Phi4 - pi;
  m_FissioningSystem->SetKineticEnergy(E4,Theta4,Phi4);

  m_TwoBodyReaction->SetParticle3(E3,Theta3,Phi3);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorGEFReader::GenerateEvent(G4Event*){
  m_FissionConditions->Clear();
  for(auto& par : m_Parameters) {
    int fileParticleMultiplicity=-1;

    vector<double>                ELab;
    vector<double>                CosThetaLab;
    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;

    TVector3 LightFragmentDirection(0,0,1.);
    TVector3 HeavyFragmentDirection(0,0,1.);
    TLorentzVector ImpulsionFrag[2];

    NPL::Particle *Neutron=new NPL::Particle("1n");
    NPL::Particle *Proton=new NPL::Particle("1H");
    NPL::Particle *Gamma=new NPL::Particle("gamma");
    NPL::Particle *Fragment;
    NPL::Particle *Ejectile;
    TVector3 FissioningSystemBoost;
    int it=0; 

    int index_Z_CN  =  0 + version_shift;                    // Zsad
    int index_A_CN  =  1 + version_shift;                    // Asad
    int index_Ex    = 24 + version_shift + version_ff_shift; // E@fission
    int data_length = 25 + version_shift + version_ff_shift; // data length of CN and FF, without particle-list info prior to fission
    
    int index_Zf[2];
    int index_Af[2];
    int index_ELabf[2];
    int index_nf[2];
    for(int frag=0; frag<2; frag++){
      index_Zf[frag] = 2 + frag + version_shift + version_ff_shift;
      index_Af[frag] = 6 + frag + version_shift + version_ff_shift;
      index_ELabf[frag] = 12 + 3*frag + version_shift + version_ff_shift;
      index_nf[frag] = 20 + frag + version_shift + version_ff_shift;
    }

    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))
      {
        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;
        //}


        // FIRST LINE OF THE FISSION EVENT : 
        // data on CN, FFl, FFh and possibly pre-fission emission
        if(DataLine.size()>=data_length && DataLine.at(index_Z_CN)>50)
        {
          double Ex  = DataLine.at(index_Ex); // Excitation energy at fission
          if(m_isTwoBody){
            GetBoostFromTwoBodyReaction(Ex);
            Ejectile = m_TwoBodyReaction->GetParticle3();
          }
          FissioningSystemBoost = m_FissioningSystem->GetEnergyImpulsion().BoostVector();

          // Compound nucleus (after pre-fission emission)
          int Z_CN = DataLine.at(index_Z_CN);
          int A_CN = DataLine.at(index_A_CN);

          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(index_Zf[frag]);
            int            Af = DataLine.at(index_Af[frag]); // Post-scission
            double      ELabf = DataLine.at(index_ELabf[frag]);
            double cos_thetaf = DataLine.at(index_ELabf[frag]+1);
            double       Phif = DataLine.at(index_ELabf[frag]+2) *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);

            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(index_nf[frag]));
            //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);
          } // end of Fragment emission
          CN_ExcitationEnergy = DataLine.at(index_Ex);

          // Particle emission before fission
          if(DataLine.size()>data_length)
          {
            string ParticleList = to_string(DataLine.at(data_length));
            int idx=0;
            for(char type: ParticleList) // read one digit of the particle list 
              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(data_length + 1 + idx*3);
                double cos_thetan = DataLine.at(data_length + 2 + 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(data_length + idx*3);
                double cos_thetap = DataLine.at(data_length + 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; // if digit is 0, end of particle emission
          }
        }
        // IF THE OPTIONS OF PROMPT-NEUTRON AND/OR PROMPT GAMMA EMISSION BY FF IS SWITCH ON IN GEF INPUT
        // NEUTRONS: first line has 0 at DataLine.(0) and gives the list of the post-scission neutrons  
        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++)
          {// Promt fission neutron treatment

            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.;

            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(FissioningSystemBoost);
            ELabn      = ImpulsionNeutron.E()-Neutron->Mass();
            cos_thetan = ImpulsionNeutron.CosTheta();
            Phin       = NeutronLabAngle.Phi();

            ELab       .push_back(ELabn);
            CosThetaLab.push_back(cos_thetan);
            PhiLab     .push_back(Phin);
            vParticle  .push_back(G4ParticleTable::GetParticleTable()->FindParticle("neutron"));

            fileParticleMultiplicity++;
          }
        // GAMMA-RAYS: gamma from light (DataLine.at(0) in [3..5]) and heavy (DataLine.at(0) in[6..8])
        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;
      }
    }
    else
    {
      cout << "Problem with Input data file ? no data to put in" << endl;
      return;
    }

    G4double x0, y0, z0;
    TVector3 Momentum;
    
    if(par.m_BeamProfile=="Flat"){
      double rand_r     = par.m_SigmaX*sqrt(RandFlat::shoot());
      double rand_theta = RandFlat::shoot() * 2. * pi;
      x0 = par.m_x0 + rand_r * cos(rand_theta);
      y0 = par.m_y0 + rand_r * sin(rand_theta);
      z0 = RandFlat::shoot(par.m_z0*mm-par.m_SigmaZ*mm, par.m_z0*mm+par.m_SigmaZ*mm);
      //cout << "Circular Flat beam profile: x0 = " << x0 << ", y0 = " << y0 << ", z0 = " << z0 << endl;
    }
    else{ // by default gaussian beam profile is considered 
      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);
    }
    
    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;
	  NPS::Particle particle(GeneratedParticle, theta,particle_energy,G4ThreeVector(Momentum.x(), Momentum.y(), Momentum.z()),G4ThreeVector(x0, y0, z0));
	  m_ParticleStack->AddParticleToStack(particle);
	}

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

	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(m_isTwoBody){
      G4double theta                          = Ejectile->GetImpulsion().Theta();
      G4double phi                            = Ejectile->GetImpulsion().Phi();
      G4double particle_energy                = Ejectile->GetEnergy();
      G4ParticleDefinition* GeneratedParticle = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(Ejectile->GetZ(),Ejectile->GetA(),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);

    }
  }
}

//....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;
}