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

modified: Calibration.C corrected the p energy being used to the 241Am peak

Browse Source 
modified:   GainMatchQQQ.C made all output files .dat instead of .txt
	modified:   TrackRecon.C made ome more plots for dissertation
	modified:   qqq_Calib.dat updated with 241Am peak position
	modified:   run_17F.sh
	modified:   run_tr.sh
	modified:   scratch/make_prettyplots.C made change so now we can input multiple runs at once and plots them together
This commit is contained in:
Vignesh Sitaraman 2026-06-07 17:18:44 -04:00
parent 9ad4e9d5d4
commit 702399be84
7 changed files with 1149 additions and 971 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
Calibration.C
View File

@ -37,7 +37,7 @@ void Calibration::Begin(TTree * /*tree*/)
plotter = new HistPlotter("Calib.root", "TFILE"); plotter = new HistPlotter("Calib.root", "TFILE");
// ----------------------- Load QQQ Gains // ----------------------- Load QQQ Gains
{ {
std::string filename = "qqq_GainMatch.txt"; std::string filename = "qqq_GainMatch.dat";
std::ifstream infile(filename); std::ifstream infile(filename);
if (!infile.is_open()) if (!infile.is_open())
{ {
@ -164,10 +164,10 @@ void Calibration::Terminate()
double calibArray[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}}; double calibArray[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
bool calibValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}}; bool calibValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
std::ofstream outFile("qqq_Calib.txt"); std::ofstream outFile("qqq_Calib.dat");
if (!outFile.is_open()) if (!outFile.is_open())
{ {
std::cerr << "Error opening qqq_Calib.txt!" << std::endl; std::cerr << "Error opening qqq_Calib.dat!" << std::endl;
return; return;
} }
@ -216,8 +216,8 @@ void Calibration::Terminate()
if (adcPeak <= 0) if (adcPeak <= 0)
continue; continue;
// double slope_keV = AM241_PEAK / adcPeak; // keV per ADC double slope_keV = AM241_PEAK / adcPeak; // keV per ADC
double slope_keV = P_PEAK / adcPeak; // keV per ADC // double slope_keV = P_PEAK / adcPeak; // keV per ADC
calibArray[det][ring][wedge] = slope_keV; calibArray[det][ring][wedge] = slope_keV;
calibValid[det][ring][wedge] = true; calibValid[det][ring][wedge] = true;
@ -229,7 +229,7 @@ void Calibration::Terminate()
} }
outFile.close(); outFile.close();
std::cout << "Wrote QQQ calibration file qqq_Calib.txt\n"; std::cout << "Wrote QQQ calibration file qqq_Calib.dat\n";
//---------------------------------------------------------------------- //----------------------------------------------------------------------
// 3. Build fully calibrated 2D combined histogram // 3. Build fully calibrated 2D combined histogram

4
GainMatchQQQ.C
View File

@ -154,10 +154,10 @@ void GainMatchQQQ::Terminate()
bool gainValid[MAX_DET][MAX_RING][MAX_WEDGE] = {{{false}}}; bool gainValid[MAX_DET][MAX_RING][MAX_WEDGE] = {{{false}}};
// Output file for the Minimizer // Output file for the Minimizer
std::ofstream outFile("qqq_GainMatch.txt"); std::ofstream outFile("qqq_GainMatch.dat");
// Benchmark/Debug file // Benchmark/Debug file
std::ofstream benchFile("benchmark_diff.txt"); std::ofstream benchFile("benchmark_diff.dat");
benchFile << "ID Wedge Ring Chi2NDF Slope SlopeErr" << std::endl; benchFile << "ID Wedge Ring Chi2NDF Slope SlopeErr" << std::endl;
if (!outFile.is_open()) { std::cerr << "Error opening output file!" << std::endl; return; } if (!outFile.is_open()) { std::cerr << "Error opening output file!" << std::endl; return; }

320
TrackRecon.C
View File

