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

modified: TrackRecon.C inlcuded timing plots for the PC

Browse Source
This commit is contained in:
Vignesh Sitaraman 2026-01-13 12:04:14 -05:00
parent 82c2127b4d
commit 19286055ea
1 changed files with 178 additions and 58 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

150
TrackRecon.C
View File

@ -30,6 +30,7 @@ double qqqGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
bool qqqGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}}; bool qqqGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
double qqqCalib[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}}; double qqqCalib[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
bool qqqCalibValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}}; bool qqqCalibValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
// TCutg *cutQQQ;
// PC Arrays // PC Arrays
double pcSlope[48]; double pcSlope[48];
@ -62,11 +63,11 @@ void TrackRecon::Begin(TTree * /*tree*/)
TVector3 a, c, diff; TVector3 a, c, diff;
double a2, ac, c2, adiff, cdiff, denom, alpha; double a2, ac, c2, adiff, cdiff, denom, alpha;
for (int i = 0; i < pwinstance.An.size(); i++) for (size_t i = 0; i < pwinstance.An.size(); i++)
{ {
a = pwinstance.An[i].first - pwinstance.An[i].second; a = pwinstance.An[i].first - pwinstance.An[i].second;
for (int j = 0; j < pwinstance.Ca.size(); j++) for (size_t j = 0; j < pwinstance.Ca.size(); j++)
{ {
c = pwinstance.Ca[j].first - pwinstance.Ca[j].second; c = pwinstance.Ca[j].first - pwinstance.Ca[j].second;
diff = pwinstance.An[i].first - pwinstance.Ca[j].first; diff = pwinstance.An[i].first - pwinstance.Ca[j].first;
@ -116,6 +117,25 @@ void TrackRecon::Begin(TTree * /*tree*/)
std::cerr << "Error opening slope_intercept.txt" << std::endl; 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) ... // ... (Load QQQ Gains and Calibs - same as before) ...
{ {
std::string filename = "qqq_GainMatch.txt"; std::string filename = "qqq_GainMatch.txt";
@ -202,7 +222,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
if (pc.index[k] < 24 && pc.e[k] > 50) if (pc.index[k] < 24 && pc.e[k] > 50)
{ {
plotter->Fill2D("QQQ_Vs_PC_Energy", 400, 0, 4000, 1000, 0, 16000, qqq.e[i], pc.e[k]); plotter->Fill2D("QQQ_Vs_PC_Energy", 400, 0, 4000, 1000, 0, 16000, qqq.e[i], pc.e[k]);
plotter->Fill2D("QQQ_Index_Vs_PC_Index", 16*4, 0, 16*4, 24, 0, 24, qqq.index[i], pc.index[k]); plotter->Fill2D("QQQ_Index_Vs_PC_Index", 16 * 8, 0, 16 * 8, 24, 0, 24, qqq.index[i], pc.index[k]);
} }
} }
@ -212,14 +232,14 @@ Bool_t TrackRecon::Process(Long64_t entry)
{ {
qqqCount++; qqqCount++;
if (qqq.id[i] == qqq.id[j])
{
int chWedge = -1; int chWedge = -1;
int chRing = -1; int chRing = -1;
float eWedge = 0.0; float eWedge = 0.0;
float eWedgeMeV = 0.0; float eWedgeMeV = 0.0;
float eRing = 0.0; float eRing = 0.0;
float eRingMeV = 0.0; float 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]) if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && qqqGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16])
{ {
@ -227,6 +247,8 @@ Bool_t TrackRecon::Process(Long64_t entry)
eWedge = qqq.e[i] * qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]; eWedge = qqq.e[i] * qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
chRing = qqq.ch[j] - 16; chRing = qqq.ch[j] - 16;
eRing = qqq.e[j]; 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]) else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && qqqGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16])
{ {
@ -234,10 +256,14 @@ Bool_t TrackRecon::Process(Long64_t entry)
eWedge = qqq.e[j] * qqqGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16]; eWedge = qqq.e[j] * qqqGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
chRing = qqq.ch[i] - 16; chRing = qqq.ch[i] - 16;
eRing = qqq.e[i]; eRing = qqq.e[i];
tRing = static_cast<double>(qqq.t[i]);
tWedge = static_cast<double>(qqq.t[j]);
} }
else else
continue; continue;
plotter->Fill1D("Wedgetime_Vs_Ringtime", 2000, -1000, 1000, tWedge - tRing, "hCalQQQ");
if (qqqCalibValid[qqq.id[i]][chRing][chWedge]) if (qqqCalibValid[qqq.id[i]][chRing][chWedge])
{ {
eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000; eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000;
@ -252,8 +278,19 @@ Bool_t TrackRecon::Process(Long64_t entry)
{ {
if (pc.index[k] < 24 && pc.e[k] > 50) if (pc.index[k] < 24 && pc.e[k] > 50)
{ {
plotter->Fill2D("QQQ_Calib_Vs_PC_Energy", 1000, 0, 16, 2000, 0, 30000, eWedgeMeV, pc.e[k], "hCalQQQ"); plotter->Fill2D("QQQ_CalibW_Vs_PC_Energy", 1000, 0, 16, 2000, 0, 30000, eWedgeMeV, pc.e[k], "hCalQQQ");
plotter->Fill2D("QQQ_Calib_Vs_PC_Energy", 1000, 0, 16, 2000, 0, 30000, eRingMeV, 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", 20000, -1000, 1000, 16, 0, 16, tRing - static_cast<double>(pc.t[k]), chRing, "hCalQQQ");
} }
} }
@ -272,11 +309,14 @@ Bool_t TrackRecon::Process(Long64_t entry)
} }
} }
} }
}
plotter->Fill1D("QQQ_Multiplicity", 10, 0, 10, qqqCount, "hCalQQQ"); plotter->Fill1D("QQQ_Multiplicity", 10, 0, 10, qqqCount, "hCalQQQ");
// PC Gain Matching and Filling // 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++) for (int i = 0; i < pc.multi; i++)
{ {
if (pc.e[i] > 10) if (pc.e[i] > 10)
@ -284,15 +324,44 @@ Bool_t TrackRecon::Process(Long64_t entry)
plotter->Fill2D("PC_Index_Vs_Energy", 48, 0, 48, 2000, 0, 30000, pc.index[i], pc.e[i], "hRawPC"); plotter->Fill2D("PC_Index_Vs_Energy", 48, 0, 48, 2000, 0, 30000, pc.index[i], pc.e[i], "hRawPC");
} }
if (pc.index[i] >= 0 && pc.index[i] < 48) if (pc.index[i] < 48)
{ {
// FIX: pcSlope defaults to 1.0 now, so this won't zero out data if file entry is missing // FIX: pcSlope defaults to 1.0 now, so this won't zero out data if file entry is missing
pc.e[i] = pcSlope[pc.index[i]] * pc.e[i] + pcIntercept[pc.index[i]]; pc.e[i] = pcSlope[pc.index[i]] * pc.e[i] + pcIntercept[pc.index[i]];
plotter->Fill2D("PC_Index_VS_GainMatched_Energy", 24, 0, 24, 2000, 0, 30000, pc.index[i], pc.e[i], "hGMPC"); plotter->Fill2D("PC_Index_VS_GainMatched_Energy", 24, 0, 24, 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", 400, -1000, 1000, anodeT - cathodeT, "hGMPC");
plotter->Fill2D("PC_Time_Vs_QQQ_ch", 400, -1000, 1000, 16*8, 0, 16*8, anodeT - cathodeT, qqq.ch[j], "hGMPC");
plotter->Fill2D("PC_Time_vs_AIndex", 400, -1000, 1000, 24, 0, 24, anodeT - cathodeT, anodeIndex, "hGMPC");
plotter->Fill2D("PC_Time_vs_CIndex", 400, -1000, 1000, 24, 0, 24, anodeT - cathodeT, cathodeIndex, "hGMPC");
plotter->Fill1D("PC_Time_A"+std::to_string(anodeIndex)+"_C"+std::to_string(cathodeIndex), 400, -1000, 1000, anodeT - cathodeT, "TimingPC");
}
for (int j = 0; j < sx3.multi; j++)
{
plotter->Fill1D("PC_Time_sx3", 400, -1000, 1000, anodeT - cathodeT, "hGMPC");
}
plotter->Fill1D("PC_Time", 400, -1000, 1000, anodeT - cathodeT, "hGMPC");
}
for (int j = i + 1; j < pc.multi; j++) for (int j = i + 1; j < pc.multi; j++)
{ {
plotter->Fill2D("PC_Coincidence_Matrix", 24, 0, 24, 24, 24, 48, pc.index[i], pc.index[j], "hRawPC"); 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");
} }
} }
@ -354,7 +423,6 @@ Bool_t TrackRecon::Process(Long64_t entry)
plotter->Fill2D("AnodeMax_Vs_Cathode_Coincidence_Matrix", 24, 0, 24, 24, 0, 24, aIDMax, cID, "hGMPC"); plotter->Fill2D("AnodeMax_Vs_Cathode_Coincidence_Matrix", 24, 0, 24, 24, 0, 24, aIDMax, cID, "hGMPC");
plotter->Fill2D("Anode_Vs_Cathode_Coincidence_Matrix", 24, 0, 24, 24, 0, 24, aID, cID, "hGMPC"); plotter->Fill2D("Anode_Vs_Cathode_Coincidence_Matrix", 24, 0, 24, 24, 0, 24, aID, cID, "hGMPC");
plotter->Fill2D("Anode_vs_CathodeE", 2000, 0, 30000, 2000, 0, 30000, aE, cE, "hGMPC"); plotter->Fill2D("Anode_vs_CathodeE", 2000, 0, 30000, 2000, 0, 30000, aE, cE, "hGMPC");
for (int j = -4; j < 3; j++) for (int j = -4; j < 3; j++)
{ {
if ((aIDMax + 24 + j) % 24 == 23 - cID) if ((aIDMax + 24 + j) % 24 == 23 - cID)
@ -389,26 +457,78 @@ Bool_t TrackRecon::Process(Long64_t entry)
if (anodeIntersection.Z() != 0) if (anodeIntersection.Z() != 0)
{ {
plotter->Fill1D("PC_Z_Projection", 600, -300, 300, anodeIntersection.Z(), "hGMPC"); plotter->Fill1D("PC_Z_Projection", 600, -300, 300, anodeIntersection.Z(), "hGMPC");
plotter->Fill2D("Z_Proj_VsDelTime", 600, -300, 300, 400, -1000, 1000, anodeIntersection.Z(), anodeT - cathodeT, "hGMPC");
} }
if(anodeIntersection.Z() != 0 && cathodeHits.size()==1) if (anodeIntersection.Z() != 0 && cathodeHits.size() == 1)
{ {
plotter->Fill1D("PC_Z_proj_1C", 600, -300, 300, anodeIntersection.Z(), "hGMPC"); plotter->Fill1D("PC_Z_proj_1C", 600, -300, 300, anodeIntersection.Z(), "hGMPC");
} }
if(anodeIntersection.Z() != 0 && cathodeHits.size()==2) if (anodeIntersection.Z() != 0 && cathodeHits.size() == 2)
{ {
plotter->Fill1D("PC_Z_proj_2C", 600, -300, 300, anodeIntersection.Z(), "hGMPC"); plotter->Fill1D("PC_Z_proj_2C", 600, -300, 300, anodeIntersection.Z(), "hGMPC");
} }
if(anodeIntersection.Z() != 0 && cathodeHits.size()==3)
for (int i = 0; i < qqq.multi; i++)
{
for (int j = i + 1; j < qqq.multi; j++)
{
if (qqq.id[i] == qqq.id[j])
{
int chWedge = -1;
int chRing = -1;
float eWedge = 0.0;
float eWedgeMeV = 0.0;
float eRing = 0.0;
float eRingMeV = 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];
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;
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, "hGMPC");
}
plotter->Fill2D("PC_Z_vs_QQQRing_Det" + std::to_string(qqqID), 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hGMPC");
}
}
}
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 3)
{ {
plotter->Fill1D("PC_Z_proj_3C", 600, -300, 300, anodeIntersection.Z(), "hGMPC"); plotter->Fill1D("PC_Z_proj_3C", 600, -300, 300, anodeIntersection.Z(), "hGMPC");
} }
plotter->Fill2D("AnodeMaxE_Vs_Cathode_Sum_Energy", 2000, 0, 30000, 2000, 0, 30000, aEMax, cESum, "hGMPC"); 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->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->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->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) if (anodeHits.size() < 1)
{ {