From 9e36743d432eb3dc39d10d4fed385d2c711634b8 Mon Sep 17 00:00:00 2001
From: "owen.syrett" <owen.syrett@cea.fr>
Date: Wed, 6 Nov 2024 09:39:01 +0100
Subject: [PATCH] correction of neutron boost and miscellaneous
---
.../EventGenerator/EventGeneratorGEFReader.cc | 117 +++++++++++++++---
1 file changed, 101 insertions(+), 16 deletions(-)
diff --git a/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc b/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc
index e06b22bbd..21b514ed8 100644
--- a/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorGEFReader.cc
@@ -109,6 +109,7 @@ void EventGeneratorGEFReader::ReadConfiguration(NPL::InputParser parser){
{
getline(fInputDataFile,line);
istringstream iss(line);
+ //cout << "line : " << line << endl;
iss >> data;
comment=to_string(data);
@@ -116,6 +117,7 @@ void EventGeneratorGEFReader::ReadConfiguration(NPL::InputParser parser){
{// Storing the first line to be read by GenerateEvent
LastLine.push_back(data);
while(iss >> data) LastLine.push_back(data);
+ break;
}
}
}
@@ -182,10 +184,12 @@ 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 FissioningSystemBoost=FissioningSystem->GetEnergyImpulsion().BoostVector();
TVector3 LightFragmentDirection(0,0,1.);
+ TVector3 HeavyFragmentDirection(0,0,1.);
TLorentzVector ImpulsionFrag[2];
NPL::Particle *Neutron=new NPL::Particle("1n");
@@ -199,7 +203,7 @@ void EventGeneratorGEFReader::GenerateEvent(G4Event*){
{
fileParticleMultiplicity=0;
// for(int GEFline=0;GEFline<4;GEFline++)
- while(LastLine.size()==1 || LastLine.at(1)<50 || it==0)
+ while(it==0 || LastLine.size()==1 || (LastLine.size()>0 && LastLine.at(1)<50))
{
it++;
string line;
@@ -217,6 +221,11 @@ void EventGeneratorGEFReader::GenerateEvent(G4Event*){
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;
+ //}
+
if(DataLine.size()>=26 && DataLine.at(1)>50)
{
// Fragment emission
@@ -227,14 +236,17 @@ void EventGeneratorGEFReader::GenerateEvent(G4Event*){
double ELabf = DataLine.at(13 + 3*frag);
double cos_thetaf = DataLine.at(14 + 3*frag);
double Phif = DataLine.at(15 + 3*frag) *M_PI/180;
-
- if(frag==0) LightFragmentDirection.SetMagThetaPhi(1.,acos(cos_thetaf),Phif);
Fragment=new NPL::Particle(Zf,Af);
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);
+ 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();
@@ -244,6 +256,7 @@ 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;
}
CN_ExcitationEnergy = DataLine.at(25);
@@ -304,24 +317,15 @@ 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)/4;it++)
- {// Promt fission neutron treatment
+ for(int it=0;it<(DataLine.size()-1)/3;it++)
+ {// Prompt fission neutron treatment (just get the angles)
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();
-
+ //std::cout << " Elab 0 " << ELabn << std::endl;
+
ELab .push_back(ELabn);
CosThetaLab.push_back(cos_thetan);
PhiLab .push_back(Phin);
@@ -329,6 +333,87 @@ void EventGeneratorGEFReader::GenerateEvent(G4Event*){
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());
+
+ TLorentzVector ImpulsionNeutron = Neutron->GetEnergyImpulsion();
+ ImpulsionNeutron.Boost(FragmentBoostforN.at(1));
+ 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)>=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++)
--
GitLab