@ -3,7 +3,7 @@
#define RAW_HISTOS #define RAW_HISTOS
// #define VTX_GATES // #define VTX_GATES
// #define AL_BEAM // #define AL_BEAM
#define F_BEAM // #define F_BEAM
// #define nA_Analysis // #define nA_Analysis
Int_t colors[40] = { Int_t colors[40] = {
@ -40,9 +40,9 @@ Int_t colors[40] = {
bool process_alpha_proton_scattering = false; bool process_alpha_proton_scattering = false;
bool doMiscHistograms = true; bool doMiscHistograms = true;
bool doPCSX3ClusterAnalysis = false; bool doPCSX3ClusterAnalysis = true;
bool doPCQQQClusterAnalysis = false; bool doPCQQQClusterAnalysis = true;
bool doOldAnalysis = true; bool doOldAnalysis = false;
bool do27AlapAnalysis = false; bool do27AlapAnalysis = false;
double source_vertex = 53; // 53 double source_vertex = 53; // 53
const double qqq_z = 100.0; const double qqq_z = 100.0;
@ -140,8 +140,8 @@ int anodeIndex = -1, cathodeIndex = -1;
void protonAlphaHistograms(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events); void protonAlphaHistograms(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events);
void miscHistograms_oneWire(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<std::vector<std::tuple<int, double, double>>> aClusters); void miscHistograms_oneWire(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<std::vector<std::tuple<int, double, double>>> aClusters);
void protonMiscHistograms(HistPlotter* plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events); void protonMiscHistograms(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events);
void protonMiscHistograms_sx3(HistPlotter* plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events); void protonMiscHistograms_sx3(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events);
void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events); void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events);
void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events); void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events);
@ -415,11 +415,15 @@ Bool_t TrackRecon::Process(Long64_t entry)
{ {
// 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", 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_id_" + std::to_string(id), 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_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)
{ {
@ -447,6 +451,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
Event sx3ev(TVector3(rho_at_strip * TMath::Cos(phi_n), rho_at_strip * TMath::Sin(phi_n), z), backE * 0.001, -1, det.ts, -1, det.stripB + 4 * id, det.stripF + 4 * id); Event sx3ev(TVector3(rho_at_strip * TMath::Cos(phi_n), rho_at_strip * TMath::Sin(phi_n), z), backE * 0.001, -1, det.ts, -1, det.stripB + 4 * id, det.stripF + 4 * id);
SX3_Events.push_back(sx3ev); SX3_Events.push_back(sx3ev);
plotter->Fill2D("sx3backs_gm", 100, 0, 100, 800, 0, 8192, det.stripB + 4 * id, backE); plotter->Fill2D("sx3backs_gm", 100, 0, 100, 800, 0, 8192, det.stripB + 4 * id, backE);
plotter->Fill1D("sx3backs_calib", 800, 0, 8192, backE);
// 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, 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, 88.0 * TMath::Cos(phi_n), 88.0 * TMath::Sin(phi_n), "hCalSX3");
@ -564,7 +569,6 @@ Bool_t TrackRecon::Process(Long64_t entry)
plotter->Fill2D("RingE_vs_Index", 16 * 4, 0, 16 * 4, 1000, 0, 16000, chRing + qqq.id[i] * 16, eRing, "hRawQQQ"); plotter->Fill2D("RingE_vs_Index", 16 * 4, 0, 16 * 4, 1000, 0, 16000, chRing + qqq.id[i] * 16, eRing, "hRawQQQ");
plotter->Fill2D("WedgeE_vs_Index", 16 * 4, 0, 16 * 4, 1000, 0, 16000, chWedge + qqq.id[i] * 16, eWedge, "hRawQQQ"); plotter->Fill2D("WedgeE_vs_Index", 16 * 4, 0, 16 * 4, 1000, 0, 16000, chWedge + qqq.id[i] * 16, eWedge, "hRawQQQ");
#endif #endif
plotter->Fill2D("WedgeE_Vs_RingECal", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ");
if (qqqCalibValid[qqq.id[i]][chWedge][chRing]) if (qqqCalibValid[qqq.id[i]][chWedge][chRing])
{ {
@ -586,6 +590,9 @@ Bool_t TrackRecon::Process(Long64_t entry)
QQQ_Events_Raw.push_back(qqqeventr); QQQ_Events_Raw.push_back(qqqeventr);
plotter->Fill2D("WedgeE_Vs_RingECal_selected", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ"); plotter->Fill2D("WedgeE_Vs_RingECal_selected", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ");
plotter->Fill1D("QQQECal", 2048, 0, 10, eRingMeV);
plotter->Fill1D("QQQECal", 2048, 0, 10, eWedgeMeV);
const int channelsPerDetector = MAX_RING + MAX_WEDGE; const int channelsPerDetector = MAX_RING + MAX_WEDGE;
int globalRingChannel = chRing + (qqq.id[i] * channelsPerDetector); int globalRingChannel = chRing + (qqq.id[i] * channelsPerDetector);
int globalWedgeChannel = chWedge + (qqq.id[i] * channelsPerDetector) + MAX_RING; int globalWedgeChannel = chWedge + (qqq.id[i] * channelsPerDetector) + MAX_RING;
@ -801,7 +808,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
bool found_rf = false; bool found_rf = false;
bool found_mcp = false; bool found_mcp = false;
bool found_needle = false; bool found_needle = false;
bool qqq_inner_ring = (qqqevent.ch1 % 16) < 6; bool qqq_inner_ring = (qqqevent.ch1 % 16) < 8;
for (int j = 0; j < misc.multi; j++) for (int j = 0; j < misc.multi; j++)
{ {
plotter->Fill1D("channels_misc_qqq", 20, 0, 20, misc.ch[j], "misc"); plotter->Fill1D("channels_misc_qqq", 20, 0, 20, misc.ch[j], "misc");
@ -816,25 +823,25 @@ Bool_t TrackRecon::Process(Long64_t entry)
{ // 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]);
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), 800, -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]);
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), 800, -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), 500, -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), 640, -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), 640, -2000, 2000, 640, -2000, 2000, dt_qqq_rf, dt_rf_mcp, "misc");
plotter->Fill2D("dt_(qqq,mcp)_vs_(qqq,rf)_innerring"+std::to_string(qqq_inner_ring), 640, -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), 640, -1400, 2000, 640, -2000, 2000, dt_qqq_mcp, dt_qqq_rf, "misc");
plotter->Fill2D("dt_(qqq,mcp)_vs_(rf,mcp)_innerring"+std::to_string(qqq_inner_ring), 640, -1400, -600, 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), 640, -1400, -600, 640, -2000, 2000, dt_qqq_mcp, dt_rf_mcp, "misc");
} }
ctr += 1; ctr += 1;
} }
@ -1078,16 +1085,21 @@ Bool_t TrackRecon::Process(Long64_t entry)
double theta_recon = (sx3event.pos - TVector3(0, 0, vertex_recon)).Theta(); double theta_recon = (sx3event.pos - TVector3(0, 0, vertex_recon)).Theta();
double sinTheta = TMath::Sin(theta_recon); double sinTheta = TMath::Sin(theta_recon);
double Ex_from_proton = apkin_p.getExc(sx3Efix, theta_recon * 180. / M_PI);
double Ex_from_alpha = apkin_a.getExc(sx3Efixalpha, theta_recon * 180. / M_PI);
// Fill standard un-gated plots // Fill standard un-gated plots
plotter->Fill1D(nA_label + "_vertex_recon_SX3", 400, -200, 200, vertex_recon, nA_label);
plotter->Fill1D(nA_label + "_vertex_recon", 400, -200, 200, vertex_recon, nA_label);
plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_SX3", 800, 0, 30, 800, 0, 30000, sx3Efix, apSumE * sinTheta, nA_label); plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_SX3", 800, 0, 30, 800, 0, 30000, sx3Efix, apSumE * sinTheta, nA_label);
plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_SX3_alpha", 800, 0, 30, 800, 0, 30000, sx3Efixalpha, apSumE * sinTheta, nA_label); plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_SX3_alpha", 800, 0, 30, 800, 0, 30000, sx3Efixalpha, apSumE * sinTheta, nA_label);
plotter->Fill2D(nA_label + "_sx3_E_vs_theta_raw_SX3", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, sx3event.Energy1, nA_label); plotter->Fill2D(nA_label + "_sx3_E_vs_theta_raw_SX3", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, sx3event.Energy1, nA_label);
plotter->Fill2D(nA_label + "_sx3_E_vs_theta_corr_SX3", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, sx3Efix, nA_label); plotter->Fill2D(nA_label + "_sx3_E_vs_theta_corr_SX3", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, sx3Efix, nA_label);
#if defined(F_BEAM) || defined(AL_BEAM)
double Ex_from_proton = apkin_p.getExc(sx3Efix, theta_recon * 180. / M_PI);
double Ex_from_alpha = apkin_a.getExc(sx3Efixalpha, theta_recon * 180. / M_PI);
plotter->Fill1D(nA_label + "_Ex_from_protons_SX3", 1200, -30, 30, Ex_from_proton, nA_label); plotter->Fill1D(nA_label + "_Ex_from_protons_SX3", 1200, -30, 30, Ex_from_proton, nA_label);
plotter->Fill1D(nA_label + "_Ex_from_alphas_SX3", 1200, -30, 30, Ex_from_alpha, nA_label); plotter->Fill1D(nA_label + "_Ex_from_alphas_SX3", 1200, -30, 30, Ex_from_alpha, nA_label);
#endif
// Fill Gated Plots // Fill Gated Plots
// if (vtx_gate != "") // if (vtx_gate != "")
@ -1172,16 +1184,19 @@ Bool_t TrackRecon::Process(Long64_t entry)
double theta_recon = (qqqevent.pos - TVector3(0, 0, vertex_recon)).Theta(); double theta_recon = (qqqevent.pos - TVector3(0, 0, vertex_recon)).Theta();
double sinTheta = TMath::Sin(theta_recon); double sinTheta = TMath::Sin(theta_recon);
double Ex_from_proton = apkin_p.getExc(qqqEfix, theta_recon * 180. / M_PI);
double Ex_from_alpha = apkin_a.getExc(qqqEfixalpha, theta_recon * 180. / M_PI);
// Fill standard un-gated plots // Fill standard un-gated plots
plotter->Fill1D(nA_label + "_vertex_recon_QQQ", 400, -200, 200, vertex_recon, nA_label);
plotter->Fill1D(nA_label + "_vertex_recon", 400, -200, 200, vertex_recon, nA_label);
plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_QQQ", 800, 0, 30, 800, 0, 30000, qqqEfix, apSumE * sinTheta, nA_label); plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_QQQ", 800, 0, 30, 800, 0, 30000, qqqEfix, apSumE * sinTheta, nA_label);
plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_QQQ_alpha", 800, 0, 30, 800, 0, 30000, qqqEfixalpha, apSumE * sinTheta, nA_label); plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_QQQ_alpha", 800, 0, 30, 800, 0, 30000, qqqEfixalpha, apSumE * sinTheta, nA_label);
plotter->Fill2D(nA_label + "_qqq_E_vs_theta_raw_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqevent.Energy1, nA_label); plotter->Fill2D(nA_label + "_qqq_E_vs_theta_raw_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqevent.Energy1, nA_label);
plotter->Fill2D(nA_label + "_qqq_E_vs_theta_corr_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqEfix, nA_label); plotter->Fill2D(nA_label + "_qqq_E_vs_theta_corr_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqEfix, nA_label);
#if defined(F_BEAM) || defined(AL_BEAM)
double Ex_from_alpha = apkin_a.getExc(qqqEfixalpha, theta_recon * 180. / M_PI);
plotter->Fill1D(nA_label + "_Ex_from_alphas_QQQ", 1200, -30, 30, Ex_from_alpha, nA_label); plotter->Fill1D(nA_label + "_Ex_from_alphas_QQQ", 1200, -30, 30, Ex_from_alpha, nA_label);
plotter->Fill1D(nA_label + "_Ex_from_protons_QQQ", 1200, -30, 30, Ex_from_proton, nA_label); plotter->Fill1D(nA_label + "_Ex_from_protons_QQQ", 1200, -30, 30, Ex_from_proton, nA_label);
#endif
// Fill Gated Plots // Fill Gated Plots
// if (vtx_gate != "") // if (vtx_gate != "")
@ -1400,6 +1415,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
double t_minimum = -1.0 * (x1.X() * v.X() + x1.Y() * v.Y()) / (v.X() * v.X() + v.Y() * v.Y()); double t_minimum = -1.0 * (x1.X() * v.X() + x1.Y() * v.Y()) / (v.X() * v.X() + v.Y() * v.Y());
TVector3 r_rhoMin_fix = x1 + t_minimum * v; TVector3 r_rhoMin_fix = x1 + t_minimum * v;
plotter->Fill1D("VertexRecon_pczfix_sx3", 800, -300, 300, r_rhoMin_fix.Z()); plotter->Fill1D("VertexRecon_pczfix_sx3", 800, -300, 300, r_rhoMin_fix.Z());
plotter->Fill1D("VertexRecon_pczfix", 800, -300, 300, r_rhoMin_fix.Z());
plotter->Fill1D("pczfix_A1C2_1d_sx3", 600, -200, 200, pcz_fix); plotter->Fill1D("pczfix_A1C2_1d_sx3", 600, -200, 200, pcz_fix);
plotter->Fill2D("pczfix_vs_sx3pczguess_A1C2", 600, -200, 200, 600, -200, 200, pczguess, pcz_fix); plotter->Fill2D("pczfix_vs_sx3pczguess_A1C2", 600, -200, 200, 600, -200, 200, pczguess, pcz_fix);
plotter->Fill2D("pcz_vs_sx3pczguess_A1C2_strip" + std::to_string(sx3event.ch2), 300, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z()); plotter->Fill2D("pcz_vs_sx3pczguess_A1C2_strip" + std::to_string(sx3event.ch2), 300, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z());
@ -1549,6 +1565,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
plotter->Fill2D("dE3_E_AnodeQQQR_A1C2_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy1 * sinTheta_customV); plotter->Fill2D("dE3_E_AnodeQQQR_A1C2_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy1 * sinTheta_customV);
plotter->Fill1D("VertexRecon_pczfix_qqq", 800, -400, 400, r_rhoMin_fix.Z()); plotter->Fill1D("VertexRecon_pczfix_qqq", 800, -400, 400, r_rhoMin_fix.Z());
plotter->Fill1D("VertexRecon_pczfix", 800, -400, 400, r_rhoMin_fix.Z());
plotter->Fill1D("VertexRecon_pczfix_qqq_PC" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath), 800, -400, 400, r_rhoMin_fix.Z()); plotter->Fill1D("VertexRecon_pczfix_qqq_PC" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath), 800, -400, 400, r_rhoMin_fix.Z());
if (TMath::Abs(r_rhoMin_fix.Z()) < 200.0) if (TMath::Abs(r_rhoMin_fix.Z()) < 200.0)
@ -2098,174 +2115,191 @@ void miscHistograms_oneWire(HistPlotter *plotter, std::vector<Event> QQQ_Events,
} // end QQQEvents loop } // end QQQEvents loop
} }
void protonMiscHistograms(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events)
void protonMiscHistograms(HistPlotter* plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events) { {
//consider the 'proton-like' QQQ branch seen in a,p data // consider the 'proton-like' QQQ branch seen in a,p data
TRandom3 rand; TRandom3 rand;
rand.SetSeed();//random seed set rand.SetSeed(); // random seed set
Kinematics apkin_a(1.008664916,4.002603254,4.002603254,1.008664916,7.1); //m3 is alpha, 6.79 MeV is 7.0 MeV proton energy after kapton+100mm 4He gas (molar mass 5.2, 250 torr) Kinematics apkin_a(1.008664916, 4.002603254, 4.002603254, 1.008664916, 7.1); // m3 is alpha, 6.79 MeV is 7.0 MeV proton energy after kapton+100mm 4He gas (molar mass 5.2, 250 torr)
for(auto qqqevent: QQQ_Events) { for (auto qqqevent : QQQ_Events)
if(qqqevent.Energy1 < 0.6) continue; //coarse gating {
//if(qqqevent.Energy1 > 5.0) continue; //coarse gating if (qqqevent.Energy1 < 0.6)
for(auto pcevent: PC_Events) { continue; // coarse gating
if(!(pcevent.multi1==1 && pcevent.multi2<=2)) continue; // if(qqqevent.Energy1 > 5.0) continue; //coarse gating
//if(pcevent.Energy1 > 11000) continue; //coarse gating for (auto pcevent : PC_Events)
bool phicut = qqqevent.pos.Phi() <= pcevent.pos.Phi()+TMath::Pi()/4. && qqqevent.pos.Phi() >= pcevent.pos.Phi()-TMath::Pi()/4.; {
if(!phicut) continue; if (!(pcevent.multi1 == 1 && pcevent.multi2 <= 2))
//if(pcevent.Time1-qqqevent.Time1<-150 || pcevent.Time1-qqqevent.Time1 >850) continue; continue;
// if(pcevent.Energy1 > 11000) continue; //coarse gating
bool phicut = qqqevent.pos.Phi() <= pcevent.pos.Phi() + TMath::Pi() / 4. && qqqevent.pos.Phi() >= pcevent.pos.Phi() - TMath::Pi() / 4.;
if (!phicut)
continue;
// if(pcevent.Time1-qqqevent.Time1<-150 || pcevent.Time1-qqqevent.Time1 >850) continue;
double pcz_fix, pcz_dith=pcevent.pos.Z(); double pcz_fix, pcz_dith = pcevent.pos.Z();
if(pcevent.multi2==2) if (pcevent.multi2 == 2)
pcz_fix = pcfix_func.Eval(pcevent.pos.Z()); pcz_fix = pcfix_func.Eval(pcevent.pos.Z());
else { else
pcz_fix = rand.Gaus(pcevent.pos.Z(),8.0);//dither for a1c1 events {
pcz_fix = rand.Gaus(pcevent.pos.Z(), 8.0); // dither for a1c1 events
pcz_dith = pcz_fix; pcz_dith = pcz_fix;
} }
TVector3 x2f(pcevent.pos.X(),pcevent.pos.Y(),pcz_fix); TVector3 x2f(pcevent.pos.X(), pcevent.pos.Y(), pcz_fix);
TVector3 x1(qqqevent.pos); TVector3 x1(qqqevent.pos);
TVector3 v = x2f-x1; TVector3 v = x2f - x1;
double t_minimum = -1.0*(x1.X()*v.X()+x1.Y()*v.Y())/(v.X()*v.X()+v.Y()*v.Y()); double t_minimum = -1.0 * (x1.X() * v.X() + x1.Y() * v.Y()) / (v.X() * v.X() + v.Y() * v.Y());
TVector3 r_rhoMin_fix = x1 + t_minimum*v; TVector3 r_rhoMin_fix = x1 + t_minimum * v;
double vertex_z = r_rhoMin_fix.Z(); double vertex_z = r_rhoMin_fix.Z();
//double theta_q = (qqqevent.pos - TVector3(0,0,vertex_z)).Theta(); // double theta_q = (qqqevent.pos - TVector3(0,0,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);
//if(r_rhoMin_fix.Perp()>6) continue; // if(r_rhoMin_fix.Perp()>6) continue;
bool cathode_alpha_select = (pcevent.Energy2 > 1400); bool cathode_alpha_select = (pcevent.Energy2 > 1400);
if(vertex_z < -173.6 || vertex_z > 100) continue; if (vertex_z < -173.6 || vertex_z > 100)
continue;
// What's below: radial cut, time coincident, phi-correlated events with possible energy selection applied to both E_si and dE_Anodes
//What's below: radial cut, time coincident, phi-correlated events with possible energy selection applied to both E_si and dE_Anodes auto plot_with_tag = [&](std::string tag = "")
auto plot_with_tag = [&](std::string tag="") { {
std::string pmlabel = "proton+misc"+tag; std::string pmlabel = "proton+misc" + tag;
plotter->Fill2D("pmisc_dE_E_AnodeQQQ"+tag,400,0,10,800,0,40000,qqqevent.Energy1,pcevent.Energy1,pmlabel); plotter->Fill2D("pmisc_dE_E_AnodeQQQ" + tag, 400, 0, 10, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, pmlabel);
plotter->Fill2D("pmisc_dE_E_CathodeQQQ"+tag,400,0,10,800,0,10000,qqqevent.Energy1,pcevent.Energy2,pmlabel); plotter->Fill2D("pmisc_dE_E_CathodeQQQ" + tag, 400, 0, 10, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2, pmlabel);
plotter->Fill2D("pmisc_dE3_E_AnodeQQQ"+tag,400,0,10,400,0,40000,qqqevent.Energy1,pcevent.Energy1*sinTheta_customV*3.,pmlabel); plotter->Fill2D("pmisc_dE3_E_AnodeQQQ" + tag, 400, 0, 10, 400, 0, 40000, qqqevent.Energy1, pcevent.Energy1 * sinTheta_customV * 3., pmlabel);
plotter->Fill2D("pmisc_dE3_E_CathodeQQQ"+tag,400,0,10,400,0,10000,qqqevent.Energy1,pcevent.Energy2*sinTheta_customV,pmlabel); plotter->Fill2D("pmisc_dE3_E_CathodeQQQ" + tag, 400, 0, 10, 400, 0, 10000, qqqevent.Energy1, pcevent.Energy2 * sinTheta_customV, pmlabel);
plotter->Fill2D("pmisc_dPhi_QQQ_PC"+tag,180,-360,360,180,-360,360,pcevent.pos.Phi()*180/M_PI,qqqevent.pos.Phi()*180/M_PI,pmlabel); plotter->Fill2D("pmisc_dPhi_QQQ_PC" + tag, 180, -360, 360, 180, -360, 360, pcevent.pos.Phi() * 180 / M_PI, qqqevent.pos.Phi() * 180 / M_PI, pmlabel);
plotter->Fill1D("pmisc_dt_Anode_QQQ_PC"+std::to_string(phicut)+tag,600,-2000,2000,pcevent.Time1-qqqevent.Time1,pmlabel); plotter->Fill1D("pmisc_dt_Anode_QQQ_PC" + std::to_string(phicut) + tag, 600, -2000, 2000, pcevent.Time1 - qqqevent.Time1, pmlabel);
plotter->Fill1D("pmisc_dt_Cathode_QQQ"+tag,600,-2000,2000,pcevent.Time2-qqqevent.Time1,pmlabel); plotter->Fill1D("pmisc_dt_Cathode_QQQ" + tag, 600, -2000, 2000, pcevent.Time2 - qqqevent.Time1, pmlabel);
plotter->Fill2D("pmisc_dt_Anode_E_QQQ_PC"+std::to_string(phicut)+tag,600,-2000,2000,400,0,10,pcevent.Time1-qqqevent.Time1,qqqevent.Energy1,pmlabel); plotter->Fill2D("pmisc_dt_Anode_E_QQQ_PC" + std::to_string(phicut) + tag, 600, -2000, 2000, 400, 0, 10, pcevent.Time1 - qqqevent.Time1, qqqevent.Energy1, pmlabel);
plotter->Fill2D("pmisc_dt_AnodeQQQ_vsPCPhi"+tag,600,-2000,2000,180,-360,360,pcevent.Time1-qqqevent.Time1,pcevent.pos.Phi()*180./M_PI,pmlabel); plotter->Fill2D("pmisc_dt_AnodeQQQ_vsPCPhi" + tag, 600, -2000, 2000, 180, -360, 360, pcevent.Time1 - qqqevent.Time1, pcevent.pos.Phi() * 180. / M_PI, pmlabel);
plotter->Fill2D("pmisc_dt_Cathode_E_QQQ"+tag,600,-2000,2000,400,0,10,pcevent.Time2-qqqevent.Time1,qqqevent.Energy1,pmlabel); plotter->Fill2D("pmisc_dt_Cathode_E_QQQ" + tag, 600, -2000, 2000, 400, 0, 10, pcevent.Time2 - qqqevent.Time1, qqqevent.Energy1, pmlabel);
plotter->Fill2D("pmisc_dt_CathodeQQQ_vsPCPhi"+tag,600,-2000,2000,180,-360,360,pcevent.Time2-qqqevent.Time1,pcevent.pos.Phi()*180./M_PI,pmlabel); plotter->Fill2D("pmisc_dt_CathodeQQQ_vsPCPhi" + tag, 600, -2000, 2000, 180, -360, 360, pcevent.Time2 - qqqevent.Time1, pcevent.pos.Phi() * 180. / M_PI, pmlabel);
plotter->Fill1D("pmisc_pczfix"+tag,600,-300,300,pcz_fix,pmlabel); plotter->Fill1D("pmisc_pczfix" + tag, 600, -300, 300, pcz_fix, pmlabel);
if(pcevent.multi2==2) { if (pcevent.multi2 == 2)
plotter->Fill1D("pmisc_pcz"+tag,600,-300,300,pcevent.pos.Z(),pmlabel); {
plotter->Fill1D("pmisc_pcz2"+tag,600,-300,300,pcevent.pos.Z(),pmlabel); plotter->Fill1D("pmisc_pcz" + tag, 600, -300, 300, pcevent.pos.Z(), pmlabel);
plotter->Fill1D("pmisc_pcz2" + tag, 600, -300, 300, pcevent.pos.Z(), pmlabel);
} }
if(pcevent.multi2==1) { if (pcevent.multi2 == 1)
plotter->Fill1D("pmisc_pcz"+tag,600,-300,300,pcz_fix,pmlabel); {
plotter->Fill1D("pmisc_pcz1"+tag,600,-300,300,pcevent.pos.Z(),pmlabel); plotter->Fill1D("pmisc_pcz" + tag, 600, -300, 300, pcz_fix, pmlabel);
plotter->Fill1D("pmisc_pcz_dith"+tag,600,-300,300,pcz_dith,pmlabel); plotter->Fill1D("pmisc_pcz1" + tag, 600, -300, 300, pcevent.pos.Z(), pmlabel);
plotter->Fill1D("pmisc_pcz_dith" + tag, 600, -300, 300, pcz_dith, pmlabel);
} }
//double path_length_q = (qqqevent.pos-TVector3(0,0,vertex_z)).Mag()*0.1; // double path_length_q = (qqqevent.pos-TVector3(0,0,vertex_z)).Mag()*0.1;
//double path_length_s = (sx3event.pos-TVector3(0,0,vertex_z)).Mag()*0.1; // double path_length_s = (sx3event.pos-TVector3(0,0,vertex_z)).Mag()*0.1;
double path_length_q = (qqqevent.pos-r_rhoMin_fix).Mag()*0.1; double path_length_q = (qqqevent.pos - r_rhoMin_fix).Mag() * 0.1;
double qqqEfix; double qqqEfix;
if(tag == "_cathode_alphas") {//satisfied when find succeeds if (tag == "_cathode_alphas")
qqqEfix = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy1)-path_length_q); { // satisfied when find succeeds
plotter->Fill1D("pmisc_Ex_from_alpha",200,-10,10,apkin_a.getExc(qqqEfix,theta_q*180/M_PI),pmlabel); qqqEfix = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy1) - path_length_q);
plotter->Fill1D("pmisc_Ex_from_alpha", 200, -10, 10, apkin_a.getExc(qqqEfix, theta_q * 180 / M_PI), pmlabel);
} }
else else
qqqEfix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(qqqevent.Energy1)-path_length_q); qqqEfix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(qqqevent.Energy1) - path_length_q);
//plotter->Fill2D("qqqEf_sx3E_matrix_all"+tag,400,0,10,400,0,10,qqqEfix,sx3event.Energy1,pmlabel); // plotter->Fill2D("qqqEf_sx3E_matrix_all"+tag,400,0,10,400,0,10,qqqEfix,sx3event.Energy1,pmlabel);
plotter->Fill2D("pmisc_dE3_Ef_AnodeQQQ"+tag,400,0,10,400,0,40000,qqqEfix,pcevent.Energy1*sinTheta_customV*3,pmlabel); plotter->Fill2D("pmisc_dE3_Ef_AnodeQQQ" + tag, 400, 0, 10, 400, 0, 40000, qqqEfix, pcevent.Energy1 * sinTheta_customV * 3, pmlabel);
plotter->Fill2D("pmisc_dE3_Ef_CathodeQQQ"+tag,400,0,10,400,0,10000,qqqEfix,pcevent.Energy2*sinTheta_customV,pmlabel); plotter->Fill2D("pmisc_dE3_Ef_CathodeQQQ" + tag, 400, 0, 10, 400, 0, 10000, qqqEfix, pcevent.Energy2 * sinTheta_customV, pmlabel);
plotter->Fill1D("pmisc_VertexReconZ"+tag,800,-400,400,vertex_z,pmlabel); plotter->Fill1D("pmisc_VertexReconZ" + tag, 800, -400, 400, vertex_z, pmlabel);
plotter->Fill2D("pmisc_VertexReconXY"+tag,200,-100,100,200,-100,100,r_rhoMin_fix.X(),r_rhoMin_fix.Y(),pmlabel); plotter->Fill2D("pmisc_VertexReconXY" + tag, 200, -100, 100, 200, -100, 100, r_rhoMin_fix.X(), r_rhoMin_fix.Y(), pmlabel);
plotter->Fill2D("pmisc_VertexReconZ_vs_Ef"+tag,800,-400,400,800,0,20,vertex_z,qqqEfix,pmlabel); plotter->Fill2D("pmisc_VertexReconZ_vs_Ef" + tag, 800, -400, 400, 800, 0, 20, vertex_z, qqqEfix, pmlabel);
plotter->Fill2D("pmisc_VertexReconZ_vs_Ef"+tag+"_a"+std::to_string(pcevent.multi1),800,-400,400,800,0,20,vertex_z,qqqEfix,pmlabel); plotter->Fill2D("pmisc_VertexReconZ_vs_Ef" + tag + "_a" + std::to_string(pcevent.multi1), 800, -400, 400, 800, 0, 20, vertex_z, qqqEfix, pmlabel);
plotter->Fill2D("pmisc_Ef_vs_theta_qqq"+tag,100,0,180,800,0,20,theta_q*180/M_PI,qqqEfix,pmlabel); plotter->Fill2D("pmisc_Ef_vs_theta_qqq" + tag, 100, 0, 180, 800, 0, 20, theta_q * 180 / M_PI, qqqEfix, pmlabel);
if(pcevent.multi2==1) { if (pcevent.multi2 == 1)
plotter->Fill2D("pmisc_Ef_vs_theta_qqq_a1c1"+tag,100,0,180,800,0,20,theta_q*180/M_PI,qqqEfix,pmlabel); {
plotter->Fill2D("pmisc_VertexReconZ_vs_Ef_a1c1"+tag,800,-400,400,800,0,20,vertex_z,qqqEfix,pmlabel); plotter->Fill2D("pmisc_Ef_vs_theta_qqq_a1c1" + tag, 100, 0, 180, 800, 0, 20, theta_q * 180 / M_PI, qqqEfix, pmlabel);
plotter->Fill2D("pmisc_VertexReconZ_vs_Ef_a1c1" + tag, 800, -400, 400, 800, 0, 20, vertex_z, qqqEfix, pmlabel);
} }
}; };
if(cathode_alpha_select) if (cathode_alpha_select)
plot_with_tag("_cathode_alphas"); plot_with_tag("_cathode_alphas");
//else // else
// plot_with_tag("_cathode_protons"); // plot_with_tag("_cathode_protons");
//plot_with_tag(); // plot_with_tag();
//plotter->Fill1D("pmisc_Ex_from_protons",200,-10,10,apkin_p.getExc(qqqEfix,theta_s*180/M_PI),pmlabel); // plotter->Fill1D("pmisc_Ex_from_protons",200,-10,10,apkin_p.getExc(qqqEfix,theta_s*180/M_PI),pmlabel);
}//end PCEvents loop } // end PCEvents loop
}//end QQQEvents loop } // end QQQEvents loop
} }
void protonMiscHistograms_sx3(HistPlotter* plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events) { void protonMiscHistograms_sx3(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events)
//consider the 'proton-like' QQQ branch seen in a,p data {
for(auto sx3event: SX3_Events) { // consider the 'proton-like' QQQ branch seen in a,p data
if(sx3event.Energy1 < 1.2) continue; //coarse gating for (auto sx3event : SX3_Events)
//if(sx3event.Energy1 > 5.0) continue; //coarse gating {
for(auto pcevent: PC_Events) { if (sx3event.Energy1 < 1.2)
if(!(pcevent.multi1==1 && pcevent.multi2==2)) continue; continue; // coarse gating
//if(pcevent.Energy1 > 11000) continue; //coarse gating // if(sx3event.Energy1 > 5.0) continue; //coarse gating
bool phicut = sx3event.pos.Phi() <= pcevent.pos.Phi()+TMath::Pi()/3. && sx3event.pos.Phi() >= pcevent.pos.Phi()-TMath::Pi()/3.; for (auto pcevent : PC_Events)
if(!phicut) continue; {
//if(pcevent.Time1-sx3event.Time1<-150 || pcevent.Time1-sx3event.Time1 >850) continue; if (!(pcevent.multi1 == 1 && pcevent.multi2 == 2))
continue;
// if(pcevent.Energy1 > 11000) continue; //coarse gating
bool phicut = sx3event.pos.Phi() <= pcevent.pos.Phi() + TMath::Pi() / 3. && sx3event.pos.Phi() >= pcevent.pos.Phi() - TMath::Pi() / 3.;
if (!phicut)
continue;
// if(pcevent.Time1-sx3event.Time1<-150 || pcevent.Time1-sx3event.Time1 >850) continue;
double pcz_fix = pcfix_func.Eval(pcevent.pos.Z()); double pcz_fix = pcfix_func.Eval(pcevent.pos.Z());
TVector3 x2f(pcevent.pos.X(),pcevent.pos.Y(),pcz_fix); TVector3 x2f(pcevent.pos.X(), pcevent.pos.Y(), pcz_fix);
TVector3 x1(sx3event.pos); TVector3 x1(sx3event.pos);
TVector3 v = x2f-x1; TVector3 v = x2f - x1;
double t_minimum = -1.0*(x1.X()*v.X()+x1.Y()*v.Y())/(v.X()*v.X()+v.Y()*v.Y()); double t_minimum = -1.0 * (x1.X() * v.X() + x1.Y() * v.Y()) / (v.X() * v.X() + v.Y() * v.Y());
TVector3 r_rhoMin_fix = x1 + t_minimum*v; TVector3 r_rhoMin_fix = x1 + t_minimum * v;
double vertex_z = r_rhoMin_fix.Z(); double vertex_z = r_rhoMin_fix.Z();
//double theta_q = (sx3event.pos - TVector3(0,0,vertex_z)).Theta(); // double theta_q = (sx3event.pos - TVector3(0,0,vertex_z)).Theta();
if(r_rhoMin_fix.Perp()>10.0) continue; if (r_rhoMin_fix.Perp() > 10.0)
continue;
double theta_s = (sx3event.pos - r_rhoMin_fix).Theta(); double theta_s = (sx3event.pos - r_rhoMin_fix).Theta();
double sinTheta_customV = TMath::Sin(theta_s); double sinTheta_customV = TMath::Sin(theta_s);
bool cathode_alpha_select = (pcevent.Energy2 > 1400); bool cathode_alpha_select = (pcevent.Energy2 > 1400);
// What's below: radial cut, time coincident, phi-correlated events with possible energy selection applied to both E_si and dE_Anodes
auto plot_with_tag = [&](std::string tag = "")
{
std::string pmlabel = "proton+miscsx3" + tag;
plotter->Fill2D("pmiscs_dE_E_Anodesx3" + tag, 400, 0, 10, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1, pmlabel);
plotter->Fill2D("pmiscs_dE_E_Cathodesx3" + tag, 400, 0, 10, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2, pmlabel);
plotter->Fill2D("pmiscs_dE3_E_Anodesx3" + tag, 400, 0, 10, 400, 0, 40000, sx3event.Energy1, pcevent.Energy1 * sinTheta_customV * 3., pmlabel);
plotter->Fill2D("pmiscs_dE3_E_Cathodesx3" + tag, 400, 0, 10, 400, 0, 10000, sx3event.Energy1, pcevent.Energy2 * sinTheta_customV, pmlabel);
plotter->Fill2D("pmiscs_dPhi_sx3_PC" + tag, 180, -360, 360, 180, -360, 360, pcevent.pos.Phi() * 180 / M_PI, sx3event.pos.Phi() * 180 / M_PI, pmlabel);
plotter->Fill1D("pmiscs_dt_Anode_sx3_PC" + std::to_string(phicut) + tag, 600, -2000, 2000, pcevent.Time1 - sx3event.Time1, pmlabel);
plotter->Fill1D("pmiscs_dt_Cathode_sx3" + tag, 600, -2000, 2000, pcevent.Time2 - sx3event.Time1, pmlabel);
plotter->Fill2D("pmiscs_dt_Anode_E_sx3_PC" + std::to_string(phicut) + tag, 600, -2000, 2000, 400, 0, 10, pcevent.Time1 - sx3event.Time1, sx3event.Energy1, pmlabel);
plotter->Fill2D("pmiscs_dt_Cathode_E_sx3" + tag, 600, -2000, 2000, 400, 0, 10, pcevent.Time2 - sx3event.Time1, sx3event.Energy1, pmlabel);
plotter->Fill2D("pmiscs_dt_Cathodesx3_vsPCPhi" + tag, 600, -2000, 2000, 180, -360, 360, pcevent.Time2 - sx3event.Time1, pcevent.pos.Phi() * 180. / M_PI, pmlabel);
plotter->Fill1D("pmiscs_pczfix" + tag, 600, -300, 300, pcz_fix, pmlabel);
plotter->Fill1D("pmiscs_pcz" + tag, 600, -300, 300, pcevent.pos.Z(), pmlabel);
// double path_length_q = (sx3event.pos-TVector3(0,0,vertex_z)).Mag()*0.1;
// double path_length_s = (sx3event.pos-TVector3(0,0,vertex_z)).Mag()*0.1;
double path_length_s = (sx3event.pos - r_rhoMin_fix).Mag() * 0.1;
double sx3Efix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(sx3event.Energy1) - path_length_s);
//What's below: radial cut, time coincident, phi-correlated events with possible energy selection applied to both E_si and dE_Anodes // plotter->Fill2D("sx3Ef_sx3E_matrix_all"+tag,400,0,10,400,0,10,sx3Efix,sx3event.Energy1,pmlabel);
auto plot_with_tag = [&](std::string tag="") { plotter->Fill2D("pmiscs_dE3_Ef_Anodesx3" + tag, 400, 0, 10, 400, 0, 40000, sx3Efix, pcevent.Energy1 * sinTheta_customV * 3, pmlabel);
std::string pmlabel = "proton+miscsx3"+tag; plotter->Fill2D("pmiscs_dE3_Ef_Cathodesx3" + tag, 400, 0, 10, 400, 0, 10000, sx3Efix, pcevent.Energy2 * sinTheta_customV, pmlabel);
plotter->Fill2D("pmiscs_dE_E_Anodesx3"+tag,400,0,10,800,0,40000,sx3event.Energy1,pcevent.Energy1,pmlabel);
plotter->Fill2D("pmiscs_dE_E_Cathodesx3"+tag,400,0,10,800,0,10000,sx3event.Energy1,pcevent.Energy2,pmlabel);
plotter->Fill2D("pmiscs_dE3_E_Anodesx3"+tag,400,0,10,400,0,40000,sx3event.Energy1,pcevent.Energy1*sinTheta_customV*3.,pmlabel);
plotter->Fill2D("pmiscs_dE3_E_Cathodesx3"+tag,400,0,10,400,0,10000,sx3event.Energy1,pcevent.Energy2*sinTheta_customV,pmlabel);
plotter->Fill2D("pmiscs_dPhi_sx3_PC"+tag,180,-360,360,180,-360,360,pcevent.pos.Phi()*180/M_PI,sx3event.pos.Phi()*180/M_PI,pmlabel);
plotter->Fill1D("pmiscs_dt_Anode_sx3_PC"+std::to_string(phicut)+tag,600,-2000,2000,pcevent.Time1-sx3event.Time1,pmlabel);
plotter->Fill1D("pmiscs_dt_Cathode_sx3"+tag,600,-2000,2000,pcevent.Time2-sx3event.Time1,pmlabel);
plotter->Fill2D("pmiscs_dt_Anode_E_sx3_PC"+std::to_string(phicut)+tag,600,-2000,2000,400,0,10,pcevent.Time1-sx3event.Time1,sx3event.Energy1,pmlabel);
plotter->Fill2D("pmiscs_dt_Cathode_E_sx3"+tag,600,-2000,2000,400,0,10,pcevent.Time2-sx3event.Time1,sx3event.Energy1,pmlabel);
plotter->Fill2D("pmiscs_dt_Cathodesx3_vsPCPhi"+tag,600,-2000,2000,180,-360,360,pcevent.Time2-sx3event.Time1,pcevent.pos.Phi()*180./M_PI,pmlabel);
plotter->Fill1D("pmiscs_pczfix"+tag,600,-300,300,pcz_fix,pmlabel);
plotter->Fill1D("pmiscs_pcz"+tag,600,-300,300,pcevent.pos.Z(),pmlabel);
//double path_length_q = (sx3event.pos-TVector3(0,0,vertex_z)).Mag()*0.1; plotter->Fill2D("pmiscs_Ef_vs_theta_sx3" + tag, 100, 0, 180, 800, 0, 20, theta_s * 180 / M_PI, sx3Efix, pmlabel);
//double path_length_s = (sx3event.pos-TVector3(0,0,vertex_z)).Mag()*0.1; plotter->Fill1D("pmiscs_VertexReconZ" + tag, 800, -400, 400, vertex_z, pmlabel);
double path_length_s = (sx3event.pos-r_rhoMin_fix).Mag()*0.1; plotter->Fill2D("pmiscs_VertexReconXY" + tag, 200, -100, 100, 200, -100, 100, r_rhoMin_fix.X(), r_rhoMin_fix.Y(), pmlabel);
double sx3Efix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(sx3event.Energy1)-path_length_s); plotter->Fill2D("pmiscs_VertexReconZ_vs_Ef" + tag, 800, -400, 400, 800, 0, 20, vertex_z, sx3Efix, pmlabel);
plotter->Fill2D("pmiscs_VertexReconZ_vs_Ef" + tag + "_a" + std::to_string(pcevent.multi1), 800, -400, 400, 800, 0, 20, vertex_z, sx3Efix, pmlabel);
//plotter->Fill2D("sx3Ef_sx3E_matrix_all"+tag,400,0,10,400,0,10,sx3Efix,sx3event.Energy1,pmlabel);
plotter->Fill2D("pmiscs_dE3_Ef_Anodesx3"+tag,400,0,10,400,0,40000,sx3Efix,pcevent.Energy1*sinTheta_customV*3,pmlabel);
plotter->Fill2D("pmiscs_dE3_Ef_Cathodesx3"+tag,400,0,10,400,0,10000,sx3Efix,pcevent.Energy2*sinTheta_customV,pmlabel);
plotter->Fill2D("pmiscs_Ef_vs_theta_sx3"+tag,100,0,180,800,0,20,theta_s*180/M_PI,sx3Efix,pmlabel);
plotter->Fill1D("pmiscs_VertexReconZ"+tag,800,-400,400,vertex_z,pmlabel);
plotter->Fill2D("pmiscs_VertexReconXY"+tag,200,-100,100,200,-100,100,r_rhoMin_fix.X(),r_rhoMin_fix.Y(),pmlabel);
plotter->Fill2D("pmiscs_VertexReconZ_vs_Ef"+tag,800,-400,400,800,0,20,vertex_z,sx3Efix,pmlabel);
plotter->Fill2D("pmiscs_VertexReconZ_vs_Ef"+tag+"_a"+std::to_string(pcevent.multi1),800,-400,400,800,0,20,vertex_z,sx3Efix,pmlabel);
}; };
plot_with_tag(); plot_with_tag();
if(cathode_alpha_select) if (cathode_alpha_select)
plot_with_tag("_cathode_alphas"); plot_with_tag("_cathode_alphas");
else else
plot_with_tag("_cathode_protons"); plot_with_tag("_cathode_protons");
//plotter->Fill1D("pmisc_Ex_from_protons",200,-10,10,apkin_p.getExc(sx3Efix,theta_s*180/M_PI),pmlabel); // plotter->Fill1D("pmisc_Ex_from_protons",200,-10,10,apkin_p.getExc(sx3Efix,theta_s*180/M_PI),pmlabel);
}//end PCEvents loop } // end PCEvents loop
}//end sx3Events loop } // end sx3Events loop
} }

