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

new sx3 gainmatch for 17F data's first part

Browse Source
This commit is contained in:
Sudarsan Balakrishnan 2026-03-24 15:55:40 -04:00
parent 623e72a197
commit 103ddb3958
40 changed files with 6528 additions and 1076 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
Analyzer.C
View File

@ -856,22 +856,21 @@ Bool_t Analyzer::Process(Long64_t entry)
// // hZProj->Draw(); // // hZProj->Draw();
// hanVScatsum->Draw("colz"); // hanVScatsum->Draw("colz");
// // TFile *outRoot = new TFile("Histograms.root", "RECREATE"); TFile *outRoot = new TFile("Histograms.root", "RECREATE");
// // if (!outRoot->IsOpen()) if (!outRoot->IsOpen())
// // { {
// // std::cerr << "Error opening file for writing!" << std::endl; std::cerr << "Error opening file for writing!" << std::endl;
// // return; return;
// // } }
// // // Loop through histograms and write them to the ROOT file // // // Loop through histograms and write them to the ROOT file
// // for (int i = 0; i < 48; i++) for (int i = 0; i < 48; i++)
// // { {
// // if (hPC_E[i] != nullptr) if (hPC_E[i] != nullptr)
// // { {
// // hPC_E[i]->Write(); // Write histogram to file hPC_E[i]->Write(); // Write histogram to file
// // } }
// // } }
outRoot->Close();
// // outRoot->Close();
} }

2
Armory/ANASEN_model.C
View File

@ -42,7 +42,7 @@ void ANASEN_model(int anodeID1 = -1, int anodeID2 = -1, int cathodeID1 = -1, int
//--- making ANASEN //--- making ANASEN
const int nWire = 24; const int nWire = 24;
const int wireShift = 3; const int wireShift = 3;
const int zLen = 350; //mm const int zLen = 300; //mm
const int radiusA = 38; const int radiusA = 38;
const int radiusC = 43; const int radiusC = 43;

95
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,72 @@ 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
std::cerr << "Here!" << std::endl;
#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 +317,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 +337,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 +358,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);

491
Armory/ClassPW.h.Questionable.Feb2026 Normal file
View File

@ -0,0 +1,491 @@
#ifndef ClassPW_h
#define ClassPW_h
#include <cstdio>
#include <iostream>
#include <TMath.h>
#include <TVector3.h>
#include <TRandom.h>
struct PWHitInfo
{
std::pair<short, short> nearestWire; // anode, cathode
std::pair<double, double> nearestDist; // anode, cathode
std::pair<short, short> nextNearestWire; // anode, cathode
std::pair<double, double> nextNearestDist; // anode, cathode
void Clear()
{
nearestWire.first = -1;
nearestWire.second = -1;
nearestDist.first = 999999999;
nearestDist.second = 999999999;
nextNearestWire.first = -1;
nextNearestWire.second = -1;
nextNearestDist.first = 999999999;
nextNearestDist.second = 999999999;
}
};
struct Coord
{
float x, y, z;
Coord() : x(0), y(0), z(0) {}
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:
PW() { ClearHitInfo(); };
~PW() {};
PWHitInfo GetHitInfo() const { return hitInfo; }
std::pair<short, short> GetNearestID() const { return hitInfo.nearestWire; }
std::pair<double, double> GetNearestDistance() const { return hitInfo.nearestDist; }
std::pair<short, short> Get2ndNearestID() const { return hitInfo.nextNearestWire; }
std::pair<double, double> Get2ndNearestDistance() const { return hitInfo.nextNearestDist; }
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
TVector3 GetTrackPos() const { return trackPos; }
TVector3 GetTrackVec() const { return trackVec; }
double GetTrackTheta() const { return trackVec.Theta(); }
double GetTrackPhi() const { return trackVec.Phi(); }
double GetZ0();
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>
FindCrossoverProperties(const std::vector<std::tuple<int,double,double>>& a_cluster, const std::vector<std::tuple<int,double,double>>& c_cluster);
inline std::vector<std::vector<std::tuple<int,double,double>>>
Make_Clusters(std::unordered_map<int,std::tuple<int,double,double>> wireEvents);
int GetNumWire() const { return nWire; }
double GetDeltaAngle() const { return dAngle; }
double GetAnodeLength() const { return anodeLength; }
double GetCathodeLength() const { return cathodeLength; }
TVector3 GetAnodeDn(short id) const { return An[id].first; }
TVector3 GetAnodeUp(short id) const { return An[id].second; }
TVector3 GetCathodeDn(short id) const { return Ca[id].first; }
TVector3 GetCathodeUp(short id) const { return Ca[id].second; }
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 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; }
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(); }
void ClearHitInfo();
void ConstructGeo();
void FindWireID(TVector3 pos, TVector3 direction, bool verbose = false);
void CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose = false);
void CalTrack2(TVector3 sx3Pos, TVector3 anodeInt, bool verbose = false);
void Print()
{
printf(" The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nearestWire.first,
hitInfo.nearestDist.first,
hitInfo.nearestWire.second,
hitInfo.nearestDist.second);
printf(" The 2nd nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nextNearestWire.first,
hitInfo.nextNearestDist.first,
hitInfo.nextNearestWire.second,
hitInfo.nextNearestDist.second);
}
private:
PWHitInfo hitInfo;
TVector3 trackPos;
TVector3 trackVec;
const int nWire = 24;
const int wireShift = 3;
//const float zLen = 380; // mm
const float zLen = 348.6; // mm
const float radiusA = 37;
const float radiusC = 43;
double dAngle;
double anodeLength;
double cathodeLength;
// 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
double Distance(TVector3 a1, TVector3 a2, TVector3 b1, TVector3 b2)
{
TVector3 na = a1 - a2;
TVector3 nb = b1 - b2;
TVector3 nd = (na.Cross(nb)).Unit();
return TMath::Abs(nd.Dot(a1 - b2));
}
};
inline void PW::ClearHitInfo()
{
hitInfo.Clear();
}
inline void PW::ConstructGeo()
{
An.clear();
Ca.clear();
std::pair<TVector3, TVector3> p1; // anode
std::pair<TVector3, TVector3> q1; // cathode
// anode and cathode start at pos-Y axis and count in right-Hand
// anode wire shift is right-hand.
// cathode wire shift is left-hand.
for (int i = 0; i < nWire; i++)
{
// Anode rotate right-hand
//updated Feb 2026, Sudarsan B
p1.first.SetXYZ(radiusA * TMath::Cos(TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
radiusA * TMath::Sin(TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
zLen / 2);
p1.second.SetXYZ(radiusA * TMath::Cos(TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()),
radiusA * TMath::Sin(TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()),
-zLen / 2);
An.push_back(p1);
// Cathod rotate left-hand with the 3 wire offset accounted for (+1 from the calculated offset from the PC coincidence spectrum)
q1.first.SetXYZ(radiusC * TMath::Cos(TMath::TwoPi() / nWire * (i + wireShift + 1) + TMath::PiOver2()),
radiusC * TMath::Sin(TMath::TwoPi() / nWire * (i + wireShift + 1) + TMath::PiOver2()),
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);
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
// NOT NECESSARY ANY MORE, HAS BEEN IMCORPORATED INTO THE WIREOFFSET IN THE BEGINNING
// std::rotate(Ca.rbegin(), Ca.rbegin() + 4, Ca.rend());
dAngle = wireShift * TMath::TwoPi() / nWire;
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)); //chord length subtending an angle alpha is 2rsin(alpha/2)
}
inline std::vector<std::vector<std::tuple<int,double,double>>>
PW::Make_Clusters(std::unordered_map<int,std::tuple<int,double,double>> wireEvents) {
std::vector<std::vector<std::tuple<int,double,double>>> wireClusters;
std::vector<std::tuple<int,double,double>> wireCluster;
//TODO: Write a macro once, call it twice
int wirecount=0;
while(wirecount < 24) {
if(wireEvents.find(wirecount)==wireEvents.end()) {
wirecount++;
continue;
}
wireCluster.clear();
int ctr2=wirecount;
do {
wireCluster.emplace_back(wireEvents[ctr2]);
ctr2+=1;
if(ctr2==24 || ctr2-wirecount == 7) break; //loose logic, needs to be looked at.
} while(wireEvents.find(ctr2)!=wireEvents.end());
wireClusters.push_back(std::move(wireCluster));
wirecount = ctr2; //we already dealt with wires until the last value of ctr2
}
if(wireClusters.size() > 1) { //Deal with wraparound if required
auto first_cluster = wireClusters.front(); //front and back provide references to the elements themselves. less copy, can modify etc
auto last_cluster = wireClusters.back();
if(std::get<0>(last_cluster.back())==23 && std::get<0>(first_cluster.front())==0) {
last_cluster.insert(last_cluster.end(),first_cluster.begin(),first_cluster.end());
}
wireClusters.erase(wireClusters.begin()); //canonically, erase() needs an iterator, hence begin() not front()
//TODO: Can also deal with 'gaps' of missing wires similarly. end of one segment and beginning of another segment will be separated by missing wire --> combine the two
//TODO: Also needs some development regarding the time-correlation. Don't put wires in the same cluster if they aren't time coincident
}
return wireClusters;
/*if(aClusters.size()>1 || cClusters.size() > 1) {
std::cout << " ============== " << std::endl;
}
if(aClusters.size()>1 && cClusters.size() >=1) {
std::cout << aClusters.size() << " new anode clusters ----> " << std::endl;
int cc=1;
for(auto ac : aClusters) {
std::cout << " Cluster " << cc << std::endl;
double first_ts = std::get<2>(ac.at(0));
for(auto item : ac) {
std::cout << " \t" << std::get<0>(item) << " " << std::get<1>(item) << " " << std::get<2>(item)-first_ts << std::endl;
}
std::cout << " ------" << std::endl;
cc++;
}
}
if(cClusters.size()>=1 ) {
std::cout << cClusters.size() << " new cathode clusters ----> " << std::endl;
int cc=1;
for(auto ac : cClusters) {
std::cout << " Cluster " << cc << std::endl;
double first_ts = std::get<2>(ac.at(0));
for(auto item : ac) {
std::cout << " \t" << std::get<0>(item) << " " << std::get<1>(item) << " " << std::get<2>(item)-first_ts << std::endl;
}
std::cout << " ------" << std::endl;
cc++;
}
} */
}
inline std::tuple<std::pair<TVector3, TVector3>, double, double, double>
PW::GetPseudoWire(const std::vector<std::tuple<int,double,double>>& cluster, std::string type) {
std::pair<TVector3,TVector3> avgvec = std::pair(TVector3(0,0,0),TVector3(0,0,0));
double sumEnergy = 0;
double maxEnergy = 0;
double tsMaxEnergy = 0;
if(type=="ANODE") {
//if(cluster.size()>1) std::cout << " -------anodes" << std::endl;
for( auto wire : cluster) {
avgvec.first += std::get<1>(wire)*TVector3(An.at(std::get<0>(wire)).first.X(), An.at(std::get<0>(wire)).first.Y(), 0) ;
avgvec.second += std::get<1>(wire)*TVector3(An.at(std::get<0>(wire)).second.X(), An.at(std::get<0>(wire)).second.Y(), 0);
sumEnergy += std::get<1>(wire);
if(std::get<1>(wire) > maxEnergy) {
maxEnergy = std::get<1>(wire);
tsMaxEnergy = std::get<2>(wire);
}
/*if(cluster.size()>1) {
std::cout << "\t\t ch:" << std::get<0>(wire) << " " << std::get<1>(wire) << " " << std::get<2>(wire) << std::endl;
std::cout << "\t\t w1(r,phi,z):" << An.at(std::get<0>(wire)).first.Perp() << " " << An.at(std::get<0>(wire)).first.Phi()*180/M_PI << " " << An.at(std::get<0>(wire)).first.Z() << std::endl;
std::cout << "\t\t w2(r,phi,z):" << An.at(std::get<0>(wire)).second.Perp() << " " << An.at(std::get<0>(wire)).second.Phi()*180/M_PI << " " << An.at(std::get<0>(wire)).second.Z() << std::endl;
}*/
}
avgvec.first = avgvec.first*(1.0/sumEnergy);
avgvec.second = avgvec.second*(1.0/sumEnergy);
double phi1 = avgvec.first.Phi();
double phi2 = avgvec.second.Phi();
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);
/*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 avg2(r,phi,z):" << avgvec.second.Perp() << " " << avgvec.second.Phi()*180/M_PI << " " << avgvec.second.Z() << std::endl;
}*/
} else if(type =="CATHODE") {
for( auto wire : cluster) {
avgvec.first += std::get<1>(wire)*TVector3(Ca.at(std::get<0>(wire)).first.X(), Ca.at(std::get<0>(wire)).first.Y(), 0) ;
avgvec.second += std::get<1>(wire)*TVector3(Ca.at(std::get<0>(wire)).second.X(), Ca.at(std::get<0>(wire)).second.Y(), 0);
sumEnergy += std::get<1>(wire);
if(std::get<1>(wire) > maxEnergy) {
maxEnergy = std::get<1>(wire);
tsMaxEnergy = std::get<2>(wire);
}
}
avgvec.first = avgvec.first*(1.0/sumEnergy);
avgvec.second = avgvec.second*(1.0/sumEnergy);
double phi1 = avgvec.first.Phi();
double phi2 = avgvec.second.Phi();
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);
}
return std::tuple(avgvec, sumEnergy, maxEnergy, tsMaxEnergy);
}
inline std::tuple<TVector3,double,double,double,double,double,double,double> PW::FindCrossoverProperties(const std::vector<std::tuple<int,double,double>>& a_cluster,
const std::vector<std::tuple<int,double,double>>& c_cluster) {
//std::pair<TVector3, TVector3> apwire = GetPseudoWire(a_cluster,"ANODE",anodeSumE);
//std::pair<TVector3, TVector3> cpwire = GetPseudoWire(c_cluster,"CATHODE",cathodeSumE);
auto [apwire, apSumE, apMaxE, apTSMaxE] = GetPseudoWire(a_cluster,"ANODE");
auto [cpwire, cpSumE, cpMaxE, cpTSMaxE] = GetPseudoWire(c_cluster,"CATHODE");
TVector3 crossover;
crossover.Clear();
TVector3 a, c, diff;
double a2, ac, c2, adiff, cdiff, denom, alpha=0;
if(apSumE && cpSumE) {
a = apwire.first - apwire.second;
c = cpwire.first - cpwire.second;
diff = apwire.first - cpwire.first;
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 = apwire.first + alpha*a;
if(crossover.z() < -190 || crossover.Z() > 190 ) {
alpha = 9999999;
apSumE=-1; cpSumE=-1;
apMaxE=-1; cpMaxE=-1;
apTSMaxE=-1; cpTSMaxE=-1;
}
}
//std::cout << apSumE << " " << cpSumE << " " << " " << crossover.Perp() << std::endl;
return std::tuple(crossover,alpha,apSumE,cpSumE,apMaxE,cpMaxE,apTSMaxE,cpTSMaxE);
}
inline void PW::FindWireID(TVector3 pos, TVector3 direction, bool verbose)
{
hitInfo.Clear();
double phi = direction.Phi();
for (int i = 0; i < nWire; i++)
{
double disA = 99999999;
double phiS = An[i].first.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());
if (phi > 0 && phiS > phiL)
phiL = phiL + TMath::TwoPi();
if (phi < 0 && phiS > phiL)
phiS = phiS - TMath::TwoPi();
if (phiS < phi && phi < phiL)
{
disA = Distance(pos, pos + direction, An[i].first, An[i].second);
if (disA < hitInfo.nearestDist.first)
{
hitInfo.nearestDist.first = disA;
hitInfo.nearestWire.first = i;
}
}
double disC = 99999999;
phiS = Ca[i].second.Phi() - TMath::PiOver4();
phiL = Ca[i].first.Phi() + TMath::PiOver4();
// 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)
phiS = phiS - TMath::TwoPi();
if (phiS < phi && phi < phiL)
{
disC = Distance(pos, pos + direction, Ca[i].first, Ca[i].second);
if (disC < hitInfo.nearestDist.second)
{
hitInfo.nearestDist.second = disC;
hitInfo.nearestWire.second = i;
}
}
if (verbose)
printf(" %2d | %8.2f, %8.2f\n", i, disA, disC);
}
//==== find the 2nd nearest wire
short anode1 = hitInfo.nearestWire.first;
short aaa1 = anode1 - 1;
if (aaa1 < 0)
aaa1 += nWire;
short aaa2 = (anode1 + 1) % nWire;
double haha1 = Distance(pos, pos + direction, An[aaa1].first, An[aaa1].second);
double haha2 = Distance(pos, pos + direction, An[aaa2].first, An[aaa2].second);
if (haha1 < haha2)
{
hitInfo.nextNearestWire.first = aaa1;
hitInfo.nextNearestDist.first = haha1;
}
else
{
hitInfo.nextNearestWire.first = aaa2;
hitInfo.nextNearestDist.first = haha2;
}
short cathode1 = hitInfo.nearestWire.second;
short ccc1 = cathode1 - 1;
if (ccc1 < 0)
ccc1 += nWire;
short ccc2 = (cathode1 + 1) % nWire;
haha1 = Distance(pos, pos + direction, Ca[ccc1].first, Ca[ccc1].second);
haha2 = Distance(pos, pos + direction, Ca[ccc2].first, Ca[ccc2].second);
if (haha1 < haha2)
{
hitInfo.nextNearestWire.second = ccc1;
hitInfo.nextNearestDist.second = haha1;
}
else
{
hitInfo.nextNearestWire.second = ccc2;
hitInfo.nextNearestDist.second = haha2;
}
if (verbose)
Print();
}
inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose)
{
trackPos = sx3Pos;
TVector3 n1 = (An[anodeID].first - An[anodeID].second).Cross((sx3Pos - An[anodeID].second)).Unit();
TVector3 n2 = (Ca[cathodeID].first - Ca[cathodeID].second).Cross((sx3Pos - Ca[cathodeID].second)).Unit();
// if the handiness of anode and cathode revered, it should be n2 cross n1
trackVec = (n2.Cross(n1)).Unit();
if (verbose)
printf("Theta, Phi = %f, %f \n", trackVec.Theta() * TMath::RadToDeg(), trackVec.Phi() * TMath::RadToDeg());
}
inline void PW::CalTrack2(TVector3 siPos, TVector3 anodeInt, bool verbose)
{
double mx, my;
double z;
mx = siPos.X() / (siPos.X() - anodeInt.X());
my = siPos.Y() / (siPos.Y() - anodeInt.Y());
z=siPos.Z() + mx * (anodeInt.Z() - siPos.Z());
// if (mx == my)
{
trackVec=TVector3(0,0,z);
}
if (verbose)
printf("X slope = %f and Y slope = %f \n", mx, my);
}
/*inline TVector3 PW::CalTrack3(TVector3 siPos, TVector3 anodeInt, bool verbose)
{
TVector3 v = anodeInt-siPos;
double t_minimum = -1.0*(siPos.X()*v.X()+siPos.Y()*v.Y())/(v.X()*v.X()+v.Y()*v.Y());
TVector3 vector_closest_to_z = siPos + t_minimum*v;
return vector_closest_to_z;
if (verbose)
printf("X slope = %f and Y slope = %f \n", mx, my);
}*/
inline double PW::GetZ0()
{
double x = trackPos.X();
double y = trackPos.Y();
double rho = TMath::Sqrt(x * x + y * y);
double theta = trackVec.Theta();
return trackVec.Z();
}
#endif

489
Armory/ClassPW.h.modified Normal file
View File

