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

Fix up detector geometry, make names clearer for dets and constexpr more of the geometry

Browse Source
This commit is contained in:
Gordon McCann 2022-08-30 15:01:11 -04:00
parent 45a80ab76f
commit bb67011bb5
11 changed files with 123 additions and 163 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/Detectors/AnasenEfficiency.cpp
View File

@ -11,16 +11,16 @@ AnasenEfficiency::AnasenEfficiency() :
{ {
for(int i=0; i<s_nSX3PerBarrel; i++) 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.emplace_back(s_barrelPhiList[i], s_barrel1Z, s_barrelRhoList[i]);
m_Ring1[i].SetPixelSmearing(true); m_Ring1[i].SetPixelSmearing(true);
m_Ring2.emplace_back(4, s_sx3Length, s_sx3Width, s_barrelPhiList[i], s_barrel2Z, s_barrelRhoList[i]); m_Ring2.emplace_back(s_barrelPhiList[i], s_barrel2Z, s_barrelRhoList[i]);
m_Ring2[i].SetPixelSmearing(true); m_Ring2[i].SetPixelSmearing(true);
} }
for(int i=0; i<s_nQQQ; i++) 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.emplace_back(s_qqqPhiList[i], s_qqqZList[i]);
m_forwardQQQs[i].SetSmearing(true); 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.emplace_back(s_qqqPhiList[i], (-1.0)*s_qqqZList[i]);
m_backwardQQQs[i].SetSmearing(true); m_backwardQQQs[i].SetSmearing(true);
} }
} }

6
src/Detectors/AnasenEfficiency.h
View File

@ -4,7 +4,7 @@
#include <string> #include <string>
#include "DetectorEfficiency.h" #include "DetectorEfficiency.h"
#include "StripDetector.h" #include "SX3Detector.h"
#include "QQQDetector.h" #include "QQQDetector.h"
#include "Target.h" #include "Target.h"
#include "Nucleus.h" #include "Nucleus.h"
@ -39,15 +39,11 @@ private:
static constexpr int s_nSX3PerBarrel = 12; static constexpr int s_nSX3PerBarrel = 12;
static constexpr int s_nQQQ = 4; static constexpr int s_nQQQ = 4;
static constexpr double s_sx3Length = 0.075; 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_barrelGap = 0.0254;
static constexpr double s_sx3FrameGap = 0.049; //0.049 is base gap due to frames 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_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_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_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_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_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, static constexpr double s_barrelRhoList[12] = {0.0890601, 0.0889871, 0.0890354, 0.0890247, 0.0890354, 0.0890354, 0.0890247,

4
src/Detectors/CMakeLists.txt
View File

