diff --git a/CLAUDE.md b/CLAUDE.md deleted file mode 100644 index 9b1aa4e..0000000 --- a/CLAUDE.md +++ /dev/null @@ -1,22 +0,0 @@ -# CLAUDE.md - -This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository. - -## Git policy - -- **Never push.** Only commit when explicitly asked. Remember only the last commit made. - -## Build - -ROOT macros (`.C`) compile via ACLiC inside ROOT — not linked into binaries. Standalone binaries (`EventBuilder`, `Mapper`) are built with `make` in `Armory/`. - -```bash -cd Armory && make # build EventBuilder, Mapper -root -q -l -b -e '.L TrackRecon.C++O' # pre-compile before parallel runs -``` - -Batch runs: `./run_17F.sh`, `./run_27Al.sh`, `./run_tr.sh`. No test suite — validation is visual via ROOT `TBrowser`. - -## Runtime configuration - -`TrackRecon.C::Begin()` reads all config from environment variables (`DATASET`, `reactiondata`, `source_vertex`, `Gain`, `CO2percent`, `A1C1_Z_SCALE`, `A1C1_Z_OFF_QQQ`, `A1C1_Z_OFF_SX3`, etc.). The batch scripts `export` these before launching ROOT. See `README.md` for the full pipeline and file reference. diff --git a/scratch/fit_a1c0_zcorr.C b/scratch/fit_a1c0_zcorr.C new file mode 100644 index 0000000..881a460 --- /dev/null +++ b/scratch/fit_a1c0_zcorr.C @@ -0,0 +1,188 @@ +// fit_a1c0_zcorr.C +// --------------------------------------------------------------------------- +// Linear-regression fit of the anode-only (A1C0) PC-z against the reference z, +// to extract the per-detector scale + offset corrections used by a1c1_zcorr() +// in TrackRecon.C. +// +// WORKFLOW +// 1. Set the 4 correction params to 0 (the default in TrackRecon.C now) and +// run all SOURCE runs so the A1C0 benchmark histograms hold the RAW, +// uncorrected anode-only z vs the reference z. +// 2. hadd the per-run outputs into one file (e.g. Output_av/all.root). +// 3. root -l -b -q 'fit_a1c0_zcorr.C("Output_av/all.root")' +// +// MODEL +// The benchmark TH2s are filled as (x = reference/truth z, y = z_a1c0 raw). +// Profiling y over x and fitting y = A*x + B gives the linear trend. +// a1c1_zcorr does z_corr = z_a1c0*(1-scale) - off, and we want z_corr = x, +// so: scale = 1 - 1/A off = B/A +// (Equivalent to OLS on the scatter; the profile gives clean per-bin means +// and avoids the per-event spread destabilising the fit.) +// +// CAVEAT +// Smaller residual width alone is NOT proof of a correct scale: as scale->1 +// the estimate collapses to a constant. The fit below pins the slope, so the +// corrected A1C0-vs-ref slope should come out ~1 by construction. Always +// re-run with the printed params and confirm the corrected slope is ~1. +// --------------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Recursively search a directory tree for the first TH2 whose name matches. +static TH2 *findTH2(TDirectory *dir, const std::string &name) +{ + // direct hit in this directory + if (TObject *o = dir->Get(name.c_str())) + if (o->InheritsFrom(TH2::Class())) + return static_cast(o); + + TIter next(dir->GetListOfKeys()); + while (TKey *key = static_cast(next())) + { + TObject *obj = key->ReadObj(); + if (obj->InheritsFrom(TDirectory::Class())) + { + if (TH2 *h = findTH2(static_cast(obj), name)) + return h; + } + else if (obj->InheritsFrom(TH2::Class()) && name == obj->GetName()) + { + return static_cast(obj); + } + } + return nullptr; +} + +// Try a list of candidate histogram names, fit the first one found. +// Returns true and sets scale/off on success. +static bool fitOne(TFile *f, const char *label, + const std::vector &candidates, + double fitlo, double fithi, + double &scale, double &off) +{ + TH2 *h2 = nullptr; + std::string used; + for (const auto &nm : candidates) + { + h2 = findTH2(f, nm); + if (h2) + { + used = nm; + break; + } + } + if (!h2) + { + std::cout << "[" << label << "] none of the candidate histograms were found:\n"; + for (const auto &nm : candidates) + std::cout << " " << nm << "\n"; + return false; + } + + // ProfileX: mean of y (z_a1c0) in each x (reference) bin. + TProfile *prof = h2->ProfileX((used + "_pfx").c_str()); + if (prof->GetEntries() < 1) + { + std::cout << "[" << label << "] profile from '" << used << "' is empty.\n"; + return false; + } + + TF1 fpol("fpol", "pol1", fitlo, fithi); + prof->Fit(&fpol, "QR"); // Quiet, restrict to range + double B = fpol.GetParameter(0); // intercept + double A = fpol.GetParameter(1); // slope + double Aerr = fpol.GetParError(1); + + if (A == 0.0) + { + std::cout << "[" << label << "] degenerate slope (A=0), cannot invert.\n"; + return false; + } + + scale = 1.0 - 1.0 / A; + off = B / A; + + std::cout << "------------------------------------------------------------\n"; + std::cout << "[" << label << "] hist: " << used << "\n"; + std::cout << " entries (2D) : " << (Long64_t)h2->GetEntries() << "\n"; + std::cout << " fit range : [" << fitlo << ", " << fithi << "] mm\n"; + std::cout << " z_a1c0 = A*z_ref + B\n"; + std::cout << " A (slope) : " << A << " +/- " << Aerr << "\n"; + std::cout << " B (intercept) : " << B << " mm\n"; + std::cout << " => scale = 1-1/A : " << scale << "\n"; + std::cout << " off = B/A : " << off << " mm\n"; + return true; +} + +void fit_a1c0_zcorr(const char *filename, + double fitlo = -150.0, double fithi = 150.0) +{ + TFile *f = TFile::Open(filename, "READ"); + if (!f || f->IsZombie()) + { + std::cerr << "Cannot open file: " << filename << std::endl; + return; + } + + std::cout << "\n==================== A1C0 z-correction fit ====================\n"; + std::cout << "file: " << filename << "\n"; + + double scale_qqq = 0, off_qqq = 0, scale_sx3 = 0, off_sx3 = 0; + + // QQQ: only a _vs_ref benchmark is filled (A1C2 crossover is the truth). + bool okQ = fitOne(f, "QQQ", + {"Benchmark_QQQ_PCZ_A1C0_Hyb_vs_ref", + "Benchmark_QQQ_PCZ_A1C0_vs_ref", + "Benchmark_QQQ_PCZ_A1C0_Hyb_vs_qqqpczguess"}, + fitlo, fithi, scale_qqq, off_qqq); + + // SX3: prefer the source-run geometric guess (fixed-vertex truth), fall + // back to the A1C2 reference. + bool okS = fitOne(f, "SX3", + {"Benchmark_SX3_PCZ_A1C0_Hyb_vs_sx3pczguess", + "Benchmark_SX3_PCZ_A1C0_Hyb_vs_ref", + "Benchmark_SX3_PCZ_A1C0_vs_sx3pczguess"}, + fitlo, fithi, scale_sx3, off_sx3); + + std::cout << "\n==================== paste-ready output ====================\n"; + if (okQ || okS) + { + std::cout << "// --- TrackRecon.C Begin() dataset block ---\n"; + if (okQ) + { + std::cout << "double a1c1_z_scale_qqq = " << scale_qqq << ";\n"; + std::cout << "double a1c1_z_off_qqq = " << off_qqq << ";\n"; + } + if (okS) + { + std::cout << "double a1c1_z_scale_sx3 = " << scale_sx3 << ";\n"; + std::cout << "double a1c1_z_off_sx3 = " << off_sx3 << ";\n"; + } + std::cout << "\n// --- or as environment overrides ---\n"; + if (okQ) + { + std::cout << "export A1C1_Z_SCALE_QQQ=" << scale_qqq << "\n"; + std::cout << "export A1C1_Z_OFF_QQQ=" << off_qqq << "\n"; + } + if (okS) + { + std::cout << "export A1C1_Z_SCALE_SX3=" << scale_sx3 << "\n"; + std::cout << "export A1C1_Z_OFF_SX3=" << off_sx3 << "\n"; + } + } + else + { + std::cout << "No histograms fitted -- check the file and run with source data.\n"; + } + std::cout << "============================================================\n\n"; + + f->Close(); +} diff --git a/scratch/fit_xy_offset.C b/scratch/fit_xy_offset.C new file mode 100644 index 0000000..1cbe0c1 --- /dev/null +++ b/scratch/fit_xy_offset.C @@ -0,0 +1,148 @@ +// fit_xy_offset.C +// --------------------------------------------------------------------------- +// Tests the "beam/source XY offset (tilted rod)" hypothesis by fitting the +// phi-dependence of the A1C2 vertex-z residual and the A1C0 z residual. +// +// z_reco - z_true = -(d/tan(theta)) * cos(phi - phi_d) +// +// where (d, phi_d) is the transverse source offset. So a cos(phi - phi_d) +// modulation of the residual is the signature of an XY offset; a flat line +// means no offset. The PHASE phi_d must be the SAME for SX3 and QQQ if the +// offset is a real, shared source/beam displacement (the falsifiable test). +// The AMPLITUDE scales as 1/tan(theta), so QQQ (forward) > SX3 (transverse) +// for the same physical d. +// +// USAGE (run ONE source position at a time -- combining z-positions and +// per-run offsets washes out the modulation): +// root -l -b -q 'fit_xy_offset.C("Output_a/results_run018.root")' +// +// Reads the Diag_XYoffset histograms filled by TrackRecon.C: +// Diag_{SX3,QQQ}_A1C2_vtxZ_resid_vs_phi (x=phi deg, y=vtxZ - source_vertex) +// Diag_{SX3,QQQ}_A1C0_zresid_vs_phi (x=phi deg, y=pcz_a1c0 - pcz_ref) +// --------------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +static TH2 *findTH2(TDirectory *dir, const std::string &name) +{ + if (TObject *o = dir->Get(name.c_str())) + if (o->InheritsFrom(TH2::Class())) + return static_cast(o); + TIter next(dir->GetListOfKeys()); + while (TKey *key = static_cast(next())) + { + TObject *obj = key->ReadObj(); + if (obj->InheritsFrom(TDirectory::Class())) + { + if (TH2 *h = findTH2(static_cast(obj), name)) + return h; + } + else if (obj->InheritsFrom(TH2::Class()) && name == obj->GetName()) + return static_cast(obj); + } + return nullptr; +} + +// Fit A + B*cos(x*pi/180 - C) to the ProfileX of a (phi, residual) TH2. +// Returns true and reports the offset constant A, amplitude B, phase C(deg). +static bool fitCos(TFile *f, const char *label, const std::string &hname, + double &A, double &B, double &Cdeg) +{ + TH2 *h2 = findTH2(f, hname); + if (!h2) + { + std::cout << "[" << label << "] histogram not found: " << hname << "\n"; + return false; + } + TProfile *prof = h2->ProfileX((hname + "_pfx").c_str()); + if (prof->GetEntries() < 10) + { + std::cout << "[" << label << "] too few entries to fit (" << prof->GetEntries() << ")\n"; + return false; + } + + // A + B*cos(phi - C), phi in radians (x stored in degrees) + TF1 fc("fc", "[0]+[1]*cos(x*TMath::DegToRad()-[2])", -180, 180); + fc.SetParameters(prof->GetMean(2), prof->GetRMS(2), 0.0); + fc.SetParLimits(1, 0.0, 1e4); // amplitude >= 0 + prof->Fit(&fc, "QR"); + + A = fc.GetParameter(0); + B = fc.GetParameter(1); + Cdeg = fc.GetParameter(2) * TMath::RadToDeg(); + // fold phase into (-180,180] + while (Cdeg > 180.0) Cdeg -= 360.0; + while (Cdeg <= -180.0) Cdeg += 360.0; + + double chi2ndf = (fc.GetNDF() > 0) ? fc.GetChisquare() / fc.GetNDF() : -1.0; + std::cout << "------------------------------------------------------------\n"; + std::cout << "[" << label << "] " << hname << "\n"; + std::cout << " constant A : " << A << " mm (mean residual)\n"; + std::cout << " amplitude B : " << B << " mm +/- " << fc.GetParError(1) << "\n"; + std::cout << " phase phi_d : " << Cdeg << " deg (offset direction)\n"; + std::cout << " chi2/ndf : " << chi2ndf << "\n"; + return true; +} + +void fit_xy_offset(const char *filename) +{ + TFile *f = TFile::Open(filename, "READ"); + if (!f || f->IsZombie()) + { + std::cerr << "Cannot open file: " << filename << std::endl; + return; + } + + std::cout << "\n================= XY-offset (rod-tilt) test =================\n"; + std::cout << "file: " << filename << "\n"; + std::cout << "Model: residual = A - (d/tan(theta))*cos(phi - phi_d)\n"; + std::cout << " -> amplitude B ~ d/tan(theta), phase phi_d = offset azimuth.\n"; + std::cout << " -> phi_d must MATCH between SX3 and QQQ for a real offset.\n"; + + double A, B, Cdeg; + double phiSX3_vtx = 0, phiQQQ_vtx = 0; + bool okSvtx = false, okQvtx = false; + + std::cout << "\n--- A1C2 vertex-z residual (cleanest: theta-projection) ---\n"; + if ((okSvtx = fitCos(f, "SX3 vtxZ", "Diag_SX3_A1C2_vtxZ_resid_vs_phi", A, B, Cdeg))) + phiSX3_vtx = Cdeg; + if ((okQvtx = fitCos(f, "QQQ vtxZ", "Diag_QQQ_A1C2_vtxZ_resid_vs_phi", A, B, Cdeg))) + phiQQQ_vtx = Cdeg; + + std::cout << "\n--- A1C0 z residual (anode-only minus A1C2 ref) ---\n"; + fitCos(f, "SX3 A1C0", "Diag_SX3_A1C0_zresid_vs_phi", A, B, Cdeg); + fitCos(f, "QQQ A1C0", "Diag_QQQ_A1C0_zresid_vs_phi", A, B, Cdeg); + + std::cout << "\n==================== verdict ====================\n"; + if (okSvtx && okQvtx) + { + double dphi = phiSX3_vtx - phiQQQ_vtx; + while (dphi > 180.0) dphi -= 360.0; + while (dphi <= -180.0) dphi += 360.0; + std::cout << "vertex-z phase: SX3 = " << phiSX3_vtx + << " deg, QQQ = " << phiQQQ_vtx << " deg\n"; + std::cout << "phase difference = " << dphi << " deg\n"; + if (TMath::Abs(dphi) < 30.0) + std::cout << "=> phases AGREE -> consistent with a real shared XY offset" + " (rod tilt). Offset direction ~ " << phiSX3_vtx << " deg.\n"; + else + std::cout << "=> phases DISAGREE -> not a single shared XY offset;" + " likely per-detector calibration instead.\n"; + } + else + { + std::cout << "Could not fit both detectors -- check the run has A1C2 stats.\n"; + } + std::cout << "=================================================\n\n"; + + f->Close(); +}