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

Style corrections, modified some plotting tools to handle detected energy effects.

Browse Source
This commit is contained in:
Gordon McCann 2021-12-03 09:28:08 -05:00
parent 2adf7c3da1
commit a19e17deb9
12 changed files with 58 additions and 70 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
include/EnergyLoss.h
View File

@ -28,7 +28,7 @@ namespace Mask {
double GetEnergyLoss(int zp, int ap, double e_initial, double thickness);
double GetReverseEnergyLoss(int zp, int ap, double e_final, double thickness);
double GetRange(double energy);
void SetTargetComponents(std::vector<int>& Zt, std::vector<int>& At, std::vector<int>& Stoich);
void SetTargetComponents(const std::vector<int>& Zt, const std::vector<int>& At, const std::vector<int>& Stoich);
private:
double GetElectronicStoppingPower(double energy);

10
include/MaskApp.h
View File

@ -17,23 +17,19 @@ namespace Mask {
bool LoadConfig(const std::string& filename);
bool SaveConfig(const std::string& filename);
inline int GetNumberOfSamples() const { return m_nsamples; };
inline const std::string GetSystemName() const { return sys == nullptr ? "" : sys->GetSystemEquation(); };
inline const std::string GetSystemName() const { return m_sys == nullptr ? "" : m_sys->GetSystemEquation(); };
inline const std::string GetOutputName() const { return m_outfile_name; };
inline const RxnType GetReactionType() const { return m_rxn_type; };
void Run();
private:
void RunOneStepRxn();
void RunOneStepDecay();
void RunTwoStep();
void RunThreeStep();
ReactionSystem* sys;
ReactionSystem* m_sys;
std::string m_outfile_name;
RxnType m_rxn_type;
uint32_t m_nsamples;
int m_nsamples;
};

4
include/MaskFile.h
View File

@ -49,8 +49,8 @@ namespace Mask {
FileType file_type;
std::string filename;
uint32_t buffer_position;
uint32_t buffer_end;
uint64_t buffer_position;
uint64_t buffer_end;
uint32_t data_size;
RxnType m_rxn_type;
uint32_t buffersize_bytes;

2
include/RootPlotter.h
View File

@ -28,7 +28,7 @@ public:
GenerateGraphs();
return table;
};
void FillData(const Mask::Nucleus& nuc, const std::string& modifier = "");
void FillData(const Mask::Nucleus& nuc, double detKE = 0.0, const std::string& modifier = "");
private:
THashTable* table;

3
include/RxnType.h
View File

@ -36,6 +36,7 @@ namespace Mask {
case RxnType::OneStepRxn: return "OneStepRxn";
case RxnType::TwoStepRxn: return "TwoStepRxn";
case RxnType::ThreeStepRxn: return "ThreeStepRxn";
case RxnType::None : return "None";
}
return "None";
@ -51,8 +52,6 @@ namespace Mask {
return static_cast<uint32_t>(type);
}
}
#endif

10
premake5.lua
View File

