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

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.

Browse Source
This commit is contained in:
Gordon McCann 2021-09-03 17:03:41 -04:00
parent 46eb87e65b
commit 170cc7afb0
39 changed files with 1917 additions and 2116 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
include/AnasenEfficiency.h
View File

@ -2,7 +2,6 @@
#define ANASEN_EFFICIENCY_H
#include <string>
#include <THashTable.h>
#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);

8
include/AngularDistribution.h
View File

@ -1,9 +1,9 @@
#ifndef ANGULARDISTRIBUTION_H
#define ANGULARDISTRIBUTION_H
#include <TRandom3.h>
#include <string>
#include <vector>
#include <random>
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<double> uniform_cosine_dist;
std::uniform_real_distribution<double> uniform_prob_dist;
double branchingRatio;
int L;

42
include/DecaySystem.h
View File

@ -6,36 +6,36 @@
namespace Mask {
class DecaySystem: public ReactionSystem {
public:
DecaySystem();
DecaySystem(std::vector<int>& z, std::vector<int>& a);
~DecaySystem();
class DecaySystem: public ReactionSystem {
public:
DecaySystem();
DecaySystem(std::vector<int>& z, std::vector<int>& a);
~DecaySystem();
bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
void RunSystem() override;
void SetRandomGenerator(TRandom3* gen) override;
bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
void RunSystem() override;
void SetRandomGenerator(std::mt19937* gen) override;
inline void SetDecay1Distribution(const std::string& filename) { decay1dist.ReadDistributionFile(filename); };
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 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(); };
inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); }
inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); }
private:
void LinkTarget() override;
void SetSystemEquation() override;
private:
void LinkTarget() override;
void SetSystemEquation() override;
Reaction step1;
Reaction step1;
AngularDistribution decay1dist;
AngularDistribution decay1dist;
};
};
};
}
#endif

36
include/DetectorEfficiency.h
View File

@ -1,52 +1,22 @@
#ifndef DETECTOREFFICIENCY_H
#define DETECTOREFFICIENCY_H
#include <THashTable.h>
#include <TH1.h>
#include <TH2.h>
#include <string>
#include <cmath>
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

87
include/Kinematics.h
View File

@ -6,66 +6,45 @@
#include "OneStepSystem.h"
#include "TwoStepSystem.h"
#include "ThreeStepSystem.h"
#include "Plotter.h"
#include <TFile.h>
#include <TTree.h>
#include <random>
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();
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_RXN,
ONESTEP_DECAY,
TWOSTEP,
THREESTEP
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

80
include/MaskFile.h
View File

@ -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<double> E, KE, p, theta, phi; //ordered: target, (if not decay)projectile, ejectile, residual, break1...
std::vector<int> Z, A;
std::vector<bool> detect_flag;
bool eof = false; //flag on end of file
};
void WriteHeader(std::vector<Nucleus>& data);
void WriteData(std::vector<Nucleus>& data);
void ReadHeader();
std::vector<Nucleus> ReadData();
class MaskFile {
public:
enum class FileType {
read,
write,
append,
none
};
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<Nucleus>& data);
void WriteData(MaskFileData& data);
MaskFileHeader ReadHeader();
MaskFileData ReadData();
private:
void FlushBuffer();
void FillBuffer();
private:
FileType file_type;
std::string filename;
unsigned int buffer_position;
int count;
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::fstream file;
std::vector<std::vector<Nucleus>> data_buffer;
std::vector<char> data_buffer;
static constexpr int buffersize = 1000;
static constexpr int width = 0;
static constexpr int precision = 3;
};
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;
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);
};
};

13
include/MassLookup.h
View File

@ -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<std::string, double> massTable;
std::unordered_map<int, std::string> 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

5
include/Nucleus.h
View File

@ -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;
};
};

36
include/OneStepSystem.h
View File

@ -5,29 +5,29 @@
namespace Mask {
class OneStepSystem: public ReactionSystem {
public:
OneStepSystem();
OneStepSystem(std::vector<int>& z, std::vector<int>& a);
~OneStepSystem();
class OneStepSystem: public ReactionSystem {
public:
OneStepSystem();
OneStepSystem(std::vector<int>& z, std::vector<int>& a);
~OneStepSystem();
bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
void RunSystem() override;
bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
void RunSystem() override;
inline void SetReactionThetaType(int type) { step1.SetEjectileThetaType(type); };
inline const Nucleus& GetTarget() const { return step1.GetTarget(); };
inline const Nucleus& GetProjectile() const { return step1.GetProjectile(); };
inline const Nucleus& GetEjectile() const { return step1.GetEjectile(); };
inline const Nucleus& GetResidual() const { return step1.GetResidual(); };
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;
private:
void LinkTarget() override;
void SetSystemEquation() override;
Reaction step1;
Reaction step1;
};
};
};
}
#endif

49
include/Plotter.h
View File

@ -1,49 +0,0 @@
#ifndef PLOTTER_H
#define PLOTTER_H
#include <TROOT.h>
#include <TH1.h>
#include <TH2.h>
#include <TGraph.h>
#include <THashTable.h>
#include "Nucleus.h"
namespace Mask {
struct GraphData {
std::string name;
std::string title;
std::vector<double> xvec;
std::vector<double> 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<TGraph*> garbage_collection;
std::vector<GraphData> graphs;
static constexpr double rad2deg = 180.0/M_PI;
};
};
#endif

2
include/QQQDetector.h
View File

@ -8,8 +8,8 @@
#ifndef QQQDETECTOR_H
#define QQQDETECTOR_H
#include <TRandom3.h>
#include <cmath>
#include <vector>
#include "Vec3.h"
#include "Rotation.h"

76
include/ReactionSystem.h
View File

