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 trying to benchmeark A1Cn events against a1c2 and failing due to flaww in logic
	modified:   run_tr.sh
	modified:   scratch/plot_slope_scan.C
This commit is contained in:
Vignesh Sitaraman 2026-06-11 16:51:41 -04:00
parent 20f8ad227c
commit f7248f6a4e
4 changed files with 222 additions and 51 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
Armory/PC_StepLadder_Correction.h
View File

@ -46,16 +46,16 @@ double model_a1c1(double* x, double* p) {
return result; return result;
} }
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 = 0.40; // /* 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;
else if(TMath::Abs(y[0]) < 85.2/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor*2; // else if(TMath::Abs(y[0]) < 85.2/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor*2;
else result=y[0]/slope-TMath::Sign(1.0,y[0])*factor**/; // else result=y[0]/slope-TMath::Sign(1.0,y[0])*factor**/;
return result+40; // return result+40;
} // }
/*void testmodel() { /*void testmodel() {

218
TrackRecon.C
View File

@ -44,6 +44,7 @@ bool doPCSX3ClusterAnalysis = true;
bool doPCQQQClusterAnalysis = true; bool doPCQQQClusterAnalysis = true;
bool doOldAnalysis = false; bool doOldAnalysis = false;
bool do27AlapAnalysis = false; bool do27AlapAnalysis = false;
bool BenchMark = true;
double source_vertex = 53; // 53 double source_vertex = 53; // 53
const double qqq_z = 105.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;
@ -143,13 +144,14 @@ void protonAlphaHistograms(HistPlotter *plotter, std::vector<Event> QQQ_Events,
void miscHistograms_oneWire(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<std::vector<std::tuple<int, double, double>>> aClusters); void miscHistograms_oneWire(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<std::vector<std::tuple<int, double, double>>> aClusters);
void protonMiscHistograms(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events); void protonMiscHistograms(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events);
void protonMiscHistograms_sx3(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events); void protonMiscHistograms_sx3(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events);
void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events); void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events,
void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events); const std::vector<std::vector<std::tuple<int, double, double>>> &aClusters, const std::vector<std::vector<std::tuple<int, double, double>>> &cClusters);
void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events,
const std::vector<std::vector<std::tuple<int, double, double>>> &aClusters, const std::vector<std::vector<std::tuple<int, double, double>>> &cClusters);
void TrackRecon::Begin(TTree * /*tree*/) void TrackRecon::Begin(TTree * /*tree*/)
{ {
pcfix_func.SetNpx(100000); pcfix_func.SetNpx(100000);
///// ---------Set Environment Variables--------- ///// ///// ---------Set Environment Variables--------- /////
TString option = GetOption(); TString option = GetOption();
if (option != "") if (option != "")
@ -598,8 +600,8 @@ Bool_t TrackRecon::Process(Long64_t entry)
// 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 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
@ -1024,11 +1026,11 @@ Bool_t TrackRecon::Process(Long64_t entry)
if (doPCSX3ClusterAnalysis) if (doPCSX3ClusterAnalysis)
{ {
PCSX3ClusterAnalysis(plotter, QQQ_Events, SX3_Events, PC_Events); PCSX3ClusterAnalysis(plotter, QQQ_Events, SX3_Events, PC_Events, aClusters, cClusters);
} }
if (doPCQQQClusterAnalysis) if (doPCQQQClusterAnalysis)
{ {
PCQQQClusterAnalysis(plotter, QQQ_Events, SX3_Events, PC_Events); PCQQQClusterAnalysis(plotter, QQQ_Events, SX3_Events, PC_Events, aClusters, cClusters);
} }
///////////////////nA analysis using pseudo-wire (GetPseudoWire + getClosestWirePosAtWirePhi)/////////////////// ///////////////////nA analysis using pseudo-wire (GetPseudoWire + getClosestWirePosAtWirePhi)///////////////////
@ -1362,7 +1364,9 @@ 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,
const std::vector<std::vector<std::tuple<int, double, double>>> &aClusters, const std::vector<std::vector<std::tuple<int, double, double>>> &cClusters)
{
for (auto pcevent : PC_Events) for (auto pcevent : PC_Events)
{ {
bool PCSX3TimeCut = false; bool PCSX3TimeCut = false;
@ -1403,8 +1407,8 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
plotter->Fill2D("pcz_vs_pcphi_TimeCut", 600, -200, 200, 120, -360, 360, pcevent.pos.Z(), pcevent.pos.Phi() * 180 / M_PI, "Z_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; double sx3rho = 88.0;
double sx3z = sx3event.pos.Z(); 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), "PID_dE_E"); plotter->Fill2D("dE2_E_Anodesx3B", 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1 * TMath::Sin(calcsx3theta), "PID_dE_E");
@ -1423,11 +1427,12 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
plotter->Fill1D("VertexReconZ_SX3", 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("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(), "Timing"); 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("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("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(), "Z_Reconstruction");
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)
{ {
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(), "Z_Reconstruction");
@ -1466,7 +1471,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
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_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(), "Z_Reconstruction"); 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"); 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.);
@ -1491,11 +1496,95 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
{ {
plotter->Fill1D("PCZ_sx3", 800, -200, 200, pcevent.pos.Z(), "Z_Reconstruction"); plotter->Fill1D("PCZ_sx3", 800, -200, 200, pcevent.pos.Z(), "Z_Reconstruction");
} }
}
//-----------------------Benchmarking Method for Source Runs------------------------//
if (BenchMark && aClusters.size() == 1 && cClusters.size() == 1)
{
const auto &aCl = aClusters.front();
const auto &cCl = cClusters.front();
auto vertexFromPCPoint = [&x1](const TVector3 &x2f)
{
TVector3 vf = x2f - x1;
double tm = -1.0 * (x1.X() * vf.X() + x1.Y() * vf.Y()) / (vf.X() * vf.X() + vf.Y() * vf.Y());
return TVector3(x1 + tm * vf);
};
auto fillSuite = [&](const std::string &tag, double pcz_method, const TVector3 &vtx)
{
plotter->Fill1D("Benchmark_SX3_VertexZ_" + tag, 800, -400, 400, vtx.Z(), "Benchmark_SX3");
plotter->Fill1D("Benchmark_SX3_VertexZ_" + tag + "_TC" + std::to_string(PCSX3TimeCut) + "_PC" + std::to_string(phicut), 800, -400, 400, vtx.Z(), "Benchmark_SX3");
plotter->Fill2D("Benchmark_SX3_VertexXY_" + tag, 200, -100, 100, 200, -100, 100, vtx.X(), vtx.Y(), "Benchmark_SX3");
plotter->Fill1D("Benchmark_SX3_PCZ_" + tag, 600, -200, 200, pcz_method, "Benchmark_SX3");
plotter->Fill2D("Benchmark_SX3_PCZ_vs_sx3pczguess_" + tag, 600, -200, 200, 600, -200, 200, pcz_guess_int, pcz_method, "Benchmark_SX3");
plotter->Fill2D("Benchmark_SX3_PCZresidual_vs_pczguess_" + tag, 600, -200, 200, 200, -100, 100, pczguess, pcz_method - pczguess, "Benchmark_SX3");
plotter->Fill1D("Benchmark_SX3_PCZ-sx3pczguess_" + tag, 200, -100, 100, pcz_method - pczguess, "Benchmark_SX3");
plotter->Fill1D("Benchmark_SX3_PCZ-sx3pczint_" + tag, 200, -100, 100, pcz_method - pcz_guess_int, "Benchmark_SX3");
};
auto fillVsRef = [&](const std::string &tag, double pcz_method, const TVector3 &vtx, double pcz_ref, const TVector3 &vtx_ref)
{
plotter->Fill2D("Benchmark_SX3_PCZ_" + tag + "_vs_ref", 400, -200, 200, 400, -200, 200, pcz_ref, pcz_method, "Benchmark_SX3_ref");
plotter->Fill1D("Benchmark_SX3_PCZ_" + tag + "_minus_ref", 400, -100, 100, pcz_method - pcz_ref, "Benchmark_SX3_ref");
plotter->Fill2D("Benchmark_SX3_VertexZ_" + tag + "_vs_ref", 400, -200, 200, 400, -200, 200, vtx_ref.Z(), vtx.Z(), "Benchmark_SX3_ref");
plotter->Fill1D("Benchmark_SX3_VertexZ_" + tag + "_minus_ref", 400, -100, 100, vtx.Z() - vtx_ref.Z(), "Benchmark_SX3_ref");
};
// cathode charge-division reference on this event (the A1C2 recipe)
double pcz_ref = pcfix_func.Eval(pcevent.pos.Z());
TVector3 vtx_ref = vertexFromPCPoint(TVector3(pcevent.pos.X(), pcevent.pos.Y(), pcz_ref));
// anode-only PC point, oneWire method: pseudo-wire interpolated to the
// SX3 hit's phi, with the same quality cuts as miscHistograms_oneWire
auto [apwire_bm, apSumE_bm, apMaxE_bm, apTSMaxE_bm] = pwinstance.GetPseudoWire(aCl, "ANODE");
TVector3 pc_anodeOnly = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, sx3event.pos.Phi());
TVector3 vtx_anodeOnly = vertexFromPCPoint(pc_anodeOnly);
bool anodeOnlyGood = vtx_anodeOnly.Perp() <= 6.0 && vtx_anodeOnly.Z() >= -173.6 && vtx_anodeOnly.Z() <= 100 && phicut;
if (pcevent.multi1 == 1 && pcevent.multi2 == 2)
{
fillSuite("A1C2", pcz_ref, vtx_ref); // baseline with identical binning
// ---- A1C1 emulation: crossover with only the max-E cathode wire (Cmax),
// used directly with no z-model ----
auto cMaxWire = *std::max_element(cCl.begin(), cCl.end(),
[](const std::tuple<int, double, double> &w1, const std::tuple<int, double, double> &w2)
{ return std::get<1>(w1) < std::get<1>(w2); });
std::vector<std::tuple<int, double, double>> cOne = {cMaxWire};
auto [xo_a1c1, alpha_a1c1, apSumE1, cpSumE1, apMaxE1, cpMaxE1, apTSMaxE1, cpTSMaxE1] = pwinstance.FindCrossoverProperties(aCl, cOne);
if (alpha_a1c1 != 9999999 && apSumE1 != -1)
{
TVector3 vtx_a1c1 = vertexFromPCPoint(xo_a1c1);
fillSuite("A1C1", xo_a1c1.Z(), vtx_a1c1);
fillVsRef("A1C1", xo_a1c1.Z(), vtx_a1c1, pcz_ref, vtx_ref);
}
// ---- A1C0 emulation: oneWire method, cathodes ignored ----
if (anodeOnlyGood)
{
fillSuite("A1C0", pc_anodeOnly.Z(), vtx_anodeOnly);
fillVsRef("A1C0", pc_anodeOnly.Z(), vtx_anodeOnly, pcz_ref, vtx_ref);
}
}
// ---- A2C0 emulation: 2-wire anode pseudo-wire, oneWire method ----
if (pcevent.multi1 == 2 && pcevent.multi2 == 2)
{
fillSuite("A2C2", pcz_ref, vtx_ref); // reference population for A2C0
if (anodeOnlyGood)
{
fillSuite("A2C0", pc_anodeOnly.Z(), vtx_anodeOnly);
fillVsRef("A2C0", pc_anodeOnly.Z(), vtx_anodeOnly, pcz_ref, vtx_ref);
}
}
}
}
} }
} }
void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events) void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events,
const std::vector<std::vector<std::tuple<int, double, double>>> &aClusters, const std::vector<std::vector<std::tuple<int, double, double>>> &cClusters)
{ {
for (auto pcevent : PC_Events) for (auto pcevent : PC_Events)
{ {
@ -1505,7 +1594,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
plotter->Fill2D("dt_pcA_qqqR_vs_qqqRE", 640, -2000, 2000, 400, 0, 30, qqqevent.Time1 - pcevent.Time1, qqqevent.Energy1, "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->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"); 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()); 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);
@ -1524,7 +1613,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
plotter->Fill1D("dt_pcC_qqqW_pidlow_PC1", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2, "Timing"); plotter->Fill1D("dt_pcC_qqqW_pidlow_PC1", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2, "Timing");
} }
if (timecut) if (timecut)
{ {
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, "PID_dE_E");
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, "PID_dE_E");
if (pcevent.multi1 == 1 && pcevent.multi2 == 2) if (pcevent.multi1 == 1 && pcevent.multi2 == 2)
@ -1546,7 +1635,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
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_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("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("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->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, "Timing"); 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->Fill1D("dt_pcC_qqqW_timecut", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2, "Timing");
@ -1622,10 +1711,89 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
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("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) //-----------------------Benchmarking Method for Source Runs------------------------//
{
if (BenchMark && aClusters.size() == 1 && cClusters.size() == 1)
} {
const auto &aCl = aClusters.front();
const auto &cCl = cClusters.front();
auto vertexFromPCPoint = [&x1](const TVector3 &x2f)
{
TVector3 vf = x2f - x1;
double tm = -1.0 * (x1.X() * vf.X() + x1.Y() * vf.Y()) / (vf.X() * vf.X() + vf.Y() * vf.Y());
return TVector3(x1 + tm * vf);
};
auto fillSuite = [&](const std::string &tag, double pcz_method, const TVector3 &vtx)
{
plotter->Fill1D("Benchmark_QQQ_VertexZ_" + tag, 800, -400, 400, vtx.Z(), "Benchmark_QQQ");
plotter->Fill1D("Benchmark_QQQ_VertexZ_" + tag + "_TC" + std::to_string(timecut) + "_PC" + std::to_string(phicut), 800, -400, 400, vtx.Z(), "Benchmark_QQQ");
plotter->Fill2D("Benchmark_QQQ_VertexXY_" + tag, 200, -100, 100, 200, -100, 100, vtx.X(), vtx.Y(), "Benchmark_QQQ");
plotter->Fill1D("Benchmark_QQQ_PCZ_" + tag, 600, -200, 200, pcz_method, "Benchmark_QQQ");
plotter->Fill2D("Benchmark_QQQ_PCZ_vs_qqqpczguess_" + tag, 600, -200, 200, 600, -200, 200, pcz_guess_int, pcz_method, "Benchmark_QQQ");
plotter->Fill2D("Benchmark_QQQ_PCZresidual_vs_pczguess_" + tag, 600, -200, 200, 200, -100, 100, pcz_guess_37, pcz_method - pcz_guess_37, "Benchmark_QQQ");
plotter->Fill1D("Benchmark_QQQ_PCZ-qqqpczguess_" + tag, 200, -100, 100, pcz_method - pcz_guess_37, "Benchmark_QQQ");
plotter->Fill1D("Benchmark_QQQ_PCZ-qqqpczint_" + tag, 200, -100, 100, pcz_method - pcz_guess_int, "Benchmark_QQQ");
};
auto fillVsRef = [&](const std::string &tag, double pcz_method, const TVector3 &vtx, double pcz_ref, const TVector3 &vtx_ref)
{
plotter->Fill2D("Benchmark_QQQ_PCZ_" + tag + "_vs_ref", 400, -200, 200, 400, -200, 200, pcz_ref, pcz_method, "Benchmark_QQQ_ref");
plotter->Fill1D("Benchmark_QQQ_PCZ_" + tag + "_minus_ref", 400, -100, 100, pcz_method - pcz_ref, "Benchmark_QQQ_ref");
plotter->Fill2D("Benchmark_QQQ_VertexZ_" + tag + "_vs_ref", 400, -200, 200, 400, -200, 200, vtx_ref.Z(), vtx.Z(), "Benchmark_QQQ_ref");
plotter->Fill1D("Benchmark_QQQ_VertexZ_" + tag + "_minus_ref", 400, -100, 100, vtx.Z() - vtx_ref.Z(), "Benchmark_QQQ_ref");
};
// cathode charge-division reference on this event (the A1C2 recipe)
double pcz_ref = pcfix_func.Eval(pcevent.pos.Z());
TVector3 vtx_ref = vertexFromPCPoint(TVector3(pcevent.pos.X(), pcevent.pos.Y(), pcz_ref));
// anode-only PC point, oneWire method: pseudo-wire interpolated to the
// QQQ hit's phi, with the same quality cuts as miscHistograms_oneWire
auto [apwire_bm, apSumE_bm, apMaxE_bm, apTSMaxE_bm] = pwinstance.GetPseudoWire(aCl, "ANODE");
TVector3 pc_anodeOnly = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, qqqevent.pos.Phi());
TVector3 vtx_anodeOnly = vertexFromPCPoint(pc_anodeOnly);
bool anodeOnlyGood = vtx_anodeOnly.Perp() <= 6.0 && vtx_anodeOnly.Z() >= -173.6 && vtx_anodeOnly.Z() <= 100 && phicut;
if (pcevent.multi1 == 1 && pcevent.multi2 == 2)
{
fillSuite("A1C2", pcz_ref, vtx_ref); // baseline with identical binning
// ---- A1C1 emulation: crossover with only the max-E cathode wire (Cmax),
// used directly with no z-model ----
auto cMaxWire = *std::max_element(cCl.begin(), cCl.end(),
[](const std::tuple<int, double, double> &w1, const std::tuple<int, double, double> &w2)
{ return std::get<1>(w1) < std::get<1>(w2); });
std::vector<std::tuple<int, double, double>> cOne = {cMaxWire};
auto [xo_a1c1, alpha_a1c1, apSumE1, cpSumE1, apMaxE1, cpMaxE1, apTSMaxE1, cpTSMaxE1] = pwinstance.FindCrossoverProperties(aCl, cOne);
if (alpha_a1c1 != 9999999 && apSumE1 != -1)
{
TVector3 vtx_a1c1 = vertexFromPCPoint(xo_a1c1);
fillSuite("A1C1", xo_a1c1.Z(), vtx_a1c1);
fillVsRef("A1C1", xo_a1c1.Z(), vtx_a1c1, pcz_ref, vtx_ref);
}
// ---- A1C0 emulation: oneWire method, cathodes ignored ----
if (anodeOnlyGood)
{
fillSuite("A1C0", pc_anodeOnly.Z(), vtx_anodeOnly);
fillVsRef("A1C0", pc_anodeOnly.Z(), vtx_anodeOnly, pcz_ref, vtx_ref);
}
}
// ---- A2C0 emulation: 2-wire anode pseudo-wire, oneWire method ----
if (pcevent.multi1 == 2 && pcevent.multi2 == 2)
{
fillSuite("A2C2", pcz_ref, vtx_ref); // reference population for A2C0
if (anodeOnlyGood)
{
fillSuite("A2C0", pc_anodeOnly.Z(), vtx_anodeOnly);
fillVsRef("A2C0", pc_anodeOnly.Z(), vtx_anodeOnly, pcz_ref, vtx_ref);
}
}
}
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;
@ -1644,7 +1812,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
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"); 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");
} }
} }
} }
} }
// 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.

