diff --git a/Armory/PC_StepLadder_Correction.h b/Armory/PC_StepLadder_Correction.h index 5fd75a7..0f5440a 100644 --- a/Armory/PC_StepLadder_Correction.h +++ b/Armory/PC_StepLadder_Correction.h @@ -1,88 +1,82 @@ #include -/*double model(double* x, double* p) { - double result = x[0]; - double factor = 29.0; - double slope = 0.40; - if(TMath::Abs(x[0]) < 16.2) result=x[0]*slope; - else if(TMath::Abs(x[0]) < 49.8 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor; - else if(TMath::Abs(x[0]) < 85.2 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*2; - else result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*3; - return result; -} -double model_invert(double *y, double *q) { - double result=y[0]; - double slope = 0.40; - double factor = 0.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*3; - return result; -}*/ - -double model_invert(double* y, double* p) { +double model_invert(double *y, double *p) +{ double result = y[0]; double slope = 0.52; - double z_grid[8] = {147.998,101.946,59.7634,19.6965,-19.6965,-59.7634,-101.946,-147.998}; - for(int i=0;i<7;i++) { - if(y[0] <= z_grid[i] && y[0] > z_grid[i+1]) { - double zavg = (z_grid[i] + z_grid[i+1])*0.5; //midpoint about which we pivot - result = (y[0]-zavg)/slope + zavg; - break; - } - } + double z_grid[8] = {147.998, 101.946, 59.7634, 19.6965, -19.6965, -59.7634, -101.946, -147.998}; + for (int i = 0; i < 7; i++) + { + if (y[0] <= z_grid[i] && y[0] > z_grid[i + 1]) + { + double zavg = (z_grid[i] + z_grid[i + 1]) * 0.5; // midpoint about which we pivot + result = (y[0] - zavg) / slope + zavg; + break; + } + } return result; } +// ---- Parabola variant (now takes optional slope; out-of-band -> Si, no clamp-snap) ---- +inline double model_invert_a1c1(double cfrac, double z_si, double z_pivot, + double C, const double* c0_table, bool& ok, + double slope = 1.0) +{ + ok = false; + static const double z_grid[8] = + {147.998, 101.946, 59.7634, 19.6965, -19.6965, -59.7634, -101.946, -147.998}; + const double cfrac_hi = 0.90, cfrac_lo = 0.12; + if (cfrac > cfrac_hi || cfrac < cfrac_lo) return z_si; -double model_a1c1(double* x, double* p) { - double result = x[0]; - double factor = 29.0; - double slope = 0.52; - if(TMath::Abs(x[0]) < 16.2) result=x[0]*slope; - else if(TMath::Abs(x[0]) < 49.8 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor; - else if(TMath::Abs(x[0]) < 85.2 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*2; - else result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*3; - return result; + const double sgn = (z_si >= z_pivot) ? +1.0 : -1.0; + int gnb = -1; double z_nb = (sgn > 0) ? +1.0e30 : -1.0e30; + for (int i = 0; i < 8; ++i) { + if (sgn > 0 && z_grid[i] > z_pivot + 1e-6 && z_grid[i] < z_nb) { z_nb = z_grid[i]; gnb = i; } + if (sgn < 0 && z_grid[i] < z_pivot - 1e-6 && z_grid[i] > z_nb) { z_nb = z_grid[i]; gnb = i; } + } + if (gnb < 0) return z_si; + const double W = TMath::Abs(z_nb - z_pivot); + if (W <= 0.0) return z_si; + int cell = (sgn > 0) ? gnb : gnb - 1; + if (cell < 0) cell = 0; if (cell > 6) cell = 6; + + double s2 = (c0_table[cell] - cfrac) / C; + if (s2 < 0.0 || s2 > 1.0) return z_si; // out-of-band -> Si (no boundary snap) + double s = TMath::Sqrt(s2); + double z_rec = z_pivot + sgn * s * W; + if (slope > 0.0 && slope != 1.0) { + const double zc = 0.5 * (z_grid[cell] + z_grid[cell + 1]); + z_rec = (z_rec - zc) / slope + zc; + } + if (TMath::Abs(z_rec - z_si) > W) return z_si; + ok = true; + return z_rec; } -// 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; -// } +// ---- Linear centre-fold variant ---- +inline double model_invert_a1c1_linear(double cfrac, double z_si, double z_pivot, + double k, const double* cfmin_table, bool& ok) +{ + ok = false; + static const double z_grid[8] = + {147.998, 101.946, 59.7634, 19.6965, -19.6965, -59.7634, -101.946, -147.998}; + const double cfrac_hi = 0.90, cfrac_lo = 0.12; + if (cfrac > cfrac_hi || cfrac < cfrac_lo) return z_si; + int cell = -1; + for (int i = 0; i < 7; ++i) + if (z_si <= z_grid[i] && z_si > z_grid[i + 1]) { cell = i; break; } + if (cell < 0) return z_si; + const double zc = 0.5 * (z_grid[cell] + z_grid[cell + 1]); + const double half = 0.5 * (z_grid[cell] - z_grid[cell + 1]); + if (half <= 0.0 || k <= 0.0) return z_si; -/*void testmodel() { - TF1 eqline("x","x",-200,200); - eqline.Draw(""); - eqline.SetLineStyle(kDashed); - - //TF1 f1("model",model,-200,200,2); - TF1 f1a("model_inv",model_a1c1,-200,200,2); - eqline.SetNpx(10000); - f1a.SetNpx(10000); - std::vector pars = {0.0,1.}; - f1a.SetParameters(pars.data()); - f1a.SetLineColor(kGreen+2); - f1a.SetLineStyle(kLine); - f1a.Draw("L SAME"); - - TF1 f1("model",model,-200,200,2); - //TF1 f1("model_inv",model_invert,-200,200,2); - eqline.SetNpx(10000); - f1.SetNpx(10000); - //std::vector pars = {0.0,1.}; - f1.SetParameters(pars.data()); - f1.SetLineColor(kGreen+2); - f1.SetLineStyle(kLine); - f1.Draw("L SAME"); - - gPad->Modified(); gPad->Update(); - while(gPad->WaitPrimitive()); -}*/ + double f = (cfrac - cfmin_table[cell]) / k; + if (f < 0.0) f = 0.0; + if (f > 1.0) return z_si; + const double d = f * half; + const double sgn = (z_pivot >= zc) ? +1.0 : -1.0; + const double z_rec = zc + sgn * d; + if (TMath::Abs(z_rec - z_si) > half) return z_si; + ok = true; + return z_rec; +} \ No newline at end of file diff --git a/TrackRecon.C b/TrackRecon.C index afe457c..efd7200 100644 --- a/TrackRecon.C +++ b/TrackRecon.C @@ -42,7 +42,7 @@ bool process_alpha_proton_scattering = false; bool doMiscHistograms = false; bool doPCSX3ClusterAnalysis = true; bool doPCQQQClusterAnalysis = true; -bool doOldAnalysis = false; +bool doOldAnalysis = true; bool do27AlapAnalysis = false; bool BenchMark = true; double source_vertex = 53; // 53 @@ -50,12 +50,34 @@ const double qqq_z = 105.0; double z_entrance = -174.3 - 9.7 - 100.0; const double anode_gain = 1.5146e-5; // channels --> MeV double dither_sigma = 8.0; -double dither_sigma_c0 = dither_sigma; +double dither_sigma_c0 = 16.0; +double Gain = 1; std::string dataset; int CO2percent; bool reactiondata = false; TF1 pcfix_func("func", model_invert, -200, 200); + +// --- A1C1 parabola calibration (cfrac = c0[cell] - C*s^2), gain-selected --- +const double a1c1_C_1 = 0.22; +const double a1c1_c0_1[7] = {0.46, 0.46, 0.46, 0.46, 0.45, 0.45, 0.45}; +const double a1c1_C_2 = 0.20; +const double a1c1_c0_2[7] = {0.24, 0.24, 0.24, 0.24, 0.24, 0.24, 0.24}; + +double a1c1_C = a1c1_C_1; +double a1c1_c0[7] = {0.46, 0.46, 0.46, 0.46, 0.45, 0.45, 0.45}; +double a1c1_slope = 1.0; // parabola staircase stretch (1.0 = off) + +// --- A1C1 linear centre-fold calibration (cfrac = cfmin[cell] + k*f) --- +const double a1c1_k_1 = 0.075; // from cfrac_vs_fold candle medians +const double a1c1_cfmin_1[7] = {0.40, 0.40, 0.40, 0.40, 0.40, 0.40, 0.40}; +const double a1c1_k_2= 0.06; // confirm with a 350-gain candle +const double a1c1_cfmin_2[7] = {0.18, 0.18, 0.18, 0.18, 0.18, 0.18, 0.18}; + +double a1c1_k = a1c1_k_1; +double a1c1_cfmin[7] = {0.40, 0.40, 0.40, 0.40, 0.40, 0.40, 0.40}; +bool a1c1_use_linear = true; // linear vs parabola toggle + TGraph *MeV_to_cm = NULL, *cm_to_MeV = NULL; TGraph *MeV_to_cm_p = NULL, *cm_to_MeVp = NULL; TGraph *MeV_to_cm_27Al = NULL, *cm_to_MeV_27Al = NULL; @@ -186,6 +208,30 @@ void TrackRecon::Begin(TTree * /*tree*/) std::cout << "Dither Sigma set to " << dither_sigma << " mm" << std::endl; } + if (getenv("Gain")) + Gain = std::atof(getenv("Gain")); + bool highP = (Gain >= 2); + + a1c1_C = highP ? a1c1_C_2 : a1c1_C_1; // parabola set + for (int i = 0; i < 7; ++i) + a1c1_c0[i] = highP ? a1c1_c0_2[i] : a1c1_c0_1[i]; + if (getenv("a1c1_C")) + a1c1_C = std::atof(getenv("a1c1_C")); + if (getenv("a1c1_slope")) + a1c1_slope = std::atof(getenv("a1c1_slope")); + + a1c1_k = highP ? a1c1_k_2 : a1c1_k_1; // linear set + for (int i = 0; i < 7; ++i) + a1c1_cfmin[i] = highP ? a1c1_cfmin_2[i] : a1c1_cfmin_1[i]; + if (getenv("a1c1_k")) + a1c1_k = std::atof(getenv("a1c1_k")); + if (getenv("a1c1_use_linear")) + a1c1_use_linear = std::atoi(getenv("a1c1_use_linear")); + + std::cout << "A1C1 cfrac model: gain=" << Gain << (a1c1_use_linear ? "LINEAR centre-fold (k=" : "parabola (C=") + << (a1c1_use_linear ? a1c1_k : a1c1_C) << ", cf/c0[3]=" + << (a1c1_use_linear ? a1c1_cfmin[3] : a1c1_c0[3]) << "), slope=" << a1c1_slope << std::endl; + pwinstance.ConstructGeo(); for (int i = 0; i < 48; i++) @@ -1556,16 +1602,14 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s double alpha_a1c1 = std::get<1>(xo_tuple); bool a1c1Good = (alpha_a1c1 != 9999999 && std::get<2>(xo_tuple) != -1); - // --- A1C1 charge-ratio diagnostic --- - // Only one cathode in A1C1, so the z-sensitive variable is cathode-vs-anode - // charge, not cathode/cathode. Use pseudo-wire sums for consistent gains. - double aSumE_bm = std::get<1>(pw_tuple); // anode sum (reuse pw_tuple) + // --- A1C1 charge fraction (single max-E cathode vs anode, pseudo-wire sums) --- + double aSumE_bm = std::get<1>(pw_tuple); auto cpw_tuple = pwinstance.GetPseudoWire(cOne, "CATHODE"); - double cSumE_bm = std::get<1>(cpw_tuple); // single-cathode sum + double cSumE_bm = std::get<1>(cpw_tuple); double ac_sum = aSumE_bm + cSumE_bm; - double cfrac = (ac_sum > 0.0) ? cSumE_bm / ac_sum : -1.0; // bounded [0,1] + double cfrac = (ac_sum > 0.0) ? cSumE_bm / ac_sum : -1.0; - // Smearing Setup + // Smearing setup double sx3_phi_pitch = 6.5 * (M_PI / 180.0); double smeared_phi = sx3event.pos.Phi() + rand.Uniform(-sx3_phi_pitch / 2.0, sx3_phi_pitch / 2.0); TVector3 smeared_sx3_pos(sx3event.pos.Perp() * TMath::Cos(smeared_phi), sx3event.pos.Perp() * TMath::Sin(smeared_phi), sx3event.pos.Z()); @@ -1574,10 +1618,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s { if (!a1c1Good) return; - // double sigma = hybrid ? (dither_sigma_c0 / 2.0) : dither_sigma_c0; - // double pcz = dither ? rand.Gaus(xo_a1c1.Z(), sigma) : xo_a1c1.Z(); double pcz = dither ? rand.Gaus(xo_a1c1.Z(), dither_sigma) : xo_a1c1.Z(); - // double pcz = dither ? rand.Uniform(xo_a1c1.Z(), dither_sigma) : xo_a1c1.Z(); TVector3 vtx = vertexFrom(si_point, TVector3(xo_a1c1.X(), xo_a1c1.Y(), pcz)); fillSuite(tag, pcz, vtx); fillVsRef(tag, pcz, vtx, pcz_ref, vtx_ref); @@ -1589,13 +1630,29 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s TVector3 vtx0 = vertexFrom(si_point, pc); if (!(vtx0.Perp() <= 6.0 && vtx0.Z() >= -173.6)) return; - // double sigma = hybrid ? (dither_sigma_c0 / 2.0) : dither_sigma_c0; double pcz = dither ? rand.Gaus(pc.Z(), dither_sigma) : pc.Z(); TVector3 vtx = vertexFrom(si_point, TVector3(pc.X(), pc.Y(), pcz)); fillSuite(tag, pcz, vtx); fillVsRef(tag, pcz, vtx, pcz_ref, vtx_ref); }; + // --- A1C1 with the cfrac sub-cell model (toggle linear/parabola) --- + // pivot = xo_a1c1.Z() (fired max-E cathode crossover). LINEAR: side from + // the fired cathode, |d| from cfrac about the cell centre. PARABOLA: side + // from the Si guess about the pivot. Both fall back to pczguess. + auto doA1C1Model = [&](const std::string &tag, const TVector3 &si_point) + { + if (!a1c1Good || cfrac < 0.0) + return; + bool ok = false; + double pcz = a1c1_use_linear + ? model_invert_a1c1_linear(cfrac, pczguess, xo_a1c1.Z(), a1c1_k, a1c1_cfmin, ok) + : model_invert_a1c1(cfrac, pczguess, xo_a1c1.Z(), a1c1_C, a1c1_c0, ok, a1c1_slope); + TVector3 vtx = vertexFrom(si_point, TVector3(xo_a1c1.X(), xo_a1c1.Y(), pcz)); + fillSuite(tag, pcz, vtx); + fillVsRef(tag, pcz, vtx, pcz_ref, vtx_ref); + }; + if (phicut && PCSX3TimeCut) { if (pcevent.multi1 == 1 && pcevent.multi2 == 2) @@ -1608,8 +1665,10 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s doA1C1("A1C1_Hyb", smeared_sx3_pos, true, true); doAnodeOnly("A1C0_Hyb", smeared_phi, smeared_sx3_pos, true, true); + // cfrac-model-corrected A1C1, directly comparable to A1C1 / A1C2 + doA1C1Model("A1C1_Cfrac", sx3event.pos); + // --- A1C1 charge-fraction diagnostics --- - // Decides whether the single cathode carries any sub-cell z information. if (a1c1Good && cfrac >= 0.0) { plotter->Fill1D("Benchmark_SX3_A1C1_cfrac", 220, -0.05, 1.05, cfrac, "Benchmark_SX3_ref"); @@ -1617,6 +1676,40 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s pcz_ref, cfrac, "Benchmark_SX3_ref"); plotter->Fill2D("Benchmark_SX3_A1C1_cfrac_vs_sx3pczguess", 400, -200, 200, 220, -0.05, 1.05, pczguess, cfrac, "Benchmark_SX3_ref"); + + static const double zg[8] = {147.998, 101.946, 59.7634, 19.6965, -19.6965, -59.7634, -101.946, -147.998}; + double zp = xo_a1c1.Z(); + auto fillCfracS = [&](const char *name, double truth) + { + double sgn = (truth >= zp) ? 1.0 : -1.0; + double znb = (sgn > 0) ? 1.0e30 : -1.0e30; + for (int i = 0; i < 8; ++i) + { + if (sgn > 0 && zg[i] > zp + 1e-6 && zg[i] < znb) + znb = zg[i]; + if (sgn < 0 && zg[i] < zp - 1e-6 && zg[i] > znb) + znb = zg[i]; + } + if (TMath::Abs(znb) < 1e8 && TMath::Abs(znb - zp) > 0.0) + plotter->Fill2D(name, 240, -1.2, 1.2, 220, -0.05, 1.05, + (truth - zp) / TMath::Abs(znb - zp), cfrac, "Benchmark_SX3_ref"); + }; + fillCfracS("Benchmark_SX3_A1C1_cfrac_vs_s", pcz_ref); + fillCfracS("Benchmark_SX3_A1C1_cfrac_vs_s_sx3pczguess", pczguess); + + // folded about the cell CENTRE: f = |ref - z_center|/halfcell in [0,1] + for (int i = 0; i < 7; ++i) + { + if (pcz_ref <= zg[i] && pcz_ref > zg[i + 1]) + { + double zc = 0.5 * (zg[i] + zg[i + 1]); + double half = 0.5 * (zg[i] - zg[i + 1]); + if (half > 0.0) + plotter->Fill2D("Benchmark_SX3_A1C1_cfrac_vs_fold", 120, 0, 1.2, 220, -0.05, 1.05, + TMath::Abs(pcz_ref - zc) / half, cfrac, "Benchmark_SX3_ref"); + break; + } + } } } } @@ -1647,7 +1740,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s TVector3 r_rhoMin = x1 + t_minimum * v; bool timecut = (qqqevent.Time1 - pcevent.Time1 < 150); - bool lowercut_cath = pcevent.Energy2 * sinTheta < 250 && (qqqevent.Energy2 < 5.0 || qqqevent.Energy1 < 5.0); + bool lowercut_cath = pcevent.Energy2 * sinTheta < 1 && (qqqevent.Energy2 < 5.0 || qqqevent.Energy1 < 5.0); bool phicut = qqqevent.pos.Phi() <= pcevent.pos.Phi() + TMath::Pi() / 4. && qqqevent.pos.Phi() >= pcevent.pos.Phi() - TMath::Pi() / 4.; if (lowercut_cath && phicut) @@ -1797,6 +1890,16 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s double alpha_a1c1 = std::get<1>(xo_tuple); bool a1c1Good = (alpha_a1c1 != 9999999 && std::get<2>(xo_tuple) != -1); + // --- A1C1 charge-fraction diagnostic --- + // One cathode in A1C1, so the z-sensitive variable is cathode-vs-anode + // charge (the V across each cell), not cathode/cathode. Use pseudo-wire + // sums for gain-consistent energies. + double aSumE_bm = std::get<1>(pw_tuple); // anode sum (reuse pw_tuple) + auto cpw_tuple = pwinstance.GetPseudoWire(cOne, "CATHODE"); + double cSumE_bm = std::get<1>(cpw_tuple); // single-cathode sum + double ac_sum = aSumE_bm + cSumE_bm; + double cfrac = (ac_sum > 0.0) ? cSumE_bm / ac_sum : -1.0; // bounded [0,1] + // Smearing Setup double qqq_wedge_pitch = (87.0 / 16.0) * (M_PI / 180.0); double qqq_ring_pitch = 48.0 / 16.0; @@ -1809,7 +1912,6 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s { if (!a1c1Good) return; - // double sigma = hybrid ? (dither_sigma / 2.0) : dither_sigma; double pcz = dither ? rand.Gaus(xo_a1c1.Z(), dither_sigma) : xo_a1c1.Z(); TVector3 vtx = vertexFrom(si_point, TVector3(xo_a1c1.X(), xo_a1c1.Y(), pcz)); fillSuite(tag, pcz, vtx); @@ -1829,6 +1931,25 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s fillVsRef(tag, pcz, vtx, pcz_ref, vtx_ref); }; + // --- A1C1 with the cfrac sub-cell model (QQQ) --- + // pivot = xo_a1c1.Z() (fired max-E cathode crossover). The active model + // is chosen by a1c1_use_linear: LINEAR centre-fold (side from the fired + // cathode, |d| from cfrac, scaled by the local cell) or the PARABOLA + // (side from the Si guess about the pivot). Both fall back to the Si + // guess on rails / out-of-band / inconsistency. + auto doA1C1Model = [&](const std::string &tag, const TVector3 &si_point) + { + if (!a1c1Good || cfrac < 0.0) + return; + bool ok = false; + double pcz = a1c1_use_linear + ? model_invert_a1c1_linear(cfrac, pcz_guess_int, xo_a1c1.Z(), a1c1_k, a1c1_cfmin, ok) + : model_invert_a1c1(cfrac, pcz_guess_int, xo_a1c1.Z(), a1c1_C, a1c1_c0, ok, a1c1_slope); + TVector3 vtx = vertexFrom(si_point, TVector3(xo_a1c1.X(), xo_a1c1.Y(), pcz)); + fillSuite(tag, pcz, vtx); + fillVsRef(tag, pcz, vtx, pcz_ref, vtx_ref); + }; + if (phicut && timecut) { if (pcevent.multi1 == 1 && pcevent.multi2 == 2) @@ -1840,6 +1961,54 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s doAnodeOnly("A1C0_Si", smeared_phi, smeared_qqq_pos, false, false); doA1C1("A1C1_Hyb", smeared_qqq_pos, true, true); doAnodeOnly("A1C0_Hyb", smeared_phi, smeared_qqq_pos, true, true); + + // --- Execute the cfrac-model for QQQ --- + doA1C1Model("A1C1_Cfrac", qqqevent.pos); + + // --- A1C1 charge-fraction diagnostics --- + if (a1c1Good && cfrac >= 0.0) + { + plotter->Fill1D("Benchmark_QQQ_A1C1_cfrac", 220, -0.05, 1.05, cfrac, "Benchmark_QQQ_ref"); + plotter->Fill2D("Benchmark_QQQ_A1C1_cfrac_vs_ref", 400, -200, 200, 220, -0.05, 1.05, + pcz_ref, cfrac, "Benchmark_QQQ_ref"); + plotter->Fill2D("Benchmark_QQQ_A1C1_cfrac_vs_qqqpczguess", 400, -200, 200, 220, -0.05, 1.05, + pcz_guess_37, cfrac, "Benchmark_QQQ_ref"); + + // --- cfrac vs normalized cell position (collapses all cells) --- + static const double zg[8] = {147.998, 101.946, 59.7634, 19.6965, -19.6965, -59.7634, -101.946, -147.998}; + double zp = xo_a1c1.Z(); + auto fillCfracS = [&](const char *name, double truth) + { + double sgn = (truth >= zp) ? 1.0 : -1.0; + double znb = (sgn > 0) ? 1.0e30 : -1.0e30; + for (int i = 0; i < 8; ++i) + { + if (sgn > 0 && zg[i] > zp + 1e-6 && zg[i] < znb) + znb = zg[i]; + if (sgn < 0 && zg[i] < zp - 1e-6 && zg[i] > znb) + znb = zg[i]; + } + if (TMath::Abs(znb) < 1e8 && TMath::Abs(znb - zp) > 0.0) + plotter->Fill2D(name, 240, -1.2, 1.2, 220, -0.05, 1.05, + (truth - zp) / TMath::Abs(znb - zp), cfrac, "Benchmark_QQQ_ref"); + }; + fillCfracS("Benchmark_QQQ_A1C1_cfrac_vs_s", pcz_ref); + fillCfracS("Benchmark_QQQ_A1C1_cfrac_vs_s_qqqpczguess", pcz_guess_37); + + // --- cfrac vs cell-centre fold: f = |ref - z_center|/halfcell in [0,1] --- + for (int i = 0; i < 7; ++i) + { + if (pcz_ref <= zg[i] && pcz_ref > zg[i + 1]) + { + double zc = 0.5 * (zg[i] + zg[i + 1]); + double half = 0.5 * (zg[i] - zg[i + 1]); + if (half > 0.0) + plotter->Fill2D("Benchmark_QQQ_A1C1_cfrac_vs_fold", 120, 0, 1.2, 220, -0.05, 1.05, + TMath::Abs(pcz_ref - zc) / half, cfrac, "Benchmark_QQQ_ref"); + break; + } + } + } } } } @@ -2235,6 +2404,24 @@ void TrackRecon::OldAnalysis() plotter->Fill2D("AnodeSumE_Vs_Cathode_Max_Energy_TC" + std::to_string(PCQQQTimeCut) + "_PC" + std::to_string(PCQQQPhiCut), 800, 0, 20000, 800, 0, 10000, aESum, cEMax, "hGMPC"); // plotter->Fill2D("AnodeSumE_Vs_Cathode_Max_Energy_path_corrected"+std::to_string(PCQQQTimeCut)+"_PC"+std::to_string(PCQQQPhiCut), 800, 0, 20000, 800, 0, 10000, aESum*sinTheta, cEMax*sinTheta, "hGMPC"); // plotter->Fill2D("AnodeSumE_Vs_Cathode_Max_Energy_path_corrected", 800, 0, 20000, 800, 0, 10000, aESum*sinTheta, cEMax*sinTheta, "hGMPC"); + if (aEMax > 0) + { + double ratio = cEMax / aEMax; + std::string folder = "Diagnostics_CMax"; + + // 1. Summary 2D Plots + plotter->Fill2D("CMax_over_Anode_vs_Z", 600, -300, 300, 200, 0, 2.0, anodeIntersection.Z(), ratio, folder); + plotter->Fill2D("CMax_over_Anode_vs_AnodeID", 24, 0, 24, 200, 0, 2.0, aIDMax, ratio, folder); + plotter->Fill2D("CMax_over_Anode_vs_CathodeID", 24, 0, 24, 200, 0, 2.0, cIDMax, ratio, folder); + + // 2. Individual 1D Histogram for this SPECIFIC Anode-Cathode Pair + std::string pairName = "Ratio_A" + std::to_string(aIDMax) + "_C" + std::to_string(cIDMax); + plotter->Fill1D(pairName, 200, 0, 2.0, ratio, folder + "/Pairs"); + + // (Optional) If you also still want the independent ones: + plotter->Fill1D("Ratio_A" + std::to_string(aIDMax), 200, 0, 2.0, ratio, folder + "/PerAnode"); + plotter->Fill1D("Ratio_C" + std::to_string(cIDMax), 200, 0, 2.0, ratio, folder + "/PerCathode"); + } if (PCQQQTimeCut && PCQQQPhiCut) { @@ -2274,11 +2461,11 @@ void miscHistograms_oneWire(HistPlotter *plotter, std::vector QQQ_Events, TRandom3 rand; rand.SetSeed(); // random seed setW double initial_energy = 7.0; - if (dataset == "27Al") /// m3 is alpha, 6.79 MeV is 7.0 MeV proton energy after kapton+100mm 4He gas (molar mass 5.6, 250 torr) + if (dataset == "27Al") /// m3 is alpha, 6.79 MeV is 7.0 MeV proton energy after kapton+100mm 4He gas (molar mass 5.6, 1 gain) initial_energy = 6.79; if (dataset == "17F") - initial_energy = 6.78; // m3 is alpha, 6.79 MeV is 7.0 MeV proton energy after kapton+100mm 4He gas (molar mass 5.6, 350 torr) - // initial_energy = 6.32; // m3 is alpha, 6.411 MeV is 7.0 MeV proton energy after havar+mylar+kapton+100mm 4He gas (molar mass 5.3, 250 torr) + initial_energy = 6.78; // m3 is alpha, 6.79 MeV is 7.0 MeV proton energy after kapton+100mm 4He gas (molar mass 5.6, 350 gain) + // initial_energy = 6.32; // m3 is alpha, 6.411 MeV is 7.0 MeV proton energy after havar+mylar+kapton+100mm 4He gas (molar mass 5.3, 1 gain) Kinematics apkin_a(1.008664916, 4.002603254, 4.002603254, 1.008664916, initial_energy); for (auto qqqevent : QQQ_Events) @@ -2368,9 +2555,9 @@ void protonMiscHistograms(HistPlotter *plotter, std::vector QQQ_Events, s rand.SetSeed(); // random seed set double initial_energy = 7.0; 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, 250 torr) + 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") - initial_energy = 6.32; // m3 is alpha, 6.411 MeV is 7.0 MeV proton energy after Havar+kapton+100mm 4He gas (molar mass 5.2, 250 torr) + initial_energy = 6.32; // m3 is alpha, 6.411 MeV is 7.0 MeV proton energy after Havar+kapton+100mm 4He gas (molar mass 5.2, 1 gain) Kinematics apkin_a(1.008664916, 4.002603254, 4.002603254, 1.008664916, initial_energy); // m3 is alpha @@ -2555,4 +2742,4 @@ void protonMiscHistograms_sx3(HistPlotter *plotter, std::vector QQQ_Event } // end PCEvents loop } // end sx3Events loop -} +} \ No newline at end of file diff --git a/run_17F.sh b/run_17F.sh index 7e62469..731bbe0 100644 --- a/run_17F.sh +++ b/run_17F.sh @@ -1,34 +1,51 @@ -rm results_run*.root +#!/bin/bash export DATASET="17F" export reactiondata=1 -#17F reaction data -#declare -i run=231 #49 -#while [[ $run -lt 258 ]]; do #392 -# wrun=$(printf "%03d" $run) -# file_exists=$(test -f ../ANASEN_analysis/data/17F_Data/Run_"$wrun"_mapped.root) -# if [[ $file_exists -ne 0 ]]; then continue; fi -# root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Run_"$wrun"_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run$wrun.root; -# run=run+1 -#done +run_once() { + local wrun=$(printf "%03d" "$1") + local prefix="${PREFIX:-Run_}" + local outdir="${OUT_DIR:-Output_default}" + local infile="../ANASEN_analysis/data/${DATASET}_Data/${prefix}${wrun}_mapped.root" + local out="${outdir}/results_run${wrun}.root" -function run_once() { - wrun=$(printf "%03d" $1) - file_exists=$(test -f ../ANASEN_analysis/data/17F_Data/Run_"$wrun"_mapped.root) - if [[ $file_exists -ne 0 ]]; then continue; fi - root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Run_"$wrun"_mapped.root -e 'tree->Process("TrackRecon.C+O","Analyzer_17F.root")'; mv Analyzer_17F.root Output_17F/results_run$wrun.root; + # Skip if the input file doesn't exist (return, not continue — this is a function) + if [ ! -f "$infile" ]; then + echo "SKIP: input $infile not found" + return + fi + + # Ensure the directory exists so ROOT doesn't fail silently + mkdir -p "$outdir" + + root -q -l -b -x "$infile" \ + -e "tree->Process(\"TrackRecon.C+\", \"${out}\")" > /dev/null 2>&1 + + if [ -f "$out" ]; then + echo "Run $wrun completed successfully in $outdir." + else + echo "ERROR: Run $wrun failed to generate $out" + fi } export -f run_once -# run_once 350 -# parallel -j 6 --ctag run_once {1} ::: {325..400} -parallel -j 6 --ctag run_once {1} ::: {351,353,355,358,359,360,362,367} -rm output_17F.root -hadd -j 4 -k output_17F.root Output_17F/results_run3*.root -unset souce_vertex +export DATASET="17F" +export PREFIX="Run_" +export OUT_DIR="Output_17F" +export source_vertex=-57.28 +export Gain=1 +rm -f ${OUT_DIR}/*.root + +# Pre-compile TrackRecon.C ONCE on a single core so parallel jobs don't race on ACLiC +echo "Pre-compiling TrackRecon.C..." +root -q -l -b -e '.L TrackRecon.C++O' + +# 3% CO2 +# parallel --bar -j 6 run_once ::: {325..400} +parallel --bar -j 6 run_once ::: 351 353 355 358 359 360 362 367 +hadd -j 4 -k ${OUT_DIR}/Output_17F.root ${OUT_DIR}/results_run*.root + +unset source_vertex unset DATASET -unset flip180 -unset flipa -unset anode_offset -unset reactiondata +unset reactiondata \ No newline at end of file diff --git a/run_27Al.sh b/run_27Al.sh index 5d7178c..85bef8c 100644 --- a/run_27Al.sh +++ b/run_27Al.sh @@ -1,22 +1,31 @@ #!/bin/bash export DATASET="27Al" +export PREFIX="Run_" +export OUT_DIR="Output_27Al" +export Gain=2 # Clean up previous runs -rm -f 27Al_output/results_run*.root output_27Al.root +rm -f ${OUT_DIR}/results_run*.root output_27Al.root echo "Pre-compiling TrackRecon.C safely on a single core..." root -q -l -b -e '.L TrackRecon.C++O' process_run() { - local wrun=$(printf "%03d" $1) - local out="Output_27Al/results_run${wrun}.root" + local wrun=$(printf "%03d" "$1") + local prefix="${PREFIX:-Run_}" + local outdir="${OUT_DIR:-Output_default}" + local infile="../ANASEN_analysis/data/${DATASET}_Data/${prefix}${wrun}_mapped.root" + local out="${outdir}/results_run${wrun}.root" - root -q -l -b -x "../ANASEN_analysis/data/27Al_Data/Run_${wrun}_mapped.root" \ + # Ensure the directory exists so ROOT doesn't fail silently + mkdir -p "$outdir" + + root -q -l -b -x "$infile" \ -e "tree->Process(\"TrackRecon.C+\", \"${out}\")" > /dev/null 2>&1 - + if [ -f "$out" ]; then - echo "Run $wrun completed successfully." + echo "Run $wrun completed successfully in $outdir." else echo "ERROR: Run $wrun failed to generate $out" fi @@ -28,6 +37,6 @@ echo "Starting parallel processing..." parallel --bar -j 4 process_run ::: {50..52} echo "Merging files..." -hadd -k -j 4 output_27Al.root Output_27Al/results_run*.root +hadd -k -j 4 ${OUT_DIR}/output_27Al.root ${OUT_DIR}/results_run*.root unset DATASET \ No newline at end of file diff --git a/run_tr.sh b/run_tr.sh index f71630e..71b4cab 100644 --- a/run_tr.sh +++ b/run_tr.sh @@ -46,25 +46,31 @@ if [[ 1 -eq 1 ]]; then export DATASET="27Al" export PREFIX="Run_" export OUT_DIR="Output_a" + export pressure_torr=350 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 pressuer_torr unset timecut_low fi # --- Block 3: 27Al Protons+Gas Runs (15, 17-22) --- -if [[ 1 -eq 0 ]]; then +if [[ 1 -eq 1 ]]; then export DATASET="27Al" export PREFIX="Run_" - export OUT_DIR="Output_p" + export OUT_DIR="Output_p" + export pressure_torr=350 + + rm -f ${OUT_DIR}/Al_protons.root export source_vertex=-200.0 # Source on the entrance window echo "Starting parallel processing for 27Al proton runs..." - process_run 17 - # 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 + unset pressuer_torr + fi # --- Block 4: 17F Source Runs (5-14) --- @@ -96,10 +102,11 @@ if [[ 1 -eq 0 ]]; then export DATASET="17F" export PREFIX="ProtonRun_" export OUT_DIR="Output_p" + rm -f ${OUT_DIR}/*pc*.root export source_vertex=-57.28 - # parallel --bar -j 6 process_run ::: {44..48} #3% CO2 - # hadd -j 4 -k ${OUT_DIR}/3pc.root ${OUT_DIR}/results_run0{44..48}.root + parallel --bar -j 6 process_run ::: {44..48} #3% CO2 + hadd -j 4 -k ${OUT_DIR}/3pc.root ${OUT_DIR}/results_run0{44..48}.root export CO2percent=4 parallel --bar -j 6 process_run ::: {38..42} #4% CO2