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

Overhaul input and reaction system construction. Much more streamlined, now can have more sampling parameters (chained decay to excited states, etc.)

Browse Source
This commit is contained in:
Gordon McCann 2022-09-01 09:52:09 -04:00
parent 5a8f8730f5
commit 76ae86d595
16 changed files with 469 additions and 394 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
input.txt
View File

@ -1,12 +1,38 @@
----------Data Information----------
OutputFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root
----------Reaction Information----------
begin_nuclei(Z,A)
NumberOfSamples: 100000
begin_chain
begin_step
Type: Reaction
begin_nuclei
5 10
2 3
2 4
end_nuclei
BeamEnergyMean(MeV): 24.0
BeamEnergySigma(MeV): 0.0
ThetaType: Lab
ThetaMin(deg): 15.0
ThetaMax(deg): 15.0
PhiMin(deg): 0.0
PhMax(deg): 0.0
ResidualExcitationMean(MeV): 16.8
ResidualExcitationSigma(MeV): 0.023
end_step
begin_step
Type: Decay
begin_nuclei
5 9
2 4
end_nuclei(Z,A)
end_nuclei
PhiMin(deg): 0.0
PhMax(deg): 360.0
ResidualExcitationMean(MeV): 0.0
ResidualExcitationSigma(MeV): 0.0
AngularDistributionFile: ./etc/isotropic_dist.txt
end_step
end_chain
----------Target Information----------
NumberOfLayers: 1
begin_layer
@ -15,13 +41,4 @@ begin_layer
element 5 10 1
end_elements
end_layer
----------Sampling Information----------
NumberOfSamples: 1000000
BeamMeanEnergy(MeV): 24.0 BeamEnergySigma(MeV): 0.0
EjectileThetaType(0=Lab,1=CM): 0
EjectileThetaMin(deg): 15.0 EjectileThetaMax(deg): 15.0
EjectilePhiMin(deg): 0.0 EjectilePhiMax(deg): 0.0
ResidualExMean(MeV): 16.8 ResidualExSigma(MeV): 0.023
Decay1_DistributionFile: ./etc/isotropic_dist.txt
Decay2_DistributionFile: ./etc/isotropic_dist.txt
--------------------------------------

2
src/Mask/AngularDistribution.cpp
View File

@ -13,7 +13,7 @@ namespace Mask {
}
AngularDistribution::AngularDistribution(const std::string& file) :
m_branchingRatio(1.0), m_L(0), m_isIsotropic(true)
m_uniformCosineDist(-1.0, 1.0), m_branchingRatio(1.0), m_L(0), m_isIsotropic(true)
{
ReadDistributionFile(file);
}

59
src/Mask/DecaySystem.cpp
View File