13
run_tr.sh
View File

@ -55,14 +55,15 @@ if [[ 1 -eq 1 ]]; then
fi fi
# --- Block 3: 27Al Protons+Gas Runs (15, 17-22) --- # --- Block 3: 27Al Protons+Gas Runs (15, 17-22) ---
if [[ 1 -eq 1 ]]; then if [[ 1 -eq 0 ]]; then
export DATASET="27Al" export DATASET="27Al"
export PREFIX="Run_" export PREFIX="Run_"
export OUT_DIR="Output_p" export OUT_DIR="Output_p"
export source_vertex=-200.0 # Source on the entrance window export source_vertex=-200.0 # Source on the entrance window
echo "Starting parallel processing for 27Al proton runs..." echo "Starting parallel processing for 27Al proton runs..."
parallel --bar -j 6 process_run ::: 15 {17..22} process_run 17
# parallel --bar -j 6 process_run ::: 15 {17..22}
hadd -j 4 -k ${OUT_DIR}/Al_protons.root ${OUT_DIR}/results_run0{17..22}.root hadd -j 4 -k ${OUT_DIR}/Al_protons.root ${OUT_DIR}/results_run0{17..22}.root
fi fi
@ -77,17 +78,17 @@ 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 0 ]]; then if [[ 1 -eq 1 ]]; then
export DATASET="17F" export DATASET="17F"
export PREFIX="SourceRun_" export PREFIX="SourceRun_"
export OUT_DIR="Output_a" export OUT_DIR="Output_a"
echo "Processing 17F alpha runs with dynamic source vertices..." echo "Processing 17F alpha runs with dynamic source vertices..."
# Running sequentially since the source_vertex variable changes per run # Running sequentially since the source_vertex variable changes per run
# export source_vertex=53.44; process_run 18 export source_vertex=53.44; process_run 18
export source_vertex=14.24; process_run 19 export source_vertex=14.24; process_run 19
# export source_vertex=-24.96; process_run 20 export source_vertex=-24.96; process_run 20
# export source_vertex=-73.96; process_run 21 export source_vertex=-73.96; process_run 21
fi fi
# --- Block 6: 17F Proton Data --- # --- Block 6: 17F Proton Data ---

