diff --git a/Armory/PC_StepLadder_Correction.h b/Armory/PC_StepLadder_Correction.h index 439fac3..79c3ea5 100644 --- a/Armory/PC_StepLadder_Correction.h +++ b/Armory/PC_StepLadder_Correction.h @@ -2,7 +2,7 @@ /*double model(double* x, double* p) { double result = x[0]; double factor = 29.0; - double slope = 0.7; + double slope = 0.40; 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]) < 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 result=y[0]; - double slope = 0.7; + double slope = 0.40; double factor = 0.0; 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; @@ -23,7 +23,7 @@ double model_invert(double *y, double *q) { double model_invert(double* y, double* p) { 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}; for(int i=0;i<7;i++) { 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 result = x[0]; double factor = 29.0; - double slope = 0.0; + double slope = 0.52; 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]) < 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 result=y[0]; -/* double slope = 1.0; +/* double slope = 0.40; double factor = 5.0; 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; diff --git a/TrackRecon.C b/TrackRecon.C index 6d3c278..913048c 100644 --- a/TrackRecon.C +++ b/TrackRecon.C @@ -45,7 +45,7 @@ bool doPCQQQClusterAnalysis = true; bool doOldAnalysis = false; bool do27AlapAnalysis = false; 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; const double anode_gain = 1.5146e-5; // channels --> MeV std::string dataset; @@ -597,7 +597,9 @@ Bool_t TrackRecon::Process(Long64_t entry) 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); //"?" // z used to be 75+30+23=128 // we found a 12mm shift towards the vertex later --> 116 @@ -1360,9 +1362,7 @@ void protonAlphaHistograms(HistPlotter *plotter, std::vector QQQ_Events, return; } -void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events) - -{ +void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events){ for (auto pcevent : PC_Events) { bool PCSX3TimeCut = false; @@ -1370,8 +1370,8 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s bool PCCSX3TimeCut = false; 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_pcC_sx3B" + std::to_string(sx3event.ch2), 640, -2000, 2000, sx3event.Time1 - pcevent.Time2, "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, "Timing"); if (sx3event.Time1 - pcevent.Time1 < 0) //-150 for alphas PCASX3TimeCut = 1; if (sx3event.Time1 - pcevent.Time2 < 0) //-200 for alphas @@ -1380,35 +1380,35 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s 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_pcC_sx3B", 640, -2000, 2000, sx3event.Time1 - pcevent.Time2); - plotter->Fill2D("dt_pcA_vs_sx3RE", 640, -2000, 2000, 400, 0, 30, sx3event.Time1 - pcevent.Time1, sx3event.Energy1); - plotter->Fill2D("dE_E_Anodesx3B", 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1); - plotter->Fill2D("dE_E_Cathodesx3B", 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2); + 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, "Timing"); + 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, "PID_dE_E"); + 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) - 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) - 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) - 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) - plotter->Fill2D("dE_E_Cathodesx3B_a2c1", 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2); - 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("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, "Kinematics_Angles"); if (PCSX3TimeCut) { - plotter->Fill1D("dt_pcA_sx3B_timecut", 640, -2000, 2000, sx3event.Time1 - pcevent.Time1); - plotter->Fill1D("dt_pcC_sx3B_timecut", 640, -2000, 2000, sx3event.Time1 - pcevent.Time2); - 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, pcevent.pos.X(), pcevent.pos.Y()); - 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->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, "Timing"); + 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(), "Vertex_Reconstruction"); + 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 sx3z = sx3event.pos.Z(); // w.r.t target origin at 90 for run12 + double sx3rho = 88.0; + double sx3z = sx3event.pos.Z(); double pcz = pcevent.pos.Z(); 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_Cathodesx3B", 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2 * 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), "PID_dE_E"); double sx3theta = TMath::ATan2(sx3rho, sx3z - source_vertex); double pczguess = 37.0 / TMath::Tan(sx3theta) + source_vertex; @@ -1419,78 +1419,79 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s 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()); 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", 600, -1300, 1300, vector_closest_to_z_sx3.Z()); - 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->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(), "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(), "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("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("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); // 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, "Timing"); + plotter->Fill2D("sx3phi_vs_time", 2000, 0, 2000, 180, -360, 360, pcevent.Time1 * 1e-9, sx3event.pos.Phi() * 180. / M_PI, "Timing"); - 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) { - plotter->Fill2D("pcz_vs_sx3pczguess_A1C2", 600, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z()); + plotter->Fill2D("pcz_vs_sx3pczguess_A1C2", 600, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z(), "Z_Reconstruction"); double pcz_fix = pcfix_func.Eval(pcevent.pos.Z()); TVector3 x2f(pcevent.pos.X(), pcevent.pos.Y(), pcz_fix); 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()); TVector3 r_rhoMin_fix = x1 + t_minimum * v; - 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->Fill2D("pczfix_vs_sx3pczguess_A1C2", 600, -200, 200, 600, -200, 200, pczguess, pcz_fix); - plotter->Fill2D("pczfix_residual_vs_pczguess_A1C2", 600, -200, 200, 200, -100, 100, pczguess, pcz_fix - pczguess); - plotter->Fill2D("pczfix_residual_vs_phi_A1C2", 200, 0, 6.28, 200, -100, 100, r_rhoMin_fix.Phi(), pcz_fix - pczguess); - plotter->Fill1D("pczfix-sx3pczguess_A1C2", 200, -100, 100, pcz_fix - pczguess); - plotter->Fill2D("pczfix_vs_sx3pczguess_A1C2_strip" + std::to_string(sx3event.ch2), 300, -200, 200, 600, -200, 200, pczguess, pcevent.pos.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(), "Vertex_Reconstruction"); + 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, "Z_Reconstruction"); + plotter->Fill2D("pczfix_vs_sx3pczguess_int_A1C2", 600, -200, 200, 600, -200, 200, pcz_guess_int, pcz_fix, "Z_Reconstruction"); + plotter->Fill2D("pczguess_vs_int", 600, -200, 200, 600, -200, 200, pcz_guess_int, pczguess, "Z_Reconstruction"); + plotter->Fill1D("pczguess_vs_int_residualsx3", 200, -50, 50, pcz_guess_int - pczguess, "Residuals"); + 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()); - 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_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_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, "PID_dE_E"); 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_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_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, "PID_dE_E"); } 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_strip" + std::to_string(sx3event.ch2), 300, -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(), "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_strip" + std::to_string(sx3event.ch2), 300, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z()); + 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(), "Z_Reconstruction"); 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->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->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(), "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); + plotter->Fill2D("sx3E_vs_sx3z", 400, 0, 30, 300, 0, 200, sx3event.Energy1, sx3z, "Kinematics_Angles"); - plotter->Fill2D("pcdEA_vs_sx3z", 800, 0, 20000, 300, 0, 200, pcevent.Energy1, sx3z); - plotter->Fill2D("pcdEC_vs_sx3z", 800, 0, 20000, 300, 0, 200, pcevent.Energy2, 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, "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("pcdEC_vs_sx3z" + std::to_string(sx3event.ch2), 800, 0, 20000, 300, 0, 200, pcevent.Energy2, 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, "Kinematics_Angles"); - plotter->Fill2D("pcdE2A_vs_sx3z", 800, 0, 20000, 300, 0, 200, pcevent.Energy1 * sinTheta, sx3z); - plotter->Fill2D("pcdE2C_vs_sx3z", 800, 0, 20000, 300, 0, 200, pcevent.Energy2 * sinTheta, sx3z); - plotter->Fill2D("phi_vs_stripnum", 180, -180, 180, 48, 0, 48, pcevent.pos.Phi() * 180. / M_PI, sx3event.ch2); - plotter->Fill2D("E_theta_AnodeSX3", 400, -20, 180, 300, 0, 15, sx3theta * 180 / M_PI, sx3event.Energy1); + 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, "Kinematics_Angles"); + 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, "Kinematics_Angles"); } 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 QQQ_Events, s { for (auto qqqevent : QQQ_Events) { - plotter->Fill1D("dt_pcA_qqqR", 640, -2000, 2000, qqqevent.Time1 - pcevent.Time1); - plotter->Fill2D("dt_pcA_qqqR_vs_qqqRE", 640, -2000, 2000, 400, 0, 30, qqqevent.Time1 - pcevent.Time1, qqqevent.Energy1); - plotter->Fill1D("dt_pcC_qqqW", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2); - plotter->Fill2D("phiPC_vs_phiQQQ", 180, -360, 360, 180, -360, 360, qqqevent.pos.Phi() * 180 / M_PI, pcevent.pos.Phi() * 180 / M_PI); - double sinTheta = TMath::Sin((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()); /// TMath::Sin((TVector3(51.5,0,128.) - TVector3(0,0,85)).Theta()); + 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, "Timing"); + 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, "Kinematics_Angles"); + double sinTheta = TMath::Sin((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()); TVector3 x2(pcevent.pos); TVector3 x1(qqqevent.pos); @@ -1512,74 +1513,67 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s 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; - // 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 phicut = qqqevent.pos.Phi() <= pcevent.pos.Phi() + TMath::Pi() / 4. && qqqevent.pos.Phi() >= pcevent.pos.Phi() - TMath::Pi() / 4.; if (lowercut_cath && phicut) { - plotter->Fill1D("dt_pcA_qqqR_pidlow_PC1", 640, -2000, 2000, qqqevent.Time1 - pcevent.Time1); - plotter->Fill2D("dt_pcA_qqqR_vs_qqqRE_pidlow_PC1", 640, -2000, 2000, 400, 0, 30, qqqevent.Time1 - pcevent.Time1, qqqevent.Energy1); - plotter->Fill1D("dt_pcC_qqqW_pidlow_PC1", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2); + 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, "Timing"); + plotter->Fill1D("dt_pcC_qqqW_pidlow_PC1", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2, "Timing"); } 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); - plotter->Fill2D("dE_E_CathodeQQQR", 400, 0, 30, 800, 0, 10000, qqqevent.Energy2, pcevent.Energy2); + { + plotter->Fill2D("dE_E_AnodeQQQR", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E"); + plotter->Fill2D("dE_E_CathodeQQQR", 400, 0, 30, 800, 0, 10000, qqqevent.Energy2, pcevent.Energy2, "PID_dE_E"); 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) - 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) - 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) - 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) { - 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_CathodeQQQW_TC1PC1_pidlow" + std::to_string(lowercut_cath), 400, 0, 30, 800, 0, 1000, qqqevent.Energy2, pcevent.Energy2 * sinTheta); - // 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); + 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, "PID_dE_E"); + 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("dE2_E_AnodeQQQR_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 4000, qqqevent.Energy1, pcevent.Energy1 * sinTheta); - plotter->Fill2D("dE2_E_CathodeQQQR_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 1000, qqqevent.Energy2, pcevent.Energy2 * sinTheta); - plotter->Fill2D("dEC_vs_dEA_TC1_PC" + std::to_string(phicut), 800, 0, 40000, 800, 0, 10000, pcevent.Energy1, pcevent.Energy2); - 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("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, "PID_dE_E"); + 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, "Timing"); - plotter->Fill1D("dt_pcA_qqqR_timecut", 640, -2000, 2000, qqqevent.Time1 - pcevent.Time1); - plotter->Fill1D("dt_pcC_qqqW_timecut", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2); - 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("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", 90, 0, 90, 400, 0, 20000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy1 * sinTheta); - 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->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, "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, "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, "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, "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, "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("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("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("dE_theta_CathodeQQQR", 75, 0, 90, 800, 0, 10000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy2); - 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_zoomin", 60, 0, 30, 800, 0, 3000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy2 * 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, "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, "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, "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_CathodeQQQR", 100, -180, 180, 800, 0, 10000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Phi() * 180 / M_PI, pcevent.Energy2); - plotter->Fill1D("PCZ", 800, -200, 200, pcevent.pos.Z(), "phicut"); - // 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"); + 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, "Kinematics_Angles"); + plotter->Fill1D("PCZ", 800, -200, 200, pcevent.pos.Z(), "Z_Reconstruction"); 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; - 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; - // 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(), "phicut"); + plotter->Fill2D("pczguess_vs_pc_int", 400, -200, 200, 600, -400, 400, pcz_guess_int, pcevent.pos.Z(), "Z_Reconstruction"); + if (pcevent.multi1 == 1 && pcevent.multi2 == 2) { double pcz_fix = pcfix_func.Eval(pcevent.pos.Z()); @@ -1589,69 +1583,69 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector QQQ_Events, s TVector3 r_rhoMin_fix = x1 + t_minimum * v; 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_AnodeQQQR_A1C2_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy1 * 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, "PID_dE_E"); - 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", 800, -400, 400, r_rhoMin_fix.Z(), "Vertex_Reconstruction"); + 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(), "Vertex_Reconstruction"); 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_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_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, "PID_dE_E"); } - plotter->Fill1D("pczfix_A1C2_1d_qqq", 600, -200, 200, pcz_fix); - plotter->Fill2D("pczfix_vs_qqqpczguess_A1C2", 600, -200, 200, 600, -200, 200, pcz_guess_int, pcz_fix); - plotter->Fill2D("pczguess_vs_pc_int_A1C2", 400, -200, 200, 600, -400, 400, pcz_guess_int, pcevent.pos.Z(), "phicut"); + 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, "Z_Reconstruction"); + 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; - // std::cout << path_length << std::endl; 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); - 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) - 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 { 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("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("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); + 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("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"); } + + if(BenchMark) + { + + } double qqqrho = qqqevent.pos.Perp(); double qqqz = (qqqevent.pos - TVector3(0, 0, source_vertex)).Z(); double tan_theta = qqqrho / qqqz; 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 tan_theta2 = qqqrho / qqqz2; 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_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"); + plotter->Fill2D("pczguess_vs_pc_int3", 180, 0, 200, 150, 0, 200, pcz_guess_int3, pcevent.pos.Z(), "Z_Reconstruction"); 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_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_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("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("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, "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, "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(), "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), // and the Kinematics object as a reference. void FillExcitationHistograms(HistPlotter *plotter, diff --git a/run_tr.sh b/run_tr.sh index 687009c..bb868c8 100644 --- a/run_tr.sh +++ b/run_tr.sh @@ -42,20 +42,20 @@ if [[ 1 -eq 0 ]]; then fi # --- Block 2: 27Al Alpha+Gas Runs (9, 12) --- -if [[ 1 -eq 0 ]]; then +if [[ 1 -eq 1 ]]; then export DATASET="27Al" export PREFIX="Run_" export OUT_DIR="Output_a" 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" + unset timecut_high export source_vertex=53.44; export timecut_low=400.0; process_run 12 "$slope" unset timecut_low - unset timecut_high fi # --- Block 3: 27Al Protons+Gas Runs (15, 17-22) --- -if [[ 1 -eq 0 ]]; then +if [[ 1 -eq 1 ]]; then export DATASET="27Al" export PREFIX="Run_" export OUT_DIR="Output_p" @@ -77,7 +77,7 @@ if [[ 1 -eq 0 ]]; then fi # --- Block 5: 17F Alpha Run with Gas (18-21) --- -if [[ 1 -eq 1 ]]; then +if [[ 1 -eq 0 ]]; then export DATASET="17F" export PREFIX="SourceRun_" export OUT_DIR="Output_a" diff --git a/scan_slope_runs.sh b/scan_slope_runs.sh index 12f505a..b5feb16 100644 --- a/scan_slope_runs.sh +++ b/scan_slope_runs.sh @@ -20,8 +20,8 @@ process_run() { export -f process_run 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 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 68 70 72 74 76 78 80 82 84 86 88 90; do slope=$(awk "BEGIN{printf \"%.2f\", $slope_x100/100}") 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) export DATASET="27Al" PREFIX="Run_" 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" unset timecut_low unset timecut_high diff --git a/scratch/plot_slope_scan.C b/scratch/plot_slope_scan.C index 018d549..f726801 100644 --- a/scratch/plot_slope_scan.C +++ b/scratch/plot_slope_scan.C @@ -8,339 +8,516 @@ #include #include #include +#include +#include void plot_slope_scan() { // ========================================================================= // --- Configuration --- // ========================================================================= - - // Set to true to inspect each Gaussian fit. Double-click the canvas to advance. bool checkFits = false; 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}; - // const int NRUNS = 6; + // 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, 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}; // const int NRUNS = 2; - int runs[] = {18, 19, 20, 21}; - const int NRUNS = 4; + // int runs[] = {18, 19, 20, 21}; + // const int NRUNS = 4; Int_t runColors[6] = {kRed + 1, kBlue + 2, kGreen + 2, kOrange + 1, kViolet + 2, kCyan + 2}; - // --- Initialize Graphs --- - TGraph *gStdDevRun[NRUNS]; - TGraph *gMeanRun[NRUNS]; + // --- Initialize Graphs for both SX3 and QQQ --- + TGraph *gStdDevRun_sx3[NRUNS], *gMeanRun_sx3[NRUNS]; + TGraph *gStdDevRun_qqq[NRUNS], *gMeanRun_qqq[NRUNS]; + for (int ir = 0; ir < NRUNS; ir++) { - gStdDevRun[ir] = new TGraph(N); - gMeanRun[ir] = new TGraph(N); + gStdDevRun_sx3[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 --- for (int ir = 0; ir < NRUNS; ir++) { 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]); TFile *f = TFile::Open(path); - if (!f || f->IsZombie()) continue; - TH1 *t1 = (TH1 *)f->Get("pczfix-sx3pczguess_A1C2"); - if (t1) + TH1 *t1_s = (TH1 *)f->Get("pczfix-sx3pczguess_A1C2"); + if (t1_s) { - t1->SetDirectory(nullptr); - h1[ir][i] = t1; + t1_s->SetDirectory(nullptr); + 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(); } } // ========================================================================= - // --- Calculate Mean & StdDev with Bounded Gaussian Fit --- + // --- Helper Lambda for Gaussian Fitting --- // ========================================================================= - std::vector nPtsRun(NRUNS, 0); - int nPtsSum = 0; - double fitWindow = 15.0; // Window around the peak to fit (mm) - - TCanvas *cDebug = nullptr; - if (checkFits) + auto fitPeak = [](TH1 *h, double fitWindow, bool debug, TCanvas *c, double &outMean, double &outStdDev) { - 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 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++) { - TH1 *hsum = nullptr; + TH1 *hsum_sx3 = nullptr, *hsum_qqq = nullptr; 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; - 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); - gMeanRun[ir]->SetPoint(nPtsRun[ir]++, slopes[i], mean); - - 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) + // Fit SX3 + if (fitPeak(h1_sx3[ir][i], fitWin, checkFits, cDebug, mean, stddev)) { - hsum->SetTitle(TString::Format("SUMMED RUNS | Slope %.2f", slopes[i])); - hsum->GetXaxis()->SetRangeUser(peakXSum - 50, peakXSum + 50); - - // Scale Y-axis for the sum + 20% headroom - hsum->SetMaximum(hsum->GetMaximum() * 1.2); - - hsum->Draw(); - cDebug->Modified(); - cDebug->Update(); - printf("Visual Check: SUMMED | Slope %.2f. Double-click canvas to continue...\n", slopes[i]); - while (cDebug->WaitPrimitive()) - ; + gStdDevRun_sx3[ir]->SetPoint(nPtsRun_sx3[ir], slopes[i], stddev); + gMeanRun_sx3[ir]->SetPoint(nPtsRun_sx3[ir]++, slopes[i], mean); + if (!hsum_sx3) + hsum_sx3 = (TH1 *)h1_sx3[ir][i]->Clone(TString::Format("hsum_sx3_%.2f", slopes[i])); + else + hsum_sx3->Add(h1_sx3[ir][i]); + } + if (runs[ir] == 20 || runs[ir] == 21) + continue; + // Fit QQQ + 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) - { - delete cDebug; // Clean up debug canvas before plotting final results - } + if (cDebug) + delete cDebug; // ========================================================================= - // ---- 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); - c1->Divide(2, 1); + TCanvas *c1 = new TCanvas("c_stats", "1D Residual Stats vs Slope", 0, 0, 1800, 800); - c1->cd(1); - gPad->SetGrid(1, 1); - TLegend *leg1 = new TLegend(0.78, 0.45, 0.97, 0.88); - leg1->SetBorderSize(1); - TMultiGraph *mg1 = new TMultiGraph("mg_stddev", "StdDev(pczfix - sx3pczguess) vs Slope;Slope;StdDev (mm)"); + // 1. Create a Master Legend in the top 8% of the canvas + c1->cd(); + TLegend *leg1_master = new TLegend(0.1, 0.92, 0.9, 1.0); + leg1_master->SetNColumns(6); + leg1_master->SetBorderSize(0); + leg1_master->SetFillStyle(0); + leg1_master->SetTextSize(0.03); for (int ir = 0; ir < NRUNS; ir++) { - if (nPtsRun[ir] == 0) - continue; - gStdDevRun[ir]->SetLineColor(runColors[ir]); - gStdDevRun[ir]->SetMarkerColor(runColors[ir]); - gStdDevRun[ir]->SetLineWidth(2); - gStdDevRun[ir]->SetMarkerStyle(20 + ir); - mg1->Add(gStdDevRun[ir], "PL"); - leg1->AddEntry(gStdDevRun[ir], TString::Format("run %d", runs[ir]), "lp"); + TLine *dSX3 = new TLine(); + dSX3->SetLineColor(runColors[ir]); + dSX3->SetLineWidth(2); + dSX3->SetLineStyle(kSolid); + leg1_master->AddEntry(dSX3, TString::Format("Run %d (SX3)", runs[ir]), "l"); } - if (nPtsSum > 0) + for (int ir = 0; ir < NRUNS; ir++) { - gStdDevSum->SetLineColor(kBlack); - gStdDevSum->SetLineWidth(3); - gStdDevSum->SetLineStyle(kDashed); - gStdDevSum->SetMarkerStyle(29); - mg1->Add(gStdDevSum, "PL"); - leg1->AddEntry(gStdDevSum, "sum", "lp"); + TLine *dQQQ = new TLine(); + dQQQ->SetLineColor(runColors[ir]); + dQQQ->SetLineWidth(2); + dQQQ->SetLineStyle(kDashed); + leg1_master->AddEntry(dQQQ, TString::Format("Run %d (QQQ)", runs[ir]), "l"); + } + 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"); - leg1->Draw(); double bestSlope = 0, bestStdDev = 1e9; - for (int i = 0; i < gStdDevRun[0]->GetN(); i++) + if (gStdDevRun_sx3[0]->GetN() > 0) { - double x, y; - gStdDevRun[0]->GetPoint(i, x, y); - if (y > 0 && y < bestStdDev) + for (int i = 0; i < gStdDevRun_sx3[0]->GetN(); i++) { - bestStdDev = y; - bestSlope = x; + double x, y; + 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(), - 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); - - c1->cd(2); + // --- Pad 2: Mean --- + pad1_grid->cd(2); gPad->SetGrid(1, 1); - TLegend *leg2 = new TLegend(0.78, 0.45, 0.97, 0.88); - leg2->SetBorderSize(1); - TMultiGraph *mg2 = new TMultiGraph("mg_mean", "Mean(pczfix - sx3pczguess) vs Slope;Slope;Mean Offset (mm)"); + TMultiGraph *mg2 = new TMultiGraph("mg_mean", "Mean of Residuals vs Slope;Slope;Mean Offset (mm)"); 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; - gMeanRun[ir]->SetLineColor(runColors[ir]); - gMeanRun[ir]->SetMarkerColor(runColors[ir]); - gMeanRun[ir]->SetLineWidth(2); - gMeanRun[ir]->SetMarkerStyle(20 + ir); - mg2->Add(gMeanRun[ir], "PL"); - leg2->AddEntry(gMeanRun[ir], TString::Format("run %d", runs[ir]), "lp"); - } - if (nPtsSum > 0) - { - gMeanSum->SetLineColor(kBlack); - gMeanSum->SetLineWidth(3); - gMeanSum->SetLineStyle(kDashed); - gMeanSum->SetMarkerStyle(29); - mg2->Add(gMeanSum, "PL"); - leg2->AddEntry(gMeanSum, "sum", "lp"); + if (nPtsRun_qqq[ir] > 0) + { + gMeanRun_qqq[ir]->SetLineColor(runColors[ir]); + gMeanRun_qqq[ir]->SetLineWidth(2); + gMeanRun_qqq[ir]->SetLineStyle(kDashed); + gMeanRun_qqq[ir]->SetMarkerStyle(24); + gMeanRun_qqq[ir]->SetMarkerColor(runColors[ir]); + mg2->Add(gMeanRun_qqq[ir], "PL"); + } } mg2->Draw("A"); - TLine *lz = new TLine(slopes[0], 0, slopes[N - 1], 0); lz->SetLineColor(kGray + 2); lz->SetLineStyle(kDotted); lz->Draw(); - leg2->Draw(); 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); - c2->Divide(4, 4); + TCanvas *c2 = new TCanvas("c_perslope_1d", "Per-slope 1D Residuals (SX3 & QQQ)", 0, 50, 2000, 1600); + + 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++) { - c2->cd(i + 1); + pad2_grid->cd(i + 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; for (int ir = 0; ir < NRUNS; ir++) { - if (!h1[ir][i] || h1[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) + if (h1_sx3[ir][i] && h1_sx3[ir][i]->GetEntries() > 0) { - 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; for (int ir = 0; ir < NRUNS; ir++) { - if (!h1[ir][i] || h1[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) + if (h1_sx3[ir][i] && h1_sx3[ir][i]->GetEntries() > 0) { - // Set the pad's Y-axis using the global max + 15% buffer so the legend fits cleanly - h1[ir][i]->SetMaximum(maxY * 1.15); - h1[ir][i]->Draw("HIST"); - first = false; + h1_sx3[ir][i]->SetLineColor(runColors[ir]); + h1_sx3[ir][i]->SetLineWidth(2); + h1_sx3[ir][i]->SetLineStyle(kSolid); + 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"); + + // ========================================================================= + // ---- 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->Update(); - while (gPad->WaitPrimitive()) - ; + c3->Modified(); + c3->Update(); } \ No newline at end of file