gwm17/Mask
1
0
Fork 0
mirror of https://github.com/gwm17/Mask.git synced 2024-11-22 18:28:51 -05:00
Code Issues Projects Releases Wiki Activity

Fixed the ReactionSystem inheritance. ReactionSystem is now abstract, added DeacySystem and OneStepSystem.

Browse Source
This commit is contained in:
Gordon McCann 2021-08-24 11:55:32 -04:00
parent 87aa3dba8a
commit a30bfa2555
13 changed files with 369 additions and 146 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
include/DecaySystem.h Normal file
View File

@ -0,0 +1,41 @@
#ifndef DECAYSYSTEM_H
#define DECAYSYSTEM_H
#include "ReactionSystem.h"
#include "AngularDistribution.h"
namespace Mask {
class DecaySystem: public ReactionSystem {
public:
DecaySystem();
DecaySystem(std::vector<int>& z, std::vector<int>& a);
~DecaySystem();
bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
void RunSystem() override;
void SetRandomGenerator(TRandom3* gen) override;
inline void SetDecay1Distribution(const std::string& filename) { decay1dist.ReadDistributionFile(filename); };
inline const Nucleus& GetTarget() const { return step1.GetTarget(); };
inline const Nucleus& GetEjectile() const { return step1.GetEjectile(); };
inline const Nucleus& GetResidual() const { return step1.GetResidual(); };
inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); };
inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); };
private:
void LinkTarget() override;
void SetSystemEquation() override;
Reaction step1;
AngularDistribution decay1dist;
};
};
#endif

4
include/Kinematics.h
View File

@ -2,6 +2,8 @@
#define KINEMATICS_H
#include "ReactionSystem.h"
#include "DecaySystem.h"
#include "OneStepSystem.h"
#include "TwoStepSystem.h"
#include "ThreeStepSystem.h"
#include "Plotter.h"
@ -32,7 +34,7 @@ public:
inline void SetTreeFlag() { save_tree_flag = true; };
inline void SetPlotterFlag() { do_plotter_flag = true; };
inline int GetNumberOfSamples() { return m_nsamples; };
inline const char* GetSystemName() { return sys == nullptr ? "" : sys->GetSystemEquation(); };
inline const char* GetSystemName() { return sys == nullptr ? "" : sys->GetSystemEquation().c_str(); };
inline const char* GetOutputName() { return m_outfile_name.c_str(); };
inline int GetReactionType() { return m_rxn_type; };
void Run();

33
include/OneStepSystem.h Normal file
View File

@ -0,0 +1,33 @@
#ifndef ONESTEPSYSTEM_H
#define ONESTEPSYSTEM_H
#include "ReactionSystem.h"
namespace Mask {
class OneStepSystem: public ReactionSystem {
public:
OneStepSystem();
OneStepSystem(std::vector<int>& z, std::vector<int>& a);
~OneStepSystem();
bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
void RunSystem() override;
inline void SetReactionThetaType(int type) { step1.SetEjectileThetaType(type); };
inline const Nucleus& GetTarget() const { return step1.GetTarget(); };
inline const Nucleus& GetProjectile() const { return step1.GetProjectile(); };
inline const Nucleus& GetEjectile() const { return step1.GetEjectile(); };
inline const Nucleus& GetResidual() const { return step1.GetResidual(); };
private:
void LinkTarget() override;
void SetSystemEquation() override;
Reaction step1;
};
};
#endif

31
include/ReactionSystem.h
View File