1516
qqq_Calib.dat
View File

File diff suppressed because it is too large Load Diff

1
run_17F.sh
View File

@ -22,6 +22,7 @@ function run_once() {
export -f run_once export -f run_once
# run_once 350 # run_once 350
# parallel -j 6 --ctag run_once {1} ::: {325..400} # parallel -j 6 --ctag run_once {1} ::: {325..400}
parallel -j 6 --ctag run_once {1} ::: {351,353,355,358,359,360,362,367}
rm output_17F.root rm output_17F.root
hadd -j 4 -k output_17F.root 17F_output/results_run3*.root hadd -j 4 -k output_17F.root 17F_output/results_run3*.root

20
run_tr.sh
View File

@ -4,7 +4,9 @@ export flipa=0
export anode_offset=0 export anode_offset=0
export DATASET="27Al" export DATASET="27Al"
if [[ 1 -eq 0 ]]; then if [[ 1 -eq 0 ]]; then
root -b -q -l -x ../ANASEN_analysis/source.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root result.root;
#root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_009_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run09.root; #root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_009_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run09.root;
# exit
root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_001_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run01.root; root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_001_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run01.root;
root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_002_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run02.root; root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_002_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run02.root;
root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_003_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run03.root; root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_003_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run03.root;
@ -26,7 +28,7 @@ export timecut_low=400.0;
#export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_009_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run09.root; #export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_009_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run09.root;
#export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_010_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run10.root; #export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_010_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run10.root;
#export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_011_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run11.root; #export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_011_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run11.root;
export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_012_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run12.root; export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_012_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root source_results_run12.root;
# exit # exit
#export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_013_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run13.root; #export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_013_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run13.root;
#exit #exit
@ -36,14 +38,14 @@ unset timecut_low
#protons+gas, 27Al #protons+gas, 27Al
#export flip180="1" #export flip180="1"
#export flip180="0" #export flip180="0"
if [[ 1 -eq 1 ]]; then if [[ 1 -eq 0 ]]; then
export flipa=0 export flipa=0
export anode_offset=0 export anode_offset=0
export source_vertex=-200.0; #put the 'source' on the entrance window export source_vertex=-200.0; #put the 'source' on the entrance window
root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_018_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run18.root;
exit
root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_015_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run15.root;
root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_017_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run17.root; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_017_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run17.root;
exit
root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_018_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run18.root;
root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_015_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run15.root;
root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_019_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run19.root; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_019_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run19.root;
root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_020_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run20.root; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_020_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run20.root;
root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_021_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run21.root; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_021_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run21.root;
@ -79,10 +81,10 @@ fi
#17F alpha run with gas #17F alpha run with gas
if [[ 1 -eq 1 ]]; then if [[ 1 -eq 1 ]]; then
export source_vertex=53.44; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_018_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run18.root; export source_vertex=53.44; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_018_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root source_results_run18.root;
export source_vertex=14.24; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_019_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run19.root; export source_vertex=14.24; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_019_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root source_results_run19.root;
export source_vertex=-24.96; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_020_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run20.root; export source_vertex=-24.96; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_020_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root source_results_run20.root;
export source_vertex=-73.96; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_021_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run21.root; export source_vertex=-73.96; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_021_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root source_results_run21.root;
fi fi
#17F reaction data #17F reaction data
#export flip180="0" #export flip180="0"

171
scratch/make_prettyplots.C
View File

@ -3,12 +3,12 @@
// //
// Feed it a ROOT file and a histogram name, get a publication-quality PNG. // Feed it a ROOT file and a histogram name, get a publication-quality PNG.
// //
// Usage: // Usage (Single File):
// root -l -b -q 'make_pretty.C("myfile.root", "histName")' // root -l -b -q 'make_pretty.C("myfile.root", "histName")'
// root -l -b -q 'make_pretty.C("myfile.root", "histName", "x label", "y label")'
// //
// Supports TH1F, TH1D, TH2F, TH2D — type is detected automatically. // Usage (Multi-File Overlay):
// Output PNG is written to the same directory as the input file. // root -l -b -q 'make_pretty.C("file1.root,file2.root", "Run 1,Run 2", "histName", "X-Axis", "Y-Axis")'
//
// ============================================================================= // =============================================================================
#include "TFile.h" #include "TFile.h"
@ -18,11 +18,15 @@
#include "TStyle.h" #include "TStyle.h"
#include "TGaxis.h" #include "TGaxis.h"
#include "TLatex.h" #include "TLatex.h"
#include "TLegend.h"
#include "TSystem.h" #include "TSystem.h"
#include "TROOT.h" #include "TROOT.h"
#include "TObjArray.h"
#include "TObjString.h"
#include <iostream> #include <iostream>
#include <string> #include <string>
#include <vector>
// --------------------------------------------------------------------------- // ---------------------------------------------------------------------------
// Style — called once before anything is drawn // Style — called once before anything is drawn
@ -63,11 +67,18 @@ void SetStyle()
gStyle->SetNumberContours(255); gStyle->SetNumberContours(255);
} }
// ---------------------------------------------------------------------------
// Single Plot Function
// ---------------------------------------------------------------------------
void make_prettyplots(const char *rootFile, void make_prettyplots(const char *rootFile,
const char *histName, const char *histName,
const char *xlabel = "", const char *xlabel = "",
const char *ylabel = "", const char *ylabel = "",
double minZ = -9999.0) double minZ = -9999.0,
double minY = -9999.0,
double minX = -9999.0,
double maxX = -9999.0,
double maxY = -9999.0)
{ {
SetStyle(); SetStyle();
@ -121,8 +132,20 @@ void make_prettyplots(const char *rootFile,
// --- Apply Minimum Z if provided --- // --- Apply Minimum Z if provided ---
if (minZ != -9999.0) if (minZ != -9999.0)
{ {
h->GetXaxis()->SetRangeUser(minX, h->GetXaxis()->GetXmax());
h->GetYaxis()->SetRangeUser(minY, h->GetYaxis()->GetXmax());
h->SetMinimum(minZ); h->SetMinimum(minZ);
} }
if (minY != -9999.0)
h->GetYaxis()->SetRangeUser(minY, h->GetYaxis()->GetXmax());
if (minX != -9999.0)
h->GetXaxis()->SetRangeUser(minX, h->GetXaxis()->GetXmax());
if (maxX != -9999.0 && maxY != -9999.0)
{
h->GetXaxis()->SetRangeUser(minX, maxX);
h->GetYaxis()->SetRangeUser(minY, maxY);
}
h->GetXaxis()->SetTitle(xlabel); h->GetXaxis()->SetTitle(xlabel);
h->GetYaxis()->SetTitle(ylabel); h->GetYaxis()->SetTitle(ylabel);
@ -135,6 +158,11 @@ void make_prettyplots(const char *rootFile,
else else
{ {
TH1 *h = (TH1 *)clone; TH1 *h = (TH1 *)clone;
if (minX != -9999.0)
h->GetXaxis()->SetRangeUser(minX, h->GetXaxis()->GetXmax());
if (maxX != -9999.0)
h->GetXaxis()->SetRangeUser(minX, maxX);
h->SetStats(0); h->SetStats(0);
h->SetLineColor(kAzure - 5); h->SetLineColor(kAzure - 5);
h->SetLineWidth(3); h->SetLineWidth(3);
@ -148,21 +176,138 @@ void make_prettyplots(const char *rootFile,
} }
// --- Save --------------------------------------------------------------- // --- Save ---------------------------------------------------------------
// Build output path: same directory as input file, named after the histogram
std::string inPath(rootFile); std::string inPath(rootFile);
std::string dir = inPath.substr(0, inPath.find_last_of("/\\")); std::string dir = inPath.substr(0, inPath.find_last_of("/\\"));
if (dir == inPath) if (dir == inPath)
dir = "."; // no directory component — use cwd dir = ".";
std::string outPath = dir + "/" + std::string(histName) + ".png"; std::string safeHistName = histName;
// Replace any "/" inside histName (e.g. "folder/hist") with "_" for (char &ch : safeHistName)
for (char &ch : outPath)
if (ch == '/') if (ch == '/')
ch = '_'; ch = '_';
c.Modified(); std::string outPath = dir + "/" + safeHistName + ".png";
c.Update();
c.SaveAs(outPath.c_str()); c.SaveAs(outPath.c_str());
std::cout << "Saved: " << outPath << "\n"; std::cout << "Saved: " << outPath << "\n";
} }
// ---------------------------------------------------------------------------
// Multi-File Overlay Function
// ---------------------------------------------------------------------------
void make_prettyplots(TString filesCSV, TString labelsCSV, TString histName, TString xAxisLabel="", TString yAxisLabel="", double yMin=-9999, double yMax=-9999, double xMin=-9999, double xMax=-9999)
{
SetStyle();
// Parse the comma-separated lists natively using ROOT's TString
TObjArray* fileArr = filesCSV.Tokenize(",");
TObjArray* labelArr = labelsCSV.Tokenize(",");
if (fileArr->GetEntriesFast() != labelArr->GetEntriesFast())
{
std::cerr << "ERROR: Number of files does not match number of labels!\n";
return;
}
// Array of distinct colors to cycle through
int colors[] = {kBlack, kRed + 1, kAzure + 2, kGreen + 2, kMagenta + 1, kOrange + 7, kTeal - 1, kPink + 2};
int nColors = 8;
std::vector<TH1 *> hists;
double globalMaxY = -999999.0;
// Load all histograms and determine the maximum peak
for (int i = 0; i < fileArr->GetEntriesFast(); i++)
{
TString currentFile = ((TObjString*)fileArr->At(i))->GetString();
TString currentLabel = ((TObjString*)labelArr->At(i))->GetString();
TFile *f = TFile::Open(currentFile, "READ");
if (!f || f->IsZombie())
continue;
TH1 *h = (TH1 *)f->Get(histName);
if (!h)
{
std::cerr << "WARNING: '" << histName << "' not found in " << currentFile << ". Skipping...\n";
f->Close();
continue;
}
TH1 *clone = (TH1 *)h->Clone(Form("h_%d", i));
clone->SetDirectory(0); // Detach from file
f->Close();
// Apply X-range BEFORE checking the max Y-height, otherwise peaks
// outside the viewing range might scale the Y-axis unnecessarily!
if (xMin != -9999.0 && xMax != -9999.0)
clone->GetXaxis()->SetRangeUser(xMin, xMax);
else if (xMin != -9999.0)
clone->GetXaxis()->SetRangeUser(xMin, clone->GetXaxis()->GetXmax());
if (clone->GetMaximum() > globalMaxY)
{
globalMaxY = clone->GetMaximum();
}
hists.push_back(clone);
}
if (hists.empty())
{
std::cerr << "ERROR: No valid histograms found to plot.\n";
return;
}
// --- Draw ---
TCanvas c("c", "", 0, 0, 2100, 1575);
c.cd();
// Create Legend (Positioned top right)
TLegend *leg = new TLegend(0.65, 0.70, 0.92, 0.90);
leg->SetBorderSize(0);
leg->SetFillStyle(0); // Transparent background
for (size_t i = 0; i < hists.size(); i++)
{
hists[i]->SetStats(0);
hists[i]->SetLineColor(colors[i % nColors]);
hists[i]->SetLineWidth(3);
if (i == 0)
{
// First histogram dictates the frame layout
if (yMin != -9999.0 && yMax != -9999.0) {
hists[i]->GetYaxis()->SetRangeUser(yMin, yMax);
} else {
hists[i]->SetMaximum(globalMaxY * 1.15); // Add 15% headroom automatically
}
if (xAxisLabel != "") hists[i]->GetXaxis()->SetTitle(xAxisLabel);
if (yAxisLabel != "") hists[i]->GetYaxis()->SetTitle(yAxisLabel);
hists[i]->GetXaxis()->SetTitleOffset(1.1);
hists[i]->GetYaxis()->SetTitleOffset(1.4);
hists[i]->GetXaxis()->CenterTitle(true);
hists[i]->GetYaxis()->CenterTitle(true);
hists[i]->Draw("hist");
}
else
{
// Subsequent histograms stack on top
hists[i]->Draw("hist same");
}
TString currentLabel = ((TObjString*)labelArr->At(i))->GetString();
leg->AddEntry(hists[i], currentLabel.Data(), "l");
}
leg->Draw();
// --- Save ---
TString safeHistName = histName;
safeHistName.ReplaceAll("/", "_");
TString outPath = safeHistName + "_overlay.png";
c.SaveAs(outPath.Data());
std::cout << "Saved: " << outPath.Data() << "\n";
}