@ -0,0 +1,489 @@
#ifndef ClassPW_h
#define ClassPW_h
#include <cstdio>
#include <iostream>
#include <TMath.h>
#include <TVector3.h>
#include <TRandom.h>
struct PWHitInfo
{
std::pair<short, short> nearestWire; // anode, cathode
std::pair<double, double> nearestDist; // anode, cathode
std::pair<short, short> nextNearestWire; // anode, cathode
std::pair<double, double> nextNearestDist; // anode, cathode
void Clear()
{
nearestWire.first = -1;
nearestWire.second = -1;
nearestDist.first = 999999999;
nearestDist.second = 999999999;
nextNearestWire.first = -1;
nextNearestWire.second = -1;
nextNearestDist.first = 999999999;
nextNearestDist.second = 999999999;
}
};
struct Coord
{
float x, y, z;
Coord() : x(0), y(0), z(0) {}
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:
PW() { ClearHitInfo(); };
~PW() {};
PWHitInfo GetHitInfo() const { return hitInfo; }
std::pair<short, short> GetNearestID() const { return hitInfo.nearestWire; }
std::pair<double, double> GetNearestDistance() const { return hitInfo.nearestDist; }
std::pair<short, short> Get2ndNearestID() const { return hitInfo.nextNearestWire; }
std::pair<double, double> Get2ndNearestDistance() const { return hitInfo.nextNearestDist; }
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
TVector3 GetTrackPos() const { return trackPos; }
TVector3 GetTrackVec() const { return trackVec; }
double GetTrackTheta() const { return trackVec.Theta(); }
double GetTrackPhi() const { return trackVec.Phi(); }
double GetZ0();
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>
FindCrossoverProperties(const std::vector<std::tuple<int,double,double>>& a_cluster, const std::vector<std::tuple<int,double,double>>& c_cluster);
inline std::vector<std::vector<std::tuple<int,double,double>>>
Make_Clusters(std::unordered_map<int,std::tuple<int,double,double>> wireEvents);
int GetNumWire() const { return nWire; }
double GetDeltaAngle() const { return dAngle; }
double GetAnodeLength() const { return anodeLength; }
double GetCathodeLength() const { return cathodeLength; }
TVector3 GetAnodeDn(short id) const { return An[id].first; }
TVector3 GetAnodeUp(short id) const { return An[id].second; }
TVector3 GetCathodeDn(short id) const { return Ca[id].first; }
TVector3 GetCathodeUp(short id) const { return Ca[id].second; }
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 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; }
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(); }
void ClearHitInfo();
void ConstructGeo();
void FindWireID(TVector3 pos, TVector3 direction, bool verbose = false);
void CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose = false);
void CalTrack2(TVector3 sx3Pos, TVector3 anodeInt, bool verbose = false);
void Print()
{
printf(" The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nearestWire.first,
hitInfo.nearestDist.first,
hitInfo.nearestWire.second,
hitInfo.nearestDist.second);
printf(" The 2nd nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nextNearestWire.first,
hitInfo.nextNearestDist.first,
hitInfo.nextNearestWire.second,
hitInfo.nextNearestDist.second);
}
private:
PWHitInfo hitInfo;
TVector3 trackPos;
TVector3 trackVec;
const int nWire = 24;
const int wireShift = 3;
//const float zLen = 380; // mm
const float zLen = 348.6; // mm
const float radiusA = 37;
const float radiusC = 43;
double dAngle;
double anodeLength;
double cathodeLength;
// 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
double Distance(TVector3 a1, TVector3 a2, TVector3 b1, TVector3 b2)
{
TVector3 na = a1 - a2;
TVector3 nb = b1 - b2;
TVector3 nd = (na.Cross(nb)).Unit();
return TMath::Abs(nd.Dot(a1 - b2));
}
};
inline void PW::ClearHitInfo()
{
hitInfo.Clear();
}
inline void PW::ConstructGeo()
{
An.clear();
Ca.clear();
std::pair<TVector3, TVector3> p1; // anode
std::pair<TVector3, TVector3> q1; // cathode
double k = TMath::TwoPi()/24.; //48 solder thru holes, wires in every other one
double offset_a1 = -6*k-3*k;
double offset_c1 = -3*k - TMath::TwoPi()/48; //correct for a half-turn
double offset_a2 = offset_a1+3*k;
double offset_c2 = offset_c1-3*k;
for (int i = 0; i < nWire; i++)
{
// Anode rotate right-hand coming in towards +z riding with the beam. In this frame, +x is to the right, and +y down
//updated Feb 2026, Sudarsan B. Photographs indicate that anode wires twist right handed, as one moves from -z to +z with the convention above
//'First' is -z locus, 'second' is +z locus
p1.first.SetXYZ(radiusA * TMath::Cos(-k*i + offset_a1),
radiusA * TMath::Sin(-k*i + offset_a1),
-zLen / 2);
p1.second.SetXYZ(radiusA * TMath::Cos(-k*i + offset_a2),
radiusA * TMath::Sin(-k*i + offset_a2),
+zLen / 2);
// Cathodes rotate left-hand, under the same system. k is positive
q1.first.SetXYZ(radiusC * TMath::Cos(k*i + offset_c1),
radiusC * TMath::Sin(k*i + offset_c1),
-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);
}
dAngle = wireShift * TMath::TwoPi() / nWire;
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)); //chord length subtending an angle alpha is 2rsin(alpha/2)
}
inline std::vector<std::vector<std::tuple<int,double,double>>>
PW::Make_Clusters(std::unordered_map<int,std::tuple<int,double,double>> wireEvents) {
std::vector<std::vector<std::tuple<int,double,double>>> wireClusters;
std::vector<std::tuple<int,double,double>> wireCluster;
//TODO: Write a macro once, call it twice
int wirecount=0;
while(wirecount < 24) {
if(wireEvents.find(wirecount)==wireEvents.end()) {
wirecount++;
continue;
}
wireCluster.clear();
int ctr2=wirecount;
do {
wireCluster.emplace_back(wireEvents[ctr2]);
ctr2+=1;
if(ctr2==24 || ctr2-wirecount == 7) break; //loose logic, needs to be looked at.
} while(wireEvents.find(ctr2)!=wireEvents.end());
wireClusters.push_back(std::move(wireCluster));
wirecount = ctr2; //we already dealt with wires until the last value of ctr2
}
if(wireClusters.size() > 1) { //Deal with wraparound if required
auto first_cluster = wireClusters.front(); //front and back provide references to the elements themselves. less copy, can modify etc
auto last_cluster = wireClusters.back();
if(std::get<0>(last_cluster.back())==23 && std::get<0>(first_cluster.front())==0) {
last_cluster.insert(last_cluster.end(),first_cluster.begin(),first_cluster.end());
}
wireClusters.erase(wireClusters.begin()); //canonically, erase() needs an iterator, hence begin() not front()
//TODO: Can also deal with 'gaps' of missing wires similarly. end of one segment and beginning of another segment will be separated by missing wire --> combine the two
//TODO: Also needs some development regarding the time-correlation. Don't put wires in the same cluster if they aren't time coincident
}
return wireClusters;
/*if(aClusters.size()>1 || cClusters.size() > 1) {
std::cout << " ============== " << std::endl;
}
if(aClusters.size()>1 && cClusters.size() >=1) {
std::cout << aClusters.size() << " new anode clusters ----> " << std::endl;
int cc=1;
for(auto ac : aClusters) {
std::cout << " Cluster " << cc << std::endl;
double first_ts = std::get<2>(ac.at(0));
for(auto item : ac) {
std::cout << " \t" << std::get<0>(item) << " " << std::get<1>(item) << " " << std::get<2>(item)-first_ts << std::endl;
}
std::cout << " ------" << std::endl;
cc++;
}
}
if(cClusters.size()>=1 ) {
std::cout << cClusters.size() << " new cathode clusters ----> " << std::endl;
int cc=1;
for(auto ac : cClusters) {
std::cout << " Cluster " << cc << std::endl;
double first_ts = std::get<2>(ac.at(0));
for(auto item : ac) {
std::cout << " \t" << std::get<0>(item) << " " << std::get<1>(item) << " " << std::get<2>(item)-first_ts << std::endl;
}
std::cout << " ------" << std::endl;
cc++;
}
} */
}
inline std::tuple<std::pair<TVector3, TVector3>, double, double, double>
PW::GetPseudoWire(const std::vector<std::tuple<int,double,double>>& cluster, std::string type) {
std::pair<TVector3,TVector3> avgvec = std::pair(TVector3(0,0,0),TVector3(0,0,0));
double sumEnergy = 0;
double maxEnergy = 0;
double tsMaxEnergy = 0;
if(type=="ANODE") {
//if(cluster.size()>1) std::cout << " -------anodes" << std::endl;
for( auto wire : cluster) {
avgvec.first += std::get<1>(wire)*TVector3(An.at(std::get<0>(wire)).first.X(), An.at(std::get<0>(wire)).first.Y(), 0) ;
avgvec.second += std::get<1>(wire)*TVector3(An.at(std::get<0>(wire)).second.X(), An.at(std::get<0>(wire)).second.Y(), 0);
sumEnergy += std::get<1>(wire);
if(std::get<1>(wire) > maxEnergy) {
maxEnergy = std::get<1>(wire);
tsMaxEnergy = std::get<2>(wire);
}
/*if(cluster.size()>1) {
std::cout << "\t\t ch:" << std::get<0>(wire) << " " << std::get<1>(wire) << " " << std::get<2>(wire) << std::endl;
std::cout << "\t\t w1(r,phi,z):" << An.at(std::get<0>(wire)).first.Perp() << " " << An.at(std::get<0>(wire)).first.Phi()*180/M_PI << " " << An.at(std::get<0>(wire)).first.Z() << std::endl;
std::cout << "\t\t w2(r,phi,z):" << An.at(std::get<0>(wire)).second.Perp() << " " << An.at(std::get<0>(wire)).second.Phi()*180/M_PI << " " << An.at(std::get<0>(wire)).second.Z() << std::endl;
}*/
}
avgvec.first = avgvec.first*(1.0/sumEnergy);
avgvec.second = avgvec.second*(1.0/sumEnergy);
double phi1 = avgvec.first.Phi();
double phi2 = avgvec.second.Phi();
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);
/*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 avg2(r,phi,z):" << avgvec.second.Perp() << " " << avgvec.second.Phi()*180/M_PI << " " << avgvec.second.Z() << std::endl;
}*/
} else if(type =="CATHODE") {
for( auto wire : cluster) {
avgvec.first += std::get<1>(wire)*TVector3(Ca.at(std::get<0>(wire)).first.X(), Ca.at(std::get<0>(wire)).first.Y(), 0) ;
avgvec.second += std::get<1>(wire)*TVector3(Ca.at(std::get<0>(wire)).second.X(), Ca.at(std::get<0>(wire)).second.Y(), 0);
sumEnergy += std::get<1>(wire);
if(std::get<1>(wire) > maxEnergy) {
maxEnergy = std::get<1>(wire);
tsMaxEnergy = std::get<2>(wire);
}
}
avgvec.first = avgvec.first*(1.0/sumEnergy);
avgvec.second = avgvec.second*(1.0/sumEnergy);
double phi1 = avgvec.first.Phi();
double phi2 = avgvec.second.Phi();
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);
}
return std::tuple(avgvec, sumEnergy, maxEnergy, tsMaxEnergy);
}
inline std::tuple<TVector3,double,double,double,double,double,double,double> PW::FindCrossoverProperties(const std::vector<std::tuple<int,double,double>>& a_cluster,
const std::vector<std::tuple<int,double,double>>& c_cluster) {
//std::pair<TVector3, TVector3> apwire = GetPseudoWire(a_cluster,"ANODE",anodeSumE);
//std::pair<TVector3, TVector3> cpwire = GetPseudoWire(c_cluster,"CATHODE",cathodeSumE);
auto [apwire, apSumE, apMaxE, apTSMaxE] = GetPseudoWire(a_cluster,"ANODE");
auto [cpwire, cpSumE, cpMaxE, cpTSMaxE] = GetPseudoWire(c_cluster,"CATHODE");
TVector3 crossover;
crossover.Clear();
TVector3 a, c, diff;
double a2, ac, c2, adiff, cdiff, denom, alpha=0;
if(apSumE && cpSumE) {
a = apwire.first - apwire.second;
c = cpwire.first - cpwire.second;
diff = apwire.first - cpwire.first;
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 = apwire.first + alpha*a;
if(crossover.z() < -190 || crossover.Z() > 190 ) {
alpha = 9999999;
apSumE=-1; cpSumE=-1;
apMaxE=-1; cpMaxE=-1;
apTSMaxE=-1; cpTSMaxE=-1;
}
}
//std::cout << apSumE << " " << cpSumE << " " << " " << crossover.Perp() << std::endl;
return std::tuple(crossover,alpha,apSumE,cpSumE,apMaxE,cpMaxE,apTSMaxE,cpTSMaxE);
}
inline void PW::FindWireID(TVector3 pos, TVector3 direction, bool verbose)
{
hitInfo.Clear();
double phi = direction.Phi();
for (int i = 0; i < nWire; i++)
{
double disA = 99999999;
double phiS = An[i].first.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());
if (phi > 0 && phiS > phiL)
phiL = phiL + TMath::TwoPi();
if (phi < 0 && phiS > phiL)
phiS = phiS - TMath::TwoPi();
if (phiS < phi && phi < phiL)
{
disA = Distance(pos, pos + direction, An[i].first, An[i].second);
if (disA < hitInfo.nearestDist.first)
{
hitInfo.nearestDist.first = disA;
hitInfo.nearestWire.first = i;
}
}
double disC = 99999999;
phiS = Ca[i].second.Phi() - TMath::PiOver4();
phiL = Ca[i].first.Phi() + TMath::PiOver4();
// 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)
phiS = phiS - TMath::TwoPi();
if (phiS < phi && phi < phiL)
{
disC = Distance(pos, pos + direction, Ca[i].first, Ca[i].second);
if (disC < hitInfo.nearestDist.second)
{
hitInfo.nearestDist.second = disC;
hitInfo.nearestWire.second = i;
}
}
if (verbose)
printf(" %2d | %8.2f, %8.2f\n", i, disA, disC);
}
//==== find the 2nd nearest wire
short anode1 = hitInfo.nearestWire.first;
short aaa1 = anode1 - 1;
if (aaa1 < 0)
aaa1 += nWire;
short aaa2 = (anode1 + 1) % nWire;
double haha1 = Distance(pos, pos + direction, An[aaa1].first, An[aaa1].second);
double haha2 = Distance(pos, pos + direction, An[aaa2].first, An[aaa2].second);
if (haha1 < haha2)
{
hitInfo.nextNearestWire.first = aaa1;
hitInfo.nextNearestDist.first = haha1;
}
else
{
hitInfo.nextNearestWire.first = aaa2;
hitInfo.nextNearestDist.first = haha2;
}
short cathode1 = hitInfo.nearestWire.second;
short ccc1 = cathode1 - 1;
if (ccc1 < 0)
ccc1 += nWire;
short ccc2 = (cathode1 + 1) % nWire;
haha1 = Distance(pos, pos + direction, Ca[ccc1].first, Ca[ccc1].second);
haha2 = Distance(pos, pos + direction, Ca[ccc2].first, Ca[ccc2].second);
if (haha1 < haha2)
{
hitInfo.nextNearestWire.second = ccc1;
hitInfo.nextNearestDist.second = haha1;
}
else
{
hitInfo.nextNearestWire.second = ccc2;
hitInfo.nextNearestDist.second = haha2;
}
if (verbose)
Print();
}
inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose)
{
trackPos = sx3Pos;
TVector3 n1 = (An[anodeID].first - An[anodeID].second).Cross((sx3Pos - An[anodeID].second)).Unit();
TVector3 n2 = (Ca[cathodeID].first - Ca[cathodeID].second).Cross((sx3Pos - Ca[cathodeID].second)).Unit();
// if the handiness of anode and cathode revered, it should be n2 cross n1
trackVec = (n2.Cross(n1)).Unit();
if (verbose)
printf("Theta, Phi = %f, %f \n", trackVec.Theta() * TMath::RadToDeg(), trackVec.Phi() * TMath::RadToDeg());
}
inline void PW::CalTrack2(TVector3 siPos, TVector3 anodeInt, bool verbose)
{
double mx, my;
double z;
mx = siPos.X() / (siPos.X() - anodeInt.X());
my = siPos.Y() / (siPos.Y() - anodeInt.Y());
z=siPos.Z() + mx * (anodeInt.Z() - siPos.Z());
// if (mx == my)
{
trackVec=TVector3(0,0,z);
}
if (verbose)
printf("X slope = %f and Y slope = %f \n", mx, my);
}
/*inline TVector3 PW::CalTrack3(TVector3 siPos, TVector3 anodeInt, bool verbose)
{
TVector3 v = anodeInt-siPos;
double t_minimum = -1.0*(siPos.X()*v.X()+siPos.Y()*v.Y())/(v.X()*v.X()+v.Y()*v.Y());
TVector3 vector_closest_to_z = siPos + t_minimum*v;
return vector_closest_to_z;
if (verbose)
printf("X slope = %f and Y slope = %f \n", mx, my);
}*/
inline double PW::GetZ0()
{
double x = trackPos.X();
double y = trackPos.Y();
double rho = TMath::Sqrt(x * x + y * y);
double theta = trackVec.Theta();
return trackVec.Z();
}
#endif

1
Armory/SX3Geom.h
View File

