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

Fixed energy loss bugs, including case of backwards travel through a non-uniform target. Modified premake file to be more universally compatible for ROOT. Requires user specification of ROOT directories in the premake file.

Browse Source
This commit is contained in:
Gordon McCann 2021-09-23 11:42:14 -04:00
parent 850daa87d8
commit 51b64e0c97
16 changed files with 472 additions and 140 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
etc/convertedAnasenCoords.txt Normal file
View File

@ -0,0 +1,30 @@
For SX3's: det id | mid rho | mid z | mid phi | name
For QQQ's: det id | mid phi | z | name
0 5.49779 0 QQQ_0
1 0.785398 0.8028 QQQ_1
2 2.35619 0 QQQ_2
3 3.92699 0 QQQ_3
6 8.90601 6.25 0.785795 R1a_A
7 8.89871 6.25 0.262014 R1a_B
8 8.90354 6.25 6.02132 R1a_C
9 8.90247 6.25 5.49779 R1a_D
10 8.90354 6.25 4.97426 R1a_E
11 8.90354 6.25 4.45052 R1a_F
12 8.90247 6.25 3.92699 R1b_A
13 8.90354 6.25 3.40346 R1b_B
14 8.90354 6.25 2.87972 R1b_C
15 8.90247 6.25 2.35619 R1b_D
16 8.90354 6.25 1.83266 R1b_E
17 8.90354 6.25 1.30893 R1b_F
18 8.90601 18.65 0.785795 R2a_A
19 8.89871 18.65 0.262014 R2a_B
20 8.90354 18.65 6.02132 R2a_C
21 8.90247 18.65 5.49779 R2a_D
22 8.90354 18.65 4.97426 R2a_E
23 8.90354 18.65 4.45052 R2a_F
24 8.90247 18.65 3.92699 R2b_A
25 8.90354 18.65 3.40346 R2b_B
26 8.90354 18.65 2.87972 R2b_C
27 8.90247 18.65 2.35619 R2b_D
28 8.90354 18.65 1.83266 R2b_E
29 8.90354 18.65 1.30893 R2b_F

23
include/AnasenEfficiency.h
View File

@ -6,6 +6,16 @@
#include "DetectorEfficiency.h"
#include "StripDetector.h"
#include "QQQDetector.h"
#include "Target.h"
#include "Nucleus.h"
#include "MaskFile.h"
struct DetectorResult {
bool detectFlag = false;
Mask::Vec3 direction;
double energy_deposited = 0.0;
std::string det_name = "";
};
class AnasenEfficiency : public DetectorEfficiency {
public:
@ -16,14 +26,18 @@ public:
double RunConsistencyCheck() override;
private:
bool IsRing1(double theta, double phi);
bool IsRing2(double theta, double phi);
bool IsQQQ(double theta, double phi);
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, int rxn_type);
std::vector<StripDetector> m_Ring1, m_Ring2;
std::vector<QQQDetector> m_forwardQQQs;
std::vector<QQQDetector> m_backwardQQQs;
Mask::Target det_silicon;
/**** ANASEN geometry constants *****/
const int n_sx3_per_ring = 12;
const int n_qqq = 4;
@ -42,8 +56,9 @@ private:
const double ring_phi[12] = {0.785795, 0.262014, 6.02132, 5.49779, 4.97426, 4.45052, 3.92699, 3.40346, 2.87972, 2.35619, 1.83266, 1.30893};
/*************************/
static constexpr double threshold = 0.2; //MeV
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

3
include/LayeredTarget.h
View File

@ -15,6 +15,8 @@ Written by G.W. McCann Aug. 2020
#include <vector>
#include <string>
#include <random>
#include "RandomGenerator.h"
#include "Target.h"
namespace Mask {
@ -37,6 +39,7 @@ namespace Mask {
private:
std::vector<Target> layers;
std::string name;
std::uniform_real_distribution<double> m_fractional_depth_dist;
};
}

