Commit d903cb15 authored by Morfouace's avatar Morfouace
Browse files

Adding a fission decay generator based on a simplified GEF code.

parent 96a69eca
......@@ -55,7 +55,7 @@ TFissionChamberPhysics::TFissionChamberPhysics()
///////////////////////////////////////////////////////////////////////////
/// A usefull method to bundle all operation to add a detector
void TFissionChamberPhysics::AddDetector(TVector3 , string ){
void TFissionChamberPhysics::AddDetector(TVector3 ){
// In That simple case nothing is done
// Typically for more complex detector one would calculate the relevant
// positions (stripped silicon) or angles (gamma array)
......@@ -63,11 +63,11 @@ void TFissionChamberPhysics::AddDetector(TVector3 , string ){
}
///////////////////////////////////////////////////////////////////////////
void TFissionChamberPhysics::AddDetector(double R, double Theta, double Phi, string shape){
void TFissionChamberPhysics::AddDetector(double R, double Theta, double Phi){
// Compute the TVector3 corresponding
TVector3 Pos(R*sin(Theta)*cos(Phi),R*sin(Theta)*sin(Phi),R*cos(Theta));
// Call the cartesian method
AddDetector(Pos,shape);
AddDetector(Pos);
}
///////////////////////////////////////////////////////////////////////////
......@@ -207,8 +207,8 @@ void TFissionChamberPhysics::ReadConfiguration(NPL::InputParser parser) {
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " detectors found " << endl;
vector<string> cart = {"POS","Shape"};
vector<string> sphe = {"R","Theta","Phi","Shape"};
vector<string> cart = {"POS","GasMaterial","Pressure"};
vector<string> sphe = {"R","Theta","Phi","GasMaterial","Pressure"};
for(unsigned int i = 0 ; i < blocks.size() ; i++){
if(blocks[i]->HasTokenList(cart)){
......@@ -216,8 +216,9 @@ void TFissionChamberPhysics::ReadConfiguration(NPL::InputParser parser) {
cout << endl << "//// FissionChamber " << i+1 << endl;
TVector3 Pos = blocks[i]->GetTVector3("POS","mm");
string Shape = blocks[i]->GetString("Shape");
AddDetector(Pos,Shape);
string gas = blocks[i]->GetString("GasMaterial");
double pressure = blocks[i]->GetDouble("Pressure","bar");
AddDetector(Pos);
}
else if(blocks[i]->HasTokenList(sphe)){
if(NPOptionManager::getInstance()->GetVerboseLevel())
......@@ -225,8 +226,9 @@ void TFissionChamberPhysics::ReadConfiguration(NPL::InputParser parser) {
double R = blocks[i]->GetDouble("R","mm");
double Theta = blocks[i]->GetDouble("Theta","deg");
double Phi = blocks[i]->GetDouble("Phi","deg");
string Shape = blocks[i]->GetString("Shape");
AddDetector(R,Theta,Phi,Shape);
string gas = blocks[i]->GetString("GasMaterial");
double pressure = blocks[i]->GetDouble("Pressure","bar");
AddDetector(R,Theta,Phi);
}
else{
cout << "ERROR: check your input file formatting " << endl;
......
......@@ -67,8 +67,8 @@ class TFissionChamberPhysics : public TObject, public NPL::VDetector {
vector<double> Time;
/// A usefull method to bundle all operation to add a detector
void AddDetector(TVector3 POS, string shape);
void AddDetector(double R, double Theta, double Phi, string shape);
void AddDetector(TVector3 POS);
void AddDetector(double R, double Theta, double Phi);
//////////////////////////////////////////////////////////////
// methods inherited from the VDetector ABC class
......
......@@ -1614,17 +1614,11 @@ void GEF::Eva(int Ilh, float Z_CN, float A_CN, float E_INIT, float T, float J_fr
// Shell effects excluded
// SN = LDMass(Zi,Ai-1.,0.0) + LyPair(Zi,Ai-1.) - (LDMass(Zi,Ai,0.0) + LyPair(Zi,Ai)) ;
// SNeff = LDMass(Zi,Ai-1.,0.0) - LDMass(Zi,Ai,0.0) ;
Zf = Zi;
Af = Ai;
Ef = Ei;
while(Ei-SN>E_MIN)
{
// Treat gamma competition
Tm = U_Temp(Zi,Ai,Ei,1,1,Tscale,Econd); // emitting nucleus
Td = U_Temp(Zi,Ai-1,Ei-SNeff,1,1,Tscale,Econd);
if(Ilh > 0) // Emission from fragments
{
Gamma_g = 0.624 * pow(Ai,1.6) * pow(Tm,5); // in meV (Ignatyuk, Bologna)
......@@ -1652,9 +1646,7 @@ void GEF::Eva(int Ilh, float Z_CN, float A_CN, float E_INIT, float T, float J_fr
}
// Reduces the even-odd effect in neutron number
// due to low level density below the pairing gap
Pgamma = Gamma_g / (Gamma_g + Gamma_n);
Pgamma = Gamma_g / (Gamma_g + Gamma_n);
if(rn->Rndm()<Pgamma) //gamma will be emitted
{
In_gamma = In_gamma + 1;
......@@ -4207,7 +4199,7 @@ void GEF::N_Saddle_Scission(void)
J_Frag_light = 0;
// Eva(0,I_Z_light_sad,I_A_light_sad,E_intr_light,TlightFF,J_Frag_light,R_Z_light_sci,
// R_A_light_sci,E_Final_Light,Array_E_n1_ss,Array_Tn,Array_Eg0_light);
Eva(0,I_Z_light_sad,I_A_light_sad,E_intr_light,TlightFF,J_Frag_light,R_Z_light_sci,
Eva(1,I_Z_light_sad,I_A_light_sad,E_intr_light,TlightFF,J_Frag_light,R_Z_light_sci,
R_A_light_sci,E_Final_Light,Array_E_n1_frag1,Array_Tn,Array_Eg0_light);
for(int i = 1;i<= 50;i++)
......@@ -4232,7 +4224,7 @@ void GEF::N_Saddle_Scission(void)
J_Frag_heavy = 0;
// Eva(0,I_Z_heavy_sad,I_A_heavy_sad,E_intr_heavy,TheavyFF,J_Frag_heavy,R_Z_heavy_sci,
// R_A_heavy_sci,E_Final_Heavy,Array_E_n2_ss,Array_Tn,Array_Eg0_heavy);
Eva(0,I_Z_heavy_sad,I_A_heavy_sad,E_intr_heavy,TheavyFF,J_Frag_heavy,R_Z_heavy_sci,
Eva(1,I_Z_heavy_sad,I_A_heavy_sad,E_intr_heavy,TheavyFF,J_Frag_heavy,R_Z_heavy_sci,
R_A_heavy_sci,E_Final_Heavy,Array_E_n2_frag2,Array_Tn,Array_Eg0_heavy);
for(int i = 1;i<=50;i++)
......
#ifndef GEF_CLASS
#define GEF_CLASS
/***********************************************************************
* *
* Original Author : Diego Ramos -> diego.ramos@ganil.fr *
* Adapted for NPTool : Pierre Morfouace -> pierre.morfouace2@cea.fr *
* Creation Date : 25/09/2020 *
* *
*----------------------------------------------------------------------*
* *
* Description: *
* This class generates fission fragments based on a simplified *
* model of the GEF code. *
* This is not the full GEF code. User should use it with caution. *
* *
*----------------------------------------------------------------------*
* Comment: *
* *
* *
************************************************************************/
// C++ header
#include <iostream>
#include <iomanip>
#include <sstream>
......@@ -17,11 +35,11 @@
#include <ctime>
#include <fcntl.h>
#include <unistd.h>
#include <sys/stat.h>
#include <time.h>
#include <typeinfo>
// ROOT header
#include <TROOT.h>
#include <TFile.h>
#include <TDirectory.h>
......@@ -33,15 +51,11 @@
#include <TKey.h>
#include <TRandom3.h>
using namespace std;
//#endif
//#define GEF_SIMPLE
// NPTool headers
// NPTool header
#include "NPParticle.h"
using namespace std;
class GEF
{
public:
......@@ -65,6 +79,8 @@ public:
float* GetNeutronEnergyFrag1() {return Array_E_n1_frag1;}
float* GetNeutronEnergyFrag2() {return Array_E_n2_frag2;}
float* GetGammaEnergyFrag1() {return Array_Eg0_light;}
float* GetGammaEnergyFrag2() {return Array_Eg0_heavy;}
float GetTKE() {return TKE;}
float GetKE1() {return Ekinlight_sci;}
float GetKE2() {return Ekinheavy_sci;}
......
......@@ -44,8 +44,11 @@ NPL::FissionDecay::FissionDecay(std::string compound, std::string fission_model)
void NPL::FissionDecay::ReadConfiguration(NPL::InputParser parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("FissionDecay");
if(blocks.size()>0 && NPOptionManager::getInstance()->GetVerboseLevel())
if(blocks.size()>0 && NPOptionManager::getInstance()->GetVerboseLevel()){
cout << endl << "\033[1;35m//// Fission decay found" << std::endl;
}
if(blocks.size()==0) HasFissionToken=0;
else HasFissionToken=1;
std::vector<std::string> token =
{"CompoundNucleus","FissionModel","Shoot_FF","Shoot_neutron","Shoot_gamma"};
......@@ -133,12 +136,6 @@ bool NPL::FissionDecay::GenerateEvent(string CompoundName, double MEx,double MEK
double Phil = m_FissionModel->GetPhffl();
double Phih = m_FissionModel->GetPhffh();
/*cout << endl;
cout << "A-> " << Al << " / " << Ah << endl;
cout << "Z-> " << Zl << " / " << Zh << endl;
cout << "Phi-> " << Phil << " / " << Phih << endl;
cout << "Theta-> " << Thetal << " / " << Thetah << endl;
cout << "KE-> " << KEl << " / " << KEh << endl;*/
TVector3 Momentuml = Pl * TVector3(sin(Thetal)*cos(Phil),
sin(Thetal)*sin(Phil),
cos(Thetal));
......@@ -158,13 +155,20 @@ bool NPL::FissionDecay::GenerateEvent(string CompoundName, double MEx,double MEK
KE1 = m_FissionModel->GetKE1();
KE2 = m_FissionModel->GetKE2();
// Neutron emission
// Neutron and gamma emission
float* En1;
float* Eg1;
En1 = m_FissionModel->GetNeutronEnergyFrag1();
/*cout << "----- Neutron energy: " << endl;
for(int i=0; i<51; i++) {
cout << "En= " << En1[i] << endl;
}*/
//cout << "----- Neutron energy: " << endl;
//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++){
// cout << "Eg= " << Eg1[i] << endl;
//}
}
}
return worked;
......
......@@ -46,7 +46,7 @@ namespace NPL{
public:
void ReadConfiguration(std::string Path);
void ReadConfiguration(NPL::InputParser parser);
int GetFissionToken() {return HasFissionToken;}
private:
GEF* m_FissionModel;
std::string m_FissionModelName;
......@@ -56,7 +56,7 @@ namespace NPL{
std::vector<std::string> m_FissionFragmentName;
std::vector<NPL::Particle> m_FissionFragment;
std::vector<double> m_FissionFragmentMasses;
int HasFissionToken;
bool m_shoot_FF;
bool m_shoot_neutron;
bool m_shoot_gamma;
......
......@@ -591,7 +591,7 @@ void Scone::ReadSensitive(const G4Event* ){
m_Event->SetCapture(2);
m_Event->SetCaptureTime(GdCaptureTime);
}
//else m_Event->SetCapture(0);
else m_Event->SetCapture(0);
GdProcess_scorer->clear();
///////////
......
......@@ -43,6 +43,10 @@ FissionDecay::FissionDecay(G4String modelName,G4Region* envelope) :
ReadConfiguration();
m_PreviousEnergy=0 ;
m_PreviousLength=0 ;
m_FissionConditions = new TFissionConditions();
//AttachFissionConditions();
//if(!RootOutput::getInstance()->GetTree()->FindBranch("FissionConditions"))
// RootOutput::getInstance()->GetTree()->Branch("FissionConditions", "TFissionConditions", &m_FissionConditions);
}
......@@ -65,22 +69,22 @@ void FissionDecay::AttachFissionConditions(){
void FissionDecay::ReadConfiguration(){
NPL::InputParser input(NPOptionManager::getInstance()->GetReactionFile());
m_FissionDecay.ReadConfiguration(input);
std::string Mother = m_FissionDecay.GetCompoundName();
m_CompoundParticle = NPL::Particle(Mother);
m_CompoundName = NPL::ChangeNameToG4Standard(Mother,true);
m_FissionConditions = new TFissionConditions();
AttachFissionConditions();
if(!RootOutput::getInstance()->GetTree()->FindBranch("FissionConditions"))
RootOutput::getInstance()->GetTree()->Branch("FissionConditions", "TFissionConditions", &m_FissionConditions);
if(m_FissionDecay.GetFissionToken()>0){
std::string Mother = m_FissionDecay.GetCompoundName();
m_CompoundParticle = NPL::Particle(Mother);
m_CompoundName = NPL::ChangeNameToG4Standard(Mother,true);
AttachFissionConditions();
if(!RootOutput::getInstance()->GetTree()->FindBranch("FissionConditions"))
RootOutput::getInstance()->GetTree()->Branch("FissionConditions", "TFissionConditions", &m_FissionConditions);
}
}
////////////////////////////////////////////////////////////////////////////////
G4bool FissionDecay::IsApplicable( const G4ParticleDefinition& particleType) {
m_CurrentName = particleType.GetParticleName();
// Extract Ex from name
if(m_CurrentName.find("[")!=std::string::npos)
......@@ -99,7 +103,7 @@ G4bool FissionDecay::IsApplicable( const G4ParticleDefinition& particleType) {
////////////////////////////////////////////////////////////////////////////////
G4bool FissionDecay::ModelTrigger(const G4FastTrack& fastTrack) {
bool Trigger = true;
m_FissionConditions->Clear();
m_PreviousEnergy=fastTrack.GetPrimaryTrack()->GetKineticEnergy();
......@@ -146,7 +150,7 @@ void FissionDecay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
pdirection.x(),pdirection.y(),pdirection.z(),
FissionFragment, Ex,DEK,DPx,DPy,DPz,
TKE, KE1, KE2);
/////////////////////////////////////////////////
// Fillion the attached Fission condition Tree //
/////////////////////////////////////////////////
......@@ -166,7 +170,7 @@ void FissionDecay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
if(size == 0)
return;
// Get Neutron Multiplicity
int Zsum = FissionFragment[0].GetZ() + FissionFragment[1].GetZ();
int Asum = FissionFragment[0].GetA() + FissionFragment[1].GetA();
......@@ -184,7 +188,7 @@ void FissionDecay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
Momentum=Momentum.unit();
double Brho = FissionFragment[i].GetBrho();
m_FissionConditions->SetFragmentZ(FFZ);
m_FissionConditions->SetFragmentA(FFA);
m_FissionConditions->SetFragmentKineticEnergy(DEK[i]);
......@@ -199,10 +203,10 @@ void FissionDecay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
if(FFZ==0){
if(FFA==1)
FissionFragmentDef=G4ParticleTable::GetParticleTable()->FindParticle("neutron");
else if(FFA==0){
FissionFragmentDef=G4ParticleTable::GetParticleTable()->FindParticle("gamma");
}
}
}
// proton
......@@ -212,11 +216,11 @@ void FissionDecay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
else
FissionFragmentDef=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(FFZ, FFA, Ex[i]);
G4DynamicParticle DynamicFissionFragment(FissionFragmentDef,Momentum,DEK[i]);
fastStep.CreateSecondaryTrack(DynamicFissionFragment, localPosition, time);
}
if(size){
// Set the end of the step conditions
fastStep.SetPrimaryTrackFinalKineticEnergyAndDirection(0,pdirection);
......@@ -225,6 +229,6 @@ void FissionDecay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
fastStep.SetPrimaryTrackFinalTime (time);
fastStep.KillPrimaryTrack();
fastStep.SetPrimaryTrackPathLength(0.0);
}
}
}
......@@ -55,7 +55,8 @@ void Analysis::TreatEvent(){
Einit = InitialConditions->GetKineticEnergy(0);
double init_ThetaLab = ReactionConditions->GetTheta(0)*deg;
double init_BeamEnergy = ReactionConditions->GetBeamEnergy();
neutron->SetKineticEnergy(init_BeamEnergy);
//neutron->SetKineticEnergy(init_BeamEnergy);
neutron->SetKineticEnergy(Einit);
double beam_TOF = neutron->GetTimeOfFlight();
double Xtarget = InitialConditions->GetIncidentPositionX();
......@@ -71,7 +72,6 @@ void Analysis::TreatEvent(){
TVector3 HitPos = DetPos-TargetPos;
//R= HitPos.Mag()*1e-3;
R= Rdet*mm;
Distance.push_back(R);
Det.push_back(m_ChiNu->DetectorNumber[i]);
T.push_back(m_ChiNu->Time[i]);
......
......@@ -49,3 +49,89 @@ PISTA
PHI= 0 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Vamos Configuration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Vamos
% R= 0.3 m
% Theta= 0 deg
%
%Vamos DC
% Z= 60 mm
% Gas= iC4H10
% Pressure= 0.01 bar
% Temperature= 295 kelvin
%
%Vamos DC
% Z= 200 mm
% Gas= iC4H10
% Pressure= 0.01 bar
% Temperature= 295 kelvin
%
%Vamos BeamCatcher
% Material= BC400
% Width= 30 mm
% Length= 80 mm
% Thickness= 30 mm
% Pos= 0 0 600 mm
%
%Vamos MWPPAC
% Z= 750 mm
% Gas= iC4H10
% Pressure= 0.01 bar
% Temperature= 295 kelvin
%
%Vamos DC
% Z= 940 mm
% Gas= iC4H10
% Pressure= 0.01 bar
% Temperature= 295 kelvin
%
%Vamos DC
% Z= 1060 mm
% Gas= iC4H10
% Pressure= 0.01 bar
% Temperature= 295 kelvin
%
%Vamos IC
% Z= 1200 mm
% Thickness= 50 mm
% Gas= CF4
% Pressure= 0.2 bar
% Temperature= 295 kelvin
%
%Vamos IC
% Z= 1250 mm
% Thickness= 50 mm
% Gas= CF4
% Pressure= 0.2 bar
% Temperature= 295 kelvin
%
%Vamos IC
% Z= 1300 mm
% Thickness= 50 mm
% Gas= CF4
% Pressure= 0.2 bar
% Temperature= 295 kelvin
%
%Vamos IC
% Z= 1375 mm
% Thickness= 100 mm
% Gas= CF4
% Pressure= 0.2 bar
% Temperature= 295 kelvin
%
%Vamos IC
% Z= 1525 mm
% Thickness= 200 mm
% Gas= CF4
% Pressure= 0.2 bar
% Temperature= 295 kelvin
%
%Vamos IC
% Z= 1725 mm
% Thickness= 200 mm
% Gas= CF4
% Pressure= 0.2 bar
% Temperature= 295 kelvin
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
......@@ -58,7 +58,6 @@ void Analysis::Init(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){
ReInitValue();
if(InitialConditions->GetParticleMultiplicity()>0){
E_init = InitialConditions->GetKineticEnergy(0);
}
......@@ -70,7 +69,6 @@ void Analysis::TreatEvent(){
new_energy = true;
}
if(new_energy == true){
if(m_DetectedNeutron!=0) vDetectedNeutron.push_back(m_DetectedNeutron);
m_DetectedNeutron = 0;
......
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Target
% THICKNESS= 10 micrometer
% RADIUS= 10 mm
% MATERIAL= CD2
% ANGLE= 0 deg
% X= 0 mm
% Y= 0 mm
% Z= 0 mm
Target
THICKNESS= 0.74 micrometer
RADIUS= 1.5 cm
MATERIAL= 238U
ANGLE= 0 deg
X= 0 mm
Y= 0 mm
Z= 0 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Scone
POS= 0 0 0 mm
Ring1= 1
Ring2= 1
FissionChamber= 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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