rtang/ANASEN_analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

modified: .vscode/settings.json

Browse Source 
modified:   Armory/ClassPW.h
This commit is contained in:
vs19g 2025-01-16 09:56:39 -05:00
parent d8b11bdcf3
commit 9bf83f8028
2 changed files with 162 additions and 110 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
.vscode/settings.json vendored
View File

@ -90,6 +90,17 @@
"processRun.C": "cpp", "processRun.C": "cpp",
"TrackRecon.C": "cpp", "TrackRecon.C": "cpp",
"processRuns.C": "cpp", "processRuns.C": "cpp",
"Analysis.C": "cpp" "Analysis.C": "cpp",
} "datastructs.h": "c",
"ANASENPlotEdit.C": "cpp",
"GetMean_Q3_new.C": "cpp",
"AlphaCal_new.C": "cpp",
"f1.C": "cpp",
"GeoCal_Maria_new.C": "cpp",
"PCPulser_All_new.C": "cpp",
"PosCal_2.C": "cpp",
"AutoFit.C": "cpp",
"Fitting.C": "cpp"
},
"github-enterprise.uri": "https://fsunuc.physics.fsu.edu"
} }

241
Armory/ClassPW.h
View File

@ -6,14 +6,16 @@
#include <TVector3.h> #include <TVector3.h>
#include <TRandom.h> #include <TRandom.h>
struct PWHitInfo{ struct PWHitInfo
std::pair<short, short> nearestWire; // anode, cathode {
std::pair<short, short> nearestWire; // anode, cathode
std::pair<double, double> nearestDist; // anode, cathode std::pair<double, double> nearestDist; // anode, cathode
std::pair<short, short> nextNearestWire; // anode, cathode std::pair<short, short> nextNearestWire; // anode, cathode
std::pair<double, double> nextNearestDist; // anode, cathode std::pair<double, double> nextNearestDist; // anode, cathode
void Clear(){ void Clear()
{
nearestWire.first = -1; nearestWire.first = -1;
nearestWire.second = -1; nearestWire.second = -1;
nearestDist.first = 999999999; nearestDist.first = 999999999;
@ -25,40 +27,52 @@ struct PWHitInfo{
} }
}; };
//!######################################################## struct Coord
class PW{ // proportional wire {
float x, y, z;
Coord(const TVector3 &vec)
{
x = vec.X(); // TVector3's X() returns the x-coordinate
y = vec.Y(); // TVector3's Y() returns the y-coordinate
z = vec.Z(); // TVector3's Z() returns the z-coordinate
}
};
//! ########################################################
class PW
{ // proportional wire
public: public:
PW(){ ClearHitInfo();}; PW() { ClearHitInfo(); };
~PW(){}; ~PW() {};
PWHitInfo GetHitInfo() const {return hitInfo;} PWHitInfo GetHitInfo() const { return hitInfo; }
std::pair<short, short> GetNearestID() const {return hitInfo.nearestWire;} std::pair<short, short> GetNearestID() const { return hitInfo.nearestWire; }
std::pair<double, double> GetNearestDistance() const {return hitInfo.nearestDist;} std::pair<double, double> GetNearestDistance() const { return hitInfo.nearestDist; }
std::pair<short, short> Get2ndNearestID() const {return hitInfo.nextNearestWire;} std::pair<short, short> Get2ndNearestID() const { return hitInfo.nextNearestWire; }
std::pair<double, double> Get2ndNearestDistance() const {return hitInfo.nextNearestDist;} std::pair<double, double> Get2ndNearestDistance() const { return hitInfo.nextNearestDist; }
TVector3 GetTrackPos() const {return trackPos;} TVector3 GetTrackPos() const { return trackPos; }
TVector3 GetTrackVec() const {return trackVec;} TVector3 GetTrackVec() const { return trackVec; }
double GetTrackTheta() const {return trackVec.Theta();} double GetTrackTheta() const { return trackVec.Theta(); }
double GetTrackPhi() const {return trackVec.Phi();} double GetTrackPhi() const { return trackVec.Phi(); }
double GetZ0(); double GetZ0();
int GetNumWire() const {return nWire;} int GetNumWire() const { return nWire; }
double GetDeltaAngle() const {return dAngle;} double GetDeltaAngle() const { return dAngle; }
double GetAnodeLength() const {return anodeLength;} double GetAnodeLength() const { return anodeLength; }
double GetCathodeLength() const {return cathodeLength;} double GetCathodeLength() const { return cathodeLength; }
TVector3 GetAnodeDn(short id) const {return An[id].first;} TVector3 GetAnodeDn(short id) const { return An[id].first; }
TVector3 GetAnodeUp(short id) const {return An[id].second;} TVector3 GetAnodeUp(short id) const { return An[id].second; }
TVector3 GetCathodeDn(short id) const {return Ca[id].first;} TVector3 GetCathodeDn(short id) const { return Ca[id].first; }
TVector3 GetCathodeUp(short id) const {return Ca[id].second;} TVector3 GetCathodeUp(short id) const { return Ca[id].second; }
TVector3 GetAnodneMid(short id) const {return (An[id].first + An[id].second) * 0.5; } TVector3 GetAnodneMid(short id) const { return (An[id].first + An[id].second) * 0.5; }
double GetAnodeTheta(short id) const {return (An[id].first - An[id].second).Theta();} double GetAnodeTheta(short id) const { return (An[id].first - An[id].second).Theta(); }
double GetAnodePhi(short id) const {return (An[id].first - An[id].second).Phi();} double GetAnodePhi(short id) const { return (An[id].first - An[id].second).Phi(); }
TVector3 GetCathodneMid(short id) const {return (Ca[id].first + Ca[id].second) * 0.5; } TVector3 GetCathodneMid(short id) const { return (Ca[id].first + Ca[id].second) * 0.5; }
double GetCathodeTheta(short id) const {return (Ca[id].first - Ca[id].second).Theta();} double GetCathodeTheta(short id) const { return (Ca[id].first - Ca[id].second).Theta(); }
double GetCathodePhi(short id) const {return (Ca[id].first - Ca[id].second).Phi();} double GetCathodePhi(short id) const { return (Ca[id].first - Ca[id].second).Phi(); }
void ClearHitInfo(); void ClearHitInfo();
void ConstructGeo(); void ConstructGeo();
@ -66,20 +80,20 @@ public:
void CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose = false); void CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose = false);
void CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA = 0, double sigmaC = 0, bool verbose = false); void CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA = 0, double sigmaC = 0, bool verbose = false);
void Print(){ void Print()
{
printf(" The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nearestWire.first, printf(" The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nearestWire.first,
hitInfo.nearestDist.first, hitInfo.nearestDist.first,
hitInfo.nearestWire.second, hitInfo.nearestWire.second,
hitInfo.nearestDist.second); hitInfo.nearestDist.second);
printf(" The 2nd nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nextNearestWire.first, printf(" The 2nd nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nextNearestWire.first,
hitInfo.nextNearestDist.first, hitInfo.nextNearestDist.first,
hitInfo.nextNearestWire.second, hitInfo.nextNearestWire.second,
hitInfo.nextNearestDist.second); hitInfo.nextNearestDist.second);
} }
private: private:
PWHitInfo hitInfo; PWHitInfo hitInfo;
TVector3 trackPos; TVector3 trackPos;
@ -87,7 +101,7 @@ private:
const int nWire = 24; const int nWire = 24;
const int wireShift = 3; const int wireShift = 3;
const float zLen = 380; //mm const float zLen = 380; // mm
const float radiusA = 37; const float radiusA = 37;
const float radiusC = 43; const float radiusC = 43;
@ -95,23 +109,25 @@ private:
double anodeLength; double anodeLength;
double cathodeLength; double cathodeLength;
std::vector<std::pair<TVector3,TVector3>> An; // the anode wire position vector in space std::vector<std::pair<TVector3, TVector3>> An; // the anode wire position vector in space
std::vector<std::pair<TVector3,TVector3>> Ca; // the cathode wire position vector in space std::vector<std::pair<TVector3, TVector3>> Ca; // the cathode wire position vector in space
double Distance(TVector3 a1, TVector3 a2, TVector3 b1, TVector3 b2){ double Distance(TVector3 a1, TVector3 a2, TVector3 b1, TVector3 b2)
{
TVector3 na = a1 - a2; TVector3 na = a1 - a2;
TVector3 nb = b1 - b2; TVector3 nb = b1 - b2;
TVector3 nd = (na.Cross(nb)).Unit(); TVector3 nd = (na.Cross(nb)).Unit();
return TMath::Abs(nd.Dot(a1-b2)); return TMath::Abs(nd.Dot(a1 - b2));
} }
}; };
inline void PW::ClearHitInfo(){ inline void PW::ClearHitInfo()
{
hitInfo.Clear(); hitInfo.Clear();
} }
inline void PW::ConstructGeo(){ inline void PW::ConstructGeo()
{
An.clear(); An.clear();
Ca.clear(); Ca.clear();
@ -119,108 +135,132 @@ inline void PW::ConstructGeo(){
std::pair<TVector3, TVector3> p1; // anode std::pair<TVector3, TVector3> p1; // anode
std::pair<TVector3, TVector3> q1; // cathode std::pair<TVector3, TVector3> q1; // cathode
//anode and cathode start at pos-Y axis and count in right-Hand // anode and cathode start at pos-Y axis and count in right-Hand
//anode wire shift is right-hand. // anode wire shift is right-hand.
//cathode wire shift is left-hand. // cathode wire shift is left-hand.
for(int i = 0; i < nWire; i++ ){ for (int i = 0; i < nWire; i++)
{
// Anode rotate right-hand // Anode rotate right-hand
p1.first.SetXYZ( radiusA * TMath::Cos( TMath::TwoPi() / nWire * (i) + TMath::PiOver2()), p1.first.SetXYZ(radiusA * TMath::Cos(TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
radiusA * TMath::Sin( TMath::TwoPi() / nWire * (i) + TMath::PiOver2()), radiusA * TMath::Sin(TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
zLen/2); zLen / 2);
p1.second.SetXYZ( radiusA * TMath::Cos( TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()), p1.second.SetXYZ(radiusA * TMath::Cos(TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()),
radiusA * TMath::Sin( TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()), radiusA * TMath::Sin(TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()),
-zLen/2); -zLen / 2);
An.push_back(p1); An.push_back(p1);
// Cathod rotate left-hand // Cathod rotate left-hand
q1.first.SetXYZ( radiusC * TMath::Cos( TMath::TwoPi() / nWire * (i) + TMath::PiOver2()), q1.first.SetXYZ(radiusC * TMath::Cos(TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
radiusC * TMath::Sin( TMath::TwoPi() / nWire * (i) + TMath::PiOver2()), radiusC * TMath::Sin(TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
zLen/2); zLen / 2);
q1.second.SetXYZ( radiusC * TMath::Cos( TMath::TwoPi() / nWire * (i - wireShift) + TMath::PiOver2()), q1.second.SetXYZ(radiusC * TMath::Cos(TMath::TwoPi() / nWire * (i - wireShift) + TMath::PiOver2()),
radiusC * TMath::Sin( TMath::TwoPi() / nWire * (i - wireShift) + TMath::PiOver2()), radiusC * TMath::Sin(TMath::TwoPi() / nWire * (i - wireShift) + TMath::PiOver2()),
-zLen/2); -zLen / 2);
Ca.push_back(q1); Ca.push_back(q1);
} }
dAngle = wireShift * TMath::TwoPi() / nWire; dAngle = wireShift * TMath::TwoPi() / nWire;
anodeLength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusA * TMath::Sin(dAngle/2),2) ); anodeLength = TMath::Sqrt(zLen * zLen + TMath::Power(2 * radiusA * TMath::Sin(dAngle / 2), 2));
cathodeLength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusC * TMath::Sin(dAngle/2),2) ); cathodeLength = TMath::Sqrt(zLen * zLen + TMath::Power(2 * radiusC * TMath::Sin(dAngle / 2), 2));
} }
inline void PW::FindWireID(TVector3 pos, TVector3 direction, bool verbose ){ inline void PW::FindWireID(TVector3 pos, TVector3 direction, bool verbose)
{
hitInfo.Clear(); hitInfo.Clear();
double phi = direction.Phi(); double phi = direction.Phi();
for( int i = 0; i < nWire; i++){ for (int i = 0; i < nWire; i++)
{
double disA = 99999999; double disA = 99999999;
double phiS = An[i].first.Phi() - TMath::PiOver4(); double phiS = An[i].first.Phi() - TMath::PiOver4();
double phiL = An[i].second.Phi() + TMath::PiOver4(); double phiL = An[i].second.Phi() + TMath::PiOver4();
// printf("A%2d: %f %f | %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg(), phi * TMath::RadToDeg()); // printf("A%2d: %f %f | %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg(), phi * TMath::RadToDeg());
if( phi > 0 && phiS > phiL ) phiL = phiL + TMath::TwoPi(); if (phi > 0 && phiS > phiL)
if( phi < 0 && phiS > phiL ) phiS = phiS - TMath::TwoPi(); phiL = phiL + TMath::TwoPi();
if (phi < 0 && phiS > phiL)
phiS = phiS - TMath::TwoPi();
if( phiS < phi && phi < phiL) { if (phiS < phi && phi < phiL)
disA = Distance( pos, pos + direction, An[i].first, An[i].second); {
if( disA < hitInfo.nearestDist.first ){ disA = Distance(pos, pos + direction, An[i].first, An[i].second);
if (disA < hitInfo.nearestDist.first)
{
hitInfo.nearestDist.first = disA; hitInfo.nearestDist.first = disA;
hitInfo.nearestWire.first = i; hitInfo.nearestWire.first = i;
} }
} }
double disC = 99999999; double disC = 99999999;
phiS = Ca[i].second.Phi()- TMath::PiOver4(); phiS = Ca[i].second.Phi() - TMath::PiOver4();
phiL = Ca[i].first.Phi() + TMath::PiOver4(); phiL = Ca[i].first.Phi() + TMath::PiOver4();
// printf("C%2d: %f %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg()); // printf("C%2d: %f %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg());
if( phi > 0 && phiS > phiL ) phiL = phiL + TMath::TwoPi(); if (phi > 0 && phiS > phiL)
if( phi < 0 && phiS > phiL ) phiS = phiS - TMath::TwoPi(); phiL = phiL + TMath::TwoPi();
if (phi < 0 && phiS > phiL)
phiS = phiS - TMath::TwoPi();
if(phiS < phi && phi < phiL) { if (phiS < phi && phi < phiL)
disC = Distance( pos, pos + direction, Ca[i].first, Ca[i].second); {
if( disC < hitInfo.nearestDist.second ){ disC = Distance(pos, pos + direction, Ca[i].first, Ca[i].second);
if (disC < hitInfo.nearestDist.second)
{
hitInfo.nearestDist.second = disC; hitInfo.nearestDist.second = disC;
hitInfo.nearestWire.second = i; hitInfo.nearestWire.second = i;
} }
} }
if(verbose) printf(" %2d | %8.2f, %8.2f\n", i, disA, disC); if (verbose)
printf(" %2d | %8.2f, %8.2f\n", i, disA, disC);
} }
//==== find the 2nd nearest wire //==== find the 2nd nearest wire
short anode1 = hitInfo.nearestWire.first; short anode1 = hitInfo.nearestWire.first;
short aaa1 = anode1 - 1; if( aaa1 < 0 ) aaa1 += nWire; short aaa1 = anode1 - 1;
if (aaa1 < 0)
aaa1 += nWire;
short aaa2 = (anode1 + 1) % nWire; short aaa2 = (anode1 + 1) % nWire;
double haha1 = Distance( pos, pos + direction, An[aaa1].first, An[aaa1].second); double haha1 = Distance(pos, pos + direction, An[aaa1].first, An[aaa1].second);
double haha2 = Distance( pos, pos + direction, An[aaa2].first, An[aaa2].second); double haha2 = Distance(pos, pos + direction, An[aaa2].first, An[aaa2].second);
if( haha1 < haha2){ if (haha1 < haha2)
{
hitInfo.nextNearestWire.first = aaa1; hitInfo.nextNearestWire.first = aaa1;
hitInfo.nextNearestDist.first = haha1; hitInfo.nextNearestDist.first = haha1;
}else{ }
else
{
hitInfo.nextNearestWire.first = aaa2; hitInfo.nextNearestWire.first = aaa2;
hitInfo.nextNearestDist.first = haha2; hitInfo.nextNearestDist.first = haha2;
} }
short cathode1 = hitInfo.nearestWire.second; short cathode1 = hitInfo.nearestWire.second;
short ccc1 = cathode1 - 1; if( ccc1 < 0 ) ccc1 += nWire; short ccc1 = cathode1 - 1;
if (ccc1 < 0)
ccc1 += nWire;
short ccc2 = (cathode1 + 1) % nWire; short ccc2 = (cathode1 + 1) % nWire;
haha1 = Distance( pos, pos + direction, Ca[ccc1].first, Ca[ccc1].second); haha1 = Distance(pos, pos + direction, Ca[ccc1].first, Ca[ccc1].second);
haha2 = Distance( pos, pos + direction, Ca[ccc2].first, Ca[ccc2].second); haha2 = Distance(pos, pos + direction, Ca[ccc2].first, Ca[ccc2].second);
if( haha1 < haha2){ if (haha1 < haha2)
{
hitInfo.nextNearestWire.second = ccc1; hitInfo.nextNearestWire.second = ccc1;
hitInfo.nextNearestDist.second = haha1; hitInfo.nextNearestDist.second = haha1;
}else{ }
else
{
hitInfo.nextNearestWire.second = ccc2; hitInfo.nextNearestWire.second = ccc2;
hitInfo.nextNearestDist.second = haha2; hitInfo.nextNearestDist.second = haha2;
} }
if( verbose ) Print(); if (verbose)
Print();
} }
inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose){ inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose)
{
trackPos = sx3Pos; trackPos = sx3Pos;
@ -230,11 +270,12 @@ inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbo
// if the handiness of anode and cathode revered, it should be n2 cross n1 // if the handiness of anode and cathode revered, it should be n2 cross n1
trackVec = (n2.Cross(n1)).Unit(); trackVec = (n2.Cross(n1)).Unit();
if( verbose ) printf("Theta, Phi = %f, %f \n", trackVec.Theta() *TMath::RadToDeg(), trackVec.Phi()*TMath::RadToDeg()); if (verbose)
printf("Theta, Phi = %f, %f \n", trackVec.Theta() * TMath::RadToDeg(), trackVec.Phi() * TMath::RadToDeg());
} }
inline void PW::CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA, double sigmaC, bool verbose){ inline void PW::CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA, double sigmaC, bool verbose)
{
trackPos = sx3Pos; trackPos = sx3Pos;
@ -266,19 +307,19 @@ inline void PW::CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA, dou
// if the handiness of anode and cathode revered, it should be n2 cross n1 // if the handiness of anode and cathode revered, it should be n2 cross n1
trackVec = (n2.Cross(n1)).Unit(); trackVec = (n2.Cross(n1)).Unit();
if( verbose ) printf("Theta, Phi = %f, %f \n", trackVec.Theta() *TMath::RadToDeg(), trackVec.Phi()*TMath::RadToDeg()); if (verbose)
printf("Theta, Phi = %f, %f \n", trackVec.Theta() * TMath::RadToDeg(), trackVec.Phi() * TMath::RadToDeg());
} }
inline double PW::GetZ0(){ inline double PW::GetZ0()
{
double x = trackPos.X(); double x = trackPos.X();
double y = trackPos.Y(); double y = trackPos.Y();
double rho = TMath::Sqrt(x*x + y*y); double rho = TMath::Sqrt(x * x + y * y);
double theta = trackVec.Theta(); double theta = trackVec.Theta();
return trackPos.Z() - rho / TMath::Tan(theta); return trackPos.Z() - rho / TMath::Tan(theta);
} }
#endif #endif