From 51e8a25ce8ca1d3d36a740e14f473d59a267b99e Mon Sep 17 00:00:00 2001
From: Morfouace <pierre.morfouace@gmail.com>
Date: Thu, 28 Jan 2021 14:05:13 +0100
Subject: [PATCH] * Updating fission event generator
---
NPLib/Physics/GEF.cxx | 4 ++--
NPLib/Physics/NPFissionDecay.cxx | 26 ++++++++++++++++++--------
NPLib/Physics/NPFissionDecay.h | 1 +
NPLib/Physics/NPParticle.cxx | 7 +++++--
NPSimulation/Process/FissionDecay.cc | 1 +
Projects/PISTA/Analysis.cxx | 9 +++++----
6 files changed, 32 insertions(+), 16 deletions(-)
diff --git a/NPLib/Physics/GEF.cxx b/NPLib/Physics/GEF.cxx
index 91fe62db9..1d13bb447 100644
--- a/NPLib/Physics/GEF.cxx
+++ b/NPLib/Physics/GEF.cxx
@@ -4786,7 +4786,7 @@ void GEF::FromCMtoLAB(void)
Bcm_light = sqrt(1.-1./pow(Gcm_light,2));
Bcm_heavy = sqrt(1.-1./pow(Gcm_heavy,2));
- //cout << Ekinlight_sci << " " << Ekinheavy_sci << " " << Ekinlight_sci+ Ekinheavy_sci << endl;
+ //cout << Ekinlight_sci << " " << Ekinheavy_sci << " " << Ekinlight_sci+ Ekinheavy_sci << endl;
KElab_light = KElab_heavy= 0;
KElab_light = (Gfis*Gcm_light*M_light +Gfis*Bfis*Gcm_light*Bcm_light*M_light*cos(Thcm_light))-M_light;
@@ -4801,7 +4801,7 @@ void GEF::FromCMtoLAB(void)
Blab_heavy = sqrt(1.-1./pow(Glab_heavy,2));
vlab_light = Blab_light*29.9792458;
vlab_heavy=Blab_heavy*29.9792458;
-
+
//Note that I am not taking account the angle of the fissioning system
Thlab_light = Thlab_heavy=0;
if(Blab_light>0)
diff --git a/NPLib/Physics/NPFissionDecay.cxx b/NPLib/Physics/NPFissionDecay.cxx
index 10a2a7d8c..e476e2292 100644
--- a/NPLib/Physics/NPFissionDecay.cxx
+++ b/NPLib/Physics/NPFissionDecay.cxx
@@ -51,17 +51,18 @@ void NPL::FissionDecay::ReadConfiguration(NPL::InputParser parser){
else HasFissionToken=1;
std::vector<std::string> token =
- {"CompoundNucleus","FissionModel","Shoot_FF","Shoot_neutron","Shoot_gamma"};
+ {"CompoundNucleus","FissionModel","VamosChargeStates","Shoot_FF","Shoot_neutron","Shoot_gamma"};
unsigned int size = blocks.size();
for(unsigned int i = 0 ; i < size ; i++){
if(blocks[i]->HasTokenList(token)){
m_CompoundName = blocks[i]->GetString("CompoundNucleus");
m_FissionModelName = blocks[i]->GetString("FissionModel");
+ m_VamosChargeStates = blocks[i]->GetInt("VamosChargeStates");
m_shoot_FF = blocks[i]->GetInt("Shoot_FF");
m_shoot_neutron = blocks[i]->GetInt("Shoot_neutron");
m_shoot_gamma = blocks[i]->GetInt("Shoot_gamma");
-
+
m_Compound = NPL::Particle(m_CompoundName);
if(m_FissionModelName=="GEF") m_FissionModel = new GEF(m_Compound);
}
@@ -93,6 +94,7 @@ bool NPL::FissionDecay::GenerateEvent(string CompoundName, double MEx,double MEK
Momentum.Unit();
double Theta = Momentum.Theta();
double Phi = Momentum.Phi();
+ double Lfis = 0;
m_Compound = NPL::Particle(CompoundName);
m_Compound.SetExcitationEnergy(MEx);
@@ -100,7 +102,7 @@ bool NPL::FissionDecay::GenerateEvent(string CompoundName, double MEx,double MEK
if(m_FissionModelName=="GEF"){
if(m_FissionModel->IsValid(m_Compound.GetZ(), m_Compound.GetA())){
worked=true;
- m_FissionModel->InitCompound(MEx,MEK);
+ m_FissionModel->InitCompound(MEx,MEK,Lfis,Theta,Phi);
m_FissionModel->Treat();
int Ah = m_FissionModel->GetAffh();
@@ -110,12 +112,20 @@ bool NPL::FissionDecay::GenerateEvent(string CompoundName, double MEx,double MEK
double KEl = m_FissionModel->GetKEffl();
double KEh = m_FissionModel->GetKEffh();
+ double Brhol = m_FissionModel->GetBrhoffl();
+ double Brhoh = m_FissionModel->GetBrhoffh();
NPL::Particle FFl = NPL::Particle(Zl,Al);
NPL::Particle FFh = NPL::Particle(Zh,Ah);
- FFl.SetKineticEnergy(KEl);
- FFh.SetKineticEnergy(KEh);
-
+ if(m_VamosChargeStates==1){
+ // Include Charge states distribtuion from Baron et al. NIM 328 (1993) 177-182
+ FFl.SetBrho(Brhol);
+ FFh.SetBrho(Brhoh);
+ }
+ else{
+ FFl.SetKineticEnergy(KEl);
+ FFh.SetKineticEnergy(KEh);
+ }
FissionFragments.push_back(FFl);
FissionFragments.push_back(FFh);
Ex.push_back(0);
@@ -150,7 +160,7 @@ bool NPL::FissionDecay::GenerateEvent(string CompoundName, double MEx,double MEK
DPy.push_back(Momentumh.Y());
DPz.push_back(Momentuml.Z());
DPz.push_back(Momentumh.Z());
-
+
TKE = m_FissionModel->GetTKE();
KE1 = m_FissionModel->GetKE1();
KE2 = m_FissionModel->GetKE2();
@@ -163,7 +173,7 @@ bool NPL::FissionDecay::GenerateEvent(string CompoundName, double MEx,double MEK
//for(int i=0; i<51; i++) {
// cout << "En= " << En1[i] << endl;
//}
-
+
Eg1 = m_FissionModel->GetGammaEnergyFrag1();
//cout << "----- Gamma energy: " << endl;
//for(int i=0; i<101; i++){
diff --git a/NPLib/Physics/NPFissionDecay.h b/NPLib/Physics/NPFissionDecay.h
index e11d1a901..9be1dda30 100644
--- a/NPLib/Physics/NPFissionDecay.h
+++ b/NPLib/Physics/NPFissionDecay.h
@@ -57,6 +57,7 @@ namespace NPL{
std::vector<NPL::Particle> m_FissionFragment;
std::vector<double> m_FissionFragmentMasses;
int HasFissionToken;
+ bool m_VamosChargeStates;
bool m_shoot_FF;
bool m_shoot_neutron;
bool m_shoot_gamma;
diff --git a/NPLib/Physics/NPParticle.cxx b/NPLib/Physics/NPParticle.cxx
index 3a1e91d44..66125ce2a 100644
--- a/NPLib/Physics/NPParticle.cxx
+++ b/NPLib/Physics/NPParticle.cxx
@@ -464,8 +464,11 @@ void Particle::GetParticleName() {
void Particle::EnergyToBrho(double Q){
if(Q==-1000)
Q=GetZ();
-
- fBrho = sqrt(pow(fKineticEnergy,2) + 2*fKineticEnergy*Mass()) * 1e6 * e_SI / (c_light*1e6) / (Q * e_SI);
+
+ EnergyToBeta();
+ BetaToGamma();
+ //fBrho = sqrt(pow(fKineticEnergy,2) + 2*fKineticEnergy*Mass()) * 1e6 * e_SI / (c_light*1e6) / (Q * e_SI);
+ fBrho = 3.107*GetA()/Q*fBeta*fGamma;
}
diff --git a/NPSimulation/Process/FissionDecay.cc b/NPSimulation/Process/FissionDecay.cc
index 55a4688d1..5ecbeb96b 100644
--- a/NPSimulation/Process/FissionDecay.cc
+++ b/NPSimulation/Process/FissionDecay.cc
@@ -160,6 +160,7 @@ void FissionDecay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
m_FissionConditions->SetA_CN(m_CompoundParticle.GetA());
m_FissionConditions->SetEx_CN(m_ExcitationEnergy);
m_FissionConditions->SetELab_CN(energy);
+ m_FissionConditions->SetThetaLab_CN(pdirection.theta()*180./3.1415);
// Fission Process
m_FissionConditions->Set_TKE(TKE);
diff --git a/Projects/PISTA/Analysis.cxx b/Projects/PISTA/Analysis.cxx
index 44cb37e3a..2f69bdf94 100644
--- a/Projects/PISTA/Analysis.cxx
+++ b/Projects/PISTA/Analysis.cxx
@@ -58,10 +58,10 @@ void Analysis::Init(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){
ReInitValue();
- //OriginalThetaLab = ReactionConditions->GetTheta(0);
- //OriginalElab = ReactionConditions->GetKineticEnergy(0);
- //OriginalBeamEnergy = ReactionConditions->GetBeamEnergy();
- //OriginalEx = ReactionConditions->GetExcitation4();
+ OriginalThetaLab = ReactionConditions->GetTheta(0);
+ OriginalElab = ReactionConditions->GetKineticEnergy(0);
+ OriginalBeamEnergy = ReactionConditions->GetBeamEnergy();
+ OriginalEx = ReactionConditions->GetExcitation4();
int mult = InteractionCoordinates->GetDetectedMultiplicity();
if(mult>0){
@@ -85,6 +85,7 @@ void Analysis::TreatEvent(){
BeamEnergy = 1428.;//InitialConditions->GetIncidentInitialKineticEnergy();
BeamEnergy = U238C.Slow(BeamEnergy,TargetThickness*0.5,0);
Transfer->SetBeamEnergy(BeamEnergy);
+ //Transfer->SetBeamEnergy(OriginalBeamEnergy);
if(PISTA->EventMultiplicity==1){
for(unsigned int i = 0; i<PISTA->EventMultiplicity; i++){
double Energy = PISTA->DE[i] + PISTA->E[i];
--
GitLab