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

modified: Armory/ClassPW.h

Browse Source
This commit is contained in:
Vignesh Sitaraman 2026-03-09 14:24:44 -04:00
parent 9af7e70f26
commit 65af941b39
1 changed files with 69 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

94
Armory/ClassPW.h
View File

@ -62,6 +62,9 @@ public:
double GetTrackPhi() const { return trackVec.Phi(); } double GetTrackPhi() const { return trackVec.Phi(); }
double GetZ0(); double GetZ0();
Coord Crossover[24][24][2];
inline std::tuple<std::pair<TVector3, TVector3>, double, double, double> GetPseudoWire(const std::vector<std::tuple<int,double,double>>& cluster, std::string type); inline std::tuple<std::pair<TVector3, TVector3>, double, double, double> GetPseudoWire(const std::vector<std::tuple<int,double,double>>& cluster, std::string type);
inline std::tuple<TVector3,double,double,double,double,double,double,double> inline std::tuple<TVector3,double,double,double,double,double,double,double>
@ -115,7 +118,8 @@ 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 zLen = 348.6; // mm // const float zLen = 348.6; // mm
const float zLen = 174.3*2; // mm
const float radiusA = 37; const float radiusA = 37;
const float radiusC = 43; const float radiusC = 43;
@ -149,35 +153,71 @@ 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 double k = TMath::TwoPi()/24.; //48 solder thru holes, wires in every other one
// anode wire shift is right-hand. double offset_a1 = -6*k-3*k;
// cathode wire shift is left-hand. double offset_c1 = -4*k -2*k - TMath::TwoPi()/48; //correct for a half-turn
#include "../scratch/testing.h"
double offset_a2 = offset_a1+wireShift*k;
double offset_c2 = offset_c1-wireShift*k;
for (int i = 0; i < nWire; i++) for (int i = 0; i < nWire; i++)
{ {
// Anode rotate right-hand // Anode rotate right-hand coming in towards +z riding with the beam. In this frame, +x is to the right, and +y down
p1.first.SetXYZ(radiusA * TMath::Cos(TMath::TwoPi() / nWire * (i) + TMath::PiOver2()), //updated Feb 2026, Sudarsan B. Photographs indicate that anode wires twist right handed, as one moves from -z to +z with the convention above
radiusA * TMath::Sin(TMath::TwoPi() / nWire * (i) + TMath::PiOver2()), //wire indices increase leftward as one moves to +z (hence -k factor), but wires themselves twist rightward - as indicated by offset_a2 being more +ve w.r.t offset_a1
zLen / 2); //'First' is -z locus, 'second' is +z locus
p1.second.SetXYZ(radiusA * TMath::Cos(TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()), p1.first.SetXYZ(radiusA * TMath::Cos(-k*i + offset_a1),
radiusA * TMath::Sin(TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()), radiusA * TMath::Sin(-k*i + offset_a1),
-zLen / 2); -zLen / 2);
An.push_back(p1); p1.second.SetXYZ(radiusA * TMath::Cos(-k*i + offset_a2),
radiusA * TMath::Sin(-k*i + offset_a2),
+zLen / 2);
// Cathod rotate left-hand with the 3 wire offset accounted for (+1 from the calculated offset from the PC coincidence spectrum) // Cathodes twist left-hand as indicated by offset_c2 being more negative than offset_c1, under the same system, while wires increase rightward (hence +k factor)
q1.first.SetXYZ(radiusC * TMath::Cos(TMath::TwoPi() / nWire * (i + wireShift + 1) + TMath::PiOver2()), q1.first.SetXYZ(radiusC * TMath::Cos(k*i + offset_c1),
radiusC * TMath::Sin(TMath::TwoPi() / nWire * (i + wireShift + 1) + TMath::PiOver2()), radiusC * TMath::Sin(k*i + offset_c1),
zLen / 2);
q1.second.SetXYZ(radiusC * TMath::Cos(TMath::TwoPi() / nWire * (i + 1) + TMath::PiOver2()),
radiusC * TMath::Sin(TMath::TwoPi() / nWire * (i + 1) + TMath::PiOver2()),
-zLen / 2); -zLen / 2);
q1.second.SetXYZ(radiusC * TMath::Cos(k*i + offset_c2),
radiusC * TMath::Sin(k*i + offset_c2),
zLen / 2);
An.push_back(p1);
Ca.push_back(q1); Ca.push_back(q1);
} }
// correcting for the fact that the order of the cathode wires is reversed
std::reverse(Ca.begin(), Ca.end());
// adjusting for the 3 wire offset, the rbegin and rend are used as the rotation of the wires is done in the opposite direction i.e. 1,2,3 -> 3,1,2 // Calculate Crossover Geometry ONCE
// NOT NECESSARY ANY MORE, HAS BEEN IMCORPORATED INTO THE WIREOFFSET IN THE BEGINNING TVector3 a, c, diff;
// std::rotate(Ca.rbegin(), Ca.rbegin() + 4, Ca.rend()); double a2, ac, c2, adiff, cdiff, denom, alpha;
for (size_t i = 0; i < An.size(); i++)
{
//a = An[i].first - An[i].second;
a = An[i].second - An[i].first;
for (size_t j = 0; j < Ca.size(); j++)
{
c = Ca[j].second- Ca[j].first;
diff = An[i].second - Ca[j].second;
a2 = a.Dot(a);
c2 = c.Dot(c);
ac = a.Dot(c);
adiff = a.Dot(diff);
cdiff = c.Dot(diff);
denom = a2 * c2 - ac * ac;
alpha = (ac * cdiff - c2 * adiff) / denom;
Crossover[i][j][0].x = An[i].second.X() + alpha * a.X();
Crossover[i][j][0].y = An[i].second.Y() + alpha * a.Y();
Crossover[i][j][0].z = An[i].second.Z() + alpha * a.Z();
if (Crossover[i][j][0].z < -190 || Crossover[i][j][0].z > 190 || (i+j)%24 == 12) {
Crossover[i][j][0].z = 9999999;
}
Crossover[i][j][1].x = alpha;
Crossover[i][j][1].y = 0;
}
}
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));
@ -276,8 +316,8 @@ PW::GetPseudoWire(const std::vector<std::tuple<int,double,double>>& cluster, std
avgvec.second = avgvec.second*(1.0/sumEnergy); avgvec.second = avgvec.second*(1.0/sumEnergy);
double phi1 = avgvec.first.Phi(); double phi1 = avgvec.first.Phi();
double phi2 = avgvec.second.Phi(); double phi2 = avgvec.second.Phi();
avgvec.first.SetXYZ(radiusA*TMath::Cos(phi1), radiusA*TMath::Sin(phi1), zLen/2); avgvec.first.SetXYZ(radiusA*TMath::Cos(phi1), radiusA*TMath::Sin(phi1), -zLen/2);
avgvec.second.SetXYZ(radiusA*TMath::Cos(phi2), radiusA*TMath::Sin(phi2), -zLen/2); avgvec.second.SetXYZ(radiusA*TMath::Cos(phi2), radiusA*TMath::Sin(phi2), zLen/2);
/*if(cluster.size()>1) { /*if(cluster.size()>1) {
std::cout << "\t\t avg1(r,phi,z):" << avgvec.first.Perp() << " " << avgvec.first.Phi()*180/M_PI << " " << avgvec.first.Z() << std::endl; std::cout << "\t\t avg1(r,phi,z):" << avgvec.first.Perp() << " " << avgvec.first.Phi()*180/M_PI << " " << avgvec.first.Z() << std::endl;
std::cout << "\t\t avg2(r,phi,z):" << avgvec.second.Perp() << " " << avgvec.second.Phi()*180/M_PI << " " << avgvec.second.Z() << std::endl; std::cout << "\t\t avg2(r,phi,z):" << avgvec.second.Perp() << " " << avgvec.second.Phi()*180/M_PI << " " << avgvec.second.Z() << std::endl;
@ -296,8 +336,8 @@ PW::GetPseudoWire(const std::vector<std::tuple<int,double,double>>& cluster, std
avgvec.second = avgvec.second*(1.0/sumEnergy); avgvec.second = avgvec.second*(1.0/sumEnergy);
double phi1 = avgvec.first.Phi(); double phi1 = avgvec.first.Phi();
double phi2 = avgvec.second.Phi(); double phi2 = avgvec.second.Phi();
avgvec.first.SetXYZ(radiusC*TMath::Cos(phi1), radiusC*TMath::Sin(phi1), zLen/2); avgvec.first.SetXYZ(radiusC*TMath::Cos(phi1), radiusC*TMath::Sin(phi1), -zLen/2);
avgvec.second.SetXYZ(radiusC*TMath::Cos(phi2), radiusC*TMath::Sin(phi2), -zLen/2); avgvec.second.SetXYZ(radiusC*TMath::Cos(phi2), radiusC*TMath::Sin(phi2), zLen/2);
} }
return std::tuple(avgvec, sumEnergy, maxEnergy, tsMaxEnergy); return std::tuple(avgvec, sumEnergy, maxEnergy, tsMaxEnergy);
} }
@ -317,7 +357,7 @@ inline std::tuple<TVector3,double,double,double,double,double,double,double> PW:
if(apSumE && cpSumE) { if(apSumE && cpSumE) {
a = apwire.first - apwire.second; a = apwire.first - apwire.second;
c = cpwire.first - cpwire.second; c = cpwire.first - cpwire.second;
diff = apwire.first - cpwire.first; diff = apwire.first- cpwire.first;
a2 = a.Dot(a); a2 = a.Dot(a);
c2 = c.Dot(c); c2 = c.Dot(c);
ac = a.Dot(c); ac = a.Dot(c);