23
include/QQQDetector.h
View File

@ -10,36 +10,45 @@
#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 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; };
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); };
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; };
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;

23
include/StripDetector.h
View File

@ -17,21 +17,28 @@
#include "Vec3.h"
#include "Rotation.h"
#include "RandomGenerator.h"
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 void SetRandomNumberGenerator(std::mt19937* random) { m_random = random; };
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);
std::pair<int,double> 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;
@ -49,14 +56,14 @@ private:
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;
std::mt19937* m_random; //Not owned by StripDetector!
bool rndmFlag;
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); };
};
#endif

8
include/Target.h
View File

@ -27,10 +27,10 @@ namespace Mask {
~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 getEnergyLossHalf(int zp, int ap, double startEnergy, double angle);
double getReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double angle);
double getReverseEnergyLossHalf(int zp, int ap, double finalEnergy, double angle);
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]; }

13
input.txt
View File

@ -1,14 +1,13 @@
----------Data Information----------
OutputFile: 7Bedp_600keV_beam_centered_target_targetgap_BackQQQ_rndmCM_test.mask
OutputFile: /data1/gwm17/mask_tests/7Be12C_870keV_beam_50CD2.mask
SaveTree: yes
SavePlots: yes
----------Reaction Information----------
ReactionType: 2
ReactionType: 1
Z A (order is target, projectile, ejectile, break1, break3(if pure decay is target, ejectile))
1 2
6 12
4 7
4 7
1 1
2 4
----------Target Information----------
Name: test_targ
Layers: 1
@ -20,8 +19,8 @@ Z A Stoich
0
~
----------Sampling Information----------
NumberOfSamples: 10000
BeamMeanEnergy(MeV): 0.6 BeamEnergySigma(MeV): 0.0
NumberOfSamples: 100000
BeamMeanEnergy(MeV): 0.87 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

47
premake5.lua
View File

@ -43,36 +43,32 @@ project "RootPlot"
"include/*.h"
}
filter "system:windows"
systemversion "latest"
--ROOT cannot be located using the config tools, so we must query for the ROOTSYS env variable
--Have to use an if statement to hide this (@penguin for example doesn't have a ROOTSYS)
if os.host() == windows then
rootpath = os.getenv("ROOTSYS")
--User specified path to ROOT CERN libraries--
ROOTIncludepath = "/usr/include/root"
ROOTLibpath = "/usr/lib64/root"
includedirs {
"include",
"./",
rootpath .. "include"
}
includedirs {
"include"
}
links {
rootpath .. "lib/**.lib"
}
end
sysincludedirs {
ROOTIncludepath
}
libdirs {
ROOTLibpath
}
links {
"Gui", "Core", "Imt", "RIO", "Net", "Hist",
"Graf", "Graf3d", "Gpad", "ROOTDataFrame", "ROOTVecOps",
"Tree", "TreePlayer", "Rint", "Postscript", "Matrix",
"Physics", "MathCore", "Thread", "MultiProc", "m", "dl"
}
filter "system:macosx or linux"
includedirs {
"include",
"./"
}
buildoptions {
"`root-config --cflags`"
}
linkoptions {
"`root-config --glibs`"
"-pthread"
}
filter "configurations:Debug"
@ -98,6 +94,7 @@ project "DetectEff"
"src/Rotation.cpp",
"src/Target.cpp",
"src/EnergyLoss.cpp",
"src/RandomGenerator.cpp",
"include/*.h"
}

260
src/Detectors/AnasenEfficiency.cpp
View File

