diff --git a/CLAUDE.md b/CLAUDE.md index 05f2aaa..53b585b 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -1,21 +1,11 @@ # ANASEN Analysis — working notes for Claude ## Agent behavior -- **Architecture & Remotes:** - - `upstream` points to the canonical public repo (read-only reference). - - `origin` points to your private sandbox repo (read/write scratchpad). - - **NEVER attempt to push to `upstream`.** All upstream synchronization is strictly one-way (`git fetch upstream` -> merge). -- **Upstream Sync Routine:** - - Before applying requested features or edits, sync with `upstream` if requested. - - Never pull into a dirty working tree. Use: `git stash` -> `git fetch upstream` -> `git merge upstream/devel_vignesh` -> `git stash pop`. - - If a merge conflict occurs in analysis macros (`TrackRecon.C`, `MakeVertex.C`, kinematics scripts, or calibration tables), **STOP immediately**. Do not attempt to guess physics logic; flag the exact conflict lines and await instruction. -- **Commits & Pushes to `origin`:** - - Commit only when a coherent, discrete change is finished or when explicitly asked. Write concise, conventional commit messages. - - Push **only** to `origin` (`git push origin `). Push after completing a requested unit of work so the sandbox stays backed up. - ALWAYS output a brief, step-by-step plan before modifying files or running commands. +- Commit only when explicitly asked. Do not push to remote on your own initiative; pushing is allowed only when the session's instructions name a branch to push to (e.g. a `claude/*` working branch), and never to `devel_vignesh` or any other shared branch. - **Targeted edits only.** Never regenerate a whole file to make a local change. -- **Never read `MakeVertex.C`** -- **Never read `TrackRecon.C` end to end** (~4,250 lines; ~63k tokens each, and it stays in the prefix for the rest of the session). Grep for the symbol, then Read with offset/limit around the hit. If a question genuinely needs the whole file, say so first. Do not explore the codebase open-endedly. +- **Never read `MakeVertex.C`** (4,003 lines, 240 KB). Same rule as `TrackRecon.C` below if you must touch it: grep, then Read with offset/limit. +- **Never read `TrackRecon.C` end to end** (4,268 lines, 218 KB; ~60k tokens, and it stays in the prefix for the rest of the session). Grep for the symbol, then Read with offset/limit around the hit. If a question genuinely needs the whole file, say so first. Do not explore the codebase open-endedly. - **Do not spawn subagents** for work that can be done inline — each starts cold and re-derives context already loaded. Worth it only for wide parallel searches where the conclusion is all that's needed. - **You cannot verify your own work here.** No test suite, no CI, no data in the repo; validation is visual, via ROOT `TBrowser`, on a machine that has the data. Say explicitly when a change is unverified — never report a physics change as working because it compiled. The dangerous bugs are silent: a flipped sign in a kinematics term, an inverted spline direction, a wrong argument order, an off-by-one in a wire mask. None of them throw. The run completes, the plot renders, and the spectrum is wrong. - **Review once per coherent change, not per edit.** Scale depth to blast radius: scratch macros, shell, and docs are shallow; `TrackRecon.C`, `MakeVertex.C`, `Armory/`, or any calibration fit is deep and should tolerate uncertain findings — silent failure is what's being hunted. Run `/security-review` only for file I/O or shell interpolation in the `run_*.sh` drivers; it is not a physics tool. @@ -34,8 +24,7 @@ All config is `getenv`-driven; the run scripts `export` it. Defaults are in `Tra Physics / results-affecting: `DATASET` `reactiondata` `CO2percent` `pressure_in_torr` `CATHODE_GAIN` `PC_ENERGY_CALIBRATION` `source_vertex` `CUTLIST` `DITHER_SIGMA` `RNG_SEED` `BEAM_AXIS_X` `BEAM_AXIS_Y` -`timecut_low` `timecut_high` `DISABLE_BAD_ANODE_WIRES` `A1C1_LOWBAND_RFACTOR` `A1C1_CFRAC_SPLIT` -`A1C1_ANODEE_COFF` `A1C1_ANODEE_REF` `A1C1_MISSING_FMAX` +`timecut_low` `timecut_high` `DISABLE_BAD_ANODE_WIRES` `A1C1_LOWBAND_RFACTOR` `A1C1_Z_SCALE_QQQ` `A1C1_Z_SCALE_SX3` `A1C1_Z_OFF_QQQ` `A1C1_Z_OFF_SX3` Plumbing: `OUT_DIR` `RUN_NUMBER` `FLUSH_BARRIER` `MAX_RSS_MB` `MEMCHECK_STRIDE` @@ -44,9 +33,9 @@ Plumbing: `OUT_DIR` `RUN_NUMBER` `FLUSH_BARRIER` `MAX_RSS_MB` `MEMCHECK_STRIDE` It multiplies `catima.dedx` when the lookup tables are generated. Changing it invalidates every `.dat` table — re-run `Eloss.py` (the run scripts do this automatically before processing). Tables are keyed `eloss_calculations/_lookup_MeV_

