rtang/ANASEN_analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Compare commits

base: rtang:38c8c7e0e0bf1457f04158c145d936ebe26915ce
Branches Tags
rtang:master
rtang:devel_vignesh
rtang:devel_james
rtang:devel_sudarshan
rtang:devel_ruchi
...
compare: rtang:83114262dba32d418f21825c343d265776380e75
Branches Tags
rtang:devel_vignesh
rtang:devel_james
rtang:devel_sudarshan
rtang:master
rtang:devel_ruchi

2 Commits
38c8c7e0e0 ... 83114262db

Author SHA1 Message Date
Vignesh Sitaraman 83114262db modified: TrackRecon.C 2026-09-09 10:49:17 -04:00
Vignesh Sitaraman 9b2c20daca modified: TrackRecon.C 
modified:   scratch/make_prettyplots.C
 2026-09-08 15:35:43 -04:00
2 changed files with 105 additions and 95 deletions
Show Stats Expand all files Collapse all files

198
TrackRecon.C
View File

@ -25,6 +25,7 @@ Int_t colors[40] = {
#include <TSpline.h> #include <TSpline.h>
#include <TSystem.h> // gSystem->mkdir for the pc_calib_raw/ output directory #include <TSystem.h> // gSystem->mkdir for the pc_calib_raw/ output directory
#include <cmath> // std::floor/fabs/lround for vertexZSliceTag()
#include <fstream> #include <fstream>
#include <iomanip> #include <iomanip>
#include <iostream> #include <iostream>
@ -60,12 +61,8 @@ bool process_alpha_proton_scattering = false,
// 1175mm instead of 1105, plus some part of the window actually lies outside the chamber // 1175mm instead of 1105, plus some part of the window actually lies outside the chamber
double source_vertex = 53.0, double source_vertex = 53.0,
z_entrance = -174.3 - 9.7 - 270.0, z_entrance = -174.3 - 9.7 - 270.0,
dither_sigma = 8.0, // single pcz dither width; A1C0/A1C1/A1C2 all use this. dither_sigma = 8.0,
// dither_sigma_c0 (=16.0, always consumed as /2.0) was cathode_gain = 1.0,
// removed -- numerically identical at the default but it
// desynced from dither_sigma the moment DITHER_SIGMA was
// set, giving the QQQ and SX3 twins a silent 2x difference.
cathode_gain = 1.0,
a1c1_cfrac_split = 0.0, a1c1_cfrac_split = 0.0,
a1c1_missing_fmax = 2.0, a1c1_missing_fmax = 2.0,
a1c1_lowband_rfactor = 0.0, a1c1_lowband_rfactor = 0.0,
@ -75,9 +72,9 @@ double source_vertex = 53.0,
a1c1_z_off_sx3 = 2.52614, a1c1_z_off_sx3 = 2.52614,
beam_axis_x = 0.0, beam_axis_x = 0.0,
beam_axis_y = 0.0, beam_axis_y = 0.0,
beam_axis_z0 = 0.0, // reference z at which (beam_axis_x, beam_axis_y) is specified beam_axis_z0 = 0.0,
beam_tilt_x = 0.0, // dx/dz slope of the beam axis beam_tilt_x = 0.0,
beam_tilt_y = 0.0, // dy/dz slope of the beam axis beam_tilt_y = 0.0,
ta_foil_z_mm = 0.0, ta_foil_z_mm = 0.0,
alpha_source_mev = 5.486; alpha_source_mev = 5.486;
@ -178,60 +175,60 @@ inline double beamPerp(const TVector3 &p)
return perp.Mag(); return perp.Mag();
} }
inline void fillBeamProfile(HistPlotter *plotter, const TVector3 &vertex, // inline void fillBeamProfile(HistPlotter *plotter, const TVector3 &vertex,
const TVector3 &si, const TVector3 &dir, const std::string &tag, // const TVector3 &si, const TVector3 &dir, const std::string &tag,
bool axisSafe) // bool axisSafe)
{ // {
const std::string folder = "BeamAxis"; // const std::string folder = "BeamAxis";
const std::string all = "beamAxis_all_"; // const std::string all = "beamAxis_all_";
const std::string one = "beamAxis_" + tag + "_"; // const std::string one = "beamAxis_" + tag + "_";
double vx = vertex.X(), vy = vertex.Y(), vz = vertex.Z(); // double vx = vertex.X(), vy = vertex.Y(), vz = vertex.Z();
double pocaDist = beamPerp(vertex); // double pocaDist = beamPerp(vertex);
plotter->Fill2D(all + "vertexX_vs_Z", 250, -450, 50, 200, -50, 50, vz, vx, folder); // plotter->Fill2D(all + "vertexX_vs_Z", 250, -450, 50, 200, -50, 50, vz, vx, folder);
plotter->Fill2D(all + "vertexY_vs_Z", 250, -450, 50, 200, -50, 50, vz, vy, folder); // plotter->Fill2D(all + "vertexY_vs_Z", 250, -450, 50, 200, -50, 50, vz, vy, folder);
plotter->Fill1D(all + "pocaDist", 400, 0, 100, pocaDist, folder); // plotter->Fill1D(all + "pocaDist", 400, 0, 100, pocaDist, folder);
// Per-branch copy so anomalies can be traced to a specific reconstruction path. // // Per-branch copy so anomalies can be traced to a specific reconstruction path.
plotter->Fill2D(one + "vertexX_vs_Z", 250, -450, 50, 200, -50, 50, vz, vx, folder); // plotter->Fill2D(one + "vertexX_vs_Z", 250, -450, 50, 200, -50, 50, vz, vx, folder);
plotter->Fill2D(one + "vertexY_vs_Z", 250, -450, 50, 200, -50, 50, vz, vy, folder); // plotter->Fill2D(one + "vertexY_vs_Z", 250, -450, 50, 200, -50, 50, vz, vy, folder);
plotter->Fill1D(one + "pocaDist", 400, 0, 100, pocaDist, folder); // plotter->Fill1D(one + "pocaDist", 400, 0, 100, pocaDist, folder);
const double zLo = -440.0, zHi = 40.0; // const double zLo = -440.0, zHi = 40.0;
const int nSlice = 16; // const int nSlice = 16;
const double sliceW = (zHi - zLo) / nSlice; // 30 mm // const double sliceW = (zHi - zLo) / nSlice; // 30 mm
if (vz >= zLo && vz < zHi)
{
int is = static_cast<int>((vz - zLo) / sliceW);
if (is >= 0 && is < nSlice)
{
char buf[8];
snprintf(buf, sizeof(buf), "%02d", is);
plotter->Fill2D(all + "vertexXY_z" + buf, 200, -50, 50, 200, -50, 50, vx, vy, folder);
}
}
double dz = dir.Z(); // if (vz >= zLo && vz < zHi)
double vzBin = beamVertexNominal(si, dir).Z(); // {
if (axisSafe && TMath::Abs(dz) > 1e-6 && vzBin >= zLo && vzBin < zHi) // int is = static_cast<int>((vz - zLo) / sliceW);
{ // if (is >= 0 && is < nSlice)
int k = static_cast<int>((vzBin - zLo) / sliceW); // {
if (k >= 0 && k < nSlice) // char buf[8];
{ // snprintf(buf, sizeof(buf), "%03d", is);
double zPlane = zLo + (k + 0.5) * sliceW; // plotter->Fill2D(all + "vertexXY_z" + buf, 200, -50, 50, 200, -50, 50, vx, vy, folder);
double s = (zPlane - si.Z()) / dz; // }
double cx = si.X() + s * dir.X(); // }
double cy = si.Y() + s * dir.Y();
if (TMath::Abs(cx) < 100.0 && TMath::Abs(cy) < 100.0) // double dz = dir.Z();
{ // double vzBin = beamVertexNominal(si, dir).Z();
char kbuf[8]; // if (axisSafe && TMath::Abs(dz) > 1e-6 && vzBin >= zLo && vzBin < zHi)
snprintf(kbuf, sizeof(kbuf), "%02d", k); // {
plotter->Fill1D(all + "crossX_z" + kbuf, 800, -100, 100, cx, folder); // int k = static_cast<int>((vzBin - zLo) / sliceW);
plotter->Fill1D(all + "crossY_z" + kbuf, 800, -100, 100, cy, folder); // if (k >= 0 && k < nSlice)
plotter->Fill1D(one + "crossX_z" + kbuf, 800, -100, 100, cx, folder); // {
plotter->Fill1D(one + "crossY_z" + kbuf, 800, -100, 100, cy, folder); // double zPlane = zLo + (k + 0.5) * sliceW;
} // double s = (zPlane - si.Z()) / dz;
} // double cx = si.X() + s * dir.X();
} // double cy = si.Y() + s * dir.Y();
} // if (TMath::Abs(cx) < 100.0 && TMath::Abs(cy) < 100.0)
// {
// char kbuf[8];
// snprintf(kbuf, sizeof(kbuf), "%03d", k);
// plotter->Fill1D(all + "crossX_z" + kbuf, 800, -100, 100, cx, folder);
// plotter->Fill1D(all + "crossY_z" + kbuf, 800, -100, 100, cy, folder);
// plotter->Fill1D(one + "crossX_z" + kbuf, 800, -100, 100, cx, folder);
// plotter->Fill1D(one + "crossY_z" + kbuf, 800, -100, 100, cy, folder);
// }
// }
// }
// }
struct PCPath struct PCPath
{ {
@ -1300,7 +1297,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
if (det.valid && diagnostic_eplots) if (det.valid && diagnostic_eplots)
{ {
// std::cout << det.frontEL << " " << det.frontEL*sx3RightGain[id][det.stripF] << std::endl; // std::cout << det.frontEL << " " << det.frontEL*sx3RightGain[id][det.stripF] << std::endl;
plotter->Fill2D("be_vs_x_sx3_id_"+std::to_string(id)+"_f"+std::to_string(det.stripF)+"_b"+std::to_string(det.stripB),200,-1,1,800,0,8192,det.frontX,det.backE,"evsx"); plotter->Fill2D("be_vs_x_sx3_id_" + std::to_string(id) + "_f" + std::to_string(det.stripF) + "_b" + std::to_string(det.stripB), 200, -1, 1, 800, 0, 8192, det.frontX, det.backE, "evsx");
plotter->Fill2D("unmatched_be_vs_x_sx3_id_" + std::to_string(id), 200, -1, 1, 800, 0, 4096, det.frontX, det.backE, "evsx"); plotter->Fill2D("unmatched_be_vs_x_sx3_id_" + std::to_string(id), 200, -1, 1, 800, 0, 4096, det.frontX, det.backE, "evsx");
plotter->Fill2D("unmatched_be_vs_x_sx3", 200, -1, 1, 800, 0, 4096, det.frontX, det.backE, "evsx"); plotter->Fill2D("unmatched_be_vs_x_sx3", 200, -1, 1, 800, 0, 4096, det.frontX, det.backE, "evsx");
plotter->Fill2D("matched_be_vs_x_sx3", 200, -60, 60, 800, 0, 8192, det.frontX * sx3FrontGain[id][det.stripF] + sx3FrontOffset[id][det.stripF], det.backE * sx3BackGain[id][det.stripF][det.stripB], "evsx"); plotter->Fill2D("matched_be_vs_x_sx3", 200, -60, 60, 800, 0, 8192, det.frontX * sx3FrontGain[id][det.stripF] + sx3FrontOffset[id][det.stripF], det.backE * sx3BackGain[id][det.stripF][det.stripB], "evsx");
@ -1308,7 +1305,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
plotter->Fill2D("matched_be_vs_x_sx3_id_" + std::to_string(id) + "_f" + std::to_string(det.stripF), 200, -60, 60, 800, 0, 8192, plotter->Fill2D("matched_be_vs_x_sx3_id_" + std::to_string(id) + "_f" + std::to_string(det.stripF), 200, -60, 60, 800, 0, 8192,
det.frontX * sx3FrontGain[id][det.stripF] + sx3FrontOffset[id][det.stripF], det.backE * sx3BackGain[id][det.stripF][det.stripB], "evsx_matched"); det.frontX * sx3FrontGain[id][det.stripF] + sx3FrontOffset[id][det.stripF], det.backE * sx3BackGain[id][det.stripF][det.stripB], "evsx_matched");
plotter->Fill2D("fe_vs_x_sx3_id_"+std::to_string(id)+"_f"+std::to_string(det.stripF)+"_"+std::to_string(det.stripB),200,-1,1,800,0,4096,det.frontX,det.backE,"evsx"); plotter->Fill2D("fe_vs_x_sx3_id_" + std::to_string(id) + "_f" + std::to_string(det.stripF) + "_" + std::to_string(det.stripB), 200, -1, 1, 800, 0, 4096, det.frontX, det.backE, "evsx");
plotter->Fill2D("l_vs_r_sx3_id_" + std::to_string(id) + "_f" + std::to_string(det.stripF), 800, 0, 4096, 800, 0, 4096, det.frontEL, det.frontER, "l_vs_r"); plotter->Fill2D("l_vs_r_sx3_id_" + std::to_string(id) + "_f" + std::to_string(det.stripF), 800, 0, 4096, 800, 0, 4096, det.frontEL, det.frontER, "l_vs_r");
} }
if (det.valid && (id == 9 || id == 7 || id == 1 || id == 3) && det.stripF != DEFAULT_NULL && det.stripB != DEFAULT_NULL) if (det.valid && (id == 9 || id == 7 || id == 1 || id == 3) && det.stripF != DEFAULT_NULL && det.stripB != DEFAULT_NULL)
@ -1335,7 +1332,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
plotter->Fill2D("sx3backs_gm", 100, 0, 100, 800, 0, 8192, det.stripB + 4 * id, backE, "hCalSX3"); plotter->Fill2D("sx3backs_gm", 100, 0, 100, 800, 0, 8192, det.stripB + 4 * id, backE, "hCalSX3");
plotter->Fill1D("sx3backs_calib", 800, 0, 8192, backE, "hCalSX3"); plotter->Fill1D("sx3backs_calib", 800, 0, 8192, backE, "hCalSX3");
plotter->Fill2D("SX3CartesianPlot", 200, -100, 100, 200, -100, 100, 88.0*TMath::Cos(phi_n),88.0*TMath::Sin(phi_n), "hCalSX3"); plotter->Fill2D("SX3CartesianPlot", 200, -100, 100, 200, -100, 100, 88.0 * TMath::Cos(phi_n), 88.0 * TMath::Sin(phi_n), "hCalSX3");
plotter->Fill2D("SX3CartesianPlot" + std::to_string(id), 200, -100, 100, 200, -100, 100, rho_at_strip * TMath::Cos(phi_n), rho_at_strip * TMath::Sin(phi_n), "hCalSX3"); plotter->Fill2D("SX3CartesianPlot" + std::to_string(id), 200, -100, 100, 200, -100, 100, rho_at_strip * TMath::Cos(phi_n), rho_at_strip * TMath::Sin(phi_n), "hCalSX3");
} }
if (diagnostic_tplots) if (diagnostic_tplots)
@ -1862,12 +1859,19 @@ Bool_t TrackRecon::Process(Long64_t entry)
double dt_rf_mcp_event = -987654321; double dt_rf_mcp_event = -987654321;
double ts_rf_event = -987654321, ts_mcp_event = -987654321; double ts_rf_event = -987654321, ts_mcp_event = -987654321;
{ {
// +/-2 ns clock dither, same dequantization as tRing/tWedge (QQQ, +/-8 ns)
// and Fsx3.ts (SX3, +/-8 ns) -- misc.t/misc.tf are raw digitizer ticks, not
// a continuous CFD fraction, so left undithered this stair-steps just like
// the Si times did before their own dither was added. Applied once here,
// at construction, same as tRing/tWedge -- everything downstream (the
// miscHistograms_17Fax parameters, dt_si_mcp/dt_si_rf in reaction_ax_core's
// siMcpRfGate cut) inherits this value and must not redither it.
for (int j = 0; j < misc.multi; j++) for (int j = 0; j < misc.multi; j++)
{ {
if (misc.ch[j] == 3) if (misc.ch[j] == 3)
ts_rf_event = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]); 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) if (misc.ch[j] == 4)
ts_mcp_event = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]); ts_mcp_event = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) + (rnd.Uniform(4.0) - 2.0);
} }
if (ts_rf_event > -987654321 && ts_mcp_event > -987654321) if (ts_rf_event > -987654321 && ts_mcp_event > -987654321)
dt_rf_mcp_event = ts_rf_event - ts_mcp_event; dt_rf_mcp_event = ts_rf_event - ts_mcp_event;
@ -1903,33 +1907,35 @@ Bool_t TrackRecon::Process(Long64_t entry)
if (misc.ch[j] == 2) if (misc.ch[j] == 2)
{ // Needle { // Needle
plotter->Fill2D("needle_vs_qqqE", 800, 0, 16384, 800, 0, 10, misc.e[j], qqqevent.Energy1, "misc"); 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]); ts_needle = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) + (rnd.Uniform(4.0) - 2.0);
found_needle = 1; found_needle = 1;
plotter->Fill1D("dt_qqq_needle", 800, -2000, 2000, ts_qqq - ts_needle, "misc"); plotter->Fill1D("dt_qqq_needle", 1600, -2000, 2000, ts_qqq - ts_needle, "misc");
} }
if (misc.ch[j] == 3) if (misc.ch[j] == 3)
{ // RF { // RF
ts_rf = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]); ts_rf = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) + (rnd.Uniform(4.0) - 2.0);
found_rf = 1; found_rf = 1;
plotter->Fill1D("dt_qqq_rf_innerring" + std::to_string(qqq_inner_ring), 800, -2000, 2000, ts_qqq - ts_rf, "misc"); 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] == 4)
{ // mcp { // mcp
ts_mcp = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]); ts_mcp = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) + (rnd.Uniform(4.0) - 2.0);
found_mcp = 1; found_mcp = 1;
plotter->Fill1D("dt_qqq_mcp_innerring" + std::to_string(qqq_inner_ring), 800, -2000, 2000, ts_qqq - ts_mcp, "misc"); plotter->Fill1D("dt_qqq_mcp_innerring" + std::to_string(qqq_inner_ring), 1600, -2000, 2000, ts_qqq - ts_mcp, "misc");
} }
} }
if (found_rf && found_mcp) if (found_rf && found_mcp)
{ {
if (ctr == 0) if (ctr == 0)
plotter->Fill1D("dt_rf_mcp_qqq_innerring" + std::to_string(qqq_inner_ring), 500, -1000, 1000, ts_rf - ts_mcp, "misc"); plotter->Fill1D("dt_rf_mcp_qqq_innerring" + std::to_string(qqq_inner_ring), 1000, -1000, 1000, ts_rf - ts_mcp, "misc");
double dt_rf_mcp = ts_rf - ts_mcp; double dt_rf_mcp = ts_rf - ts_mcp;
double dt_qqq_rf = ts_qqq - ts_rf; double dt_qqq_rf = ts_qqq - ts_rf;
double dt_qqq_mcp = ts_qqq - ts_mcp; double dt_qqq_mcp = ts_qqq - ts_mcp;
plotter->Fill2D("dt(qqq,rf)_vs_(rf,mcp)_innerring" + std::to_string(qqq_inner_ring), 800, -2000, 2000, 640, -2000, 2000, dt_qqq_rf, dt_rf_mcp, "misc"); plotter->Fill2D("dt(qqq,rf)_vs_(rf,mcp)_innerring" + std::to_string(qqq_inner_ring), 1600, -2000, 2000, 1280, -2000, 2000, dt_qqq_rf, dt_rf_mcp, "misc");
plotter->Fill2D("dt_(qqq,mcp)_vs_(qqq,rf)_innerring" + std::to_string(qqq_inner_ring), 800, -1400, 2000, 640, -2000, 2000, dt_qqq_mcp, dt_qqq_rf, "misc"); plotter->Fill2D("dt_(qqq,mcp)_vs_(qqq,rf)_innerring" + std::to_string(qqq_inner_ring), 1600, -1400, 2000, 1280, -2000, 2000, dt_qqq_mcp, dt_qqq_rf, "misc");
plotter->Fill2D("dt_(qqq,mcp)_vs_(rf,mcp)_innerring" + std::to_string(qqq_inner_ring), 1000, -1400, 1000, 640, -2000, 2000, dt_qqq_mcp, dt_rf_mcp, "misc"); plotter->Fill2D("dt_(qqq,mcp)_vs_(rf,mcp)_innerring" + std::to_string(qqq_inner_ring), 2000, -1400, 1000, 1280, -2000, 2000, dt_qqq_mcp, dt_rf_mcp, "misc");
plotter->Fill2D("dt_(qqq,mcp)_vs_(rf,mcp)", 2000, -1400, 1000, 1280, -2000, 2000, dt_qqq_mcp, dt_rf_mcp, "misc");
plotter->Fill2D("dt_(si,mcp)_vs_(rf,mcp)", 2000, -1400, 1000, 1280, -2000, 2000, dt_qqq_mcp, dt_rf_mcp, "misc");
} }
ctr += 1; ctr += 1;
} }
@ -1949,33 +1955,34 @@ Bool_t TrackRecon::Process(Long64_t entry)
if (misc.ch[j] == 2) if (misc.ch[j] == 2)
{ // Needle { // Needle
plotter->Fill2D("needle_vs_sx3E", 800, 0, 16384, 800, 0, 10, misc.e[j], sx3event.Energy1, "misc"); 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]); ts_needle = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) + (rnd.Uniform(4.0) - 2.0);
found_needle = 1; found_needle = 1;
plotter->Fill1D("dt_sx3_needle", 800, -2000, 2000, ts_sx3 - ts_needle, "misc"); plotter->Fill1D("dt_sx3_needle", 1600, -2000, 2000, ts_sx3 - ts_needle, "misc");
} }
if (misc.ch[j] == 3) if (misc.ch[j] == 3)
{ // RF { // RF
ts_rf = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]); ts_rf = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) + (rnd.Uniform(4.0) - 2.0);
found_rf = 1; found_rf = 1;
plotter->Fill1D("dt_sx3_rf", 800, -2000, 2000, ts_sx3 - ts_rf, "misc"); plotter->Fill1D("dt_sx3_rf", 1600, -2000, 2000, ts_sx3 - ts_rf, "misc");
} }
if (misc.ch[j] == 4) if (misc.ch[j] == 4)
{ // mcp { // mcp
ts_mcp = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]); ts_mcp = static_cast<double>(misc.t[j]) + static_cast<double>(misc.tf[j]) + (rnd.Uniform(4.0) - 2.0);
found_mcp = 1; found_mcp = 1;
plotter->Fill1D("dt_sx3_mcp", 800, -2000, 2000, ts_sx3 - ts_mcp, "misc"); plotter->Fill1D("dt_sx3_mcp", 1600, -2000, 2000, ts_sx3 - ts_mcp, "misc");
} }
} }
if (found_rf && found_mcp) if (found_rf && found_mcp)
{ {
if (ctr == 0) if (ctr == 0)
plotter->Fill1D("dt_rf_mcp_sx3", 500, -1000, 1000, ts_rf - ts_mcp, "misc"); plotter->Fill1D("dt_rf_mcp_sx3", 1000, -1000, 1000, ts_rf - ts_mcp, "misc");
double dt_rf_mcp = ts_rf - ts_mcp; double dt_rf_mcp = ts_rf - ts_mcp;
double dt_sx3_rf = ts_sx3 - ts_rf; double dt_sx3_rf = ts_sx3 - ts_rf;
double dt_sx3_mcp = ts_sx3 - ts_mcp; double dt_sx3_mcp = ts_sx3 - ts_mcp;
plotter->Fill2D("dt(sx3,rf)_vs_(rf,mcp)", 800, -2000, 2000, 640, -2000, 2000, dt_sx3_rf, dt_rf_mcp, "misc"); plotter->Fill2D("dt(sx3,rf)_vs_(rf,mcp)", 1600, -2000, 2000, 1280, -2000, 2000, dt_sx3_rf, dt_rf_mcp, "misc");
plotter->Fill2D("dt_(sx3,mcp)_vs_(sx3,rf)", 800, -1400, 2000, 640, -2000, 2000, dt_sx3_mcp, dt_sx3_rf, "misc"); plotter->Fill2D("dt_(sx3,mcp)_vs_(sx3,rf)", 1600, -1400, 2000, 1280, -2000, 2000, dt_sx3_mcp, dt_sx3_rf, "misc");
plotter->Fill2D("dt_(sx3,mcp)_vs_(rf,mcp)", 1000, -1400, 1000, 640, -2000, 2000, dt_sx3_mcp, dt_rf_mcp, "misc"); plotter->Fill2D("dt_(sx3,mcp)_vs_(rf,mcp)", 2000, -1400, 1000, 1280, -2000, 2000, dt_sx3_mcp, dt_rf_mcp, "misc");
plotter->Fill2D("dt_(si,mcp)_vs_(rf,mcp)", 2000, -1400, 1000, 1280, -2000, 2000, dt_sx3_mcp, dt_rf_mcp, "misc");
} }
ctr += 1; ctr += 1;
} }
@ -2146,8 +2153,7 @@ void protonAlphaHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_E
// purely longitudinal direction. // purely longitudinal direction.
TVector3 r_rhoMin_fix = beamVertex(x1, v); TVector3 r_rhoMin_fix = beamVertex(x1, v);
double vertex_z = r_rhoMin_fix.Z(); double vertex_z = r_rhoMin_fix.Z();
fillBeamProfile(plotter, r_rhoMin_fix, x1, v, "apCoinc", // fillBeamProfile(plotter, r_rhoMin_fix, x1, v, "apCoinc", pcevent.multi1 == 1 && pcevent.multi2 == 2);
pcevent.multi1 == 1 && pcevent.multi2 == 2);
double theta_q = (qqqevent.pos - beamAxisPoint(vertex_z)).Theta(); double theta_q = (qqqevent.pos - beamAxisPoint(vertex_z)).Theta();
// double theta_q = (qqqevent.pos - r_rhoMin_fix).Theta(); // double theta_q = (qqqevent.pos - r_rhoMin_fix).Theta();
double sinTheta_customV = TMath::Sin(theta_q); double sinTheta_customV = TMath::Sin(theta_q);
@ -2588,7 +2594,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
// beam axis at (0,0), silently ignoring BEAM_AXIS_X/Y, and had no guard for a // beam axis at (0,0), silently ignoring BEAM_AXIS_X/Y, and had no guard for a
// purely longitudinal direction. // purely longitudinal direction.
TVector3 r_rhoMin_fix = beamVertex(x1, v); TVector3 r_rhoMin_fix = beamVertex(x1, v);
fillBeamProfile(plotter, r_rhoMin_fix, x1, v, "sx3a1c2", true); // a1c2 only, gated above // fillBeamProfile(plotter, r_rhoMin_fix, x1, v, "sx3a1c2", true); // a1c2 only, gated above
plotter->Fill1D("VertexRecon_pczfix_sx3", 800, -300, 300, r_rhoMin_fix.Z(), "Vertex_Reconstruction"); plotter->Fill1D("VertexRecon_pczfix_sx3", 800, -300, 300, r_rhoMin_fix.Z(), "Vertex_Reconstruction");
plotter->Fill1D("VertexRecon_pczfix", 800, -300, 300, r_rhoMin_fix.Z(), "Vertex_Reconstruction"); plotter->Fill1D("VertexRecon_pczfix", 800, -300, 300, r_rhoMin_fix.Z(), "Vertex_Reconstruction");
plotter->Fill1D("pczfix_A1C2_1d_sx3", 600, -200, 200, pcz_fix, "PCZ_Recon"); plotter->Fill1D("pczfix_A1C2_1d_sx3", 600, -200, 200, pcz_fix, "PCZ_Recon");
@ -3140,7 +3146,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
// beam axis at (0,0), silently ignoring BEAM_AXIS_X/Y, and had no guard for a // beam axis at (0,0), silently ignoring BEAM_AXIS_X/Y, and had no guard for a
// purely longitudinal direction. // purely longitudinal direction.
TVector3 r_rhoMin_fix = beamVertex(x1, v); TVector3 r_rhoMin_fix = beamVertex(x1, v);
fillBeamProfile(plotter, r_rhoMin_fix, x1, v, "qqqa1c2", true); // a1c2 only, gated above // fillBeamProfile(plotter, r_rhoMin_fix, x1, v, "qqqa1c2", true); // a1c2 only, gated above
double sinTheta_customV = TMath::Sin((qqqevent.pos - beamAxisPoint(r_rhoMin_fix.Z())).Theta()); double sinTheta_customV = TMath::Sin((qqqevent.pos - beamAxisPoint(r_rhoMin_fix.Z())).Theta());
plotter->Fill2D("dE3_E_CathodeQQQR_A1C2_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2 * sinTheta_customV, "PID_dE_E"); plotter->Fill2D("dE3_E_CathodeQQQR_A1C2_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2 * sinTheta_customV, "PID_dE_E");
@ -3699,7 +3705,7 @@ void protonAlphaElastic_core(HistPlotter *plotter, const std::vector<Event> &Si_
TVector3 r_rhoMin_fix = beamVertex(sievent.pos, x2f - sievent.pos); TVector3 r_rhoMin_fix = beamVertex(sievent.pos, x2f - sievent.pos);
double vertex_z = r_rhoMin_fix.Z(); double vertex_z = r_rhoMin_fix.Z();
const bool axisSafe = (multi2 == 2) || (multi1 == 2 && multi2 == 0); const bool axisSafe = (multi2 == 2) || (multi1 == 2 && multi2 == 0);
fillBeamProfile(plotter, r_rhoMin_fix, sievent.pos, x2f - sievent.pos, "elastic_" + det, axisSafe); // fillBeamProfile(plotter, r_rhoMin_fix, sievent.pos, x2f - sievent.pos, "elastic_" + det, axisSafe);
if (vertex_z < z_entrance) if (vertex_z < z_entrance)
return; return;
double theta = (sievent.pos - r_rhoMin_fix).Theta(); double theta = (sievent.pos - r_rhoMin_fix).Theta();
@ -3738,6 +3744,8 @@ void protonAlphaElastic_core(HistPlotter *plotter, const std::vector<Event> &Si_
plotter->Fill2D(rx + "_Ef_vs_theta" + ejtag + sfx, 100, 0, 180, 800, 0, 10, theta * 180 / M_PI, Efix, pmlabel); plotter->Fill2D(rx + "_Ef_vs_theta" + ejtag + sfx, 100, 0, 180, 800, 0, 10, theta * 180 / M_PI, Efix, pmlabel);
plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + sfx, 720, 0, 180, 800, -10, 10, theta * 180 / M_PI, Ex, pmlabel); plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + sfx, 720, 0, 180, 800, -10, 10, theta * 180 / M_PI, Ex, pmlabel);
plotter->Fill2D(rx + "_Ex_vs_phi" + ejtag + sfx, 180, -180, 180, 800, -10, 10, sievent.pos.Phi() * 180 / M_PI, Ex, pmlabel); plotter->Fill2D(rx + "_Ex_vs_phi" + ejtag + sfx, 180, -180, 180, 800, -10, 10, sievent.pos.Phi() * 180 / M_PI, Ex, pmlabel);
plotter->Fill2D(rx + "_Ex_vs_X" + ejtag + sfx, 800, -200, 200, 800, -10, 10, sievent.pos.X(), Ex, pmlabel);
plotter->Fill2D(rx + "_Ex_vs_Y" + ejtag + sfx, 800, -200, 200, 800, -10, 10, sievent.pos.Y(), Ex, pmlabel);
for (const auto &pcevent : PC_Events) for (const auto &pcevent : PC_Events)
{ {
@ -4078,8 +4086,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
TVector3 r_rhoMin_fix = beamVertex(sievent.pos, x2f - sievent.pos); TVector3 r_rhoMin_fix = beamVertex(sievent.pos, x2f - sievent.pos);
double vertex_z = r_rhoMin_fix.Z(); double vertex_z = r_rhoMin_fix.Z();
const bool axisSafe = (topo1 == "a1c2fix") || (topo1 == "a2c0") || (topo2 == "a1c1_inband"); const bool axisSafe = (topo1 == "a1c2fix") || (topo1 == "a2c0") || (topo2 == "a1c1_inband");
fillBeamProfile(plotter, r_rhoMin_fix, sievent.pos, x2f - sievent.pos, // fillBeamProfile(plotter, r_rhoMin_fix, sievent.pos, x2f - sievent.pos, "reaction_" + rx + "_" + det, axisSafe);
"reaction_" + rx + "_" + det, axisSafe);
if (beamPerp(r_rhoMin_fix) > perp_cut || vertex_z < z_entrance || vertex_z > 30.0) if (beamPerp(r_rhoMin_fix) > perp_cut || vertex_z < z_entrance || vertex_z > 30.0)
return; return;
@ -4138,6 +4145,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
std::string t = topo.empty() ? "" : ("_" + topo); std::string t = topo.empty() ? "" : ("_" + topo);
plotter->Fill1D(rx + "_Ex_from" + ejtag + t + sfx, 600, -10, 20, Ex, pmlabel); plotter->Fill1D(rx + "_Ex_from" + ejtag + t + sfx, 600, -10, 20, Ex, pmlabel);
plotter->Fill2D(rx + "_VertexReconZ_vs_Ef" + ejtag + t + sfx, 800, -400, 400, 800, 0, ef_max, vertex_z, Efix, pmlabel); plotter->Fill2D(rx + "_VertexReconZ_vs_Ef" + ejtag + t + sfx, 800, -400, 400, 800, 0, ef_max, vertex_z, Efix, pmlabel);
plotter->Fill2D(rx + "_VertexReconZ_vs_E" + ejtag + sfx, 800, -400, 400, 800, 0, ef_max, vertex_z, sievent.Energy1, pmlabel);
plotter->Fill2D(rx + "_VertexReconZ_vs_Ex" + ejtag + t + sfx, 800, -400, 400, 600, -10, 20, vertex_z, Ex, pmlabel); plotter->Fill2D(rx + "_VertexReconZ_vs_Ex" + ejtag + t + sfx, 800, -400, 400, 600, -10, 20, vertex_z, Ex, pmlabel);
plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + t + sfx, 720, 0, 180, 600, -10, 20, theta * 180 / M_PI, Ex, pmlabel); plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + t + sfx, 720, 0, 180, 600, -10, 20, theta * 180 / M_PI, Ex, pmlabel);
@ -4181,6 +4189,8 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
plotter->Fill2D(rx + "_Ef_vs_theta" + ejtag + sfx, 100, 0, 180, 800, 0, ef_max, theta * 180 / M_PI, Efix, 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_theta" + ejtag + sfx, 720, 0, 180, 800, -20, 20, theta * 180 / M_PI, Ex, pmlabel); plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + sfx, 720, 0, 180, 800, -20, 20, theta * 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); plotter->Fill2D(rx + "_Ex_vs_phi" + ejtag + sfx, 45, -180, 180, 600, -10, 20, phi * 180 / M_PI, Ex, pmlabel);
plotter->Fill2D(rx + "_Ex_vs_X" + ejtag + sfx, 100, -100, 120, 800, -6, 15, sievent.pos.X(), Ex, pmlabel);
plotter->Fill2D(rx + "_Ex_vs_Y" + ejtag + sfx, 100, -100, 120, 800, -6, 15, sievent.pos.Y(), Ex, pmlabel);
for (const auto &pcevent : PC_Events) for (const auto &pcevent : PC_Events)
{ {

2
scratch/make_prettyplots.C
View File

@ -167,7 +167,7 @@ void make_prettyplots(const char *rootFile,
h->GetXaxis()->CenterTitle(true); h->GetXaxis()->CenterTitle(true);
h->GetYaxis()->CenterTitle(true); h->GetYaxis()->CenterTitle(true);
gPad->SetLogz(); // gPad->SetLogz();
h->Draw("colz"); h->Draw("colz");
} }