adding ejectile to particlestack for GEFReader generator

NPSimulation/EventGenerator/EventGeneratorGEFReader.cc

@@ -160,7 +160,7 @@ void EventGeneratorGEFReader::ReadConfiguration(NPL::InputParser parser){
       {
         it->m_FissioningSystemName=blocks[i]->GetString("FissioningSystem");
         m_FissioningSystem=new NPL::Particle(it->m_FissioningSystemName);
-      	m_FissioningSystem->SetKineticEnergy(0,0,0);
+        m_FissioningSystem->SetKineticEnergy(0,0,0);
       }
       if(blocks[i]->HasToken("SigmaX"))
         it->m_SigmaX=blocks[i]->GetDouble("SigmaX","mm");
@@ -198,8 +198,10 @@ void EventGeneratorGEFReader::GetBoostFromTwoBodyReaction(double Ex){
   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......
@@ -229,6 +231,7 @@ void EventGeneratorGEFReader::GenerateEvent(G4Event*){
     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;
 
@@ -264,12 +267,13 @@ void EventGeneratorGEFReader::GenerateEvent(G4Event*){
           double Ex  = DataLine.at(25); // Excitation energy at fission
           if(m_isTwoBody){
             GetBoostFromTwoBodyReaction(Ex);
+            Ejectile = m_TwoBodyReaction->GetParticle3();
           }
           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);
@@ -371,29 +375,29 @@ void EventGeneratorGEFReader::GenerateEvent(G4Event*){
         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
+          {// 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.;
+            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());
+            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();
+            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"));
+            ELab       .push_back(ELabn);
+            CosThetaLab.push_back(cos_thetan);
+            PhiLab     .push_back(Phin);
+            vParticle  .push_back(G4ParticleTable::GetParticleTable()->FindParticle("neutron"));
 
-          	fileParticleMultiplicity++;
+            fileParticleMultiplicity++;
 
           }
 
@@ -476,6 +480,19 @@ void EventGeneratorGEFReader::GenerateEvent(G4Event*){
         NPS::Particle particle(GeneratedParticle, theta,particle_energy,G4ThreeVector(Momentum.x(), Momentum.y(), Momentum.z()),G4ThreeVector(x0, y0, z0));
         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);
+
+    }
   }
 }