diff --git a/TrackRecon.C b/TrackRecon.C index 06b7d8b..7183b88 100644 --- a/TrackRecon.C +++ b/TrackRecon.C @@ -74,6 +74,46 @@ static const double a1c1_k_27Al[7] = {0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06}; double a1c1_cfmin_cell[7] = {0.40, 0.40, 0.40, 0.40, 0.40, 0.40, 0.40}; double a1c1_k_cell[7] = {0.075, 0.075, 0.075, 0.075, 0.075, 0.075, 0.075}; +// --- Dead / missing PC wires ------------------------------------------------- +// Some channels are unresponsive in a run (the white vertical gaps in +// PC_Index_Vs_Energy). A genuine two-wire track that straddles a dead wire +// collapses to a single fired wire: a "pseudo-1-wire" event. We still treat it +// as A1C1, but flag it so the excitation function can be split three ways: +// _all : every A1C1 event (unchanged from before) +// _true1w : neither neighbour of the fired anode/cathode is dead (genuine single) +// _missingw : a neighbouring wire is dead (suspected masked two-wire event) +// Index 0-23 within EACH plane; 1 = dead. Read the dead channels off the gaps +// in PC_Index_Vs_Energy (anode = index 0-23, cathode = index 24-47 -> w = idx-24) +// and fill the per-dataset arrays. Default all-alive => everything is _true1w +// and _missingw stays empty (no behaviour change until channels are entered). + +// static std::vector a1c1_dead_anode_17F = {}; // 1 can be recovered +// static std::vector a1c1_dead_cathode_17F = {}; // 0,13,15 gain-matched to 0 +// static std::vector a1c1_dead_anode_27Al = {}; +// static std::vector a1c1_dead_cathode_27Al = {}; + +static std::vector a1c1_dead_anode_17F = {9, 12}; // 1 can be recovered +static std::vector a1c1_dead_cathode_17F = {}; // 0,13,15 can be recovered +static std::vector a1c1_dead_anode_27Al = {0, 9, 12, 19}; +static std::vector a1c1_dead_cathode_27Al = {13}; +std::vector *a1c1_dead_anode = &a1c1_dead_anode_17F; // active set, chosen in Begin() +std::vector *a1c1_dead_cathode = &a1c1_dead_cathode_17F; + +// True if a neighbouring wire (index +/-1) of the fired anode OR cathode is +// dead -- i.e. this single-wire event may actually be a masked two-wire one. +inline bool a1c1_missing_neighbor(int awire, int cwire) +{ + auto deadAdj = [](const std::vector &dead, int w) + { + if (w < 0 || w >= 24) + return false; + auto isDead = [&](int x) + { return std::find(dead.begin(), dead.end(), x) != dead.end(); }; + return isDead(w - 1) || isDead(w + 1); + }; + return deadAdj(*a1c1_dead_anode, awire) || deadAdj(*a1c1_dead_cathode, cwire); +} + // --- A1C1 LOW BAND (incomplete charge integration) ----------------------------- // 17F shows a second, parallel cfrac band (~0.10 -> 0.15 over the fold) where the // cathode charge is only partially integrated. It still tracks z (real physics), @@ -87,7 +127,8 @@ static const double a1c1_k2_27Al[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; double a1c1_cfmin2_cell[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; double a1c1_k2_cell[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; -double a1c1_cfrac_split = 0.0; // cfrac below this uses the low band; <=0 disables +double a1c1_cfrac_split = 0.0; // cfrac below this uses the low band; <=0 disables +double a1c1_missing_fmax = 2.0; // f-ceiling when a neighbouring wire is dead (1.0 = no extension) // Sub-cell A1C1 z from cfrac (linear centre-fold). zf = crossover z (fired // cathode, on the grid), z_a1c0 = anode-only z (side reference). Anchors on the @@ -98,7 +139,8 @@ double a1c1_cfrac_split = 0.0; // cfrac below this uses the low band; <=0 disabl // pitchok : |pcz - zf| <= pitch (consistent with the fired wire) struct A1C1Sol { - double pcz; + double pcz; // raw inversion (can fall outside the cell if f<0 or f>1) + double pcz_clamped; // f clamped to [0,1] -> always inside the cell (no events lost) double f; int cell; double pitch; @@ -106,10 +148,14 @@ struct A1C1Sol bool pitchok; int band; }; -inline A1C1Sol a1c1_solve(double cfrac, double zf, double z_a1c0) +inline A1C1Sol a1c1_solve(double cfrac, double zf, double z_a1c0, int cwire = -1, int awire = -1) { - A1C1Sol s{zf, 0.0, 0, 0.0, false, false, 0}; - // band select: incomplete-integration low band uses its own per-cell cfmin/k + // Result structure. + // By default return the fired-wire position (zf) and mark the solution invalid. + A1C1Sol s{zf, zf, 0.0, 0, 0.0, false, false, 0}; + // Two independent calibrations are supported: + // band 0 : nominal integration + // band 1 : low-cfrac / incomplete-integration events const double *cfmin = a1c1_cfmin_cell; const double *kk = a1c1_k_cell; if (a1c1_cfrac_split > 0.0 && cfrac >= 0.0 && cfrac < a1c1_cfrac_split) @@ -118,25 +164,52 @@ inline A1C1Sol a1c1_solve(double cfrac, double zf, double z_a1c0) kk = a1c1_k2_cell; s.band = 1; } + + // Identify the fired cathode wire from the measured cathode position zf + // which is assumed to lie on one side of the calibrated cathode wire positions a1c1_zg[]. + int wf = 0; for (int i = 1; i < 8; ++i) if (TMath::Abs(a1c1_zg[i] - zf) < TMath::Abs(a1c1_zg[wf] - zf)) wf = i; + + // Determine the side of the fired cathode wire the event occurred onfrom teh anode only recon z_a1c0 + // Cells are indexed: wire0 --- cell0 --- wire1 --- cell1 --- ... --- wire7 giving seven physical drift cells. + int cell = (z_a1c0 >= a1c1_zg[wf]) ? (wf - 1) : wf; if (cell < 0) cell = 0; if (cell > 6) cell = 6; s.cell = cell; + + // Cell geometry. + // zc : cell centre half : half-cell width pitch: full wire-to-wire spacing + double zc = 0.5 * (a1c1_zg[cell] + a1c1_zg[cell + 1]); double half = 0.5 * (a1c1_zg[cell] - a1c1_zg[cell + 1]); s.pitch = a1c1_zg[cell] - a1c1_zg[cell + 1]; if (half <= 0.0 || kk[cell] <= 0.0) return s; + + // Convert cfrac into a normalized position coordinate f, with f = 0 -> cell centre and f = 1 -> fired wire + // Values outside [0,1] indicate an event outside the calibrated cfrac band. + s.f = (cfrac - cfmin[cell]) / kk[cell]; + // Determine whether the fired wire is above or below the cell centre. + // The sign maps increasing f toward the appropriate cathode wire. double sgn = (a1c1_zg[wf] >= zc) ? +1.0 : -1.0; s.pcz = zc + sgn * s.f * half; - s.inband = (s.f >= 0.0 && s.f <= 1.0); + // A dead neighbouring wire means the fired wire also collects the charge it + // would normally share, so cfrac (hence f) runs past the usual cfmin+k ceiling. + // Lift the upper limit toward the dead wire's cell so these events land there + // instead of being pinned to the fired-wire edge (f=1). + double fmax = a1c1_missing_neighbor(awire, cwire) ? a1c1_missing_fmax : 1.0; + double fc = (s.f < 0.0) ? 0.0 : (s.f > fmax ? fmax : s.f); + s.pcz_clamped = zc + sgn * fc * half; + s.inband = (s.f >= 0.0 && s.f <= fmax); + // Consistency check:reconstructed position should remain within one cell pitch of + // originally fired cathode wire. s.pitchok = (TMath::Abs(s.pcz - zf) <= s.pitch); return s; } @@ -283,6 +356,8 @@ void TrackRecon::Begin(TTree * /*tree*/) const double *cfmin2_src = a1c1_cfmin2_17F; const double *k2_src = a1c1_k2_17F; a1c1_cfrac_split = 0.25; // 17F: split off the incomplete-integration low band + a1c1_dead_anode = &a1c1_dead_anode_17F; + a1c1_dead_cathode = &a1c1_dead_cathode_17F; if (dataset == "27Al") { cfmin_src = a1c1_cfmin_27Al; @@ -290,6 +365,8 @@ void TrackRecon::Begin(TTree * /*tree*/) cfmin2_src = a1c1_cfmin2_27Al; k2_src = a1c1_k2_27Al; a1c1_cfrac_split = 0.0; // 27Al: no second band, low band disabled + a1c1_dead_anode = &a1c1_dead_anode_27Al; + a1c1_dead_cathode = &a1c1_dead_cathode_27Al; } for (int i = 0; i < 7; ++i) { @@ -1721,7 +1798,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s if (!a1c1Good || cfrac < 0.0) return; double z_a1c0 = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, si_point.Phi()).Z(); - A1C1Sol s = a1c1_solve(cfrac, xo_a1c1.Z(), z_a1c0); + A1C1Sol s = a1c1_solve(cfrac, xo_a1c1.Z(), z_a1c0, std::get<0>(cMaxWire)); if (!(s.inband && s.pitchok)) return; TVector3 vtx = vertexFrom(si_point, TVector3(xo_a1c1.X(), xo_a1c1.Y(), s.pcz)); @@ -1807,7 +1884,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s if (cfrac >= 0.0) { double z_a1c0 = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, sx3event.pos.Phi()).Z(); - A1C1Sol s = a1c1_solve(cfrac, xo_a1c1.Z(), z_a1c0); + A1C1Sol s = a1c1_solve(cfrac, xo_a1c1.Z(), z_a1c0, std::get<0>(cMaxWire)); int cell = s.cell; double f = s.f; double pcz_cf = s.pcz; @@ -2102,7 +2179,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s if (!a1c1Good || cfrac < 0.0) return; double z_a1c0 = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, si_point.Phi()).Z(); - A1C1Sol s = a1c1_solve(cfrac, xo_a1c1.Z(), z_a1c0); + A1C1Sol s = a1c1_solve(cfrac, xo_a1c1.Z(), z_a1c0, std::get<0>(cMaxWire)); if (!(s.inband && s.pitchok)) return; TVector3 vtx = vertexFrom(si_point, TVector3(xo_a1c1.X(), xo_a1c1.Y(), s.pcz)); @@ -2190,7 +2267,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s if (cfrac >= 0.0) { double z_a1c0 = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, qqqevent.pos.Phi()).Z(); - A1C1Sol s = a1c1_solve(cfrac, xo_a1c1.Z(), z_a1c0); + A1C1Sol s = a1c1_solve(cfrac, xo_a1c1.Z(), z_a1c0, std::get<0>(cMaxWire)); int cell = s.cell; double f = s.f; double pcz_cf = s.pcz; @@ -2348,7 +2425,7 @@ void TrackRecon::OldAnalysis() cESum += cE; } }*/ - if ((aIDMax + cID) % 24 == 22 || (aIDMax + cID) % 24 == 23 || (aIDMax + cID) % 24 >= 0 || (aIDMax + cID) % 24 <= 3) + if (((aIDMax + cID) % 24) >= 20 || ((aIDMax + cID) % 24) <= 3) { corrcatMax.push_back(std::pair(cID, cE)); cESum += cE; @@ -2568,10 +2645,10 @@ void TrackRecon::OldAnalysis() else continue; - if (qqqCalibValid[qqq.id[i]][chRing][chWedge]) + if (qqqCalibValid[qqq.id[i]][chWedge][chRing]) { - eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000; - eRingMeV = eRing * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000; + eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chWedge][chRing] / 1000; + eRingMeV = eRing * qqqCalib[qqq.id[i]][chWedge][chRing] / 1000; } else continue; @@ -2835,6 +2912,9 @@ void protonMiscHistograms(HistPlotter *plotter, std::vector QQQ_Events, s // cathode, side from the anode-only z, offset from the per-cell autocal. if (pcevent.multi2 == 1 && pcevent.Energy2 > 1400) { + // wire-topology category: _true1w (no dead neighbour) vs _missingw + // (a neighbouring wire is dead -> possible masked two-wire event). + const std::string wcat = a1c1_missing_neighbor(pcevent.Anodech, pcevent.Cathodech) ? "_missingw" : "_true1w"; auto fillCmp = [&](double pcz, const std::string &m) { TVector3 x2(pcevent.pos.X(), pcevent.pos.Y(), pcz); @@ -2846,12 +2926,17 @@ void protonMiscHistograms(HistPlotter *plotter, std::vector QQQ_Events, s double th = (qqqevent.pos - rv).Theta(); double pl = (qqqevent.pos - rv).Mag() * 0.1; double Ef = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy1) - pl); + double Ex = apkin_a.getExc(Ef, th * 180 / M_PI); std::string lbl = "proton+misc_a1c1cmp"; - plotter->Fill1D("pmisc_a1c1cmp_pcz_" + m, 600, -300, 300, pcz, lbl); - plotter->Fill1D("pmisc_a1c1cmp_Ex_" + m, 200, -10, 10, apkin_a.getExc(Ef, th * 180 / M_PI), lbl); - plotter->Fill1D("pmisc_a1c1cmp_VertexZ_" + m, 800, -400, 400, rv.Z(), lbl); - plotter->Fill2D("pmisc_a1c1cmp_VertexZ_vs_Ef_" + m, 800, -400, 400, 800, 0, 20, rv.Z(), Ef, lbl); - plotter->Fill2D("pmisc_a1c1cmp_Ef_vs_theta_" + m, 100, 0, 180, 800, 0, 20, th * 180 / M_PI, Ef, lbl); + // fill "all" (existing names) plus the wire-topology split (_true1w/_missingw) + for (const std::string &w : {std::string(""), wcat}) + { + plotter->Fill1D("pmisc_a1c1cmp_pcz_" + m + w, 600, -300, 300, pcz, lbl); + plotter->Fill1D("pmisc_a1c1cmp_Ex_" + m + w, 200, -10, 10, Ex, lbl); + plotter->Fill1D("pmisc_a1c1cmp_VertexZ_" + m + w, 800, -400, 400, rv.Z(), lbl); + plotter->Fill2D("pmisc_a1c1cmp_VertexZ_vs_Ef_" + m + w, 800, -400, 400, 800, 0, 20, rv.Z(), Ef, lbl); + plotter->Fill2D("pmisc_a1c1cmp_Ef_vs_theta_" + m + w, 100, 0, 180, 800, 0, 20, th * 180 / M_PI, Ef, lbl); + } }; fillCmp(pcz_dith, "dither"); // method 1: Gaussian dither (main-flow value) @@ -2865,7 +2950,10 @@ void protonMiscHistograms(HistPlotter *plotter, std::vector QQQ_Events, s std::vector> aOne = {std::make_tuple(pcevent.Anodech, 1.0, 0.0)}; auto apw = pwinstance.GetPseudoWire(aOne, "ANODE"); double z_a1c0 = pwinstance.getClosestWirePosAtWirePhi(std::get<0>(apw), qqqevent.pos.Phi()).Z(); - A1C1Sol s = a1c1_solve(cfrac, pcevent.pos.Z(), z_a1c0); + A1C1Sol s = a1c1_solve(cfrac, pcevent.pos.Z(), z_a1c0, pcevent.Cathodech); + // ungated: every alpha A1C1 event, f clamped to the cell (same stats as + // dither) -- shows the cfrac method without the acceptance cut. + fillCmp(s.pcz_clamped, "cfrac_all"); if (s.inband && s.pitchok) { fillCmp(s.pcz, "cfrac"); @@ -3051,6 +3139,9 @@ void protonMiscHistograms_sx3(HistPlotter *plotter, std::vector QQQ_Event if (!(pcevent.Energy2 > 1400)) continue; + // wire-topology category: _true1w (no dead neighbour) vs _missingw + // (a neighbouring wire is dead -> possible masked two-wire event). + const std::string wcat = a1c1_missing_neighbor(pcevent.Anodech, pcevent.Cathodech) ? "_missingw" : "_true1w"; auto fillCmp = [&](double pcz, const std::string &m) { TVector3 x2(pcevent.pos.X(), pcevent.pos.Y(), pcz); @@ -3063,10 +3154,14 @@ void protonMiscHistograms_sx3(HistPlotter *plotter, std::vector QQQ_Event double pl = (sx3event.pos - rv).Mag() * 0.1; double Ef = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(sx3event.Energy1) - pl); std::string lbl = "proton+miscsx3_a1c1cmp"; - plotter->Fill1D("pmiscs_a1c1cmp_pcz_" + m, 600, -300, 300, pcz, lbl); - plotter->Fill1D("pmiscs_a1c1cmp_VertexZ_" + m, 800, -400, 400, rv.Z(), lbl); - plotter->Fill2D("pmiscs_a1c1cmp_VertexZ_vs_Ef_" + m, 800, -400, 400, 800, 0, 20, rv.Z(), Ef, lbl); - plotter->Fill2D("pmiscs_a1c1cmp_Ef_vs_theta_" + m, 100, 0, 180, 800, 0, 20, th * 180 / M_PI, Ef, lbl); + // fill "all" (existing names) plus the wire-topology split (_true1w/_missingw) + for (const std::string &w : {std::string(""), wcat}) + { + plotter->Fill1D("pmiscs_a1c1cmp_pcz_" + m + w, 600, -300, 300, pcz, lbl); + plotter->Fill1D("pmiscs_a1c1cmp_VertexZ_" + m + w, 800, -400, 400, rv.Z(), lbl); + plotter->Fill2D("pmiscs_a1c1cmp_VertexZ_vs_Ef_" + m + w, 800, -400, 400, 800, 0, 20, rv.Z(), Ef, lbl); + plotter->Fill2D("pmiscs_a1c1cmp_Ef_vs_theta_" + m + w, 100, 0, 180, 800, 0, 20, th * 180 / M_PI, Ef, lbl); + } }; fillCmp(rand.Gaus(pcevent.pos.Z(), 8.0), "dither"); // method 1: Gaussian dither baseline @@ -3080,7 +3175,10 @@ void protonMiscHistograms_sx3(HistPlotter *plotter, std::vector QQQ_Event std::vector> aOne = {std::make_tuple(pcevent.Anodech, 1.0, 0.0)}; auto apw = pwinstance.GetPseudoWire(aOne, "ANODE"); double z_a1c0 = pwinstance.getClosestWirePosAtWirePhi(std::get<0>(apw), sx3event.pos.Phi()).Z(); - A1C1Sol s = a1c1_solve(cfrac, pcevent.pos.Z(), z_a1c0); + A1C1Sol s = a1c1_solve(cfrac, pcevent.pos.Z(), z_a1c0, pcevent.Cathodech); + // ungated: every alpha A1C1 event, f clamped to the cell (same stats as + // dither) -- shows the cfrac method without the acceptance cut. + fillCmp(s.pcz_clamped, "cfrac_all"); if (s.inband && s.pitchok) fillCmp(s.pcz, "cfrac"); } diff --git a/fit_a1c1_cfrac.C b/fit_a1c1_cfrac.C index 86a2ce2..eb3e4ec 100644 --- a/fit_a1c1_cfrac.C +++ b/fit_a1c1_cfrac.C @@ -1,3 +1,36 @@ +// fit_a1c1_cfrac.C +// ============================================================================= +// Offline per-cell "gain match" for the A1C1 linear centre-fold model. +// +// The raw cfrac = cpMax/(apSum+cpMax) pedestal floats with the arbitrary +// anode/cathode gain, so the model cfrac = cfmin[cell] + k[cell]*fold has to be +// re-fit per dataset. This macro does that fit OFFLINE from an already-processed +// results file (no reprocessing of raw data), replacing the in-line autocal that +// used to live in TrackRecon.C. +// +// It reads the reference histograms +// Benchmark_QQQ_A1C1_cfrac_vs_ref (X = true z = pcz_ref, Y = cfrac) +// Benchmark_SX3_A1C1_cfrac_vs_ref +// filled from A1C2 events (where the true z is known from the crossover), folds +// each cell about its centre (fold = |z - z_centre| / halfcell, in [0,1]) and +// does a weighted least-squares fit of cfrac vs fold per cathode cell. QQQ and +// SX3 are combined (same proportional-counter cells). +// +// Output: a per-cell table, paste-ready C++ array literals for TrackRecon.C +// (a1c1_cfmin_[7] / a1c1_k_[7]), and optionally a .dat file. +// +// Usage (compile once with + then call): +// root -l -b -q 'fit_a1c1_cfrac.C+("Output_17F/Output_17F.root")' // build +// MAIN band: fit_a1c1_cfrac("Output_17F/Output_17F.root","17F",0.25,1.05) +// LOW band: fit_a1c1_cfrac("Output_17F/Output_17F.root","17F",0.00,0.25,50,false,"2") +// 27Al: fit_a1c1_cfrac("Output_27Al/output_27Al.root","27Al",0.10,1.05) +// Args: cfrac_min, cfrac_max (band window), min_counts (per-cell fit threshold), +// write_dat (dump a1c1_cfrac_.dat), band ("" main, "2" low band -> +// prints a1c1_cfmin2_), all_cells (pool all 7 cells into one fit +// and broadcast -- robust when per-cell stats are thin / pedestals uniform). +// ALL-CELLS: fit_a1c1_cfrac("Output_cal/17F.root","17F",0.25,1.05,50,false,"",true) +// ============================================================================= + #include #include #include @@ -9,6 +42,7 @@ #include #include #include +#include namespace { // cathode wire grid; cell i spans zg[i] (high) .. zg[i+1] (low) @@ -21,7 +55,7 @@ int cellOf(double z) { return -1; } -// recursively search a directory tree for a TH2 by name +// recursively search a directory tree for a TH2 by name (folders vary by group) TH2 *findTH2(TDirectory *dir, const std::string &name) { if (TH2 *h = dynamic_cast(dir->Get(name.c_str()))) return h; TIter next(dir->GetListOfKeys()); @@ -36,151 +70,227 @@ TH2 *findTH2(TDirectory *dir, const std::string &name) { } } // namespace -// ============================================================================= -// GLOBAL WITH REFERENCE (Least Squares Fit) -// ============================================================================= +// cfrac_min/cfrac_max bracket the band to fit: [0.25, 1.05] for the MAIN band, +// [0.0, 0.25] for the 17F low (incomplete-integration) band. band="" prints the +// a1c1_cfmin_ arrays; band="2" prints a1c1_cfmin2_ (the low band). void fit_a1c1_cfrac(const char *filename, const char *dataset = "DATA", - double cfrac_min = 0.0, double cfrac_max = 1.05, - long min_counts = 50, bool write_dat = false, - const char *band = "") { + double cfrac_min = 0.0, double cfrac_max = 1.05, + long min_counts = 50, bool write_dat = false, + const char *band = "", bool all_cells = false) { TFile *f = TFile::Open(filename, "READ"); if (!f || f->IsZombie()) { - std::cerr << "Cannot open " << filename << std::endl; + std::cerr << "fit_a1c1_cfrac: cannot open " << filename << std::endl; return; } const char *names[2] = {"Benchmark_QQQ_A1C1_cfrac_vs_ref", "Benchmark_SX3_A1C1_cfrac_vs_ref"}; - // Global weighted OLS accumulators across ALL cells - double n_tot = 0, sf_tot = 0, sff_tot = 0, sc_tot = 0, sfc_tot = 0; + // weighted OLS accumulators per cell: cfrac = cfmin + k*fold + double n[7] = {0}, sf[7] = {0}, sff[7] = {0}, sc[7] = {0}, sfc[7] = {0}; for (int s = 0; s < 2; ++s) { TH2 *h = findTH2(f, names[s]); - if (!h) continue; - + if (!h) { + std::cout << " (note: " << names[s] << " not found -- skipping)" << std::endl; + continue; + } + std::cout << "Using " << names[s] << " (" << h->GetEntries() << " entries)" << std::endl; const int nx = h->GetNbinsX(), ny = h->GetNbinsY(); for (int ix = 1; ix <= nx; ++ix) { const double z = h->GetXaxis()->GetBinCenter(ix); const int c = cellOf(z); if (c < 0) continue; - const double zc = 0.5 * (zg[c] + zg[c + 1]); const double half = 0.5 * (zg[c] - zg[c + 1]); if (half <= 0.0) continue; - const double fold = TMath::Abs(z - zc) / half; // [0,1] - for (int iy = 1; iy <= ny; ++iy) { const double w = h->GetBinContent(ix, iy); if (w <= 0.0) continue; const double cfrac = h->GetYaxis()->GetBinCenter(iy); if (cfrac < cfrac_min || cfrac > cfrac_max) continue; - - // Accumulate globally - n_tot += w; - sf_tot += w * fold; - sff_tot += w * fold * fold; - sc_tot += w * cfrac; - sfc_tot += w * fold * cfrac; + n[c] += w; + sf[c] += w * fold; + sff[c] += w * fold * fold; + sc[c] += w * cfrac; + sfc[c] += w * fold * cfrac; } } } - // Global least-squares calculation - double cfmin_global = 0.40, k_global = 0.075; - const char *status = "FALLBACK"; + // per-cell weighted least-squares; keep a fallback for starved cells + const double cfmin_fallback = 0.40, k_fallback = 0.075; + double cfmin[7], k[7]; - if (n_tot > min_counts) { - const double denom = n_tot * sff_tot - sf_tot * sf_tot; - if (TMath::Abs(denom) > 1e-9) { - const double kfit = (n_tot * sfc_tot - sf_tot * sc_tot) / denom; - const double cffit = (sc_tot - kfit * sf_tot) / n_tot; - if (kfit > 0.01) { // reject degenerate/flat fit - cfmin_global = cffit; - k_global = kfit; - status = "FIT"; + if (all_cells) { + // pool every cell into one fit. fold is already cell-relative (|z-zc|/half), + // so pooling the (fold, cfrac) pairs is valid -- robust when per-cell stats + // are thin and the pedestals are ~uniform. Same value written to all 7 cells. + double N = 0, SF = 0, SFF = 0, SC = 0, SFC = 0; + for (int c = 0; c < 7; ++c) { N += n[c]; SF += sf[c]; SFF += sff[c]; SC += sc[c]; SFC += sfc[c]; } + double cf = cfmin_fallback, kk = k_fallback; + const char *status = "FALLBACK (starved)"; + if (N > min_counts) { + const double denom = N * SFF - SF * SF; + if (TMath::Abs(denom) > 1e-9) { + const double kfit = (N * SFC - SF * SC) / denom; + const double cffit = (SC - kfit * SF) / N; + if (kfit > 0.01) { cf = cffit; kk = kfit; status = "fit"; } + else status = "FALLBACK (flat k)"; } } + for (int c = 0; c < 7; ++c) { cfmin[c] = cf; k[c] = kk; } + printf("\n ALL-CELLS pooled fit: cfmin=%.6f k=%.6f n(eff)=%.0f %s\n", cf, kk, N, status); + } else { + std::cout << "\n cell cfmin k n(eff) status" << std::endl; + std::cout << " ---- --------- --------- -------- ------" << std::endl; + for (int c = 0; c < 7; ++c) { + double cf = cfmin_fallback, kk = k_fallback; + const char *status = "FALLBACK (starved)"; + if (n[c] > min_counts) { + const double denom = n[c] * sff[c] - sf[c] * sf[c]; + if (TMath::Abs(denom) > 1e-9) { + const double kfit = (n[c] * sfc[c] - sf[c] * sc[c]) / denom; + const double cffit = (sc[c] - kfit * sf[c]) / n[c]; + if (kfit > 0.01) { // reject degenerate/flat fit + cf = cffit; + kk = kfit; + status = "fit"; + } else { + status = "FALLBACK (flat k)"; + } + } + } + cfmin[c] = cf; + k[c] = kk; + printf(" %d %9.6f %9.6f %8.0f %s\n", c, cf, kk, n[c], status); + } } - printf("\n=== GLOBAL FIT RESULT ===\n"); - printf(" cfmin: %9.6f | k: %9.6f | N: %8.0f | Status: %s\n", - cfmin_global, k_global, n_tot, status); - - // Print 7-element array with the identical global value for backward compatibility - printf("\n// ---- paste into TrackRecon.C (Global Fit) ----\n"); + // paste-ready literals for TrackRecon.C (band="" -> main set, "2" -> low band) + printf("\n// ---- paste into TrackRecon.C (cfrac in [%.2f, %.2f]) ----\n", cfrac_min, cfrac_max); printf("static const double a1c1_cfmin%s_%s[7] = {", band, dataset); - for (int c = 0; c < 7; ++c) printf("%.6f%s", cfmin_global, c < 6 ? ", " : "};\n"); - + for (int c = 0; c < 7; ++c) printf("%.6f%s", cfmin[c], c < 6 ? ", " : "};\n"); printf("static const double a1c1_k%s_%s[7] = {", band, dataset); - for (int c = 0; c < 7; ++c) printf("%.6f%s", k_global, c < 6 ? ", " : "};\n"); + for (int c = 0; c < 7; ++c) printf("%.6f%s", k[c], c < 6 ? ", " : "};\n"); if (write_dat) { - const std::string fn = std::string("a1c1_cfrac_global_") + dataset + ".dat"; + const std::string fn = std::string("a1c1_cfrac_") + dataset + ".dat"; std::ofstream out(fn); - out << "# GLOBAL cfmin k (offline fit cfrac = cfmin + k*fold; dataset=" << dataset << ")\n"; - for (int c = 0; c < 7; ++c) out << c << " " << cfmin_global << " " << k_global << "\n"; + out << "# cell cfmin k (offline fit cfrac = cfmin + k*fold; dataset=" << dataset << ")\n"; + for (int c = 0; c < 7; ++c) out << c << " " << cfmin[c] << " " << k[c] << "\n"; out.close(); + std::cout << "\nwrote " << fn << std::endl; } - + f->Close(); } // ============================================================================= -// GLOBAL REFERENCE-FREE (Percentile Fit) +// REFERENCE-FREE gain match. +// +// Does NOT use the A1C2 true-z reference. Reads the per-cell A1C1 cfrac +// histogram Benchmark_{QQQ,SX3}_trueA1C1_cfrac_vs_cell (cell assignment is pure +// geometry) and recovers the constants from the EDGES of each cell's cfrac +// distribution: within a cell cfrac = cfmin + k*fold with fold in [0,1], so the +// low percentile ~ cfmin (fold->0, cell centre) and the high percentile ~ +// cfmin+k (fold->1, at the wire). Robust percentiles (plo/phi) are used instead +// of min/max. +// +// This aligns the per-cell pedestals (the actual gain match) WITHOUT a reference, +// but it cannot fix the absolute z scale -- it assumes the linear model and that +// each cell is illuminated across its length. If shared_k=true, a single +// geometric k (median of the per-cell spans) is used for every cell and only the +// per-cell cfmin pedestal is taken from the data -- the most robust mode. +// +// Usage: +// root -l -b -q 'fit_a1c1_cfrac.C+("Output_17F/Output_17F.root")' // compile +// root -l -b -q 'fit_a1c1_cfrac_reffree("Output_17F/Output_17F.root","17F")' // ============================================================================= -void fit_a1c1_cfrac_reffree_global(const char *filename, const char *dataset = "DATA", - double plo = 0.05, double phi = 0.95, - double min_counts = 50, bool write_dat = false) { +void fit_a1c1_cfrac_reffree(const char *filename, const char *dataset = "DATA", + double plo = 0.05, double phi = 0.95, + bool shared_k = false, double min_counts = 50, + bool write_dat = false) { TFile *f = TFile::Open(filename, "READ"); - if (!f || f->IsZombie()) return; + if (!f || f->IsZombie()) { + std::cerr << "fit_a1c1_cfrac_reffree: cannot open " << filename << std::endl; + return; + } const char *names[2] = {"Benchmark_QQQ_trueA1C1_cfrac_vs_cell", "Benchmark_SX3_trueA1C1_cfrac_vs_cell"}; + // sum the QQQ + SX3 per-cell cfrac histograms into one (cell x cfrac) TH2 *acc = nullptr; for (int s = 0; s < 2; ++s) { TH2 *h = findTH2(f, names[s]); - if (!h) continue; + if (!h) { + std::cout << " (note: " << names[s] << " not found -- skipping)" << std::endl; + continue; + } + std::cout << "Using " << names[s] << " (" << h->GetEntries() << " entries)" << std::endl; if (!acc) { acc = (TH2 *)h->Clone("acc_cfrac_vs_cell"); acc->SetDirectory(nullptr); } else acc->Add(h); } - if (!acc) { std::cerr << "no cfrac_vs_cell histograms found\n"; f->Close(); return; } + if (!acc) { std::cerr << "no cfrac_vs_cell histograms found" << std::endl; f->Close(); return; } - // Project the ENTIRE detector onto Y to get a global distribution - TH1D *proj = acc->ProjectionY("proj_global"); - const double n_tot = proj->Integral(); + const double cfmin_fallback = 0.40, k_fallback = 0.075; + double cfmin[7], k[7], span[7]; + bool ok[7]; - double cfmin_global = 0.40, k_global = 0.075; - const char *status = "FALLBACK"; - - if (n_tot > min_counts) { - double q[2], p[2] = {plo, phi}; - proj->GetQuantiles(2, q, p); - double lo = q[0]; - double hi = q[1]; - - if (hi - lo > 0.01) { - cfmin_global = lo; - k_global = hi - lo; - status = "FIT"; + std::cout << "\n cell cfmin cfmin+k k n status" << std::endl; + std::cout << " ---- --------- --------- --------- ------ ------" << std::endl; + for (int c = 0; c < 7; ++c) { + TH1D *proj = acc->ProjectionY(Form("proj_c%d", c), c + 1, c + 1); + const double n = proj->Integral(); + double lo = cfmin_fallback, hi = cfmin_fallback + k_fallback; + ok[c] = false; + if (n > min_counts) { + double q[2], p[2] = {plo, phi}; + proj->GetQuantiles(2, q, p); + lo = q[0]; + hi = q[1]; + ok[c] = (hi - lo > 0.01); } + cfmin[c] = ok[c] ? lo : cfmin_fallback; + span[c] = ok[c] ? (hi - lo) : k_fallback; + k[c] = span[c]; + printf(" %d %9.6f %9.6f %9.6f %6.0f %s\n", c, cfmin[c], + cfmin[c] + span[c], span[c], n, ok[c] ? "fit" : "FALLBACK"); + delete proj; } - printf("\n=== GLOBAL REFERENCE-FREE FIT ===\n"); - printf(" cfmin: %9.6f | k: %9.6f | N: %8.0f | Status: %s\n", - cfmin_global, k_global, n_tot, status); + if (shared_k) { + // median of the well-determined per-cell spans -> one geometric k for all cells + double v[7]; + int m = 0; + for (int c = 0; c < 7; ++c) if (ok[c]) v[m++] = span[c]; + double kg = k_fallback; + if (m > 0) { + std::sort(v, v + m); + kg = (m % 2) ? v[m / 2] : 0.5 * (v[m / 2 - 1] + v[m / 2]); + } + for (int c = 0; c < 7; ++c) k[c] = kg; + std::cout << "\nshared geometric k = " << kg << " (median of " << m << " fitted cells)" << std::endl; + } - // Print arrays - printf("\n// ---- paste into TrackRecon.C (Global Ref-Free) ----\n"); + printf("\n// ---- paste into TrackRecon.C (reference-free gain match) ----\n"); printf("static const double a1c1_cfmin_%s[7] = {", dataset); - for (int c = 0; c < 7; ++c) printf("%.6f%s", cfmin_global, c < 6 ? ", " : "};\n"); - + for (int c = 0; c < 7; ++c) printf("%.6f%s", cfmin[c], c < 6 ? ", " : "};\n"); printf("static const double a1c1_k_%s[7] = {", dataset); - for (int c = 0; c < 7; ++c) printf("%.6f%s", k_global, c < 6 ? ", " : "};\n"); + for (int c = 0; c < 7; ++c) printf("%.6f%s", k[c], c < 6 ? ", " : "};\n"); + + if (write_dat) { + const std::string fn = std::string("a1c1_cfrac_") + dataset + ".dat"; + std::ofstream out(fn); + out << "# cell cfmin k (reference-free gain match, percentiles " << plo << "/" << phi + << "; dataset=" << dataset << ")\n"; + for (int c = 0; c < 7; ++c) out << c << " " << cfmin[c] << " " << k[c] << "\n"; + out.close(); + std::cout << "\nwrote " << fn << std::endl; + } - delete proj; delete acc; f->Close(); } \ No newline at end of file diff --git a/slope_intercept_results_17F.dat b/slope_intercept_results_17F.dat index 92d002f..db5b54b 100644 --- a/slope_intercept_results_17F.dat +++ b/slope_intercept_results_17F.dat @@ -1,6 +1,6 @@ #Histogram Number Slope Intercept 0 0.937314 -16.871 -1 0 0 +1 1 -1.87356e-10 2 0.965461 -1.54376 3 0.926501 -3.27662 4 0.905634 2.54577 @@ -8,10 +8,10 @@ 6 0.853919 6.23079 7 0.945588 -9.54044 8 0.884454 -11.8262 -9 0.922501 -3.42538 +9 0 0 10 0.903053 9.28069 11 0.914653 9.87642 -12 0.965332 13.2526 +12 0 0 13 0.923847 -3.41775 14 0.93845 25.9901 15 0.955424 12.324 @@ -23,7 +23,7 @@ 21 0.969834 -45.001 22 0.89304 -31.5635 23 0.933226 4.02193 -24 0 0 +24 1 -2.89219e-10 25 0.941896 6.16135 26 0.980284 2.86886 27 0.983166 -3.82952 @@ -36,14 +36,14 @@ 34 0.928892 7.61384 35 0.947376 -0.644223 36 0.875342 6.066 -37 0 0 +37 0.918408 -3.27891 38 0.970953 6.262 -39 0 0 +39 0.918408 -3.27891 40 0.918408 -3.27891 41 0.913619 4.11288 42 0.954083 2.21261 43 0.993037 5.48924 -44 0 0 +44 0.954083 2.21261 45 0.926406 -19.719 46 1.00459 5.14574 47 0.942483 5.54183