@ -10,7 +10,7 @@
#define REACTIONSYSTEM_H
#include "Reaction.h"
#include "AngularDistribution.h"
#include "KinematicsExceptions.h"
#include <vector>
#include <TRandom3.h>
@ -19,47 +19,32 @@ namespace Mask {
class ReactionSystem {
public:
ReactionSystem();
ReactionSystem(std::vector<int>& z, std::vector<int>& a);
virtual ~ReactionSystem();
virtual bool SetNuclei(std::vector<int>& z, std::vector<int>& a);
virtual bool SetNuclei(std::vector<int>& z, std::vector<int>& a) = 0;
virtual void RunSystem() = 0;
void AddTargetLayer(std::vector<int>& zt, std::vector<int>& at, std::vector<int>& stoich, double thickness);
/*Set sampling parameters*/
void SetRandomGenerator(TRandom3* gen);
virtual void SetRandomGenerator(TRandom3* gen);
inline void SetBeamDistro(double mean, double sigma) { m_beamDist = std::make_pair(mean, sigma); };
inline void SetReactionThetaType(int type) { step1.SetEjectileThetaType(type); };
inline void SetTheta1Range(double min, double max) { m_theta1Range = std::make_pair(min*deg2rad, max*deg2rad); };
inline void SetPhi1Range(double min, double max) { m_phi1Range = std::make_pair(min*deg2rad, max*deg2rad); };
inline void SetExcitationDistro(double mean, double sigma) { m_exDist = std::make_pair(mean, sigma); };
inline void SetDecay1Distribution(const std::string& file) { decay1dist.ReadDistributionFile(file); };
inline void SetDecay2Distribution (const std::string& file) { decay2dist.ReadDistributionFile(file); };
virtual void RunSystem();
inline const Nucleus& GetTarget() const { return step1.GetTarget(); };
inline const Nucleus& GetProjectile() const { return step1.GetProjectile(); };
inline const Nucleus& GetEjectile() const { return step1.GetEjectile(); };
inline const Nucleus& GetResidual() const { return step1.GetResidual(); };
inline const char* GetSystemEquation() const { return m_sys_equation.c_str(); };
inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); };
inline int GetDecay2AngularMomentum(){ return decay2dist.GetL(); };
inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); };
inline double GetDecay2BranchingRatio(){ return decay2dist.GetBranchingRatio(); };
inline const std::string& GetSystemEquation() const { return m_sys_equation; };
protected:
virtual void LinkTarget();
virtual void SetSystemEquation();
virtual void LinkTarget() = 0;
virtual void SetSystemEquation() = 0;
Reaction step1;
LayeredTarget target;
//Sampling information
std::pair<double, double> m_beamDist, m_theta1Range, m_phi1Range, m_exDist;
TRandom3* generator; //not owned by ReactionSystem
AngularDistribution decay1dist, decay2dist;
bool target_set_flag, gen_set_flag;
int rxnLayer;
std::string m_sys_equation;

29
include/ThreeStepSystem.h
View File

@ -2,6 +2,7 @@
#define THREESTEPSYSTEM_H
#include "ReactionSystem.h"
#include "AngularDistribution.h"
namespace Mask {
@ -10,19 +11,35 @@ public:
ThreeStepSystem();
ThreeStepSystem(std::vector<int>& z, std::vector<int>& a);
~ThreeStepSystem();
bool SetNuclei(std::vector<int>& z, std::vector<int>& a);
void RunSystem();
bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
void RunSystem() override;
void SetRandomGenerator(TRandom3* gen) override;
inline void SetDecay1Distribution(const std::string& filename) { decay1dist.ReadDistributionFile(filename); };
inline void SetDecay2Distribution(const std::string& filename) { decay2dist.ReadDistributionFile(filename); };
inline void SetReactionThetaType(int type) { step1.SetEjectileThetaType(type); };
inline const Nucleus& GetTarget() const { return step1.GetTarget(); };
inline const Nucleus& GetProjectile() const { return step1.GetProjectile(); };
inline const Nucleus& GetEjectile() const { return step1.GetEjectile(); };
inline const Nucleus& GetResidual() const { return step1.GetResidual(); };
inline const Nucleus& GetBreakup1() const { return step2.GetEjectile(); };
inline const Nucleus& GetBreakup2() const { return step2.GetResidual(); };
inline const Nucleus& GetBreakup3() const { return step3.GetEjectile(); };
inline const Nucleus& GetBreakup4() const { return step3.GetResidual(); };
protected:
void LinkTarget();
void SetSystemEquation();
inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); };
inline int GetDecay2AngularMomentum(){ return decay2dist.GetL(); };
inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); };
inline double GetDecay2BranchingRatio(){ return decay2dist.GetBranchingRatio(); };
Reaction step2, step3;
protected:
void LinkTarget() override;
void SetSystemEquation() override;
Reaction step1, step2, step3;
AngularDistribution decay1dist, decay2dist;
};

25
include/TwoStepSystem.h
View File