@ -12,45 +12,65 @@
#include "Reaction.h"
#include "KinematicsExceptions.h"
#include <vector>
#include <TRandom3.h>
#include <random>
namespace Mask {
class ReactionSystem {
public:
ReactionSystem();
virtual ~ReactionSystem();
class ReactionSystem {
public:
ReactionSystem();
virtual ~ReactionSystem();
virtual bool SetNuclei(std::vector<int>& z, std::vector<int>& a) = 0;
virtual void RunSystem() = 0;
virtual bool SetNuclei(std::vector<int>& z, std::vector<int>& a) = 0;
virtual void RunSystem() = 0;
void AddTargetLayer(std::vector<int>& zt, std::vector<int>& at, std::vector<int>& stoich, double thickness);
void AddTargetLayer(std::vector<int>& zt, std::vector<int>& at, std::vector<int>& 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); };
/*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<double>(mean, sigma);
}
inline const std::string& GetSystemEquation() const { return m_sys_equation; };
inline void SetTheta1Range(double min, double max) {
if(m_theta1Range)
delete m_theta1Range;
m_theta1Range = new std::uniform_real_distribution<double>(std::cos(min*deg2rad), std::cos(max*deg2rad));
}
protected:
virtual void LinkTarget() = 0;
virtual void SetSystemEquation() = 0;
inline void SetPhi1Range(double min, double max) {
if(m_phi1Range)
delete m_phi1Range;
m_phi1Range = new std::uniform_real_distribution<double>(min*deg2rad, max*deg2rad);
}
LayeredTarget target;
inline void SetExcitationDistro(double mean, double sigma) {
if(m_exDist)
delete m_exDist;
m_exDist = new std::normal_distribution<double>(mean, sigma);
}
//Sampling information
std::pair<double, double> m_beamDist, m_theta1Range, m_phi1Range, m_exDist;
TRandom3* generator; //not owned by ReactionSystem
inline const std::string& GetSystemEquation() const { return m_sys_equation; }
bool target_set_flag, gen_set_flag;
int rxnLayer;
std::string m_sys_equation;
static constexpr double deg2rad = M_PI/180.0;
};
protected:
virtual void LinkTarget() = 0;
virtual void SetSystemEquation() = 0;
};
LayeredTarget target;
//Sampling information
std::normal_distribution<double> *m_beamDist, *m_exDist;
std::uniform_real_distribution<double> *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;
};
}
#endif

48
include/RootPlotter.h Normal file
View File

@ -0,0 +1,48 @@
#ifndef ROOTPLOTTER_H
#define ROOTPLOTTER_H
#include <vector>
#include <string>
#include "Nucleus.h"
#include <THashTable.h>
#include <TH1.h>
#include <TH2.h>
#include <TGraph.h>
struct GraphData {
std::string name;
std::string title;
std::vector<double> xvec;
std::vector<double> 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<TGraph*> garbage_collection;
std::vector<GraphData> graphs;
static constexpr double rad2deg = 180.0/M_PI;
};
#endif

11
include/SabreEfficiency.h
View File

@ -5,24 +5,19 @@
#include "SabreDetector.h"
#include "Target.h"
#include "DeadChannelMap.h"
#include "Kinematics.h"
#include <THashTable.h>
#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<bool,double> 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<bool,double> IsSabre(Mask::Nucleus& nucleus);
std::vector<SabreDetector> detectors;

30
include/Stopwatch.h Normal file
View File

@ -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 <chrono>
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

8
include/StripDetector.h
View File

@ -11,9 +11,9 @@
//Back strips from lowest z to highest z
#include <TRandom3.h>
#include <cmath>
#include <vector>
#include <random>
#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<int,double> GetChannelRatio(double theta, double phi);
@ -51,7 +51,9 @@ private:
Mask::ZRotation zRot;
TRandom3* m_random; //Not owned by StripDetector!
std::uniform_real_distribution<double> 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); };

62
include/ThreeStepSystem.h
View File

@ -6,43 +6,45 @@
namespace Mask {
class ThreeStepSystem : public ReactionSystem {
public:
ThreeStepSystem();
ThreeStepSystem(std::vector<int>& z, std::vector<int>& a);
~ThreeStepSystem();
bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
void RunSystem() override;
void SetRandomGenerator(TRandom3* gen) override;
class ThreeStepSystem : public ReactionSystem {
public:
ThreeStepSystem();
ThreeStepSystem(std::vector<int>& z, std::vector<int>& a);
~ThreeStepSystem();
bool SetNuclei(std::vector<int>& z, std::vector<int>& 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 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 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(); };
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;
protected:
void LinkTarget() override;
void SetSystemEquation() override;
Reaction step1, step2, step3;
std::uniform_real_distribution<double> m_phi2Range;
AngularDistribution decay1dist, decay2dist;
Reaction step1, step2, step3;
};
AngularDistribution decay1dist, decay2dist;
};
};
}
#endif

52
include/TwoStepSystem.h
View File

@ -6,39 +6,41 @@
namespace Mask {
class TwoStepSystem : public ReactionSystem {
public:
TwoStepSystem();
TwoStepSystem(std::vector<int>& z, std::vector<int>& a);
~TwoStepSystem();
bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
void RunSystem() override;
class TwoStepSystem : public ReactionSystem {
public:
TwoStepSystem();
TwoStepSystem(std::vector<int>& z, std::vector<int>& a);
~TwoStepSystem();
bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
void RunSystem() override;
void SetRandomGenerator(TRandom3* gen) override;
void SetRandomGenerator(std::mt19937* gen) override;
inline void SetDecay1Distribution(const std::string& filename) { decay1dist.ReadDistributionFile(filename); };
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 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(); };
inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); };
inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); };
private:
void LinkTarget() override;
void SetSystemEquation() override;
private:
void LinkTarget() override;
void SetSystemEquation() override;
Reaction step1, step2;
std::uniform_real_distribution<double> m_phi2Range;
AngularDistribution decay1dist;
Reaction step1, step2;
};
AngularDistribution decay1dist;
};
};
}
#endif

6
input.txt
View File

@ -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

63
premake5.lua
View File

@ -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"
@ -59,3 +80,33 @@ project "Mask"
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"
filter "configurations:Release"
optimize "On"

10
src/AngularDistribution.cpp
View File

@ -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<constants.size(); i++)
probability += constants[i]*P_l(i*2, costheta);
} while(test > probability);

126
src/DecaySystem.cpp
View File

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

398
src/Detectors/AnasenEfficiency.cpp
View File

@ -1,11 +1,8 @@
#include "AnasenEfficiency.h"
#include "Kinematics.h"
#include <TH2.h>
#include <TH1.h>
#include <TGraph2D.h>
#include <TGraph.h>
#include <TCanvas.h>
#include <TParameter.h>
#include "MaskFile.h"
#include <fstream>
#include <iomanip>
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<double> x, y, z;
std::vector<double> 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"<<std::endl;
output<<"Edges: x y z"<<std::endl;
for(unsigned int i=0; i<x.size(); i++) {
output<<x[i]<<" "<<y[i]<<" "<<z[i]<<std::endl;
}
output<<"Centers: x y z"<<std::endl;
for(unsigned int i=0; i<cx.size(); i++) {
output<<cx[i]<<" "<<cy[i]<<" "<<cz[i]<<std::endl;
}
TCanvas* canvas = new TCanvas();
canvas->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----------"<<std::endl;
std::cout<<"Loading in output from kinematics simulation: "<<file<<std::endl;
std::cout<<"Loading in output from kinematics simulation: "<<inputname<<std::endl;
std::cout<<"Running efficiency calculation..."<<std::endl;
switch(m_rxn_type) {
case Mask::Kinematics::ONESTEP_DECAY:
{
RunDecay(file);
Mask::MaskFile input(inputname, Mask::MaskFile::FileType::read);
Mask::MaskFile output(outputname, Mask::MaskFile::FileType::write);
std::ofstream stats(statsname);
stats<<std::setprecision(5);
Mask::MaskFileHeader header = input.ReadHeader();
output.WriteHeader(header.rxn_type, header.nsamples);
stats<<"Efficiency statistics for data from "<<inputname<<" using the ANASEN geometry"<<std::endl;
stats<<"Given in order of target=0, projectile=1, ejectile=2, residual=3, .... etc."<<std::endl;
Mask::MaskFileData data;
std::vector<int> 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: "<<header.rxn_type<<"). Quiting..."<<std::endl;
input.Close();
output.Close();
stats.close();
return;
}
}
int percent5 = header.nsamples*0.05;
int count = 0;
int npercent = 0;
while(true) {
if(++count == percent5) {//Show progress every 5%
npercent++;
count = 0;
std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
}
data = input.ReadData();
if(data.eof)
break;
for(unsigned int i=0; i<data.Z.size(); i++) {
if(data.KE[i] >= 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<<std::setw(10)<<"Index"<<"|"<<std::setw(10)<<"Efficiency"<<std::endl;
stats<<"---------------------"<<std::endl;
for(unsigned int i=0; i<counts.size(); i++) {
stats<<std::setw(10)<<i<<"|"<<std::setw(10)<<((double)counts[i]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
}
stats.close();
std::cout<<std::endl;
std::cout<<"Complete."<<std::endl;
std::cout<<"---------------------------------------------"<<std::endl;
}
void AnasenEfficiency::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;
int detected_eject=0;
int detected_resid=0;
for(int i=0; i<tree->GetEntries(); i++) {
if(++count == percent5) {//Show progress every 5%
npercent++;
count = 0;
std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
}
tree->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<double> eject_eff("Light Breakup Efficiency", ejecteff);
TParameter<double> 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!"<<std::endl;
return;
}
void AnasenEfficiency::Run2Step(const std::string& filename) {
TFile* input = TFile::Open(filename.c_str(), "UPDATE");
TTree* tree = (TTree*) input->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; i<tree->GetEntries(); i++) {
if(++count == percent5) {//Show progress every 5%
npercent++;
count = 0;
std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
}
break1_flag = false;
break2_flag = false;
eject_flag = false;
tree->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<double> eject_eff("Ejectile Efficiency", eeff);
TParameter<double> break1_eff("Breakup 1 Efficiency", b1eff);
TParameter<double> break2_eff("Breakup 2 Efficiency", b2eff);
TParameter<double> break1break2_eff("Breakup1 Coincident with Breakup2", b1b2eff);
TParameter<double> break1eject_eff("Breakup1 Coincident with Ejectile", b1eeff);
TParameter<double> break2eject_eff("Breakup2 Coincident with Ejectile", b2eeff);
TParameter<double> 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; i<tree->GetEntries(); i++) {
if(++count == percent5) {//Show progress every 5%
npercent++;
count = 0;
std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
}
break1_flag = false;
break3_flag = false;
break4_flag = false;
tree->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<double> break1_eff("Breakup1 Efficiency", b1eff);
TParameter<double> break3_eff("Breakup3 Efficiency", b3eff);
TParameter<double> break4_eff("Breakup4 Efficiency", b4eff);
TParameter<double> b1b3_eff("Break1 Coincident with Break3", b1b3eff);
TParameter<double> b1b4_eff("Break1 Coincident with Break4", b1b4eff);
TParameter<double> b3b4_eff("Break3 Coincident with Break4", b3b4eff);
TParameter<double> 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();
}

33
src/Detectors/DetectorEfficiency.cpp Normal file
View File

@ -0,0 +1,33 @@
#include "SabreEfficiency.h"
#include "AnasenEfficiency.h"
#include <iostream>
#include <string>
int main(int argc, char** argv) {
if(argc != 4) {
std::cerr<<"Incorrect number of commandline arguments! Returning."<<std::endl;
return 1;
}
std::string inputname = argv[1];
std::string outputname = argv[2];
std::string statsname = argv[3];
/*
SabreEfficiency sabre;
std::string mapfile = "./etc/DeadChannels.txt";
sabre.SetDeadChannelMap(mapfile);
sabre.CalculateEfficiency(inputname, outputname, statsname);
//std::cout<<"Running consistency check(1=success): "<<sabre.RunConsistencyCheck()<<std::endl;
//sabre.DrawDetectorSystem("/data1/gwm17/10B3He/Feb2021/simulation/SABREGeo.txt");
*/
AnasenEfficiency anasen;
anasen.CalculateEfficiency(inputname, outputname, statsname);
//std::cout<<"Running consistency check(1=success): "<<anasen.RunConsistencyCheck()<<std::endl;
//anasen.DrawDetectorSystem("/data1/gwm17/TRIUMF_7Bed/simulation/ANASENGeo_centered_target_targetGap_BackQQQ_test.root");
return 0;
}

385
src/Detectors/SabreEfficiency.cpp
View File

@ -1,12 +1,9 @@
#include "SabreEfficiency.h"
#include "MaskFile.h"
#include <fstream>
#include <iostream>
#include <TFile.h>
#include <TTree.h>
#include <TParameter.h>
#include <TGraph.h>
#include <TGraph2D.h>
#include <TCanvas.h>
#include <iomanip>
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----------"<<std::endl;
std::cout<<"Loading in output from kinematics simulation: "<<file<<std::endl;
std::cout<<"Loading in output from kinematics simulation: "<<inputname<<std::endl;
std::cout<<"Running efficiency calculation..."<<std::endl;
if(!dmap.IsValid()) {
@ -40,35 +37,92 @@ void SabreEfficiency::CalculateEfficiency(const std::string& file) {
std::cout<<"---------------------------------------------"<<std::endl;
}
switch(m_rxn_type) {
case Mask::Kinematics::ONESTEP_DECAY:
{
RunDecay(file);
Mask::MaskFile input(inputname, Mask::MaskFile::FileType::read);
Mask::MaskFile output(outputname, Mask::MaskFile::FileType::read);
std::ofstream stats(statsname);
stats<<std::setprecision(5);
Mask::MaskFileHeader header = input.ReadHeader();
output.WriteHeader(header.rxn_type, header.nsamples);
stats<<"Efficiency statistics for data from "<<inputname<<" using the SPS-SABRE geometry"<<std::endl;
stats<<"Given in order of target=0, projectile=1, ejectile=2, residual=3, .... etc."<<std::endl;
Mask::MaskFileData data;
Mask::Nucleus nucleus;
std::vector<int> 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: "<<header.rxn_type<<"). Quiting..."<<std::endl;
input.Close();
output.Close();
stats.close();
return;
}
}
int percent5 = header.nsamples*0.05;
int count = 0;
int npercent = 0;
while(true) {
if(++count == percent5) {//Show progress every 5%
npercent++;
count = 0;
std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
}
data = input.ReadData();
if(data.eof)
break;
for(unsigned int i=0; i<data.Z.size(); i++) {
nucleus.SetIsotope(data.Z[i], data.A[i]);
nucleus.SetVectorSpherical(data.theta[i], data.phi[i], data.p[i], data.E[i]);
auto result = IsSabre(nucleus);
if(result.first) {
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<<std::setw(10)<<"Index"<<"|"<<std::setw(10)<<"Efficiency"<<std::endl;
stats<<"---------------------"<<std::endl;
for(unsigned int i=0; i<counts.size(); i++) {
stats<<std::setw(10)<<i<<"|"<<std::setw(10)<<((double)counts[i]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
}
stats.close();
std::cout<<std::endl;
std::cout<<"Complete."<<std::endl;
std::cout<<"---------------------------------------------"<<std::endl;
}
void SabreEfficiency::DrawDetectorSystem(const std::string& filename) {
TFile* output = TFile::Open(filename.c_str(), "RECREATE");
std::ofstream output(filename);
std::vector<double> ringxs, ringys, ringzs;
std::vector<double> 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"<<std::endl;
output<<"Edges: x y z"<<std::endl;
for(unsigned int i=0; i<ringxs.size(); i++) {
output<<ringxs[i]<<" "<<ringys[i]<<" "<<ringzs[i]<<std::endl;
}
for(unsigned int i=0; i<wedgexs.size(); i++) {
output<<wedgexs[i]<<" "<<wedgeys[i]<<" "<<wedgezs[i]<<std::endl;
}
TGraph2D* gw = new TGraph2D(wedgexs.size(), &(wedgexs[0]), &(wedgeys[0]), &(wedgezs[0]));
gw->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<bool,double> SabreEfficiency::IsSabre(Mask::NucData* nucleus) {
if(nucleus->KE <= ENERGY_THRESHOLD) {
std::pair<bool,double> 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<bool, double> eject_result, resid_result;
int detected_eject=0, detected_resid=0;
for(int i=0; i<tree->GetEntries(); i++) {
if(++count == percent5) {//Show progress every 5%
npercent++;
count = 0;
std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
}
tree->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<double> eject_eff("Ejectile Efficiency", ejecteff);
TParameter<double> 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!"<<std::endl;
return;
}
void SabreEfficiency::Run2Step(const std::string& filename) {
TFile* input = TFile::Open(filename.c_str(), "UPDATE");
TTree* tree = (TTree*) input->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<bool, double> break1_result, break2_result;
int detected_b1=0, detected_b2=0;
for(int i=0; i<tree->GetEntries(); i++) {
if(++count == percent5) {//Show progress every 5%
npercent++;
count = 0;
std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
}
tree->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<double> break1_eff("Breakup1 Efficiency", b1eff);
TParameter<double> 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<bool, double> 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; i<tree->GetEntries(); i++) {
if(++count == percent5) {//Show progress every 5%
npercent++;
count = 0;
std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
}
break1_result = IsSabre(break1);
break3_result = IsSabre(break3);
break4_result = IsSabre(break4);
tree->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<double> break1_eff("Breakup1 Efficiency", b1eff);
TParameter<double> break3_eff("Breakup3 Efficiency", b3eff);
TParameter<double> break4_eff("Breakup4 Efficiency", b4eff);
TParameter<double> b1b3_eff("Break1 Coincident with Break3", b1b3eff);
TParameter<double> b1b4_eff("Break1 Coincident with Break4", b1b4eff);
TParameter<double> b3b4_eff("Break3 Coincident with Break4", b3b4eff);
TParameter<double> 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();
}

4
src/Detectors/StripDetector.cpp
View File

@ -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;
}

4
src/EnergyLoss.cpp
View File

@ -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;

758
src/Kinematics.cpp
View File

@ -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----------"<<std::endl;
}
Kinematics::Kinematics() :
sys(nullptr)
{
std::cout<<"----------GWM Kinematics Simulation----------"<<std::endl;
Kinematics::~Kinematics() {
delete global_generator;
if(sys) delete sys;
}
bool Kinematics::LoadConfig(const char* filename) {
std::cout<<"Loading configuration in "<<filename<<"..."<<std::endl;
std::ifstream input(filename);
if(!input.is_open()) {
std::cerr<<"Unable to load configuration in "<<filename<<", check that it exists"<<std::endl;
return false;
std::random_device rd;
global_generator = new std::mt19937(rd());
}
std::string junk;
getline(input, junk);
input>>junk>>m_outfile_name;
input>>junk>>junk;
if(junk == "yes") save_tree_flag = true;
input>>junk>>junk;
if(junk == "yes") do_plotter_flag = true;
Kinematics::~Kinematics() {
delete global_generator;
if(sys) delete sys;
}
std::vector<int> 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:
bool Kinematics::LoadConfig(const std::string& filename) {
std::cout<<"Loading configuration in "<<filename<<"..."<<std::endl;
std::ifstream input(filename);
if(!input.is_open()) {
std::cerr<<"Unable to load configuration in "<<filename<<", check that it exists"<<std::endl;
return false;
}
sys->SetNuclei(zvec, avec);
}
int nlayers;
double thickness;
getline(input, junk);
getline(input, junk);
input>>junk>>junk;
input>>junk>>nlayers;
for(int i=0; i<nlayers; i++) {
input>>junk>>junk>>thickness;
std::string junk;
getline(input, junk);
input>>junk>>m_outfile_name;
input>>junk>>junk;
input>>junk>>junk;
std::vector<int> 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: "<<GetSystemName()<<std::endl;
sys->SetNuclei(zvec, avec);
double par1, par2;
std::string dfile1, dfile2;
getline(input, junk);
getline(input, junk);
int nlayers;
double thickness;
getline(input, junk);
getline(input, junk);
input>>junk>>junk;
input>>junk>>nlayers;
for(int i=0; i<nlayers; i++) {
input>>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: "<<GetSystemName()<<std::endl;
input>>junk>>m_nsamples;
input>>junk>>par1>>junk>>par2;
sys->SetBeamDistro(par1, par2);
input>>junk>>par1;
switch(m_rxn_type) {
case ONESTEP_RXN :
{
dynamic_cast<OneStepSystem*>(sys)->SetReactionThetaType(par1);
break;
}
case TWOSTEP:
{
dynamic_cast<TwoStepSystem*>(sys)->SetReactionThetaType(par1);
break;
}
case THREESTEP:
{
dynamic_cast<ThreeStepSystem*>(sys)->SetReactionThetaType(par1);
break;
}
}
input>>junk>>par1>>junk>>par2;
sys->SetTheta1Range(par1, par2);
input>>junk>>par1>>junk>>par2;
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<DecaySystem*>(sys);
this_sys->SetDecay1Distribution(dfile1);
std::cout<<"Decay1 angular momentum: "<<this_sys->GetDecay1AngularMomentum()<<std::endl;
std::cout<<"Decay1 total branching ratio: "<<this_sys->GetDecay1BranchingRatio()<<std::endl;
break;
}
case TWOSTEP:
{
TwoStepSystem* this_sys = dynamic_cast<TwoStepSystem*>(sys);
this_sys->SetDecay1Distribution(dfile1);
std::cout<<"Decay1 angular momentum: "<<this_sys->GetDecay1AngularMomentum()<<std::endl;
std::cout<<"Decay1 total branching ratio: "<<this_sys->GetDecay1BranchingRatio()<<std::endl;
break;
}
case THREESTEP:
{
ThreeStepSystem* this_sys = dynamic_cast<ThreeStepSystem*>(sys);
this_sys->SetDecay1Distribution(dfile1);
this_sys->SetDecay2Distribution(dfile2);
std::cout<<"Decay1 angular momentum: "<<this_sys->GetDecay1AngularMomentum()<<" Decay2 angular momentum: "<<this_sys->GetDecay2AngularMomentum()<<std::endl;
std::cout<<"Decay1 total branching ratio: "<<this_sys->GetDecay1BranchingRatio()<<" Decay2 total branching ratio: "<<this_sys->GetDecay2BranchingRatio()<<std::endl;
break;
}
}
sys->SetRandomGenerator(global_generator);
double par1, par2;
std::string dfile1, dfile2;
getline(input, junk);
getline(input, junk);
std::cout<<"Number of samples: "<<GetNumberOfSamples()<<std::endl;
return true;
}
bool Kinematics::SaveConfig(const char* filename) { return true; }
NucData Kinematics::ConvertNucleus(const Nucleus& nuc) {
NucData datum;
datum.E = nuc.GetE();
datum.KE = nuc.GetKE();
datum.p = nuc.GetP();
datum.theta = nuc.GetTheta();
datum.theta_cm = nuc.GetThetaCM();
datum.phi = nuc.GetPhi();
datum.Ex = nuc.GetExcitationEnergy();
datum.Z = nuc.GetZ();
datum.A = nuc.GetA();
return datum;
}
void Kinematics::Run() {
std::cout<<"Running simulation..."<<std::endl;
switch(m_rxn_type) {
case ONESTEP_DECAY:
{
RunOneStepDecay();
break;
input>>junk>>m_nsamples;
input>>junk>>par1>>junk>>par2;
sys->SetBeamDistro(par1, par2);
input>>junk>>par1;
switch(m_rxn_type) {
case ONESTEP_RXN :
{
dynamic_cast<OneStepSystem*>(sys)->SetReactionThetaType(par1);
break;
}
case TWOSTEP:
{
dynamic_cast<TwoStepSystem*>(sys)->SetReactionThetaType(par1);
break;
}
case THREESTEP:
{
dynamic_cast<ThreeStepSystem*>(sys)->SetReactionThetaType(par1);
break;
}
}
case ONESTEP_RXN:
{
RunOneStepRxn();
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<DecaySystem*>(sys);
this_sys->SetDecay1Distribution(dfile1);
std::cout<<"Decay1 angular momentum: "<<this_sys->GetDecay1AngularMomentum()<<std::endl;
std::cout<<"Decay1 total branching ratio: "<<this_sys->GetDecay1BranchingRatio()<<std::endl;
break;
}
case TWOSTEP:
{
TwoStepSystem* this_sys = dynamic_cast<TwoStepSystem*>(sys);
this_sys->SetDecay1Distribution(dfile1);
std::cout<<"Decay1 angular momentum: "<<this_sys->GetDecay1AngularMomentum()<<std::endl;
std::cout<<"Decay1 total branching ratio: "<<this_sys->GetDecay1BranchingRatio()<<std::endl;
break;
}
case THREESTEP:
{
ThreeStepSystem* this_sys = dynamic_cast<ThreeStepSystem*>(sys);
this_sys->SetDecay1Distribution(dfile1);
this_sys->SetDecay2Distribution(dfile2);
std::cout<<"Decay1 angular momentum: "<<this_sys->GetDecay1AngularMomentum()<<" Decay2 angular momentum: "<<this_sys->GetDecay2AngularMomentum()<<std::endl;
std::cout<<"Decay1 total branching ratio: "<<this_sys->GetDecay1BranchingRatio()<<" Decay2 total branching ratio: "<<this_sys->GetDecay2BranchingRatio()<<std::endl;
break;
}
}
case TWOSTEP:
{
RunTwoStep();
break;
}
case THREESTEP:
{
RunThreeStep();
break;
}
}
std::cout<<std::endl;
std::cout<<"Complete."<<std::endl;
std::cout<<"---------------------------------------------"<<std::endl;
}
sys->SetRandomGenerator(global_generator);
void Kinematics::RunOneStepRxn() {
OneStepSystem* this_sys = dynamic_cast<OneStepSystem*>(sys);
if(this_sys == nullptr) {
return;
std::cout<<"Number of samples: "<<GetNumberOfSamples()<<std::endl;
return true;
}
bool Kinematics::SaveConfig(const std::string& filename) { return true; }
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);
void Kinematics::Run() {
std::cout<<"Running simulation..."<<std::endl;
switch(m_rxn_type) {
case ONESTEP_DECAY:
{
RunOneStepDecay();
break;
}
case ONESTEP_RXN:
{
RunOneStepRxn();
break;
}
case TWOSTEP:
{
RunTwoStep();
break;
}
case THREESTEP:
{
RunThreeStep();
break;
}
}
std::cout<<std::endl;
std::cout<<"Complete."<<std::endl;
std::cout<<"---------------------------------------------"<<std::endl;
}
//For progress tracking
int percent5 = 0.05*m_nsamples;
int count = 0;
int npercent = 0;
std::vector<Nucleus> data;
data.resize(4);
for(int i=0; i<m_nsamples; i++) {
if(++count == percent5) {//Show update every 5 percent
npercent++;
count = 0;
std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
void Kinematics::RunOneStepRxn() {
OneStepSystem* this_sys = dynamic_cast<OneStepSystem*>(sys);
if(this_sys == nullptr) {
return;
}
this_sys->RunSystem();
MaskFile output(m_outfile_name, MaskFile::FileType::write);
std::vector<Nucleus> 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; i<m_nsamples; i++) {
if(++count == percent5) {//Show update every 5 percent
npercent++;
count = 0;
std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
}
this_sys->RunSystem();
data[0] = this_sys->GetTarget();
data[1] = this_sys->GetProjectile();
data[2] = this_sys->GetEjectile();
data[3] = this_sys->GetResidual();
output.WriteData(data);
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");
}
output.Close();
}
void Kinematics::RunOneStepDecay() {
DecaySystem* this_sys = dynamic_cast<DecaySystem*>(sys);
if(this_sys == nullptr) {
return;
}
output->cd();
if(save_tree_flag) tree->Write(tree->GetName(), TObject::kOverwrite);
if(do_plotter_flag) {
plotman.GetTable()->Write();
plotman.ClearTable();
MaskFile output(m_outfile_name, MaskFile::FileType::write);
std::vector<Nucleus> 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; i<m_nsamples; i++) {
if(++count == percent5) {//Show update every 5 percent
npercent++;
count = 0;
std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
}
this_sys->RunSystem();
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<TwoStepSystem*>(sys);
if(this_sys == nullptr) {
return;
}
MaskFile output(m_outfile_name, MaskFile::FileType::write);
std::vector<Nucleus> 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; i<m_nsamples; i++) {
if(++count == percent5) {//Show update every 5 percent
npercent++;
count = 0;
std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
}
this_sys->RunSystem();
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<ThreeStepSystem*>(sys);
if(this_sys == nullptr) {
return;
}
MaskFile output(m_outfile_name, MaskFile::FileType::write);
std::vector<Nucleus> 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; i<m_nsamples; i++) {
if(++count == percent5) {//Show update every 5 percent
npercent++;
count = 0;
std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
}
this_sys->RunSystem();
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->Close();
}
void Kinematics::RunOneStepDecay() {
DecaySystem* this_sys = dynamic_cast<DecaySystem*>(sys);
if(this_sys == nullptr) {
return;
}
TFile* output = TFile::Open(m_outfile_name.c_str(), "RECREATE");
TTree* tree;
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; i<m_nsamples; i++) {
if(++count == percent5) {//Show update every 5 percent
npercent++;
count = 0;
std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
}
this_sys->RunSystem();
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<TwoStepSystem*>(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);
}
//For progress tracking
int percent5 = 0.05*m_nsamples;
int count = 0;
int npercent = 0;
for(int i=0; i<m_nsamples; i++) {
if(++count == percent5) {//Show update every 5 percent
npercent++;
count = 0;
std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
}
this_sys->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");
}
}
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<ThreeStepSystem*>(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; i<m_nsamples; i++) {
if(++count == percent5) {//Show update every 5 percent
npercent++;
count = 0;
std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
}
this_sys->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());
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();
}
};

479
src/MaskFile.cpp
View File

@ -3,177 +3,336 @@
#include <iomanip>
#include <sstream>
/*
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!"<<std::endl;
return false;
MaskFile::MaskFile() :
file_type(MaskFile::FileType::none), filename(""), buffer_position(0), buffer_end(0), data_size(0), buffersize_bytes(0), file()
{
}
if(type == FileType::read) {
file.open(name, std::ios::in);
buffer_position = 0;
if(IsOpen()) ReadHeader();
} else if (type == FileType::write) {
file.open(name, std::ios::out | std::ios::trunc);
file<<std::setprecision(precision);
file<<std::scientific;
} else if (type == FileType::append) {
file.open(name, std::ios::out | std::ios::app);
file<<std::setprecision(precision);
file<<std::scientific;
} else {
std::cerr<<"Invalid FileType at MaskFile::Open()"<<std::endl;
MaskFile::MaskFile(const std::string& name, MaskFile::FileType type) :
file_type(type), filename(name), buffer_position(0), buffer_end(0), data_size(0), buffersize_bytes(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!"<<std::endl;
return false;
}
if(type == FileType::read) {
file.open(name, std::ios::in);
buffer_position = 0;
} else if (type == FileType::write) {
file.open(name, std::ios::out | std::ios::trunc);
} else if (type == FileType::append) {
file.open(name, std::ios::out | std::ios::app);
} else {
std::cerr<<"Invalid FileType at MaskFile::Open()"<<std::endl;
return IsOpen();
}
return IsOpen();
}
count = 0;
return IsOpen();
}
void MaskFile::Close() {
if(data_buffer.size() > 0) {
FlushBuffer();
}
file.close();
}
void MaskFile::WriteHeader(std::vector<Nucleus>& nuclei) {
file<<std::setw(width)<<"ZT"<<","
<<std::setw(width)<<"AT"<<","
<<std::setw(width)<<"ZP"<<","
<<std::setw(width)<<"AP"<<","
<<std::setw(width)<<"ZE"<<","
<<std::setw(width)<<"AE"<<","
<<std::setw(width)<<"ZR"<<","
<<std::setw(width)<<"AR"<<","
<<std::setw(width)<<"ZB1"<<","
<<std::setw(width)<<"AB1"<<","
<<std::setw(width)<<"ZB2"<<","
<<std::setw(width)<<"AB2"<<","
<<std::setw(width)<<"ZB3"<<","
<<std::setw(width)<<"AB3"<<","
<<std::setw(width)<<"ZB4"<<","
<<std::setw(width)<<"AB4"<<","
<<std::setw(width)<<"KEP"<<","
<<std::endl;
file<<std::setw(width)<<nuclei[0].GetZ()<<","
<<std::setw(width)<<nuclei[0].GetA()<<","
<<std::setw(width)<<nuclei[1].GetZ()<<","
<<std::setw(width)<<nuclei[1].GetA()<<","
<<std::setw(width)<<nuclei[2].GetZ()<<","
<<std::setw(width)<<nuclei[2].GetA()<<","
<<std::setw(width)<<nuclei[3].GetZ()<<","
<<std::setw(width)<<nuclei[3].GetA()<<",";
if(nuclei.size() == 6) {
file<<std::setw(width)<<nuclei[4].GetZ()<<","
<<std::setw(width)<<nuclei[4].GetA()<<","
<<std::setw(width)<<nuclei[5].GetZ()<<","
<<std::setw(width)<<nuclei[5].GetA()<<",";
} else if(nuclei.size() == 8) {
file<<std::setw(width)<<nuclei[4].GetZ()<<","
<<std::setw(width)<<nuclei[4].GetA()<<","
<<std::setw(width)<<nuclei[5].GetZ()<<","
<<std::setw(width)<<nuclei[5].GetA()<<","
<<std::setw(width)<<nuclei[6].GetZ()<<","
<<std::setw(width)<<nuclei[6].GetA()<<","
<<std::setw(width)<<nuclei[7].GetZ()<<","
<<std::setw(width)<<nuclei[7].GetA()<<",";
}
file<<std::setw(width)<<nuclei[1].GetKE()<<",";
file<<std::endl;
file<<std::setw(width)<<"Event"<<","
<<std::setw(width)<<"EE"<<","<<std::setw(width)<<"KEE"<<","<<std::setw(width)<<"PE"<<","<<std::setw(width)<<"ThetaE"<<","<<std::setw(width)<<"PhiE"<<","<<std::setw(width)<<"ThetaCME"<<","
<<std::setw(width)<<"ER"<<","<<std::setw(width)<<"KER"<<","<<std::setw(width)<<"PR"<<","<<std::setw(width)<<"ThetaR"<<","<<std::setw(width)<<"PhiR"<<","<<std::setw(width)<<"ThetaCMR"<<","
<<std::setw(width)<<"EB1"<<","<<std::setw(width)<<"KEB1"<<","<<std::setw(width)<<"PB1"<<","<<std::setw(width)<<"ThetaB1"<<","<<std::setw(width)<<"PhiB1"<<","<<std::setw(width)<<"ThetaCMB1"<<","
<<std::setw(width)<<"EB2"<<","<<std::setw(width)<<"KEB2"<<","<<std::setw(width)<<"PB2"<<","<<std::setw(width)<<"ThetaB2"<<","<<std::setw(width)<<"PhiB2"<<","<<std::setw(width)<<"ThetaCMB2"<<","
<<std::setw(width)<<"EB3"<<","<<std::setw(width)<<"KEB3"<<","<<std::setw(width)<<"PB3"<<","<<std::setw(width)<<"ThetaB3"<<","<<std::setw(width)<<"PhiB3"<<","<<std::setw(width)<<"ThetaCMB3"<<","
<<std::setw(width)<<"EB4"<<","<<std::setw(width)<<"KEB4"<<","<<std::setw(width)<<"PB4"<<","<<std::setw(width)<<"ThetaB4"<<","<<std::setw(width)<<"PhiB4"<<","<<std::setw(width)<<"ThetaCMB4"<<","
<<std::endl;
}
void MaskFile::ReadHeader() {
std::string junk;
std::getline(file, junk);
}
void MaskFile::WriteData(std::vector<Nucleus>& data) {
if(count == 0 && data_buffer.size() == 0) {
WriteHeader(data);
}
data_buffer.push_back(data);
if(data_buffer.size() == buffersize) {
FlushBuffer();
}
}
void MaskFile::FlushBuffer() {
for(auto& event : data_buffer) {
file<<std::setw(width)<<count<<",";
for(unsigned int i=2; i< event.size(); i++) {
Nucleus& nuc = event[i];
file<<std::setw(width)<<nuc.GetE()<<","
<<std::setw(width)<<nuc.GetKE()<<","
<<std::setw(width)<<nuc.GetP()<<","
<<std::setw(width)<<nuc.GetTheta()<<","
<<std::setw(width)<<nuc.GetPhi()<<","
<<std::setw(width)<<nuc.GetThetaCM()<<",";
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<<std::endl;
count++;
file.close();
}
data_buffer.clear();
}
std::vector<Nucleus> MaskFile::ReadData() {
if(buffer_position == data_buffer.size()) {
FillBuffer();
buffer_position = 0;
}
std::vector<Nucleus> data = data_buffer[buffer_position];
buffer_position++;
return data;
}
void MaskFile::FillBuffer() {
std::vector<Nucleus> 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);
void MaskFile::WriteHeader(int rxn_type, int nsamples) {
//size of a datum = data nuclei per event * (# doubles per nucleus * sizeof double + # ints per nucleus * sizeof int + sizeof bool)
if(rxn_type == 0) {
m_rxn_type = 0;
data_size = 3 * ( 5 * double_size + 2 * int_size + bool_size);
}
data_buffer.push_back(data);
data.clear();
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"<<std::endl;
return;
}
buffersize_bytes = buffersize * data_size; //buffersize_bytes = # of events * size of an event
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<char> 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 = "<<m_rxn_type<<")! Returning"<<std::endl;
return header;
}
buffersize_bytes = buffersize * data_size;//buffersize_bytes = size of a datum * # of events
header.rxn_type = m_rxn_type;
data_buffer.resize(buffersize_bytes);
return header;
}
void MaskFile::WriteData(std::vector<Nucleus>& data) {
char* data_pointer;
double datum;
int number;
bool flag;
std::size_t j;
for(unsigned int i=0; i<data.size(); i++) {
number = data[i].GetZ();
data_pointer = (char*) &number;
for(j=0; j<int_size; j++) {
data_buffer[buffer_position] = *(data_pointer + j);
buffer_position++;
}
number = data[i].GetA();
data_pointer = (char*) &number;
for(j=0; j<int_size; j++) {
data_buffer[buffer_position] = *(data_pointer + j);
buffer_position++;
}
flag = data[i].IsDetected();
data_pointer = (char*) &flag;
for(j=0; j<bool_size; j++) {
data_buffer[buffer_position] = *(data_pointer + j);
buffer_position++;
}
datum = data[i].GetE();
data_pointer = (char*) &datum;
for(j=0; j<double_size; j++) {
data_buffer[buffer_position] = *(data_pointer + j);
buffer_position++;
}
datum = data[i].GetKE();
data_pointer = (char*) &datum;
for(j=0; j<double_size; j++) {
data_buffer[buffer_position] = *(data_pointer + j);
buffer_position++;
}
datum = data[i].GetP();
data_pointer = (char*) &datum;
for(j=0; j<double_size; j++) {
data_buffer[buffer_position] = *(data_pointer + j);
buffer_position++;
}
datum = data[i].GetTheta();
data_pointer = (char*) &datum;
for(j=0; j<double_size; j++) {
data_buffer[buffer_position] = *(data_pointer + j);
buffer_position++;
}
datum = data[i].GetPhi();
data_pointer = (char*) &datum;
for(j=0; j<double_size; j++) {
data_buffer[buffer_position] = *(data_pointer + j);
buffer_position++;
}
}
//Flush the buffer when it is full, and reset the position.
if(buffer_position == buffersize_bytes) {
file.write(data_buffer.data(), data_buffer.size());
buffer_position = 0;
}
}
void MaskFile::WriteData(MaskFileData& data) {
char* data_pointer;
double datum;
int number;
bool flag;
std::size_t j;
for(unsigned int i=0; i<data.Z.size(); i++) {
number = data.Z[i];
data_pointer = (char*) &number;
for(j=0; j<int_size; j++) {
data_buffer[buffer_position] = *(data_pointer + j);
buffer_position++;
}
number = data.A[i];
data_pointer = (char*) &number;
for(j=0; j<int_size; j++) {
data_buffer[buffer_position] = *(data_pointer + j);
buffer_position++;
}
flag = data.detect_flag[i];
data_pointer = (char*) &flag;
for(j=0; j<bool_size; j++) {
data_buffer[buffer_position] = *(data_pointer + j);
buffer_position++;
}
datum = data.E[i];
data_pointer = (char*) &datum;
for(j=0; j<double_size; j++) {
data_buffer[buffer_position] = *(data_pointer + j);
buffer_position++;
}
datum = data.KE[i];
data_pointer = (char*) &datum;
for(j=0; j<double_size; j++) {
data_buffer[buffer_position] = *(data_pointer + j);
buffer_position++;
}
datum = data.p[i];
data_pointer = (char*) &datum;
for(j=0; j<double_size; j++) {
data_buffer[buffer_position] = *(data_pointer + j);
buffer_position++;
}
datum = data.theta[i];
data_pointer = (char*) &datum;
for(j=0; j<double_size; j++) {
data_buffer[buffer_position] = *(data_pointer + j);
buffer_position++;
}
datum = data.phi[i];
data_pointer = (char*) &datum;
for(j=0; j<double_size; j++) {
data_buffer[buffer_position] = *(data_pointer + j);
buffer_position++;
}
}
//Flush the buffer when it is full, and reset the position.
if(buffer_position == buffersize_bytes) {
file.write(data_buffer.data(), data_buffer.size());
buffer_position = 0;
}
}
/*
Read data from the buffer and submit it to the client side as a MaskFileData struct.
When file reaches the end of the file (no more data to read), an empty MaskFileData with
eof == true is sent out signaling that the file is finished.
Should be used like
Mask::MaskFile input(file, Mask::MaskFile::FileType::read);
Mask::MaskFileHeader header = input.ReadHeader();
Mask::MaskFileData data;
while(true) {
data = input.ReadData();
if(data.eof) break;
Do some stuff...
}
input.Close();
Good luck
*/
MaskFileData MaskFile::ReadData() {
MaskFileData data;
//Fill the buffer when needed, reset the positon, and find the end
if(buffer_position == data_buffer.size() || buffer_position == buffer_end) {
file.read(data_buffer.data(), buffersize_bytes);
buffer_position = 0;
buffer_end = file.gcount();
if(buffer_end == 0 && file.eof()) {
data.eof = true;
return data;
}
}
unsigned int local_end = buffer_position + data_size;
if(local_end > buffer_end) {
std::cerr<<"Attempting to read past end of file at MaskFile::ReadData! Returning empty"<<std::endl;
data.eof = true;
return data;
}
while(buffer_position < local_end) {
data.Z.push_back(*(int*)(&data_buffer[buffer_position]));
buffer_position += int_size;
data.A.push_back(*(int*)(&data_buffer[buffer_position]));
buffer_position += int_size;
data.detect_flag.push_back(*(bool*)(&data_buffer[buffer_position]));
buffer_position += bool_size;
data.E.push_back(*(double*)(&data_buffer[buffer_position]));
buffer_position += double_size;
data.KE.push_back(*(double*)(&data_buffer[buffer_position]));
buffer_position += double_size;
data.p.push_back(*(double*)(&data_buffer[buffer_position]));
buffer_position += double_size;
data.theta.push_back(*(double*)(&data_buffer[buffer_position]));
buffer_position += double_size;
data.phi.push_back(*(double*)(&data_buffer[buffer_position]));
buffer_position += double_size;
}
return data;
}
}
};

78
src/MassLookup.cpp
View File

@ -16,61 +16,47 @@ Written by G.W. McCann Aug. 2020
Read in AMDC mass file, preformated to remove excess info. Here assumes that by default
the file is in a local directory etc/
*/
MassLookup* MassLookup::s_instance = nullptr;
MassLookup::MassLookup() {
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();
}
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 "<<key<<" (Z,A) not found in Mass Table! Returning 1"<<std::endl;
return 1;
}*/
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;
}
//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: "<<Z<<" not found in Element Table! Returning void."<<std::endl;
std::string element = "void";
return element;
}*/
auto data = elementTable.find(Z);
if(data == elementTable.end()) {
throw MassException();
}
std::string fullsymbol = std::to_string(A) + data->second;
return fullsymbol;
}

