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

Update project structure. Tighten up code, improve syntax, enforce style. Still bugged, but now working from proper environment. Switch back to ROOT as file I/O and data management

Browse Source
This commit is contained in:
Gordon McCann 2022-08-18 10:32:48 -04:00
parent 66900ea9b4
commit 0b7b06e4f1
74 changed files with 2644 additions and 2596 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
CMakeLists.txt
View File

@ -8,8 +8,10 @@ set(MASK_INCLUDE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/include)
set(CMAKE_CXX_STANDARD 17)
find_package(ROOT REQUIRED COMPONENTS GenVector)
add_subdirectory(src/vendor/catima)
add_subdirectory(src/Mask)
add_subdirectory(src/MaskApp)
add_subdirectory(src/Detectors)
add_subdirectory(src/Plotters/ROOT)
add_subdirectory(src/Plotters)

69
include/AnasenEfficiency.h
View File

@ -1,69 +0,0 @@
#ifndef ANASEN_EFFICIENCY_H
#define ANASEN_EFFICIENCY_H
#include <string>
#include "DetectorEfficiency.h"
#include "StripDetector.h"
#include "QQQDetector.h"
#include "Target.h"
#include "Nucleus.h"
#include "MaskFile.h"
#include "AnasenDeadChannelMap.h"
struct DetectorResult {
bool detectFlag = false;
Mask::Vec3 direction;
double energy_deposited = 0.0;
std::string det_name = "";
};
class AnasenEfficiency : public DetectorEfficiency {
public:
AnasenEfficiency();
~AnasenEfficiency();
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;
inline void SetDeadChannelMap(const std::string& filename) { dmap.LoadMapfile(filename); }
private:
DetectorResult IsRing1(Mask::Nucleus& nucleus);
DetectorResult IsRing2(Mask::Nucleus& nucleus);
DetectorResult IsQQQ(Mask::Nucleus& nucleus);
DetectorResult IsAnasen(Mask::Nucleus& nucleus);
void CountCoincidences(const Mask::MaskFileData& data, std::vector<int>& counts, Mask::RxnType rxn_type);
std::vector<StripDetector> m_Ring1, m_Ring2;
std::vector<QQQDetector> m_forwardQQQs;
std::vector<QQQDetector> m_backwardQQQs;
Mask::Target det_silicon;
AnasenDeadChannelMap dmap;
/**** ANASEN geometry constants *****/
const int n_sx3_per_ring = 12;
const int n_qqq = 4;
const double sx3_length = 0.075;
const double sx3_width = 0.04;
const double barrel_gap = 0.0254;
const double sx3_frame = 0.049; //0.049 is base gap due to frames
const double ring1_z = sx3_length/2.0 + sx3_frame + barrel_gap/2.0;
const double ring2_z = (-1.0)*(barrel_gap/2.0 + sx3_length/2.0);
const double qqq_nom_z = 0.0125 + sx3_length + sx3_frame + barrel_gap/2.0;
const double qqq_rinner = 0.0501;
const double qqq_router = 0.0990;
const double qqq_deltaphi = 1.52119;
const double qqq_z[4] = {qqq_nom_z, qqq_nom_z, qqq_nom_z, qqq_nom_z};
const double qqq_phi[4] = {5.49779, 0.785398, 2.35619, 3.92699};
const double ring_rho[12] = {0.0890601, 0.0889871, 0.0890354, 0.0890247, 0.0890354, 0.0890354, 0.0890247, 0.0890354, 0.0890354, 0.0890247, 0.0890354, 0.0890354};
const double ring_phi[12] = {4.97426, 5.49739, 6.02132, 0.261868, 0.785398, 1.30893, 1.83266, 2.35619, 2.87972, 3.40346, 3.92699, 4.45052};
/*************************/
static constexpr double threshold = 0.6; //MeV
static constexpr double deg2rad = M_PI/180.0;
static constexpr double si_thickness = 1000 * 1e-4 * 2.3926 * 1e6; //thickness in um -> eff thickness in ug/cm^2 for detector
};
#endif

34
include/AngularDistribution.h
View File

@ -1,34 +0,0 @@
#ifndef ANGULARDISTRIBUTION_H
#define ANGULARDISTRIBUTION_H
#include <string>
#include <vector>
#include <random>
namespace Mask {
class AngularDistribution {
public:
AngularDistribution();
AngularDistribution(const std::string& file);
~AngularDistribution();
void ReadDistributionFile(const std::string& file);
double GetRandomCosTheta();
inline int GetL() { return L; }
inline double GetBranchingRatio() { return branchingRatio; }
private:
inline bool IsIsotropic() { return isoFlag; }
std::uniform_real_distribution<double> uniform_cosine_dist;
std::uniform_real_distribution<double> uniform_prob_dist;
double branchingRatio;
int L;
std::vector<double> constants;
bool isoFlag;
};
}
#endif

41
include/DecaySystem.h
View File

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

38
include/MaskApp.h
View File

@ -1,38 +0,0 @@
#ifndef MASKAPP_H
#define MASKAPP_H
#include "ReactionSystem.h"
#include "DecaySystem.h"
#include "OneStepSystem.h"
#include "TwoStepSystem.h"
#include "ThreeStepSystem.h"
#include "RxnType.h"
namespace Mask {
class MaskApp {
public:
MaskApp();
~MaskApp();
bool LoadConfig(const std::string& filename);
bool SaveConfig(const std::string& filename);
inline int GetNumberOfSamples() const { return m_nsamples; };
inline const std::string GetSystemName() const { return m_sys == nullptr ? "" : m_sys->GetSystemEquation(); };
inline const std::string GetOutputName() const { return m_outfile_name; };
inline const RxnType GetReactionType() const { return m_rxn_type; };
void Run();
private:
ReactionSystem* m_sys;
std::string m_outfile_name;
RxnType m_rxn_type;
int m_nsamples;
};
}
#endif

64
include/Nucleus.h
View File

@ -1,64 +0,0 @@
/*
Nucleus.h
Nucleus is a derived class of Vec4. A nucleus is the kinematics is essentially a 4 vector with the
additional properties of the number of total nucleons (A), the number of protons (Z), a ground state mass,
an exctitation energy, and an isotopic symbol.
--GWM Jan 2021
*/
#ifndef NUCLEUS_H
#define NUCLEUS_H
#include "Vec4.h"
#include <string>
namespace Mask {
class Nucleus : public Vec4 {
public:
Nucleus();
Nucleus(int Z, int A);
Nucleus(int Z, int A, double px, double py, double pz, double E);
virtual ~Nucleus() override;
bool SetIsotope(int Z, int A);
inline void SetThetaCM(double tcm) { m_theta_cm = tcm; }; //save theta in rxn CM frame
inline int GetZ() const { return m_z; };
inline int GetA() const { return m_a; };
inline double GetExcitationEnergy() const { return GetInvMass() - m_gs_mass; };
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() const { return m_detectFlag; };
inline Nucleus& operator=(const Nucleus& rhs) {
SetIsotope(rhs.GetZ(), rhs.GetA());
SetVectorCartesian(rhs.GetPx(), rhs.GetPy(), rhs.GetPz(), rhs.GetE());
SetThetaCM(rhs.GetThetaCM());
return *this;
};
//Conservation of nucleons and momentum
inline Nucleus operator+(const Nucleus& daughter) {
return Nucleus(GetZ()+daughter.GetZ(), GetA()+daughter.GetA(), GetPx()+daughter.GetPx(), GetPy()+daughter.GetPy(), GetPz()+daughter.GetPz(), GetE()+daughter.GetE());
};
inline Nucleus operator-(const Nucleus& daughter) {
return (GetZ() - daughter.GetZ()) <= 0 || (GetA() - daughter.GetA()) <= 0 ? Nucleus() :
Nucleus(GetZ()-daughter.GetZ(), GetA() - daughter.GetA(), GetPx()-daughter.GetPx(), GetPy()-daughter.GetPy(), GetPz()-daughter.GetPz(), GetE()-daughter.GetE());
};
private:
int m_z, m_a;
double m_gs_mass;
double m_theta_cm;
std::string m_symbol;
bool m_detectFlag;
};
};
#endif

34
include/OneStepSystem.h
View File

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

57
include/QQQDetector.h
View File

@ -1,57 +0,0 @@
/*
QQQDetector.h
Class implementing geometry for QQQDetector where the detector is perpendicular to the beam axis.
Detector is first generated centered on the x-axis (phi=0)
Coordinate convention : +z is downstream, -z is upstream. +y is vertically down in the lab.
*/
#ifndef QQQDETECTOR_H
#define QQQDETECTOR_H
#include <cmath>
#include <vector>
#include <random>
#include "Vec3.h"
#include "Rotation.h"
#include "RandomGenerator.h"
class QQQDetector {
public:
QQQDetector(double R_in, double R_out, double deltaPhi, double phiCentral, double z, double x=0, double y=0);
~QQQDetector();
inline Mask::Vec3 GetRingCoordinates(int ringch, int corner) { return m_ringCoords[ringch][corner]; }
inline Mask::Vec3 GetWedgeCoordinates(int wedgech, int corner) { return m_wedgeCoords[wedgech][corner]; }
inline Mask::Vec3 GetNorm() { return m_norm; }
Mask::Vec3 GetTrajectoryCoordinates(double theta, double phi);
std::pair<int, int> GetTrajectoryRingWedge(double theta, double phi);
Mask::Vec3 GetHitCoordinates(int ringch, int wedgech);
inline void TurnOnRandomizedCoordinates() { rndmFlag = true; }
inline void TurnOffRandomizedCoordinates() { rndmFlag = false; }
inline int GetNumberOfRings() { return nrings; }
inline int GetNumberOfWedges() { return nwedges; }
private:
inline bool CheckChannel(int ch) { return (ch >=0 && ch < nrings); }
inline bool CheckCorner(int corner) { return (corner >=0 && corner < 4); }
void CalculateCorners();
Mask::Vec3 TransformCoordinates(Mask::Vec3& vector) { return m_ZRot*vector + m_translation; }
double m_Rinner, m_Router, m_deltaR, m_deltaPhi, m_deltaPhi_per_wedge, m_phiCentral;
std::vector<std::vector<Mask::Vec3>> m_ringCoords, m_wedgeCoords;
Mask::Vec3 m_translation;
Mask::Vec3 m_norm;
Mask::ZRotation m_ZRot;
std::uniform_real_distribution<double> m_uniform_fraction;
bool rndmFlag;
static constexpr int nrings = 16;
static constexpr int nwedges = 16;
static constexpr double deg2rad = M_PI/180.0;
};
#endif

25
include/RandomGenerator.h
View File

@ -1,25 +0,0 @@
#ifndef RANDOMGENERATOR_H
#define RANDOMGENERATOR_H
#include <random>
namespace Mask {
class RandomGenerator {
public:
~RandomGenerator();
inline std::mt19937& GetGenerator() { return rng; }
inline static RandomGenerator& GetInstance() {
static RandomGenerator s_instance;
return s_instance;
}
private:
RandomGenerator();
std::mt19937 rng;
};
}
#endif

78
include/Reaction.h
View File

@ -1,78 +0,0 @@
/*
Reaction.h
Reaction is a class which implements either a decay or scattering reaction. As such it requires either
3 (decay) or 4 (scattering) nuclei to perform any calcualtions. I also links together the target, which provides
energy loss calculations, with the kinematics. Note that Reaction does not own the LayeredTarget.
--GWM Jan. 2021
*/
#ifndef REACTION_H
#define REACTION_H
#include "Nucleus.h"
#include "LayeredTarget.h"
namespace Mask {
class Reaction {
public:
Reaction();
Reaction(int zt, int at, int zp, int ap, int ze, int ae);
~Reaction();
bool Calculate();
void SetNuclei(int zt, int at, int zp, int ap, int ze, int ae);
void SetNuclei(const Nucleus* nucs);
void SetBeamKE(double bke);
void SetEjectileThetaType(int type);
inline void SetLayeredTarget(LayeredTarget* targ) { target = targ; };
inline void SetPolarRxnAngle(double theta) { m_theta = theta; };
inline void SetAzimRxnAngle(double phi) { m_phi = phi; };
inline void SetExcitation(double ex) { m_ex = ex; };
inline void SetTarget(const Nucleus& nuc) { reactants[0] = nuc; };
inline void SetTarget(int z, int a) { reactants[0] = Nucleus(z, a); };
inline void SetProjectile(const Nucleus& nuc) { reactants[1] = nuc; };
inline void SetProjectile(int z, int a) { reactants[1] = Nucleus(z, a); };
inline void SetEjectile(const Nucleus& nuc) { reactants[2] = nuc; };
inline void SetEjectile(int z, int a) { reactants[2] = Nucleus(z, a); };
inline void SetResidual(const Nucleus& nuc) { reactants[3] = nuc; };
inline void SetResidual(int z, int a) { reactants[3] = Nucleus(z, a); };
inline void SetRxnLayer(int layer) { rxnLayer = layer; };
inline void TurnOffResidualEloss() { resid_elossFlag = false; };
inline void TurnOnResidualEloss() { resid_elossFlag = true; };
inline bool IsDecay() { return decayFlag; };
inline const Nucleus* GetNuclei() const { return &(reactants[0]); };
inline const Nucleus& GetProjectile() const { return reactants[1]; };
inline const Nucleus& GetTarget() const { return reactants[0]; };
inline const Nucleus& GetEjectile() const { return reactants[2]; };
inline const Nucleus& GetResidual() const { return reactants[3]; };
inline int GetRxnLayer() { return rxnLayer; };
inline void ResetTarget() { reactants[0].SetVectorCartesian(0,0,0, reactants[0].GetGroundStateMass()); }
inline void ResetProjectile() { reactants[1].SetVectorCartesian(0,0,0, reactants[1].GetGroundStateMass()); }
inline void ResetEjectile() { reactants[2].SetVectorCartesian(0,0,0, reactants[2].GetGroundStateMass()); }
inline void ResetResidual() { reactants[3].SetVectorCartesian(0,0,0, reactants[3].GetGroundStateMass()); }
private:
void CalculateDecay(); //target -> light_decay (eject) + heavy_decay(resid)
void CalculateReaction(); //target + project -> eject + resid
void CalculateReactionThetaLab();
void CalculateReactionThetaCM();
Nucleus reactants[4]; //0=target, 1=projectile, 2=ejectile, 3=residual
LayeredTarget* target; //not owned by Reaction
double m_bke, m_theta, m_phi, m_ex;
int rxnLayer;
int m_eject_theta_type;
bool decayFlag, nuc_initFlag, resid_elossFlag;
static constexpr int lab = 0;
static constexpr int center_of_mass = 1;
};
}
#endif

52
include/RootPlotter.h
View File

@ -1,52 +0,0 @@
#ifndef ROOTPLOTTER_H
#define ROOTPLOTTER_H
#include <vector>
#include <string>
#include "Nucleus.h"
#include "RxnType.h"
#include "MaskFile.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, double detKE = 0.0, const std::string& modifier = "");
void FillCorrelations(const Mask::MaskFileData& data, Mask::RxnType type);
void FillCorrelationsDetected(const Mask::MaskFileData& data, Mask::RxnType type);
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

177
include/SabreDetector.h
View File

@ -1,177 +0,0 @@
/*
Class which represents a single MMM detector in the SABRE array at FSU. Origial code by KGH, re-written by
GWM.
Distances in meters, angles in radians.
The channel arrays have four points, one for each corner. The corners are
as follows, as if looking BACK along beam (i.e. from the target's pov):
0---------------------1
| |
| | x
| | <-----
| | |
| | |
3---------------------2 y
(z is hence positive along beam direction)
The channel numbers, also as looking back from target pov, are:
>> rings are 0 -- 15 from inner to outer:
15 -------------------
14 -------------------
13 -------------------
.
.
.
2 -------------------
1 -------------------
0 -------------------
>> wedges are 0 -- 7 moving counterclockwise:
7 6 ... 1 0
| | | | | |
| | | | | |
| | | | | |
| | | | | |
| | | | | |
| | | | | |
>> Note that the detector starts centered on the x-axis (central phi = 0) untilted, and then is rotated to wherever the frick
it is supposed to go; phi = 90 is centered on y axis, pointing down towards the bottom of the scattering chamber
-- GWM, Dec 2020; based on the og code from kgh
*/
#ifndef SABREDETECTOR_H
#define SABREDETECTOR_H
#include <vector>
#include <cmath>
#include "Vec3.h"
#include "Rotation.h"
class SabreDetector {
public:
SabreDetector();
SabreDetector(int detID, double Rin, double Rout, double deltaPhi_flat, double phiCentral, double tiltFromVert, double zdist, double xdist=0, double ydist=0);
~SabreDetector();
/*Return coordinates of the corners of each ring/wedge in SABRE*/
inline Mask::Vec3 GetRingFlatCoords(int ch, int corner) { return CheckRingLocation(ch, corner) ? m_ringCoords_flat[ch][corner] : Mask::Vec3(); };
inline Mask::Vec3 GetWedgeFlatCoords(int ch, int corner) { return CheckWedgeLocation(ch, corner) ? m_wedgeCoords_flat[ch][corner] : Mask::Vec3(); };
inline Mask::Vec3 GetRingTiltCoords(int ch, int corner) { return CheckRingLocation(ch, corner) ? m_ringCoords_tilt[ch][corner] : Mask::Vec3(); };
inline Mask::Vec3 GetWedgeTiltCoords(int ch, int corner) { return CheckWedgeLocation(ch, corner) ? m_wedgeCoords_tilt[ch][corner] : Mask::Vec3(); };
Mask::Vec3 GetTrajectoryCoordinates(double theta, double phi);
std::pair<int, int> GetTrajectoryRingWedge(double theta, double phi);
Mask::Vec3 GetHitCoordinates(int ringch, int wedgech);
/*Basic getters*/
inline int GetNumberOfWedges() { return m_nWedges; };
inline int GetNumberOfRings() { return m_nRings; };
inline double GetInnerRadius() { return m_Rinner; };
inline double GetOuterRadius() { return m_Router; };
inline double GetPhiCentral() { return m_phiCentral; };
inline double GetTiltAngle() { return m_tilt; };
inline Mask::Vec3 GetTranslation() { return m_translation; };
inline Mask::Vec3 GetNormTilted() { return TransformToTiltedFrame(m_norm_flat); };
int GetDetectorID() { return m_detectorID; }
private:
/*Class constants*/
static constexpr int m_nRings = 16;
static constexpr int m_nWedges = 8;
static constexpr double deg2rad = M_PI/180.0;
/*These are implicitly the width of the spacing between detector active strips*/
static constexpr double POSITION_TOL = 0.0001; //0.1 mm position tolerance
static constexpr double ANGULAR_TOL = 0.1*M_PI/180.0; // 0.1 degree angular tolerance
void CalculateCorners();
/*Performs the transformation to the tilted,rotated,translated frame of the SABRE detector*/
inline Mask::Vec3 TransformToTiltedFrame(Mask::Vec3& vector) { return (m_ZRot*(m_YRot*vector)) + m_translation; };
/*Determine if a given channel/corner combo is valid*/
inline bool CheckRingChannel(int ch) { return (ch<m_nRings && ch>=0) ? true : false; };
inline bool CheckWedgeChannel(int ch) { return (ch<m_nWedges && ch >=0) ? true : false; };
inline bool CheckCorner(int corner) { return (corner < 4 && corner >=0) ? true : false; };
inline bool CheckRingLocation(int ch, int corner) { return CheckRingChannel(ch) && CheckCorner(corner); };
inline bool CheckWedgeLocation(int ch, int corner) { return CheckWedgeChannel(ch) && CheckCorner(corner); };
/*
For all of the calculations, need a limit precision to determine if values are actually equal or not
Here the approx. size of the strip spacing is used as the precision.
*/
inline bool CheckPositionEqual(double val1,double val2) { return fabs(val1-val2) > POSITION_TOL ? false : true; };
inline bool CheckAngleEqual(double val1,double val2) { return fabs(val1-val2) > ANGULAR_TOL ? false : true; };
/*Determine if a hit is within the bulk detector*/
inline bool IsInside(double r, double phi) {
double phi_1 = m_deltaPhi_flat/2.0;
double phi_2 = M_PI*2.0 - m_deltaPhi_flat/2.0;
return (((r > m_Rinner && r < m_Router) || CheckPositionEqual(r, m_Rinner) || CheckPositionEqual(r, m_Router)) && (phi > phi_2 || phi < phi_1 || CheckAngleEqual(phi, phi_1) || CheckAngleEqual(phi, phi_2)));
};
/*
For a given radius/phi are you inside of a given ring/wedge channel,
or are you on the spacing between these channels
*/
inline bool IsRing(double r, int ringch) {
double ringtop = m_Rinner + m_deltaR_flat_ring*(ringch + 1);
double ringbottom = m_Rinner + m_deltaR_flat_ring*(ringch);
return (r>ringbottom && r<ringtop);
};
inline bool IsRingTopEdge(double r, int ringch) {
double ringtop = m_Rinner + m_deltaR_flat_ring*(ringch + 1);
return CheckPositionEqual(r, ringtop);
};
inline bool IsRingBottomEdge(double r, int ringch) {
double ringbottom = m_Rinner + m_deltaR_flat_ring*(ringch);
return CheckPositionEqual(r, ringbottom);
};
inline bool IsWedge(double phi, int wedgech) {
double wedgetop = -m_deltaPhi_flat/2.0 + m_deltaPhi_flat_wedge*(wedgech+1);
double wedgebottom = -m_deltaPhi_flat/2.0 + m_deltaPhi_flat_wedge*(wedgech);
return ((phi>wedgebottom && phi<wedgetop));
};
inline bool IsWedgeTopEdge(double phi, int wedgech) {
double wedgetop = -m_deltaPhi_flat/2.0 + m_deltaPhi_flat_wedge*(wedgech+1);
return CheckAngleEqual(phi, wedgetop);
}
inline bool IsWedgeBottomEdge(double phi, int wedgech) {
double wedgebottom = -m_deltaPhi_flat/2.0 + m_deltaPhi_flat_wedge*(wedgech);
return CheckAngleEqual(phi, wedgebottom);
}
/*Class data*/
double m_Router, m_Rinner, m_deltaPhi_flat, m_phiCentral, m_tilt;
Mask::Vec3 m_translation;
Mask::YRotation m_YRot;
Mask::ZRotation m_ZRot;
double m_deltaR_flat, m_deltaR_flat_ring, m_deltaPhi_flat_wedge;
Mask::Vec3 m_norm_flat;
int m_detectorID;
std::vector<std::vector<Mask::Vec3>> m_ringCoords_flat, m_wedgeCoords_flat;
std::vector<std::vector<Mask::Vec3>> m_ringCoords_tilt, m_wedgeCoords_tilt;
};
#endif

51
include/SabreEfficiency.h
View File

@ -1,51 +0,0 @@
#ifndef SABREEFFICIENCY_H
#define SABREEFFICIENCY_H
#include "MaskFile.h"
#include "DetectorEfficiency.h"
#include "SabreDetector.h"
#include "Target.h"
#include "SabreDeadChannelMap.h"
#include "Nucleus.h"
class SabreEfficiency : public DetectorEfficiency {
public:
SabreEfficiency();
~SabreEfficiency();
void SetDeadChannelMap(std::string& filename) { dmap.LoadMapfile(filename); };
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::Nucleus& nucleus);
void CountCoincidences(const Mask::MaskFileData& data, std::vector<int>& counts, Mask::RxnType rxn_type);
std::vector<SabreDetector> detectors;
Mask::Target deadlayer;
Mask::Target sabre_eloss;
Mask::Target degrader;
SabreDeadChannelMap dmap;
//Sabre constants
const double INNER_R = 0.0326;
const double OUTER_R = 0.1351;
const double TILT = 40.0;
const double DIST_2_TARG = -0.1245;
const double PHI_COVERAGE = 54.4; //delta phi for each det
const double PHI0 = 306.0; //center phi values for each det in array
const double PHI1 = 18.0; //# is equal to detID in channel map
const double PHI2 = 234.0;
const double PHI3 = 162.0;
const double PHI4 = 90.0;
const double DEG2RAD = M_PI/180.0;
static constexpr double DEADLAYER_THIN = 50 * 1e-7 * 2.3296 * 1e6; // ug/cm^2 (50 nm thick * density)
static constexpr double SABRE_THICKNESS = 500 * 1e-4 * 2.3926 * 1e6; // ug/cm^2 (500 um thick * density)
static constexpr double DEGRADER_THICKNESS = 70.0 * 1.0e-4 * 16.69 * 1e6; //tantalum degrader (70 um thick)
const double ENERGY_THRESHOLD = 0.25; //in MeV
};
#endif

30
include/Stopwatch.h
View File

