modified: TrackRecon.C consolidated the a1c1 linear functions and fits
new file: fit_a1c1_cfrac.C modified: run_tr.sh
This commit is contained in:
parent
9e2023a9db
commit
ff21920bae
460
TrackRecon.C
460
TrackRecon.C
|
|
@ -58,25 +58,76 @@ bool reactiondata = false;
|
||||||
|
|
||||||
TF1 pcfix_func("func", model_invert, -200, 200);
|
TF1 pcfix_func("func", model_invert, -200, 200);
|
||||||
|
|
||||||
// --- A1C1 parabola calibration (cfrac = c0[cell] - C*s^2), gain-selected ---
|
// --- A1C1 linear centre-fold model: cfrac = cfmin[cell] + k[cell]*fold ----------
|
||||||
const double a1c1_C_1 = 0.22;
|
// fold = |z - z_cellcentre| / halfcell in [0,1]. The raw cfrac pedestal floats
|
||||||
const double a1c1_c0_1[7] = {0.46, 0.46, 0.46, 0.46, 0.45, 0.45, 0.45};
|
// with the arbitrary anode/cathode gain, so cfmin/k are dataset-specific: fit
|
||||||
const double a1c1_C_2 = 0.20;
|
// them OFFLINE with fit_a1c1_cfrac.C on prebuilt data and paste the per-cell
|
||||||
const double a1c1_c0_2[7] = {0.24, 0.24, 0.24, 0.24, 0.24, 0.24, 0.24};
|
// results below; Begin() selects the active set by dataset.
|
||||||
|
const double a1c1_zg[8] = {147.998, 101.946, 59.7634, 19.6965, -19.6965, -59.7634, -101.946, -147.998};
|
||||||
|
|
||||||
double a1c1_C = a1c1_C_1;
|
static const double a1c1_cfmin_17F[7] = {0.557612, 0.400000, 0.400000, 0.342771, 0.400000, 0.682479, 0.400000};
|
||||||
double a1c1_c0[7] = {0.46, 0.46, 0.46, 0.46, 0.45, 0.45, 0.45};
|
static const double a1c1_k_17F[7] = {0.0688469, 0.075000, 0.075000, 0.160131, 0.075000, 0.273423, 0.075000};
|
||||||
double a1c1_slope = 1.0; // parabola staircase stretch (1.0 = off)
|
static const double a1c1_cfmin_27Al[7] = {0.18, 0.18, 0.18, 0.18, 0.18, 0.18, 0.18}; // TODO: optimise on 27Al data
|
||||||
|
static const double a1c1_k_27Al[7] = {0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06};
|
||||||
|
|
||||||
// --- A1C1 linear centre-fold calibration (cfrac = cfmin[cell] + k*f) ---
|
// active per-cell set, populated by dataset in Begin()
|
||||||
const double a1c1_k_1 = 0.075; // from cfrac_vs_fold candle medians
|
double a1c1_cfmin_cell[7] = {0.40, 0.40, 0.40, 0.40, 0.40, 0.40, 0.40};
|
||||||
const double a1c1_cfmin_1[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};
|
||||||
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;
|
// --- A1C1 LOW BAND (incomplete charge integration) -----------------------------
|
||||||
double a1c1_cfmin[7] = {0.40, 0.40, 0.40, 0.40, 0.40, 0.40, 0.40};
|
// 17F shows a second, parallel cfrac band (~0.10 -> 0.15 over the fold) where the
|
||||||
bool a1c1_use_linear = true; // linear vs parabola toggle
|
// cathode charge is only partially integrated. It still tracks z (real physics),
|
||||||
|
// so it gets its own per-cell cfmin/k. Events with cfrac below a1c1_cfrac_split
|
||||||
|
// are reconstructed with this set instead of being rejected. a1c1_cfrac_split<=0
|
||||||
|
// disables the low band (e.g. 27Al, which shows no second band).
|
||||||
|
static const double a1c1_cfmin2_17F[7] = {0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10};
|
||||||
|
static const double a1c1_k2_17F[7] = {0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05};
|
||||||
|
static const double a1c1_cfmin2_27Al[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; // no low band
|
||||||
|
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
|
||||||
|
|
||||||
|
// 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
|
||||||
|
// fired wire, folds about the adjacent cell centre, inverts cfrac = cfmin+k*fold
|
||||||
|
// per cell. Picks the LOW band (band=1) when cfrac < a1c1_cfrac_split, else the
|
||||||
|
// MAIN band (band=0). Callers apply their own acceptance via the returned flags:
|
||||||
|
// inband : f in [0,1] (cfrac within the calibrated band)
|
||||||
|
// pitchok : |pcz - zf| <= pitch (consistent with the fired wire)
|
||||||
|
struct A1C1Sol { double pcz; double f; int cell; double pitch; bool inband; bool pitchok; int band; };
|
||||||
|
inline A1C1Sol a1c1_solve(double cfrac, double zf, double z_a1c0)
|
||||||
|
{
|
||||||
|
A1C1Sol s{zf, 0.0, 0, 0.0, false, false, 0};
|
||||||
|
// band select: incomplete-integration low band uses its own per-cell cfmin/k
|
||||||
|
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)
|
||||||
|
{
|
||||||
|
cfmin = a1c1_cfmin2_cell;
|
||||||
|
kk = a1c1_k2_cell;
|
||||||
|
s.band = 1;
|
||||||
|
}
|
||||||
|
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;
|
||||||
|
int cell = (z_a1c0 >= a1c1_zg[wf]) ? (wf - 1) : wf;
|
||||||
|
if (cell < 0) cell = 0;
|
||||||
|
if (cell > 6) cell = 6;
|
||||||
|
s.cell = cell;
|
||||||
|
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;
|
||||||
|
s.f = (cfrac - cfmin[cell]) / kk[cell];
|
||||||
|
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);
|
||||||
|
s.pitchok = (TMath::Abs(s.pcz - zf) <= s.pitch);
|
||||||
|
return s;
|
||||||
|
}
|
||||||
|
|
||||||
TGraph *MeV_to_cm = NULL, *cm_to_MeV = NULL;
|
TGraph *MeV_to_cm = NULL, *cm_to_MeV = NULL;
|
||||||
TGraph *MeV_to_cm_p = NULL, *cm_to_MeVp = NULL;
|
TGraph *MeV_to_cm_p = NULL, *cm_to_MeVp = NULL;
|
||||||
|
|
@ -210,27 +261,33 @@ void TrackRecon::Begin(TTree * /*tree*/)
|
||||||
|
|
||||||
if (getenv("Gain"))
|
if (getenv("Gain"))
|
||||||
Gain = std::atof(getenv("Gain"));
|
Gain = std::atof(getenv("Gain"));
|
||||||
bool highP = (Gain >= 2);
|
|
||||||
|
|
||||||
a1c1_C = highP ? a1c1_C_2 : a1c1_C_1; // parabola set
|
// --- A1C1 per-cell linear centre-fold constants: select the static, offline-
|
||||||
|
// optimised set for this dataset (cfmin floats with the arbitrary gain, so the
|
||||||
|
// values are dataset-specific). Re-fit them with fit_a1c1_cfrac.C on prebuilt
|
||||||
|
// data and paste into the a1c1_{cfmin,k}_<dataset> arrays above. ---
|
||||||
|
const double *cfmin_src = a1c1_cfmin_17F;
|
||||||
|
const double *k_src = a1c1_k_17F;
|
||||||
|
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
|
||||||
|
if (dataset == "27Al")
|
||||||
|
{
|
||||||
|
cfmin_src = a1c1_cfmin_27Al;
|
||||||
|
k_src = a1c1_k_27Al;
|
||||||
|
cfmin2_src = a1c1_cfmin2_27Al;
|
||||||
|
k2_src = a1c1_k2_27Al;
|
||||||
|
a1c1_cfrac_split = 0.0; // 27Al: no second band, low band disabled
|
||||||
|
}
|
||||||
for (int i = 0; i < 7; ++i)
|
for (int i = 0; i < 7; ++i)
|
||||||
a1c1_c0[i] = highP ? a1c1_c0_2[i] : a1c1_c0_1[i];
|
{
|
||||||
if (getenv("a1c1_C"))
|
a1c1_cfmin_cell[i] = cfmin_src[i];
|
||||||
a1c1_C = std::atof(getenv("a1c1_C"));
|
a1c1_k_cell[i] = k_src[i];
|
||||||
if (getenv("a1c1_slope"))
|
a1c1_cfmin2_cell[i] = cfmin2_src[i];
|
||||||
a1c1_slope = std::atof(getenv("a1c1_slope"));
|
a1c1_k2_cell[i] = k2_src[i];
|
||||||
|
}
|
||||||
a1c1_k = highP ? a1c1_k_2 : a1c1_k_1; // linear set
|
std::cout << "A1C1 per-cell constants: using static " << (dataset.empty() ? "(default 17F)" : dataset)
|
||||||
for (int i = 0; i < 7; ++i)
|
<< " set; low-band split cfrac<" << a1c1_cfrac_split << std::endl;
|
||||||
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();
|
pwinstance.ConstructGeo();
|
||||||
|
|
||||||
|
|
@ -771,7 +828,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
|
||||||
{
|
{
|
||||||
cathodeT = static_cast<double>(pc.t[i]);
|
cathodeT = static_cast<double>(pc.t[i]);
|
||||||
cathodeIndex = pc.index[i] - 24;
|
cathodeIndex = pc.index[i] - 24;
|
||||||
cWireEvents[pc.index[i] - 24] = std::tuple(pc.index[i] - 24, 4*pc.e[i], static_cast<double>(pc.t[i]));
|
cWireEvents[pc.index[i] - 24] = std::tuple(pc.index[i] - 24, pc.e[i], static_cast<double>(pc.t[i]));
|
||||||
}
|
}
|
||||||
|
|
||||||
if (anodeT != -99999 && cathodeT != 99999)
|
if (anodeT != -99999 && cathodeT != 99999)
|
||||||
|
|
@ -1604,11 +1661,17 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
|
||||||
|
|
||||||
// --- A1C1 charge fraction (single max-E cathode vs anode, pseudo-wire sums) ---
|
// --- A1C1 charge fraction (single max-E cathode vs anode, pseudo-wire sums) ---
|
||||||
double aSumE_bm = std::get<1>(pw_tuple);
|
double aSumE_bm = std::get<1>(pw_tuple);
|
||||||
auto cpw_tuple = pwinstance.GetPseudoWire(cOne, "CATHODE");
|
double cSumE_bm = std::get<1>(cMaxWire); // Extract the energy directly!
|
||||||
double cSumE_bm = std::get<1>(cpw_tuple);
|
|
||||||
double ac_sum = aSumE_bm + cSumE_bm;
|
double ac_sum = aSumE_bm + cSumE_bm;
|
||||||
double cfrac = (ac_sum > 0.0) ? cSumE_bm / ac_sum : -1.0;
|
double cfrac = (ac_sum > 0.0) ? cSumE_bm / ac_sum : -1.0;
|
||||||
|
|
||||||
|
// Cmax/anode vs Si phi for ALL events reaching this block (not just the
|
||||||
|
// A1C2 subset). Cmax = single max-E cathode wire (cSumE_bm), anode =
|
||||||
|
// anode pseudo-wire sum. Ratio is unbounded above (arbitrary cathode gains).
|
||||||
|
if (aSumE_bm > 0.0)
|
||||||
|
plotter->Fill2D("Benchmark_SX3_CmaxOverAnode_vs_phi", 180, -180, 180, 250, 0, 5,
|
||||||
|
sx3event.pos.Phi() * 180. / M_PI, cSumE_bm / aSumE_bm, "Benchmark_SX3_ref");
|
||||||
|
|
||||||
// Smearing setup
|
// Smearing setup
|
||||||
double sx3_phi_pitch = 6.5 * (M_PI / 180.0);
|
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);
|
double smeared_phi = sx3event.pos.Phi() + rand.Uniform(-sx3_phi_pitch / 2.0, sx3_phi_pitch / 2.0);
|
||||||
|
|
@ -1636,26 +1699,25 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
|
||||||
fillVsRef(tag, pcz, vtx, pcz_ref, vtx_ref);
|
fillVsRef(tag, pcz, vtx, pcz_ref, vtx_ref);
|
||||||
};
|
};
|
||||||
|
|
||||||
// --- A1C1 with the cfrac sub-cell model (toggle linear/parabola) ---
|
// --- A1C1 with the cfrac sub-cell model (linear centre-fold) ---
|
||||||
// pivot = xo_a1c1.Z() (fired max-E cathode crossover). LINEAR: side from
|
// Anchored on the fired cathode (xo_a1c1.Z()); cell from the anode z;
|
||||||
// the fired cathode, |d| from cfrac about the cell centre. PARABOLA: side
|
// |offset| from cfrac. REJECT (no fill) when out of band or inconsistent.
|
||||||
// from the Si guess about the pivot. Both fall back to pczguess.
|
|
||||||
auto doA1C1Model = [&](const std::string &tag, const TVector3 &si_point)
|
auto doA1C1Model = [&](const std::string &tag, const TVector3 &si_point)
|
||||||
{
|
{
|
||||||
if (!a1c1Good || cfrac < 0.0)
|
if (!a1c1Good || cfrac < 0.0)
|
||||||
return;
|
return;
|
||||||
bool ok = false;
|
double z_a1c0 = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, si_point.Phi()).Z();
|
||||||
double pcz = a1c1_use_linear
|
A1C1Sol s = a1c1_solve(cfrac, xo_a1c1.Z(), z_a1c0);
|
||||||
? model_invert_a1c1_linear(cfrac, pczguess, xo_a1c1.Z(), a1c1_k, a1c1_cfmin, ok)
|
if (!(s.inband && s.pitchok))
|
||||||
: model_invert_a1c1(cfrac, pczguess, xo_a1c1.Z(), a1c1_C, a1c1_c0, ok, a1c1_slope);
|
return;
|
||||||
TVector3 vtx = vertexFrom(si_point, TVector3(xo_a1c1.X(), xo_a1c1.Y(), pcz));
|
TVector3 vtx = vertexFrom(si_point, TVector3(xo_a1c1.X(), xo_a1c1.Y(), s.pcz));
|
||||||
fillSuite(tag, pcz, vtx);
|
fillSuite(tag, s.pcz, vtx);
|
||||||
fillVsRef(tag, pcz, vtx, pcz_ref, vtx_ref);
|
fillVsRef(tag, s.pcz, vtx, pcz_ref, vtx_ref);
|
||||||
};
|
};
|
||||||
|
|
||||||
if (phicut && PCSX3TimeCut)
|
if (phicut && PCSX3TimeCut)
|
||||||
{
|
{
|
||||||
if (pcevent.multi1 == 1 && pcevent.multi2 >= 1)
|
if (pcevent.multi1 == 1 && pcevent.multi2 == 2)
|
||||||
{
|
{
|
||||||
fillSuite("A1C2", pcz_ref, vtx_ref);
|
fillSuite("A1C2", pcz_ref, vtx_ref);
|
||||||
doA1C1("A1C1", sx3event.pos, true, false);
|
doA1C1("A1C1", sx3event.pos, true, false);
|
||||||
|
|
@ -1676,17 +1738,12 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
|
||||||
pcz_ref, cfrac, "Benchmark_SX3_ref");
|
pcz_ref, cfrac, "Benchmark_SX3_ref");
|
||||||
plotter->Fill2D("Benchmark_SX3_A1C1_cfrac_vs_sx3pczguess", 400, -200, 200, 220, -0.05, 1.05,
|
plotter->Fill2D("Benchmark_SX3_A1C1_cfrac_vs_sx3pczguess", 400, -200, 200, 220, -0.05, 1.05,
|
||||||
pczguess, cfrac, "Benchmark_SX3_ref");
|
pczguess, cfrac, "Benchmark_SX3_ref");
|
||||||
if (pcevent.multi2 == 1)
|
|
||||||
{
|
|
||||||
plotter->Fill2D("Benchmark_SX3_A1C1_cfrac_vs_sx3pczguess_a1c1", 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};
|
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();
|
double zp = xo_a1c1.Z();
|
||||||
auto fillCfracS = [&](const char *name, double val)
|
auto fillCfracS = [&](const char *name, double truth)
|
||||||
{
|
{
|
||||||
double sgn = (val >= zp) ? 1.0 : -1.0;
|
double sgn = (truth >= zp) ? 1.0 : -1.0;
|
||||||
double znb = (sgn > 0) ? 1.0e30 : -1.0e30;
|
double znb = (sgn > 0) ? 1.0e30 : -1.0e30;
|
||||||
for (int i = 0; i < 8; ++i)
|
for (int i = 0; i < 8; ++i)
|
||||||
{
|
{
|
||||||
|
|
@ -1697,7 +1754,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
|
||||||
}
|
}
|
||||||
if (TMath::Abs(znb) < 1e8 && TMath::Abs(znb - zp) > 0.0)
|
if (TMath::Abs(znb) < 1e8 && TMath::Abs(znb - zp) > 0.0)
|
||||||
plotter->Fill2D(name, 240, -1.2, 1.2, 220, -0.05, 1.05,
|
plotter->Fill2D(name, 240, -1.2, 1.2, 220, -0.05, 1.05,
|
||||||
(val - zp) / TMath::Abs(znb - zp), cfrac, "Benchmark_SX3_ref");
|
(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", pcz_ref);
|
||||||
fillCfracS("Benchmark_SX3_A1C1_cfrac_vs_s_sx3pczguess", pczguess);
|
fillCfracS("Benchmark_SX3_A1C1_cfrac_vs_s_sx3pczguess", pczguess);
|
||||||
|
|
@ -1715,6 +1772,81 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
|
||||||
break;
|
break;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// ---- TRUE A1C1: genuine one-anode/one-cathode events. There is no
|
||||||
|
// A1C2 crossover reference for these, so compare only to the Si
|
||||||
|
// geometric guess (pczguess). xo_a1c1 is the real measured crossover. ----
|
||||||
|
else if (pcevent.multi1 == 1 && pcevent.multi2 == 1 && a1c1Good)
|
||||||
|
{
|
||||||
|
double pcz_raw = xo_a1c1.Z();
|
||||||
|
TVector3 vtx_raw = vertexFrom(sx3event.pos, TVector3(xo_a1c1.X(), xo_a1c1.Y(), pcz_raw));
|
||||||
|
fillSuite("trueA1C1", pcz_raw, vtx_raw);
|
||||||
|
plotter->Fill2D("Benchmark_SX3_PCZ_trueA1C1_vs_sx3pczguess", 400, -200, 200, 400, -200, 200, pczguess, pcz_raw, "Benchmark_SX3_trueA1C1");
|
||||||
|
plotter->Fill1D("Benchmark_SX3_PCZ_trueA1C1_minus_sx3pczguess", 400, -100, 100, pcz_raw - pczguess, "Benchmark_SX3_trueA1C1");
|
||||||
|
|
||||||
|
// cfrac sub-cell reconstruction, anchored on the FIRED CATHODE.
|
||||||
|
// Vertex-independent (no pczguess): the cell is the one adjacent to
|
||||||
|
// the cathode that fired (zf), the side is taken from the anode z,
|
||||||
|
// and the hit must land within one cell pitch of zf to be valid.
|
||||||
|
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);
|
||||||
|
int cell = s.cell;
|
||||||
|
double f = s.f;
|
||||||
|
double pcz_cf = s.pcz;
|
||||||
|
bool valid = s.pitchok; // accept if consistent with the fired wire
|
||||||
|
|
||||||
|
TVector3 vtx_cf = vertexFrom(sx3event.pos, TVector3(xo_a1c1.X(), xo_a1c1.Y(), pcz_cf));
|
||||||
|
fillSuite(valid ? "trueA1C1_Cfrac" : "trueA1C1_Cfrac_invalid", pcz_cf, vtx_cf);
|
||||||
|
plotter->Fill1D("Benchmark_SX3_trueA1C1_cfrac", 220, -0.05, 1.05, cfrac, "Benchmark_SX3_trueA1C1");
|
||||||
|
// reference-free per-cell cfrac (cell from geometry, no A1C2 ref): the
|
||||||
|
// offline fitter reads per-cell edges/percentiles to gain-match cfmin/k.
|
||||||
|
plotter->Fill2D("Benchmark_SX3_trueA1C1_cfrac_vs_cell", 7, 0, 7, 220, -0.05, 1.05, cell + 0.5, cfrac, "Benchmark_SX3_trueA1C1");
|
||||||
|
plotter->Fill1D("Benchmark_SX3_trueA1C1_f", 260, -1.5, 2.5, f, "Benchmark_SX3_trueA1C1");
|
||||||
|
plotter->Fill1D("Benchmark_SX3_trueA1C1_valid", 2, 0, 2, valid ? 1.0 : 0.0, "Benchmark_SX3_trueA1C1");
|
||||||
|
// failure-reason breakdown (why the cfrac estimate is / isn't usable):
|
||||||
|
// 0 valid & f in [0,1] ideal, inside the calibrated band
|
||||||
|
// 1 valid & f in [-1,0) below band but within the pitch tolerance
|
||||||
|
// 2 valid & f in (1,3] above band but within the pitch tolerance
|
||||||
|
// 3 invalid: f < -1 cfrac far below band (cfmin too high / cathode low?)
|
||||||
|
// 4 invalid: f > 3 cfrac far above band (cfmin too low / cathode high?)
|
||||||
|
// 5 cell k<=0 cell not autocalibrated
|
||||||
|
int reason;
|
||||||
|
if (a1c1_k_cell[cell] <= 0.0) reason = 5;
|
||||||
|
else if (!valid) reason = (f < 0.0) ? 3 : 4;
|
||||||
|
else if (f < 0.0) reason = 1;
|
||||||
|
else if (f > 1.0) reason = 2;
|
||||||
|
else reason = 0;
|
||||||
|
plotter->Fill1D("Benchmark_SX3_trueA1C1_failreason", 6, 0, 6, reason + 0.5, "Benchmark_SX3_trueA1C1");
|
||||||
|
// valid-only breakdown (reason is 0/1/2 here): how clean the accepted
|
||||||
|
// sample is -- 0 ideal (f in [0,1]), 1 below-band, 2 above-band marginal.
|
||||||
|
if (valid)
|
||||||
|
plotter->Fill1D("Benchmark_SX3_trueA1C1_validreason", 3, 0, 3, reason + 0.5, "Benchmark_SX3_trueA1C1");
|
||||||
|
// which cfrac band reconstructed this event: 0 = main, 1 = low (incomplete integration)
|
||||||
|
plotter->Fill1D("Benchmark_SX3_trueA1C1_band", 2, 0, 2, s.band + 0.5, "Benchmark_SX3_trueA1C1");
|
||||||
|
if (valid)
|
||||||
|
plotter->Fill1D("Benchmark_SX3_trueA1C1_band_valid", 2, 0, 2, s.band + 0.5, "Benchmark_SX3_trueA1C1");
|
||||||
|
// source-run-only diagnostics: pczguess is meaningful only when the vertex is fixed
|
||||||
|
if (valid)
|
||||||
|
{
|
||||||
|
plotter->Fill1D("Benchmark_SX3_PCZ_trueA1C1_Cfrac_minus_sx3pczguess_DIAG", 400, -100, 100, pcz_cf - pczguess, "Benchmark_SX3_trueA1C1");
|
||||||
|
plotter->Fill2D("Benchmark_SX3_PCZ_trueA1C1_Cfrac_vs_sx3pczguess_DIAG", 400, -200, 200, 400, -200, 200, pczguess, pcz_cf, "Benchmark_SX3_trueA1C1");
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// anode-only (A1C0) estimate for the same events
|
||||||
|
{
|
||||||
|
TVector3 pc = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, sx3event.pos.Phi());
|
||||||
|
TVector3 vtx0 = vertexFrom(sx3event.pos, pc);
|
||||||
|
if (vtx0.Perp() <= 6.0 && vtx0.Z() >= -173.6)
|
||||||
|
{
|
||||||
|
fillSuite("trueA1C0", pc.Z(), vtx0);
|
||||||
|
plotter->Fill2D("Benchmark_SX3_PCZ_trueA1C0_vs_sx3pczguess", 400, -200, 200, 400, -200, 200, pczguess, pc.Z(), "Benchmark_SX3_trueA1C1");
|
||||||
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
@ -1905,6 +2037,13 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
|
||||||
double ac_sum = aSumE_bm + cSumE_bm;
|
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; // bounded [0,1]
|
||||||
|
|
||||||
|
// Cmax/anode vs Si phi for ALL events reaching this block (not just the
|
||||||
|
// A1C2 subset). Cmax = single max-E cathode wire (cSumE_bm), anode =
|
||||||
|
// anode pseudo-wire sum. Ratio is unbounded above (arbitrary cathode gains).
|
||||||
|
if (aSumE_bm > 0.0)
|
||||||
|
plotter->Fill2D("Benchmark_QQQ_CmaxOverAnode_vs_phi", 180, -180, 180, 250, 0, 5,
|
||||||
|
qqqevent.pos.Phi() * 180. / M_PI, cSumE_bm / aSumE_bm, "Benchmark_QQQ_ref");
|
||||||
|
|
||||||
// Smearing Setup
|
// Smearing Setup
|
||||||
double qqq_wedge_pitch = (87.0 / 16.0) * (M_PI / 180.0);
|
double qqq_wedge_pitch = (87.0 / 16.0) * (M_PI / 180.0);
|
||||||
double qqq_ring_pitch = 48.0 / 16.0;
|
double qqq_ring_pitch = 48.0 / 16.0;
|
||||||
|
|
@ -1936,25 +2075,21 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
|
||||||
fillVsRef(tag, pcz, vtx, pcz_ref, vtx_ref);
|
fillVsRef(tag, pcz, vtx, pcz_ref, vtx_ref);
|
||||||
};
|
};
|
||||||
|
|
||||||
// --- A1C1 with the cfrac sub-cell model (QQQ) ---
|
// --- A1C1 with the cfrac sub-cell model (QQQ, linear centre-fold) ---
|
||||||
// pivot = xo_a1c1.Z() (fired max-E cathode crossover). The active model
|
// Anchored on the fired cathode (xo_a1c1.Z()); cell from the anode-only
|
||||||
// is chosen by a1c1_use_linear: LINEAR centre-fold (side from the fired
|
// z (z_a1c0, independent of the Si guess); |offset| from cfrac. REJECT
|
||||||
// cathode, |d| from cfrac, scaled by the local cell) or the PARABOLA
|
// when out of band or inconsistent with the fired wire.
|
||||||
// (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)
|
auto doA1C1Model = [&](const std::string &tag, const TVector3 &si_point)
|
||||||
{
|
{
|
||||||
if (!a1c1Good || cfrac < 0.0)
|
if (!a1c1Good || cfrac < 0.0)
|
||||||
return;
|
return;
|
||||||
bool ok = false;
|
|
||||||
// anchor on the A1C0 anode-only z (independent of the Si geometric guess) to avoid circularity
|
|
||||||
double z_a1c0 = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, si_point.Phi()).Z();
|
double z_a1c0 = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, si_point.Phi()).Z();
|
||||||
double pcz = a1c1_use_linear
|
A1C1Sol s = a1c1_solve(cfrac, xo_a1c1.Z(), z_a1c0);
|
||||||
? model_invert_a1c1_linear(cfrac, z_a1c0, xo_a1c1.Z(), a1c1_k, a1c1_cfmin, ok)
|
if (!(s.inband && s.pitchok))
|
||||||
: model_invert_a1c1(cfrac, z_a1c0, xo_a1c1.Z(), a1c1_C, a1c1_c0, ok, a1c1_slope);
|
return;
|
||||||
TVector3 vtx = vertexFrom(si_point, TVector3(xo_a1c1.X(), xo_a1c1.Y(), pcz));
|
TVector3 vtx = vertexFrom(si_point, TVector3(xo_a1c1.X(), xo_a1c1.Y(), s.pcz));
|
||||||
fillSuite(tag, pcz, vtx);
|
fillSuite(tag, s.pcz, vtx);
|
||||||
fillVsRef(tag, pcz, vtx, pcz_ref, vtx_ref);
|
fillVsRef(tag, s.pcz, vtx, pcz_ref, vtx_ref);
|
||||||
};
|
};
|
||||||
|
|
||||||
if (phicut && timecut)
|
if (phicut && timecut)
|
||||||
|
|
@ -2015,6 +2150,81 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
|
||||||
break;
|
break;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// ---- TRUE A1C1 (QQQ): genuine one-anode/one-cathode events. No A1C2
|
||||||
|
// crossover reference exists, so compare only to the geometric guess
|
||||||
|
// (pcz_guess_int). xo_a1c1 is the real measured crossover. ----
|
||||||
|
else if (pcevent.multi1 >= 1 && pcevent.multi2 == 1 && a1c1Good)
|
||||||
|
{
|
||||||
|
double pcz_raw = xo_a1c1.Z();
|
||||||
|
TVector3 vtx_raw = vertexFrom(qqqevent.pos, TVector3(xo_a1c1.X(), xo_a1c1.Y(), pcz_raw));
|
||||||
|
fillSuite("trueA1C1", pcz_raw, vtx_raw);
|
||||||
|
plotter->Fill2D("Benchmark_QQQ_PCZ_trueA1C1_vs_qqqpczguess", 400, -200, 200, 400, -200, 200, pcz_guess_int, pcz_raw, "Benchmark_QQQ_trueA1C1");
|
||||||
|
plotter->Fill1D("Benchmark_QQQ_PCZ_trueA1C1_minus_qqqpczguess", 400, -100, 100, pcz_raw - pcz_guess_int, "Benchmark_QQQ_trueA1C1");
|
||||||
|
|
||||||
|
// cfrac sub-cell reconstruction, anchored on the FIRED CATHODE.
|
||||||
|
// Vertex-independent (no pcz_guess_int): the cell is the one adjacent
|
||||||
|
// to the cathode that fired (zf), the side is taken from the anode z,
|
||||||
|
// and the hit must land within one cell pitch of zf to be valid.
|
||||||
|
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);
|
||||||
|
int cell = s.cell;
|
||||||
|
double f = s.f;
|
||||||
|
double pcz_cf = s.pcz;
|
||||||
|
bool valid = s.pitchok; // accept if consistent with the fired wire
|
||||||
|
|
||||||
|
TVector3 vtx_cf = vertexFrom(qqqevent.pos, TVector3(xo_a1c1.X(), xo_a1c1.Y(), pcz_cf));
|
||||||
|
fillSuite(valid ? "trueA1C1_Cfrac" : "trueA1C1_Cfrac_invalid", pcz_cf, vtx_cf);
|
||||||
|
plotter->Fill1D("Benchmark_QQQ_trueA1C1_cfrac", 220, -0.05, 1.05, cfrac, "Benchmark_QQQ_trueA1C1");
|
||||||
|
// reference-free per-cell cfrac (cell from geometry, no A1C2 ref): the
|
||||||
|
// offline fitter reads per-cell edges/percentiles to gain-match cfmin/k.
|
||||||
|
plotter->Fill2D("Benchmark_QQQ_trueA1C1_cfrac_vs_cell", 7, 0, 7, 220, -0.05, 1.05, cell + 0.5, cfrac, "Benchmark_QQQ_trueA1C1");
|
||||||
|
plotter->Fill1D("Benchmark_QQQ_trueA1C1_f", 260, -1.5, 2.5, f, "Benchmark_QQQ_trueA1C1");
|
||||||
|
plotter->Fill1D("Benchmark_QQQ_trueA1C1_valid", 2, 0, 2, valid ? 1.0 : 0.0, "Benchmark_QQQ_trueA1C1");
|
||||||
|
// failure-reason breakdown (why the cfrac estimate is / isn't usable):
|
||||||
|
// 0 valid & f in [0,1] ideal, inside the calibrated band
|
||||||
|
// 1 valid & f in [-1,0) below band but within the pitch tolerance
|
||||||
|
// 2 valid & f in (1,3] above band but within the pitch tolerance
|
||||||
|
// 3 invalid: f < -1 cfrac far below band (cfmin too high / cathode low?)
|
||||||
|
// 4 invalid: f > 3 cfrac far above band (cfmin too low / cathode high?)
|
||||||
|
// 5 cell k<=0 cell not autocalibrated
|
||||||
|
int reason;
|
||||||
|
if (a1c1_k_cell[cell] <= 0.0) reason = 5;
|
||||||
|
else if (!valid) reason = (f < 0.0) ? 3 : 4;
|
||||||
|
else if (f < 0.0) reason = 1;
|
||||||
|
else if (f > 1.0) reason = 2;
|
||||||
|
else reason = 0;
|
||||||
|
plotter->Fill1D("Benchmark_QQQ_trueA1C1_failreason", 6, 0, 6, reason + 0.5, "Benchmark_QQQ_trueA1C1");
|
||||||
|
// valid-only breakdown (reason is 0/1/2 here): how clean the accepted
|
||||||
|
// sample is -- 0 ideal (f in [0,1]), 1 below-band, 2 above-band marginal.
|
||||||
|
if (valid)
|
||||||
|
plotter->Fill1D("Benchmark_QQQ_trueA1C1_validreason", 3, 0, 3, reason + 0.5, "Benchmark_QQQ_trueA1C1");
|
||||||
|
// which cfrac band reconstructed this event: 0 = main, 1 = low (incomplete integration)
|
||||||
|
plotter->Fill1D("Benchmark_QQQ_trueA1C1_band", 2, 0, 2, s.band + 0.5, "Benchmark_QQQ_trueA1C1");
|
||||||
|
if (valid)
|
||||||
|
plotter->Fill1D("Benchmark_QQQ_trueA1C1_band_valid", 2, 0, 2, s.band + 0.5, "Benchmark_QQQ_trueA1C1");
|
||||||
|
// source-run-only diagnostics: pcz_guess_int is meaningful only when the vertex is fixed
|
||||||
|
if (valid)
|
||||||
|
{
|
||||||
|
plotter->Fill1D("Benchmark_QQQ_PCZ_trueA1C1_Cfrac_minus_qqqpczguess_DIAG", 400, -100, 100, pcz_cf - pcz_guess_int, "Benchmark_QQQ_trueA1C1");
|
||||||
|
plotter->Fill2D("Benchmark_QQQ_PCZ_trueA1C1_Cfrac_vs_qqqpczguess_DIAG", 400, -200, 200, 400, -200, 200, pcz_guess_int, pcz_cf, "Benchmark_QQQ_trueA1C1");
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// anode-only (A1C0) estimate for the same events
|
||||||
|
{
|
||||||
|
TVector3 pc = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, qqqevent.pos.Phi());
|
||||||
|
TVector3 vtx0 = vertexFrom(qqqevent.pos, pc);
|
||||||
|
if (vtx0.Perp() <= 6.0 && vtx0.Z() >= -173.6)
|
||||||
|
{
|
||||||
|
fillSuite("trueA1C0", pc.Z(), vtx0);
|
||||||
|
plotter->Fill2D("Benchmark_QQQ_PCZ_trueA1C0_vs_qqqpczguess", 400, -200, 200, 400, -200, 200, pcz_guess_int, pc.Z(), "Benchmark_QQQ_trueA1C1");
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
@ -2591,6 +2801,56 @@ void protonMiscHistograms(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
|
||||||
pcz_fix = rand.Gaus(pcevent.pos.Z(), 8.0); // dither for a1c1 events
|
pcz_fix = rand.Gaus(pcevent.pos.Z(), 8.0); // dither for a1c1 events
|
||||||
pcz_dith = pcz_fix;
|
pcz_dith = pcz_fix;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
// --- A1C1 dither vs cfrac comparison (keep BOTH methods) ---
|
||||||
|
// For genuine A1C1 events, reconstruct the PC z two ways and fill a parallel
|
||||||
|
// set of excitation observables (suffix _dither / _cfrac) so the methods can
|
||||||
|
// be overlaid directly. This runs independently of the main-flow vertex cut
|
||||||
|
// below; each method applies its own vertex gate. The dither value also
|
||||||
|
// feeds the main (untagged) histograms exactly as before. cfrac =
|
||||||
|
// cpMax/(apSum+cpMax) = Energy2/(Energy1+Energy2): cell anchored on the fired
|
||||||
|
// cathode, side from the anode-only z, offset from the per-cell autocal.
|
||||||
|
if (pcevent.multi2 == 1 && pcevent.Energy2 > 1400)
|
||||||
|
{
|
||||||
|
auto fillCmp = [&](double pcz, const std::string &m)
|
||||||
|
{
|
||||||
|
TVector3 x2(pcevent.pos.X(), pcevent.pos.Y(), pcz);
|
||||||
|
TVector3 vv = x2 - qqqevent.pos;
|
||||||
|
double tm = -1.0 * (qqqevent.pos.X() * vv.X() + qqqevent.pos.Y() * vv.Y()) / (vv.X() * vv.X() + vv.Y() * vv.Y());
|
||||||
|
TVector3 rv = qqqevent.pos + tm * vv;
|
||||||
|
if (rv.Perp() > 6.0 || rv.Z() < -173.6 || rv.Z() > 100)
|
||||||
|
return;
|
||||||
|
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);
|
||||||
|
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);
|
||||||
|
};
|
||||||
|
|
||||||
|
fillCmp(pcz_dith, "dither"); // method 1: Gaussian dither (main-flow value)
|
||||||
|
|
||||||
|
// method 2: cfrac sub-cell linear centre-fold (side ref = anode-only z
|
||||||
|
// rebuilt from the fired anode wire)
|
||||||
|
double ac = pcevent.Energy1 + pcevent.Energy2;
|
||||||
|
double cfrac = (ac > 0.0) ? pcevent.Energy2 / ac : -1.0;
|
||||||
|
if (cfrac >= 0.0)
|
||||||
|
{
|
||||||
|
std::vector<std::tuple<int, double, double>> 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);
|
||||||
|
if (s.inband && s.pitchok)
|
||||||
|
{
|
||||||
|
fillCmp(s.pcz, "cfrac");
|
||||||
|
plotter->Fill2D("pmisc_a1c1cmp_pcz_cfrac_vs_dither", 600, -300, 300, 600, -300, 300, pcz_dith, s.pcz, "proton+misc_a1c1cmp");
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
TVector3 x2f(pcevent.pos.X(), pcevent.pos.Y(), pcz_fix);
|
TVector3 x2f(pcevent.pos.X(), pcevent.pos.Y(), pcz_fix);
|
||||||
TVector3 x1(qqqevent.pos);
|
TVector3 x1(qqqevent.pos);
|
||||||
TVector3 v = x2f - x1;
|
TVector3 v = x2f - x1;
|
||||||
|
|
@ -2676,6 +2936,8 @@ void protonMiscHistograms(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
|
||||||
void protonMiscHistograms_sx3(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events)
|
void protonMiscHistograms_sx3(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events)
|
||||||
{
|
{
|
||||||
// consider the 'proton-like' QQQ branch seen in a,p data
|
// consider the 'proton-like' QQQ branch seen in a,p data
|
||||||
|
TRandom3 rand;
|
||||||
|
rand.SetSeed(); // for the A1C1 dither baseline in the dither-vs-cfrac comparison
|
||||||
for (auto sx3event : SX3_Events)
|
for (auto sx3event : SX3_Events)
|
||||||
{
|
{
|
||||||
if (sx3event.Energy1 < 1.2)
|
if (sx3event.Energy1 < 1.2)
|
||||||
|
|
@ -2747,6 +3009,58 @@ void protonMiscHistograms_sx3(HistPlotter *plotter, std::vector<Event> QQQ_Event
|
||||||
|
|
||||||
// plotter->Fill1D("pmisc_Ex_from_protons",200,-10,10,apkin_p.getExc(sx3Efix,theta_s*180/M_PI),pmlabel);
|
// plotter->Fill1D("pmisc_Ex_from_protons",200,-10,10,apkin_p.getExc(sx3Efix,theta_s*180/M_PI),pmlabel);
|
||||||
|
|
||||||
} // end PCEvents loop
|
} // end PCEvents loop (A1C2 main flow)
|
||||||
|
|
||||||
|
// --- A1C1 dither vs cfrac comparison (keep BOTH methods) ---
|
||||||
|
// The A1C2 loop above does not touch A1C1 events. Here we reconstruct the PC z
|
||||||
|
// for genuine A1C1 events two ways and fill a parallel set of excitation
|
||||||
|
// observables (suffix _dither / _cfrac) so the methods can be overlaid. The
|
||||||
|
// dither is a Gaussian-smeared crossover z baseline; cfrac is the linear
|
||||||
|
// centre-fold sub-cell model (cell anchored on the fired cathode, side from
|
||||||
|
// the anode-only z at the SX3 phi, offset from the per-cell autocal).
|
||||||
|
for (auto pcevent : PC_Events)
|
||||||
|
{
|
||||||
|
if (!(pcevent.multi1 == 1 && pcevent.multi2 == 1))
|
||||||
|
continue;
|
||||||
|
bool phicut = sx3event.pos.Phi() <= pcevent.pos.Phi() + TMath::Pi() / 3. && sx3event.pos.Phi() >= pcevent.pos.Phi() - TMath::Pi() / 3.;
|
||||||
|
if (!phicut)
|
||||||
|
continue;
|
||||||
|
if (!(pcevent.Energy2 > 1400))
|
||||||
|
continue;
|
||||||
|
|
||||||
|
auto fillCmp = [&](double pcz, const std::string &m)
|
||||||
|
{
|
||||||
|
TVector3 x2(pcevent.pos.X(), pcevent.pos.Y(), pcz);
|
||||||
|
TVector3 vv = x2 - sx3event.pos;
|
||||||
|
double tm = -1.0 * (sx3event.pos.X() * vv.X() + sx3event.pos.Y() * vv.Y()) / (vv.X() * vv.X() + vv.Y() * vv.Y());
|
||||||
|
TVector3 rv = sx3event.pos + tm * vv;
|
||||||
|
if (rv.Perp() > 10.0)
|
||||||
|
return;
|
||||||
|
double th = (sx3event.pos - rv).Theta();
|
||||||
|
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);
|
||||||
|
};
|
||||||
|
|
||||||
|
fillCmp(rand.Gaus(pcevent.pos.Z(), 8.0), "dither"); // method 1: Gaussian dither baseline
|
||||||
|
|
||||||
|
// method 2: cfrac sub-cell linear centre-fold (side ref = anode-only z
|
||||||
|
// rebuilt from the fired anode wire)
|
||||||
|
double ac = pcevent.Energy1 + pcevent.Energy2;
|
||||||
|
double cfrac = (ac > 0.0) ? pcevent.Energy2 / ac : -1.0;
|
||||||
|
if (cfrac >= 0.0)
|
||||||
|
{
|
||||||
|
std::vector<std::tuple<int, double, double>> 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);
|
||||||
|
if (s.inband && s.pitchok)
|
||||||
|
fillCmp(s.pcz, "cfrac");
|
||||||
|
}
|
||||||
|
} // end A1C1 comparison loop
|
||||||
} // end sx3Events loop
|
} // end sx3Events loop
|
||||||
}
|
}
|
||||||
186
fit_a1c1_cfrac.C
Normal file
186
fit_a1c1_cfrac.C
Normal file
|
|
@ -0,0 +1,186 @@
|
||||||
|
#include <TFile.h>
|
||||||
|
#include <TH1.h>
|
||||||
|
#include <TH2.h>
|
||||||
|
#include <TKey.h>
|
||||||
|
#include <TDirectory.h>
|
||||||
|
#include <TMath.h>
|
||||||
|
#include <TString.h>
|
||||||
|
#include <algorithm>
|
||||||
|
#include <fstream>
|
||||||
|
#include <iostream>
|
||||||
|
#include <string>
|
||||||
|
|
||||||
|
namespace {
|
||||||
|
// cathode wire grid; cell i spans zg[i] (high) .. zg[i+1] (low)
|
||||||
|
const double zg[8] = {147.998, 101.946, 59.7634, 19.6965,
|
||||||
|
-19.6965, -59.7634, -101.946, -147.998};
|
||||||
|
|
||||||
|
int cellOf(double z) {
|
||||||
|
for (int i = 0; i < 7; ++i)
|
||||||
|
if (z <= zg[i] && z > zg[i + 1]) return i;
|
||||||
|
return -1;
|
||||||
|
}
|
||||||
|
|
||||||
|
// recursively search a directory tree for a TH2 by name
|
||||||
|
TH2 *findTH2(TDirectory *dir, const std::string &name) {
|
||||||
|
if (TH2 *h = dynamic_cast<TH2 *>(dir->Get(name.c_str()))) return h;
|
||||||
|
TIter next(dir->GetListOfKeys());
|
||||||
|
TKey *key;
|
||||||
|
while ((key = (TKey *)next())) {
|
||||||
|
TObject *obj = key->ReadObj();
|
||||||
|
if (TDirectory *sub = dynamic_cast<TDirectory *>(obj)) {
|
||||||
|
if (TH2 *h = findTH2(sub, name)) return h;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return nullptr;
|
||||||
|
}
|
||||||
|
} // namespace
|
||||||
|
|
||||||
|
// =============================================================================
|
||||||
|
// GLOBAL WITH REFERENCE (Least Squares Fit)
|
||||||
|
// =============================================================================
|
||||||
|
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 = "") {
|
||||||
|
TFile *f = TFile::Open(filename, "READ");
|
||||||
|
if (!f || f->IsZombie()) {
|
||||||
|
std::cerr << "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;
|
||||||
|
|
||||||
|
for (int s = 0; s < 2; ++s) {
|
||||||
|
TH2 *h = findTH2(f, names[s]);
|
||||||
|
if (!h) continue;
|
||||||
|
|
||||||
|
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;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Global least-squares calculation
|
||||||
|
double cfmin_global = 0.40, k_global = 0.075;
|
||||||
|
const char *status = "FALLBACK";
|
||||||
|
|
||||||
|
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";
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
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");
|
||||||
|
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");
|
||||||
|
|
||||||
|
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");
|
||||||
|
|
||||||
|
if (write_dat) {
|
||||||
|
const std::string fn = std::string("a1c1_cfrac_global_") + 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.close();
|
||||||
|
}
|
||||||
|
|
||||||
|
f->Close();
|
||||||
|
}
|
||||||
|
|
||||||
|
// =============================================================================
|
||||||
|
// GLOBAL REFERENCE-FREE (Percentile Fit)
|
||||||
|
// =============================================================================
|
||||||
|
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) {
|
||||||
|
TFile *f = TFile::Open(filename, "READ");
|
||||||
|
if (!f || f->IsZombie()) return;
|
||||||
|
|
||||||
|
const char *names[2] = {"Benchmark_QQQ_trueA1C1_cfrac_vs_cell",
|
||||||
|
"Benchmark_SX3_trueA1C1_cfrac_vs_cell"};
|
||||||
|
|
||||||
|
TH2 *acc = nullptr;
|
||||||
|
for (int s = 0; s < 2; ++s) {
|
||||||
|
TH2 *h = findTH2(f, names[s]);
|
||||||
|
if (!h) continue;
|
||||||
|
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; }
|
||||||
|
|
||||||
|
// Project the ENTIRE detector onto Y to get a global distribution
|
||||||
|
TH1D *proj = acc->ProjectionY("proj_global");
|
||||||
|
const double n_tot = proj->Integral();
|
||||||
|
|
||||||
|
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";
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
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);
|
||||||
|
|
||||||
|
// Print arrays
|
||||||
|
printf("\n// ---- paste into TrackRecon.C (Global Ref-Free) ----\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");
|
||||||
|
|
||||||
|
printf("static const double a1c1_k_%s[7] = {", dataset);
|
||||||
|
for (int c = 0; c < 7; ++c) printf("%.6f%s", k_global, c < 6 ? ", " : "};\n");
|
||||||
|
|
||||||
|
delete proj;
|
||||||
|
delete acc;
|
||||||
|
f->Close();
|
||||||
|
}
|
||||||
59
run_tr.sh
59
run_tr.sh
|
|
@ -30,6 +30,47 @@ export -f process_run
|
||||||
|
|
||||||
export CO2percent=3
|
export CO2percent=3
|
||||||
|
|
||||||
|
if [[ 1 -eq 0 ]]; then
|
||||||
|
export DATASET=27Al Gain=2 a1c1_calibrate=1
|
||||||
|
# root -q -l -b -e '
|
||||||
|
# TChain *t = new TChain("tree");
|
||||||
|
# t->Add("../ANASEN_analysis/data/27Al_Data/Run_015_mapped.root");
|
||||||
|
# t->Add("../ANASEN_analysis/data/27Al_Data/Run_017_mapped.root");
|
||||||
|
# t->Add("../ANASEN_analysis/data/27Al_Data/Run_018_mapped.root");
|
||||||
|
# t->Add("../ANASEN_analysis/data/27Al_Data/Run_019_mapped.root");
|
||||||
|
# t->Add("../ANASEN_analysis/data/27Al_Data/Run_020_mapped.root");
|
||||||
|
# t->Add("../ANASEN_analysis/data/27Al_Data/Run_021_mapped.root");
|
||||||
|
# t->Add("../ANASEN_analysis/data/27Al_Data/Run_022_mapped.root");
|
||||||
|
# t->Process("TrackRecon.C+", "Output_cal/calib_27Al.root");
|
||||||
|
# '
|
||||||
|
export DATASET=17F Gain=1
|
||||||
|
root -q -l -b -e '
|
||||||
|
TChain *t = new TChain("tree");
|
||||||
|
// proton runs, 3% CO2
|
||||||
|
t->Add("../ANASEN_analysis/data/17F_Data/ProtonRun_044_mapped.root");
|
||||||
|
t->Add("../ANASEN_analysis/data/17F_Data/ProtonRun_045_mapped.root");
|
||||||
|
t->Add("../ANASEN_analysis/data/17F_Data/ProtonRun_046_mapped.root");
|
||||||
|
t->Add("../ANASEN_analysis/data/17F_Data/ProtonRun_047_mapped.root");
|
||||||
|
t->Add("../ANASEN_analysis/data/17F_Data/ProtonRun_048_mapped.root");
|
||||||
|
// proton runs, 4% CO2
|
||||||
|
t->Add("../ANASEN_analysis/data/17F_Data/ProtonRun_038_mapped.root");
|
||||||
|
t->Add("../ANASEN_analysis/data/17F_Data/ProtonRun_039_mapped.root");
|
||||||
|
t->Add("../ANASEN_analysis/data/17F_Data/ProtonRun_040_mapped.root");
|
||||||
|
t->Add("../ANASEN_analysis/data/17F_Data/ProtonRun_041_mapped.root");
|
||||||
|
t->Add("../ANASEN_analysis/data/17F_Data/ProtonRun_042_mapped.root");
|
||||||
|
// source runs (extra stats; same Gain=1)
|
||||||
|
t->Add("../ANASEN_analysis/data/17F_Data/SourceRun_018_mapped.root");
|
||||||
|
t->Add("../ANASEN_analysis/data/17F_Data/SourceRun_019_mapped.root");
|
||||||
|
t->Add("../ANASEN_analysis/data/17F_Data/SourceRun_020_mapped.root");
|
||||||
|
t->Add("../ANASEN_analysis/data/17F_Data/SourceRun_021_mapped.root");
|
||||||
|
t->Process("TrackRecon.C+", "Output_cal/calib_17F.root");
|
||||||
|
'
|
||||||
|
unset a1c1_calibrate
|
||||||
|
unset DATASET
|
||||||
|
unset Gain
|
||||||
|
exit
|
||||||
|
fi
|
||||||
|
|
||||||
# --- Block 1: 27Al Source Runs No Gas (1-8) ---
|
# --- Block 1: 27Al Source Runs No Gas (1-8) ---
|
||||||
if [[ 1 -eq 0 ]]; then
|
if [[ 1 -eq 0 ]]; then
|
||||||
export DATASET="27Al"
|
export DATASET="27Al"
|
||||||
|
|
@ -46,12 +87,12 @@ if [[ 1 -eq 1 ]]; then
|
||||||
export DATASET="27Al"
|
export DATASET="27Al"
|
||||||
export PREFIX="Run_"
|
export PREFIX="Run_"
|
||||||
export OUT_DIR="Output_a"
|
export OUT_DIR="Output_a"
|
||||||
export pressure_torr=350
|
export Gain=2
|
||||||
echo "Processing 27Al alpha+gas runs..."
|
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"
|
export source_vertex=-5.36; export timecut_low=12.0; export timecut_high=119.0; process_run 9 "$slope"
|
||||||
unset timecut_high
|
unset timecut_high
|
||||||
export source_vertex=53.44; export timecut_low=400.0; process_run 12 "$slope"
|
export source_vertex=53.44; export timecut_low=400.0; process_run 12 "$slope"
|
||||||
unset pressuer_torr
|
unset Gain
|
||||||
unset timecut_low
|
unset timecut_low
|
||||||
fi
|
fi
|
||||||
|
|
||||||
|
|
@ -59,8 +100,8 @@ fi
|
||||||
if [[ 1 -eq 0 ]]; then
|
if [[ 1 -eq 0 ]]; then
|
||||||
export DATASET="27Al"
|
export DATASET="27Al"
|
||||||
export PREFIX="Run_"
|
export PREFIX="Run_"
|
||||||
export OUT_DIR="Output_p"
|
export OUT_DIR="Output_p"
|
||||||
export pressure_torr=350
|
export Gain=2
|
||||||
|
|
||||||
rm -f ${OUT_DIR}/Al_protons.root
|
rm -f ${OUT_DIR}/Al_protons.root
|
||||||
export source_vertex=-200.0 # Source on the entrance window
|
export source_vertex=-200.0 # Source on the entrance window
|
||||||
|
|
@ -68,8 +109,8 @@ if [[ 1 -eq 0 ]]; then
|
||||||
|
|
||||||
# process_run 17
|
# process_run 17
|
||||||
parallel --bar -j 6 process_run ::: 15 {17..22}
|
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
|
hadd -j 4 -k ${OUT_DIR}/Al_protons.root ${OUT_DIR}/results_run0{15..22}.root
|
||||||
unset pressuer_torr
|
unset Gain
|
||||||
|
|
||||||
fi
|
fi
|
||||||
|
|
||||||
|
|
@ -84,7 +125,7 @@ if [[ 1 -eq 0 ]]; then
|
||||||
fi
|
fi
|
||||||
|
|
||||||
# --- Block 5: 17F Alpha Run with Gas (18-21) ---
|
# --- Block 5: 17F Alpha Run with Gas (18-21) ---
|
||||||
if [[ 1 -eq 0 ]]; then
|
if [[ 1 -eq 1 ]]; then
|
||||||
export DATASET="17F"
|
export DATASET="17F"
|
||||||
export PREFIX="SourceRun_"
|
export PREFIX="SourceRun_"
|
||||||
export OUT_DIR="Output_a"
|
export OUT_DIR="Output_a"
|
||||||
|
|
@ -98,12 +139,12 @@ if [[ 1 -eq 0 ]]; then
|
||||||
fi
|
fi
|
||||||
|
|
||||||
# --- Block 6: 17F Proton Data ---
|
# --- Block 6: 17F Proton Data ---
|
||||||
if [[ 1 -eq 0 ]]; then
|
if [[ 1 -eq 1 ]]; then
|
||||||
export DATASET="17F"
|
export DATASET="17F"
|
||||||
export PREFIX="ProtonRun_"
|
export PREFIX="ProtonRun_"
|
||||||
export OUT_DIR="Output_p"
|
export OUT_DIR="Output_p"
|
||||||
rm -f ${OUT_DIR}/*pc*.root
|
rm -f ${OUT_DIR}/*pc*.root
|
||||||
export source_vertex=-57.28
|
export source_vertex=-200.0
|
||||||
|
|
||||||
parallel --bar -j 6 process_run ::: {44..48} #3% CO2
|
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
|
hadd -j 4 -k ${OUT_DIR}/3pc.root ${OUT_DIR}/results_run0{44..48}.root
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue
Block a user