16
src/Nucleus.cpp
View File

@ -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;
}

128
src/OneStepSystem.cpp
View File

@ -2,76 +2,76 @@
namespace Mask {
OneStepSystem::OneStepSystem() :
ReactionSystem()
{
}
OneStepSystem::OneStepSystem(std::vector<int>& z, std::vector<int>& a) :
ReactionSystem()
{
SetNuclei(z, a);
}
OneStepSystem::~OneStepSystem() {}
bool OneStepSystem::SetNuclei(std::vector<int>& z, std::vector<int>& a) {
if(z.size() != a.size() || z.size() != 3) {
return false;
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<int>& z, std::vector<int>& 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);
OneStepSystem::~OneStepSystem() {}
step1.SetBeamKE(bke);
step1.SetPolarRxnAngle(rxnTheta);
step1.SetAzimRxnAngle(rxnPhi);
step1.SetExcitation(residEx);
bool OneStepSystem::SetNuclei(std::vector<int>& z, std::vector<int>& a) {
if(z.size() != a.size() || z.size() != 3) {
return false;
}
step1.TurnOnResidualEloss();
step1.Calculate();
} else {
return;
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;
}
}
}
}

89
src/Plotter.cpp
View File

@ -1,89 +0,0 @@
#include "Plotter.h"
#include <iostream>
namespace Mask {
Plotter::Plotter() :
table(new THashTable())
{
}
Plotter::~Plotter() {
for(unsigned int i=0; i<garbage_collection.size(); i++) {
delete garbage_collection[i];
}
garbage_collection.clear();
delete table;
}
void Plotter::FillData(const 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 Plotter::MyFill(const char* name, const char* title, int bins, float min, float max, double val) {
TH1F* h = (TH1F*) table->FindObject(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);
}
}
};

