Commit 51e8a25c authored by Morfouace's avatar Morfouace
Browse files

* Updating fission event generator

parent c6c793e9
Pipeline #101374 passed with stages
in 12 minutes and 36 seconds
......@@ -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)
......
......@@ -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++){
......
......@@ -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;
......
......@@ -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;
}
......
......@@ -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);
......
......@@ -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];
......
Markdown is supported
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