torr_pc.dat` (e.g. `alpha_lookup_50MeV_250torr_3pc.dat`). ## Footguns -- **`HistPlotter` keys `oMap` by histogram NAME only.** The folder argument does not namespace. Two `Fill` calls with the same name in different folders merge silently — no warning, no error, wrong plot. Always make names unique unless required. -- **Si–PC coincidence sits at negative dt.** The real band is roughly `-450 < t_si - t_pc < -200`, not near zero. Do not "correct" a timing gate toward zero. -- **RNG**: one shared `anasenRandom` (`TRandom3`, seed 4357, `RNG_SEED` to override) — do not introduce fresh `TRandom3(0)` instances; it breaks reproducibility. +- **`HistPlotter` keys `oMap` by histogram NAME only.** The folder argument does not namespace. Two `Fill` calls with the same name and the *same* folder merge silently — no warning, no error, wrong plot. A same-name/different-folder (or toplevel-vs-folder) clash does print to `stderr` and then merges anyway, but the run scripts send stdout+stderr to `/dev/null`, so in a batch run it is invisible either way. Always make names unique. +- **Si–PC coincidence sits at negative dt.** The band is centred well below zero, not near it. The one gate of record is `siPcCoincident()` in `TrackRecon.C` (`kSiPcDtMin = -500`, `kSiPcDtMax = 150` ns); anything else is a plotting range, not a cut. Do not "correct" a timing gate toward zero. +- **RNG**: in `TrackRecon.C`, one shared `anasenRandom` (`TRandom3`, seed 4357, `RNG_SEED` to override), aliased by reference wherever a local `rnd`/`rand` appears — do not introduce fresh `TRandom3(0)` instances; it breaks reproducibility. `MakeVertex.C` predates this and still has ~12 independent default-seeded `TRandom3` objects; do not assume its output is reproducible against `TrackRecon.C`. - **Pass By Reference in Event Loops:** Always iterate via `const auto&` over cluster/hit vectors to prevent expensive struct copies in multi-million event runs. ## Invariants to check on any physics change diff --git a/TrackRecon.C b/TrackRecon.C index 466d88c..ada6dbf 100644 --- a/TrackRecon.C +++ b/TrackRecon.C @@ -42,8 +42,8 @@ Int_t colors[40] = { bool process_alpha_proton_scattering = false, doMiscHistograms = true, doRawHistos = false, - doPCSX3ClusterAnalysis = false, - doPCQQQClusterAnalysis = false, + doPCSX3ClusterAnalysis = true, + doPCQQQClusterAnalysis = true, doOldAnalysis = false, BenchMark = false, onewire_analysis = true, @@ -75,6 +75,9 @@ double source_vertex = 53.0, a1c1_z_off_sx3 = 2.52614, beam_axis_x = 0.0, beam_axis_y = 0.0, + beam_axis_z0 = 0.0, // reference z at which (beam_axis_x, beam_axis_y) is specified + beam_tilt_x = 0.0, // dx/dz slope of the beam axis + beam_tilt_y = 0.0, // dy/dz slope of the beam axis ta_foil_z_mm = 0.0, alpha_source_mev = 5.486; @@ -122,23 +125,111 @@ inline SiPcPid classifyByAnodeDe(double anodeE_MeV) return (anodeE_MeV < kProtonAlphaAnodeDeGate_MeV) ? SiPcPid::kProton : SiPcPid::kAlpha; } -inline TVector3 beamVertex(const TVector3 &si, const TVector3 &dir) +inline TVector3 beamAxisOrigin() { - double d = dir.X() * dir.X() + dir.Y() * dir.Y(); - double t = (d > 0.0) ? -((si.X() - beam_axis_x) * dir.X() + (si.Y() - beam_axis_y) * dir.Y()) / d : 0.0; - return si + t * dir; + return TVector3(beam_axis_x, beam_axis_y, beam_axis_z0); } -inline double beamPerp(const TVector3 &p) +inline TVector3 beamAxisDirUnit() { - return TMath::Sqrt((p.X() - beam_axis_x) * (p.X() - beam_axis_x) + (p.Y() - beam_axis_y) * (p.Y() - beam_axis_y)); + return TVector3(beam_tilt_x, beam_tilt_y, 1.0).Unit(); } -// A point on the beam axis at height z. Every theta/phi reference point used to be -// spelled beamAxisPoint(z), which silently ignored BEAM_AXIS_X/Y even though -// Begin() prints them as configured parameters and pcEnergyCalibrationAccumulate -// already built its source_pos the correct way. + inline TVector3 beamAxisPoint(double z) { - return TVector3(beam_axis_x, beam_axis_y, z); + double dz = z - beam_axis_z0; + return TVector3(beam_axis_x + beam_tilt_x * dz, + beam_axis_y + beam_tilt_y * dz, + z); +} + +inline TVector3 pocaOnTrack(const TVector3 &si, const TVector3 &dir, + const TVector3 &p, const TVector3 &d) +{ + TVector3 e = dir; + double ee = e.Dot(e); + double de = d.Dot(e); + double denom = ee - de * de; // = |e|^2 * sin^2(angle between beam and track) + double t = 0.0; + if (denom > 1e-18 * (ee > 0.0 ? ee : 1.0)) + { + TVector3 w = si - p; + double dw = d.Dot(w); + double ew = e.Dot(w); + t = (de * dw - ew) / denom; + } + return si + t * e; +} + +inline TVector3 beamVertex(const TVector3 &si, const TVector3 &dir) +{ + return pocaOnTrack(si, dir, beamAxisOrigin(), beamAxisDirUnit()); +} +inline TVector3 beamVertexNominal(const TVector3 &si, const TVector3 &dir) +{ + return pocaOnTrack(si, dir, TVector3(0.0, 0.0, 0.0), TVector3(0.0, 0.0, 1.0)); +} + +inline double beamPerp(const TVector3 &p) +{ + TVector3 p0 = beamAxisOrigin(); + TVector3 d = beamAxisDirUnit(); + TVector3 w = p - p0; + TVector3 perp = w - d * w.Dot(d); + return perp.Mag(); +} + +inline void fillBeamProfile(HistPlotter *plotter, const TVector3 &vertex, + const TVector3 &si, const TVector3 &dir, const std::string &tag, + bool axisSafe) +{ + const std::string folder = "BeamAxis"; + const std::string all = "beamAxis_all_"; + const std::string one = "beamAxis_" + tag + "_"; + double vx = vertex.X(), vy = vertex.Y(), vz = vertex.Z(); + double pocaDist = beamPerp(vertex); + plotter->Fill2D(all + "vertexX_vs_Z", 250, -450, 50, 200, -50, 50, vz, vx, folder); + plotter->Fill2D(all + "vertexY_vs_Z", 250, -450, 50, 200, -50, 50, vz, vy, folder); + plotter->Fill1D(all + "pocaDist", 400, 0, 100, pocaDist, folder); + // Per-branch copy so anomalies can be traced to a specific reconstruction path. + plotter->Fill2D(one + "vertexX_vs_Z", 250, -450, 50, 200, -50, 50, vz, vx, folder); + plotter->Fill2D(one + "vertexY_vs_Z", 250, -450, 50, 200, -50, 50, vz, vy, folder); + plotter->Fill1D(one + "pocaDist", 400, 0, 100, pocaDist, folder); + const double zLo = -440.0, zHi = 40.0; + const int nSlice = 16; + const double sliceW = (zHi - zLo) / nSlice; // 40 mm + if (vz >= zLo && vz < zHi) + { + int is = static_cast((vz - zLo) / sliceW); + if (is >= 0 && is < nSlice) + { + char buf[8]; + snprintf(buf, sizeof(buf), "%02d", is); + plotter->Fill2D(all + "vertexXY_z" + buf, 200, -50, 50, 200, -50, 50, vx, vy, folder); + } + } + + double dz = dir.Z(); + double vzBin = beamVertexNominal(si, dir).Z(); + if (axisSafe && TMath::Abs(dz) > 1e-6 && vzBin >= zLo && vzBin < zHi) + { + int k = static_cast((vzBin - zLo) / sliceW); + if (k >= 0 && k < nSlice) + { + double zPlane = zLo + (k + 0.5) * sliceW; + double s = (zPlane - si.Z()) / dz; + double cx = si.X() + s * dir.X(); + double cy = si.Y() + s * dir.Y(); + if (TMath::Abs(cx) < 100.0 && TMath::Abs(cy) < 100.0) + { + char kbuf[8]; + snprintf(kbuf, sizeof(kbuf), "%02d", k); + plotter->Fill1D(all + "crossX_z" + kbuf, 800, -100, 100, cx, folder); + plotter->Fill1D(all + "crossY_z" + kbuf, 800, -100, 100, cy, folder); + plotter->Fill1D(one + "crossX_z" + kbuf, 800, -100, 100, cx, folder); + plotter->Fill1D(one + "crossY_z" + kbuf, 800, -100, 100, cy, folder); + } + } + } } struct PCPath @@ -587,7 +678,16 @@ void TrackRecon::Begin(TTree * /*tree*/) beam_axis_x = std::atof(getenv("BEAM_AXIS_X")); if (getenv("BEAM_AXIS_Y")) beam_axis_y = std::atof(getenv("BEAM_AXIS_Y")); - std::cout << "Beam-axis origin (x,y) = (" << beam_axis_x << ", " << beam_axis_y << ") mm" << std::endl; + if (getenv("BEAM_AXIS_Z0")) + beam_axis_z0 = std::atof(getenv("BEAM_AXIS_Z0")); + if (getenv("BEAM_TILT_X")) + beam_tilt_x = std::atof(getenv("BEAM_TILT_X")); + if (getenv("BEAM_TILT_Y")) + beam_tilt_y = std::atof(getenv("BEAM_TILT_Y")); + std::cout << "Beam axis line: (x,y) = (" << beam_axis_x << ", " << beam_axis_y + << ") mm at z0 = " << beam_axis_z0 + << " mm, tilts (dx/dz, dy/dz) = (" << beam_tilt_x << ", " << beam_tilt_y << ")" + << std::endl; if (doPCEnergyCalibration) std::cout << "PC energy calibration ON: alpha source = " << alpha_source_mev << " MeV, source position = (" << beam_axis_x << ", " << beam_axis_y << ", " << source_vertex @@ -2045,6 +2145,8 @@ void protonAlphaHistograms(HistPlotter *plotter, const std::vector &QQQ_E // purely longitudinal direction. TVector3 r_rhoMin_fix = beamVertex(x1, v); double vertex_z = r_rhoMin_fix.Z(); + fillBeamProfile(plotter, r_rhoMin_fix, x1, v, "apCoinc", + pcevent.multi1 == 1 && pcevent.multi2 == 2); double theta_q = (qqqevent.pos - beamAxisPoint(vertex_z)).Theta(); // double theta_q = (qqqevent.pos - r_rhoMin_fix).Theta(); double sinTheta_customV = TMath::Sin(theta_q); @@ -2485,6 +2587,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector &QQQ_Ev // beam axis at (0,0), silently ignoring BEAM_AXIS_X/Y, and had no guard for a // purely longitudinal direction. TVector3 r_rhoMin_fix = beamVertex(x1, v); + fillBeamProfile(plotter, r_rhoMin_fix, x1, v, "sx3a1c2", true); // a1c2 only, gated above plotter->Fill1D("VertexRecon_pczfix_sx3", 800, -300, 300, r_rhoMin_fix.Z(), "Vertex_Reconstruction"); plotter->Fill1D("VertexRecon_pczfix", 800, -300, 300, r_rhoMin_fix.Z(), "Vertex_Reconstruction"); plotter->Fill1D("pczfix_A1C2_1d_sx3", 600, -200, 200, pcz_fix, "PCZ_Recon"); @@ -3036,6 +3139,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector &QQQ_Ev // beam axis at (0,0), silently ignoring BEAM_AXIS_X/Y, and had no guard for a // purely longitudinal direction. TVector3 r_rhoMin_fix = beamVertex(x1, v); + fillBeamProfile(plotter, r_rhoMin_fix, x1, v, "qqqa1c2", true); // a1c2 only, gated above double sinTheta_customV = TMath::Sin((qqqevent.pos - beamAxisPoint(r_rhoMin_fix.Z())).Theta()); plotter->Fill2D("dE3_E_CathodeQQQR_A1C2_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2 * sinTheta_customV, "PID_dE_E"); @@ -3593,6 +3697,8 @@ void protonAlphaElastic_core(HistPlotter *plotter, const std::vector &Si_ TVector3 x2f(pcXY.X(), pcXY.Y(), pcz_fix); TVector3 r_rhoMin_fix = beamVertex(sievent.pos, x2f - sievent.pos); double vertex_z = r_rhoMin_fix.Z(); + const bool axisSafe = (multi2 == 2) || (multi1 == 2 && multi2 == 0); + fillBeamProfile(plotter, r_rhoMin_fix, sievent.pos, x2f - sievent.pos, "elastic_" + det, axisSafe); if (vertex_z < z_entrance) return; double theta = (sievent.pos - r_rhoMin_fix).Theta(); @@ -3629,7 +3735,7 @@ void protonAlphaElastic_core(HistPlotter *plotter, const std::vector &Si_ plotter->Fill1D(rx + "_VertexReconZ" + ejtag + sfx, 800, -400, 400, vertex_z, pmlabel); plotter->Fill2D(rx + "_VertexReconXY" + ejtag + sfx, 200, -100, 100, 200, -100, 100, r_rhoMin_fix.X(), r_rhoMin_fix.Y(), pmlabel); plotter->Fill2D(rx + "_Ef_vs_theta" + ejtag + sfx, 100, 0, 180, 800, 0, 10, theta * 180 / M_PI, Efix, pmlabel); - plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + sfx, 360, 0, 180, 800, -10, 10, theta * 180 / M_PI, Ex, pmlabel); + plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + sfx, 720, 0, 180, 800, -10, 10, theta * 180 / M_PI, Ex, pmlabel); plotter->Fill2D(rx + "_Ex_vs_phi" + ejtag + sfx, 180, -180, 180, 800, -10, 10, sievent.pos.Phi() * 180 / M_PI, Ex, pmlabel); for (const auto &pcevent : PC_Events) @@ -3970,6 +4076,9 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector &Si_ TVector3 x2f(pcXY.X(), pcXY.Y(), pcz_fix); TVector3 r_rhoMin_fix = beamVertex(sievent.pos, x2f - sievent.pos); double vertex_z = r_rhoMin_fix.Z(); + const bool axisSafe = (topo1 == "a1c2fix") || (topo1 == "a2c0") || (topo2 == "a1c1_inband"); + fillBeamProfile(plotter, r_rhoMin_fix, sievent.pos, x2f - sievent.pos, + "reaction_" + globaltag + "_" + det, axisSafe); if (beamPerp(r_rhoMin_fix) > perp_cut || vertex_z < z_entrance) return; @@ -4029,7 +4138,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector &Si_ plotter->Fill1D(rx + "_Ex_from" + ejtag + t + sfx, 600, -10, 20, Ex, pmlabel); plotter->Fill2D(rx + "_VertexReconZ_vs_Ef" + ejtag + t + sfx, 800, -400, 400, 800, 0, ef_max, vertex_z, Efix, pmlabel); plotter->Fill2D(rx + "_VertexReconZ_vs_Ex" + ejtag + t + sfx, 800, -400, 400, 600, -10, 20, vertex_z, Ex, pmlabel); - plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + t + sfx, 360, 0, 180, 600, -10, 20, theta * 180 / M_PI, Ex, pmlabel); + plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + t + sfx, 720, 0, 180, 600, -10, 20, theta * 180 / M_PI, Ex, pmlabel); if (ebeam_kin_MeV > 0.0) plotter->Fill2D(rx + "_BeamEnergy_ETrack_vs_EKin" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5, @@ -4069,6 +4178,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector &Si_ } plotter->Fill1D(rx + "_pczfix" + sfx, 600, -300, 300, pcz_fix, pmlabel); plotter->Fill2D(rx + "_Ef_vs_theta" + ejtag + sfx, 100, 0, 180, 800, 0, ef_max, theta * 180 / M_PI, Efix, pmlabel); + plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + sfx, 720, 0, 180, 800, -20, 20, theta * 180 / M_PI, Ex, pmlabel); plotter->Fill2D(rx + "_Ex_vs_phi" + ejtag + sfx, 45, -180, 180, 600, -10, 20, phi * 180 / M_PI, Ex, pmlabel); for (const auto &pcevent : PC_Events) @@ -4107,7 +4217,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector &Si_ plotter->Fill2D(rx + "_dEgasCalib_vs_VertexZ" + ejtag + sfx, 800, -400, 400, 800, 0, 0.6, vertex_z, anodeE_MeV, pmlabel); plotter->Fill2D(rx + "_dEgasRaw_vs_VertexZ" + ejtag + sfx, 800, -400, 400, 800, 0, 20000, vertex_z, anodeE, pmlabel); plotter->Fill2D(rx + "_dEgasRaw_vs_theta" + ejtag + sfx, 180, 0, 180, 800, 0, 20000, theta * 180 / M_PI, anodeE, pmlabel); - plotter->Fill2D(rx + "_dEgasCalib_vs_theta" + ejtag + sfx, 360, 0, 180, 800, 0, 0.6, theta * 180 / M_PI, anodeE_MeV, pmlabel); + plotter->Fill2D(rx + "_dEgasCalib_vs_theta" + ejtag + sfx, 720, 0, 180, 800, 0, 0.6, theta * 180 / M_PI, anodeE_MeV, pmlabel); plotter->Fill2D(rx + "_dEgasCalib_vs_phi" + ejtag + sfx, 90, -180, 180, 800, 0, 0.6, phi * 180 / M_PI, anodeE_MeV, pmlabel); // if (anodeCh >= 0) // plotter->Fill2D(rx + "_dEgasCalib_vs_E" + ejtag + sfx + "_anode" + pad2(anodeCh), 400, 0, ef_max, 800, 0, 0.6, sievent.Energy1, anodeE_MeV, pmlabel); diff --git a/run_27Al.sh b/run_27Al.sh index f8af870..2f9b383 100644 --- a/run_27Al.sh +++ b/run_27Al.sh @@ -11,14 +11,12 @@ export CATHODE_GAIN=3.0 export source_vertex=-200.0 export DEDX_SCALE=0.89 export CUTLIST=cuts_list.txt -export BEAM_AXIS_X=0 -export BEAM_AXIS_Y=0 -# Tilted-beam parameters: (x, y) at BEAM_AXIS_Z0 (default 0), plus dx/dz and dy/dz slopes. -# Default (zeros) reproduces the original vertical beam bit-for-bit. Fit them -# with scratch/iterate_beam.sh + scratch/FitBeamAxis.C. -export BEAM_AXIS_Z0=${BEAM_AXIS_Z0:-0} -export BEAM_TILT_X=${BEAM_TILT_X:-0} -export BEAM_TILT_Y=${BEAM_TILT_Y:-0} + +export BEAM_AXIS_Z0=-200 +export BEAM_AXIS_X=-0.136398 +export BEAM_AXIS_Y=0.034824 +export BEAM_TILT_X=-0.00083798 +export BEAM_TILT_Y=0.00191377 echo "Pre-compiling TrackRecon.C safely on a single core..." root -q -l -b -e '.L TrackRecon.C++O' @@ -64,11 +62,11 @@ export -f process_run echo "Starting parallel processing..." # time parallel --bar -j 12 process_run ::: {24..41} 44 45 46 {50..59} - time parallel --bar -j 12 process_run ::: {24..41} + # time parallel --bar -j 12 process_run ::: {24..41} # time parallel --bar -j 10 process_run ::: 44 45 46 {50..59} # time parallel --bar -j 1 process_run ::: 48 # pc without coincidence # mv "${CURRENT_OUT_DIR}/results_run048.root" "Output_27Al_run48/." - # time parallel --bar -j 10 process_run ::: 62 63 66 67 68 73 74 {78..89} + time parallel --bar -j 10 process_run ::: {24..41} 44 45 46 {50..59} 62 63 66 67 68 73 74 {78..89} # time parallel --bar -j 4 process_run ::: 62 63 66 67 68 # time parallel --bar -j 1 process_run ::: 73 # time parallel --bar -j 1 process_run ::: 74 diff --git a/scratch/FitBeamAxis.C b/scratch/FitBeamAxis.C new file mode 100644 index 0000000..2dd920c --- /dev/null +++ b/scratch/FitBeamAxis.C @@ -0,0 +1,261 @@ +// FitBeamAxis.C +// --------------------------------------------------------------------------- +// Single-pass beam-axis line fit. Replaces the iterative +// FitBeamAxis.C + iterate_beam.sh scheme, which diverges. +// +// WHY THE ITERATIVE SCHEME FAILS +// ------------------------------ +// FitBeamAxis.C fits beamAxis_all_vertexX_vs_Z and feeds the residual back as a +// correction to BEAM_AXIS_X / BEAM_TILT_X. But that vertex comes from +// beamVertex(), which returns the track's point of closest approach to the +// ASSUMED beam line -- so the reconstructed vertex is pulled onto whatever axis +// the run was configured with, by construction. Its residual therefore measures +// acceptance asymmetry and geometry artifacts far more than it measures "how far +// off is my assumed axis". Feeding that back adds a roughly fixed bias each pass +// with no restoring force, so the parameters drift monotonically and never +// converge. Measured behaviour: step L2 grew 1.48 -> 2.65 over 6 passes at full +// step, and still grew (1.27 -> 1.50 over 10 passes) at relax=0.3 -- damping only +// scales the drift rate, which is the signature of a wrong-direction update +// rather than overshoot. +// +// WHAT THIS DOES INSTEAD +// ---------------------- +// TrackRecon.C's fillBeamProfile() also histograms where each track CROSSES a set +// of fixed z-planes (beamAxis__crossX_z**, crossY_z**). A track crossing is a +// property of the track alone -- it never references the assumed beam axis. Tracks +// originate on the true beam line and fan out with roughly uniform azimuth, so at +// each plane the crossing distribution centres on the true beam position at that z. +// Take a robust centre per plane, line-fit those centres against z, and the axis +// falls out directly. One pass, no iteration, no feedback. +// +// WHICH TRACKS FEED THE PLANES +// ---------------------------- +// fillBeamProfile() only sends a track to the crossing planes when its PC z is a +// real measurement: a1c2 (cathode charge division), a2c0 (two wires, geometry +// fixes z), and a1c1 whose a1c1_cfrac_pcz() solve landed in band. Bare a1c1 and +// a1c0 are excluded -- their z is a Gaussian dither about the raw wire z, which +// puts them on the assumed axis by construction and would make this fit return +// whatever axis the run was configured with. The vertexX_vs_Z and pocaDist +// diagnostics in the same folder are still filled for every topology, so their +// entry counts will not match the crossing planes'. Do not fit those. +// +// CAVEAT -- azimuthal symmetry. The per-plane centre is only unbiased if +// acceptance is uniform in phi. Dead channels break that (this analysis already +// excludes QQQ wedge 48/49, ring 63, and several anode wires). The median is used +// rather than the mean to blunt the resulting asymmetric tails, but a strongly +// lopsided acceptance will still pull the centre. If the fitted axis looks +// implausible, check the crossX/crossY distributions at a few planes for +// visible asymmetry before trusting the number. +// +// USAGE: +// root -l -b -q 'scratch/FitBeamAxisCrossings.C("Output_27Al/Output_27Al.root")' +// root -l -b -q 'scratch/FitBeamAxisCrossings.C("out.root", -440, 40, "all")' +// root -l -b -q 'scratch/FitBeamAxisCrossings.C("out.root", -440, 40, "elastic")' +// +// Arguments: +// filename TrackRecon.C output ROOT file +// zmin/zmax only planes whose z falls in this window are used in the line fit. +// Defaults span every plane; narrow it to exclude planes you have +// reason to distrust, not to chase a better chi2. +// tag which branch's crossings to use. "all" pools everything; the +// per-branch tags are "reaction__" (e.g. +// "reaction_m27Alax_qqq", "..._sx3") and "elastic_" ("elastic_QQQ", +// "elastic_SX3"). Fitting the qqq and sx3 tags separately is the check +// that a tilt is the beam and not one branch's acceptance: same beam, +// different geometry and phi coverage, so they must agree. +// Which tags exist at all depends on run configuration -- on a dataset +// where doPCSX3ClusterAnalysis / doPCQQQClusterAnalysis / +// process_alpha_proton_scattering are false and OUT_DIR is not +// Output_p, the sx3a1c2, qqqa1c2, apCoinc and elastic_* tags are never +// filled and "all" is exactly the reaction_* tags pooled. +// +// NOTE on chi2/ndf: medianAndError reports a STATISTICS-ONLY uncertainty +// (~8 um at a 4e5-entry plane) while observed plane-to-plane scatter in a +// well-behaved region is ~0.4 mm, i.e. ~50x larger. chi2/ndf is therefore +// inflated by a missing systematic term and will not approach 1 even for a +// correct axis. Judge the fit on whether the medians lie on a straight line, +// not on chi2/ndf alone. +// +// Prints the fitted BEAM_AXIS_X/Y and BEAM_TILT_X/Y as ready-to-paste export +// lines. These are ABSOLUTE values, not increments -- set them directly in +// run_27Al.sh; do NOT add them to the current values. +// --------------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Must match fillBeamProfile() in TrackRecon.C. The range follows the physical +// vertex acceptance (z_entrance = -454 mm up to the +100 mm cut), not a +// symmetric window about 0, and stops well below the QQQ plane at z = 105 -- +// see the comment in fillBeamProfile() for why planes near a detector are +// unusable. +static const double kZLo = -440.0; +static const double kZHi = 40.0; +static const int kNSlice = 16; + +static TH1 *findTH1(TDirectory *dir, const std::string &name) +{ + if (TObject *o = dir->Get(name.c_str())) + if (o && o->InheritsFrom(TH1::Class()) && !o->InheritsFrom("TH2")) + return static_cast(o); + TIter next(dir->GetListOfKeys()); + while (TKey *key = static_cast(next())) + { + TObject *obj = key->ReadObj(); + if (!obj) + continue; + if (obj->InheritsFrom(TDirectory::Class())) + { + if (TH1 *h = findTH1(static_cast(obj), name)) + return h; + } + else if (obj->InheritsFrom(TH1::Class()) && !obj->InheritsFrom("TH2") && name == obj->GetName()) + return static_cast(obj); + } + return nullptr; +} + +// Robust centre of a crossing distribution: the median. Chosen over the mean +// because the crossing distribution has long tails (tracks originating far from +// this plane) that are not necessarily symmetric once dead channels bite. +// Uncertainty on the median ~ 1.253 * sigma / sqrt(N), with sigma taken from the +// interquartile range (IQR/1.349) so a few wild tails don't inflate it. +static bool medianAndError(TH1 *h, double &med, double &err, double &nEff, int minEntries = 500) +{ + if (!h) + return false; + double n = h->Integral(); + if (n < minEntries) + return false; + double probs[3] = {0.25, 0.5, 0.75}; + double q[3] = {0, 0, 0}; + h->GetQuantiles(3, q, probs); + med = q[1]; + double sigma = (q[2] - q[0]) / 1.349; // IQR -> Gaussian-equivalent sigma + if (sigma <= 0.0) + return false; + err = 1.253 * sigma / TMath::Sqrt(n); + nEff = n; + return true; +} + +static bool fitOneProjection(TFile *f, const std::string &base, const char *coord, + double zmin, double zmax, double z0, + double &axis0, double &tilt, double &axis0_err, double &tilt_err, + double &chi2ndf, int &nUsed) +{ + std::vector zs, cs, ces, zes; + std::cout << " plane z (mm) median (mm) +/- entries\n"; + for (int k = 0; k < kNSlice; k++) + { + double sliceW = (kZHi - kZLo) / kNSlice; + double zPlane = kZLo + (k + 0.5) * sliceW; + if (zPlane < zmin || zPlane > zmax) + continue; + char kbuf[8]; + snprintf(kbuf, sizeof(kbuf), "%02d", k); + TH1 *h = findTH1(f, base + "cross" + coord + "_z" + kbuf); + double med = 0, err = 0, nEff = 0; + if (!medianAndError(h, med, err, nEff)) + { + std::cout << " " << zPlane << "\t(skipped -- missing or too few entries)\n"; + continue; + } + std::cout << " " << zPlane << "\t\t" << med << "\t" << err << "\t" << (long)nEff << "\n"; + zs.push_back(zPlane); + zes.push_back(0.0); + cs.push_back(med); + ces.push_back(err); + } + nUsed = static_cast(zs.size()); + if (nUsed < 3) + { + std::cout << " [" << coord << "] only " << nUsed << " usable planes -- need >= 3 for a line fit\n"; + return false; + } + + TGraphErrors g(nUsed, zs.data(), cs.data(), zes.data(), ces.data()); + TF1 fL("fL", Form("[0]+[1]*(x-%.8g)", z0), zmin, zmax); + fL.SetParameters(cs[nUsed / 2], 0.0); + if (g.Fit(&fL, "QRN") != 0) + { + std::cout << " [" << coord << "] line fit failed\n"; + return false; + } + axis0 = fL.GetParameter(0); + tilt = fL.GetParameter(1); + axis0_err = fL.GetParError(0); + tilt_err = fL.GetParError(1); + chi2ndf = (fL.GetNDF() > 0) ? fL.GetChisquare() / fL.GetNDF() : -1.0; + return true; +} + +void FitBeamAxis(const char *filename, double zmin = -440.0, double zmax = 40.0, + const char *tag = "all") +{ + TFile *f = TFile::Open(filename, "READ"); + if (!f || f->IsZombie()) + { + std::cerr << "Cannot open file: " << filename << std::endl; + return; + } + + double z0 = 0.0; + if (const char *s = std::getenv("BEAM_AXIS_Z0")) + z0 = std::atof(s); + + std::string base = std::string("beamAxis_") + tag + "_"; + + std::cout << "\n============ Beam-axis fit (track crossings) ============\n"; + std::cout << "file: " << filename << "\n"; + std::cout << "tag: " << tag << " (histogram prefix: " << base << ")\n"; + std::cout << "range: " << zmin << " <= z <= " << zmax << " mm\n"; + std::cout << "z0: " << z0 << " mm (BEAM_AXIS_Z0)\n\n"; + + double bx = 0, tx = 0, bxe = 0, txe = 0, chi2x = -1; + double by = 0, ty = 0, bye = 0, tye = 0, chi2y = -1; + int nx = 0, ny = 0; + + std::cout << "--- X crossings ---\n"; + bool okx = fitOneProjection(f, base, "X", zmin, zmax, z0, bx, tx, bxe, txe, chi2x, nx); + std::cout << "\n--- Y crossings ---\n"; + bool oky = fitOneProjection(f, base, "Y", zmin, zmax, z0, by, ty, bye, tye, chi2y, ny); + + std::cout << "\n--- fitted beam line ---\n"; + if (okx) + std::cout << " x(z0) = " << bx << " +/- " << bxe << " mm, dx/dz = " << tx << " +/- " << txe + << " (" << tx * 1000.0 << " mrad, chi2/ndf = " << chi2x << ", " << nx << " planes)\n"; + if (oky) + std::cout << " y(z0) = " << by << " +/- " << bye << " mm, dy/dz = " << ty << " +/- " << tye + << " (" << ty * 1000.0 << " mrad, chi2/ndf = " << chi2y << ", " << ny << " planes)\n"; + + if (!okx || !oky) + { + std::cout << "\nOne or both projections failed -- not emitting exports.\n"; + std::cout << "=========================================================\n\n"; + f->Close(); + return; + } + + std::cout << "\n--- set these directly (ABSOLUTE values, do not add to current) ---\n"; + std::cout << "export BEAM_AXIS_X=" << bx << "\n"; + std::cout << "export BEAM_AXIS_Y=" << by << "\n"; + std::cout << "export BEAM_AXIS_Z0=" << z0 << "\n"; + std::cout << "export BEAM_TILT_X=" << tx << "\n"; + std::cout << "export BEAM_TILT_Y=" << ty << "\n"; + std::cout << "=========================================================\n\n"; + + f->Close(); +} \ No newline at end of file