137
src/Plotters/ROOT/RootPlotter.cpp Normal file
View File

@ -0,0 +1,137 @@
#include "RootPlotter.h"
#include "MaskFile.h"
#include <TFile.h>
#include <iostream>
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: "<<argc<<" Exiting."<<std::endl;
return 1;
}
std::string inputname = argv[1];
std::string outputname = argv[2];
Mask::MaskFile input(inputname, Mask::MaskFile::FileType::read);
TFile* root_out = TFile::Open(outputname.c_str(), "RECREATE");
RootPlotter plotter;
Mask::MaskFileHeader header = input.ReadHeader();
std::cout<<"File Header -- rxn type: "<<header.rxn_type<<" nsamples: "<<header.nsamples<<std::endl;
Mask::MaskFileData data;
Mask::Nucleus nucleus;
double flush_frac = 0.05;
int count = 0, flush_val = flush_frac*header.nsamples, flush_count = 0;
while(true) {
if(count == flush_val) {
count = 0;
flush_count++;
std::cout<<"\rPercent of file processed: "<<flush_frac*flush_count*100<<"%"<<std::flush;
}
data = input.ReadData();
if(data.eof)
break;
for(unsigned int i=0; i<data.Z.size(); i++) {
nucleus.SetIsotope(data.Z[i], data.A[i]);
nucleus.SetVectorSpherical(data.theta[i], data.phi[i], data.p[i], data.E[i]);
plotter.FillData(nucleus);
if(data.detect_flag[i] == true) {
plotter.FillData(nucleus, "detected");
}
}
count++;
}
std::cout<<std::endl;
input.Close();
root_out->cd();
plotter.GetTable()->Write();
root_out->Close();
return 0;
}

