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

modified: TrackRecon.C

Browse Source 
modified:   run_27Al.sh
	modified:   run_tr.sh
This commit is contained in:
Vignesh Sitaraman 2026-05-26 15:43:21 -04:00
parent 07fe2a5d38
commit 15773d4606
3 changed files with 285 additions and 313 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

535
TrackRecon.C
View File

@ -1,5 +1,10 @@
#define TrackRecon_cxx #define TrackRecon_cxx
// #define RAW_HISTOS
// #define VTX_GATES
#define AL_BEAM
// #define F_BEAM
Int_t colors[40] = { Int_t colors[40] = {
kBlack, kRed, kGreen, kBlue, kYellow, kMagenta, kCyan, kOrange, kBlack, kRed, kGreen, kBlue, kYellow, kMagenta, kCyan, kOrange,
kSpring, kTeal, kAzure, kViolet, kPink, kGray, kWhite, kSpring, kTeal, kAzure, kViolet, kPink, kGray, kWhite,
@ -33,8 +38,10 @@ Int_t colors[40] = {
bool process_alpha_proton_scattering = true; bool process_alpha_proton_scattering = true;
bool doPCSX3ClusterAnalysis = true; bool doPCSX3ClusterAnalysis = true;
bool doPCQQQClusterAnalysis = true; bool doPCQQQClusterAnalysis = true;
bool do27AlapAnalysis = true;
double source_vertex = 53; // 53 double source_vertex = 53; // 53
const double qqq_z = 100.0; const double qqq_z = 100.0;
double z_entrance = -174.3 - 9.7 - 100.0;
const double anode_gain = 1.5146e-5; // channels --> MeV const double anode_gain = 1.5146e-5; // channels --> MeV
std::string dataset; std::string dataset;
bool productionrun = false; bool productionrun = false;
@ -42,7 +49,16 @@ bool productionrun = false;
TF1 pcfix_func("func", model_invert, -200, 200); TF1 pcfix_func("func", model_invert, -200, 200);
TGraph *MeV_to_cm = NULL, *cm_to_MeV = NULL; TGraph *MeV_to_cm = NULL, *cm_to_MeV = NULL;
TGraph *MeV_to_cm_p = NULL, *cm_to_MeVp = NULL; TGraph *MeV_to_cm_p = NULL, *cm_to_MeVp = NULL;
TGraph *MeV_to_cm_beam = NULL, *cm_to_MeV_beam = NULL; TGraph *MeV_to_cm_27Al = NULL, *cm_to_MeV_27Al = NULL;
TGraph *MeV_to_cm_17F = NULL, *cm_to_MeV_17F = NULL;
// declaring masses for kinematics calculations
double mass_27Al = 26.981538;
double mass_4He = 4.002603;
double mass_1H = 1.007825;
double mass_30Si = 29.973770;
double mass_17F = 17.002095;
double mass_20Ne = 19.992440;
/* /*
double z_to_crossover_rho(double z) double z_to_crossover_rho(double z)
@ -127,20 +143,17 @@ void TrackRecon::Begin(TTree * /*tree*/)
pwinstance.ConstructGeo(); pwinstance.ConstructGeo();
// ---------------------------------------------------------
// 1. CRITICAL FIX: Initialize PC Arrays to Default (Raw)
// ---------------------------------------------------------
for (int i = 0; i < 48; i++) for (int i = 0; i < 48; i++)
{ {
pcSlope[i] = 1.0; // Default slope = 1 (preserves Raw energy) pcSlope[i] = 1.0;
pcIntercept[i] = 0.0; // Default intercept = 0 pcIntercept[i] = 0.0;
} }
if (getenv("DATASET")) if (getenv("DATASET"))
dataset = std::string(getenv("DATASET")); dataset = std::string(getenv("DATASET"));
if (getenv("source_vertex")) if (getenv("source_vertex"))
source_vertex = (double)std::atof(std::string(getenv("source_vertex")).c_str()); source_vertex = (double)std::atof(std::string(getenv("source_vertex")).c_str());
if(dataset == "17F" && getenv("PRODUCTION")) if (dataset == "17F" && getenv("PRODUCTION"))
productionrun = true; productionrun = true;
std::cout << "Dataset set to " << dataset << std::endl; std::cout << "Dataset set to " << dataset << std::endl;
std::cout << "source_vertex set to " << source_vertex << std::endl; std::cout << "source_vertex set to " << source_vertex << std::endl;
@ -238,8 +251,11 @@ void TrackRecon::Begin(TTree * /*tree*/)
cm_to_MeVp = new TGraph(MeV_to_cm_p->GetN(), MeV_to_cm_p->GetY(), MeV_to_cm_p->GetX()); cm_to_MeVp = new TGraph(MeV_to_cm_p->GetN(), MeV_to_cm_p->GetY(), MeV_to_cm_p->GetX());
// Add these alongside your existing proton and alpha tables // Add these alongside your existing proton and alpha tables
MeV_to_cm_beam = new TGraph("eloss_calculations/aluminum_lookup_80MeV.dat", "%lf %*lf %lf"); MeV_to_cm_27Al = new TGraph("eloss_calculations/aluminum_lookup_80MeV.dat", "%lf %*lf %lf");
cm_to_MeV_beam = new TGraph(MeV_to_cm_beam->GetN(), MeV_to_cm_beam->GetY(), MeV_to_cm_beam->GetX()); cm_to_MeV_27Al = new TGraph(MeV_to_cm_27Al->GetN(), MeV_to_cm_27Al->GetY(), MeV_to_cm_27Al->GetX());
MeV_to_cm_17F = new TGraph("eloss_calculations/fluorine_lookup_70MeV.dat", "%lf %*lf %lf");
cm_to_MeV_17F = new TGraph(MeV_to_cm_17F->GetN(), MeV_to_cm_17F->GetY(), MeV_to_cm_17F->GetX());
// cm_to_MeV.Eval(MeV_to_cm.Eval(detectedE)-PathLength) gives energy of particle before it traversed 'path length' // cm_to_MeV.Eval(MeV_to_cm.Eval(detectedE)-PathLength) gives energy of particle before it traversed 'path length'
} }
@ -323,7 +339,9 @@ Bool_t TrackRecon::Process(Long64_t entry)
if (id < 12) if (id < 12)
Fsx3.at(id).fillevent("BACK", sx3ch, value); Fsx3.at(id).fillevent("BACK", sx3ch, value);
Fsx3.at(id).ts = static_cast<double>(sx3.t[i]); Fsx3.at(id).ts = static_cast<double>(sx3.t[i]);
#ifdef RAW_HISTOS
plotter->Fill2D("sx3backs_all_raw", 100, 0, 100, 800, 0, 4096, gch, sx3.e[i]); plotter->Fill2D("sx3backs_all_raw", 100, 0, 100, 800, 0, 4096, gch, sx3.e[i]);
#endif
} }
else else
{ {
@ -400,10 +418,12 @@ Bool_t TrackRecon::Process(Long64_t entry)
// plotter->Fill2D("SX3CartesianPlot", 200, -100, 100, 200, -100, 100, 88.0*TMath::Cos(phi_n),88.0*TMath::Sin(phi_n), "hCalSX3"); // plotter->Fill2D("SX3CartesianPlot", 200, -100, 100, 200, -100, 100, 88.0*TMath::Cos(phi_n),88.0*TMath::Sin(phi_n), "hCalSX3");
plotter->Fill2D("SX3CartesianPlot" + std::to_string(id), 200, -100, 100, 200, -100, 100, 88.0 * TMath::Cos(phi_n), 88.0 * TMath::Sin(phi_n), "hCalSX3"); plotter->Fill2D("SX3CartesianPlot" + std::to_string(id), 200, -100, 100, 200, -100, 100, 88.0 * TMath::Cos(phi_n), 88.0 * TMath::Sin(phi_n), "hCalSX3");
} }
#ifdef RAW_HISTOS
if (det.valid && det.stripF != DEFAULT_NULL && det.stripB != DEFAULT_NULL) if (det.valid && det.stripF != DEFAULT_NULL && det.stripB != DEFAULT_NULL)
{ {
plotter->Fill2D("sx3backs_raw", 100, 0, 100, 800, 0, 8192, det.stripB + 4 * id, det.backE); plotter->Fill2D("sx3backs_raw", 100, 0, 100, 800, 0, 8192, det.stripB + 4 * id, det.backE);
} }
#endif
} }
} }
// return kTRUE; // return kTRUE;
@ -452,6 +472,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
bool PCCQQQTimeCut = false; bool PCCQQQTimeCut = false;
for (int i = 0; i < qqq.multi; i++) for (int i = 0; i < qqq.multi; i++)
{ {
#ifdef RAW_HISTOS
plotter->Fill2D("QQQ_Index_Vs_Energy", 16 * 8, 0, 16 * 8, 2000, 0, 16000, qqq.index[i], qqq.e[i], "hRawQQQ"); plotter->Fill2D("QQQ_Index_Vs_Energy", 16 * 8, 0, 16 * 8, 2000, 0, 16000, qqq.index[i], qqq.e[i], "hRawQQQ");
for (int j = 0; j < qqq.multi; j++) for (int j = 0; j < qqq.multi; j++)
@ -473,7 +494,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
plotter->Fill2D("QQQ_Vs_Cathode_Energy", 400, 0, 4000, 1000, 0, 16000, qqq.e[i], pc.e[k], "hRawQQQ"); plotter->Fill2D("QQQ_Vs_Cathode_Energy", 400, 0, 4000, 1000, 0, 16000, qqq.e[i], pc.e[k], "hRawQQQ");
} }
} }
#endif
for (int j = i + 1; j < qqq.multi; j++) for (int j = i + 1; j < qqq.multi; j++)
{ {
if (qqq.id[i] == qqq.id[j]) if (qqq.id[i] == qqq.id[j])
@ -511,8 +532,10 @@ Bool_t TrackRecon::Process(Long64_t entry)
continue; continue;
plotter->Fill1D("Wedgetime_Vs_Ringtime", 100, -1000, 1000, tWedge - tRing, "hTiming"); plotter->Fill1D("Wedgetime_Vs_Ringtime", 100, -1000, 1000, tWedge - tRing, "hTiming");
#ifdef RAW_HISTOS
plotter->Fill2D("RingE_vs_Index", 16 * 4, 0, 16 * 4, 1000, 0, 16000, chRing + qqq.id[i] * 16, eRing, "hRawQQQ"); plotter->Fill2D("RingE_vs_Index", 16 * 4, 0, 16 * 4, 1000, 0, 16000, chRing + qqq.id[i] * 16, eRing, "hRawQQQ");
plotter->Fill2D("WedgeE_vs_Index", 16 * 4, 0, 16 * 4, 1000, 0, 16000, chWedge + qqq.id[i] * 16, eWedge, "hRawQQQ"); plotter->Fill2D("WedgeE_vs_Index", 16 * 4, 0, 16 * 4, 1000, 0, 16000, chWedge + qqq.id[i] * 16, eWedge, "hRawQQQ");
#endif
plotter->Fill2D("WedgeE_Vs_RingECal", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ"); plotter->Fill2D("WedgeE_Vs_RingECal", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ");
if (qqqCalibValid[qqq.id[i]][chWedge][chRing]) if (qqqCalibValid[qqq.id[i]][chWedge][chRing])
@ -530,12 +553,11 @@ Bool_t TrackRecon::Process(Long64_t entry)
// we found a 12mm shift towards the vertex later --> 116 // we found a 12mm shift towards the vertex later --> 116
Event qqqevent(TVector3(rho * TMath::Cos(theta), rho * TMath::Sin(theta), qqq_z), eRingMeV, eWedgeMeV, tRing, tWedge, chRing + qqq.id[i] * 16, chWedge + qqq.id[i] * 16); Event qqqevent(TVector3(rho * TMath::Cos(theta), rho * TMath::Sin(theta), qqq_z), eRingMeV, eWedgeMeV, tRing, tWedge, chRing + qqq.id[i] * 16, chWedge + qqq.id[i] * 16);
Event qqqeventr(TVector3(rho * TMath::Cos(theta), rho * TMath::Sin(theta), qqq_z), eRing, eWedge, tRing, tWedge, chRing + qqq.id[i] * 16, chWedge + qqq.id[i] * 16); Event qqqeventr(TVector3(rho * TMath::Cos(theta), rho * TMath::Sin(theta), qqq_z), eRing, eWedge, tRing, tWedge, chRing + qqq.id[i] * 16, chWedge + qqq.id[i] * 16);
if (qqq.id[i] >= 0)
{ QQQ_Events.push_back(qqqevent);
QQQ_Events.push_back(qqqevent); QQQ_Events_Raw.push_back(qqqeventr);
QQQ_Events_Raw.push_back(qqqeventr); plotter->Fill2D("WedgeE_Vs_RingECal_selected", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ");
plotter->Fill2D("WedgeE_Vs_RingECal_selected", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ");
}
plotter->Fill2D("QQQCartesianPlot", 200, -100, 100, 200, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hCalQQQ"); plotter->Fill2D("QQQCartesianPlot", 200, -100, 100, 200, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hCalQQQ");
plotter->Fill2D("QQQCartesianPlot" + std::to_string(qqq.id[i]), 200, -100, 100, 200, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hCalQQQ"); plotter->Fill2D("QQQCartesianPlot" + std::to_string(qqq.id[i]), 200, -100, 100, 200, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hCalQQQ");
plotter->Fill2D("PC_XY_Projection_QQQ" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hPCQQQ"); plotter->Fill2D("PC_XY_Projection_QQQ" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hPCQQQ");
@ -545,12 +567,13 @@ Bool_t TrackRecon::Process(Long64_t entry)
for (int k = 0; k < pc.multi; k++) for (int k = 0; k < pc.multi; k++)
{ {
#ifdef RAW_HISTOS
plotter->Fill2D("RingCh_vs_Anode_Index", 16 * 4, 0, 16 * 4, 24, 0, 24, chRing + qqq.id[i] * 16, pc.index[k], "hRawQQQ"); plotter->Fill2D("RingCh_vs_Anode_Index", 16 * 4, 0, 16 * 4, 24, 0, 24, chRing + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
plotter->Fill2D("WedgeCh_vs_Anode_Index", 16 * 4, 0, 16 * 4, 24, 0, 24, chWedge + qqq.id[i] * 16, pc.index[k], "hRawQQQ"); plotter->Fill2D("WedgeCh_vs_Anode_Index", 16 * 4, 0, 16 * 4, 24, 0, 24, chWedge + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
plotter->Fill2D("WedgeCh_vs_Anode_Index" + std::to_string(qqq.id[i]), 16 * 4, 0, 16 * 4, 24, 0, 24, chWedge + qqq.id[i] * 16, pc.index[k]); plotter->Fill2D("WedgeCh_vs_Anode_Index" + std::to_string(qqq.id[i]), 16 * 4, 0, 16 * 4, 24, 0, 24, chWedge + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
plotter->Fill2D("RingCh_vs_Cathode_Index", 16 * 4, 0, 16 * 4, 24, 24, 48, chRing + qqq.id[i] * 16, pc.index[k], "hRawQQQ"); plotter->Fill2D("RingCh_vs_Cathode_Index", 16 * 4, 0, 16 * 4, 24, 24, 48, chRing + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
plotter->Fill2D("WedgeCh_vs_Cathode_Index", 16 * 4, 0, 16 * 4, 24, 24, 48, chWedge + qqq.id[i] * 16, pc.index[k], "hRawQQQ"); plotter->Fill2D("WedgeCh_vs_Cathode_Index", 16 * 4, 0, 16 * 4, 24, 24, 48, chWedge + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
#endif
if (pc.index[k] < 24 && pc.e[k] > 10) if (pc.index[k] < 24 && pc.e[k] > 10)
{ {
plotter->Fill2D("Timing_Difference_QQQ_PC", 500, -2000, 2000, 16, 0, 16, tRing - static_cast<double>(pc.t[k]), chRing, "hTiming"); plotter->Fill2D("Timing_Difference_QQQ_PC", 500, -2000, 2000, 16, 0, 16, tRing - static_cast<double>(pc.t[k]), chRing, "hTiming");
@ -591,8 +614,9 @@ Bool_t TrackRecon::Process(Long64_t entry)
} // i loop end } // i loop end
PCQQQTimeCut = PCAQQQTimeCut && PCCQQQTimeCut; PCQQQTimeCut = PCAQQQTimeCut && PCCQQQTimeCut;
#ifdef RAW_HISTOS
plotter->Fill1D("QQQ_Multiplicity", 10, 0, 10, qqqCount, "hRawQQQ"); plotter->Fill1D("QQQ_Multiplicity", 10, 0, 10, qqqCount, "hRawQQQ");
#endif
typedef std::unordered_map<int, std::tuple<int, double, double>> WireEvent; // this stores nearest neighbour wire events, or a 'cluster' typedef std::unordered_map<int, std::tuple<int, double, double>> WireEvent; // this stores nearest neighbour wire events, or a 'cluster'
WireEvent aWireEvents, cWireEvents; // naming for book keeping WireEvent aWireEvents, cWireEvents; // naming for book keeping
@ -608,14 +632,12 @@ Bool_t TrackRecon::Process(Long64_t entry)
for (int i = 0; i < pc.multi; i++) for (int i = 0; i < pc.multi; i++)
{ {
// std::cout << pc.index[i] << " " << pc.e[i] << " " << std::endl; // std::cout << pc.index[i] << " " << pc.e[i] << " " << std::endl;
#ifdef RAW_HISTOS
if (pc.e[i] > 10) if (pc.e[i] > 10)
{ {
plotter->Fill2D("PC_Index_Vs_Energy", 48, 0, 48, 2000, 0, 30000, pc.index[i], static_cast<double>(pc.e[i]), "hRawPC"); plotter->Fill2D("PC_Index_Vs_Energy", 48, 0, 48, 2000, 0, 30000, pc.index[i], static_cast<double>(pc.e[i]), "hRawPC");
} }
else #endif
{
continue;
}
if (pc.index[i] < 48) if (pc.index[i] < 48)
{ {
@ -667,8 +689,10 @@ Bool_t TrackRecon::Process(Long64_t entry)
for (int j = i + 1; j < pc.multi; j++) for (int j = i + 1; j < pc.multi; j++)
{ {
#ifdef RAW_HISTOS
plotter->Fill2D("PC_Coincidence_Matrix", 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC"); plotter->Fill2D("PC_Coincidence_Matrix", 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC");
plotter->Fill2D("PC_Coincidence_Matrix_anodeMinusCathode_lt_-200_" + std::to_string(anodeT - cathodeT < -200), 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC"); plotter->Fill2D("PC_Coincidence_Matrix_anodeMinusCathode_lt_-200_" + std::to_string(anodeT - cathodeT < -200), 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC");
#endif
plotter->Fill2D("Anode_V_Anode", 24, 0, 24, 24, 0, 24, pc.index[i], pc.index[j], "hGMPC"); plotter->Fill2D("Anode_V_Anode", 24, 0, 24, 24, 0, 24, pc.index[i], pc.index[j], "hGMPC");
} }
} }
@ -748,7 +772,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
protonAlphaHistograms(plotter, QQQ_Events, SX3_Events, PC_Events); protonAlphaHistograms(plotter, QQQ_Events, SX3_Events, PC_Events);
// return kTRUE; // return kTRUE;
} // end if(process_alpha_proton_scattering) } // end if(process_alpha_proton_scattering)
#ifdef RAW_HISTOS
if (QQQ_Events.size() && PC_Events.size()) if (QQQ_Events.size() && PC_Events.size())
plotter->Fill2D("PCEv_vs_QQQEv", 20, 0, 20, 20, 0, 20, QQQ_Events.size(), PC_Events.size()); plotter->Fill2D("PCEv_vs_QQQEv", 20, 0, 20, 20, 0, 20, QQQ_Events.size(), PC_Events.size());
@ -766,6 +790,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
{ {
plotter->Fill2D("ac_vs_cc_ign0", 20, 0, 20, 20, 0, 20, aClusters.size(), cClusters.size(), "wiremult"); plotter->Fill2D("ac_vs_cc_ign0", 20, 0, 20, 20, 0, 20, aClusters.size(), cClusters.size(), "wiremult");
} }
#endif
for (auto sx3event : SX3_Events) for (auto sx3event : SX3_Events)
{ {
@ -850,282 +875,195 @@ Bool_t TrackRecon::Process(Long64_t entry)
{ {
PCQQQClusterAnalysis(plotter, QQQ_Events, SX3_Events, PC_Events); PCQQQClusterAnalysis(plotter, QQQ_Events, SX3_Events, PC_Events);
} }
return kTRUE;
///////////////////Single wire analysis for the anodes/////////////////// ///////////////////nA analysis using pseudo-wire (GetPseudoWire + getClosestWirePosAtWirePhi)///////////////////
if (aClusters.size() > 0)
if (aClusters.size() == 1 && cClusters.size() == 0)
{ {
// Extract the primary anode hit properties // ---------------------------------------------------------
auto anodeHit = aClusters.front().front(); // PROTON LOOP (SX3 BARREL)
int aWireID = std::get<0>(anodeHit); // ---------------------------------------------------------
double aEnergy = std::get<1>(anodeHit);
double aTime = std::get<2>(anodeHit);
// Get the 3D endpoints of the fired twisted anode wire from your geometry class
TVector3 a1 = pwinstance.An[aWireID].first;
TVector3 wireVec = pwinstance.An[aWireID].first - pwinstance.An[aWireID].second;
// Loop over SX3_Events directly
for (auto sx3event : SX3_Events) for (auto sx3event : SX3_Events)
{ {
// Pick the anode cluster closest in phi to this SX3 hit
const std::vector<std::tuple<int, double, double>> *bestCluster = &aClusters[0];
double bestDphi = 9999.0;
if (sx3event.Time1 - aTime < -150) // Time cut for protons for (const auto &acluster : aClusters)
{ {
// 1. Define the plane of the track (Z-axis to SX3 hit) auto [pw, sumE, maxE, tsMax] = pwinstance.GetPseudoWire(acluster, "ANODE");
TVector3 planeNormal(-TMath::Sin(sx3event.pos.Phi()), TMath::Cos(sx3event.pos.Phi()), 0.0); TVector3 pos = pwinstance.getClosestWirePosAtWirePhi(pw, sx3event.pos.Phi());
double dphi = TMath::Abs(TVector2::Phi_mpi_pi(sx3event.pos.Phi() - pos.Phi()));
// 2. Find intersection of the twisted wire with this track plane if (dphi < bestDphi)
double dot_wireVec = wireVec.Dot(planeNormal);
// Prevent divide-by-zero if wire is perfectly parallel to the track plane
if (TMath::Abs(dot_wireVec) < 1e-6)
continue;
double t_intersect = -(a1.Dot(planeNormal)) / dot_wireVec;
// Calculate the exact 3D coordinate on the wire that matches the SX3 phi
TVector3 pcz_intersect = a1 + t_intersect * wireVec;
// 3. Reconstruct Vertex Z
double deltaRho = sx3event.pos.Perp() - pcz_intersect.Perp();
double deltaZ = sx3event.pos.Z() - pcz_intersect.Z();
double vertex_recon = sx3event.pos.Z() - sx3event.pos.Perp() * (deltaZ / deltaRho);
double z_entrance = -274.3;
double beam_path_length = TMath::Abs(vertex_recon - z_entrance) * 0.1;
double initial_beam_energy = 72.0;
double beam_energy_at_vertex = cm_to_MeV_beam->Eval(
MeV_to_cm_beam->Eval(initial_beam_energy) - beam_path_length);
Kinematics apkin_p(26.981538, 4.002603, 1.007825, 29.973770, beam_energy_at_vertex);
Kinematics apkin_a(26.981538, 4.002603, 4.002603, 26.981538, beam_energy_at_vertex);
std::string vtx_gate = "";
if (vertex_recon >= -176.0 && vertex_recon < -100.0)
{ {
vtx_gate = "_Z[-176_to_-100]"; bestDphi = dphi;
} bestCluster = &acluster;
else if (vertex_recon >= -100.0 && vertex_recon < -50.0)
{
vtx_gate = "_Z[-100_to_-50]";
}
else if (vertex_recon >= -50.0 && vertex_recon < 0.0)
{
vtx_gate = "_Z[-50_to_0]";
}
else if (vertex_recon >= 0.0 && vertex_recon < 50.0)
{
vtx_gate = "_Z[0_to_50]";
}
else if (vertex_recon >= 50.0 && vertex_recon < 100.0)
{
vtx_gate = "_Z[50_to_100]";
}
else if (vertex_recon >= 100.0 && vertex_recon < 176.0)
{
vtx_gate = "_Z[100_to_176]";
}
// 4. Energy Loss Correction in Silicon
double path_length = (sx3event.pos - TVector3(0, 0, vertex_recon)).Mag() * 0.1;
double sx3Efix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(sx3event.Energy1) - path_length);
double sx3Efixalpha = cm_to_MeV->Eval(MeV_to_cm->Eval(sx3event.Energy1) - path_length);
double theta_recon = (sx3event.pos - TVector3(0, 0, vertex_recon)).Theta();
double sinTheta = TMath::Sin(theta_recon);
// Now these functions will use the correct, event-specific beam energy!
double Ex_from_proton = apkin_p.getExc(sx3Efix, theta_recon * 180. / M_PI);
double Ex_from_alpha = apkin_a.getExc(sx3Efixalpha, theta_recon * 180. / M_PI);
// 5. Fill Diagnostics
plotter->Fill2D("1A0C_dE_Ecorr_Anode_SX3", 400, 0, 30, 800, 0, 40000, sx3Efix, aEnergy * sinTheta, "1A0C");
plotter->Fill1D("1A0C_Ex_from_protons_SX3", 200, -10, 10, Ex_from_proton, "1A0C");
plotter->Fill1D("1A0C_Ex_from_alphas_SX3", 200, -10, 10, Ex_from_alpha, "1A0C");
plotter->Fill2D("1A0C_sx3_E_vs_theta_raw_SX3", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, sx3event.Energy1, "1A0C");
plotter->Fill2D("1A0C_sx3_E_vs_theta_corr_SX3", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, sx3Efix, "1A0C");
if (vtx_gate != "")
{
plotter->Fill1D("1A0C_twisted_vertex_recon_SX3" + vtx_gate, 600, -300, 300, vertex_recon, "1A0C");
plotter->Fill1D("1A0C_twisted_pcz_recon_SX3" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), "1A0C");
plotter->Fill2D("1A0C_dE_Ecorr_Anode_SX3" + vtx_gate, 400, 0, 30, 800, 0, 40000, sx3Efix, aEnergy * sinTheta, "1A0C");
plotter->Fill1D("1A0C_Ex_from_protons_SX3" + vtx_gate, 200, -10, 10, Ex_from_proton, "1A0C");
plotter->Fill1D("1A0C_Ex_from_alphas_SX3" + vtx_gate, 200, -10, 10, Ex_from_alpha, "1A0C");
// Track where on the wire the hit occurred (0 to 1 is inside the physical PC)
// plotter->Fill1D("1A0C_wire_t_parameter" + vtx_gate, 200, -0.5, 1.5, t_intersect, "1A0C");
} }
} }
// Extract the virtual wire specifically for the best cluster
auto [apwire, apSumE, apMaxE, apTSMaxE] = pwinstance.GetPseudoWire(*bestCluster, "ANODE");
std::string nA_label = std::to_string(bestCluster->size()) + "A";
TVector3 pcz_intersect = pwinstance.getClosestWirePosAtWirePhi(apwire, sx3event.pos.Phi());
double deltaRho = sx3event.pos.Perp() - pcz_intersect.Perp();
double deltaZ = sx3event.pos.Z() - pcz_intersect.Z();
double vertex_recon = sx3event.pos.Z() - sx3event.pos.Perp() * (deltaZ / deltaRho);
std::string vtx_gate = "";
#ifdef VTX_GATE
if (vertex_recon >= -176.0 && vertex_recon < -100.0)
vtx_gate = "_Z[-176_to_-100]";
else if (vertex_recon >= -100.0 && vertex_recon < -50.0)
vtx_gate = "_Z[-100_to_-50]";
else if (vertex_recon >= -50.0 && vertex_recon < 0.0)
vtx_gate = "_Z[-50_to_0]";
else if (vertex_recon >= 0.0 && vertex_recon < 50.0)
vtx_gate = "_Z[0_to_50]";
else if (vertex_recon >= 50.0 && vertex_recon < 100.0)
vtx_gate = "_Z[50_to_100]";
else if (vertex_recon >= 100.0 && vertex_recon < 176.0)
vtx_gate = "_Z[100_to_176]";
#endif
// *0.1 converts mm to cm for the Eloss
double beam_path_length = TMath::Abs(vertex_recon - z_entrance) * 0.1;
#ifdef AL_BEAM
double beam_energy_at_vertex = cm_to_MeV_27Al->Eval(MeV_to_cm_27Al->Eval(72.0) - beam_path_length);
Kinematics apkin_p(mass_27Al, mass_4He, mass_1H, mass_30Si, beam_energy_at_vertex / mass_27Al);
Kinematics apkin_a(mass_27Al, mass_4He, mass_4He, mass_30Si, beam_energy_at_vertex / mass_27Al);
#endif
#ifdef F_BEAM
double beam_energy_at_vertex = cm_to_MeV_17F->Eval(MeV_to_cm_17F->Eval(65.0) - beam_path_length);
Kinematics apkin_p(mass_17F, mass_4He, mass_1H, mass_20Ne, beam_energy_at_vertex / mass_17F);
Kinematics apkin_a(mass_17F, mass_4He, mass_4He, mass_20Ne, beam_energy_at_vertex / mass_17F);
#endif
double path_length = (sx3event.pos - TVector3(0, 0, vertex_recon)).Mag() * 0.1;
double sx3Efix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(sx3event.Energy1) - path_length);
double sx3Efixalpha = cm_to_MeV->Eval(MeV_to_cm->Eval(sx3event.Energy1) - path_length);
double theta_recon = (sx3event.pos - TVector3(0, 0, vertex_recon)).Theta();
double sinTheta = TMath::Sin(theta_recon);
double Ex_from_proton = apkin_p.getExc(sx3Efix, theta_recon * 180. / M_PI);
double Ex_from_alpha = apkin_a.getExc(sx3Efixalpha, theta_recon * 180. / M_PI);
// Fill standard un-gated plots
plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_SX3", 800, 0, 30, 800, 0, 30000, sx3Efix, apSumE * sinTheta, nA_label);
plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_SX3_alpha", 800, 0, 30, 800, 0, 30000, sx3Efixalpha, apSumE * sinTheta, nA_label);
plotter->Fill2D(nA_label + "_sx3_E_vs_theta_raw_SX3", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, sx3event.Energy1, nA_label);
plotter->Fill2D(nA_label + "_sx3_E_vs_theta_corr_SX3", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, sx3Efix, nA_label);
plotter->Fill1D(nA_label + "_Ex_from_protons_SX3", 1200, -30, 30, Ex_from_proton, nA_label);
plotter->Fill1D(nA_label + "_Ex_from_alphas_SX3", 1200, -30, 30, Ex_from_alpha, nA_label);
// Fill Gated Plots
// if (vtx_gate != "")
// {
// plotter->Fill2D("dE_Ecorr_Anode_SX3" + vtx_gate, 400, 0, 30, 800, 0, 40000, sx3Efix, apSumE * sinTheta, nA_label);
// plotter->Fill2D("dE_Ecorr_Anode_SX3_alpha" + vtx_gate, 400, 0, 30, 800, 0, 40000, sx3Efixalpha, apSumE * sinTheta, nA_label);
// plotter->Fill1D("twisted_pcz_recon_SX3" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA_label);
// plotter->Fill1D("twisted_vertex_recon_SX3" + vtx_gate, 600, -300, 300, vertex_recon, nA_label);
// plotter->Fill1D("Ex_from_protons_SX3" + vtx_gate, 1200, -30, 30, Ex_from_proton, nA_label);
// plotter->Fill1D("Ex_from_alphas_SX3" + vtx_gate, 1200, -30, 30, Ex_from_alpha, nA_label);
// }
} }
// Loop over QQQ_Events directly // ---------------------------------------------------------
// PROTON LOOP (QQQ ENDCAP)
// ---------------------------------------------------------
for (auto qqqevent : QQQ_Events) for (auto qqqevent : QQQ_Events)
{ {
if (qqqevent.Time1 - aTime < -150) // Time cut for protons const std::vector<std::tuple<int, double, double>> *bestCluster = nullptr;
double bestDphi = 9999.0;
for (const auto &acluster : aClusters)
{ {
// 1. Define the plane of the track (Z-axis to SX3 hit) auto [apw, sumE, maxE, tsMax] = pwinstance.GetPseudoWire(acluster, "ANODE");
TVector3 planeNormal(-TMath::Sin(qqqevent.pos.Phi()), TMath::Cos(qqqevent.pos.Phi()), 0.0); TVector3 pcPos = pwinstance.getClosestWirePosAtWirePhi(apw, qqqevent.pos.Phi());
double dphi = TMath::Abs(TVector2::Phi_mpi_pi(qqqevent.pos.Phi() - pcPos.Phi()));
// 2. Find intersection of the twisted wire with this track plane if (dphi < bestDphi)
double dot_wireVec = wireVec.Dot(planeNormal);
// Prevent divide-by-zero if wire is perfectly parallel to the track plane
if (TMath::Abs(dot_wireVec) < 1e-6)
continue;
double t_intersect_QQQ = -(a1.Dot(planeNormal)) / dot_wireVec;
// Calculate the exact 3D coordinate on the wire that matches the SX3 phi
TVector3 pcz_intersect = a1 + t_intersect_QQQ * wireVec;
// 3. Reconstruct Vertex Z
double deltaRho = qqqevent.pos.Perp() - pcz_intersect.Perp();
double deltaZ = qqqevent.pos.Z() - pcz_intersect.Z();
double vertex_recon = qqqevent.pos.Z() - qqqevent.pos.Perp() * (deltaZ / deltaRho);
double z_entrance = -274.3;
double beam_path_length = TMath::Abs(vertex_recon - z_entrance) * 0.1;
double initial_beam_energy = 72.0;
double beam_energy_at_vertex = cm_to_MeV_beam->Eval(
MeV_to_cm_beam->Eval(initial_beam_energy) - beam_path_length);
Kinematics apkin_p(26.981538, 4.002603, 1.007825, 29.973770, beam_energy_at_vertex);
Kinematics apkin_a(26.981538, 4.002603, 4.002603, 26.981538, beam_energy_at_vertex);
// ==========================================================
std::string vtx_gate = "";
if (vertex_recon >= -176.0 && vertex_recon < -100.0)
{ {
vtx_gate = "_Z[-176_to_-100]"; bestDphi = dphi;
} bestCluster = &acluster;
else if (vertex_recon >= -100.0 && vertex_recon < -50.0)
{
vtx_gate = "_Z[-100_to_-50]";
}
else if (vertex_recon >= -50.0 && vertex_recon < 0.0)
{
vtx_gate = "_Z[-50_to_0]";
}
else if (vertex_recon >= 0.0 && vertex_recon < 50.0)
{
vtx_gate = "_Z[0_to_50]";
}
else if (vertex_recon >= 50.0 && vertex_recon < 100.0)
{
vtx_gate = "_Z[50_to_100]";
}
else if (vertex_recon >= 100.0 && vertex_recon < 176.0)
{
vtx_gate = "_Z[100_to_176]";
}
// 4. Energy Loss Correction in Silicon
double path_length = (qqqevent.pos - TVector3(0, 0, vertex_recon)).Mag() * 0.1;
double qqqEfix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(qqqevent.Energy1) - path_length);
double qqqEfixalpha = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy2) - path_length);
double theta_recon = (qqqevent.pos - TVector3(0, 0, vertex_recon)).Theta();
double sinTheta = TMath::Sin(theta_recon);
double Ex_from_proton = apkin_p.getExc(qqqEfix, theta_recon * 180. / M_PI);
double Ex_from_alpha = apkin_a.getExc(qqqEfixalpha, theta_recon * 180. / M_PI);
// 5. Fill Diagnostics
// plotter->Fill2D("1A0C_dE_Ecorr_Anode_QQQ", 400, 0, 30, 800, 0, 40000, qqqEfix, aEnergy * sinTheta, "1A0C");
// plotter->Fill1D("1A0C_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, "1A0C");
// plotter->Fill1D("1A0C_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, "1A0C");
// plotter->Fill2D("1A0C_qqq_E_vs_theta_raw_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqevent.Energy1, "1A0C");
// plotter->Fill2D("1A0C_qqq_E_vs_theta_corr_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqEfix, "1A0C");
if (vtx_gate != "")
{
plotter->Fill1D("1A0C_twisted_pcz_recon_QQQ" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), "1A0C");
plotter->Fill1D("1A0C_twisted_vertex_recon_QQQ" + vtx_gate, 600, -300, 300, vertex_recon, "1A0C");
plotter->Fill2D("1A0C_dE_Ecorr_Anode_QQQ" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfix, aEnergy * sinTheta, "1A0C");
plotter->Fill2D("1A0C_dE_Ecorr_Anode_QQQ_alpha" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfixalpha, aEnergy * sinTheta, "1A0C");
plotter->Fill1D("1A0C_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, "1A0C");
plotter->Fill1D("1A0C_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, "1A0C");
// Track where on the wire the hit occurred (0 to 1 is inside the physical PC)
plotter->Fill1D("1A0C_wire_t_parameter_QQQ" + vtx_gate, 200, -0.5, 1.5, t_intersect_QQQ, "1A0C");
} }
} }
if (!bestCluster)
continue;
///////////////////nA0C analysis using pseudo-wire (GetPseudoWire + getClosestWirePosAtWirePhi)/////////////////// // Extract the virtual wire specifically for the best cluster
if (cClusters.size() == 0 && aClusters.size() > 0) auto [apwire, apSumE, apMaxE, apTSMaxE] = pwinstance.GetPseudoWire(*bestCluster, "ANODE");
{ std::string nA_label = std::to_string(bestCluster->size()) + "A";
std::string nA0C_label = std::to_string(aClusters.size()) + "A0C";
// Flatten all anode clusters into a combined hit list for the pseudo-wire TVector3 pcz_intersect = pwinstance.getClosestWirePosAtWirePhi(apwire, qqqevent.pos.Phi());
std::vector<std::tuple<int, double, double>> allAnodeHits;
for (const auto &acluster : aClusters) double deltaRho = qqqevent.pos.Perp() - pcz_intersect.Perp();
for (const auto &hit : acluster) double deltaZ = qqqevent.pos.Z() - pcz_intersect.Z();
allAnodeHits.push_back(hit); double vertex_recon = qqqevent.pos.Z() - qqqevent.pos.Perp() * (deltaZ / deltaRho);
auto [apwire, apSumE, apMaxE, apTSMaxE] = pwinstance.GetPseudoWire(allAnodeHits, "ANODE");
for (auto qqqevent : QQQ_Events) std::string vtx_gate = "";
{
if (qqqevent.Time1 - apTSMaxE < -150) #ifdef VTX_GATE
{ if (vertex_recon >= -176.0 && vertex_recon < -100.0)
// Use pseudo-wire to find the PC hit position at the QQQ phi vtx_gate = "_Z[-176_to_-100]";
TVector3 pcz_intersect = pwinstance.getClosestWirePosAtWirePhi(apwire, qqqevent.pos.Phi()); else if (vertex_recon >= -100.0 && vertex_recon < -50.0)
// Reconstruct vertex Z vtx_gate = "_Z[-100_to_-50]";
double deltaRho = qqqevent.pos.Perp() - pcz_intersect.Perp(); else if (vertex_recon >= -50.0 && vertex_recon < 0.0)
double deltaZ = qqqevent.pos.Z() - pcz_intersect.Z(); vtx_gate = "_Z[-50_to_0]";
double vertex_recon = qqqevent.pos.Z() - qqqevent.pos.Perp() * (deltaZ / deltaRho); else if (vertex_recon >= 0.0 && vertex_recon < 50.0)
double z_entrance = -274.3; vtx_gate = "_Z[0_to_50]";
double beam_path_length = TMath::Abs(vertex_recon - z_entrance) * 0.1; else if (vertex_recon >= 50.0 && vertex_recon < 100.0)
double initial_beam_energy = 72.0; vtx_gate = "_Z[50_to_100]";
double beam_energy_at_vertex = cm_to_MeV_beam->Eval( else if (vertex_recon >= 100.0 && vertex_recon < 176.0)
MeV_to_cm_beam->Eval(initial_beam_energy) - beam_path_length); vtx_gate = "_Z[100_to_176]";
Kinematics apkin_p(26.981538, 4.002603, 1.007825, 29.973770, beam_energy_at_vertex); #endif
Kinematics apkin_a(26.981538, 4.002603, 4.002603, 26.981538, beam_energy_at_vertex);
std::string vtx_gate = ""; double beam_path_length = TMath::Abs(vertex_recon - z_entrance) * 0.1;
if (vertex_recon >= -176.0 && vertex_recon < -100.0)
vtx_gate = "_Z[-176_to_-100]"; #ifdef AL_BEAM
else if (vertex_recon >= -100.0 && vertex_recon < -50.0) double beam_energy_at_vertex = cm_to_MeV_27Al->Eval(MeV_to_cm_27Al->Eval(72.0) - beam_path_length);
vtx_gate = "_Z[-100_to_-50]"; Kinematics apkin_p(mass_27Al, mass_4He, mass_1H, mass_30Si, beam_energy_at_vertex / mass_27Al);
else if (vertex_recon >= -50.0 && vertex_recon < 0.0) Kinematics apkin_a(mass_27Al, mass_4He, mass_4He, mass_30Si, beam_energy_at_vertex / mass_27Al);
vtx_gate = "_Z[-50_to_0]"; #endif
else if (vertex_recon >= 0.0 && vertex_recon < 50.0)
vtx_gate = "_Z[0_to_50]"; #ifdef F_BEAM
else if (vertex_recon >= 50.0 && vertex_recon < 100.0) double beam_energy_at_vertex = cm_to_MeV_17F->Eval(MeV_to_cm_17F->Eval(65.0) - beam_path_length);
vtx_gate = "_Z[50_to_100]"; Kinematics apkin_p(mass_17F, mass_4He, mass_1H, mass_20Ne, beam_energy_at_vertex / mass_17F);
else if (vertex_recon >= 100.0 && vertex_recon < 176.0) Kinematics apkin_a(mass_17F, mass_4He, mass_4He, mass_20Ne, beam_energy_at_vertex / mass_17F);
vtx_gate = "_Z[100_to_176]"; #endif
double path_length = (qqqevent.pos - TVector3(0, 0, vertex_recon)).Mag() * 0.1;
double qqqEfix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(qqqevent.Energy1) - path_length); // Energy Loss Correction
double qqqEfixalpha = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy2) - path_length); double path_length = (qqqevent.pos - TVector3(0, 0, vertex_recon)).Mag() * 0.1;
double theta_recon = (qqqevent.pos - TVector3(0, 0, vertex_recon)).Theta(); double qqqEfix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(qqqevent.Energy1) - path_length);
double sinTheta = TMath::Sin(theta_recon); double qqqEfixalpha = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy2) - path_length);
double Ex_from_proton = apkin_p.getExc(qqqEfix, theta_recon * 180. / M_PI);
double Ex_from_alpha = apkin_a.getExc(qqqEfixalpha, theta_recon * 180. / M_PI); double theta_recon = (qqqevent.pos - TVector3(0, 0, vertex_recon)).Theta();
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_QQQ", 400, 0, 30, 800, 0, 40000, qqqEfix, apSumE * sinTheta, nA0C_label); double sinTheta = TMath::Sin(theta_recon);
plotter->Fill1D(nA0C_label + "_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label);
plotter->Fill1D(nA0C_label + "_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label); double Ex_from_proton = apkin_p.getExc(qqqEfix, theta_recon * 180. / M_PI);
plotter->Fill2D(nA0C_label + "_qqq_E_vs_theta_raw_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqevent.Energy1, nA0C_label); double Ex_from_alpha = apkin_a.getExc(qqqEfixalpha, theta_recon * 180. / M_PI);
plotter->Fill2D(nA0C_label + "_qqq_E_vs_theta_corr_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqEfix, nA0C_label);
if (vtx_gate != "") // Fill standard un-gated plots
{ plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_QQQ", 800, 0, 30, 800, 0, 30000, qqqEfix, apSumE * sinTheta, nA_label);
plotter->Fill1D(nA0C_label + "_twisted_pcz_recon_QQQ" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA0C_label); plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_QQQ_alpha", 800, 0, 30, 800, 0, 30000, qqqEfixalpha, apSumE * sinTheta, nA_label);
plotter->Fill1D(nA0C_label + "_twisted_vertex_recon_QQQ" + vtx_gate, 600, -300, 300, vertex_recon, nA0C_label); plotter->Fill2D(nA_label + "_qqq_E_vs_theta_raw_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqevent.Energy1, nA_label);
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_QQQ" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfix, apSumE * sinTheta, nA0C_label); plotter->Fill2D(nA_label + "_qqq_E_vs_theta_corr_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqEfix, nA_label);
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_QQQ_alpha" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfixalpha, apSumE * sinTheta, nA0C_label); plotter->Fill1D(nA_label + "_Ex_from_alphas_QQQ", 1200, -30, 30, Ex_from_alpha, nA_label);
plotter->Fill1D(nA0C_label + "_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label); plotter->Fill1D(nA_label + "_Ex_from_protons_QQQ", 1200, -30, 30, Ex_from_proton, nA_label);
plotter->Fill1D(nA0C_label + "_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label);
} // Fill Gated Plots
} // if (vtx_gate != "")
} // {
} // plotter->Fill1D(nA_label + "_twisted_pcz_recon_QQQ" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA_label);
// plotter->Fill1D(nA_label + "_twisted_vertex_recon_QQQ" + vtx_gate, 600, -300, 300, vertex_recon, nA_label);
// plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_QQQ" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfix, apSumE * sinTheta, nA_label);
// plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_QQQ_alpha" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfixalpha, apSumE * sinTheta, nA_label);
// plotter->Fill1D(nA_label + "_Ex_from_alphas_QQQ" + vtx_gate, 1200, -30, 30, Ex_from_alpha, nA_label);
// plotter->Fill1D(nA_label + "_Ex_from_protons_QQQ" + vtx_gate, 1200, -30, 30, Ex_from_proton, nA_label);
// }
} }
} }
@ -1290,6 +1228,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
// --------------------------------------------------------- // ---------------------------------------------------------
if (anodeHits.size() > 0 && cathodeHits.size() > 0) if (anodeHits.size() > 0 && cathodeHits.size() > 0)
{ {
#ifdef RAW_HISTOS
plotter->Fill2D("AHits_vs_CHits_NA" + std::to_string(hasNeighbourAnodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC"); plotter->Fill2D("AHits_vs_CHits_NA" + std::to_string(hasNeighbourAnodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
plotter->Fill2D("AHits_vs_CHits_NC" + std::to_string(hasNeighbourCathodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC"); plotter->Fill2D("AHits_vs_CHits_NC" + std::to_string(hasNeighbourCathodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
@ -1299,6 +1238,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
{ {
plotter->Fill2D("AHits_vs_CHits_NN", 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC"); plotter->Fill2D("AHits_vs_CHits_NN", 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
} }
#endif
} }
if (HitNonZero && anodeIntersection.Z() != 0) if (HitNonZero && anodeIntersection.Z() != 0)
@ -1912,4 +1852,35 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, std::vector<Event> QQQ_Events, s
} // end PC QQQ coincidence } // end PC QQQ coincidence
// HALFTIME! Can stop here in future versions // HALFTIME! Can stop here in future versions
// return kTRUE; // 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,
const std::string &targetName,
const std::string &anodeLabel,
Kinematics &kin,
double energyMeV,
double thetaRad)
{
// 1. Convert angle to degrees for your kinematics class
double thetaDeg = thetaRad * 180.0 / M_PI;
// 2. Calculate the Excitation Energy (Ex)
double Ex = kin.getExc(energyMeV, thetaDeg);
// 3. Dynamically build the histogram names
// Example: "27Al_Ex_a1" or "17F_Kinematics_aN"
std::string hNameEx = targetName + "_Ex_" + anodeLabel;
std::string hNameKin = targetName + "_Kinematics_" + anodeLabel;
// 4. Group them in a specific folder in the ROOT file
std::string folderName = "Excitation_" + targetName;
// 5. Fill the Histograms
// Fill 1D Excitation Energy plot
plotter->Fill1D(hNameEx, 400, -10, 20, Ex, folderName);
// Fill 2D Kinematics plot (Angle vs Energy)
plotter->Fill2D(hNameKin, 180, 0, 180, 500, 0, 25, thetaDeg, energyMeV, folderName);
} }

55
run_27Al.sh
View File

@ -1,32 +1,33 @@
rm results_run*.root #!/bin/bash
export DATASET="27Al" export DATASET="27Al"
export flip180="0"
export flipa=0
export flipc=0
export anode_offset=0
export cathode_offset=0
#declare -i run=28
#while [[ $run -lt 34 ]]; do #runs 1 to 84
# wrun=$(printf "%03d" $run)
# root -q -l -b -x ../ANASEN_analysis/data/27Al_Data/Run_"$wrun"_mapped.root -e 'tree->Process("Make#Vertex.C+O")'; mv Analyzer_SX3.root 27Al_output/results_run$wrun.root;
# run=run+1
#done
declare -i run=50 # Clean up previous runs
while [[ $run -lt 51 ]]; do #runs 1 to 84 rm -f 27Al_output/results_run*.root output_27Al.root
wrun=$(printf "%03d" $run)
root -q -l -b -x ../ANASEN_analysis/data/27Al_Data/Run_"$wrun"_mapped.root -e 'tree->Process("MakeVertex.C+O","Analyzer_27Al.root")'; mv Analyzer_27Al.root 27Al_output/results_run$wrun.root;
run=run+1
done
rm output.root echo "Pre-compiling TrackRecon.C safely on a single core..."
hadd -k -j 4 output.root 27Al_output/results_run*.root root -q -l -b -e '.L TrackRecon.C++O'
mv output.root output_27Al.root
#root -q -l -b -x -e '.L MakeVertex.C+O'; process_run() {
#halfproc=3 local wrun=$(printf "%03d" $1)
#parallel --ctag --bar -j $halfproc ./run.sh ::: {028..034} ::: 27Al #color-tag, linebuffer, then run run.sh in parallel local out="27Al_output/results_run${wrun}.root"
unset souce_vertex root -q -l -b -x "../ANASEN_analysis/data/27Al_Data/Run_${wrun}_mapped.root" \
unset DATASET -e "tree->Process(\"TrackRecon.C+\", \"${out}\")" > /dev/null 2>&1
unset flip180
if [ -f "$out" ]; then
echo "Run $wrun completed successfully."
else
echo "ERROR: Run $wrun failed to generate $out"
fi
}
export -f process_run
echo "Starting parallel processing..."
parallel --bar -j 4 process_run ::: {50..52}
echo "Merging files..."
hadd -k -j 4 output_27Al.root 27Al_output/results_run*.root
unset DATASET

8
run_tr.sh
View File

@ -27,7 +27,7 @@ export timecut_low=400.0;
#export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_010_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run10.root; #export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_010_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run10.root;
#export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_011_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run11.root; #export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_011_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run11.root;
export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_012_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run12.root; export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_012_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run12.root;
# exit exit
#export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_013_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run13.root; #export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_013_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run13.root;
#exit #exit
fi fi
@ -80,9 +80,9 @@ fi
#17F alpha run with gas #17F alpha run with gas
if [[ 1 -eq 1 ]]; then if [[ 1 -eq 1 ]]; then
export source_vertex=53.44; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_018_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run18.root; export source_vertex=53.44; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_018_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run18.root;
# export source_vertex=14.24; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_019_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run19.root; export source_vertex=14.24; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_019_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run19.root;
# export source_vertex=-24.96; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_020_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run20.root; export source_vertex=-24.96; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_020_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run20.root;
# export source_vertex=-73.96; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_021_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run21.root; export source_vertex=-73.96; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_021_mapped.root -e 'tree->Process("TrackRecon.C+O")'; mv Analyzer_SX3.root results_run21.root;
fi fi
#17F reaction data #17F reaction data
#export flip180="0" #export flip180="0"