@ -2,6 +2,7 @@
#define TWOSTEPSYSTEM_H
#include "ReactionSystem.h"
#include "AngularDistribution.h"
namespace Mask {
@ -10,17 +11,31 @@ public:
TwoStepSystem();
TwoStepSystem(std::vector<int>& z, std::vector<int>& a);
~TwoStepSystem();
bool SetNuclei(std::vector<int>& z, std::vector<int>& a);
void RunSystem();
bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
void RunSystem() override;
void SetRandomGenerator(TRandom3* gen) override;
inline void SetDecay1Distribution(const std::string& filename) { decay1dist.ReadDistributionFile(filename); };
inline void SetReactionThetaType(int type) { step1.SetEjectileThetaType(type); };
inline const Nucleus& GetTarget() const { return step1.GetTarget(); };
inline const Nucleus& GetProjectile() const { return step1.GetProjectile(); };
inline const Nucleus& GetEjectile() const { return step1.GetEjectile(); };
inline const Nucleus& GetResidual() const { return step1.GetResidual(); };
inline const Nucleus& GetBreakup1() const { return step2.GetEjectile(); };
inline const Nucleus& GetBreakup2() const { return step2.GetResidual(); };
inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); };
inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); };
private:
void LinkTarget();
void SetSystemEquation();
void LinkTarget() override;
void SetSystemEquation() override;
Reaction step2;
Reaction step1, step2;
AngularDistribution decay1dist;
};

4
input.txt
View File

@ -1,10 +1,10 @@
----------Data Information----------
OutputFile: /data1/gwm17/TRIUMF_7Bed/simulation/7Bedp_600keV_beam_centered_target_targetgap_BackQQQ_rndmCM_test.root
OutputFile: /data1/gwm17/mask_tests/7Bedp_600keV_beam_centered_target_targetgap_BackQQQ_rndmCM_test.root
SaveTree: yes
SavePlots: yes
----------Reaction Information----------
ReactionType: 2
Z A (order is target, projectile, ejectile, break1, break3)
Z A (order is target, projectile, ejectile, break1, break3(if pure decay is target, ejectile))
1 2
4 7
1 1

75
src/DecaySystem.cpp Normal file
View File

@ -0,0 +1,75 @@
#include "DecaySystem.h"
namespace Mask {
DecaySystem::DecaySystem() :
ReactionSystem()
{
}
DecaySystem::DecaySystem(std::vector<int>& z, std::vector<int>& a) :
ReactionSystem()
{
SetNuclei(z, a);
}
DecaySystem::~DecaySystem() {}
void DecaySystem::SetRandomGenerator(TRandom3* gen) {
generator = gen;
decay1dist.AttachRandomNumberGenerator(gen);
gen_set_flag = true;
}
bool DecaySystem::SetNuclei(std::vector<int>& z, std::vector<int>& a) {
if(z.size() != a.size() || z.size() != 2) {
return false;
}
step1.SetNuclei(z[0], a[0], 0, 0, z[1], a[1]);
SetSystemEquation();
return true;
}
void DecaySystem::LinkTarget() {
step1.SetLayeredTarget(&target);
rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA());
if(rxnLayer != -1) {
step1.SetRxnLayer(rxnLayer);
target_set_flag = true;
} else {
throw ReactionLayerException();
}
}
void DecaySystem::SetSystemEquation() {
m_sys_equation = step1.GetTarget().GetIsotopicSymbol();
m_sys_equation += "-> ";
m_sys_equation += step1.GetEjectile().GetIsotopicSymbol();
m_sys_equation += "+";
m_sys_equation += step1.GetResidual().GetIsotopicSymbol();
}
void DecaySystem::RunSystem() {
if(!gen_set_flag) return;
//Link up the target if it hasn't been done yet
if(!target_set_flag) {
LinkTarget();
}
if(step1.IsDecay()) {
double rxnTheta = std::acos(decay1dist.GetRandomCosTheta());
double rxnPhi = generator->Uniform(0, 2.0*M_PI);
step1.SetPolarRxnAngle(rxnTheta);
step1.SetAzimRxnAngle(rxnPhi);
step1.TurnOnResidualEloss();
step1.Calculate();
} else {
return;
}
}
}

97
src/Kinematics.cpp
View File