38
src/ReactionSystem.cpp
View File

@ -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<int>& zt, std::vector<int>& at, std::vector<int>& 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<int>& zt, std::vector<int>& at, std::vector<int>& stoich, double thickness) {
target.AddLayer(zt, at, stoich, thickness);
}
};

31
src/Stopwatch.cpp Normal file
View File

@ -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<std::chrono::duration<double>>(stop_time-start_time).count();
}
double Stopwatch::GetElapsedMilliseconds() {
return std::chrono::duration_cast<std::chrono::duration<double>>(stop_time-start_time).count()*1000.0;
}

240
src/ThreeStepSystem.cpp
View File

@ -3,131 +3,129 @@
namespace Mask {
ThreeStepSystem::ThreeStepSystem() :
ReactionSystem()
{
}
ThreeStepSystem::ThreeStepSystem(std::vector<int>& z, std::vector<int>& 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<int>&z, std::vector<int>& 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();
ThreeStepSystem::ThreeStepSystem() :
ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
{
}
//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;
ThreeStepSystem::ThreeStepSystem(std::vector<int>& z, std::vector<int>& a) :
ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
{
SetNuclei(z, a);
}
step2.Calculate();
step3.SetTarget(step2.GetResidual());
if(decay2costheta == -10) {
step3.ResetEjectile();
step3.ResetResidual();
return;
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<int>&z, std::vector<int>& 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();
}
step3.TurnOnResidualEloss();
step3.Calculate();
}
};