@ -1,20 +1,29 @@
#include "AnasenEfficiency.h"
#include "Kinematics.h"
#include "MaskFile.h"
#include <fstream>
#include <iomanip>
AnasenEfficiency::AnasenEfficiency() :
DetectorEfficiency()
DetectorEfficiency(), det_silicon(si_thickness)
{
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], -1.0*ring1_z, ring_rho[i]);
m_Ring2[i].TurnOnRandomizedCoordinates();
}
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();
}
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() {}
@ -196,47 +205,199 @@ double AnasenEfficiency::RunConsistencyCheck() {
}
bool AnasenEfficiency::IsRing1(double theta, double phi) {
DetectorResult AnasenEfficiency::IsRing1(Mask::Nucleus& nucleus) {
DetectorResult observation;
//Mask::Vec3 coords;
double thetaIncident, eloss, e_dep;
for(auto& sx3 : m_Ring1) {
auto result = sx3.GetChannelRatio(theta, phi);
auto result = sx3.GetChannelRatio(nucleus.GetTheta(), nucleus.GetPhi());
if(result.first != -1) {
return true;
//coords = sx3.GetHitCoordinates(result.first, result.second);
observation.detectFlag = true;
observation.direction = sx3.GetHitCoordinates(result.first, result.second);
thetaIncident = std::acos(observation.direction.Dot(sx3.GetNormRotated())/observation.direction.GetR());
if(thetaIncident > M_PI/2.0)
thetaIncident = M_PI - thetaIncident;
//e_dep = det_silicon.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
observation.det_name = "R1";
//return std::make_pair(true, e_dep);
return observation;
}
}
return false;
//return std::make_pair(false, 0.0);
return observation;
}
bool AnasenEfficiency::IsRing2(double theta, double phi) {
DetectorResult AnasenEfficiency::IsRing2(Mask::Nucleus& nucleus) {
DetectorResult observation;
//Mask::Vec3 coords;
double thetaIncident, eloss, e_dep;
for(auto& sx3 : m_Ring2) {
auto result = sx3.GetChannelRatio(theta, phi);
auto result = sx3.GetChannelRatio(nucleus.GetTheta(), nucleus.GetPhi());
if(result.first != -1) {
return true;
//coords = sx3.GetHitCoordinates(result.first, result.second);
observation.detectFlag = true;
observation.direction = sx3.GetHitCoordinates(result.first, result.second);
thetaIncident = std::acos(observation.direction.Dot(sx3.GetNormRotated())/observation.direction.GetR());
if(thetaIncident > M_PI/2.0)
thetaIncident = M_PI - thetaIncident;
//e_dep = det_silicon.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
observation.det_name = "R2";
//return std::make_pair(true, e_dep);
return observation;
}
}
return false;
return observation;
}
bool AnasenEfficiency::IsQQQ(double theta, double phi) {
DetectorResult AnasenEfficiency::IsQQQ(Mask::Nucleus& nucleus) {
DetectorResult observation;
//Mask::Vec3 coords;
double thetaIncident, eloss, e_dep;
for(auto& qqq : m_forwardQQQs) {
auto result = qqq.GetTrajectoryRingWedge(theta, phi);
auto result = qqq.GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi());
if(result.first != -1) {
return true;
//coords = qqq.GetHitCoordinates(result.first, result.second);
observation.detectFlag = true;
observation.direction = qqq.GetHitCoordinates(result.first, result.second);
thetaIncident = std::acos(observation.direction.Dot(qqq.GetNorm())/observation.direction.GetR());
if(thetaIncident > M_PI/2.0)
thetaIncident = M_PI - thetaIncident;
//e_dep = det_silicon.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
//return std::make_pair(true, e_dep);
observation.det_name = "FQQQ";
return observation;
}
}
for(auto& qqq : m_backwardQQQs) {
auto result = qqq.GetTrajectoryRingWedge(theta, phi);
auto result = qqq.GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi());
if(result.first != -1) {
return true;
//coords = qqq.GetHitCoordinates(result.first, result.second);
observation.detectFlag = true;
observation.direction = qqq.GetHitCoordinates(result.first, result.second);
thetaIncident = std::acos(observation.direction.Dot(qqq.GetNorm())/observation.direction.GetR());
if(thetaIncident > M_PI/2.0)
thetaIncident = M_PI - thetaIncident;
//e_dep = det_silicon.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
//return std::make_pair(true, e_dep);
observation.det_name = "BQQQ";
return observation;
}
}
//return std::make_pair(false, 0.0);
return observation;
}
DetectorResult AnasenEfficiency::IsAnasen(Mask::Nucleus& nucleus) {
DetectorResult result;
if(nucleus.GetKE() <= threshold)
return result;
return false;
if(!result.detectFlag) {
result = IsRing1(nucleus);
}
if(!result.detectFlag) {
result = IsRing2(nucleus);
}
if(!result.detectFlag) {
result = IsQQQ(nucleus);
}
return result;
}
void AnasenEfficiency::CountCoincidences(const Mask::MaskFileData& data, std::vector<int>& counts, int rxn_type) {
if (rxn_type == 0 && data.detect_flag[1] && data.detect_flag[2])
{
counts[0]++;
}
else if (rxn_type == 1 && data.detect_flag[2] && data.detect_flag[3])
{
counts[0]++;
}
else if(rxn_type == 2)
{
if(data.detect_flag[2] && data.detect_flag[4])
{
counts[0]++;
}
if(data.detect_flag[2] && data.detect_flag[5])
{
counts[1]++;
}
if(data.detect_flag[4] && data.detect_flag[5])
{
counts[2]++;
}
if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[5])
{
counts[3]++;
}
}
else if(rxn_type == 3)
{
if(data.detect_flag[2] && data.detect_flag[4])
{
counts[0]++;
}
if(data.detect_flag[2] && data.detect_flag[6])
{
counts[1]++;
}
if(data.detect_flag[2] && data.detect_flag[7])
{
counts[2]++;
}
if(data.detect_flag[4] && data.detect_flag[6])
{
counts[3]++;
}
if(data.detect_flag[4] && data.detect_flag[7])
{
counts[4]++;
}
if(data.detect_flag[6] && data.detect_flag[7])
{
counts[5]++;
}
if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[6])
{
counts[6]++;
}
if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[7])
{
counts[7]++;
}
if(data.detect_flag[2] && data.detect_flag[6] && data.detect_flag[7])
{
counts[8]++;
}
if(data.detect_flag[4] && data.detect_flag[6] && data.detect_flag[7])
{
counts[9]++;
}
if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[6] && data.detect_flag[7])
{
counts[10]++;
}
}
}
void AnasenEfficiency::CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) {
@ -257,18 +418,23 @@ void AnasenEfficiency::CalculateEfficiency(const std::string& inputname, const s
Mask::MaskFileData data;
std::vector<int> counts;
std::vector<int> coinc_counts;
switch(header.rxn_type) {
case 0:
counts.resize(3, 0);
coinc_counts.resize(1, 0);
break;
case 1:
counts.resize(4, 0);
coinc_counts.resize(1, 0);
break;
case 2:
counts.resize(6, 0);
coinc_counts.resize(4, 0);
break;
case 3:
counts.resize(8, 0);
coinc_counts.resize(11, 0);
break;
default:
{
@ -284,6 +450,8 @@ void AnasenEfficiency::CalculateEfficiency(const std::string& inputname, const s
int count = 0;
int npercent = 0;
Mask::Nucleus nucleus;
int index=0;
while(true) {
if(++count == percent5) {//Show progress every 5%
npercent++;
@ -296,15 +464,25 @@ void AnasenEfficiency::CalculateEfficiency(const std::string& inputname, const s
break;
for(unsigned int i=0; i<data.Z.size(); i++) {
if(data.KE[i] >= threshold && (IsRing1(data.theta[i], data.phi[i]) || IsRing2(data.theta[i], data.phi[i]) || IsQQQ(data.theta[i], data.phi[i]))) {
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;
}
}
CountCoincidences(data, coinc_counts, header.rxn_type);
output.WriteData(data);
index++;
}
input.Close();
@ -316,6 +494,54 @@ void AnasenEfficiency::CalculateEfficiency(const std::string& inputname, const s
stats<<std::setw(10)<<i<<"|"<<std::setw(10)<<((double)counts[i]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
}
stats<<"Coincidence Efficiency"<<std::endl;
stats<<"---------------------"<<std::endl;
if(header.rxn_type == 0)
{
stats<<std::setw(10)<<"1 + 2"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
}
else if(header.rxn_type == 1)
{
stats<<std::setw(10)<<"2 + 3"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
}
else if(header.rxn_type == 2)
{
stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/header.nsamples)<<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::endl;
}
else if(header.rxn_type == 3)
{
stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"2 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[4]/header.nsamples)<<std::endl;
stats<<"---------------------"<<std::endl;
stats<<std::setw(10)<<"6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[5]/header.nsamples)<<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::endl;
stats<<std::setw(10)<<"2 + 4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[7]/header.nsamples)<<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::endl;
stats<<std::setw(10)<<"4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[9]/header.nsamples)<<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::endl;
}
stats.close();
std::cout<<std::endl;

9
src/Detectors/DetectorEfficiency.cpp
View File

@ -1,5 +1,6 @@
#include "SabreEfficiency.h"
#include "AnasenEfficiency.h"
#include "KinematicsExceptions.h"
#include <iostream>
#include <string>
@ -23,11 +24,15 @@ int main(int argc, char** argv) {
//sabre.DrawDetectorSystem("/data1/gwm17/10B3He/Feb2021/simulation/SABREGeo.txt");
*/
try {
AnasenEfficiency anasen;
anasen.CalculateEfficiency(inputname, outputname, statsname);
//std::cout<<"Running consistency check(1=success): "<<anasen.RunConsistencyCheck()<<std::endl;
//anasen.DrawDetectorSystem("/data1/gwm17/TRIUMF_7Bed/simulation/ANASENGeo_centered_target_targetGap_BackQQQ_test.root");
} catch(const std::exception& e) {
std::cerr<<"Error: "<<e.what()<<std::endl;
std::cerr<<"Terminating."<<std::endl;
return 1;
}
return 0;
}

18
src/Detectors/QQQDetector.cpp
View File

@ -1,7 +1,7 @@
#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_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_deltaR = (m_Router - m_Rinner)/nrings;
m_deltaPhi_per_wedge = m_deltaPhi/nwedges;
@ -139,12 +139,20 @@ std::pair<int,int> QQQDetector::GetTrajectoryRingWedge(double theta, double phi)
}
Mask::Vec3 QQQDetector::GetHitCoordinates(int ringch, int wedgech) {
if(!CheckChannel(ringch) || !CheckChannel(wedgech)) {
if(!CheckChannel(ringch) || !CheckChannel(wedgech))
return Mask::Vec3();
}
double r_center = m_Rinner + (0.5+ringch)*m_deltaR;
double phi_center = -m_deltaPhi/2.0 + (0.5+wedgech)*m_deltaPhi_per_wedge;
double r_center, phi_center;
if(rndmFlag)
{
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;
}
else
{
r_center = m_Rinner + (0.5+ringch)*m_deltaR;
phi_center = -m_deltaPhi/2.0 + (0.5+wedgech)*m_deltaPhi_per_wedge;
}
double x = r_center*std::cos(phi_center);
double y = r_center*std::sin(phi_center);
double z = 0;

4
src/Detectors/SabreEfficiency.cpp
View File

@ -197,9 +197,9 @@ std::pair<bool,double> SabreEfficiency::IsSabre(Mask::Nucleus& nucleus) {
if(dmap.IsDead(i, chan.first, 0) || dmap.IsDead(i, chan.second, 1)) break; //dead channel check
coords = detectors[i].GetTrajectoryCoordinates(nucleus.GetTheta(), nucleus.GetPhi());
thetaIncident = std::acos(coords.Dot(detectors[i].GetNormTilted())/(coords.GetR()));
eloss = deadlayer.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), M_PI - thetaIncident);
eloss = deadlayer.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), M_PI - thetaIncident);
if((nucleus.GetKE() - eloss) <= ENERGY_THRESHOLD) break; //deadlayer check
e_deposited = sabre_eloss.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE() - eloss, M_PI - thetaIncident);
e_deposited = sabre_eloss.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE() - eloss, M_PI - thetaIncident);
return std::make_pair(true, e_deposited);
}
}

