From 170cc7afb043923fd2dd602db15298c5e2615520 Mon Sep 17 00:00:00 2001 From: Gordon McCann Date: Fri, 3 Sep 2021 17:03:41 -0400 Subject: [PATCH] Major update. Now uses a homeade file type MaskFile to save data. Separated kinematics simulation (Mask), plotting (RootPlot/python), and detector efficiency (DetectEff). Mask itself no longer depends on ROOT, only RootPlot does. --- include/AnasenEfficiency.h | 8 +- include/AngularDistribution.h | 8 +- include/DecaySystem.h | 58 +- include/DetectorEfficiency.h | 36 +- include/Kinematics.h | 89 ++- include/MaskFile.h | 88 +-- include/MassLookup.h | 13 +- include/Nucleus.h | 5 + include/OneStepSystem.h | 46 +- include/Plotter.h | 49 -- include/QQQDetector.h | 2 +- include/ReactionSystem.h | 90 +-- include/RootPlotter.h | 48 ++ include/SabreEfficiency.h | 11 +- include/Stopwatch.h | 30 + include/StripDetector.h | 8 +- include/ThreeStepSystem.h | 74 +-- include/TwoStepSystem.h | 64 +-- input.txt | 6 +- premake5.lua | 63 ++- src/AngularDistribution.cpp | 10 +- src/DecaySystem.cpp | 130 ++--- src/Detectors/AnasenEfficiency.cpp | 398 +++----------- src/Detectors/DetectorEfficiency.cpp | 33 ++ src/Detectors/SabreEfficiency.cpp | 385 ++++--------- src/Detectors/StripDetector.cpp | 4 +- src/EnergyLoss.cpp | 4 +- src/Kinematics.cpp | 786 ++++++++++++--------------- src/MaskFile.cpp | 483 ++++++++++------ src/MassLookup.cpp | 78 ++- src/Nucleus.cpp | 16 +- src/OneStepSystem.cpp | 130 ++--- src/Plotter.cpp | 89 --- src/Plotters/ROOT/RootPlotter.cpp | 137 +++++ src/ReactionSystem.cpp | 38 +- src/Stopwatch.cpp | 31 ++ src/ThreeStepSystem.cpp | 246 +++++---- src/TwoStepSystem.cpp | 192 +++---- src/main.cpp | 47 +- 39 files changed, 1917 insertions(+), 2116 deletions(-) delete mode 100644 include/Plotter.h create mode 100644 include/RootPlotter.h create mode 100644 include/Stopwatch.h create mode 100644 src/Detectors/DetectorEfficiency.cpp delete mode 100644 src/Plotter.cpp create mode 100644 src/Plotters/ROOT/RootPlotter.cpp create mode 100644 src/Stopwatch.cpp diff --git a/include/AnasenEfficiency.h b/include/AnasenEfficiency.h index de2f276..fbe1ec8 100644 --- a/include/AnasenEfficiency.h +++ b/include/AnasenEfficiency.h @@ -2,7 +2,6 @@ #define ANASEN_EFFICIENCY_H #include -#include #include "DetectorEfficiency.h" #include "StripDetector.h" @@ -12,16 +11,11 @@ class AnasenEfficiency : public DetectorEfficiency { public: AnasenEfficiency(); ~AnasenEfficiency(); - void CalculateEfficiency(const std::string& filename) override; + void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) override; void DrawDetectorSystem(const std::string& filename) override; double RunConsistencyCheck() override; private: - void RunDecay(const std::string& filename) override; - void Run2Step(const std::string& filename) override; - void Run3Step(const std::string& filename) override; - void Run1Step(const std::string& filename) override; - bool IsRing1(double theta, double phi); bool IsRing2(double theta, double phi); bool IsQQQ(double theta, double phi); diff --git a/include/AngularDistribution.h b/include/AngularDistribution.h index ef9f3a8..393699f 100644 --- a/include/AngularDistribution.h +++ b/include/AngularDistribution.h @@ -1,9 +1,9 @@ #ifndef ANGULARDISTRIBUTION_H #define ANGULARDISTRIBUTION_H -#include #include #include +#include class AngularDistribution { public: @@ -11,7 +11,7 @@ public: AngularDistribution(const std::string& file); ~AngularDistribution(); void ReadDistributionFile(const std::string& file); - void AttachRandomNumberGenerator(TRandom3* random) { generator = random; }; + void AttachRandomNumberGenerator(std::mt19937* random) { generator = random; }; double GetRandomCosTheta(); int GetL() { return L; }; double GetBranchingRatio() { return branchingRatio; }; @@ -26,7 +26,9 @@ private: } } - TRandom3* generator; //NOT OWNED BY ANGULAR DISTRIBUTION + std::mt19937* generator; //NOT OWNED BY ANGULAR DISTRIBUTION + std::uniform_real_distribution uniform_cosine_dist; + std::uniform_real_distribution uniform_prob_dist; double branchingRatio; int L; diff --git a/include/DecaySystem.h b/include/DecaySystem.h index d435561..8ab0bfd 100644 --- a/include/DecaySystem.h +++ b/include/DecaySystem.h @@ -6,36 +6,36 @@ namespace Mask { -class DecaySystem: public ReactionSystem { -public: - DecaySystem(); - DecaySystem(std::vector& z, std::vector& a); - ~DecaySystem(); - - bool SetNuclei(std::vector& z, std::vector& 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; + class DecaySystem: public ReactionSystem { + public: + DecaySystem(); + DecaySystem(std::vector& z, std::vector& a); + ~DecaySystem(); -}; + bool SetNuclei(std::vector& z, std::vector& a) override; + void RunSystem() override; + void SetRandomGenerator(std::mt19937* 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 \ No newline at end of file diff --git a/include/DetectorEfficiency.h b/include/DetectorEfficiency.h index 7c24e0c..e0a1385 100644 --- a/include/DetectorEfficiency.h +++ b/include/DetectorEfficiency.h @@ -1,52 +1,22 @@ #ifndef DETECTOREFFICIENCY_H #define DETECTOREFFICIENCY_H -#include -#include -#include #include +#include class DetectorEfficiency { public: - DetectorEfficiency() { m_rxn_type = -1; }; + DetectorEfficiency() {}; virtual ~DetectorEfficiency() {}; - inline void SetReactionType(int rxntype) { m_rxn_type = rxntype; }; - virtual void CalculateEfficiency(const std::string& filename) = 0; + virtual void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) = 0; virtual void DrawDetectorSystem(const std::string& filename) = 0; virtual double RunConsistencyCheck() = 0; protected: - void MyFill(THashTable* table, const std::string& name, const std::string& title, int bins, float min, float max, double val) { - TH1F* h = (TH1F*) table->FindObject(name.c_str()); - if(h) { - h->Fill(val); - } else { - h = new TH1F(name.c_str(), title.c_str(), bins, min, max); - h->Fill(val); - table->Add(h); - } - } - void MyFill(THashTable* table, const std::string& name, const std::string& title, int binsx, float minx, float maxx, double valx, int binsy, float miny, float maxy, double valy) { - TH2F* h = (TH2F*) table->FindObject(name.c_str()); - if(h) { - h->Fill(valx, valy); - } else { - h = new TH2F(name.c_str(), title.c_str(), binsx, minx, maxx, binsy, miny, maxy); - h->Fill(valx, valy); - table->Add(h); - } - } inline bool IsDoubleEqual(double x, double y) { return std::fabs(x-y) < epsilon ? true : false; }; - - virtual void Run2Step(const std::string& filename) = 0; - virtual void Run3Step(const std::string& filename) = 0; - virtual void RunDecay(const std::string& filename) = 0; - virtual void Run1Step(const std::string& filename) = 0; - static constexpr double epsilon = 1.0e-6; - int m_rxn_type; }; #endif \ No newline at end of file diff --git a/include/Kinematics.h b/include/Kinematics.h index a015cee..574a59d 100644 --- a/include/Kinematics.h +++ b/include/Kinematics.h @@ -6,66 +6,45 @@ #include "OneStepSystem.h" #include "TwoStepSystem.h" #include "ThreeStepSystem.h" -#include "Plotter.h" -#include -#include +#include namespace Mask { -//For tree -struct NucData { - double KE = -1; - double E = -1; - double Ex = -1; - double p = -1; - double theta = -1; - double theta_cm = -1; - double phi = -1; - int Z = -1; - int A = -1; -}; -class Kinematics { -public: - Kinematics(); - ~Kinematics(); - bool LoadConfig(const char* filename); - bool SaveConfig(const char* filename); - 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().c_str(); }; - inline const char* GetOutputName() { return m_outfile_name.c_str(); }; - inline int GetReactionType() { return m_rxn_type; }; - void Run(); - - enum RxnType { - ONESTEP_RXN, - ONESTEP_DECAY, - TWOSTEP, - THREESTEP + class Kinematics { + public: + Kinematics(); + ~Kinematics(); + bool LoadConfig(const std::string& filename); + bool SaveConfig(const std::string& filename); + inline int GetNumberOfSamples() { return m_nsamples; }; + inline const std::string GetSystemName() { return sys == nullptr ? "" : sys->GetSystemEquation(); }; + inline const std::string GetOutputName() { return m_outfile_name; }; + inline const int GetReactionType() { return m_rxn_type; }; + void Run(); + + enum RxnType { + ONESTEP_DECAY, + ONESTEP_RXN, + TWOSTEP, + THREESTEP + }; + + private: + void RunOneStepRxn(); + void RunOneStepDecay(); + void RunTwoStep(); + void RunThreeStep(); + + ReactionSystem* sys; + + std::string m_outfile_name; + + int m_rxn_type, m_nsamples; + + std::mt19937* global_generator; }; -private: - void RunOneStepRxn(); - void RunOneStepDecay(); - void RunTwoStep(); - void RunThreeStep(); - - ReactionSystem* sys; - Plotter plotman; - - NucData ConvertNucleus(const Nucleus& nuc); - - bool save_tree_flag, do_plotter_flag; - - std::string m_outfile_name; - - int m_rxn_type, m_nsamples; - - TRandom3* global_generator; -}; - -}; +} #endif diff --git a/include/MaskFile.h b/include/MaskFile.h index b584f89..e0a8bae 100644 --- a/include/MaskFile.h +++ b/include/MaskFile.h @@ -9,45 +9,67 @@ namespace Mask { -class MaskFile { -public: - enum class FileType { - read, - write, - append, - none + struct MaskFileHeader { + int rxn_type = -1; + int nsamples = -1; }; - MaskFile(); - MaskFile(const std::string& name, MaskFile::FileType type); - bool Open(const std::string& name, MaskFile::FileType type); - inline bool IsOpen() { return file.is_open(); } - void Close(); + struct MaskFileData { + std::vector E, KE, p, theta, phi; //ordered: target, (if not decay)projectile, ejectile, residual, break1... + std::vector Z, A; + std::vector detect_flag; + bool eof = false; //flag on end of file + }; + + class MaskFile { + public: + enum class FileType { + read, + write, + append, + none + }; - void WriteHeader(std::vector& data); - void WriteData(std::vector& data); - void ReadHeader(); - std::vector ReadData(); - + MaskFile(); + MaskFile(const std::string& name, MaskFile::FileType type); + bool Open(const std::string& name, MaskFile::FileType type); + inline bool IsOpen() { return file.is_open(); } + void Close(); + + void WriteHeader(int rxn_type, int nsamples); + void WriteData(std::vector& data); + void WriteData(MaskFileData& data); + MaskFileHeader ReadHeader(); + MaskFileData ReadData(); + + + private: -private: - void FlushBuffer(); - void FillBuffer(); + FileType file_type; + std::string filename; + unsigned int buffer_position; + unsigned int buffer_end; + unsigned int data_size; + int m_rxn_type; + int buffersize_bytes; + + std::fstream file; + + std::vector data_buffer; + + static constexpr int onestep_rxn_n = 2; + static constexpr int twostep_rxn_n = 4; + static constexpr int threestep_rxn_n = 6; + + static constexpr int buffersize = 10000; //number of events + static constexpr int width = 0; + static constexpr int precision = 3; - FileType file_type; - std::string filename; - unsigned int buffer_position; - int count; - - std::fstream file; - - std::vector> data_buffer; - - static constexpr int buffersize = 1000; - static constexpr int width = 0; - static constexpr int precision = 3; -}; + static constexpr std::size_t double_size = sizeof(double); + static constexpr std::size_t int_size = sizeof(int); + static constexpr std::size_t bool_size = sizeof(bool); + }; }; diff --git a/include/MassLookup.h b/include/MassLookup.h index 404116f..ba7bd52 100644 --- a/include/MassLookup.h +++ b/include/MassLookup.h @@ -7,6 +7,7 @@ of this map (MASS) for use throughout code it is included into. Written by G.W. McCann Aug. 2020 +Converted to true singleton to simplify usage -- Aug. 2021 GWM */ #ifndef MASS_LOOKUP_H #define MASS_LOOKUP_H @@ -20,21 +21,27 @@ Written by G.W. McCann Aug. 2020 class MassLookup { public: - MassLookup(); ~MassLookup(); double FindMass(int Z, int A); std::string FindSymbol(int Z, int A); + static MassLookup* GetInstance() { + if(s_instance == nullptr) { + s_instance = new MassLookup(); + } + return s_instance; + } private: + MassLookup(); std::unordered_map massTable; std::unordered_map elementTable; + static MassLookup* s_instance; + //constants static constexpr double u_to_mev = 931.4940954; static constexpr double electron_mass = 0.000548579909; }; -//static instance for use throught program -static MassLookup MASS; #endif diff --git a/include/Nucleus.h b/include/Nucleus.h index 8a42f52..70cca6f 100644 --- a/include/Nucleus.h +++ b/include/Nucleus.h @@ -28,6 +28,9 @@ public: inline double GetGroundStateMass() const { return m_gs_mass; }; inline const char* GetIsotopicSymbol() const { return m_symbol.c_str(); }; inline double GetThetaCM() const { return m_theta_cm; }; + inline void SetDetected() { m_detectFlag = true; }; + inline void SetNotDetected() { m_detectFlag = false; }; + inline bool IsDetected() { return m_detectFlag; }; inline Nucleus& operator=(const Nucleus& rhs) { SetIsotope(rhs.GetZ(), rhs.GetA()); @@ -52,6 +55,8 @@ private: double m_theta_cm; std::string m_symbol; + bool m_detectFlag; + }; }; diff --git a/include/OneStepSystem.h b/include/OneStepSystem.h index 6b58fc4..06a276b 100644 --- a/include/OneStepSystem.h +++ b/include/OneStepSystem.h @@ -5,29 +5,29 @@ namespace Mask { -class OneStepSystem: public ReactionSystem { -public: - OneStepSystem(); - OneStepSystem(std::vector& z, std::vector& a); - ~OneStepSystem(); + class OneStepSystem: public ReactionSystem { + public: + OneStepSystem(); + OneStepSystem(std::vector& z, std::vector& a); + ~OneStepSystem(); + + bool SetNuclei(std::vector& z, std::vector& 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; + + }; - bool SetNuclei(std::vector& z, std::vector& 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 \ No newline at end of file diff --git a/include/Plotter.h b/include/Plotter.h deleted file mode 100644 index 5e4bb16..0000000 --- a/include/Plotter.h +++ /dev/null @@ -1,49 +0,0 @@ -#ifndef PLOTTER_H -#define PLOTTER_H - -#include -#include -#include -#include -#include -#include "Nucleus.h" - -namespace Mask { - -struct GraphData { - std::string name; - std::string title; - std::vector xvec; - std::vector yvec; - int color; -}; - -class Plotter { -public: - Plotter(); - ~Plotter(); - inline void ClearTable() { table->Clear(); }; - inline THashTable* GetTable() { - GenerateGraphs(); - return table; - }; - void FillData(const Nucleus& nuc, const std::string& modifier = ""); - -private: - THashTable* table; - - void GenerateGraphs(); - void MyFill(const char* name, const char* title, int bins, float min, float max, double val); - void MyFill(const char* name, const char* title, int binsx, float minx, float maxx, int binsy, float miny, float maxy, double valx, double valy); - void MyFill(const char* name, const char* title, double valx, double valy, int color); //TGraph - - std::vector garbage_collection; - std::vector graphs; - - static constexpr double rad2deg = 180.0/M_PI; - -}; - -}; - -#endif \ No newline at end of file diff --git a/include/QQQDetector.h b/include/QQQDetector.h index 3ffee56..9ee6d90 100644 --- a/include/QQQDetector.h +++ b/include/QQQDetector.h @@ -8,8 +8,8 @@ #ifndef QQQDETECTOR_H #define QQQDETECTOR_H -#include #include +#include #include "Vec3.h" #include "Rotation.h" diff --git a/include/ReactionSystem.h b/include/ReactionSystem.h index 4aa885f..7e64138 100644 --- a/include/ReactionSystem.h +++ b/include/ReactionSystem.h @@ -12,45 +12,65 @@ #include "Reaction.h" #include "KinematicsExceptions.h" #include -#include +#include namespace Mask { -class ReactionSystem { -public: - ReactionSystem(); - virtual ~ReactionSystem(); - - virtual bool SetNuclei(std::vector& z, std::vector& a) = 0; - virtual void RunSystem() = 0; - - void AddTargetLayer(std::vector& zt, std::vector& at, std::vector& stoich, double thickness); - - /*Set sampling parameters*/ - virtual void SetRandomGenerator(TRandom3* gen); - inline void SetBeamDistro(double mean, double sigma) { m_beamDist = std::make_pair(mean, sigma); }; - 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 const std::string& GetSystemEquation() const { return m_sys_equation; }; - -protected: - virtual void LinkTarget() = 0; - virtual void SetSystemEquation() = 0; + class ReactionSystem { + public: + ReactionSystem(); + virtual ~ReactionSystem(); - LayeredTarget target; + virtual bool SetNuclei(std::vector& z, std::vector& a) = 0; + virtual void RunSystem() = 0; + + void AddTargetLayer(std::vector& zt, std::vector& at, std::vector& stoich, double thickness); + + /*Set sampling parameters*/ + virtual void SetRandomGenerator(std::mt19937* gen); + inline void SetBeamDistro(double mean, double sigma) { + if(m_beamDist) + delete m_beamDist; + m_beamDist = new std::normal_distribution(mean, sigma); + } + + inline void SetTheta1Range(double min, double max) { + if(m_theta1Range) + delete m_theta1Range; + m_theta1Range = new std::uniform_real_distribution(std::cos(min*deg2rad), std::cos(max*deg2rad)); + } + + inline void SetPhi1Range(double min, double max) { + if(m_phi1Range) + delete m_phi1Range; + m_phi1Range = new std::uniform_real_distribution(min*deg2rad, max*deg2rad); + } + + inline void SetExcitationDistro(double mean, double sigma) { + if(m_exDist) + delete m_exDist; + m_exDist = new std::normal_distribution(mean, sigma); + } + + inline const std::string& GetSystemEquation() const { return m_sys_equation; } + + protected: + virtual void LinkTarget() = 0; + virtual void SetSystemEquation() = 0; + + LayeredTarget target; + + //Sampling information + std::normal_distribution *m_beamDist, *m_exDist; + std::uniform_real_distribution *m_theta1Range, *m_phi1Range; + std::mt19937* generator; //not owned by ReactionSystem + + bool target_set_flag, gen_set_flag; + int rxnLayer; + std::string m_sys_equation; + static constexpr double deg2rad = M_PI/180.0; + }; - //Sampling information - std::pair m_beamDist, m_theta1Range, m_phi1Range, m_exDist; - TRandom3* generator; //not owned by ReactionSystem - - bool target_set_flag, gen_set_flag; - int rxnLayer; - std::string m_sys_equation; - static constexpr double deg2rad = M_PI/180.0; -}; - -}; +} #endif \ No newline at end of file diff --git a/include/RootPlotter.h b/include/RootPlotter.h new file mode 100644 index 0000000..2ea5559 --- /dev/null +++ b/include/RootPlotter.h @@ -0,0 +1,48 @@ +#ifndef ROOTPLOTTER_H +#define ROOTPLOTTER_H + +#include +#include + +#include "Nucleus.h" + +#include +#include +#include +#include + +struct GraphData { + std::string name; + std::string title; + std::vector xvec; + std::vector yvec; + int color; +}; + +class RootPlotter { +public: + RootPlotter(); + ~RootPlotter(); + inline void ClearTable() { table->Clear(); }; + inline THashTable* GetTable() { + GenerateGraphs(); + return table; + }; + void FillData(const Mask::Nucleus& nuc, const std::string& modifier = ""); + +private: + THashTable* table; + + void GenerateGraphs(); + void MyFill(const std::string& name, const std::string& title, int bins, float min, float max, double val); + void MyFill(const std::string& name, const std::string& title, int binsx, float minx, float maxx, int binsy, float miny, float maxy, double valx, double valy); + void MyFill(const std::string& name, const std::string& title, double valx, double valy, int color); //TGraph + + std::vector garbage_collection; + std::vector graphs; + + static constexpr double rad2deg = 180.0/M_PI; + +}; + +#endif \ No newline at end of file diff --git a/include/SabreEfficiency.h b/include/SabreEfficiency.h index 6bee112..b795acf 100644 --- a/include/SabreEfficiency.h +++ b/include/SabreEfficiency.h @@ -5,24 +5,19 @@ #include "SabreDetector.h" #include "Target.h" #include "DeadChannelMap.h" -#include "Kinematics.h" -#include +#include "Nucleus.h" class SabreEfficiency : public DetectorEfficiency { public: SabreEfficiency(); ~SabreEfficiency(); void SetDeadChannelMap(std::string& filename) { dmap.LoadMapfile(filename); }; - void CalculateEfficiency(const std::string& file) override; + void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) override; void DrawDetectorSystem(const std::string& filename) override; double RunConsistencyCheck() override; private: - std::pair IsSabre(Mask::NucData* nucleus); - void Run2Step(const std::string& filename) override; - void Run3Step(const std::string& filename) override; - void RunDecay(const std::string& filename) override; - void Run1Step(const std::string& filename) override; + std::pair IsSabre(Mask::Nucleus& nucleus); std::vector detectors; diff --git a/include/Stopwatch.h b/include/Stopwatch.h new file mode 100644 index 0000000..af217c1 --- /dev/null +++ b/include/Stopwatch.h @@ -0,0 +1,30 @@ +/* + Stopwatch.h + Simple class designed to provide timing info on parts of the process. + Only for use in development. + + Written by G.W. McCann Oct. 2020 +*/ +#ifndef STOPWATCH_H +#define STOPWATCH_H + +#include + +class Stopwatch { + +public: + Stopwatch(); + ~Stopwatch(); + void Start(); + void Stop(); + double GetElapsedSeconds(); + double GetElapsedMilliseconds(); + +private: + using Time = std::chrono::high_resolution_clock::time_point; + using Clock = std::chrono::high_resolution_clock; + + Time start_time, stop_time; +}; + +#endif \ No newline at end of file diff --git a/include/StripDetector.h b/include/StripDetector.h index 2aaef13..7479cf0 100644 --- a/include/StripDetector.h +++ b/include/StripDetector.h @@ -11,9 +11,9 @@ //Back strips from lowest z to highest z -#include #include #include +#include #include "Vec3.h" #include "Rotation.h" @@ -27,7 +27,7 @@ public: inline Mask::Vec3 GetBackStripCoordinates(int stripch, int corner) { return back_strip_coords[stripch][corner]; }; inline Mask::Vec3 GetRotatedFrontStripCoordinates(int stripch, int corner) { return rotated_front_strip_coords[stripch][corner]; }; inline Mask::Vec3 GetRotatedBackStripCoordinates(int stripch, int corner) { return rotated_back_strip_coords[stripch][corner]; }; - inline void SetRandomNumberGenerator(TRandom3* random) { m_random = random; }; + inline void SetRandomNumberGenerator(std::mt19937* random) { m_random = random; }; Mask::Vec3 GetHitCoordinates(int front_stripch, double front_strip_ratio); std::pair GetChannelRatio(double theta, double phi); @@ -51,7 +51,9 @@ private: Mask::ZRotation zRot; - TRandom3* m_random; //Not owned by StripDetector! + std::uniform_real_distribution m_uniform_fraction; + + std::mt19937* m_random; //Not owned by StripDetector! inline bool ValidChannel(int f) { return ((f >= 0 && f < num_strips) ? true : false); }; inline bool ValidRatio(double r) { return ((r >= -1 && r <= 1) ? true : false); }; diff --git a/include/ThreeStepSystem.h b/include/ThreeStepSystem.h index 4d95356..76447c8 100644 --- a/include/ThreeStepSystem.h +++ b/include/ThreeStepSystem.h @@ -6,43 +6,45 @@ namespace Mask { -class ThreeStepSystem : public ReactionSystem { -public: - ThreeStepSystem(); - ThreeStepSystem(std::vector& z, std::vector& a); - ~ThreeStepSystem(); - bool SetNuclei(std::vector& z, std::vector& a) override; - void RunSystem() override; - void SetRandomGenerator(TRandom3* gen) override; + class ThreeStepSystem : public ReactionSystem { + public: + ThreeStepSystem(); + ThreeStepSystem(std::vector& z, std::vector& a); + ~ThreeStepSystem(); + bool SetNuclei(std::vector& z, std::vector& a) override; + void RunSystem() override; + void SetRandomGenerator(std::mt19937* 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(); }; + + inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); }; + inline int GetDecay2AngularMomentum(){ return decay2dist.GetL(); }; + inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); }; + inline double GetDecay2BranchingRatio(){ return decay2dist.GetBranchingRatio(); }; + + protected: + void LinkTarget() override; + void SetSystemEquation() override; - inline void SetDecay1Distribution(const std::string& filename) { decay1dist.ReadDistributionFile(filename); }; - inline void SetDecay2Distribution(const std::string& filename) { decay2dist.ReadDistributionFile(filename); }; + std::uniform_real_distribution m_phi2Range; + + Reaction step1, step2, step3; + + AngularDistribution decay1dist, decay2dist; + + }; - 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(); }; - - inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); }; - inline int GetDecay2AngularMomentum(){ return decay2dist.GetL(); }; - inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); }; - inline double GetDecay2BranchingRatio(){ return decay2dist.GetBranchingRatio(); }; - -protected: - void LinkTarget() override; - void SetSystemEquation() override; - - Reaction step1, step2, step3; - - AngularDistribution decay1dist, decay2dist; - -}; - -}; +} #endif \ No newline at end of file diff --git a/include/TwoStepSystem.h b/include/TwoStepSystem.h index 44c927a..c985069 100644 --- a/include/TwoStepSystem.h +++ b/include/TwoStepSystem.h @@ -6,39 +6,41 @@ namespace Mask { -class TwoStepSystem : public ReactionSystem { -public: - TwoStepSystem(); - TwoStepSystem(std::vector& z, std::vector& a); - ~TwoStepSystem(); - bool SetNuclei(std::vector& z, std::vector& 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() override; - void SetSystemEquation() override; + class TwoStepSystem : public ReactionSystem { + public: + TwoStepSystem(); + TwoStepSystem(std::vector& z, std::vector& a); + ~TwoStepSystem(); + bool SetNuclei(std::vector& z, std::vector& a) override; + void RunSystem() override; - Reaction step1, step2; + void SetRandomGenerator(std::mt19937* 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() override; + void SetSystemEquation() override; - AngularDistribution decay1dist; + std::uniform_real_distribution m_phi2Range; + + Reaction step1, step2; + + AngularDistribution decay1dist; + + }; -}; - -}; +} #endif \ No newline at end of file diff --git a/input.txt b/input.txt index 69ced3c..bc9234c 100644 --- a/input.txt +++ b/input.txt @@ -1,5 +1,5 @@ ----------Data Information---------- -OutputFile: /data1/gwm17/mask_tests/7Bedp_600keV_beam_centered_target_targetgap_BackQQQ_rndmCM_test.root +OutputFile: /data1/gwm17/mask_tests/7Bedp_600keV_beam_centered_target_targetgap_BackQQQ_rndmCM_test.mask SaveTree: yes SavePlots: yes ----------Reaction Information---------- @@ -20,8 +20,8 @@ Z A Stoich 0 ~ ----------Sampling Information---------- -NumberOfSamples: 100000 -BeamMeanEnergy(MeV): 0.6 BeamEnergySigma(MeV): 0.001 +NumberOfSamples: 10000 +BeamMeanEnergy(MeV): 0.6 BeamEnergySigma(MeV): 0.0 EjectileThetaType(0=Lab,1=CM): 1 EjectileThetaMin(deg): 0.0 EjectileThetaMax(deg): 180.0 EjectilePhiMin(deg): 0.0 EjectilePhiMax(deg): 360.0 diff --git a/premake5.lua b/premake5.lua index 4cd98ee..1eb4b76 100644 --- a/premake5.lua +++ b/premake5.lua @@ -12,14 +12,35 @@ project "Mask" cppdialect "C++11" files { - "src/**.cpp", - "include/**.h", - "src/**.cxx" + "src/*.cpp", + "include/*.h" } - prebuildcommands { - "rootcint -f src/kinematics_dict.cxx include/Kinematics.h include/LinkDef_Kinematics.h", - "{MOVE} src/kinematics_dict_rdict.pcm bin/" + includedirs { + "include" + } + + filter "configurations:Debug" + symbols "On" + + filter "configurations:Release" + optimize "On" + +project "RootPlot" + kind "ConsoleApp" + language "C++" + targetdir "bin" + objdir "objs" + cppdialect "c++11" + + files { + "src/Plotters/ROOT/RootPlotter.cpp", + "src/MaskFile.cpp", + "src/Nucleus.cpp", + "src/Vec4.cpp", + "src/Vec3.cpp", + "src/MassLookup.cpp", + "include/*.h" } filter "system:windows" @@ -54,6 +75,36 @@ project "Mask" "`root-config --glibs`" } + filter "configurations:Debug" + symbols "On" + + filter "configurations:Release" + optimize "On" + +project "DetectEff" + kind "ConsoleApp" + language "C++" + targetdir "bin" + objdir "objs" + cppdialect "c++11" + + files { + "src/Detectors/*.cpp", + "src/MaskFile.cpp", + "src/Nucleus.cpp", + "src/Vec4.cpp", + "src/Vec3.cpp", + "src/MassLookup.cpp", + "src/Rotation.cpp", + "src/Target.cpp", + "src/EnergyLoss.cpp", + "include/*.h" + } + + includedirs { + "include" + } + filter "configurations:Debug" symbols "On" diff --git a/src/AngularDistribution.cpp b/src/AngularDistribution.cpp index 453b862..1065c71 100644 --- a/src/AngularDistribution.cpp +++ b/src/AngularDistribution.cpp @@ -5,7 +5,7 @@ #include "LegendrePoly.h" AngularDistribution::AngularDistribution() : - generator(nullptr), branchingRatio(1.0), L(0), isoFlag(true) + generator(nullptr), uniform_cosine_dist(-1.0, 1.0), uniform_prob_dist(0.0, 1.0), branchingRatio(1.0), L(0), isoFlag(true) { } @@ -84,18 +84,18 @@ double AngularDistribution::GetRandomCosTheta() { return 0.0; } - if(isoFlag) return generator->Uniform(-1.0, 1.0); + if(isoFlag) return uniform_cosine_dist(*generator); double test, probability; double costheta; - test = generator->Uniform(0.0, 1.0); + test = uniform_prob_dist(*generator); if(test > branchingRatio) return -10; do { probability = 0.0; - costheta = generator->Uniform(-1.0, 1.0); - test = generator->Uniform(0.0, 1.0); + costheta = uniform_cosine_dist(*generator); + test = uniform_prob_dist(*generator); for(unsigned int i=0; i probability); diff --git a/src/DecaySystem.cpp b/src/DecaySystem.cpp index b18165a..eaa0ba0 100644 --- a/src/DecaySystem.cpp +++ b/src/DecaySystem.cpp @@ -2,74 +2,74 @@ namespace Mask { -DecaySystem::DecaySystem() : - ReactionSystem() -{ -} - -DecaySystem::DecaySystem(std::vector& z, std::vector& 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& z, std::vector& a) { - if(z.size() != a.size() || z.size() != 2) { - return false; + DecaySystem::DecaySystem() : + ReactionSystem() + { } - - 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(); + DecaySystem::DecaySystem(std::vector& z, std::vector& a) : + ReactionSystem() + { + SetNuclei(z, a); } - - 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; + + DecaySystem::~DecaySystem() {} + + void DecaySystem::SetRandomGenerator(std::mt19937* gen) { + generator = gen; + decay1dist.AttachRandomNumberGenerator(gen); + gen_set_flag = true; + } + + bool DecaySystem::SetNuclei(std::vector& z, std::vector& 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 = (*m_phi1Range)(*generator); + step1.SetPolarRxnAngle(rxnTheta); + step1.SetAzimRxnAngle(rxnPhi); + + step1.TurnOnResidualEloss(); + step1.Calculate(); + } else { + return; + } } -} } \ No newline at end of file diff --git a/src/Detectors/AnasenEfficiency.cpp b/src/Detectors/AnasenEfficiency.cpp index 2c43f19..117644a 100644 --- a/src/Detectors/AnasenEfficiency.cpp +++ b/src/Detectors/AnasenEfficiency.cpp @@ -1,11 +1,8 @@ #include "AnasenEfficiency.h" #include "Kinematics.h" -#include -#include -#include -#include -#include -#include +#include "MaskFile.h" +#include +#include AnasenEfficiency::AnasenEfficiency() : DetectorEfficiency() @@ -24,7 +21,7 @@ AnasenEfficiency::~AnasenEfficiency() {} void AnasenEfficiency::DrawDetectorSystem(const std::string& filename) { - TFile* file = TFile::Open(filename.c_str(), "RECREATE"); + std::ofstream output(filename); std::vector x, y, z; std::vector cx, cy, cz; @@ -105,24 +102,18 @@ void AnasenEfficiency::DrawDetectorSystem(const std::string& filename) { } } } - TGraph2D* graph = new TGraph2D(x.size(), &(x[0]), &(y[0]), &(z[0])); - graph->SetName("CornerGraph"); - graph->SetMarkerStyle(2); - graph->SetLineColor(1); - TGraph2D* graph2 = new TGraph2D(cx.size(), &(cx[0]), &(cy[0]), &(cz[0])); - graph2->SetName("CenterGraph"); - graph2->SetMarkerStyle(2); - graph2->SetMarkerColor(4); + output<<"ANASEN Geometry File -- Coordinates for Detectors"<SetName("ANASEN Detector"); - graph->Draw("A|P"); - graph2->Draw("same|P"); - canvas->Write(); - graph->Write(); - graph2->Write(); - file->Close(); + output.close(); } double AnasenEfficiency::RunConsistencyCheck() { @@ -248,305 +239,86 @@ bool AnasenEfficiency::IsQQQ(double theta, double phi) { return false; } -void AnasenEfficiency::CalculateEfficiency(const std::string& file) { +void AnasenEfficiency::CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) { std::cout<<"----------ANASEN Efficiency Calculation----------"< counts; + switch(header.rxn_type) { + case 0: + counts.resize(3, 0); break; - } - case Mask::Kinematics::ONESTEP_RXN: - { - Run1Step(file); - } - case Mask::Kinematics::TWOSTEP: - { - Run2Step(file); + case 1: + counts.resize(4, 0); break; - } - case Mask::Kinematics::THREESTEP: - { - Run3Step(file); + case 2: + counts.resize(6, 0); break; + case 3: + counts.resize(8, 0); + break; + default: + { + std::cerr<<"Bad reaction type at AnasenEfficiency::CalculateEfficiency (given value: "<= threshold && (IsRing1(data.theta[i], data.phi[i]) || IsRing2(data.theta[i], data.phi[i]) || IsQQQ(data.theta[i], data.phi[i]))) { + data.detect_flag[i] = true; + counts[i]++; + } else if(data.detect_flag[i] == true) { + data.detect_flag[i] = false; + } + } + + output.WriteData(data); + } + + input.Close(); + output.Close(); + + stats<Get("DataTree"); - - THashTable* table = new THashTable(); - Mask::NucData* eject = new Mask::NucData(); - Mask::NucData* resid = new Mask::NucData(); - - tree->SetBranchAddress("ejectile", &eject); - tree->SetBranchAddress("residual", &resid); - - double nevents = tree->GetEntries(); - - //Progress tracking - int percent5 = nevents*0.05; - int count = 0; - int npercent = 0; - - int detected_eject=0; - int detected_resid=0; - for(int i=0; iGetEntries(); i++) { - if(++count == percent5) {//Show progress every 5% - npercent++; - count = 0; - std::cout<<"\rPercent completed: "<GetEntry(i); - - if(eject->KE >= threshold && (IsRing1(eject->theta, eject->phi) || IsRing2(eject->theta, eject->phi) || IsQQQ(eject->theta, eject->phi))) { - detected_eject++; - MyFill(table, "detected_eject_theta_ke","detected ejectile theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,eject->theta/deg2rad,300,0.0,30.0,eject->KE); - } - - if(resid->KE > threshold && (IsRing1(resid->theta, resid->phi) || IsRing2(resid->theta, resid->phi) || IsQQQ(resid->theta, resid->phi))) { - detected_resid++; - MyFill(table, "detected_resid_theta_ke","detected residual theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,resid->theta/deg2rad,300,0.0,30.0,resid->KE); - } - - } - - double ejecteff = ((double) detected_eject)/nevents; - double resideff = ((double) detected_resid)/nevents; - TParameter eject_eff("Light Breakup Efficiency", ejecteff); - TParameter resid_eff("Heavy Breakup Efficiency", resideff); - - input->cd(); - eject_eff.Write(); - resid_eff.Write(); - table->Write(); - input->Close(); -} - -void AnasenEfficiency::Run1Step(const std::string& filename) { - std::cout<<"SabreEfficiency::Run1Step Not Implemented!"<Get("DataTree"); - - Mask::NucData* eject = new Mask::NucData(); - Mask::NucData* break1 = new Mask::NucData(); - Mask::NucData* break2 = new Mask::NucData(); - - THashTable* table = new THashTable(); - - tree->SetBranchAddress("ejectile", &eject); - tree->SetBranchAddress("breakup1", &break1); - tree->SetBranchAddress("breakup2", &break2); - - double nevents = tree->GetEntries(); - - //Progress tracking - int percent5 = nevents*0.05; - int count = 0; - int npercent = 0; - - double costheta_cm =0; - - bool break1_flag, break2_flag, eject_flag; - int b1_count=0, b2_count=0, eject_count=0; - int b1b2_count=0, b1e_count=0, b2e_count=0, b1b2e_count=0; - for(int i=0; iGetEntries(); i++) { - if(++count == percent5) {//Show progress every 5% - npercent++; - count = 0; - std::cout<<"\rPercent completed: "<GetEntry(i); - - if(eject->KE >= threshold && (IsRing1(eject->theta, eject->phi) || IsRing2(eject->theta, eject->phi) || IsQQQ(eject->theta, eject->phi))) { - eject_count++; - eject_flag = true; - MyFill(table, "detected_eject_theta_ke","detected ejectile theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,eject->theta/deg2rad,300,0.0,30.0,eject->KE); - } - - if(break1->KE >= threshold && (IsRing1(break1->theta, break1->phi) || IsRing2(break1->theta, break1->phi) || IsQQQ(break1->theta, break1->phi))) { - b1_count++; - break1_flag = true; - costheta_cm = std::cos(break1->theta_cm); - MyFill(table,"detected_break1_cm_theta","detected breakup1 CM theta;cos(#theta_{CM})",20,-1.0,1.0,costheta_cm); - MyFill(table, "detected_break1_theta_ke","detected breakup1 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break1->theta/deg2rad,300,0.0,30.0,break1->KE); - } - - if(break2->KE >= threshold && (IsRing1(break2->theta, break2->phi) || IsRing2(break2->theta, break2->phi) || IsQQQ(break2->theta, break2->phi))) { - b2_count++; - break2_flag = true; - MyFill(table, "detected_break2_theta_ke","detected breakup2 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break2->theta/deg2rad,300,0.0,30.0,break2->KE); - } - - if(break1_flag && break2_flag && eject_flag) { - b1b2e_count++; - b1b2_count++; - b1e_count++; - b2e_count++; - MyFill(table,"b1b2_theta_theta","b1b2_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break1->theta/deg2rad,180,0.0,180.0,break2->theta/deg2rad); - MyFill(table,"b1b2_ke_ke","b1b2_ke_ke;KE_{3};KE_{4}",300,0.0,30.0,break1->KE,300,0.0,30.0,break2->KE); - MyFill(table,"b1b2_phi_phi","b1b2_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break1->phi/deg2rad,360,0.0,360.0,break2->phi/deg2rad); - } else if(break1_flag && eject_flag) { - b1e_count++; - } else if(break2_flag && eject_flag) { - b2e_count++; - } else if(break1_flag && break2_flag) { - b1b2_count++; - MyFill(table,"b1b2_theta_theta","b1b2_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break1->theta/deg2rad,180,0.0,180.0,break2->theta/deg2rad); - MyFill(table,"b1b2_ke_ke","b1b2_ke_ke;KE_{3};KE_{4}",300,0.0,30.0,break1->KE,300,0.0,30.0,break2->KE); - MyFill(table,"b1b2_phi_phi","b1b2_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break1->phi/deg2rad,360,0.0,360.0,break2->phi/deg2rad); - } - } - - double eeff = ((double) eject_count)/nevents; - double b1eff = ((double) b1_count)/nevents; - double b2eff = ((double) b2_count)/nevents; - double b1b2eff = b1b2_count/nevents; - double b1eeff = b1e_count/nevents; - double b2eeff = b2e_count/nevents; - double b1b2eeff = b1b2e_count/nevents; - TParameter eject_eff("Ejectile Efficiency", eeff); - TParameter break1_eff("Breakup 1 Efficiency", b1eff); - TParameter break2_eff("Breakup 2 Efficiency", b2eff); - TParameter break1break2_eff("Breakup1 Coincident with Breakup2", b1b2eff); - TParameter break1eject_eff("Breakup1 Coincident with Ejectile", b1eeff); - TParameter break2eject_eff("Breakup2 Coincident with Ejectile", b2eeff); - TParameter break1break2eject_eff("All Three Particles Coincident", b1b2eeff); - - - input->cd(); - table->Write(); - eject_eff.Write(); - break1_eff.Write(); - break2_eff.Write(); - break1break2_eff.Write(); - break1eject_eff.Write(); - break2eject_eff.Write(); - break1break2eject_eff.Write(); - - input->Close(); -} - -void AnasenEfficiency::Run3Step(const std::string& filename) { - TFile* input = TFile::Open(filename.c_str(), "UPDATE"); - TTree* tree = (TTree*) input->Get("DataTree"); - - THashTable* table = new THashTable(); - - Mask::NucData* break1 = new Mask::NucData(); - Mask::NucData* break3 = new Mask::NucData(); - Mask::NucData* break4 = new Mask::NucData(); - - - tree->SetBranchAddress("breakup1", &break1); - tree->SetBranchAddress("breakup3", &break3); - tree->SetBranchAddress("breakup4", &break4); - - double nevents = tree->GetEntries(); - - //Progress tracking - int percent5 = nevents*0.05; - int count = 0; - int npercent = 0; - - bool break1_flag, break3_flag, break4_flag; - int b1_count=0, b3_count=0, b4_count=0; - int b1b3_count=0, b1b4_count=0, b3b4_count=0, b1b3b4_count=0; - - for(int i=0; iGetEntries(); i++) { - if(++count == percent5) {//Show progress every 5% - npercent++; - count = 0; - std::cout<<"\rPercent completed: "<GetEntry(i); - - if(break1->KE >= threshold && (IsRing1(break1->theta, break1->phi) || IsRing2(break1->theta, break1->phi) || IsQQQ(break1->theta, break1->phi))) { - b1_count++; - break1_flag = true; - MyFill(table, "detected_break1_theta_ke","detected breakup1 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break1->theta/deg2rad,300,0.0,30.0,break1->KE); - } - - - if(break3->KE >= threshold && (IsRing1(break3->theta, break3->phi) || IsRing2(break3->theta, break3->phi) || IsQQQ(break3->theta, break3->phi))) { - b3_count++; - break3_flag = true; - MyFill(table, "detected_break3_theta_ke","detected breakup3 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break3->theta/deg2rad,300,0.0,30.0,break3->KE); - } - - if(break4->KE >= threshold && (IsRing1(break4->theta, break4->phi) || IsRing2(break4->theta, break4->phi) || IsQQQ(break4->theta, break4->phi))) { - b4_count++; - break4_flag = true; - MyFill(table, "detected_break4_theta_ke","detected breakup4 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break4->theta/deg2rad,300,0.0,30.0,break4->KE); - } - - if(break1_flag && break3_flag && break4_flag) { - b1b3b4_count++; - b1b3_count++; - b1b4_count++; - b3b4_count++; - MyFill(table,"b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break3->theta/deg2rad,180,0.0,180.0,break4->theta/deg2rad); - MyFill(table,"b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",300,0.0,30.0,break3->KE,300,0.0,30.0,break4->KE); - MyFill(table,"b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break3->phi/deg2rad,360,0.0,360.0,break4->phi/deg2rad); - } else if(break1_flag && break3_flag) { - b1b3_count++; - } else if(break1_flag && break4_flag) { - b1b4_count++; - } else if(break3_flag && break4_flag) { - b3b4_count++; - MyFill(table,"b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break3->theta/deg2rad,180,0.0,180.0,break4->theta/deg2rad); - MyFill(table,"b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",300,0.0,30.0,break3->KE,300,0.0,30.0,break4->KE); - MyFill(table,"b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break3->phi/deg2rad,360,0.0,360.0,break4->phi/deg2rad); - } - } - - double b1eff = ((double) b1_count)/nevents; - double b3eff = ((double) b3_count)/nevents; - double b4eff = ((double) b4_count)/nevents; - double b1b3eff = b1b3_count/nevents; - double b1b4eff = b1b4_count/nevents; - double b3b4eff = b3b4_count/nevents; - double b1b3b4eff = b1b3b4_count/nevents; - TParameter break1_eff("Breakup1 Efficiency", b1eff); - TParameter break3_eff("Breakup3 Efficiency", b3eff); - TParameter break4_eff("Breakup4 Efficiency", b4eff); - TParameter b1b3_eff("Break1 Coincident with Break3", b1b3eff); - TParameter b1b4_eff("Break1 Coincident with Break4", b1b4eff); - TParameter b3b4_eff("Break3 Coincident with Break4", b3b4eff); - TParameter b1b3b4_eff("All Breakups Coincident", b1b3b4eff); - - input->cd(); - break1_eff.Write(); - break3_eff.Write(); - break4_eff.Write(); - b1b3_eff.Write(); - b1b4_eff.Write(); - b3b4_eff.Write(); - b1b3b4_eff.Write(); - table->Write(); - input->Close(); } \ No newline at end of file diff --git a/src/Detectors/DetectorEfficiency.cpp b/src/Detectors/DetectorEfficiency.cpp new file mode 100644 index 0000000..5d84933 --- /dev/null +++ b/src/Detectors/DetectorEfficiency.cpp @@ -0,0 +1,33 @@ +#include "SabreEfficiency.h" +#include "AnasenEfficiency.h" +#include +#include + +int main(int argc, char** argv) { + + if(argc != 4) { + std::cerr<<"Incorrect number of commandline arguments! Returning."< #include -#include -#include -#include -#include -#include -#include +#include + SabreEfficiency::SabreEfficiency() : DetectorEfficiency(), deadlayer(DEADLAYER_THIN), sabre_eloss(SABRE_THICKNESS) @@ -28,9 +25,9 @@ SabreEfficiency::SabreEfficiency() : SabreEfficiency::~SabreEfficiency() {} -void SabreEfficiency::CalculateEfficiency(const std::string& file) { +void SabreEfficiency::CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) { std::cout<<"----------SABRE Efficiency Calculation----------"< counts; + switch(header.rxn_type) { + case 0: + counts.resize(3, 0); break; - } - case Mask::Kinematics::ONESTEP_RXN: - { - Run1Step(file); + case 1: + counts.resize(4, 0); break; - } - case Mask::Kinematics::TWOSTEP: - { - Run2Step(file); + case 2: + counts.resize(6, 0); break; - } - case Mask::Kinematics::THREESTEP: - { - Run3Step(file); + case 3: + counts.resize(8, 0); break; + default: + { + std::cerr<<"Bad reaction type at AnasenEfficiency::CalculateEfficiency (given value: "< ringxs, ringys, ringzs; std::vector wedgexs, wedgeys, wedgezs; @@ -95,27 +149,16 @@ void SabreEfficiency::DrawDetectorSystem(const std::string& filename) { } } - TGraph2D* gr = new TGraph2D(ringxs.size(), &(ringxs[0]), &(ringys[0]), &(ringzs[0])); - gr->SetName("ring_detector_edges"); - gr->SetTitle("SABRE Detector; x(m); y(m); z(m)"); - gr->SetMarkerStyle(1); + output<<"SABRE Geometry File -- Coordinates for Detectors"<SetName("wedge_detector_edges"); - gw->SetTitle("SABRE Detector Wedges; x(m); y(m); z(m)"); - gw->SetMarkerStyle(1); - - TCanvas* canvas = new TCanvas(); - canvas->SetName("detector_system"); - canvas->cd(); - gr->Draw("AP"); - gw->Draw("same P"); - - canvas->Write(); - gr->Write(); - gw->Write(); - - output->Close(); + output.close(); } double SabreEfficiency::RunConsistencyCheck() { @@ -141,263 +184,25 @@ double SabreEfficiency::RunConsistencyCheck() { } /*Returns if detected, as well as total energy deposited in SABRE*/ -std::pair SabreEfficiency::IsSabre(Mask::NucData* nucleus) { - if(nucleus->KE <= ENERGY_THRESHOLD) { +std::pair SabreEfficiency::IsSabre(Mask::Nucleus& nucleus) { + if(nucleus.GetKE() <= ENERGY_THRESHOLD) { return std::make_pair(false, 0.0); } Mask::Vec3 coords; double thetaIncident, eloss, e_deposited; for(int i=0; i<5; i++) { - auto chan = detectors[i].GetTrajectoryRingWedge(nucleus->theta, nucleus->phi); + auto chan = detectors[i].GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi()); if(chan.first != -1 && chan.second != -1) { if(dmap.IsDead(i, chan.first, 0) || dmap.IsDead(i, chan.second, 1)) break; //dead channel check - coords = detectors[i].GetTrajectoryCoordinates(nucleus->theta, nucleus->phi); + coords = detectors[i].GetTrajectoryCoordinates(nucleus.GetTheta(), nucleus.GetPhi()); thetaIncident = std::acos(coords.Dot(detectors[i].GetNormTilted())/(coords.GetR())); - eloss = deadlayer.getEnergyLossTotal(nucleus->Z, nucleus->A, nucleus->KE, M_PI - thetaIncident); - if((nucleus->KE - eloss) <= ENERGY_THRESHOLD) break; //deadlayer check - e_deposited = sabre_eloss.getEnergyLossTotal(nucleus->Z, nucleus->A, nucleus->KE - eloss, M_PI - thetaIncident); + eloss = deadlayer.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), M_PI - thetaIncident); + if((nucleus.GetKE() - eloss) <= ENERGY_THRESHOLD) break; //deadlayer check + e_deposited = sabre_eloss.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE() - eloss, M_PI - thetaIncident); return std::make_pair(true, e_deposited); } } return std::make_pair(false,0.0); } - -void SabreEfficiency::RunDecay(const std::string& filename) { - TFile* input = TFile::Open(filename.c_str(), "UPDATE"); - TTree* tree = (TTree*) input->Get("DataTree"); - - THashTable* table = new THashTable(); - - Mask::NucData* eject = new Mask::NucData(); - Mask::NucData* resid = new Mask::NucData(); - - - tree->SetBranchAddress("ejectile", &eject); - tree->SetBranchAddress("residual", &resid); - - double nevents = tree->GetEntries(); - - //Progress tracking - int percent5 = nevents*0.05; - int count = 0; - int npercent = 0; - - std::pair eject_result, resid_result; - int detected_eject=0, detected_resid=0; - for(int i=0; iGetEntries(); i++) { - if(++count == percent5) {//Show progress every 5% - npercent++; - count = 0; - std::cout<<"\rPercent completed: "<GetEntry(i); - - eject_result = IsSabre(eject); - resid_result = IsSabre(resid); - - if(eject_result.first) { - detected_eject++; - MyFill(table, "detected_eject_theta_ke","detected ejectile theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,eject->theta/DEG2RAD,300,0.0,30.0,eject->KE); - MyFill(table, "detected_eject_theta_edep","detected ejectile theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,eject->theta/DEG2RAD,300,0.0,30.0,eject_result.second); - } - - if(resid_result.first) { - detected_resid++; - MyFill(table, "detected_resid_theta_ke","detected residual theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,resid->theta/DEG2RAD,300,0.0,30.0,resid->KE); - MyFill(table, "detected_resid_theta_edep","detected residual theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,resid->theta/DEG2RAD,300,0.0,30.0,resid_result.second); - } - } - - double ejecteff = ((double) detected_eject)/nevents; - double resideff = ((double) detected_resid)/nevents; - TParameter eject_eff("Ejectile Efficiency", ejecteff); - TParameter resid_eff("Residual Efficiency", resideff); - - input->cd(); - eject_eff.Write(); - resid_eff.Write(); - table->Write(); - input->Close(); - -} - -void SabreEfficiency::Run1Step(const std::string& filename) { - std::cout<<"SabreEfficiency::Run1Step Not Implemented!"<Get("DataTree"); - - Mask::NucData* break1 = new Mask::NucData(); - Mask::NucData* break2 = new Mask::NucData(); - - THashTable* table = new THashTable(); - - tree->SetBranchAddress("breakup1", &break1); - tree->SetBranchAddress("breakup2", &break2); - - double nevents = tree->GetEntries(); - - //Progress tracking - int percent5 = nevents*0.05; - int count = 0; - int npercent = 0; - - double costheta_cm =0; - - std::pair break1_result, break2_result; - int detected_b1=0, detected_b2=0; - - for(int i=0; iGetEntries(); i++) { - if(++count == percent5) {//Show progress every 5% - npercent++; - count = 0; - std::cout<<"\rPercent completed: "<GetEntry(i); - - break1_result = IsSabre(break1); - break2_result = IsSabre(break2); - - if(break1_result.first) { - detected_b1++; - costheta_cm = std::cos(break1->theta_cm); - MyFill(table,"detected_break1_cm_theta","cos(#theta_{CM})",20,-1.0,1.0,costheta_cm); - MyFill(table, "detected_break1_theta_ke","detected break1 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break1->theta/DEG2RAD,300,0.0,30.0,break1->KE); - MyFill(table, "detected_break1_theta_edep","detected break1 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break1->theta/DEG2RAD,300,0.0,30.0,break1_result.second); - } - - if(break2_result.first) { - detected_b2++; - MyFill(table, "detected_break2_theta_ke","detected break2 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break2->theta/DEG2RAD,300,0.0,30.0,break2->KE); - MyFill(table, "detected_break2_theta_edep","detected break2 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break2->theta/DEG2RAD,300,0.0,30.0,break2_result.second); - - } - } - - double b1eff = ((double) detected_b1)/nevents; - double b2eff = ((double) detected_b2)/nevents; - TParameter break1_eff("Breakup1 Efficiency", b1eff); - TParameter break2_eff("Breakup2 Efficiency", b2eff); - - input->cd(); - table->Write(); - break1_eff.Write(); - break2_eff.Write(); - input->Close(); - -} - -void SabreEfficiency::Run3Step(const std::string& filename) { - TFile* input = TFile::Open(filename.c_str(), "UPDATE"); - TTree* tree = (TTree*) input->Get("DataTree"); - - THashTable* table = new THashTable(); - - Mask::NucData* break1 = new Mask::NucData(); - Mask::NucData* break3 = new Mask::NucData(); - Mask::NucData* break4 = new Mask::NucData(); - - - tree->SetBranchAddress("breakup1", &break1); - tree->SetBranchAddress("breakup3", &break3); - tree->SetBranchAddress("breakup4", &break4); - - double nevents = tree->GetEntries(); - - //Progress tracking - int percent5 = nevents*0.05; - int count = 0; - int npercent = 0; - - std::pair break1_result, break3_result, break4_result; - int detected_b1=0, detected_b3=0, detected_b4=0; - int b1b3_count=0, b1b4_count=0, b3b4_count=0, b1b3b4_count=0; - - - for(int i=0; iGetEntries(); i++) { - if(++count == percent5) {//Show progress every 5% - npercent++; - count = 0; - std::cout<<"\rPercent completed: "<GetEntry(i); - - if(break1_result.first) { - detected_b1++; - MyFill(table, "detected_break1_theta_ke","detected break1 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break1->theta/DEG2RAD,300,0.0,30.0,break1->KE); - MyFill(table, "detected_break1_theta_edep","detected break1 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break1->theta/DEG2RAD,300,0.0,30.0,break1_result.second); - } - - if(break3_result.first) { - detected_b3++; - MyFill(table, "detected_break3_theta_ke","detected break3 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break3->theta/DEG2RAD,300,0.0,30.0,break3->KE); - MyFill(table, "detected_break3_theta_edep","detected break3 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break3->theta/DEG2RAD,300,0.0,30.0,break3_result.second); - } - - if(break4_result.first) { - detected_b1++; - MyFill(table, "detected_break4_theta_ke","detected break4 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break4->theta/DEG2RAD,300,0.0,30.0,break4->KE); - MyFill(table, "detected_break4_theta_edep","detected break4 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break4->theta/DEG2RAD,300,0.0,30.0,break4_result.second); - } - - if(break1_result.first && break3_result.first && break4_result.first) { - b1b3b4_count++; - b1b3_count++; - b1b4_count++; - b3b4_count++; - MyFill(table,"b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break3->theta/DEG2RAD,180,0,180,break4->theta/DEG2RAD); - MyFill(table,"b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",150,0.0,15.0,break3->KE,150,0.0,15.0,break4->KE); - MyFill(table,"b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break3->phi/DEG2RAD,360,0.0,360.0,break4->phi/DEG2RAD); - } else if(break1_result.first && break3_result.first) { - b1b3_count++; - } else if(break1_result.first && break4_result.first) { - b1b4_count++; - } else if(break3_result.first && break4_result.first) { - b3b4_count++; - MyFill(table,"b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break3->theta/DEG2RAD,180,0,180,break4->theta/DEG2RAD); - MyFill(table,"b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",150,0.0,15.0,break3->KE,150,0.0,15.0,break4->KE); - MyFill(table,"b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break3->phi/DEG2RAD,360,0.0,360.0,break4->phi/DEG2RAD); - } - } - - double b1eff = ((double) detected_b1)/nevents; - double b3eff = ((double) detected_b3)/nevents; - double b4eff = ((double) detected_b4)/nevents; - double b1b3eff = b1b3_count/nevents; - double b1b4eff = b1b4_count/nevents; - double b3b4eff = b3b4_count/nevents; - double b1b3b4eff = b1b3b4_count/nevents; - TParameter break1_eff("Breakup1 Efficiency", b1eff); - TParameter break3_eff("Breakup3 Efficiency", b3eff); - TParameter break4_eff("Breakup4 Efficiency", b4eff); - TParameter b1b3_eff("Break1 Coincident with Break3", b1b3eff); - TParameter b1b4_eff("Break1 Coincident with Break4", b1b4eff); - TParameter b3b4_eff("Break3 Coincident with Break4", b3b4eff); - TParameter b1b3b4_eff("All Breakups Coincident", b1b3b4eff); - - input->cd(); - break1_eff.Write(); - break3_eff.Write(); - break4_eff.Write(); - b1b3_eff.Write(); - b1b4_eff.Write(); - b3b4_eff.Write(); - b1b3b4_eff.Write(); - table->Write(); - - input->Close(); -} diff --git a/src/Detectors/StripDetector.cpp b/src/Detectors/StripDetector.cpp index 80ad21b..b08a346 100644 --- a/src/Detectors/StripDetector.cpp +++ b/src/Detectors/StripDetector.cpp @@ -13,7 +13,7 @@ */ StripDetector::StripDetector(int ns, double len, double wid, double cphi, double cz, double crho) : - m_random(nullptr) + m_random(nullptr), m_uniform_fraction(0.0, 1.0) { num_strips = ns; @@ -88,7 +88,7 @@ Mask::Vec3 StripDetector::GetHitCoordinates(int front_stripch, double front_stri double y; //If we have a random number generator, randomize y position within pixel. Otherwise take halfway. if(m_random) { - y = -total_width/2.0 + (front_stripch + m_random->Uniform(0.0, 1.0))*front_strip_width; + y = -total_width/2.0 + (front_stripch + m_uniform_fraction(*m_random))*front_strip_width; } else { y = -total_width/2.0 + (front_stripch+0.5)*front_strip_width; } diff --git a/src/EnergyLoss.cpp b/src/EnergyLoss.cpp index 266d942..c332ecf 100644 --- a/src/EnergyLoss.cpp +++ b/src/EnergyLoss.cpp @@ -46,7 +46,7 @@ double EnergyLoss::GetEnergyLoss(int zp, int ap, double e_initial, double thickn if( ZP != zp) { ZP = zp; AP = ap; - MP = MASS.FindMass(ZP, AP)*MEV2U; + MP = MassLookup::GetInstance()->FindMass(ZP, AP)*MEV2U; } double e_final = e_initial; @@ -98,7 +98,7 @@ double EnergyLoss::GetReverseEnergyLoss(int zp, int ap, double e_final, double t if( ZP != zp) { ZP = zp; AP = ap; - MP = MASS.FindMass(ZP, AP)*MEV2U; + MP = MassLookup::GetInstance()->FindMass(ZP, AP)*MEV2U; } double e_initial = e_final; diff --git a/src/Kinematics.cpp b/src/Kinematics.cpp index 50d25a6..55b737e 100644 --- a/src/Kinematics.cpp +++ b/src/Kinematics.cpp @@ -5,471 +5,355 @@ namespace Mask { -Kinematics::Kinematics() : - sys(nullptr), save_tree_flag(false), do_plotter_flag(false), global_generator(new TRandom3(0)) -{ - std::cout<<"----------GWM Kinematics Simulation----------"<>junk>>m_outfile_name; - input>>junk>>junk; - if(junk == "yes") save_tree_flag = true; - input>>junk>>junk; - if(junk == "yes") do_plotter_flag = true; - - std::vector avec, zvec, svec; - int z, a, s; - getline(input, junk); - getline(input, junk); - input>>junk>>m_rxn_type; - getline(input, junk); - getline(input, junk); - switch(m_rxn_type) { - case 0: - { - sys = new DecaySystem(); - m_rxn_type = ONESTEP_DECAY; - for(int i=0; i<2; i++) { - input>>z>>a; - avec.push_back(a); - zvec.push_back(z); - } - break; - } - case 1: - { - sys = new OneStepSystem(); - m_rxn_type = ONESTEP_RXN; - for(int i=0; i<3; i++) { - input>>z>>a; - avec.push_back(a); - zvec.push_back(z); - } - break; - } - case 2: - { - sys = new TwoStepSystem(); - m_rxn_type = TWOSTEP; - for(int i=0; i<4; i++) { - input>>z>>a; - avec.push_back(a); - zvec.push_back(z); - } - break; - } - case 3: - { - sys = new ThreeStepSystem(); - m_rxn_type = THREESTEP; - for(int i=0; i<5; i++) { - input>>z>>a; - avec.push_back(a); - zvec.push_back(z); - } - break; - } - default: + + Kinematics::~Kinematics() { + delete global_generator; + if(sys) delete sys; + } + + bool Kinematics::LoadConfig(const std::string& filename) { + std::cout<<"Loading configuration in "<SetNuclei(zvec, avec); - - int nlayers; - double thickness; - getline(input, junk); - getline(input, junk); - input>>junk>>junk; - input>>junk>>nlayers; - for(int i=0; i>junk>>junk>>thickness; + } + + std::string junk; + getline(input, junk); + input>>junk>>m_outfile_name; + input>>junk>>junk; + input>>junk>>junk; + + std::vector avec, zvec, svec; + int z, a, s; getline(input, junk); getline(input, junk); - avec.clear(); zvec.clear(); svec.clear(); - while(input>>z) { - if(z == 0) break; - input>>a>>s; - zvec.push_back(z); avec.push_back(a); svec.push_back(s); + input>>junk>>m_rxn_type; + getline(input, junk); + getline(input, junk); + switch(m_rxn_type) { + case 0: + { + sys = new DecaySystem(); + m_rxn_type = ONESTEP_DECAY; + for(int i=0; i<2; i++) { + input>>z>>a; + avec.push_back(a); + zvec.push_back(z); + } + break; + } + case 1: + { + sys = new OneStepSystem(); + m_rxn_type = ONESTEP_RXN; + for(int i=0; i<3; i++) { + input>>z>>a; + avec.push_back(a); + zvec.push_back(z); + } + break; + } + case 2: + { + sys = new TwoStepSystem(); + m_rxn_type = TWOSTEP; + for(int i=0; i<4; i++) { + input>>z>>a; + avec.push_back(a); + zvec.push_back(z); + } + break; + } + case 3: + { + sys = new ThreeStepSystem(); + m_rxn_type = THREESTEP; + for(int i=0; i<5; i++) { + input>>z>>a; + avec.push_back(a); + zvec.push_back(z); + } + break; + } + default: + return false; } - sys->AddTargetLayer(zvec, avec, svec, thickness); - input>>junk; - } - std::cout<<"Reaction equation: "<>junk>>m_nsamples; - input>>junk>>par1>>junk>>par2; - sys->SetBeamDistro(par1, par2); - input>>junk>>par1; - switch(m_rxn_type) { - case ONESTEP_RXN : - { - dynamic_cast(sys)->SetReactionThetaType(par1); - break; - } - case TWOSTEP: - { - dynamic_cast(sys)->SetReactionThetaType(par1); - break; - } - case THREESTEP: - { - dynamic_cast(sys)->SetReactionThetaType(par1); - break; - } - } - input>>junk>>par1>>junk>>par2; - sys->SetTheta1Range(par1, par2); - input>>junk>>par1>>junk>>par2; - sys->SetPhi1Range(par1, par2); - input>>junk>>par1>>junk>>par2; - sys->SetExcitationDistro(par1, par2); - input>>junk>>dfile1; - input>>junk>>dfile2; - switch(m_rxn_type) { - case ONESTEP_DECAY: - { - DecaySystem* this_sys = dynamic_cast(sys); - this_sys->SetDecay1Distribution(dfile1); - std::cout<<"Decay1 angular momentum: "<GetDecay1AngularMomentum()<GetDecay1BranchingRatio()<(sys); - this_sys->SetDecay1Distribution(dfile1); - std::cout<<"Decay1 angular momentum: "<GetDecay1AngularMomentum()<GetDecay1BranchingRatio()<(sys); - this_sys->SetDecay1Distribution(dfile1); - this_sys->SetDecay2Distribution(dfile2); - std::cout<<"Decay1 angular momentum: "<GetDecay1AngularMomentum()<<" Decay2 angular momentum: "<GetDecay2AngularMomentum()<GetDecay1BranchingRatio()<<" Decay2 total branching ratio: "<GetDecay2BranchingRatio()<SetRandomGenerator(global_generator); - - std::cout<<"Number of samples: "<(sys); - if(this_sys == nullptr) { - return; - } - + sys->SetNuclei(zvec, avec); - TFile* output = TFile::Open(m_outfile_name.c_str(), "RECREATE"); - TTree* tree; - NucData targ, proj, eject, residual; - if(save_tree_flag) { - tree = new TTree("DataTree","DataTree"); - tree->Branch("target", &targ); - tree->Branch("projectile", &proj); - tree->Branch("ejectile", &eject); - tree->Branch("residual", &residual); + int nlayers; + double thickness; + getline(input, junk); + getline(input, junk); + input>>junk>>junk; + input>>junk>>nlayers; + for(int i=0; i>junk>>junk>>thickness; + getline(input, junk); + getline(input, junk); + avec.clear(); zvec.clear(); svec.clear(); + while(input>>z) { + if(z == 0) break; + input>>a>>s; + zvec.push_back(z); avec.push_back(a); svec.push_back(s); + } + sys->AddTargetLayer(zvec, avec, svec, thickness); + input>>junk; + } + std::cout<<"Reaction equation: "<>junk>>m_nsamples; + input>>junk>>par1>>junk>>par2; + sys->SetBeamDistro(par1, par2); + input>>junk>>par1; + switch(m_rxn_type) { + case ONESTEP_RXN : + { + dynamic_cast(sys)->SetReactionThetaType(par1); + break; + } + case TWOSTEP: + { + dynamic_cast(sys)->SetReactionThetaType(par1); + break; + } + case THREESTEP: + { + dynamic_cast(sys)->SetReactionThetaType(par1); + break; + } + } + input>>junk>>par1>>junk>>par2; + sys->SetTheta1Range(par1, par2); + input>>junk>>par1>>junk>>par2; + sys->SetPhi1Range(par1, par2); + input>>junk>>par1>>junk>>par2; + sys->SetExcitationDistro(par1, par2); + input>>junk>>dfile1; + input>>junk>>dfile2; + switch(m_rxn_type) { + case ONESTEP_DECAY: + { + DecaySystem* this_sys = dynamic_cast(sys); + this_sys->SetDecay1Distribution(dfile1); + std::cout<<"Decay1 angular momentum: "<GetDecay1AngularMomentum()<GetDecay1BranchingRatio()<(sys); + this_sys->SetDecay1Distribution(dfile1); + std::cout<<"Decay1 angular momentum: "<GetDecay1AngularMomentum()<GetDecay1BranchingRatio()<(sys); + this_sys->SetDecay1Distribution(dfile1); + this_sys->SetDecay2Distribution(dfile2); + std::cout<<"Decay1 angular momentum: "<GetDecay1AngularMomentum()<<" Decay2 angular momentum: "<GetDecay2AngularMomentum()<GetDecay1BranchingRatio()<<" Decay2 total branching ratio: "<GetDecay2BranchingRatio()<SetRandomGenerator(global_generator); + + std::cout<<"Number of samples: "< data; - data.resize(4); - for(int i=0; iRunSystem(); - - if(save_tree_flag) { - 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(this_sys->GetTarget(), "targ"); - plotman.FillData(this_sys->GetProjectile(), "proj"); - plotman.FillData(this_sys->GetEjectile(), "eject"); - plotman.FillData(this_sys->GetResidual(), "resid"); - } - - } - + bool Kinematics::SaveConfig(const std::string& filename) { return true; } - output->cd(); - if(save_tree_flag) tree->Write(tree->GetName(), TObject::kOverwrite); - if(do_plotter_flag) { - plotman.GetTable()->Write(); - plotman.ClearTable(); - } - output->Close(); - -} - -void Kinematics::RunOneStepDecay() { - DecaySystem* this_sys = dynamic_cast(sys); - if(this_sys == nullptr) { - return; - } - - TFile* output = TFile::Open(m_outfile_name.c_str(), "RECREATE"); - TTree* tree; - NucData targ, eject, residual; - if(save_tree_flag) { - tree = new TTree("DataTree","DataTree"); - tree->Branch("target", &targ); - tree->Branch("ejectile", &eject); - tree->Branch("residual", &residual); - } - - //For progress tracking - int percent5 = 0.05*m_nsamples; - int count = 0; - int npercent = 0; - - for(int i=0; iRunSystem(); - if(save_tree_flag) { - targ = ConvertNucleus(this_sys->GetTarget()); - eject = ConvertNucleus(this_sys->GetEjectile()); - residual = ConvertNucleus(this_sys->GetResidual()); - tree->Fill(); - } - if(do_plotter_flag) { - plotman.FillData(this_sys->GetTarget(), "targ"); - plotman.FillData(this_sys->GetEjectile(), "eject"); - plotman.FillData(this_sys->GetResidual(), "resid"); - } - - } - - output->cd(); - if(save_tree_flag) tree->Write(tree->GetName(), TObject::kOverwrite); - if(do_plotter_flag) { - plotman.GetTable()->Write(); - plotman.ClearTable(); - } - output->Close(); -} - -void Kinematics::RunTwoStep() { - - TwoStepSystem* this_sys = dynamic_cast(sys); - if(this_sys == nullptr) { - return; - } - - - TFile* output = TFile::Open(m_outfile_name.c_str(), "RECREATE"); - TTree* tree; - NucData targ, proj, eject, residual, break1, break2; - if(save_tree_flag) { - tree = new TTree("DataTree","DataTree"); - tree->Branch("target", &targ); - tree->Branch("projectile", &proj); - tree->Branch("ejectile", &eject); - tree->Branch("residual", &residual); - tree->Branch("breakup1", &break1); - tree->Branch("breakup2", &break2); + std::cout<RunSystem(); - if(save_tree_flag) { - targ = ConvertNucleus(this_sys->GetTarget()); - proj = ConvertNucleus(this_sys->GetProjectile()); - eject = ConvertNucleus(this_sys->GetEjectile()); - residual = ConvertNucleus(this_sys->GetResidual()); - break1 = ConvertNucleus(this_sys->GetBreakup1()); - break2 = ConvertNucleus(this_sys->GetBreakup2()); - tree->Fill(); - } - if(do_plotter_flag) { - plotman.FillData(this_sys->GetTarget(), "targ"); - plotman.FillData(this_sys->GetProjectile(), "proj"); - plotman.FillData(this_sys->GetEjectile(), "eject"); - plotman.FillData(this_sys->GetResidual(), "resid"); - plotman.FillData(this_sys->GetBreakup1(), "break1"); - plotman.FillData(this_sys->GetBreakup2(), "break2"); + void Kinematics::RunOneStepRxn() { + OneStepSystem* this_sys = dynamic_cast(sys); + if(this_sys == nullptr) { + return; } + MaskFile output(m_outfile_name, MaskFile::FileType::write); + std::vector data; + data.resize(4); + output.WriteHeader(m_rxn_type, m_nsamples); + + + //For progress tracking + int percent5 = 0.05*m_nsamples; + int count = 0; + int npercent = 0; + + for(int i=0; iRunSystem(); + data[0] = this_sys->GetTarget(); + data[1] = this_sys->GetProjectile(); + data[2] = this_sys->GetEjectile(); + data[3] = this_sys->GetResidual(); + output.WriteData(data); + + } + + output.Close(); + } + + void Kinematics::RunOneStepDecay() { + DecaySystem* this_sys = dynamic_cast(sys); + if(this_sys == nullptr) { + return; + } + + MaskFile output(m_outfile_name, MaskFile::FileType::write); + std::vector data; + data.resize(3); + output.WriteHeader(m_rxn_type, m_nsamples); + + //For progress tracking + int percent5 = 0.05*m_nsamples; + int count = 0; + int npercent = 0; + + for(int i=0; iRunSystem(); + data[0] = this_sys->GetTarget(); + data[1] = this_sys->GetEjectile(); + data[2] = this_sys->GetResidual(); + output.WriteData(data); + } + + output.Close(); + } + + void Kinematics::RunTwoStep() { + + TwoStepSystem* this_sys = dynamic_cast(sys); + if(this_sys == nullptr) { + return; + } + + MaskFile output(m_outfile_name, MaskFile::FileType::write); + std::vector data; + data.resize(6); + output.WriteHeader(m_rxn_type, m_nsamples); + + //For progress tracking + int percent5 = 0.05*m_nsamples; + int count = 0; + int npercent = 0; + + for(int i=0; iRunSystem(); + data[0] = this_sys->GetTarget(); + data[1] = this_sys->GetProjectile(); + data[2] = this_sys->GetEjectile(); + data[3] = this_sys->GetResidual(); + data[4] = this_sys->GetBreakup1(); + data[5] = this_sys->GetBreakup2(); + output.WriteData(data); + } + + output.Close(); + } + + void Kinematics::RunThreeStep() { + + ThreeStepSystem* this_sys = dynamic_cast(sys); + if(this_sys == nullptr) { + return; + } + + MaskFile output(m_outfile_name, MaskFile::FileType::write); + std::vector data; + data.resize(8); + output.WriteHeader(m_rxn_type, m_nsamples); + + //For progress updating + int percent5 = 0.05*m_nsamples; + int count = 0; + int npercent = 0; + + for(int i=0; iRunSystem(); + data[0] = this_sys->GetTarget(); + data[1] = this_sys->GetProjectile(); + data[2] = this_sys->GetEjectile(); + data[3] = this_sys->GetResidual(); + data[4] = this_sys->GetBreakup1(); + data[5] = this_sys->GetBreakup2(); + data[6] = this_sys->GetBreakup3(); + data[7] = this_sys->GetBreakup4(); + output.WriteData(data); + } + + output.Close(); } - - output->cd(); - if(save_tree_flag) tree->Write(tree->GetName(), TObject::kOverwrite); - if(do_plotter_flag) { - plotman.GetTable()->Write(); - plotman.ClearTable(); - } - output->Close(); - } - -void Kinematics::RunThreeStep() { - - ThreeStepSystem* this_sys = dynamic_cast(sys); - if(this_sys == nullptr) { - return; - } - - TFile* output = TFile::Open(m_outfile_name.c_str(), "RECREATE"); - TTree* tree; - NucData targ, proj, eject, residual, break1, break2, break3, break4; - if(save_tree_flag) { - tree = new TTree("DataTree","DataTree"); - tree->Branch("target", &targ); - tree->Branch("projectile", &proj); - tree->Branch("ejectile", &eject); - tree->Branch("residual", &residual); - tree->Branch("breakup1", &break1); - tree->Branch("breakup2", &break2); - tree->Branch("breakup3", &break3); - tree->Branch("breakup4", &break4); - } - - //For progress updating - int percent5 = 0.05*m_nsamples; - int count = 0; - int npercent = 0; - - for(int i=0; iRunSystem(); - if(save_tree_flag) { - targ = ConvertNucleus(this_sys->GetTarget()); - proj = ConvertNucleus(this_sys->GetProjectile()); - eject = ConvertNucleus(this_sys->GetEjectile()); - residual = ConvertNucleus(this_sys->GetResidual()); - break1 = ConvertNucleus(this_sys->GetBreakup1()); - break2 = ConvertNucleus(this_sys->GetBreakup2()); - break3 = ConvertNucleus(this_sys->GetBreakup3()); - break4 = ConvertNucleus(this_sys->GetBreakup4()); - tree->Fill(); - } - if(do_plotter_flag) { - plotman.FillData(this_sys->GetTarget(), "targ"); - plotman.FillData(this_sys->GetProjectile(), "proj"); - plotman.FillData(this_sys->GetEjectile(), "eject"); - plotman.FillData(this_sys->GetResidual(), "resid"); - plotman.FillData(this_sys->GetBreakup1(), "break1"); - plotman.FillData(this_sys->GetBreakup2(), "break2"); - plotman.FillData(this_sys->GetBreakup3(), "break3"); - plotman.FillData(this_sys->GetBreakup4(), "break4"); - } - } - - output->cd(); - if(save_tree_flag) tree->Write(tree->GetName(), TObject::kOverwrite); - if(do_plotter_flag) { - plotman.GetTable()->Write(); - plotman.ClearTable(); - } - output->Close(); -} - -}; diff --git a/src/MaskFile.cpp b/src/MaskFile.cpp index af9b279..5af094a 100644 --- a/src/MaskFile.cpp +++ b/src/MaskFile.cpp @@ -3,177 +3,336 @@ #include #include +/* + + FORMAT + + HEADER (contains rxntype & nuclei numbers & beam kinetic energy) (64 bits, 8 bytes) + NSAMPLES(32bit) RXNTYPE(32bit) + end HEADER + + There are NSAMPLES * (number of saved nuclei) data in the file. The number of nuclei saved is related to the + RXNTYPE. All nuclei (target, projectile, ejectile, residual, break1, etc...) are saved. A datum is as follows: + + DATUM (contains kinematic data for a nucleus) (384 bits, 48 bytes) + Z(32bit) A(32bit) DETECTFLAG(32bit) E(64bit) KE(64bit) P(64bit) THETA(64bit) PHI(64bit) + end DATUM + + +*/ + namespace Mask { -MaskFile::MaskFile() : - file_type(MaskFile::FileType::none), filename(""), buffer_position(0), count(0), file() -{ -} - -MaskFile::MaskFile(const std::string& name, MaskFile::FileType type) : - file_type(type), filename(name), buffer_position(0), count(0), file() -{ - Open(filename, type); -} - -bool MaskFile::Open(const std::string& name, MaskFile::FileType type) { - if(IsOpen()) { - std::cerr<<"Attempted to open file that is already open!"< 0) { - FlushBuffer(); - } - file.close(); -} - -void MaskFile::WriteHeader(std::vector& nuclei) { - file<& data) { - if(count == 0 && data_buffer.size() == 0) { - WriteHeader(data); + + void MaskFile::Close() { + //Final flush (if necessary) + if(buffer_position > 0 && buffer_position < buffersize_bytes && (file_type == FileType::write || file_type == FileType::append)) { + file.write(data_buffer.data(), buffer_position); + } + file.close(); } - data_buffer.push_back(data); - if(data_buffer.size() == buffersize) { - FlushBuffer(); - } -} - -void MaskFile::FlushBuffer() { - for(auto& event : data_buffer) { - file< MaskFile::ReadData() { - if(buffer_position == data_buffer.size()) { - FillBuffer(); - buffer_position = 0; - } - std::vector data = data_buffer[buffer_position]; - buffer_position++; + else if (rxn_type == 1) { + m_rxn_type = 1; + data_size = 4 * ( 5 * double_size + 2 * int_size + bool_size); + } + else if (rxn_type == 2) { + m_rxn_type = 2; + data_size = 6 * ( 5 * double_size + 2 * int_size + bool_size); + } + else if (rxn_type == 3) { + m_rxn_type = 3; + data_size = 8 * ( 5 * double_size + 2 * int_size + bool_size); + } else { + std::cerr<<"Invalid number of nuclei at MaskFile::WriteHeader! Returning"< data; - std::string line; - std::stringstream linebuf; - std::string junk; - int Z, A; - double E, p, theta, phi, tcm; - data_buffer.clear(); - - while(data_buffer.size() <= buffersize) { - if(!std::getline(file, line)) break; - linebuf.str(line); - while(linebuf>>Z) { - linebuf>>A; - linebuf>>E>>junk>>junk>>p>>theta>>phi>>tcm; - data.emplace_back(Z, A); - data[data.size()-1].SetVectorSpherical(theta, phi, p, E); - data[data.size()-1].SetThetaCM(tcm); - } - data_buffer.push_back(data); - data.clear(); + data_buffer.resize(buffersize_bytes); + + file.write((char*) &nsamples, int_size); + file.write((char*) &m_rxn_type, int_size); + } + + MaskFileHeader MaskFile::ReadHeader() { + MaskFileHeader header; + std::vector temp_buffer(4); + file.read(temp_buffer.data(), 4); + header.nsamples = *(int*)(&temp_buffer[0]); + file.read(temp_buffer.data(), 4); + m_rxn_type = *(int*)(&temp_buffer[0]); + + //size of a datum = data nuclei per event * (# doubles per nucleus * sizeof double + # ints per nucleus * sizeof int + sizeof bool) + if(m_rxn_type == 0) { + data_size = 3 * ( 5 * double_size + 2 * int_size + bool_size); + } + else if (m_rxn_type == 1) { + data_size = 4 * ( 5 * double_size + 2 * int_size + bool_size); + } + else if (m_rxn_type == 2) { + data_size = 6 * ( 5 * double_size + 2 * int_size + bool_size); + } + else if (m_rxn_type == 3) { + data_size = 8 * ( 5 * double_size + 2 * int_size + bool_size); + } else { + std::cerr<<"Unexpected reaction type at MaskFile::ReadHeader (rxn type = "<& data) { + + char* data_pointer; + double datum; + int number; + bool flag; + std::size_t j; + for(unsigned int i=0; i buffer_end) { + std::cerr<<"Attempting to read past end of file at MaskFile::ReadData! Returning empty"<>junk) { - massfile>>Z>>A>>element>>atomicMassBig>>atomicMassSmall; - isotopicMass = (atomicMassBig + atomicMassSmall*1e-6 - Z*electron_mass)*u_to_mev; - std::string key = "("+std::to_string(Z)+","+A+")"; - massTable[key] = isotopicMass; - elementTable[Z] = element; - } - } else { - throw MassFileException(); - } + std::ifstream massfile("./etc/mass.txt"); + if(massfile.is_open()) { + std::string junk, A, element; + int Z; + double atomicMassBig, atomicMassSmall, isotopicMass; + getline(massfile,junk); + getline(massfile,junk); + while(massfile>>junk) { + massfile>>Z>>A>>element>>atomicMassBig>>atomicMassSmall; + isotopicMass = (atomicMassBig + atomicMassSmall*1e-6 - Z*electron_mass)*u_to_mev; + std::string key = "("+std::to_string(Z)+","+A+")"; + massTable[key] = isotopicMass; + elementTable[Z] = element; + } + } else { + throw MassFileException(); + } } MassLookup::~MassLookup() {} //Returns nuclear mass in MeV double MassLookup::FindMass(int Z, int A) { - std::string key = "("+std::to_string(Z)+","+std::to_string(A)+")"; - auto data = massTable.find(key); - if(data == massTable.end()) { - throw MassException(); - } - return data->second; - /*try { - double mass = massTable.at(key); - return mass; - } catch (std::out_of_range& oor) { - std::cerr<<"Mass of "<second; } //returns element symbol std::string MassLookup::FindSymbol(int Z, int A) { - auto data = elementTable.find(Z); - if(data == elementTable.end()) { - throw MassException(); - } - std::string fullsymbol = std::to_string(A) + data->second; - return fullsymbol; - - /*try { - std::string element = elementTable.at(Z); - std::string fullsymbol = std::to_string(A) + element; - return fullsymbol; - } catch (std::out_of_range& oor) { - std::cerr<<"Atomic number: "<second; + return fullsymbol; } diff --git a/src/Nucleus.cpp b/src/Nucleus.cpp index f98d241..d5d621f 100644 --- a/src/Nucleus.cpp +++ b/src/Nucleus.cpp @@ -12,23 +12,23 @@ namespace Mask { Nucleus::Nucleus () : - Vec4(), m_z(0), m_a(0), m_gs_mass(0), m_theta_cm(0), m_symbol("") + Vec4(), m_z(0), m_a(0), m_gs_mass(0), m_theta_cm(0), m_symbol(""), m_detectFlag(false) { } Nucleus::Nucleus(int Z, int A) : - Vec4(), m_z(Z), m_a(A), m_theta_cm(0) + Vec4(), m_z(Z), m_a(A), m_theta_cm(0), m_detectFlag(false) { - m_gs_mass = MASS.FindMass(Z, A); - m_symbol = MASS.FindSymbol(Z, A); + m_gs_mass = MassLookup::GetInstance()->FindMass(Z, A); + m_symbol = MassLookup::GetInstance()->FindSymbol(Z, A); SetVectorCartesian(0,0,0,m_gs_mass); //by defualt a nucleus has mass given by the g.s. } Nucleus::Nucleus(int Z, int A, double px, double py, double pz, double E) : Vec4(px, py, pz, E), m_z(Z), m_a(A) { - m_gs_mass = MASS.FindMass(Z, A); - m_symbol = MASS.FindSymbol(Z, A); + m_gs_mass = MassLookup::GetInstance()->FindMass(Z, A); + m_symbol = MassLookup::GetInstance()->FindSymbol(Z, A); } Nucleus::~Nucleus() {} @@ -38,8 +38,8 @@ bool Nucleus::SetIsotope(int Z, int A) { m_z = Z; m_a = A; - m_gs_mass = MASS.FindMass(Z, A); - m_symbol = MASS.FindSymbol(Z, A); + m_gs_mass = MassLookup::GetInstance()->FindMass(Z, A); + m_symbol = MassLookup::GetInstance()->FindSymbol(Z, A); SetVectorCartesian(0,0,0,m_gs_mass); return true; } diff --git a/src/OneStepSystem.cpp b/src/OneStepSystem.cpp index d9b2b69..ad2e257 100644 --- a/src/OneStepSystem.cpp +++ b/src/OneStepSystem.cpp @@ -2,76 +2,76 @@ namespace Mask { -OneStepSystem::OneStepSystem() : - ReactionSystem() -{ -} - -OneStepSystem::OneStepSystem(std::vector& z, std::vector& a) : - ReactionSystem() -{ - SetNuclei(z, a); -} - -OneStepSystem::~OneStepSystem() {} - -bool OneStepSystem::SetNuclei(std::vector& z, std::vector& a) { - if(z.size() != a.size() || z.size() != 3) { - return false; + OneStepSystem::OneStepSystem() : + ReactionSystem() + { } - - 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(); + OneStepSystem::OneStepSystem(std::vector& z, std::vector& a) : + ReactionSystem() + { + SetNuclei(z, a); } - - 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); + OneStepSystem::~OneStepSystem() {} - step1.TurnOnResidualEloss(); - step1.Calculate(); - } else { - return; + bool OneStepSystem::SetNuclei(std::vector& z, std::vector& 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 = (*m_beamDist)(*generator); + double rxnTheta = std::acos((*m_theta1Range)(*generator)); + double rxnPhi = (*m_phi1Range)(*generator); + double residEx = (*m_exDist)(*generator); + + step1.SetBeamKE(bke); + step1.SetPolarRxnAngle(rxnTheta); + step1.SetAzimRxnAngle(rxnPhi); + step1.SetExcitation(residEx); + + step1.TurnOnResidualEloss(); + step1.Calculate(); + } else { + return; + } } -} } \ No newline at end of file diff --git a/src/Plotter.cpp b/src/Plotter.cpp deleted file mode 100644 index f9a5753..0000000 --- a/src/Plotter.cpp +++ /dev/null @@ -1,89 +0,0 @@ -#include "Plotter.h" -#include - -namespace Mask { - -Plotter::Plotter() : - table(new THashTable()) -{ -} - -Plotter::~Plotter() { - for(unsigned int i=0; iFindObject(name); - if(h) { - h->Fill(val); - } else { - h = new TH1F(name, title, bins, min, max); - h->Fill(val); - table->Add(h); - } -} - -void Plotter::MyFill(const char* name, const char* title, int binsx, float minx, float maxx, int binsy, float miny, float maxy, double valx, double valy) { - TH2F* h = (TH2F*) table->FindObject(name); - if(h) { - h->Fill(valx, valy); - } else { - h = new TH2F(name, title, binsx, minx, maxx, binsy, miny, maxy); - h->Fill(valx, valy); - table->Add(h); - } -} - -void Plotter::MyFill(const char* name, const char* title, double valx, double valy, int color) { - for(auto& g : graphs) { - if(g.name == name) { - g.xvec.push_back(valx); - g.yvec.push_back(valy); - return; - } - } - - GraphData new_g; - new_g.name = name; - new_g.title = title; - new_g.xvec.push_back(valx); - new_g.yvec.push_back(valy); - new_g.color = color; - - graphs.push_back(new_g); -} - -void Plotter::GenerateGraphs() { - for(auto& g : graphs) { - TGraph* graph = new TGraph(g.xvec.size(), &(g.xvec[0]), &(g.yvec[0])); - graph->SetName(g.name.c_str()); - graph->SetTitle(g.title.c_str()); - graph->SetMarkerColor(g.color); - table->Add(graph); - garbage_collection.push_back(graph); - } -} - -}; \ No newline at end of file diff --git a/src/Plotters/ROOT/RootPlotter.cpp b/src/Plotters/ROOT/RootPlotter.cpp new file mode 100644 index 0000000..1996863 --- /dev/null +++ b/src/Plotters/ROOT/RootPlotter.cpp @@ -0,0 +1,137 @@ +#include "RootPlotter.h" +#include "MaskFile.h" +#include + +#include + +RootPlotter::RootPlotter() : + table(new THashTable()) +{ +} + +RootPlotter::~RootPlotter() {} + +void RootPlotter::FillData(const Mask::Nucleus& nuc, 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)"; + std::string ke_vs_ph_name = sym + modifier + "_ke_vs_phi"; + std::string ke_vs_ph_title = ke_vs_ph_name + ";#phi_{lab} (degrees);Kinetic Energy (MeV)"; + std::string ex_name = sym + modifier + "_ex"; + std::string ex_title = ex_name + ";E_{ex} (MeV);counts"; + 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())); + +} + +void RootPlotter::MyFill(const std::string& name, const std::string& title, int bins, float min, float max, double val) { + TH1F* h = (TH1F*) table->FindObject(name.c_str()); + if(h) { + h->Fill(val); + } else { + h = new TH1F(name.c_str(), title.c_str(), bins, min, max); + h->Fill(val); + table->Add(h); + } +} + +void RootPlotter::MyFill(const std::string& name, const std::string& title, int binsx, float minx, float maxx, int binsy, float miny, float maxy, double valx, double valy) { + TH2F* h = (TH2F*) table->FindObject(name.c_str()); + if(h) { + h->Fill(valx, valy); + } else { + h = new TH2F(name.c_str(), title.c_str(), binsx, minx, maxx, binsy, miny, maxy); + h->Fill(valx, valy); + table->Add(h); + } +} + +void RootPlotter::MyFill(const std::string& name, const std::string& title, double valx, double valy, int color) { + for(auto& g : graphs) { + if(g.name == name) { + g.xvec.push_back(valx); + g.yvec.push_back(valy); + return; + } + } + + GraphData new_g; + new_g.name = name; + new_g.title = title; + new_g.xvec.push_back(valx); + new_g.yvec.push_back(valy); + new_g.color = color; + + graphs.push_back(new_g); +} + +void RootPlotter::GenerateGraphs() { + for(auto& g : graphs) { + TGraph* graph = new TGraph(g.xvec.size(), &(g.xvec[0]), &(g.yvec[0])); + graph->SetName(g.name.c_str()); + graph->SetTitle(g.title.c_str()); + graph->SetMarkerColor(g.color); + table->Add(graph); + garbage_collection.push_back(graph); + } +} + +int main(int argc, char** argv) { + if(argc != 3) { + std::cout<<"Unable to run ROOT plotting tool with incorrect number of arguments! Expected 2 args, given: "<cd(); + plotter.GetTable()->Write(); + root_out->Close(); + + return 0; +} \ No newline at end of file diff --git a/src/ReactionSystem.cpp b/src/ReactionSystem.cpp index 4034dff..d282832 100644 --- a/src/ReactionSystem.cpp +++ b/src/ReactionSystem.cpp @@ -4,20 +4,26 @@ namespace Mask { -ReactionSystem::ReactionSystem() : - 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("") -{ -} + ReactionSystem::ReactionSystem() : + m_beamDist(nullptr), m_theta1Range(nullptr), m_phi1Range(nullptr), m_exDist(nullptr), generator(nullptr), target_set_flag(false), + gen_set_flag(false), rxnLayer(0), m_sys_equation("") + { + } + + ReactionSystem::~ReactionSystem() { + delete m_beamDist; + delete m_theta1Range; + delete m_phi1Range; + delete m_exDist; + } + + void ReactionSystem::SetRandomGenerator(std::mt19937* gen) { + generator = gen; + gen_set_flag = true; + } + + void ReactionSystem::AddTargetLayer(std::vector& zt, std::vector& at, std::vector& stoich, double thickness) { + target.AddLayer(zt, at, stoich, thickness); + } -ReactionSystem::~ReactionSystem() {} - -void ReactionSystem::SetRandomGenerator(TRandom3* gen) { - generator = gen; - gen_set_flag = true; -} - -void ReactionSystem::AddTargetLayer(std::vector& zt, std::vector& at, std::vector& stoich, double thickness) { - target.AddLayer(zt, at, stoich, thickness); -} - -}; \ No newline at end of file +} \ No newline at end of file diff --git a/src/Stopwatch.cpp b/src/Stopwatch.cpp new file mode 100644 index 0000000..337098e --- /dev/null +++ b/src/Stopwatch.cpp @@ -0,0 +1,31 @@ +/* + Stopwatch.cpp + Simple class designed to provide timing info on parts of the process. + Only for use in development. + + Written by G.W. McCann Oct. 2020 +*/ +#include "Stopwatch.h" + +Stopwatch::Stopwatch() { + start_time = Clock::now(); + stop_time = start_time; +} + +Stopwatch::~Stopwatch() {} + +void Stopwatch::Start() { + start_time = Clock::now(); +} + +void Stopwatch::Stop() { + stop_time = Clock::now(); +} + +double Stopwatch::GetElapsedSeconds() { + return std::chrono::duration_cast>(stop_time-start_time).count(); +} + +double Stopwatch::GetElapsedMilliseconds() { + return std::chrono::duration_cast>(stop_time-start_time).count()*1000.0; +} \ No newline at end of file diff --git a/src/ThreeStepSystem.cpp b/src/ThreeStepSystem.cpp index 88c5256..83b5113 100644 --- a/src/ThreeStepSystem.cpp +++ b/src/ThreeStepSystem.cpp @@ -3,131 +3,129 @@ namespace Mask { -ThreeStepSystem::ThreeStepSystem() : - ReactionSystem() -{ -} - -ThreeStepSystem::ThreeStepSystem(std::vector& z, std::vector& a) : - ReactionSystem() -{ - SetNuclei(z, a); -} - -ThreeStepSystem::~ThreeStepSystem() { - -} - -void ThreeStepSystem::SetRandomGenerator(TRandom3* gen) { - generator = gen; - decay1dist.AttachRandomNumberGenerator(gen); - decay2dist.AttachRandomNumberGenerator(gen); - gen_set_flag = true; -} - -bool ThreeStepSystem::SetNuclei(std::vector&z, std::vector& a) { - if(z.size() != a.size() || z.size() < 5) { - return false; + ThreeStepSystem::ThreeStepSystem() : + ReactionSystem(), m_phi2Range(0, 2.0*M_PI) + { } - int zr = z[0] + z[1] - z[2]; - int zb2 = zr - z[3]; - int ar = a[0] + a[1] - a[2]; - int ab2 = ar - a[3]; - - step1.SetNuclei(z[0], a[0], z[1], a[1], z[2], a[2]); - step2.SetNuclei(zr, ar, 0, 0, z[3], a[3]); - step3.SetNuclei(zb2, ab2, 0, 0, z[4], a[4]); - SetSystemEquation(); - return true; -} - -void ThreeStepSystem::LinkTarget() { - step1.SetLayeredTarget(&target); - step2.SetLayeredTarget(&target); - step3.SetLayeredTarget(&target); - - rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA()); - if(rxnLayer != -1) { - step1.SetRxnLayer(rxnLayer); - step2.SetRxnLayer(rxnLayer); - step3.SetRxnLayer(rxnLayer); - target_set_flag = true; - } else { - throw ReactionLayerException(); - } -} - -void ThreeStepSystem::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(); - m_sys_equation += "-> "; - m_sys_equation += step2.GetEjectile().GetIsotopicSymbol(); - m_sys_equation += " + "; - m_sys_equation += step2.GetResidual().GetIsotopicSymbol(); - m_sys_equation += "-> "; - m_sys_equation += step3.GetEjectile().GetIsotopicSymbol(); - m_sys_equation += " + "; - m_sys_equation += step3.GetResidual().GetIsotopicSymbol(); -} - -void ThreeStepSystem::RunSystem() { - if(!gen_set_flag) return; - //Link up the target if it hasn't been done yet - if(!target_set_flag) { - LinkTarget(); + ThreeStepSystem::ThreeStepSystem(std::vector& z, std::vector& a) : + ReactionSystem(), m_phi2Range(0, 2.0*M_PI) + { + SetNuclei(z, a); + } + + ThreeStepSystem::~ThreeStepSystem() {} + + void ThreeStepSystem::SetRandomGenerator(std::mt19937* gen) { + generator = gen; + decay1dist.AttachRandomNumberGenerator(gen); + decay2dist.AttachRandomNumberGenerator(gen); + gen_set_flag = true; + } + + bool ThreeStepSystem::SetNuclei(std::vector&z, std::vector& a) { + if(z.size() != a.size() || z.size() < 5) { + return false; + } + int zr = z[0] + z[1] - z[2]; + int zb2 = zr - z[3]; + int ar = a[0] + a[1] - a[2]; + int ab2 = ar - a[3]; + + step1.SetNuclei(z[0], a[0], z[1], a[1], z[2], a[2]); + step2.SetNuclei(zr, ar, 0, 0, z[3], a[3]); + step3.SetNuclei(zb2, ab2, 0, 0, z[4], a[4]); + SetSystemEquation(); + return true; + } + + void ThreeStepSystem::LinkTarget() { + step1.SetLayeredTarget(&target); + step2.SetLayeredTarget(&target); + step3.SetLayeredTarget(&target); + + rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA()); + if(rxnLayer != -1) { + step1.SetRxnLayer(rxnLayer); + step2.SetRxnLayer(rxnLayer); + step3.SetRxnLayer(rxnLayer); + target_set_flag = true; + } else { + throw ReactionLayerException(); + } + } + + void ThreeStepSystem::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(); + m_sys_equation += "-> "; + m_sys_equation += step2.GetEjectile().GetIsotopicSymbol(); + m_sys_equation += " + "; + m_sys_equation += step2.GetResidual().GetIsotopicSymbol(); + m_sys_equation += "-> "; + m_sys_equation += step3.GetEjectile().GetIsotopicSymbol(); + m_sys_equation += " + "; + m_sys_equation += step3.GetResidual().GetIsotopicSymbol(); + } + + void ThreeStepSystem::RunSystem() { + if(!gen_set_flag) return; + + //Link up the target if it hasn't been done yet + if(!target_set_flag) { + LinkTarget(); + } + + //Sample parameters + double bke = (*m_beamDist)(*generator); + double rxnTheta = acos((*m_theta1Range)(*generator)); + double rxnPhi = (*m_phi1Range)(*generator); + double decay1costheta = decay1dist.GetRandomCosTheta(); + double decay1Theta = std::acos(decay1costheta); + double decay1Phi = m_phi2Range(*generator); + double decay2costheta = decay2dist.GetRandomCosTheta(); + double decay2Theta = std::acos(decay2costheta); + double decay2Phi = m_phi2Range(*generator); + double residEx = (*m_exDist)(*generator); + + step1.SetBeamKE(bke); + step1.SetPolarRxnAngle(rxnTheta); + step1.SetAzimRxnAngle(rxnPhi); + step1.SetExcitation(residEx); + + step2.SetPolarRxnAngle(decay1Theta); + step2.SetAzimRxnAngle(decay1Phi); + + step3.SetPolarRxnAngle(decay2Theta); + step3.SetAzimRxnAngle(decay2Phi); + + step1.Calculate(); + + step2.SetTarget(step1.GetResidual()); + if(decay1costheta == -10) { + step2.ResetEjectile(); + step2.ResetResidual(); + step3.ResetTarget(); + step3.ResetEjectile(); + step3.ResetResidual(); + return; + } + step2.Calculate(); + + step3.SetTarget(step2.GetResidual()); + if(decay2costheta == -10) { + step3.ResetEjectile(); + step3.ResetResidual(); + return; + } + step3.TurnOnResidualEloss(); + step3.Calculate(); + } - //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 decay1costheta = decay1dist.GetRandomCosTheta(); - double decay1Theta = std::acos(decay1costheta); - double decay1Phi = generator->Uniform(0, M_PI*2.0); - double decay2costheta = decay2dist.GetRandomCosTheta(); - double decay2Theta = std::acos(decay2costheta); - double decay2Phi = generator->Uniform(0, M_PI*2.0); - double residEx = generator->Gaus(m_exDist.first, m_exDist.second); - - step1.SetBeamKE(bke); - step1.SetPolarRxnAngle(rxnTheta); - step1.SetAzimRxnAngle(rxnPhi); - step1.SetExcitation(residEx); - - step2.SetPolarRxnAngle(decay1Theta); - step2.SetAzimRxnAngle(decay1Phi); - - step3.SetPolarRxnAngle(decay2Theta); - step3.SetAzimRxnAngle(decay2Phi); - - step1.Calculate(); - - step2.SetTarget(step1.GetResidual()); - if(decay1costheta == -10) { - step2.ResetEjectile(); - step2.ResetResidual(); - step3.ResetTarget(); - step3.ResetEjectile(); - step3.ResetResidual(); - return; - } - step2.Calculate(); - - step3.SetTarget(step2.GetResidual()); - if(decay2costheta == -10) { - step3.ResetEjectile(); - step3.ResetResidual(); - return; - } - step3.TurnOnResidualEloss(); - step3.Calculate(); - -} - -}; \ No newline at end of file +} \ No newline at end of file diff --git a/src/TwoStepSystem.cpp b/src/TwoStepSystem.cpp index c4f69d9..575cbc6 100644 --- a/src/TwoStepSystem.cpp +++ b/src/TwoStepSystem.cpp @@ -3,105 +3,105 @@ namespace Mask { -TwoStepSystem::TwoStepSystem() : - ReactionSystem() -{ -} - -TwoStepSystem::TwoStepSystem(std::vector& z, std::vector& a) : - ReactionSystem() -{ - SetNuclei(z, a); -} - -TwoStepSystem::~TwoStepSystem() { - -} - -void TwoStepSystem::SetRandomGenerator(TRandom3* gen) { - generator = gen; - decay1dist.AttachRandomNumberGenerator(gen); - gen_set_flag = true; -} - -bool TwoStepSystem::SetNuclei(std::vector&z, std::vector& a) { - if(z.size() != a.size() || z.size() != 4) { - return false; + TwoStepSystem::TwoStepSystem() : + ReactionSystem(), m_phi2Range(0, 2.0*M_PI) + { } - int zr = z[0] + z[1] - z[2]; - int ar = a[0] + a[1] - a[2]; - - step1.SetNuclei(z[0], a[0], z[1], a[1], z[2], a[2]); - step2.SetNuclei(zr, ar, 0, 0, z[3], a[3]); - SetSystemEquation(); - return true; -} - -void TwoStepSystem::LinkTarget() { - step1.SetLayeredTarget(&target); - step2.SetLayeredTarget(&target); - - rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA()); - if(rxnLayer != -1) { - step1.SetRxnLayer(rxnLayer); - step2.SetRxnLayer(rxnLayer); - target_set_flag = true; - } else { - throw ReactionLayerException(); - } -} - -void TwoStepSystem::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(); - m_sys_equation += "-> "; - m_sys_equation += step2.GetEjectile().GetIsotopicSymbol(); - m_sys_equation += "+"; - m_sys_equation += step2.GetResidual().GetIsotopicSymbol(); -} - -void TwoStepSystem::RunSystem() { - if(!gen_set_flag) return; - //Link up the target if it hasn't been done yet - if(!target_set_flag) { - LinkTarget(); + TwoStepSystem::TwoStepSystem(std::vector& z, std::vector& a) : + ReactionSystem(), m_phi2Range(0, 2.0*M_PI) + { + SetNuclei(z, a); } - - //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 decay1costheta = decay1dist.GetRandomCosTheta(); - double decay1Theta = std::acos(decay1costheta); - double decay1Phi = generator->Uniform(0, M_PI*2.0); - double residEx = generator->Gaus(m_exDist.first, m_exDist.second); - - step1.SetBeamKE(bke); - step1.SetPolarRxnAngle(rxnTheta); - step1.SetAzimRxnAngle(rxnPhi); - step1.SetExcitation(residEx); - - step2.SetPolarRxnAngle(decay1Theta); - step2.SetAzimRxnAngle(decay1Phi); - step1.Calculate(); - - step2.SetTarget(step1.GetResidual()); - if(decay1costheta == -10) { - step2.ResetEjectile(); - step2.ResetResidual(); - return; + TwoStepSystem::~TwoStepSystem() { + + } + + void TwoStepSystem::SetRandomGenerator(std::mt19937* gen) { + generator = gen; + decay1dist.AttachRandomNumberGenerator(gen); + gen_set_flag = true; + } + + bool TwoStepSystem::SetNuclei(std::vector&z, std::vector& a) { + if(z.size() != a.size() || z.size() != 4) { + return false; + } + int zr = z[0] + z[1] - z[2]; + int ar = a[0] + a[1] - a[2]; + + step1.SetNuclei(z[0], a[0], z[1], a[1], z[2], a[2]); + step2.SetNuclei(zr, ar, 0, 0, z[3], a[3]); + SetSystemEquation(); + return true; + } + + void TwoStepSystem::LinkTarget() { + step1.SetLayeredTarget(&target); + step2.SetLayeredTarget(&target); + + rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA()); + if(rxnLayer != -1) { + step1.SetRxnLayer(rxnLayer); + step2.SetRxnLayer(rxnLayer); + target_set_flag = true; + } else { + throw ReactionLayerException(); + } + } + + void TwoStepSystem::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(); + m_sys_equation += "-> "; + m_sys_equation += step2.GetEjectile().GetIsotopicSymbol(); + m_sys_equation += "+"; + m_sys_equation += step2.GetResidual().GetIsotopicSymbol(); + } + + void TwoStepSystem::RunSystem() { + if(!gen_set_flag) return; + + //Link up the target if it hasn't been done yet + if(!target_set_flag) { + LinkTarget(); + } + + //Sample parameters + double bke = (*m_beamDist)(*generator); + double rxnTheta = acos((*m_theta1Range)(*generator)); + double rxnPhi = (*m_phi1Range)(*generator); + double decay1costheta = decay1dist.GetRandomCosTheta(); + double decay1Theta = std::acos(decay1costheta); + double decay1Phi = m_phi2Range(*generator); + double residEx = (*m_beamDist)(*generator); + + step1.SetBeamKE(bke); + step1.SetPolarRxnAngle(rxnTheta); + step1.SetAzimRxnAngle(rxnPhi); + step1.SetExcitation(residEx); + + step2.SetPolarRxnAngle(decay1Theta); + step2.SetAzimRxnAngle(decay1Phi); + + step1.Calculate(); + + step2.SetTarget(step1.GetResidual()); + if(decay1costheta == -10) { + step2.ResetEjectile(); + step2.ResetResidual(); + return; + } + step2.TurnOnResidualEloss(); + step2.Calculate(); + + } - step2.TurnOnResidualEloss(); - step2.Calculate(); - -} - -}; \ No newline at end of file +} \ No newline at end of file diff --git a/src/main.cpp b/src/main.cpp index 197e8af..3fd5446 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,21 +1,20 @@ #include #include +#include "Stopwatch.h" +#include "MaskFile.h" #include "Kinematics.h" -#include "SabreEfficiency.h" -#include "AnasenEfficiency.h" #include "KinematicsExceptions.h" -#include -#include - int main(int argc, char** argv) { if(argc<2) { std::cerr<<"Incorrect number of arguments!"<