22
scratch/plot_slope_scan.C
View File

@ -25,9 +25,11 @@ void plot_slope_scan()
// 0.82, 0.84, 0.86, 0.88, 0.90}; // 0.82, 0.84, 0.86, 0.88, 0.90};
// const int N = 25; // 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}; // 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 // const int N = 12; // Updated to 14 to match the size of your slopes array
const int N = 1; // Updated to 14 to match the size of your slopes array
double slopes[] = {0.52};
int runs[] = {9, 12, 18, 19, 20, 21}; int runs[] = {9, 12, 18, 19, 20, 21};
const int NRUNS = 6; const int NRUNS = 6;
// int runs[] = {9, 12}; // int runs[] = {9, 12};
@ -272,12 +274,12 @@ void plot_slope_scan()
mg1->Draw("A"); mg1->Draw("A");
double bestSlope = 0, bestStdDev = 1e9; double bestSlope = 0, bestStdDev = 1e9;
if (gStdDevRun_sx3[0]->GetN() > 0) if (gStdDevRun_sx3[3]->GetN() > 0)
{ {
for (int i = 0; i < gStdDevRun_sx3[0]->GetN(); i++) for (int i = 0; i < gStdDevRun_sx3[3]->GetN(); i++)
{ {
double x, y; double x, y;
gStdDevRun_sx3[0]->GetPoint(i, x, y); gStdDevRun_sx3[3]->GetPoint(i, x, y);
if (y > 0 && y < bestStdDev) if (y > 0 && y < bestStdDev)
{ {
bestStdDev = y; bestStdDev = y;
@ -288,12 +290,12 @@ void plot_slope_scan()
lb->SetLineColor(kRed); lb->SetLineColor(kRed);
lb->SetLineStyle(kDashed); lb->SetLineStyle(kDashed);
lb->Draw(); lb->Draw();
printf(">>> Best slope (Run %d) = %.2f StdDev=%.2f mm\n", runs[0], bestSlope, bestStdDev); printf(">>> Best slope (Run %d) = %.2f StdDev=%.2f mm\n", runs[3], bestSlope, bestStdDev);
} }
// --- Pad 2: Mean --- // --- Pad 2: Mean ---
pad1_grid->cd(2); pad1_grid->cd(2);
gPad->SetGrid(1, 1); // gPad->SetGrid(1, 1);
TMultiGraph *mg2 = new TMultiGraph("mg_mean", "Mean of Residuals 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++) for (int ir = 0; ir < NRUNS; ir++)
@ -358,7 +360,7 @@ void plot_slope_scan()
TPad *pad2_grid = new TPad("pad2_grid", "Grid", 0.0, 0.0, 1.0, 0.92); TPad *pad2_grid = new TPad("pad2_grid", "Grid", 0.0, 0.0, 1.0, 0.92);
pad2_grid->Draw(); pad2_grid->Draw();
pad2_grid->Divide(5, 5); pad2_grid->Divide(3,4);
for (int i = 0; i < N; i++) for (int i = 0; i < N; i++)
{ {
@ -453,7 +455,7 @@ void plot_slope_scan()
TPad *pad3_grid = new TPad("pad3_grid", "Grid", 0.0, 0.0, 1.0, 0.92); TPad *pad3_grid = new TPad("pad3_grid", "Grid", 0.0, 0.0, 1.0, 0.92);
pad3_grid->Draw(); pad3_grid->Draw();
pad3_grid->Divide(5, 5); pad3_grid->Divide(3,4);
for (int i = 0; i < N; i++) for (int i = 0; i < N; i++)
{ {