From 3cfecd4c29886f23c248f835861b8e8dfb86c881 Mon Sep 17 00:00:00 2001 From: vsitaraman Date: Mon, 6 Jul 2026 10:31:52 -0400 Subject: [PATCH] modified: Armory/ClassPW.h modified: TrackRecon.C modified: run_tr.sh --- Armory/ClassPW.h | 37 ++++++++++++++ TrackRecon.C | 122 +++++++++++++++++++++++++++-------------------- run_tr.sh | 41 ++++++++-------- 3 files changed, 129 insertions(+), 71 deletions(-) diff --git a/Armory/ClassPW.h b/Armory/ClassPW.h index a56c6d5..0d05cb8 100755 --- a/Armory/ClassPW.h +++ b/Armory/ClassPW.h @@ -96,6 +96,8 @@ public: 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); + static constexpr double kPCPathFail = 54321.0; + inline double PCPathLength(const TVector3 &x1, const TVector3 &x2) const; void Print() { @@ -617,4 +619,39 @@ inline double PW::GetZ0() return trackVec.Z(); } +inline double PW::PCPathLength(const TVector3 &x1, const TVector3 &x2) const +{ + TVector3 dx = x2 - x1; + double t2 = 1.0; // parametric endpoint: x2 corresponds to t=1 + + // Returns the point where the ray x1+t*dx intersects the one-sheet + // hyperboloid (x²+y²)/a² - z²/c² = 1, for t in [0, t2]. + // Returns a sentinel with Z==kPCPathFail on no intersection. + auto intersect = [&](double a, double c) -> TVector3 + { + double A = dx.Perp2() / (a * a) - dx.Z() * dx.Z() / (c * c); + double B = 2.0 * (dx.X() * x1.X() + dx.Y() * x1.Y()) / (a * a) - 2.0 * (dx.Z() * x1.Z()) / (c * c); + double C = x1.Perp2() / (a * a) - x1.Z() * x1.Z() / (c * c) - 1.0; + double disc = B * B - 4.0 * A * C; + if (disc < 0.0) + return TVector3(0, 0, kPCPathFail); + double t1s = (-B + TMath::Sqrt(disc)) / (2.0 * A); + double t2s = (-B - TMath::Sqrt(disc)) / (2.0 * A); + if (t1s >= 0.0 && t1s <= t2) + return x1 + t1s * dx; + else if (t2s >= 0.0 && t2s <= t2) + return x1 + t2s * dx; + else + return TVector3(0, 0, kPCPathFail); + }; + + // Hyperboloid parameters (mm) from fits to anode/cathode crossover points. + // Cathode waist is anode waist scaled by the outermost-radius ratio 43/37. + TVector3 an = intersect(32.0429, 301.895); + TVector3 ca = intersect(37.239045, 301.895); + if (an.Z() != kPCPathFail && ca.Z() != kPCPathFail) + return (ca - an).Mag(); + return kPCPathFail; +} + #endif diff --git a/TrackRecon.C b/TrackRecon.C index 843498b..a2e83e2 100644 --- a/TrackRecon.C +++ b/TrackRecon.C @@ -41,7 +41,7 @@ Int_t colors[40] = { // --- Analysis Control Flags --- bool process_alpha_proton_scattering = false, - doMiscHistograms = false, + doMiscHistograms = true, doPCSX3ClusterAnalysis = true, doPCQQQClusterAnalysis = true, doOldAnalysis = false, @@ -581,31 +581,27 @@ void TrackRecon::Begin(TTree * /*tree*/) if (dataset == "17F" && CO2percent == 3) { - if (!MeV_to_cm) - MeV_to_cm = new TGraph("eloss_calculations/alpha_lookup_20MeV_3pc.dat", "%lf %*lf %lf"); - if (!MeV_to_cm_p) - MeV_to_cm_p = new TGraph("eloss_calculations/proton_lookup_20MeV_3pc.dat", "%lf %*lf %lf"); - if (!MeV_to_cm_27Al) - MeV_to_cm_27Al = new TGraph("eloss_calculations/aluminum_lookup_80MeV_3pc.dat", "%lf %*lf %lf"); - if (!MeV_to_cm_17F) - MeV_to_cm_17F = new TGraph("eloss_calculations/fluorine_lookup_70MeV_3pc.dat", "%lf %*lf %lf"); + // MeV_to_cm = new TGraph("eloss_calculations/alpha_lookup_20MeV_3pc.dat", "%lf %*lf %lf"); + // MeV_to_cm_p = new TGraph("eloss_calculations/proton_lookup_20MeV_3pc.dat", "%lf %*lf %lf"); + // MeV_to_cm_27Al = new TGraph("eloss_calculations/aluminum_lookup_80MeV_3pc.dat", "%lf %*lf %lf"); + // MeV_to_cm_17F = new TGraph("eloss_calculations/fluorine_lookup_70MeV_3pc.dat", "%lf %*lf %lf"); - // MeV_to_cm = new TGraph("eloss_calculations/alpha_lookup_20MeV_3pc_350Torr.dat", "%lf %*lf %lf"); - // MeV_to_cm_p = new TGraph("eloss_calculations/proton_lookup_20MeV_3pc_350Torr.dat", "%lf %*lf %lf"); - // MeV_to_cm_17F = new TGraph("eloss_calculations/fluorine_lookup_70MeV_3pc_350Torr.dat", "%lf %*lf %lf"); - // MeV_to_cm_27Al = new TGraph("eloss_calculations/aluminum_lookup_80MeV_3pc_350Torr.dat", "%lf %*lf %lf"); + MeV_to_cm = new TGraph("eloss_calculations/alpha_lookup_20MeV_3pc_350Torr.dat", "%lf %*lf %lf"); + MeV_to_cm_p = new TGraph("eloss_calculations/proton_lookup_20MeV_3pc_350Torr.dat", "%lf %*lf %lf"); + MeV_to_cm_17F = new TGraph("eloss_calculations/fluorine_lookup_70MeV_3pc_350Torr.dat", "%lf %*lf %lf"); + MeV_to_cm_27Al = new TGraph("eloss_calculations/aluminum_lookup_80MeV_3pc_350Torr.dat", "%lf %*lf %lf"); } else if (dataset == "17F" && CO2percent == 4) { - MeV_to_cm = new TGraph("eloss_calculations/alpha_lookup_20MeV_4pc.dat", "%lf %*lf %lf"); - MeV_to_cm_p = new TGraph("eloss_calculations/proton_lookup_20MeV_4pc.dat", "%lf %*lf %lf"); - MeV_to_cm_27Al = new TGraph("eloss_calculations/aluminum_lookup_80MeV_4pc.dat", "%lf %*lf %lf"); - MeV_to_cm_17F = new TGraph("eloss_calculations/fluorine_lookup_70MeV_4pc.dat", "%lf %*lf %lf"); + // MeV_to_cm = new TGraph("eloss_calculations/alpha_lookup_20MeV_4pc.dat", "%lf %*lf %lf"); + // MeV_to_cm_p = new TGraph("eloss_calculations/proton_lookup_20MeV_4pc.dat", "%lf %*lf %lf"); + // MeV_to_cm_27Al = new TGraph("eloss_calculations/aluminum_lookup_80MeV_4pc.dat", "%lf %*lf %lf"); + // MeV_to_cm_17F = new TGraph("eloss_calculations/fluorine_lookup_70MeV_4pc.dat", "%lf %*lf %lf"); - // MeV_to_cm = new TGraph("eloss_calculations/alpha_lookup_20MeV_4pc_350Torr.dat", "%lf %*lf %lf"); - // MeV_to_cm_p = new TGraph("eloss_calculations/proton_lookup_20MeV_4pc_350Torr.dat", "%lf %*lf %lf"); - // MeV_to_cm_17F = new TGraph("eloss_calculations/fluorine_lookup_70MeV_4pc_350Torr.dat", "%lf %*lf %lf"); - // MeV_to_cm_27Al = new TGraph("eloss_calculations/aluminum_lookup_80MeV_4pc_350Torr.dat", "%lf %*lf %lf"); + MeV_to_cm = new TGraph("eloss_calculations/alpha_lookup_20MeV_4pc_350Torr.dat", "%lf %*lf %lf"); + MeV_to_cm_p = new TGraph("eloss_calculations/proton_lookup_20MeV_4pc_350Torr.dat", "%lf %*lf %lf"); + MeV_to_cm_17F = new TGraph("eloss_calculations/fluorine_lookup_70MeV_4pc_350Torr.dat", "%lf %*lf %lf"); + MeV_to_cm_27Al = new TGraph("eloss_calculations/aluminum_lookup_80MeV_4pc_350Torr.dat", "%lf %*lf %lf"); } else { @@ -1422,6 +1418,8 @@ Bool_t TrackRecon::Process(Long64_t entry) #endif double path_length = (sx3event.pos - TVector3(0, 0, vertex_recon)).Mag() * 0.1; + double _pc_sx3 = pwinstance.PCPathLength(TVector3(0, 0, vertex_recon), sx3event.pos); + double path_length = (_pc_sx3 < PW::kPCPathFail ? _pc_sx3 : (sx3event.pos - TVector3(0, 0, vertex_recon)).Mag()) * 0.1; double sx3Efix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(sx3event.Energy1) - path_length); double sx3Efixalpha = cm_to_MeV->Eval(MeV_to_cm->Eval(sx3event.Energy1) - path_length); @@ -1521,6 +1519,8 @@ Bool_t TrackRecon::Process(Long64_t entry) // Energy Loss Correction double path_length = (qqqevent.pos - TVector3(0, 0, vertex_recon)).Mag() * 0.1; + double _pc_qqq = pwinstance.PCPathLength(TVector3(0, 0, vertex_recon), qqqevent.pos); + double path_length = (_pc_qqq < PW::kPCPathFail ? _pc_qqq : (qqqevent.pos - TVector3(0, 0, vertex_recon)).Mag()) * 0.1; double qqqEfix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(qqqevent.Energy1) - path_length); double qqqEfixalpha = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy1) - path_length); @@ -1628,8 +1628,13 @@ void protonAlphaHistograms(HistPlotter *plotter, const std::vector &QQQ_E plotter->Fill1D("ap_pczfix", 600, -300, 300, pcz_fix, aplabel); plotter->Fill1D("ap_pcz", 600, -300, 300, pcevent.pos.Z(), aplabel); - double path_length_q = (qqqevent.pos - TVector3(0, 0, vertex_z)).Mag() * 0.1; - double path_length_s = (sx3event.pos - TVector3(0, 0, vertex_z)).Mag() * 0.1; + double _pc_q0 = pwinstance.PCPathLength(TVector3(0, 0, vertex_z), qqqevent.pos); + double _pc_s0 = pwinstance.PCPathLength(TVector3(0, 0, vertex_z), sx3event.pos); + double path_length_q = (_pc_q0 < PW::kPCPathFail ? _pc_q0 : (qqqevent.pos - TVector3(0, 0, vertex_z)).Mag()) * 0.1; + double path_length_s = (_pc_s0 < PW::kPCPathFail ? _pc_s0 : (sx3event.pos - TVector3(0, 0, vertex_z)).Mag()) * 0.1; + + // double path_length_q = (qqqevent.pos - TVector3(0, 0, vertex_z)).Mag() * 0.1; + // double path_length_s = (sx3event.pos - TVector3(0, 0, vertex_z)).Mag() * 0.1; // double path_length_q = (qqqevent.pos-r_rhoMin_fix).Mag()*0.1; // double path_length_s = (sx3event.pos-r_rhoMin_fix).Mag()*0.1; @@ -1858,7 +1863,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector &QQQ_Ev // plotter->Fill2D("pcdEA_vs_sx3z", 300, 0, 200, 800, 0, 20000, sx3z, pcevent.Energy1, "Kinematics_Angles"); // plotter->Fill2D("pcdEA_vs_sx3pczguess", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy1, "Kinematics_Angles"); - plotter->Fill2D("pcdEA_vs_pczfix", 600, -200, 200, 800, 0, 20000, pcz_fix, pcevent.Energy1, "Kinematics_Angles"); + plotter->Fill2D("pcdEA_vs_pczfix", 600, -200, 200, 800, 0, 20000, pcz_fix, pcevent.Energy1, "PCdE_vs_Z"); // plotter->Fill2D("pcdEC_vs_sx3z", 300, 0, 200, 800, 0, 20000, sx3z, pcevent.Energy2, "Kinematics_Angles"); // plotter->Fill2D("pcdEC_vs_sx3pczguess", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy2, "Kinematics_Angles"); plotter->Fill2D("pcdEC_vs_pczfix", 600, -200, 200, 800, 0, 20000, pcz_fix, pcevent.Energy2, "PCdE_vs_Z"); @@ -1877,7 +1882,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector &QQQ_Ev // plotter->Fill2D("pcdEA_vs_sx3pczguess_A" + std::to_string(pcevent.multi1) + "C" + std::to_string(pcevent.multi2), 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy1, "PCdE_vs_Z"); plotter->Fill2D("pcdEA_vs_sx3pczguess", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy1, "PCdE_vs_Z"); - // plotter->Fill2D("pcdEA_vs_pczguess", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy1, "PCdE_vs_Z"); + plotter->Fill2D("pcdEA_vs_pczguess", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy1, "PCdE_vs_Z"); // plotter->Fill2D("pcdEA_vs_pczfix" + std::to_string(pcevent.multi1) + "A" + std::to_string(pcevent.multi1) + "C", 600, -200, 200, 800, 0, 20000, pcz_fix, pcevent.Energy1, "PCdE_vs_Z"); // plotter->Fill2D("pcdEC_vs_sx3pczguess_A" + std::to_string(pcevent.multi1) + "C" + std::to_string(pcevent.multi2), 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy2, "PCdE_vs_Z"); if (pcevent.multi1 == 1) @@ -1888,6 +1893,10 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector &QQQ_Ev { plotter->Fill2D("pcdEC_vs_sx3pczguess_C1", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy2, "PCdE_vs_Z"); } + if (pcevent.multi2 == 2) + { + plotter->Fill2D("pcdEC_vs_sx3pczguess_C2", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy2, "PCdE_vs_Z"); + } plotter->Fill2D("pcdEC_vs_sx3pczguess", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy2, "PCdE_vs_Z"); plotter->Fill2D("pcdEC_vs_pczguess", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy2, "PCdE_vs_Z"); plotter->Fill2D("pcdEACSum_vs_pczguess", 600, -200, 200, 800, 0, 20000, pczguess, (pcevent.Energy2 + pcevent.Energy1) / 2, "PCdE_vs_Z"); @@ -2130,10 +2139,10 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector &QQQ_Ev double zc = 0.5 * (a1c1_zg[cell_truth] + a1c1_zg[cell_truth + 1]); double half = 0.5 * (a1c1_zg[cell_truth] - a1c1_zg[cell_truth + 1]); - plotter->Fill2D("AnodeEnergy_vs_CellSX3", 120, 0, 1.2, 100, 0, 40000, 1 - TMath::Abs(pcz_ref - zc) / half, pcevent.Energy1); + plotter->Fill2D("AnodeEnergy_vs_CellSX3", 120, 0, 1.2, 800, 0, 40000, 1 - TMath::Abs(pcz_ref - zc) / half, pcevent.Energy1); plotter->Fill2D("CathodeEnergy_vs_CellSX3", 120, 0, 1.2, 800, 0, 40000, TMath::Abs(pcz_ref - zc) / half, pcevent.Energy2); - plotter->Fill2D("FracEnergy_vs_CellSX3", 120, 0, 1.2, 1200, 0, 20, TMath::Abs(pcz_ref - zc) / half, pcevent.Energy2 / pcevent.Energy1); - plotter->Fill2D("SumEnergy_vs_CellSX3", 120, 0, 1.2, 1200, 0, 20, TMath::Abs(pcz_ref - zc) / half, (pcevent.Energy1 + pcevent.Energy2) / 2); + plotter->Fill2D("FracEnergy_vs_CellSX3", 120, 0, 1.2, 800, 0, 10, TMath::Abs(pcz_ref - zc) / half, pcevent.Energy2 / pcevent.Energy1); + plotter->Fill2D("SumEnergy_vs_CellSX3", 120, 0, 1.2, 800, 0, 10, TMath::Abs(pcz_ref - zc) / half, (pcevent.Energy1 + pcevent.Energy2) / 2); if (half > 0.0) { @@ -2300,8 +2309,8 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector &QQQ_Ev for (const auto &pcevent : PC_Events) { - plotter->Fill2D("pcdEACSum_vs_anodechannel", 24, 0, 23, 800, 0, 20000, pcevent.ch1, (pcevent.Energy1 + pcevent.Energy2) / 2); - plotter->Fill2D("pcdEACSum_vs_cathodechannel", 24, 0, 23, 800, 0, 20000, pcevent.ch2, (pcevent.Energy1 + pcevent.Energy2) / 2); + plotter->Fill2D("pcdEACAvg_vs_anodechannel", 24, 0, 23, 800, 0, 20000, pcevent.Anodech, (pcevent.Energy1 + pcevent.Energy2) / 2); + plotter->Fill2D("pcdEACAvg_vs_cathodechannel", 24, 0, 23, 800, 0, 20000, pcevent.Cathodech, (pcevent.Energy1 + pcevent.Energy2) / 2); for (const auto &qqqevent : QQQ_Events) { @@ -2432,7 +2441,9 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector &QQQ_Ev plotter->Fill2D("pcdEACSum_vs_pczfix", 600, -200, 200, 800, 0, 20000, pcz_fix, (pcevent.Energy1 + pcevent.Energy2) / 2, "PCdE_vs_Z"); - double path_length = (qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Mag() * 0.1; + // double path_length = (qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Mag() * 0.1; + double _pc_bq = pwinstance.PCPathLength(r_rhoMin_fix, qqqevent.pos); + double path_length = (_pc_bq < PW::kPCPathFail ? _pc_bq : (qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Mag()) * 0.1; double qqqEfix = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy1) - path_length); double qqqEfix_p = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(qqqevent.Energy1) - path_length); @@ -2447,6 +2458,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector &QQQ_Ev plotter->Fill2D("dE3_Ef_AnodeQQQR_TC1" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath), 600, 0, 15, 800, 0, 40000, qqqEfix, pcevent.Energy1 * sinTheta_customV, "PID_dE_E"); 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"); + plotter->Fill2D("pcdEACAvg_A1C2_vs_qqqE", 400, 0, 10, 1600, 0, 40000, qqqevent.Energy1, (pcevent.Energy1 + pcevent.Energy2) / 2); } // plotter->Fill2D("pcdEA_vs_qqqpczguess_A" + std::to_string(pcevent.multi1) + "C" + std::to_string(pcevent.multi2), 600, -200, 200, 800, 0, 20000, pcz_guess_37, pcevent.Energy1, "PCdE_vs_Z"); @@ -2457,24 +2469,10 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector &QQQ_Ev plotter->Fill2D("pcdEC_vs_qqqpczguess", 600, -200, 200, 800, 0, 20000, pcz_guess_37, pcevent.Energy2, "PCdE_vs_Z"); plotter->Fill2D("pcdEC_vs_pczguess", 600, -200, 200, 800, 0, 20000, pcz_guess_37, pcevent.Energy2, "PCdE_vs_Z"); - plotter->Fill2D("pcdEACSum_vs_qqqpczguess", 600, -200, 200, 800, 0, 20000, pcz_guess_37, (pcevent.Energy1 + pcevent.Energy2) / 2, "PCdE_vs_Z"); - plotter->Fill2D("pcdEACSum_vs_pczguess", 600, -200, 200, 800, 0, 20000, pcz_guess_37, (pcevent.Energy1 + pcevent.Energy2) / 2, "PCdE_vs_Z"); + plotter->Fill2D("pcdEACAvg_vs_qqqpczguess", 600, -200, 200, 800, 0, 20000, pcz_guess_37, (pcevent.Energy1 + pcevent.Energy2) / 2, "PCdE_vs_Z"); + plotter->Fill2D("pcdEACAvg_vs_pczguess", 600, -200, 200, 800, 0, 20000, pcz_guess_37, (pcevent.Energy1 + pcevent.Energy2) / 2, "PCdE_vs_Z"); - // A1CmaxC low-band corrected PC energy sum vs QQQ energy. - { - double ac_sum_ev = pcevent.Energy1 + pcevent.Energy2; - if (ac_sum_ev > 0.0) - { - double cfrac_ev = pcevent.Energy2 / ac_sum_ev; - bool is_lowband = (cfrac_ev < a1c1_cfrac_split); - double pcEsum_corr_cmaxc = (is_lowband && a1c1_lowband_rfactor > 0.0) - ? pcevent.Energy1 + pcevent.Energy2 * a1c1_lowband_rfactor - : ac_sum_ev; - plotter->Fill2D("pcdEACSum_cmaxc_corr_vs_qqqE", 400, 0, 10, 800, 0, 40000, qqqevent.Energy1, pcEsum_corr_cmaxc); - plotter->Fill2D("pcdEA_vs_qqqE", 400, 0, 10, 800, 0, 40000, qqqevent.Energy1, pcEsum_corr_cmaxc); - plotter->Fill2D("pcdEACSum_cmaxc_raw_vs_qqqE", 400, 0, 10, 800, 0, 40000, qqqevent.Energy1, ac_sum_ev); - } - } + plotter->Fill2D("pcdEACAvg_vs_qqqE", 400, 0, 10, 1600, 0, 40000, qqqevent.Energy1, (pcevent.Energy1 + pcevent.Energy2) / 2); // plotter->Fill2D("pcdEC_vs_pczfix" + std::to_string(pcevent.multi1) + "A" + std::to_string(pcevent.multi1) + "C", 800, 0, 20000, 600, -200, 200, pcevent.Energy2, pcz_fix, "PCdE_vs_Z"); @@ -3268,7 +3266,9 @@ void miscHistograms_oneWire(HistPlotter *plotter, const std::vector &QQQ_ plotter->Fill1D("vertexZ1d_ainterp_qqqZ_TC1_ignC_a" + std::to_string(acluster.size()), 800, -400, 400, r_rhoMin_fix.Z(), "ainterp_noc"); plotter->Fill2D("vertexXY_ainterp_TC1_ignC_a" + std::to_string(acluster.size()), 200, -100, 100, 200, -100, 100, r_rhoMin_fix.X(), r_rhoMin_fix.Y(), "ainterp_noc"); - double path_length_q = (qqqevent.pos - r_rhoMin_fix).Mag() * 0.1; + // double path_length_q = (qqqevent.pos - r_rhoMin_fix).Mag() * 0.1; + double _pc_q1 = pwinstance.PCPathLength(r_rhoMin_fix, qqqevent.pos); + double path_length_q = (_pc_q1 < PW::kPCPathFail ? _pc_q1 : (qqqevent.pos - r_rhoMin_fix).Mag()) * 0.1; double qqqEfix; qqqEfix = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy1) - path_length_q); plotter->Fill1D("pmisc_ow_Ex_from_alpha", 200, -10, 10, apkin_a.getExc(qqqEfix, theta_q * 180 / M_PI), "ainterp_noc"); @@ -3283,7 +3283,8 @@ void protonMiscHistograms(HistPlotter *plotter, const std::vector &QQQ_Ev { // consider the 'proton-like' QQQ branch seen in a,p data static TRandom3 rand(0); // seeded once (random seed via TUUID), not per call - double initial_energy = 6.78; + double initial_energy = 6.89; //only Kapton + // double initial_energy = 6.78; // Kapton plus 100mm He-CO2 // if (dataset == "27Al") // initial_energy = 6.79; // m3 is alpha, 6.79 MeV is 7.0 MeV proton energy after kapton+100mm 4He gas (molar mass 5.2, 1 gain) // if (dataset == "17F") @@ -3395,6 +3396,12 @@ void protonMiscHistograms(HistPlotter *plotter, const std::vector &QQQ_Ev if (vertex_z < -173.6 || vertex_z > 100) continue; + double beam_path_length_q = TMath::Abs(vertex_z - z_entrance) * 0.1; + double beam_energy_at_vertex_q = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(initial_energy) - beam_path_length_q); + if (beam_energy_at_vertex_q <= 0.0) + beam_energy_at_vertex_q = 0.001; + Kinematics apkin_a(mass_1H, mass_4He, mass_4He, mass_1H, beam_energy_at_vertex_q); + // What's below: radial cut, time coincident, phi-correlated events with possible energy selection applied to both E_si and dE_Anodes auto plot_with_tag = [&](std::string tag = "") { @@ -3425,7 +3432,9 @@ void protonMiscHistograms(HistPlotter *plotter, const std::vector &QQQ_Ev // double path_length_q = (qqqevent.pos-TVector3(0,0,vertex_z)).Mag()*0.1; // double path_length_s = (sx3event.pos-TVector3(0,0,vertex_z)).Mag()*0.1; - double path_length_q = (qqqevent.pos - r_rhoMin_fix).Mag() * 0.1; + // double path_length_q = (qqqevent.pos - r_rhoMin_fix).Mag() * 0.1; + double _pc_q2 = pwinstance.PCPathLength(r_rhoMin_fix, qqqevent.pos); + double path_length_q = (_pc_q2 < PW::kPCPathFail ? _pc_q2 : (qqqevent.pos - r_rhoMin_fix).Mag()) * 0.1; double qqqEfix; if (tag == "_cathode_alphas") { // satisfied when find succeeds @@ -3500,6 +3509,13 @@ void protonMiscHistograms_sx3(HistPlotter *plotter, const std::vector &QQ bool cathode_alpha_select = (pcevent.Energy2 > 1400); // What's below: radial cut, time coincident, phi-correlated events with possible energy selection applied to both E_si and dE_Anodes + + double beam_path_length_s = TMath::Abs(vertex_z - z_entrance) * 0.1; + double beam_energy_at_vertex_s = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(6.89) - beam_path_length_s); + if (beam_energy_at_vertex_s <= 0.0) + beam_energy_at_vertex_s = 0.001; + Kinematics apkin_a_s(mass_1H, mass_4He, mass_4He, mass_1H, beam_energy_at_vertex_s); + auto plot_with_tag = [&](std::string tag = "") { std::string pmlabel = "proton+miscsx3" + tag; @@ -3518,7 +3534,9 @@ void protonMiscHistograms_sx3(HistPlotter *plotter, const std::vector &QQ // double path_length_q = (sx3event.pos-TVector3(0,0,vertex_z)).Mag()*0.1; // double path_length_s = (sx3event.pos-TVector3(0,0,vertex_z)).Mag()*0.1; - double path_length_s = (sx3event.pos - r_rhoMin_fix).Mag() * 0.1; + // double path_length_s = (sx3event.pos - r_rhoMin_fix).Mag() * 0.1; + double _pc_s1 = pwinstance.PCPathLength(r_rhoMin_fix, sx3event.pos); + double path_length_s = (_pc_s1 < PW::kPCPathFail ? _pc_s1 : (sx3event.pos - r_rhoMin_fix).Mag()) * 0.1; double sx3Efix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(sx3event.Energy1) - path_length_s); // plotter->Fill2D("sx3Ef_sx3E_matrix_all"+tag,400,0,10,400,0,10,sx3Efix,sx3event.Energy1,pmlabel); diff --git a/run_tr.sh b/run_tr.sh index 5437da4..aef0e61 100644 --- a/run_tr.sh +++ b/run_tr.sh @@ -41,22 +41,6 @@ if [[ 1 -eq 0 ]]; then parallel --bar -j 6 process_run ::: {1..8} fi -# --- Block 2: 27Al Alpha+Gas Runs (9, 12) --- -if [[ 1 -eq 1 ]]; then - export DATASET="27Al" - export PREFIX="Run_" - export OUT_DIR="Output_a" - export Gain=2 - export CATHODE_GAIN=3.0 - echo "Processing 27Al alpha+gas runs..." - export source_vertex=-5.36; export timecut_low=12.0; export timecut_high=119.0; process_run 9 "$slope" - unset timecut_high - export source_vertex=53.44; export timecut_low=400.0; process_run 12 "$slope" - unset Gain - unset CATHODE_GAIN - unset timecut_low -fi - # --- Block 4: 17F Source Runs (5-14) --- if [[ 1 -eq 0 ]]; then export DATASET="17F" @@ -67,8 +51,25 @@ if [[ 1 -eq 0 ]]; then parallel --bar -j 6 process_run ::: {5..13} fi -# --- Block 5: 17F Alpha Run with Gas (18-21) --- -if [[ 1 -eq 1 ]]; then +# --- Block 2: 27Al Alpha+Gas Runs (9, 12) --- +if [[ 1 -eq 0 ]]; then + export DATASET="27Al" + export PREFIX="Run_" + export OUT_DIR="Output_a" + export Gain=2 + export CATHODE_GAIN=3.0 + rm -f ${OUT_DIR}/all.root + echo "Processing 27Al alpha+gas runs..." + export source_vertex=-5.36; export timecut_low=12.0; export timecut_high=119.0; process_run 9 "$slope" + unset timecut_high + export source_vertex=53.44; export timecut_low=400.0; process_run 12 "$slope" + unset Gain + unset CATHODE_GAIN + unset timecut_low +fi + +# --- Block 5: 17F Alpha+Gas Runs (18-21) --- +if [[ 1 -eq 0 ]]; then export DATASET="17F" export PREFIX="SourceRun_" export OUT_DIR="Output_a" @@ -79,11 +80,12 @@ if [[ 1 -eq 1 ]]; then export source_vertex=14.24; process_run 19 export source_vertex=-24.96; process_run 20 export source_vertex=-73.96; process_run 21 + hadd -j 4 -k ${OUT_DIR}/all.root ${OUT_DIR}/results_run*.root # exit fi # --- Block 3: 27Al Protons+Gas Runs (15, 17-22) --- -if [[ 1 -eq 1 ]]; then +if [[ 1 -eq 0 ]]; then # export CO2percent=4 export DATASET="27Al" @@ -101,6 +103,7 @@ if [[ 1 -eq 1 ]]; then hadd -j 4 -k ${OUT_DIR}/Al_protons.root ${OUT_DIR}/results_run0{15..22}.root unset Gain unset CATHODE_GAIN + # exit fi # --- Block 6: 17F Proton Data ---