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

modified: Armory/PC_StepLadder_Correction.h

Browse Source 
modified:   TrackRecon.C
	modified:   run_tr.sh
	modified:   scan_slope_runs.sh
	modified:   scratch/plot_slope_scan.C
This commit is contained in:
Vignesh Sitaraman 2026-06-11 09:29:14 -04:00
parent 8dd6526221
commit 20f8ad227c
5 changed files with 554 additions and 383 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
Armory/PC_StepLadder_Correction.h
View File

@ -2,7 +2,7 @@
/*double model(double* x, double* p) { /*double model(double* x, double* p) {
double result = x[0]; double result = x[0];
double factor = 29.0; double factor = 29.0;
double slope = 0.7; double slope = 0.40;
if(TMath::Abs(x[0]) < 16.2) result=x[0]*slope; if(TMath::Abs(x[0]) < 16.2) result=x[0]*slope;
else if(TMath::Abs(x[0]) < 49.8 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor; else if(TMath::Abs(x[0]) < 49.8 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor;
else if(TMath::Abs(x[0]) < 85.2 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*2; else if(TMath::Abs(x[0]) < 85.2 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*2;
@ -12,7 +12,7 @@
double model_invert(double *y, double *q) { double model_invert(double *y, double *q) {
double result=y[0]; double result=y[0];
double slope = 0.7; double slope = 0.40;
double factor = 0.0; double factor = 0.0;
if(TMath::Abs(y[0]) < 16.2/slope) result = y[0]/slope; if(TMath::Abs(y[0]) < 16.2/slope) result = y[0]/slope;
else if(TMath::Abs(y[0]) < 49.8/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor; else if(TMath::Abs(y[0]) < 49.8/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor;
@ -23,7 +23,7 @@ double model_invert(double *y, double *q) {
double model_invert(double* y, double* p) { double model_invert(double* y, double* p) {
double result = y[0]; double result = y[0];
double slope = 0.6; double slope = 0.52;
double z_grid[8] = {147.998,101.946,59.7634,19.6965,-19.6965,-59.7634,-101.946,-147.998}; double z_grid[8] = {147.998,101.946,59.7634,19.6965,-19.6965,-59.7634,-101.946,-147.998};
for(int i=0;i<7;i++) { for(int i=0;i<7;i++) {
if(y[0] <= z_grid[i] && y[0] > z_grid[i+1]) { if(y[0] <= z_grid[i] && y[0] > z_grid[i+1]) {
@ -38,7 +38,7 @@ double model_invert(double* y, double* p) {
double model_a1c1(double* x, double* p) { double model_a1c1(double* x, double* p) {
double result = x[0]; double result = x[0];
double factor = 29.0; double factor = 29.0;
double slope = 0.0; double slope = 0.52;
if(TMath::Abs(x[0]) < 16.2) result=x[0]*slope; if(TMath::Abs(x[0]) < 16.2) result=x[0]*slope;
else if(TMath::Abs(x[0]) < 49.8 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor; else if(TMath::Abs(x[0]) < 49.8 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor;
else if(TMath::Abs(x[0]) < 85.2 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*2; else if(TMath::Abs(x[0]) < 85.2 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*2;
@ -48,7 +48,7 @@ double model_a1c1(double* x, double* p) {
double model_invert_a1c1(double *y, double *q) { double model_invert_a1c1(double *y, double *q) {
double result=y[0]; double result=y[0];
/* double slope = 1.0; /* double slope = 0.40;
double factor = 5.0; double factor = 5.0;
if(TMath::Abs(y[0]) < 16.2/slope) result = y[0]/slope; if(TMath::Abs(y[0]) < 16.2/slope) result = y[0]/slope;
else if(TMath::Abs(y[0]) < 49.8/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor; else if(TMath::Abs(y[0]) < 49.8/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor;

282
TrackRecon.C
View File

@ -45,7 +45,7 @@ bool doPCQQQClusterAnalysis = true;
bool doOldAnalysis = false; 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 = 105.0;
double z_entrance = -174.3 - 9.7 - 100.0; double z_entrance = -174.3 - 9.7 - 100.0;
const double anode_gain = 1.5146e-5; // channels --> MeV const double anode_gain = 1.5146e-5; // channels --> MeV
std::string dataset; std::string dataset;
@ -597,7 +597,9 @@ Bool_t TrackRecon::Process(Long64_t entry)
continue; continue;
// if(eRingMeV<1.2 || eWedgeMeV<1.2) continue; // if(eRingMeV<1.2 || eWedgeMeV<1.2) continue;
double theta = 2 * TMath::Pi() * (-qqq.id[i] * 16 + (15 - chWedge) + 0.5) / (16 * 4); // double theta = 2 * TMath::Pi() * (-qqq.id[i] * 16 + (15 - chWedge) + 0.5) / (16 * 4);
double theta = (M_PI/180.)*(-90*qqq.id[i]+(87./16.)*((15-chWedge)+0.5)+3.0);
double rho = 50. + (50. / 16.) * (chRing + 0.5); //"?" double rho = 50. + (50. / 16.) * (chRing + 0.5); //"?"
// z used to be 75+30+23=128 // z used to be 75+30+23=128
// we found a 12mm shift towards the vertex later --> 116 // we found a 12mm shift towards the vertex later --> 116
@ -1360,9 +1362,7 @@ void protonAlphaHistograms(HistPlotter *plotter, std::vector<Event> QQQ_Events,
return; return;
} }
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){
{
for (auto pcevent : PC_Events) for (auto pcevent : PC_Events)
{ {
bool PCSX3TimeCut = false; bool PCSX3TimeCut = false;
@ -1370,8 +1370,8 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
bool PCCSX3TimeCut = false; bool PCCSX3TimeCut = false;
for (auto sx3event : SX3_Events) for (auto sx3event : SX3_Events)
{ {
plotter->Fill1D("dt_pcA_sx3B" + std::to_string(sx3event.ch2), 640, -2000, 2000, sx3event.Time1 - pcevent.Time1, "hTiming"); plotter->Fill1D("dt_pcA_sx3B" + std::to_string(sx3event.ch2), 640, -2000, 2000, sx3event.Time1 - pcevent.Time1, "Timing");
plotter->Fill1D("dt_pcC_sx3B" + std::to_string(sx3event.ch2), 640, -2000, 2000, sx3event.Time1 - pcevent.Time2, "hTiming"); plotter->Fill1D("dt_pcC_sx3B" + std::to_string(sx3event.ch2), 640, -2000, 2000, sx3event.Time1 - pcevent.Time2, "Timing");
if (sx3event.Time1 - pcevent.Time1 < 0) //-150 for alphas if (sx3event.Time1 - pcevent.Time1 < 0) //-150 for alphas
PCASX3TimeCut = 1; PCASX3TimeCut = 1;
if (sx3event.Time1 - pcevent.Time2 < 0) //-200 for alphas if (sx3event.Time1 - pcevent.Time2 < 0) //-200 for alphas
@ -1380,35 +1380,35 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
bool phicut = sx3event.pos.Phi() <= pcevent.pos.Phi() + TMath::Pi() / 4. && sx3event.pos.Phi() >= pcevent.pos.Phi() - TMath::Pi() / 4.; bool phicut = sx3event.pos.Phi() <= pcevent.pos.Phi() + TMath::Pi() / 4. && sx3event.pos.Phi() >= pcevent.pos.Phi() - TMath::Pi() / 4.;
plotter->Fill1D("dt_pcA_sx3B", 640, -2000, 2000, sx3event.Time1 - pcevent.Time1); plotter->Fill1D("dt_pcA_sx3B", 640, -2000, 2000, sx3event.Time1 - pcevent.Time1, "Timing");
plotter->Fill1D("dt_pcC_sx3B", 640, -2000, 2000, sx3event.Time1 - pcevent.Time2); plotter->Fill1D("dt_pcC_sx3B", 640, -2000, 2000, sx3event.Time1 - pcevent.Time2, "Timing");
plotter->Fill2D("dt_pcA_vs_sx3RE", 640, -2000, 2000, 400, 0, 30, sx3event.Time1 - pcevent.Time1, sx3event.Energy1); plotter->Fill2D("dt_pcA_vs_sx3RE", 640, -2000, 2000, 400, 0, 30, sx3event.Time1 - pcevent.Time1, sx3event.Energy1, "Timing");
plotter->Fill2D("dE_E_Anodesx3B", 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1); plotter->Fill2D("dE_E_Anodesx3B", 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1, "PID_dE_E");
plotter->Fill2D("dE_E_Cathodesx3B", 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2); plotter->Fill2D("dE_E_Cathodesx3B", 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2, "PID_dE_E");
if (pcevent.multi1 == 1 && pcevent.multi2 == 2) if (pcevent.multi1 == 1 && pcevent.multi2 == 2)
plotter->Fill2D("dE_E_Anodesx3B_a1c2", 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1); plotter->Fill2D("dE_E_Anodesx3B_a1c2", 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1, "PID_dE_E");
if (pcevent.multi1 == 1 && pcevent.multi2 == 2) if (pcevent.multi1 == 1 && pcevent.multi2 == 2)
plotter->Fill2D("dE_E_Cathodesx3B_a1c2", 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2); plotter->Fill2D("dE_E_Cathodesx3B_a1c2", 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2, "PID_dE_E");
if (pcevent.multi1 == 2 && pcevent.multi2 == 1) if (pcevent.multi1 == 2 && pcevent.multi2 == 1)
plotter->Fill2D("dE_E_Anodesx3B_a2c1", 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1); plotter->Fill2D("dE_E_Anodesx3B_a2c1", 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1, "PID_dE_E");
if (pcevent.multi1 == 2 && pcevent.multi2 == 1) if (pcevent.multi1 == 2 && pcevent.multi2 == 1)
plotter->Fill2D("dE_E_Cathodesx3B_a2c1", 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2); plotter->Fill2D("dE_E_Cathodesx3B_a2c1", 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2, "PID_dE_E");
plotter->Fill2D("sx3phi_vs_pcphi" + std::to_string(sx3event.Time1 - pcevent.Time1 < -150), 100, -360, 360, 100, -360, 360, sx3event.pos.Phi() * 180 / M_PI, pcevent.pos.Phi() * 180 / M_PI); plotter->Fill2D("sx3phi_vs_pcphi" + std::to_string(sx3event.Time1 - pcevent.Time1 < -150), 100, -360, 360, 100, -360, 360, sx3event.pos.Phi() * 180 / M_PI, pcevent.pos.Phi() * 180 / M_PI, "Kinematics_Angles");
if (PCSX3TimeCut) if (PCSX3TimeCut)
{ {
plotter->Fill1D("dt_pcA_sx3B_timecut", 640, -2000, 2000, sx3event.Time1 - pcevent.Time1); plotter->Fill1D("dt_pcA_sx3B_timecut", 640, -2000, 2000, sx3event.Time1 - pcevent.Time1, "Timing");
plotter->Fill1D("dt_pcC_sx3B_timecut", 640, -2000, 2000, sx3event.Time1 - pcevent.Time2); plotter->Fill1D("dt_pcC_sx3B_timecut", 640, -2000, 2000, sx3event.Time1 - pcevent.Time2, "Timing");
plotter->Fill2D("xyplot_sx3" + std::to_string(sx3event.ch2 / 4), 100, -100, 100, 100, -100, 100, sx3event.pos.X(), sx3event.pos.Y()); plotter->Fill2D("xyplot_sx3" + std::to_string(sx3event.ch2 / 4), 100, -100, 100, 100, -100, 100, sx3event.pos.X(), sx3event.pos.Y(), "Vertex_Reconstruction");
plotter->Fill2D("xyplot_sx3" + std::to_string(sx3event.ch2 / 4), 100, -100, 100, 100, -100, 100, pcevent.pos.X(), pcevent.pos.Y()); plotter->Fill2D("xyplot_sx3" + std::to_string(sx3event.ch2 / 4), 100, -100, 100, 100, -100, 100, pcevent.pos.X(), pcevent.pos.Y(), "Vertex_Reconstruction");
plotter->Fill2D("pcz_vs_pcphi_TimeCut", 600, -200, 200, 120, -360, 360, pcevent.pos.Z(), pcevent.pos.Phi() * 180 / M_PI); // x-axis is all Si det, y-axis is PC anode+cathode only plotter->Fill2D("pcz_vs_pcphi_TimeCut", 600, -200, 200, 120, -360, 360, pcevent.pos.Z(), pcevent.pos.Phi() * 180 / M_PI, "Z_Reconstruction");
} }
double sx3rho = 88.0; // approximate barrel radius double sx3rho = 88.0;
double sx3z = sx3event.pos.Z(); // w.r.t target origin at 90 for run12 double sx3z = sx3event.pos.Z();
double pcz = pcevent.pos.Z(); double pcz = pcevent.pos.Z();
double calcsx3theta = TMath::ATan2(sx3rho - z_to_crossover_rho(pcz), sx3z - pcz); double calcsx3theta = TMath::ATan2(sx3rho - z_to_crossover_rho(pcz), sx3z - pcz);
plotter->Fill2D("dE2_E_Anodesx3B", 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1 * TMath::Sin(calcsx3theta)); plotter->Fill2D("dE2_E_Anodesx3B", 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1 * TMath::Sin(calcsx3theta), "PID_dE_E");
plotter->Fill2D("dE2_E_Cathodesx3B", 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2 * TMath::Sin(calcsx3theta)); plotter->Fill2D("dE2_E_Cathodesx3B", 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2 * TMath::Sin(calcsx3theta), "PID_dE_E");
double sx3theta = TMath::ATan2(sx3rho, sx3z - source_vertex); double sx3theta = TMath::ATan2(sx3rho, sx3z - source_vertex);
double pczguess = 37.0 / TMath::Tan(sx3theta) + source_vertex; double pczguess = 37.0 / TMath::Tan(sx3theta) + source_vertex;
@ -1419,78 +1419,79 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
TVector3 v = x2 - x1; TVector3 v = x2 - 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 vector_closest_to_z_sx3 = x1 + t_minimum * v; TVector3 vector_closest_to_z_sx3 = x1 + t_minimum * v;
plotter->Fill1D("VertexReconZ_SX3" + std::to_string(PCSX3TimeCut), 600, -1300, 1300, vector_closest_to_z_sx3.Z(), "hPCZSX3"); plotter->Fill1D("VertexReconZ_SX3" + std::to_string(PCSX3TimeCut), 600, -1300, 1300, vector_closest_to_z_sx3.Z(), "Vertex_Reconstruction");
plotter->Fill1D("VertexReconZ_SX3", 600, -1300, 1300, vector_closest_to_z_sx3.Z()); plotter->Fill1D("VertexReconZ_SX3", 600, -1300, 1300, vector_closest_to_z_sx3.Z(), "Vertex_Reconstruction");
plotter->Fill2D("VertexReconXY_SX3" + std::to_string(PCSX3TimeCut), 100, -100, 100, 100, -100, 100, vector_closest_to_z_sx3.X(), vector_closest_to_z_sx3.Y(), "hPCZSX3"); plotter->Fill2D("VertexReconXY_SX3" + std::to_string(PCSX3TimeCut), 100, -100, 100, 100, -100, 100, vector_closest_to_z_sx3.X(), vector_closest_to_z_sx3.Y(), "Vertex_Reconstruction");
plotter->Fill2D("pcz_vs_time", 2000, 0, 2000, 600, -200, 200, pcevent.Time1 * 1e-9, pcevent.pos.Z()); // x-axis is all Si det, y-axis is PC anode+cathode only plotter->Fill2D("pcz_vs_time", 2000, 0, 2000, 600, -200, 200, pcevent.Time1 * 1e-9, pcevent.pos.Z(), "Timing");
plotter->Fill2D("pcphi_vs_time", 2000, 0, 2000, 180, -360, 360, pcevent.Time1 * 1e-9, pcevent.pos.Phi() * 180. / M_PI); // x-axis is all Si det, y-axis is PC anode+cathode only plotter->Fill2D("pcphi_vs_time", 2000, 0, 2000, 180, -360, 360, pcevent.Time1 * 1e-9, pcevent.pos.Phi() * 180. / M_PI, "Timing");
// plotter->Fill2D("pcz_vs_time_strip"+std::to_string(sx3event.ch2),2000,0,2000,600,-200,200,pcevent.Time1*1e-9,pcevent.pos.Z()); //x-axis is all Si det, y-axis is PC anode+cathode only plotter->Fill2D("sx3phi_vs_time", 2000, 0, 2000, 180, -360, 360, pcevent.Time1 * 1e-9, sx3event.pos.Phi() * 180. / M_PI, "Timing");
plotter->Fill2D("sx3phi_vs_time", 2000, 0, 2000, 180, -360, 360, pcevent.Time1 * 1e-9, sx3event.pos.Phi() * 180. / M_PI); // x-axis is all Si det, y-axis is PC anode+cathode only
plotter->Fill2D("pcz_vs_sx3pczguess", 600, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z()); // x-axis is all Si det, y-axis is PC anode+cathode only plotter->Fill2D("pcz_vs_sx3pczguess", 600, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z(), "Z_Reconstruction");
if (pcevent.multi1 == 1 && pcevent.multi2 == 2) if (pcevent.multi1 == 1 && pcevent.multi2 == 2)
{ {
// if(pcevent.multi1==1) { plotter->Fill2D("pcz_vs_sx3pczguess_A1C2", 600, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z(), "Z_Reconstruction");
plotter->Fill2D("pcz_vs_sx3pczguess_A1C2", 600, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z());
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 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;
plotter->Fill1D("VertexRecon_pczfix_sx3", 800, -300, 300, r_rhoMin_fix.Z()); 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()); plotter->Fill1D("VertexRecon_pczfix", 800, -300, 300, r_rhoMin_fix.Z(), "Vertex_Reconstruction");
plotter->Fill1D("pczfix_A1C2_1d_sx3", 600, -200, 200, pcz_fix); plotter->Fill1D("pczfix_A1C2_1d_sx3", 600, -200, 200, pcz_fix, "Z_Reconstruction");
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, "Z_Reconstruction");
plotter->Fill2D("pczfix_residual_vs_pczguess_A1C2", 600, -200, 200, 200, -100, 100, pczguess, pcz_fix - pczguess); plotter->Fill2D("pczfix_vs_sx3pczguess_int_A1C2", 600, -200, 200, 600, -200, 200, pcz_guess_int, pcz_fix, "Z_Reconstruction");
plotter->Fill2D("pczfix_residual_vs_phi_A1C2", 200, 0, 6.28, 200, -100, 100, r_rhoMin_fix.Phi(), pcz_fix - pczguess); plotter->Fill2D("pczguess_vs_int", 600, -200, 200, 600, -200, 200, pcz_guess_int, pczguess, "Z_Reconstruction");
plotter->Fill1D("pczfix-sx3pczguess_A1C2", 200, -100, 100, pcz_fix - pczguess); plotter->Fill1D("pczguess_vs_int_residualsx3", 200, -50, 50, pcz_guess_int - pczguess, "Residuals");
plotter->Fill2D("pczfix_vs_sx3pczguess_A1C2_strip" + std::to_string(sx3event.ch2), 300, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z()); plotter->Fill2D("pczfix_residual_vs_pczguess_A1C2", 600, -200, 200, 200, -100, 100, pczguess, pcz_fix - pczguess, "Residuals");
plotter->Fill2D("pczfix_residual_vs_phi_A1C2", 200, 0, 6.28, 200, -100, 100, r_rhoMin_fix.Phi(), pcz_fix - pczguess, "Residuals");
plotter->Fill2D("pczguess_vs_int_residual_vs_phi_A1C2", 200, 0, 6.28, 200, -100, 100, r_rhoMin_fix.Phi(), pcz_guess_int - pczguess, "Residuals");
plotter->Fill1D("pczfix-sx3pczguess_A1C2", 200, -100, 100, pcz_fix - pczguess, "Residuals");
plotter->Fill2D("pczfix_vs_sx3pczguess_A1C2_strip" + std::to_string(sx3event.ch2), 300, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z(), "Z_Reconstruction");
double sinTheta_customV = TMath::Sin((sx3event.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Theta()); double sinTheta_customV = TMath::Sin((sx3event.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Theta());
plotter->Fill2D("dE3_E_CathodeSX3_A1C2_TC" + std::to_string(PCSX3TimeCut) + "_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2 * sinTheta_customV); plotter->Fill2D("dE3_E_CathodeSX3_A1C2_TC" + std::to_string(PCSX3TimeCut) + "_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2 * sinTheta_customV, "PID_dE_E");
plotter->Fill2D("dE3_E_AnodeSX3_A1C2_TC" + std::to_string(PCSX3TimeCut) + "_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1 * sinTheta_customV); plotter->Fill2D("dE3_E_AnodeSX3_A1C2_TC" + std::to_string(PCSX3TimeCut) + "_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1 * sinTheta_customV, "PID_dE_E");
if (TMath::Abs(r_rhoMin_fix.Z()) < 200.0) if (TMath::Abs(r_rhoMin_fix.Z()) < 200.0)
{ {
plotter->Fill2D("dE3_E_AnodeSX3B_A1C2_(vertex_fix_z/100)=" + std::to_string(floor(r_rhoMin_fix.Z() / 100.0)), 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1 * sinTheta_customV); plotter->Fill2D("dE3_E_AnodeSX3B_A1C2_(vertex_fix_z/100)=" + std::to_string(floor(r_rhoMin_fix.Z() / 100.0)), 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1 * sinTheta_customV, "PID_dE_E");
plotter->Fill2D("dE3_E_CathodeSX3B_A1C2_(vertex_fix_z/100)=" + std::to_string(floor(r_rhoMin_fix.Z() / 100.0)), 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2 * sinTheta_customV); plotter->Fill2D("dE3_E_CathodeSX3B_A1C2_(vertex_fix_z/100)=" + std::to_string(floor(r_rhoMin_fix.Z() / 100.0)), 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2 * sinTheta_customV, "PID_dE_E");
} }
if (pcevent.multi1 == 1 && pcevent.multi2 == 3) if (pcevent.multi1 == 1 && pcevent.multi2 == 3)
{ {
plotter->Fill2D("pcz_vs_sx3pczguess_A1C3", 600, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z()); plotter->Fill2D("pcz_vs_sx3pczguess_A1C3", 600, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z(), "Z_Reconstruction");
plotter->Fill2D("pcz_vs_sx3pczguess_A1C3_strip" + std::to_string(sx3event.ch2), 300, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z()); plotter->Fill2D("pcz_vs_sx3pczguess_A1C3_strip" + std::to_string(sx3event.ch2), 300, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z(), "Z_Reconstruction");
} }
plotter->Fill2D("pcz_vs_sx3pczguess_int", 600, -200, 200, 600, -200, 200, pcz_guess_int, pcevent.pos.Z()); // x-axis is all Si det, y-axis is PC anode+cathode only plotter->Fill2D("pcz_vs_sx3pczguess_int", 600, -200, 200, 600, -200, 200, pcz_guess_int, pcevent.pos.Z(), "Z_Reconstruction");
plotter->Fill2D("pcz_vs_sx3pczguess_strip" + std::to_string(sx3event.ch2), 300, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z()); plotter->Fill2D("pcz_vs_sx3pczguess_strip" + std::to_string(sx3event.ch2), 300, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z(), "Z_Reconstruction");
bool sx3PhiCut = (TMath::Abs(sx3event.pos.Phi() - pcevent.pos.Phi()) < 45.0 * M_PI / 180.); bool sx3PhiCut = (TMath::Abs(sx3event.pos.Phi() - pcevent.pos.Phi()) < 45.0 * M_PI / 180.);
plotter->Fill1D("pcz_sx3Coinc_phiCut" + std::to_string(sx3PhiCut) + "_TC" + std::to_string(PCSX3TimeCut), 300, 0, 200, sx3z); plotter->Fill1D("pcz_sx3Coinc_phiCut" + std::to_string(sx3PhiCut) + "_TC" + std::to_string(PCSX3TimeCut), 300, 0, 200, sx3z, "Z_Reconstruction");
plotter->Fill2D("pcz_vs_sx3z_phiCut" + std::to_string(sx3PhiCut) + "_TC" + std::to_string(PCSX3TimeCut), 300, 0, 200, 600, -400, 400, sx3z, pcevent.pos.Z()); plotter->Fill2D("pcz_vs_sx3z_phiCut" + std::to_string(sx3PhiCut) + "_TC" + std::to_string(PCSX3TimeCut), 300, 0, 200, 600, -400, 400, sx3z, pcevent.pos.Z(), "Z_Reconstruction");
// plotter->Fill2D("sx3E_vs_sx3z"+std::to_string(sx3event.ch2),400,0,30,300,0,200,sx3event.Energy1,sx3z); plotter->Fill2D("sx3E_vs_sx3z", 400, 0, 30, 300, 0, 200, sx3event.Energy1, sx3z, "Kinematics_Angles");
plotter->Fill2D("sx3E_vs_sx3z", 400, 0, 30, 300, 0, 200, sx3event.Energy1, sx3z);
plotter->Fill2D("pcdEA_vs_sx3z", 800, 0, 20000, 300, 0, 200, pcevent.Energy1, sx3z); plotter->Fill2D("pcdEA_vs_sx3z", 800, 0, 20000, 300, 0, 200, pcevent.Energy1, sx3z, "Kinematics_Angles");
plotter->Fill2D("pcdEC_vs_sx3z", 800, 0, 20000, 300, 0, 200, pcevent.Energy2, sx3z); plotter->Fill2D("pcdEC_vs_sx3z", 800, 0, 20000, 300, 0, 200, pcevent.Energy2, sx3z, "Kinematics_Angles");
plotter->Fill2D("pcdEA_vs_sx3z" + std::to_string(sx3event.ch2), 800, 0, 20000, 300, 0, 200, pcevent.Energy1, sx3z, "pcE_vs_sx3pos"); plotter->Fill2D("pcdEA_vs_sx3z" + std::to_string(sx3event.ch2), 800, 0, 20000, 300, 0, 200, pcevent.Energy1, sx3z, "Kinematics_Angles");
plotter->Fill2D("pcdEC_vs_sx3z" + std::to_string(sx3event.ch2), 800, 0, 20000, 300, 0, 200, pcevent.Energy2, sx3z, "pcE_vs_sx3pos"); plotter->Fill2D("pcdEC_vs_sx3z" + std::to_string(sx3event.ch2), 800, 0, 20000, 300, 0, 200, pcevent.Energy2, sx3z, "Kinematics_Angles");
plotter->Fill2D("pcdE2A_vs_sx3z", 800, 0, 20000, 300, 0, 200, pcevent.Energy1 * sinTheta, sx3z); plotter->Fill2D("pcdE2A_vs_sx3z", 800, 0, 20000, 300, 0, 200, pcevent.Energy1 * sinTheta, sx3z, "Kinematics_Angles");
plotter->Fill2D("pcdE2C_vs_sx3z", 800, 0, 20000, 300, 0, 200, pcevent.Energy2 * sinTheta, sx3z); plotter->Fill2D("pcdE2C_vs_sx3z", 800, 0, 20000, 300, 0, 200, pcevent.Energy2 * sinTheta, sx3z, "Kinematics_Angles");
plotter->Fill2D("phi_vs_stripnum", 180, -180, 180, 48, 0, 48, pcevent.pos.Phi() * 180. / M_PI, sx3event.ch2); plotter->Fill2D("phi_vs_stripnum", 180, -180, 180, 48, 0, 48, pcevent.pos.Phi() * 180. / M_PI, sx3event.ch2, "Kinematics_Angles");
plotter->Fill2D("E_theta_AnodeSX3", 400, -20, 180, 300, 0, 15, sx3theta * 180 / M_PI, sx3event.Energy1); plotter->Fill2D("E_theta_AnodeSX3", 400, -20, 180, 300, 0, 15, sx3theta * 180 / M_PI, sx3event.Energy1, "Kinematics_Angles");
} }
if (PCSX3TimeCut) if (PCSX3TimeCut)
{ {
plotter->Fill1D("PCZ_sx3", 800, -200, 200, pcevent.pos.Z(), "hPCZSX3"); plotter->Fill1D("PCZ_sx3", 800, -200, 200, pcevent.pos.Z(), "Z_Reconstruction");
} }
} // end PC-SX3 coincidence }
} }
} }
@ -1500,11 +1501,11 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
{ {
for (auto qqqevent : QQQ_Events) for (auto qqqevent : QQQ_Events)
{ {
plotter->Fill1D("dt_pcA_qqqR", 640, -2000, 2000, qqqevent.Time1 - pcevent.Time1); plotter->Fill1D("dt_pcA_qqqR", 640, -2000, 2000, qqqevent.Time1 - pcevent.Time1, "Timing");
plotter->Fill2D("dt_pcA_qqqR_vs_qqqRE", 640, -2000, 2000, 400, 0, 30, qqqevent.Time1 - pcevent.Time1, qqqevent.Energy1); plotter->Fill2D("dt_pcA_qqqR_vs_qqqRE", 640, -2000, 2000, 400, 0, 30, qqqevent.Time1 - pcevent.Time1, qqqevent.Energy1, "Timing");
plotter->Fill1D("dt_pcC_qqqW", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2); plotter->Fill1D("dt_pcC_qqqW", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2, "Timing");
plotter->Fill2D("phiPC_vs_phiQQQ", 180, -360, 360, 180, -360, 360, qqqevent.pos.Phi() * 180 / M_PI, pcevent.pos.Phi() * 180 / M_PI); plotter->Fill2D("phiPC_vs_phiQQQ", 180, -360, 360, 180, -360, 360, qqqevent.pos.Phi() * 180 / M_PI, pcevent.pos.Phi() * 180 / M_PI, "Kinematics_Angles");
double sinTheta = TMath::Sin((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()); /// TMath::Sin((TVector3(51.5,0,128.) - TVector3(0,0,85)).Theta()); double sinTheta = TMath::Sin((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta());
TVector3 x2(pcevent.pos); TVector3 x2(pcevent.pos);
TVector3 x1(qqqevent.pos); TVector3 x1(qqqevent.pos);
@ -1512,74 +1513,67 @@ void PCQQQClusterAnalysis(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 = x1 + t_minimum * v; TVector3 r_rhoMin = x1 + t_minimum * v;
// bool timecut = (qqqevent.Time1 - pcevent.Time1 < -150);
bool timecut = (qqqevent.Time1 - pcevent.Time1 < 150); bool timecut = (qqqevent.Time1 - pcevent.Time1 < 150);
bool lowercut_cath = pcevent.Energy2 * sinTheta < 250 && (qqqevent.Energy2 < 5.0 || qqqevent.Energy1 < 5.0); bool lowercut_cath = pcevent.Energy2 * sinTheta < 250 && (qqqevent.Energy2 < 5.0 || qqqevent.Energy1 < 5.0);
bool phicut = qqqevent.pos.Phi() <= pcevent.pos.Phi() + TMath::Pi() / 4. && qqqevent.pos.Phi() >= pcevent.pos.Phi() - TMath::Pi() / 4.; bool phicut = qqqevent.pos.Phi() <= pcevent.pos.Phi() + TMath::Pi() / 4. && qqqevent.pos.Phi() >= pcevent.pos.Phi() - TMath::Pi() / 4.;
if (lowercut_cath && phicut) if (lowercut_cath && phicut)
{ {
plotter->Fill1D("dt_pcA_qqqR_pidlow_PC1", 640, -2000, 2000, qqqevent.Time1 - pcevent.Time1); plotter->Fill1D("dt_pcA_qqqR_pidlow_PC1", 640, -2000, 2000, qqqevent.Time1 - pcevent.Time1, "Timing");
plotter->Fill2D("dt_pcA_qqqR_vs_qqqRE_pidlow_PC1", 640, -2000, 2000, 400, 0, 30, qqqevent.Time1 - pcevent.Time1, qqqevent.Energy1); plotter->Fill2D("dt_pcA_qqqR_vs_qqqRE_pidlow_PC1", 640, -2000, 2000, 400, 0, 30, qqqevent.Time1 - pcevent.Time1, qqqevent.Energy1, "Timing");
plotter->Fill1D("dt_pcC_qqqW_pidlow_PC1", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2); plotter->Fill1D("dt_pcC_qqqW_pidlow_PC1", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2, "Timing");
} }
if (timecut) if (timecut)
{ // && qqqevent.pos.Phi() <= pcevent.pos.Phi()+TMath::Pi()/4. && qqqevent.pos.Phi() >= pcevent.pos.Phi()-TMath::Pi()/4. ) { {
plotter->Fill2D("dE_E_AnodeQQQR", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E");
plotter->Fill2D("dE_E_AnodeQQQR", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1); plotter->Fill2D("dE_E_CathodeQQQR", 400, 0, 30, 800, 0, 10000, qqqevent.Energy2, pcevent.Energy2, "PID_dE_E");
plotter->Fill2D("dE_E_CathodeQQQR", 400, 0, 30, 800, 0, 10000, qqqevent.Energy2, pcevent.Energy2);
if (pcevent.multi1 == 1 && pcevent.multi2 == 2) if (pcevent.multi1 == 1 && pcevent.multi2 == 2)
plotter->Fill2D("dE_E_AnodeQQQR_a1c2", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1); plotter->Fill2D("dE_E_AnodeQQQR_a1c2", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E");
if (pcevent.multi1 == 1 && pcevent.multi2 == 2) if (pcevent.multi1 == 1 && pcevent.multi2 == 2)
plotter->Fill2D("dE_E_CathodeQQQR_a1c2", 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2); plotter->Fill2D("dE_E_CathodeQQQR_a1c2", 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2, "PID_dE_E");
if (pcevent.multi1 == 2 && pcevent.multi2 == 1) if (pcevent.multi1 == 2 && pcevent.multi2 == 1)
plotter->Fill2D("dE_E_AnodeQQQR_a2c1", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1); plotter->Fill2D("dE_E_AnodeQQQR_a2c1", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E");
if (pcevent.multi1 == 2 && pcevent.multi2 == 1) if (pcevent.multi1 == 2 && pcevent.multi2 == 1)
plotter->Fill2D("dE_E_CathodeQQQR_a2c1", 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2); plotter->Fill2D("dE_E_CathodeQQQR_a2c1", 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2, "PID_dE_E");
if (phicut) if (phicut)
{ {
plotter->Fill2D("dE2_E_AnodeQQQR_TC1PC1_pidlow" + std::to_string(lowercut_cath), 400, 0, 30, 800, 0, 4000, qqqevent.Energy1, pcevent.Energy1 * sinTheta); plotter->Fill2D("dE2_E_AnodeQQQR_TC1PC1_pidlow" + std::to_string(lowercut_cath), 400, 0, 30, 800, 0, 4000, qqqevent.Energy1, pcevent.Energy1 * sinTheta, "PID_dE_E");
plotter->Fill2D("dE2_E_CathodeQQQW_TC1PC1_pidlow" + std::to_string(lowercut_cath), 400, 0, 30, 800, 0, 1000, qqqevent.Energy2, pcevent.Energy2 * sinTheta); plotter->Fill2D("dE2_E_CathodeQQQW_TC1PC1_pidlow" + std::to_string(lowercut_cath), 400, 0, 30, 800, 0, 1000, qqqevent.Energy2, pcevent.Energy2 * sinTheta, "PID_dE_E");
// plotter->Fill2D("E_theta_AnodeQQQR_TC1PC1_pidlow"+std::to_string(lowercut_cath),75,0,90,300,0,15,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,qqqevent.Energy1); plotter->Fill2D("E_theta_zoomin_AnodeQQQR_TC1PC1_pidlow" + std::to_string(lowercut_cath), 60, 0, 30, 300, 0, 15, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, qqqevent.Energy1, "Kinematics_Angles");
plotter->Fill2D("E_theta_zoomin_AnodeQQQR_TC1PC1_pidlow" + std::to_string(lowercut_cath), 60, 0, 30, 300, 0, 15, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, qqqevent.Energy1);
} }
plotter->Fill2D("dE2_E_AnodeQQQR_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 4000, qqqevent.Energy1, pcevent.Energy1 * sinTheta); plotter->Fill2D("dE2_E_AnodeQQQR_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 4000, qqqevent.Energy1, pcevent.Energy1 * sinTheta, "PID_dE_E");
plotter->Fill2D("dE2_E_CathodeQQQR_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 1000, qqqevent.Energy2, pcevent.Energy2 * sinTheta); plotter->Fill2D("dE2_E_CathodeQQQR_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 1000, qqqevent.Energy2, pcevent.Energy2 * sinTheta, "PID_dE_E");
plotter->Fill2D("dEC_vs_dEA_TC1_PC" + std::to_string(phicut), 800, 0, 40000, 800, 0, 10000, pcevent.Energy1, pcevent.Energy2); plotter->Fill2D("dEC_vs_dEA_TC1_PC" + std::to_string(phicut), 800, 0, 40000, 800, 0, 10000, pcevent.Energy1, pcevent.Energy2, "PID_dE_E");
plotter->Fill2D("qqqphi_vs_time", 2000, 0, 2000, 180, -360, 360, pcevent.Time1 * 1e-9, qqqevent.pos.Phi() * 180. / M_PI); // x-axis is all Si det, y-axis is PC anode+cathode only plotter->Fill2D("qqqphi_vs_time", 2000, 0, 2000, 180, -360, 360, pcevent.Time1 * 1e-9, qqqevent.pos.Phi() * 180. / M_PI, "Timing");
plotter->Fill1D("dt_pcA_qqqR_timecut", 640, -2000, 2000, qqqevent.Time1 - pcevent.Time1); plotter->Fill1D("dt_pcA_qqqR_timecut", 640, -2000, 2000, qqqevent.Time1 - pcevent.Time1, "Timing");
plotter->Fill1D("dt_pcC_qqqW_timecut", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2); plotter->Fill1D("dt_pcC_qqqW_timecut", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2, "Timing");
plotter->Fill2D("dE_theta_AnodeQQQR", 90, 0, 90, 400, 0, 20000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy1); plotter->Fill2D("dE_theta_AnodeQQQR", 90, 0, 90, 400, 0, 20000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy1, "Kinematics_Angles");
plotter->Fill2D("dE2_theta_AnodeQQQR_zoomin", 60, 0, 30, 400, 0, 5000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy1 * sinTheta); plotter->Fill2D("dE2_theta_AnodeQQQR_zoomin", 60, 0, 30, 400, 0, 5000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy1 * sinTheta, "Kinematics_Angles");
plotter->Fill2D("dE2_theta_AnodeQQQR", 90, 0, 90, 400, 0, 20000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy1 * sinTheta); plotter->Fill2D("dE2_theta_AnodeQQQR", 90, 0, 90, 400, 0, 20000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy1 * sinTheta, "Kinematics_Angles");
plotter->Fill2D("phiPC_vs_phiQQQ_TimeCut", 180, -360, 360, 180, -360, 360, qqqevent.pos.Phi() * 180 / M_PI, pcevent.pos.Phi() * 180 / M_PI); plotter->Fill2D("phiPC_vs_phiQQQ_TimeCut", 180, -360, 360, 180, -360, 360, qqqevent.pos.Phi() * 180 / M_PI, pcevent.pos.Phi() * 180 / M_PI, "Kinematics_Angles");
// plotter->Fill2D("E_theta_AnodeQQQR_TC1_PC"+std::to_string(phicut),75,0,90,300,0,15,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,qqqevent.Energy1); plotter->Fill2D("Etot2_theta_AnodeQQQR", 75, 0, 90, 300, 0, 15, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, qqqevent.Energy1 + pcevent.Energy1 * anode_gain * sinTheta, "Kinematics_Angles");
// plotter->Fill2D("E_theta_zoomin_AnodeQQQR_TC1_PC"+std::to_string(phicut),60,0,30,300,0,15,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,qqqevent.Energy1);
// plotter->Fill2D("E2_theta_AnodeQQQR",75,0,90,300,0,15,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,qqqevent.Energy1);
plotter->Fill2D("Etot2_theta_AnodeQQQR", 75, 0, 90, 300, 0, 15, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, qqqevent.Energy1 + pcevent.Energy1 * anode_gain * sinTheta);
plotter->Fill2D("dE_theta_CathodeQQQR", 75, 0, 90, 800, 0, 10000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy2); plotter->Fill2D("dE_theta_CathodeQQQR", 75, 0, 90, 800, 0, 10000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy2, "Kinematics_Angles");
plotter->Fill2D("dE2_theta_CathodeQQQR", 75, 0, 90, 800, 0, 10000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy2 * sinTheta); plotter->Fill2D("dE2_theta_CathodeQQQR", 75, 0, 90, 800, 0, 10000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy2 * sinTheta, "Kinematics_Angles");
plotter->Fill2D("dE2_theta_CathodeQQQR_zoomin", 60, 0, 30, 800, 0, 3000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy2 * sinTheta); plotter->Fill2D("dE2_theta_CathodeQQQR_zoomin", 60, 0, 30, 800, 0, 3000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy2 * sinTheta, "Kinematics_Angles");
plotter->Fill2D("dE_phi_AnodeQQQR", 100, -180, 180, 800, 0, 40000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Phi() * 180 / M_PI, pcevent.Energy1); plotter->Fill2D("dE_phi_AnodeQQQR", 100, -180, 180, 800, 0, 40000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Phi() * 180 / M_PI, pcevent.Energy1, "Kinematics_Angles");
plotter->Fill2D("dE_phi_CathodeQQQR", 100, -180, 180, 800, 0, 10000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Phi() * 180 / M_PI, pcevent.Energy2); plotter->Fill2D("dE_phi_CathodeQQQR", 100, -180, 180, 800, 0, 10000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Phi() * 180 / M_PI, pcevent.Energy2, "Kinematics_Angles");
plotter->Fill1D("PCZ", 800, -200, 200, pcevent.pos.Z(), "phicut"); plotter->Fill1D("PCZ", 800, -200, 200, pcevent.pos.Z(), "Z_Reconstruction");
// plotter->Fill1D("PCZ_phicut_a"+std::to_string(aClusters.at(0).size())+"_c"+std::to_string(cClusters.at(0).size()),800,-200,200,pcevent.pos.Z(),"wiremult");
double pcz_guess_37 = 37. / TMath::Tan((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()) + source_vertex; double pcz_guess_37 = 37. / TMath::Tan((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()) + source_vertex;
plotter->Fill2D("pczguess_vs_pc_37", 180, 0, 200, 150, 0, 200, pcz_guess_37, pcevent.pos.Z(), "phicut"); plotter->Fill2D("pczguess_vs_pc_37", 180, 0, 200, 150, 0, 200, pcz_guess_37, pcevent.pos.Z(), "Z_Reconstruction");
double pcz_guess_42 = 42. / TMath::Tan((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()) + source_vertex; double pcz_guess_42 = 42. / TMath::Tan((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()) + source_vertex;
plotter->Fill2D("pczguess_vs_pc_42", 180, 0, 200, 150, 0, 200, pcz_guess_42, pcevent.pos.Z(), "phicut"); plotter->Fill2D("pczguess_vs_pc_42", 180, 0, 200, 150, 0, 200, pcz_guess_42, pcevent.pos.Z(), "Z_Reconstruction");
double pcz_guess_int = z_to_crossover_rho(pcevent.pos.Z()) / TMath::Tan((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()) + source_vertex; double pcz_guess_int = z_to_crossover_rho(pcevent.pos.Z()) / TMath::Tan((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()) + source_vertex;
// plotter->Fill2D("pczguess_vs_pc_int",180,0,200,150,0,200,pcz_guess_int,pcevent.pos.Z(),"phicut"); plotter->Fill2D("pczguess_vs_pc_int", 400, -200, 200, 600, -400, 400, pcz_guess_int, pcevent.pos.Z(), "Z_Reconstruction");
plotter->Fill2D("pczguess_vs_pc_int", 400, -200, 200, 600, -400, 400, pcz_guess_int, pcevent.pos.Z(), "phicut");
if (pcevent.multi1 == 1 && pcevent.multi2 == 2) if (pcevent.multi1 == 1 && pcevent.multi2 == 2)
{ {
double pcz_fix = pcfix_func.Eval(pcevent.pos.Z()); double pcz_fix = pcfix_func.Eval(pcevent.pos.Z());
@ -1589,69 +1583,69 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
TVector3 r_rhoMin_fix = x1 + t_minimum * v; TVector3 r_rhoMin_fix = x1 + t_minimum * v;
double sinTheta_customV = TMath::Sin((qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Theta()); double sinTheta_customV = TMath::Sin((qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Theta());
plotter->Fill2D("dE3_E_CathodeQQQW_A1C2_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 10000, qqqevent.Energy2, pcevent.Energy2 * sinTheta_customV); plotter->Fill2D("dE3_E_CathodeQQQW_A1C2_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 10000, qqqevent.Energy2, pcevent.Energy2 * sinTheta_customV, "PID_dE_E");
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, "PID_dE_E");
plotter->Fill1D("VertexRecon_pczfix_qqq", 800, -400, 400, r_rhoMin_fix.Z()); plotter->Fill1D("VertexRecon_pczfix_qqq", 800, -400, 400, r_rhoMin_fix.Z(), "Vertex_Reconstruction");
plotter->Fill1D("VertexRecon_pczfix", 800, -400, 400, r_rhoMin_fix.Z()); plotter->Fill1D("VertexRecon_pczfix", 800, -400, 400, r_rhoMin_fix.Z(), "Vertex_Reconstruction");
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(), "Vertex_Reconstruction");
if (TMath::Abs(r_rhoMin_fix.Z()) < 200.0) if (TMath::Abs(r_rhoMin_fix.Z()) < 200.0)
{ {
plotter->Fill2D("dE3_E_AnodeQQQR_A1C2_(vertex_fix_z/100)=" + std::to_string(floor(r_rhoMin_fix.Z() / 100.0)), 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1 * sinTheta_customV); plotter->Fill2D("dE3_E_AnodeQQQR_A1C2_(vertex_fix_z/100)=" + std::to_string(floor(r_rhoMin_fix.Z() / 100.0)), 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1 * sinTheta_customV, "PID_dE_E");
plotter->Fill2D("dE3_E_CathodeQQQR_A1C2_(vertex_fix_z/100)=" + std::to_string(floor(r_rhoMin_fix.Z() / 100.0)), 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2 * sinTheta_customV); plotter->Fill2D("dE3_E_CathodeQQQR_A1C2_(vertex_fix_z/100)=" + std::to_string(floor(r_rhoMin_fix.Z() / 100.0)), 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2 * sinTheta_customV, "PID_dE_E");
} }
plotter->Fill1D("pczfix_A1C2_1d_qqq", 600, -200, 200, pcz_fix); plotter->Fill1D("pczfix_A1C2_1d_qqq", 600, -200, 200, pcz_fix, "Z_Reconstruction");
plotter->Fill2D("pczfix_vs_qqqpczguess_A1C2", 600, -200, 200, 600, -200, 200, pcz_guess_int, pcz_fix); plotter->Fill2D("pczfix_vs_qqqpczguess_A1C2", 600, -200, 200, 600, -200, 200, pcz_guess_int, pcz_fix, "Z_Reconstruction");
plotter->Fill2D("pczguess_vs_pc_int_A1C2", 400, -200, 200, 600, -400, 400, pcz_guess_int, pcevent.pos.Z(), "phicut"); plotter->Fill2D("pczguess_vs_pc_int_A1C2", 400, -200, 200, 600, -400, 400, pcz_guess_int, pcevent.pos.Z(), "Z_Reconstruction");
plotter->Fill2D("pczfix_residual_vs_pczguess_A1C2", 600, -200, 200, 200, -100, 100, pcz_guess_37, pcz_fix - pcz_guess_37, "Residuals");
plotter->Fill2D("pczfix_residual_vs_phi_A1C2", 200, 0, 6.28, 200, -100, 100, r_rhoMin_fix.Phi(), pcz_fix - pcz_guess_37, "Residuals");
plotter->Fill1D("pczfix-qqqpczguess_A1C2", 200, -100, 100, pcz_fix - pcz_guess_37, "Residuals");
plotter->Fill1D("pczfix-qqqpczint_A1C2", 200, -100, 100, pcz_fix - pcz_guess_int, "Residuals");
plotter->Fill1D("pczguess_vs_int_residualsqqq", 200, -50, 50, pcz_guess_int - pcz_guess_37, "Residuals");
double path_length = (qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Mag() * 0.1; double path_length = (qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Mag() * 0.1;
// std::cout << path_length << std::endl;
double qqqEfix = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy1) - path_length); double qqqEfix = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy1) - path_length);
double qqqEfix_p = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(qqqevent.Energy1) - path_length); double qqqEfix_p = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(qqqevent.Energy1) - path_length);
plotter->Fill2D("E_thetaf_AnodeQQQR_TC1_PC" + std::to_string(phicut), 180, 0, 180, 600, 0, 15, (qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Theta() * 180 / M_PI, qqqevent.Energy1); plotter->Fill2D("E_thetaf_AnodeQQQR_TC1_PC" + std::to_string(phicut), 180, 0, 180, 600, 0, 15, (qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Theta() * 180 / M_PI, qqqevent.Energy1, "Kinematics_Angles");
if (lowercut_cath) if (lowercut_cath)
plotter->Fill2D("Ef_thetaf_AnodeQQQR_TC1_PC" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath), 180, 0, 180, 600, 0, 15, (qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Theta() * 180 / M_PI, qqqEfix_p); plotter->Fill2D("Ef_thetaf_AnodeQQQR_TC1_PC" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath), 180, 0, 180, 600, 0, 15, (qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Theta() * 180 / M_PI, qqqEfix_p, "Kinematics_Angles");
else else
{ {
std::string zcut = "_" + std::to_string((TMath::Abs(r_rhoMin_fix.Z()) < 180)); std::string zcut = "_" + std::to_string((TMath::Abs(r_rhoMin_fix.Z()) < 180));
plotter->Fill2D("Ef_thetaf_AnodeQQQR_TC1_PC" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath) + zcut, 180, 0, 180, 600, 0, 15, (qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Theta() * 180 / M_PI, qqqEfix); plotter->Fill2D("Ef_thetaf_AnodeQQQR_TC1_PC" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath) + zcut, 180, 0, 180, 600, 0, 15, (qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Theta() * 180 / M_PI, qqqEfix, "Kinematics_Angles");
} }
std::string morecuts = "_pidlow" + std::to_string(lowercut_cath) + "_vertexfix=" + std::to_string(floor(r_rhoMin_fix.Z() / 20) * 20 + 10); plotter->Fill2D("dE3_Ef_AnodeQQQR_TC1" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath), 600, 0, 15, 800, 0, 40000, qqqEfix, pcevent.Energy1 * sinTheta_customV, "PID_dE_E");
// plotter->Fill2D("E_thetaf_AnodeQQQR_TC1_PC"+std::to_string(phicut)+morecuts,180,0,180,800,0,8,(qqqevent.pos - TVector3(0,0,r_rhoMin_fix.Z())).Theta()*180/M_PI,qqqevent.Energy1,"morecuts"); plotter->Fill2D("dE3_Ef_CathodeQQQR_TC1PC" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath), 600, 0, 15, 800, 0, 10000, qqqEfix, pcevent.Energy2 * sinTheta_customV, "PID_dE_E");
// plotter->Fill2D("Ef_thetaf_AnodeQQQR_TC1_PC"+std::to_string(phicut)+morecuts,180,0,180,800,0,8,(qqqevent.pos - TVector3(0,0,r_rhoMin_fix.Z())).Theta()*180/M_PI,qqqEfix,"morecuts");
plotter->Fill2D("dE3_Ef_AnodeQQQR_TC1" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath), 600, 0, 15, 800, 0, 40000, qqqEfix, pcevent.Energy1 * sinTheta_customV);
plotter->Fill2D("dE3_Ef_CathodeQQQR_TC1PC" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath), 600, 0, 15, 800, 0, 10000, qqqEfix, pcevent.Energy2 * sinTheta_customV);
} }
if(BenchMark)
{
}
double qqqrho = qqqevent.pos.Perp(); double qqqrho = qqqevent.pos.Perp();
double qqqz = (qqqevent.pos - TVector3(0, 0, source_vertex)).Z(); double qqqz = (qqqevent.pos - TVector3(0, 0, source_vertex)).Z();
double tan_theta = qqqrho / qqqz; double tan_theta = qqqrho / qqqz;
double pcz_guess_int2 = z_to_crossover_rho(pcevent.pos.Z()) / tan_theta + source_vertex; double pcz_guess_int2 = z_to_crossover_rho(pcevent.pos.Z()) / tan_theta + source_vertex;
plotter->Fill2D("pczguess_vs_pc_int2", 180, 0, 200, 150, 0, 200, pcz_guess_int2, pcevent.pos.Z(), "phicut"); plotter->Fill2D("pczguess_vs_pc_int2", 180, 0, 200, 150, 0, 200, pcz_guess_int2, pcevent.pos.Z(), "Z_Reconstruction");
double qqqz2 = (qqqevent.pos - r_rhoMin).Z(); double qqqz2 = (qqqevent.pos - r_rhoMin).Z();
double tan_theta2 = qqqrho / qqqz2; double tan_theta2 = qqqrho / qqqz2;
double pcz_guess_int3 = z_to_crossover_rho(pcevent.pos.Z()) / tan_theta2 + r_rhoMin.Z(); double pcz_guess_int3 = z_to_crossover_rho(pcevent.pos.Z()) / tan_theta2 + r_rhoMin.Z();
plotter->Fill2D("pczguess_vs_pc_int3", 180, 0, 200, 150, 0, 200, pcz_guess_int3, pcevent.pos.Z(), "phicut"); plotter->Fill2D("pczguess_vs_pc_int3", 180, 0, 200, 150, 0, 200, pcz_guess_int3, pcevent.pos.Z(), "Z_Reconstruction");
// plotter->Fill2D("pczguess_vs_pc_int2_a"+std::to_string(pcevent.multi1)+"_c"+std::to_string(pcevent.multi2),180,0,200,150,0,200,pcz_guess_int2,pcevent.pos.Z(),"phicut");
double pcz_guess = pcz_guess_int; double pcz_guess = pcz_guess_int;
plotter->Fill2D("pctheta_vs_qqqtheta_sv", 180, -360, 360, 180, -360, 360, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, (pcevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, "phicut"); plotter->Fill2D("pctheta_vs_qqqtheta_sv", 180, -360, 360, 180, -360, 360, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, (pcevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, "Kinematics_Angles");
plotter->Fill2D("pctheta_vs_qqqtheta_rmz", 180, -360, 360, 180, -360, 360, (qqqevent.pos - TVector3(0, 0, r_rhoMin.Z())).Theta() * 180 / M_PI, (pcevent.pos - TVector3(0, 0, r_rhoMin.Z())).Theta() * 180 / M_PI, "phicut"); plotter->Fill2D("pctheta_vs_qqqtheta_rmz", 180, -360, 360, 180, -360, 360, (qqqevent.pos - TVector3(0, 0, r_rhoMin.Z())).Theta() * 180 / M_PI, (pcevent.pos - TVector3(0, 0, r_rhoMin.Z())).Theta() * 180 / M_PI, "Kinematics_Angles");
plotter->Fill2D("pctheta_vs_qqqtheta_rm", 180, -360, 360, 180, -360, 360, (qqqevent.pos - r_rhoMin).Theta() * 180 / M_PI, (pcevent.pos - r_rhoMin).Theta() * 180 / M_PI, "phicut"); plotter->Fill2D("pctheta_vs_qqqtheta_rm", 180, -360, 360, 180, -360, 360, (qqqevent.pos - r_rhoMin).Theta() * 180 / M_PI, (pcevent.pos - r_rhoMin).Theta() * 180 / M_PI, "Kinematics_Angles");
plotter->Fill2D("pczguess_vs_pc_phi=" + std::to_string(qqqevent.pos.Phi() * 180. / M_PI), 300, 0, 200, 150, 0, 200, pcz_guess, pcevent.pos.Z(), "phicut"); plotter->Fill2D("pczguess_vs_pc_phi=" + std::to_string(qqqevent.pos.Phi() * 180. / M_PI), 300, 0, 200, 150, 0, 200, pcz_guess, pcevent.pos.Z(), "Z_Reconstruction");
} }
} }
} // end PC QQQ coincidence }
// HALFTIME! Can stop here in future versions
// return kTRUE;
} }
// We pass plotter as a pointer, strings as const references (for speed), // We pass plotter as a pointer, strings as const references (for speed),
// and the Kinematics object as a reference. // and the Kinematics object as a reference.
void FillExcitationHistograms(HistPlotter *plotter, void FillExcitationHistograms(HistPlotter *plotter,

8
run_tr.sh
View File

@ -42,20 +42,20 @@ if [[ 1 -eq 0 ]]; then
fi fi
# --- Block 2: 27Al Alpha+Gas Runs (9, 12) --- # --- Block 2: 27Al Alpha+Gas Runs (9, 12) ---
if [[ 1 -eq 0 ]]; then if [[ 1 -eq 1 ]]; then
export DATASET="27Al" export DATASET="27Al"
export PREFIX="Run_" export PREFIX="Run_"
export OUT_DIR="Output_a" export OUT_DIR="Output_a"
echo "Processing 27Al alpha+gas runs..." echo "Processing 27Al alpha+gas runs..."
export source_vertex=-5.36; export timecut_low=12.0; export timecut_high=119.0; process_run 9 "$slope" export source_vertex=-5.36; export timecut_low=12.0; export timecut_high=119.0; process_run 9 "$slope"
unset timecut_high
export source_vertex=53.44; export timecut_low=400.0; process_run 12 "$slope" export source_vertex=53.44; export timecut_low=400.0; process_run 12 "$slope"
unset timecut_low unset timecut_low
unset timecut_high
fi fi
# --- Block 3: 27Al Protons+Gas Runs (15, 17-22) --- # --- Block 3: 27Al Protons+Gas Runs (15, 17-22) ---
if [[ 1 -eq 0 ]]; then if [[ 1 -eq 1 ]]; then
export DATASET="27Al" export DATASET="27Al"
export PREFIX="Run_" export PREFIX="Run_"
export OUT_DIR="Output_p" export OUT_DIR="Output_p"
@ -77,7 +77,7 @@ if [[ 1 -eq 0 ]]; then
fi fi
# --- Block 5: 17F Alpha Run with Gas (18-21) --- # --- Block 5: 17F Alpha Run with Gas (18-21) ---
if [[ 1 -eq 1 ]]; then if [[ 1 -eq 0 ]]; then
export DATASET="17F" export DATASET="17F"
export PREFIX="SourceRun_" export PREFIX="SourceRun_"
export OUT_DIR="Output_a" export OUT_DIR="Output_a"

6
scan_slope_runs.sh
View File

@ -20,8 +20,8 @@ process_run() {
export -f process_run export -f process_run
export SCAN_DIR export SCAN_DIR
for slope_x100 in 40 42 44 46 48 50 52 54 56 58 60 62 64 66; do for slope_x100 in 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 72 74 76 78 80 82 84 86 88 90; do
# for slope_x100 in 58 60 62 64 66; do # for slope_x100 in 68 70 72 74 76 78 80 82 84 86 88 90; do
slope=$(awk "BEGIN{printf \"%.2f\", $slope_x100/100}") slope=$(awk "BEGIN{printf \"%.2f\", $slope_x100/100}")
echo "=== slope=${slope} ===" echo "=== slope=${slope} ==="
@ -33,7 +33,7 @@ for slope_x100 in 40 42 44 46 48 50 52 54 56 58 60 62 64 66; do
# 27Al alpha+gas runs (9, 12, 13) # 27Al alpha+gas runs (9, 12, 13)
export DATASET="27Al" PREFIX="Run_" export DATASET="27Al" PREFIX="Run_"
echo " 27Al runs 9 12 13..." echo " 27Al runs 9 12 13..."
# export source_vertex=53.44; export timecut_low=400.0; process_run 12 "$slope" export source_vertex=53.44; export timecut_low=400.0; process_run 12 "$slope"
export source_vertex=-5.36; export timecut_low=12.0; export timecut_high=120.0; process_run 9 "$slope" export source_vertex=-5.36; export timecut_low=12.0; export timecut_high=120.0; process_run 9 "$slope"
unset timecut_low unset timecut_low
unset timecut_high unset timecut_high

631
scratch/plot_slope_scan.C
View File

@ -8,339 +8,516 @@
#include <TLegend.h> #include <TLegend.h>
#include <TLine.h> #include <TLine.h>
#include <TString.h> #include <TString.h>
#include <TStyle.h>
#include <TApplication.h>
void plot_slope_scan() void plot_slope_scan()
{ {
// ========================================================================= // =========================================================================
// --- Configuration --- // --- Configuration ---
// ========================================================================= // =========================================================================
// Set to true to inspect each Gaussian fit. Double-click the canvas to advance.
bool checkFits = false; bool checkFits = false;
gStyle->SetOptStat(0); gStyle->SetOptStat(0);
double slopes[] = {0.40, 0.42, 0.44, 0.46, 0.48, 0.50, 0.52, 0.54, 0.56, 0.58, 0.60, 0.62, 0.64, 0.66};
const int N = 14; // Updated to 14 to match the size of your slopes array
// int runs[] = {9, 12, 18, 19, 20, 21}; // double slopes[] = {
// const int NRUNS = 6; // 0.42, 0.44, 0.46, 0.48, 0.50, 0.52, 0.54, 0.56, 0.58, 0.60,
// 0.62, 0.64, 0.66, 0.68, 0.70, 0.72, 0.74, 0.76, 0.78, 0.80,
// 0.82, 0.84, 0.86, 0.88, 0.90};
// const int N = 25;
double slopes[] = {0.42, 0.44, 0.46, 0.48, 0.50, 0.52, 0.54, 0.56, 0.58, 0.60, 0.62, 0.64};
const int N = 12; // Updated to 14 to match the size of your slopes array
int runs[] = {9, 12, 18, 19, 20, 21};
const int NRUNS = 6;
// int runs[] = {9, 12}; // int runs[] = {9, 12};
// const int NRUNS = 2; // const int NRUNS = 2;
int runs[] = {18, 19, 20, 21}; // int runs[] = {18, 19, 20, 21};
const int NRUNS = 4; // const int NRUNS = 4;
Int_t runColors[6] = {kRed + 1, kBlue + 2, kGreen + 2, kOrange + 1, kViolet + 2, kCyan + 2}; Int_t runColors[6] = {kRed + 1, kBlue + 2, kGreen + 2, kOrange + 1, kViolet + 2, kCyan + 2};
// --- Initialize Graphs --- // --- Initialize Graphs for both SX3 and QQQ ---
TGraph *gStdDevRun[NRUNS]; TGraph *gStdDevRun_sx3[NRUNS], *gMeanRun_sx3[NRUNS];
TGraph *gMeanRun[NRUNS]; TGraph *gStdDevRun_qqq[NRUNS], *gMeanRun_qqq[NRUNS];
for (int ir = 0; ir < NRUNS; ir++) for (int ir = 0; ir < NRUNS; ir++)
{ {
gStdDevRun[ir] = new TGraph(N); gStdDevRun_sx3[ir] = new TGraph(N);
gMeanRun[ir] = new TGraph(N); gMeanRun_sx3[ir] = new TGraph(N);
gStdDevRun_qqq[ir] = new TGraph(N);
gMeanRun_qqq[ir] = new TGraph(N);
} }
TGraph *gStdDevSum = new TGraph(N);
TGraph *gMeanSum = new TGraph(N);
TH1 *h1[NRUNS][N]; TGraph *gStdDevSum_sx3 = new TGraph(N), *gMeanSum_sx3 = new TGraph(N);
TGraph *gStdDevSum_qqq = new TGraph(N), *gMeanSum_qqq = new TGraph(N);
// Arrays to hold histograms
TH1 *h1_sx3[NRUNS][N];
TH1 *h1_qqq[NRUNS][N];
TH1 *h_vtx_sx3[NRUNS][N];
TH1 *h_vtx_qqq[NRUNS][N];
// --- Load Data --- // --- Load Data ---
for (int ir = 0; ir < NRUNS; ir++) for (int ir = 0; ir < NRUNS; ir++)
{ {
for (int i = 0; i < N; i++) for (int i = 0; i < N; i++)
{ {
h1[ir][i] = nullptr; h1_sx3[ir][i] = nullptr;
h1_qqq[ir][i] = nullptr;
h_vtx_sx3[ir][i] = nullptr;
h_vtx_qqq[ir][i] = nullptr;
TString path = TString::Format("slope_scan/run%03d/slope_%.2f.root", runs[ir], slopes[i]); TString path = TString::Format("slope_scan/run%03d/slope_%.2f.root", runs[ir], slopes[i]);
TFile *f = TFile::Open(path); TFile *f = TFile::Open(path);
if (!f || f->IsZombie()) if (!f || f->IsZombie())
continue; continue;
TH1 *t1 = (TH1 *)f->Get("pczfix-sx3pczguess_A1C2"); TH1 *t1_s = (TH1 *)f->Get("pczfix-sx3pczguess_A1C2");
if (t1) if (t1_s)
{ {
t1->SetDirectory(nullptr); t1_s->SetDirectory(nullptr);
h1[ir][i] = t1; h1_sx3[ir][i] = t1_s;
} }
TH1 *t1_q = (TH1 *)f->Get("pczfix-qqqpczguess_A1C2");
if (t1_q)
{
t1_q->SetDirectory(nullptr);
h1_qqq[ir][i] = t1_q;
}
TH1 *t_sx3 = (TH1 *)f->Get("VertexRecon_pczfix_sx3");
if (t_sx3)
{
t_sx3->SetDirectory(nullptr);
h_vtx_sx3[ir][i] = t_sx3;
}
TH1 *t_qqq = (TH1 *)f->Get("VertexRecon_pczfix_qqq");
if (t_qqq)
{
t_qqq->SetDirectory(nullptr);
h_vtx_qqq[ir][i] = t_qqq;
}
f->Close(); f->Close();
} }
} }
// ========================================================================= // =========================================================================
// --- Calculate Mean & StdDev with Bounded Gaussian Fit --- // --- Helper Lambda for Gaussian Fitting ---
// ========================================================================= // =========================================================================
std::vector<int> nPtsRun(NRUNS, 0); auto fitPeak = [](TH1 *h, double fitWindow, bool debug, TCanvas *c, double &outMean, double &outStdDev)
int nPtsSum = 0;
double fitWindow = 15.0; // Window around the peak to fit (mm)
TCanvas *cDebug = nullptr;
if (checkFits)
{ {
cDebug = new TCanvas("cDebug", "Fit Diagnostic Debugger", 0, 0, 800, 600); if (!h || h->GetEntries() <= 0)
} return false;
int maxBin = h->GetMaximumBin();
double peakX = h->GetXaxis()->GetBinCenter(maxBin);
TF1 *gfit = new TF1("gfit", "gaus", peakX - fitWindow, peakX + fitWindow);
TString fitOpt = debug ? "RQS" : "RQS0";
if (debug && c)
c->cd();
TFitResultPtr fitRes = h->Fit(gfit, fitOpt);
if (debug && c)
{
h->GetXaxis()->SetRangeUser(peakX - 50, peakX + 50);
h->SetMaximum(h->GetMaximum() * 1.2);
h->Draw();
c->Modified();
c->Update();
while (c->WaitPrimitive())
;
}
if (fitRes >= 0)
{
outMean = gfit->GetParameter(1);
outStdDev = gfit->GetParameter(2);
}
else
{
outMean = h->GetMean();
outStdDev = h->GetStdDev();
}
delete gfit;
return true;
};
// =========================================================================
// --- Calculate Mean & StdDev for SX3 & QQQ ---
// =========================================================================
std::vector<int> nPtsRun_sx3(NRUNS, 0), nPtsRun_qqq(NRUNS, 0);
int nPtsSum_sx3 = 0, nPtsSum_qqq = 0;
double fitWin = 15.0;
TCanvas *cDebug = checkFits ? new TCanvas("cDebug", "Debugger", 0, 0, 800, 600) : nullptr;
for (int i = 0; i < N; i++) for (int i = 0; i < N; i++)
{ {
TH1 *hsum = nullptr; TH1 *hsum_sx3 = nullptr, *hsum_qqq = nullptr;
for (int ir = 0; ir < NRUNS; ir++) for (int ir = 0; ir < NRUNS; ir++)
{ {
if (!h1[ir][i] || h1[ir][i]->GetEntries() <= 0)
continue;
// 1. Find center of the highest peak
int maxBin = h1[ir][i]->GetMaximumBin();
double peakX = h1[ir][i]->GetXaxis()->GetBinCenter(maxBin);
// 2. Define the fit
TF1 *gfit = new TF1("gfit", "gaus", peakX - fitWindow, peakX + fitWindow);
// 3. Fit Conditionally (Draw if checking, otherwise run quietly in background)
TString fitOpt = checkFits ? "RQS" : "RQS0";
if (checkFits)
cDebug->cd(); // Ensure we draw on the debug canvas
TFitResultPtr fitRes = h1[ir][i]->Fit(gfit, fitOpt);
// --- DIAGNOSTIC PAUSE ---
if (checkFits)
{
h1[ir][i]->SetTitle(TString::Format("Run %d | Slope %.2f", runs[ir], slopes[i]));
h1[ir][i]->GetXaxis()->SetRangeUser(peakX - 50, peakX + 50); // Zoom in X
// Scale Y-axis to the maximum of this specific zoomed window + 20% headroom
h1[ir][i]->SetMaximum(h1[ir][i]->GetMaximum() * 1.2);
h1[ir][i]->Draw();
cDebug->Modified();
cDebug->Update();
printf("Visual Check: Run %d | Slope %.2f. Double-click canvas to continue...\n", runs[ir], slopes[i]);
while (cDebug->WaitPrimitive())
;
}
double mean, stddev; double mean, stddev;
if (fitRes >= 0)
{ // If fit succeeds
mean = gfit->GetParameter(1);
stddev = gfit->GetParameter(2);
}
else
{ // Fallback to raw stats if fit fails
mean = h1[ir][i]->GetMean();
stddev = h1[ir][i]->GetStdDev();
}
gStdDevRun[ir]->SetPoint(nPtsRun[ir], slopes[i], stddev); // Fit SX3
gMeanRun[ir]->SetPoint(nPtsRun[ir]++, slopes[i], mean); if (fitPeak(h1_sx3[ir][i], fitWin, checkFits, cDebug, mean, stddev))
delete gfit;
// Build the cumulative sum histogram
if (!hsum)
hsum = (TH1 *)h1[ir][i]->Clone(TString::Format("hsum_%.2f", slopes[i]));
else
hsum->Add(h1[ir][i]);
}
// --- Process Summed Histogram ---
if (hsum && hsum->GetEntries() > 0)
{
int maxBinSum = hsum->GetMaximumBin();
double peakXSum = hsum->GetXaxis()->GetBinCenter(maxBinSum);
TF1 *gfitSum = new TF1("gfitSum", "gaus", peakXSum - fitWindow, peakXSum + fitWindow);
TString fitOptSum = checkFits ? "RQS" : "RQS0";
if (checkFits)
cDebug->cd();
TFitResultPtr fitResSum = hsum->Fit(gfitSum, fitOptSum);
// --- DIAGNOSTIC PAUSE FOR SUM ---
if (checkFits)
{ {
hsum->SetTitle(TString::Format("SUMMED RUNS | Slope %.2f", slopes[i])); gStdDevRun_sx3[ir]->SetPoint(nPtsRun_sx3[ir], slopes[i], stddev);
hsum->GetXaxis()->SetRangeUser(peakXSum - 50, peakXSum + 50); gMeanRun_sx3[ir]->SetPoint(nPtsRun_sx3[ir]++, slopes[i], mean);
if (!hsum_sx3)
// Scale Y-axis for the sum + 20% headroom hsum_sx3 = (TH1 *)h1_sx3[ir][i]->Clone(TString::Format("hsum_sx3_%.2f", slopes[i]));
hsum->SetMaximum(hsum->GetMaximum() * 1.2); else
hsum_sx3->Add(h1_sx3[ir][i]);
hsum->Draw(); }
cDebug->Modified(); if (runs[ir] == 20 || runs[ir] == 21)
cDebug->Update(); continue;
printf("Visual Check: SUMMED | Slope %.2f. Double-click canvas to continue...\n", slopes[i]); // Fit QQQ
while (cDebug->WaitPrimitive()) if (fitPeak(h1_qqq[ir][i], fitWin, checkFits, cDebug, mean, stddev))
; {
gStdDevRun_qqq[ir]->SetPoint(nPtsRun_qqq[ir], slopes[i], stddev);
gMeanRun_qqq[ir]->SetPoint(nPtsRun_qqq[ir]++, slopes[i], mean);
if (!hsum_qqq)
hsum_qqq = (TH1 *)h1_qqq[ir][i]->Clone(TString::Format("hsum_qqq_%.2f", slopes[i]));
else
hsum_qqq->Add(h1_qqq[ir][i]);
} }
// if (fitResSum >= 0)
// {
// gStdDevSum->SetPoint(nPtsSum, slopes[i], gfitSum->GetParameter(2));
// gMeanSum->SetPoint(nPtsSum++, slopes[i], gfitSum->GetParameter(1));
// }
// else
// {
// gStdDevSum->SetPoint(nPtsSum, slopes[i], hsum->GetStdDev());
// gMeanSum->SetPoint(nPtsSum++, slopes[i], hsum->GetMean());
// }
// delete gfitSum;
// delete hsum;
} }
// Fit Sums
double sMean, sStdDev;
if (fitPeak(hsum_sx3, fitWin, checkFits, cDebug, sMean, sStdDev))
{
gStdDevSum_sx3->SetPoint(nPtsSum_sx3, slopes[i], sStdDev);
gMeanSum_sx3->SetPoint(nPtsSum_sx3++, slopes[i], sMean);
}
if (fitPeak(hsum_qqq, fitWin, checkFits, cDebug, sMean, sStdDev))
{
gStdDevSum_qqq->SetPoint(nPtsSum_qqq, slopes[i], sStdDev);
gMeanSum_qqq->SetPoint(nPtsSum_qqq++, slopes[i], sMean);
}
delete hsum_sx3;
delete hsum_qqq;
} }
if (checkFits && cDebug) if (cDebug)
{ delete cDebug;
delete cDebug; // Clean up debug canvas before plotting final results
}
// ========================================================================= // =========================================================================
// ---- Canvas 1: StdDev and Mean vs Slope // ---- Canvas 1: StdDev and Mean vs Slope (Both SX3 and QQQ) ----
// ========================================================================= // =========================================================================
TCanvas *c1 = new TCanvas("c_stats", "1D Residual Stats vs Slope", 0, 0, 1600, 600); TCanvas *c1 = new TCanvas("c_stats", "1D Residual Stats vs Slope", 0, 0, 1800, 800);
c1->Divide(2, 1);
c1->cd(1); // 1. Create a Master Legend in the top 8% of the canvas
gPad->SetGrid(1, 1); c1->cd();
TLegend *leg1 = new TLegend(0.78, 0.45, 0.97, 0.88); TLegend *leg1_master = new TLegend(0.1, 0.92, 0.9, 1.0);
leg1->SetBorderSize(1); leg1_master->SetNColumns(6);
TMultiGraph *mg1 = new TMultiGraph("mg_stddev", "StdDev(pczfix - sx3pczguess) vs Slope;Slope;StdDev (mm)"); leg1_master->SetBorderSize(0);
leg1_master->SetFillStyle(0);
leg1_master->SetTextSize(0.03);
for (int ir = 0; ir < NRUNS; ir++) for (int ir = 0; ir < NRUNS; ir++)
{ {
if (nPtsRun[ir] == 0) TLine *dSX3 = new TLine();
continue; dSX3->SetLineColor(runColors[ir]);
gStdDevRun[ir]->SetLineColor(runColors[ir]); dSX3->SetLineWidth(2);
gStdDevRun[ir]->SetMarkerColor(runColors[ir]); dSX3->SetLineStyle(kSolid);
gStdDevRun[ir]->SetLineWidth(2); leg1_master->AddEntry(dSX3, TString::Format("Run %d (SX3)", runs[ir]), "l");
gStdDevRun[ir]->SetMarkerStyle(20 + ir);
mg1->Add(gStdDevRun[ir], "PL");
leg1->AddEntry(gStdDevRun[ir], TString::Format("run %d", runs[ir]), "lp");
} }
if (nPtsSum > 0) for (int ir = 0; ir < NRUNS; ir++)
{ {
gStdDevSum->SetLineColor(kBlack); TLine *dQQQ = new TLine();
gStdDevSum->SetLineWidth(3); dQQQ->SetLineColor(runColors[ir]);
gStdDevSum->SetLineStyle(kDashed); dQQQ->SetLineWidth(2);
gStdDevSum->SetMarkerStyle(29); dQQQ->SetLineStyle(kDashed);
mg1->Add(gStdDevSum, "PL"); leg1_master->AddEntry(dQQQ, TString::Format("Run %d (QQQ)", runs[ir]), "l");
leg1->AddEntry(gStdDevSum, "sum", "lp"); }
leg1_master->Draw();
// 2. Create the Grid Pad for Canvas 1
TPad *pad1_grid = new TPad("pad1_grid", "Grid", 0.0, 0.0, 1.0, 0.92);
pad1_grid->Draw();
pad1_grid->Divide(2, 1);
// --- Pad 1: StdDev ---
pad1_grid->cd(1);
gPad->SetGrid(1, 1);
TMultiGraph *mg1 = new TMultiGraph("mg_stddev", "StdDev of Residuals vs Slope;Slope;StdDev (mm)");
for (int ir = 0; ir < NRUNS; ir++)
{
if (nPtsRun_sx3[ir] > 0)
{
gStdDevRun_sx3[ir]->SetLineColor(runColors[ir]);
gStdDevRun_sx3[ir]->SetLineWidth(2);
gStdDevRun_sx3[ir]->SetMarkerStyle(20);
gStdDevRun_sx3[ir]->SetMarkerColor(runColors[ir]);
mg1->Add(gStdDevRun_sx3[ir], "PL");
}
if (runs[ir] == 20 || runs[ir] == 21)
continue;
if (nPtsRun_qqq[ir] > 0)
{
gStdDevRun_qqq[ir]->SetLineColor(runColors[ir]);
gStdDevRun_qqq[ir]->SetLineWidth(2);
gStdDevRun_qqq[ir]->SetLineStyle(kDashed);
gStdDevRun_qqq[ir]->SetMarkerStyle(24);
gStdDevRun_qqq[ir]->SetMarkerColor(runColors[ir]);
mg1->Add(gStdDevRun_qqq[ir], "PL");
}
} }
mg1->Draw("A"); mg1->Draw("A");
leg1->Draw();
double bestSlope = 0, bestStdDev = 1e9; double bestSlope = 0, bestStdDev = 1e9;
for (int i = 0; i < gStdDevRun[0]->GetN(); i++) if (gStdDevRun_sx3[0]->GetN() > 0)
{ {
double x, y; for (int i = 0; i < gStdDevRun_sx3[0]->GetN(); i++)
gStdDevRun[0]->GetPoint(i, x, y);
if (y > 0 && y < bestStdDev)
{ {
bestStdDev = y; double x, y;
bestSlope = x; gStdDevRun_sx3[0]->GetPoint(i, x, y);
if (y > 0 && y < bestStdDev)
{
bestStdDev = y;
bestSlope = x;
}
} }
TLine *lb = new TLine(bestSlope, mg1->GetHistogram()->GetMinimum(), bestSlope, mg1->GetHistogram()->GetMaximum());
lb->SetLineColor(kRed);
lb->SetLineStyle(kDashed);
lb->Draw();
printf(">>> Best slope (Run %d) = %.2f StdDev=%.2f mm\n", runs[0], bestSlope, bestStdDev);
} }
TLine *lb = new TLine(bestSlope, mg1->GetHistogram()->GetMinimum(), // --- Pad 2: Mean ---
bestSlope, mg1->GetHistogram()->GetMaximum()); pad1_grid->cd(2);
lb->SetLineColor(kRed);
lb->SetLineStyle(kDashed);
lb->Draw();
printf(">>> Best slope (Run %d) = %.2f StdDev=%.2f mm\n", runs[0], bestSlope, bestStdDev);
c1->cd(2);
gPad->SetGrid(1, 1); gPad->SetGrid(1, 1);
TLegend *leg2 = new TLegend(0.78, 0.45, 0.97, 0.88); TMultiGraph *mg2 = new TMultiGraph("mg_mean", "Mean of Residuals vs Slope;Slope;Mean Offset (mm)");
leg2->SetBorderSize(1);
TMultiGraph *mg2 = new TMultiGraph("mg_mean", "Mean(pczfix - sx3pczguess) vs Slope;Slope;Mean Offset (mm)");
for (int ir = 0; ir < NRUNS; ir++) for (int ir = 0; ir < NRUNS; ir++)
{ {
if (nPtsRun[ir] == 0) if (nPtsRun_sx3[ir] > 0)
{
gMeanRun_sx3[ir]->SetLineColor(runColors[ir]);
gMeanRun_sx3[ir]->SetLineWidth(2);
gMeanRun_sx3[ir]->SetMarkerStyle(20);
gMeanRun_sx3[ir]->SetMarkerColor(runColors[ir]);
mg2->Add(gMeanRun_sx3[ir], "PL");
}
if (runs[ir] == 20 || runs[ir] == 21)
continue; continue;
gMeanRun[ir]->SetLineColor(runColors[ir]); if (nPtsRun_qqq[ir] > 0)
gMeanRun[ir]->SetMarkerColor(runColors[ir]); {
gMeanRun[ir]->SetLineWidth(2); gMeanRun_qqq[ir]->SetLineColor(runColors[ir]);
gMeanRun[ir]->SetMarkerStyle(20 + ir); gMeanRun_qqq[ir]->SetLineWidth(2);
mg2->Add(gMeanRun[ir], "PL"); gMeanRun_qqq[ir]->SetLineStyle(kDashed);
leg2->AddEntry(gMeanRun[ir], TString::Format("run %d", runs[ir]), "lp"); gMeanRun_qqq[ir]->SetMarkerStyle(24);
} gMeanRun_qqq[ir]->SetMarkerColor(runColors[ir]);
if (nPtsSum > 0) mg2->Add(gMeanRun_qqq[ir], "PL");
{ }
gMeanSum->SetLineColor(kBlack);
gMeanSum->SetLineWidth(3);
gMeanSum->SetLineStyle(kDashed);
gMeanSum->SetMarkerStyle(29);
mg2->Add(gMeanSum, "PL");
leg2->AddEntry(gMeanSum, "sum", "lp");
} }
mg2->Draw("A"); mg2->Draw("A");
TLine *lz = new TLine(slopes[0], 0, slopes[N - 1], 0); TLine *lz = new TLine(slopes[0], 0, slopes[N - 1], 0);
lz->SetLineColor(kGray + 2); lz->SetLineColor(kGray + 2);
lz->SetLineStyle(kDotted); lz->SetLineStyle(kDotted);
lz->Draw(); lz->Draw();
leg2->Draw();
c1->SaveAs("slope_scan_1d_metric.png"); c1->SaveAs("slope_scan_1d_metric.png");
// ========================================================================= // =========================================================================
// ---- Canvas 2: Per-slope pads, overlaying the 1D Residual distributions // ---- Canvas 2: Per-slope pads, 1D Residuals (SX3 + QQQ) ----
// ========================================================================= // =========================================================================
TCanvas *c2 = new TCanvas("c_perslope_1d", "Per-slope 1D Residuals: all runs overlaid", 0, 100, 1800, 1200); TCanvas *c2 = new TCanvas("c_perslope_1d", "Per-slope 1D Residuals (SX3 & QQQ)", 0, 50, 2000, 1600);
c2->Divide(4, 4);
c2->cd();
TLegend *leg2_master = new TLegend(0.1, 0.92, 0.9, 1.0);
leg2_master->SetNColumns(6);
leg2_master->SetBorderSize(0);
leg2_master->SetFillStyle(0);
leg2_master->SetTextSize(0.025);
for (int ir = 0; ir < NRUNS; ir++)
{
TLine *dSX3 = new TLine();
dSX3->SetLineColor(runColors[ir]);
dSX3->SetLineWidth(2);
dSX3->SetLineStyle(kSolid);
leg2_master->AddEntry(dSX3, TString::Format("Run %d (SX3)", runs[ir]), "l");
}
for (int ir = 0; ir < NRUNS; ir++)
{
TLine *dQQQ = new TLine();
dQQQ->SetLineColor(runColors[ir]);
dQQQ->SetLineWidth(2);
dQQQ->SetLineStyle(kDashed);
leg2_master->AddEntry(dQQQ, TString::Format("Run %d (QQQ)", runs[ir]), "l");
}
leg2_master->Draw();
TPad *pad2_grid = new TPad("pad2_grid", "Grid", 0.0, 0.0, 1.0, 0.92);
pad2_grid->Draw();
pad2_grid->Divide(5, 5);
for (int i = 0; i < N; i++) for (int i = 0; i < N; i++)
{ {
c2->cd(i + 1); pad2_grid->cd(i + 1);
gPad->SetGrid(1, 1); gPad->SetGrid(1, 1);
TLegend *legS = new TLegend(0.62, 0.55, 0.88, 0.88);
legS->SetBorderSize(0);
// --- STEP A: Pre-scan to find the highest Y value across all runs for this slope ---
double maxY = 0; double maxY = 0;
for (int ir = 0; ir < NRUNS; ir++) for (int ir = 0; ir < NRUNS; ir++)
{ {
if (!h1[ir][i] || h1[ir][i]->GetEntries() <= 0) if (h1_sx3[ir][i] && h1_sx3[ir][i]->GetEntries() > 0)
continue;
// Set X range first so GetMaximum() doesn't accidentally grab a tail far away
h1[ir][i]->GetXaxis()->SetRangeUser(-50, 50);
if (h1[ir][i]->GetMaximum() > maxY)
{ {
maxY = h1[ir][i]->GetMaximum(); h1_sx3[ir][i]->GetXaxis()->SetRangeUser(-50, 50);
if (h1_sx3[ir][i]->GetMaximum() > maxY)
maxY = h1_sx3[ir][i]->GetMaximum();
}
if (runs[ir] == 20 || runs[ir] == 21)
continue;
if (h1_qqq[ir][i] && h1_qqq[ir][i]->GetEntries() > 0)
{
h1_qqq[ir][i]->GetXaxis()->SetRangeUser(-50, 50);
if (h1_qqq[ir][i]->GetMaximum() > maxY)
maxY = h1_qqq[ir][i]->GetMaximum();
} }
} }
// --- STEP B: Draw with the dynamically scaled Y-axis ---
bool first = true; bool first = true;
for (int ir = 0; ir < NRUNS; ir++) for (int ir = 0; ir < NRUNS; ir++)
{ {
if (!h1[ir][i] || h1[ir][i]->GetEntries() <= 0) if (h1_sx3[ir][i] && h1_sx3[ir][i]->GetEntries() > 0)
continue;
h1[ir][i]->SetLineColor(runColors[ir]);
h1[ir][i]->SetLineWidth(2);
h1[ir][i]->SetTitle(TString::Format("slope=%.2f;Residual (mm);Counts", slopes[i]));
if (first)
{ {
// Set the pad's Y-axis using the global max + 15% buffer so the legend fits cleanly h1_sx3[ir][i]->SetLineColor(runColors[ir]);
h1[ir][i]->SetMaximum(maxY * 1.15); h1_sx3[ir][i]->SetLineWidth(2);
h1[ir][i]->Draw("HIST"); h1_sx3[ir][i]->SetLineStyle(kSolid);
first = false; h1_sx3[ir][i]->SetTitle(TString::Format("slope=%.2f;Residual (mm);Counts", slopes[i]));
if (first)
{
h1_sx3[ir][i]->SetMaximum(maxY * 1.15);
h1_sx3[ir][i]->Draw("HIST");
first = false;
}
else
h1_sx3[ir][i]->Draw("HIST SAME");
} }
else if (h1_qqq[ir][i] && h1_qqq[ir][i]->GetEntries() > 0)
{ {
h1[ir][i]->Draw("HIST SAME"); h1_qqq[ir][i]->SetLineColor(runColors[ir]);
h1_qqq[ir][i]->SetLineWidth(2);
h1_qqq[ir][i]->SetLineStyle(kDashed);
if (first)
{
h1_qqq[ir][i]->SetTitle(TString::Format("slope=%.2f;Residual (mm);Counts", slopes[i]));
h1_qqq[ir][i]->SetMaximum(maxY * 1.15);
h1_qqq[ir][i]->Draw("HIST");
first = false;
}
else
h1_qqq[ir][i]->Draw("HIST SAME");
} }
legS->AddEntry(h1[ir][i], TString::Format("run %d", runs[ir]), "l");
} }
if (!first)
legS->Draw();
} }
c2->SaveAs("slope_scan_perslope_1d.png"); c2->SaveAs("slope_scan_perslope_1d.png");
// =========================================================================
// ---- Canvas 3: Vertex Recon (SX3 & QQQ Overlaid) ----
// =========================================================================
TCanvas *c3 = new TCanvas("c_vtx_recon", "Vertex Recon: SX3 vs QQQ", 0, 100, 2000, 1600);
c3->cd();
TLegend *leg3_master = new TLegend(0.1, 0.92, 0.9, 1.0);
leg3_master->SetNColumns(6);
leg3_master->SetBorderSize(0);
leg3_master->SetFillStyle(0);
leg3_master->SetTextSize(0.025);
for (int ir = 0; ir < NRUNS; ir++)
{
TLine *dSX3 = new TLine();
dSX3->SetLineColor(runColors[ir]);
dSX3->SetLineWidth(2);
dSX3->SetLineStyle(kSolid);
leg3_master->AddEntry(dSX3, TString::Format("Run %d (SX3)", runs[ir]), "l");
}
for (int ir = 0; ir < NRUNS; ir++)
{
TLine *dQQQ = new TLine();
dQQQ->SetLineColor(runColors[ir]);
dQQQ->SetLineWidth(2);
dQQQ->SetLineStyle(kDashed);
leg3_master->AddEntry(dQQQ, TString::Format("Run %d (QQQ)", runs[ir]), "l");
}
leg3_master->Draw();
TPad *pad3_grid = new TPad("pad3_grid", "Grid", 0.0, 0.0, 1.0, 0.92);
pad3_grid->Draw();
pad3_grid->Divide(5, 5);
for (int i = 0; i < N; i++)
{
pad3_grid->cd(i + 1);
gPad->SetGrid(1, 1);
double maxY3 = 0;
for (int ir = 0; ir < NRUNS; ir++)
{
if (h_vtx_sx3[ir][i] && h_vtx_sx3[ir][i]->GetEntries() > 0)
{
h_vtx_sx3[ir][i]->GetXaxis()->SetRangeUser(-200, 200);
if (h_vtx_sx3[ir][i]->GetMaximum() > maxY3)
maxY3 = h_vtx_sx3[ir][i]->GetMaximum();
}
if (h_vtx_qqq[ir][i] && h_vtx_qqq[ir][i]->GetEntries() > 0)
{
h_vtx_qqq[ir][i]->GetXaxis()->SetRangeUser(-200, 200);
if (h_vtx_qqq[ir][i]->GetMaximum() > maxY3)
maxY3 = h_vtx_qqq[ir][i]->GetMaximum();
}
}
bool first3 = true;
for (int ir = 0; ir < NRUNS; ir++)
{
if (h_vtx_sx3[ir][i] && h_vtx_sx3[ir][i]->GetEntries() > 0)
{
h_vtx_sx3[ir][i]->SetLineColor(runColors[ir]);
h_vtx_sx3[ir][i]->SetLineWidth(2);
h_vtx_sx3[ir][i]->SetLineStyle(kSolid);
h_vtx_sx3[ir][i]->SetTitle(TString::Format("slope=%.2f;Z_{Vertex} (mm);Counts", slopes[i]));
if (first3)
{
h_vtx_sx3[ir][i]->SetMaximum(maxY3 * 1.15);
h_vtx_sx3[ir][i]->Draw("HIST");
first3 = false;
}
else
h_vtx_sx3[ir][i]->Draw("HIST SAME");
}
if (h_vtx_qqq[ir][i] && h_vtx_qqq[ir][i]->GetEntries() > 0)
{
h_vtx_qqq[ir][i]->SetLineColor(runColors[ir]);
h_vtx_qqq[ir][i]->SetLineWidth(2);
h_vtx_qqq[ir][i]->SetLineStyle(kDashed);
if (first3)
{
h_vtx_qqq[ir][i]->SetTitle(TString::Format("slope=%.2f;Z_{Vertex} (mm);Counts", slopes[i]));
h_vtx_qqq[ir][i]->SetMaximum(maxY3 * 1.15);
h_vtx_qqq[ir][i]->Draw("HIST");
first3 = false;
}
else
h_vtx_qqq[ir][i]->Draw("HIST SAME");
}
}
}
c3->SaveAs("slope_scan_vtx_recon.png");
c1->Modified();
c1->Update();
c2->Modified(); c2->Modified();
c2->Update(); c2->Update();
while (gPad->WaitPrimitive()) c3->Modified();
; c3->Update();
} }