modified: TrackRecon.C added fine timing correction

modified:   run_17F.sh updating to new format and including some other changes also updated the eloss files for 4%
This commit is contained in:
Vignesh Sitaraman 2026-10-23 10:32:40 -04:00
parent d4856783e7
commit 8be662d0d1
8 changed files with 163538 additions and 189237 deletions

View File

@ -429,6 +429,7 @@ double mass_4He = 4.002603254;
double mass_1H = 1.007825032;
double mass_30Si = 29.973770;
double mass_17F = 17.002095;
double mass_16O = 15.994915;
double mass_20Ne = 19.992440;
double mass_2H = 2.014101778; // deuteron, for (a,d) ejectile kinematics
// Recoil masses for the (a,X) ejectile channels (from the MakeVertex branch).
@ -512,7 +513,21 @@ bool pcA1C0PerCluster = true; // set in Begin() from PC_A1C0_PER_CLUSTER
std::vector<char> aClusterHasCathode;
bool fillPooledTier = false;
int qqqPairRule = 0;
int qqqSprayNch = 0;
int sx3SprayNch = 0;
double qqqFbTol = 0.15;
double qqqRingWedgeDtMax = 0.0;
double siEcut17F_qqq = 2.0;
double siEcut17F_sx3 = 1.5;
bool qqqPcGate17F = true;
bool pcGmNoWrap = false;
int miscChRF = 3;
int miscChMCP = 4;
inline bool isAnodeWireExcluded(int wire)
{
return excludeBadAnodeWires && badAnodeWires.count(wire) > 0;
@ -731,6 +746,25 @@ void TrackRecon::Begin(TTree * /*tree*/)
<< " (ratio " << (1 - qqqFbTol) / (1 + qqqFbTol) << " to "
<< (1 + qqqFbTol) / (1 - qqqFbTol) << ")" << std::endl;
}
if (getenv("QQQ_RW_DTMAX"))
qqqRingWedgeDtMax = std::atof(getenv("QQQ_RW_DTMAX"));
if (getenv("SI_ECUT_17F_QQQ"))
siEcut17F_qqq = std::atof(getenv("SI_ECUT_17F_QQQ"));
if (getenv("SI_ECUT_17F_SX3"))
siEcut17F_sx3 = std::atof(getenv("SI_ECUT_17F_SX3"));
if (getenv("QQQ_PC_GATE_17F"))
qqqPcGate17F = std::atoi(getenv("QQQ_PC_GATE_17F")) != 0;
if (getenv("PC_GM_NOWRAP"))
pcGmNoWrap = std::atoi(getenv("PC_GM_NOWRAP")) != 0;
if (getenv("MISC_SWAP_RF_MCP") && std::atoi(getenv("MISC_SWAP_RF_MCP")) != 0)
{
miscChRF = 4;
miscChMCP = 3;
}
std::cout << "QQQ ring/wedge dt max " << qqqRingWedgeDtMax << " ns (0 = off); 17F Si thresholds "
<< siEcut17F_qqq << "/" << siEcut17F_sx3 << " MeV; 17F QQQ-PC pre-gate " << qqqPcGate17F
<< "; PC GM no-wrap " << pcGmNoWrap << "; misc RF/MCP = ch" << miscChRF << "/ch" << miscChMCP
<< std::endl;
if (getenv("QQQ_PAIR_RULE"))
{
std::string r = getenv("QQQ_PAIR_RULE");
@ -738,6 +772,13 @@ void TrackRecon::Begin(TTree * /*tree*/)
: 0;
std::cout << "QQQ pairing rule: " << r << " (" << qqqPairRule << ")" << std::endl;
}
if (getenv("QQQ_SPRAY_NCH"))
qqqSprayNch = std::atoi(getenv("QQQ_SPRAY_NCH"));
if (getenv("SX3_SPRAY_NCH"))
sx3SprayNch = std::atoi(getenv("SX3_SPRAY_NCH"));
if (qqqSprayNch > 0 || sx3SprayNch > 0)
std::cout << "Spray veto: drop a detector with >= " << qqqSprayNch << " QQQ / "
<< sx3SprayNch << " SX3 hits in an event (0 = off)" << std::endl;
if (getenv("A1C1_LOWBAND_RFACTOR"))
a1c1_lowband_rfactor = std::atof(getenv("A1C1_LOWBAND_RFACTOR"));
if (getenv("A1C1_Z_SCALE_QQQ"))
@ -1319,6 +1360,30 @@ Bool_t TrackRecon::Process(Long64_t entry)
TRandom3 &rnd_qqq = anasenRandom; // shared stream (see anasenRandom)
TRandom3 &rnd_sx3 = anasenRandom;
// Spray veto (see qqqSprayNch): hits per detector this event, before any pairing.
bool qqqSpray[8] = {};
bool sx3Spray[24] = {};
if (qqqSprayNch > 0)
{
int n[8] = {};
for (int i = 0; i < qqq.multi; i++)
if (qqq.id[i] < 8 && qqq.index[i] != 112)
++n[qqq.id[i]];
for (int d = 0; d < 8; d++)
if ((qqqSpray[d] = (n[d] >= qqqSprayNch)))
plotter->Fill1D("SprayVeto_qqq_det", 8, 0, 8, d, "SprayVeto");
}
if (sx3SprayNch > 0)
{
int n[24] = {};
for (int i = 0; i < sx3.multi; i++)
if (sx3.id[i] < 24)
++n[sx3.id[i]];
for (int d = 0; d < 24; d++)
if ((sx3Spray[d] = (n[d] >= sx3SprayNch)))
plotter->Fill1D("SprayVeto_sx3_det", 24, 0, 24, d, "SprayVeto");
}
std::vector<Event> SX3_Events;
if (sx3.multi > 1)
{
@ -1330,6 +1395,8 @@ Bool_t TrackRecon::Process(Long64_t entry)
int id = sx3.id[i];
if (id >= 12)
continue;
if (sx3Spray[id])
continue;
if (sx3.ch[i] >= 8)
{
int sx3ch = sx3.ch[i] - 8;
@ -1485,6 +1552,8 @@ Bool_t TrackRecon::Process(Long64_t entry)
}
}
}
if (qqq.id[i] < 8 && qqqSpray[qqq.id[i]])
continue; // spray veto: raw histograms above still see the hit, pairing does not
for (int j = i + 1; j < qqq.multi; j++)
{
if (qqq.id[i] == qqq.id[j])
@ -1526,6 +1595,8 @@ Bool_t TrackRecon::Process(Long64_t entry)
continue;
if (chRing + qqq.id[i] * 16 == 63)
continue;
if (qqqRingWedgeDtMax > 0.0 && TMath::Abs(tRing - tWedge) > qqqRingWedgeDtMax)
continue; // ring and wedge not from the same particle (QQQ_RW_DTMAX)
if (diagnostic_tplots)
{
@ -1767,6 +1838,12 @@ Bool_t TrackRecon::Process(Long64_t entry)
}
}
if (pcGmNoWrap)
{
const double eGM = pcSlope[pc.index[i]] * pc.e[i] + pcIntercept[pc.index[i]];
pc.e[i] = eGM <= 0.0 ? 0 : (eGM >= 65535.0 ? 65535 : static_cast<unsigned short>(eGM));
}
else
pc.e[i] = pcSlope[pc.index[i]] * pc.e[i] + pcIntercept[pc.index[i]];
if (diagnostic_eplots)
{
@ -2143,10 +2220,10 @@ Bool_t TrackRecon::Process(Long64_t entry)
// siMcpRfGate cut) inherits this value and must not redither it.
for (int j = 0; j < misc.multi; j++)
{
if (misc.ch[j] == 3)
ts_rf_event = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) + (rnd.Uniform(4.0) - 2.0);
if (misc.ch[j] == 4)
ts_mcp_event = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) + (rnd.Uniform(4.0) - 2.0);
if (misc.ch[j] == miscChRF)
ts_rf_event = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) * 4.0 / 1024.0 + (rnd.Uniform(4.0) - 2.0);
if (misc.ch[j] == miscChMCP)
ts_mcp_event = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) * 4.0 / 1024.0 + (rnd.Uniform(4.0) - 2.0);
}
if (ts_rf_event > -987654321 && ts_mcp_event > -987654321)
dt_rf_mcp_event = ts_rf_event - ts_mcp_event;
@ -2182,19 +2259,19 @@ Bool_t TrackRecon::Process(Long64_t entry)
if (misc.ch[j] == 2)
{ // Needle
plotter->Fill2D("needle_vs_qqqE", 800, 0, 16384, 800, 0, 10, misc.e[j], qqqevent.Energy1, "misc");
ts_needle = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) + (rnd.Uniform(4.0) - 2.0);
ts_needle = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) * 4.0 / 1024.0 + (rnd.Uniform(4.0) - 2.0);
found_needle = 1;
plotter->Fill1D("dt_qqq_needle", 1600, -2000, 2000, ts_qqq - ts_needle, "misc");
}
if (misc.ch[j] == 3)
if (misc.ch[j] == miscChRF)
{ // RF
ts_rf = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) + (rnd.Uniform(4.0) - 2.0);
ts_rf = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) * 4.0 / 1024.0 + (rnd.Uniform(4.0) - 2.0);
found_rf = 1;
plotter->Fill1D("dt_qqq_rf_innerring" + std::to_string(qqq_inner_ring), 1600, -2000, 2000, ts_qqq - ts_rf, "misc");
}
if (misc.ch[j] == 4)
if (misc.ch[j] == miscChMCP)
{ // mcp
ts_mcp = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) + (rnd.Uniform(4.0) - 2.0);
ts_mcp = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) * 4.0 / 1024.0 + (rnd.Uniform(4.0) - 2.0);
found_mcp = 1;
plotter->Fill1D("dt_qqq_mcp_innerring" + std::to_string(qqq_inner_ring), 1600, -2000, 2000, ts_qqq - ts_mcp, "misc");
}
@ -2230,19 +2307,19 @@ Bool_t TrackRecon::Process(Long64_t entry)
if (misc.ch[j] == 2)
{ // Needle
plotter->Fill2D("needle_vs_sx3E", 800, 0, 16384, 800, 0, 10, misc.e[j], sx3event.Energy1, "misc");
ts_needle = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) + (rnd.Uniform(4.0) - 2.0);
ts_needle = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) * 4.0 / 1024.0 + (rnd.Uniform(4.0) - 2.0);
found_needle = 1;
plotter->Fill1D("dt_sx3_needle", 1600, -2000, 2000, ts_sx3 - ts_needle, "misc");
}
if (misc.ch[j] == 3)
if (misc.ch[j] == miscChRF)
{ // RF
ts_rf = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) + (rnd.Uniform(4.0) - 2.0);
ts_rf = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) * 4.0 / 1024.0 + (rnd.Uniform(4.0) - 2.0);
found_rf = 1;
plotter->Fill1D("dt_sx3_rf", 1600, -2000, 2000, ts_sx3 - ts_rf, "misc");
}
if (misc.ch[j] == 4)
if (misc.ch[j] == miscChMCP)
{ // mcp
ts_mcp = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) + (rnd.Uniform(4.0) - 2.0);
ts_mcp = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) * 4.0 / 1024.0 + (rnd.Uniform(4.0) - 2.0);
found_mcp = 1;
plotter->Fill1D("dt_sx3_mcp", 1600, -2000, 2000, ts_sx3 - ts_mcp, "misc");
}
@ -4591,7 +4668,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
if (sievent.Energy1 < si_ecut)
continue; // coarse Si energy cut
if (isQQQ && rx == "m17Fax")
if (isQQQ && rx == "m17Fax" && qqqPcGate17F)
{
bool hasCoincidentPC = false;
for (const auto &pcevent : PC_Events)
@ -4630,8 +4707,9 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
double vertex_z = r_rhoMin_fix.Z();
const bool axisSafe = (topo1 == "a1c2fix") || (topo1 == "a2c0") || (topo2 == "a1c1_inband");
// fillBeamProfile(plotter, r_rhoMin_fix, si_pos, x2f - si_pos, "reaction_" + rx + "_" + det, axisSafe);
// if (beamPerp(r_rhoMin_fix) > perp_cut || vertex_z < z_entrance || vertex_z > 30.0)
if (vertex_z < z_entrance || vertex_z > 30.0)
if (beamPerp(r_rhoMin_fix) > perp_cut || vertex_z < z_entrance)
return;
if (rx == "m27Alax" && vertex_z > 30.0)
return;
double theta = (si_pos - r_rhoMin_fix).Theta();
@ -4943,29 +5021,30 @@ void miscHistograms_17Fax(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
double ebeam_17F_MeV = 56.72;
// 17F(a,a)/(a,d)/(a,p): ejectile + recoil masses per channel.
AAEjectileMasses ej17F{mass_4He, mass_17F, mass_2H, mass_19Ne_rec, mass_1H, mass_20Ne};
reaction_ax_core(plotter, QQQ_Events, PC_Events, aClusters, true, "m17Fax", "qqq", 2.0, 6.0, TMath::Pi() / 4.0,
// AAEjectileMasses ej16O{mass_4He, mass_16O, mass_2H, mass_19Ne_rec, mass_1H, mass_20Ne};
reaction_ax_core(plotter, QQQ_Events, PC_Events, aClusters, true, "m17Fax", "qqq", siEcut17F_qqq, 6.0, TMath::Pi() / 4.0,
30.0, 40000.0, 30.0, ebeam_17F_MeV, MeV_to_cm_17F_spl, cm_to_MeV_17F_spl, mass_17F, ej17F, globaltag, dt_rf_mcp);
reaction_ax_core(plotter, SX3_Events, PC_Events, aClusters, false, "m17Fax", "sx3", 1.5, 10.0, TMath::Pi() / 3.0,
reaction_ax_core(plotter, SX3_Events, PC_Events, aClusters, false, "m17Fax", "sx3", siEcut17F_sx3, 10.0, TMath::Pi() / 3.0,
30.0, 40000.0, 30.0, ebeam_17F_MeV, MeV_to_cm_17F_spl, cm_to_MeV_17F_spl, mass_17F, ej17F, globaltag, dt_rf_mcp);
if (siMcpRfGate1Cut)
{
reaction_ax_core(plotter, QQQ_Events, PC_Events, aClusters, true, "m17Fax", "qqq", 2.0, 6.0, TMath::Pi() / 4.0,
30.0, 40000.0, 30.0, ebeam_17F_MeV, MeV_to_cm_17F_spl, cm_to_MeV_17F_spl, mass_17F, ej17F, globaltag, dt_rf_mcp,
siMcpRfGate1Cut, ts_rf_event, ts_mcp_event, "_Gate1", "Gate1_");
reaction_ax_core(plotter, SX3_Events, PC_Events, aClusters, false, "m17Fax", "sx3", 1.5, 10.0, TMath::Pi() / 3.0,
30.0, 40000.0, 30.0, ebeam_17F_MeV, MeV_to_cm_17F_spl, cm_to_MeV_17F_spl, mass_17F, ej17F, globaltag, dt_rf_mcp,
siMcpRfGate1Cut, ts_rf_event, ts_mcp_event, "_Gate1", "Gate1_");
}
if (siMcpRfGate2Cut)
{
reaction_ax_core(plotter, QQQ_Events, PC_Events, aClusters, true, "m17Fax", "qqq", 2.0, 6.0, TMath::Pi() / 4.0,
30.0, 40000.0, 30.0, ebeam_17F_MeV, MeV_to_cm_17F_spl, cm_to_MeV_17F_spl, mass_17F, ej17F, globaltag, dt_rf_mcp,
siMcpRfGate2Cut, ts_rf_event, ts_mcp_event, "_Gate2", "Gate2_");
reaction_ax_core(plotter, SX3_Events, PC_Events, aClusters, false, "m17Fax", "sx3", 1.5, 10.0, TMath::Pi() / 3.0,
30.0, 40000.0, 30.0, ebeam_17F_MeV, MeV_to_cm_17F_spl, cm_to_MeV_17F_spl, mass_17F, ej17F, globaltag, dt_rf_mcp,
siMcpRfGate2Cut, ts_rf_event, ts_mcp_event, "_Gate2", "Gate2_");
}
// if (siMcpRfGate1Cut)
// {
// reaction_ax_core(plotter, QQQ_Events, PC_Events, aClusters, true, "m17Fax", "qqq", siEcut17F_qqq, 6.0, TMath::Pi() / 4.0,
// 30.0, 40000.0, 30.0, ebeam_17F_MeV, MeV_to_cm_17F_spl, cm_to_MeV_17F_spl, mass_17F, ej17F, globaltag, dt_rf_mcp,
// siMcpRfGate1Cut, ts_rf_event, ts_mcp_event, "_Gate1", "Gate1_");
// reaction_ax_core(plotter, SX3_Events, PC_Events, aClusters, false, "m17Fax", "sx3", siEcut17F_sx3, 10.0, TMath::Pi() / 3.0,
// 30.0, 40000.0, 30.0, ebeam_17F_MeV, MeV_to_cm_17F_spl, cm_to_MeV_17F_spl, mass_17F, ej17F, globaltag, dt_rf_mcp,
// siMcpRfGate1Cut, ts_rf_event, ts_mcp_event, "_Gate1", "Gate1_");
// }
// if (siMcpRfGate2Cut)
// {
// reaction_ax_core(plotter, QQQ_Events, PC_Events, aClusters, true, "m17Fax", "qqq", siEcut17F_qqq, 6.0, TMath::Pi() / 4.0,
// 30.0, 40000.0, 30.0, ebeam_17F_MeV, MeV_to_cm_17F_spl, cm_to_MeV_17F_spl, mass_17F, ej17F, globaltag, dt_rf_mcp,
// siMcpRfGate2Cut, ts_rf_event, ts_mcp_event, "_Gate2", "Gate2_");
// reaction_ax_core(plotter, SX3_Events, PC_Events, aClusters, false, "m17Fax", "sx3", siEcut17F_sx3, 10.0, TMath::Pi() / 3.0,
// 30.0, 40000.0, 30.0, ebeam_17F_MeV, MeV_to_cm_17F_spl, cm_to_MeV_17F_spl, mass_17F, ej17F, globaltag, dt_rf_mcp,
// siMcpRfGate2Cut, ts_rf_event, ts_mcp_event, "_Gate2", "Gate2_");
// }
}
// 27Al(a,a) excitation functions for BOTH silicon branches (QQQ + SX3), with the

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -7,8 +7,20 @@ export CO2percent=4
export pressure_in_torr=250
export CATHODE_GAIN=1
export source_vertex=-200.0
export DEDX_SCALE=0.89
export DEDX_SCALE=0.962
export CUTLIST=cuts_list.txt
# Spray veto: drop a QQQ/SX3 detector that fires this many hits in one event
# (whole-detector bursts, mostly saturated). SX3 6 is the tighter option.
export QQQ_SPRAY_NCH=7
export SX3_SPRAY_NCH=9
# Cleanups (TrackRecon defaults reproduce the old behaviour when these are unset)
export QQQ_RW_DTMAX=100 # QQQ ring/wedge time coincidence [ns]
export PC_GM_NOWRAP=1 # negative gain-matched PC energies -> 0 instead of ~65000
# To test, one at a time against a reference output:
export SI_ECUT_17F_QQQ=0.45 # was 2.0 MeV
export SI_ECUT_17F_SX3=0.6 # was 1.5 MeV
export QQQ_PC_GATE_17F=0 # drop the QQQ -450..-200 ns pre-gate
export MISC_SWAP_RF_MCP=1 # 405 ch14 = MCP, ch10 = RF; transpose simcprfgate1/2 first
run_once() {
local wrun=$(printf "%03d" "$1")
@ -42,7 +54,7 @@ rm -f ${OUT_DIR}/*.root
echo "Running Eloss.py with a scaling parameter of $DEDX_SCALE"
echo "running with a beam offset of $BEAM_AXIS_X $BEAM_AXIS_Y"
python3 eloss_calculations/Eloss.py
myenv/bin/python eloss_calculations/Eloss.py || { echo "Eloss.py failed"; exit 1; }
# Pre-compile TrackRecon.C ONCE on a single core so parallel jobs don't race on ACLiC
echo "Pre-compiling TrackRecon.C..."
@ -62,6 +74,14 @@ unset pressure_in_torr
unset CATHODE_GAIN
unset A1C1_LOWBAND_RFACTOR
unset CUTLIST
unset QQQ_SPRAY_NCH
unset SX3_SPRAY_NCH
unset QQQ_RW_DTMAX
unset PC_GM_NOWRAP
unset SI_ECUT_17F_QQQ
unset SI_ECUT_17F_SX3
unset QQQ_PC_GATE_17F
unset MISC_SWAP_RF_MCP
unset A1C1_Z_SCALE_QQQ
unset A1C1_Z_OFF_QQQ
unset A1C1_Z_OFF_SX3