@ -43,9 +43,9 @@ bool Kinematics::LoadConfig(const char* filename) {
switch(m_rxn_type) {
case 0:
{
sys = new ReactionSystem();
sys = new DecaySystem();
m_rxn_type = ONESTEP_DECAY;
for(int i=0; i<3; i++) {
for(int i=0; i<2; i++) {
input>>z>>a;
avec.push_back(a);
zvec.push_back(z);
@ -54,7 +54,7 @@ bool Kinematics::LoadConfig(const char* filename) {
}
case 1:
{
sys = new ReactionSystem();
sys = new OneStepSystem();
m_rxn_type = ONESTEP_RXN;
for(int i=0; i<3; i++) {
input>>z>>a;
@ -109,6 +109,7 @@ bool Kinematics::LoadConfig(const char* filename) {
sys->AddTargetLayer(zvec, avec, svec, thickness);
input>>junk;
}
std::cout<<"Reaction equation: "<<GetSystemName()<<std::endl;
double par1, par2;
std::string dfile1, dfile2;
@ -119,7 +120,23 @@ bool Kinematics::LoadConfig(const char* filename) {
input>>junk>>par1>>junk>>par2;
sys->SetBeamDistro(par1, par2);
input>>junk>>par1;
sys->SetReactionThetaType(par1);
switch(m_rxn_type) {
case ONESTEP_RXN :
{
dynamic_cast<OneStepSystem*>(sys)->SetReactionThetaType(par1);
break;
}
case TWOSTEP:
{
dynamic_cast<TwoStepSystem*>(sys)->SetReactionThetaType(par1);
break;
}
case THREESTEP:
{
dynamic_cast<ThreeStepSystem*>(sys)->SetReactionThetaType(par1);
break;
}
}
input>>junk>>par1>>junk>>par2;
sys->SetTheta1Range(par1, par2);
input>>junk>>par1>>junk>>par2;
@ -128,13 +145,35 @@ bool Kinematics::LoadConfig(const char* filename) {
sys->SetExcitationDistro(par1, par2);
input>>junk>>dfile1;
input>>junk>>dfile2;
sys->SetDecay1Distribution(dfile1);
sys->SetDecay2Distribution(dfile2);
switch(m_rxn_type) {
case ONESTEP_DECAY:
{
DecaySystem* this_sys = dynamic_cast<DecaySystem*>(sys);
this_sys->SetDecay1Distribution(dfile1);
std::cout<<"Decay1 angular momentum: "<<this_sys->GetDecay1AngularMomentum()<<std::endl;
std::cout<<"Decay1 total branching ratio: "<<this_sys->GetDecay1BranchingRatio()<<std::endl;
break;
}
case TWOSTEP:
{
TwoStepSystem* this_sys = dynamic_cast<TwoStepSystem*>(sys);
this_sys->SetDecay1Distribution(dfile1);
std::cout<<"Decay1 angular momentum: "<<this_sys->GetDecay1AngularMomentum()<<std::endl;
std::cout<<"Decay1 total branching ratio: "<<this_sys->GetDecay1BranchingRatio()<<std::endl;
break;
}
case THREESTEP:
{
ThreeStepSystem* this_sys = dynamic_cast<ThreeStepSystem*>(sys);
this_sys->SetDecay1Distribution(dfile1);
this_sys->SetDecay2Distribution(dfile2);
std::cout<<"Decay1 angular momentum: "<<this_sys->GetDecay1AngularMomentum()<<" Decay2 angular momentum: "<<this_sys->GetDecay2AngularMomentum()<<std::endl;
std::cout<<"Decay1 total branching ratio: "<<this_sys->GetDecay1BranchingRatio()<<" Decay2 total branching ratio: "<<this_sys->GetDecay2BranchingRatio()<<std::endl;
break;
}
}
sys->SetRandomGenerator(global_generator);
std::cout<<"Reaction equation: "<<GetSystemName()<<std::endl;
std::cout<<"Decay1 angular momentum: "<<sys->GetDecay1AngularMomentum()<<" Decay2 angular momentum: "<<sys->GetDecay2AngularMomentum()<<std::endl;
std::cout<<"Decay1 total branching ratio: "<<sys->GetDecay1BranchingRatio()<<" Decay2 total branching ratio: "<<sys->GetDecay2BranchingRatio()<<std::endl;
std::cout<<"Number of samples: "<<GetNumberOfSamples()<<std::endl;
return true;
@ -186,6 +225,10 @@ void Kinematics::Run() {
}
void Kinematics::RunOneStepRxn() {
OneStepSystem* this_sys = dynamic_cast<OneStepSystem*>(sys);
if(this_sys == nullptr) {
return;
}
TFile* output = TFile::Open(m_outfile_name.c_str(), "RECREATE");
@ -214,20 +257,20 @@ void Kinematics::RunOneStepRxn() {
std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
}
sys->RunSystem();
this_sys->RunSystem();
if(save_tree_flag) {
targ = ConvertNucleus(sys->GetTarget());
proj = ConvertNucleus(sys->GetProjectile());
eject = ConvertNucleus(sys->GetEjectile());
residual = ConvertNucleus(sys->GetResidual());
targ = ConvertNucleus(this_sys->GetTarget());
proj = ConvertNucleus(this_sys->GetProjectile());
eject = ConvertNucleus(this_sys->GetEjectile());
residual = ConvertNucleus(this_sys->GetResidual());
tree->Fill();
}
if(do_plotter_flag) {
plotman.FillData(sys->GetTarget());
plotman.FillData(sys->GetProjectile());
plotman.FillData(sys->GetEjectile());
plotman.FillData(sys->GetResidual());
plotman.FillData(this_sys->GetTarget());
plotman.FillData(this_sys->GetProjectile());
plotman.FillData(this_sys->GetEjectile());
plotman.FillData(this_sys->GetResidual());
}
}
@ -244,6 +287,10 @@ void Kinematics::RunOneStepRxn() {
}
void Kinematics::RunOneStepDecay() {
DecaySystem* this_sys = dynamic_cast<DecaySystem*>(sys);
if(this_sys == nullptr) {
return;
}
TFile* output = TFile::Open(m_outfile_name.c_str(), "RECREATE");
TTree* tree;
@ -267,17 +314,17 @@ void Kinematics::RunOneStepDecay() {
std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
}
sys->RunSystem();
this_sys->RunSystem();
if(save_tree_flag) {
targ = ConvertNucleus(sys->GetTarget());
eject = ConvertNucleus(sys->GetEjectile());
residual = ConvertNucleus(sys->GetResidual());
targ = ConvertNucleus(this_sys->GetTarget());
eject = ConvertNucleus(this_sys->GetEjectile());
residual = ConvertNucleus(this_sys->GetResidual());
tree->Fill();
}
if(do_plotter_flag) {
plotman.FillData(sys->GetTarget());
plotman.FillData(sys->GetEjectile());
plotman.FillData(sys->GetResidual());
plotman.FillData(this_sys->GetTarget());
plotman.FillData(this_sys->GetEjectile());
plotman.FillData(this_sys->GetResidual());
}
}

77
src/OneStepSystem.cpp Normal file
View File

@ -0,0 +1,77 @@
#include "OneStepSystem.h"
namespace Mask {
OneStepSystem::OneStepSystem() :
ReactionSystem()
{
}
OneStepSystem::OneStepSystem(std::vector<int>& z, std::vector<int>& a) :
ReactionSystem()
{
SetNuclei(z, a);
}
OneStepSystem::~OneStepSystem() {}
bool OneStepSystem::SetNuclei(std::vector<int>& z, std::vector<int>& a) {
if(z.size() != a.size() || z.size() != 3) {
return false;
}
step1.SetNuclei(z[0], a[0], z[1], a[1], z[2], a[2]);
SetSystemEquation();
return true;
}
void OneStepSystem::LinkTarget() {
step1.SetLayeredTarget(&target);
rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA());
if(rxnLayer != -1) {
step1.SetRxnLayer(rxnLayer);
target_set_flag = true;
} else {
throw ReactionLayerException();
}
}
void OneStepSystem::SetSystemEquation() {
m_sys_equation = step1.GetTarget().GetIsotopicSymbol();
m_sys_equation += "(";
m_sys_equation += step1.GetProjectile().GetIsotopicSymbol();
m_sys_equation += ", ";
m_sys_equation += step1.GetEjectile().GetIsotopicSymbol();
m_sys_equation += ")";
m_sys_equation += step1.GetResidual().GetIsotopicSymbol();
}
void OneStepSystem::RunSystem() {
if(!gen_set_flag) return;
//Link up the target if it hasn't been done yet
if(!target_set_flag) {
LinkTarget();
}
if(!step1.IsDecay()) {
//Sample parameters
double bke = generator->Gaus(m_beamDist.first, m_beamDist.second);
double rxnTheta = acos(generator->Uniform(cos(m_theta1Range.first), cos(m_theta1Range.second)));
double rxnPhi = generator->Uniform(m_phi1Range.first, m_phi1Range.second);
double residEx = generator->Gaus(m_exDist.first, m_exDist.second);
step1.SetBeamKE(bke);
step1.SetPolarRxnAngle(rxnTheta);
step1.SetAzimRxnAngle(rxnPhi);
step1.SetExcitation(residEx);
step1.TurnOnResidualEloss();
step1.Calculate();
} else {
return;
}
}
}

84
src/ReactionSystem.cpp
View File

@ -9,30 +9,10 @@ ReactionSystem::ReactionSystem() :
{
}
ReactionSystem::ReactionSystem(std::vector<int>& z, std::vector<int>& a) :
m_beamDist(0,0), m_theta1Range(0,0), m_phi1Range(0,0), m_exDist(0,0), generator(nullptr), target_set_flag(false), gen_set_flag(false), rxnLayer(0), m_sys_equation("")
{
SetNuclei(z, a);
}
ReactionSystem::~ReactionSystem() {
}
bool ReactionSystem::SetNuclei(std::vector<int>&z, std::vector<int>& a) {
if(z.size() != a.size() || z.size() < 3) {
return false;
}
step1.SetNuclei(z[0], a[0], z[1], a[1], z[2], a[2]);
SetSystemEquation();
return true;
}
ReactionSystem::~ReactionSystem() {}
void ReactionSystem::SetRandomGenerator(TRandom3* gen) {
generator = gen;
decay1dist.AttachRandomNumberGenerator(gen);
decay2dist.AttachRandomNumberGenerator(gen);
gen_set_flag = true;
}
@ -40,66 +20,4 @@ void ReactionSystem::AddTargetLayer(std::vector<int>& zt, std::vector<int>& at,
target.AddLayer(zt, at, stoich, thickness);
}
void ReactionSystem::LinkTarget() {
step1.SetLayeredTarget(&target);
rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA());
if(rxnLayer != -1) {
step1.SetRxnLayer(rxnLayer);
target_set_flag = true;
} else {
throw ReactionLayerException();
}
}
void ReactionSystem::SetSystemEquation() {
m_sys_equation = step1.GetTarget().GetIsotopicSymbol();
m_sys_equation += "(";
m_sys_equation += step1.GetProjectile().GetIsotopicSymbol();
m_sys_equation += ", ";
m_sys_equation += step1.GetEjectile().GetIsotopicSymbol();
m_sys_equation += ")";
m_sys_equation += step1.GetResidual().GetIsotopicSymbol();
}
void ReactionSystem::RunSystem() {
if(!gen_set_flag) return;
//Link up the target if it hasn't been done yet
if(!target_set_flag) {
LinkTarget();
}
if(step1.IsDecay()) {
double decaycostheta = decay1dist.GetRandomCosTheta();
if(decaycostheta == -10.0) {
step1.ResetEjectile();
step1.ResetResidual();
return;
}
double rxnTheta = std::acos(decay1dist.GetRandomCosTheta());
double rxnPhi = generator->Uniform(0, 2.0*M_PI);
step1.SetPolarRxnAngle(rxnTheta);
step1.SetAzimRxnAngle(rxnPhi);
step1.TurnOnResidualEloss();
step1.Calculate();
} else {
//Sample parameters
double bke = generator->Gaus(m_beamDist.first, m_beamDist.second);
double rxnTheta = acos(generator->Uniform(cos(m_theta1Range.first), cos(m_theta1Range.second)));
double rxnPhi = generator->Uniform(m_phi1Range.first, m_phi1Range.second);
double residEx = generator->Gaus(m_exDist.first, m_exDist.second);
step1.SetBeamKE(bke);
step1.SetPolarRxnAngle(rxnTheta);
step1.SetAzimRxnAngle(rxnPhi);
step1.SetExcitation(residEx);
step1.TurnOnResidualEloss();
step1.Calculate();
}
}
};

7
src/ThreeStepSystem.cpp
View File

@ -18,6 +18,13 @@ ThreeStepSystem::~ThreeStepSystem() {
}
void ThreeStepSystem::SetRandomGenerator(TRandom3* gen) {
generator = gen;
decay1dist.AttachRandomNumberGenerator(gen);
decay2dist.AttachRandomNumberGenerator(gen);
gen_set_flag = true;
}
bool ThreeStepSystem::SetNuclei(std::vector<int>&z, std::vector<int>& a) {
if(z.size() != a.size() || z.size() < 5) {
return false;

8
src/TwoStepSystem.cpp
View File

@ -18,8 +18,14 @@ TwoStepSystem::~TwoStepSystem() {
}
void TwoStepSystem::SetRandomGenerator(TRandom3* gen) {
generator = gen;
decay1dist.AttachRandomNumberGenerator(gen);
gen_set_flag = true;
}
bool TwoStepSystem::SetNuclei(std::vector<int>&z, std::vector<int>& a) {
if(z.size() != a.size() || z.size() < 4) {
if(z.size() != a.size() || z.size() != 4) {
return false;
}
int zr = z[0] + z[1] - z[2];