@ -1,30 +0,0 @@
/*
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

76
include/StripDetector.h
View File

@ -1,76 +0,0 @@
#ifndef __STRIP_DETECTOR_H
#define __STRIP_DETECTOR_H
// +z is along beam axis
// +y is vertically "downward" in the lab frame
//angles must be in radians, but distances can be whatever
//PROVIDED all input distances are the same
//Front strips from largest y to smallest y
//Back strips from lowest z to highest z
#include <cmath>
#include <vector>
#include <random>
#include "Vec3.h"
#include "Rotation.h"
#include "RandomGenerator.h"
struct StripHit
{
int front_strip_index=-1;
int back_strip_index=-1;
double front_ratio=0.0;
};
class StripDetector {
public:
StripDetector(int ns, double len, double wid, double cphi, double cz, double crho);
~StripDetector();
inline Mask::Vec3 GetFrontStripCoordinates(int stripch, int corner) { return front_strip_coords[stripch][corner]; }
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 Mask::Vec3 GetNormRotated() { return zRot*m_norm_unrot; }
inline void TurnOnRandomizedCoordinates() { rndmFlag = true; }
inline void TurnOffRandomizedCoordinates() { rndmFlag = false; }
Mask::Vec3 GetHitCoordinates(int front_stripch, double front_strip_ratio);
StripHit GetChannelRatio(double theta, double phi);
private:
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); };
void CalculateCorners();
int num_strips;
static constexpr int num_corners = 4;
double length; //common to all strips, hence total
double total_width;
double front_strip_width; //assuming equal widths
double back_strip_length; //assuming equal widths
double center_phi; //assuming det centered above x-axis (corresponds to zero phi)
double center_z;
double center_rho; //perpendicular radius from axis
std::vector<std::vector<Mask::Vec3>> front_strip_coords, back_strip_coords;
std::vector<std::vector<Mask::Vec3>> rotated_front_strip_coords, rotated_back_strip_coords;
Mask::Vec3 m_norm_unrot;
Mask::ZRotation zRot;
std::uniform_real_distribution<double> m_uniform_fraction;
bool rndmFlag;
};
#endif

50
include/ThreeStepSystem.h
View File

@ -1,50 +0,0 @@
#ifndef THREESTEPSYSTEM_H
#define THREESTEPSYSTEM_H
#include "ReactionSystem.h"
#include "AngularDistribution.h"
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;
const std::vector<Nucleus>& GetNuclei() override;
inline void SetDecay1Distribution(const std::string& filename) { decay1dist.ReadDistributionFile(filename); };
inline void SetDecay2Distribution(const std::string& filename) { decay2dist.ReadDistributionFile(filename); };
inline void SetReactionThetaType(int type) { step1.SetEjectileThetaType(type); };
inline const Nucleus& GetTarget() const { return step1.GetTarget(); };
inline const Nucleus& GetProjectile() const { return step1.GetProjectile(); };
inline const Nucleus& GetEjectile() const { return step1.GetEjectile(); };
inline const Nucleus& GetResidual() const { return step1.GetResidual(); };
inline const Nucleus& GetBreakup1() const { return step2.GetEjectile(); };
inline const Nucleus& GetBreakup2() const { return step2.GetResidual(); };
inline const Nucleus& GetBreakup3() const { return step3.GetEjectile(); };
inline const Nucleus& GetBreakup4() const { return step3.GetResidual(); };
inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); };
inline int GetDecay2AngularMomentum(){ return decay2dist.GetL(); };
inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); };
inline double GetDecay2BranchingRatio(){ return decay2dist.GetBranchingRatio(); };
protected:
void LinkTarget() override;
void SetSystemEquation() override;
std::uniform_real_distribution<double> m_phi2Range;
Reaction step1, step2, step3;
AngularDistribution decay1dist, decay2dist;
};
}
#endif

45
include/TwoStepSystem.h
View File

@ -1,45 +0,0 @@
#ifndef TWOSTEPSYSTEM_H
#define TWOSTEPSYSTEM_H
#include "ReactionSystem.h"
#include "AngularDistribution.h"
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;
const std::vector<Nucleus>& GetNuclei() override;
inline void SetDecay1Distribution(const std::string& filename) { decay1dist.ReadDistributionFile(filename); };
inline void SetReactionThetaType(int type) { step1.SetEjectileThetaType(type); };
inline const Nucleus& GetTarget() const { return step1.GetTarget(); };
inline const Nucleus& GetProjectile() const { return step1.GetProjectile(); };
inline const Nucleus& GetEjectile() const { return step1.GetEjectile(); };
inline const Nucleus& GetResidual() const { return step1.GetResidual(); };
inline const Nucleus& GetBreakup1() const { return step2.GetEjectile(); };
inline const Nucleus& GetBreakup2() const { return step2.GetResidual(); };
inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); };
inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); };
private:
void LinkTarget() override;
void SetSystemEquation() override;
std::uniform_real_distribution<double> m_phi2Range;
Reaction step1, step2;
AngularDistribution decay1dist;
};
}
#endif

18
input.txt
View File

@ -1,28 +1,26 @@
----------Data Information----------
OutputFile: /data1/gwm17/7BeNov2021/sim/7Bedp_1100keV_beam_50CD2.mask
OutputFile: /media/data/gwm17/mask_tests/10B3Hea_test.root
----------Reaction Information----------
ReactionType: TwoStepRxn
ReactionType: OneStepRxn
Z A (order is target, projectile, ejectile, break1, break3(if pure decay is target, ejectile))
1 2
4 7
1 1
5 10
2 3
2 4
----------Target Information----------
NumberOfLayers: 1
begin_layer
Thickness(ug/cm^2): 50
begin_elements (Z, A, Stoich.)
element 1 2 2
element 6 12 1
element 5 10 1
end_elements
end_layer
----------Sampling Information----------
NumberOfSamples: 100000
BeamMeanEnergy(MeV): 1.1 BeamEnergySigma(MeV): 0.0
NumberOfSamples: 1000
BeamMeanEnergy(MeV): 24.0 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
ResidualExMean(MeV): 0.0 ResidualExSigma(MeV): 0.0
ResidualExMean(MeV): 16.8 ResidualExSigma(MeV): 0.023
Decay1_DistributionFile: ./etc/isotropic_dist.txt
Decay2_DistributionFile: ./etc/isotropic_dist.txt
--------------------------------------

2
include/AnasenDeadChannelMap.h → src/Detectors/AnasenDeadChannelMap.h
View File

@ -25,7 +25,7 @@ public:
AnasenDeadChannelMap(const std::string& filename);
~AnasenDeadChannelMap();
void LoadMapfile(const std::string& filename);
inline const bool IsValid() const { return valid_flag; }
const bool IsValid() const { return valid_flag; }
const bool IsDead(AnasenDetectorType type, int detIndex, int channel, AnasenDetectorSide side) const;
private:
void InitMap();

504
src/Detectors/AnasenEfficiency.cpp
View File

@ -3,194 +3,221 @@
#include <iomanip>
#include <iostream>
#include "TFile.h"
#include "TTree.h"
AnasenEfficiency::AnasenEfficiency() :
DetectorEfficiency(), det_silicon(si_thickness)
DetectorEfficiency(), m_detectorEloss({14}, {28}, {1}, s_detectorThickness)
{
for(int i=0; i<n_sx3_per_ring; i++) {
m_Ring1.emplace_back(4, sx3_length, sx3_width, ring_phi[i], ring1_z, ring_rho[i]);
m_Ring1[i].TurnOnRandomizedCoordinates();
m_Ring2.emplace_back(4, sx3_length, sx3_width, ring_phi[i], ring2_z, ring_rho[i]);
m_Ring2[i].TurnOnRandomizedCoordinates();
for(int i=0; i<s_nSX3PerBarrel; i++)
{
m_Ring1.emplace_back(4, s_sx3Length, s_sx3Width, s_barrelPhiList[i], s_barrel1Z, s_barrelRhoList[i]);
m_Ring1[i].SetPixelSmearing(true);
m_Ring2.emplace_back(4, s_sx3Length, s_sx3Width, s_barrelPhiList[i], s_barrel2Z, s_barrelRhoList[i]);
m_Ring2[i].SetPixelSmearing(true);
}
for(int i=0; i<n_qqq; i++) {
m_forwardQQQs.emplace_back(qqq_rinner, qqq_router, qqq_deltaphi, qqq_phi[i], qqq_z[i]);
m_forwardQQQs[i].TurnOnRandomizedCoordinates();
m_backwardQQQs.emplace_back(qqq_rinner, qqq_router, qqq_deltaphi, qqq_phi[i], (-1.0)*qqq_z[i]);
m_backwardQQQs[i].TurnOnRandomizedCoordinates();
for(int i=0; i<s_nQQQ; i++)
{
m_forwardQQQs.emplace_back(s_qqqInnerR, s_qqqOuterR, s_qqqDeltaPhi, s_qqqPhiList[i], s_qqqZList[i]);
m_forwardQQQs[i].SetSmearing(true);
m_backwardQQQs.emplace_back(s_qqqInnerR, s_qqqOuterR, s_qqqDeltaPhi, s_qqqPhiList[i], (-1.0)*s_qqqZList[i]);
m_backwardQQQs[i].SetSmearing(true);
}
std::vector<int> det_z = {14};
std::vector<int> det_a = {28};
std::vector<int> det_stoich = {1};
det_silicon.SetElements(det_z, det_a, det_stoich);
}
AnasenEfficiency::~AnasenEfficiency() {}
void AnasenEfficiency::DrawDetectorSystem(const std::string& filename) {
void AnasenEfficiency::DrawDetectorSystem(const std::string& filename)
{
std::ofstream output(filename);
std::vector<double> x, y, z;
std::vector<double> cx, cy, cz;
Mask::Vec3 coords;
for(int i=0; i<n_sx3_per_ring; i++) {
for(int j=0; j<4; j++) {
for(int k=0; k<4; k++) {
ROOT::Math::XYZPoint coords;
for(int i=0; i<s_nSX3PerBarrel; i++)
{
for(int j=0; j<4; j++)
{
for(int k=0; k<4; k++)
{
coords = m_Ring1[i].GetRotatedFrontStripCoordinates(j, k);
x.push_back(coords.GetX());
y.push_back(coords.GetY());
z.push_back(coords.GetZ());
x.push_back(coords.X());
y.push_back(coords.Y());
z.push_back(coords.Z());
coords = m_Ring1[i].GetRotatedBackStripCoordinates(j, k);
x.push_back(coords.GetX());
y.push_back(coords.GetY());
z.push_back(coords.GetZ());
x.push_back(coords.X());
y.push_back(coords.Y());
z.push_back(coords.Z());
}
coords = m_Ring1[i].GetHitCoordinates(j, 0);
cx.push_back(coords.GetX());
cy.push_back(coords.GetY());
cz.push_back(coords.GetZ());
cx.push_back(coords.X());
cy.push_back(coords.Y());
cz.push_back(coords.Z());
}
}
for(int i=0; i<n_sx3_per_ring; i++) {
for(int j=0; j<4; j++) {
for(int k=0; k<4; k++) {
for(int i=0; i<s_nSX3PerBarrel; i++)
{
for(int j=0; j<4; j++)
{
for(int k=0; k<4; k++)
{
coords = m_Ring2[i].GetRotatedFrontStripCoordinates(j, k);
x.push_back(coords.GetX());
y.push_back(coords.GetY());
z.push_back(coords.GetZ());
x.push_back(coords.X());
y.push_back(coords.Y());
z.push_back(coords.Z());
coords = m_Ring2[i].GetRotatedBackStripCoordinates(j, k);
x.push_back(coords.GetX());
y.push_back(coords.GetY());
z.push_back(coords.GetZ());
x.push_back(coords.X());
y.push_back(coords.Y());
z.push_back(coords.Z());
}
coords = m_Ring2[i].GetHitCoordinates(j, 0);
cx.push_back(coords.GetX());
cy.push_back(coords.GetY());
cz.push_back(coords.GetZ());
cx.push_back(coords.X());
cy.push_back(coords.Y());
cz.push_back(coords.Z());
}
}
for(int i=0; i<n_qqq; i++) {
for(int j=0; j<16; j++) {
for(int k=0; k<4; k++) {
for(int i=0; i<s_nQQQ; i++)
{
for(int j=0; j<16; j++)
{
for(int k=0; k<4; k++)
{
coords = m_forwardQQQs[i].GetRingCoordinates(j, k);
x.push_back(coords.GetX());
y.push_back(coords.GetY());
z.push_back(coords.GetZ());
x.push_back(coords.X());
y.push_back(coords.Y());
z.push_back(coords.Z());
coords = m_forwardQQQs[i].GetWedgeCoordinates(j, k);
x.push_back(coords.GetX());
y.push_back(coords.GetY());
z.push_back(coords.GetZ());
x.push_back(coords.X());
y.push_back(coords.Y());
z.push_back(coords.Z());
}
for(int k=0; k<16; k++) {
for(int k=0; k<16; k++)
{
coords = m_forwardQQQs[i].GetHitCoordinates(j, k);
cx.push_back(coords.GetX());
cy.push_back(coords.GetY());
cz.push_back(coords.GetZ());
cx.push_back(coords.X());
cy.push_back(coords.Y());
cz.push_back(coords.Z());
}
}
}
for(int i=0; i<n_qqq; i++) {
for(int j=0; j<16; j++) {
for(int k=0; k<4; k++) {
for(int i=0; i<s_nQQQ; i++)
{
for(int j=0; j<16; j++)
{
for(int k=0; k<4; k++)
{
coords = m_backwardQQQs[i].GetRingCoordinates(j, k);
x.push_back(coords.GetX());
y.push_back(coords.GetY());
z.push_back(coords.GetZ());
x.push_back(coords.X());
y.push_back(coords.Y());
z.push_back(coords.Z());
coords = m_backwardQQQs[i].GetWedgeCoordinates(j, k);
x.push_back(coords.GetX());
y.push_back(coords.GetY());
z.push_back(coords.GetZ());
x.push_back(coords.X());
y.push_back(coords.Y());
z.push_back(coords.Z());
}
for(int k=0; k<16; k++) {
for(int k=0; k<16; k++)
{
coords = m_backwardQQQs[i].GetHitCoordinates(j, k);
cx.push_back(coords.GetX());
cy.push_back(coords.GetY());
cz.push_back(coords.GetZ());
cx.push_back(coords.X());
cy.push_back(coords.Y());
cz.push_back(coords.Z());
}
}
}
output<<"ANASEN Geometry File -- Coordinates for Detectors"<<std::endl;
for(unsigned int i=0; i<x.size(); i++) {
for(std::size_t i=0; i<x.size(); i++)
output<<x[i]<<" "<<y[i]<<" "<<z[i]<<std::endl;
}
for(unsigned int i=0; i<cx.size(); i++) {
for(std::size_t i=0; i<cx.size(); i++)
output<<cx[i]<<" "<<cy[i]<<" "<<cz[i]<<std::endl;
}
output.close();
}
double AnasenEfficiency::RunConsistencyCheck() {
std::vector<Mask::Vec3> r1_points;
std::vector<Mask::Vec3> r2_points;
std::vector<Mask::Vec3> fqqq_points;
std::vector<Mask::Vec3> bqqq_points;
for(int i=0; i<n_sx3_per_ring; i++) {
for(int j=0; j<4; j++) {
double AnasenEfficiency::RunConsistencyCheck()
{
std::vector<ROOT::Math::XYZPoint> r1_points;
std::vector<ROOT::Math::XYZPoint> r2_points;
std::vector<ROOT::Math::XYZPoint> fqqq_points;
std::vector<ROOT::Math::XYZPoint> bqqq_points;
for(int i=0; i<s_nSX3PerBarrel; i++)
{
for(int j=0; j<4; j++)
r1_points.push_back(m_Ring1[i].GetHitCoordinates(j, 0));
}
}
for(int i=0; i<n_sx3_per_ring; i++) {
for(int j=0; j<4; j++) {
for(int i=0; i<s_nSX3PerBarrel; i++)
{
for(int j=0; j<4; j++)
r2_points.push_back(m_Ring2[i].GetHitCoordinates(j, 0));
}
}
for(int i=0; i<n_qqq; i++) {
for(int j=0; j<16; j++) {
for(int k=0; k<16; k++) {
for(int i=0; i<s_nQQQ; i++)
{
for(int j=0; j<16; j++)
{
for(int k=0; k<16; k++)
fqqq_points.push_back(m_forwardQQQs[i].GetHitCoordinates(j, k));
}
}
}
for(int i=0; i<n_qqq; i++) {
for(int j=0; j<16; j++) {
for(int k=0; k<16; k++) {
for(int i=0; i<s_nQQQ; i++)
{
for(int j=0; j<16; j++)
{
for(int k=0; k<16; k++)
bqqq_points.push_back(m_backwardQQQs[i].GetHitCoordinates(j, k));
}
}
}
int npoints = r1_points.size() + r2_points.size() + fqqq_points.size() + bqqq_points.size();
int count = 0;
Mask::Vec3 coords;
for(auto& point : r1_points) {
for(auto& sx3 : m_Ring1) {
auto result = sx3.GetChannelRatio(point.GetTheta(), point.GetPhi());
std::size_t npoints = r1_points.size() + r2_points.size() + fqqq_points.size() + bqqq_points.size();
std::size_t count = 0;
ROOT::Math::XYZPoint coords;
for(auto& point : r1_points)
{
for(auto& sx3 : m_Ring1)
{
auto result = sx3.GetChannelRatio(point.Theta(), point.Phi());
coords = sx3.GetHitCoordinates(result.front_strip_index, result.front_ratio);
if(IsDoubleEqual(point.GetX(), coords.GetX()) && IsDoubleEqual(point.GetY(), coords.GetY()) && IsDoubleEqual(point.GetZ(), coords.GetZ())) {
if(IsDoubleEqual(point.X(), coords.X()) && IsDoubleEqual(point.Y(), coords.Y()) && IsDoubleEqual(point.Z(), coords.Z()))
{
count++;
break;
}
}
}
for(auto& point : r2_points) {
for(auto& sx3 : m_Ring2) {
auto result = sx3.GetChannelRatio(point.GetTheta(), point.GetPhi());
for(auto& point : r2_points)
{
for(auto& sx3 : m_Ring2)
{
auto result = sx3.GetChannelRatio(point.Theta(), point.Phi());
coords = sx3.GetHitCoordinates(result.front_strip_index, result.front_ratio);
if(IsDoubleEqual(point.GetX(), coords.GetX()) && IsDoubleEqual(point.GetY(), coords.GetY()) && IsDoubleEqual(point.GetZ(), coords.GetZ())) {
if(IsDoubleEqual(point.X(), coords.X()) && IsDoubleEqual(point.Y(), coords.Y()) && IsDoubleEqual(point.Z(), coords.Z()))
{
count++;
break;
}
}
}
for(auto& point : fqqq_points) {
for(auto& qqq : m_forwardQQQs) {
auto result = qqq.GetTrajectoryRingWedge(point.GetTheta(), point.GetPhi());
for(auto& point : fqqq_points)
{
for(auto& qqq : m_forwardQQQs)
{
auto result = qqq.GetTrajectoryRingWedge(point.Theta(), point.Phi());
coords = qqq.GetHitCoordinates(result.first, result.second);
if(IsDoubleEqual(point.GetX(), coords.GetX()) && IsDoubleEqual(point.GetY(), coords.GetY()) && IsDoubleEqual(point.GetZ(), coords.GetZ())) {
if(IsDoubleEqual(point.X(), coords.X()) && IsDoubleEqual(point.Y(), coords.Y()) && IsDoubleEqual(point.Z(), coords.Z()))
{
count++;
break;
}
}
}
for(auto& point : bqqq_points) {
for(auto& qqq : m_backwardQQQs) {
auto result = qqq.GetTrajectoryRingWedge(point.GetTheta(), point.GetPhi());
for(auto& point : bqqq_points)
{
for(auto& qqq : m_backwardQQQs)
{
auto result = qqq.GetTrajectoryRingWedge(point.Theta(), point.Phi());
coords = qqq.GetHitCoordinates(result.first, result.second);
if(IsDoubleEqual(point.GetX(), coords.GetX()) && IsDoubleEqual(point.GetY(), coords.GetY()) && IsDoubleEqual(point.GetZ(), coords.GetZ())) {
if(IsDoubleEqual(point.X(), coords.X()) && IsDoubleEqual(point.Y(), coords.Y()) && IsDoubleEqual(point.Z(), coords.Z()))
{
count++;
break;
}
@ -203,22 +230,23 @@ double AnasenEfficiency::RunConsistencyCheck() {
}
DetectorResult AnasenEfficiency::IsRing1(Mask::Nucleus& nucleus) {
DetectorResult AnasenEfficiency::IsRing1(Mask::Nucleus& nucleus)
{
DetectorResult observation;
double thetaIncident;
for(int i=0; i<n_sx3_per_ring; i++) {
auto result = m_Ring1[i].GetChannelRatio(nucleus.GetTheta(), nucleus.GetPhi());
for(int i=0; i<s_nSX3PerBarrel; i++)
{
auto result = m_Ring1[i].GetChannelRatio(nucleus.vec4.Theta(), nucleus.vec4.Phi());
if(result.front_strip_index != -1 /*&& !dmap.IsDead(AnasenDetectorType::Barrel1, i, result.front_strip_index, AnasenDetectorSide::Front)*/
&& !dmap.IsDead(AnasenDetectorType::Barrel1, i, result.back_strip_index, AnasenDetectorSide::Back))
{
observation.detectFlag = true;
observation.direction = m_Ring1[i].GetHitCoordinates(result.front_strip_index, result.front_ratio);
thetaIncident = std::acos(observation.direction.Dot(m_Ring1[i].GetNormRotated())/observation.direction.GetR());
thetaIncident = std::acos(observation.direction.Dot(m_Ring1[i].GetNormRotated())/observation.direction.R());
if(thetaIncident > M_PI/2.0)
thetaIncident = M_PI - thetaIncident;
observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
observation.energy_deposited = m_detectorEloss.GetEnergyLossTotal(nucleus.Z, nucleus.A, nucleus.GetKE(), thetaIncident);
observation.det_name = "R1";
return observation;
}
@ -227,22 +255,23 @@ DetectorResult AnasenEfficiency::IsRing1(Mask::Nucleus& nucleus) {
return observation;
}
DetectorResult AnasenEfficiency::IsRing2(Mask::Nucleus& nucleus) {
DetectorResult AnasenEfficiency::IsRing2(Mask::Nucleus& nucleus)
{
DetectorResult observation;
double thetaIncident;
for(int i=0; i<n_sx3_per_ring; i++) {
auto result = m_Ring2[i].GetChannelRatio(nucleus.GetTheta(), nucleus.GetPhi());
for(int i=0; i<s_nSX3PerBarrel; i++)
{
auto result = m_Ring2[i].GetChannelRatio(nucleus.vec4.Theta(), nucleus.vec4.Phi());
if(result.front_strip_index != -1 /*&& !dmap.IsDead(AnasenDetectorType::Barrel2, i, result.front_strip_index, AnasenDetectorSide::Front)*/
&& !dmap.IsDead(AnasenDetectorType::Barrel2, i, result.back_strip_index, AnasenDetectorSide::Back))
{
observation.detectFlag = true;
observation.direction = m_Ring2[i].GetHitCoordinates(result.front_strip_index, result.front_ratio);
thetaIncident = std::acos(observation.direction.Dot(m_Ring2[i].GetNormRotated())/observation.direction.GetR());
thetaIncident = std::acos(observation.direction.Dot(m_Ring2[i].GetNormRotated())/observation.direction.R());
if(thetaIncident > M_PI/2.0)
thetaIncident = M_PI - thetaIncident;
observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
observation.energy_deposited = m_detectorEloss.GetEnergyLossTotal(nucleus.Z, nucleus.A, nucleus.GetKE(), thetaIncident);
observation.det_name = "R2";
return observation;
}
@ -251,40 +280,42 @@ DetectorResult AnasenEfficiency::IsRing2(Mask::Nucleus& nucleus) {
return observation;
}
DetectorResult AnasenEfficiency::IsQQQ(Mask::Nucleus& nucleus) {
DetectorResult AnasenEfficiency::IsQQQ(Mask::Nucleus& nucleus)
{
DetectorResult observation;
double thetaIncident;
for(int i=0; i<n_qqq; i++) {
auto result = m_forwardQQQs[i].GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi());
for(int i=0; i<s_nQQQ; i++)
{
auto result = m_forwardQQQs[i].GetTrajectoryRingWedge(nucleus.vec4.Theta(), nucleus.vec4.Phi());
if(result.first != -1 /*&& !dmap.IsDead(AnasenDetectorType::FQQQ, i, result.first, AnasenDetectorSide::Front)*/ &&
!dmap.IsDead(AnasenDetectorType::FQQQ, i, result.second, AnasenDetectorSide::Back))
{
observation.detectFlag = true;
observation.direction = m_forwardQQQs[i].GetHitCoordinates(result.first, result.second);
thetaIncident = std::acos(observation.direction.Dot(m_forwardQQQs[i].GetNorm())/observation.direction.GetR());
thetaIncident = std::acos(observation.direction.Dot(m_forwardQQQs[i].GetNorm())/observation.direction.R());
if(thetaIncident > M_PI/2.0)
thetaIncident = M_PI - thetaIncident;
observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
observation.energy_deposited = m_detectorEloss.GetEnergyLossTotal(nucleus.Z, nucleus.A, nucleus.GetKE(), thetaIncident);
observation.det_name = "FQQQ";
return observation;
}
}
for(int i=0; i<n_qqq; i++) {
auto result = m_backwardQQQs[i].GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi());
for(int i=0; i<s_nQQQ; i++)
{
auto result = m_backwardQQQs[i].GetTrajectoryRingWedge(nucleus.vec4.Theta(), nucleus.vec4.Phi());
if(result.first != -1 /*&& !dmap.IsDead(AnasenDetectorType::BQQQ, i, result.first, AnasenDetectorSide::Front)*/ &&
!dmap.IsDead(AnasenDetectorType::BQQQ, i, result.second, AnasenDetectorSide::Back))
{
observation.detectFlag = true;
observation.direction = m_backwardQQQs[i].GetHitCoordinates(result.first, result.second);
thetaIncident = std::acos(observation.direction.Dot(m_backwardQQQs[i].GetNorm())/observation.direction.GetR());
thetaIncident = std::acos(observation.direction.Dot(m_backwardQQQs[i].GetNorm())/observation.direction.R());
if(thetaIncident > M_PI/2.0)
thetaIncident = M_PI - thetaIncident;
observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
observation.energy_deposited = m_detectorEloss.GetEnergyLossTotal(nucleus.Z, nucleus.A, nucleus.GetKE(), thetaIncident);
observation.det_name = "BQQQ";
return observation;
}
@ -293,95 +324,93 @@ DetectorResult AnasenEfficiency::IsQQQ(Mask::Nucleus& nucleus) {
return observation;
}
DetectorResult AnasenEfficiency::IsAnasen(Mask::Nucleus& nucleus) {
DetectorResult AnasenEfficiency::IsAnasen(Mask::Nucleus& nucleus)
{
DetectorResult result;
if(nucleus.GetKE() <= threshold)
if(nucleus.GetKE() <= s_energyThreshold)
return result;
if(!result.detectFlag) {
if(!result.detectFlag)
result = IsRing1(nucleus);
}
if(!result.detectFlag) {
if(!result.detectFlag)
result = IsRing2(nucleus);
}
if(!result.detectFlag) {
if(!result.detectFlag)
result = IsQQQ(nucleus);
}
return result;
}
void AnasenEfficiency::CountCoincidences(const Mask::MaskFileData& data, std::vector<int>& counts, Mask::RxnType rxn_type) {
if (rxn_type == Mask::RxnType::PureDecay && data.detect_flag[1] && data.detect_flag[2])
void AnasenEfficiency::CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts)
{
if (data.size() == 3 && data[1].isDetected && data[2].isDetected)
{
counts[0]++;
}
else if (rxn_type == Mask::RxnType::OneStepRxn && data.detect_flag[2] && data.detect_flag[3])
else if (data.size() == 4 && data[2].isDetected && data[3].isDetected)
{
counts[0]++;
}
else if(rxn_type == Mask::RxnType::TwoStepRxn)
else if(data.size() == 6)
{
if(data.detect_flag[2] && data.detect_flag[4])
if(data[2].isDetected && data[4].isDetected)
{
counts[0]++;
}
if(data.detect_flag[2] && data.detect_flag[5])
if(data[2].isDetected && data[5].isDetected)
{
counts[1]++;
}
if(data.detect_flag[4] && data.detect_flag[5])
if(data[4].isDetected && data[5].isDetected)
{
counts[2]++;
}
if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[5])
if(data[2].isDetected && data[4].isDetected && data[5].isDetected)
{
counts[3]++;
}
}
else if(rxn_type == Mask::RxnType::ThreeStepRxn)
else if(data.size() == 8)
{
if(data.detect_flag[2] && data.detect_flag[4])
if(data[2].isDetected && data[4].isDetected)
{
counts[0]++;
}
if(data.detect_flag[2] && data.detect_flag[6])
if(data[2].isDetected && data[6].isDetected)
{
counts[1]++;
}
if(data.detect_flag[2] && data.detect_flag[7])
if(data[2].isDetected && data[7].isDetected)
{
counts[2]++;
}
if(data.detect_flag[4] && data.detect_flag[6])
if(data[4].isDetected && data[6].isDetected)
{
counts[3]++;
}
if(data.detect_flag[4] && data.detect_flag[7])
if(data[4].isDetected && data[7].isDetected)
{
counts[4]++;
}
if(data.detect_flag[6] && data.detect_flag[7])
if(data[6].isDetected && data[7].isDetected)
{
counts[5]++;
}
if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[6])
if(data[2].isDetected && data[4].isDetected && data[6].isDetected)
{
counts[6]++;
}
if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[7])
if(data[2].isDetected && data[4].isDetected && data[7].isDetected)
{
counts[7]++;
}
if(data.detect_flag[2] && data.detect_flag[6] && data.detect_flag[7])
if(data[2].isDetected && data[6].isDetected && data[7].isDetected)
{
counts[8]++;
}
if(data.detect_flag[4] && data.detect_flag[6] && data.detect_flag[7])
if(data[4].isDetected && data[6].isDetected && data[7].isDetected)
{
counts[9]++;
}
if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[6] && data.detect_flag[7])
if(data[2].isDetected && data[4].isDetected && data[6].isDetected && data[7].isDetected)
{
counts[10]++;
}
@ -401,141 +430,138 @@ void AnasenEfficiency::CalculateEfficiency(const std::string& inputname, const s
std::cout<<"Loading in output from kinematics simulation: "<<inputname<<std::endl;
std::cout<<"Running efficiency calculation..."<<std::endl;
Mask::MaskFile input(inputname, Mask::MaskFile::FileType::read);
Mask::MaskFile output(outputname, Mask::MaskFile::FileType::write);
TFile* input = TFile::Open(inputname.c_str(), "READ");
TFile* output = TFile::Open(outputname.c_str(), "RECREATE");
std::ofstream stats(statsname);
stats<<std::setprecision(5);
Mask::MaskFileHeader header = input.ReadHeader();
output.WriteHeader(header.rxn_type, header.nsamples);
TTree* intree = (TTree*) input->Get("SimTree");
std::vector<Mask::Nucleus>* dataHandle = new std::vector<Mask::Nucleus>();
intree->SetBranchAddress("nuclei", &dataHandle);
output->cd();
TTree* outtree = new TTree("SimTree", "SimTree");
outtree->Branch("nuclei", dataHandle);
input->cd();
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;
intree->GetEntry(1);
std::vector<int> counts;
std::vector<int> coinc_counts;
switch(header.rxn_type) {
case Mask::RxnType::PureDecay:
counts.resize(3, 0);
coinc_counts.resize(1, 0);
break;
case Mask::RxnType::OneStepRxn:
counts.resize(4, 0);
coinc_counts.resize(1, 0);
break;
case Mask::RxnType::TwoStepRxn:
counts.resize(6, 0);
coinc_counts.resize(4, 0);
break;
case Mask::RxnType::ThreeStepRxn:
counts.resize(8, 0);
coinc_counts.resize(11, 0);
break;
counts.resize(dataHandle->size());
switch(counts.size())
{
case 3: coinc_counts.resize(1, 0); break;
case 4: coinc_counts.resize(1, 0); break;
case 6: coinc_counts.resize(4, 0); break;
case 8: coinc_counts.resize(11, 0); break;
default:
{
std::cerr<<"Bad reaction type at AnasenEfficiency::CalculateEfficiency (given value: "<<GetStringFromRxnType(header.rxn_type)<<"). Quiting..."<<std::endl;
input.Close();
output.Close();
std::cerr<<"Bad reaction type at AnasenEfficiency::CalculateEfficiency (given value: "<<counts.size()<<"). Quiting..."<<std::endl;
input->Close();
output->Close();
stats.close();
return;
}
}
int percent5 = header.nsamples*0.05;
int count = 0;
int npercent = 0;
uint64_t nentries = intree->GetEntries();
uint64_t percent5 = nentries*0.05;
uint64_t count = 0;
uint64_t npercent = 0;
Mask::Nucleus nucleus;
int index=0;
while(true) {
if(++count == percent5) {//Show progress every 5%
for(uint64_t i=0; i<nentries; i++)
{
intree->GetEntry(i);
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 = IsAnasen(nucleus);
if(result.detectFlag) {
data.detect_flag[i] = true;
data.KE[i] = result.energy_deposited;
data.theta[i] = result.direction.GetTheta();
data.phi[i] = result.direction.GetPhi();
counts[i]++;
} else if(data.detect_flag[i] == true) {
data.detect_flag[i] = false;
for(std::size_t j=0; j<dataHandle->size(); j++)
{
Mask::Nucleus& nucleus = (*dataHandle)[j];
DetectorResult result = IsAnasen(nucleus);
if(result.detectFlag)
{
nucleus.isDetected = true;
nucleus.detectedKE = result.energy_deposited;
nucleus.detectedTheta = result.direction.Theta();
nucleus.detectedPhi = result.direction.Phi();
counts[j]++;
}
}
CountCoincidences(data, coinc_counts, header.rxn_type);
CountCoincidences(*dataHandle, coinc_counts);
output.WriteData(data);
index++;
outtree->Fill();
}
input->Close();
output->cd();
outtree->Write(outtree->GetName(), TObject::kOverwrite);
output->Close();
input.Close();
output.Close();
delete dataHandle;
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;
for(unsigned int i=0; i<counts.size(); i++)
{
stats<<std::setw(10)<<i<<"|"<<std::setw(10)<<((double)counts[i]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
}
stats<<"Coincidence Efficiency"<<std::endl;
stats<<"---------------------"<<std::endl;
if(header.rxn_type == Mask::RxnType::PureDecay)
if(dataHandle->size() == 3)
{
stats<<std::setw(10)<<"1 + 2"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"1 + 2"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
}
else if(header.rxn_type == Mask::RxnType::OneStepRxn)
else if(dataHandle->size() == 4)
{
stats<<std::setw(10)<<"2 + 3"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 3"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
}
else if(header.rxn_type == Mask::RxnType::TwoStepRxn)
else if(dataHandle->size() == 6)
{
stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
}
else if(header.rxn_type == Mask::RxnType::ThreeStepRxn)
else if(dataHandle->size() == 8)
{
stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[4]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[4]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[5]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[5]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[6]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[6]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[7]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[7]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[8]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[8]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[9]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[9]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[10]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[10]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
}
stats.close();

64
src/Detectors/AnasenEfficiency.h Normal file
View File

@ -0,0 +1,64 @@
#ifndef ANASEN_EFFICIENCY_H
#define ANASEN_EFFICIENCY_H
#include <string>
#include "DetectorEfficiency.h"
#include "StripDetector.h"
#include "QQQDetector.h"
#include "Target.h"
#include "Nucleus.h"
#include "AnasenDeadChannelMap.h"
class AnasenEfficiency : public DetectorEfficiency
{
public:
AnasenEfficiency();
~AnasenEfficiency();
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;
inline void SetDeadChannelMap(const std::string& filename) { dmap.LoadMapfile(filename); }
private:
DetectorResult IsRing1(Mask::Nucleus& nucleus);
DetectorResult IsRing2(Mask::Nucleus& nucleus);
DetectorResult IsQQQ(Mask::Nucleus& nucleus);
DetectorResult IsAnasen(Mask::Nucleus& nucleus);
void CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts);
std::vector<StripDetector> m_Ring1, m_Ring2;
std::vector<QQQDetector> m_forwardQQQs;
std::vector<QQQDetector> m_backwardQQQs;
Mask::Target m_detectorEloss;
AnasenDeadChannelMap dmap;
/**** ANASEN geometry constants *****/
static constexpr int s_nSX3PerBarrel = 12;
static constexpr int s_nQQQ = 4;
static constexpr double s_sx3Length = 0.075;
static constexpr double s_sx3Width = 0.04;
static constexpr double s_barrelGap = 0.0254;
static constexpr double s_sx3FrameGap = 0.049; //0.049 is base gap due to frames
static constexpr double s_barrel1Z = s_sx3Length/2.0 + s_sx3FrameGap + s_barrelGap/2.0;
static constexpr double s_barrel2Z = (-1.0)*(s_barrelGap/2.0 + s_sx3Length/2.0);
static constexpr double s_qqqZ = 0.0125 + s_sx3Length + s_sx3FrameGap + s_barrelGap/2.0;
static constexpr double s_qqqInnerR = 0.0501;
static constexpr double s_qqqOuterR = 0.0990;
static constexpr double s_qqqDeltaPhi = 1.52119;
static constexpr double s_qqqZList[4] = {s_qqqZ, s_qqqZ, s_qqqZ, s_qqqZ};
static constexpr double s_qqqPhiList[4] = {5.49779, 0.785398, 2.35619, 3.92699};
static constexpr double s_barrelRhoList[12] = {0.0890601, 0.0889871, 0.0890354, 0.0890247, 0.0890354, 0.0890354, 0.0890247,
0.0890354, 0.0890354, 0.0890247, 0.0890354, 0.0890354};
static constexpr double s_barrelPhiList[12] = {4.97426, 5.49739, 6.02132, 0.261868, 0.785398, 1.30893, 1.83266, 2.35619, 2.87972,
3.40346, 3.92699, 4.45052};
/*************************/
static constexpr double s_energyThreshold = 0.6; //MeV
static constexpr double s_deg2rad = M_PI/180.0;
static constexpr double s_detectorThickness = 1000 * 1e-4 * 2.3926 * 1e6; //thickness in um -> eff thickness in ug/cm^2 for detector
};
#endif

13
src/Detectors/CMakeLists.txt
View File

@ -1,18 +1,27 @@
add_executable(DetectEff)
target_include_directories(DetectEff PUBLIC ${MASK_INCLUDE_DIR})
target_include_directories(DetectEff PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..)
target_sources(DetectEff PUBLIC
AnasenDeadChannelMap.cpp
AnasenDeadChannelMap.h
AnasenEfficiency.cpp
DetectorEfficiency.cpp
AnasenEfficiency.h
main.cpp
DetectorEfficiency.h
QQQDetector.cpp
QQQDetector.h
SabreDeadChannelMap.cpp
SabreDeadChannelMap.h
SabreDetector.cpp
SabreDetector.h
SabreEfficiency.cpp
SabreEfficiency.h
StripDetector.cpp
StripDetector.h
)
target_link_libraries(DetectEff
MaskDict
Mask
catima
)

17
include/DetectorEfficiency.h → src/Detectors/DetectorEfficiency.h
View File

@ -4,7 +4,18 @@
#include <string>
#include <cmath>
class DetectorEfficiency {
#include "Math/Point3D.h"
struct DetectorResult
{
bool detectFlag = false;
ROOT::Math::XYZPoint direction;
double energy_deposited = 0.0;
std::string det_name = "";
};
class DetectorEfficiency
{
public:
DetectorEfficiency() {};
virtual ~DetectorEfficiency() {};
@ -14,9 +25,9 @@ public:
virtual double RunConsistencyCheck() = 0;
protected:
inline bool IsDoubleEqual(double x, double y) { return std::fabs(x-y) < epsilon ? true : false; };
inline bool IsDoubleEqual(double x, double y) { return std::fabs(x-y) < s_epsilon ? true : false; };
static constexpr double epsilon = 1.0e-6;
static constexpr double s_epsilon = 1.0e-6;
};
#endif

163
src/Detectors/QQQDetector.cpp
View File

@ -1,107 +1,116 @@
#include "QQQDetector.h"
QQQDetector::QQQDetector(double R_in, double R_out, double deltaPhi, double phiCentral, double z, double x, double y) :
m_Rinner(R_in), m_Router(R_out), m_deltaPhi(deltaPhi), m_phiCentral(phiCentral), m_translation(x,y,z), m_norm(0.0,0.0,1.0), m_uniform_fraction(0.0, 1.0), rndmFlag(false)
m_innerR(R_in), m_outerR(R_out), m_deltaPhi(deltaPhi), m_centralPhi(phiCentral), m_translation(x,y,z), m_norm(0.0,0.0,1.0),
m_uniformFraction(0.0, 1.0), m_isSmearing(false)
{
m_deltaR = (m_Router - m_Rinner)/nrings;
m_deltaPhi_per_wedge = m_deltaPhi/nwedges;
m_ZRot.SetAngle(m_phiCentral);
m_ringCoords.resize(nrings);
m_wedgeCoords.resize(nwedges);
for(auto& ring : m_ringCoords) {
m_deltaR = (m_outerR - m_innerR)/s_nRings;
m_deltaPhiWedge = m_deltaPhi/s_nWedges;
m_zRotation.SetAngle(m_centralPhi);
m_ringCoords.resize(s_nRings);
m_wedgeCoords.resize(s_nWedges);
for(auto& ring : m_ringCoords)
ring.resize(4);
}
for(auto& wedge : m_wedgeCoords) {
for(auto& wedge : m_wedgeCoords)
wedge.resize(4);
}
CalculateCorners();
}
QQQDetector::~QQQDetector() {}
void QQQDetector::CalculateCorners() {
void QQQDetector::CalculateCorners()
{
double x0, x1, x2, x3;
double y0, y1, y2, y3;
double z0, z1, z2, z3;
//Generate flat ring corner coordinates
for(int i=0; i<nrings; i++) {
x0 = (m_Rinner + m_deltaR*(i+1))*std::cos(-m_deltaPhi/2.0);
y0 = (m_Rinner + m_deltaR*(i+1))*std::sin(-m_deltaPhi/2.0);
for(int i=0; i<s_nRings; i++)
{
x0 = (m_innerR + m_deltaR*(i+1))*std::cos(-m_deltaPhi/2.0);
y0 = (m_innerR + m_deltaR*(i+1))*std::sin(-m_deltaPhi/2.0);
z0 = 0.0;
m_ringCoords[i][0].SetVectorCartesian(x0, y0, z0);
m_ringCoords[i][0].SetXYZ(x0, y0, z0);
x1 = (m_Rinner + m_deltaR*(i))*std::cos(-m_deltaPhi/2.0);
y1 = (m_Rinner + m_deltaR*(i))*std::sin(-m_deltaPhi/2.0);
x1 = (m_innerR + m_deltaR*(i))*std::cos(-m_deltaPhi/2.0);
y1 = (m_innerR + m_deltaR*(i))*std::sin(-m_deltaPhi/2.0);
z1 = 0.0;
m_ringCoords[i][1].SetVectorCartesian(x1, y1, z1);
m_ringCoords[i][1].SetXYZ(x1, y1, z1);
x2 = (m_Rinner + m_deltaR*(i))*std::cos(m_deltaPhi/2.0);
y2 = (m_Rinner + m_deltaR*(i))*std::sin(m_deltaPhi/2.0);
x2 = (m_innerR + m_deltaR*(i))*std::cos(m_deltaPhi/2.0);
y2 = (m_innerR + m_deltaR*(i))*std::sin(m_deltaPhi/2.0);
z2 = 0.0;
m_ringCoords[i][2].SetVectorCartesian(x2, y2, z2);
m_ringCoords[i][2].SetXYZ(x2, y2, z2);
x3 = (m_Rinner + m_deltaR*(i+1))*std::cos(m_deltaPhi/2.0);
y3 = (m_Rinner + m_deltaR*(i+1))*std::sin(m_deltaPhi/2.0);
x3 = (m_innerR + m_deltaR*(i+1))*std::cos(m_deltaPhi/2.0);
y3 = (m_innerR + m_deltaR*(i+1))*std::sin(m_deltaPhi/2.0);
z3 = 0.0;
m_ringCoords[i][3].SetVectorCartesian(x3, y3, z3);
m_ringCoords[i][3].SetXYZ(x3, y3, z3);
}
//Generate flat wedge corner coordinates
for(int i=0; i<nwedges; i++) {
x0 = m_Router * std::cos(-m_deltaPhi/2.0 + i*m_deltaPhi_per_wedge);
y0 = m_Router * std::sin(-m_deltaPhi/2.0 + i*m_deltaPhi_per_wedge);
for(int i=0; i<s_nWedges; i++)
{
x0 = m_outerR * std::cos(-m_deltaPhi/2.0 + i*m_deltaPhiWedge);
y0 = m_outerR * std::sin(-m_deltaPhi/2.0 + i*m_deltaPhiWedge);
z0 = 0.0;
m_wedgeCoords[i][0].SetVectorCartesian(x0, y0, z0);
m_wedgeCoords[i][0].SetXYZ(x0, y0, z0);
x1 = m_Rinner * std::cos(-m_deltaPhi/2.0 + i*m_deltaPhi_per_wedge);
y1 = m_Rinner * std::sin(-m_deltaPhi/2.0 + i*m_deltaPhi_per_wedge);
x1 = m_innerR * std::cos(-m_deltaPhi/2.0 + i*m_deltaPhiWedge);
y1 = m_innerR * std::sin(-m_deltaPhi/2.0 + i*m_deltaPhiWedge);
z1 = 0.0;
m_wedgeCoords[i][1].SetVectorCartesian(x1, y1, z1);
m_wedgeCoords[i][1].SetXYZ(x1, y1, z1);
x2 = m_Rinner * std::cos(-m_deltaPhi/2.0 + (i+1)*m_deltaPhi_per_wedge);
y2 = m_Rinner * std::sin(-m_deltaPhi/2.0 + (i+1)*m_deltaPhi_per_wedge);
x2 = m_innerR * std::cos(-m_deltaPhi/2.0 + (i+1)*m_deltaPhiWedge);
y2 = m_innerR * std::sin(-m_deltaPhi/2.0 + (i+1)*m_deltaPhiWedge);
z2 = 0.0;
m_wedgeCoords[i][2].SetVectorCartesian(x2, y2, z2);
m_wedgeCoords[i][2].SetXYZ(x2, y2, z2);
x3 = m_Router * std::cos(-m_deltaPhi/2.0 + (i+1)*m_deltaPhi_per_wedge);
y3 = m_Router * std::sin(-m_deltaPhi/2.0 + (i+1)*m_deltaPhi_per_wedge);
x3 = m_outerR * std::cos(-m_deltaPhi/2.0 + (i+1)*m_deltaPhiWedge);
y3 = m_outerR * std::sin(-m_deltaPhi/2.0 + (i+1)*m_deltaPhiWedge);
z3 = 0.0;
m_wedgeCoords[i][3].SetVectorCartesian(x3, y3, z3);
m_wedgeCoords[i][3].SetXYZ(x3, y3, z3);
}
for(int i=0; i<nrings; i++) {
for(int j=0; j<4; j++) {
for(int i=0; i<s_nRings; i++)
{
for(int j=0; j<4; j++)
m_ringCoords[i][j] = TransformCoordinates(m_ringCoords[i][j]);
}
}
for(int i=0; i<nwedges; i++) {
for(int j=0; j<4; j++) {
for(int i=0; i<s_nWedges; i++)
{
for(int j=0; j<4; j++)
m_wedgeCoords[i][j] = TransformCoordinates(m_wedgeCoords[i][j]);
}
}
}
Mask::Vec3 QQQDetector::GetTrajectoryCoordinates(double theta, double phi) {
double z_to_detector = m_translation.GetZ();
ROOT::Math::XYZPoint QQQDetector::GetTrajectoryCoordinates(double theta, double phi)
{
double z_to_detector = m_translation.Vect().Z();
double rho_traj = z_to_detector*std::tan(theta);
double r_traj = std::sqrt(rho_traj*rho_traj + z_to_detector*z_to_detector);
double min_rho, max_rho, min_phi, max_phi;
Mask::Vec3 result;
ROOT::Math::XYZPoint result;
for(auto& ring : m_ringCoords) {
min_rho = ring[1].GetRho();
max_rho = ring[0].GetRho();
if(rho_traj < max_rho && rho_traj > min_rho) {
for(auto& wedge : m_wedgeCoords) {
min_phi = wedge[0].GetPhi();
max_phi = wedge[3].GetPhi();
if(phi < min_phi && phi < max_phi) {
result.SetVectorSpherical(r_traj, theta, phi);
for(auto& ring : m_ringCoords)
{
min_rho = ring[1].Rho();
max_rho = ring[0].Rho();
if(rho_traj < max_rho && rho_traj > min_rho)
{
for(auto& wedge : m_wedgeCoords)
{
min_phi = wedge[0].Phi();
max_phi = wedge[3].Phi();
if(phi < min_phi && phi < max_phi)
{
result.SetXYZ(std::sin(theta)*std::cos(phi)*r_traj,
std::sin(theta)*std::sin(phi)*r_traj,
std::cos(theta)*r_traj);
break;
}
}
@ -109,55 +118,57 @@ Mask::Vec3 QQQDetector::GetTrajectoryCoordinates(double theta, double phi) {
}
return result;
}
std::pair<int,int> QQQDetector::GetTrajectoryRingWedge(double theta, double phi) {
double z_to_detector = m_translation.GetZ();
std::pair<int,int> QQQDetector::GetTrajectoryRingWedge(double theta, double phi)
{
double z_to_detector = m_translation.Vect().Z();
double rho_traj = z_to_detector*std::tan(theta);
double min_rho, max_rho, min_phi, max_phi;
for(int r=0; r<nrings; r++) {
for(int r=0; r<s_nRings; r++)
{
auto& ring = m_ringCoords[r];
min_rho = ring[1].GetRho();
max_rho = ring[0].GetRho();
if(rho_traj < max_rho && rho_traj > min_rho) {
for(int w=0; w<nwedges; w++) {
min_rho = ring[1].Rho();
max_rho = ring[0].Rho();
if(rho_traj < max_rho && rho_traj > min_rho)
{
for(int w=0; w<s_nWedges; w++)
{
auto& wedge = m_wedgeCoords[w];
min_phi = wedge[0].GetPhi();
max_phi = wedge[3].GetPhi();
if(phi > min_phi && phi < max_phi) {
min_phi = wedge[0].Phi();
max_phi = wedge[3].Phi();
if(phi > min_phi && phi < max_phi)
return std::make_pair(r, w);
}
}
}
}
return std::make_pair(-1, -1);
}
Mask::Vec3 QQQDetector::GetHitCoordinates(int ringch, int wedgech) {
ROOT::Math::XYZPoint QQQDetector::GetHitCoordinates(int ringch, int wedgech)
{
if(!CheckChannel(ringch) || !CheckChannel(wedgech))
return Mask::Vec3();
return ROOT::Math::XYZPoint();
double r_center, phi_center;
if(rndmFlag)
if(m_isSmearing)
{
r_center = m_Rinner + (m_uniform_fraction(Mask::RandomGenerator::GetInstance().GetGenerator())+ringch)*m_deltaR;
phi_center = -m_deltaPhi/2.0 + (m_uniform_fraction(Mask::RandomGenerator::GetInstance().GetGenerator())+wedgech)*m_deltaPhi_per_wedge;
r_center = m_innerR + (m_uniformFraction(Mask::RandomGenerator::GetInstance().GetGenerator())+ringch)*m_deltaR;
phi_center = -m_deltaPhi/2.0 + (m_uniformFraction(Mask::RandomGenerator::GetInstance().GetGenerator())+wedgech)*m_deltaPhiWedge;
}
else
{
r_center = m_Rinner + (0.5+ringch)*m_deltaR;
phi_center = -m_deltaPhi/2.0 + (0.5+wedgech)*m_deltaPhi_per_wedge;
r_center = m_innerR + (0.5+ringch)*m_deltaR;
phi_center = -m_deltaPhi/2.0 + (0.5+wedgech)*m_deltaPhiWedge;
}
double x = r_center*std::cos(phi_center);
double y = r_center*std::sin(phi_center);
double z = 0;
Mask::Vec3 hit(x, y, z);
ROOT::Math::XYZPoint hit(x, y, z);
return TransformCoordinates(hit);
}

65
src/Detectors/QQQDetector.h Normal file
View File

@ -0,0 +1,65 @@
/*
QQQDetector.h
Class implementing geometry for QQQDetector where the detector is perpendicular to the beam axis.
Detector is first generated centered on the x-axis (phi=0)
Coordinate convention : +z is downstream, -z is upstream. +y is vertically down in the lab.
*/
#ifndef QQQDETECTOR_H
#define QQQDETECTOR_H
#include <cmath>
#include <vector>
#include "RandomGenerator.h"
#include "Math/Point3D.h"
#include "Math/Vector3D.h"
#include "Math/RotationZ.h"
#include "Math/Translation3D.h"
class QQQDetector
{
public:
QQQDetector(double R_in, double R_out, double deltaPhi, double phiCentral, double z, double x=0, double y=0);
~QQQDetector();
const ROOT::Math::XYZPoint& GetRingCoordinates(int ringch, int corner) { return m_ringCoords[ringch][corner]; }
const ROOT::Math::XYZPoint& GetWedgeCoordinates(int wedgech, int corner) { return m_wedgeCoords[wedgech][corner]; }
const ROOT::Math::XYZVector& GetNorm() { return m_norm; }
ROOT::Math::XYZPoint GetTrajectoryCoordinates(double theta, double phi);
std::pair<int, int> GetTrajectoryRingWedge(double theta, double phi);
ROOT::Math::XYZPoint GetHitCoordinates(int ringch, int wedgech);
void SetSmearing(bool isSmearing) { m_isSmearing = isSmearing; }
int GetNumberOfRings() { return s_nRings; }
int GetNumberOfWedges() { return s_nWedges; }
private:
bool CheckChannel(int ch) { return (ch >=0 && ch < s_nRings); }
bool CheckCorner(int corner) { return (corner >=0 && corner < 4); }
void CalculateCorners();
ROOT::Math::XYZPoint TransformCoordinates(ROOT::Math::XYZPoint& vector) { return m_translation * (m_zRotation * vector) ; }
double m_innerR;
double m_outerR;
double m_deltaR;
double m_deltaPhi;
double m_deltaPhiWedge;
double m_centralPhi;
std::vector<std::vector<ROOT::Math::XYZPoint>> m_ringCoords, m_wedgeCoords;
ROOT::Math::Translation3D m_translation;
ROOT::Math::XYZVector m_norm;
ROOT::Math::RotationZ m_zRotation;
std::uniform_real_distribution<double> m_uniformFraction;
bool m_isSmearing;
static constexpr int s_nRings = 16;
static constexpr int s_nWedges = 16;
static constexpr double s_deg2rad = M_PI/180.0;
};
#endif

0
include/SabreDeadChannelMap.h → src/Detectors/SabreDeadChannelMap.h
View File

242
src/Detectors/SabreDetector.cpp
View File

@ -53,126 +53,133 @@
#include "SabreDetector.h"
SabreDetector::SabreDetector() :
m_Router(0.1351), m_Rinner(0.0326), m_deltaPhi_flat(54.4*deg2rad), m_phiCentral(0.0), m_tilt(0.0), m_translation(0.,0.,0.), m_norm_flat(0,0,1.0),
m_detectorID(-1)
m_outerR(0.1351), m_innerR(0.0326), m_deltaPhiFlat(54.4*s_deg2rad), m_centerPhi(0.0), m_tilt(0.0),
m_translation(0.,0.,0.), m_norm(0,0,1.0), m_detectorID(-1)
{
m_YRot.SetAngle(m_tilt);
m_ZRot.SetAngle(m_phiCentral);
m_yRotation.SetAngle(m_tilt);
m_zRotation.SetAngle(m_centerPhi);
//Initialize the coordinate arrays
m_ringCoords_flat.resize(m_nRings);
m_ringCoords_tilt.resize(m_nRings);
m_wedgeCoords_flat.resize(m_nWedges);
m_wedgeCoords_tilt.resize(m_nWedges);
for(int i=0; i<m_nRings; i++) {
m_ringCoords_flat[i].resize(4);
m_ringCoords_tilt[i].resize(4);
m_flatRingCoords.resize(s_nRings);
m_tiltRingCoords.resize(s_nRings);
m_flatWedgeCoords.resize(s_nWedges);
m_tiltWedgeCoords.resize(s_nWedges);
for(int i=0; i<s_nRings; i++) {
m_flatRingCoords[i].resize(4);
m_tiltRingCoords[i].resize(4);
}
for(int i=0; i<m_nWedges; i++) {
m_wedgeCoords_flat[i].resize(4);
m_wedgeCoords_tilt[i].resize(4);
for(int i=0; i<s_nWedges; i++) {
m_flatWedgeCoords[i].resize(4);
m_tiltWedgeCoords[i].resize(4);
}
m_deltaR_flat = m_Router - m_Rinner;
m_deltaR_flat_ring = m_deltaR_flat/m_nRings;
m_deltaRFlat = m_outerR - m_innerR;
m_deltaRFlatRing = m_deltaRFlat/s_nRings;
m_deltaPhiFlatWedge = m_deltaPhiFlat/s_nWedges;
CalculateCorners();
}
SabreDetector::SabreDetector(int detID, double Rin, double Rout, double deltaPhi_flat, double phiCentral, double tiltFromVert, double zdist, double xdist, double ydist) :
m_Router(Rout), m_Rinner(Rin), m_deltaPhi_flat(deltaPhi_flat), m_phiCentral(phiCentral), m_tilt(tiltFromVert), m_translation(xdist, ydist, zdist), m_norm_flat(0,0,1.0),
m_detectorID(detID)
SabreDetector::SabreDetector(int detID, double Rin, double Rout, double deltaPhi_flat, double phiCentral,
double tiltFromVert, double zdist, double xdist, double ydist) :
m_outerR(Rout), m_innerR(Rin), m_deltaPhiFlat(deltaPhi_flat), m_centerPhi(phiCentral), m_tilt(tiltFromVert),
m_translation(xdist, ydist, zdist), m_norm(0,0,1.0), m_detectorID(detID)
{
m_YRot.SetAngle(m_tilt);
m_ZRot.SetAngle(m_phiCentral);
m_yRotation.SetAngle(m_tilt);
m_zRotation.SetAngle(m_centerPhi);
//Initialize coordinate arrays
m_ringCoords_flat.resize(m_nRings);
m_ringCoords_tilt.resize(m_nRings);
m_wedgeCoords_flat.resize(m_nWedges);
m_wedgeCoords_tilt.resize(m_nWedges);
for(int i=0; i<m_nRings; i++) {
m_ringCoords_flat[i].resize(4);
m_ringCoords_tilt[i].resize(4);
m_flatRingCoords.resize(s_nRings);
m_tiltRingCoords.resize(s_nRings);
m_flatWedgeCoords.resize(s_nWedges);
m_tiltWedgeCoords.resize(s_nWedges);
for(int i=0; i<s_nRings; i++)
{
m_flatRingCoords[i].resize(4);
m_tiltRingCoords[i].resize(4);
}
for(int i=0; i<m_nWedges; i++) {
m_wedgeCoords_flat[i].resize(4);
m_wedgeCoords_tilt[i].resize(4);
for(int i=0; i<s_nWedges; i++)
{
m_flatWedgeCoords[i].resize(4);
m_tiltWedgeCoords[i].resize(4);
}
m_deltaR_flat = m_Router - m_Rinner;
m_deltaR_flat_ring = m_deltaR_flat/m_nRings;
m_deltaPhi_flat_wedge = m_deltaPhi_flat/m_nWedges;
m_deltaRFlat = m_outerR - m_innerR;
m_deltaRFlatRing = m_deltaRFlat/s_nRings;
m_deltaPhiFlatWedge = m_deltaPhiFlat/s_nWedges;
CalculateCorners();
}
SabreDetector::~SabreDetector() {}
void SabreDetector::CalculateCorners() {
void SabreDetector::CalculateCorners()
{
double x0, x1, x2, x3;
double y0, y1, y2, y3;
double z0, z1, z2, z3;
//Generate flat ring corner coordinates
for(int i=0; i<m_nRings; i++) {
x0 = (m_Rinner + m_deltaR_flat_ring*(i+1))*std::cos(-m_deltaPhi_flat/2.0);
y0 = (m_Rinner + m_deltaR_flat_ring*(i+1))*std::sin(-m_deltaPhi_flat/2.0);
for(int i=0; i<s_nRings; i++)
{
x0 = (m_innerR + m_deltaRFlatRing*(i+1))*std::cos(-m_deltaPhiFlat/2.0);
y0 = (m_innerR + m_deltaRFlatRing*(i+1))*std::sin(-m_deltaPhiFlat/2.0);
z0 = 0.0;
m_ringCoords_flat[i][0].SetVectorCartesian(x0, y0, z0);
m_flatRingCoords[i][0].SetXYZ(x0, y0, z0);
x1 = (m_Rinner + m_deltaR_flat_ring*(i))*std::cos(-m_deltaPhi_flat/2.0);
y1 = (m_Rinner + m_deltaR_flat_ring*(i))*std::sin(-m_deltaPhi_flat/2.0);
x1 = (m_innerR + m_deltaRFlatRing*(i))*std::cos(-m_deltaPhiFlat/2.0);
y1 = (m_innerR + m_deltaRFlatRing*(i))*std::sin(-m_deltaPhiFlat/2.0);
z1 = 0.0;
m_ringCoords_flat[i][1].SetVectorCartesian(x1, y1, z1);
m_flatRingCoords[i][1].SetXYZ(x1, y1, z1);
x2 = (m_Rinner + m_deltaR_flat_ring*(i))*std::cos(m_deltaPhi_flat/2.0);
y2 = (m_Rinner + m_deltaR_flat_ring*(i))*std::sin(m_deltaPhi_flat/2.0);
x2 = (m_innerR + m_deltaRFlatRing*(i))*std::cos(m_deltaPhiFlat/2.0);
y2 = (m_innerR + m_deltaRFlatRing*(i))*std::sin(m_deltaPhiFlat/2.0);
z2 = 0.0;
m_ringCoords_flat[i][2].SetVectorCartesian(x2, y2, z2);
m_flatRingCoords[i][2].SetXYZ(x2, y2, z2);
x3 = (m_Rinner + m_deltaR_flat_ring*(i+1))*std::cos(m_deltaPhi_flat/2.0);
y3 = (m_Rinner + m_deltaR_flat_ring*(i+1))*std::sin(m_deltaPhi_flat/2.0);
x3 = (m_innerR + m_deltaRFlatRing*(i+1))*std::cos(m_deltaPhiFlat/2.0);
y3 = (m_innerR + m_deltaRFlatRing*(i+1))*std::sin(m_deltaPhiFlat/2.0);
z3 = 0.0;
m_ringCoords_flat[i][3].SetVectorCartesian(x3, y3, z3);
m_flatRingCoords[i][3].SetXYZ(x3, y3, z3);
}
//Generate flat wedge corner coordinates
for(int i=0; i<m_nWedges; i++) {
x0 = m_Router * std::cos(-m_deltaPhi_flat/2.0 + i*m_deltaPhi_flat_wedge);
y0 = m_Router * std::sin(-m_deltaPhi_flat/2.0 + i*m_deltaPhi_flat_wedge);
for(int i=0; i<s_nWedges; i++)
{
x0 = m_outerR * std::cos(-m_deltaPhiFlat/2.0 + i*m_deltaPhiFlatWedge);
y0 = m_outerR * std::sin(-m_deltaPhiFlat/2.0 + i*m_deltaPhiFlatWedge);
z0 = 0.0;
m_wedgeCoords_flat[i][0].SetVectorCartesian(x0, y0, z0);
m_flatWedgeCoords[i][0].SetXYZ(x0, y0, z0);
x1 = m_Rinner * std::cos(-m_deltaPhi_flat/2.0 + i*m_deltaPhi_flat_wedge);
y1 = m_Rinner * std::sin(-m_deltaPhi_flat/2.0 + i*m_deltaPhi_flat_wedge);
x1 = m_innerR * std::cos(-m_deltaPhiFlat/2.0 + i*m_deltaPhiFlatWedge);
y1 = m_innerR * std::sin(-m_deltaPhiFlat/2.0 + i*m_deltaPhiFlatWedge);
z1 = 0.0;
m_wedgeCoords_flat[i][1].SetVectorCartesian(x1, y1, z1);
m_flatWedgeCoords[i][1].SetXYZ(x1, y1, z1);
x2 = m_Rinner * std::cos(-m_deltaPhi_flat/2.0 + (i+1)*m_deltaPhi_flat_wedge);
y2 = m_Rinner * std::sin(-m_deltaPhi_flat/2.0 + (i+1)*m_deltaPhi_flat_wedge);
x2 = m_innerR * std::cos(-m_deltaPhiFlat/2.0 + (i+1)*m_deltaPhiFlatWedge);
y2 = m_innerR * std::sin(-m_deltaPhiFlat/2.0 + (i+1)*m_deltaPhiFlatWedge);
z2 = 0.0;
m_wedgeCoords_flat[i][2].SetVectorCartesian(x2, y2, z2);
m_flatWedgeCoords[i][2].SetXYZ(x2, y2, z2);
x3 = m_Router * std::cos(-m_deltaPhi_flat/2.0 + (i+1)*m_deltaPhi_flat_wedge);
y3 = m_Router * std::sin(-m_deltaPhi_flat/2.0 + (i+1)*m_deltaPhi_flat_wedge);
x3 = m_outerR * std::cos(-m_deltaPhiFlat/2.0 + (i+1)*m_deltaPhiFlatWedge);
y3 = m_outerR * std::sin(-m_deltaPhiFlat/2.0 + (i+1)*m_deltaPhiFlatWedge);
z3 = 0.0;
m_wedgeCoords_flat[i][3].SetVectorCartesian(x3, y3, z3);
m_flatWedgeCoords[i][3].SetXYZ(x3, y3, z3);
}
//Generate tilted rings
for(int i=0; i<m_nRings; i++) {
for(int j=0; j<4; j++) {
m_ringCoords_tilt[i][j] = TransformToTiltedFrame(m_ringCoords_flat[i][j]);
}
for(int i=0; i<s_nRings; i++)
{
for(int j=0; j<4; j++)
m_tiltRingCoords[i][j] = Transform(m_flatRingCoords[i][j]);
}
//Generate tilted wedges
for(int i=0; i<m_nWedges; i++) {
for(int j=0; j<4; j++) {
m_wedgeCoords_tilt[i][j] = TransformToTiltedFrame(m_wedgeCoords_flat[i][j]);
}
for(int i=0; i<s_nWedges; i++)
{
for(int j=0; j<4; j++)
m_tiltWedgeCoords[i][j] = Transform(m_flatWedgeCoords[i][j]);
}
}
@ -193,35 +200,37 @@ void SabreDetector::CalculateCorners() {
!NOTE: This currently only applies to a configuration where there is no translation in x & y. The math becomes significantly messier in these cases.
Also, don't use tan(). It's behavior near PI/2 makes it basically useless for these.
*/
Mask::Vec3 SabreDetector::GetTrajectoryCoordinates(double theta, double phi) {
if(m_translation.GetX() != 0.0 || m_translation.GetY() != 0.0) {
return Mask::Vec3();
}
ROOT::Math::XYZPoint SabreDetector::GetTrajectoryCoordinates(double theta, double phi)
{
if(m_translation.Vect().X() != 0.0 || m_translation.Vect().Y() != 0.0)
return ROOT::Math::XYZPoint();
//Calculate the *potential* phi in the flat detector
double phi_numerator = std::cos(m_tilt)*(std::sin(phi)*std::cos(m_phiCentral) - std::sin(m_phiCentral)*std::cos(phi));
double phi_denominator = std::cos(m_phiCentral)*std::cos(phi) + std::sin(m_phiCentral)*std::sin(phi);
double phi_numerator = std::cos(m_tilt)*(std::sin(phi)*std::cos(m_centerPhi) - std::sin(m_centerPhi)*std::cos(phi));
double phi_denominator = std::cos(m_centerPhi)*std::cos(phi) + std::sin(m_centerPhi)*std::sin(phi);
double phi_flat = std::atan2(phi_numerator, phi_denominator);
if(phi_flat < 0) phi_flat += M_PI*2.0;
if(phi_flat < 0)
phi_flat += M_PI*2.0;
//Calculate the *potential* R in the flat detector
double r_numerator = m_translation.GetZ()*std::cos(phi)*std::sin(theta);
double r_denominator = std::cos(phi_flat)*std::cos(m_phiCentral)*std::cos(m_tilt)*std::cos(theta) - std::sin(phi_flat)*std::sin(m_phiCentral)*std::cos(theta) - std::cos(phi_flat)*std::sin(m_tilt)*std::cos(phi)*std::sin(theta);
double r_numerator = m_translation.Vect().Z()*std::cos(phi)*std::sin(theta);
double r_denominator = std::cos(phi_flat)*std::cos(m_centerPhi)*std::cos(m_tilt)*std::cos(theta) -
std::sin(phi_flat)*std::sin(m_centerPhi)*std::cos(theta) -
std::cos(phi_flat)*std::sin(m_tilt)*std::cos(phi)*std::sin(theta);
double r_flat = r_numerator/r_denominator;
//Calculate the distance from the origin to the hit on the detector
double R_to_detector = (r_flat*std::cos(phi_flat)*std::sin(m_tilt) + m_translation.GetZ())/std::cos(theta);
double R_to_detector = (r_flat*std::cos(phi_flat)*std::sin(m_tilt) + m_translation.Vect().Z())/std::cos(theta);
double xhit = R_to_detector*std::sin(theta)*std::cos(phi);
double yhit = R_to_detector*std::sin(theta)*std::sin(phi);
double zhit = R_to_detector*std::cos(theta);
//Check to see if our flat coords fall inside the flat detector
if(IsInside(r_flat, phi_flat)) {
return Mask::Vec3(xhit, yhit, zhit);
} else {
return Mask::Vec3();
}
if(IsInside(r_flat, phi_flat))
return ROOT::Math::XYZPoint(xhit, yhit, zhit);
else
return ROOT::Math::XYZPoint();
}
/*
@ -241,20 +250,23 @@ Mask::Vec3 SabreDetector::GetTrajectoryCoordinates(double theta, double phi) {
!NOTE: This currently only applies to a configuration where there is no translation in x & y. The math becomes significantly messier in these cases.
Also, don't use tan(). It's behavior near PI/2 makes it basically useless for these.
*/
std::pair<int, int> SabreDetector::GetTrajectoryRingWedge(double theta, double phi) {
if(m_translation.GetX() != 0.0 || m_translation.GetY() != 0.0) {
std::pair<int, int> SabreDetector::GetTrajectoryRingWedge(double theta, double phi)
{
if(m_translation.Vect().X() != 0.0 || m_translation.Vect().Y() != 0.0)
return std::make_pair(-1, -1);
}
//Calculate the *potential* phi in the flat detector
double phi_numerator = std::cos(m_tilt)*(std::sin(phi)*std::cos(m_phiCentral) - std::sin(m_phiCentral)*std::cos(phi));
double phi_denominator = std::cos(m_phiCentral)*std::cos(phi) + std::sin(m_phiCentral)*std::sin(phi);
double phi_numerator = std::cos(m_tilt)*(std::sin(phi)*std::cos(m_centerPhi) - std::sin(m_centerPhi)*std::cos(phi));
double phi_denominator = std::cos(m_centerPhi)*std::cos(phi) + std::sin(m_centerPhi)*std::sin(phi);
double phi_flat = std::atan2(phi_numerator, phi_denominator);
if(phi_flat < 0) phi_flat += M_PI*2.0;
if(phi_flat < 0)
phi_flat += M_PI*2.0;
//Calculate the *potential* R in the flat detector
double r_numerator = m_translation.GetZ()*std::cos(phi)*std::sin(theta);
double r_denominator = std::cos(phi_flat)*std::cos(m_phiCentral)*std::cos(m_tilt)*std::cos(theta) - std::sin(phi_flat)*std::sin(m_phiCentral)*std::cos(theta) - std::cos(phi_flat)*std::sin(m_tilt)*std::cos(phi)*std::sin(theta);
double r_numerator = m_translation.Vect().Z()*std::cos(phi)*std::sin(theta);
double r_denominator = std::cos(phi_flat)*std::cos(m_centerPhi)*std::cos(m_tilt)*std::cos(theta) -
std::sin(phi_flat)*std::sin(m_centerPhi)*std::cos(theta) -
std::cos(phi_flat)*std::sin(m_tilt)*std::cos(phi)*std::sin(theta);
double r_flat = r_numerator/r_denominator;
//Calculate the distance from the origin to the hit on the detector
@ -262,31 +274,41 @@ std::pair<int, int> SabreDetector::GetTrajectoryRingWedge(double theta, double p
//Check to see if our flat coords fall inside the flat detector
if(IsInside(r_flat, phi_flat)) {
if(IsInside(r_flat, phi_flat))
{
int ringch, wedgech;
if(phi_flat > M_PI) phi_flat -= 2.0*M_PI; //Need phi in terms of [-deltaPhi_flat/2, deltaPhi_flat/2]
for(int i=0; i<m_nRings; i++) {
if(IsRingTopEdge(r_flat, i) || IsRingBottomEdge(r_flat, i)) { //If it falls in the interstrip spacing, kill it
if(phi_flat > M_PI)
phi_flat -= 2.0*M_PI; //Need phi in terms of [-deltaPhi_flat/2, deltaPhi_flat/2]
for(int i=0; i<s_nRings; i++)
{
if(IsRingTopEdge(r_flat, i) || IsRingBottomEdge(r_flat, i))
{ //If it falls in the interstrip spacing, kill it
ringch = -1;
break;
} else if(IsRing(r_flat, i)) {
}
else if(IsRing(r_flat, i))
{
ringch = i;
break;
}
}
for(int i=0; i<m_nWedges; i++) {
if(IsWedgeTopEdge(phi_flat, i) || IsWedgeBottomEdge(phi_flat, i)){ //If it falls in the interstrip spacing, kill it
for(int i=0; i<s_nWedges; i++)
{
if(IsWedgeTopEdge(phi_flat, i) || IsWedgeBottomEdge(phi_flat, i))
{ //If it falls in the interstrip spacing, kill it
wedgech = -1;
break;
} else if(IsWedge(phi_flat, i)) {
}
else if(IsWedge(phi_flat, i))
{
wedgech = i;
break;
}
}
return std::make_pair(ringch, wedgech);
} else {
return std::make_pair(-1,-1);
}
else
return std::make_pair(-1,-1);
}
/*
@ -295,18 +317,18 @@ std::pair<int, int> SabreDetector::GetTrajectoryRingWedge(double theta, double p
randomly wiggle the point within the pixel. Method intended for use with data, or
to smear out simulated data to mimic real data.
*/
Mask::Vec3 SabreDetector::GetHitCoordinates(int ringch, int wedgech) {
if(!CheckRingChannel(ringch) || !CheckWedgeChannel(wedgech)) {
return Mask::Vec3();
}
ROOT::Math::XYZPoint SabreDetector::GetHitCoordinates(int ringch, int wedgech)
{
if(!CheckRingChannel(ringch) || !CheckWedgeChannel(wedgech))
return ROOT::Math::XYZPoint();
double r_center = m_Rinner + (0.5+ringch)*m_deltaR_flat_ring;
double phi_center = -m_deltaPhi_flat/2.0 + (0.5+wedgech)*m_deltaPhi_flat_wedge;
double r_center = m_innerR + (0.5+ringch)*m_deltaRFlatRing;
double phi_center = -m_deltaPhiFlat/2.0 + (0.5+wedgech)*m_deltaPhiFlatWedge;
double x = r_center*std::cos(phi_center);
double y = r_center*std::sin(phi_center);
double z = 0;
Mask::Vec3 hit(x, y, z);
ROOT::Math::XYZPoint hit(x, y, z);
return TransformToTiltedFrame(hit);
return Transform(hit);
}

189
src/Detectors/SabreDetector.h Normal file
View File

@ -0,0 +1,189 @@
/*
Class which represents a single MMM detector in the SABRE array at FSU. Origial code by KGH, re-written by
GWM.
Distances in meters, angles in radians.
The channel arrays have four points, one for each corner. The corners are
as follows, as if looking BACK along beam (i.e. from the target's pov):
0---------------------1
| |
| | x
| | <-----
| | |
| | |
3---------------------2 y
(z is hence positive along beam direction)
The channel numbers, also as looking back from target pov, are:
>> rings are 0 -- 15 from inner to outer:
15 -------------------
14 -------------------
13 -------------------
.
.
.
2 -------------------
1 -------------------
0 -------------------
>> wedges are 0 -- 7 moving counterclockwise:
7 6 ... 1 0
| | | | | |
| | | | | |
| | | | | |
| | | | | |
| | | | | |
| | | | | |
>> Note that the detector starts centered on the x-axis (central phi = 0) untilted, and then is rotated to wherever the frick
it is supposed to go; phi = 90 is centered on y axis, pointing down towards the bottom of the scattering chamber
-- GWM, Dec 2020; based on the og code from kgh
*/
#ifndef SABREDETECTOR_H
#define SABREDETECTOR_H
#include <vector>
#include <cmath>
#include "Math/Point3D.h"
#include "Math/Vector3D.h"
#include "Math/RotationZ.h"
#include "Math/RotationY.h"
#include "Math/Translation3D.h"
class SabreDetector {
public:
SabreDetector();
SabreDetector(int detID, double Rin, double Rout, double deltaPhi_flat, double phiCentral, double tiltFromVert, double zdist, double xdist=0, double ydist=0);
~SabreDetector();
/*Return coordinates of the corners of each ring/wedge in SABRE*/
const ROOT::Math::XYZPoint& GetRingFlatCoords(int ch, int corner) const { return m_flatRingCoords[ch][corner]; }
const ROOT::Math::XYZPoint& GetWedgeFlatCoords(int ch, int corner) const { return m_flatWedgeCoords[ch][corner]; }
const ROOT::Math::XYZPoint& GetRingTiltCoords(int ch, int corner) const { return m_tiltRingCoords[ch][corner]; }
const ROOT::Math::XYZPoint& GetWedgeTiltCoords(int ch, int corner) const { return m_tiltWedgeCoords[ch][corner]; }
ROOT::Math::XYZPoint GetTrajectoryCoordinates(double theta, double phi);
std::pair<int, int> GetTrajectoryRingWedge(double theta, double phi);
ROOT::Math::XYZPoint GetHitCoordinates(int ringch, int wedgech);
int GetNumberOfWedges() { return s_nWedges; }
int GetNumberOfRings() { return s_nRings; }
ROOT::Math::XYZVector GetNormTilted() { return m_zRotation*(m_yRotation*m_norm); }
int GetDetectorID() { return m_detectorID; }
private:
void CalculateCorners();
/*Performs the transformation to the tilted,rotated,translated frame of the SABRE detector*/
ROOT::Math::XYZPoint Transform(const ROOT::Math::XYZPoint& vector) { return m_translation*(m_zRotation*(m_yRotation*vector)); };
/*Determine if a given channel/corner combo is valid*/
bool CheckRingChannel(int ch) { return (ch<s_nRings && ch>=0) ? true : false; };
bool CheckWedgeChannel(int ch) { return (ch<s_nWedges && ch >=0) ? true : false; };
bool CheckCorner(int corner) { return (corner < 4 && corner >=0) ? true : false; };
bool CheckRingLocation(int ch, int corner) { return CheckRingChannel(ch) && CheckCorner(corner); };
bool CheckWedgeLocation(int ch, int corner) { return CheckWedgeChannel(ch) && CheckCorner(corner); };
/*
For all of the calculations, need a limit precision to determine if values are actually equal or not
Here the approx. size of the strip spacing is used as the precision.
*/
bool CheckPositionEqual(double val1,double val2) { return fabs(val1-val2) > s_positionTol ? false : true; };
bool CheckAngleEqual(double val1,double val2) { return fabs(val1-val2) > s_angularTol ? false : true; };
/*Determine if a hit is within the bulk detector*/
bool IsInside(double r, double phi)
{
double phi_1 = m_deltaPhiFlat/2.0;
double phi_2 = M_PI*2.0 - m_deltaPhiFlat/2.0;
return (((r > m_innerR && r < m_outerR) || CheckPositionEqual(r, m_innerR) || CheckPositionEqual(r, m_outerR))
&& (phi > phi_2 || phi < phi_1 || CheckAngleEqual(phi, phi_1) || CheckAngleEqual(phi, phi_2)));
};
/*
For a given radius/phi are you inside of a given ring/wedge channel,
or are you on the spacing between these channels
*/
bool IsRing(double r, int ringch)
{
double ringtop = m_innerR + m_deltaRFlatRing*(ringch + 1);
double ringbottom = m_innerR + m_deltaRFlatRing*(ringch);
return (r>ringbottom && r<ringtop);
};
inline bool IsRingTopEdge(double r, int ringch)
{
double ringtop = m_innerR + m_deltaRFlatRing*(ringch + 1);
return CheckPositionEqual(r, ringtop);
};
inline bool IsRingBottomEdge(double r, int ringch)
{
double ringbottom = m_innerR + m_deltaRFlatRing*(ringch);
return CheckPositionEqual(r, ringbottom);
};
inline bool IsWedge(double phi, int wedgech)
{
double wedgetop = -m_deltaPhiFlat/2.0 + m_deltaPhiFlatWedge*(wedgech+1);
double wedgebottom = -m_deltaPhiFlat/2.0 + m_deltaPhiFlatWedge*(wedgech);
return ((phi>wedgebottom && phi<wedgetop));
};
inline bool IsWedgeTopEdge(double phi, int wedgech)
{
double wedgetop = -m_deltaPhiFlat/2.0 + m_deltaPhiFlatWedge*(wedgech+1);
return CheckAngleEqual(phi, wedgetop);
}
inline bool IsWedgeBottomEdge(double phi, int wedgech)
{
double wedgebottom = -m_deltaPhiFlat/2.0 + m_deltaPhiFlatWedge*(wedgech);
return CheckAngleEqual(phi, wedgebottom);
}
/*Class constants*/
static constexpr int s_nRings = 16;
static constexpr int s_nWedges = 8;
static constexpr double s_deg2rad = M_PI/180.0;
/*These are implicitly the width of the spacing between detector active strips*/
static constexpr double s_positionTol = 0.0001; //0.1 mm position tolerance
static constexpr double s_angularTol = 0.1*M_PI/180.0; // 0.1 degree angular tolerance
/*Class data*/
double m_outerR;
double m_innerR;
double m_deltaRFlat;
double m_deltaRFlatRing;
double m_deltaPhiFlat;
double m_deltaPhiFlatWedge;
double m_centerPhi;
double m_tilt;
ROOT::Math::Translation3D m_translation;
ROOT::Math::RotationY m_yRotation;
ROOT::Math::RotationZ m_zRotation;
ROOT::Math::XYZVector m_norm;
int m_detectorID;
std::vector<std::vector<ROOT::Math::XYZPoint>> m_flatRingCoords, m_flatWedgeCoords;
std::vector<std::vector<ROOT::Math::XYZPoint>> m_tiltRingCoords, m_tiltWedgeCoords;
};
#endif

358
src/Detectors/SabreEfficiency.cpp
View File

@ -1,40 +1,44 @@
#include "SabreEfficiency.h"
#include "MaskFile.h"
#include <fstream>
#include <iostream>
#include <iomanip>
#include "TFile.h"
#include "TTree.h"
SabreEfficiency::SabreEfficiency() :
DetectorEfficiency(), deadlayer(DEADLAYER_THIN), sabre_eloss(SABRE_THICKNESS), degrader(DEGRADER_THICKNESS)
DetectorEfficiency(), m_deadlayerEloss({14}, {28}, {1}, s_deadlayerThickness),
m_detectorEloss({14}, {28}, {1}, s_detectorThickness), m_degraderEloss({73}, {181}, {1}, s_degraderThickness)
{
//Only 0, 1, 4 valid in degrader land
//detectors.emplace_back(0, INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI0*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
//detectors.emplace_back(1, INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI1*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
detectors.emplace_back(2, INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI2*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
detectors.emplace_back(3, INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI3*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
//detectors.emplace_back(4, INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI4*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
std::vector<int> dead_z = {14};
std::vector<int> dead_a = {28};
std::vector<int> dead_stoich = {1};
deadlayer.SetElements(dead_z, dead_a, dead_stoich);
sabre_eloss.SetElements(dead_z, dead_a, dead_stoich);
std::vector<int> deg_z = {73};
std::vector<int> deg_a = {181};
std::vector<int> deg_s = {1};
degrader.SetElements(deg_z, deg_a, deg_s);
for(int i=0; i<s_nDets; i++)
m_detectors.emplace_back(i, s_innerR, s_outerR, s_deltaPhi*s_deg2rad, s_centerPhiList[i]*s_deg2rad, s_tilt*s_deg2rad, s_zOffset);
//Only 0, 1, 4 valid in degrader land
m_degradedDetectors[0] = true;
m_degradedDetectors[1] = true;
m_degradedDetectors[2] = false;
m_degradedDetectors[3] = false;
m_degradedDetectors[4] = true;
//Choose who to look at right now. Usually switch on or off degraded/non-degraded.
m_activeDetectors[0] = false;
m_activeDetectors[1] = false;
m_activeDetectors[2] = true;
m_activeDetectors[3] = true;
m_activeDetectors[4] = false;
}
SabreEfficiency::~SabreEfficiency() {}
void SabreEfficiency::CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) {
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: "<<inputname<<std::endl;
std::cout<<"Running efficiency calculation..."<<std::endl;
if(!dmap.IsValid()) {
if(!m_deadMap.IsValid())
{
std::cerr<<"Unable to run SABRE Efficiency without a dead channel map."<<std::endl;
std::cerr<<"If you have no dead channels, simply make a file that's empty"<<std::endl;
std::cerr<<"Exiting."<<std::endl;
@ -42,137 +46,137 @@ void SabreEfficiency::CalculateEfficiency(const std::string& inputname, const st
}
Mask::MaskFile input(inputname, Mask::MaskFile::FileType::read);
Mask::MaskFile output(outputname, Mask::MaskFile::FileType::write);
TFile* input = TFile::Open(inputname.c_str(), "READ");
TFile* output = TFile::Open(outputname.c_str(), "RECREATE");
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;
TTree* intree = (TTree*) input->Get("SimTree");
std::vector<Mask::Nucleus>* dataHandle = new std::vector<Mask::Nucleus>();
intree->SetBranchAddress("nuclei", &dataHandle);
output->cd();
TTree* outtree = new TTree("SimTree", "SimTree");
outtree->Branch("nuclei", dataHandle);
input->cd();
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;
Mask::Nucleus nucleus;
intree->GetEntry(1);
std::vector<int> counts;
std::vector<int> coinc_counts;
switch(header.rxn_type) {
case Mask::RxnType::PureDecay:
counts.resize(3, 0);
coinc_counts.resize(1, 0);
break;
case Mask::RxnType::OneStepRxn:
counts.resize(4, 0);
coinc_counts.resize(1, 0);
break;
case Mask::RxnType::TwoStepRxn:
counts.resize(6, 0);
coinc_counts.resize(4, 0);
break;
case Mask::RxnType::ThreeStepRxn:
counts.resize(8, 0);
coinc_counts.resize(11, 0);
break;
counts.resize(dataHandle->size());
switch(counts.size())
{
case 3: coinc_counts.resize(1, 0); break;
case 4: coinc_counts.resize(1, 0); break;
case 6: coinc_counts.resize(4, 0); break;
case 8: coinc_counts.resize(11, 0); break;
default:
{
std::cerr<<"Bad reaction type at AnasenEfficiency::CalculateEfficiency (given value: "<<GetStringFromRxnType(header.rxn_type)<<"). Quiting..."<<std::endl;
input.Close();
output.Close();
std::cerr<<"Bad reaction type at AnasenEfficiency::CalculateEfficiency (given value: "<<counts.size()<<"). Quiting..."<<std::endl;
input->Close();
output->Close();
stats.close();
return;
}
}
int percent5 = header.nsamples*0.05;
int count = 0;
int npercent = 0;
uint64_t nentries = intree->GetEntries();
uint64_t percent5 = nentries*0.05;
uint64_t count = 0;
uint64_t npercent = 0;
while(true) {
if(++count == percent5) {//Show progress every 5%
Mask::Nucleus nucleus;
for(uint64_t i=0; i<nentries; i++)
{
intree->GetEntry(i);
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;
data.KE[i] = result.second;
counts[i]++;
} else if(data.detect_flag[i] == true) {
data.detect_flag[i] = false;
for(std::size_t j=0; j<dataHandle->size(); j++)
{
Mask::Nucleus& nucleus = (*dataHandle)[j];
DetectorResult result = IsSabre(nucleus);
if(result.detectFlag)
{
nucleus.isDetected = true;
nucleus.detectedKE = result.energy_deposited;
nucleus.detectedTheta = result.direction.Theta();
nucleus.detectedPhi = result.direction.Phi();
counts[j]++;
}
}
CountCoincidences(data, coinc_counts, header.rxn_type);
CountCoincidences(*dataHandle, coinc_counts);
output.WriteData(data);
outtree->Fill();
}
input->Close();
output->cd();
outtree->Write(outtree->GetName(), TObject::kOverwrite);
output->Close();
input.Close();
output.Close();
delete dataHandle;
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::setw(10)<<i<<"|"<<std::setw(10)<<((double)counts[i]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
}
stats<<"Coincidence Efficiency"<<std::endl;
stats<<"---------------------"<<std::endl;
if(header.rxn_type == Mask::RxnType::PureDecay)
if(counts.size() == 3)
{
stats<<std::setw(10)<<"1 + 2"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"1 + 2"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
}
else if(header.rxn_type == Mask::RxnType::OneStepRxn)
else if(counts.size() == 4)
{
stats<<std::setw(10)<<"2 + 3"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 3"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
}
else if(header.rxn_type == Mask::RxnType::TwoStepRxn)
else if(counts.size() == 6)
{
stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
}
else if(header.rxn_type == Mask::RxnType::ThreeStepRxn)
else if(counts.size() == 8)
{
stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[4]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[4]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[5]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[5]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[6]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[6]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[7]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[7]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[8]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[8]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[9]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[9]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[10]/header.nsamples)<<std::endl;
stats<<std::setw(10)<<"2 + 4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[10]/nentries)<<std::endl;
stats<<"---------------------"<<std::endl;
}
stats.close();
@ -182,58 +186,69 @@ void SabreEfficiency::CalculateEfficiency(const std::string& inputname, const st
std::cout<<"---------------------------------------------"<<std::endl;
}
void SabreEfficiency::DrawDetectorSystem(const std::string& filename) {
void SabreEfficiency::DrawDetectorSystem(const std::string& filename)
{
std::ofstream output(filename);
std::vector<double> ringxs, ringys, ringzs;
std::vector<double> wedgexs, wedgeys, wedgezs;
Mask::Vec3 coords;
for(int i=0; i<5; i++) {
for(int j=0; j<detectors[i].GetNumberOfRings(); j++) {
for(int k=0; k<4; k++) {
coords = detectors[i].GetRingTiltCoords(j, k);
ringxs.push_back(coords.GetX());
ringys.push_back(coords.GetY());
ringzs.push_back(coords.GetZ());
ROOT::Math::XYZPoint coords;
for(int i=0; i<s_nDets; i++)
{
for(int j=0; j<m_detectors[i].GetNumberOfRings(); j++)
{
for(int k=0; k<4; k++)
{
coords = m_detectors[i].GetRingTiltCoords(j, k);
ringxs.push_back(coords.X());
ringys.push_back(coords.Y());
ringzs.push_back(coords.Z());
}
}
}
for(int i=0; i<5; i++) {
for(int j=0; j<detectors[i].GetNumberOfWedges(); j++) {
for(int k=0; k<4; k++) {
coords = detectors[i].GetWedgeTiltCoords(j, k);
wedgexs.push_back(coords.GetX());
wedgeys.push_back(coords.GetY());
wedgezs.push_back(coords.GetZ());
for(int i=0; i<s_nDets; i++)
{
for(int j=0; j<m_detectors[i].GetNumberOfWedges(); j++)
{
for(int k=0; k<4; k++)
{
coords = m_detectors[i].GetWedgeTiltCoords(j, k);
wedgexs.push_back(coords.X());
wedgeys.push_back(coords.Y());
wedgezs.push_back(coords.Z());
}
}
}
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++) {
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++) {
for(unsigned int i=0; i<wedgexs.size(); i++)
output<<wedgexs[i]<<" "<<wedgeys[i]<<" "<<wedgezs[i]<<std::endl;
}
output.close();
}
double SabreEfficiency::RunConsistencyCheck() {
double SabreEfficiency::RunConsistencyCheck()
{
double theta, phi;
double npoints = 5.0*16.0*4.0;
int count=0;
for(int h=0; h<5; h++) {
for(int j=0; j<16; j++) {
for(int k=0; k<4; k ++) {
theta = detectors[h].GetRingTiltCoords(j, k).GetTheta();
phi = detectors[h].GetRingTiltCoords(j, k).GetPhi();
for(int i=0; i<5; i++) {
auto channels = detectors[i].GetTrajectoryRingWedge(theta, phi);
if(channels.first != -1) {
ROOT::Math::XYZPoint corner;
for(auto& detector : m_detectors)
{
for(int j=0; j<detector.GetNumberOfRings(); j++)
{
for(int k=0; k<4; k ++) //Check corners
{
corner = detector.GetRingTiltCoords(j, k);
for(int i=0; i<5; i++)
{
auto channels = m_detectors[i].GetTrajectoryRingWedge(corner.Theta(), corner.Phi());
if(channels.first != -1)
{
count++;
}
}
@ -245,105 +260,118 @@ double SabreEfficiency::RunConsistencyCheck() {
}
/*Returns if detected, as well as total energy deposited in SABRE*/
std::pair<bool,double> SabreEfficiency::IsSabre(Mask::Nucleus& nucleus) {
if(nucleus.GetKE() <= ENERGY_THRESHOLD) {
return std::make_pair(false, 0.0);
}
DetectorResult SabreEfficiency::IsSabre(Mask::Nucleus& nucleus)
{
DetectorResult observation;
if(nucleus.GetKE() <= s_energyThreshold)
return observation;
Mask::Vec3 coords;
double thetaIncident, eloss, e_deposited;
double thetaIncident;
double ke = 0.0;
for(auto& detector : detectors) {
auto chan = detector.GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi());
if(chan.first != -1 && chan.second != -1) {
if(dmap.IsDead(detector.GetDetectorID(), chan.first, 0) || dmap.IsDead(detector.GetDetectorID(), chan.second, 1)) break; //dead channel check
coords = detector.GetTrajectoryCoordinates(nucleus.GetTheta(), nucleus.GetPhi());
thetaIncident = std::acos(coords.Dot(detector.GetNormTilted())/(coords.GetR()));
//eloss = degrader.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
//ke = nucleus.GetKE() - eloss;
eloss = deadlayer.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), ke, M_PI - thetaIncident);
ke = nucleus.GetKE() - eloss;
//ke -= eloss;
if(ke <= ENERGY_THRESHOLD) break; //deadlayer check
e_deposited = sabre_eloss.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), ke, M_PI - thetaIncident);
return std::make_pair(true, e_deposited);
}
for(auto& detector : m_detectors)
{
if(!m_activeDetectors[detector.GetDetectorID()])
continue;
auto channel = detector.GetTrajectoryRingWedge(nucleus.vec4.Theta(), nucleus.vec4.Phi());
if(channel.first == -1 || channel.second == -1)
continue;
if(m_deadMap.IsDead(detector.GetDetectorID(), channel.first, 0) || m_deadMap.IsDead(detector.GetDetectorID(), channel.second, 1))
break; //dead channel check
observation.detectFlag = true;
observation.direction = detector.GetTrajectoryCoordinates(nucleus.vec4.Theta(), nucleus.vec4.Phi());
thetaIncident = std::acos(observation.direction.Dot(detector.GetNormTilted())/(observation.direction.R()));
//Energy loss
ke = nucleus.GetKE();
if(m_degradedDetectors[detector.GetDetectorID()])
ke -= m_degraderEloss.GetEnergyLossTotal(nucleus.Z, nucleus.A, ke, M_PI - thetaIncident);
ke -= m_deadlayerEloss.GetEnergyLossTotal(nucleus.Z, nucleus.A, ke, M_PI - thetaIncident);
if(ke <= s_energyThreshold)
break;
observation.det_name = "SABRE"+std::to_string(detector.GetDetectorID());
observation.energy_deposited = m_detectorEloss.GetEnergyLossTotal(nucleus.Z, nucleus.A, ke, M_PI - thetaIncident);
return observation;
}
return std::make_pair(false,0.0);
observation.detectFlag = false;
return observation;
}
void SabreEfficiency::CountCoincidences(const Mask::MaskFileData& data, std::vector<int>& counts, Mask::RxnType rxn_type) {
if (rxn_type == Mask::RxnType::PureDecay && data.detect_flag[1] && data.detect_flag[2])
void SabreEfficiency::CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts)
{
if (data.size() == 3 && data[1].isDetected && data[2].isDetected)
{
counts[0]++;
}
else if (rxn_type == Mask::RxnType::OneStepRxn && data.detect_flag[2] && data.detect_flag[3])
else if (data.size() == 4 && data[2].isDetected && data[3].isDetected)
{
counts[0]++;
}
else if(rxn_type == Mask::RxnType::TwoStepRxn)
else if(data.size() == 6)
{
if(data.detect_flag[2] && data.detect_flag[4])
if(data[2].isDetected && data[4].isDetected)
{
counts[0]++;
}
if(data.detect_flag[2] && data.detect_flag[5])
if(data[2].isDetected && data[5].isDetected)
{
counts[1]++;
}
if(data.detect_flag[4] && data.detect_flag[5])
if(data[4].isDetected && data[5].isDetected)
{
counts[2]++;
}
if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[5])
if(data[2].isDetected && data[4].isDetected && data[5].isDetected)
{
counts[3]++;
}
}
else if(rxn_type == Mask::RxnType::ThreeStepRxn)
else if(data.size() == 8)
{
if(data.detect_flag[2] && data.detect_flag[4])
if(data[2].isDetected && data[4].isDetected)
{
counts[0]++;
}
if(data.detect_flag[2] && data.detect_flag[6])
if(data[2].isDetected && data[6].isDetected)
{
counts[1]++;
}
if(data.detect_flag[2] && data.detect_flag[7])
if(data[2].isDetected && data[7].isDetected)
{
counts[2]++;
}
if(data.detect_flag[4] && data.detect_flag[6])
if(data[4].isDetected && data[6].isDetected)
{
counts[3]++;
}
if(data.detect_flag[4] && data.detect_flag[7])
if(data[4].isDetected && data[7].isDetected)
{
counts[4]++;
}
if(data.detect_flag[6] && data.detect_flag[7])
if(data[6].isDetected && data[7].isDetected)
{
counts[5]++;
}
if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[6])
if(data[2].isDetected && data[4].isDetected && data[6].isDetected)
{
counts[6]++;
}
if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[7])
if(data[2].isDetected && data[4].isDetected && data[7].isDetected)
{
counts[7]++;
}
if(data.detect_flag[2] && data.detect_flag[6] && data.detect_flag[7])
if(data[2].isDetected && data[6].isDetected && data[7].isDetected)
{
counts[8]++;
}
if(data.detect_flag[4] && data.detect_flag[6] && data.detect_flag[7])
if(data[4].isDetected && data[6].isDetected && data[7].isDetected)
{
counts[9]++;
}
if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[6] && data.detect_flag[7])
if(data[2].isDetected && data[4].isDetected && data[6].isDetected && data[7].isDetected)
{
counts[10]++;
}

51
src/Detectors/SabreEfficiency.h Normal file
View File

@ -0,0 +1,51 @@
#ifndef SABREEFFICIENCY_H
#define SABREEFFICIENCY_H
#include "DetectorEfficiency.h"
#include "SabreDetector.h"
#include "Target.h"
#include "SabreDeadChannelMap.h"
#include "Mask/Nucleus.h"
class SabreEfficiency : public DetectorEfficiency
{
public:
SabreEfficiency();
~SabreEfficiency();
void SetDeadChannelMap(const std::string& filename) { m_deadMap.LoadMapfile(filename); };
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:
DetectorResult IsSabre(Mask::Nucleus& nucleus);
void CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts);
std::vector<SabreDetector> m_detectors;
Mask::Target m_deadlayerEloss;
Mask::Target m_detectorEloss;
Mask::Target m_degraderEloss;
SabreDeadChannelMap m_deadMap;
bool m_activeDetectors[5];
bool m_degradedDetectors[5];
//Sabre constants
static constexpr double s_innerR = 0.0326;
static constexpr double s_outerR = 0.1351;
static constexpr double s_tilt = 40.0;
static constexpr double s_zOffset = -0.1245;
static constexpr double s_deltaPhi = 54.4; //delta phi for each det
static constexpr int s_nDets = 5;
static constexpr double s_centerPhiList[s_nDets] = {306.0, 18.0, 234.0, 162.0, 90.0};
static constexpr double s_deg2rad = M_PI/180.0;
static constexpr double s_deadlayerThickness = 50 * 1e-7 * 2.3296 * 1e6; // ug/cm^2 (50 nm thick * density)
static constexpr double s_detectorThickness = 500 * 1e-4 * 2.3926 * 1e6; // ug/cm^2 (500 um thick * density)
static constexpr double s_degraderThickness = 70.0 * 1.0e-4 * 16.69 * 1e6; //tantalum degrader (70 um thick)
static constexpr double s_energyThreshold = 0.25; //in MeV
};
#endif

154
src/Detectors/StripDetector.cpp
View File

@ -17,32 +17,34 @@
*/
StripDetector::StripDetector(int ns, double len, double wid, double cphi, double cz, double crho) :
m_norm_unrot(1.0,0.0,0.0), m_uniform_fraction(0.0, 1.0), rndmFlag(false)
m_norm(1.0,0.0,0.0), m_uniformFraction(0.0, 1.0), m_isSmearing(false)
{
num_strips = ns;
m_nStrips = ns;
length = std::fabs(len);
total_width = std::fabs(wid);
front_strip_width = total_width/num_strips;
back_strip_length = length/num_strips;
m_stripLength = std::fabs(len);
m_totalWidth = std::fabs(wid);
m_frontStripWidth = m_totalWidth/m_nStrips;
m_backStripLength = m_stripLength/m_nStrips;
while (cphi < 0) cphi += 2*M_PI;
center_phi = cphi;
center_z = cz;
center_rho = std::fabs(crho);
while (cphi < 0)
cphi += 2*M_PI;
m_centerPhi = cphi;
m_centerZ = cz;
m_centerRho = std::fabs(crho);
zRot.SetAngle(center_phi);
m_zRotation.SetAngle(m_centerPhi);
front_strip_coords.resize(num_strips);
back_strip_coords.resize(num_strips);
rotated_front_strip_coords.resize(num_strips);
rotated_back_strip_coords.resize(num_strips);
for(int i=0; i<num_strips; i++) {
front_strip_coords[i].resize(num_corners);
back_strip_coords[i].resize(num_corners);
rotated_front_strip_coords[i].resize(num_corners);
rotated_back_strip_coords[i].resize(num_corners);
m_frontStripCoords.resize(m_nStrips);
m_backStripCoords.resize(m_nStrips);
m_rotFrontStripCoords.resize(m_nStrips);
m_rotBackStripCoords.resize(m_nStrips);
for(int i=0; i<m_nStrips; i++)
{
m_frontStripCoords[i].resize(s_nCorners);
m_backStripCoords[i].resize(s_nCorners);
m_rotFrontStripCoords[i].resize(s_nCorners);
m_rotBackStripCoords[i].resize(s_nCorners);
}
CalculateCorners();
@ -50,99 +52,107 @@ StripDetector::StripDetector(int ns, double len, double wid, double cphi, double
StripDetector::~StripDetector() {}
void StripDetector::CalculateCorners() {
void StripDetector::CalculateCorners()
{
double y_min, y_max, z_min, z_max;
for (int s=0; s<num_strips; s++) {
y_max = total_width/2.0 - front_strip_width*s;
y_min = total_width/2.0 - front_strip_width*(s+1);
z_max = center_z + length/2.0;
z_min = center_z - length/2.0;
front_strip_coords[s][2] = Mask::Vec3(center_rho, y_max, z_max);
front_strip_coords[s][3] = Mask::Vec3(center_rho, y_max, z_min);
front_strip_coords[s][0] = Mask::Vec3(center_rho, y_min, z_max);
front_strip_coords[s][1] = Mask::Vec3(center_rho, y_min, z_min);
for (int s=0; s<m_nStrips; s++)
{
y_max = m_totalWidth/2.0 - m_frontStripWidth*s;
y_min = m_totalWidth/2.0 - m_frontStripWidth*(s+1);
z_max = m_centerZ + m_stripLength/2.0;
z_min = m_centerZ - m_stripLength/2.0;
m_frontStripCoords[s][2].SetXYZ(m_centerRho, y_max, z_max);
m_frontStripCoords[s][3].SetXYZ(m_centerRho, y_max, z_min);
m_frontStripCoords[s][0].SetXYZ(m_centerRho, y_min, z_max);
m_frontStripCoords[s][1].SetXYZ(m_centerRho, y_min, z_min);
z_max = (center_z - length/2.0) + (s+1)*back_strip_length;
z_min = (center_z - length/2.0) + (s)*back_strip_length;
y_max = total_width/2.0;
y_min = -total_width/2.0;
back_strip_coords[s][2] = Mask::Vec3(center_rho, y_max, z_max);
back_strip_coords[s][3] = Mask::Vec3(center_rho, y_max, z_min);
back_strip_coords[s][0] = Mask::Vec3(center_rho, y_min, z_max);
back_strip_coords[s][1] = Mask::Vec3(center_rho, y_min, z_min);
z_max = (m_centerZ - m_stripLength/2.0) + (s+1)*m_backStripLength;
z_min = (m_centerZ - m_stripLength/2.0) + (s)*m_backStripLength;
y_max = m_totalWidth/2.0;
y_min = -m_totalWidth/2.0;
m_backStripCoords[s][2].SetXYZ(m_centerRho, y_max, z_max);
m_backStripCoords[s][3].SetXYZ(m_centerRho, y_max, z_min);
m_backStripCoords[s][0].SetXYZ(m_centerRho, y_min, z_max);
m_backStripCoords[s][1].SetXYZ(m_centerRho, y_min, z_min);
}
for(int s=0; s<num_strips; s++) {
rotated_front_strip_coords[s][0] = zRot*front_strip_coords[s][0];
rotated_front_strip_coords[s][1] = zRot*front_strip_coords[s][1];
rotated_front_strip_coords[s][2] = zRot*front_strip_coords[s][2];
rotated_front_strip_coords[s][3] = zRot*front_strip_coords[s][3];
for(int s=0; s<m_nStrips; s++)
{
m_rotFrontStripCoords[s][0] = m_zRotation*m_frontStripCoords[s][0];
m_rotFrontStripCoords[s][1] = m_zRotation*m_frontStripCoords[s][1];
m_rotFrontStripCoords[s][2] = m_zRotation*m_frontStripCoords[s][2];
m_rotFrontStripCoords[s][3] = m_zRotation*m_frontStripCoords[s][3];
rotated_back_strip_coords[s][0] = zRot*back_strip_coords[s][0];
rotated_back_strip_coords[s][1] = zRot*back_strip_coords[s][1];
rotated_back_strip_coords[s][2] = zRot*back_strip_coords[s][2];
rotated_back_strip_coords[s][3] = zRot*back_strip_coords[s][3];
m_rotBackStripCoords[s][0] = m_zRotation*m_backStripCoords[s][0];
m_rotBackStripCoords[s][1] = m_zRotation*m_backStripCoords[s][1];
m_rotBackStripCoords[s][2] = m_zRotation*m_backStripCoords[s][2];
m_rotBackStripCoords[s][3] = m_zRotation*m_backStripCoords[s][3];
}
}
Mask::Vec3 StripDetector::GetHitCoordinates(int front_stripch, double front_strip_ratio) {
ROOT::Math::XYZPoint StripDetector::GetHitCoordinates(int front_stripch, double front_strip_ratio)
{
if (!ValidChannel(front_stripch) || !ValidRatio(front_strip_ratio)) return Mask::Vec3(0,0,0);
if (!ValidChannel(front_stripch) || !ValidRatio(front_strip_ratio))
return ROOT::Math::XYZPoint(0,0,0);
double y;
if(rndmFlag) {
y = -total_width/2.0 + (front_stripch + m_uniform_fraction(Mask::RandomGenerator::GetInstance().GetGenerator()))*front_strip_width;
} else {
y = -total_width/2.0 + (front_stripch+0.5)*front_strip_width;
}
if(m_isSmearing)
y = -m_totalWidth/2.0 + (front_stripch + m_uniformFraction(Mask::RandomGenerator::GetInstance().GetGenerator()))*m_frontStripWidth;
else
y = -m_totalWidth/2.0 + (front_stripch+0.5)*m_frontStripWidth;
//recall we're still assuming phi=0 det:
Mask::Vec3 coords(center_rho, y, front_strip_ratio*(length/2) + center_z);
ROOT::Math::XYZPoint coords(m_centerRho, y, front_strip_ratio*(m_stripLength/2) + m_centerZ);
//NOW rotate by appropriate phi
return zRot*coords;
return m_zRotation*coords;
}
StripHit StripDetector::GetChannelRatio(double theta, double phi) {
StripHit StripDetector::GetChannelRatio(double theta, double phi)
{
StripHit hit;
while (phi < 0) phi += 2*M_PI;
while (phi < 0)
phi += 2*M_PI;
//to make the math easier (and the same for each det), rotate the input phi
//BACKWARD by the phi of the det, s.t. we are looking along x-axis
phi -= center_phi;
phi -= m_centerPhi;
if (phi > M_PI) phi -= 2*M_PI;
if (phi > M_PI)
phi -= 2*M_PI;
//then we can check easily whether it even hit the detector in phi
double det_max_phi = atan2(total_width/2,center_rho);
double det_max_phi = std::atan2(m_totalWidth/2,m_centerRho);
double det_min_phi = -det_max_phi;
if (phi < det_min_phi || phi > det_max_phi) return hit;
if (phi < det_min_phi || phi > det_max_phi)
return hit;
//for theta it's not so simple, so we have to go through the typical plane-intersect method
//first thing's first: we have a fixed x for the entire detector plane:
double xhit = center_rho;
double xhit = m_centerRho;
//thus we find the corresponding y and z for that fixed x, given the input theta and phi:
double yhit = xhit*tan(phi);
double zhit = sqrt(xhit*xhit+yhit*yhit)/tan(theta);
for (int s=0; s<num_strips; s++) {
if (xhit >= front_strip_coords[s][0].GetX() && xhit <= front_strip_coords[s][0].GetX() && //Check min and max x (constant in flat)
yhit >= front_strip_coords[s][1].GetY() && yhit <= front_strip_coords[s][2].GetY() && //Check min and max y
zhit >= front_strip_coords[s][1].GetZ() && zhit <= front_strip_coords[s][0].GetZ()) //Check min and max z
for (int s=0; s<m_nStrips; s++) {
if (xhit >=m_frontStripCoords[s][0].X() && xhit <=m_frontStripCoords[s][0].X() && //Check min and max x (constant in flat)
yhit >=m_frontStripCoords[s][1].Y() && yhit <=m_frontStripCoords[s][2].Y() && //Check min and max y
zhit >=m_frontStripCoords[s][1].Z() && zhit <=m_frontStripCoords[s][0].Z()) //Check min and max z
{
hit.front_strip_index = s;
hit.front_ratio = (zhit-center_z)/(length/2);
hit.front_ratio = (zhit-m_centerZ)/(m_stripLength/2);
break;
}
}
for (int s=0; s<num_strips; s++) {
if (xhit >= back_strip_coords[s][0].GetX() && xhit <= back_strip_coords[s][0].GetX() && //Check min and max x (constant in flat)
yhit >= back_strip_coords[s][1].GetY() && yhit <= back_strip_coords[s][2].GetY() && //Check min and max y
zhit >= back_strip_coords[s][1].GetZ() && zhit <= back_strip_coords[s][0].GetZ()) //Check min and max z
for (int s=0; s<m_nStrips; s++) {
if (xhit >= m_backStripCoords[s][0].X() && xhit <= m_backStripCoords[s][0].X() && //Check min and max x (constant in flat)
yhit >= m_backStripCoords[s][1].Y() && yhit <= m_backStripCoords[s][2].Y() && //Check min and max y
zhit >= m_backStripCoords[s][1].Z() && zhit <= m_backStripCoords[s][0].Z()) //Check min and max z
{
hit.back_strip_index = s;
break;

82
src/Detectors/StripDetector.h Normal file
View File

@ -0,0 +1,82 @@
#ifndef STRIP_DETECTOR_H
#define STRIP_DETECTOR_H
// +z is along beam axis
// +y is vertically "downward" in the lab frame
//angles must be in radians, but distances can be whatever
//PROVIDED all input distances are the same
//Front strips from largest y to smallest y
//Back strips from lowest z to highest z
#include <cmath>
#include <vector>
#include "Math/Point3D.h"
#include "Math/Vector3D.h"
#include "Math/RotationZ.h"
#include "Mask/RandomGenerator.h"
struct StripHit
{
int front_strip_index=-1;
int back_strip_index=-1;
double front_ratio=0.0;
};
class StripDetector
{
public:
StripDetector(int ns, double len, double wid, double cphi, double cz, double crho);
~StripDetector();
const ROOT::Math::XYZPoint& GetFrontStripCoordinates(int stripch, int corner) const { return m_frontStripCoords[stripch][corner]; }
const ROOT::Math::XYZPoint& GetBackStripCoordinates(int stripch, int corner) const { return m_backStripCoords[stripch][corner]; }
const ROOT::Math::XYZPoint& GetRotatedFrontStripCoordinates(int stripch, int corner) const
{
return m_rotFrontStripCoords[stripch][corner];
}
const ROOT::Math::XYZPoint& GetRotatedBackStripCoordinates(int stripch, int corner) const
{
return m_rotBackStripCoords[stripch][corner];
}
ROOT::Math::XYZVector GetNormRotated() const { return m_zRotation * m_norm; }
void SetPixelSmearing(bool isSmearing) { m_isSmearing = isSmearing; }
ROOT::Math::XYZPoint GetHitCoordinates(int front_stripch, double front_strip_ratio);
StripHit GetChannelRatio(double theta, double phi);
private:
bool ValidChannel(int f) { return ((f >= 0 && f < m_nStrips) ? true : false); };
bool ValidRatio(double r) { return ((r >= -1 && r <= 1) ? true : false); };
void CalculateCorners();
int m_nStrips;
static constexpr int s_nCorners = 4;
double m_stripLength; //common to all strips, hence total
double m_totalWidth;
double m_frontStripWidth; //assuming equal widths
double m_backStripLength; //assuming equal widths
double m_centerPhi; //assuming det centered above x-axis (corresponds to zero phi)
double m_centerZ;
double m_centerRho; //perpendicular radius from axis
std::vector<std::vector<ROOT::Math::XYZPoint>> m_frontStripCoords, m_backStripCoords;
std::vector<std::vector<ROOT::Math::XYZPoint>> m_rotFrontStripCoords, m_rotBackStripCoords;
ROOT::Math::XYZVector m_norm;
ROOT::Math::RotationZ m_zRotation;
std::uniform_real_distribution<double> m_uniformFraction;
bool m_isSmearing;
};
#endif

20
src/Detectors/DetectorEfficiency.cpp → src/Detectors/main.cpp
View File

@ -4,9 +4,11 @@
#include <iostream>
#include <string>
int main(int argc, char** argv) {
int main(int argc, char** argv)
{
if(argc != 4) {
if(argc != 4)
{
std::cerr<<"Incorrect number of commandline arguments! Returning."<<std::endl;
return 1;
}
@ -25,18 +27,12 @@ int main(int argc, char** argv) {
/*
try {
AnasenEfficiency anasen;
std::string mapfile = "./etc/AnasenDeadChannels.txt";
anasen.SetDeadChannelMap(mapfile);
anasen.CalculateEfficiency(inputname, outputname, statsname);
AnasenEfficiency anasen;
std::string mapfile = "./etc/AnasenDeadChannels.txt";
anasen.SetDeadChannelMap(mapfile);
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_fixedZ.txt");
} catch(const std::exception& e) {
std::cerr<<"Error: "<<e.what()<<std::endl;
std::cerr<<"Terminating."<<std::endl;
return 1;
}
*/
return 0;
}

91
src/Mask/AngularDistribution.cpp
View File

@ -8,26 +8,28 @@
namespace Mask {
AngularDistribution::AngularDistribution() :
uniform_cosine_dist(-1.0, 1.0), uniform_prob_dist(0.0, 1.0), branchingRatio(1.0), L(0), isoFlag(true)
m_uniformCosineDist(-1.0, 1.0), m_uniformProbDist(0.0, 1.0), m_branchingRatio(1.0), m_L(0), m_isIsotropic(true)
{
}
AngularDistribution::AngularDistribution(const std::string& file) :
branchingRatio(1.0), L(0), isoFlag(true)
m_branchingRatio(1.0), m_L(0), m_isIsotropic(true)
{
ReadDistributionFile(file);
}
AngularDistribution::~AngularDistribution() {}
void AngularDistribution::ReadDistributionFile(const std::string& file) {
void AngularDistribution::ReadDistributionFile(const std::string& file)
{
if(file == "none" || file == "") {
L=0;
branchingRatio=1.0;
constants.clear();
constants.push_back(0.5);
isoFlag = true;
if(file == "none" || file == "")
{
m_L=0;
m_branchingRatio=1.0;
m_constants.clear();
m_constants.push_back(0.5);
m_isIsotropic = true;
return;
}
@ -36,67 +38,76 @@ namespace Mask {
int l;
double par;
if(!input.is_open()) {
if(!input.is_open())
{
std::cerr<<"Unable to open distribution file. All values reset to default."<<std::endl;
L=0;
branchingRatio=1.0;
constants.clear();
constants.push_back(0.5);
isoFlag = true;
m_L=0;
m_branchingRatio=1.0;
m_constants.clear();
m_constants.push_back(0.5);
m_isIsotropic = true;
return;
}
input>>junk>>l;
while(input>>junk) {
while(input>>junk)
{
input>>par;
constants.push_back(par);
m_constants.push_back(par);
}
input.close();
if(constants.size() != ((unsigned int) l+1)) {
std::cerr<<"Unexpected number of constants for given angular momentum! Expected "<<l+1<<" and given "<<constants.size()<<std::endl;
if(m_constants.size() != ((std::size_t) l+1))
{
std::cerr<<"Unexpected number of constants for given angular momentum! Expected "<<l+1<<" and given "<<m_constants.size()<<std::endl;
std::cerr<<"Setting all values to default."<<std::endl;
branchingRatio=1.0;
constants.clear();
constants.push_back(0.5);
isoFlag = true;
m_L=0;
m_branchingRatio=1.0;
m_constants.clear();
m_constants.push_back(0.5);
m_isIsotropic = true;
return;
}
//Total branching ratio
branchingRatio = constants[0]*2.0;
L = l;
m_branchingRatio = m_constants[0]*2.0;
m_L = l;
//Renormalize distribution such that total prob is 1.0.
//Test branching ratio to see if we "make" a decay particle,
//then use re-normalized distribution to pick an angle.
if(constants[0] < 0.5) {
double norm = 0.5/constants[0];
for(auto& value : constants)
if(m_constants[0] < 0.5)
{
double norm = 0.5/m_constants[0];
for(auto& value : m_constants)
value *= norm;
}
isoFlag = false;
m_isIsotropic = false;
}
double AngularDistribution::GetRandomCosTheta() {
if(isoFlag)
return uniform_cosine_dist(RandomGenerator::GetInstance().GetGenerator());
double AngularDistribution::GetRandomCosTheta()
{
if(m_isIsotropic)
return m_uniformCosineDist(RandomGenerator::GetInstance().GetGenerator());
double test, probability;
double costheta;
test = uniform_prob_dist(RandomGenerator::GetInstance().GetGenerator());
if(test > branchingRatio) return -10;
test = m_uniformProbDist(RandomGenerator::GetInstance().GetGenerator());
if(test > m_branchingRatio)
return -10;
do {
do
{
probability = 0.0;
costheta = uniform_cosine_dist(RandomGenerator::GetInstance().GetGenerator());
test = uniform_prob_dist(RandomGenerator::GetInstance().GetGenerator());
for(unsigned int i=0; i<constants.size(); i++)
probability += constants[i]*P_l(i*2, costheta);
} while(test > probability);
costheta = m_uniformCosineDist(RandomGenerator::GetInstance().GetGenerator());
test = m_uniformProbDist(RandomGenerator::GetInstance().GetGenerator());
for(std::size_t i=0; i<m_constants.size(); i++)
probability += m_constants[i]*P_l(i*2, costheta);
}
while(test > probability);
return costheta;
}

35
src/Mask/AngularDistribution.h Normal file
View File

@ -0,0 +1,35 @@
#ifndef ANGULARDISTRIBUTION_H
#define ANGULARDISTRIBUTION_H
#include <string>
#include <vector>
#include <random>
namespace Mask {
class AngularDistribution
{
public:
AngularDistribution();
AngularDistribution(const std::string& file);
~AngularDistribution();
void ReadDistributionFile(const std::string& file);
double GetRandomCosTheta();
int GetL() { return m_L; }
double GetBranchingRatio() { return m_branchingRatio; }
private:
bool IsIsotropic() { return m_isIsotropic; }
std::uniform_real_distribution<double> m_uniformCosineDist;
std::uniform_real_distribution<double> m_uniformProbDist;
double m_branchingRatio;
int m_L;
std::vector<double> m_constants;
bool m_isIsotropic;
};
}
#endif

42
src/Mask/CMakeLists.txt
View File

@ -1,30 +1,56 @@
add_library(MaskDict SHARED)
target_include_directories(MaskDict
PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
SYSTEM PUBLIC ${ROOT_INCLUDE_DIRS}
)
ROOT_GENERATE_DICTIONARY(mask_dict Nucleus.h LINKDEF LinkDef_Nucleus.h MODULE MaskDict)
target_sources(MaskDict PRIVATE Nucleus.h Nucleus.cpp MassLookup.h MassLookup.cpp)
target_link_libraries(MaskDict ${ROOT_LIBRARIES})
set_target_properties(MaskDict PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${MASK_LIBRARY_DIR})
add_custom_command(TARGET MaskDict POST_BUILD
COMMAND ${CMAKE_COMMAND} -E copy
${CMAKE_CURRENT_BINARY_DIR}/libMaskDict_rdict.pcm
${MASK_LIBRARY_DIR}/libMaskDict_rdict.pcm
)
add_library(Mask STATIC)
target_include_directories(Mask
PUBLIC ${MASK_INCLUDE_DIR}
${ROOT_INCLUDE_DIRS}
)
target_sources(Mask PRIVATE
AngularDistribution.cpp
AngularDistribution.h
DecaySystem.cpp
EnergyLoss.cpp
DecaySystem.h
LayeredTarget.cpp
LayeredTarget.h
LegendrePoly.cpp
LegendrePoly.h
MaskApp.cpp
MaskFile.cpp
MaskApp.h
MassLookup.cpp
Nucleus.cpp
MassLookup.h
OneStepSystem.cpp
OneStepSystem.h
RandomGenerator.cpp
RandomGenerator.h
Reaction.cpp
Reaction.h
ReactionSystem.cpp
Rotation.cpp
ReactionSystem.h
Stopwatch.cpp
Stopwatch.h
Target.cpp
Target.h
ThreeStepSystem.cpp
ThreeStepSystem.h
TwoStepSystem.cpp
Vec3.cpp
Vec4.cpp
TwoStepSystem.h
)
target_link_libraries(Mask catima)
set_target_properties(Mask PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${MASK_LIBRARY_DIR})
target_link_libraries(Mask catima MaskDict ${ROOT_LIBRARIES})
set_target_properties(Mask PROPERTIES ARCHIVE_OUTPUT_DIRECTORY ${MASK_LIBRARY_DIR})

77
src/Mask/DecaySystem.cpp
View File

@ -1,76 +1,81 @@
#include "DecaySystem.h"
#include "RandomGenerator.h"
#include <sstream>
namespace Mask {
DecaySystem::DecaySystem() :
ReactionSystem()
{
nuclei.resize(3);
m_nuclei.resize(3);
}
DecaySystem::DecaySystem(std::vector<int>& z, std::vector<int>& a) :
DecaySystem::DecaySystem(const std::vector<int>& z, const std::vector<int>& a) :
ReactionSystem()
{
nuclei.resize(3);
m_nuclei.resize(3);
SetNuclei(z, a);
}
DecaySystem::~DecaySystem() {}
bool DecaySystem::SetNuclei(std::vector<int>& z, std::vector<int>& a) {
if(z.size() != a.size() || z.size() != 2) {
bool DecaySystem::SetNuclei(const std::vector<int>& z, const 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]);
int zr = z[0] - z[1];
int ar = a[0] - a[1];
m_nuclei[0] = CreateNucleus(z[0], a[0]); //target
m_nuclei[1] = CreateNucleus(z[1], a[1]); //breakup1
m_nuclei[2] = CreateNucleus(zr, ar); //breakup2
m_step1.BindNuclei(&(m_nuclei[0]), nullptr, &(m_nuclei[1]), &(m_nuclei[2]));
SetSystemEquation();
return true;
}
const std::vector<Nucleus>& DecaySystem::GetNuclei()
std::vector<Nucleus>* DecaySystem::GetNuclei()
{
nuclei[0] = step1.GetTarget();
nuclei[1] = step1.GetEjectile();
nuclei[2] = step1.GetResidual();
return nuclei;
return &m_nuclei;
}
void DecaySystem::LinkTarget() {
step1.SetLayeredTarget(&target);
m_step1.SetLayeredTarget(&m_target);
rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA());
if(rxnLayer != -1) {
step1.SetRxnLayer(rxnLayer);
target_set_flag = true;
} else {
m_rxnLayer = m_target.FindLayerContaining(m_nuclei[0].Z, m_nuclei[0].A);
if(m_rxnLayer != m_target.GetNumberOfLayers())
{
m_step1.SetRxnLayer(m_rxnLayer);
m_isTargetSet = 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::SetSystemEquation()
{
std::stringstream stream;
stream << m_nuclei[0].isotopicSymbol << "->"
<< m_nuclei[1].isotopicSymbol << "+"
<< m_nuclei[2].isotopicSymbol;
m_sysEquation = stream.str();
}
void DecaySystem::RunSystem() {
void DecaySystem::RunSystem()
{
//Link up the target if it hasn't been done yet
if(!target_set_flag) {
if(!m_isTargetSet)
LinkTarget();
}
double rxnTheta = std::acos(decay1dist.GetRandomCosTheta());
double rxnTheta = std::acos(m_step1Distribution.GetRandomCosTheta());
double rxnPhi = (*m_phi1Range)(RandomGenerator::GetInstance().GetGenerator());
step1.SetPolarRxnAngle(rxnTheta);
step1.SetAzimRxnAngle(rxnPhi);
step1.TurnOnResidualEloss();
step1.Calculate();
m_step1.SetPolarRxnAngle(rxnTheta);
m_step1.SetAzimRxnAngle(rxnPhi);
m_step1.SetResidualEnergyLoss(true);
m_step1.Calculate();
}
}

36
src/Mask/DecaySystem.h Normal file
View File

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

0
include/KinematicsExceptions.h → src/Mask/KinematicsExceptions.h
View File

135
src/Mask/LayeredTarget.cpp
View File

@ -16,17 +16,16 @@ Written by G.W. McCann Aug. 2020
namespace Mask {
LayeredTarget::LayeredTarget() :
name(""), m_fractional_depth_dist(0.0, 1.0)
m_name(""), m_fractionalDepthDistribution(0.0, 1.0)
{
}
LayeredTarget::~LayeredTarget() {}
/*Add in a Target made of a compound defined by a set of Zs, As, Ss, and a thickness*/
void LayeredTarget::AddLayer(std::vector<int>& Z, std::vector<int>& A, std::vector<int>& stoich, double thickness) {
Target t(thickness);
t.SetElements(Z, A, stoich);
layers.push_back(t);
void LayeredTarget::AddLayer(const std::vector<int>& Z, const std::vector<int>& A, const std::vector<int>& stoich, double thickness)
{
m_layers.emplace_back(Z, A, stoich, thickness);
}
/*
@ -34,9 +33,10 @@ namespace Mask {
Calculates energy loss assuming that the reaction occurs in the middle of the target layer
Note that the layer order can matter!
*/
double LayeredTarget::GetProjectileEnergyLoss(int zp, int ap, double startEnergy, int rxnLayer, double angle) {
if(rxnLayer < 0 || ((unsigned int) rxnLayer) > layers.size()) {
double LayeredTarget::GetProjectileEnergyLoss(int zp, int ap, double startEnergy, std::size_t rxnLayer, double angle)
{
if(rxnLayer > m_layers.size())
{
std::cerr<<"Reaction layer in eloss calculation is not in range! Returning 0"<<std::endl;
return 0.0;
}
@ -44,13 +44,17 @@ namespace Mask {
double eloss = 0.0;
double newEnergy = startEnergy;
double frac;
for(int i=0; i<=rxnLayer; i++) {
if(i == rxnLayer) {
frac = m_fractional_depth_dist(RandomGenerator::GetInstance().GetGenerator());
eloss += layers[i].GetEnergyLossFractionalDepth(zp, ap, newEnergy, angle, frac);
for(std::size_t i=0; i<=rxnLayer; i++)
{
if(i == rxnLayer)
{
frac = m_fractionalDepthDistribution(RandomGenerator::GetInstance().GetGenerator());
eloss += m_layers[i].GetEnergyLossFractionalDepth(zp, ap, newEnergy, angle, frac);
newEnergy = startEnergy - eloss;
} else {
eloss += layers[i].GetEnergyLossTotal(zp, ap, newEnergy, angle);
}
else
{
eloss += m_layers[i].GetEnergyLossTotal(zp, ap, newEnergy, angle);
newEnergy = startEnergy-eloss;
}
}
@ -63,83 +67,96 @@ namespace Mask {
Calculates energy loss assuming that the reaction occurs in the middle of the target
Note that the layer order can matter!
*/
double LayeredTarget::GetEjectileEnergyLoss(int ze, int ae, double startEnergy, int rxnLayer, double angle) {
double LayeredTarget::GetEjectileEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle) {
if(rxnLayer < 0 || ((unsigned int) rxnLayer) > layers.size()) {
if(rxnLayer > m_layers.size())
{
std::cerr<<"Reaction layer in eloss calculation is not in range! Returning 0"<<std::endl;
return 0.0;
}
double eloss = 0.0;
if(angle < M_PI/2.0) {
RandomGenerator& gen = RandomGenerator::GetInstance();
if(angle < M_PI/2.0)
{
double newEnergy = startEnergy;
for(unsigned int i=rxnLayer; i<layers.size(); i++) {
if(i == ((unsigned int)rxnLayer)) {
eloss += layers[i].GetEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractional_depth_dist(RandomGenerator::GetInstance().GetGenerator()));
newEnergy = startEnergy - eloss;
} else {
eloss += layers[i].GetEnergyLossTotal(ze, ae, newEnergy, angle);
newEnergy = startEnergy - eloss;
}
eloss += m_layers[rxnLayer].GetEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractionalDepthDistribution(gen.GetGenerator()));
newEnergy = startEnergy - eloss;
if(rxnLayer == m_layers.size())
return eloss;
for(std::size_t i=rxnLayer; i<m_layers.size(); i++)
{
eloss += m_layers[i].GetEnergyLossTotal(ze, ae, newEnergy, angle);
newEnergy = startEnergy - eloss;
}
} else { //Travelling backwards through target
}
else
{ //Travelling backwards through target
double newEnergy = startEnergy;
for(int i=rxnLayer; i>=0; i--) {
if(i == ((unsigned int)rxnLayer)) {
eloss += layers[i].GetEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractional_depth_dist(RandomGenerator::GetInstance().GetGenerator()));
newEnergy = startEnergy - eloss;
} else {
eloss += layers[i].GetEnergyLossTotal(ze, ae, newEnergy, angle);
newEnergy = startEnergy - eloss;
}
eloss += m_layers[rxnLayer].GetEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractionalDepthDistribution(gen.GetGenerator()));
newEnergy = startEnergy - eloss;
if(rxnLayer == 0)
return eloss;
for(std::size_t i=rxnLayer-1; i>0; i--) //unsigned ints cant be less than 0
{
eloss += m_layers[i].GetEnergyLossTotal(ze, ae, newEnergy, angle);
newEnergy = startEnergy - eloss;
}
eloss += m_layers[0].GetEnergyLossTotal(ze, ae, newEnergy, angle);
newEnergy = startEnergy - eloss;
}
return eloss;
}
/*ReverseEnergyLoss version of GetEjectileEnergyLoss*/
double LayeredTarget::GetEjectileReverseEnergyLoss(int ze, int ae, double startEnergy, int rxnLayer, double angle) {
if(rxnLayer < 0 || ((unsigned int) rxnLayer) > layers.size()) {
double LayeredTarget::GetEjectileReverseEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle)
{
if(rxnLayer > m_layers.size())
{
std::cerr<<"Reaction layer in eloss calculation is not in range! Returning 0"<<std::endl;
return 0.0;
}
double eloss = 0.0;
if(angle < M_PI/2.0) {
RandomGenerator& gen = RandomGenerator::GetInstance();
if(angle < M_PI/2.0)
{
double newEnergy = startEnergy;
for(int i=(layers.size()-1); i>=rxnLayer; i--) {
if(i == rxnLayer) {
eloss += layers[i].GetReverseEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractional_depth_dist(RandomGenerator::GetInstance().GetGenerator()));
newEnergy = startEnergy + eloss;
} else {
eloss += layers[i].GetReverseEnergyLossTotal(ze, ae, newEnergy, angle);
newEnergy = startEnergy + eloss;
}
for(std::size_t i=(m_layers.size()-1); i>rxnLayer; i--)
{
eloss += m_layers[i].GetReverseEnergyLossTotal(ze, ae, newEnergy, angle);
newEnergy = startEnergy + eloss;
}
} else {
eloss += m_layers[rxnLayer].GetReverseEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractionalDepthDistribution(gen.GetGenerator()));
newEnergy = startEnergy + eloss;
}
else
{
double newEnergy = startEnergy;
for(int i=0; i <= rxnLayer; i++) {
if(i == rxnLayer) {
eloss += layers[i].GetReverseEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractional_depth_dist(RandomGenerator::GetInstance().GetGenerator()));
newEnergy = startEnergy + eloss;
} else {
eloss += layers[i].GetReverseEnergyLossTotal(ze, ae, newEnergy, angle);
newEnergy = startEnergy + eloss;
}
for(std::size_t i=0; i < rxnLayer; i++)
{
eloss += m_layers[i].GetReverseEnergyLossTotal(ze, ae, newEnergy, angle);
newEnergy = startEnergy + eloss;
}
eloss += m_layers[rxnLayer].GetReverseEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractionalDepthDistribution(gen.GetGenerator()));
newEnergy = startEnergy + eloss;
}
return eloss;
}
int LayeredTarget::FindLayerContaining(int Z, int A) {
for(unsigned int i=0; i<layers.size(); i++)
if(layers[i].ContainsElement(Z, A))
std::size_t LayeredTarget::FindLayerContaining(int Z, int A)
{
for(std::size_t i=0; i<m_layers.size(); i++)
{
if(m_layers[i].ContainsElement(Z, A))
return i;
}
return -1;
return m_layers.size();
}

28
include/LayeredTarget.h → src/Mask/LayeredTarget.h
View File

@ -21,25 +21,25 @@ Written by G.W. McCann Aug. 2020
namespace Mask {
class LayeredTarget {
class LayeredTarget
{
public:
LayeredTarget();
~LayeredTarget();
void AddLayer(std::vector<int>& Z, std::vector<int>& A, std::vector<int>& stoich, double thickness);
double GetProjectileEnergyLoss(int zp, int ap, double startEnergy, int rxnLayer, double angle);
double GetEjectileEnergyLoss(int ze, int ae, double startEnergy, int rxnLayer, double angle);
double GetEjectileReverseEnergyLoss(int ze, int ae, double startEnergy, int rxnLayer, double angle);
int FindLayerContaining(int Z, int A);
inline int GetNumberOfLayers() { return layers.size(); }
inline void SetName(std::string& n) { name = n; }
inline const Target& GetLayerInfo(int index) { return layers[index]; }
inline const std::string& GetName() { return name; }
void AddLayer(const std::vector<int>& Z, const std::vector<int>& A, const std::vector<int>& stoich, double thickness);
double GetProjectileEnergyLoss(int zp, int ap, double startEnergy, std::size_t rxnLayer, double angle);
double GetEjectileEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle);
double GetEjectileReverseEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle);
std::size_t FindLayerContaining(int Z, int A);
std::size_t GetNumberOfLayers() { return m_layers.size(); }
void SetName(std::string& n) { m_name = n; }
const Target& GetLayerInfo(int index) { return m_layers[index]; }
const std::string& GetName() { return m_name; }
private:
std::vector<Target> layers;
std::string name;
std::uniform_real_distribution<double> m_fractional_depth_dist;
std::vector<Target> m_layers;
std::string m_name;
std::uniform_real_distribution<double> m_fractionalDepthDistribution;
};
}

25
src/Mask/LegendrePoly.cpp
View File

@ -3,33 +3,38 @@
namespace Mask {
double P_l(int l, double x) {
if(l == 0) {
double P_l(int l, double x)
{
if(l == 0)
return 1.0;
} else if (l == 1) {
else if (l == 1)
return x;
} else {
else
return (2.0*l - 1.0)/l*x*P_l(l-1, x) - (l-1.0)/l*P_l(l-2, x);
}
}
double Normed_P_l_sq(int l, double x) {
double Normed_P_l_sq(int l, double x)
{
return (2.0*l+1.0)/2.0*std::pow(P_l(l, x), 2.0);
}
double P_0(double x) {
double P_0(double x)
{
return 1.0;
}
double P_1(double x) {
double P_1(double x)
{
return x;
}
double P_2(double x) {
double P_2(double x)
{
return 0.5*(3.0*x*x -1.0);
}
double P_l_ROOT(double* x, double* pars) {
double P_l_ROOT(double* x, double* pars)
{
return P_l(pars[0], x[0]);
}

2
include/LegendrePoly.h → src/Mask/LegendrePoly.h
View File

@ -4,7 +4,7 @@
namespace Mask {
double P_l(int l, double x);
double P_l_ROOT(double* x, double* pars);
double P_l_ROOT(double* x, double* pars); //Only for use with ROOT cern in plotting
double Normed_P_l_sq(int l, double x);
double P_0(double x);

7
src/Mask/LinkDef_Nucleus.h Normal file
View File

@ -0,0 +1,7 @@
#ifdef __CLING__
#pragma link C++ struct Mask::Nucleus+;
#pragma link C++ class std::vector<Mask::Nucleus>+;
#pragma link C++ class std::string+;
#endif

87
src/Mask/MaskApp.cpp
View File

@ -1,19 +1,21 @@
#include "MaskApp.h"
#include "MaskFile.h"
#include <fstream>
#include <iostream>
#include "TFile.h"
#include "TTree.h"
namespace Mask {
MaskApp::MaskApp() :
m_sys(nullptr)
m_system(nullptr)
{
std::cout<<"----------GWM Kinematics Simulation----------"<<std::endl;
}
MaskApp::~MaskApp()
{
if(m_sys) delete m_sys;
delete m_system;
}
bool MaskApp::LoadConfig(const std::string& filename)
@ -28,23 +30,23 @@ namespace Mask {
std::string junk;
getline(input, junk);
input>>junk>>m_outfile_name;
input>>junk>>m_outputName;
std::vector<int> avec, zvec, svec;
int z, a, s;
getline(input, junk);
getline(input, junk);
input>>junk>>junk;
m_rxn_type = GetRxnTypeFromString(junk);
m_rxnType = GetRxnTypeFromString(junk);
getline(input, junk);
getline(input, junk);
switch(m_rxn_type)
switch(m_rxnType)
{
case RxnType::PureDecay:
{
m_sys = new DecaySystem();
m_rxn_type = RxnType::PureDecay;
for(int i=0; i<2; i++) {
m_system = new DecaySystem();
for(int i=0; i<2; i++)
{
input>>z>>a;
avec.push_back(a);
zvec.push_back(z);
@ -53,20 +55,21 @@ namespace Mask {
}
case RxnType::OneStepRxn:
{
m_sys = new OneStepSystem();
m_rxn_type = RxnType::OneStepRxn;
for(int i=0; i<3; i++) {
m_system = new OneStepSystem();
for(int i=0; i<3; i++)
{
input>>z>>a;
avec.push_back(a);
zvec.push_back(z);
}
std::cout<<"here"<<std::endl;
break;
}
case RxnType::TwoStepRxn:
{
m_sys = new TwoStepSystem();
m_rxn_type = RxnType::TwoStepRxn;
for(int i=0; i<4; i++) {
m_system = new TwoStepSystem();
for(int i=0; i<4; i++)
{
input>>z>>a;
avec.push_back(a);
zvec.push_back(z);
@ -75,9 +78,9 @@ namespace Mask {
}
case RxnType::ThreeStepRxn:
{
m_sys = new ThreeStepSystem();
m_rxn_type = RxnType::ThreeStepRxn;
for(int i=0; i<5; i++) {
m_system = new ThreeStepSystem();
for(int i=0; i<5; i++)
{
input>>z>>a;
avec.push_back(a);
zvec.push_back(z);
@ -87,7 +90,7 @@ namespace Mask {
default:
return false;
}
m_sys->SetNuclei(zvec, avec);
m_system->SetNuclei(zvec, avec);
int nlayers;
double thickness;
@ -110,7 +113,7 @@ namespace Mask {
input>>z>>a>>s;
zvec.push_back(z); avec.push_back(a); svec.push_back(s);
}
m_sys->AddTargetLayer(zvec, avec, svec, thickness);
m_system->AddTargetLayer(zvec, avec, svec, thickness);
input>>junk;
}
std::cout<<"Reaction equation: "<<GetSystemName()<<std::endl;
@ -122,51 +125,51 @@ namespace Mask {
input>>junk>>m_nsamples;
input>>junk>>par1>>junk>>par2;
m_sys->SetBeamDistro(par1, par2);
m_system->SetBeamDistro(par1, par2);
input>>junk>>par1;
switch(m_rxn_type)
switch(m_rxnType)
{
case RxnType::PureDecay : break;
case RxnType::None : break;
case RxnType::OneStepRxn :
{
dynamic_cast<OneStepSystem*>(m_sys)->SetReactionThetaType(par1);
dynamic_cast<OneStepSystem*>(m_system)->SetReactionThetaType(par1);
break;
}
case RxnType::TwoStepRxn :
{
dynamic_cast<TwoStepSystem*>(m_sys)->SetReactionThetaType(par1);
dynamic_cast<TwoStepSystem*>(m_system)->SetReactionThetaType(par1);
break;
}
case RxnType::ThreeStepRxn :
{
dynamic_cast<ThreeStepSystem*>(m_sys)->SetReactionThetaType(par1);
dynamic_cast<ThreeStepSystem*>(m_system)->SetReactionThetaType(par1);
break;
}
}
input>>junk>>par1>>junk>>par2;
m_sys->SetTheta1Range(par1, par2);
m_system->SetTheta1Range(par1, par2);
input>>junk>>par1>>junk>>par2;
m_sys->SetPhi1Range(par1, par2);
m_system->SetPhi1Range(par1, par2);
input>>junk>>par1>>junk>>par2;
m_sys->SetExcitationDistro(par1, par2);
m_system->SetExcitationDistro(par1, par2);
input>>junk>>dfile1;
input>>junk>>dfile2;
switch(m_rxn_type)
switch(m_rxnType)
{
case RxnType::PureDecay : break;
case RxnType::None : break;
case RxnType::OneStepRxn :
case RxnType::PureDecay :
{
DecaySystem* this_sys = dynamic_cast<DecaySystem*>(m_sys);
DecaySystem* this_sys = dynamic_cast<DecaySystem*>(m_system);
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 RxnType::None : break;
case RxnType::OneStepRxn : break;
case RxnType::TwoStepRxn :
{
TwoStepSystem* this_sys = dynamic_cast<TwoStepSystem*>(m_sys);
TwoStepSystem* this_sys = dynamic_cast<TwoStepSystem*>(m_system);
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;
@ -174,7 +177,7 @@ namespace Mask {
}
case RxnType::ThreeStepRxn :
{
ThreeStepSystem* this_sys = dynamic_cast<ThreeStepSystem*>(m_sys);
ThreeStepSystem* this_sys = dynamic_cast<ThreeStepSystem*>(m_system);
this_sys->SetDecay1Distribution(dfile1);
this_sys->SetDecay2Distribution(dfile2);
std::cout<<"Decay1 angular momentum: "<<this_sys->GetDecay1AngularMomentum()<<" Decay2 angular momentum: "<<this_sys->GetDecay2AngularMomentum()<<std::endl;
@ -192,13 +195,14 @@ namespace Mask {
void MaskApp::Run() {
std::cout<<"Running simulation..."<<std::endl;
if(m_sys == nullptr)
if(m_system == nullptr)
{
return;
}
MaskFile output(m_outfile_name, MaskFile::FileType::write);
output.WriteHeader(m_rxn_type, m_nsamples);
TFile* output = TFile::Open(m_outputName.c_str(), "RECREATE");
TTree* tree = new TTree("SimTree", "SimTree");
tree->Branch("nuclei", m_system->GetNuclei());
//For progress tracking
uint32_t percent5 = 0.05*m_nsamples;
@ -214,11 +218,12 @@ namespace Mask {
std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
}
m_sys->RunSystem();
output.WriteData(m_sys->GetNuclei());
m_system->RunSystem();
tree->Fill();
}
output.Close();
tree->Write(tree->GetName(), TObject::kOverwrite);
output->Close();
std::cout<<std::endl;
std::cout<<"Complete."<<std::endl;

37
src/Mask/MaskApp.h Normal file
View File

@ -0,0 +1,37 @@
#ifndef MASKAPP_H
#define MASKAPP_H
#include "ReactionSystem.h"
#include "DecaySystem.h"
#include "OneStepSystem.h"
#include "TwoStepSystem.h"
#include "ThreeStepSystem.h"
#include "RxnType.h"
namespace Mask {
class MaskApp
{
public:
MaskApp();
~MaskApp();
bool LoadConfig(const std::string& filename);
bool SaveConfig(const std::string& filename);
int GetNumberOfSamples() const { return m_nsamples; }
const std::string GetSystemName() const { return m_system == nullptr ? "Invalid System" : m_system->GetSystemEquation(); }
const std::string GetOutputName() const { return m_outputName; }
const RxnType GetReactionType() const { return m_rxnType; }
void Run();
private:
ReactionSystem* m_system;
std::string m_outputName;
RxnType m_rxnType;
int m_nsamples;
};
}
#endif

50
src/Mask/MassLookup.cpp
View File

@ -10,38 +10,41 @@ Written by G.W. McCann Aug. 2020
*/
#include "MassLookup.h"
#include "KinematicsExceptions.h"
#include <sstream>
namespace Mask {
MassLookup::MassLookup() {
MassLookup* MassLookup::s_instance = new MassLookup();
MassLookup::MassLookup()
{
std::ifstream massfile("etc/mass.txt");
if(massfile.is_open()) {
std::string junk, A, element;
int Z;
KeyPair key;
if(massfile.is_open())
{
std::string junk, element;
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;
while(massfile>>junk)
{
massfile>>key.Z>>key.A>>element>>atomicMassBig>>atomicMassSmall;
isotopicMass = (atomicMassBig + atomicMassSmall*1e-6 - key.Z*electron_mass)*u_to_mev;
massTable[key.GetID()] = isotopicMass;
elementTable[key.GetID()] = std::to_string(key.A) + element;
}
} else {
throw MassFileException();
}
else
throw MassFileException();
}
MassLookup::~MassLookup() {
}
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);
double MassLookup::FindMass(int Z, int A)
{
KeyPair key({Z, A});
auto data = massTable.find(key.GetID());
if(data == massTable.end())
throw MassException();
@ -49,13 +52,14 @@ namespace Mask {
}
//returns element symbol
std::string MassLookup::FindSymbol(int Z, int A) {
auto data = elementTable.find(Z);
std::string MassLookup::FindSymbol(int Z, int A)
{
KeyPair key({Z, A});
auto data = elementTable.find(key.GetID());
if(data == elementTable.end())
throw MassException();
std::string fullsymbol = std::to_string(A) + data->second;
return fullsymbol;
return data->second;
}
}

27
include/MassLookup.h → src/Mask/MassLookup.h
View File

@ -18,22 +18,35 @@ Converted to true singleton to simplify usage -- Aug. 2021 GWM
namespace Mask {
class MassLookup {
class MassLookup
{
public:
struct KeyPair
{
uint32_t Z;
uint32_t A;
//Use szudzik pairing method to make unqiue key out of two unsigned ints. Use size_t as extra safety.
std::size_t GetID()
{
return Z >= A ? Z*Z + Z + A : A*A + Z;
}
};
~MassLookup();
double FindMass(int Z, int A);
double FindMassU(int Z, int A) { return FindMass(Z, A)/u_to_mev; }
std::string FindSymbol(int Z, int A);
static MassLookup& GetInstance() {
static MassLookup s_instance;
return s_instance;
}
static MassLookup& GetInstance() { return *s_instance; }
private:
MassLookup();
std::unordered_map<std::string, double> massTable;
std::unordered_map<int, std::string> elementTable;
static MassLookup* s_instance;
std::unordered_map<std::size_t, double> massTable;
std::unordered_map<std::size_t, std::string> elementTable;
//constants
static constexpr double u_to_mev = 931.4940954;

17
src/Mask/Nucleus.cpp Normal file
View File

@ -0,0 +1,17 @@
#include "Nucleus.h"
namespace Mask {
Nucleus CreateNucleus(int z, int a)
{
Nucleus nuc;
nuc.Z = z;
nuc.A = a;
nuc.groundStateMass = MassLookup::GetInstance().FindMass(z, a);
nuc.isotopicSymbol = MassLookup::GetInstance().FindSymbol(z, a);
nuc.vec4 = ROOT::Math::PxPyPzEVector(0., 0., 0., nuc.groundStateMass);
return nuc;
}
bool EnforceDictionaryLinked() { return true; }
}

59
src/Mask/Nucleus.h Normal file
View File

@ -0,0 +1,59 @@
/*
Nucleus.h
Nucleus is a derived class of Vec4. A nucleus is the kinematics is essentially a 4 vector with the
additional properties of the number of total nucleons (A), the number of protons (Z), a ground state mass,
an exctitation energy, and an isotopic symbol.
--GWM Jan 2021
*/
#ifndef NUCLEUS_H
#define NUCLEUS_H
#include <string>
#include <vector>
#include "Math/Vector4D.h"
#include "MassLookup.h"
namespace Mask {
struct Nucleus
{
void SetVec4Spherical(double theta, double phi, double p, double E)
{
vec4.SetPxPyPzE(std::sin(theta)*std::cos(phi)*p,
std::sin(theta)*std::sin(phi)*p,
std::cos(theta)*p,
E
);
}
double GetKE() const
{
return vec4.E() - vec4.M();
}
double GetExcitationEnergy() const
{
return vec4.M() - groundStateMass;
}
int Z = 0;
int A = 0;
double groundStateMass = 0.0;
std::string isotopicSymbol = "";
double thetaCM = 0.0;
ROOT::Math::PxPyPzEVector vec4;
bool isDetected = false;
double detectedKE = 0.0;
double detectedTheta = 0.0;
double detectedPhi = 0.0;
};
Nucleus CreateNucleus(int z, int a);
bool EnforceDictionaryLinked();
};
#endif

84
src/Mask/OneStepSystem.cpp
View File

@ -1,68 +1,76 @@
#include "OneStepSystem.h"
#include "RandomGenerator.h"
#include <sstream>
namespace Mask {
OneStepSystem::OneStepSystem() :
ReactionSystem()
{
nuclei.resize(4);
m_nuclei.resize(4);
}
OneStepSystem::OneStepSystem(std::vector<int>& z, std::vector<int>& a) :
OneStepSystem::OneStepSystem(const std::vector<int>& z, const std::vector<int>& a) :
ReactionSystem()
{
nuclei.resize(4);
m_nuclei.resize(4);
SetNuclei(z, a);
}
OneStepSystem::~OneStepSystem() {}
bool OneStepSystem::SetNuclei(std::vector<int>& z, std::vector<int>& a) {
if(z.size() != a.size() || z.size() != 3) {
bool OneStepSystem::SetNuclei(const std::vector<int>& z, const std::vector<int>& a)
{
if(z.size() != a.size() || z.size() != 3)
return false;
}
step1.SetNuclei(z[0], a[0], z[1], a[1], z[2], a[2]);
int zr = z[0] + z[1] - z[2];
int ar = a[0] + a[1] - a[2];
m_nuclei[0] = CreateNucleus(z[0], a[0]); //target
m_nuclei[1] = CreateNucleus(z[1], a[1]); //projectile
m_nuclei[2] = CreateNucleus(z[2], a[2]); //ejectile
m_nuclei[3] = CreateNucleus(zr, ar); //residual
m_step1.BindNuclei(&(m_nuclei[0]), &(m_nuclei[1]), &(m_nuclei[2]), &(m_nuclei[3]));
SetSystemEquation();
return true;
}
const std::vector<Nucleus>& OneStepSystem::GetNuclei()
std::vector<Nucleus>* OneStepSystem::GetNuclei()
{
nuclei[0] = step1.GetTarget();
nuclei[1] = step1.GetProjectile();
nuclei[2] = step1.GetEjectile();
nuclei[3] = step1.GetResidual();
return nuclei;
return &m_nuclei;
}
void OneStepSystem::LinkTarget() {
step1.SetLayeredTarget(&target);
void OneStepSystem::LinkTarget()
{
m_step1.SetLayeredTarget(&m_target);
rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA());
if(rxnLayer != -1) {
step1.SetRxnLayer(rxnLayer);
target_set_flag = true;
} else {
throw ReactionLayerException();
m_rxnLayer = m_target.FindLayerContaining(m_nuclei[0].Z, m_nuclei[0].A);
if(m_rxnLayer != m_target.GetNumberOfLayers())
{
m_step1.SetRxnLayer(m_rxnLayer);
m_isTargetSet = 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::SetSystemEquation()
{
std::stringstream stream;
stream << m_nuclei[0].isotopicSymbol << "("
<< m_nuclei[1].isotopicSymbol << ", "
<< m_nuclei[2].isotopicSymbol << ")"
<< m_nuclei[3].isotopicSymbol;
m_sysEquation = stream.str();
}
void OneStepSystem::RunSystem() {
//Link up the target if it hasn't been done yet
if(!target_set_flag) {
if(!m_isTargetSet)
{
LinkTarget();
}
@ -72,13 +80,13 @@ namespace Mask {
double rxnPhi = (*m_phi1Range)(RandomGenerator::GetInstance().GetGenerator());
double residEx = (*m_exDist)(RandomGenerator::GetInstance().GetGenerator());
step1.SetBeamKE(bke);
step1.SetPolarRxnAngle(rxnTheta);
step1.SetAzimRxnAngle(rxnPhi);
step1.SetExcitation(residEx);
m_step1.SetBeamKE(bke);
m_step1.SetPolarRxnAngle(rxnTheta);
m_step1.SetAzimRxnAngle(rxnPhi);
m_step1.SetExcitation(residEx);
step1.TurnOnResidualEloss();
step1.Calculate();
m_step1.SetResidualEnergyLoss(true);
m_step1.Calculate();
}
}

30
src/Mask/OneStepSystem.h Normal file
View File

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

8
src/Mask/RandomGenerator.cpp
View File

@ -1,11 +1,13 @@
#include "RandomGenerator.h"
namespace Mask {
RandomGenerator::RandomGenerator() {
RandomGenerator* RandomGenerator::s_instance = new RandomGenerator();
RandomGenerator::RandomGenerator()
{
std::random_device rd;
rng.seed(rd());
}
RandomGenerator::~RandomGenerator() {
}
RandomGenerator::~RandomGenerator() {}
}

24
src/Mask/RandomGenerator.h Normal file
View File

@ -0,0 +1,24 @@
#ifndef RANDOMGENERATOR_H
#define RANDOMGENERATOR_H
#include <random>
namespace Mask {
class RandomGenerator
{
public:
~RandomGenerator();
std::mt19937& GetGenerator() { return rng; }
static RandomGenerator& GetInstance() { return *s_instance; }
private:
RandomGenerator();
static RandomGenerator* s_instance;
std::mt19937 rng;
};
}
#endif

295
src/Mask/Reaction.cpp
View File

@ -1,4 +1,4 @@
/*
/*ayer
Reaction.cpp
Reaction is a class which implements either a decay or scattering reaction. As such it requires either
3 (decay) or 4 (scattering) nuclei to perform any calcualtions. I also links together the target, which provides
@ -9,29 +9,33 @@
#include "Reaction.h"
#include "KinematicsExceptions.h"
#include "Math/Boost.h"
namespace Mask {
Reaction::Reaction() :
target(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), rxnLayer(0), m_eject_theta_type(lab), nuc_initFlag(false), resid_elossFlag(false)
m_target(nullptr), m_projectile(nullptr), m_ejectile(nullptr), m_residual(nullptr), m_layeredTarget(nullptr),
m_bke(0), m_theta(0), m_phi(0), m_ex(0), m_rxnLayer(0), m_ejectThetaType(s_lab), m_isInit(false), m_isResidEloss(false)
{
}
Reaction::Reaction(int zt, int at, int zp, int ap, int ze, int ae) :
target(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), rxnLayer(0), m_eject_theta_type(lab), resid_elossFlag(false)
Reaction::Reaction(Nucleus* target, Nucleus* projectile, Nucleus* ejectile, Nucleus* residual) :
m_target(nullptr), m_projectile(nullptr), m_ejectile(nullptr), m_residual(nullptr),
m_layeredTarget(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), m_rxnLayer(0), m_ejectThetaType(s_lab), m_isResidEloss(false)
{
SetNuclei(zt, at, zp, ap, ze, ae);
BindNuclei(target, projectile, ejectile, residual);
}
Reaction::~Reaction()
{
}
bool Reaction::Calculate() {
if(!nuc_initFlag)
bool Reaction::Calculate()
{
if(!m_isInit)
return false;
if(decayFlag) {
if(m_isDecay) {
CalculateDecay();
return true;
} else {
@ -40,203 +44,184 @@ namespace Mask {
}
}
//Deep copy of nucleus array
void Reaction::SetNuclei(const Nucleus* nucs) {
reactants[0] = nucs[0];
reactants[1] = nucs[1];
reactants[2] = nucs[2];
reactants[3] = nucs[3];
nuc_initFlag = true;
void Reaction::BindNuclei(Nucleus* target, Nucleus* projectile, Nucleus* ejectile, Nucleus* residual)
{
m_target = target;
m_projectile = projectile;
m_ejectile = ejectile;
m_residual = residual;
if(m_projectile == nullptr)
m_isDecay = true;
else
m_isDecay = false;
if(m_target == nullptr || m_ejectile == nullptr || m_residual == nullptr)
m_isInit = false;
else
m_isInit = true;
}
void Reaction::SetNuclei(int zt, int at, int zp, int ap, int ze, int ae) {
int zr, ar;
reactants[0] = Nucleus(zt, at);
reactants[2] = Nucleus(ze, ae);
if(ap == 0) {
decayFlag = true;
zr = zt - ze;
ar = at - ae;
} else {
reactants[1] = Nucleus(zp, ap);
decayFlag = false;
zr = zt + zp - ze;
ar = at + ap - ae;
}
if(zr < 0 || ar <= 0) {
nuc_initFlag = false;
} else {
reactants[3] = Nucleus(zr, ar);
nuc_initFlag = true;
}
}
void Reaction::SetBeamKE(double bke) {
if(!nuc_initFlag || decayFlag)
void Reaction::SetBeamKE(double bke)
{
if(!m_isInit || m_isDecay)
return;
m_bke = bke - target->GetProjectileEnergyLoss(reactants[1].GetZ(), reactants[1].GetA(), bke, rxnLayer, 0);
m_bke = bke - m_layeredTarget->GetProjectileEnergyLoss(m_projectile->Z, m_projectile->A, bke, m_rxnLayer, 0);
}
void Reaction::SetEjectileThetaType(int type) {
if(decayFlag) return;
if(type != center_of_mass && type != lab) return;
void Reaction::SetEjectileThetaType(int type)
{
if(m_isDecay)
return;
if(type != s_centerOfMass && type != s_lab)
return;
m_eject_theta_type = type;
m_ejectThetaType = type;
}
//Methods given by Iliadis in Nuclear Physics of Stars, Appendix C
//For use with lab frame restricted angles. May not give appropriate disribution for ejectile
void Reaction::CalculateReactionThetaLab() {
reactants[0].SetVectorCartesian(0.,0.,0.,reactants[0].GetGroundStateMass());
double beam_pz = std::sqrt(m_bke*(m_bke + 2.0 * reactants[1].GetGroundStateMass()));
double beam_E = m_bke + reactants[1].GetGroundStateMass();
reactants[1].SetVectorCartesian(0.,0.,beam_pz,beam_E);
void Reaction::CalculateReactionThetaLab()
{
m_target->vec4.SetPxPyPzE(0.,0.,0.,m_target->groundStateMass);
double beam_pz = std::sqrt(m_bke*(m_bke + 2.0 * m_projectile->groundStateMass));
double beam_E = m_bke + m_projectile->groundStateMass;
m_projectile->vec4.SetPxPyPzE(0.,0.,beam_pz,beam_E);
double Q = reactants[0].GetGroundStateMass() + reactants[1].GetGroundStateMass() - (reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() + m_ex);
double Q = m_target->groundStateMass + m_projectile->groundStateMass - (m_ejectile->groundStateMass + m_residual->groundStateMass + m_ex);
double Ethresh = -Q*(reactants[2].GetGroundStateMass()+reactants[3].GetGroundStateMass())/(reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() - reactants[1].GetGroundStateMass());
if(m_bke < Ethresh) {
double Ethresh = -Q*(m_ejectile->groundStateMass+m_residual->groundStateMass) /
(m_ejectile->groundStateMass + m_residual->groundStateMass - m_projectile->groundStateMass);
if(m_bke < Ethresh)
throw EnergyThresholdException();
}
double term1 = sqrt(reactants[1].GetGroundStateMass()*reactants[2].GetGroundStateMass()*m_bke)/(reactants[2].GetGroundStateMass()+reactants[3].GetGroundStateMass())*cos(m_theta);
double term2 = (m_bke*(reactants[3].GetGroundStateMass() - reactants[1].GetGroundStateMass()) + reactants[3].GetGroundStateMass()*Q)/(reactants[3].GetGroundStateMass() + reactants[2].GetGroundStateMass());
double sqrt_pos_ejectKE = term1 + sqrt(term1*term1 + term2);
double sqrt_neg_ejectKE = term1 - sqrt(term1*term1 + term2);
double term1 = std::sqrt(m_projectile->groundStateMass * m_ejectile->groundStateMass * m_bke)/
(m_ejectile->groundStateMass + m_residual->groundStateMass) * std::cos(m_theta);
double term2 = (m_bke * (m_residual->groundStateMass - m_projectile->groundStateMass) + m_residual->groundStateMass*Q) /
(m_residual->groundStateMass + m_ejectile->groundStateMass);
double sqrt_pos_ejectKE = term1 + std::sqrt(term1*term1 + term2);
double sqrt_neg_ejectKE = term1 - std::sqrt(term1*term1 + term2);
double ejectKE;
if(sqrt_pos_ejectKE > 0) {
if(sqrt_pos_ejectKE > 0)
ejectKE = sqrt_pos_ejectKE*sqrt_pos_ejectKE;
} else {
else
ejectKE = sqrt_neg_ejectKE*sqrt_neg_ejectKE;
}
double ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * reactants[2].GetGroundStateMass()));
double ejectE = ejectKE + reactants[2].GetGroundStateMass();
double ejectP = std::sqrt(ejectKE * (ejectKE + 2.0 * m_ejectile->groundStateMass));
double ejectE = ejectKE + m_ejectile->groundStateMass;
reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP, ejectE);
m_ejectile->SetVec4Spherical(m_theta, m_phi, ejectP, ejectE);
reactants[3] = reactants[0] + reactants[1] - reactants[2];
m_residual->vec4 = m_target->vec4 + m_projectile->vec4 - m_ejectile->vec4;
ejectKE -= target->GetEjectileEnergyLoss(reactants[2].GetZ(), reactants[2].GetA(), ejectKE, rxnLayer, m_theta);
ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * reactants[2].GetGroundStateMass()));
ejectE = ejectKE + reactants[2].GetGroundStateMass();
reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP, ejectE);
ejectKE -= m_layeredTarget->GetEjectileEnergyLoss(m_ejectile->Z, m_ejectile->A, ejectKE, m_rxnLayer, m_theta);
ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * m_ejectile->groundStateMass));
ejectE = ejectKE + m_ejectile->groundStateMass;
m_ejectile->SetVec4Spherical(m_theta, m_phi, ejectP, ejectE);
if(resid_elossFlag) {
double residKE = reactants[3].GetKE() - target->GetEjectileEnergyLoss(reactants[3].GetZ(), reactants[3].GetA(), reactants[3].GetKE(), rxnLayer, reactants[3].GetTheta());
double residP = std::sqrt(residKE*(residKE + 2.0*reactants[3].GetInvMass()));
double residE = residKE + reactants[3].GetInvMass();
reactants[3].SetVectorSpherical(reactants[3].GetTheta(), reactants[3].GetPhi(), residP, residE);
if(m_isResidEloss) {
double residKE = m_residual->GetKE() - m_layeredTarget->GetEjectileEnergyLoss(m_residual->Z, m_residual->A, m_residual->GetKE(),
m_rxnLayer, m_residual->vec4.Theta());
double residP = std::sqrt(residKE*(residKE + 2.0*m_residual->vec4.M()));
double residE = residKE + m_residual->vec4.M();
m_residual->SetVec4Spherical(m_residual->vec4.Theta(), m_residual->vec4.Phi(), residP, residE);
}
}
//Methods from original ANASEN. Gives proper distribution for inverse kinematics.
void Reaction::CalculateReactionThetaCM() {
void Reaction::CalculateReactionThetaCM()
{
//Target assumed at rest, with 0 excitation energy
reactants[0].SetVectorCartesian(0.,0.,0.,reactants[0].GetGroundStateMass());
double beam_pz = std::sqrt(m_bke*(m_bke + 2.0 * reactants[1].GetGroundStateMass()));
double beam_E = m_bke + reactants[1].GetGroundStateMass();
reactants[1].SetVectorCartesian(0.,0.,beam_pz,beam_E);
m_target->vec4.SetPxPyPzE(0.,0.,0.,m_target->groundStateMass);
double beam_pz = std::sqrt(m_bke*(m_bke + 2.0 * m_projectile->groundStateMass));
double beam_E = m_bke + m_projectile->groundStateMass;
m_projectile->vec4.SetPxPyPzE(0.,0.,beam_pz,beam_E);
double Q = m_target->groundStateMass + m_projectile->groundStateMass - (m_ejectile->groundStateMass + m_residual->groundStateMass + m_ex);
double Q = reactants[0].GetGroundStateMass() + reactants[1].GetGroundStateMass() - (reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() + m_ex);
double Ethresh = -Q*(reactants[2].GetGroundStateMass()+reactants[3].GetGroundStateMass())/(reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() - reactants[1].GetGroundStateMass());
if(m_bke < Ethresh) {
double Ethresh = -Q*(m_ejectile->groundStateMass + m_residual->groundStateMass) /
(m_ejectile->groundStateMass + m_residual->groundStateMass - m_projectile->groundStateMass);
if(m_bke < Ethresh)
throw EnergyThresholdException();
}
auto parent = reactants[0] + reactants[1];
double boost2lab[3];
double boost2cm[3];
const double* boost = parent.GetBoost();
for(int i=0; i<3; i++) {
boost2lab[i] = boost[i];
boost2cm[i] = boost[i]*(-1.0);
}
parent.ApplyBoost(boost2cm);
double ejectE_cm = (std::pow(reactants[2].GetGroundStateMass(), 2.0) - std::pow(reactants[3].GetGroundStateMass() + m_ex, 2.0) + std::pow(parent.GetE(),2.0))/(2.0*parent.GetE());
double ejectP_cm = std::sqrt(ejectE_cm*ejectE_cm - std::pow(reactants[2].GetGroundStateMass(), 2.0));
reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP_cm, ejectE_cm);
reactants[2].ApplyBoost(boost2lab);
reactants[3] = reactants[0] + reactants[1] - reactants[2];
ROOT::Math::PxPyPzEVector parent = m_target->vec4 + m_projectile->vec4;
ROOT::Math::Boost boost(parent.BoostToCM());
parent = boost*parent;
double ejectE_cm = (std::pow(m_ejectile->groundStateMass, 2.0) -
std::pow(m_residual->groundStateMass + m_ex, 2.0) + std::pow(parent.E(),2.0))/
(2.0*parent.E());
double ejectP_cm = std::sqrt(ejectE_cm*ejectE_cm - std::pow(m_ejectile->groundStateMass, 2.0));
m_ejectile->SetVec4Spherical(m_theta, m_phi, ejectP_cm, ejectE_cm);
m_ejectile->vec4 = boost.Inverse() * m_ejectile->vec4;
m_residual->vec4 = m_target->vec4 + m_projectile->vec4 - m_ejectile->vec4;
double ejectKE = reactants[2].GetKE();
double ejectP = reactants[2].GetP();
double ejectE = reactants[2].GetE();
double ejectKE = m_ejectile->GetKE();
double ejectP = m_ejectile->vec4.P();
double ejectE = m_ejectile->vec4.E();
//energy loss for ejectile (after reaction!)
ejectKE -= target->GetEjectileEnergyLoss(reactants[2].GetZ(), reactants[2].GetA(), ejectKE, rxnLayer, reactants[2].GetTheta());
ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * reactants[2].GetGroundStateMass()));
ejectE = ejectKE + reactants[2].GetGroundStateMass();
reactants[2].SetVectorSpherical(reactants[2].GetTheta(), reactants[2].GetPhi(), ejectP, ejectE);
ejectKE -= m_layeredTarget->GetEjectileEnergyLoss(m_ejectile->Z, m_ejectile->A, ejectKE, m_rxnLayer, m_ejectile->vec4.Theta());
ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * m_ejectile->groundStateMass));
ejectE = ejectKE + m_ejectile->groundStateMass;
m_ejectile->SetVec4Spherical(m_ejectile->vec4.Theta(), m_ejectile->vec4.Phi(), ejectP, ejectE);
//if on, get eloss for residual (after reaction!)
if(resid_elossFlag) {
double residKE = reactants[3].GetKE() - target->GetEjectileEnergyLoss(reactants[3].GetZ(), reactants[3].GetA(), reactants[3].GetKE(), rxnLayer, reactants[3].GetTheta());
double residP = std::sqrt(residKE*(residKE + 2.0*reactants[3].GetInvMass()));
double residE = residKE + reactants[3].GetInvMass();
reactants[3].SetVectorSpherical(reactants[3].GetTheta(), reactants[3].GetPhi(), residP, residE);
if(m_isResidEloss)
{
double residKE = m_residual->GetKE() -
m_layeredTarget->GetEjectileEnergyLoss(m_residual->Z, m_residual->A, m_residual->GetKE(), m_rxnLayer, m_residual->vec4.Theta());
double residP = std::sqrt(residKE*(residKE + 2.0*m_residual->vec4.M()));
double residE = residKE + m_residual->vec4.M();
m_residual->SetVec4Spherical(m_residual->vec4.Theta(), m_residual->vec4.Phi(), residP, residE);
}
}
void Reaction::CalculateReaction() {
switch(m_eject_theta_type) {
case center_of_mass:
{
CalculateReactionThetaCM();
break;
}
case lab:
{
CalculateReactionThetaLab();
break;
}
void Reaction::CalculateReaction()
{
switch(m_ejectThetaType)
{
case s_centerOfMass: CalculateReactionThetaCM(); break;
case s_lab: CalculateReactionThetaLab(); break;
}
}
//Calculate in CM, where decay is isotropic
void Reaction::CalculateDecay() {
double Q = reactants[0].GetInvMass() - reactants[2].GetGroundStateMass() - reactants[3].GetGroundStateMass();
if(Q < 0) {
void Reaction::CalculateDecay()
{
double Q = m_target->vec4.M() - m_ejectile->groundStateMass - m_residual->groundStateMass;
if(Q < 0)
throw QValueException();
}
const double* boost = reactants[0].GetBoost();
double boost2cm[3];
double boost2lab[3];
for(int i=0; i<3; i++) {
boost2lab[i] = boost[i];
boost2cm[i] = boost[i]*(-1.0);
}
ROOT::Math::Boost boost(m_target->vec4.BoostToCM());
m_target->vec4 = boost*m_target->vec4;
double ejectE_cm = (m_ejectile->groundStateMass*m_ejectile->groundStateMass -
m_residual->groundStateMass*m_residual->groundStateMass + m_target->vec4.E()*m_target->vec4.E()) /
(2.0*m_target->vec4.E());
double ejectP_cm = std::sqrt(ejectE_cm*ejectE_cm - m_ejectile->groundStateMass*m_ejectile->groundStateMass);
reactants[0].ApplyBoost(&(boost2cm[0]));
double ejectE_cm = (reactants[2].GetGroundStateMass()*reactants[2].GetGroundStateMass() - reactants[3].GetGroundStateMass()*reactants[3].GetGroundStateMass() + reactants[0].GetE()*reactants[0].GetE())/
(2.0*reactants[0].GetE());
double ejectP_cm = std::sqrt(ejectE_cm*ejectE_cm - reactants[2].GetGroundStateMass()*reactants[2].GetGroundStateMass());
m_ejectile->SetVec4Spherical(m_theta, m_phi, ejectP_cm, ejectE_cm);
m_ejectile->thetaCM = m_theta;
reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP_cm, ejectE_cm);
reactants[2].SetThetaCM(m_theta);
m_target->vec4 = boost.Inverse() * m_target->vec4;
m_ejectile->vec4 = boost.Inverse() * m_ejectile->vec4;
reactants[0].ApplyBoost(boost2lab);
reactants[2].ApplyBoost(boost2lab);
reactants[3] = reactants[0] - reactants[2];
m_residual->vec4 = m_target->vec4 - m_ejectile->vec4;
//energy loss for the *light* break up nucleus
double ejectKE = reactants[2].GetKE() - target->GetEjectileEnergyLoss(reactants[2].GetZ(), reactants[2].GetA(), reactants[2].GetKE(), rxnLayer, reactants[2].GetTheta());
double ejectP = std::sqrt(ejectKE*(ejectKE + 2.0*reactants[2].GetGroundStateMass()));
double ejectE = ejectKE + reactants[2].GetGroundStateMass();
reactants[2].SetVectorSpherical(reactants[2].GetTheta(), reactants[2].GetPhi(), ejectP, ejectE);
double ejectKE = m_ejectile->GetKE() -
m_layeredTarget->GetEjectileEnergyLoss(m_ejectile->Z, m_ejectile->A, m_ejectile->GetKE(), m_rxnLayer, m_ejectile->vec4.Theta());
double ejectP = std::sqrt(ejectKE*(ejectKE + 2.0*m_ejectile->groundStateMass));
double ejectE = ejectKE + m_ejectile->groundStateMass;
m_ejectile->SetVec4Spherical(m_ejectile->vec4.Theta(), m_ejectile->vec4.Phi(), ejectP, ejectE);
//if on, get eloss for *heavy* break up nucleus
if(resid_elossFlag) {
double residKE = reactants[3].GetKE() - target->GetEjectileEnergyLoss(reactants[3].GetZ(), reactants[3].GetA(), reactants[3].GetKE(), rxnLayer, reactants[3].GetTheta());
double residP = std::sqrt(residKE*(residKE + 2.0*reactants[3].GetInvMass()));
double residE = residKE + reactants[3].GetInvMass();
reactants[3].SetVectorSpherical(reactants[3].GetTheta(), reactants[3].GetPhi(), residP, residE);
if(m_isResidEloss)
{
double residKE = m_residual->GetKE() -
m_layeredTarget->GetEjectileEnergyLoss(m_residual->Z, m_residual->A, m_residual->GetKE(), m_rxnLayer, m_residual->vec4.Theta());
double residP = std::sqrt(residKE*(residKE + 2.0*m_residual->vec4.M()));
double residE = residKE + m_residual->vec4.M();
m_residual->SetVec4Spherical(m_residual->vec4.Theta(), m_residual->vec4.Phi(), residP, residE);
}
}

74
src/Mask/Reaction.h Normal file
View File

@ -0,0 +1,74 @@
/*
Reaction.h
Reaction is a class which implements either a decay or scattering reaction. As such it requires either
3 (decay) or 4 (scattering) nuclei to perform any calcualtions. I also links together the target, which provides
energy loss calculations, with the kinematics. Note that Reaction does not own the LayeredTarget.
--GWM Jan. 2021
*/
#ifndef REACTION_H
#define REACTION_H
#include "Nucleus.h"
#include "LayeredTarget.h"
namespace Mask {
class Reaction
{
public:
Reaction();
Reaction(Nucleus* target, Nucleus* projectile, Nucleus* ejectile, Nucleus* residual);
~Reaction();
bool Calculate(); //do sim
void BindNuclei(Nucleus* target, Nucleus* projectile, Nucleus* ejectile, Nucleus* residual);
void SetBeamKE(double bke);
void SetEjectileThetaType(int type);
void SetLayeredTarget(LayeredTarget* targ) { m_layeredTarget = targ; };
void SetPolarRxnAngle(double theta) { m_theta = theta; };
void SetAzimRxnAngle(double phi) { m_phi = phi; };
void SetExcitation(double ex) { m_ex = ex; };
void BindTarget(Nucleus* nuc) { m_target = nuc; };
void BindProjectile(Nucleus* nuc) { m_projectile = nuc; };
void BindEjectile(Nucleus* nuc) { m_ejectile = nuc; };
void BindResidual(Nucleus* nuc) { m_residual = nuc; };
void SetRxnLayer(std::size_t layer) { m_rxnLayer = layer; };
void SetResidualEnergyLoss(bool isEloss) { m_isResidEloss = isEloss; };
bool IsDecay() const { return m_isDecay; };
std::size_t GetRxnLayer() const { return m_rxnLayer; };
private:
void CalculateDecay(); //target -> light_decay (eject) + heavy_decay(resid)
void CalculateReaction(); //target + project -> eject + resid
void CalculateReactionThetaLab();
void CalculateReactionThetaCM();
//Reactants -> NOT OWNED BY RXN
Nucleus* m_target;
Nucleus* m_projectile;
Nucleus* m_ejectile;
Nucleus* m_residual;
LayeredTarget* m_layeredTarget; //not owned by Reaction
double m_bke, m_theta, m_phi, m_ex;
int m_rxnLayer;
int m_ejectThetaType;
bool m_isDecay, m_isInit, m_isResidEloss;
static constexpr int s_lab = 0;
static constexpr int s_centerOfMass = 1;
};
}
#endif

12
src/Mask/ReactionSystem.cpp
View File

@ -5,20 +5,22 @@
namespace Mask {
ReactionSystem::ReactionSystem() :
m_beamDist(nullptr), m_theta1Range(nullptr), m_phi1Range(nullptr), m_exDist(nullptr), target_set_flag(false),
rxnLayer(0), m_sys_equation("")
m_beamDist(nullptr), m_theta1Range(nullptr), m_phi1Range(nullptr), m_exDist(nullptr), m_isTargetSet(false),
m_rxnLayer(0), m_sysEquation("")
{
}
ReactionSystem::~ReactionSystem() {
ReactionSystem::~ReactionSystem()
{
delete m_beamDist;
delete m_theta1Range;
delete m_phi1Range;
delete m_exDist;
}
void ReactionSystem::AddTargetLayer(std::vector<int>& zt, std::vector<int>& at, std::vector<int>& stoich, double thickness) {
target.AddLayer(zt, at, stoich, thickness);
void ReactionSystem::AddTargetLayer(const std::vector<int>& zt, const std::vector<int>& at, const std::vector<int>& stoich, double thickness)
{
m_target.AddLayer(zt, at, stoich, thickness);
}
}

39
include/ReactionSystem.h → src/Mask/ReactionSystem.h
View File

@ -16,59 +16,64 @@
namespace Mask {
class ReactionSystem {
class ReactionSystem
{
public:
ReactionSystem();
virtual ~ReactionSystem();
virtual bool SetNuclei(std::vector<int>& z, std::vector<int>& a) = 0;
virtual bool SetNuclei(const std::vector<int>& z, const std::vector<int>& a) = 0;
virtual void RunSystem() = 0;
virtual const std::vector<Nucleus>& GetNuclei() = 0;
virtual std::vector<Nucleus>* GetNuclei() = 0;
void AddTargetLayer(std::vector<int>& zt, std::vector<int>& at, std::vector<int>& stoich, double thickness);
void AddTargetLayer(const std::vector<int>& zt, const std::vector<int>& at, const std::vector<int>& stoich, double thickness);
/*Set sampling parameters*/
inline void SetBeamDistro(double mean, double sigma) {
void SetBeamDistro(double mean, double sigma)
{
if(m_beamDist)
delete m_beamDist;
m_beamDist = new std::normal_distribution<double>(mean, sigma);
}
inline void SetTheta1Range(double min, double max) {
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));
m_theta1Range = new std::uniform_real_distribution<double>(std::cos(min*s_deg2rad), std::cos(max*s_deg2rad));
}
inline void SetPhi1Range(double min, double max) {
void SetPhi1Range(double min, double max)
{
if(m_phi1Range)
delete m_phi1Range;
m_phi1Range = new std::uniform_real_distribution<double>(min*deg2rad, max*deg2rad);
m_phi1Range = new std::uniform_real_distribution<double>(min*s_deg2rad, max*s_deg2rad);
}
inline void SetExcitationDistro(double mean, double sigma) {
void SetExcitationDistro(double mean, double sigma)
{
if(m_exDist)
delete m_exDist;
m_exDist = new std::normal_distribution<double>(mean, sigma);
}
inline const std::string& GetSystemEquation() const { return m_sys_equation; }
const std::string& GetSystemEquation() const { return m_sysEquation; }
protected:
virtual void LinkTarget() = 0;
virtual void SetSystemEquation() = 0;
LayeredTarget target;
LayeredTarget m_target;
//Sampling information
std::normal_distribution<double> *m_beamDist, *m_exDist;
std::uniform_real_distribution<double> *m_theta1Range, *m_phi1Range;
bool target_set_flag, gen_set_flag;
int rxnLayer;
std::string m_sys_equation;
std::vector<Nucleus> nuclei;
static constexpr double deg2rad = M_PI/180.0;
bool m_isTargetSet;
std::size_t m_rxnLayer;
std::string m_sysEquation;
std::vector<Nucleus> m_nuclei;
static constexpr double s_deg2rad = M_PI/180.0;
};
}

0
include/RxnType.h → src/Mask/RxnType.h
View File

41
src/Mask/Stopwatch.cpp
View File

@ -7,25 +7,34 @@
*/
#include "Stopwatch.h"
Stopwatch::Stopwatch() {
start_time = Clock::now();
stop_time = start_time;
}
namespace Mask {
Stopwatch::~Stopwatch() {}
Stopwatch::Stopwatch()
{
start_time = Clock::now();
stop_time = start_time;
}
void Stopwatch::Start() {
start_time = Clock::now();
}
Stopwatch::~Stopwatch() {}
void Stopwatch::Stop() {
stop_time = Clock::now();
}
void Stopwatch::Start()
{
start_time = Clock::now();
}
double Stopwatch::GetElapsedSeconds() {
return std::chrono::duration_cast<std::chrono::duration<double>>(stop_time-start_time).count();
}
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;
}
double Stopwatch::GetElapsedMilliseconds() {
return std::chrono::duration_cast<std::chrono::duration<double>>(stop_time-start_time).count()*1000.0;
}

32
src/Mask/Stopwatch.h Normal file
View File

@ -0,0 +1,32 @@
/*
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>
namespace Mask {
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

57
src/Mask/Target.cpp
View File

@ -16,45 +16,50 @@ Written by G.W. McCann Aug. 2020
namespace Mask {
/*Targets must be of known thickness*/
Target::Target(double thick) {
thickness = thick;
thickness_gcm2 = thickness*1.0e-6;
Target::Target(const std::vector<int>& z, const std::vector<int>& a, const std::vector<int>& stoich, double thick) :
m_thickness(thick), m_thickness_gcm2(thick*1.0e-6)
{
Init(z, a, stoich);
}
Target::~Target() {}
/*Set target elements of given Z, A, S*/
void Target::SetElements(std::vector<int>& z, std::vector<int>& a, std::vector<int>& stoich) {
Z = z;
A = a;
Stoich = stoich;
void Target::Init(const std::vector<int>& z, const std::vector<int>& a, const std::vector<int>& stoich)
{
m_Z = z;
m_A = a;
m_stoich = stoich;
MassLookup& masses = MassLookup::GetInstance();
eloss.SetTargetComponents(Z, A, Stoich);
for(size_t i=0; i<Z.size(); i++)
for(size_t i=0; i<m_Z.size(); i++)
{
target_material.add_element(masses.FindMassU(Z[i], A[i]), Z[i], Stoich[i]);
m_material.add_element(masses.FindMassU(m_Z[i], m_A[i]), m_Z[i], m_stoich[i]);
}
}
/*Element verification*/
bool Target::ContainsElement(int z, int a) {
for(unsigned int i=0; i<Z.size(); i++)
if( z == Z[i] && a == A[i])
bool Target::ContainsElement(int z, int a)
{
for(std::size_t i=0; i<m_Z.size(); i++)
{
if( z == m_Z[i] && a == m_A[i])
return true;
}
return false;
}
/*Calculates energy loss for travelling all the way through the target*/
double Target::GetEnergyLossTotal(int zp, int ap, double startEnergy, double theta) {
double Target::GetEnergyLossTotal(int zp, int ap, double startEnergy, double theta)
{
if(theta == M_PI/2.)
return startEnergy;
else if (theta > M_PI/2.)
theta = M_PI - theta;
catima::Projectile proj(MassLookup::GetInstance().FindMassU(zp, ap), zp, 0.0, 0.0);
proj.T = startEnergy/proj.A;
target_material.thickness(thickness_gcm2/fabs(cos(theta)));
return catima::integrate_energyloss(proj, target_material);
//return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/fabs(cos(theta)));
m_material.thickness(m_thickness_gcm2/fabs(cos(theta)));
return catima::integrate_energyloss(proj, m_material);
}
/*Calculates the energy loss for traveling some fraction through the target*/
@ -67,13 +72,13 @@ namespace Mask {
catima::Projectile proj(MassLookup::GetInstance().FindMassU(zp, ap), zp, 0.0, 0.0);
proj.T = finalEnergy/proj.A;
target_material.thickness(thickness_gcm2*percent_depth/fabs(cos(theta)));
return catima::integrate_energyloss(proj, target_material);
//return eloss.GetEnergyLoss(zp, ap, finalEnergy, thickness*percent_depth/(std::fabs(std::cos(theta))));
m_material.thickness(m_thickness_gcm2*percent_depth/fabs(cos(theta)));
return catima::integrate_energyloss(proj, m_material);
}
/*Calculates reverse energy loss for travelling all the way through the target*/
double Target::GetReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double theta) {
double Target::GetReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double theta)
{
if(theta == M_PI/2.)
return finalEnergy;
else if (theta > M_PI/2.)
@ -81,9 +86,8 @@ namespace Mask {
catima::Projectile proj(MassLookup::GetInstance().FindMassU(zp, ap), zp, 0.0, 0.0);
proj.T = finalEnergy/proj.A;
target_material.thickness(thickness_gcm2/fabs(cos(theta)));
return catima::reverse_integrate_energyloss(proj, target_material);
//return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/fabs(cos(theta)));
m_material.thickness(m_thickness_gcm2/fabs(cos(theta)));
return catima::reverse_integrate_energyloss(proj, m_material);
}
/*Calculates the reverse energy loss for traveling some fraction through the target*/
@ -95,9 +99,8 @@ namespace Mask {
theta = M_PI-theta;
catima::Projectile proj(MassLookup::GetInstance().FindMassU(zp, ap), zp, 0.0, 0.0);
proj.T = finalEnergy/proj.A;
target_material.thickness(thickness_gcm2*percent_depth/fabs(cos(theta)));
return catima::reverse_integrate_energyloss(proj, target_material);
//return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness*percent_depth/(std::fabs(std::cos(theta))));
m_material.thickness(m_thickness_gcm2*percent_depth/fabs(cos(theta)));
return catima::reverse_integrate_energyloss(proj, m_material);
}
}

25
include/Target.h → src/Mask/Target.h
View File

@ -16,7 +16,6 @@ Written by G.W. McCann Aug. 2020
#include <string>
#include <vector>
#include <cmath>
#include "EnergyLoss.h"
#include "catima/gwm_integrators.h"
#include "MassLookup.h"
@ -25,26 +24,26 @@ namespace Mask {
class Target {
public:
Target(double thick);
Target(const std::vector<int>& z, const std::vector<int>& a, const std::vector<int>& stoich, double thick);
~Target();
void SetElements(std::vector<int>& z, std::vector<int>& a, std::vector<int>& stoich);
bool ContainsElement(int z, int a);
double GetEnergyLossTotal(int zp, int ap, double startEnergy, double angle);
double GetReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double angle);
double GetEnergyLossFractionalDepth(int zp, int ap, double startEnergy, double angle, double percent_depth);
double GetReverseEnergyLossFractionalDepth(int zp, int ap, double finalEnergy, double angle, double percent_depth);
inline const double& GetThickness() { return thickness; }
inline int GetNumberOfElements() { return Z.size(); }
inline int GetElementZ(int index) { return Z[index]; }
inline int GetElementA(int index) { return A[index]; }
inline int GetElementStoich(int index) { return Stoich[index]; }
inline const double& GetThickness() { return m_thickness; }
inline int GetNumberOfElements() { return m_Z.size(); }
inline int GetElementZ(int index) { return m_Z[index]; }
inline int GetElementA(int index) { return m_A[index]; }
inline int GetElementStoich(int index) { return m_stoich[index]; }
private:
EnergyLoss eloss;
catima::Material target_material;
double thickness;
double thickness_gcm2;
std::vector<int> Z, A, Stoich;
void Init(const std::vector<int>& z, const std::vector<int>& a, const std::vector<int>& stoich);
catima::Material m_material;
double m_thickness;
double m_thickness_gcm2;
std::vector<int> m_Z, m_A, m_stoich;
};

150
src/Mask/ThreeStepSystem.cpp
View File

@ -7,133 +7,131 @@ namespace Mask {
ThreeStepSystem::ThreeStepSystem() :
ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
{
nuclei.resize(8);
m_nuclei.resize(8);
}
ThreeStepSystem::ThreeStepSystem(std::vector<int>& z, std::vector<int>& a) :
ThreeStepSystem::ThreeStepSystem(const std::vector<int>& z, const std::vector<int>& a) :
ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
{
nuclei.resize(8);
m_nuclei.resize(8);
SetNuclei(z, a);
}
ThreeStepSystem::~ThreeStepSystem() {}
bool ThreeStepSystem::SetNuclei(std::vector<int>&z, std::vector<int>& a) {
if(z.size() != a.size() || z.size() < 5) {
bool ThreeStepSystem::SetNuclei(const std::vector<int>& z, const 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 zb4 = zb2 - z[4];
int ar = a[0] + a[1] - a[2];
int ab2 = ar - a[3];
int ab4 = ab2 - a[4];
m_nuclei[0] = CreateNucleus(z[0], a[0]); //target
m_nuclei[1] = CreateNucleus(z[1], a[1]); //projectile
m_nuclei[2] = CreateNucleus(z[2], a[2]); //ejectile
m_nuclei[3] = CreateNucleus(zr, ar); //residual
m_nuclei[4] = CreateNucleus(z[3], a[3]); //breakup1
m_nuclei[5] = CreateNucleus(zb2, ab2); //breakup2
m_nuclei[5] = CreateNucleus(z[4], a[4]); //breakup3
m_nuclei[5] = CreateNucleus(zb4, ab4); //breakup4
m_step1.BindNuclei(&(m_nuclei[0]), &(m_nuclei[1]), &(m_nuclei[2]), &(m_nuclei[3]));
m_step2.BindNuclei(&(m_nuclei[3]), nullptr, &(m_nuclei[4]), &(m_nuclei[5]));
m_step3.BindNuclei(&(m_nuclei[5]), nullptr, &(m_nuclei[6]), &(m_nuclei[7]));
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;
}
const std::vector<Nucleus>& ThreeStepSystem::GetNuclei()
std::vector<Nucleus>* ThreeStepSystem::GetNuclei()
{
nuclei[0] = step1.GetTarget();
nuclei[1] = step1.GetProjectile();
nuclei[2] = step1.GetEjectile();
nuclei[3] = step1.GetResidual();
nuclei[4] = step2.GetEjectile();
nuclei[5] = step2.GetResidual();
nuclei[6] = step3.GetEjectile();
nuclei[7] = step3.GetResidual();
return nuclei;
return &m_nuclei;
}
void ThreeStepSystem::LinkTarget() {
step1.SetLayeredTarget(&target);
step2.SetLayeredTarget(&target);
step3.SetLayeredTarget(&target);
void ThreeStepSystem::LinkTarget()
{
m_step1.SetLayeredTarget(&m_target);
m_step2.SetLayeredTarget(&m_target);
m_step3.SetLayeredTarget(&m_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 {
m_rxnLayer = m_target.FindLayerContaining(m_nuclei[0].Z, m_nuclei[0].A);
if(m_rxnLayer != m_target.GetNumberOfLayers())
{
m_step1.SetRxnLayer(m_rxnLayer);
m_step2.SetRxnLayer(m_rxnLayer);
m_step3.SetRxnLayer(m_rxnLayer);
m_isTargetSet = true;
} else
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::SetSystemEquation()
{
std::stringstream stream;
stream << m_nuclei[0].isotopicSymbol << "("
<< m_nuclei[1].isotopicSymbol << ", "
<< m_nuclei[2].isotopicSymbol << ")"
<< m_nuclei[3].isotopicSymbol << "->"
<< m_nuclei[4].isotopicSymbol << "+"
<< m_nuclei[5].isotopicSymbol << "->"
<< m_nuclei[6].isotopicSymbol << "+"
<< m_nuclei[7].isotopicSymbol;
m_sysEquation = stream.str();
}
void ThreeStepSystem::RunSystem() {
//Link up the target if it hasn't been done yet
if(!target_set_flag) {
if(!m_isTargetSet)
LinkTarget();
}
//Sample parameters
double bke = (*m_beamDist)(RandomGenerator::GetInstance().GetGenerator());
double rxnTheta = acos((*m_theta1Range)(RandomGenerator::GetInstance().GetGenerator()));
double rxnPhi = (*m_phi1Range)(RandomGenerator::GetInstance().GetGenerator());
double decay1costheta = decay1dist.GetRandomCosTheta();
double decay1costheta = m_step2Distribution.GetRandomCosTheta();
double decay1Theta = std::acos(decay1costheta);
double decay1Phi = m_phi2Range(RandomGenerator::GetInstance().GetGenerator());
double decay2costheta = decay2dist.GetRandomCosTheta();
double decay2costheta = m_step3Distribution.GetRandomCosTheta();
double decay2Theta = std::acos(decay2costheta);
double decay2Phi = m_phi2Range(RandomGenerator::GetInstance().GetGenerator());
double residEx = (*m_exDist)(RandomGenerator::GetInstance().GetGenerator());
step1.SetBeamKE(bke);
step1.SetPolarRxnAngle(rxnTheta);
step1.SetAzimRxnAngle(rxnPhi);
step1.SetExcitation(residEx);
m_step1.SetBeamKE(bke);
m_step1.SetPolarRxnAngle(rxnTheta);
m_step1.SetAzimRxnAngle(rxnPhi);
m_step1.SetExcitation(residEx);
step2.SetPolarRxnAngle(decay1Theta);
step2.SetAzimRxnAngle(decay1Phi);
m_step2.SetPolarRxnAngle(decay1Theta);
m_step2.SetAzimRxnAngle(decay1Phi);
step3.SetPolarRxnAngle(decay2Theta);
step3.SetAzimRxnAngle(decay2Phi);
m_step3.SetPolarRxnAngle(decay2Theta);
m_step3.SetAzimRxnAngle(decay2Phi);
step1.Calculate();
m_step1.Calculate();
step2.SetTarget(step1.GetResidual());
if(decay1costheta == -10) {
step2.ResetEjectile();
step2.ResetResidual();
step3.ResetTarget();
step3.ResetEjectile();
step3.ResetResidual();
if(decay1costheta == -10)
{
m_nuclei[4].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[4].groundStateMass);
m_nuclei[5].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[5].groundStateMass);
m_nuclei[6].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[6].groundStateMass);
m_nuclei[7].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[7].groundStateMass);
return;
}
step2.Calculate();
m_step2.Calculate();
step3.SetTarget(step2.GetResidual());
if(decay2costheta == -10) {
step3.ResetEjectile();
step3.ResetResidual();
if(decay2costheta == -10)
{
m_nuclei[6].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[6].groundStateMass);
m_nuclei[7].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[7].groundStateMass);
return;
}
step3.TurnOnResidualEloss();
step3.Calculate();
m_step3.SetResidualEnergyLoss(true);
m_step3.Calculate();
}

42
src/Mask/ThreeStepSystem.h Normal file
View File

@ -0,0 +1,42 @@
#ifndef THREESTEPSYSTEM_H
#define THREESTEPSYSTEM_H
#include "ReactionSystem.h"
#include "AngularDistribution.h"
namespace Mask {
class ThreeStepSystem : public ReactionSystem {
public:
ThreeStepSystem();
ThreeStepSystem(const std::vector<int>& z, const std::vector<int>& a);
~ThreeStepSystem();
bool SetNuclei(const std::vector<int>& z, const std::vector<int>& a) override;
void RunSystem() override;
std::vector<Nucleus>* GetNuclei() override;
inline void SetDecay1Distribution(const std::string& filename) { m_step2Distribution.ReadDistributionFile(filename); };
inline void SetDecay2Distribution(const std::string& filename) { m_step3Distribution.ReadDistributionFile(filename); };
void SetReactionThetaType(int type) { m_step1.SetEjectileThetaType(type); };
int GetDecay1AngularMomentum() { return m_step2Distribution.GetL(); };
int GetDecay2AngularMomentum(){ return m_step3Distribution.GetL(); };
double GetDecay1BranchingRatio() { return m_step2Distribution.GetBranchingRatio(); };
double GetDecay2BranchingRatio(){ return m_step3Distribution.GetBranchingRatio(); };
protected:
void LinkTarget() override;
void SetSystemEquation() override;
std::uniform_real_distribution<double> m_phi2Range;
Reaction m_step1, m_step2, m_step3;
AngularDistribution m_step2Distribution, m_step3Distribution;
};
}
#endif

121
src/Mask/TwoStepSystem.cpp
View File

@ -2,111 +2,114 @@
#include "RandomGenerator.h"
#include "KinematicsExceptions.h"
#include <sstream>
namespace Mask {
TwoStepSystem::TwoStepSystem() :
ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
{
nuclei.resize(6);
m_nuclei.resize(6);
}
TwoStepSystem::TwoStepSystem(std::vector<int>& z, std::vector<int>& a) :
TwoStepSystem::TwoStepSystem(const std::vector<int>& z, const std::vector<int>& a) :
ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
{
nuclei.resize(6);
m_nuclei.resize(6);
SetNuclei(z, a);
}
TwoStepSystem::~TwoStepSystem() {
TwoStepSystem::~TwoStepSystem() {}
}
bool TwoStepSystem::SetNuclei(std::vector<int>&z, std::vector<int>& a) {
if(z.size() != a.size() || z.size() != 4) {
bool TwoStepSystem::SetNuclei(const std::vector<int>&z, const 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];
int zb = zr - z[3];
int ab = ar - a[3];
m_nuclei[0] = CreateNucleus(z[0], a[0]); //target
m_nuclei[1] = CreateNucleus(z[1], a[1]); //projectile
m_nuclei[2] = CreateNucleus(z[2], a[2]); //ejectile
m_nuclei[3] = CreateNucleus(zr, ar); //residual
m_nuclei[4] = CreateNucleus(z[3], a[3]); //breakup1
m_nuclei[5] = CreateNucleus(zb, ab); //breakup2
step1.SetNuclei(z[0], a[0], z[1], a[1], z[2], a[2]);
step2.SetNuclei(zr, ar, 0, 0, z[3], a[3]);
m_step1.BindNuclei(&(m_nuclei[0]), &(m_nuclei[1]), &(m_nuclei[2]), &(m_nuclei[3]));
m_step2.BindNuclei(&(m_nuclei[3]), nullptr, &(m_nuclei[4]), &(m_nuclei[5]));
SetSystemEquation();
return true;
}
const std::vector<Nucleus>& TwoStepSystem::GetNuclei()
std::vector<Nucleus>* TwoStepSystem::GetNuclei()
{
nuclei[0] = step1.GetTarget();
nuclei[1] = step1.GetProjectile();
nuclei[2] = step1.GetEjectile();
nuclei[3] = step1.GetResidual();
nuclei[4] = step2.GetEjectile();
nuclei[5] = step2.GetResidual();
return nuclei;
return &m_nuclei;
}
void TwoStepSystem::LinkTarget() {
step1.SetLayeredTarget(&target);
step2.SetLayeredTarget(&target);
void TwoStepSystem::LinkTarget()
{
m_step1.SetLayeredTarget(&m_target);
m_step2.SetLayeredTarget(&m_target);
rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA());
if(rxnLayer != -1) {
step1.SetRxnLayer(rxnLayer);
step2.SetRxnLayer(rxnLayer);
target_set_flag = true;
} else {
m_rxnLayer = m_target.FindLayerContaining(m_nuclei[0].Z, m_nuclei[0].A);
if(m_rxnLayer != m_target.GetNumberOfLayers())
{
m_step1.SetRxnLayer(m_rxnLayer);
m_step2.SetRxnLayer(m_rxnLayer);
m_isTargetSet = 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::SetSystemEquation()
{
std::stringstream stream;
stream << m_nuclei[0].isotopicSymbol << "("
<< m_nuclei[1].isotopicSymbol << ", "
<< m_nuclei[2].isotopicSymbol << ")"
<< m_nuclei[3].isotopicSymbol << "->"
<< m_nuclei[4].isotopicSymbol << "+"
<< m_nuclei[5].isotopicSymbol;
m_sysEquation = stream.str();
}
void TwoStepSystem::RunSystem() {
void TwoStepSystem::RunSystem()
{
//Link up the target if it hasn't been done yet
if(!target_set_flag) {
if(!m_isTargetSet)
LinkTarget();
}
//Sample parameters
double bke = (*m_beamDist)(RandomGenerator::GetInstance().GetGenerator());
double rxnTheta = acos((*m_theta1Range)(RandomGenerator::GetInstance().GetGenerator()));
double rxnPhi = (*m_phi1Range)(RandomGenerator::GetInstance().GetGenerator());
double decay1costheta = decay1dist.GetRandomCosTheta();
double decay1costheta = m_step2Distribution.GetRandomCosTheta();
double decay1Theta = std::acos(decay1costheta);
double decay1Phi = m_phi2Range(RandomGenerator::GetInstance().GetGenerator());
double residEx = (*m_exDist)(RandomGenerator::GetInstance().GetGenerator());
step1.SetBeamKE(bke);
step1.SetPolarRxnAngle(rxnTheta);
step1.SetAzimRxnAngle(rxnPhi);
step1.SetExcitation(residEx);
m_step1.SetBeamKE(bke);
m_step1.SetPolarRxnAngle(rxnTheta);
m_step1.SetAzimRxnAngle(rxnPhi);
m_step1.SetExcitation(residEx);
step2.SetPolarRxnAngle(decay1Theta);
step2.SetAzimRxnAngle(decay1Phi);
m_step2.SetPolarRxnAngle(decay1Theta);
m_step2.SetAzimRxnAngle(decay1Phi);
step1.Calculate();
m_step1.Calculate();
step2.SetTarget(step1.GetResidual());
if(decay1costheta == -10) {
step2.ResetEjectile();
step2.ResetResidual();
if(decay1costheta == -10)
{
m_nuclei[4].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[4].groundStateMass);
m_nuclei[5].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[5].groundStateMass);
return;
}
step2.TurnOnResidualEloss();
step2.Calculate();
m_step2.SetResidualEnergyLoss(true);
m_step2.Calculate();
}

40
src/Mask/TwoStepSystem.h Normal file
View File

@ -0,0 +1,40 @@
#ifndef TWOSTEPSYSTEM_H
#define TWOSTEPSYSTEM_H
#include "ReactionSystem.h"
#include "AngularDistribution.h"
namespace Mask {
class TwoStepSystem : public ReactionSystem
{
public:
TwoStepSystem();
TwoStepSystem(const std::vector<int>& z, const std::vector<int>& a);
~TwoStepSystem();
bool SetNuclei(const std::vector<int>& z, const std::vector<int>& a) override;
void RunSystem() override;
std::vector<Nucleus>* GetNuclei() override;
void SetDecay1Distribution(const std::string& filename) { m_step2Distribution.ReadDistributionFile(filename); };
void SetReactionThetaType(int type) { m_step1.SetEjectileThetaType(type); };
int GetDecay1AngularMomentum() { return m_step2Distribution.GetL(); };
double GetDecay1BranchingRatio() { return m_step2Distribution.GetBranchingRatio(); };
private:
void LinkTarget() override;
void SetSystemEquation() override;
std::uniform_real_distribution<double> m_phi2Range;
Reaction m_step1, m_step2;
AngularDistribution m_step2Distribution;
};
}
#endif

2
src/MaskApp/CMakeLists.txt
View File

@ -1,5 +1,5 @@
add_executable(MaskApp)
target_include_directories(MaskApp PUBLIC ${MASK_INCLUDE_DIR})
target_include_directories(MaskApp PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/..)
target_sources(MaskApp PUBLIC
main.cpp

25
src/MaskApp/main.cpp
View File

@ -1,27 +1,32 @@
#include <iostream>
#include <string>
#include "Stopwatch.h"
#include "MaskFile.h"
#include "MaskApp.h"
#include "KinematicsExceptions.h"
#include "Mask/Stopwatch.h"
#include "Mask/MaskApp.h"
#include "Mask/KinematicsExceptions.h"
int main(int argc, char** argv) {
if(argc<2) {
int main(int argc, char** argv)
{
if(argc<2)
{
std::cerr<<"Incorrect number of arguments!"<<std::endl;
return 1;
}
Stopwatch sw;
Mask::Stopwatch sw;
Mask::MaskApp calculator;
sw.Start();
try {
if(!calculator.LoadConfig(argv[1])) {
try
{
if(!calculator.LoadConfig(argv[1]))
{
std::cerr<<"Unable to read input file!"<<std::endl;
return 1;
}
calculator.Run();
} catch(const std::exception& e) {
}
catch(const std::exception& e)
{
std::cerr<<"Exception caught! Information: "<<e.what()<<std::endl;
std::cerr<<"Terminating process."<<std::endl;
return 1;

7
src/Plotters/ROOT/CMakeLists.txt → src/Plotters/CMakeLists.txt
View File

@ -1,15 +1,18 @@
find_package(ROOT REQUIRED)
add_executable(RootPlot)
target_include_directories(RootPlot
PUBLIC ${MASK_INCLUDE_DIR}
PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
${CMAKE_CURRENT_SOURCE_DIR}/..
SYSTEM PUBLIC ${ROOT_INCLUDE_DIRS}
)
target_sources(RootPlot PUBLIC
RootPlotter.cpp
RootPlotter.h
main.cpp
)
target_link_libraries(RootPlot
MaskDict
Mask
${ROOT_LIBRARIES}
)

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

@ -1,323 +0,0 @@
#include "RootPlotter.h"
#include <TFile.h>
#include <TVector3.h>
#include <iostream>
RootPlotter::RootPlotter() :
table(new THashTable())
{
}
RootPlotter::~RootPlotter() {}
void RootPlotter::FillData(const Mask::Nucleus& nuc, double detKE, const std::string& modifier) {
std::string sym = nuc.GetIsotopicSymbol();
std::string ke_vs_th_name = sym + modifier + "_ke_vs_theta";
std::string ke_vs_th_title = ke_vs_th_name + ";#theta_{lab} (degrees);Kinetic Energy (MeV)";
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 th_vs_ph_name = sym + modifier + "_theta_vs_phi";
std::string th_vs_ph_title = th_vs_ph_name + ";#theta_{lab};#phi_{lab}";
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";
if(detKE == 0.0)
{
MyFill(ke_vs_th_name.c_str(), ke_vs_th_title.c_str(), nuc.GetTheta()*rad2deg, nuc.GetKE(), 2);
MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), nuc.GetPhi()*rad2deg, nuc.GetKE(), 4);
MyFill(th_vs_ph_name.c_str(), th_vs_ph_title.c_str(), nuc.GetTheta()*rad2deg, nuc.GetPhi()*rad2deg, 2);
MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy());
MyFill(angdist_name.c_str(), angdist_title.c_str(),20,-1.0,1.0,std::cos(nuc.GetThetaCM()));
}
else
{
MyFill(ke_vs_th_name.c_str(), ke_vs_th_title.c_str(), nuc.GetTheta()*rad2deg, detKE, 2);
MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), nuc.GetPhi()*rad2deg, detKE, 4);
MyFill(th_vs_ph_name.c_str(), th_vs_ph_title.c_str(), nuc.GetTheta()*rad2deg, nuc.GetPhi()*rad2deg, 2);
MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy());
MyFill(angdist_name.c_str(), angdist_title.c_str(),20,-1.0,1.0,std::cos(nuc.GetThetaCM()));
}
}
void RootPlotter::FillCorrelations(const Mask::MaskFileData& data, Mask::RxnType type)
{
std::string theta_eject_theta_resid_name = "theta_eject_theta_resid_cor";
std::string theta_eject_theta_resid_title = theta_eject_theta_resid_name + ";#theta_{lab} Ejectile (deg);#theta_{lab} Residual";
if(type == Mask::RxnType::PureDecay)
{
MyFill(theta_eject_theta_resid_name, theta_eject_theta_resid_title, data.theta[1]*rad2deg, data.theta[2]*rad2deg, 4);
}
else
{
MyFill(theta_eject_theta_resid_name, theta_eject_theta_resid_title, data.theta[2]*rad2deg, data.theta[3]*rad2deg, 4);
}
if(type == Mask::RxnType::TwoStepRxn || type == Mask::RxnType::ThreeStepRxn)
{
TVector3 p1, p2;
p1.SetMagThetaPhi(1.0, data.theta[3], data.phi[3]);
p2.SetMagThetaPhi(1.0, data.theta[4], data.phi[4]);
double theta_resid_break1 = std::acos(p1.Dot(p2));
std::string theta_break1_theta_break2_name = "theta_break1_theta_break2_cor";
std::string theta_break1_theta_break2_title = theta_break1_theta_break2_name + ";#theta_{lab} Breakup1 (deg);#theta_{lab} Breakup2 (deg)";
MyFill(theta_break1_theta_break2_name, theta_break1_theta_break2_title, data.theta[4]*rad2deg, data.theta[5]*rad2deg, 4);
std::string theta_resid_theta_break1_name = "theta_resid_theta_break1_cor";
std::string theta_resid_theta_break1_title = theta_resid_theta_break1_name + ";#theta_{lab} Residual (deg);#theta_{lab} Breakup1 (deg)";
MyFill(theta_resid_theta_break1_name, theta_resid_theta_break1_title, data.theta[3]*rad2deg, data.theta[4]*rad2deg, 4);
std::string ke_break1_theta_rel_name = "ke_break1_theta_rel";
std::string ke_break1_theta_rel_title = ke_break1_theta_rel_name + ";#theta_{resid-break1};KE_{break1} (MeV)";
MyFill(ke_break1_theta_rel_name, ke_break1_theta_rel_title, theta_resid_break1*rad2deg, data.KE[4], 4);
}
if(type == Mask::RxnType::ThreeStepRxn)
{
std::string theta_break3_theta_break4_name = "theta_break3_theta_break4_cor";
std::string theta_break3_theta_break4_title = theta_break3_theta_break4_name + ";#theta_{lab} Breakup3 (deg);#theta_{lab} Breakup4 (deg)";
MyFill(theta_break3_theta_break4_name, theta_break3_theta_break4_title, data.theta[6]*rad2deg, data.theta[7]*rad2deg, 4);
}
}
void RootPlotter::FillCorrelationsDetected(const Mask::MaskFileData& data, Mask::RxnType type)
{
std::string theta_eject_theta_resid_name = "theta_eject_theta_resid_cor_detected";
std::string theta_eject_theta_resid_title = theta_eject_theta_resid_name + ";#theta_{lab} Ejectile (deg);#theta_{lab} Residual";
if(type == Mask::RxnType::PureDecay && data.detect_flag[1] && data.detect_flag[2])
{
MyFill(theta_eject_theta_resid_name, theta_eject_theta_resid_title, data.theta[1]*rad2deg, data.theta[2]*rad2deg, 4);
}
else if(data.detect_flag[2] && data.detect_flag[3])
{
MyFill(theta_eject_theta_resid_name, theta_eject_theta_resid_title, data.theta[2]*rad2deg, data.theta[3]*rad2deg, 4);
}
if((type == Mask::RxnType::TwoStepRxn || type == Mask::RxnType::ThreeStepRxn) && data.detect_flag[4])
{
TVector3 p1, p2;
p1.SetMagThetaPhi(1.0, data.theta[3], data.phi[3]);
p2.SetMagThetaPhi(1.0, data.theta[4], data.phi[4]);
double theta_resid_break1 = std::acos(p1.Dot(p2));
std::string theta_resid_theta_break1_name = "theta_resid_theta_break1_cor_detected";
std::string theta_resid_theta_break1_title = theta_resid_theta_break1_name + ";#theta_{lab} Residual (deg);#theta_{lab} Breakup1 (deg)";
MyFill(theta_resid_theta_break1_name, theta_resid_theta_break1_title, data.theta[3]*rad2deg, data.theta[4]*rad2deg, 4);
std::string ke_break1_theta_rel_name = "ke_break1_theta_rel_detected";
std::string ke_break1_theta_rel_title = ke_break1_theta_rel_name + ";#theta_{resid-break1};KE_{break1} (MeV)";
MyFill(ke_break1_theta_rel_name, ke_break1_theta_rel_title, theta_resid_break1*rad2deg, data.KE[4], 4);
}
if((type == Mask::RxnType::TwoStepRxn || type == Mask::RxnType::ThreeStepRxn) && data.detect_flag[4] && data.detect_flag[5])
{
std::string theta_break1_theta_break2_name = "theta_break1_theta_break2_cor_detected";
std::string theta_break1_theta_break2_title = theta_break1_theta_break2_name + ";#theta_{lab} Breakup1 (deg);#theta_{lab} Breakup2 (deg)";
MyFill(theta_break1_theta_break2_name, theta_break1_theta_break2_title, data.theta[4]*rad2deg, data.theta[5]*rad2deg, 4);
}
if(type == Mask::RxnType::ThreeStepRxn && data.detect_flag[6] && data.detect_flag[7])
{
std::string theta_break3_theta_break4_name = "theta_break3_theta_break4_cor_detected";
std::string theta_break3_theta_break4_title = theta_break3_theta_break4_name + ";#theta_{lab} Breakup3 (deg);#theta_{lab} Breakup4 (deg)";
MyFill(theta_break3_theta_break4_name, theta_break3_theta_break4_title, data.theta[6]*rad2deg, data.theta[7]*rad2deg, 4);
}
}
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);
}
}
std::vector<Mask::Nucleus> GetParents(const Mask::MaskFileData& data, Mask::RxnType rxn_type)
{
std::vector<Mask::Nucleus> parents;
Mask::Nucleus temp1, temp2, temp3;
switch(rxn_type)
{
case Mask::RxnType::None : break;
case Mask::RxnType::PureDecay :
{
temp1.SetIsotope(data.Z[0], data.A[0]);
temp1.SetVectorSpherical(data.theta[0], data.phi[0], data.p[0], data.E[0]);
parents.push_back(temp1);
return parents;
}
case Mask::RxnType::OneStepRxn :
{
temp1.SetIsotope(data.Z[0], data.A[0]);
temp1.SetVectorSpherical(data.theta[0], data.phi[0], data.p[0], data.E[0]);
temp2.SetIsotope(data.Z[1], data.A[1]);
temp2.SetVectorSpherical(data.theta[1], data.phi[1], data.p[1], data.E[1]);
temp3 = temp1 + temp2;
parents.push_back(temp3);
return parents;
}
case Mask::RxnType::TwoStepRxn :
{
temp1.SetIsotope(data.Z[0], data.A[0]);
temp1.SetVectorSpherical(data.theta[0], data.phi[0], data.p[0], data.E[0]);
temp2.SetIsotope(data.Z[1], data.A[1]);
temp2.SetVectorSpherical(data.theta[1], data.phi[1], data.p[1], data.E[1]);
temp3 = temp1 + temp2;
parents.push_back(temp3);
temp3.SetIsotope(data.Z[3], data.A[3]);
temp3.SetVectorSpherical(data.theta[3], data.phi[3], data.p[3], data.E[3]);
parents.push_back(temp3);
return parents;
}
case Mask::RxnType::ThreeStepRxn :
{
temp1.SetIsotope(data.Z[0], data.A[0]);
temp1.SetVectorSpherical(data.theta[0], data.phi[0], data.p[0], data.E[0]);
temp2.SetIsotope(data.Z[1], data.A[1]);
temp2.SetVectorSpherical(data.theta[1], data.phi[1], data.p[1], data.E[1]);
temp3 = temp1 + temp2;
parents.push_back(temp3);
temp3.SetIsotope(data.Z[3], data.A[3]);
temp3.SetVectorSpherical(data.theta[3], data.phi[3], data.p[3], data.E[3]);
parents.push_back(temp3);
temp3.SetIsotope(data.Z[5], data.A[5]);
temp3.SetVectorSpherical(data.theta[5], data.phi[5], data.p[5], data.E[5]);
parents.push_back(temp3);
return parents;
}
}
return parents;
}
void SetThetaCM(Mask::Nucleus& daughter, const Mask::Nucleus& parent)
{
const double* boost = parent.GetBoost();
double boost2cm[3];
double boost2lab[3];
for(int i=0; i<3; i++)
{
boost2lab[i] = boost[i];
boost2cm[i] = -1.0*boost[i];
}
daughter.ApplyBoost(boost2cm);
daughter.SetThetaCM(daughter.GetTheta());
daughter.ApplyBoost(boost2lab);
}
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: "<<GetStringFromRxnType(header.rxn_type)<<" nsamples: "<<header.nsamples<<std::endl;
Mask::MaskFileData data;
Mask::Nucleus nucleus;
std::vector<Mask::Nucleus> parents; //for use with CM theta calc
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;
parents = GetParents(data, header.rxn_type);
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]);
/*
Little dirty, but works well. Save theta CM using parent from specific rxn step.
I.e. ejectile calculated from first parent, break1 from second parent, break3 from third...
*/
if(i==1 || i==2 || i==3)
{
SetThetaCM(nucleus, parents[0]);
}
else if (i==4 || i==5)
{
SetThetaCM(nucleus, parents[1]);
}
else if(i==6 || i==7)
{
SetThetaCM(nucleus, parents[2]);
}
plotter.FillData(nucleus);
if(data.detect_flag[i] == true) {
plotter.FillData(nucleus, data.KE[i], "detected");
}
}
plotter.FillCorrelations(data, header.rxn_type);
plotter.FillCorrelationsDetected(data, header.rxn_type);
count++;
}
std::cout<<std::endl;
input.Close();
root_out->cd();
plotter.GetTable()->Write();
root_out->Close();
return 0;
}

139
src/Plotters/RootPlotter.cpp Normal file
View File

@ -0,0 +1,139 @@
#include "RootPlotter.h"
#include <TFile.h>
#include <TTree.h>
#include <TVector3.h>
#include <iostream>
RootPlotter::RootPlotter()
{
TH1::AddDirectory(kFALSE);
//Enforce dictionary linking
if(Mask::EnforceDictionaryLinked())
{
std::cout<<"Dictionary Linked"<<std::endl;
}
}
RootPlotter::~RootPlotter() {}
void RootPlotter::Run(const std::string& inputname, const std::string& outputname)
{
TFile* input = TFile::Open(inputname.c_str(), "READ");
TTree* tree = (TTree*) input->Get("SimTree");
std::vector<Mask::Nucleus>* dataHandle = new std::vector<Mask::Nucleus>();
tree->SetBranchAddress("nuclei", &dataHandle);
TFile* output = TFile::Open(outputname.c_str(), "RECREATE");
double flushFrac = 0.05;
uint64_t nentries = tree->GetEntries();
uint64_t flushVal = flushFrac*nentries;
uint64_t count=0;
uint64_t flushCount = 0;
for(uint64_t i=0; i<nentries; i++)
{
tree->GetEntry(i);
for(Mask::Nucleus& nuc : *(dataHandle))
{
FillData(nuc);
}
}
input->Close();
delete dataHandle;
output->cd();
for(auto& obj : m_map)
obj.second->Write(obj.second->GetName(), TObject::kOverwrite);
output->Close();
}
void RootPlotter::FillData(const Mask::Nucleus& nuc)
{
std::string modifier = "";
if(nuc.isDetected)
modifier = "detected";
std::string sym = nuc.isotopicSymbol;
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 th_vs_ph_name = sym + modifier + "_theta_vs_phi";
std::string th_vs_ph_title = th_vs_ph_name + ";#theta_{lab};#phi_{lab}";
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";
if(nuc.isDetected)
{
MyFill(ke_vs_th_name.c_str(), ke_vs_th_title.c_str(), nuc.vec4.Theta()*s_rad2deg, nuc.GetKE(), 2);
MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), nuc.vec4.Phi()*s_rad2deg, nuc.GetKE(), 4);
MyFill(th_vs_ph_name.c_str(), th_vs_ph_title.c_str(), nuc.vec4.Theta()*s_rad2deg, nuc.vec4.Phi()*s_rad2deg, 2);
MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy());
MyFill(angdist_name.c_str(), angdist_title.c_str(),20,-1.0,1.0,std::cos(nuc.thetaCM));
}
else
{
MyFill(ke_vs_th_name.c_str(), ke_vs_th_title.c_str(), nuc.vec4.Theta()*s_rad2deg, nuc.detectedKE, 2);
MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), nuc.vec4.Phi()*s_rad2deg, nuc.detectedKE, 4);
MyFill(th_vs_ph_name.c_str(), th_vs_ph_title.c_str(), nuc.vec4.Theta()*s_rad2deg, nuc.vec4.Phi()*s_rad2deg, 2);
MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy());
MyFill(angdist_name.c_str(), angdist_title.c_str(),20,-1.0,1.0,std::cos(nuc.thetaCM));
}
}
void RootPlotter::MyFill(const std::string& name, const std::string& title, int bins, float min, float max, double val)
{
auto iter = m_map.find(name);
if(iter != m_map.end())
{
std::shared_ptr<TH1> h = std::static_pointer_cast<TH1>(iter->second);
h->Fill(val);
}
else
{
std::shared_ptr<TH1F> h = std::make_shared<TH1F>(name.c_str(), title.c_str(), bins, min, max);
h->Fill(val);
m_map[name] = 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)
{
auto iter = m_map.find(name);
if(iter != m_map.end())
{
std::shared_ptr<TH2> h = std::static_pointer_cast<TH2>(iter->second);
h->Fill(valx, valy);
}
else
{
std::shared_ptr<TH2F> h = std::make_shared<TH2F>(name.c_str(), title.c_str(), binsx, minx, maxx, binsy, miny, maxy);
h->Fill(valx, valy);
m_map[name] = h;
}
}
void RootPlotter::MyFill(const std::string& name, const std::string& title, double valx, double valy, int color)
{
auto iter = m_map.find(name);
if(iter != m_map.end())
{
std::shared_ptr<TGraph> g = std::static_pointer_cast<TGraph>(iter->second);
g->SetPoint(g->GetN(), valx, valy);
}
else
{
std::shared_ptr<TGraph> g = std::make_shared<TGraph>(1, &valx, &valy);
g->SetName(name.c_str());
g->SetTitle(title.c_str());
g->SetMarkerColor(color);
m_map[name] = g;
}
}

34
src/Plotters/RootPlotter.h Normal file
View File

@ -0,0 +1,34 @@
#ifndef ROOTPLOTTER_H
#define ROOTPLOTTER_H
#include <vector>
#include <string>
#include <unordered_map>
#include "Nucleus.h"
#include <TH1.h>
#include <TH2.h>
#include <TGraph.h>
class RootPlotter
{
public:
RootPlotter();
~RootPlotter();
void Run(const std::string& inputname, const std::string& outputname);
private:
void FillData(const Mask::Nucleus& nuc);
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::unordered_map<std::string, std::shared_ptr<TObject>> m_map;
static constexpr double s_rad2deg = 180.0/M_PI;
};
#endif

14
src/Plotters/main.cpp Normal file
View File

@ -0,0 +1,14 @@
#include "RootPlotter.h"
#include <iostream>
int main(int argc, char** argv)
{
if(argc != 3)
{
std::cerr<<"Root plotter requires two commandline arguments: path to input file and path to outputfile"<<std::endl;
return 1;
}
RootPlotter plotter;
plotter.Run(argv[1], argv[2]);
}