new file: scratch/fit_a1c0_zcorr.C not sure if this is the right apporach, doing this for now as a stopgap
This commit is contained in:
parent
8c10235be9
commit
18b948a335
22
CLAUDE.md
22
CLAUDE.md
|
|
@ -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.
|
||||
188
scratch/fit_a1c0_zcorr.C
Normal file
188
scratch/fit_a1c0_zcorr.C
Normal file
|
|
@ -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 <TFile.h>
|
||||
#include <TH2.h>
|
||||
#include <TProfile.h>
|
||||
#include <TF1.h>
|
||||
#include <TKey.h>
|
||||
#include <TDirectory.h>
|
||||
#include <iostream>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
// 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<TH2 *>(o);
|
||||
|
||||
TIter next(dir->GetListOfKeys());
|
||||
while (TKey *key = static_cast<TKey *>(next()))
|
||||
{
|
||||
TObject *obj = key->ReadObj();
|
||||
if (obj->InheritsFrom(TDirectory::Class()))
|
||||
{
|
||||
if (TH2 *h = findTH2(static_cast<TDirectory *>(obj), name))
|
||||
return h;
|
||||
}
|
||||
else if (obj->InheritsFrom(TH2::Class()) && name == obj->GetName())
|
||||
{
|
||||
return static_cast<TH2 *>(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<std::string> &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();
|
||||
}
|
||||
148
scratch/fit_xy_offset.C
Normal file
148
scratch/fit_xy_offset.C
Normal file
|
|
@ -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 <TFile.h>
|
||||
#include <TH2.h>
|
||||
#include <TProfile.h>
|
||||
#include <TF1.h>
|
||||
#include <TKey.h>
|
||||
#include <TDirectory.h>
|
||||
#include <TMath.h>
|
||||
#include <iostream>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
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<TH2 *>(o);
|
||||
TIter next(dir->GetListOfKeys());
|
||||
while (TKey *key = static_cast<TKey *>(next()))
|
||||
{
|
||||
TObject *obj = key->ReadObj();
|
||||
if (obj->InheritsFrom(TDirectory::Class()))
|
||||
{
|
||||
if (TH2 *h = findTH2(static_cast<TDirectory *>(obj), name))
|
||||
return h;
|
||||
}
|
||||
else if (obj->InheritsFrom(TH2::Class()) && name == obj->GetName())
|
||||
return static_cast<TH2 *>(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();
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user