@ -5,49 +5,52 @@
namespace Mask {
DecaySystem::DecaySystem() :
DecaySystem::DecaySystem(const std::vector<StepParameters>& params) :
ReactionSystem()
{
m_nuclei.resize(3);
}
DecaySystem::DecaySystem(const std::vector<int>& z, const std::vector<int>& a) :
ReactionSystem()
{
m_nuclei.resize(3);
SetNuclei(z, a);
Init(params);
}
DecaySystem::~DecaySystem() {}
bool DecaySystem::SetNuclei(const std::vector<int>& z, const std::vector<int>& a)
void DecaySystem::Init(const std::vector<StepParameters>& params)
{
if(z.size() != a.size() || z.size() != 2)
return false;
if(params.size() != 1 || params[0].rxnType != RxnType::Decay ||
params[0].Z.size() != 2 || params[0].A.size() != 2)
{
m_isValid = false;
std::cerr << "Invalid parameters at DecaySystem::Init(), does not match Decay signature!" << std::endl;
return;
}
int zr = z[0] - z[1];
int ar = a[0] - a[1];
const StepParameters& step1Params = params[0];
m_nuclei[0] = CreateNucleus(z[0], a[0]); //target
m_nuclei[1] = CreateNucleus(z[1], a[1]); //breakup1
int zr = step1Params.Z[0] - step1Params.Z[1];
int ar = step1Params.A[0] - step1Params.A[1];
m_nuclei[0] = CreateNucleus(step1Params.Z[0], step1Params.A[0]); //target
m_nuclei[1] = CreateNucleus(step1Params.Z[1], step1Params.A[1]); //breakup1
m_nuclei[2] = CreateNucleus(zr, ar); //breakup2
m_step1.BindNuclei(&(m_nuclei[0]), nullptr, &(m_nuclei[1]), &(m_nuclei[2]));
SetSystemEquation();
return true;
AddPhiRange(step1Params.phiMin, step1Params.phiMax);
AddDecayAngularDistribution(step1Params.angularDistFile);
AddExcitationDistribution(step1Params.meanResidualEx, step1Params.sigmaResidualEx);
return;
}
std::vector<Nucleus>* DecaySystem::GetNuclei()
void DecaySystem::SetLayeredTarget(const LayeredTarget& target)
{
return &m_nuclei;
}
void DecaySystem::LinkTarget() {
m_step1.SetLayeredTarget(&m_target);
m_target = target;
m_rxnLayer = m_target.FindLayerContaining(m_nuclei[0].Z, m_nuclei[0].A);
if(m_rxnLayer != m_target.GetNumberOfLayers())
{
m_step1.SetLayeredTarget(&m_target);
m_step1.SetRxnLayer(m_rxnLayer);
m_isTargetSet = true;
}
@ -66,14 +69,14 @@ namespace Mask {
void DecaySystem::RunSystem()
{
//Link up the target if it hasn't been done yet
if(!m_isTargetSet)
LinkTarget();
std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
double rxnTheta = m_decayAngularDistributions[0].GetRandomCosTheta();
double rxnPhi = m_phiRanges[0](gen);
double ex = m_exDistributions[0](gen);
double rxnTheta = std::acos(m_step1Distribution.GetRandomCosTheta());
double rxnPhi = (*m_phi1Range)(RandomGenerator::GetInstance().GetGenerator());
m_step1.SetPolarRxnAngle(rxnTheta);
m_step1.SetAzimRxnAngle(rxnPhi);
m_step1.SetExcitation(ex);
m_step1.SetResidualEnergyLoss(true);
m_step1.Calculate();
}

21
src/Mask/DecaySystem.h
View File

@ -6,29 +6,20 @@
namespace Mask {
class DecaySystem: public ReactionSystem {
class DecaySystem: public ReactionSystem
{
public:
DecaySystem();
DecaySystem(const std::vector<int>& z, const std::vector<int>& a);
DecaySystem(const std::vector<StepParameters>& params);
~DecaySystem();
bool SetNuclei(const std::vector<int>& z, const std::vector<int>& a) override;
void RunSystem() override;
std::vector<Nucleus>*GetNuclei() override;
virtual void SetDecay1Distribution(const std::string& filename) override { m_step1Distribution.ReadDistributionFile(filename); }
int GetDecay1AngularMomentum() { return m_step1Distribution.GetL(); }
double GetDecay1BranchingRatio() { return m_step1Distribution.GetBranchingRatio(); }
virtual void SetLayeredTarget(const LayeredTarget& target) override;
virtual void RunSystem() override;
private:
void LinkTarget() override;
void Init(const std::vector<StepParameters>& params);
void SetSystemEquation() override;
Reaction m_step1;
AngularDistribution m_step1Distribution;
};
}

119
src/Mask/MaskApp.cpp
View File

@ -10,7 +10,7 @@ namespace Mask {
MaskApp::MaskApp() :
m_system(nullptr)
{
std::cout<<"----------GWM Kinematics Simulation----------"<<std::endl;
std::cout<<"----------Monte Carlo Simulation of Kinematics----------"<<std::endl;
}
MaskApp::~MaskApp()
@ -29,40 +29,87 @@ namespace Mask {
}
std::string junk;
getline(input, junk);
std::getline(input, junk);
input>>junk>>m_outputName;
std::getline(input, junk);
std::getline(input, junk);
input>>junk>>m_nsamples;
std::vector<int> avec, zvec, svec;
int z, a, s;
getline(input, junk);
getline(input, junk);
std::vector<StepParameters> params;
int z, a;
while(input>>junk)
{
if(junk == "begin_nuclei(Z,A)")
if(junk == "begin_chain")
continue;
else if (junk == "end_nuclei(Z,A)")
else if (junk == "end_chain")
break;
else if(junk == "begin_step")
{
StepParameters currentParams;
input >> junk >> junk;
currentParams.rxnType = StringToRxnType(junk);
if(currentParams.rxnType == RxnType::Reaction)
{
input >> junk;
for(int i=0; i<3; i++)
{
input >> z >> a;
currentParams.Z.push_back(z);
currentParams.A.push_back(a);
}
input >> junk;
input >> junk >> currentParams.meanBeamEnergy;
input >> junk >> currentParams.sigmaBeamEnergy;
input >> junk >> junk;
currentParams.thetaType = StringToRxnThetaType(junk);
input >> junk >> currentParams.thetaMin;
input >> junk >> currentParams.thetaMax;
input >> junk >> currentParams.phiMin;
input >> junk >> currentParams.phiMax;
input >> junk >> currentParams.meanResidualEx;
input >> junk >> currentParams.sigmaResidualEx;
params.push_back(currentParams);
}
else if(currentParams.rxnType == RxnType::Decay)
{
input >> junk;
for(int i=0; i<2; i++)
{
input >> z >> a;
currentParams.Z.push_back(z);
currentParams.A.push_back(a);
}
input >> junk;
input >> junk >> currentParams.phiMin;
input >> junk >> currentParams.phiMax;
input >> junk >> currentParams.meanResidualEx;
input >> junk >> currentParams.sigmaResidualEx;
input >> junk >> currentParams.angularDistFile;
params.push_back(currentParams);
}
else
{
z = std::stoi(junk);
input>>a;
zvec.push_back(z);
avec.push_back(a);
std::cerr << "Invalid reaction information at MaskApp::LoadConfig!" << std::endl;
return false;
}
}
}
m_system = CreateSystem(zvec, avec);
if(m_system == nullptr)
m_system = CreateSystem(params);
if(m_system == nullptr || !m_system->IsValid())
{
std::cerr<<"Param size: "<<params.size()<<std::endl;
std::cerr<<"Failure to parse reaction system... configuration not loaded."<<std::endl;
return false;
}
std::getline(input, junk);
std::getline(input, junk);
LayeredTarget target;
int nlayers;
double thickness;
getline(input, junk);
getline(input, junk);
std::vector<int> avec, zvec, svec;
int s;
input>>junk>>nlayers;
for(int i=0; i<nlayers; i++)
{
@ -80,42 +127,26 @@ namespace Mask {
input>>z>>a>>s;
zvec.push_back(z); avec.push_back(a); svec.push_back(s);
}
m_system->AddTargetLayer(zvec, avec, svec, thickness);
target.AddLayer(zvec, avec, svec, thickness);
input>>junk;
}
m_system->SetLayeredTarget(target);
std::cout<<"Outputing data to file: "<<m_outputName<<std::endl;
std::cout<<"Reaction equation: "<<GetSystemName()<<std::endl;
double par1, par2;
std::string dfile1, dfile2;
std::string thetaTypeString;
getline(input, junk);
getline(input, junk);
input>>junk>>m_nsamples;
input>>junk>>par1>>junk>>par2;
m_system->SetBeamDistro(par1, par2);
input>>junk>>thetaTypeString;
m_system->SetReactionThetaType(StringToRxnThetaType(thetaTypeString));
input>>junk>>par1>>junk>>par2;
m_system->SetTheta1Range(par1, par2);
input>>junk>>par1>>junk>>par2;
m_system->SetPhi1Range(par1, par2);
input>>junk>>par1>>junk>>par2;
m_system->SetExcitationDistro(par1, par2);
input>>junk>>dfile1;
input>>junk>>dfile2;
m_system->SetDecay1Distribution(dfile1);
m_system->SetDecay2Distribution(dfile2);
std::cout<<"Number of samples: "<<GetNumberOfSamples()<<std::endl;
return true;
}
bool MaskApp::SaveConfig(const std::string& filename) { return true; }
//Not implemented... yet!
bool MaskApp::SaveConfig(const std::string& filename)
{
return true;
}
void MaskApp::Run() {
void MaskApp::Run()
{
std::cout<<"Running simulation..."<<std::endl;
if(m_system == nullptr)
{

71
src/Mask/OneStepSystem.cpp
View File

@ -5,51 +5,56 @@
namespace Mask {
OneStepSystem::OneStepSystem() :
OneStepSystem::OneStepSystem(const std::vector<StepParameters>& params) :
ReactionSystem()
{
m_nuclei.resize(4);
}
OneStepSystem::OneStepSystem(const std::vector<int>& z, const std::vector<int>& a) :
ReactionSystem()
{
m_nuclei.resize(4);
SetNuclei(z, a);
Init(params);
}
OneStepSystem::~OneStepSystem() {}
bool OneStepSystem::SetNuclei(const std::vector<int>& z, const std::vector<int>& a)
void OneStepSystem::Init(const std::vector<StepParameters>& params)
{
if(z.size() != a.size() || z.size() != 3)
return false;
if(params.size() != 1 || params[0].rxnType != RxnType::Reaction ||
params[0].Z.size() != 3 || params[0].A.size() != 3)
{
m_isValid = false;
std::cerr << "Invalid parameters at OneStepSystem::Init(), does not match OneStep signature!" << std::endl;
return;
}
int zr = z[0] + z[1] - z[2];
int ar = a[0] + a[1] - a[2];
const StepParameters& step1Params = params[0];
m_nuclei[0] = CreateNucleus(z[0], a[0]); //target
m_nuclei[1] = CreateNucleus(z[1], a[1]); //projectile
m_nuclei[2] = CreateNucleus(z[2], a[2]); //ejectile
//Set nuclei
int zr = step1Params.Z[0] + step1Params.Z[1] - step1Params.Z[2];
int ar = step1Params.A[0] + step1Params.A[1] - step1Params.A[2];
m_nuclei[0] = CreateNucleus(step1Params.Z[0], step1Params.A[0]); //target
m_nuclei[1] = CreateNucleus(step1Params.Z[1], step1Params.A[1]); //projectile
m_nuclei[2] = CreateNucleus(step1Params.Z[2], step1Params.A[2]); //ejectile
m_nuclei[3] = CreateNucleus(zr, ar); //residual
m_step1.BindNuclei(&(m_nuclei[0]), &(m_nuclei[1]), &(m_nuclei[2]), &(m_nuclei[3]));
SetSystemEquation();
return true;
//Set sampling parameters
AddBeamDistribution(step1Params.meanBeamEnergy, step1Params.sigmaBeamEnergy);
m_step1.SetEjectileThetaType(step1Params.thetaType);
AddThetaRange(step1Params.thetaMin, step1Params.thetaMax);
AddPhiRange(step1Params.phiMin, step1Params.phiMax);
AddExcitationDistribution(step1Params.meanResidualEx, step1Params.sigmaResidualEx);
}
std::vector<Nucleus>* OneStepSystem::GetNuclei()
void OneStepSystem::SetLayeredTarget(const LayeredTarget& target)
{
return &m_nuclei;
}
void OneStepSystem::LinkTarget()
{
m_step1.SetLayeredTarget(&m_target);
m_target = target;
m_rxnLayer = m_target.FindLayerContaining(m_nuclei[0].Z, m_nuclei[0].A);
if(m_rxnLayer != m_target.GetNumberOfLayers())
{
m_step1.SetLayeredTarget(&m_target);
m_step1.SetRxnLayer(m_rxnLayer);
m_isTargetSet = true;
}
@ -67,18 +72,14 @@ namespace Mask {
m_sysEquation = stream.str();
}
void OneStepSystem::RunSystem() {
//Link up the target if it hasn't been done yet
if(!m_isTargetSet)
void OneStepSystem::RunSystem()
{
LinkTarget();
}
//Sample parameters
double bke = (*m_beamDist)(RandomGenerator::GetInstance().GetGenerator());
double rxnTheta = std::acos((*m_theta1Range)(RandomGenerator::GetInstance().GetGenerator()));
double rxnPhi = (*m_phi1Range)(RandomGenerator::GetInstance().GetGenerator());
double residEx = (*m_exDist)(RandomGenerator::GetInstance().GetGenerator());
std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
double bke = (m_beamDistributions[0])(gen);
double rxnTheta = std::acos((m_thetaRanges[0])(gen));
double rxnPhi = (m_phiRanges[0])(gen);
double residEx = (m_exDistributions[0])(gen);
m_step1.SetBeamKE(bke);
m_step1.SetPolarRxnAngle(rxnTheta);

14
src/Mask/OneStepSystem.h
View File

@ -5,24 +5,20 @@
namespace Mask {
class OneStepSystem: public ReactionSystem {
class OneStepSystem: public ReactionSystem
{
public:
OneStepSystem();
OneStepSystem(const std::vector<int>& z, const std::vector<int>& a);
OneStepSystem(const std::vector<StepParameters>& params);
~OneStepSystem();
bool SetNuclei(const std::vector<int>& z, const std::vector<int>& a) override;
virtual void SetLayeredTarget(const LayeredTarget& target) override;
void RunSystem() override;
std::vector<Nucleus>* GetNuclei() override;
virtual void SetReactionThetaType(RxnThetaType type) override { m_step1.SetEjectileThetaType(type); };
private:
void LinkTarget() override;
void Init(const std::vector<StepParameters>& params);
void SetSystemEquation() override;
Reaction m_step1;
};
}

5
src/Mask/Reaction.cpp
View File

@ -187,14 +187,15 @@ namespace Mask {
//Calculate in CM, where decay is isotropic
void Reaction::CalculateDecay()
{
double Q = m_target->vec4.M() - m_ejectile->groundStateMass - m_residual->groundStateMass;
double residualMass = m_residual->groundStateMass + m_ex;
double Q = m_target->vec4.M() - m_ejectile->groundStateMass - residualMass;
if(Q < 0)
throw QValueException();
ROOT::Math::Boost boost(m_target->vec4.BoostToCM());
m_target->vec4 = boost*m_target->vec4;
double ejectE_cm = (m_ejectile->groundStateMass*m_ejectile->groundStateMass -
m_residual->groundStateMass*m_residual->groundStateMass + m_target->vec4.E()*m_target->vec4.E()) /
residualMass*residualMass + m_target->vec4.E()*m_target->vec4.E()) /
(2.0*m_target->vec4.E());
double ejectP_cm = std::sqrt(ejectE_cm*ejectE_cm - m_ejectile->groundStateMass*m_ejectile->groundStateMass);

96
src/Mask/ReactionSystem.cpp
View File

@ -10,41 +10,91 @@
namespace Mask {
ReactionSystem::ReactionSystem() :
m_beamDist(nullptr), m_theta1Range(nullptr), m_phi1Range(nullptr), m_exDist(nullptr), m_isTargetSet(false),
m_rxnLayer(0), m_sysEquation("")
m_isTargetSet(false), m_isValid(true), m_rxnLayer(0), m_sysEquation("")
{
}
ReactionSystem::~ReactionSystem()
{
delete m_beamDist;
delete m_theta1Range;
delete m_phi1Range;
delete m_exDist;
}
void ReactionSystem::AddTargetLayer(const std::vector<int>& zt, const std::vector<int>& at, const std::vector<int>& stoich, double thickness)
ReactionSystem* CreateSystem(const std::vector<StepParameters>& params)
{
m_target.AddLayer(zt, at, stoich, thickness);
switch(params.size())
{
case 1:
{
if(params[0].rxnType == RxnType::Decay)
return new DecaySystem(params);
else if (params[0].rxnType == RxnType::Reaction)
return new OneStepSystem(params);
}
ReactionSystem* CreateSystem(const std::vector<int>& z, const std::vector<int>& a)
{
if(z.size() != a.size())
{
std::cerr<<"Size of Z list does not equal size of A list!"<<std::endl;
return nullptr;
}
switch(a.size())
{
case RxnSize::DecaySize: return new DecaySystem(z, a);
case RxnSize::OneStepSize: return new OneStepSystem(z, a);
case RxnSize::TwoStepSize: return new TwoStepSystem(z, a);
case RxnSize::ThreeStepSize: return new ThreeStepSystem(z, a);
case 2: return new TwoStepSystem(params);
case 3: return new ThreeStepSystem(params);
}
return nullptr;
}
/*Set sampling parameters*/
//Should only ever really be one of these, but they cannot be set a priori
void ReactionSystem::AddBeamDistribution(double mean, double sigma)
{
if(mean == -1.0 || sigma == -1.0)
{
m_isValid = false;
std::cerr << "Invalid parameters at ReactionSystem::AddBeamDistribution(), distribution is invalid -> mean: " << mean
<< " sigma: " << sigma << std::endl;
return;
}
m_beamDistributions.emplace_back(mean, sigma);
}
//Again only really one of these, but same issue as beam. This is the sampling of the primary reaction/decay (step 1)
void ReactionSystem::AddThetaRange(double min, double max)
{
if(min == -1.0 || max == -1.0 || min > max)
{
m_isValid = false;
std::cerr << "Invalid parameters at ReactionSystem::AddThetaRange(), range is invalid -> min: " << min
<< " max: " << max << std::endl;
return;
}
m_thetaRanges.emplace_back(std::cos(min*s_deg2rad), std::cos(max*s_deg2rad));
}
void ReactionSystem::AddPhiRange(double min, double max)
{
if(min == -1.0 || max == -1.0 || min > max)
{
m_isValid = false;
std::cerr << "Invalid parameters at ReactionSystem::AddPhiRange(), range is invalid -> min: " << min
<< " max: " << max << std::endl;
return;
}
m_phiRanges.emplace_back(min*s_deg2rad, max*s_deg2rad);
}
//Each reaction step can generate an excited nucleus (for a reaction can make an excited residual, decay can make an excited
//breakup2 or "heavy")
void ReactionSystem::AddExcitationDistribution(double mean, double sigma)
{
if(mean == -1.0 || sigma == -1.0)
{
m_isValid = false;
std::cerr << "Invalid parameters at ReactionSystem::AddExcitationDistribution(), distribution is invalid -> mean: " << mean
<< " sigma: " << sigma << std::endl;
return;
}
m_exDistributions.emplace_back(mean, sigma);
}
//Apply angular distributions to decay products. Here apply them to the breakup1 or "light" fragment
void ReactionSystem::AddDecayAngularDistribution(const std::string& filename)
{
if(filename.empty())
m_decayAngularDistributions.emplace_back();
else
m_decayAngularDistributions.emplace_back(filename);
}
}

74
src/Mask/ReactionSystem.h
View File

@ -12,75 +12,69 @@
#include "Reaction.h"
#include "KinematicsExceptions.h"
#include "RxnType.h"
#include "AngularDistribution.h"
#include <vector>
#include <random>
namespace Mask {
struct StepParameters
{
RxnType rxnType = RxnType::None;
std::vector<int> Z;
std::vector<int> A;
double meanBeamEnergy = -1.0;
double sigmaBeamEnergy = -1.0;
RxnThetaType thetaType = RxnThetaType::None;
double thetaMin = -1.0;
double thetaMax = -1.0;
double phiMin = -1.0;
double phiMax = -1.0;
double meanResidualEx = -1.0;
double sigmaResidualEx = -1.0;
std::string angularDistFile;
};
class ReactionSystem
{
public:
ReactionSystem();
virtual ~ReactionSystem();
virtual bool SetNuclei(const std::vector<int>& z, const std::vector<int>& a) = 0;
virtual void SetLayeredTarget(const LayeredTarget& target) = 0;
virtual void RunSystem() = 0;
virtual std::vector<Nucleus>* GetNuclei() = 0;
virtual void SetReactionThetaType(RxnThetaType type) {}
virtual void SetDecay1Distribution(const std::string& filename) {}
virtual void SetDecay2Distribution(const std::string& filename) {}
void AddTargetLayer(const std::vector<int>& zt, const std::vector<int>& at, const std::vector<int>& stoich, double thickness);
/*Set sampling parameters*/
void SetBeamDistro(double mean, double sigma)
{
if(m_beamDist)
delete m_beamDist;
m_beamDist = new std::normal_distribution<double>(mean, sigma);
}
void SetTheta1Range(double min, double max)
{
if(m_theta1Range)
delete m_theta1Range;
m_theta1Range = new std::uniform_real_distribution<double>(std::cos(min*s_deg2rad), std::cos(max*s_deg2rad));
}
void SetPhi1Range(double min, double max)
{
if(m_phi1Range)
delete m_phi1Range;
m_phi1Range = new std::uniform_real_distribution<double>(min*s_deg2rad, max*s_deg2rad);
}
void SetExcitationDistro(double mean, double sigma)
{
if(m_exDist)
delete m_exDist;
m_exDist = new std::normal_distribution<double>(mean, sigma);
}
std::vector<Nucleus>* GetNuclei() { return &m_nuclei; }
const std::string& GetSystemEquation() const { return m_sysEquation; }
bool IsValid() const { return m_isValid; }
protected:
virtual void LinkTarget() = 0;
virtual void SetSystemEquation() = 0;
void AddBeamDistribution(double mean, double sigma);
void AddThetaRange(double min, double max);
void AddPhiRange(double min, double max);
void AddExcitationDistribution(double mean, double sigma);
void AddDecayAngularDistribution(const std::string& filename);
LayeredTarget m_target;
//Sampling information
std::normal_distribution<double> *m_beamDist, *m_exDist;
std::uniform_real_distribution<double> *m_theta1Range, *m_phi1Range;
std::vector<std::normal_distribution<double>> m_beamDistributions, m_exDistributions;
std::vector<std::uniform_real_distribution<double>> m_thetaRanges, m_phiRanges;
std::vector<AngularDistribution> m_decayAngularDistributions;
bool m_isTargetSet;
bool m_isValid;
std::size_t m_rxnLayer;
std::string m_sysEquation;
std::vector<Nucleus> m_nuclei;
static constexpr double s_deg2rad = M_PI/180.0;
};
ReactionSystem* CreateSystem(const std::vector<int>& z, const std::vector<int>& a);
ReactionSystem* CreateSystem(const std::vector<StepParameters>& params);
}
#endif

46
src/Mask/RxnType.h
View File

@ -7,11 +7,9 @@ namespace Mask {
enum class RxnType
{
PureDecay=0,
OneStepRxn=1,
TwoStepRxn=2,
ThreeStepRxn=3,
None=4
Decay,
Reaction,
None
};
enum class RxnThetaType
@ -29,44 +27,16 @@ namespace Mask {
ThreeStepSize = 5
};
static RxnType GetRxnTypeFromString(const std::string& type_str)
static RxnType StringToRxnType(const std::string& type)
{
if (type_str == "PureDecay")
return RxnType::PureDecay;
else if (type_str == "OneStepRxn")
return RxnType::OneStepRxn;
else if (type_str == "TwoStepRxn")
return RxnType::TwoStepRxn;
else if (type_str == "ThreeStepRxn")
return RxnType::ThreeStepRxn;
if (type == "Decay")
return RxnType::Decay;
else if (type == "Reaction")
return RxnType::Reaction;
else
return RxnType::None;
}
static std::string GetStringFromRxnType(RxnType type)
{
switch(type)
{
case RxnType::PureDecay: return "PureDecay";
case RxnType::OneStepRxn: return "OneStepRxn";
case RxnType::TwoStepRxn: return "TwoStepRxn";
case RxnType::ThreeStepRxn: return "ThreeStepRxn";
case RxnType::None : return "None";
}
return "None";
}
static RxnType GetRxnTypeFromInt(uint32_t value)
{
return static_cast<RxnType>(value);
}
static uint32_t GetIntFromRxnType(RxnType type)
{
return static_cast<uint32_t>(type);
}
static RxnThetaType StringToRxnThetaType(const std::string& type)
{
if(type == "CenterOfMass")

127
src/Mask/ThreeStepSystem.cpp
View File

@ -4,40 +4,59 @@
namespace Mask {
ThreeStepSystem::ThreeStepSystem() :
ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
ThreeStepSystem::ThreeStepSystem(const std::vector<StepParameters>& params) :
ReactionSystem()
{
m_nuclei.resize(8);
}
ThreeStepSystem::ThreeStepSystem(const std::vector<int>& z, const std::vector<int>& a) :
ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
{
m_nuclei.resize(8);
SetNuclei(z, a);
Init(params);
}
ThreeStepSystem::~ThreeStepSystem() {}
bool ThreeStepSystem::SetNuclei(const std::vector<int>& z, const std::vector<int>& a)
void ThreeStepSystem::Init(const std::vector<StepParameters>& params)
{
if(z.size() != a.size() || z.size() < 5)
return false;
if(params.size() != 3 || params[0].rxnType != RxnType::Reaction || params[1].rxnType != RxnType::Decay ||
params[2].rxnType != RxnType::Decay || params[0].Z.size() != 3 || params[0].A.size() != 3 || params[1].Z.size() != 2 ||
params[1].A.size() != 2 || params[2].Z.size() != 2 || params[2].A.size() != 2)
{
m_isValid = false;
std::cerr << "Invalid parameters at ThreeStepSystem::Init(), does not match ThreeStep signature!" << std::endl;
return;
}
int zr = z[0] + z[1] - z[2];
int zb2 = zr - z[3];
int zb4 = zb2 - z[4];
int ar = a[0] + a[1] - a[2];
int ab2 = ar - a[3];
int ab4 = ab2 - a[4];
const StepParameters& step1Params = params[0];
const StepParameters& step2Params = params[1];
const StepParameters& step3Params = params[2];
m_nuclei[0] = CreateNucleus(z[0], a[0]); //target
m_nuclei[1] = CreateNucleus(z[1], a[1]); //projectile
m_nuclei[2] = CreateNucleus(z[2], a[2]); //ejectile
//Setup nuclei
int zr = step1Params.Z[0] + step1Params.Z[1] - step1Params.Z[2];
int ar = step1Params.A[0] + step1Params.A[1] - step1Params.A[2];
if(zr != step2Params.Z[0] || ar != step2Params.A[0])
{
m_isValid = false;
std::cerr << "Invalid parameters at ThreeStepSystem::Init(), step one and step two are not sequential! Step one recoil (Z,A): ("
<< zr << "," << ar << ") Step two target (Z,A): (" << step2Params.Z[0] << "," << step2Params.A[0] << ")" <<std::endl;
return;
}
int zb2 = step2Params.Z[0] - step2Params.Z[1];
int ab2 = step2Params.A[0] - step2Params.A[1];
if(zb2 != step3Params.Z[0] || ab2 != step3Params.A[0])
{
m_isValid = false;
std::cerr << "Invalid parameters at ThreeStepSystem::Init(), step two and step three are not sequential! Step two heavy (Z,A): ("
<< zb2 << "," << ab2 << ") Step three target (Z,A): (" << step3Params.Z[0] << "," << step3Params.A[0] << ")" <<std::endl;
return;
}
int zb4 = step3Params.Z[0] - step3Params.Z[1];
int ab4 = step3Params.A[0] - step3Params.A[1];
m_nuclei[0] = CreateNucleus(step1Params.Z[0], step1Params.A[0]); //target
m_nuclei[1] = CreateNucleus(step1Params.Z[1], step1Params.A[1]); //projectile
m_nuclei[2] = CreateNucleus(step1Params.Z[2], step1Params.A[2]); //ejectile
m_nuclei[3] = CreateNucleus(zr, ar); //residual
m_nuclei[4] = CreateNucleus(z[3], a[3]); //breakup1
m_nuclei[4] = CreateNucleus(step2Params.Z[1], step2Params.A[1]); //breakup1
m_nuclei[5] = CreateNucleus(zb2, ab2); //breakup2
m_nuclei[5] = CreateNucleus(z[4], a[4]); //breakup3
m_nuclei[5] = CreateNucleus(step3Params.Z[1], step3Params.A[1]); //breakup3
m_nuclei[5] = CreateNucleus(zb4, ab4); //breakup4
m_step1.BindNuclei(&(m_nuclei[0]), &(m_nuclei[1]), &(m_nuclei[2]), &(m_nuclei[3]));
@ -45,28 +64,40 @@ namespace Mask {
m_step3.BindNuclei(&(m_nuclei[5]), nullptr, &(m_nuclei[6]), &(m_nuclei[7]));
SetSystemEquation();
return true;
//Step one sampling parameters
AddBeamDistribution(step1Params.meanBeamEnergy, step1Params.sigmaBeamEnergy);
m_step1.SetEjectileThetaType(step1Params.thetaType);
AddThetaRange(step1Params.thetaMin, step1Params.thetaMax);
AddPhiRange(step1Params.phiMin, step1Params.phiMax);
AddExcitationDistribution(step1Params.meanResidualEx, step1Params.sigmaResidualEx);
//Step two sampling parameters
AddPhiRange(step2Params.phiMin, step2Params.phiMax);
AddDecayAngularDistribution(step2Params.angularDistFile);
AddExcitationDistribution(step2Params.meanResidualEx, step2Params.sigmaResidualEx);
//Step three sampling parameters
AddPhiRange(step3Params.phiMin, step3Params.phiMax);
AddDecayAngularDistribution(step3Params.angularDistFile);
AddExcitationDistribution(step3Params.meanResidualEx, step3Params.sigmaResidualEx);
}
std::vector<Nucleus>* ThreeStepSystem::GetNuclei()
void ThreeStepSystem::SetLayeredTarget(const LayeredTarget& target)
{
return &m_nuclei;
}
void ThreeStepSystem::LinkTarget()
m_target = target;
m_rxnLayer = m_target.FindLayerContaining(m_nuclei[0].Z, m_nuclei[0].A);
if(m_rxnLayer != m_target.GetNumberOfLayers())
{
m_step1.SetLayeredTarget(&m_target);
m_step2.SetLayeredTarget(&m_target);
m_step3.SetLayeredTarget(&m_target);
m_rxnLayer = m_target.FindLayerContaining(m_nuclei[0].Z, m_nuclei[0].A);
if(m_rxnLayer != m_target.GetNumberOfLayers())
{
m_step1.SetRxnLayer(m_rxnLayer);
m_step2.SetRxnLayer(m_rxnLayer);
m_step3.SetRxnLayer(m_rxnLayer);
m_isTargetSet = true;
} else
}
else
throw ReactionLayerException();
}
@ -84,22 +115,22 @@ namespace Mask {
m_sysEquation = stream.str();
}
void ThreeStepSystem::RunSystem() {
//Link up the target if it hasn't been done yet
if(!m_isTargetSet)
LinkTarget();
void ThreeStepSystem::RunSystem()
{
//Sample parameters
double bke = (*m_beamDist)(RandomGenerator::GetInstance().GetGenerator());
double rxnTheta = acos((*m_theta1Range)(RandomGenerator::GetInstance().GetGenerator()));
double rxnPhi = (*m_phi1Range)(RandomGenerator::GetInstance().GetGenerator());
double decay1costheta = m_step2Distribution.GetRandomCosTheta();
std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
double bke = m_beamDistributions[0](gen);
double rxnTheta = std::acos((m_thetaRanges[0])(gen));
double rxnPhi = m_phiRanges[0](gen);
double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta();
double decay1Theta = std::acos(decay1costheta);
double decay1Phi = m_phi2Range(RandomGenerator::GetInstance().GetGenerator());
double decay2costheta = m_step3Distribution.GetRandomCosTheta();
double decay1Phi = m_phiRanges[1](gen);
double decay2costheta = m_decayAngularDistributions[1].GetRandomCosTheta();
double decay2Theta = std::acos(decay2costheta);
double decay2Phi = m_phi2Range(RandomGenerator::GetInstance().GetGenerator());
double residEx = (*m_exDist)(RandomGenerator::GetInstance().GetGenerator());
double decay2Phi = m_phiRanges[2](gen);
double residEx = m_exDistributions[0](gen);
double decay1Ex = m_exDistributions[1](gen);
double decay2Ex = m_exDistributions[2](gen);
m_step1.SetBeamKE(bke);
m_step1.SetPolarRxnAngle(rxnTheta);
@ -108,9 +139,11 @@ namespace Mask {
m_step2.SetPolarRxnAngle(decay1Theta);
m_step2.SetAzimRxnAngle(decay1Phi);
m_step2.SetExcitation(decay1Ex);
m_step3.SetPolarRxnAngle(decay2Theta);
m_step3.SetAzimRxnAngle(decay2Phi);
m_step3.SetExcitation(decay2Ex);
m_step1.Calculate();

26
src/Mask/ThreeStepSystem.h
View File

@ -6,35 +6,21 @@
namespace Mask {
class ThreeStepSystem : public ReactionSystem {
class ThreeStepSystem : public ReactionSystem
{
public:
ThreeStepSystem();
ThreeStepSystem(const std::vector<int>& z, const std::vector<int>& a);
ThreeStepSystem(const std::vector<StepParameters>& params);
~ThreeStepSystem();
bool SetNuclei(const std::vector<int>& z, const std::vector<int>& a) override;
void RunSystem() override;
std::vector<Nucleus>* GetNuclei() override;
virtual void SetDecay1Distribution(const std::string& filename) override { m_step2Distribution.ReadDistributionFile(filename); };
virtual void SetDecay2Distribution(const std::string& filename) override { m_step3Distribution.ReadDistributionFile(filename); };
virtual void SetReactionThetaType(RxnThetaType type) override { m_step1.SetEjectileThetaType(type); };
int GetDecay1AngularMomentum() { return m_step2Distribution.GetL(); };
int GetDecay2AngularMomentum(){ return m_step3Distribution.GetL(); };
double GetDecay1BranchingRatio() { return m_step2Distribution.GetBranchingRatio(); };
double GetDecay2BranchingRatio(){ return m_step3Distribution.GetBranchingRatio(); };
virtual void SetLayeredTarget(const LayeredTarget& target) override;
virtual void RunSystem() override;
protected:
void LinkTarget() override;
void Init(const std::vector<StepParameters>& params);
void SetSystemEquation() override;
std::uniform_real_distribution<double> m_phi2Range;
Reaction m_step1, m_step2, m_step3;
AngularDistribution m_step2Distribution, m_step3Distribution;
};
}

95
src/Mask/TwoStepSystem.cpp
View File

@ -6,57 +6,73 @@
namespace Mask {
TwoStepSystem::TwoStepSystem() :
ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
TwoStepSystem::TwoStepSystem(const std::vector<StepParameters>& params) :
ReactionSystem()
{
m_nuclei.resize(6);
}
TwoStepSystem::TwoStepSystem(const std::vector<int>& z, const std::vector<int>& a) :
ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
{
m_nuclei.resize(6);
SetNuclei(z, a);
Init(params);
}
TwoStepSystem::~TwoStepSystem() {}
bool TwoStepSystem::SetNuclei(const std::vector<int>&z, const std::vector<int>& a)
void TwoStepSystem::Init(const std::vector<StepParameters>& params)
{
if(z.size() != a.size() || z.size() != 4)
return false;
if(params.size() != 2 || params[0].rxnType != RxnType::Reaction || params[1].rxnType != RxnType::Decay ||
params[0].Z.size() != 3 || params[0].A.size() != 3 || params[1].Z.size() != 2 || params[1].A.size() != 2)
{
m_isValid = false;
std::cerr << "Invalid parameters at TwoStepSystem::Init(), does not match TwoStep signature!" << std::endl;
return;
}
int zr = z[0] + z[1] - z[2];
int ar = a[0] + a[1] - a[2];
int zb = zr - z[3];
int ab = ar - a[3];
const StepParameters& step1Params = params[0];
const StepParameters& step2Params = params[1];
m_nuclei[0] = CreateNucleus(z[0], a[0]); //target
m_nuclei[1] = CreateNucleus(z[1], a[1]); //projectile
m_nuclei[2] = CreateNucleus(z[2], a[2]); //ejectile
//Setup nuclei
int zr = step1Params.Z[0] + step1Params.Z[1] - step1Params.Z[2];
int ar = step1Params.A[0] + step1Params.A[1] - step1Params.A[2];
if(zr != step2Params.Z[0] || ar != step2Params.A[0])
{
m_isValid = false;
std::cerr << "Invalid parameters at TwoStepSystem::Init(), step one and step two are not sequential! Step one recoil (Z,A): ("
<< zr << "," << ar << ") Step two target (Z,A): (" << step2Params.Z[0] << "," << step2Params.A[0] << ")" <<std::endl;
return;
}
int zb = step2Params.Z[0] - step2Params.Z[1];
int ab = step2Params.A[0] - step2Params.A[1];
m_nuclei[0] = CreateNucleus(step1Params.Z[0], step1Params.A[0]); //target
m_nuclei[1] = CreateNucleus(step1Params.Z[1], step1Params.A[1]); //projectile
m_nuclei[2] = CreateNucleus(step1Params.Z[2], step1Params.A[2]); //ejectile
m_nuclei[3] = CreateNucleus(zr, ar); //residual
m_nuclei[4] = CreateNucleus(z[3], a[3]); //breakup1
m_nuclei[4] = CreateNucleus(step2Params.Z[1], step2Params.A[1]); //breakup1
m_nuclei[5] = CreateNucleus(zb, ab); //breakup2
m_step1.BindNuclei(&(m_nuclei[0]), &(m_nuclei[1]), &(m_nuclei[2]), &(m_nuclei[3]));
m_step2.BindNuclei(&(m_nuclei[3]), nullptr, &(m_nuclei[4]), &(m_nuclei[5]));
SetSystemEquation();
return true;
//Step one sampling parameters
AddBeamDistribution(step1Params.meanBeamEnergy, step1Params.sigmaBeamEnergy);
m_step1.SetEjectileThetaType(step1Params.thetaType);
AddThetaRange(step1Params.thetaMin, step1Params.thetaMax);
AddPhiRange(step1Params.phiMin, step1Params.phiMax);
AddExcitationDistribution(step1Params.meanResidualEx, step1Params.sigmaResidualEx);
//Step two sampling parameters
AddPhiRange(step2Params.phiMin, step2Params.phiMax);
AddDecayAngularDistribution(step2Params.angularDistFile);
AddExcitationDistribution(step2Params.meanResidualEx, step2Params.sigmaResidualEx);
}
std::vector<Nucleus>* TwoStepSystem::GetNuclei()
void TwoStepSystem::SetLayeredTarget(const LayeredTarget& target)
{
return &m_nuclei;
}
void TwoStepSystem::LinkTarget()
{
m_step1.SetLayeredTarget(&m_target);
m_step2.SetLayeredTarget(&m_target);
m_target = target;
m_rxnLayer = m_target.FindLayerContaining(m_nuclei[0].Z, m_nuclei[0].A);
if(m_rxnLayer != m_target.GetNumberOfLayers())
{
m_step1.SetLayeredTarget(&m_target);
m_step2.SetLayeredTarget(&m_target);
m_step1.SetRxnLayer(m_rxnLayer);
m_step2.SetRxnLayer(m_rxnLayer);
m_isTargetSet = true;
@ -79,18 +95,16 @@ namespace Mask {
void TwoStepSystem::RunSystem()
{
//Link up the target if it hasn't been done yet
if(!m_isTargetSet)
LinkTarget();
//Sample parameters
double bke = (*m_beamDist)(RandomGenerator::GetInstance().GetGenerator());
double rxnTheta = acos((*m_theta1Range)(RandomGenerator::GetInstance().GetGenerator()));
double rxnPhi = (*m_phi1Range)(RandomGenerator::GetInstance().GetGenerator());
double decay1costheta = m_step2Distribution.GetRandomCosTheta();
std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
double bke = (m_beamDistributions[0])(gen);
double rxnTheta = std::acos((m_thetaRanges[0])(gen));
double rxnPhi = (m_phiRanges[0])(gen);
double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta();
double decay1Theta = std::acos(decay1costheta);
double decay1Phi = m_phi2Range(RandomGenerator::GetInstance().GetGenerator());
double residEx = (*m_exDist)(RandomGenerator::GetInstance().GetGenerator());
double decay1Phi = m_phiRanges[1](gen);
double residEx = (m_exDistributions[0])(gen);
double decay2Ex = m_exDistributions[1](gen);
m_step1.SetBeamKE(bke);
m_step1.SetPolarRxnAngle(rxnTheta);
@ -99,6 +113,7 @@ namespace Mask {
m_step2.SetPolarRxnAngle(decay1Theta);
m_step2.SetAzimRxnAngle(decay1Phi);
m_step2.SetExcitation(decay2Ex);
m_step1.Calculate();

21
src/Mask/TwoStepSystem.h
View File

@ -9,30 +9,17 @@ namespace Mask {
class TwoStepSystem : public ReactionSystem
{
public:
TwoStepSystem();
TwoStepSystem(const std::vector<int>& z, const std::vector<int>& a);
TwoStepSystem(const std::vector<StepParameters>& params);
~TwoStepSystem();
bool SetNuclei(const std::vector<int>& z, const std::vector<int>& a) override;
void RunSystem() override;
std::vector<Nucleus>* GetNuclei() override;
virtual void SetDecay1Distribution(const std::string& filename) override { m_step2Distribution.ReadDistributionFile(filename); };
virtual void SetReactionThetaType(RxnThetaType type) override { m_step1.SetEjectileThetaType(type); };
int GetDecay1AngularMomentum() { return m_step2Distribution.GetL(); };
double GetDecay1BranchingRatio() { return m_step2Distribution.GetBranchingRatio(); };
virtual void SetLayeredTarget(const LayeredTarget& target) override;
virtual void RunSystem() override;
private:
void LinkTarget() override;
void Init(const std::vector<StepParameters>& params);
void SetSystemEquation() override;
std::uniform_real_distribution<double> m_phi2Range;
Reaction m_step1, m_step2;
AngularDistribution m_step2Distribution;
};
}

20
src/Plotters/RootPlotter.cpp
View File

@ -75,18 +75,18 @@ void RootPlotter::FillData(const Mask::Nucleus& nuc)
std::string angdist_name = sym + modifier +"_angDist";
std::string angdist_title = angdist_name+";cos#right(#theta_{CM}#left);counts";
MyFill(ke_vs_th_name.c_str(), ke_vs_th_title.c_str(), nuc.vec4.Theta()*s_rad2deg, nuc.GetKE(), 2);
MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.GetKE(), 4);
MyFill(th_vs_ph_name.c_str(), th_vs_ph_title.c_str(), nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2);
MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy());
MyFill(angdist_name.c_str(), angdist_title.c_str(),20,-1.0,1.0,std::cos(nuc.thetaCM));
MyFill(ke_vs_th_name, ke_vs_th_title, nuc.vec4.Theta()*s_rad2deg, nuc.GetKE(), 2);
MyFill(ke_vs_ph_name, ke_vs_ph_title, FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.GetKE(), 4);
MyFill(th_vs_ph_name, th_vs_ph_title, nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2);
MyFill(ex_name, ex_title, 260, -1.0, 25, nuc.GetExcitationEnergy());
MyFill(angdist_name, angdist_title, 20, -1.0, 1.0, std::cos(nuc.thetaCM));
if(nuc.isDetected)
{
MyFill(ke_vs_th_name.c_str(), ke_vs_th_title.c_str(), nuc.vec4.Theta()*s_rad2deg, nuc.detectedKE, 2);
MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.detectedKE, 4);
MyFill(th_vs_ph_name.c_str(), th_vs_ph_title.c_str(), nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2);
MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy());
MyFill(angdist_name.c_str(), angdist_title.c_str(),20,-1.0,1.0,std::cos(nuc.thetaCM));
MyFill(ke_vs_th_name, ke_vs_th_title, nuc.vec4.Theta()*s_rad2deg, nuc.detectedKE, 2);
MyFill(ke_vs_ph_name, ke_vs_ph_title, FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.detectedKE, 4);
MyFill(th_vs_ph_name, th_vs_ph_title, nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2);
MyFill(ex_name, ex_title, 260, -1.0, 25, nuc.GetExcitationEnergy());
MyFill(angdist_name, angdist_title, 20, -1.0, 1.0, std::cos(nuc.thetaCM));
}
}