41
src/Detectors/StripDetector.cpp
View File

@ -1,19 +1,23 @@
#include "StripDetector.h"
#include <iostream>
/*
Corner layout for each strip in the un-rotated frame
0--------------------------1
| | y
| | ^
| | |
| | |
| | |
| |
| |
| |
| |
| | x
2--------------------------3 z<--------X
x
|
|
|
y
*/
StripDetector::StripDetector(int ns, double len, double wid, double cphi, double cz, double crho) :
m_random(nullptr), m_uniform_fraction(0.0, 1.0)
m_norm_unrot(1.0,0.0,0.0), m_uniform_fraction(0.0, 1.0), rndmFlag(false)
{
num_strips = ns;
@ -53,19 +57,19 @@ void StripDetector::CalculateCorners() {
y_min = -total_width/2.0 + front_strip_width*s;
z_max = center_z + length/2.0;
z_min = center_z - length/2.0;
front_strip_coords[s][0] = Mask::Vec3(center_rho, y_max, z_max);
front_strip_coords[s][1] = Mask::Vec3(center_rho, y_max, z_min);
front_strip_coords[s][2] = Mask::Vec3(center_rho, y_min, z_max);
front_strip_coords[s][3] = Mask::Vec3(center_rho, y_min, z_min);
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);
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][0] = Mask::Vec3(center_rho, y_max, z_max);
back_strip_coords[s][1] = Mask::Vec3(center_rho, y_max, z_min);
back_strip_coords[s][2] = Mask::Vec3(center_rho, y_min, z_max);
back_strip_coords[s][3] = Mask::Vec3(center_rho, y_min, z_min);
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);
}
for(int s=0; s<num_strips; s++) {
@ -86,9 +90,8 @@ Mask::Vec3 StripDetector::GetHitCoordinates(int front_stripch, double front_stri
if (!ValidChannel(front_stripch) || !ValidRatio(front_strip_ratio)) return Mask::Vec3(0,0,0);
double y;
//If we have a random number generator, randomize y position within pixel. Otherwise take halfway.
if(m_random) {
y = -total_width/2.0 + (front_stripch + m_uniform_fraction(*m_random))*front_strip_width;
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;
}
@ -128,7 +131,7 @@ std::pair<int,double> StripDetector::GetChannelRatio(double theta, double phi) {
double ratio=0;
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][2].GetY() && yhit <= front_strip_coords[s][0].GetY() && //Check min and max y
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
{
front_ch_hit = s;

68
src/LayeredTarget.cpp
View File

@ -16,7 +16,7 @@ Written by G.W. McCann Aug. 2020
namespace Mask {
LayeredTarget::LayeredTarget() :
name("")
name(""), m_fractional_depth_dist(0.0, 1.0)
{
}
@ -31,7 +31,7 @@ namespace Mask {
/*
Here projectile refers to the incoming reactant particle (i.e. the beam)
Calculates energy loss assuming that the reaction occurs in the middle of the target
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) {
@ -43,12 +43,14 @@ namespace Mask {
double eloss = 0.0;
double newEnergy = startEnergy;
double frac;
for(int i=0; i<=rxnLayer; i++) {
if(i == rxnLayer) {
eloss += layers[i].getEnergyLossHalf(zp, ap, newEnergy, angle);
frac = m_fractional_depth_dist(RandomGenerator::GetInstance().GetGenerator());
eloss += layers[i].GetEnergyLossFractionalDepth(zp, ap, newEnergy, angle, frac);
newEnergy = startEnergy - eloss;
} else {
eloss += layers[i].getEnergyLossTotal(zp, ap, newEnergy, angle);
eloss += layers[i].GetEnergyLossTotal(zp, ap, newEnergy, angle);
newEnergy = startEnergy-eloss;
}
}
@ -69,14 +71,27 @@ namespace Mask {
}
double eloss = 0.0;
double newEnergy = startEnergy;
for(unsigned int i=rxnLayer; i<layers.size(); i++) {
if(i == ((unsigned int)rxnLayer)) {
eloss += layers[i].getEnergyLossHalf(ze, ae, newEnergy, angle);
newEnergy = startEnergy - eloss;
} else {
eloss += layers[i].getEnergyLossTotal(ze, ae, newEnergy, angle);
newEnergy = startEnergy - eloss;
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;
}
}
} 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;
}
}
}
@ -92,14 +107,27 @@ namespace Mask {
}
double eloss = 0.0;
double newEnergy = startEnergy;
for(int i=(layers.size()-1); i>=rxnLayer; i--) {
if(i == rxnLayer) {
eloss += layers[i].getReverseEnergyLossHalf(ze, ae, newEnergy, angle);
newEnergy = startEnergy + eloss;
} else {
eloss += layers[i].getReverseEnergyLossTotal(ze, ae, newEnergy, angle);
newEnergy = startEnergy + eloss;
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;
}
}
} 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;
}
}
}

38
src/Target.cpp
View File

@ -39,7 +39,7 @@ namespace Mask {
}
/*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.)
@ -47,19 +47,20 @@ namespace Mask {
return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/fabs(cos(theta)));
}
/*Calculates energy loss for travelling halfway through the target*/
double Target::getEnergyLossHalf(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;
return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/(2.0*fabs(cos(theta))));
/*Calculates the energy loss for traveling some fraction through the target*/
double Target::GetEnergyLossFractionalDepth(int zp, int ap, double finalEnergy, double theta, double percent_depth)
{
if(theta == M_PI/2.)
return finalEnergy;
else if (theta > M_PI/2.)
theta = M_PI-theta;
return eloss.GetEnergyLoss(zp, ap, finalEnergy, thickness*percent_depth/(std::fabs(std::cos(theta))));
}
/*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.)
@ -67,15 +68,16 @@ namespace Mask {
return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/fabs(cos(theta)));
}
/*Calculates reverse energy loss for travelling half way through the target*/
double Target::getReverseEnergyLossHalf(int zp, int ap, double finalEnergy, double theta) {
if(theta == M_PI/2.)
return finalEnergy;
else if (theta > M_PI/2.)
theta = M_PI - theta;
return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/(2.0*fabs(cos(theta))));
/*Calculates the reverse energy loss for traveling some fraction through the target*/
double Target::GetReverseEnergyLossFractionalDepth(int zp, int ap, double finalEnergy, double theta, double percent_depth)
{
if(theta == M_PI/2.)
return finalEnergy;
else if (theta > M_PI/2.)
theta = M_PI-theta;
return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness*percent_depth/(std::fabs(std::cos(theta))));
}
}

4
src/TwoStepSystem.cpp
View File

@ -73,7 +73,7 @@ namespace Mask {
double decay1costheta = decay1dist.GetRandomCosTheta();
double decay1Theta = std::acos(decay1costheta);
double decay1Phi = m_phi2Range(RandomGenerator::GetInstance().GetGenerator());
double residEx = (*m_beamDist)(RandomGenerator::GetInstance().GetGenerator());
double residEx = (*m_exDist)(RandomGenerator::GetInstance().GetGenerator());
step1.SetBeamKE(bke);
step1.SetPolarRxnAngle(rxnTheta);
@ -97,4 +97,4 @@ namespace Mask {
}
}
}