deleted: Armory/mass20.txt

modified:   TrackRecon.C paap gating
	modified:   run_17F.sh
	modified:   run_27Al.sh
	modified:   run_tr.sh
	modified:   scratch/overlay_2d.C
This commit is contained in:
Vignesh Sitaraman 2026-09-04 10:16:34 -04:00
parent 78add5432a
commit 80b079b81d
6 changed files with 80 additions and 3639 deletions

File diff suppressed because it is too large Load Diff

View File

@ -45,7 +45,7 @@ bool process_alpha_proton_scattering = false,
doPCSX3ClusterAnalysis = true,
doPCQQQClusterAnalysis = true,
doOldAnalysis = false,
BenchMark = false,
BenchMark = true,
onewire_analysis = true,
diagnostic_eplots = false,
diagnostic_tplots = true,
@ -3605,7 +3605,15 @@ void protonAlphaElastic_core(HistPlotter *plotter, const std::vector<Event> &Si_
auto fillHypothesis = [&](bool alphaHyp)
{
const std::string ejtag = (alphaHyp && sievent.Energy1 < 5) ? "_a" : "_p";
// const std::string ejtag = (alphaHyp && sievent.Energy1 < 5) ? "_a" : "_p";
std::string ejtag;
if (alphaHyp && sievent.Energy1 < 5)
ejtag = "_a";
else if (sievent.Energy1 >= 6.4 && sievent.Energy1 < 7.0)
ejtag = "_p";
else
ejtag = "_maybep";
std::string pmlabel = misclabel + ejtag;
TSpline3 *ej_fwd = alphaHyp ? MeV_to_cm_spl : MeV_to_cm_p_spl;
TSpline3 *ej_inv = alphaHyp ? cm_to_MeV_spl : cm_to_MeVp_spl;
@ -3684,9 +3692,9 @@ void protonAlphaElastic_core(HistPlotter *plotter, const std::vector<Event> &Si_
}
};
if (pid != SiPcPid::kAlpha)
if (pid != SiPcPid::kAlpha && sievent.Energy1 >= 6.4 && sievent.Energy1 < 7.0)
fillHypothesis(false); // proton, or PID unavailable (legacy default)
if (pid == SiPcPid::kAlpha)
if (pid == SiPcPid::kAlpha && sievent.Energy1 < 6.2)
fillHypothesis(true);
};
@ -3964,8 +3972,6 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
plotter->Fill2D(rx + "_BeamEnergy_vs_VertexZ" + sfx, 800, -400, 400, 400, 0, beamE0, vertex_z, beam_energy_at_vertex, folderPrefix + globaltag + "_" + rx + "+misc_" + det);
bool trueProton = (beam_energy_at_vertex < 10.0);
double ex_as_proton = 0.0, ex_as_alpha = 0.0;
auto fillHypothesis = [&](double m3, double m4, TSpline3 *ej_fwd, TSpline3 *ej_inv, const std::string &ejtag)
{
@ -4016,28 +4022,25 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
if (ebeam_kin_MeV > 0.0)
plotter->Fill2D(rx + "_BeamEnergy_ETrack_vs_EKin" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5,
beam_energy_at_vertex, ebeam_kin_MeV, pmlabel);
if (ejtag == "_p")
if (beam_energy_at_vertex > 4.0 && beam_energy_at_vertex <= 12.0)
{
plotter->Fill2D(rx + "_EKin-ETrack2235keV_vs_phi" + ejtag + t + sfx, 45, -180, 180, 600, -20, 40, phi * 180 / M_PI, ebeam_kin_2235keV - beam_energy_at_vertex, pmlabel);
}
if (ejtag == "_p" && rx == "m27Alax")
{
// if (beam_energy_at_vertex < 4.0 && beam_energy_at_vertex > 0.1)
plotter->Fill2D(rx + "_ETrack_vs_EKinGS" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5,
beam_energy_at_vertex, ebeam_kin_GS, folderPrefix + "ETrackvsKin_assumed");
// else if (beam_energy_at_vertex <= 12.0)
plotter->Fill2D(rx + "_ETrack_vs_EKin2235keV" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5,
beam_energy_at_vertex, ebeam_kin_2235keV, folderPrefix + "ETrackvsKin_assumed");
// else if (beam_energy_at_vertex < 24.0)
{
plotter->Fill2D(rx + "_ETrack_vs_EKin3498keV" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5,
beam_energy_at_vertex, ebeam_kin_3498keV, folderPrefix + "ETrackvsKin_assumed");
plotter->Fill2D(rx + "_ETrack_vs_EKin3774keV" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5,
beam_energy_at_vertex, ebeam_kin_3774keV, folderPrefix + "ETrackvsKin_assumed");
}
// else if (beam_energy_at_vertex <= 36.00)
plotter->Fill2D(rx + "_ETrack_vs_EKin4809keV" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5,
beam_energy_at_vertex, ebeam_kin_4809keV, folderPrefix + "ETrackvsKin_assumed");
// else if (beam_energy_at_vertex < 42.0)
plotter->Fill2D(rx + "_ETrack_vs_EKin5614keV" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5,
beam_energy_at_vertex, ebeam_kin_5614keV, folderPrefix + "ETrackvsKin_assumed");
// else
plotter->Fill2D(rx + "_ETrack_vs_EKin6550keV" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5,
beam_energy_at_vertex, ebeam_kin_6550keV, folderPrefix + "ETrackvsKin_assumed");
}
@ -4045,7 +4048,6 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
plotter->Fill2D(rx + "_dE_E_Anode" + sfx, 400, 0, dEa_max, 800, 0, 40000, sievent.Energy1, anodeE, pmlabel);
plotter->Fill2D(rx + "_dE_E_Anode" + sfx + "_E<10MeV" + std::to_string(beam_energy_at_vertex < 10), 400, 0, dEa_max, 800, 0, 40000, sievent.Energy1, anodeE, pmlabel);
if (cathodeE >= 0.0)
{
plotter->Fill2D(rx + "_dE_E_Cathode" + sfx, 400, 0, dEa_max, 800, 0, dEc_max, sievent.Energy1, cathodeE, pmlabel);
@ -4055,7 +4057,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
}
plotter->Fill1D(rx + "_pczfix" + sfx, 600, -300, 300, pcz_fix, pmlabel);
plotter->Fill2D(rx + "_Ef_vs_theta" + ejtag + sfx, 100, 0, 180, 800, 0, ef_max, theta * 180 / M_PI, Efix, pmlabel);
plotter->Fill2D(rx + "_Ex_vs_phi" + ejtag + sfx, 180, -180, 180, 600, -10, 20, phi * 180 / M_PI, Ex, pmlabel);
plotter->Fill2D(rx + "_Ex_vs_phi" + ejtag + sfx, 45, -180, 180, 600, -10, 20, phi * 180 / M_PI, Ex, pmlabel);
for (const auto &pcevent : PC_Events)
{
@ -4065,6 +4067,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
continue;
plotter->Fill2D(rx + "_Ex_vs_dT" + ejtag + sfx, 500, -2000, 2000, 600, -10, 20, (sievent.Time1 - pcevent.Time1), Ex, pmlabel);
plotter->Fill2D(rx + "_dEgasCalib_vs_dT" + ejtag + sfx, 500, -2000, 2000, 800, 0, 0.6, (sievent.Time1 - pcevent.Time1), anodeE_MeV, pmlabel);
plotter->Fill2D(rx + "_dEgasCalibCathode_vs_dT" + ejtag + sfx, 500, -2000, 2000, 800, 0, 0.6, (sievent.Time1 - pcevent.Time1), cathodeE_MeV, pmlabel);
}
if (dt_rf_mcp > -900000000)
@ -4072,8 +4075,6 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
plotter->Fill1D(rx + "_VertexReconZ" + sfx, 800, -400, 400, vertex_z, pmlabel);
forEachTier(topo1, topo2, methodGroup, plot_with_tag);
if (trueProton)
plot_with_tag("trueProton"); // clean, alpha-free proton sub-sample
// Gas segmentation validation
PCCollect pcc = pcCollectionPath(r_rhoMin_fix, sievent.pos);
@ -4095,7 +4096,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
plotter->Fill2D(rx + "_dEgasRaw_vs_VertexZ" + ejtag + sfx, 800, -400, 400, 800, 0, 20000, vertex_z, anodeE, pmlabel);
plotter->Fill2D(rx + "_dEgasRaw_vs_theta" + ejtag + sfx, 180, 0, 180, 800, 0, 20000, theta * 180 / M_PI, anodeE, pmlabel);
plotter->Fill2D(rx + "_dEgasCalib_vs_theta" + ejtag + sfx, 360, 0, 180, 800, 0, 0.6, theta * 180 / M_PI, anodeE_MeV, pmlabel);
plotter->Fill2D(rx + "_dEgasCalib_vs_phi" + ejtag + sfx, 90, -200, 200, 800, 0, 0.6, phi * 180 / M_PI, anodeE_MeV, pmlabel);
plotter->Fill2D(rx + "_dEgasCalib_vs_phi" + ejtag + sfx, 90, -180, 180, 800, 0, 0.6, phi * 180 / M_PI, anodeE_MeV, pmlabel);
// if (anodeCh >= 0)
// plotter->Fill2D(rx + "_dEgasCalib_vs_E" + ejtag + sfx + "_anode" + pad2(anodeCh), 400, 0, ef_max, 800, 0, 0.6, sievent.Energy1, anodeE_MeV, pmlabel);
plotter->Fill2D(rx + "_dEgasCalib_vs_Ex" + ejtag + sfx, 600, -10, 20, 800, 0, 0.6, Ex, anodeE_MeV, pmlabel);

View File

@ -35,7 +35,7 @@ export -f run_once
export DATASET="17F"
export PREFIX="Run_"
export OUT_DIR="Output_17F"
export source_vertex=-57.28
export source_vertex=-200.0
export CATHODE_GAIN=1
# A1C1 cfrac low-band fold + z autocal (see TrackRecon.C Begin()). Defaults are the
# 17F values; override here to re-tune without recompiling.

View File

@ -23,8 +23,8 @@ process_run() {
local infile="../ANASEN_analysis/data/${DATASET}_Data/${prefix}${wrun}_mapped.root"
# Dynamically point to the correct output directory for this X/Y iteration
# local current_out_dir="Output_27Al_X${BEAM_AXIS_X}_Y${BEAM_AXIS_Y}"
local current_out_dir="Output_27Al"
# local current_out_dir="Output_27Al_X${BEAM_AXIS_X}_Y${BEAM_AXIS_Y}"
local out="${current_out_dir}/results_run${wrun}.root"
root -q -l -b -x "$infile" \
@ -47,9 +47,9 @@ export -f process_run
# BEAM_AXIS_Y=$y
# Define and create a clean directory name BEFORE running parallel tasks
# CURRENT_OUT_DIR="Output_27Al_X${BEAM_AXIS_X}_Y${BEAM_AXIS_Y}"
CURRENT_OUT_DIR="Output_27Al"
# CURRENT_OUT_DIR="Output_27Al_run_48"
# CURRENT_OUT_DIR="Output_27Al_X${BEAM_AXIS_X}_Y${BEAM_AXIS_Y}"
rm -f ${OUT_DIR}/*.root
mkdir -p "$CURRENT_OUT_DIR"
echo "Running Eloss.py with a scaling parameter of $DEDX_SCALE"
@ -57,9 +57,11 @@ export -f process_run
python3 eloss_calculations/Eloss.py
echo "Starting parallel processing..."
time parallel --bar -j 10 process_run ::: {24..41}
time parallel --bar -j 10 process_run ::: 44 45 46 {50..59}
# time parallel --bar -j 10 process_run ::: 48 # pc without coincidence
time parallel --bar -j 12 process_run ::: {24..41} 44 45 46 {50..59}
# time parallel --bar -j 10 process_run ::: {24..41}
# time parallel --bar -j 10 process_run ::: 44 45 46 {50..59}
time parallel --bar -j 1 process_run ::: 48 # pc without coincidence
mv "${CURRENT_OUT_DIR}/results_run048.root" "Output_27Al_run48/."
# time parallel --bar -j 4 process_run ::: 62 63 66 67 68
# time parallel --bar -j 1 process_run ::: 73
# time parallel --bar -j 1 process_run ::: 74
@ -67,7 +69,7 @@ export -f process_run
echo "Merging files..."
# Fixed: Safely merge using the clean directory variable (added -f to overwrite if re-running)
hadd -k -f -j 4 "${CURRENT_OUT_DIR}/output_27Al.root" "${CURRENT_OUT_DIR}/results_run"*.root
hadd -k -f -j 4 "${CURRENT_OUT_DIR}/Output_27Al.root" "${CURRENT_OUT_DIR}/results_run"*.root
# done
# done

View File

@ -70,7 +70,7 @@ if [[ 1 -eq 1 ]]; then
unset CATHODE_GAIN
unset timecut_low
unset pressure_in_torr
exit
# exit
fi
# --- Block 4: 17F Alpha+Gas Runs (18-21) ---
@ -102,8 +102,8 @@ if [[ 1 -eq 1 ]]; then
echo "Starting parallel processing for 27Al proton runs..."
# process_run 18
parallel --bar -j 8 process_run ::: 15 {17..22}
# parallel --bar -j 8 process_run ::: {17..22}
# parallel --bar -j 8 process_run ::: 15 {17..22}
parallel --bar -j 8 process_run ::: {17..22}
hadd -j 4 -k ${OUT_DIR}/Al_protons.root ${OUT_DIR}/results_run0{15..22}.root
unset CATHODE_GAIN
# exit

View File

@ -48,6 +48,7 @@
// Call, now just: rootFile, axis labels, and the drawing knobs --
//
// root -l -b -q 'scratch/overlay_2d.C("Output_27Al/output_27Al.root", "Tracked Beam Energy (MeV)", "Kinematic Energy (MeV)", 2, 2, 2.0)'
// root -l -b -q 'scratch/overlay_2d.C("Output_27Al/output_27Al.root", "Tracked Beam Energy (MeV)", "Kinematic Energy (MeV)")'
#include "TFile.h"
#include "TH2.h"
@ -59,7 +60,7 @@
#include <vector>
void overlay_2d(TString rootFile, TString xAxisLabel, TString yAxisLabel,
int rebinX = 8, int rebinY = 8, double minz = 5.0)
int rebinX = 4, int rebinY = 4, double minz = 8.0)
{
gROOT->SetStyle("Plain");
gStyle->SetOptStat(0);
@ -67,7 +68,7 @@ void overlay_2d(TString rootFile, TString xAxisLabel, TString yAxisLabel,
// ---- Hardcoded selection: edit these to choose what to plot ----
TString reaction = "m27Alax"; // reaction tag, e.g. "m27Alax", "m17Fax"
TString ejtag = "p"; // ejectile: "p" (proton) or "a" (alpha)
TString topology = "a2c0"; // e.g. "a1c0", "a1c1", "a1c2fix", "a2c0", "a1c1c2"
TString topology = "a1c1_inband"; // e.g. "a1c0", "a1c1", "a1c2fix", "a2c0", "a1c1c2"
TString detector = "sx3"; // "sx3" or "qqq"
// Excited states to overlay, in the same order as labels[]/xlo[]/xhi[]
@ -75,14 +76,16 @@ void overlay_2d(TString rootFile, TString xAxisLabel, TString yAxisLabel,
// note 3498's tag is lowercase "kev" (TrackRecon.C:4028), unlike every
// other state's "keV". 3774 keV is deliberately left out (too close
// to 3498, spin-parity predicts it's weaker anyway).
TString states[] = {"GS", "2235keV", "3498kev", "4809keV", "5614keV", "6550keV"};
TString states[] = {"GS", "2235keV", "3498keV", "4810keV", "6550keV"};
TString labels[] = {"Ground State", "2.235 MeV", "3.498 MeV",
"4.809 MeV", "5.614 MeV", "6.550 MeV"};
"4.810 MeV", "6.550 MeV"};
// Per-state display window (MeV), same order as states[] above.
// See the header comment for why these are currently a no-op.
double xlo[] = {0, 4, 12, 24, 36, 42};
double xhi[] = {4, 12, 24, 36, 42, 90};
// double xlo[] = {0, 4, 12, 24, 36};
// double xhi[] = {4, 12, 24, 36, 90};
double xlo[] = {0, 4, 16, 28, 36};
double xhi[] = {4, 16, 28, 36, 90};
const size_t nStates = sizeof(states) / sizeof(states[0]);
if (sizeof(labels) / sizeof(labels[0]) != nStates ||
@ -94,8 +97,18 @@ void overlay_2d(TString rootFile, TString xAxisLabel, TString yAxisLabel,
return;
}
// Standard ANASEN color sequence
int colors[] = {kBlack, kRed + 1, kAzure + 2, kGreen + 2, kMagenta + 1};
// Custom 12-color set (named ROOT colors + hex via TColor::GetColor), not the
// ANASEN colors[40] sequence in TrackRecon.C -- 12 distinct colors is plenty of
// headroom for however many states get overlaid here without collisions.
int colors[] = {kBlack, kRed + 1, kAzure + 2, kGreen + 2, kMagenta + 1,
TColor::GetColor("#3f90da"),
TColor::GetColor("#bc1e00"),
kGreen + 2,
TColor::GetColor("#832db5"),
kPink - 7,
TColor::GetColor("#e76300"),
TColor::GetColor("#b8ab6f")};
const int nColors = sizeof(colors) / sizeof(colors[0]);
std::vector<TH2 *> hists;
std::vector<TString> usedLabels;
@ -144,6 +157,12 @@ void overlay_2d(TString rootFile, TString xAxisLabel, TString yAxisLabel,
hists[i]->SetMinimum(minz);
}
double globalMax = 0;
for (size_t i = 0; i < hists.size(); i++)
globalMax = std::max(globalMax, hists[i]->GetMaximum());
for (size_t i = 0; i < hists.size(); i++)
hists[i]->SetMaximum(globalMax);
TCanvas *c = new TCanvas("c", "", 1800, 1800);
c->SetLeftMargin(0.15);
c->SetBottomMargin(0.15);
@ -152,7 +171,15 @@ void overlay_2d(TString rootFile, TString xAxisLabel, TString yAxisLabel,
c->SetGridx(1);
c->SetGridy(1);
TLegend *leg = new TLegend(0.65, 0.75, 0.9, 0.9);
// Box height scales with the number of series instead of a fixed constant --
// at SetTextSize(0.03) a fixed 0.15 NDC-tall box only gives each of 6 rows
// 0.025 NDC, less than the text height itself, which is what made the legend
// look crowded. 0.045 NDC/row leaves real padding and stays correctly sized
// for however many states are overlaid.
const double legEntryHeight = 0.045;
const double legY2 = 0.9;
double legY1 = legY2 - legEntryHeight * hists.size();
TLegend *leg = new TLegend(0.65, legY1, 0.9, legY2);
leg->SetBorderSize(1);
leg->SetFillStyle(1001); // Solid background so grid doesn't bleed through
leg->SetFillColor(kWhite);
@ -160,7 +187,12 @@ void overlay_2d(TString rootFile, TString xAxisLabel, TString yAxisLabel,
for (size_t i = 0; i < hists.size(); i++)
{
hists[i]->SetFillColor(colors[i % 5]);
int color = colors[i % nColors];
hists[i]->SetMarkerColor(color);
hists[i]->SetLineColor(color);
hists[i]->SetFillColor(color);
hists[i]->SetFillStyle(1001); // solid fill for the "box" boxes below --
// "Plain" style's default fill is hollow
if (i == 0)
{
hists[i]->SetTitle("");