diff --git a/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc b/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc index 8b72743c198db8020b7662e2fbe338a450423d91..c4e64349249c92cc71176129623eaa79360ce272 100644 --- a/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc +++ b/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); + + } } }