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

Fix some dead channel related bugs. Preparing for a code base update.

Browse Source
This commit is contained in:
Gordon McCann 2022-08-15 11:03:57 -04:00
parent 0b63196077
commit 59c0dc99a7
7 changed files with 36 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
etc/sabreDeadChannels_May2022.txt Normal file
View File

@ -0,0 +1,4 @@
DetectorID RING/WEDGE Channel
4 WEDGE 4
3 RING 2
1 RING 13

4
include/SabreDetector.h
View File

@ -62,7 +62,7 @@ class SabreDetector {
public:
SabreDetector();
SabreDetector(double Rin, double Rout, double deltaPhi_flat, double phiCentral, double tiltFromVert, double zdist, double xdist=0, double ydist=0);
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*/
@ -84,6 +84,7 @@ public:
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:
@ -165,6 +166,7 @@ private:
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;

10
include/SabreEfficiency.h
View File

@ -34,17 +34,17 @@ private:
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 = 234.0; //center phi values for each det in array
const double PHI1 = 162.0; //# is equal to detID in channel map
const double PHI2 = 306.0;
const double PHI3 = 18.0;
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.2; //in MeV
const double ENERGY_THRESHOLD = 0.25; //in MeV
};

2
src/Detectors/DetectorEfficiency.cpp
View File

@ -17,7 +17,7 @@ int main(int argc, char** argv) {
SabreEfficiency sabre;
std::string mapfile = "./etc/SabreDeadChannels.txt";
std::string mapfile = "./etc/sabreDeadChannels_May2022.txt";
sabre.SetDeadChannelMap(mapfile);
sabre.CalculateEfficiency(inputname, outputname, statsname);
//std::cout<<"Running consistency check(1=success): "<<sabre.RunConsistencyCheck()<<std::endl;

8
src/Detectors/SabreDetector.cpp
View File

@ -53,7 +53,8 @@
#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_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_YRot.SetAngle(m_tilt);
m_ZRot.SetAngle(m_phiCentral);
@ -78,8 +79,9 @@ m_Router(0.1351), m_Rinner(0.0326), m_deltaPhi_flat(54.4*deg2rad), m_phiCentral(
CalculateCorners();
}
SabreDetector::SabreDetector(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)
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)
{
m_YRot.SetAngle(m_tilt);
m_ZRot.SetAngle(m_phiCentral);

31
src/Detectors/SabreEfficiency.cpp
View File

@ -8,13 +8,12 @@
SabreEfficiency::SabreEfficiency() :
DetectorEfficiency(), deadlayer(DEADLAYER_THIN), sabre_eloss(SABRE_THICKNESS), degrader(DEGRADER_THICKNESS)
{
detectors.reserve(5);
detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI0*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI1*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
//Only 0,1,4 valid in degrader land
//detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI2*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
//detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI3*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI4*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
//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};
@ -103,6 +102,7 @@ void SabreEfficiency::CalculateEfficiency(const std::string& inputname, const st
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;
@ -253,16 +253,17 @@ std::pair<bool,double> SabreEfficiency::IsSabre(Mask::Nucleus& nucleus) {
Mask::Vec3 coords;
double thetaIncident, eloss, e_deposited;
double ke = 0.0;
for(int i=0; i<5; i++) {
auto chan = detectors[i].GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi());
for(auto& detector : detectors) {
auto chan = detector.GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi());
if(chan.first != -1 && chan.second != -1) {
if(dmap.IsDead(i, chan.first, 0) || dmap.IsDead(i, chan.second, 1)) break; //dead channel check
coords = detectors[i].GetTrajectoryCoordinates(nucleus.GetTheta(), nucleus.GetPhi());
thetaIncident = std::acos(coords.Dot(detectors[i].GetNormTilted())/(coords.GetR()));
eloss = degrader.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident);
ke = nucleus.GetKE() - eloss;
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 -= eloss;
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);

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

@ -30,7 +30,7 @@ void RootPlotter::FillData(const Mask::Nucleus& nuc, double detKE, const std::st
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(),100,-1.0,1.0,std::cos(nuc.GetThetaCM()));
MyFill(angdist_name.c_str(), angdist_title.c_str(),20,-1.0,1.0,std::cos(nuc.GetThetaCM()));
}
else
{
@ -38,7 +38,7 @@ void RootPlotter::FillData(const Mask::Nucleus& nuc, double detKE, const std::st
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(),100,-1.0,1.0,std::cos(nuc.GetThetaCM()));
MyFill(angdist_name.c_str(), angdist_title.c_str(),20,-1.0,1.0,std::cos(nuc.GetThetaCM()));
}
}