@ -72,6 +72,7 @@ public:
void sx3::fillevent(const std::string& positionstring, const int subchannel, const float value) { void sx3::fillevent(const std::string& positionstring, const int subchannel, const float value) {
assert(subchannel>=0 && subchannel<4); assert(subchannel>=0 && subchannel<4);
foundevent=1;
if(positionstring=="FRONT_L") { if(positionstring=="FRONT_L") {
frontL[subchannel].push_back(value); frontL[subchannel].push_back(value);
unmatched_front_chans.insert(subchannel); unmatched_front_chans.insert(subchannel);

639
DataDump.C Normal file
View File

@ -0,0 +1,639 @@
#define DataDump_cxx
#include "DataDump.h"
#include "Armory/ClassPW.h"
#include "Armory/HistPlotter.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TBranch.h>
#include "TVector3.h"
#include <fstream>
#include <iostream>
#include <sstream>
#include <map>
#include <utility>
#include <vector>
#include <utility>
#include <algorithm>
// Global instances
PW pw_contr;
PW pwinstance;
TVector3 hitPos;
long long int gcount=0;
class Event {
public:
TVector3 pos;
int ch1=-1; //ring# for QQQ, anode# for PC
int ch2=-1; //wedge# for QQQ, cathode# for PC
double Energy1=-1; //Front for QQQ, Anode for PC
double Energy2=-1; //Back for QQQ, Cathode for PC
double Time1=-1;
double Time2=-1;
Event(TVector3 p, double E1, double E2, double T1, double T2): pos(p), Energy1(E1), Energy2(E2), Time1(T1), Time2(T2) {}
};
/*
void testfunction()
{
for(auto cathode: cathodes) {
std::unordered_set<int> chans;
chans.insert(cathode.ch);
}
}
*/
// Calibration globals
const int MAX_QQQ = 4;
const int MAX_RING = 16;
const int MAX_WEDGE = 16;
double qqqGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
bool qqqGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
double qqqCalib[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
bool qqqCalibValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
// TCutg *cutQQQ;
// PC Arrays
double pcSlope[48];
double pcIntercept[48];
HistPlotter *plotter;
bool HitNonZero;
bool sx3ecut;
bool qqqEcut;
void DataDump::Begin(TTree * /*tree*/)
{
TString option = GetOption();
plotter = new HistPlotter("Analyzer_QQQ.root", "TFILE");
pw_contr.ConstructGeo();
pwinstance.ConstructGeo();
// ---------------------------------------------------------
// 1. CRITICAL FIX: Initialize PC Arrays to Default (Raw)
// ---------------------------------------------------------
for (int i = 0; i < 48; i++) {
pcSlope[i] = 1.0; // Default slope = 1 (preserves Raw energy)
pcIntercept[i] = 0.0; // Default intercept = 0
}
// Calculate Crossover Geometry ONCE
TVector3 a, c, diff;
double a2, ac, c2, adiff, cdiff, denom, alpha;
for (size_t i = 0; i < pwinstance.An.size(); i++)
{
a = pwinstance.An[i].first - pwinstance.An[i].second;
for (size_t j = 0; j < pwinstance.Ca.size(); j++)
{
c = pwinstance.Ca[j].first - pwinstance.Ca[j].second;
diff = pwinstance.An[i].first - pwinstance.Ca[j].first;
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 = pwinstance.An[i].first.X() + alpha * a.X();
Crossover[i][j][0].y = pwinstance.An[i].first.Y() + alpha * a.Y();
Crossover[i][j][0].z = pwinstance.An[i].first.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;
}
}
// Load PC Calibrations
std::ifstream inputFile("slope_intercept_results.txt");
if (inputFile.is_open())
{
std::string line;
int index;
double slope, intercept;
while (std::getline(inputFile, line))
{
std::stringstream ss(line);
ss >> index >> slope >> intercept;
if (index >= 0 && index <= 47)
{
pcSlope[index] = slope;
pcIntercept[index] = intercept;
}
}
inputFile.close();
}
else
{
std::cerr << "Error opening slope_intercept.txt" << std::endl;
}
// ... (Load QQQ Gains and Calibs - same as before) ...
{
std::string filename = "qqq_GainMatch.dat";
std::ifstream infile(filename);
if (infile.is_open())
{
int det, ring, wedge;
double gainw, gainr;
while (infile >> det >> wedge >> ring >> gainw >> gainr)
{
qqqGain[det][wedge][ring] = gainw;
qqqGainValid[det][wedge][ring] = (gainw > 0);
}
infile.close();
}
}
{
std::string filename = "qqq_Calib.dat";
std::ifstream infile(filename);
if (infile.is_open())
{
int det, ring, wedge;
double slope;
while (infile >> det >> wedge >> ring >> slope)
{
qqqCalib[det][wedge][ring] = slope;
qqqCalibValid[det][wedge][ring] = (slope > 0);
}
infile.close();
}
}
}
Bool_t DataDump::Process(Long64_t entry)
{
hitPos.Clear();
HitNonZero = false;
bool qqq1000cut = false;
b_sx3Multi->GetEntry(entry);
b_sx3ID->GetEntry(entry);
b_sx3Ch->GetEntry(entry);
b_sx3E->GetEntry(entry);
b_sx3T->GetEntry(entry);
b_qqqMulti->GetEntry(entry);
b_qqqID->GetEntry(entry);
b_qqqCh->GetEntry(entry);
b_qqqE->GetEntry(entry);
b_qqqT->GetEntry(entry);
b_pcMulti->GetEntry(entry);
b_pcID->GetEntry(entry);
b_pcCh->GetEntry(entry);
b_pcE->GetEntry(entry);
b_pcT->GetEntry(entry);
sx3.CalIndex();
qqq.CalIndex();
pc.CalIndex();
// QQQ Processing
int qqqCount = 0;
int qqqAdjCh = 0;
// REMOVE WHEN RERUNNING USING THE NEW CALIBRATION FILE
for (int i = 0; i < qqq.multi; i++)
{
if ((qqq.id[i] == 3 || qqq.id[i] == 1) && qqq.ch[i] < 16)
{
qqq.ch[i] = 16 - qqq.ch[i];
}
}
for (int i = 0; i < qqq.multi; i++)
{
if (qqq.id[i] == 0 && qqq.ch[i] >= 16)
{
qqq.ch[i] = 31 - qqq.ch[i] + 16;
}
}
std::vector<std::tuple<int,int,double,int,double>> qqqlist;
std::vector<Event> QQQ_Events, PC_Events;
std::vector<Event> QQQ_Events_Raw, PC_Events_Raw;
bool PCQQQTimeCut = false;
for (int i = 0; i < qqq.multi; i++)
{
for (int j = i + 1; j < qqq.multi; j++)
{
if (qqq.id[i] == qqq.id[j])
{
qqqCount++;
int chWedge = -1;
int chRing = -1;
double eWedge = 0.0;
double eWedgeMeV = 0.0;
double eRing = 0.0;
double eRingMeV = 0.0;
double tRing = 0.0;
double tWedge = 0.0;
if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && qqqGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16])
{
chWedge = qqq.ch[i];
eWedge = qqq.e[i] * qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
chRing = qqq.ch[j] - 16;
eRing = qqq.e[j];
tRing = static_cast<double>(qqq.t[j]);
tWedge = static_cast<double>(qqq.t[i]);
}
else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && qqqGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16])
{
chWedge = qqq.ch[j];
eWedge = qqq.e[j] * qqqGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
chRing = qqq.ch[i] - 16;
eRing = qqq.e[i];
tRing = static_cast<double>(qqq.t[i]);
tWedge = static_cast<double>(qqq.t[j]);
}
else
continue;
if (qqqCalibValid[qqq.id[i]][chWedge][chRing]) {
eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chWedge][chRing] / 1000;
eRingMeV = eRing * qqqCalib[qqq.id[i]][chWedge][chRing] / 1000;
}
else
continue;
double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
double rho = 50. + 40. / 16. * (chRing + 0.5);
Event qqqevent(TVector3(rho*TMath::Cos(theta),rho*TMath::Sin(theta),23+75+30), eRingMeV, eWedgeMeV, tRing, tWedge);
Event qqqeventr(TVector3(rho*TMath::Cos(theta),rho*TMath::Sin(theta),23+75+30), eRing, eWedge, tRing, tWedge);
QQQ_Events.push_back(qqqevent);
QQQ_Events_Raw.push_back(qqqeventr);
qqqlist.push_back(std::tuple(qqq.id[i],chRing,eRingMeV,chWedge,eWedgeMeV));
} //end if qqq.id[i] == qqq.id[j]
}//inner qqq loop, j
}//outer qqq loop, i
// PC Gain Matching and Filling
double anodeT = -99999;
double cathodeT = 99999;
int anodeIndex = -1;
int cathodeIndex = -1;
int aID = 0;
int cID = 0;
double aE = 0;
double cE = 0;
double aESum = 0;
double cESum = 0;
double aEMax = 0;
int aIDMax = 0;
anodeHits.clear();
cathodeHits.clear();
anodeTimes.clear();
cathodeTimes.clear();
corrcatMax.clear();
corranoMax.clear();
std::array<bool,24> caths_seen{0}, anos_seen{0};
std::vector<std::tuple<int,double,double>> anodeChunks, cathodeChunks;
for (int i = 0; i < pc.multi; i++)
{
if (pc.e[i] > 4)
{
;//plotter->Fill2D("PC_Index_Vs_Energy", 48, 0, 48, 2000, 0, 30000, pc.index[i], static_cast<double>(pc.e[i]), "hRawPC");
} else
continue;
if (pc.index[i] < 48) {
pc.e[i] = pcSlope[pc.index[i]] * pc.e[i] + pcIntercept[pc.index[i]];
//plotter->Fill2D("PC_Index_VS_GainMatched_Energy", 48, 0, 48, 2000, 0, 30000, pc.index[i], pc.e[i], "hGMPC");
}
if (pc.index[i] < 24) {
anodeT = static_cast<double>(pc.t[i]);
anodeIndex = pc.index[i];
anos_seen[anodeIndex] = 1;
anodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
anodeTimes.push_back(anodeT);
anodeChunks.push_back(std::tuple<int,double,double>(pc.index[i],anodeT,pc.e[i]));
} else {
cathodeT = static_cast<double>(pc.t[i]);
cathodeIndex = pc.index[i] - 24;
caths_seen[cathodeIndex] = 1;
cathodeChunks.push_back(std::tuple<int,double,double>(pc.index[i]-24,cathodeT,pc.e[i]));
cathodeHits.push_back(std::pair<int, double>(pc.index[i] - 24, pc.e[i]));
cathodeTimes.push_back(cathodeT);
}
}//end of pc.multi loop
if(anodeHits.size() && cathodeHits.size()) {
for(size_t ii=0; ii<anodeHits.size(); ii++) {
const auto& an = anodeHits.at(ii);
const auto& at = anodeTimes.at(ii);
for(size_t jj=0; jj<cathodeHits.size(); jj++) {
const auto& ca = cathodeHits.at(jj);
const auto& ct = cathodeTimes.at(jj);
plotter->Fill2D("ach_minus_cch_vs_ach",60,-30,30,24,0,24,an.first-ca.first,an.first);
plotter->Fill2D("ach_minus_cch_vs_dt",60,-30,30,400,-1000,1000,an.first-ca.first,at-ct);
plotter->Fill2D("ach_vs_cch",24,0,24,24,0,24,an.first,ca.first);
}
}
gcount++;
}
bool all_three = anodeHits.size() > 0 && cathodeHits.size() > 0 && qqqlist.size() > 0;
if(all_three) std::cout << "---" << std::endl;
for(size_t ii=0; ii<anodeHits.size() && all_three; ii++) {
const auto& an = anodeHits.at(ii);
const auto& at = anodeTimes.at(ii);
std::cout << "an" << std::setprecision(16) << ", " << an.first << ", " << an.second << ", " << at << " ,-1, -1" << std::endl;
}
for(size_t jj=0; jj<cathodeHits.size() && all_three; jj++) {
const auto& ca = cathodeHits.at(jj);
const auto& ct = cathodeTimes.at(jj);
std::cout << "ca" << std::setprecision(16) << ", " << ca.first << ", " << ca.second << ", " << ct << " ,-1, -1" << std::endl;
}
for(size_t jj=0; jj<qqqlist.size() && all_three; jj++) {
const auto[id,chr,er,chw,ew] = qqqlist.at(jj);
std::cout << "q" << std::setprecision(16) << ", " << id << ", " << chr << ", " << er << ", " << chw << ", " << ew << std::endl;
}
if(all_three) std::cout << " --end-- " << std::endl;
if(gcount == 100)
return -1;
// std::sort(anodeChunks.begin(),anodeChunks.end(),);
if (anodeHits.size() >= 1 && cathodeHits.size() >= 1)
{
// 2. CRITICAL FIX: Define reference vector 'a'
// In Analyzer.cxx, 'a' was left over from the loop. We use the first anode wire as reference here.
// (Assuming pwinstance.An is populated and wires are generally parallel).
TVector3 refAnode = pwinstance.An[0].first - pwinstance.An[0].second;
for (const auto &anode : anodeHits)
{
aID = anode.first;
aE = anode.second;
aESum += aE;
if (aE > aEMax)
{
aEMax = aE;
aIDMax = aID;
}
}
for (const auto &cathode : cathodeHits)
{
cID = cathode.first;
cE = cathode.second;
for (int j = -4; j < 3; j++)
{
if ((aIDMax + 24 + j) % 24 == 23 - cID)
{
corrcatMax.push_back(std::pair<int, double>(cID, cE));
cESum += cE;
}
}
}
}
TVector3 anodeIntersection;
anodeIntersection.Clear();
if (corrcatMax.size() > 0)
{
double x = 0, y = 0, z = 0;
for (const auto &corr : corrcatMax)
{
if (Crossover[aIDMax][corr.first][0].z > 9000000)
continue;
if (cESum > 0)
{
x += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].x;
y += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].y;
z += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].z;
}
}
if (x == 0 && y == 0 && z == 0)
;
// to ignore events with no valid crossover points
else
anodeIntersection = TVector3(x, y, z);
}
bool PCQQQPhiCut = false;
// flip the algorithm for cathode 1 multi anode events
if ((hitPos.Phi() > (anodeIntersection.Phi() - TMath::PiOver4())) && (hitPos.Phi() < (anodeIntersection.Phi() + TMath::PiOver4()))) {
PCQQQPhiCut = true;
}
if (anodeIntersection.Z() != 0)
{
plotter->Fill1D("PC_Z_Projection", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
plotter->Fill2D("Z_Proj_VsDelTime", 600, -300, 300, 200, -2000, 2000, anodeIntersection.Z(), anodeT - cathodeT, "hPCzQQQ");
plotter->Fill2D("IntPhi_vs_QQQphi", 100, -200, 200, 80, -200, 200, anodeIntersection.Phi() * 180. / TMath::Pi(), hitPos.Phi() * 180. / TMath::Pi(), "hPCQQQ");
plotter->Fill2D("Inttheta_vs_QQQtheta", 90, 0, 180, 20, 0, 45, anodeIntersection.Theta() * 180. / TMath::Pi(), hitPos.Theta() * 180. / TMath::Pi(), "hPCQQQ");
plotter->Fill2D("Inttheta_vs_QQQtheta_TC" + std::to_string(PCQQQTimeCut), 90, 0, 180, 20, 0, 45, anodeIntersection.Theta() * 180. / TMath::Pi(), hitPos.Theta() * 180. / TMath::Pi(), "hPCQQQ");
plotter->Fill2D("IntPhi_vs_QQQphi_TC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 100, -200, 200, 80, -200, 200, anodeIntersection.Phi() * 180. / TMath::Pi(), hitPos.Phi() * 180. / TMath::Pi(), "hPCQQQ");
}
if (anodeIntersection.Z() != 0 && cathodeHits.size() >= 2)
plotter->Fill1D("PC_Z_Projection_TC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 1)
{
plotter->Fill1D("PC_Z_proj_1C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_1C", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hPCzQQQ");
}
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 2)
{
plotter->Fill1D("PC_Z_proj_2C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_2C", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hGMPC");
}
if (anodeIntersection.Z() != 0 && cathodeHits.size() > 2)
{
plotter->Fill1D("PC_Z_proj_nC", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_nC", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hGMPC");
}
if (anodeHits.size() > 0 && cathodeHits.size() > 0)
plotter->Fill2D("AHits_vs_CHits", 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
// make another plot with nearest neighbour constraint
bool hasNeighbourAnodes = false;
bool hasNeighbourCathodes = false;
// 1. Check Anodes for neighbours (including wrap-around 0-23)
for (size_t i = 0; i < anodeHits.size(); i++)
{
for (size_t j = i + 1; j < anodeHits.size(); j++)
{
int diff = std::abs(anodeHits[i].first - anodeHits[j].first);
if (diff == 1 || diff == 23)
{ // 23 handles the cylindrical wrap
hasNeighbourAnodes = true;
break;
}
}
if (hasNeighbourAnodes)
break;
}
// 2. Check Cathodes for neighbours (including wrap-around 0-23)
for (size_t i = 0; i < cathodeHits.size(); i++)
{
for (size_t j = i + 1; j < cathodeHits.size(); j++)
{
int diff = std::abs(cathodeHits[i].first - cathodeHits[j].first);
if (diff == 1 || diff == 23)
{
hasNeighbourCathodes = true;
break;
}
}
if (hasNeighbourCathodes)
break;
}
// ---------------------------------------------------------
// FILL PLOTS
// ---------------------------------------------------------
if (anodeHits.size() > 0 && cathodeHits.size() > 0)
{
plotter->Fill2D("AHits_vs_CHits_NA" + std::to_string(hasNeighbourAnodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
plotter->Fill2D("AHits_vs_CHits_NC" + std::to_string(hasNeighbourCathodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
// Constraint Plot: Only fill if BOTH planes have adjacent hits
// This effectively removes events with only isolated single-wire hits (noise)
if (hasNeighbourAnodes && hasNeighbourCathodes)
{
plotter->Fill2D("AHits_vs_CHits_NN", 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
}
}
if (HitNonZero && anodeIntersection.Z() != 0)
{
pw_contr.CalTrack2(hitPos, anodeIntersection);
plotter->Fill1D("VertexRecon", 600, -300, 300, pw_contr.GetZ0());
plotter->Fill1D("VertexRecon_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -300, 300, pw_contr.GetZ0());
if (cathodeHits.size() == 2)
plotter->Fill1D("VertexRecon_2c_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -300, 300, pw_contr.GetZ0());
}
for (int i = 0; i < qqq.multi; i++)
{
if (PCQQQTimeCut)
{
plotter->Fill2D("PC_XY_Projection_QQQ_TimeCut" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, anodeIntersection.X(), anodeIntersection.Y(), "hPCQQQ");
}
plotter->Fill2D("PC_XY_Projection_QQQ" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, anodeIntersection.X(), anodeIntersection.Y(), "hPCQQQ");
for (int j = i + 1; j < qqq.multi; j++)
{
if (qqq.id[i] == qqq.id[j])
{
int chWedge = -1;
int chRing = -1;
double eWedge = 0.0;
double eWedgeMeV = 0.0;
double eRing = 0.0;
double eRingMeV = 0.0;
double tRing = 0.0;
int qqqID = -1;
if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && qqqGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16])
{
chWedge = qqq.ch[i];
eWedge = qqq.e[i] * qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
chRing = qqq.ch[j] - 16;
eRing = qqq.e[j];
tRing = static_cast<double>(qqq.t[j]);
qqqID = qqq.id[i];
}
else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && qqqGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16])
{
chWedge = qqq.ch[j];
eWedge = qqq.e[j] * qqqGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
chRing = qqq.ch[i] - 16;
tRing = static_cast<double>(qqq.t[i]);
eRing = qqq.e[i];
qqqID = qqq.id[i];
}
else
continue;
if (qqqCalibValid[qqq.id[i]][chRing][chWedge])
{
eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000;
eRingMeV = eRing * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000;
}
else
continue;
// if (anodeIntersection.Z() != 0)
{
plotter->Fill2D("PC_Z_vs_QQQRing", 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCzQQQ");
}
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 2)
{
plotter->Fill2D("PC_Z_vs_QQQRing_2C", 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCzQQQ");
plotter->Fill2D("PC_Z_vs_QQQRing_2C" + std::to_string(qqq.id[i]), 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCzQQQ");
plotter->Fill2D("PC_Z_vs_QQQWedge_2C", 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chWedge, "hPCzQQQ");
}
plotter->Fill2D("Vertex_V_QQQRingTC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 600, -300, 300, 16, 0, 16, pw_contr.GetZ0(), chRing, "hPCQQQ");
double phi = TMath::ATan2(anodeIntersection.Y(), anodeIntersection.X()) * 180. / TMath::Pi();
plotter->Fill2D("PolarAngle_Vs_QQQWedge" + std::to_string(qqqID), 360, -360, 360, 16, 0, 16, phi, chWedge, "hPCQQQ");
// plotter->Fill2D("EdE_PC_vs_QQQ_timegate_ls1000"+std::to_string())
plotter->Fill2D("PC_Z_vs_QQQRing_Det" + std::to_string(qqqID), 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCQQQ");
//double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
//double rho = 50. + 40. / 16. * (chRing + 0.5);
for (int k = 0; k < pc.multi; k++)
{
if(pc.index[k] >= 24)
continue;
double sinTheta = TMath::Sin(hitPos.Theta());
plotter->Fill2D("CalibratedQQQE_RvsPCE_TC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 1000, 0, 10, 2000, 0, 30000, eRingMeV, pc.e[k]*sinTheta, "hPCQQQ");
plotter->Fill2D("CalibratedQQQE_WvsPCE_TC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 1000, 0, 10, 2000, 0, 30000, eWedgeMeV, pc.e[k]*sinTheta, "hPCQQQ");
plotter->Fill2D("PCQQQ_dTimevsdPhi", 200, -2000, 2000, 80, -200, 200, tRing - static_cast<double>(pc.t[k]), (hitPos.Phi()-anodeIntersection.Phi()) * 180. / TMath::Pi(), "hTiming");
}
}
}
}
for (int i = 0; i < sx3.multi; i++)
{
// plotting sx3 strip hits vs anode phi
if (sx3.ch[i] < 8)
plotter->Fill2D("AnodePhi_vs_SX3Strip", 100, -200, 200, 8 * 24, 0, 8 * 24, anodeIntersection.Phi() * 180. / TMath::Pi(), sx3.id[i] * 8 + sx3.ch[i]);
}
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 3)
{
plotter->Fill1D("PC_Z_proj_3C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
}
plotter->Fill2D("AnodeMaxE_Vs_Cathode_Sum_Energy", 2000, 0, 30000, 2000, 0, 30000, aEMax, cESum, "hGMPC");
plotter->Fill1D("Correlated_Cathode_MaxAnode", 6, 0, 5, corrcatMax.size(), "hGMPC");
plotter->Fill2D("Correlated_Cathode_VS_MaxAnodeEnergy", 6, 0, 5, 2000, 0, 30000, corrcatMax.size(), aEMax, "hGMPC");
plotter->Fill1D("AnodeHits", 12, 0, 11, anodeHits.size(), "hGMPC");
plotter->Fill2D("AnodeMaxE_vs_AnodeHits", 12, 0, 11, 2000, 0, 30000, anodeHits.size(), aEMax, "hGMPC");
if (anodeHits.size() < 1)
{
plotter->Fill1D("NoAnodeHits_CathodeHits", 6, 0, 5, cathodeHits.size(), "hGMPC");
}
return kTRUE;
}
void DataDump::Terminate()
{
plotter->FlushToDisk();
}

132
DataDump.h Normal file
View File

@ -0,0 +1,132 @@
#ifndef DataDump_h
#define DataDump_h
#include <TROOT.h>
#include <TChain.h>
#include <TApplication.h>
#include <TFile.h>
#include <TSelector.h>
#include <iomanip>
#include <vector> // Required for vectors
#include <utility> // Required for std::pair
#include "Armory/ClassDet.h"
#include "Armory/ClassPW.h" // YOU ADDED THIS (Correct! Defines Coord)
class DataDump : public TSelector {
public :
TTree *fChain; //!pointer to the analyzed TTree or TChain
// Declaration of leaf types
Det sx3;
Det qqq;
Det pc ;
Det misc;
ULong64_t evID;
UInt_t run;
// List of branches
TBranch *b_eventID; //!
TBranch *b_run; //!
TBranch *b_sx3Multi; //!
TBranch *b_sx3ID; //!
TBranch *b_sx3Ch; //!
TBranch *b_sx3E; //!
TBranch *b_sx3T; //!
TBranch *b_qqqMulti; //!
TBranch *b_qqqID; //!
TBranch *b_qqqCh; //!
TBranch *b_qqqE; //!
TBranch *b_qqqT; //!
TBranch *b_pcMulti; //!
TBranch *b_pcID; //!
TBranch *b_pcCh; //!
TBranch *b_pcE; //!
TBranch *b_pcT; //!
TBranch *b_miscMulti; //!
TBranch *b_miscID; //!
TBranch *b_miscCh; //!
TBranch *b_miscE; //!
TBranch *b_miscT; //!
TBranch *b_miscTf; //!
// 1. Geometry Cache
Coord Crossover[24][24][2];
// 2. Persistent Vectors (REQUIRED for the optimized .cxx to work)
std::vector<std::pair<int, double>> anodeHits;
std::vector<std::pair<int, double>> cathodeHits;
std::vector<std::pair<int, double>> corrcatMax;
std::vector<std::pair<int, double>> corranoMax;
std::vector<double> cathodeTimes;
std::vector<double> anodeTimes;
DataDump(TTree * /*tree*/ =0) : fChain(0) { }
virtual ~DataDump() { }
virtual Int_t Version() const { return 2; }
virtual void Begin(TTree *tree);
virtual void SlaveBegin(TTree *tree);
virtual void Init(TTree *tree);
virtual Bool_t Notify();
virtual Bool_t Process(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
virtual void SetOption(const char *option) { fOption = option; }
virtual void SetObject(TObject *obj) { fObject = obj; }
virtual void SetInputList(TList *input) { fInput = input; }
virtual TList *GetOutputList() const { return fOutput; }
virtual void SlaveTerminate();
virtual void Terminate();
ClassDef(DataDump,0);
};
#endif
#ifdef DataDump_cxx
void DataDump::Init(TTree *tree){
if (!tree) return;
fChain = tree;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("evID", &evID, &b_eventID);
fChain->SetBranchAddress("run", &run, &b_run);
sx3.SetDetDimension(24,12);
qqq.SetDetDimension(4,32);
pc.SetDetDimension(2,24);
fChain->SetBranchAddress("sx3Multi", &sx3.multi, &b_sx3Multi);
fChain->SetBranchAddress("sx3ID", &sx3.id, &b_sx3ID);
fChain->SetBranchAddress("sx3Ch", &sx3.ch, &b_sx3Ch);
fChain->SetBranchAddress("sx3E", &sx3.e, &b_sx3E);
fChain->SetBranchAddress("sx3T", &sx3.t, &b_sx3T);
fChain->SetBranchAddress("qqqMulti", &qqq.multi, &b_qqqMulti);
fChain->SetBranchAddress("qqqID", &qqq.id, &b_qqqID);
fChain->SetBranchAddress("qqqCh", &qqq.ch, &b_qqqCh);
fChain->SetBranchAddress("qqqE", &qqq.e, &b_qqqE);
fChain->SetBranchAddress("qqqT", &qqq.t, &b_qqqT);
fChain->SetBranchAddress("pcMulti", &pc.multi, &b_pcMulti);
fChain->SetBranchAddress("pcID", &pc.id, &b_pcID);
fChain->SetBranchAddress("pcCh", &pc.ch, &b_pcCh);
fChain->SetBranchAddress("pcE", &pc.e, &b_pcE);
fChain->SetBranchAddress("pcT", &pc.t, &b_pcT);
fChain->SetBranchAddress("miscMulti", &misc.multi, &b_miscMulti);
fChain->SetBranchAddress("miscID", &misc.id, &b_miscID);
fChain->SetBranchAddress("miscCh", &misc.ch, &b_miscCh);
fChain->SetBranchAddress("miscE", &misc.e, &b_miscE);
fChain->SetBranchAddress("miscT", &misc.t, &b_miscT);
}
Bool_t DataDump::Notify(){
return kTRUE;
}
void DataDump::SlaveBegin(TTree * /*tree*/){
// TString option = GetOption();
}
void DataDump::SlaveTerminate(){
}
#endif // #ifdef DataDump_cxx

979
MakeVertex.C
View File

File diff suppressed because it is too large Load Diff

38
MakeVertex.h
View File

@ -51,8 +51,6 @@ public :
TBranch *b_miscT; //! TBranch *b_miscT; //!
TBranch *b_miscTf; //! TBranch *b_miscTf; //!
// 1. Geometry Cache
Coord Crossover[24][24][2];
// 2. Persistent Vectors (REQUIRED for the optimized .cxx to work) // 2. Persistent Vectors (REQUIRED for the optimized .cxx to work)
std::vector<std::pair<int, double>> anodeHits; std::vector<std::pair<int, double>> anodeHits;
@ -64,17 +62,29 @@ public :
MakeVertex(TTree * /*tree*/ =0) : fChain(0) { } MakeVertex(TTree * /*tree*/ =0) : fChain(0) { }
virtual ~MakeVertex() { } virtual ~MakeVertex() { }
virtual Int_t Version() const { return 2; } virtual Int_t Version() const {
return 2;
}
virtual void Begin(TTree *tree); virtual void Begin(TTree *tree);
virtual void SlaveBegin(TTree *tree); virtual void SlaveBegin(TTree *tree);
virtual void Init(TTree *tree); virtual void Init(TTree *tree);
virtual Bool_t Notify(); virtual Bool_t Notify();
virtual Bool_t Process(Long64_t entry); virtual Bool_t Process(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; } virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) {
virtual void SetOption(const char *option) { fOption = option; } return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0;
virtual void SetObject(TObject *obj) { fObject = obj; } }
virtual void SetInputList(TList *input) { fInput = input; } virtual void SetOption(const char *option) {
virtual TList *GetOutputList() const { return fOutput; } fOption = option;
}
virtual void SetObject(TObject *obj) {
fObject = obj;
}
virtual void SetInputList(TList *input) {
fInput = input;
}
virtual TList *GetOutputList() const {
return fOutput;
}
virtual void SlaveTerminate(); virtual void SlaveTerminate();
virtual void Terminate(); virtual void Terminate();
@ -84,7 +94,7 @@ public :
#endif #endif
#ifdef MakeVertex_cxx #ifdef MakeVertex_cxx
void MakeVertex::Init(TTree *tree){ void MakeVertex::Init(TTree *tree) {
if (!tree) return; if (!tree) return;
fChain = tree; fChain = tree;
@ -112,21 +122,21 @@ void MakeVertex::Init(TTree *tree){
fChain->SetBranchAddress("pcCh", &pc.ch, &b_pcCh); fChain->SetBranchAddress("pcCh", &pc.ch, &b_pcCh);
fChain->SetBranchAddress("pcE", &pc.e, &b_pcE); fChain->SetBranchAddress("pcE", &pc.e, &b_pcE);
fChain->SetBranchAddress("pcT", &pc.t, &b_pcT); fChain->SetBranchAddress("pcT", &pc.t, &b_pcT);
fChain->SetBranchAddress("miscMulti", &misc.multi, &b_miscMulti); /*fChain->SetBranchAddress("miscMulti", &misc.multi, &b_miscMulti);
fChain->SetBranchAddress("miscID", &misc.id, &b_miscID); fChain->SetBranchAddress("miscID", &misc.id, &b_miscID);
fChain->SetBranchAddress("miscCh", &misc.ch, &b_miscCh); fChain->SetBranchAddress("miscCh", &misc.ch, &b_miscCh);
fChain->SetBranchAddress("miscE", &misc.e, &b_miscE); fChain->SetBranchAddress("miscE", &misc.e, &b_miscE);
fChain->SetBranchAddress("miscT", &misc.t, &b_miscT); fChain->SetBranchAddress("miscT", &misc.t, &b_miscT);*/
} }
Bool_t MakeVertex::Notify(){ Bool_t MakeVertex::Notify() {
return kTRUE; return kTRUE;
} }
void MakeVertex::SlaveBegin(TTree * /*tree*/){ void MakeVertex::SlaveBegin(TTree * /*tree*/) {
// TString option = GetOption(); // TString option = GetOption();
} }
void MakeVertex::SlaveTerminate(){ void MakeVertex::SlaveTerminate() {
} }
#endif // #ifdef MakeVertex_cxx #endif // #ifdef MakeVertex_cxx

958
MakeVertexSX3.C Normal file
View File

@ -0,0 +1,958 @@
#define MakeVertexSX3_cxx
#include "MakeVertexSX3.h"
#include "Armory/ClassPW.h"
#include "Armory/HistPlotter.h"
#include "Armory/SX3Geom.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TBranch.h>
#include "TVector3.h"
#include <fstream>
#include <iostream>
#include <sstream>
#include <vector>
#include <array>
#include <map>
#include <utility>
#include <algorithm>
// Global instances
PW pw_contr;
PW pwinstance;
TVector3 hitPos;
double qqqenergy, qqqtimestamp;
class Event {
public:
Event(TVector3 p, double e1, double e2, double t1, double t2) : pos(p), Energy1(e1), Energy2(e2), Time1(t1), Time2(t2) {}
Event(TVector3 p, double e1, double e2, double t1, double t2, int c1, int c2) : pos(p), Energy1(e1), Energy2(e2), Time1(t1), Time2(t2), ch1(c1), ch2(c2) {}
TVector3 pos;
int ch1=-1; //int(ch1/16) gives qqq id, ch1%16 gives ring#
int ch2=-1; //int(ch2/16) gives qqq id, ch2%16 gives wedge#
double Energy1=-1; //Front for QQQ, Anode for PC
double Energy2=-1; //Back for QQQ, Cathode for PC
double Time1=-1;
double Time2=-1;
};
// Calibration globals
const int MAX_QQQ = 4;
const int MAX_RING = 16;
const int MAX_WEDGE = 16;
double qqqGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
bool qqqGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
double qqqCalib[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
bool qqqCalibValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
// TCutg *cutQQQ;
double sx3BackGain[24][4][4] = {{{1.}}};
double sx3FrontGain[24][4] = {{1.}};
double sx3FrontOffset[24][4] = {{0.}};
double sx3RightGain[24][4] = {{1.}};
// PC Arrays
double pcSlope[48];
double pcIntercept[48];
HistPlotter *plotter;
bool HitNonZero;
bool sx3ecut;
bool qqqEcut;
void MakeVertexSX3::Begin(TTree * /*tree*/)
{
TString option = GetOption();
plotter = new HistPlotter("Analyzer_SX3.root", "TFILE");
pw_contr.ConstructGeo();
pwinstance.ConstructGeo();
// ---------------------------------------------------------
// 1. CRITICAL FIX: Initialize PC Arrays to Default (Raw)
// ---------------------------------------------------------
for (int i = 0; i < 48; i++)
{
pcSlope[i] = 1.0; // Default slope = 1 (preserves Raw energy)
pcIntercept[i] = 0.0; // Default intercept = 0
}
// Calculate Crossover Geometry ONCE
TVector3 a, c, diff;
double a2, ac, c2, adiff, cdiff, denom, alpha;
for (size_t i = 0; i < pwinstance.An.size(); i++)
{
a = pwinstance.An[i].first - pwinstance.An[i].second;
for (size_t j = 0; j < pwinstance.Ca.size(); j++)
{
c = pwinstance.Ca[j].first - pwinstance.Ca[j].second;
diff = pwinstance.An[i].first - pwinstance.Ca[j].first;
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 = pwinstance.An[i].first.X() + alpha * a.X();
Crossover[i][j][0].y = pwinstance.An[i].first.Y() + alpha * a.Y();
Crossover[i][j][0].z = pwinstance.An[i].first.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;
}
}
// Load PC Calibrations
std::ifstream inputFile("slope_intercept_results.txt");
if (inputFile.is_open())
{
std::string line;
int index;
double slope, intercept;
while (std::getline(inputFile, line))
{
std::stringstream ss(line);
ss >> index >> slope >> intercept;
if (index >= 0 && index <= 47)
{
pcSlope[index] = slope;
pcIntercept[index] = intercept;
}
}
inputFile.close();
}
else
{
std::cerr << "Error opening slope_intercept.txt" << std::endl;
}
// Load QQQ Cuts from file
// {
// std::string filename = "QQQ_PCCut.root";
// TFile *cutFile = TFile::Open(filename.c_str(), "READ");
// if (cutFile && !cutFile->IsZombie())
// {
// cutQQQ = (TCutg *)cutFile->Get("cutQQQPC");
// if (cutQQQ)
// {
// std::cout << "Loaded QQQ PC cut from " << filename << std::endl;
// }
// else
// {
// std::cerr << "Error: cutQQQPC not found in " << filename << std::endl;
// }
// cutFile->Close();
// }
// }
// ... (Load QQQ Gains and Calibs - same as before) ...
{
std::string filename = "qqq_GainMatch.dat";
std::ifstream infile(filename);
if (infile.is_open())
{
int det, ring, wedge;
double gainw, gainr;
while (infile >> det >> wedge >> ring >> gainw >> gainr)
{
qqqGain[det][wedge][ring] = gainw;
qqqGainValid[det][wedge][ring] = (gainw > 0);
// std::cout << "QQQ Gain Loaded: Det " << det << " Ring " << ring << " Wedge " << wedge << " GainW " << gainw << " GainR " << gainr << std::endl;
}
infile.close();
}
}
{
std::string filename = "qqq_Calib.dat";
std::ifstream infile(filename);
if (infile.is_open())
{
int det, ring, wedge;
double slope;
while (infile >> det >> wedge >> ring >> slope)
{
qqqCalib[det][wedge][ring] = slope;
qqqCalibValid[det][wedge][ring] = (slope > 0);
// std::cout << "QQQ Calib Loaded: Det " << det << " Ring " << ring << " Wedge " << wedge << " Slope " << slope << std::endl;
}
infile.close();
}
}
{
std::ifstream infile("sx3cal/backgains.dat");
std::string temp;
int backpos, frontpos, clkpos;
std::cout << "foo" << std::endl;
if (infile.is_open())
while(infile>>clkpos>>temp>>frontpos>>temp>>backpos>>sx3BackGain[clkpos][frontpos][backpos])
std::cout << sx3BackGain[clkpos][frontpos][backpos] << std::endl;
infile.close();
infile.open("sx3cal/frontgains.dat");
if (infile.is_open())
while(infile>>clkpos>>temp>>temp>>frontpos>>sx3FrontOffset[clkpos][frontpos]>>sx3FrontGain[clkpos][frontpos])
std::cout << sx3FrontOffset[clkpos][frontpos] << " " << sx3FrontGain[clkpos][frontpos] << std::endl;
infile.close();
infile.open("sx3cal/rightgains.dat");
if (infile.is_open())
while(infile>>clkpos>>frontpos>>temp>>sx3RightGain[clkpos][frontpos]) {
sx3RightGain[clkpos][frontpos]=TMath::Abs(sx3RightGain[clkpos][frontpos]);
}
infile.close();
}
std::cout << "aaa" << std::endl;
}
Bool_t MakeVertexSX3::Process(Long64_t entry)
{
hitPos.Clear();
qqqenergy = -1;
qqqtimestamp=-1;
HitNonZero = false;
bool qqq1000cut = false;
b_sx3Multi->GetEntry(entry);
b_sx3ID->GetEntry(entry);
b_sx3Ch->GetEntry(entry);
b_sx3E->GetEntry(entry);
b_sx3T->GetEntry(entry);
b_qqqMulti->GetEntry(entry);
b_qqqID->GetEntry(entry);
b_qqqCh->GetEntry(entry);
b_qqqE->GetEntry(entry);
b_qqqT->GetEntry(entry);
b_pcMulti->GetEntry(entry);
b_pcID->GetEntry(entry);
b_pcCh->GetEntry(entry);
b_pcE->GetEntry(entry);
b_pcT->GetEntry(entry);
sx3.CalIndex();
qqq.CalIndex();
pc.CalIndex();
std::vector<Event> sx3Events;
if(sx3.multi>1) {
std::array<sx3det,12> Fsx3;
//std::cout << "-----" << std::endl;
for(int i=0; i<sx3.multi; i++) {
if(sx3.id[i]>=12) continue;
int id = sx3.id[i];
if(sx3.ch[i]>=8) {
int sx3ch=sx3.ch[i]-8;
sx3ch=(sx3ch+3)%4;
if(sx3ch==0 || sx3ch==3) continue;
float value=sx3.e[i];
int gch = sx3.id[i]*4+(sx3.ch[i]-8);
Fsx3.at(id).fillevent("BACK",sx3ch,value);
Fsx3.at(id).ts = static_cast<double>(sx3.t[i]);
plotter->Fill2D("sx3backs_raw",100,0,100,800,0,4096,gch,sx3.e[i]);
} else {
int sx3ch=sx3.ch[i]/2;
double value=sx3.e[i];
if(sx3.ch[i]%2==0) {
Fsx3.at(id).fillevent("FRONT_L",sx3ch,value*sx3RightGain[id][sx3ch]);
} else {
Fsx3.at(id).fillevent("FRONT_R",sx3ch,value);
}
}
}
for(int id=0; id<12; id++) {
Fsx3.at(id).validate();
auto det = Fsx3.at(id);
bool no_charge_sharing_strict = det.valid_front_chans.size()==1 && det.valid_back_chans.size()==1;
if(det.valid) {
//std::cout << det.frontEL << " " << det.frontEL*sx3RightGain[id][det.stripF] << std::endl;
plotter->Fill2D("be_vs_x_sx3_id_"+std::to_string(id)+"_f"+std::to_string(det.stripF)+"_b"+std::to_string(det.stripB),200,-1,1,800,0,8192,
det.frontX,det.backE,"evsx");
//std::cout << sx3BackGain[id][det.stripF][det.stripB] << " " << sx3FrontGain[id][det.stripF] << std::endl;
plotter->Fill2D("matched_be_vs_x_sx3_id_"+std::to_string(id)+"_f"+std::to_string(det.stripF),200,-30,30,800,0,8192,
det.frontX*sx3FrontGain[id][det.stripF]+sx3FrontOffset[id][det.stripF],det.backE*sx3BackGain[id][det.stripF][det.stripB],"evsx_matched");
//plotter->Fill2D("fe_vs_x_sx3_id_"+std::to_string(id)+"_f"+std::to_string(det.stripF)+"_"+std::to_string(det.stripB),200,-1,1,800,0,4096,det.frontX,det.backE,"evsx");
plotter->Fill2D("l_vs_r_sx3_id_"+std::to_string(id)+"_f"+std::to_string(det.stripF),800,0,4096,800,0,4096,det.frontEL,det.frontER,"l_vs_r");
}
if(det.valid && (id ==9 || id==7 || id == 1 || id==3) && det.stripF!=DEFAULT_NULL && det.stripB!=DEFAULT_NULL) {
double z = det.frontX*sx3FrontGain[id][det.stripF]+sx3FrontOffset[id][det.stripF];
double backE = det.backE*sx3BackGain[id][det.stripF][det.stripB];
Event sx3ev(TVector3(0,0,z),backE,-1,det.ts,-1,det.stripB+4*id,det.stripF+4*id);
sx3Events.push_back(sx3ev);
}
}
}
//return kTRUE;
// QQQ Processing
int qqqCount = 0;
int qqqAdjCh = 0;
// REMOVE WHEN RERUNNING USING THE NEW CALIBRATION FILE
for (int i = 0; i < qqq.multi; i++)
{
//if ((qqq.id[i] == 3 || qqq.id[i] == 1) && qqq.ch[i] < 16)
if (qqq.id[i] == 1 && qqq.ch[i] < 16) //for run 12, 26Al
{
qqq.ch[i] = 16 - qqq.ch[i];
}
}
for (int i = 0; i < qqq.multi; i++)
{
if (qqq.id[i] == 0 && qqq.ch[i] >= 16)
{
qqq.ch[i] = 31 - qqq.ch[i] + 16;
}
}
std::vector<Event> QQQ_Events, PC_Events;
std::vector<Event> QQQ_Events_Raw, PC_Events_Raw;
std::vector<Event> QQQ_Events2; //clustering done
std::unordered_map<int,std::tuple<int,int,double,double>> qvecr[4], qvecw[4];
if(qqq.multi>1) {
//if(qqq.multi>=3) std::cout << "-----" << std::endl;
for(int i=0; i<qqq.multi; i++) {
//if(qqq.multi>=3) std::cout << std::setprecision(16) << "qqq"<< qqq.id[i] << " " << std::string(qqq.ch[i]/16?"ring":"wedge") << qqq.ch[i]%16 << " " << qqq.e[i] << " " << qqq.t[i] - qqq.t[0] << std::endl;
if(qqq.ch[i]/16) {
if(qvecr[qqq.id[i]].find(qqq.ch[i])!=qvecr[qqq.id[i]].end()) std::cout << "mayday!" << std::endl;
qvecr[qqq.id[i]][qqq.ch[i]] = std::tuple(qqq.id[i],qqq.ch[i],qqq.e[i],qqq.t[i]);
} else {
if(qvecw[qqq.id[i]].find(qqq.ch[i])!=qvecw[qqq.id[i]].end()) std::cout << "mayday!" << std::endl;
qvecw[qqq.id[i]][qqq.ch[i]] = std::tuple(qqq.id[i],qqq.ch[i],qqq.e[i],qqq.t[i]);
}
}
}
bool PCQQQTimeCut = false;
for (int i = 0; i < qqq.multi; i++) {
plotter->Fill2D("QQQ_Index_Vs_Energy", 16 * 8, 0, 16 * 8, 2000, 0, 16000, qqq.index[i], qqq.e[i], "hRawQQQ");
for (int j = 0; j < qqq.multi; j++) {
if (j == i)
continue;
plotter->Fill2D("QQQ_Coincidence_Matrix", 16 * 8, 0, 16 * 8, 16 * 8, 0, 16 * 8, qqq.index[i], qqq.index[j], "hRawQQQ");
}
for (int k = 0; k < pc.multi; k++) {
if (pc.index[k] < 24 && pc.e[k] > 50) {
plotter->Fill2D("QQQ_Vs_Anode_Energy", 400, 0, 4000, 1000, 0, 16000, qqq.e[i], pc.e[k], "hRawQQQ");
plotter->Fill2D("QQQ_Vs_PC_Index", 16 * 8, 0, 16 * 8, 24, 0, 24, qqq.index[i], pc.index[k], "hRawQQQ");
}
else if (pc.index[k] >= 24 && pc.e[k] > 50) {
plotter->Fill2D("QQQ_Vs_Cathode_Energy", 400, 0, 4000, 1000, 0, 16000, qqq.e[i], pc.e[k], "hRawQQQ");
}
}
for (int j = i + 1; j < qqq.multi; j++) {
if (qqq.id[i] == qqq.id[j]) {
qqqCount++;
int chWedge = -1;
int chRing = -1;
double eWedge = 0.0;
double eWedgeMeV = 0.0;
double eRing = 0.0;
double eRingMeV = 0.0;
double tRing = 0.0;
double tWedge = 0.0;
if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && qqqGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]) {
chWedge = qqq.ch[i];
eWedge = qqq.e[i] * qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
chRing = qqq.ch[j] - 16;
eRing = qqq.e[j];
tRing = static_cast<double>(qqq.t[j]);
tWedge = static_cast<double>(qqq.t[i]);
}
else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && qqqGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16]) {
chWedge = qqq.ch[j];
eWedge = qqq.e[j] * qqqGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
chRing = qqq.ch[i] - 16;
eRing = qqq.e[i];
tRing = static_cast<double>(qqq.t[i]);
tWedge = static_cast<double>(qqq.t[j]);
}
else
continue;
plotter->Fill1D("Wedgetime_Vs_Ringtime", 100, -1000, 1000, tWedge - tRing, "hTiming");
plotter->Fill2D("RingE_vs_Index", 16 * 4, 0, 16 * 4, 1000, 0, 16000, chRing + qqq.id[i] * 16, eRing, "hRawQQQ");
plotter->Fill2D("WedgeE_vs_Index", 16 * 4, 0, 16 * 4, 1000, 0, 16000, chWedge + qqq.id[i] * 16, eWedge, "hRawQQQ");
if (qqqCalibValid[qqq.id[i]][chWedge][chRing]) {
eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chWedge][chRing] / 1000;
eRingMeV = eRing * qqqCalib[qqq.id[i]][chWedge][chRing] / 1000;
if(eRingMeV/eWedgeMeV > 3.0 || eRingMeV/eWedgeMeV<1.0/3.0) continue;
//if(eRingMeV<4.0 || eWedgeMeV<4.0) continue;
double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
double rho = 50. + (50. / 16.) * (chRing + 0.5); //"?"
//z used to be 75+30+23=128
//we found a 12mm shift towards the vertex later --> 116
Event qqqevent(TVector3(rho*TMath::Cos(theta),rho*TMath::Sin(theta),116), eRingMeV, eWedgeMeV, tRing, tWedge,chRing+qqq.id[i]*16, chWedge+qqq.id[i]*16);
Event qqqeventr(TVector3(rho*TMath::Cos(theta),rho*TMath::Sin(theta),116), eRing, eWedge, tRing, tWedge,chRing+qqq.id[i]*16, chWedge+qqq.id[i]*16);
QQQ_Events.push_back(qqqevent);
QQQ_Events_Raw.push_back(qqqeventr);
plotter->Fill2D("QQQCartesianPlot", 200, -100, 100, 200, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hCalQQQ");
plotter->Fill2D("QQQCartesianPlot" + std::to_string(qqq.id[i]), 200, -100, 100, 200, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hCalQQQ");
if (PCQQQTimeCut) {
plotter->Fill2D("PC_XY_Projection_QQQ_TimeCut" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hPCQQQ");
}
plotter->Fill2D("PC_XY_Projection_QQQ" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hPCQQQ");
}
else
continue;
plotter->Fill2D("WedgeE_Vs_RingECal", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ");
plotter->Fill2D("WedgeE_Vs_RingECal_selected", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ");
for (int k = 0; k < pc.multi; k++)
{
plotter->Fill2D("RingCh_vs_Anode_Index", 16 * 4, 0, 16 * 4, 24, 0, 24, chRing + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
plotter->Fill2D("WedgeCh_vs_Anode_Index", 16 * 4, 0, 16 * 4, 24, 0, 24, chWedge + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
plotter->Fill2D("WedgeCh_vs_Anode_Index" + std::to_string(qqq.id[i]), 16 * 4, 0, 16 * 4, 24, 0, 24, chWedge + qqq.id[i] * 16, pc.index[k]);
plotter->Fill2D("RingCh_vs_Cathode_Index", 16 * 4, 0, 16 * 4, 24, 24, 48, chRing + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
plotter->Fill2D("WedgeCh_vs_Cathode_Index", 16 * 4, 0, 16 * 4, 24, 24, 48, chWedge + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
if (pc.index[k] < 24 && pc.e[k] > 50)
{
plotter->Fill2D("Timing_Difference_QQQ_PC", 500, -2000, 2000, 16, 0, 16, tRing - static_cast<double>(pc.t[k]), chRing, "hTiming");
plotter->Fill2D("DelT_Vs_QQQRingECal", 500, -2000, 2000, 1000, 0, 10, tRing - static_cast<double>(pc.t[k]), eRingMeV, "hTiming");
plotter->Fill2D("CalibratedQQQEvsPCE_R", 1000, 0, 10, 2000, 0, 30000, eRingMeV, pc.e[k], "hPCQQQ");
plotter->Fill2D("CalibratedQQQEvsPCE_W", 1000, 0, 10, 2000, 0, 30000, eWedgeMeV, pc.e[k], "hPCQQQ");
if (tRing - static_cast<double>(pc.t[k]) < -150) // proton tests, 27Al
//if (tRing - static_cast<double>(pc.t[k]) < -150 && tRing - static_cast<double>(pc.t[k]) > -450) // 27Al
//if (tRing - static_cast<double>(pc.t[k]) < -70 && tRing - static_cast<double>(pc.t[k]) > -150) // 17F
{
PCQQQTimeCut = true;
}
}
if (pc.index[k] >= 24 && pc.e[k] > 50) {
plotter->Fill2D("Timing_Difference_QQQ_PC_Cathode", 500, -2000, 2000, 16, 0, 16, tRing - static_cast<double>(pc.t[k]), chRing, "hTiming");
}
} //end of pc k loop
if (!HitNonZero) {
double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
double rho = 50. + (50. / 16.) * (chRing + 0.5); //"?"
double x = rho * TMath::Cos(theta);
double y = rho * TMath::Sin(theta);
hitPos.SetXYZ(x, y, (23 + 75 + 30));
qqqenergy = eRingMeV;
qqqtimestamp = tRing;
HitNonZero = true;
}
} // if j==i
} //j loop end
} //i loop end
plotter->Fill1D("QQQ_Multiplicity", 10, 0, 10, qqqCount, "hRawQQQ");
/*if(QQQ_Events.size()>=1) {
std::cout<< " ---->" << std::endl;
for(auto qe: QQQ_Events) {
std::cout << qe.ch1/16 << " " <<qe.ch2/16 << " " << qe.ch1%16 << " "<< qe.ch2%16 << " " << qe.Energy1 << " " << qe.Energy2 << " " << std::endl;
}
}*/
typedef std::unordered_map<int,std::tuple<int,double,double>> WireEvent; //this stores nearest neighbour wire events, or a 'cluster'
WireEvent aWireEvents, cWireEvents; //naming for book keeping
aWireEvents.clear();
aWireEvents.reserve(24);
// PC Gain Matching and Filling
double anodeT = -99999;
double cathodeT = 99999;
int anodeIndex = -1;
int cathodeIndex = -1;
for (int i = 0; i < pc.multi; i++)
{
if (pc.e[i] > 50)
{
plotter->Fill2D("PC_Index_Vs_Energy", 48, 0, 48, 2000, 0, 30000, pc.index[i], static_cast<double>(pc.e[i]), "hRawPC");
} else {
continue;
}
if (pc.index[i] < 48)
{
pc.e[i] = pcSlope[pc.index[i]] * pc.e[i] + pcIntercept[pc.index[i]];
plotter->Fill2D("PC_Index_VS_GainMatched_Energy", 48, 0, 48, 2000, 0, 30000, pc.index[i], pc.e[i], "hGMPC");
}
if (pc.index[i] < 24)
{
anodeT = static_cast<double>(pc.t[i]);
anodeIndex = pc.index[i];
aWireEvents[pc.index[i]] = std::tuple(pc.index[i],pc.e[i],static_cast<double>(pc.t[i]));
}
else
{
cathodeT = static_cast<double>(pc.t[i]);
cathodeIndex = pc.index[i] - 24;
cWireEvents[pc.index[i]-24] = std::tuple(pc.index[i]-24,pc.e[i],static_cast<double>(pc.t[i]));
}
if (anodeT != -99999 && cathodeT != 99999)
{
for (int j = 0; j < qqq.multi; j++)
{
plotter->Fill1D("PC_Time_qqq", 200, -2000, 2000, anodeT - cathodeT, "hTiming");
plotter->Fill2D("PC_Time_Vs_QQQ_ch", 200, -2000, 2000, 16 * 8, 0, 16 * 8, anodeT - cathodeT, qqq.ch[j], "hTiming");
plotter->Fill2D("PC_Time_vs_AIndex", 200, -2000, 2000, 24, 0, 24, anodeT - cathodeT, anodeIndex, "hTiming");
plotter->Fill2D("PC_Time_vs_CIndex", 200, -2000, 2000, 24, 0, 24, anodeT - cathodeT, cathodeIndex, "hTiming");
// plotter->Fill1D("PC_Time_A" + std::to_string(anodeIndex) + "_C" + std::to_string(cathodeIndex), 200, -1000, 1000, anodeT - cathodeT, "TimingPC");
}
for (int j = 0; j < sx3.multi; j++)
{
plotter->Fill1D("PC_Time_sx3", 200, -2000, 2000, anodeT - cathodeT, "hTiming");
}
plotter->Fill1D("PC_Time", 200, -2000, 2000, anodeT - cathodeT, "hTiming");
}
for (int j = i + 1; j < pc.multi; j++)
{
plotter->Fill2D("PC_Coincidence_Matrix", 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC");
plotter->Fill2D("PC_Coincidence_Matrix_anodeMinusCathode_lt_-200_" + std::to_string(anodeT - cathodeT < -200), 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC");
plotter->Fill2D("Anode_V_Anode", 24, 0, 24, 24, 0, 24, pc.index[i], pc.index[j], "hGMPC");
}
}
anodeHits.clear();
cathodeHits.clear();
corrcatMax.clear();
int aID = 0;
int cID = 0;
double aE = 0;
double cE = 0;
double aESum = 0;
double cESum = 0;
double aEMax = 0;
int aIDMax = 0;
for (int i = 0; i < pc.multi; i++) {
// if (pc.e[i] > 100)
{
if (pc.index[i] < 24) {
anodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
}
else if (pc.index[i] >= 24) {
cathodeHits.push_back(std::pair<int, double>(pc.index[i] - 24, pc.e[i]));
}
}
}
std::sort(anodeHits.begin(),anodeHits.end(),[](std::pair<int,double> a, std::pair<int,double> b){ return a.first < b.first;});
std::sort(cathodeHits.begin(),cathodeHits.end(),[](std::pair<int,double> a, std::pair<int,double> b){ return a.first < b.first;});
//clusters = collection of (collection of wires) where each wire is (index, energy, timestamp)
std::vector<std::vector<std::tuple<int,double,double>>> aClusters = pwinstance.Make_Clusters(aWireEvents);
std::vector<std::vector<std::tuple<int,double,double>>> cClusters = pwinstance.Make_Clusters(cWireEvents);
std::vector<std::pair<double,double>> sumE_AC;
for(auto aCluster: aClusters) {
for(auto cCluster: cClusters) {
if(aCluster.size()<=1 && cCluster.size()<=1) continue;
auto [crossover,alpha,apSumE,cpSumE,apMaxE,cpMaxE,apTSMaxE,cpTSMaxE] = pwinstance.FindCrossoverProperties(aCluster, cCluster);
if(alpha!=9999999 && apSumE!=-1) {
//Event PCEvent(crossover,apMaxE,cpMaxE,apTSMaxE,cpTSMaxE);
//Event PCEvent(crossover,apSumE,cpSumE,apTSMaxE,cpTSMaxE);
Event PCEvent(crossover,apSumE,cpMaxE,apTSMaxE,cpTSMaxE); //run12 shows cathode-max and anode-sum provide best dE signals.
//std::cout << apSumE << " " << crossover.Perp() << " " << apMaxE << " " << apTSMaxE << std::endl;
PC_Events.push_back(PCEvent);
sumE_AC.push_back(std::pair(apSumE,cpSumE));
}
}
}
if(QQQ_Events.size() && PC_Events.size())
plotter->Fill2D("PCEv_vs_QQQEv",20,0,20,20,0,20,QQQ_Events.size(),PC_Events.size());
for(auto pcevent:PC_Events) {
for(auto sx3event:sx3Events) {
plotter->Fill1D("dt_pcA_sx3B"+std::to_string(sx3event.ch2),640,-2000,2000,sx3event.Time1 - pcevent.Time1);
plotter->Fill1D("dt_pcC_sx3B"+std::to_string(sx3event.ch2),640,-2000,2000,sx3event.Time1 - pcevent.Time2);
plotter->Fill2D("dE_E_Anodesx3B",400,0,10,800,0,40000,sx3event.Energy1*0.001,pcevent.Energy1);
plotter->Fill2D("dE_E_Cathodesx3B",400,0,10,800,0,10000,sx3event.Energy1*0.001,pcevent.Energy2);
double sx3z = sx3event.pos.Z()+(75.0/2.0)+23.0-90.0; //w.r.t target origin at 90 for run12
double sx3rho = 88.0;
double sx3theta = TMath::ATan2(sx3rho,sx3z);
double pczguess = 40.0/TMath::Tan(sx3theta) + 90.0;
plotter->Fill2D("pcz_vs_sx3pczguess",300,0,200,150,0,200,pczguess,pcevent.pos.Z());
plotter->Fill2D("pcz_vs_sx3pczguess"+std::to_string(sx3event.ch2),300,0,200,150,0,200,pczguess,pcevent.pos.Z());
plotter->Fill2D("pcz_vs_sx3z",300,0,200,150,0,200,sx3z+90,pcevent.pos.Z());
}
}
for(auto pcevent: PC_Events) {
for(auto qqqevent: QQQ_Events) {
plotter->Fill1D("dt_pcA_qqqR",640,-2000,2000,qqqevent.Time1 - pcevent.Time1);
plotter->Fill1D("dt_pcC_qqqW",640,-2000,2000,qqqevent.Time2 - pcevent.Time2);
plotter->Fill2D("dE_E_AnodeQQQR",400,0,10,800,0,40000,qqqevent.Energy1,pcevent.Energy1);
plotter->Fill2D("dE_E_CathodeQQQR",400,0,10,800,0,10000,qqqevent.Energy2,pcevent.Energy2);
double sinTheta = TMath::Sin((qqqevent.pos - TVector3(0,0,90)).Theta())/TMath::Sin((TVector3(51.5,0,128.) - TVector3(0,0,90)).Theta());
plotter->Fill2D("dE2_E_AnodeQQQR",400,0,10,800,0,40000,qqqevent.Energy1,pcevent.Energy1*sinTheta);
plotter->Fill2D("dE2_E_CathodeQQQR",400,0,10,800,0,10000,qqqevent.Energy2,pcevent.Energy2*sinTheta);
if(qqqevent.pos.Phi() <= pcevent.pos.Phi()+TMath::Pi()/4. && qqqevent.pos.Phi() >= pcevent.pos.Phi()-TMath::Pi()/4.) {
plotter->Fill1D("PCZ",800,-200,200,pcevent.pos.Z(),"phicut");
double pcz_guess = 40.0/TMath::Tan((qqqevent.pos-TVector3(0,0,90)).Theta()) + 90; //this is ideally kept to be all QQQ+userinput for calibration of pcz
plotter->Fill2D("pczguess_vs_pc",300,0,200,150,0,200,pcz_guess,pcevent.pos.Z(),"phicut");
plotter->Fill2D("pczguess_vs_pc_phi="+std::to_string(qqqevent.pos.Phi()*180./M_PI),300,0,200,150,0,200,pcz_guess,pcevent.pos.Z(),"phicut");
//plotter->Fill1D("PCZ",800,-200,200,pcevent.pos.Z(),"phicut");
}
}
}
//HALFTIME! Can stop here in future versions
//return kTRUE;
if (anodeHits.size() >= 1 && cathodeHits.size() >= 1)
{
// 2. CRITICAL FIX: Define reference vector 'a'
// In Analyzer.cxx, 'a' was left over from the loop. We use the first anode wire as reference here.
// (Assuming pwinstance.An is populated and wires are generally parallel).
TVector3 refAnode = pwinstance.An[0].first - pwinstance.An[0].second;
{
for (const auto &anode : anodeHits)
{
aID = anode.first;
aE = anode.second;
aESum += aE;
if (aE > aEMax)
{
aEMax = aE;
aIDMax = aID;
}
}
for (const auto &cathode : cathodeHits)
{
cID = cathode.first;
cE = cathode.second;
plotter->Fill2D("AnodeMax_Vs_Cathode_Coincidence_Matrix", 24, 0, 24, 24, 0, 24, aIDMax, cID, "hRawPC");
plotter->Fill2D("Anode_Vs_Cathode_Coincidence_Matrix", 24, 0, 24, 24, 0, 24, aID, cID, "hRawPC");
plotter->Fill2D("Anode_vs_CathodeE", 2000, 0, 30000, 2000, 0, 30000, aE, cE, "hGMPC");
plotter->Fill2D("CathodeMult_V_CathodeE", 6, 0, 6, 2000, 0, 30000, cathodeHits.size(), cE, "hGMPC");
for (int j = -4; j < 3; j++)
{
if ((aIDMax + 24 + j) % 24 == 23 - cID)
{
corrcatMax.push_back(std::pair<int, double>(cID, cE));
cESum += cE;
}
}
}
}
}
TVector3 anodeIntersection,vector_closest_to_z;
anodeIntersection.Clear();
vector_closest_to_z.Clear();
if (corrcatMax.size() > 0)
{
double x = 0, y = 0, z = 0;
for (const auto &corr : corrcatMax)
{
if (Crossover[aIDMax][corr.first][0].z > 9000000)
continue;
if (cESum > 0)
{
x += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].x;
y += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].y;
z += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].z;
}
}
if (x == 0 && y == 0 && z == 0)
;
// to ignore events with no valid crossover points
else
anodeIntersection = TVector3(x, y, z);
// << "Anode Intersection: " << anodeIntersection.X() << ", " << anodeIntersection.Y() << ", " << anodeIntersection.Z() << std::endl;
}
bool PCQQQPhiCut = false;
// flip the algorithm for cathode 1 multi anode events
if ((hitPos.Phi() > (anodeIntersection.Phi() - TMath::PiOver4())) && (hitPos.Phi() < (anodeIntersection.Phi() + TMath::PiOver4()))) {
PCQQQPhiCut = true;
}
for (double Tz = 60; Tz <= 100; Tz += 1.0)
{
TVector3 TargetPos(0, 0, Tz);
if(PCQQQPhiCut && anodeIntersection.Perp()>0 && anodeIntersection.Z()!=0 && cathodeHits.size()>=2) {
plotter->Fill2D("Inttheta_vs_QQQtheta_TC" + std::to_string(PCQQQTimeCut) + "_TZ" + std::to_string(Tz), 400, 0, 180, 90, 0, 90, (anodeIntersection - TargetPos).Theta() * 180. / TMath::Pi(), (hitPos - TargetPos).Theta() * 180. / TMath::Pi(), "TPosVariation");
//plotter->Fill2D("R_ratio_to_Z_ratio" + std::to_string(PCQQQTimeCut) + "_TZ" + std::to_string(Tz), 100, -2, 2, 100, -2, 2, (anodeIntersection - TargetPos).Z()/(hitPos-TargetPos).Z(), ((anodeIntersection - TargetPos).Perp()+2.5)/(hitPos-TargetPos).Perp(), "TPosVariation");
}
}
if (anodeIntersection.Z() != 0 && anodeIntersection.Perp()>0 && HitNonZero)
{
plotter->Fill1D("PC_Z_Projection", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
plotter->Fill2D("Z_Proj_VsDelTime", 600, -300, 300, 200, -2000, 2000, anodeIntersection.Z(), anodeT - cathodeT, "hPCzQQQ");
plotter->Fill2D("IntPhi_vs_QQQphi", 100, -200, 200, 80, -200, 200, anodeIntersection.Phi() * 180. / TMath::Pi(), hitPos.Phi() * 180. / TMath::Pi(), "hPCQQQ");
//plotter->Fill2D("Inttheta_vs_QQQtheta", 90, 0, 180, 20, 0, 45, anodeIntersection.Theta() * 180. / TMath::Pi(), hitPos.Theta() * 180. / TMath::Pi(), "hPCQQQ");
//plotter->Fill2D("Inttheta_vs_QQQtheta_TC" + std::to_string(PCQQQTimeCut)+ "_PC"+std::to_string(PCQQQPhiCut), 90, 0, 180, 20, 0, 45, anodeIntersection.Theta() * 180. / TMath::Pi(), hitPos.Theta() * 180. / TMath::Pi(), "hPCQQQ");
plotter->Fill2D("IntPhi_vs_QQQphi_TC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 100, -200, 200, 80, -200, 200, anodeIntersection.Phi() * 180. / TMath::Pi(), hitPos.Phi() * 180. / TMath::Pi(), "hPCQQQ");
}
if (anodeIntersection.Z() != 0 && cathodeHits.size() >= 2)
plotter->Fill1D("PC_Z_Projection_TC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 1)
{
plotter->Fill1D("PC_Z_proj_1C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_1C", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hPCzQQQ");
}
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 2)
{
plotter->Fill1D("PC_Z_proj_2C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_2C", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hGMPC");
}
if (anodeIntersection.Z() != 0 && cathodeHits.size() > 2)
{
plotter->Fill1D("PC_Z_proj_nC", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_nC", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hGMPC");
}
if (anodeHits.size() > 0 && cathodeHits.size() > 0)
plotter->Fill2D("AHits_vs_CHits", 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
// make another plot with nearest neighbour constraint
bool hasNeighbourAnodes = false;
bool hasNeighbourCathodes = false;
// 1. Check Anodes for neighbours (including wrap-around 0-23)
for (size_t i = 0; i < anodeHits.size(); i++)
{
for (size_t j = i + 1; j < anodeHits.size(); j++)
{
int diff = std::abs(anodeHits[i].first - anodeHits[j].first);
if (diff == 1 || diff == 23)
{ // 23 handles the cylindrical wrap
hasNeighbourAnodes = true;
break;
}
}
if (hasNeighbourAnodes)
break;
}
// 2. Check Cathodes for neighbours (including wrap-around 0-23)
for (size_t i = 0; i < cathodeHits.size(); i++)
{
for (size_t j = i + 1; j < cathodeHits.size(); j++)
{
int diff = std::abs(cathodeHits[i].first - cathodeHits[j].first);
if (diff == 1 || diff == 23)
{
hasNeighbourCathodes = true;
break;
}
}
if (hasNeighbourCathodes)
break;
}
// ---------------------------------------------------------
// FILL PLOTS
// ---------------------------------------------------------
if (anodeHits.size() > 0 && cathodeHits.size() > 0)
{
plotter->Fill2D("AHits_vs_CHits_NA" + std::to_string(hasNeighbourAnodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
plotter->Fill2D("AHits_vs_CHits_NC" + std::to_string(hasNeighbourCathodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
// Constraint Plot: Only fill if BOTH planes have adjacent hits
// This effectively removes events with only isolated single-wire hits (noise)
if (hasNeighbourAnodes && hasNeighbourCathodes)
{
plotter->Fill2D("AHits_vs_CHits_NN", 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
}
}
if (HitNonZero && anodeIntersection.Z() != 0)
{
pw_contr.CalTrack2(hitPos, anodeIntersection);
plotter->Fill1D("VertexRecon", 600, -1300, 1300, pw_contr.GetZ0());
plotter->Fill1D("VertexRecon_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, pw_contr.GetZ0());
if (cathodeHits.size() == 2)
plotter->Fill1D("VertexRecon_2c_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, pw_contr.GetZ0());
TVector3 x2(anodeIntersection), x1(hitPos);
TVector3 v = x2-x1;
double t_minimum = -1.0*(x1.X()*v.X()+x1.Y()*v.Y())/(v.X()*v.X()+v.Y()*v.Y());
vector_closest_to_z = x1 + t_minimum*v;
plotter->Fill1D("VertexRecon_Z_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, vector_closest_to_z.Z() ,"customVertex");
if(vector_closest_to_z.Perp() < 20) {
plotter->Fill1D("VertexRecon_RadialCut_Z_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, vector_closest_to_z.Z() ,"customVertex");
}
plotter->Fill2D("VertexRecon_XY_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 100, -100, 100, 100,-100,100, vector_closest_to_z.X(), vector_closest_to_z.Y() ,"customVertex");
if(cathodeHits.size()==2) {
plotter->Fill1D("VertexRecon2C_Z_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, vector_closest_to_z.Z() ,"customVertex");
if(vector_closest_to_z.Perp() < 20) {
plotter->Fill1D("VertexRecon2C_RadialCut_Z_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, vector_closest_to_z.Z() ,"customVertex");
}
plotter->Fill2D("VertexRecon2C_XY_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 100, -100, 100, 100,-100,100, vector_closest_to_z.X(), vector_closest_to_z.Y() ,"customVertex");
plotter->Fill2D("VertexRecon2C_RhoZ_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 100, -100, 100, 600,-1300,1300, vector_closest_to_z.Perp(), vector_closest_to_z.Z() ,"customVertex");
plotter->Fill2D("VertexRecon2C_Z_vs_QQQE_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, 800,0,20000, vector_closest_to_z.Z(), qqqenergy ,"customVertex");
}
}
for (int i = 0; i < qqq.multi; i++)
{
if(anodeIntersection.Perp() > 0) { //suppress x,y=0,0 events
if (PCQQQTimeCut) {
plotter->Fill2D("PC_XY_Projection_QQQ_TimeCut" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, anodeIntersection.X(), anodeIntersection.Y(), "hPCQQQ");
}
plotter->Fill2D("PC_XY_Projection_QQQ" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, anodeIntersection.X(), anodeIntersection.Y(), "hPCQQQ");
}
for (int j = i + 1; j < qqq.multi; j++)
{
if (qqq.id[i] == qqq.id[j])
{
int chWedge = -1;
int chRing = -1;
double eWedge = 0.0;
double eWedgeMeV = 0.0;
double eRing = 0.0;
double eRingMeV = 0.0;
double tRing = 0.0;
int qqqID = -1;
if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && qqqGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16])
{
chWedge = qqq.ch[i];
eWedge = qqq.e[i] * qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
chRing = qqq.ch[j] - 16;
eRing = qqq.e[j];
tRing = static_cast<double>(qqq.t[j]);
qqqID = qqq.id[i];
}
else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && qqqGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16])
{
chWedge = qqq.ch[j];
eWedge = qqq.e[j] * qqqGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
chRing = qqq.ch[i] - 16;
tRing = static_cast<double>(qqq.t[i]);
eRing = qqq.e[i];
qqqID = qqq.id[i];
}
else
continue;
if (qqqCalibValid[qqq.id[i]][chRing][chWedge])
{
eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000;
eRingMeV = eRing * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000;
}
else
continue;
// if (anodeIntersection.Z() != 0)
{
plotter->Fill2D("PC_Z_vs_QQQRing", 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCzQQQ");
}
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 2)
{
plotter->Fill2D("PC_Z_vs_QQQRing_2C", 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCzQQQ");
plotter->Fill2D("PC_Z_vs_QQQRing_2C" + std::to_string(qqq.id[i]), 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCzQQQ");
plotter->Fill2D("PC_Z_vs_QQQWedge_2C", 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chWedge, "hPCzQQQ");
}
plotter->Fill2D("VertexRecon_QQQRingTC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 600, -1300, 1300, 16, 0, 16, vector_closest_to_z.Z(), chRing, "hPCQQQ");
double phi = TMath::ATan2(anodeIntersection.Y(), anodeIntersection.X()) * 180. / TMath::Pi();
plotter->Fill2D("PolarAngle_Vs_QQQWedge" + std::to_string(qqqID), 360, -360, 360, 16, 0, 16, phi, chWedge, "hPCQQQ");
// plotter->Fill2D("EdE_PC_vs_QQQ_timegate_ls1000"+std::to_string())
plotter->Fill2D("PC_Z_vs_QQQRing_Det" + std::to_string(qqqID), 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCQQQ");
//double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
//double rho = 50. + 40. / 16. * (chRing + 0.5);
for (int k = 0; k < pc.multi; k++)
{
if(pc.index[k] >= 24)
continue;
// double sinTheta = TMath::Sin((hitPos-vector_closest_to_z).Theta());
double sinTheta = TMath::Sin((anodeIntersection-TVector3(0,0,90.0)).Theta());
// double sinTheta = TMath::Sin((anodeIntersection-vector_closest_to_z).Theta());
// double sinTheta = TMath::Sin((hitPos-TVector3(0,0,30.0)).Theta());
// double sinTheta = TMath::Sin(hitPos.Theta());
if(cathodeHits.size()==2 && PCQQQPhiCut) {
plotter->Fill2D("CalibratedQQQE_RvsCPCE_TC" + std::to_string(PCQQQTimeCut) , 400, 0, 10, 400, 0, 30000, eRingMeV, pc.e[k]*sinTheta, "hPCQQQ");
plotter->Fill2D("CalibratedQQQE_WvsCPCE_TC" + std::to_string(PCQQQTimeCut) , 400, 0, 10, 400, 0, 30000, eWedgeMeV, pc.e[k]*sinTheta, "hPCQQQ");
plotter->Fill2D("CalibratedQQQE_RvsPCE_TC" + std::to_string(PCQQQTimeCut) , 400, 0, 10, 400, 0, 30000, eRingMeV, pc.e[k], "hPCQQQ");
plotter->Fill2D("CalibratedQQQE_WvsPCE_TC" + std::to_string(PCQQQTimeCut) , 400, 0, 10, 400, 0, 30000, eWedgeMeV, pc.e[k], "hPCQQQ");
plotter->Fill2D("PCQQQ_dTimevsdPhi", 200, -2000, 2000, 80, -200, 200, tRing - static_cast<double>(pc.t[k]), (hitPos.Phi()-anodeIntersection.Phi()) * 180. / TMath::Pi(), "hTiming");
}
}
}///qqq i==j case end
} //j loop end
} // qqq i loop end
TVector3 guessVertex(0,0,90.); //for run12, subtract anodeIntersection.Z() by ~74.0 seems to work
//rho=40.0 mm is halfway between the cathodes(rho=42) and anodes(rho=37)
double pcz_guess = 42.0/TMath::Tan((hitPos-guessVertex).Theta()) + guessVertex.Z(); //this is ideally kept to be all QQQ+userinput for calibration of pcz
if(PCQQQTimeCut && PCQQQPhiCut && hitPos.Perp()>0 && anodeIntersection.Perp()>0 && cathodeHits.size()>=2) {
plotter->Fill2D("pczguess_vs_qqqE",100,0,200,800,0,20,pcz_guess,qqqenergy,"pczguess");
double pczoffset=30.0;
//plotter->Fill2D("pczguess_vs_pcz_rad="+std::to_string(hitPos.Perp()),100,0,200,150,0,200,pcz_guess,anodeIntersection.Z(),"pczguess"); //entirely qqq-derived position vs entirely PC derived position
plotter->Fill2D("pczguess_vs_pcz_phi="+std::to_string(hitPos.Phi()*180./M_PI),100,0,200,150,0,200,pcz_guess,anodeIntersection.Z()+pczoffset,"pczguess"); //entirely qqq-derived position vs entirely PC derived position
plotter->Fill2D("pczguess_vs_pcz",100,0,200,150,0,200,pcz_guess,anodeIntersection.Z()+pczoffset);
plotter->Fill2D("pcz_vs_pcPhi_rad="+std::to_string(hitPos.Perp()),360,0,360,150,0,200,anodeIntersection.Phi()*180./M_PI,anodeIntersection.Z()+pczoffset,"pczguess");
}
for (int i = 0; i < sx3.multi; i++)
{
// plotting sx3 strip hits vs anode phi
if (sx3.ch[i] < 8 && anodeIntersection.Perp()>0)
plotter->Fill2D("PCPhi_vs_SX3Strip", 100, -200, 200, 8 * 24, 0, 8 * 24, anodeIntersection.Phi() * 180. / TMath::Pi(), sx3.id[i] * 8 + sx3.ch[i]);
}
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 3)
{
plotter->Fill1D("PC_Z_proj_3C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
}
plotter->Fill2D("AnodeMaxE_Vs_Cathode_Sum_Energy", 2000, 0, 30000, 2000, 0, 30000, aEMax, cESum, "hGMPC");
plotter->Fill1D("Correlated_Cathode_MaxAnode", 6, 0, 5, corrcatMax.size(), "hGMPC");
plotter->Fill2D("Correlated_Cathode_VS_MaxAnodeEnergy", 6, 0, 5, 2000, 0, 30000, corrcatMax.size(), aEMax, "hGMPC");
plotter->Fill1D("AnodeHits", 12, 0, 11, anodeHits.size(), "hGMPC");
plotter->Fill2D("AnodeMaxE_vs_AnodeHits", 12, 0, 11, 2000, 0, 30000, anodeHits.size(), aEMax, "hGMPC");
if (anodeHits.size() < 1)
{
plotter->Fill1D("NoAnodeHits_CathodeHits", 6, 0, 5, cathodeHits.size(), "hGMPC");
}
return kTRUE;
}
void MakeVertexSX3::Terminate()
{
plotter->FlushToDisk();
}

1114
MakeVertexSX3.C.backup Normal file
View File

File diff suppressed because it is too large Load Diff

132
MakeVertexSX3.h Normal file
View File

@ -0,0 +1,132 @@
#ifndef MakeVertexSX3_h
#define MakeVertexSX3_h
#include <TROOT.h>
#include <TChain.h>
#include <TApplication.h>
#include <TFile.h>
#include <TSelector.h>
#include <iomanip>
#include <vector> // Required for vectors
#include <utility> // Required for std::pair
#include "Armory/ClassDet.h"
#include "Armory/ClassPW.h" // YOU ADDED THIS (Correct! Defines Coord)
class MakeVertexSX3 : public TSelector {
public :
TTree *fChain; //!pointer to the analyzed TTree or TChain
// Declaration of leaf types
Det sx3;
Det qqq;
Det pc ;
Det misc;
ULong64_t evID;
UInt_t run;
// List of branches
TBranch *b_eventID; //!
TBranch *b_run; //!
TBranch *b_sx3Multi; //!
TBranch *b_sx3ID; //!
TBranch *b_sx3Ch; //!
TBranch *b_sx3E; //!
TBranch *b_sx3T; //!
TBranch *b_qqqMulti; //!
TBranch *b_qqqID; //!
TBranch *b_qqqCh; //!
TBranch *b_qqqE; //!
TBranch *b_qqqT; //!
TBranch *b_pcMulti; //!
TBranch *b_pcID; //!
TBranch *b_pcCh; //!
TBranch *b_pcE; //!
TBranch *b_pcT; //!
TBranch *b_miscMulti; //!
TBranch *b_miscID; //!
TBranch *b_miscCh; //!
TBranch *b_miscE; //!
TBranch *b_miscT; //!
TBranch *b_miscTf; //!
// 1. Geometry Cache
Coord Crossover[24][24][2];
// 2. Persistent Vectors (REQUIRED for the optimized .cxx to work)
std::vector<std::pair<int, double>> anodeHits;
std::vector<std::pair<int, double>> cathodeHits;
std::vector<std::pair<int, double>> corrcatMax;
std::vector<std::pair<int, double>> corranoMax;
std::vector<double> cathodeTimes;
std::vector<double> anodeTimes;
MakeVertexSX3(TTree * /*tree*/ =0) : fChain(0) { }
virtual ~MakeVertexSX3() { }
virtual Int_t Version() const { return 2; }
virtual void Begin(TTree *tree);
virtual void SlaveBegin(TTree *tree);
virtual void Init(TTree *tree);
virtual Bool_t Notify();
virtual Bool_t Process(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
virtual void SetOption(const char *option) { fOption = option; }
virtual void SetObject(TObject *obj) { fObject = obj; }
virtual void SetInputList(TList *input) { fInput = input; }
virtual TList *GetOutputList() const { return fOutput; }
virtual void SlaveTerminate();
virtual void Terminate();
ClassDef(MakeVertexSX3,0);
};
#endif
#ifdef MakeVertexSX3_cxx
void MakeVertexSX3::Init(TTree *tree){
if (!tree) return;
fChain = tree;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("evID", &evID, &b_eventID);
fChain->SetBranchAddress("run", &run, &b_run);
sx3.SetDetDimension(24,12);
qqq.SetDetDimension(4,32);
pc.SetDetDimension(2,24);
fChain->SetBranchAddress("sx3Multi", &sx3.multi, &b_sx3Multi);
fChain->SetBranchAddress("sx3ID", &sx3.id, &b_sx3ID);
fChain->SetBranchAddress("sx3Ch", &sx3.ch, &b_sx3Ch);
fChain->SetBranchAddress("sx3E", &sx3.e, &b_sx3E);
fChain->SetBranchAddress("sx3T", &sx3.t, &b_sx3T);
fChain->SetBranchAddress("qqqMulti", &qqq.multi, &b_qqqMulti);
fChain->SetBranchAddress("qqqID", &qqq.id, &b_qqqID);
fChain->SetBranchAddress("qqqCh", &qqq.ch, &b_qqqCh);
fChain->SetBranchAddress("qqqE", &qqq.e, &b_qqqE);
fChain->SetBranchAddress("qqqT", &qqq.t, &b_qqqT);
fChain->SetBranchAddress("pcMulti", &pc.multi, &b_pcMulti);
fChain->SetBranchAddress("pcID", &pc.id, &b_pcID);
fChain->SetBranchAddress("pcCh", &pc.ch, &b_pcCh);
fChain->SetBranchAddress("pcE", &pc.e, &b_pcE);
fChain->SetBranchAddress("pcT", &pc.t, &b_pcT);
fChain->SetBranchAddress("miscMulti", &misc.multi, &b_miscMulti);
fChain->SetBranchAddress("miscID", &misc.id, &b_miscID);
fChain->SetBranchAddress("miscCh", &misc.ch, &b_miscCh);
fChain->SetBranchAddress("miscE", &misc.e, &b_miscE);
fChain->SetBranchAddress("miscT", &misc.t, &b_miscT);
}
Bool_t MakeVertexSX3::Notify(){
return kTRUE;
}
void MakeVertexSX3::SlaveBegin(TTree * /*tree*/){
// TString option = GetOption();
}
void MakeVertexSX3::SlaveTerminate(){
}
#endif // #ifdef MakeVertexSX3_cxx

12
ProcessRun.sh
View File

@ -7,22 +7,22 @@ if [ "$#" -ne 3 ]; then
exit 1 exit 1
fi fi
runID=$(printf "%03d" $1) runID=$1
timeWindow=$2 timeWindow=$2
option=$3 option=$3
# rawFolder=/home/tandem/data1/2024_09_17Fap/data # rawFolder=/home/tandem/data1/2024_09_17Fap/data
rawFolder=/mnt/d/17F rawFolder=../Raw_data
rootFolder=/mnt/d/Remapped_files/17F_data/root_data rootFolder=../root_data
if [ $option -eq 0 ]; then if [ $option -eq 0 ]; then
# rsync -auh --info=progress2 splitpole@128.186.111.223:/media/nvmeData/2024_09_17Fap/*.fsu /home/tandem/data1/2024_09_17Fap/data # rsync -auh --info=progress2 splitpole@128.186.111.223:/media/nvmeData/2024_09_17Fap/*.fsu /home/tandem/data1/2024_09_17Fap/data
fileList=`\ls -1 ${rawFolder}/*SourceRun_${runID}_*.fsu` fileList=`\ls -1 ${rawFolder}/*Run_${runID}_*.fsu`
./EventBuilder ${timeWindow} 0 0 100000000 ${fileList} # ./EventBuilder ${timeWindow} 0 0 100000000 ${fileList}
outFile=${rawFolder}/*${runID}*${timeWindow}.root outFile=${rawFolder}/*${runID}*${timeWindow}.root
@ -31,4 +31,4 @@ if [ $option -eq 0 ]; then
./Mapper ${rootFolder}/*${runID}*${timeWindow}.root ./Mapper ${rootFolder}/*${runID}*${timeWindow}.root
fi fi
# root "processRun.C(\"${rootFolder}/Run_${runID}_mapped.root\")" root "processRun.C(\"${rootFolder}/Run_${runID}_mapped.root\")"

742
TrackRecon.C.backup Normal file
View File

@ -0,0 +1,742 @@
#define TrackRecon_cxx
#include "TrackRecon.h"
#include "Armory/ClassPW.h"
#include "Armory/HistPlotter.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TBranch.h>
#include "TVector3.h"
#include <fstream>
#include <iostream>
#include <sstream>
#include <map>
#include <utility>
#include <algorithm>
// Global instances
PW pw_contr;
PW pwinstance;
TVector3 hitPos;
struct Event {
TVector3 pos;
double Energy1=-1; //Front for QQQ, Anode for PC
double Energy2=-1; //Back for QQQ, Cathode for PC
double Time1=-1;
double Time2=-1;
};
// Calibration globals
const int MAX_QQQ = 4;
const int MAX_RING = 16;
const int MAX_WEDGE = 16;
double qqqGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
bool qqqGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
double qqqCalib[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
bool qqqCalibValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
// TCutg *cutQQQ;
// PC Arrays
double pcSlope[48];
double pcIntercept[48];
HistPlotter *plotter;
bool HitNonZero;
bool sx3ecut;
bool qqqEcut;
void TrackRecon::Begin(TTree * /*tree*/)
{
TString option = GetOption();
plotter = new HistPlotter("Analyzer_QQQ.root", "TFILE");
pw_contr.ConstructGeo();
pwinstance.ConstructGeo();
// ---------------------------------------------------------
// 1. CRITICAL FIX: Initialize PC Arrays to Default (Raw)
// ---------------------------------------------------------
for (int i = 0; i < 48; i++)
{
pcSlope[i] = 1.0; // Default slope = 1 (preserves Raw energy)
pcIntercept[i] = 0.0; // Default intercept = 0
}
// Calculate Crossover Geometry ONCE
TVector3 a, c, diff;
double a2, ac, c2, adiff, cdiff, denom, alpha;
for (size_t i = 0; i < pwinstance.An.size(); i++)
{
a = pwinstance.An[i].first - pwinstance.An[i].second;
for (size_t j = 0; j < pwinstance.Ca.size(); j++)
{
c = pwinstance.Ca[j].first - pwinstance.Ca[j].second;
diff = pwinstance.An[i].first - pwinstance.Ca[j].first;
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 = pwinstance.An[i].first.X() + alpha * a.X();
Crossover[i][j][0].y = pwinstance.An[i].first.Y() + alpha * a.Y();
Crossover[i][j][0].z = pwinstance.An[i].first.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;
}
}
// Load PC Calibrations
std::ifstream inputFile("slope_intercept_results.txt");
if (inputFile.is_open())
{
std::string line;
int index;
double slope, intercept;
while (std::getline(inputFile, line))
{
std::stringstream ss(line);
ss >> index >> slope >> intercept;
if (index >= 0 && index <= 47)
{
pcSlope[index] = slope;
pcIntercept[index] = intercept;
}
}
inputFile.close();
}
else
{
std::cerr << "Error opening slope_intercept.txt" << std::endl;
}
// Load QQQ Cuts from file
// {
// std::string filename = "QQQ_PCCut.root";
// TFile *cutFile = TFile::Open(filename.c_str(), "READ");
// if (cutFile && !cutFile->IsZombie())
// {
// cutQQQ = (TCutg *)cutFile->Get("cutQQQPC");
// if (cutQQQ)
// {
// std::cout << "Loaded QQQ PC cut from " << filename << std::endl;
// }
// else
// {
// std::cerr << "Error: cutQQQPC not found in " << filename << std::endl;
// }
// cutFile->Close();
// }
// }
// ... (Load QQQ Gains and Calibs - same as before) ...
{
std::string filename = "qqq_GainMatch.dat";
std::ifstream infile(filename);
if (infile.is_open())
{
int det, ring, wedge;
double gainw, gainr;
while (infile >> det >> wedge >> ring >> gainw >> gainr)
{
qqqGain[det][wedge][ring] = gainw;
qqqGainValid[det][wedge][ring] = (gainw > 0);
// std::cout << "QQQ Gain Loaded: Det " << det << " Ring " << ring << " Wedge " << wedge << " GainW " << gainw << " GainR " << gainr << std::endl;
}
infile.close();
}
}
{
std::string filename = "qqq_Calib.dat";
std::ifstream infile(filename);
if (infile.is_open())
{
int det, ring, wedge;
double slope;
while (infile >> det >> wedge >> ring >> slope)
{
qqqCalib[det][wedge][ring] = slope;
qqqCalibValid[det][wedge][ring] = (slope > 0);
// std::cout << "QQQ Calib Loaded: Det " << det << " Ring " << ring << " Wedge " << wedge << " Slope " << slope << std::endl;
}
infile.close();
}
}
}
Bool_t TrackRecon::Process(Long64_t entry)
{
hitPos.Clear();
HitNonZero = false;
bool qqq1000cut = false;
b_sx3Multi->GetEntry(entry);
b_sx3ID->GetEntry(entry);
b_sx3Ch->GetEntry(entry);
b_sx3E->GetEntry(entry);
b_sx3T->GetEntry(entry);
b_qqqMulti->GetEntry(entry);
b_qqqID->GetEntry(entry);
b_qqqCh->GetEntry(entry);
b_qqqE->GetEntry(entry);
b_qqqT->GetEntry(entry);
b_pcMulti->GetEntry(entry);
b_pcID->GetEntry(entry);
b_pcCh->GetEntry(entry);
b_pcE->GetEntry(entry);
b_pcT->GetEntry(entry);
sx3.CalIndex();
qqq.CalIndex();
pc.CalIndex();
// QQQ Processing
int qqqCount = 0;
int qqqAdjCh = 0;
// REMOVE WHEN RERUNNING USING THE NEW CALIBRATION FILE
for (int i = 0; i < qqq.multi; i++)
{
if ((qqq.id[i] == 3 || qqq.id[i] == 1) && qqq.ch[i] < 16)
{
qqq.ch[i] = 16 - qqq.ch[i];
}
}
for (int i = 0; i < qqq.multi; i++)
{
if (qqq.id[i] == 0 && qqq.ch[i] >= 16)
{
qqq.ch[i] = 31 - qqq.ch[i] + 16;
}
}
std::vector<Event> QQQ_Events, PC_Events;
std::vector<Event> QQQ_Events_Raw, PC_Events_Raw;
bool PCQQQTimeCut = false;
for (int i = 0; i < qqq.multi; i++)
{
plotter->Fill2D("QQQ_Index_Vs_Energy", 16 * 8, 0, 16 * 8, 2000, 0, 16000, qqq.index[i], qqq.e[i], "hRawQQQ");
for (int j = 0; j < qqq.multi; j++)
{
if (j == i)
continue;
plotter->Fill2D("QQQ_Coincidence_Matrix", 16 * 8, 0, 16 * 8, 16 * 8, 0, 16 * 8, qqq.index[i], qqq.index[j], "hRawQQQ");
}
for (int k = 0; k < pc.multi; k++)
{
if (pc.index[k] < 24 && pc.e[k] > 50)
{
plotter->Fill2D("QQQ_Vs_Anode_Energy", 400, 0, 4000, 1000, 0, 16000, qqq.e[i], pc.e[k], "hRawQQQ");
plotter->Fill2D("QQQ_Vs_PC_Index", 16 * 8, 0, 16 * 8, 24, 0, 24, qqq.index[i], pc.index[k], "hRawQQQ");
}
else if (pc.index[k] >= 24 && pc.e[k] > 50)
{
plotter->Fill2D("QQQ_Vs_Cathode_Energy", 400, 0, 4000, 1000, 0, 16000, qqq.e[i], pc.e[k], "hRawQQQ");
}
}
for (int j = i + 1; j < qqq.multi; j++)
{
if (qqq.id[i] == qqq.id[j])
{
qqqCount++;
int chWedge = -1;
int chRing = -1;
double eWedge = 0.0;
double eWedgeMeV = 0.0;
double eRing = 0.0;
double eRingMeV = 0.0;
double tRing = 0.0;
double tWedge = 0.0;
if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && qqqGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16])
{
chWedge = qqq.ch[i];
eWedge = qqq.e[i] * qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
chRing = qqq.ch[j] - 16;
eRing = qqq.e[j];
tRing = static_cast<double>(qqq.t[j]);
tWedge = static_cast<double>(qqq.t[i]);
}
else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && qqqGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16])
{
chWedge = qqq.ch[j];
eWedge = qqq.e[j] * qqqGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
chRing = qqq.ch[i] - 16;
eRing = qqq.e[i];
tRing = static_cast<double>(qqq.t[i]);
tWedge = static_cast<double>(qqq.t[j]);
}
else
continue;
plotter->Fill1D("Wedgetime_Vs_Ringtime", 100, -1000, 1000, tWedge - tRing, "hTiming");
plotter->Fill2D("RingE_vs_Index", 16 * 4, 0, 16 * 4, 1000, 0, 16000, chRing + qqq.id[i] * 16, eRing, "hRawQQQ");
plotter->Fill2D("WedgeE_vs_Index", 16 * 4, 0, 16 * 4, 1000, 0, 16000, chWedge + qqq.id[i] * 16, eWedge, "hRawQQQ");
if (qqqCalibValid[qqq.id[i]][chWedge][chRing])
{
eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chWedge][chRing] / 1000;
eRingMeV = eRing * qqqCalib[qqq.id[i]][chWedge][chRing] / 1000;
}
else
continue;
plotter->Fill2D("WedgeE_Vs_RingECal", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ");
for (int k = 0; k < pc.multi; k++)
{
plotter->Fill2D("RingCh_vs_Anode_Index", 16 * 4, 0, 16 * 4, 24, 0, 24, chRing + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
plotter->Fill2D("WedgeCh_vs_Anode_Index", 16 * 4, 0, 16 * 4, 24, 0, 24, chWedge + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
plotter->Fill2D("WedgeCh_vs_Anode_Index" + std::to_string(qqq.id[i]), 16 * 4, 0, 16 * 4, 24, 0, 24, chWedge + qqq.id[i] * 16, pc.index[k]);
plotter->Fill2D("RingCh_vs_Cathode_Index", 16 * 4, 0, 16 * 4, 24, 24, 48, chRing + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
plotter->Fill2D("WedgeCh_vs_Cathode_Index", 16 * 4, 0, 16 * 4, 24, 24, 48, chWedge + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
if (pc.index[k] < 24 && pc.e[k] > 50)
{
// plotter->Fill2D("QQQ_CalibW_Vs_PC_Energy", 1000, 0, 16, 2000, 0, 30000, eWedgeMeV, pc.e[k], "hCalQQQ");
// plotter->Fill2D("QQQ_CalibR_Vs_PC_Energy", 1000, 0, 16, 2000, 0, 30000, eRingMeV, pc.e[k], "hCalQQQ");
// if (tRing - static_cast<double>(pc.t[k]) < 0 && tRing - static_cast<double>(pc.t[k]) > -600)
// // {
// // plotter->Fill2D("QQQ_CalibW_Vs_PC_Energy_Tight", 1000, 0, 16, 2000, 0, 30000, eWedgeMeV, pc.e[k], "hCalQQQ");
// // plotter->Fill2D("QQQ_CalibR_Vs_PC_Energy_Tight", 1000, 0, 16, 2000, 0, 30000, eRingMeV, pc.e[k], "hCalQQQ");
// // }
// // else
// // {
// // plotter->Fill2D("QQQ_CalibW_Vs_PC_Energy_OffTime", 1000, 0, 16, 2000, 0, 30000, eWedgeMeV, pc.e[k], "hCalQQQ");
// // plotter->Fill2D("QQQ_CalibR_Vs_PC_Energy_OffTime", 1000, 0, 16, 2000, 0, 30000, eRingMeV, pc.e[k], "hCalQQQ");
// // }
plotter->Fill2D("Timing_Difference_QQQ_PC", 500, -1000, 1000, 16, 0, 16, tRing - static_cast<double>(pc.t[k]), chRing, "hTiming");
plotter->Fill2D("DelT_Vs_QQQRingECal", 500, -1000, 1000, 1000, 0, 10, tRing - static_cast<double>(pc.t[k]), eRingMeV, "hTiming");
plotter->Fill2D("CalibratedQQQEvsPCE_R", 1000, 0, 10, 2000, 0, 30000, eRingMeV, pc.e[k], "hPCQQQ");
plotter->Fill2D("CalibratedQQQEvsPCE_W", 1000, 0, 10, 2000, 0, 30000, eWedgeMeV, pc.e[k], "hPCQQQ");
if (tRing - static_cast<double>(pc.t[k]) < -150 && tRing - static_cast<double>(pc.t[k]) > -450) // 27Al
//if (tRing - static_cast<double>(pc.t[k]) < -70 && tRing - static_cast<double>(pc.t[k]) > -150) // 17F
{
PCQQQTimeCut = true;
}
}
if (pc.index[k] >= 24 && pc.e[k] > 50) {
plotter->Fill2D("Timing_Difference_QQQ_PC_Cathode", 500, -1000, 1000, 16, 0, 16, tRing - static_cast<double>(pc.t[k]), chRing, "hTiming");
}
} //end of pc loop
double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
double rho = 50. + 40. / 16. * (chRing + 0.5);
//Event qqqevent(TVector3(rho*TMath::Cos(theta),rho*TMath::Sin(theta),23+75+30), eRingMeV, eWedgeMeV, tRing, tWedge);
//Event qqqeventr(TVector3(rho*TMath::Cos(theta),rho*TMath::Sin(theta),23+75+30), eRing, eWedge, tRing, tWedge);
//QQQ_Events.push_back(qqqevent);
//QQQ_Events_Raw.push_back(qqqeventr);
plotter->Fill2D("QQQPolarPlot", 16 * 4, -TMath::Pi(), TMath::Pi(), 32, 40, 100, theta, rho, "hCalQQQ");
plotter->Fill2D("QQQCartesianPlot", 200, -100, 100, 200, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hCalQQQ");
plotter->Fill2D("QQQCartesianPlot" + std::to_string(qqq.id[i]), 200, -100, 100, 200, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hCalQQQ");
if (PCQQQTimeCut)
{
plotter->Fill2D("PC_XY_Projection_QQQ_TimeCut" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hPCQQQ");
}
plotter->Fill2D("PC_XY_Projection_QQQ" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hPCQQQ");
if (!HitNonZero)
{
double x = rho * TMath::Cos(theta);
double y = rho * TMath::Sin(theta);
hitPos.SetXYZ(x, y, 23 + 75 + 30);
HitNonZero = true;
}
}
}
}
plotter->Fill1D("QQQ_Multiplicity", 10, 0, 10, qqqCount, "hRawQQQ");
// PC Gain Matching and Filling
double anodeT = -99999;
double cathodeT = 99999;
int anodeIndex = -1;
int cathodeIndex = -1;
for (int i = 0; i < pc.multi; i++)
{
if (pc.e[i] > 10)
{
plotter->Fill2D("PC_Index_Vs_Energy", 48, 0, 48, 2000, 0, 30000, pc.index[i], static_cast<double>(pc.e[i]), "hRawPC");
}
if (pc.index[i] < 48)
{
pc.e[i] = pcSlope[pc.index[i]] * pc.e[i] + pcIntercept[pc.index[i]];
plotter->Fill2D("PC_Index_VS_GainMatched_Energy", 48, 0, 48, 2000, 0, 30000, pc.index[i], pc.e[i], "hGMPC");
}
if (pc.index[i] < 24)
{
anodeT = static_cast<double>(pc.t[i]);
anodeIndex = pc.index[i];
}
else
{
cathodeT = static_cast<double>(pc.t[i]);
cathodeIndex = pc.index[i] - 24;
}
if (anodeT != -99999 && cathodeT != 99999)
{
for (int j = 0; j < qqq.multi; j++)
{
plotter->Fill1D("PC_Time_qqq", 200, -2000, 2000, anodeT - cathodeT, "hTiming");
plotter->Fill2D("PC_Time_Vs_QQQ_ch", 200, -2000, 2000, 16 * 8, 0, 16 * 8, anodeT - cathodeT, qqq.ch[j], "hTiming");
plotter->Fill2D("PC_Time_vs_AIndex", 200, -2000, 2000, 24, 0, 24, anodeT - cathodeT, anodeIndex, "hTiming");
plotter->Fill2D("PC_Time_vs_CIndex", 200, -2000, 2000, 24, 0, 24, anodeT - cathodeT, cathodeIndex, "hTiming");
// plotter->Fill1D("PC_Time_A" + std::to_string(anodeIndex) + "_C" + std::to_string(cathodeIndex), 200, -1000, 1000, anodeT - cathodeT, "TimingPC");
}
for (int j = 0; j < sx3.multi; j++)
{
plotter->Fill1D("PC_Time_sx3", 200, -2000, 2000, anodeT - cathodeT, "hTiming");
}
plotter->Fill1D("PC_Time", 200, -2000, 2000, anodeT - cathodeT, "hTiming");
}
for (int j = i + 1; j < pc.multi; j++)
{
plotter->Fill2D("PC_Coincidence_Matrix", 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC");
plotter->Fill2D("PC_Coincidence_Matrix_anodeMinusCathode_lt_-200_" + std::to_string(anodeT - cathodeT < -200), 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC");
plotter->Fill2D("Anode_V_Anode", 24, 0, 24, 24, 0, 24, pc.index[i], pc.index[j], "hGMPC");
}
}
anodeHits.clear();
cathodeHits.clear();
corrcatMax.clear();
int aID = 0;
int cID = 0;
double aE = 0;
double cE = 0;
double aESum = 0;
double cESum = 0;
double aEMax = 0;
int aIDMax = 0;
for (int i = 0; i < pc.multi; i++)
{
// if (pc.e[i] > 100)
{
if (pc.index[i] < 24)
anodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
else if (pc.index[i] >= 24)
cathodeHits.push_back(std::pair<int, double>(pc.index[i] - 24, pc.e[i]));
}
}
// std::sort(anodeHits.begin(), anodeHits.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b)
// { return a.second > b.second; });
// std::sort(cathodeHits.begin(), cathodeHits.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b)
// { return a.second > b.second; });
if (anodeHits.size() >= 1 && cathodeHits.size() >= 1)
{
// 2. CRITICAL FIX: Define reference vector 'a'
// In Analyzer.cxx, 'a' was left over from the loop. We use the first anode wire as reference here.
// (Assuming pwinstance.An is populated and wires are generally parallel).
TVector3 refAnode = pwinstance.An[0].first - pwinstance.An[0].second;
{
for (const auto &anode : anodeHits)
{
aID = anode.first;
aE = anode.second;
aESum += aE;
if (aE > aEMax)
{
aEMax = aE;
aIDMax = aID;
}
}
for (const auto &cathode : cathodeHits)
{
cID = cathode.first;
cE = cathode.second;
plotter->Fill2D("AnodeMax_Vs_Cathode_Coincidence_Matrix", 24, 0, 24, 24, 0, 24, aIDMax, cID, "hRawPC");
plotter->Fill2D("Anode_Vs_Cathode_Coincidence_Matrix", 24, 0, 24, 24, 0, 24, aID, cID, "hRawPC");
plotter->Fill2D("Anode_vs_CathodeE", 2000, 0, 30000, 2000, 0, 30000, aE, cE, "hGMPC");
plotter->Fill2D("CathodeMult_V_CathodeE", 6, 0, 6, 2000, 0, 30000, cathodeHits.size(), cE, "hGMPC");
for (int j = -4; j < 3; j++)
{
if ((aIDMax + 24 + j) % 24 == 23 - cID)
{
corrcatMax.push_back(std::pair<int, double>(cID, cE));
cESum += cE;
}
}
}
}
}
TVector3 anodeIntersection;
anodeIntersection.Clear();
if (corrcatMax.size() > 0)
{
double x = 0, y = 0, z = 0;
for (const auto &corr : corrcatMax)
{
if (Crossover[aIDMax][corr.first][0].z > 9000000)
continue;
if (cESum > 0)
{
x += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].x;
y += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].y;
z += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].z;
}
}
if (x == 0 && y == 0 && z == 0)
;
// to ignore events with no valid crossover points
else
anodeIntersection = TVector3(x, -y, -z);
// std::cout << "Anode Intersection: " << anodeIntersection.X() << ", " << anodeIntersection.Y() << ", " << anodeIntersection.Z() << std::endl;
}
bool PCQQQPhiCut = false;
// flip the algorithm for cathode 1 multi anode events
if ((hitPos.Phi() > (anodeIntersection.Phi() - TMath::PiOver4())) && (hitPos.Phi() < (anodeIntersection.Phi() + TMath::PiOver4()))) {
PCQQQPhiCut = true;
}
// for (double Tz = -190; Tz <= 190; Tz += 10.0)
// {
// TVector3 TargetPos(0, 0, Tz);
// plotter->Fill2D("Inttheta_vs_QQQtheta_TC" + std::to_string(PCQQQTimeCut) + "_TZ" + std::to_string(Tz), 90, 0, 180, 120, 0, 180, (anodeIntersection - TargetPos).Theta() * 180. / TMath::Pi(), (hitPos - TargetPos).Theta() * 180. / TMath::Pi(), "TPosVariation");
// }
if (anodeIntersection.Z() != 0)
{
plotter->Fill1D("PC_Z_Projection", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
plotter->Fill2D("Z_Proj_VsDelTime", 600, -300, 300, 200, -2000, 2000, anodeIntersection.Z(), anodeT - cathodeT, "hPCzQQQ");
plotter->Fill2D("IntPhi_vs_QQQphi", 100, -200, 200, 80, -200, 200, anodeIntersection.Phi() * 180. / TMath::Pi(), hitPos.Phi() * 180. / TMath::Pi(), "hPCQQQ");
plotter->Fill2D("Inttheta_vs_QQQtheta", 90, 0, 180, 20, 0, 45, anodeIntersection.Theta() * 180. / TMath::Pi(), hitPos.Theta() * 180. / TMath::Pi(), "hPCQQQ");
plotter->Fill2D("Inttheta_vs_QQQtheta_TC" + std::to_string(PCQQQTimeCut), 90, 0, 180, 20, 0, 45, anodeIntersection.Theta() * 180. / TMath::Pi(), hitPos.Theta() * 180. / TMath::Pi(), "hPCQQQ");
plotter->Fill2D("IntPhi_vs_QQQphi_TC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 100, -200, 200, 80, -200, 200, anodeIntersection.Phi() * 180. / TMath::Pi(), hitPos.Phi() * 180. / TMath::Pi(), "hPCQQQ");
}
if (anodeIntersection.Z() != 0 && cathodeHits.size() >= 2)
plotter->Fill1D("PC_Z_Projection_TC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 1)
{
plotter->Fill1D("PC_Z_proj_1C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_1C", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hPCzQQQ");
}
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 2)
{
plotter->Fill1D("PC_Z_proj_2C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_2C", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hGMPC");
}
if (anodeIntersection.Z() != 0 && cathodeHits.size() > 2)
{
plotter->Fill1D("PC_Z_proj_nC", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_nC", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hGMPC");
}
if (anodeHits.size() > 0 && cathodeHits.size() > 0)
plotter->Fill2D("AHits_vs_CHits", 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
// make another plot with nearest neighbour constraint
bool hasNeighbourAnodes = false;
bool hasNeighbourCathodes = false;
// 1. Check Anodes for neighbours (including wrap-around 0-23)
for (size_t i = 0; i < anodeHits.size(); i++)
{
for (size_t j = i + 1; j < anodeHits.size(); j++)
{
int diff = std::abs(anodeHits[i].first - anodeHits[j].first);
if (diff == 1 || diff == 23)
{ // 23 handles the cylindrical wrap
hasNeighbourAnodes = true;
break;
}
}
if (hasNeighbourAnodes)
break;
}
// 2. Check Cathodes for neighbours (including wrap-around 0-23)
for (size_t i = 0; i < cathodeHits.size(); i++)
{
for (size_t j = i + 1; j < cathodeHits.size(); j++)
{
int diff = std::abs(cathodeHits[i].first - cathodeHits[j].first);
if (diff == 1 || diff == 23)
{
hasNeighbourCathodes = true;
break;
}
}
if (hasNeighbourCathodes)
break;
}
// ---------------------------------------------------------
// FILL PLOTS
// ---------------------------------------------------------
if (anodeHits.size() > 0 && cathodeHits.size() > 0)
{
plotter->Fill2D("AHits_vs_CHits_NA" + std::to_string(hasNeighbourAnodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
plotter->Fill2D("AHits_vs_CHits_NC" + std::to_string(hasNeighbourCathodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
// Constraint Plot: Only fill if BOTH planes have adjacent hits
// This effectively removes events with only isolated single-wire hits (noise)
if (hasNeighbourAnodes && hasNeighbourCathodes)
{
plotter->Fill2D("AHits_vs_CHits_NN", 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
}
}
if (HitNonZero && anodeIntersection.Z() != 0)
{
pw_contr.CalTrack2(hitPos, anodeIntersection);
plotter->Fill1D("VertexRecon", 600, -1300, 1300, pw_contr.GetZ0());
plotter->Fill1D("VertexRecon_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -300, 300, pw_contr.GetZ0());
if (cathodeHits.size() == 2)
plotter->Fill1D("VertexRecon_2c_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -300, 300, pw_contr.GetZ0());
}
for (int i = 0; i < qqq.multi; i++)
{
if (PCQQQTimeCut) {
plotter->Fill2D("PC_XY_Projection_QQQ_TimeCut" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, anodeIntersection.X(), anodeIntersection.Y(), "hPCQQQ");
}
plotter->Fill2D("PC_XY_Projection_QQQ" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, anodeIntersection.X(), anodeIntersection.Y(), "hPCQQQ");
for (int j = i + 1; j < qqq.multi; j++)
{
if (qqq.id[i] == qqq.id[j])
{
int chWedge = -1;
int chRing = -1;
double eWedge = 0.0;
double eWedgeMeV = 0.0;
double eRing = 0.0;
double eRingMeV = 0.0;
double tRing = 0.0;
int qqqID = -1;
if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && qqqGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16])
{
chWedge = qqq.ch[i];
eWedge = qqq.e[i] * qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
chRing = qqq.ch[j] - 16;
eRing = qqq.e[j];
tRing = static_cast<double>(qqq.t[j]);
qqqID = qqq.id[i];
}
else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && qqqGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16])
{
chWedge = qqq.ch[j];
eWedge = qqq.e[j] * qqqGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
chRing = qqq.ch[i] - 16;
tRing = static_cast<double>(qqq.t[i]);
eRing = qqq.e[i];
qqqID = qqq.id[i];
}
else
continue;
if (qqqCalibValid[qqq.id[i]][chRing][chWedge])
{
eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000;
eRingMeV = eRing * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000;
}
else
continue;
// if (anodeIntersection.Z() != 0)
{
plotter->Fill2D("PC_Z_vs_QQQRing", 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCzQQQ");
}
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 2)
{
plotter->Fill2D("PC_Z_vs_QQQRing_2C", 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCzQQQ");
plotter->Fill2D("PC_Z_vs_QQQRing_2C" + std::to_string(qqq.id[i]), 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCzQQQ");
plotter->Fill2D("PC_Z_vs_QQQWedge_2C", 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chWedge, "hPCzQQQ");
}
plotter->Fill2D("Vertex_V_QQQRingTC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 600, -300, 300, 16, 0, 16, pw_contr.GetZ0(), chRing, "hPCQQQ");
double phi = TMath::ATan2(anodeIntersection.Y(), anodeIntersection.X()) * 180. / TMath::Pi();
plotter->Fill2D("PolarAngle_Vs_QQQWedge" + std::to_string(qqqID), 360, -360, 360, 16, 0, 16, phi, chWedge, "hPCQQQ");
// plotter->Fill2D("EdE_PC_vs_QQQ_timegate_ls1000"+std::to_string())
plotter->Fill2D("PC_Z_vs_QQQRing_Det" + std::to_string(qqqID), 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCQQQ");
//double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
//double rho = 50. + 40. / 16. * (chRing + 0.5);
for (int k = 0; k < pc.multi; k++)
{
if(pc.index[k] >= 24)
continue;
double sinTheta = TMath::Sin(hitPos.Theta());
plotter->Fill2D("CalibratedQQQE_RvsPCE_TC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 1000, 0, 10, 2000, 0, 30000, eRingMeV, pc.e[k]*sinTheta, "hPCQQQ");
plotter->Fill2D("CalibratedQQQE_WvsPCE_TC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 1000, 0, 10, 2000, 0, 30000, eWedgeMeV, pc.e[k]*sinTheta, "hPCQQQ");
plotter->Fill2D("PCQQQ_dTimevsdPhi", 200, -2000, 2000, 80, -200, 200, tRing - static_cast<double>(pc.t[k]), (hitPos.Phi()-anodeIntersection.Phi()) * 180. / TMath::Pi(), "hTiming");
}
}
}
}
for (int i = 0; i < sx3.multi; i++)
{
// plotting sx3 strip hits vs anode phi
if (sx3.ch[i] < 8)
plotter->Fill2D("AnodePhi_vs_SX3Strip", 100, -200, 200, 8 * 24, 0, 8 * 24, anodeIntersection.Phi() * 180. / TMath::Pi(), sx3.id[i] * 8 + sx3.ch[i]);
}
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 3)
{
plotter->Fill1D("PC_Z_proj_3C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
}
plotter->Fill2D("AnodeMaxE_Vs_Cathode_Sum_Energy", 2000, 0, 30000, 2000, 0, 30000, aEMax, cESum, "hGMPC");
plotter->Fill1D("Correlated_Cathode_MaxAnode", 6, 0, 5, corrcatMax.size(), "hGMPC");
plotter->Fill2D("Correlated_Cathode_VS_MaxAnodeEnergy", 6, 0, 5, 2000, 0, 30000, corrcatMax.size(), aEMax, "hGMPC");
plotter->Fill1D("AnodeHits", 12, 0, 11, anodeHits.size(), "hGMPC");
plotter->Fill2D("AnodeMaxE_vs_AnodeHits", 12, 0, 11, 2000, 0, 30000, anodeHits.size(), aEMax, "hGMPC");
if (anodeHits.size() < 1)
{
plotter->Fill1D("NoAnodeHits_CathodeHits", 6, 0, 5, cathodeHits.size(), "hGMPC");
}
return kTRUE;
}
void TrackRecon::Terminate()
{
plotter->FlushToDisk();
}

24
feeder.py Normal file
View File

@ -0,0 +1,24 @@
import keyboard
import time
fi = open("/tmp/eventlist.dat","r")
lines = fi.readlines()[3:]
fo = open("/tmp/coords","w")
for line in lines:
if("--" not in line):
print(line,end='')
fo.write(line)
fo.flush()
else:
if("end" in line):
while True:
time.sleep(2)
key = "n"
#key = input()
if 'n' in key:
break
fo.seek(0)
fo.truncate(0)
print(line,end='')
fo.close()
fi.close()

3
rootlogon.C Normal file
View File

@ -0,0 +1,3 @@
{
gSystem->SetBuildDir("./obj/",1);
}

33
run_19_22.sh Normal file
View File

@ -0,0 +1,33 @@
#Alpha runs at different spacer positions
#root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_009_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_QQQ.root Analyzer_QQQ_09.root;
###root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_010_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_QQQ.root Analyzer_QQQ_10.root;
#root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_012_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_QQQ.root Analyzer_QQQ_12.root;
#root -l -x Analyzer_QQQ_09.root Analyzer_QQQ_12.root -e "new TBrowser"
#root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_009_mapped.root -e 'tree->Process("DataDump.C+")'; mv Analyzer_QQQ.root Analyzer_QQQ_09.root;
#26Al runs
#root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_027_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_QQQ.root Analyzer_QQQ_27.root;
#root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_028_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_QQQ.root Analyzer_QQQ_28.root;
# root -l -x Analyzer_QQQ_27.root Analyzer_QQQ_28.root -e "new TBrowser"
#Proton runs at different spacer positions 26Al
root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_019_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run19.root;
#root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_021_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_QQQ.root Analyzer_QQQ_21.root;
root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_022_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run22.root;
root -l -x results_run19.root results_run22.root
#root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_084_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run84.root;
#root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_078_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run78.root;
#root -l -x results_run84.root -e "new TBrowser"
#root -q -l -x ../ANASEN_analysis/data/27Al_Data/ProtonRun_032_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run32.root;
#root -l -x results_run32.root -e "new TBrowser"
#root -q -l -x ../ANASEN_analysis/data/17F_Data/ProtonRun_022_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_QQQ.root Analyzer_QQQ_22.root;
#root -q -l -x ../ANASEN_analysis/data/17F_Data/ProtonRun_025_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_QQQ.root Analyzer_QQQ_25.root;
#root -q -l -x ../ANASEN_analysis/data/17F_Data/ProtonRun_026_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_QQQ.root Analyzer_QQQ_26.root;
#root -q -l -x ../ANASEN_analysis/data/17F_Data/ProtonRun_027_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_QQQ.root Analyzer_QQQ_27.root;
#root -q -l -x ../ANASEN_analysis/data/17F_Data/ProtonRun_028_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_QQQ.root Analyzer_QQQ_28.root;

30
run_sx3.sh Normal file
View File

@ -0,0 +1,30 @@
#Alpha runs at different spacer positions
#rm results_run*.root
#root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_009_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run09.root;
#root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_001_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run01.root;
#root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_002_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run02.root;
#root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_003_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run03.root;
#root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_004_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run04.root;
#alpha+gas 27Al
#root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_012_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run12.root;
#protons+gas, 27Al
#root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_022_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run22.root;
#root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_021_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run21.root;
#27Al reaction data
#root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_051_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run51.root;
#root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_078_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run78.root;
#root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_081_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run81.root;
#root -l -x results_run19.root results_run12.root -e "new TBrowser"
#17F alpha run with gas
#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_018_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run18.root;
root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_019_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run19.root;
#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_020_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run20.root;
#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_021_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run21.root;
#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Run_104_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run104.root;
#mv Analyzer_SX3.root results_run19.root;

23
sx3cal/17F/backgains.dat Normal file
View File

@ -0,0 +1,23 @@
1 front 1 back 2 4.4094
1 front 1 back 1 4.4094
1 front 2 back 2 4.4832
1 front 3 back 2 4.52103
1 front 3 back 1 4.5210
3 front 1 back 2 3.63215
3 front 1 back 1 3.70756
3 front 2 back 2 3.68208
3 front 2 back 1 3.86817
3 front 3 back 2 3.7334
3 front 3 back 1 3.70756
7 front 1 back 2 3.60769
7 front 1 back 1 3.46759
7 front 2 back 2 3.58356
7 front 2 back 1 3.49018
7 front 3 back 2 3.60769
7 front 3 back 1 3.46759
9 front 1 back 2 3.58356
9 front 1 back 1 3.44529
9 front 2 back 2 3.58356
9 front 2 back 1 3.46759
9 front 3 back 2 3.63215
9 front 3 back 1 3.46759

96
sx3cal/17F/backgains.dat.unity Normal file
View File

@ -0,0 +1,96 @@
1 0 0 1.
1 1 0 1.
1 2 0 1.
1 3 0 1.
1 0 1 1.
1 1 1 1.
1 2 1 1.
1 3 1 1.
1 0 2 1.
1 1 2 1.
1 2 2 1.
1 3 2 1.
1 0 3 1.
1 1 3 1.
1 2 3 1.
1 3 3 1.
7 0 0 1.
7 1 0 1.
7 2 0 1.
7 3 0 1.
7 0 1 1.
7 1 1 1.
7 2 1 1.
7 3 1 1.
7 0 2 1.
7 1 2 1.
7 2 2 1.
7 3 2 1.
7 0 3 1.
7 1 3 1.
7 2 3 1.
7 3 3 1.
0 0 0 1.
0 1 0 1.
0 2 0 1.
0 3 0 1.
0 0 1 1.
0 1 1 1.
0 2 1 1.
0 3 1 1.
0 0 2 1.
0 1 2 1.
0 2 2 1.
0 3 2 1.
0 0 3 1.
0 1 3 1.
0 2 3 1.
0 3 3 1.
2 0 0 1.
2 1 0 1.
2 2 0 1.
2 3 0 1.
2 0 1 1.
2 1 1 1.
2 2 1 1.
2 3 1 1.
2 0 2 1.
2 1 2 1.
2 2 2 1.
2 3 2 1.
2 0 3 1.
2 1 3 1.
2 2 3 1.
2 3 3 1.
9 0 0 1.
9 1 0 1.
9 2 0 1.
9 3 0 1.
9 0 1 1.
9 1 1 1.
9 2 1 1.
9 3 1 1.
9 0 2 1.
9 1 2 1.
9 2 2 1.
9 3 2 1.
9 0 3 1.
9 1 3 1.
9 2 3 1.
9 3 3 1.
3 0 0 1.
3 1 0 1.
3 2 0 1.
3 3 0 1.
3 0 1 1.
3 1 1 1.
3 2 1 1.
3 3 1 1.
3 0 2 1.
3 1 2 1.
3 2 2 1.
3 3 2 1.
3 0 3 1.
3 1 3 1.
3 2 3 1.
3 3 3 1.

12
sx3cal/17F/frontgains.dat Normal file
View File

@ -0,0 +1,12 @@
1 lengthcal front 1 1.5121 60.4839
1 lengthcal front 2 -1.5625 62.5
1 lengthcal front 3 2.72177 60.4839
3 lengthcal front 1 -0.595088 59.5088
3 lengthcal front 2 -4.53935 58.5723
3 lengthcal front 3 4.08107 60.4603
7 lengthcal front 1 1.14329 45.7317
7 lengthcal front 2 0.115661 46.2646
7 lengthcal front 3 2.90179 44.6429
9 lengthcal front 1 0.115732 46.2928
9 lengthcal front 2 0.799176 45.6672
9 lengthcal front 3 1.68159 48.0453

24
sx3cal/17F/frontgains.dat.unity Normal file
View File

@ -0,0 +1,24 @@
9 temp temp 0 0. 1.
9 temp temp 1 0. 1.
9 temp temp 2 0. 1.
9 temp temp 3 0. 1.
7 temp temp 0 0. 1.
7 temp temp 1 0. 1.
7 temp temp 2 0. 1.
7 temp temp 3 0. 1.
1 temp temp 0 0. 1.
1 temp temp 1 0. 1.
1 temp temp 2 0. 1.
1 temp temp 3 0. 1.
2 temp temp 0 0. 1.
2 temp temp 1 0. 1.
2 temp temp 2 0. 1.
2 temp temp 3 0. 1.
0 temp temp 0 0. 1.
0 temp temp 1 0. 1.
0 temp temp 2 0. 1.
0 temp temp 3 0. 1.
3 temp temp 0 0. 1.
3 temp temp 1 0. 1.
3 temp temp 2 0. 1.
3 temp temp 3 0. 1.

16
sx3cal/17F/rightgains.dat Normal file
View File

@ -0,0 +1,16 @@
1 0 1678.38 1.0
1 1 1678.38 1.07163
1 2 1693.99 1.1035
1 3 1667.9 0.975015
3 0 1597.96 1.0
3 1 1597.96 1.02536
3 2 1821.48 1.29182
3 3 1607.52 0.928543
7 0 1773.34 1.
7 1 1773.34 1.14263
7 2 1573.79 1.06715
7 3 1542.41 0.956475
9 0 1555.08 1.
9 1 1555.08 0.998252
9 2 1559.36 1.00299
9 3 1585.79 1.01582

24
sx3cal/17F/rightgains.dat.unity Normal file
View File

@ -0,0 +1,24 @@
9 0 1000 1.0
9 1 1000 1.0
9 2 1000 1.0
9 3 1000 1.0
7 0 1000 1.0
7 1 1000 1.0
7 2 1000 1.0
7 3 1000 1.0
2 0 1000 1.0
2 1 1000 1.0
2 2 1000 1.0
2 3 1000 1.0
0 0 1000 1.0
0 1 1000 1.0
0 2 1000 1.0
0 3 1000 1.0
1 0 1000 1.0
1 1 1000 1.0
1 2 1000 1.0
1 3 1000 1.0
3 0 1000 1.0
3 1 1000 1.0
3 2 1000 1.0
3 3 1000 1.0

28
sx3cal/26Al/backgains.dat Normal file
View File

@ -0,0 +1,28 @@
1 front 0 back 2 4.03168
1 front 1 back 2 4.03168
1 front 2 back 2 4.11533
1 front 3 back 2 4.17315
7 front 0 back 2 4.26886
7 front 0 back 1 3.44529
7 front 1 back 2 4.26886
7 front 1 back 1 3.44529
7 front 2 back 2 4.26886
7 front 2 back 1 3.46759
7 front 3 back 2 4.26886
7 front 3 back 1 3.44529
9 front 0 back 2 3.63215
9 front 0 back 1 3.42327
9 front 1 back 2 3.63215
9 front 1 back 1 3.42327
9 front 2 back 2 3.65694
9 front 2 back 1 3.46759
9 front 3 back 2 3.68208
9 front 3 back 1 3.42327
3 front 0 back 2 3.
3 front 0 back 1 3.
3 front 1 back 2 3.65694
3 front 1 back 1 3.68208
3 front 2 back 2 3.70756
3 front 2 back 1 3.78616
3 front 3 back 2 3.7334
3 front 3 back 1 3.68208

96
sx3cal/26Al/backgains.dat.unity Normal file
View File

@ -0,0 +1,96 @@
1 0 0 1.
1 1 0 1.
1 2 0 1.
1 3 0 1.
1 0 1 1.
1 1 1 1.
1 2 1 1.
1 3 1 1.
1 0 2 1.
1 1 2 1.
1 2 2 1.
1 3 2 1.
1 0 3 1.
1 1 3 1.
1 2 3 1.
1 3 3 1.
7 0 0 1.
7 1 0 1.
7 2 0 1.
7 3 0 1.
7 0 1 1.
7 1 1 1.
7 2 1 1.
7 3 1 1.
7 0 2 1.
7 1 2 1.
7 2 2 1.
7 3 2 1.
7 0 3 1.
7 1 3 1.
7 2 3 1.
7 3 3 1.
0 0 0 1.
0 1 0 1.
0 2 0 1.
0 3 0 1.
0 0 1 1.
0 1 1 1.
0 2 1 1.
0 3 1 1.
0 0 2 1.
0 1 2 1.
0 2 2 1.
0 3 2 1.
0 0 3 1.
0 1 3 1.
0 2 3 1.
0 3 3 1.
2 0 0 1.
2 1 0 1.
2 2 0 1.
2 3 0 1.
2 0 1 1.
2 1 1 1.
2 2 1 1.
2 3 1 1.
2 0 2 1.
2 1 2 1.
2 2 2 1.
2 3 2 1.
2 0 3 1.
2 1 3 1.
2 2 3 1.
2 3 3 1.
9 0 0 1.
9 1 0 1.
9 2 0 1.
9 3 0 1.
9 0 1 1.
9 1 1 1.
9 2 1 1.
9 3 1 1.
9 0 2 1.
9 1 2 1.
9 2 2 1.
9 3 2 1.
9 0 3 1.
9 1 3 1.
9 2 3 1.
9 3 3 1.
3 0 0 1.
3 1 0 1.
3 2 0 1.
3 3 0 1.
3 0 1 1.
3 1 1 1.
3 2 1 1.
3 3 1 1.
3 0 2 1.
3 1 2 1.
3 2 2 1.
3 3 2 1.
3 0 3 1.
3 1 3 1.
3 2 3 1.
3 3 3 1.

16
sx3cal/26Al/frontgains.dat Normal file
View File

@ -0,0 +1,16 @@
1 lengthcal front 0 0.878906 58.5938
1 lengthcal front 1 1.42045 56.8182
1 lengthcal front 2 -2.55682 56.8182
1 lengthcal front 3 2.55682 56.8182
7 lengthcal front 0 0.425806 42.5806
7 lengthcal front 1 1.92004 45.1774
7 lengthcal front 2 1.11607 44.6429
7 lengthcal front 3 3.45909 44.6334
9 lengthcal front 0 1.82872 45.7181
9 lengthcal front 1 1.01649 45.1774
9 lengthcal front 2 1.46827 45.1774
9 lengthcal front 3 2.54513 46.2751
3 lengthcal front 0 0. 50.
3 lengthcal front 1 1.1713 58.5652
3 lengthcal front 2 -3.07505 58.5723
3 lengthcal front 3 4.0726 60.3348

20
sx3cal/26Al/frontgains.dat.unity Normal file
View File

@ -0,0 +1,20 @@
9 temp temp 0 0. 1.
9 temp temp 1 0. 1.
9 temp temp 2 0. 1.
9 temp temp 3 0. 1.
7 temp temp 0 0. 1.
7 temp temp 1 0. 1.
7 temp temp 2 0. 1.
7 temp temp 3 0. 1.
1 temp temp 0 0. 1.
1 temp temp 1 0. 1.
1 temp temp 2 0. 1.
1 temp temp 3 0. 1.
2 temp temp 0 0. 1.
2 temp temp 1 0. 1.
2 temp temp 2 0. 1.
2 temp temp 3 0. 1.
0 temp temp 0 0. 1.
0 temp temp 1 0. 1.
0 temp temp 2 0. 1.
0 temp temp 3 0. 1.

16
sx3cal/26Al/rightgains.dat Normal file
View File

@ -0,0 +1,16 @@
1 0 1221.23 0.648782
1 1 1819.66 1.06196
1 2 1860.02 1.11979
1 3 1825.44 0.964989
7 0 1609.63 1.04668
7 1 1734.45 1.12285
7 2 1538.97 1.0486
7 3 1524.57 0.951587
9 0 1672.38 1.11321
9 1 1542.13 1.01442
9 2 1540.38 0.967847
9 3 1560.42 0.969022
3 0 1000.0 1.
3 1 1539.42 1.0422
3 2 1720.12 1.31534
3 3 1562.16 1.00415

24
sx3cal/26Al/rightgains.dat.unity Normal file
View File

@ -0,0 +1,24 @@
9 0 1000 1.0
9 1 1000 1.0
9 2 1000 1.0
9 3 1000 1.0
7 0 1000 1.0
7 1 1000 1.0
7 2 1000 1.0
7 3 1000 1.0
2 0 1000 1.0
2 1 1000 1.0
2 2 1000 1.0
2 3 1000 1.0
0 0 1000 1.0
0 1 1000 1.0
0 2 1000 1.0
0 3 1000 1.0
1 0 1000 1.0
1 1 1000 1.0
1 2 1000 1.0
1 3 1000 1.0
3 0 1000 1.0
3 1 1000 1.0
3 2 1000 1.0
3 3 1000 1.0

8
sx3cal/EXFit.C
View File

@ -1,6 +1,6 @@
{ {
int index = 3; int index = 1;
TFile *f = new TFile("../results_SX3_run12.root"); TFile *f = new TFile("../results_run19.root");
TH2F *h2=NULL; TH2F *h2=NULL;
TH1F *h1x=NULL, *h1y=NULL; TH1F *h1x=NULL, *h1y=NULL;
//f->cd("evsx"); //f->cd("evsx");
@ -17,8 +17,8 @@
h2 = (TH2F*)(f->Get(Form("evsx/be_vs_x_sx3_id_%d_f%d_b%d",index,i,backnum))); h2 = (TH2F*)(f->Get(Form("evsx/be_vs_x_sx3_id_%d_f%d_b%d",index,i,backnum)));
auto macro = [&]() { auto macro = [&]() {
h1x = (TH1F*)(h2->ProjectionX("_px")); h1x = (TH1F*)(h2->ProjectionX("_px"));
double xleft = h1x->GetBinCenter(h1x->FindFirstBinAbove(h1x->GetMaximum()*0.25)); double xleft = h1x->GetBinCenter(h1x->FindFirstBinAbove(h1x->GetMaximum()*0.4));
double xright = h1x->GetBinCenter(h1x->FindLastBinAbove(h1x->GetMaximum()*0.25)); double xright = h1x->GetBinCenter(h1x->FindLastBinAbove(h1x->GetMaximum()*0.4));
//h1x->GetXaxis()->SetRangeUser(4*xleft, xright*4); //h1x->GetXaxis()->SetRangeUser(4*xleft, xright*4);
h1x->Draw(); h1x->Draw();
TLine L1(xleft,0,xleft,h1x->GetMaximum()); L1.SetLineColor(kRed); L1.Draw("SAME"); TLine L1(xleft,0,xleft,h1x->GetMaximum()); L1.SetLineColor(kRed); L1.Draw("SAME");

8
sx3cal/LRFit.C
View File

@ -1,8 +1,8 @@
{ {
TFile *f = new TFile("../results_SX3_run12.root"); TFile *f = new TFile("../results_run19.root");
f->cd("l_vs_r"); f->cd("l_vs_r");
f->ls(); gDirectory->ls();
int clkpos = 3; int clkpos = 13;
std::ofstream ofile(Form("rightgains%d.dat",clkpos)); std::ofstream ofile(Form("rightgains%d.dat",clkpos));
for(int i=1; i<4; i++) { for(int i=1; i<4; i++) {
TH2F h2(*(TH2F*)(f->Get(Form("l_vs_r/l_vs_r_sx3_id_%d_f%d",clkpos,i)))); TH2F h2(*(TH2F*)(f->Get(Form("l_vs_r/l_vs_r_sx3_id_%d_f%d",clkpos,i))));
@ -14,7 +14,7 @@
gPad->Update(); gPad->Update();
while(gPad->WaitPrimitive());*/ while(gPad->WaitPrimitive());*/
int leftbin = hproj.FindFirstBinAbove(hproj.GetMaximum()*0.1); int leftbin = hproj.FindFirstBinAbove(hproj.GetMaximum()*0.4);
int rightbin = hproj.FindLastBinAbove(hproj.GetMaximum()*0.1); int rightbin = hproj.FindLastBinAbove(hproj.GetMaximum()*0.1);
TH1F h1(*(TH1F*)(h2.ProfileX("_pfx",leftbin,rightbin))); TH1F h1(*(TH1F*)(h2.ProfileX("_pfx",leftbin,rightbin)));

3
sx3cal/backgains1.dat Normal file
View File

@ -0,0 +1,3 @@
1 front 1 back 2 4.4094
1 front 2 back 2 4.4832
1 front 3 back 2 4.52103

6
sx3cal/backgains9.dat Normal file
View File

@ -0,0 +1,6 @@
9 front 1 back 2 3.58356
9 front 1 back 1 3.44529
9 front 2 back 2 3.58356
9 front 2 back 1 3.46759
9 front 3 back 2 3.63215
9 front 3 back 1 3.46759

3
sx3cal/frontgains1.dat Normal file
View File

@ -0,0 +1,3 @@
1 lengthcal front 1 1.5121 60.4839
1 lengthcal front 2 -1.5625 62.5
1 lengthcal front 3 2.72177 60.4839

3
sx3cal/frontgains9.dat Normal file
View File

@ -0,0 +1,3 @@
9 lengthcal front 1 0.115732 46.2928
9 lengthcal front 2 0.799176 45.6672
9 lengthcal front 3 1.68159 48.0453

3
sx3cal/rightgains1.dat Normal file
View File

@ -0,0 +1,3 @@
1 1 1678.38 1.07163
1 2 1693.99 1.1035
1 3 1667.9 0.975015

3
sx3cal/rightgains13.dat Normal file
View File

@ -0,0 +1,3 @@
13 1 1 1
13 2 1 1
13 3 1 1

3
sx3cal/rightgains9.dat Normal file
View File

@ -0,0 +1,3 @@
9 1 1555.08 0.998252
9 2 1559.36 1.00299
9 3 1585.79 1.01582