@ -1,7 +1,7 @@
workspace "Mask"
configurations {
"Debug",
"Release"
"Release",
"Debug"
}
project "Mask"
@ -44,8 +44,8 @@ project "RootPlot"
}
--User specified path to ROOT CERN libraries--
ROOTIncludepath = "/home/gordon/cern/root-6.22.02/root-install/include"
ROOTLibpath = "/home/gordon/cern/root-6.22.02/root-install/lib"
ROOTIncludepath = "/usr/include/root"
ROOTLibpath = "/usr/lib64/root"
includedirs {
"include"
@ -106,4 +106,4 @@ project "DetectEff"
symbols "On"
filter "configurations:Release"
optimize "On"
optimize "On"

23
src/Detectors/AnasenEfficiency.cpp
View File

@ -206,48 +206,40 @@ DetectorResult AnasenEfficiency::IsRing1(Mask::Nucleus& nucleus) {
DetectorResult observation;
//Mask::Vec3 coords;
double thetaIncident, eloss, e_dep;
double thetaIncident;
for(auto& sx3 : m_Ring1) {
auto result = sx3.GetChannelRatio(nucleus.GetTheta(), nucleus.GetPhi());
if(result.first != -1) {
//coords = sx3.GetHitCoordinates(result.first, result.second);
observation.detectFlag = true;
observation.direction = sx3.GetHitCoordinates(result.first, result.second);
thetaIncident = std::acos(observation.direction.Dot(sx3.GetNormRotated())/observation.direction.GetR());
if(thetaIncident > M_PI/2.0)
thetaIncident = M_PI - thetaIncident;
//e_dep = det_silicon.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
observation.det_name = "R1";
//return std::make_pair(true, e_dep);
return observation;
}
}
//return std::make_pair(false, 0.0);
return observation;
}
DetectorResult AnasenEfficiency::IsRing2(Mask::Nucleus& nucleus) {
DetectorResult observation;
//Mask::Vec3 coords;
double thetaIncident, eloss, e_dep;
double thetaIncident;
for(auto& sx3 : m_Ring2) {
auto result = sx3.GetChannelRatio(nucleus.GetTheta(), nucleus.GetPhi());
if(result.first != -1) {
//coords = sx3.GetHitCoordinates(result.first, result.second);
observation.detectFlag = true;
observation.direction = sx3.GetHitCoordinates(result.first, result.second);
thetaIncident = std::acos(observation.direction.Dot(sx3.GetNormRotated())/observation.direction.GetR());
if(thetaIncident > M_PI/2.0)
thetaIncident = M_PI - thetaIncident;
//e_dep = det_silicon.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
observation.det_name = "R2";
//return std::make_pair(true, e_dep);
return observation;
}
}
@ -258,22 +250,17 @@ DetectorResult AnasenEfficiency::IsRing2(Mask::Nucleus& nucleus) {
DetectorResult AnasenEfficiency::IsQQQ(Mask::Nucleus& nucleus) {
DetectorResult observation;
//Mask::Vec3 coords;
double thetaIncident, eloss, e_dep;
double thetaIncident;
for(auto& qqq : m_forwardQQQs) {
auto result = qqq.GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi());
if(result.first != -1) {
//coords = qqq.GetHitCoordinates(result.first, result.second);
observation.detectFlag = true;
observation.direction = qqq.GetHitCoordinates(result.first, result.second);
thetaIncident = std::acos(observation.direction.Dot(qqq.GetNorm())/observation.direction.GetR());
if(thetaIncident > M_PI/2.0)
thetaIncident = M_PI - thetaIncident;
//e_dep = det_silicon.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
//return std::make_pair(true, e_dep);
observation.det_name = "FQQQ";
return observation;
}
@ -283,22 +270,18 @@ DetectorResult AnasenEfficiency::IsQQQ(Mask::Nucleus& nucleus) {
for(auto& qqq : m_backwardQQQs) {
auto result = qqq.GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi());
if(result.first != -1) {
//coords = qqq.GetHitCoordinates(result.first, result.second);
observation.detectFlag = true;
observation.direction = qqq.GetHitCoordinates(result.first, result.second);
thetaIncident = std::acos(observation.direction.Dot(qqq.GetNorm())/observation.direction.GetR());
if(thetaIncident > M_PI/2.0)
thetaIncident = M_PI - thetaIncident;
//e_dep = det_silicon.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
//return std::make_pair(true, e_dep);
observation.det_name = "BQQQ";
return observation;
}
}
//return std::make_pair(false, 0.0);
return observation;
}

6
src/Detectors/DetectorEfficiency.cpp
View File

@ -25,10 +25,10 @@ int main(int argc, char** argv) {
*/
try {
AnasenEfficiency anasen;
anasen.CalculateEfficiency(inputname, outputname, statsname);
AnasenEfficiency anasen;
anasen.CalculateEfficiency(inputname, outputname, statsname);
//std::cout<<"Running consistency check(1=success): "<<anasen.RunConsistencyCheck()<<std::endl;
//anasen.DrawDetectorSystem("/data1/gwm17/TRIUMF_7Bed/simulation/ANASENGeo_centered_target_targetGap_BackQQQ_test.root");
//anasen.DrawDetectorSystem("/data1/gwm17/TRIUMF_7Bed/simulation/ANASENGeo_centered_target_targetGap_BackQQQ_fixedZ.txt");
} catch(const std::exception& e) {
std::cerr<<"Error: "<<e.what()<<std::endl;
std::cerr<<"Terminating."<<std::endl;

2
src/EnergyLoss.cpp
View File

@ -24,7 +24,7 @@ namespace Mask {
EnergyLoss::~EnergyLoss() {}
/*Targets are defined by their atomic number, total number of nucleons, and their stoichiometry within the target compound*/
void EnergyLoss::SetTargetComponents(std::vector<int>& Zt, std::vector<int>& At, std::vector<int>& Stoich) {
void EnergyLoss::SetTargetComponents(const std::vector<int>& Zt, const std::vector<int>& At, const std::vector<int>& Stoich) {
comp_denom = 0;
ZT = Zt;
AT = At;

42
src/MaskApp.cpp
View File

@ -6,14 +6,14 @@
namespace Mask {
MaskApp::MaskApp() :
sys(nullptr)
m_sys(nullptr)
{
std::cout<<"----------GWM Kinematics Simulation----------"<<std::endl;
}
MaskApp::~MaskApp()
{
if(sys) delete sys;
if(m_sys) delete m_sys;
}
bool MaskApp::LoadConfig(const std::string& filename)
@ -42,7 +42,7 @@ namespace Mask {
{
case RxnType::PureDecay:
{
sys = new DecaySystem();
m_sys = new DecaySystem();
m_rxn_type = RxnType::PureDecay;
for(int i=0; i<2; i++) {
input>>z>>a;
@ -53,7 +53,7 @@ namespace Mask {
}
case RxnType::OneStepRxn:
{
sys = new OneStepSystem();
m_sys = new OneStepSystem();
m_rxn_type = RxnType::OneStepRxn;
for(int i=0; i<3; i++) {
input>>z>>a;
@ -64,7 +64,7 @@ namespace Mask {
}
case RxnType::TwoStepRxn:
{
sys = new TwoStepSystem();
m_sys = new TwoStepSystem();
m_rxn_type = RxnType::TwoStepRxn;
for(int i=0; i<4; i++) {
input>>z>>a;
@ -75,7 +75,7 @@ namespace Mask {
}
case RxnType::ThreeStepRxn:
{
sys = new ThreeStepSystem();
m_sys = new ThreeStepSystem();
m_rxn_type = RxnType::TwoStepRxn;
for(int i=0; i<5; i++) {
input>>z>>a;
@ -87,7 +87,7 @@ namespace Mask {
default:
return false;
}
sys->SetNuclei(zvec, avec);
m_sys->SetNuclei(zvec, avec);
int nlayers;
double thickness;
@ -110,7 +110,7 @@ namespace Mask {
input>>z>>a>>s;
zvec.push_back(z); avec.push_back(a); svec.push_back(s);
}
sys->AddTargetLayer(zvec, avec, svec, thickness);
m_sys->AddTargetLayer(zvec, avec, svec, thickness);
input>>junk;
}
std::cout<<"Reaction equation: "<<GetSystemName()<<std::endl;
@ -122,39 +122,39 @@ namespace Mask {
input>>junk>>m_nsamples;
input>>junk>>par1>>junk>>par2;
sys->SetBeamDistro(par1, par2);
m_sys->SetBeamDistro(par1, par2);
input>>junk>>par1;
switch(m_rxn_type)
{
case RxnType::OneStepRxn :
{
dynamic_cast<OneStepSystem*>(sys)->SetReactionThetaType(par1);
dynamic_cast<OneStepSystem*>(m_sys)->SetReactionThetaType(par1);
break;
}
case RxnType::TwoStepRxn :
{
dynamic_cast<TwoStepSystem*>(sys)->SetReactionThetaType(par1);
dynamic_cast<TwoStepSystem*>(m_sys)->SetReactionThetaType(par1);
break;
}
case RxnType::ThreeStepRxn :
{
dynamic_cast<ThreeStepSystem*>(sys)->SetReactionThetaType(par1);
dynamic_cast<ThreeStepSystem*>(m_sys)->SetReactionThetaType(par1);
break;
}
}
input>>junk>>par1>>junk>>par2;
sys->SetTheta1Range(par1, par2);
m_sys->SetTheta1Range(par1, par2);
input>>junk>>par1>>junk>>par2;
sys->SetPhi1Range(par1, par2);
m_sys->SetPhi1Range(par1, par2);
input>>junk>>par1>>junk>>par2;
sys->SetExcitationDistro(par1, par2);
m_sys->SetExcitationDistro(par1, par2);
input>>junk>>dfile1;
input>>junk>>dfile2;
switch(m_rxn_type)
{
case RxnType::OneStepRxn :
{
DecaySystem* this_sys = dynamic_cast<DecaySystem*>(sys);
DecaySystem* this_sys = dynamic_cast<DecaySystem*>(m_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;
@ -162,7 +162,7 @@ namespace Mask {
}
case RxnType::TwoStepRxn :
{
TwoStepSystem* this_sys = dynamic_cast<TwoStepSystem*>(sys);
TwoStepSystem* this_sys = dynamic_cast<TwoStepSystem*>(m_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;
@ -170,7 +170,7 @@ namespace Mask {
}
case RxnType::ThreeStepRxn :
{
ThreeStepSystem* this_sys = dynamic_cast<ThreeStepSystem*>(sys);
ThreeStepSystem* this_sys = dynamic_cast<ThreeStepSystem*>(m_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;
@ -188,7 +188,7 @@ namespace Mask {
void MaskApp::Run() {
std::cout<<"Running simulation..."<<std::endl;
if(sys == nullptr)
if(m_sys == nullptr)
{
return;
}
@ -210,8 +210,8 @@ namespace Mask {
std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
}
sys->RunSystem();
output.WriteData(sys->GetNuclei());
m_sys->RunSystem();
output.WriteData(m_sys->GetNuclei());
}
output.Close();

2
src/MaskFile.cpp
View File

@ -332,7 +332,7 @@ namespace Mask {
}
}
unsigned int local_end = buffer_position + data_size;
uint64_t local_end = buffer_position + data_size;
if(local_end > buffer_end) {
std::cerr<<"Attempting to read past end of file at MaskFile::ReadData! Returning empty"<<std::endl;
data.eof = true;

22
src/Plotters/ROOT/RootPlotter.cpp
View File

@ -11,7 +11,7 @@ RootPlotter::RootPlotter() :
RootPlotter::~RootPlotter() {}
void RootPlotter::FillData(const Mask::Nucleus& nuc, const std::string& modifier) {
void RootPlotter::FillData(const Mask::Nucleus& nuc, double detKE, const std::string& modifier) {
std::string sym = nuc.GetIsotopicSymbol();
std::string ke_vs_th_name = sym + modifier + "_ke_vs_theta";
std::string ke_vs_th_title = ke_vs_th_name + ";#theta_{lab} (degrees);Kinetic Energy (MeV)";
@ -22,10 +22,20 @@ void RootPlotter::FillData(const Mask::Nucleus& nuc, const std::string& modifier
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.GetTheta()*rad2deg, nuc.GetKE(), 2);
MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), nuc.GetPhi()*rad2deg, nuc.GetKE(), 4);
MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy());
MyFill(angdist_name.c_str(), angdist_title.c_str(),100,-1.0,1.0,std::cos(nuc.GetThetaCM()));
if(detKE == 0.0)
{
MyFill(ke_vs_th_name.c_str(), ke_vs_th_title.c_str(), nuc.GetTheta()*rad2deg, nuc.GetKE(), 2);
MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), nuc.GetPhi()*rad2deg, nuc.GetKE(), 4);
MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy());
MyFill(angdist_name.c_str(), angdist_title.c_str(),100,-1.0,1.0,std::cos(nuc.GetThetaCM()));
}
else
{
MyFill(ke_vs_th_name.c_str(), ke_vs_th_title.c_str(), nuc.GetTheta()*rad2deg, nuc.GetKE(), 2);
MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), nuc.GetPhi()*rad2deg, nuc.GetKE(), 4);
MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy());
MyFill(angdist_name.c_str(), angdist_title.c_str(),100,-1.0,1.0,std::cos(nuc.GetThetaCM()));
}
}
@ -120,7 +130,7 @@ int main(int argc, char** argv) {
nucleus.SetVectorSpherical(data.theta[i], data.phi[i], data.p[i], data.E[i]);
plotter.FillData(nucleus);
if(data.detect_flag[i] == true) {
plotter.FillData(nucleus, "detected");
plotter.FillData(nucleus, data.KE[i], "detected");
}
}
count++;