From 19286055eaf592969d0d41f99a3b3c6e7eff1a63 Mon Sep 17 00:00:00 2001 From: vsitaraman Date: Tue, 13 Jan 2026 12:04:14 -0500 Subject: [PATCH] modified: TrackRecon.C inlcuded timing plots for the PC --- TrackRecon.C | 236 ++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 178 insertions(+), 58 deletions(-) diff --git a/TrackRecon.C b/TrackRecon.C index 5cc32e6..983ac91 100644 --- a/TrackRecon.C +++ b/TrackRecon.C @@ -30,6 +30,7 @@ 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]; @@ -62,11 +63,11 @@ void TrackRecon::Begin(TTree * /*tree*/) TVector3 a, c, diff; 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; - 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; 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; } + // 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.txt"; @@ -202,7 +222,7 @@ Bool_t TrackRecon::Process(Long64_t entry) 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_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,63 +232,79 @@ Bool_t TrackRecon::Process(Long64_t entry) { qqqCount++; - 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; + 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]) { - int chWedge = -1; - int chRing = -1; - float eWedge = 0.0; - float eWedgeMeV = 0.0; - float eRing = 0.0; - float eRingMeV = 0.0; + 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(qqq.t[j]); + tWedge = static_cast(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(qqq.t[i]); + tWedge = static_cast(qqq.t[j]); + } + else + continue; - 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]; - } - 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]; - } - else - continue; + plotter->Fill1D("Wedgetime_Vs_Ringtime", 2000, -1000, 1000, tWedge - tRing, "hCalQQQ"); - 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 (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; - plotter->Fill2D("WedgeE_Vs_RingECal", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ"); + plotter->Fill2D("WedgeE_Vs_RingECal", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ"); - for (int k = 0; k < pc.multi; k++) + for (int k = 0; k < pc.multi; k++) + { + if (pc.index[k] < 24 && pc.e[k] > 50) { - 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(pc.t[k]) < 0 && tRing - static_cast(pc.t[k]) > -600) { - plotter->Fill2D("QQQ_Calib_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_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(pc.t[k]), chRing, "hCalQQQ"); } + } - 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); + 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); - plotter->Fill2D("QQQPolarPlot", 16 * 4, -TMath::Pi(), TMath::Pi(), 32, 40, 100, theta, rho, "hCalQQQ"); + plotter->Fill2D("QQQPolarPlot", 16 * 4, -TMath::Pi(), TMath::Pi(), 32, 40, 100, theta, rho, "hCalQQQ"); - if (!HitNonZero) - { - double x = rho * TMath::Cos(theta); - double y = rho * TMath::Sin(theta); - hitPos.SetXYZ(x, y, 23 + 75 + 30); - HitNonZero = true; - } + if (!HitNonZero) + { + double x = rho * TMath::Cos(theta); + double y = rho * TMath::Sin(theta); + hitPos.SetXYZ(x, y, 23 + 75 + 30); + HitNonZero = true; } } } @@ -277,6 +313,10 @@ Bool_t TrackRecon::Process(Long64_t entry) plotter->Fill1D("QQQ_Multiplicity", 10, 0, 10, qqqCount, "hCalQQQ"); // 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) @@ -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"); } - 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 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"); } + + if (pc.index[i] < 24) + { + anodeT = static_cast(pc.t[i]); + anodeIndex = pc.index[i]; + } + else + { + cathodeT = static_cast(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++) { - 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("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"); - for (int j = -4; j < 3; j++) { if ((aIDMax + 24 + j) % 24 == 23 - cID) @@ -366,7 +434,7 @@ Bool_t TrackRecon::Process(Long64_t entry) } } } - + TVector3 anodeIntersection; anodeIntersection.Clear(); if (qqq1000cut) @@ -375,7 +443,7 @@ Bool_t TrackRecon::Process(Long64_t entry) for (const auto &corr : corrcatMax) { if (Crossover[aIDMax][corr.first][0].z > 9000000) - continue; + continue; if (cESum > 0) { x += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].x; @@ -385,30 +453,82 @@ Bool_t TrackRecon::Process(Long64_t entry) } anodeIntersection = TVector3(x, y, z); } - + if (anodeIntersection.Z() != 0) { 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"); } - 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"); } - 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->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) {