192
src/TwoStepSystem.cpp
View File

@ -3,105 +3,105 @@
namespace Mask {
TwoStepSystem::TwoStepSystem() :
ReactionSystem()
{
}
TwoStepSystem::TwoStepSystem(std::vector<int>& z, std::vector<int>& 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<int>&z, std::vector<int>& 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();
TwoStepSystem::TwoStepSystem() :
ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
{
}
//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(std::vector<int>& z, std::vector<int>& a) :
ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
{
SetNuclei(z, a);
}
step2.TurnOnResidualEloss();
step2.Calculate();
TwoStepSystem::~TwoStepSystem() {
}
void TwoStepSystem::SetRandomGenerator(std::mt19937* gen) {
generator = gen;
decay1dist.AttachRandomNumberGenerator(gen);
gen_set_flag = true;
}
bool TwoStepSystem::SetNuclei(std::vector<int>&z, std::vector<int>& a) {
if(z.size() != a.size() || z.size() != 4) {
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();
}
}
};

43
src/main.cpp
View File

@ -1,21 +1,20 @@
#include <iostream>
#include <string>
#include "Stopwatch.h"
#include "MaskFile.h"
#include "Kinematics.h"
#include "SabreEfficiency.h"
#include "AnasenEfficiency.h"
#include "KinematicsExceptions.h"
#include <TGraph2D.h>
#include <TCanvas.h>
int main(int argc, char** argv) {
if(argc<2) {
std::cerr<<"Incorrect number of arguments!"<<std::endl;
return 1;
}
Stopwatch sw;
Mask::Kinematics calculator;
sw.Start();
try {
if(!calculator.LoadConfig(argv[1])) {
return 1;
@ -26,26 +25,26 @@ int main(int argc, char** argv) {
std::cerr<<"Terminating process."<<std::endl;
return 1;
}
sw.Stop();
std::cout<<"Time elapsed(seconds): "<<sw.GetElapsedSeconds()<<std::endl;
Mask::MaskFile input(calculator.GetOutputName(), Mask::MaskFile::FileType::read);
Mask::MaskFileData data;
Mask::MaskFileHeader header = input.ReadHeader();
std::cout<<"Header Found -- rxn type: "<<header.rxn_type<<" nsamples: "<<header.nsamples;
std::cout<<std::endl;
/*
SabreEfficiency sabre;
std::string mapfile = "./etc/DeadChannels.txt";
sabre.SetReactionType(calculator.GetReactionType());
sabre.SetDeadChannelMap(mapfile);
sabre.CalculateEfficiency(calculator.GetOutputName());
//std::cout<<"Running consistency check(1=success): "<<sabre.RunConsistencyCheck()<<std::endl;
//sabre.DrawDetectorSystem("/data1/gwm17/10B3He/Feb2021/simulation/SABREGeo.root");
*/
AnasenEfficiency anasen;
anasen.SetReactionType(calculator.GetReactionType());
anasen.CalculateEfficiency(calculator.GetOutputName());
//std::cout<<"Running consistency check(1=success): "<<anasen.RunConsistencyCheck()<<std::endl;
//anasen.DrawDetectorSystem("/data1/gwm17/TRIUMF_7Bed/simulation/ANASENGeo_centered_target_targetGap_BackQQQ_test.root");
int counter=0;
while(!data.eof) {
data = input.ReadData();
for(unsigned int i=0; i<data.E.size(); i++)
counter++;
//std::cout<<"Data Found -- E: "<<data.E[i]<<" KE: "<<data.KE[i]<<" p: "<<data.p[i]<<" theta: "<<data.theta[i]<<" phi: "<<data.phi[i]<<std::endl;
}
std::cout<<"events found: "<<counter<<std::endl;
input.Close();
std::cout<<"File closed."<<std::endl;
return 0;