Skip to content
Snippets Groups Projects
Commit 62df19d6 authored by Owen Syrett's avatar Owen Syrett
Browse files

Debugginf Event Generator GEF reader

parent 191caabb
No related branches found
No related tags found
1 merge request!27Draft: [Epic] Preparation of the environement for the new GaseousDetectorScorers...
Pipeline #364150 passed
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Hide whitespace changes
Inline Side-by-side
NPSimulation/EventGenerator/EventGeneratorGEFReader.cc
+ 20
− 92
View file @ 62df19d6
......@@ -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)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment