modified: Armory/ClassPW.h fixed typo
modified: TrackRecon.C found some interesting behaviour in the A1c1 plots going to implement new a1c1 model invert to acount for a/c charege variation next modified: scratch/plot_dither_scan.C modified: scratch/scan_dither_runs.sh
This commit is contained in:
parent
1bbf2ae059
commit
fcd412dc59
|
|
@ -86,7 +86,7 @@ public:
|
|||
double GetAnodeTheta(short id) const { return (An[id].first - An[id].second).Theta(); }
|
||||
double GetAnodePhi(short id) const { return (An[id].first - An[id].second).Phi(); }
|
||||
|
||||
TVector3 GetCathodneMid(short id) const { return (Ca[id].first + Ca[id].second) * 0.5; }
|
||||
TVector3 GetCathodeMid(short id) const { return (Ca[id].first + Ca[id].second) * 0.5; }
|
||||
double GetCathodeTheta(short id) const { return (Ca[id].first - Ca[id].second).Theta(); }
|
||||
double GetCathodePhi(short id) const { return (Ca[id].first - Ca[id].second).Phi(); }
|
||||
|
||||
|
|
|
|||
229
TrackRecon.C
229
TrackRecon.C
|
|
@ -50,6 +50,7 @@ const double qqq_z = 105.0;
|
|||
double z_entrance = -174.3 - 9.7 - 100.0;
|
||||
const double anode_gain = 1.5146e-5; // channels --> MeV
|
||||
double dither_sigma = 8.0;
|
||||
double dither_sigma_c0 = dither_sigma;
|
||||
std::string dataset;
|
||||
int CO2percent;
|
||||
bool reactiondata = false;
|
||||
|
|
@ -181,6 +182,7 @@ void TrackRecon::Begin(TTree * /*tree*/)
|
|||
if (getenv("DITHER_SIGMA"))
|
||||
{
|
||||
dither_sigma = std::atof(getenv("DITHER_SIGMA"));
|
||||
dither_sigma_c0 = dither_sigma;
|
||||
std::cout << "Dither Sigma set to " << dither_sigma << " mm" << std::endl;
|
||||
}
|
||||
|
||||
|
|
@ -1510,18 +1512,17 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
|
|||
plotter->Fill1D("PCZ_sx3", 800, -200, 200, pcevent.pos.Z(), "Z_Reconstruction");
|
||||
}
|
||||
|
||||
//-----------------------Benchmarking Method for Source Runs------------------------//
|
||||
|
||||
//-----------------------Benchmarking Method for Source Runs (SX3)------------------------//
|
||||
if (BenchMark && aClusters.size() == 1 && cClusters.size() == 1)
|
||||
{
|
||||
const auto &aCl = aClusters.front();
|
||||
const auto &cCl = cClusters.front();
|
||||
|
||||
auto vertexFromPCPoint = [&x1](const TVector3 &x2f)
|
||||
auto vertexFrom = [](const TVector3 &si, const TVector3 &pcpoint)
|
||||
{
|
||||
TVector3 vf = x2f - x1;
|
||||
double tm = -1.0 * (x1.X() * vf.X() + x1.Y() * vf.Y()) / (vf.X() * vf.X() + vf.Y() * vf.Y());
|
||||
return TVector3(x1 + tm * vf);
|
||||
TVector3 vf = pcpoint - si;
|
||||
double tm = -1.0 * (si.X() * vf.X() + si.Y() * vf.Y()) / (vf.X() * vf.X() + vf.Y() * vf.Y());
|
||||
return TVector3(si + tm * vf);
|
||||
};
|
||||
|
||||
auto fillSuite = [&](const std::string &tag, double pcz_method, const TVector3 &vtx)
|
||||
|
|
@ -1530,58 +1531,93 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
|
|||
plotter->Fill1D("Benchmark_SX3_VertexZ_" + tag + "_TC" + std::to_string(PCSX3TimeCut) + "_PC" + std::to_string(phicut), 800, -400, 400, vtx.Z(), "Benchmark_SX3");
|
||||
plotter->Fill2D("Benchmark_SX3_VertexXY_" + tag, 200, -100, 100, 200, -100, 100, vtx.X(), vtx.Y(), "Benchmark_SX3");
|
||||
plotter->Fill1D("Benchmark_SX3_PCZ_" + tag, 600, -200, 200, pcz_method, "Benchmark_SX3");
|
||||
plotter->Fill2D("Benchmark_SX3_PCZ_vs_sx3pczguess_" + tag, 600, -200, 200, 600, -200, 200, pcz_guess_int, pcz_method, "Benchmark_SX3");
|
||||
plotter->Fill2D("Benchmark_SX3_PCZresidual_vs_pczguess_" + tag, 600, -200, 200, 200, -100, 100, pczguess, pcz_method - pczguess, "Benchmark_SX3");
|
||||
plotter->Fill1D("Benchmark_SX3_PCZ-sx3pczguess_" + tag, 200, -100, 100, pcz_method - pczguess, "Benchmark_SX3");
|
||||
plotter->Fill1D("Benchmark_SX3_PCZ-sx3pczint_" + tag, 200, -100, 100, pcz_method - pcz_guess_int, "Benchmark_SX3");
|
||||
};
|
||||
|
||||
auto fillVsRef = [&](const std::string &tag, double pcz_method, const TVector3 &vtx, double pcz_ref, const TVector3 &vtx_ref)
|
||||
{
|
||||
plotter->Fill2D("Benchmark_SX3_PCZ_" + tag + "_vs_ref", 400, -200, 200, 400, -200, 200, pcz_ref, pcz_method, "Benchmark_SX3_ref");
|
||||
plotter->Fill1D("Benchmark_SX3_PCZ_" + tag + "_minus_ref", 400, -100, 100, pcz_method - pcz_ref, "Benchmark_SX3_ref");
|
||||
plotter->Fill2D("Benchmark_SX3_VertexZ_" + tag + "_vs_ref", 400, -200, 200, 400, -200, 200, vtx_ref.Z(), vtx.Z(), "Benchmark_SX3_ref");
|
||||
plotter->Fill1D("Benchmark_SX3_VertexZ_" + tag + "_minus_ref", 400, -100, 100, vtx.Z() - vtx_ref.Z(), "Benchmark_SX3_ref");
|
||||
plotter->Fill2D("Benchmark_SX3_PCZ_" + tag + "_vs_sx3pczguess", 400, -200, 200, 400, -200, 200, pczguess, pcz_method, "Benchmark_SX3_ref");
|
||||
plotter->Fill1D("Benchmark_SX3_PCZ_" + tag + "_minus_sx3pczguess", 400, -100, 100, pcz_method - pczguess, "Benchmark_SX3_ref");
|
||||
};
|
||||
|
||||
// cathode charge-division reference on this event (the A1C2 recipe)
|
||||
double pcz_ref = pcfix_func.Eval(pcevent.pos.Z());
|
||||
TVector3 vtx_ref = vertexFromPCPoint(TVector3(pcevent.pos.X(), pcevent.pos.Y(), pcz_ref));
|
||||
TVector3 vtx_ref = vertexFrom(sx3event.pos, TVector3(pcevent.pos.X(), pcevent.pos.Y(), pcz_ref));
|
||||
|
||||
// anode-only PC point, oneWire method: pseudo-wire interpolated to the
|
||||
// SX3 hit's phi, with the same quality cuts as miscHistograms_oneWire
|
||||
auto [apwire_bm, apSumE_bm, apMaxE_bm, apTSMaxE_bm] = pwinstance.GetPseudoWire(aCl, "ANODE");
|
||||
TVector3 pc_anodeOnly = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, sx3event.pos.Phi());
|
||||
pc_anodeOnly.SetZ(rand.Gaus(pc_anodeOnly.Z(), dither_sigma));
|
||||
TVector3 vtx_anodeOnly = vertexFromPCPoint(pc_anodeOnly);
|
||||
bool anodeOnlyGood = vtx_anodeOnly.Perp() <= 6.0 && vtx_anodeOnly.Z() >= -173.6 && vtx_anodeOnly.Z() <= 100 && phicut;
|
||||
auto pw_tuple = pwinstance.GetPseudoWire(aCl, "ANODE");
|
||||
std::pair<TVector3, TVector3> apwire_bm = std::get<0>(pw_tuple);
|
||||
|
||||
if (pcevent.multi1 == 1 && pcevent.multi2 == 2)
|
||||
auto cMaxWire = *std::max_element(cCl.begin(), cCl.end(), [](const auto &a, const auto &b)
|
||||
{ return std::get<1>(a) < std::get<1>(b); });
|
||||
std::vector<std::tuple<int, double, double>> cOne = {cMaxWire};
|
||||
|
||||
auto xo_tuple = pwinstance.FindCrossoverProperties(aCl, cOne);
|
||||
TVector3 xo_a1c1 = std::get<0>(xo_tuple);
|
||||
double alpha_a1c1 = std::get<1>(xo_tuple);
|
||||
bool a1c1Good = (alpha_a1c1 != 9999999 && std::get<2>(xo_tuple) != -1);
|
||||
|
||||
// --- A1C1 charge-ratio diagnostic ---
|
||||
// Only one cathode in A1C1, so the z-sensitive variable is cathode-vs-anode
|
||||
// charge, not cathode/cathode. Use pseudo-wire sums for consistent gains.
|
||||
double aSumE_bm = std::get<1>(pw_tuple); // anode sum (reuse pw_tuple)
|
||||
auto cpw_tuple = pwinstance.GetPseudoWire(cOne, "CATHODE");
|
||||
double cSumE_bm = std::get<1>(cpw_tuple); // single-cathode sum
|
||||
double ac_sum = aSumE_bm + cSumE_bm;
|
||||
double cfrac = (ac_sum > 0.0) ? cSumE_bm / ac_sum : -1.0; // bounded [0,1]
|
||||
|
||||
// Smearing Setup
|
||||
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);
|
||||
TVector3 smeared_sx3_pos(sx3event.pos.Perp() * TMath::Cos(smeared_phi), sx3event.pos.Perp() * TMath::Sin(smeared_phi), sx3event.pos.Z());
|
||||
|
||||
auto doA1C1 = [&](const std::string &tag, const TVector3 &si_point, bool dither, bool hybrid)
|
||||
{
|
||||
fillSuite("A1C2", pcz_ref, vtx_ref); // baseline with identical binning
|
||||
if (!a1c1Good)
|
||||
return;
|
||||
// double sigma = hybrid ? (dither_sigma_c0 / 2.0) : dither_sigma_c0;
|
||||
// double pcz = dither ? rand.Gaus(xo_a1c1.Z(), sigma) : xo_a1c1.Z();
|
||||
double pcz = dither ? rand.Gaus(xo_a1c1.Z(), dither_sigma) : xo_a1c1.Z();
|
||||
// double pcz = dither ? rand.Uniform(xo_a1c1.Z(), dither_sigma) : xo_a1c1.Z();
|
||||
TVector3 vtx = vertexFrom(si_point, TVector3(xo_a1c1.X(), xo_a1c1.Y(), pcz));
|
||||
fillSuite(tag, pcz, vtx);
|
||||
fillVsRef(tag, pcz, vtx, pcz_ref, vtx_ref);
|
||||
};
|
||||
|
||||
double pcz_dith = rand.Gaus(pcevent.pos.Z(), dither_sigma);
|
||||
|
||||
TVector3 vtx_a1c1 = vertexFromPCPoint(TVector3(pcevent.pos.X(), pcevent.pos.Y(), pcz_dith));
|
||||
fillSuite("A1C1", pcz_dith, vtx_a1c1);
|
||||
fillVsRef("A1C1", pcz_dith, vtx_a1c1, pcz_ref, vtx_ref);
|
||||
|
||||
// ---- A1C0 emulation: oneWire method, cathodes ignored ----
|
||||
if (anodeOnlyGood)
|
||||
{
|
||||
fillSuite("A1C0", pc_anodeOnly.Z(), vtx_anodeOnly);
|
||||
fillVsRef("A1C0", pc_anodeOnly.Z(), vtx_anodeOnly, pcz_ref, vtx_ref);
|
||||
}
|
||||
}
|
||||
|
||||
// ---- A2C0 emulation: 2-wire anode pseudo-wire, oneWire method ----
|
||||
if (pcevent.multi1 == 2 && pcevent.multi2 == 2)
|
||||
auto doAnodeOnly = [&](const std::string &tag, double phi_use, const TVector3 &si_point, bool dither, bool hybrid)
|
||||
{
|
||||
fillSuite("A2C2", pcz_ref, vtx_ref); // reference population for A2C0
|
||||
if (anodeOnlyGood)
|
||||
TVector3 pc = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, phi_use);
|
||||
TVector3 vtx0 = vertexFrom(si_point, pc);
|
||||
if (!(vtx0.Perp() <= 6.0 && vtx0.Z() >= -173.6))
|
||||
return;
|
||||
// double sigma = hybrid ? (dither_sigma_c0 / 2.0) : dither_sigma_c0;
|
||||
double pcz = dither ? rand.Gaus(pc.Z(), dither_sigma) : pc.Z();
|
||||
TVector3 vtx = vertexFrom(si_point, TVector3(pc.X(), pc.Y(), pcz));
|
||||
fillSuite(tag, pcz, vtx);
|
||||
fillVsRef(tag, pcz, vtx, pcz_ref, vtx_ref);
|
||||
};
|
||||
|
||||
if (phicut && PCSX3TimeCut)
|
||||
{
|
||||
if (pcevent.multi1 == 1 && pcevent.multi2 == 2)
|
||||
{
|
||||
fillSuite("A2C0", pc_anodeOnly.Z(), vtx_anodeOnly);
|
||||
fillVsRef("A2C0", pc_anodeOnly.Z(), vtx_anodeOnly, pcz_ref, vtx_ref);
|
||||
fillSuite("A1C2", pcz_ref, vtx_ref);
|
||||
doA1C1("A1C1", sx3event.pos, true, false);
|
||||
doAnodeOnly("A1C0", sx3event.pos.Phi(), sx3event.pos, true, false);
|
||||
doA1C1("A1C1_Si", smeared_sx3_pos, false, false);
|
||||
doAnodeOnly("A1C0_Si", smeared_phi, smeared_sx3_pos, false, false);
|
||||
doA1C1("A1C1_Hyb", smeared_sx3_pos, true, true);
|
||||
doAnodeOnly("A1C0_Hyb", smeared_phi, smeared_sx3_pos, true, true);
|
||||
|
||||
// --- A1C1 charge-fraction diagnostics ---
|
||||
// Decides whether the single cathode carries any sub-cell z information.
|
||||
if (a1c1Good && cfrac >= 0.0)
|
||||
{
|
||||
plotter->Fill1D("Benchmark_SX3_A1C1_cfrac", 220, -0.05, 1.05, cfrac, "Benchmark_SX3_ref");
|
||||
plotter->Fill2D("Benchmark_SX3_A1C1_cfrac_vs_ref", 400, -200, 200, 220, -0.05, 1.05,
|
||||
pcz_ref, cfrac, "Benchmark_SX3_ref");
|
||||
plotter->Fill2D("Benchmark_SX3_A1C1_cfrac_vs_sx3pczguess", 400, -200, 200, 220, -0.05, 1.05,
|
||||
pczguess, cfrac, "Benchmark_SX3_ref");
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -1719,18 +1755,17 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
|
|||
plotter->Fill2D("dE3_Ef_CathodeQQQR_TC1PC" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath), 600, 0, 15, 800, 0, 10000, qqqEfix, pcevent.Energy2 * sinTheta_customV, "PID_dE_E");
|
||||
}
|
||||
|
||||
//-----------------------Benchmarking Method for Source Runs------------------------//
|
||||
|
||||
//-----------------------Benchmarking Method for Source Runs (QQQ)------------------------//
|
||||
if (BenchMark && aClusters.size() == 1 && cClusters.size() == 1)
|
||||
{
|
||||
const auto &aCl = aClusters.front();
|
||||
const auto &cCl = cClusters.front();
|
||||
|
||||
auto vertexFromPCPoint = [&x1](const TVector3 &x2f)
|
||||
auto vertexFrom = [](const TVector3 &si, const TVector3 &pcpoint)
|
||||
{
|
||||
TVector3 vf = x2f - x1;
|
||||
double tm = -1.0 * (x1.X() * vf.X() + x1.Y() * vf.Y()) / (vf.X() * vf.X() + vf.Y() * vf.Y());
|
||||
return TVector3(x1 + tm * vf);
|
||||
TVector3 vf = pcpoint - si;
|
||||
double tm = -1.0 * (si.X() * vf.X() + si.Y() * vf.Y()) / (vf.X() * vf.X() + vf.Y() * vf.Y());
|
||||
return TVector3(si + tm * vf);
|
||||
};
|
||||
|
||||
auto fillSuite = [&](const std::string &tag, double pcz_method, const TVector3 &vtx)
|
||||
|
|
@ -1739,69 +1774,75 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
|
|||
plotter->Fill1D("Benchmark_QQQ_VertexZ_" + tag + "_TC" + std::to_string(timecut) + "_PC" + std::to_string(phicut), 800, -400, 400, vtx.Z(), "Benchmark_QQQ");
|
||||
plotter->Fill2D("Benchmark_QQQ_VertexXY_" + tag, 200, -100, 100, 200, -100, 100, vtx.X(), vtx.Y(), "Benchmark_QQQ");
|
||||
plotter->Fill1D("Benchmark_QQQ_PCZ_" + tag, 600, -200, 200, pcz_method, "Benchmark_QQQ");
|
||||
plotter->Fill2D("Benchmark_QQQ_PCZ_vs_qqqpczguess_" + tag, 600, -200, 200, 600, -200, 200, pcz_guess_int, pcz_method, "Benchmark_QQQ");
|
||||
plotter->Fill2D("Benchmark_QQQ_PCZresidual_vs_pczguess_" + tag, 600, -200, 200, 200, -100, 100, pcz_guess_37, pcz_method - pcz_guess_37, "Benchmark_QQQ");
|
||||
plotter->Fill1D("Benchmark_QQQ_PCZ-qqqpczguess_" + tag, 200, -100, 100, pcz_method - pcz_guess_37, "Benchmark_QQQ");
|
||||
plotter->Fill1D("Benchmark_QQQ_PCZ-qqqpczint_" + tag, 200, -100, 100, pcz_method - pcz_guess_int, "Benchmark_QQQ");
|
||||
};
|
||||
|
||||
auto fillVsRef = [&](const std::string &tag, double pcz_method, const TVector3 &vtx, double pcz_ref, const TVector3 &vtx_ref)
|
||||
{
|
||||
plotter->Fill2D("Benchmark_QQQ_PCZ_" + tag + "_vs_ref", 400, -200, 200, 400, -200, 200, pcz_ref, pcz_method, "Benchmark_QQQ_ref");
|
||||
plotter->Fill1D("Benchmark_QQQ_PCZ_" + tag + "_minus_ref", 400, -100, 100, pcz_method - pcz_ref, "Benchmark_QQQ_ref");
|
||||
plotter->Fill2D("Benchmark_QQQ_VertexZ_" + tag + "_vs_ref", 400, -200, 200, 400, -200, 200, vtx_ref.Z(), vtx.Z(), "Benchmark_QQQ_ref");
|
||||
plotter->Fill1D("Benchmark_QQQ_VertexZ_" + tag + "_minus_ref", 400, -100, 100, vtx.Z() - vtx_ref.Z(), "Benchmark_QQQ_ref");
|
||||
};
|
||||
|
||||
// cathode charge-division reference on this event (the A1C2 recipe)
|
||||
double pcz_ref = pcfix_func.Eval(pcevent.pos.Z());
|
||||
TVector3 vtx_ref = vertexFromPCPoint(TVector3(pcevent.pos.X(), pcevent.pos.Y(), pcz_ref));
|
||||
TVector3 vtx_ref = vertexFrom(qqqevent.pos, TVector3(pcevent.pos.X(), pcevent.pos.Y(), pcz_ref));
|
||||
|
||||
// anode-only PC point, oneWire method: pseudo-wire interpolated to the
|
||||
// QQQ hit's phi, with the same quality cuts as miscHistograms_oneWire
|
||||
auto [apwire_bm, apSumE_bm, apMaxE_bm, apTSMaxE_bm] = pwinstance.GetPseudoWire(aCl, "ANODE");
|
||||
TVector3 pc_anodeOnly = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, qqqevent.pos.Phi());
|
||||
pc_anodeOnly.SetZ(rand.Gaus(pc_anodeOnly.Z(), dither_sigma));
|
||||
TVector3 vtx_anodeOnly = vertexFromPCPoint(pc_anodeOnly);
|
||||
bool anodeOnlyGood = vtx_anodeOnly.Perp() <= 6.0 && vtx_anodeOnly.Z() >= -173.6 && vtx_anodeOnly.Z() <= 100 && phicut;
|
||||
auto pw_tuple = pwinstance.GetPseudoWire(aCl, "ANODE");
|
||||
std::pair<TVector3, TVector3> apwire_bm = std::get<0>(pw_tuple);
|
||||
|
||||
if (pcevent.multi1 == 1 && pcevent.multi2 == 2)
|
||||
auto cMaxWire = *std::max_element(cCl.begin(), cCl.end(), [](const auto &a, const auto &b)
|
||||
{ return std::get<1>(a) < std::get<1>(b); });
|
||||
std::vector<std::tuple<int, double, double>> cOne = {cMaxWire};
|
||||
|
||||
auto xo_tuple = pwinstance.FindCrossoverProperties(aCl, cOne);
|
||||
TVector3 xo_a1c1 = std::get<0>(xo_tuple);
|
||||
double alpha_a1c1 = std::get<1>(xo_tuple);
|
||||
bool a1c1Good = (alpha_a1c1 != 9999999 && std::get<2>(xo_tuple) != -1);
|
||||
|
||||
// Smearing Setup
|
||||
double qqq_wedge_pitch = (87.0 / 16.0) * (M_PI / 180.0);
|
||||
double qqq_ring_pitch = 48.0 / 16.0;
|
||||
double smeared_phi = qqqevent.pos.Phi() + rand.Uniform(-qqq_wedge_pitch / 2.0, qqq_wedge_pitch / 2.0);
|
||||
double smeared_rho = qqqevent.pos.Perp() + rand.Uniform(-qqq_ring_pitch / 2.0, qqq_ring_pitch / 2.0);
|
||||
|
||||
TVector3 smeared_qqq_pos(smeared_rho * TMath::Cos(smeared_phi), smeared_rho * TMath::Sin(smeared_phi), qqqevent.pos.Z());
|
||||
|
||||
auto doA1C1 = [&](const std::string &tag, const TVector3 &si_point, bool dither, bool hybrid)
|
||||
{
|
||||
fillSuite("A1C2", pcz_ref, vtx_ref); // baseline with identical binning
|
||||
if (!a1c1Good)
|
||||
return;
|
||||
// double sigma = hybrid ? (dither_sigma / 2.0) : dither_sigma;
|
||||
double pcz = dither ? rand.Gaus(xo_a1c1.Z(), dither_sigma) : xo_a1c1.Z();
|
||||
TVector3 vtx = vertexFrom(si_point, TVector3(xo_a1c1.X(), xo_a1c1.Y(), pcz));
|
||||
fillSuite(tag, pcz, vtx);
|
||||
fillVsRef(tag, pcz, vtx, pcz_ref, vtx_ref);
|
||||
};
|
||||
|
||||
// ---- A1C1 emulation: crossover with only the max-E cathode wire (Cmax),
|
||||
// used directly with no z-model ----
|
||||
if (!phicut)
|
||||
continue;
|
||||
if (pcevent.Time1 - qqqevent.Time1 < -150 || pcevent.Time1 - qqqevent.Time1 > 850)
|
||||
continue;
|
||||
|
||||
double pcz_dith = rand.Gaus(pcevent.pos.Z(), dither_sigma);
|
||||
|
||||
TVector3 vtx_a1c1 = vertexFromPCPoint(TVector3(pcevent.pos.X(), pcevent.pos.Y(), pcz_dith));
|
||||
fillSuite("A1C1", pcz_dith, vtx_a1c1);
|
||||
fillVsRef("A1C1", pcz_dith, vtx_a1c1, pcz_ref, vtx_ref);
|
||||
|
||||
// ---- A1C0 emulation: oneWire method, cathodes ignored ----
|
||||
if (anodeOnlyGood)
|
||||
{
|
||||
fillSuite("A1C0", pc_anodeOnly.Z(), vtx_anodeOnly);
|
||||
fillVsRef("A1C0", pc_anodeOnly.Z(), vtx_anodeOnly, pcz_ref, vtx_ref);
|
||||
}
|
||||
}
|
||||
|
||||
// ---- A2C0 emulation: 2-wire anode pseudo-wire, oneWire method ----
|
||||
if (pcevent.multi1 == 2 && pcevent.multi2 == 2)
|
||||
auto doAnodeOnly = [&](const std::string &tag, double phi_use, const TVector3 &si_point, bool dither, bool hybrid)
|
||||
{
|
||||
fillSuite("A2C2", pcz_ref, vtx_ref); // reference population for A2C0
|
||||
if (anodeOnlyGood)
|
||||
TVector3 pc = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, phi_use);
|
||||
TVector3 vtx0 = vertexFrom(si_point, pc);
|
||||
if (!(vtx0.Perp() <= 6.0 && vtx0.Z() >= -173.6))
|
||||
return;
|
||||
double sigma = hybrid ? (dither_sigma_c0 / 2.0) : dither_sigma_c0;
|
||||
double pcz = dither ? rand.Gaus(pc.Z(), sigma) : pc.Z();
|
||||
TVector3 vtx = vertexFrom(si_point, TVector3(pc.X(), pc.Y(), pcz));
|
||||
fillSuite(tag, pcz, vtx);
|
||||
fillVsRef(tag, pcz, vtx, pcz_ref, vtx_ref);
|
||||
};
|
||||
|
||||
if (phicut && timecut)
|
||||
{
|
||||
if (pcevent.multi1 == 1 && pcevent.multi2 == 2)
|
||||
{
|
||||
fillSuite("A2C0", pc_anodeOnly.Z(), vtx_anodeOnly);
|
||||
fillVsRef("A2C0", pc_anodeOnly.Z(), vtx_anodeOnly, pcz_ref, vtx_ref);
|
||||
fillSuite("A1C2", pcz_ref, vtx_ref);
|
||||
doA1C1("A1C1", qqqevent.pos, true, false);
|
||||
doAnodeOnly("A1C0", qqqevent.pos.Phi(), qqqevent.pos, true, false);
|
||||
doA1C1("A1C1_Si", smeared_qqq_pos, false, false);
|
||||
doAnodeOnly("A1C0_Si", smeared_phi, smeared_qqq_pos, false, false);
|
||||
doA1C1("A1C1_Hyb", smeared_qqq_pos, true, true);
|
||||
doAnodeOnly("A1C0_Hyb", smeared_phi, smeared_qqq_pos, true, true);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double qqqrho = qqqevent.pos.Perp();
|
||||
double qqqz = (qqqevent.pos - TVector3(0, 0, source_vertex)).Z();
|
||||
double tan_theta = qqqrho / qqqz;
|
||||
|
|
|
|||
|
|
@ -1,5 +1,6 @@
|
|||
#include <TFile.h>
|
||||
#include <TH1.h>
|
||||
#include <TF1.h>
|
||||
#include <TGraphErrors.h>
|
||||
#include <TGraph.h>
|
||||
#include <TMultiGraph.h>
|
||||
|
|
@ -17,49 +18,77 @@ void plot_dither_scan()
|
|||
// =========================================================================
|
||||
// --- Configuration ---
|
||||
// =========================================================================
|
||||
double dithers[] = {0.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};
|
||||
// double dithers[] = { 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0,
|
||||
// 19.0, 20.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0, 30.0, 32.0, 33.0, 34.0, 35.0, 36.0, 37.0, 38.0};
|
||||
double dithers[] = { 8.0 ,12.0, 16.0, 20.0, 24.0, 28.0, 32.0, 36.0, 40.0};
|
||||
const int N_DITHER = sizeof(dithers) / sizeof(dithers[0]);
|
||||
|
||||
int runs[] = {12, 18, 19, 20, 21};
|
||||
int runs[] = {9, 12, 18, 19, 20, 21};
|
||||
const int N_RUNS = sizeof(runs) / sizeof(runs[0]);
|
||||
|
||||
int runColors[] = {kBlack, kRed, kBlue, kGreen + 2, kMagenta};
|
||||
TString runNames[] = {"27Al Run 12", "17F Run 18", "17F Run 19", "17F Run 20", "17F Run 21"};
|
||||
// Safety: 6 runs defined, 6 colors and 6 names provided.
|
||||
int runColors[] = {kOrange + 1, kBlack, kRed, kBlue, kGreen + 2, kMagenta};
|
||||
TString runNames[] = {"27Al Run 9", "27Al Run 12", "17F Run 18", "17F Run 19", "17F Run 20", "17F Run 21"};
|
||||
|
||||
// Graph arrays
|
||||
TGraphErrors *g_rms_A1C0[N_RUNS];
|
||||
TGraphErrors *g_rms_A1C1[N_RUNS];
|
||||
TGraph *g_kurt_A1C0[N_RUNS];
|
||||
TGraph *g_kurt_A1C1[N_RUNS];
|
||||
// =========================================================================
|
||||
// --- Graph Arrays ---
|
||||
// =========================================================================
|
||||
// A1C0 Arrays
|
||||
TGraphErrors *g_rms_A1C0_orig[N_RUNS];
|
||||
TGraphErrors *g_rms_A1C0_hyb[N_RUNS];
|
||||
TGraphErrors *g_fit_A1C0_orig[N_RUNS];
|
||||
TGraphErrors *g_fit_A1C0_hyb[N_RUNS];
|
||||
TGraph *g_kurt_A1C0_orig[N_RUNS];
|
||||
TGraph *g_kurt_A1C0_hyb[N_RUNS];
|
||||
|
||||
// Initialize graphs
|
||||
// A1C1 Arrays
|
||||
TGraphErrors *g_rms_A1C1_orig[N_RUNS];
|
||||
TGraphErrors *g_rms_A1C1_hyb[N_RUNS];
|
||||
TGraphErrors *g_fit_A1C1_orig[N_RUNS];
|
||||
TGraphErrors *g_fit_A1C1_hyb[N_RUNS];
|
||||
TGraph *g_kurt_A1C1_orig[N_RUNS];
|
||||
TGraph *g_kurt_A1C1_hyb[N_RUNS];
|
||||
|
||||
// Initialize all graphs
|
||||
for (int ir = 0; ir < N_RUNS; ir++)
|
||||
{
|
||||
// RMS (Spread) Graphs
|
||||
g_rms_A1C0[ir] = new TGraphErrors(N_DITHER);
|
||||
g_rms_A1C0[ir]->SetTitle(runNames[ir]);
|
||||
g_rms_A1C0[ir]->SetMarkerStyle(20);
|
||||
g_rms_A1C0[ir]->SetMarkerColor(runColors[ir]);
|
||||
g_rms_A1C0[ir]->SetLineColor(runColors[ir]);
|
||||
// ------------------ A1C0 Initialization ------------------
|
||||
g_rms_A1C0_orig[ir] = new TGraphErrors(N_DITHER);
|
||||
g_rms_A1C0_orig[ir]->SetTitle(runNames[ir]); g_rms_A1C0_orig[ir]->SetMarkerStyle(20); g_rms_A1C0_orig[ir]->SetMarkerColor(runColors[ir]); g_rms_A1C0_orig[ir]->SetLineColor(runColors[ir]);
|
||||
|
||||
g_rms_A1C1[ir] = new TGraphErrors(N_DITHER);
|
||||
g_rms_A1C1[ir]->SetTitle(runNames[ir]);
|
||||
g_rms_A1C1[ir]->SetMarkerStyle(21);
|
||||
g_rms_A1C1[ir]->SetMarkerColor(runColors[ir]);
|
||||
g_rms_A1C1[ir]->SetLineColor(runColors[ir]);
|
||||
g_rms_A1C0_hyb[ir] = new TGraphErrors(N_DITHER);
|
||||
g_rms_A1C0_hyb[ir]->SetTitle(runNames[ir]); g_rms_A1C0_hyb[ir]->SetMarkerStyle(21); g_rms_A1C0_hyb[ir]->SetMarkerColor(runColors[ir]); g_rms_A1C0_hyb[ir]->SetLineColor(runColors[ir]);
|
||||
|
||||
// Kurtosis Graphs
|
||||
g_kurt_A1C0[ir] = new TGraph(N_DITHER);
|
||||
g_kurt_A1C0[ir]->SetTitle(runNames[ir]);
|
||||
g_kurt_A1C0[ir]->SetMarkerStyle(24);
|
||||
g_kurt_A1C0[ir]->SetMarkerColor(runColors[ir]);
|
||||
g_kurt_A1C0[ir]->SetLineColor(runColors[ir]);
|
||||
g_fit_A1C0_orig[ir] = new TGraphErrors(N_DITHER);
|
||||
g_fit_A1C0_orig[ir]->SetTitle(runNames[ir]); g_fit_A1C0_orig[ir]->SetMarkerStyle(22); g_fit_A1C0_orig[ir]->SetMarkerColor(runColors[ir]); g_fit_A1C0_orig[ir]->SetLineColor(runColors[ir]);
|
||||
|
||||
g_kurt_A1C1[ir] = new TGraph(N_DITHER);
|
||||
g_kurt_A1C1[ir]->SetTitle(runNames[ir]);
|
||||
g_kurt_A1C1[ir]->SetMarkerStyle(25);
|
||||
g_kurt_A1C1[ir]->SetMarkerColor(runColors[ir]);
|
||||
g_kurt_A1C1[ir]->SetLineColor(runColors[ir]);
|
||||
g_fit_A1C0_hyb[ir] = new TGraphErrors(N_DITHER);
|
||||
g_fit_A1C0_hyb[ir]->SetTitle(runNames[ir]); g_fit_A1C0_hyb[ir]->SetMarkerStyle(23); g_fit_A1C0_hyb[ir]->SetMarkerColor(runColors[ir]); g_fit_A1C0_hyb[ir]->SetLineColor(runColors[ir]);
|
||||
|
||||
g_kurt_A1C0_orig[ir] = new TGraph(N_DITHER);
|
||||
g_kurt_A1C0_orig[ir]->SetTitle(runNames[ir]); g_kurt_A1C0_orig[ir]->SetMarkerStyle(24); g_kurt_A1C0_orig[ir]->SetMarkerColor(runColors[ir]); g_kurt_A1C0_orig[ir]->SetLineColor(runColors[ir]);
|
||||
|
||||
g_kurt_A1C0_hyb[ir] = new TGraph(N_DITHER);
|
||||
g_kurt_A1C0_hyb[ir]->SetTitle(runNames[ir]); g_kurt_A1C0_hyb[ir]->SetMarkerStyle(25); g_kurt_A1C0_hyb[ir]->SetMarkerColor(runColors[ir]); g_kurt_A1C0_hyb[ir]->SetLineColor(runColors[ir]);
|
||||
|
||||
// ------------------ A1C1 Initialization ------------------
|
||||
g_rms_A1C1_orig[ir] = new TGraphErrors(N_DITHER);
|
||||
g_rms_A1C1_orig[ir]->SetTitle(runNames[ir]); g_rms_A1C1_orig[ir]->SetMarkerStyle(20); g_rms_A1C1_orig[ir]->SetMarkerColor(runColors[ir]); g_rms_A1C1_orig[ir]->SetLineColor(runColors[ir]);
|
||||
|
||||
g_rms_A1C1_hyb[ir] = new TGraphErrors(N_DITHER);
|
||||
g_rms_A1C1_hyb[ir]->SetTitle(runNames[ir]); g_rms_A1C1_hyb[ir]->SetMarkerStyle(21); g_rms_A1C1_hyb[ir]->SetMarkerColor(runColors[ir]); g_rms_A1C1_hyb[ir]->SetLineColor(runColors[ir]);
|
||||
|
||||
g_fit_A1C1_orig[ir] = new TGraphErrors(N_DITHER);
|
||||
g_fit_A1C1_orig[ir]->SetTitle(runNames[ir]); g_fit_A1C1_orig[ir]->SetMarkerStyle(22); g_fit_A1C1_orig[ir]->SetMarkerColor(runColors[ir]); g_fit_A1C1_orig[ir]->SetLineColor(runColors[ir]);
|
||||
|
||||
g_fit_A1C1_hyb[ir] = new TGraphErrors(N_DITHER);
|
||||
g_fit_A1C1_hyb[ir]->SetTitle(runNames[ir]); g_fit_A1C1_hyb[ir]->SetMarkerStyle(23); g_fit_A1C1_hyb[ir]->SetMarkerColor(runColors[ir]); g_fit_A1C1_hyb[ir]->SetLineColor(runColors[ir]);
|
||||
|
||||
g_kurt_A1C1_orig[ir] = new TGraph(N_DITHER);
|
||||
g_kurt_A1C1_orig[ir]->SetTitle(runNames[ir]); g_kurt_A1C1_orig[ir]->SetMarkerStyle(24); g_kurt_A1C1_orig[ir]->SetMarkerColor(runColors[ir]); g_kurt_A1C1_orig[ir]->SetLineColor(runColors[ir]);
|
||||
|
||||
g_kurt_A1C1_hyb[ir] = new TGraph(N_DITHER);
|
||||
g_kurt_A1C1_hyb[ir]->SetTitle(runNames[ir]); g_kurt_A1C1_hyb[ir]->SetMarkerStyle(25); g_kurt_A1C1_hyb[ir]->SetMarkerColor(runColors[ir]); g_kurt_A1C1_hyb[ir]->SetLineColor(runColors[ir]);
|
||||
}
|
||||
|
||||
// =========================================================================
|
||||
|
|
@ -77,10 +106,13 @@ void plot_dither_scan()
|
|||
continue;
|
||||
}
|
||||
|
||||
auto fetchStats = [&](const char *tag, TGraphErrors *graph_rms, TGraph *graph_kurt)
|
||||
auto fetchStats = [&](const char *tag, TGraphErrors *graph_rms, TGraphErrors *graph_fit, TGraph *graph_kurt)
|
||||
{
|
||||
TH1D *h_sx3 = (TH1D *)f->Get(TString::Format("Benchmark_SX3/Benchmark_SX3_PCZ-sx3pczguess_%s", tag));
|
||||
TH1D *h_qqq = (TH1D *)f->Get(TString::Format("Benchmark_QQQ/Benchmark_QQQ_PCZ-qqqpczguess_%s", tag));
|
||||
TString sx3_path = TString::Format("Benchmark_SX3_ref/Benchmark_SX3_PCZ_%s_minus_ref", tag);
|
||||
TString qqq_path = TString::Format("Benchmark_QQQ_ref/Benchmark_QQQ_PCZ_%s_minus_ref", tag);
|
||||
|
||||
TH1D *h_sx3 = (TH1D *)f->Get(sx3_path);
|
||||
TH1D *h_qqq = (TH1D *)f->Get(qqq_path);
|
||||
|
||||
TH1D *h_combo = nullptr;
|
||||
if (h_sx3 && h_qqq) {
|
||||
|
|
@ -94,21 +126,42 @@ void plot_dither_scan()
|
|||
|
||||
if (h_combo && h_combo->GetEntries() > 50)
|
||||
{
|
||||
// Use raw statistical properties instead of a Gaussian fit
|
||||
// 1. Raw Stats
|
||||
double rms = h_combo->GetStdDev();
|
||||
double rmsErr = h_combo->GetStdDevError();
|
||||
double kurtosis = h_combo->GetKurtosis(); // < 0 means bifurcated, ~0 means Gaussian
|
||||
double kurtosis = h_combo->GetKurtosis();
|
||||
|
||||
// 2. Gaussian Fit Stats
|
||||
h_combo->Fit("gaus", "Q0"); // Q = Quiet, 0 = Do not draw
|
||||
TF1 *gausFit = h_combo->GetFunction("gaus");
|
||||
double fitSigma = 0.0;
|
||||
double fitSigmaErr = 0.0;
|
||||
|
||||
if (gausFit) {
|
||||
fitSigma = gausFit->GetParameter(2); // Param 2 is Sigma
|
||||
fitSigmaErr = gausFit->GetParError(2);
|
||||
}
|
||||
|
||||
// 3. Fill Graphs
|
||||
graph_rms->SetPoint(id, dithers[id], rms);
|
||||
graph_rms->SetPointError(id, 0.0, rmsErr);
|
||||
|
||||
graph_fit->SetPoint(id, dithers[id], fitSigma);
|
||||
graph_fit->SetPointError(id, 0.0, fitSigmaErr);
|
||||
|
||||
graph_kurt->SetPoint(id, dithers[id], kurtosis);
|
||||
|
||||
delete h_combo;
|
||||
}
|
||||
};
|
||||
|
||||
fetchStats("A1C0", g_rms_A1C0[ir], g_kurt_A1C0[ir]);
|
||||
fetchStats("A1C1", g_rms_A1C1[ir], g_kurt_A1C1[ir]);
|
||||
// Fetch Normal A1C0 vs Hybrid A1C0
|
||||
fetchStats("A1C0", g_rms_A1C0_orig[ir], g_fit_A1C0_orig[ir], g_kurt_A1C0_orig[ir]);
|
||||
fetchStats("A1C0_Hyb", g_rms_A1C0_hyb[ir], g_fit_A1C0_hyb[ir], g_kurt_A1C0_hyb[ir]);
|
||||
|
||||
// Fetch Normal A1C1 vs Hybrid A1C1
|
||||
fetchStats("A1C1", g_rms_A1C1_orig[ir], g_fit_A1C1_orig[ir], g_kurt_A1C1_orig[ir]);
|
||||
fetchStats("A1C1_Hyb", g_rms_A1C1_hyb[ir], g_fit_A1C1_hyb[ir], g_kurt_A1C1_hyb[ir]);
|
||||
|
||||
f->Close();
|
||||
delete f;
|
||||
|
|
@ -116,73 +169,132 @@ void plot_dither_scan()
|
|||
}
|
||||
|
||||
// =========================================================================
|
||||
// --- Plotting ---
|
||||
// --- CANVAS 1: A1C0 Plotting ---
|
||||
// =========================================================================
|
||||
TCanvas *c1 = new TCanvas("c1", "Dither Optimization", 1400, 1000);
|
||||
c1->Divide(2, 2);
|
||||
TCanvas *c1 = new TCanvas("c1", "A1C0 Normal vs Hybrid", 1400, 1400);
|
||||
c1->Divide(2, 3);
|
||||
|
||||
// -- Pad 1: A1C0 RMS
|
||||
c1->cd(1);
|
||||
gPad->SetGrid(1, 1);
|
||||
TMultiGraph *mg_A1C0_rms = new TMultiGraph();
|
||||
mg_A1C0_rms->SetTitle("A1C0 Emulation Resolution; Dither Sigma (mm); Raw RMS Spread (mm)");
|
||||
TLegend *leg0 = new TLegend(0.15, 0.65, 0.45, 0.85);
|
||||
leg0->SetBorderSize(1);
|
||||
// ROW 1: RAW RMS SPREAD
|
||||
c1->cd(1); gPad->SetGrid(1, 1);
|
||||
TMultiGraph *mg_A1C0_orig_rms = new TMultiGraph();
|
||||
mg_A1C0_orig_rms->SetTitle("Original A1C0 (Pure Z-Dither); Dither Sigma (mm); PC Z Raw RMS (mm)");
|
||||
TLegend *leg0 = new TLegend(0.15, 0.65, 0.45, 0.85); leg0->SetBorderSize(1);
|
||||
for (int ir = 0; ir < N_RUNS; ir++) {
|
||||
if (g_rms_A1C0[ir]->GetN() > 0) {
|
||||
mg_A1C0_rms->Add(g_rms_A1C0[ir], "PL");
|
||||
leg0->AddEntry(g_rms_A1C0[ir], runNames[ir], "pl");
|
||||
if (g_rms_A1C0_orig[ir]->GetN() > 0) {
|
||||
mg_A1C0_orig_rms->Add(g_rms_A1C0_orig[ir], "PL");
|
||||
leg0->AddEntry(g_rms_A1C0_orig[ir], runNames[ir], "pl");
|
||||
}
|
||||
}
|
||||
mg_A1C0_rms->Draw("A");
|
||||
leg0->Draw();
|
||||
mg_A1C0_orig_rms->Draw("A"); leg0->Draw();
|
||||
|
||||
// -- Pad 2: A1C1 RMS
|
||||
c1->cd(2);
|
||||
gPad->SetGrid(1, 1);
|
||||
TMultiGraph *mg_A1C1_rms = new TMultiGraph();
|
||||
mg_A1C1_rms->SetTitle("A1C1 Emulation Resolution; Dither Sigma (mm); Raw RMS Spread (mm)");
|
||||
TLegend *leg1 = new TLegend(0.15, 0.65, 0.45, 0.85);
|
||||
leg1->SetBorderSize(1);
|
||||
c1->cd(2); gPad->SetGrid(1, 1);
|
||||
TMultiGraph *mg_A1C0_hyb_rms = new TMultiGraph();
|
||||
mg_A1C0_hyb_rms->SetTitle("Hybrid A1C0 (Si Smear + Z-Dither); Dither Sigma (mm); PC Z Raw RMS (mm)");
|
||||
TLegend *leg1 = new TLegend(0.15, 0.65, 0.45, 0.85); leg1->SetBorderSize(1);
|
||||
for (int ir = 0; ir < N_RUNS; ir++) {
|
||||
if (g_rms_A1C1[ir]->GetN() > 0) {
|
||||
mg_A1C1_rms->Add(g_rms_A1C1[ir], "PL");
|
||||
leg1->AddEntry(g_rms_A1C1[ir], runNames[ir], "pl");
|
||||
if (g_rms_A1C0_hyb[ir]->GetN() > 0) {
|
||||
mg_A1C0_hyb_rms->Add(g_rms_A1C0_hyb[ir], "PL");
|
||||
leg1->AddEntry(g_rms_A1C0_hyb[ir], runNames[ir], "pl");
|
||||
}
|
||||
}
|
||||
mg_A1C1_rms->Draw("A");
|
||||
leg1->Draw();
|
||||
mg_A1C0_hyb_rms->Draw("A"); leg1->Draw();
|
||||
|
||||
// -- Pad 3: A1C0 Kurtosis
|
||||
c1->cd(3);
|
||||
gPad->SetGrid(1, 1);
|
||||
TMultiGraph *mg_A1C0_kurt = new TMultiGraph();
|
||||
mg_A1C0_kurt->SetTitle("A1C0 Bifurcation Metric; Dither Sigma (mm); Excess Kurtosis (0 = Gaussian)");
|
||||
// ROW 2: GAUSSIAN FIT SIGMA
|
||||
c1->cd(3); gPad->SetGrid(1, 1);
|
||||
TMultiGraph *mg_A1C0_orig_fit = new TMultiGraph();
|
||||
mg_A1C0_orig_fit->SetTitle("Original A1C0 Fit; Dither Sigma (mm); Gaussian Fit #sigma (mm)");
|
||||
for (int ir = 0; ir < N_RUNS; ir++) if (g_fit_A1C0_orig[ir]->GetN() > 0) mg_A1C0_orig_fit->Add(g_fit_A1C0_orig[ir], "PL");
|
||||
mg_A1C0_orig_fit->Draw("A");
|
||||
|
||||
c1->cd(4); gPad->SetGrid(1, 1);
|
||||
TMultiGraph *mg_A1C0_hyb_fit = new TMultiGraph();
|
||||
mg_A1C0_hyb_fit->SetTitle("Hybrid A1C0 Fit; Dither Sigma (mm); Gaussian Fit #sigma (mm)");
|
||||
for (int ir = 0; ir < N_RUNS; ir++) if (g_fit_A1C0_hyb[ir]->GetN() > 0) mg_A1C0_hyb_fit->Add(g_fit_A1C0_hyb[ir], "PL");
|
||||
mg_A1C0_hyb_fit->Draw("A");
|
||||
|
||||
// ROW 3: KURTOSIS
|
||||
c1->cd(5); gPad->SetGrid(1, 1);
|
||||
TMultiGraph *mg_A1C0_orig_kurt = new TMultiGraph();
|
||||
mg_A1C0_orig_kurt->SetTitle("Original A1C0 Bifurcation Metric; Dither Sigma (mm); Excess Kurtosis");
|
||||
for (int ir = 0; ir < N_RUNS; ir++) if (g_kurt_A1C0_orig[ir]->GetN() > 0) mg_A1C0_orig_kurt->Add(g_kurt_A1C0_orig[ir], "PL");
|
||||
mg_A1C0_orig_kurt->Draw("A");
|
||||
|
||||
c1->Update(); TLine *line0 = new TLine(c1->cd(5)->GetUxmin(), 0, c1->cd(5)->GetUxmax(), 0);
|
||||
line0->SetLineColor(kBlack); line0->SetLineStyle(2); line0->SetLineWidth(2); line0->Draw();
|
||||
|
||||
c1->cd(6); gPad->SetGrid(1, 1);
|
||||
TMultiGraph *mg_A1C0_hyb_kurt = new TMultiGraph();
|
||||
mg_A1C0_hyb_kurt->SetTitle("Hybrid A1C0 Bifurcation Metric; Dither Sigma (mm); Excess Kurtosis");
|
||||
for (int ir = 0; ir < N_RUNS; ir++) if (g_kurt_A1C0_hyb[ir]->GetN() > 0) mg_A1C0_hyb_kurt->Add(g_kurt_A1C0_hyb[ir], "PL");
|
||||
mg_A1C0_hyb_kurt->Draw("A");
|
||||
|
||||
c1->Update(); TLine *line1 = new TLine(c1->cd(6)->GetUxmin(), 0, c1->cd(6)->GetUxmax(), 0);
|
||||
line1->SetLineColor(kBlack); line1->SetLineStyle(2); line1->SetLineWidth(2); line1->Draw();
|
||||
|
||||
c1->SaveAs("a1c0_normal_vs_hybrid_results.png");
|
||||
|
||||
// =========================================================================
|
||||
// --- CANVAS 2: A1C1 Plotting ---
|
||||
// =========================================================================
|
||||
TCanvas *c2 = new TCanvas("c2", "A1C1 Normal vs Hybrid", 1400, 1400);
|
||||
c2->Divide(2, 3);
|
||||
|
||||
// ROW 1: RAW RMS SPREAD
|
||||
c2->cd(1); gPad->SetGrid(1, 1);
|
||||
TMultiGraph *mg_A1C1_orig_rms = new TMultiGraph();
|
||||
mg_A1C1_orig_rms->SetTitle("Original A1C1 (Pure Z-Dither); Dither Sigma (mm); PC Z Raw RMS (mm)");
|
||||
TLegend *leg2 = new TLegend(0.15, 0.65, 0.45, 0.85); leg2->SetBorderSize(1);
|
||||
for (int ir = 0; ir < N_RUNS; ir++) {
|
||||
if (g_kurt_A1C0[ir]->GetN() > 0) mg_A1C0_kurt->Add(g_kurt_A1C0[ir], "PL");
|
||||
if (g_rms_A1C1_orig[ir]->GetN() > 0) {
|
||||
mg_A1C1_orig_rms->Add(g_rms_A1C1_orig[ir], "PL");
|
||||
leg2->AddEntry(g_rms_A1C1_orig[ir], runNames[ir], "pl");
|
||||
}
|
||||
}
|
||||
mg_A1C0_kurt->Draw("A");
|
||||
mg_A1C1_orig_rms->Draw("A"); leg2->Draw();
|
||||
|
||||
// Draw a horizontal line at y=0 to easily spot where it becomes Gaussian
|
||||
c1->Update();
|
||||
TLine *line0 = new TLine(c1->cd(3)->GetUxmin(), 0, c1->cd(3)->GetUxmax(), 0);
|
||||
line0->SetLineColor(kBlack); line0->SetLineStyle(2); line0->SetLineWidth(2);
|
||||
line0->Draw();
|
||||
|
||||
// -- Pad 4: A1C1 Kurtosis
|
||||
c1->cd(4);
|
||||
gPad->SetGrid(1, 1);
|
||||
TMultiGraph *mg_A1C1_kurt = new TMultiGraph();
|
||||
mg_A1C1_kurt->SetTitle("A1C1 Bifurcation Metric; Dither Sigma (mm); Excess Kurtosis (0 = Gaussian)");
|
||||
c2->cd(2); gPad->SetGrid(1, 1);
|
||||
TMultiGraph *mg_A1C1_hyb_rms = new TMultiGraph();
|
||||
mg_A1C1_hyb_rms->SetTitle("Hybrid A1C1 (Si Smear + Z-Dither); Dither Sigma (mm); PC Z Raw RMS (mm)");
|
||||
TLegend *leg3 = new TLegend(0.15, 0.65, 0.45, 0.85); leg3->SetBorderSize(1);
|
||||
for (int ir = 0; ir < N_RUNS; ir++) {
|
||||
if (g_kurt_A1C1[ir]->GetN() > 0) mg_A1C1_kurt->Add(g_kurt_A1C1[ir], "PL");
|
||||
if (g_rms_A1C1_hyb[ir]->GetN() > 0) {
|
||||
mg_A1C1_hyb_rms->Add(g_rms_A1C1_hyb[ir], "PL");
|
||||
leg3->AddEntry(g_rms_A1C1_hyb[ir], runNames[ir], "pl");
|
||||
}
|
||||
}
|
||||
mg_A1C1_kurt->Draw("A");
|
||||
mg_A1C1_hyb_rms->Draw("A"); leg3->Draw();
|
||||
|
||||
c1->Update();
|
||||
TLine *line1 = new TLine(c1->cd(4)->GetUxmin(), 0, c1->cd(4)->GetUxmax(), 0);
|
||||
line1->SetLineColor(kBlack); line1->SetLineStyle(2); line1->SetLineWidth(2);
|
||||
line1->Draw();
|
||||
// ROW 2: GAUSSIAN FIT SIGMA
|
||||
c2->cd(3); gPad->SetGrid(1, 1);
|
||||
TMultiGraph *mg_A1C1_orig_fit = new TMultiGraph();
|
||||
mg_A1C1_orig_fit->SetTitle("Original A1C1 Fit; Dither Sigma (mm); Gaussian Fit #sigma (mm)");
|
||||
for (int ir = 0; ir < N_RUNS; ir++) if (g_fit_A1C1_orig[ir]->GetN() > 0) mg_A1C1_orig_fit->Add(g_fit_A1C1_orig[ir], "PL");
|
||||
mg_A1C1_orig_fit->Draw("A");
|
||||
|
||||
c1->SaveAs("dither_optimization_results.png");
|
||||
c2->cd(4); gPad->SetGrid(1, 1);
|
||||
TMultiGraph *mg_A1C1_hyb_fit = new TMultiGraph();
|
||||
mg_A1C1_hyb_fit->SetTitle("Hybrid A1C1 Fit; Dither Sigma (mm); Gaussian Fit #sigma (mm)");
|
||||
for (int ir = 0; ir < N_RUNS; ir++) if (g_fit_A1C1_hyb[ir]->GetN() > 0) mg_A1C1_hyb_fit->Add(g_fit_A1C1_hyb[ir], "PL");
|
||||
mg_A1C1_hyb_fit->Draw("A");
|
||||
|
||||
// ROW 3: KURTOSIS
|
||||
c2->cd(5); gPad->SetGrid(1, 1);
|
||||
TMultiGraph *mg_A1C1_orig_kurt = new TMultiGraph();
|
||||
mg_A1C1_orig_kurt->SetTitle("Original A1C1 Bifurcation Metric; Dither Sigma (mm); Excess Kurtosis");
|
||||
for (int ir = 0; ir < N_RUNS; ir++) if (g_kurt_A1C1_orig[ir]->GetN() > 0) mg_A1C1_orig_kurt->Add(g_kurt_A1C1_orig[ir], "PL");
|
||||
mg_A1C1_orig_kurt->Draw("A");
|
||||
|
||||
c2->Update(); TLine *line2 = new TLine(c2->cd(5)->GetUxmin(), 0, c2->cd(5)->GetUxmax(), 0);
|
||||
line2->SetLineColor(kBlack); line2->SetLineStyle(2); line2->SetLineWidth(2); line2->Draw();
|
||||
|
||||
c2->cd(6); gPad->SetGrid(1, 1);
|
||||
TMultiGraph *mg_A1C1_hyb_kurt = new TMultiGraph();
|
||||
mg_A1C1_hyb_kurt->SetTitle("Hybrid A1C1 Bifurcation Metric; Dither Sigma (mm); Excess Kurtosis");
|
||||
for (int ir = 0; ir < N_RUNS; ir++) if (g_kurt_A1C1_hyb[ir]->GetN() > 0) mg_A1C1_hyb_kurt->Add(g_kurt_A1C1_hyb[ir], "PL");
|
||||
mg_A1C1_hyb_kurt->Draw("A");
|
||||
|
||||
c2->Update(); TLine *line3 = new TLine(c2->cd(6)->GetUxmin(), 0, c2->cd(6)->GetUxmax(), 0);
|
||||
line3->SetLineColor(kBlack); line3->SetLineStyle(2); line3->SetLineWidth(2); line3->Draw();
|
||||
|
||||
c2->SaveAs("a1c1_normal_vs_hybrid_results.png");
|
||||
}
|
||||
|
|
@ -7,6 +7,7 @@ SCAN_DIR="dither_scan"
|
|||
|
||||
echo "=== Compiling TrackRecon.C once... ==="
|
||||
root -q -l -b -e '.L TrackRecon.C++O' 2>/dev/null
|
||||
rm -rf dither_scan/*
|
||||
|
||||
process_run() {
|
||||
local wrun=$(printf "%03d" "$1")
|
||||
|
|
@ -19,7 +20,7 @@ process_run() {
|
|||
# Pass the variable to C++
|
||||
export DITHER_SIGMA="$dither"
|
||||
|
||||
root -q -l -b -x "../../ANASEN_analysis/data/${DATASET}_Data/${PREFIX}${wrun}_mapped.root" \
|
||||
root -q -l -b -x "../ANASEN_analysis/data/${DATASET}_Data/${PREFIX}${wrun}_mapped.root" \
|
||||
-e "tree->Process(\"TrackRecon.C+\", \"${out}\")" > /dev/null 2>&1
|
||||
|
||||
[ -f "$out" ] && echo " run $wrun dither $dither OK" || echo " run $wrun dither $dither FAILED"
|
||||
|
|
@ -29,7 +30,8 @@ export -f process_run
|
|||
export SCAN_DIR
|
||||
|
||||
# Loop through the dither amounts (in mm) you want to test
|
||||
for dither in 0.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0; do
|
||||
# for dither in 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 12.0 13.0 14.0 15.0 16.0 17.0 18.0 19.0 20.0 22.0 23.0 24.0 25.0 26.0 27.0 28.0 29.0 30.0 32.0 33.0 34.0 35.0 36.0 37.0 38.0; do
|
||||
for dither in 8.0 12.0 16.0 20.0 24.0 28.0 32.0 36.0 40.0; do
|
||||
echo "=== Scanning Dither Sigma = ${dither} mm ==="
|
||||
|
||||
# 27Al alpha+gas runs (9, 12)
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user