From f7248f6a4e3169c8f9ddf71e711c0e5de354d3a3 Mon Sep 17 00:00:00 2001 From: vsitaraman Date: Thu, 11 Jun 2026 16:51:41 -0400 Subject: [PATCH] modified: Armory/PC_StepLadder_Correction.h modified: TrackRecon.C trying to benchmeark A1Cn events against a1c2 and failing due to flaww in logic modified: run_tr.sh modified: scratch/plot_slope_scan.C --- Armory/PC_StepLadder_Correction.h | 20 +-- TrackRecon.C | 218 ++++++++++++++++++++++++++---- run_tr.sh | 13 +- scratch/plot_slope_scan.C | 22 +-- 4 files changed, 222 insertions(+), 51 deletions(-) diff --git a/Armory/PC_StepLadder_Correction.h b/Armory/PC_StepLadder_Correction.h index 79c3ea5..5fd75a7 100644 --- a/Armory/PC_StepLadder_Correction.h +++ b/Armory/PC_StepLadder_Correction.h @@ -46,16 +46,16 @@ double model_a1c1(double* x, double* p) { return result; } -double model_invert_a1c1(double *y, double *q) { - double result=y[0]; -/* double slope = 0.40; - double factor = 5.0; - if(TMath::Abs(y[0]) < 16.2/slope) result = y[0]/slope; - else if(TMath::Abs(y[0]) < 49.8/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor; - else if(TMath::Abs(y[0]) < 85.2/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor*2; - else result=y[0]/slope-TMath::Sign(1.0,y[0])*factor**/; - return result+40; -} +// double model_invert_a1c1(double *y, double *q) { +// double result=y[0]; +// /* double slope = 0.40; +// double factor = 5.0; +// if(TMath::Abs(y[0]) < 16.2/slope) result = y[0]/slope; +// else if(TMath::Abs(y[0]) < 49.8/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor; +// else if(TMath::Abs(y[0]) < 85.2/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor*2; +// else result=y[0]/slope-TMath::Sign(1.0,y[0])*factor**/; +// return result+40; +// } /*void testmodel() { diff --git a/TrackRecon.C b/TrackRecon.C index 913048c..66df707 100644 --- a/TrackRecon.C +++ b/TrackRecon.C @@ -44,6 +44,7 @@ bool doPCSX3ClusterAnalysis = true; bool doPCQQQClusterAnalysis = true; bool doOldAnalysis = false; bool do27AlapAnalysis = false; +bool BenchMark = true; double source_vertex = 53; // 53 const double qqq_z = 105.0; double z_entrance = -174.3 - 9.7 - 100.0; @@ -143,13 +144,14 @@ void protonAlphaHistograms(HistPlotter *plotter, std::vector QQQ_Events, void miscHistograms_oneWire(HistPlotter *plotter, std::vector QQQ_Events, std::vector>> aClusters); void protonMiscHistograms(HistPlotter *plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events); void protonMiscHistograms_sx3(HistPlotter *plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events); -void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events); -void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events); +void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events, + const std::vector>> &aClusters, const std::vector>> &cClusters); +void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events, + const std::vector>> &aClusters, const std::vector>> &cClusters); void TrackRecon::Begin(TTree * /*tree*/) { pcfix_func.SetNpx(100000); - ///// ---------Set Environment Variables--------- ///// TString option = GetOption(); if (option != "") @@ -598,8 +600,8 @@ Bool_t TrackRecon::Process(Long64_t entry) // if(eRingMeV<1.2 || eWedgeMeV<1.2) continue; // double theta = 2 * TMath::Pi() * (-qqq.id[i] * 16 + (15 - chWedge) + 0.5) / (16 * 4); - - double theta = (M_PI/180.)*(-90*qqq.id[i]+(87./16.)*((15-chWedge)+0.5)+3.0); + + double theta = (M_PI / 180.) * (-90 * qqq.id[i] + (87. / 16.) * ((15 - chWedge) + 0.5) + 3.0); 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 @@ -1024,11 +1026,11 @@ Bool_t TrackRecon::Process(Long64_t entry) if (doPCSX3ClusterAnalysis) { - PCSX3ClusterAnalysis(plotter, QQQ_Events, SX3_Events, PC_Events); + PCSX3ClusterAnalysis(plotter, QQQ_Events, SX3_Events, PC_Events, aClusters, cClusters); } if (doPCQQQClusterAnalysis) { - PCQQQClusterAnalysis(plotter, QQQ_Events, SX3_Events, PC_Events); + PCQQQClusterAnalysis(plotter, QQQ_Events, SX3_Events, PC_Events, aClusters, cClusters); } ///////////////////nA analysis using pseudo-wire (GetPseudoWire + getClosestWirePosAtWirePhi)/////////////////// @@ -1362,7 +1364,9 @@ void protonAlphaHistograms(HistPlotter *plotter, std::vector QQQ_Events, return; } -void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events){ +void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events, + const std::vector>> &aClusters, const std::vector>> &cClusters) +{ for (auto pcevent : PC_Events) { bool PCSX3TimeCut = false; @@ -1403,8 +1407,8 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s plotter->Fill2D("pcz_vs_pcphi_TimeCut", 600, -200, 200, 120, -360, 360, pcevent.pos.Z(), pcevent.pos.Phi() * 180 / M_PI, "Z_Reconstruction"); } - double sx3rho = 88.0; - double sx3z = sx3event.pos.Z(); + double sx3rho = 88.0; + double sx3z = sx3event.pos.Z(); double pcz = pcevent.pos.Z(); double calcsx3theta = TMath::ATan2(sx3rho - z_to_crossover_rho(pcz), sx3z - pcz); plotter->Fill2D("dE2_E_Anodesx3B", 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1 * TMath::Sin(calcsx3theta), "PID_dE_E"); @@ -1423,11 +1427,12 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s plotter->Fill1D("VertexReconZ_SX3", 600, -1300, 1300, vector_closest_to_z_sx3.Z(), "Vertex_Reconstruction"); plotter->Fill2D("VertexReconXY_SX3" + std::to_string(PCSX3TimeCut), 100, -100, 100, 100, -100, 100, vector_closest_to_z_sx3.X(), vector_closest_to_z_sx3.Y(), "Vertex_Reconstruction"); - plotter->Fill2D("pcz_vs_time", 2000, 0, 2000, 600, -200, 200, pcevent.Time1 * 1e-9, pcevent.pos.Z(), "Timing"); - plotter->Fill2D("pcphi_vs_time", 2000, 0, 2000, 180, -360, 360, pcevent.Time1 * 1e-9, pcevent.pos.Phi() * 180. / M_PI, "Timing"); - plotter->Fill2D("sx3phi_vs_time", 2000, 0, 2000, 180, -360, 360, pcevent.Time1 * 1e-9, sx3event.pos.Phi() * 180. / M_PI, "Timing"); + plotter->Fill2D("pcz_vs_time", 2000, 0, 2000, 600, -200, 200, pcevent.Time1 * 1e-9, pcevent.pos.Z(), "Timing"); + plotter->Fill2D("pcphi_vs_time", 2000, 0, 2000, 180, -360, 360, pcevent.Time1 * 1e-9, pcevent.pos.Phi() * 180. / M_PI, "Timing"); + plotter->Fill2D("sx3phi_vs_time", 2000, 0, 2000, 180, -360, 360, pcevent.Time1 * 1e-9, sx3event.pos.Phi() * 180. / M_PI, "Timing"); + + plotter->Fill2D("pcz_vs_sx3pczguess", 600, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z(), "Z_Reconstruction"); - plotter->Fill2D("pcz_vs_sx3pczguess", 600, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z(), "Z_Reconstruction"); if (pcevent.multi1 == 1 && pcevent.multi2 == 2) { plotter->Fill2D("pcz_vs_sx3pczguess_A1C2", 600, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z(), "Z_Reconstruction"); @@ -1466,7 +1471,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s plotter->Fill2D("pcz_vs_sx3pczguess_A1C3_strip" + std::to_string(sx3event.ch2), 300, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z(), "Z_Reconstruction"); } - plotter->Fill2D("pcz_vs_sx3pczguess_int", 600, -200, 200, 600, -200, 200, pcz_guess_int, pcevent.pos.Z(), "Z_Reconstruction"); + plotter->Fill2D("pcz_vs_sx3pczguess_int", 600, -200, 200, 600, -200, 200, pcz_guess_int, pcevent.pos.Z(), "Z_Reconstruction"); plotter->Fill2D("pcz_vs_sx3pczguess_strip" + std::to_string(sx3event.ch2), 300, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z(), "Z_Reconstruction"); bool sx3PhiCut = (TMath::Abs(sx3event.pos.Phi() - pcevent.pos.Phi()) < 45.0 * M_PI / 180.); @@ -1491,11 +1496,95 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s { plotter->Fill1D("PCZ_sx3", 800, -200, 200, pcevent.pos.Z(), "Z_Reconstruction"); } - } + + //-----------------------Benchmarking Method for Source Runs------------------------// + + if (BenchMark && aClusters.size() == 1 && cClusters.size() == 1) + { + const auto &aCl = aClusters.front(); + const auto &cCl = cClusters.front(); + + auto vertexFromPCPoint = [&x1](const TVector3 &x2f) + { + TVector3 vf = x2f - x1; + double tm = -1.0 * (x1.X() * vf.X() + x1.Y() * vf.Y()) / (vf.X() * vf.X() + vf.Y() * vf.Y()); + return TVector3(x1 + tm * vf); + }; + + auto fillSuite = [&](const std::string &tag, double pcz_method, const TVector3 &vtx) + { + plotter->Fill1D("Benchmark_SX3_VertexZ_" + tag, 800, -400, 400, vtx.Z(), "Benchmark_SX3"); + plotter->Fill1D("Benchmark_SX3_VertexZ_" + tag + "_TC" + std::to_string(PCSX3TimeCut) + "_PC" + std::to_string(phicut), 800, -400, 400, vtx.Z(), "Benchmark_SX3"); + plotter->Fill2D("Benchmark_SX3_VertexXY_" + tag, 200, -100, 100, 200, -100, 100, vtx.X(), vtx.Y(), "Benchmark_SX3"); + plotter->Fill1D("Benchmark_SX3_PCZ_" + tag, 600, -200, 200, pcz_method, "Benchmark_SX3"); + plotter->Fill2D("Benchmark_SX3_PCZ_vs_sx3pczguess_" + tag, 600, -200, 200, 600, -200, 200, pcz_guess_int, pcz_method, "Benchmark_SX3"); + plotter->Fill2D("Benchmark_SX3_PCZresidual_vs_pczguess_" + tag, 600, -200, 200, 200, -100, 100, pczguess, pcz_method - pczguess, "Benchmark_SX3"); + plotter->Fill1D("Benchmark_SX3_PCZ-sx3pczguess_" + tag, 200, -100, 100, pcz_method - pczguess, "Benchmark_SX3"); + plotter->Fill1D("Benchmark_SX3_PCZ-sx3pczint_" + tag, 200, -100, 100, pcz_method - pcz_guess_int, "Benchmark_SX3"); + }; + + auto fillVsRef = [&](const std::string &tag, double pcz_method, const TVector3 &vtx, double pcz_ref, const TVector3 &vtx_ref) + { + plotter->Fill2D("Benchmark_SX3_PCZ_" + tag + "_vs_ref", 400, -200, 200, 400, -200, 200, pcz_ref, pcz_method, "Benchmark_SX3_ref"); + plotter->Fill1D("Benchmark_SX3_PCZ_" + tag + "_minus_ref", 400, -100, 100, pcz_method - pcz_ref, "Benchmark_SX3_ref"); + plotter->Fill2D("Benchmark_SX3_VertexZ_" + tag + "_vs_ref", 400, -200, 200, 400, -200, 200, vtx_ref.Z(), vtx.Z(), "Benchmark_SX3_ref"); + plotter->Fill1D("Benchmark_SX3_VertexZ_" + tag + "_minus_ref", 400, -100, 100, vtx.Z() - vtx_ref.Z(), "Benchmark_SX3_ref"); + }; + + // cathode charge-division reference on this event (the A1C2 recipe) + double pcz_ref = pcfix_func.Eval(pcevent.pos.Z()); + TVector3 vtx_ref = vertexFromPCPoint(TVector3(pcevent.pos.X(), pcevent.pos.Y(), pcz_ref)); + + // anode-only PC point, oneWire method: pseudo-wire interpolated to the + // SX3 hit's phi, with the same quality cuts as miscHistograms_oneWire + auto [apwire_bm, apSumE_bm, apMaxE_bm, apTSMaxE_bm] = pwinstance.GetPseudoWire(aCl, "ANODE"); + TVector3 pc_anodeOnly = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, sx3event.pos.Phi()); + TVector3 vtx_anodeOnly = vertexFromPCPoint(pc_anodeOnly); + bool anodeOnlyGood = vtx_anodeOnly.Perp() <= 6.0 && vtx_anodeOnly.Z() >= -173.6 && vtx_anodeOnly.Z() <= 100 && phicut; + + if (pcevent.multi1 == 1 && pcevent.multi2 == 2) + { + fillSuite("A1C2", pcz_ref, vtx_ref); // baseline with identical binning + + // ---- A1C1 emulation: crossover with only the max-E cathode wire (Cmax), + // used directly with no z-model ---- + auto cMaxWire = *std::max_element(cCl.begin(), cCl.end(), + [](const std::tuple &w1, const std::tuple &w2) + { return std::get<1>(w1) < std::get<1>(w2); }); + std::vector> cOne = {cMaxWire}; + auto [xo_a1c1, alpha_a1c1, apSumE1, cpSumE1, apMaxE1, cpMaxE1, apTSMaxE1, cpTSMaxE1] = pwinstance.FindCrossoverProperties(aCl, cOne); + if (alpha_a1c1 != 9999999 && apSumE1 != -1) + { + TVector3 vtx_a1c1 = vertexFromPCPoint(xo_a1c1); + fillSuite("A1C1", xo_a1c1.Z(), vtx_a1c1); + fillVsRef("A1C1", xo_a1c1.Z(), vtx_a1c1, pcz_ref, vtx_ref); + } + + // ---- A1C0 emulation: oneWire method, cathodes ignored ---- + if (anodeOnlyGood) + { + fillSuite("A1C0", pc_anodeOnly.Z(), vtx_anodeOnly); + fillVsRef("A1C0", pc_anodeOnly.Z(), vtx_anodeOnly, pcz_ref, vtx_ref); + } + } + + // ---- A2C0 emulation: 2-wire anode pseudo-wire, oneWire method ---- + if (pcevent.multi1 == 2 && pcevent.multi2 == 2) + { + fillSuite("A2C2", pcz_ref, vtx_ref); // reference population for A2C0 + if (anodeOnlyGood) + { + fillSuite("A2C0", pc_anodeOnly.Z(), vtx_anodeOnly); + fillVsRef("A2C0", pc_anodeOnly.Z(), vtx_anodeOnly, pcz_ref, vtx_ref); + } + } + } + } } } -void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events) +void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events, + const std::vector>> &aClusters, const std::vector>> &cClusters) { for (auto pcevent : PC_Events) { @@ -1505,7 +1594,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s plotter->Fill2D("dt_pcA_qqqR_vs_qqqRE", 640, -2000, 2000, 400, 0, 30, qqqevent.Time1 - pcevent.Time1, qqqevent.Energy1, "Timing"); plotter->Fill1D("dt_pcC_qqqW", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2, "Timing"); plotter->Fill2D("phiPC_vs_phiQQQ", 180, -360, 360, 180, -360, 360, qqqevent.pos.Phi() * 180 / M_PI, pcevent.pos.Phi() * 180 / M_PI, "Kinematics_Angles"); - double sinTheta = TMath::Sin((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()); + double sinTheta = TMath::Sin((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()); TVector3 x2(pcevent.pos); TVector3 x1(qqqevent.pos); @@ -1524,7 +1613,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s plotter->Fill1D("dt_pcC_qqqW_pidlow_PC1", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2, "Timing"); } if (timecut) - { + { plotter->Fill2D("dE_E_AnodeQQQR", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E"); plotter->Fill2D("dE_E_CathodeQQQR", 400, 0, 30, 800, 0, 10000, qqqevent.Energy2, pcevent.Energy2, "PID_dE_E"); if (pcevent.multi1 == 1 && pcevent.multi2 == 2) @@ -1546,7 +1635,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s plotter->Fill2D("dE2_E_AnodeQQQR_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 4000, qqqevent.Energy1, pcevent.Energy1 * sinTheta, "PID_dE_E"); plotter->Fill2D("dE2_E_CathodeQQQR_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 1000, qqqevent.Energy2, pcevent.Energy2 * sinTheta, "PID_dE_E"); plotter->Fill2D("dEC_vs_dEA_TC1_PC" + std::to_string(phicut), 800, 0, 40000, 800, 0, 10000, pcevent.Energy1, pcevent.Energy2, "PID_dE_E"); - plotter->Fill2D("qqqphi_vs_time", 2000, 0, 2000, 180, -360, 360, pcevent.Time1 * 1e-9, qqqevent.pos.Phi() * 180. / M_PI, "Timing"); + plotter->Fill2D("qqqphi_vs_time", 2000, 0, 2000, 180, -360, 360, pcevent.Time1 * 1e-9, qqqevent.pos.Phi() * 180. / M_PI, "Timing"); plotter->Fill1D("dt_pcA_qqqR_timecut", 640, -2000, 2000, qqqevent.Time1 - pcevent.Time1, "Timing"); plotter->Fill1D("dt_pcC_qqqW_timecut", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2, "Timing"); @@ -1622,10 +1711,89 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s plotter->Fill2D("dE3_Ef_CathodeQQQR_TC1PC" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath), 600, 0, 15, 800, 0, 10000, qqqEfix, pcevent.Energy2 * sinTheta_customV, "PID_dE_E"); } - if(BenchMark) - { - - } + //-----------------------Benchmarking Method for Source Runs------------------------// + + if (BenchMark && aClusters.size() == 1 && cClusters.size() == 1) + { + const auto &aCl = aClusters.front(); + const auto &cCl = cClusters.front(); + + auto vertexFromPCPoint = [&x1](const TVector3 &x2f) + { + TVector3 vf = x2f - x1; + double tm = -1.0 * (x1.X() * vf.X() + x1.Y() * vf.Y()) / (vf.X() * vf.X() + vf.Y() * vf.Y()); + return TVector3(x1 + tm * vf); + }; + + auto fillSuite = [&](const std::string &tag, double pcz_method, const TVector3 &vtx) + { + plotter->Fill1D("Benchmark_QQQ_VertexZ_" + tag, 800, -400, 400, vtx.Z(), "Benchmark_QQQ"); + plotter->Fill1D("Benchmark_QQQ_VertexZ_" + tag + "_TC" + std::to_string(timecut) + "_PC" + std::to_string(phicut), 800, -400, 400, vtx.Z(), "Benchmark_QQQ"); + plotter->Fill2D("Benchmark_QQQ_VertexXY_" + tag, 200, -100, 100, 200, -100, 100, vtx.X(), vtx.Y(), "Benchmark_QQQ"); + plotter->Fill1D("Benchmark_QQQ_PCZ_" + tag, 600, -200, 200, pcz_method, "Benchmark_QQQ"); + plotter->Fill2D("Benchmark_QQQ_PCZ_vs_qqqpczguess_" + tag, 600, -200, 200, 600, -200, 200, pcz_guess_int, pcz_method, "Benchmark_QQQ"); + plotter->Fill2D("Benchmark_QQQ_PCZresidual_vs_pczguess_" + tag, 600, -200, 200, 200, -100, 100, pcz_guess_37, pcz_method - pcz_guess_37, "Benchmark_QQQ"); + plotter->Fill1D("Benchmark_QQQ_PCZ-qqqpczguess_" + tag, 200, -100, 100, pcz_method - pcz_guess_37, "Benchmark_QQQ"); + plotter->Fill1D("Benchmark_QQQ_PCZ-qqqpczint_" + tag, 200, -100, 100, pcz_method - pcz_guess_int, "Benchmark_QQQ"); + }; + + auto fillVsRef = [&](const std::string &tag, double pcz_method, const TVector3 &vtx, double pcz_ref, const TVector3 &vtx_ref) + { + plotter->Fill2D("Benchmark_QQQ_PCZ_" + tag + "_vs_ref", 400, -200, 200, 400, -200, 200, pcz_ref, pcz_method, "Benchmark_QQQ_ref"); + plotter->Fill1D("Benchmark_QQQ_PCZ_" + tag + "_minus_ref", 400, -100, 100, pcz_method - pcz_ref, "Benchmark_QQQ_ref"); + plotter->Fill2D("Benchmark_QQQ_VertexZ_" + tag + "_vs_ref", 400, -200, 200, 400, -200, 200, vtx_ref.Z(), vtx.Z(), "Benchmark_QQQ_ref"); + plotter->Fill1D("Benchmark_QQQ_VertexZ_" + tag + "_minus_ref", 400, -100, 100, vtx.Z() - vtx_ref.Z(), "Benchmark_QQQ_ref"); + }; + + // cathode charge-division reference on this event (the A1C2 recipe) + double pcz_ref = pcfix_func.Eval(pcevent.pos.Z()); + TVector3 vtx_ref = vertexFromPCPoint(TVector3(pcevent.pos.X(), pcevent.pos.Y(), pcz_ref)); + + // anode-only PC point, oneWire method: pseudo-wire interpolated to the + // QQQ hit's phi, with the same quality cuts as miscHistograms_oneWire + auto [apwire_bm, apSumE_bm, apMaxE_bm, apTSMaxE_bm] = pwinstance.GetPseudoWire(aCl, "ANODE"); + TVector3 pc_anodeOnly = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, qqqevent.pos.Phi()); + TVector3 vtx_anodeOnly = vertexFromPCPoint(pc_anodeOnly); + bool anodeOnlyGood = vtx_anodeOnly.Perp() <= 6.0 && vtx_anodeOnly.Z() >= -173.6 && vtx_anodeOnly.Z() <= 100 && phicut; + + if (pcevent.multi1 == 1 && pcevent.multi2 == 2) + { + fillSuite("A1C2", pcz_ref, vtx_ref); // baseline with identical binning + + // ---- A1C1 emulation: crossover with only the max-E cathode wire (Cmax), + // used directly with no z-model ---- + auto cMaxWire = *std::max_element(cCl.begin(), cCl.end(), + [](const std::tuple &w1, const std::tuple &w2) + { return std::get<1>(w1) < std::get<1>(w2); }); + std::vector> cOne = {cMaxWire}; + auto [xo_a1c1, alpha_a1c1, apSumE1, cpSumE1, apMaxE1, cpMaxE1, apTSMaxE1, cpTSMaxE1] = pwinstance.FindCrossoverProperties(aCl, cOne); + if (alpha_a1c1 != 9999999 && apSumE1 != -1) + { + TVector3 vtx_a1c1 = vertexFromPCPoint(xo_a1c1); + fillSuite("A1C1", xo_a1c1.Z(), vtx_a1c1); + fillVsRef("A1C1", xo_a1c1.Z(), vtx_a1c1, pcz_ref, vtx_ref); + } + + // ---- A1C0 emulation: oneWire method, cathodes ignored ---- + if (anodeOnlyGood) + { + fillSuite("A1C0", pc_anodeOnly.Z(), vtx_anodeOnly); + fillVsRef("A1C0", pc_anodeOnly.Z(), vtx_anodeOnly, pcz_ref, vtx_ref); + } + } + + // ---- A2C0 emulation: 2-wire anode pseudo-wire, oneWire method ---- + if (pcevent.multi1 == 2 && pcevent.multi2 == 2) + { + fillSuite("A2C2", pcz_ref, vtx_ref); // reference population for A2C0 + if (anodeOnlyGood) + { + fillSuite("A2C0", pc_anodeOnly.Z(), vtx_anodeOnly); + fillVsRef("A2C0", pc_anodeOnly.Z(), vtx_anodeOnly, pcz_ref, vtx_ref); + } + } + } + double qqqrho = qqqevent.pos.Perp(); double qqqz = (qqqevent.pos - TVector3(0, 0, source_vertex)).Z(); double tan_theta = qqqrho / qqqz; @@ -1644,7 +1812,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s 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(), "Z_Reconstruction"); } } - } + } } // We pass plotter as a pointer, strings as const references (for speed), // and the Kinematics object as a reference. diff --git a/run_tr.sh b/run_tr.sh index bb868c8..f71630e 100644 --- a/run_tr.sh +++ b/run_tr.sh @@ -55,14 +55,15 @@ if [[ 1 -eq 1 ]]; then fi # --- Block 3: 27Al Protons+Gas Runs (15, 17-22) --- -if [[ 1 -eq 1 ]]; then +if [[ 1 -eq 0 ]]; then export DATASET="27Al" export PREFIX="Run_" export OUT_DIR="Output_p" export source_vertex=-200.0 # Source on the entrance window echo "Starting parallel processing for 27Al proton runs..." - parallel --bar -j 6 process_run ::: 15 {17..22} + process_run 17 + # parallel --bar -j 6 process_run ::: 15 {17..22} hadd -j 4 -k ${OUT_DIR}/Al_protons.root ${OUT_DIR}/results_run0{17..22}.root fi @@ -77,17 +78,17 @@ if [[ 1 -eq 0 ]]; then fi # --- Block 5: 17F Alpha Run with Gas (18-21) --- -if [[ 1 -eq 0 ]]; then +if [[ 1 -eq 1 ]]; then export DATASET="17F" export PREFIX="SourceRun_" export OUT_DIR="Output_a" echo "Processing 17F alpha runs with dynamic source vertices..." # Running sequentially since the source_vertex variable changes per run - # export source_vertex=53.44; process_run 18 + export source_vertex=53.44; process_run 18 export source_vertex=14.24; process_run 19 - # export source_vertex=-24.96; process_run 20 - # export source_vertex=-73.96; process_run 21 + export source_vertex=-24.96; process_run 20 + export source_vertex=-73.96; process_run 21 fi # --- Block 6: 17F Proton Data --- diff --git a/scratch/plot_slope_scan.C b/scratch/plot_slope_scan.C index f726801..3de0fac 100644 --- a/scratch/plot_slope_scan.C +++ b/scratch/plot_slope_scan.C @@ -25,9 +25,11 @@ void plot_slope_scan() // 0.82, 0.84, 0.86, 0.88, 0.90}; // const int N = 25; - double slopes[] = {0.42, 0.44, 0.46, 0.48, 0.50, 0.52, 0.54, 0.56, 0.58, 0.60, 0.62, 0.64}; - const int N = 12; // Updated to 14 to match the size of your slopes array - + // double slopes[] = {0.42, 0.44, 0.46, 0.48, 0.50, 0.52, 0.54, 0.56, 0.58, 0.60, 0.62, 0.64}; + // const int N = 12; // Updated to 14 to match the size of your slopes array + const int N = 1; // Updated to 14 to match the size of your slopes array + double slopes[] = {0.52}; + int runs[] = {9, 12, 18, 19, 20, 21}; const int NRUNS = 6; // int runs[] = {9, 12}; @@ -272,12 +274,12 @@ void plot_slope_scan() mg1->Draw("A"); double bestSlope = 0, bestStdDev = 1e9; - if (gStdDevRun_sx3[0]->GetN() > 0) + if (gStdDevRun_sx3[3]->GetN() > 0) { - for (int i = 0; i < gStdDevRun_sx3[0]->GetN(); i++) + for (int i = 0; i < gStdDevRun_sx3[3]->GetN(); i++) { double x, y; - gStdDevRun_sx3[0]->GetPoint(i, x, y); + gStdDevRun_sx3[3]->GetPoint(i, x, y); if (y > 0 && y < bestStdDev) { bestStdDev = y; @@ -288,12 +290,12 @@ void plot_slope_scan() lb->SetLineColor(kRed); lb->SetLineStyle(kDashed); lb->Draw(); - printf(">>> Best slope (Run %d) = %.2f StdDev=%.2f mm\n", runs[0], bestSlope, bestStdDev); + printf(">>> Best slope (Run %d) = %.2f StdDev=%.2f mm\n", runs[3], bestSlope, bestStdDev); } // --- Pad 2: Mean --- pad1_grid->cd(2); - gPad->SetGrid(1, 1); + // gPad->SetGrid(1, 1); TMultiGraph *mg2 = new TMultiGraph("mg_mean", "Mean of Residuals vs Slope;Slope;Mean Offset (mm)"); for (int ir = 0; ir < NRUNS; ir++) @@ -358,7 +360,7 @@ void plot_slope_scan() TPad *pad2_grid = new TPad("pad2_grid", "Grid", 0.0, 0.0, 1.0, 0.92); pad2_grid->Draw(); - pad2_grid->Divide(5, 5); + pad2_grid->Divide(3,4); for (int i = 0; i < N; i++) { @@ -453,7 +455,7 @@ void plot_slope_scan() TPad *pad3_grid = new TPad("pad3_grid", "Grid", 0.0, 0.0, 1.0, 0.92); pad3_grid->Draw(); - pad3_grid->Divide(5, 5); + pad3_grid->Divide(3,4); for (int i = 0; i < N; i++) {