@ -16,8 +16,8 @@ target_sources(Detectors PUBLIC
SabreDetector.h SabreDetector.h
SabreEfficiency.cpp SabreEfficiency.cpp
SabreEfficiency.h SabreEfficiency.h
StripDetector.cpp SX3Detector.cpp
StripDetector.h SX3Detector.h
) )
target_link_libraries(Detectors target_link_libraries(Detectors

47
src/Detectors/QQQDetector.cpp
View File

@ -1,11 +1,8 @@
#include "QQQDetector.h" #include "QQQDetector.h"
QQQDetector::QQQDetector(double R_in, double R_out, double deltaPhi, double phiCentral, double z, double x, double y) : QQQDetector::QQQDetector(double phiCentral, double zOffset, double xOffset, double yOffset) :
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_centralPhi(phiCentral), m_translation(xOffset,yOffset,zOffset), m_norm(0.0,0.0,1.0), m_uniformFraction(0.0, 1.0), m_isSmearing(false)
m_uniformFraction(0.0, 1.0), m_isSmearing(false)
{ {
m_deltaR = (m_outerR - m_innerR)/s_nRings;
m_deltaPhiWedge = m_deltaPhi/s_nWedges;
m_zRotation.SetAngle(m_centralPhi); m_zRotation.SetAngle(m_centralPhi);
m_ringCoords.resize(s_nRings); m_ringCoords.resize(s_nRings);
m_wedgeCoords.resize(s_nWedges); m_wedgeCoords.resize(s_nWedges);
@ -28,23 +25,23 @@ void QQQDetector::CalculateCorners()
//Generate flat ring corner coordinates //Generate flat ring corner coordinates
for(int i=0; i<s_nRings; i++) for(int i=0; i<s_nRings; i++)
{ {
x0 = (m_innerR + m_deltaR*(i+1))*std::cos(-m_deltaPhi/2.0); x0 = (s_innerR + s_deltaR*(i+1))*std::cos(-s_deltaPhiTotal/2.0);
y0 = (m_innerR + m_deltaR*(i+1))*std::sin(-m_deltaPhi/2.0); y0 = (s_innerR + s_deltaR*(i+1))*std::sin(-s_deltaPhiTotal/2.0);
z0 = 0.0; z0 = 0.0;
m_ringCoords[i][0].SetXYZ(x0, y0, z0); m_ringCoords[i][0].SetXYZ(x0, y0, z0);
x1 = (m_innerR + m_deltaR*(i))*std::cos(-m_deltaPhi/2.0); x1 = (s_innerR + s_deltaR*(i))*std::cos(-s_deltaPhiTotal/2.0);
y1 = (m_innerR + m_deltaR*(i))*std::sin(-m_deltaPhi/2.0); y1 = (s_innerR + s_deltaR*(i))*std::sin(-s_deltaPhiTotal/2.0);
z1 = 0.0; z1 = 0.0;
m_ringCoords[i][1].SetXYZ(x1, y1, z1); m_ringCoords[i][1].SetXYZ(x1, y1, z1);
x2 = (m_innerR + m_deltaR*(i))*std::cos(m_deltaPhi/2.0); x2 = (s_innerR + s_deltaR*(i))*std::cos(s_deltaPhiTotal/2.0);
y2 = (m_innerR + m_deltaR*(i))*std::sin(m_deltaPhi/2.0); y2 = (s_innerR + s_deltaR*(i))*std::sin(s_deltaPhiTotal/2.0);
z2 = 0.0; z2 = 0.0;
m_ringCoords[i][2].SetXYZ(x2, y2, z2); m_ringCoords[i][2].SetXYZ(x2, y2, z2);
x3 = (m_innerR + m_deltaR*(i+1))*std::cos(m_deltaPhi/2.0); x3 = (s_innerR + s_deltaR*(i+1))*std::cos(s_deltaPhiTotal/2.0);
y3 = (m_innerR + m_deltaR*(i+1))*std::sin(m_deltaPhi/2.0); y3 = (s_innerR + s_deltaR*(i+1))*std::sin(s_deltaPhiTotal/2.0);
z3 = 0.0; z3 = 0.0;
m_ringCoords[i][3].SetXYZ(x3, y3, z3); m_ringCoords[i][3].SetXYZ(x3, y3, z3);
} }
@ -52,23 +49,23 @@ void QQQDetector::CalculateCorners()
//Generate flat wedge corner coordinates //Generate flat wedge corner coordinates
for(int i=0; i<s_nWedges; i++) for(int i=0; i<s_nWedges; i++)
{ {
x0 = m_outerR * std::cos(-m_deltaPhi/2.0 + i*m_deltaPhiWedge); x0 = s_outerR * std::cos(-s_deltaPhiTotal/2.0 + i*s_deltaPhi);
y0 = m_outerR * std::sin(-m_deltaPhi/2.0 + i*m_deltaPhiWedge); y0 = s_outerR * std::sin(-s_deltaPhiTotal/2.0 + i*s_deltaPhi);
z0 = 0.0; z0 = 0.0;
m_wedgeCoords[i][0].SetXYZ(x0, y0, z0); m_wedgeCoords[i][0].SetXYZ(x0, y0, z0);
x1 = m_innerR * std::cos(-m_deltaPhi/2.0 + i*m_deltaPhiWedge); x1 = s_innerR * std::cos(-s_deltaPhiTotal/2.0 + i*s_deltaPhi);
y1 = m_innerR * std::sin(-m_deltaPhi/2.0 + i*m_deltaPhiWedge); y1 = s_innerR * std::sin(-s_deltaPhiTotal/2.0 + i*s_deltaPhi);
z1 = 0.0; z1 = 0.0;
m_wedgeCoords[i][1].SetXYZ(x1, y1, z1); m_wedgeCoords[i][1].SetXYZ(x1, y1, z1);
x2 = m_innerR * std::cos(-m_deltaPhi/2.0 + (i+1)*m_deltaPhiWedge); x2 = s_innerR * std::cos(-s_deltaPhiTotal/2.0 + (i+1)*s_deltaPhi);
y2 = m_innerR * std::sin(-m_deltaPhi/2.0 + (i+1)*m_deltaPhiWedge); y2 = s_innerR * std::sin(-s_deltaPhiTotal/2.0 + (i+1)*s_deltaPhi);
z2 = 0.0; z2 = 0.0;
m_wedgeCoords[i][2].SetXYZ(x2, y2, z2); m_wedgeCoords[i][2].SetXYZ(x2, y2, z2);
x3 = m_outerR * std::cos(-m_deltaPhi/2.0 + (i+1)*m_deltaPhiWedge); x3 = s_outerR * std::cos(-s_deltaPhiTotal/2.0 + (i+1)*s_deltaPhi);
y3 = m_outerR * std::sin(-m_deltaPhi/2.0 + (i+1)*m_deltaPhiWedge); y3 = s_outerR * std::sin(-s_deltaPhiTotal/2.0 + (i+1)*s_deltaPhi);
z3 = 0.0; z3 = 0.0;
m_wedgeCoords[i][3].SetXYZ(x3, y3, z3); m_wedgeCoords[i][3].SetXYZ(x3, y3, z3);
} }
@ -156,13 +153,13 @@ ROOT::Math::XYZPoint QQQDetector::GetHitCoordinates(int ringch, int wedgech)
double r_center, phi_center; double r_center, phi_center;
if(m_isSmearing) if(m_isSmearing)
{ {
r_center = m_innerR + (m_uniformFraction(Mask::RandomGenerator::GetInstance().GetGenerator())+ringch)*m_deltaR; r_center = s_innerR + (m_uniformFraction(Mask::RandomGenerator::GetInstance().GetGenerator())+ringch)*s_deltaR;
phi_center = -m_deltaPhi/2.0 + (m_uniformFraction(Mask::RandomGenerator::GetInstance().GetGenerator())+wedgech)*m_deltaPhiWedge; phi_center = -s_deltaPhiTotal/2.0 + (m_uniformFraction(Mask::RandomGenerator::GetInstance().GetGenerator())+wedgech)*s_deltaPhi;
} }
else else
{ {
r_center = m_innerR + (0.5+ringch)*m_deltaR; r_center = s_innerR + (0.5+ringch)*s_deltaR;
phi_center = -m_deltaPhi/2.0 + (0.5+wedgech)*m_deltaPhiWedge; phi_center = -s_deltaPhiTotal/2.0 + (0.5+wedgech)*s_deltaPhi;
} }
double x = r_center*std::cos(phi_center); double x = r_center*std::cos(phi_center);
double y = r_center*std::sin(phi_center); double y = r_center*std::sin(phi_center);

12
src/Detectors/QQQDetector.h
View File

@ -20,7 +20,7 @@
class QQQDetector class QQQDetector
{ {
public: public:
QQQDetector(double R_in, double R_out, double deltaPhi, double phiCentral, double z, double x=0, double y=0); QQQDetector(double phiCentral, double zOffset, double xOffset=0, double yOffset=0);
~QQQDetector(); ~QQQDetector();
const ROOT::Math::XYZPoint& GetRingCoordinates(int ringch, int corner) { return m_ringCoords[ringch][corner]; } 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::XYZPoint& GetWedgeCoordinates(int wedgech, int corner) { return m_wedgeCoords[wedgech][corner]; }
@ -42,11 +42,6 @@ private:
void CalculateCorners(); void CalculateCorners();
ROOT::Math::XYZPoint TransformCoordinates(ROOT::Math::XYZPoint& vector) { return m_translation * (m_zRotation * vector) ; } 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; double m_centralPhi;
std::vector<std::vector<ROOT::Math::XYZPoint>> m_ringCoords, m_wedgeCoords; std::vector<std::vector<ROOT::Math::XYZPoint>> m_ringCoords, m_wedgeCoords;
@ -60,6 +55,11 @@ private:
static constexpr int s_nRings = 16; static constexpr int s_nRings = 16;
static constexpr int s_nWedges = 16; static constexpr int s_nWedges = 16;
static constexpr double s_deg2rad = M_PI/180.0; static constexpr double s_deg2rad = M_PI/180.0;
static constexpr double s_innerR = 0.0501;
static constexpr double s_outerR = 0.0990;
static constexpr double s_deltaR = (s_outerR - s_innerR) / s_nRings;
static constexpr double s_deltaPhiTotal = 1.52119;
static constexpr double s_deltaPhi = s_deltaPhiTotal / s_nWedges;
}; };
#endif #endif

65
src/Detectors/StripDetector.cpp → src/Detectors/SX3Detector.cpp
View File

@ -1,5 +1,4 @@
#include "StripDetector.h" #include "SX3Detector.h"
#include <iostream>
/* /*
Corner layout for each strip in the un-rotated frame Corner layout for each strip in the un-rotated frame
@ -16,30 +15,16 @@
y y
*/ */
StripDetector::StripDetector(int ns, double len, double wid, double cphi, double cz, double crho) : StripDetector::StripDetector(double centerPhi, double centerZ, double centerRho) :
m_norm(1.0,0.0,0.0), m_uniformFraction(0.0, 1.0), m_isSmearing(false) m_centerPhi(centerPhi), m_centerZ(centerZ), m_centerRho(centerRho), m_norm(1.0,0.0,0.0), m_uniformFraction(0.0, 1.0), m_isSmearing(false)
{ {
m_nStrips = ns;
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;
m_centerPhi = cphi;
m_centerZ = cz;
m_centerRho = std::fabs(crho);
m_zRotation.SetAngle(m_centerPhi); m_zRotation.SetAngle(m_centerPhi);
m_frontStripCoords.resize(m_nStrips); m_frontStripCoords.resize(s_nStrips);
m_backStripCoords.resize(m_nStrips); m_backStripCoords.resize(s_nStrips);
m_rotFrontStripCoords.resize(m_nStrips); m_rotFrontStripCoords.resize(s_nStrips);
m_rotBackStripCoords.resize(m_nStrips); m_rotBackStripCoords.resize(s_nStrips);
for(int i=0; i<m_nStrips; i++) for(int i=0; i<s_nStrips; i++)
{ {
m_frontStripCoords[i].resize(s_nCorners); m_frontStripCoords[i].resize(s_nCorners);
m_backStripCoords[i].resize(s_nCorners); m_backStripCoords[i].resize(s_nCorners);
@ -55,28 +40,28 @@ StripDetector::~StripDetector() {}
void StripDetector::CalculateCorners() void StripDetector::CalculateCorners()
{ {
double y_min, y_max, z_min, z_max; double y_min, y_max, z_min, z_max;
for (int s=0; s<m_nStrips; s++) for (int s=0; s<s_nStrips; s++)
{ {
y_max = m_totalWidth/2.0 - m_frontStripWidth*s; y_max = s_totalWidth/2.0 - s_frontStripWidth*s;
y_min = m_totalWidth/2.0 - m_frontStripWidth*(s+1); y_min = s_totalWidth/2.0 - s_frontStripWidth*(s+1);
z_max = m_centerZ + m_stripLength/2.0; z_max = m_centerZ + s_totalLength/2.0;
z_min = m_centerZ - m_stripLength/2.0; z_min = m_centerZ - s_totalLength/2.0;
m_frontStripCoords[s][2].SetXYZ(m_centerRho, y_max, z_max); 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][3].SetXYZ(m_centerRho, y_max, z_min);
m_frontStripCoords[s][0].SetXYZ(m_centerRho, y_min, z_max); m_frontStripCoords[s][0].SetXYZ(m_centerRho, y_min, z_max);
m_frontStripCoords[s][1].SetXYZ(m_centerRho, y_min, z_min); m_frontStripCoords[s][1].SetXYZ(m_centerRho, y_min, z_min);
z_max = (m_centerZ - m_stripLength/2.0) + (s+1)*m_backStripLength; z_max = (m_centerZ - s_totalLength/2.0) + (s+1)*s_backStripLength;
z_min = (m_centerZ - m_stripLength/2.0) + (s)*m_backStripLength; z_min = (m_centerZ - s_totalLength/2.0) + (s)*s_backStripLength;
y_max = m_totalWidth/2.0; y_max = s_totalWidth/2.0;
y_min = -m_totalWidth/2.0; y_min = -s_totalWidth/2.0;
m_backStripCoords[s][2].SetXYZ(m_centerRho, y_max, z_max); 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][3].SetXYZ(m_centerRho, y_max, z_min);
m_backStripCoords[s][0].SetXYZ(m_centerRho, y_min, z_max); m_backStripCoords[s][0].SetXYZ(m_centerRho, y_min, z_max);
m_backStripCoords[s][1].SetXYZ(m_centerRho, y_min, z_min); m_backStripCoords[s][1].SetXYZ(m_centerRho, y_min, z_min);
} }
for(int s=0; s<m_nStrips; s++) for(int s=0; s<s_nStrips; s++)
{ {
m_rotFrontStripCoords[s][0] = m_zRotation*m_frontStripCoords[s][0]; m_rotFrontStripCoords[s][0] = m_zRotation*m_frontStripCoords[s][0];
m_rotFrontStripCoords[s][1] = m_zRotation*m_frontStripCoords[s][1]; m_rotFrontStripCoords[s][1] = m_zRotation*m_frontStripCoords[s][1];
@ -98,12 +83,12 @@ ROOT::Math::XYZPoint StripDetector::GetHitCoordinates(int front_stripch, double
double y; double y;
if(m_isSmearing) if(m_isSmearing)
y = -m_totalWidth/2.0 + (front_stripch + m_uniformFraction(Mask::RandomGenerator::GetInstance().GetGenerator()))*m_frontStripWidth; y = -s_totalWidth/2.0 + (front_stripch + m_uniformFraction(Mask::RandomGenerator::GetInstance().GetGenerator()))*s_frontStripWidth;
else else
y = -m_totalWidth/2.0 + (front_stripch+0.5)*m_frontStripWidth; y = -s_totalWidth/2.0 + (front_stripch+0.5)*s_frontStripWidth;
//recall we're still assuming phi=0 det: //recall we're still assuming phi=0 det:
ROOT::Math::XYZPoint coords(m_centerRho, y, front_strip_ratio*(m_stripLength/2) + m_centerZ); ROOT::Math::XYZPoint coords(m_centerRho, y, front_strip_ratio*(s_totalLength/2) + m_centerZ);
//NOW rotate by appropriate phi //NOW rotate by appropriate phi
return m_zRotation*coords; return m_zRotation*coords;
@ -125,7 +110,7 @@ StripHit StripDetector::GetChannelRatio(double theta, double phi)
phi -= 2*M_PI; phi -= 2*M_PI;
//then we can check easily whether it even hit the detector in phi //then we can check easily whether it even hit the detector in phi
double det_max_phi = std::atan2(m_totalWidth/2,m_centerRho); double det_max_phi = std::atan2(s_totalWidth/2,m_centerRho);
double det_min_phi = -det_max_phi; double det_min_phi = -det_max_phi;
if (phi < det_min_phi || phi > det_max_phi) if (phi < det_min_phi || phi > det_max_phi)
@ -138,18 +123,18 @@ StripHit StripDetector::GetChannelRatio(double theta, double phi)
double yhit = xhit*tan(phi); double yhit = xhit*tan(phi);
double zhit = sqrt(xhit*xhit+yhit*yhit)/tan(theta); double zhit = sqrt(xhit*xhit+yhit*yhit)/tan(theta);
for (int s=0; s<m_nStrips; s++) { for (int s=0; s<s_nStrips; s++) {
if (xhit >=m_frontStripCoords[s][0].X() && xhit <=m_frontStripCoords[s][0].X() && //Check min and max x (constant in flat) 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 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 zhit >=m_frontStripCoords[s][1].Z() && zhit <=m_frontStripCoords[s][0].Z()) //Check min and max z
{ {
hit.front_strip_index = s; hit.front_strip_index = s;
hit.front_ratio = (zhit-m_centerZ)/(m_stripLength/2); hit.front_ratio = (zhit-m_centerZ)/(s_totalLength/2);
break; break;
} }
} }
for (int s=0; s<m_nStrips; s++) { for (int s=0; s<s_nStrips; s++) {
if (xhit >= m_backStripCoords[s][0].X() && xhit <= m_backStripCoords[s][0].X() && //Check min and max x (constant in flat) 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 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 zhit >= m_backStripCoords[s][1].Z() && zhit <= m_backStripCoords[s][0].Z()) //Check min and max z

25
src/Detectors/StripDetector.h → src/Detectors/SX3Detector.h
View File

@ -1,5 +1,5 @@
#ifndef STRIP_DETECTOR_H #ifndef SX3_DETECTOR_H
#define STRIP_DETECTOR_H #define SX3_DETECTOR_H
// +z is along beam axis // +z is along beam axis
// +y is vertically "downward" in the lab frame // +y is vertically "downward" in the lab frame
@ -30,7 +30,7 @@ class StripDetector
{ {
public: public:
StripDetector(int ns, double len, double wid, double cphi, double cz, double crho); StripDetector(double centerPhi, double centerZ, double centerRho);
~StripDetector(); ~StripDetector();
const ROOT::Math::XYZPoint& GetFrontStripCoordinates(int stripch, int corner) const { return m_frontStripCoords[stripch][corner]; } 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& GetBackStripCoordinates(int stripch, int corner) const { return m_backStripCoords[stripch][corner]; }
@ -50,18 +50,10 @@ public:
StripHit GetChannelRatio(double theta, double phi); StripHit GetChannelRatio(double theta, double phi);
private: private:
bool ValidChannel(int f) { return ((f >= 0 && f < m_nStrips) ? true : false); }; bool ValidChannel(int f) { return ((f >= 0 && f < s_nStrips) ? true : false); };
bool ValidRatio(double r) { return ((r >= -1 && r <= 1) ? true : false); }; bool ValidRatio(double r) { return ((r >= -1 && r <= 1) ? true : false); };
void CalculateCorners(); 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_centerPhi; //assuming det centered above x-axis (corresponds to zero phi)
double m_centerZ; double m_centerZ;
double m_centerRho; //perpendicular radius from axis double m_centerRho; //perpendicular radius from axis
@ -77,6 +69,15 @@ private:
bool m_isSmearing; bool m_isSmearing;
//Units in meters
static constexpr double s_nStrips = 4; //Same for front and back
static constexpr double s_nCorners = 4;
static constexpr double s_totalLength = 0.075; //length of front strips
static constexpr double s_backStripLength = s_totalLength / s_nStrips;
static constexpr double s_totalWidth = 0.04; //width of back strips
static constexpr double s_frontStripWidth = s_totalWidth / s_nStrips;
static constexpr double s_deg2rad = M_PI/180.0;
}; };
#endif #endif

53
src/Detectors/SabreDetector.cpp
View File

@ -53,8 +53,7 @@
#include "SabreDetector.h" #include "SabreDetector.h"
SabreDetector::SabreDetector() : SabreDetector::SabreDetector() :
m_outerR(0.1351), m_innerR(0.0326), m_deltaPhiFlat(54.4*s_deg2rad), m_centerPhi(0.0), m_tilt(0.0), m_centerPhi(0.0), m_tilt(0.0), m_translation(0.,0.,0.), m_norm(0,0,1.0), m_detectorID(-1)
m_translation(0.,0.,0.), m_norm(0,0,1.0), m_detectorID(-1)
{ {
m_yRotation.SetAngle(m_tilt); m_yRotation.SetAngle(m_tilt);
m_zRotation.SetAngle(m_centerPhi); m_zRotation.SetAngle(m_centerPhi);
@ -73,17 +72,11 @@ m_outerR(0.1351), m_innerR(0.0326), m_deltaPhiFlat(54.4*s_deg2rad), m_centerPhi(
m_tiltWedgeCoords[i].resize(4); m_tiltWedgeCoords[i].resize(4);
} }
m_deltaRFlat = m_outerR - m_innerR;
m_deltaRFlatRing = m_deltaRFlat/s_nRings;
m_deltaPhiFlatWedge = m_deltaPhiFlat/s_nWedges;
CalculateCorners(); CalculateCorners();
} }
SabreDetector::SabreDetector(int detID, double Rin, double Rout, double deltaPhi_flat, double phiCentral, SabreDetector::SabreDetector(int detID, double phiCentral, double tiltFromVert, double zdist, double xdist, double ydist) :
double tiltFromVert, double zdist, double xdist, double ydist) : m_centerPhi(phiCentral), m_tilt(tiltFromVert), m_translation(xdist, ydist, zdist), m_norm(0,0,1.0), m_detectorID(detID)
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_yRotation.SetAngle(m_tilt); m_yRotation.SetAngle(m_tilt);
m_zRotation.SetAngle(m_centerPhi); m_zRotation.SetAngle(m_centerPhi);
@ -104,10 +97,6 @@ SabreDetector::SabreDetector(int detID, double Rin, double Rout, double deltaPhi
m_tiltWedgeCoords[i].resize(4); m_tiltWedgeCoords[i].resize(4);
} }
m_deltaRFlat = m_outerR - m_innerR;
m_deltaRFlatRing = m_deltaRFlat/s_nRings;
m_deltaPhiFlatWedge = m_deltaPhiFlat/s_nWedges;
CalculateCorners(); CalculateCorners();
} }
@ -123,23 +112,23 @@ void SabreDetector::CalculateCorners()
//Generate flat ring corner coordinates //Generate flat ring corner coordinates
for(int i=0; i<s_nRings; i++) for(int i=0; i<s_nRings; i++)
{ {
x0 = (m_innerR + m_deltaRFlatRing*(i+1))*std::cos(-m_deltaPhiFlat/2.0); x0 = (s_innerR + s_deltaR*(i+1))*std::cos(-s_deltaPhiTotal/2.0);
y0 = (m_innerR + m_deltaRFlatRing*(i+1))*std::sin(-m_deltaPhiFlat/2.0); y0 = (s_innerR + s_deltaR*(i+1))*std::sin(-s_deltaPhiTotal/2.0);
z0 = 0.0; z0 = 0.0;
m_flatRingCoords[i][0].SetXYZ(x0, y0, z0); m_flatRingCoords[i][0].SetXYZ(x0, y0, z0);
x1 = (m_innerR + m_deltaRFlatRing*(i))*std::cos(-m_deltaPhiFlat/2.0); x1 = (s_innerR + s_deltaR*(i))*std::cos(-s_deltaPhiTotal/2.0);
y1 = (m_innerR + m_deltaRFlatRing*(i))*std::sin(-m_deltaPhiFlat/2.0); y1 = (s_innerR + s_deltaR*(i))*std::sin(-s_deltaPhiTotal/2.0);
z1 = 0.0; z1 = 0.0;
m_flatRingCoords[i][1].SetXYZ(x1, y1, z1); m_flatRingCoords[i][1].SetXYZ(x1, y1, z1);
x2 = (m_innerR + m_deltaRFlatRing*(i))*std::cos(m_deltaPhiFlat/2.0); x2 = (s_innerR + s_deltaR*(i))*std::cos(s_deltaPhiTotal/2.0);
y2 = (m_innerR + m_deltaRFlatRing*(i))*std::sin(m_deltaPhiFlat/2.0); y2 = (s_innerR + s_deltaR*(i))*std::sin(s_deltaPhiTotal/2.0);
z2 = 0.0; z2 = 0.0;
m_flatRingCoords[i][2].SetXYZ(x2, y2, z2); m_flatRingCoords[i][2].SetXYZ(x2, y2, z2);
x3 = (m_innerR + m_deltaRFlatRing*(i+1))*std::cos(m_deltaPhiFlat/2.0); x3 = (s_innerR + s_deltaR*(i+1))*std::cos(s_deltaPhiTotal/2.0);
y3 = (m_innerR + m_deltaRFlatRing*(i+1))*std::sin(m_deltaPhiFlat/2.0); y3 = (s_innerR + s_deltaR*(i+1))*std::sin(s_deltaPhiTotal/2.0);
z3 = 0.0; z3 = 0.0;
m_flatRingCoords[i][3].SetXYZ(x3, y3, z3); m_flatRingCoords[i][3].SetXYZ(x3, y3, z3);
} }
@ -147,23 +136,23 @@ void SabreDetector::CalculateCorners()
//Generate flat wedge corner coordinates //Generate flat wedge corner coordinates
for(int i=0; i<s_nWedges; i++) for(int i=0; i<s_nWedges; i++)
{ {
x0 = m_outerR * std::cos(-m_deltaPhiFlat/2.0 + i*m_deltaPhiFlatWedge); x0 = s_outerR * std::cos(-s_deltaPhiTotal/2.0 + i*s_deltaPhi);
y0 = m_outerR * std::sin(-m_deltaPhiFlat/2.0 + i*m_deltaPhiFlatWedge); y0 = s_outerR * std::sin(-s_deltaPhiTotal/2.0 + i*s_deltaPhi);
z0 = 0.0; z0 = 0.0;
m_flatWedgeCoords[i][0].SetXYZ(x0, y0, z0); m_flatWedgeCoords[i][0].SetXYZ(x0, y0, z0);
x1 = m_innerR * std::cos(-m_deltaPhiFlat/2.0 + i*m_deltaPhiFlatWedge); x1 = s_innerR * std::cos(-s_deltaPhiTotal/2.0 + i*s_deltaPhi);
y1 = m_innerR * std::sin(-m_deltaPhiFlat/2.0 + i*m_deltaPhiFlatWedge); y1 = s_innerR * std::sin(-s_deltaPhiTotal/2.0 + i*s_deltaPhi);
z1 = 0.0; z1 = 0.0;
m_flatWedgeCoords[i][1].SetXYZ(x1, y1, z1); m_flatWedgeCoords[i][1].SetXYZ(x1, y1, z1);
x2 = m_innerR * std::cos(-m_deltaPhiFlat/2.0 + (i+1)*m_deltaPhiFlatWedge); x2 = s_innerR * std::cos(-s_deltaPhiTotal/2.0 + (i+1)*s_deltaPhi);
y2 = m_innerR * std::sin(-m_deltaPhiFlat/2.0 + (i+1)*m_deltaPhiFlatWedge); y2 = s_innerR * std::sin(-s_deltaPhiTotal/2.0 + (i+1)*s_deltaPhi);
z2 = 0.0; z2 = 0.0;
m_flatWedgeCoords[i][2].SetXYZ(x2, y2, z2); m_flatWedgeCoords[i][2].SetXYZ(x2, y2, z2);
x3 = m_outerR * std::cos(-m_deltaPhiFlat/2.0 + (i+1)*m_deltaPhiFlatWedge); x3 = s_outerR * std::cos(-s_deltaPhiTotal/2.0 + (i+1)*s_deltaPhi);
y3 = m_outerR * std::sin(-m_deltaPhiFlat/2.0 + (i+1)*m_deltaPhiFlatWedge); y3 = s_outerR * std::sin(-s_deltaPhiTotal/2.0 + (i+1)*s_deltaPhi);
z3 = 0.0; z3 = 0.0;
m_flatWedgeCoords[i][3].SetXYZ(x3, y3, z3); m_flatWedgeCoords[i][3].SetXYZ(x3, y3, z3);
} }
@ -322,8 +311,8 @@ ROOT::Math::XYZPoint SabreDetector::GetHitCoordinates(int ringch, int wedgech)
if(!CheckRingChannel(ringch) || !CheckWedgeChannel(wedgech)) if(!CheckRingChannel(ringch) || !CheckWedgeChannel(wedgech))
return ROOT::Math::XYZPoint(); return ROOT::Math::XYZPoint();
double r_center = m_innerR + (0.5+ringch)*m_deltaRFlatRing; double r_center = s_innerR + (0.5+ringch)*s_deltaR;
double phi_center = -m_deltaPhiFlat/2.0 + (0.5+wedgech)*m_deltaPhiFlatWedge; double phi_center = -s_deltaPhiTotal/2.0 + (0.5+wedgech)*s_deltaPhi;
double x = r_center*std::cos(phi_center); double x = r_center*std::cos(phi_center);
double y = r_center*std::sin(phi_center); double y = r_center*std::sin(phi_center);
double z = 0; double z = 0;

53
src/Detectors/SabreDetector.h
View File

@ -65,7 +65,7 @@ class SabreDetector {
public: public:
SabreDetector(); SabreDetector();
SabreDetector(int detID, double Rin, double Rout, double deltaPhi_flat, double phiCentral, double tiltFromVert, double zdist, double xdist=0, double ydist=0); SabreDetector(int detID, double phiCentral, double tiltFromVert, double zdist, double xdist=0, double ydist=0);
~SabreDetector(); ~SabreDetector();
/*Return coordinates of the corners of each ring/wedge in SABRE*/ /*Return coordinates of the corners of each ring/wedge in SABRE*/
@ -84,9 +84,6 @@ public:
int GetDetectorID() { return m_detectorID; } int GetDetectorID() { return m_detectorID; }
private: private:
void CalculateCorners(); void CalculateCorners();
/*Performs the transformation to the tilted,rotated,translated frame of the SABRE detector*/ /*Performs the transformation to the tilted,rotated,translated frame of the SABRE detector*/
@ -109,9 +106,9 @@ private:
/*Determine if a hit is within the bulk detector*/ /*Determine if a hit is within the bulk detector*/
bool IsInside(double r, double phi) bool IsInside(double r, double phi)
{ {
double phi_1 = m_deltaPhiFlat/2.0; double phi_1 = s_deltaPhiTotal/2.0;
double phi_2 = M_PI*2.0 - m_deltaPhiFlat/2.0; double phi_2 = M_PI*2.0 - s_deltaPhiTotal/2.0;
return (((r > m_innerR && r < m_outerR) || CheckPositionEqual(r, m_innerR) || CheckPositionEqual(r, m_outerR)) return (((r > s_innerR && r < s_outerR) || CheckPositionEqual(r, s_innerR) || CheckPositionEqual(r, s_outerR))
&& (phi > phi_2 || phi < phi_1 || CheckAngleEqual(phi, phi_1) || CheckAngleEqual(phi, phi_2))); && (phi > phi_2 || phi < phi_1 || CheckAngleEqual(phi, phi_1) || CheckAngleEqual(phi, phi_2)));
}; };
@ -121,57 +118,43 @@ private:
*/ */
bool IsRing(double r, int ringch) bool IsRing(double r, int ringch)
{ {
double ringtop = m_innerR + m_deltaRFlatRing*(ringch + 1); double ringtop = s_innerR + s_deltaR*(ringch + 1);
double ringbottom = m_innerR + m_deltaRFlatRing*(ringch); double ringbottom = s_innerR + s_deltaR*(ringch);
return (r>ringbottom && r<ringtop); return (r>ringbottom && r<ringtop);
}; };
inline bool IsRingTopEdge(double r, int ringch) inline bool IsRingTopEdge(double r, int ringch)
{ {
double ringtop = m_innerR + m_deltaRFlatRing*(ringch + 1); double ringtop = s_innerR + s_deltaR*(ringch + 1);
return CheckPositionEqual(r, ringtop); return CheckPositionEqual(r, ringtop);
}; };
inline bool IsRingBottomEdge(double r, int ringch) inline bool IsRingBottomEdge(double r, int ringch)
{ {
double ringbottom = m_innerR + m_deltaRFlatRing*(ringch); double ringbottom = s_innerR + s_deltaR*(ringch);
return CheckPositionEqual(r, ringbottom); return CheckPositionEqual(r, ringbottom);
}; };
inline bool IsWedge(double phi, int wedgech) inline bool IsWedge(double phi, int wedgech)
{ {
double wedgetop = -m_deltaPhiFlat/2.0 + m_deltaPhiFlatWedge*(wedgech+1); double wedgetop = -s_deltaPhiTotal/2.0 + s_deltaPhi*(wedgech+1);
double wedgebottom = -m_deltaPhiFlat/2.0 + m_deltaPhiFlatWedge*(wedgech); double wedgebottom = -s_deltaPhiTotal/2.0 + s_deltaPhi*(wedgech);
return ((phi>wedgebottom && phi<wedgetop)); return ((phi>wedgebottom && phi<wedgetop));
}; };
inline bool IsWedgeTopEdge(double phi, int wedgech) inline bool IsWedgeTopEdge(double phi, int wedgech)
{ {
double wedgetop = -m_deltaPhiFlat/2.0 + m_deltaPhiFlatWedge*(wedgech+1); double wedgetop = -s_deltaPhiTotal/2.0 + s_deltaPhi*(wedgech+1);
return CheckAngleEqual(phi, wedgetop); return CheckAngleEqual(phi, wedgetop);
} }
inline bool IsWedgeBottomEdge(double phi, int wedgech) inline bool IsWedgeBottomEdge(double phi, int wedgech)
{ {
double wedgebottom = -m_deltaPhiFlat/2.0 + m_deltaPhiFlatWedge*(wedgech); double wedgebottom = -s_deltaPhiTotal/2.0 + s_deltaPhi*(wedgech);
return CheckAngleEqual(phi, wedgebottom); 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*/ /*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_centerPhi;
double m_tilt; double m_tilt;
ROOT::Math::Translation3D m_translation; ROOT::Math::Translation3D m_translation;
@ -183,6 +166,18 @@ private:
std::vector<std::vector<ROOT::Math::XYZPoint>> m_flatRingCoords, m_flatWedgeCoords; std::vector<std::vector<ROOT::Math::XYZPoint>> m_flatRingCoords, m_flatWedgeCoords;
std::vector<std::vector<ROOT::Math::XYZPoint>> m_tiltRingCoords, m_tiltWedgeCoords; std::vector<std::vector<ROOT::Math::XYZPoint>> m_tiltRingCoords, m_tiltWedgeCoords;
/*Class constants*/
static constexpr double s_deg2rad = M_PI/180.0;
static constexpr int s_nRings = 16;
static constexpr int s_nWedges = 8;
static constexpr double s_outerR = 0.1351;
static constexpr double s_innerR = 0.0326;
static constexpr double s_deltaR = (s_outerR - s_innerR) / s_nRings;
static constexpr double s_deltaPhiTotal = 54.4 * s_deg2rad;
static constexpr double s_deltaPhi = (s_deltaPhiTotal / s_nWedges);
/*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
}; };

2
src/Detectors/SabreEfficiency.cpp
View File

@ -13,7 +13,7 @@ SabreEfficiency::SabreEfficiency() :
{ {
for(int i=0; i<s_nDets; i++) 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); m_detectors.emplace_back(i, s_centerPhiList[i]*s_deg2rad, s_tilt*s_deg2rad, s_zOffset);
//Only 0, 1, 4 valid in degrader land //Only 0, 1, 4 valid in degrader land
m_degradedDetectors[0] = true; m_degradedDetectors[0] = true;
m_degradedDetectors[1] = true; m_degradedDetectors[1] = true;

3
src/Detectors/SabreEfficiency.h
View File

@ -32,11 +32,8 @@ private:
bool m_degradedDetectors[5]; bool m_degradedDetectors[5];
//Sabre constants //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_tilt = 40.0;
static constexpr double s_zOffset = -0.1245; 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 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_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_deg2rad = M_PI/180.0;