From 62df19d6349796ae7ae7ad11b43a556baee897d4 Mon Sep 17 00:00:00 2001
From: "owen.syrett" <owen.syrett@cea.fr>
Date: Thu, 7 Nov 2024 13:45:10 +0100
Subject: [PATCH] Debugginf Event Generator GEF reader
---
.../EventGenerator/EventGeneratorGEFReader.cc | 112 ++++--------------
1 file changed, 20 insertions(+), 92 deletions(-)
diff --git a/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc b/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc
index 70b703588..8b72743c1 100644
--- a/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc
@@ -160,6 +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);
}
if(blocks[i]->HasToken("SigmaX"))
it->m_SigmaX=blocks[i]->GetDouble("SigmaX","mm");
@@ -203,6 +204,7 @@ void EventGeneratorGEFReader::GetBoostFromTwoBodyReaction(double Ex){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorGEFReader::GenerateEvent(G4Event*){
+ m_FissionConditions->Clear();
for(auto& par : m_Parameters) {
int fileParticleMultiplicity=-1;
@@ -217,7 +219,6 @@ void EventGeneratorGEFReader::GenerateEvent(G4Event*){
vector<double> cos_thetaFrag;
vector<double> PhiFrag;
vector<int> NneutronsEmission; // may be used to correct the energy of the fragments
- vector<TVector3> FragmentBoostforN;
double CN_ExcitationEnergy;
TVector3 LightFragmentDirection(0,0,1.);
@@ -291,8 +292,6 @@ void EventGeneratorGEFReader::GenerateEvent(G4Event*){
ImpulsionFrag[frag].Boost(FissioningSystemBoost);
if(frag==0) LightFragmentDirection.SetMagThetaPhi(1.,acos(cos_thetaf),Phif);
- if(frag==1) HeavyFragmentDirection.SetMagThetaPhi(1.,acos(cos_thetaf),Phif);
- TVector3 Boost=Fragment->GetEnergyImpulsion().BoostVector();
ELabf = ImpulsionFrag[frag].E() - Fragment->Mass();
cos_thetaf = ImpulsionFrag[frag].CosTheta();
@@ -303,7 +302,6 @@ void EventGeneratorGEFReader::GenerateEvent(G4Event*){
cos_thetaFrag .push_back( cos_thetaf);
PhiFrag .push_back( Phif);
NneutronsEmission.push_back( DataLine.at(21 + frag));
- FragmentBoostforN.push_back(Boost);
//delete Fragment;
// Filling fission conditions
@@ -373,100 +371,30 @@ 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++)
- {// Prompt fission neutron treatment (just get the angles)
+ {// 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.;
- //std::cout << " Elab 0 " << ELabn << std::endl;
+ TVector3 NeutronLabAngle(0,0,1.);
+ NeutronLabAngle.SetMagThetaPhi(1.,acos(cos_thetan),Phin);
+ NeutronLabAngle.RotateUz(LightFragmentDirection);
+ Neutron->SetKineticEnergy(ELabn,NeutronLabAngle.Theta(),NeutronLabAngle.Phi());
- ELab .push_back(ELabn);
- CosThetaLab.push_back(cos_thetan);
- PhiLab .push_back(Phin);
- vParticle .push_back(G4ParticleTable::GetParticleTable()->FindParticle("neutron"));
+ TLorentzVector ImpulsionNeutron = Neutron->GetEnergyImpulsion();
+ ImpulsionNeutron.Boost(FissioningSystemBoost);
+ ELabn = ImpulsionNeutron.E()-Neutron->Mass();
+ cos_thetan = ImpulsionNeutron.CosTheta();
+ Phin = NeutronLabAngle.Phi();
- fileParticleMultiplicity++;
- }
-
- else if(DataLine.size()>0 && DataLine.at(0)==1
- && (AllowedParticles.size()==0 || std::find(AllowedParticles.begin(), AllowedParticles.end(), "neutron") != AllowedParticles.end()))
- for(int it=0;it<(DataLine.size()-1)/2;it++)
- {// LIGHT fragment prompt fission neutron treatment (apply boost by linking neutron to fragment emitting it)
-
- double ELabn_for_ID = DataLine.at(2+2*it);
-
- //std::cout << " Elab for ID " << ELabn_for_ID << std::endl;
-
- // Find the index of ELabn_for_ID in Elab
- auto iterator = std::find(ELab.begin(), ELab.end(), ELabn_for_ID);
-
- int index;
- if (iterator != ELab.end()) {
- // Calculate the index
- index = std::distance(ELab.begin(), iterator);
- //std::cout << "Index of " << ELabn_for_ID << " is " << index << std::endl;
- } else {
- //std::cout << " " << ELabn_for_ID << " not found in ELab." << std::endl;
- }
-
- double ELabn = ELab.at(index);
- double cos_thetan = CosThetaLab.at(index);
- double Phin = PhiLab.at(index);
-
- 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(FragmentBoostforN.at(0));
- ELabn = ImpulsionNeutron.E()-Neutron->Mass();
- cos_thetan = ImpulsionNeutron.CosTheta();
- Phin = NeutronLabAngle.Phi();
-
- ELab.at(index) = ELabn;
- CosThetaLab.at(index) = cos_thetan;
- PhiLab.at(index) = Phin;
- }
-
- else if(DataLine.size()>0 && DataLine.at(0)==2
- && (AllowedParticles.size()==0 || std::find(AllowedParticles.begin(), AllowedParticles.end(), "neutron") != AllowedParticles.end()))
- for(int it=0;it<(DataLine.size()-1)/2;it++)
- {// HEAVY fragment prompt fission neutron treatment (apply boost by linking neutron to fragment emitting it)
-
- double ELabn_for_ID = DataLine.at(2+2*it);
-
- // Find the index of ELabn_for_ID in Elab
- auto iterator = std::find(ELab.begin(), ELab.end(), ELabn_for_ID);
-
- int index;
- if (iterator != ELab.end()) {
- // Calculate the index
- index = std::distance(ELab.begin(), iterator);
- //std::cout << "Index of " << ELabn_for_ID << " is " << index << std::endl;
- } else {
- //std::cout << " " << ELabn_for_ID << " not found in ELab." << std::endl;
- }
-
- double ELabn = ELab.at(index);
- double cos_thetan = CosThetaLab.at(index);
- double Phin = PhiLab.at(index);
-
- TVector3 NeutronLabAngle(0,0,1.);
- NeutronLabAngle.SetMagThetaPhi(1.,acos(cos_thetan),Phin);
- NeutronLabAngle.RotateUz(HeavyFragmentDirection);
- Neutron->SetKineticEnergy(ELabn,NeutronLabAngle.Theta(),NeutronLabAngle.Phi());
+ ELab .push_back(ELabn);
+ CosThetaLab.push_back(cos_thetan);
+ PhiLab .push_back(Phin);
+ vParticle .push_back(G4ParticleTable::GetParticleTable()->FindParticle("neutron"));
- TLorentzVector ImpulsionNeutron = Neutron->GetEnergyImpulsion();
- ImpulsionNeutron.Boost(FragmentBoostforN.at(1));
- ELabn = ImpulsionNeutron.E()-Neutron->Mass();
- cos_thetan = ImpulsionNeutron.CosTheta();
- Phin = NeutronLabAngle.Phi();
+ fileParticleMultiplicity++;
- ELab.at(index) = ELabn;
- CosThetaLab.at(index) = cos_thetan;
- PhiLab.at(index) = Phin;
}
else if(DataLine.size()>0 && (DataLine.at(0)>=3 && DataLine.at(0)<=8)
--
GitLab