diff --git a/MiscEdE.C b/MiscEdE.C new file mode 100644 index 0000000..26622d1 --- /dev/null +++ b/MiscEdE.C @@ -0,0 +1,359 @@ +#include "TFile.h" +#include "TTree.h" +#include "TChain.h" +#include "TH1D.h" +#include "TH2D.h" +#include "TCutG.h" +#include "TKey.h" +#include "TString.h" +#include "TSystem.h" +#include +#include +#include +#include + +namespace +{ + constexpr double kAdcMax = 16384; // 14-bit-ish ADC range for axes + + struct HitData + { + double e; + unsigned long long t; + }; +} + +void MiscEdE(const char *f0, + const char *f1 = "", const char *f2 = "", const char *f3 = "", + const char *f4 = "", const char *f5 = "", const char *f6 = "", + const char *f7 = "", const char *f8 = "", const char *f9 = "") +{ + std::vector inputs; + for (const char *f : {f0, f1, f2, f3, f4, f5, f6, f7, f8, f9}) + if (f && f[0]) + inputs.emplace_back(f); + if (inputs.empty()) + { + printf("MiscEdE: no input files given.\n"); + return; + } + + TChain *tree = new TChain("tree"); + for (auto &in : inputs) + { + if (gSystem->AccessPathName(in)) + { + printf("MiscEdE: WARNING input not found, skipping: %s\n", in.Data()); + continue; + } + tree->Add(in); + } + if (tree->GetNtrees() == 0) + { + printf("MiscEdE: no readable inputs.\n"); + return; + } + + // --------------------------------------------------------- + // Bash Environment Variable Extraction + // --------------------------------------------------------- + int runNumber = -1; + if (gSystem->Getenv("CURRENT_RUN")) { + runNumber = std::atoi(gSystem->Getenv("CURRENT_RUN")); + } + printf("MiscEdE: Run Number loaded from Bash: %d\n", runNumber); + + TString cutDir = "~/ANASEN_analysis"; // Fallback if not run via bash + if (gSystem->Getenv("CUT_DIR")) { + cutDir = gSystem->Getenv("CUT_DIR"); + } + + // --------------------------------------------------------- + // Run-Dependent Channel Routing + // --------------------------------------------------------- + // Defaults (Run <= 180) + int snLolli = 405, chLolli = 9; + int snSi = 405, chSi = 11; + int snRF = 405, chRF = 15; + int snMCP = 405, chMCP = 14; + + if (runNumber > 180 && runNumber <= 282) { + chRF = 10; + } + else if (runNumber > 282 && runNumber <= 322) { + snRF = 89; chRF = 0; + snMCP = 89; chMCP = 1; + } + else if (runNumber > 322) { + snRF = 89; chRF = 0; + snMCP = 89; chMCP = 2; + } + + // Composite Unique IDs (UID = Board * 100 + Channel) to prevent collisions + int uidLollipop = snLolli * 100 + chLolli; + int uidSiMon = snSi * 100 + chSi; + int uidRF = snRF * 100 + chRF; + int uidMCP = snMCP * 100 + chMCP; + + printf(" -> Routed Lollipop : Board %d, Ch %d (UID: %d)\n", snLolli, chLolli, uidLollipop); + printf(" -> Routed Si Monitor: Board %d, Ch %d (UID: %d)\n", snSi, chSi, uidSiMon); + printf(" -> Routed RF : Board %d, Ch %d (UID: %d)\n", snRF, chRF, uidRF); + printf(" -> Routed MCP : Board %d, Ch %d (UID: %d)\n", snMCP, chMCP, uidMCP); + + // --------------------------------------------------------- + // Exact schema matching EventBuilder.cpp + // --------------------------------------------------------- + const int MAX_MULTI = 2000; + unsigned int multi = 0; + unsigned short sn[MAX_MULTI]; + unsigned short ch[MAX_MULTI]; + unsigned short e[MAX_MULTI]; + unsigned long long e_t[MAX_MULTI]; // Timestamp branch is "e_t" + + tree->SetBranchAddress("multi", &multi); + tree->SetBranchAddress("sn", sn); + tree->SetBranchAddress("ch", ch); + tree->SetBranchAddress("e", e); + tree->SetBranchAddress("e_t", e_t); + + TString stem = gSystem->BaseName(inputs[0].Data()); + if (stem.EndsWith(".root")) + stem.Remove(stem.Length() - 5); + TString outName = "MiscEdE_" + stem + ".root"; + TFile *out = new TFile(outName, "recreate"); + + // --------------------------------------------------------- + // Load the 17F and 16O TCutGs from file + // --------------------------------------------------------- + TCutG *cut17F = nullptr; + TString path17F = gSystem->ExpandPathName(Form("%s/17FCut.root", cutDir.Data())); + if (!gSystem->AccessPathName(path17F)) + { + TFile *fCut = TFile::Open(path17F); + if (fCut && !fCut->IsZombie()) + { + for (auto keyObj : *fCut->GetListOfKeys()) + { + TKey *key = (TKey *)keyObj; + if (TString(key->GetClassName()) == "TCutG") + { + cut17F = (TCutG *)key->ReadObj()->Clone("cut17F"); + break; + } + } + fCut->Close(); + } + } + + TCutG *cut16O = nullptr; + TString path16O = gSystem->ExpandPathName(Form("%s/16OCut.root", cutDir.Data())); + if (!gSystem->AccessPathName(path16O)) + { + TFile *fCut = TFile::Open(path16O); + if (fCut && !fCut->IsZombie()) + { + for (auto keyObj : *fCut->GetListOfKeys()) + { + TKey *key = (TKey *)keyObj; + if (TString(key->GetClassName()) == "TCutG") + { + cut16O = (TCutG *)key->ReadObj()->Clone("cut16O"); + break; + } + } + fCut->Close(); + } + } + + out->cd(); + + // --------------------------------------------------------- + // Prepare Histograms + // --------------------------------------------------------- + std::map hEvsChMap; + auto getEvsCh = [&](int boardSN) -> TH2D* { + if (hEvsChMap.count(boardSN)) return hEvsChMap[boardSN]; + TH2D *h = new TH2D(Form("h2_E_vs_ch_bd%d", boardSN), Form("Board %d: E vs channel;board channel;E [ADC]", boardSN), 16, -0.5, 15.5, 800, 0, kAdcMax); + hEvsChMap[boardSN] = h; + return h; + }; + + TH2D *h2_EvsdT_LolliSi = new TH2D("EvsdT_Lollipop_Si", "Energy Si vs dT (Lollipop - Si);dT [ticks];Si E [ADC]", 500, 0, 2000, 400, 0, kAdcMax); + TH2D *h2_TOF_SiRF = new TH2D("TOF_Si_RF", "Time of Flight: Energy Si vs (T_Si - T_RF);T_Si - T_RF [ticks];Si E [ADC]", 500, 0, 2000, 400, 0, kAdcMax); + + // Dynamic Pairwise E-dE mapped by UID + std::map, TH2D *> hEdE; + auto getEdE = [&](int uidA, int uidB) -> TH2D * { + if (uidA > uidB) std::swap(uidA, uidB); + auto key = std::make_pair(uidA, uidB); + auto it = hEdE.find(key); + if (it != hEdE.end()) return it->second; + + TString name = Form("h2_ede_bd%d_ch%d_vs_bd%d_ch%d", uidA/100, uidA%100, uidB/100, uidB%100); + TString title = Form("E-dE: Bd %d Ch %d (x) vs Bd %d Ch %d (y);X E [ADC];Y E [ADC]", uidA/100, uidA%100, uidB/100, uidB%100); + + TH2D *h = new TH2D(name, title, 400, 0, kAdcMax, 400, 0, kAdcMax); + hEdE[key] = h; + return h; + }; + + // 17F Gated Plots + TH2D *hEdE_17F = nullptr; + TH2D *hEdE_inv = nullptr; + TH2D *h2_EvsdT_LolliSi_17F = nullptr; + TH2D *h2_TOF_SiRF_17F = nullptr; + TH1D *h1_dT_17F = nullptr; + + if (cut17F) + { + hEdE_17F = new TH2D("EdE_lollipopIC_vs_SiMonitor_17FCut", "lollipop IC vs Si monitor (17F Gated);Si monitor E [ADC];lollipop IC E [ADC]", 400, 0, kAdcMax, 400, 0, kAdcMax); + hEdE_inv = new TH2D("EdE_lollipopIC_vs_SiMonitor_17FCut_Inverse", "lollipop IC vs Si monitor (Inverse 17F Gated);Si monitor E [ADC];lollipop IC E [ADC]", 400, 0, kAdcMax, 400, 0, kAdcMax); + h2_EvsdT_LolliSi_17F = new TH2D("EvsdT_Lollipop_Si_17FCut", "Energy Si vs dT (17F Gated);dT (Lollipop - Si) [ticks];Si E [ADC]", 500, 0, 2000, 400, 0, kAdcMax); + h2_TOF_SiRF_17F = new TH2D("TOF_Si_RF_17FCut", "Si TOF vs Energy (17F Gated);T_Si - T_RF [ticks];Si E [ADC]", 500, 0, 2000, 400, 0, kAdcMax); + h1_dT_17F = new TH1D("dT_Lollipop_Si_17FCut", "dT Lollipop-Si (17F Gated);T_Lollipop - T_Si [ticks]", 500, 0, 2000); + } + + // 16O Gated Plots + TH2D *hEdE_16O = nullptr; + TH2D *h2_EvsdT_LolliSi_16O = nullptr; + TH2D *h2_TOF_SiRF_16O = nullptr; + TH1D *h1_dT_16O = nullptr; + + if (cut16O) + { + hEdE_16O = new TH2D("EdE_lollipopIC_vs_SiMonitor_16OCut", "lollipop IC vs Si monitor (16O Gated);Si monitor E [ADC];lollipop IC E [ADC]", 400, 0, kAdcMax, 400, 0, kAdcMax); + h2_EvsdT_LolliSi_16O = new TH2D("EvsdT_Lollipop_Si_16OCut", "Energy Si vs dT (16O Gated);dT (Lollipop - Si) [ticks];Si E [ADC]", 500, 0, 2000, 400, 0, kAdcMax); + h2_TOF_SiRF_16O = new TH2D("TOF_Si_RF_16OCut", "Si TOF vs Energy (16O Gated);T_Si - T_RF [ticks];Si E [ADC]", 500, 0, 2000, 800, 0, kAdcMax); + h1_dT_16O = new TH1D("dT_Lollipop_Si_16OCut", "dT Lollipop-Si (16O Gated);T_Lollipop - T_Si [ticks]", 500, 0, 2000); + } + + Long64_t nEnt = tree->GetEntries(); + printf("MiscEdE: %lld entries across %d file(s) -> %s\n", nEnt, tree->GetNtrees(), outName.Data()); + + Long64_t used = 0; + for (Long64_t i = 0; i < nEnt; ++i) + { + tree->GetEntry(i); + + // We map hits by UID instead of standard channel to prevent overlaps between Bd405 and Bd89 + std::map hit; + for (unsigned int j = 0; j < multi && j < (unsigned)MAX_MULTI; ++j) + { + // Only care about MISC boards 405 and 89 + if (sn[j] != 405 && sn[j] != 89) + continue; + + int c = ch[j]; + int uid = (sn[j] * 100) + ch[j]; + double en = e[j]; + unsigned long long ts = e_t[j]; + + auto it = hit.find(uid); + if (it == hit.end() || en > it->second.e) + hit[uid] = {en, ts}; + + getEvsCh(sn[j])->Fill(c, en); + } + + if (hit.empty()) + continue; + ++used; + + // --------------------------------------------------------- + // Generic Dynamic Pairwise E-dE Logic across all UIDs + // --------------------------------------------------------- + for (auto a = hit.begin(); a != hit.end(); ++a) + { + for (auto b = std::next(a); b != hit.end(); ++b) + { + if (a->second.e > 500 && b->second.e > 500) + { + getEdE(a->first, b->first)->Fill(a->second.e, b->second.e); + } + } + } + + // --------------------------------------------------------- + // Core Timing and PID Gating Logic (Lollipop-Si specific) + // --------------------------------------------------------- + if (hit.count(uidLollipop) && hit.count(uidSiMon)) + { + double e_lol = hit[uidLollipop].e; + double e_si = hit[uidSiMon].e; + + // Calculate delta T safely across potential CAEN rollovers + double dt_lol_si = static_cast(hit[uidLollipop].t) - static_cast(hit[uidSiMon].t); + + if (e_lol > 500 && e_si > 500) + { + h2_EvsdT_LolliSi->Fill(dt_lol_si, e_si); + + bool is17F = cut17F && cut17F->IsInside(e_si, e_lol); + bool is16O = cut16O && cut16O->IsInside(e_si, e_lol); + + if (is17F) + { + hEdE_17F->Fill(e_si, e_lol); + h2_EvsdT_LolliSi_17F->Fill(dt_lol_si, e_si); + h1_dT_17F->Fill(dt_lol_si); + } + else if (cut17F) + { + hEdE_inv->Fill(e_si, e_lol); + } + + if (is16O) + { + hEdE_16O->Fill(e_si, e_lol); + h2_EvsdT_LolliSi_16O->Fill(dt_lol_si, e_si); + h1_dT_16O->Fill(dt_lol_si); + } + + // Generate absolute TOF if the RF channel fired + if (hit.count(uidRF)) + { + double tof = static_cast(hit[uidSiMon].t) - static_cast(hit[uidRF].t); + h2_TOF_SiRF->Fill(tof, e_si); + + if (is17F) + h2_TOF_SiRF_17F->Fill(tof, e_si); + if (is16O) + h2_TOF_SiRF_16O->Fill(tof, e_si); + } + } + } + } + + // --------------------------------------------------------- + // Purge Sparse E-dE Histograms (< 500 counts) + // --------------------------------------------------------- + for (auto it = hEdE.begin(); it != hEdE.end();) + { + if (it->second->GetEntries() < 500) + { + delete it->second; // Free memory and remove from ROOT file + it = hEdE.erase(it); + } + else + { + ++it; + } + } + + // Convenience alias: Transpose generic plot to correctly place stopping E on X + if (hEdE.count({uidLollipop, uidSiMon})) + { + TH2D *src = hEdE[{uidLollipop, uidSiMon}]; + TH2D *alias = new TH2D("EdE_lollipopIC_vs_SiMonitor_Ungated", + Form("lollipop IC vs Si monitor (Ungated);Bd %d Ch %d (Si) E [ADC];Bd %d Ch %d (Lolli) E [ADC]", snSi, chSi, snLolli, chLolli), + src->GetNbinsY(), src->GetYaxis()->GetXmin(), src->GetYaxis()->GetXmax(), + src->GetNbinsX(), src->GetXaxis()->GetXmin(), src->GetXaxis()->GetXmax()); + for (int ix = 1; ix <= src->GetNbinsX(); ++ix) + for (int iy = 1; iy <= src->GetNbinsY(); ++iy) + alias->SetBinContent(iy, ix, src->GetBinContent(ix, iy)); // (x,y) -> (y,x) + } + + out->Write(); + out->Close(); + printf("MiscEdE: Processed %lld valid MISC events. Wrote %s\n", used, outName.Data()); +} \ No newline at end of file diff --git a/TrackRecon.C b/TrackRecon.C index 14adcc4..89bf59e 100644 --- a/TrackRecon.C +++ b/TrackRecon.C @@ -3951,6 +3951,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector &Si_ double ebeam_kin_3774keV = invertBeamEnergyMeV(m_beam, mass_4He, m3, m4, Efix, theta * 180 / M_PI, 3.774); double ebeam_kin_4809keV = invertBeamEnergyMeV(m_beam, mass_4He, m3, m4, Efix, theta * 180 / M_PI, 4.809); double ebeam_kin_5614keV = invertBeamEnergyMeV(m_beam, mass_4He, m3, m4, Efix, theta * 180 / M_PI, 5.164); + double ebeam_kin_6550keV = invertBeamEnergyMeV(m_beam, mass_4He, m3, m4, Efix, theta * 180 / M_PI, 6.550); // Gated output: only fill when this hypothesis (proton "_p" / alpha "_a") agrees // with which side of the proton_locus gate the event fell on, so each event @@ -3962,28 +3963,38 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector &Si_ auto plot_with_tag = [&](const std::string &topo) { std::string t = topo.empty() ? "" : ("_" + topo); - plotter->Fill1D(rx + "_Ex_from" + ejtag + t + sfx, 600, -15, 15, Ex, pmlabel); + plotter->Fill1D(rx + "_Ex_from" + ejtag + t + sfx, 600, -10, 20, Ex, pmlabel); plotter->Fill2D(rx + "_VertexReconZ_vs_Ef" + ejtag + t + sfx, 800, -400, 400, 800, 0, ef_max, vertex_z, Efix, pmlabel); - plotter->Fill2D(rx + "_VertexReconZ_vs_Ex" + ejtag + t + sfx, 800, -400, 400, 800, -20, 20, vertex_z, Ex, pmlabel); - plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + t + sfx, 360, 0, 180, 800, -20, 20, theta * 180 / M_PI, Ex, pmlabel); + plotter->Fill2D(rx + "_VertexReconZ_vs_Ex" + ejtag + t + sfx, 800, -400, 400, 600, -10, 20, vertex_z, Ex, pmlabel); + plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + t + sfx, 360, 0, 180, 600, -10, 20, theta * 180 / M_PI, Ex, pmlabel); if (ebeam_kin_MeV > 0.0) plotter->Fill2D(rx + "_BeamEnergy_ETrack_vs_EKin" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5, beam_energy_at_vertex, ebeam_kin_MeV, pmlabel); if (ejtag == "_p") { - plotter->Fill2D(rx + "_ETrack_vs_EKinGS" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5, - beam_energy_at_vertex, ebeam_kin_GS, "ETrackvsKin_assumed"); - plotter->Fill2D(rx + "_ETrack_vs_EKin2235keV" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5, - beam_energy_at_vertex, ebeam_kin_2235keV, "ETrackvsKin_assumed"); - plotter->Fill2D(rx + "_ETrack_vs_EKin3498kev" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5, - beam_energy_at_vertex, ebeam_kin_3498keV, "ETrackvsKin_assumed"); - plotter->Fill2D(rx + "_ETrack_vs_EKin3774keV" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5, - beam_energy_at_vertex, ebeam_kin_3774keV, "ETrackvsKin_assumed"); - plotter->Fill2D(rx + "_ETrack_vs_EKin4809keV" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5, - beam_energy_at_vertex, ebeam_kin_4809keV, "ETrackvsKin_assumed"); - plotter->Fill2D(rx + "_ETrack_vs_EKin5614keV" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5, - beam_energy_at_vertex, ebeam_kin_5614keV, "ETrackvsKin_assumed"); + if (beam_energy_at_vertex < 4.0 && beam_energy_at_vertex > 0.1) + plotter->Fill2D(rx + "_ETrack_vs_EKinGS" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5, + beam_energy_at_vertex, ebeam_kin_GS, "ETrackvsKin_assumed"); + else if (beam_energy_at_vertex <= 12.0) + plotter->Fill2D(rx + "_ETrack_vs_EKin2235keV" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5, + beam_energy_at_vertex, ebeam_kin_2235keV, "ETrackvsKin_assumed"); + else if (beam_energy_at_vertex < 24.0) + { + plotter->Fill2D(rx + "_ETrack_vs_EKin3498kev" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5, + beam_energy_at_vertex, ebeam_kin_3498keV, "ETrackvsKin_assumed"); + plotter->Fill2D(rx + "_ETrack_vs_EKin3774keV" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5, + beam_energy_at_vertex, ebeam_kin_3774keV, "ETrackvsKin_assumed"); + } + else if (beam_energy_at_vertex >= 36.00) + plotter->Fill2D(rx + "_ETrack_vs_EKin4809keV" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5, + beam_energy_at_vertex, ebeam_kin_4809keV, "ETrackvsKin_assumed"); + else if (beam_energy_at_vertex > 42.0) + plotter->Fill2D(rx + "_ETrack_vs_EKin5614keV" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5, + beam_energy_at_vertex, ebeam_kin_5614keV, "ETrackvsKin_assumed"); + else + plotter->Fill2D(rx + "_ETrack_vs_EKin6550keV" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5, + beam_energy_at_vertex, ebeam_kin_6550keV, "ETrackvsKin_assumed"); } }; @@ -3999,7 +4010,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector &Si_ } plotter->Fill1D(rx + "_pczfix" + sfx, 600, -300, 300, pcz_fix, pmlabel); plotter->Fill2D(rx + "_Ef_vs_theta" + ejtag + sfx, 100, 0, 180, 800, 0, ef_max, theta * 180 / M_PI, Efix, pmlabel); - plotter->Fill2D(rx + "_Ex_vs_phi" + ejtag + sfx, 180, -180, 180, 800, -20, 20, phi * 180 / M_PI, Ex, pmlabel); + plotter->Fill2D(rx + "_Ex_vs_phi" + ejtag + sfx, 180, -180, 180, 600, -10, 20, phi * 180 / M_PI, Ex, pmlabel); if (dt_rf_mcp > -900000000) plotter->Fill2D(rx + "_Ex_vs_TOF_rf_mcp" + ejtag + sfx, 500, -1000, 1000, 800, -20, 20, dt_rf_mcp, Ex, pmlabel); plotter->Fill1D(rx + "_VertexReconZ" + sfx, 800, -400, 400, vertex_z, pmlabel); @@ -4029,12 +4040,10 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector &Si_ plotter->Fill2D(rx + "_dEgasRaw_vs_theta" + ejtag + sfx, 180, 0, 180, 800, 0, 20000, theta * 180 / M_PI, anodeE, pmlabel); plotter->Fill2D(rx + "_dEgasCalib_vs_theta" + ejtag + sfx, 360, 0, 180, 800, 0, 0.6, theta * 180 / M_PI, anodeE_MeV, pmlabel); plotter->Fill2D(rx + "_dEgasCalib_vs_phi" + ejtag + sfx, 90, -200, 200, 800, 0, 0.6, phi * 180 / M_PI, anodeE_MeV, pmlabel); - plotter->Fill2D(rx + "_dEgasCalib_vs_E" + ejtag + sfx + "_E<10MeV" + std::to_string(beam_energy_at_vertex < 10), 400, 0, ef_max, 800, 0, 0.6, sievent.Energy1, anodeE_MeV, pmlabel); - if (anodeCh >= 0) - plotter->Fill2D(rx + "_dEgasCalib_vs_E" + ejtag + sfx + "_anode" + pad2(anodeCh), - 400, 0, ef_max, 800, 0, 0.6, sievent.Energy1, anodeE_MeV, pmlabel); - plotter->Fill2D(rx + "_dEgasCalib_vs_Ex" + ejtag + sfx + "_E<10MeV" + std::to_string(beam_energy_at_vertex < 10), 800, -10, 10, 800, 0, 0.6, Ex, anodeE_MeV, pmlabel); - plotter->Fill2D(rx + "_dEgasCalib_vs_Z" + ejtag + sfx + "_E<10MeV" + std::to_string(beam_energy_at_vertex < 10), 800, -400, 400, 800, 0, 0.6, vertex_z, anodeE_MeV, pmlabel); + // if (anodeCh >= 0) + // plotter->Fill2D(rx + "_dEgasCalib_vs_E" + ejtag + sfx + "_anode" + pad2(anodeCh), 400, 0, ef_max, 800, 0, 0.6, sievent.Energy1, anodeE_MeV, pmlabel); + plotter->Fill2D(rx + "_dEgasCalib_vs_Ex" + ejtag + sfx, 600, -10, 20, 800, 0, 0.6, Ex, anodeE_MeV, pmlabel); + plotter->Fill2D(rx + "_dEgasCalib_vs_Z" + ejtag + sfx, 800, -400, 400, 800, 0, 0.6, vertex_z, anodeE_MeV, pmlabel); plotter->Fill2D(rx + "_dEgasPred_vs_dEgasCalib" + ejtag + sfx, 800, 0, 0.6, 800, 0, 0.6, anodeE_MeV, dE_pred, pmlabel); } } diff --git a/run_17F.sh b/run_17F.sh index 229e526..50b388c 100644 --- a/run_17F.sh +++ b/run_17F.sh @@ -48,7 +48,8 @@ root -q -l -b -e '.L TrackRecon.C++O' # 3% CO2 # parallel --bar -j 6 run_once ::: {325..400} -parallel --bar -j 7 run_once ::: 351 353 355 358 359 360 362 367 +parallel --bar -j 11 run_once ::: {351..400} +# parallel --bar -j 7 run_once ::: 351 353 355 358 359 360 362 367 hadd -j 4 -k ${OUT_DIR}/Output_17F.root ${OUT_DIR}/results_run*.root unset source_vertex diff --git a/run_27Al.sh b/run_27Al.sh index 3b9b5fe..c99a6b6 100644 --- a/run_27Al.sh +++ b/run_27Al.sh @@ -8,7 +8,7 @@ export CO2percent=3 export pressure_in_torr=250 export CATHODE_GAIN=3.0 export source_vertex=-200.0 -export DEDX_SCALE=0.89 +export DEDX_SCALE=0.9 export CUTLIST=cuts_list.txt export BEAM_AXIS_X=0 export BEAM_AXIS_Y=0 diff --git a/run_misc_ede.sh b/run_misc_ede.sh new file mode 100644 index 0000000..3d651de --- /dev/null +++ b/run_misc_ede.sh @@ -0,0 +1,93 @@ +#!/bin/bash +# run_misc_ede.sh -- run the standalone MISC E-dE telescope plots (MiscEdE.C) +# over raw, pre-mapping event-built files. Independent of the main TrackRecon +# pipeline: no mapping/calibration/reconstruction, works on beam-monitor runs. +# +# The macro reads the flat sn/ch/e event-built tree and produces per-channel +# spectra + every pairwise E-dE 2D, so you can see all correlations and confirm +# telescopes directly from the data. +# +# Point EVT_DIR / PREFIX / SUFFIX at wherever your event-built files live and +# how they're named, then list the run numbers. Each run is processed into its +# own MiscEdE_.root; set MERGE=1 to also hadd them together. +# +# Usage: +# ./run_misc_ede.sh # uses the RUNS list below +# ./run_misc_ede.sh 351 353 355 # override runs on the command line + +set -u + +# ---- configure these to your layout ---- +EVT_DIR="${EVT_DIR:-../ANASEN_analysis/data/17F_Data}" # where the event-built files are +PREFIX="${PREFIX:-Run_}" # filename prefix before the run number +SUFFIX="${SUFFIX:-_2000.root}" # filename suffix after the run number +PAD="${PAD:-3}" # zero-pad width for the run number (Run_053 -> 3) +OUT_DIR="${OUT_DIR:-Output_misc_ede}" +MERGE="${MERGE:-1}" # 1 = also hadd all per-run outputs +MACRO="${MACRO:-MiscEdE.C}" +CUT_DIR="${CUT_DIR:-/home/vsitaraman/ANASEN_analysis}" # Where 17FCut.root and 16OCut.root live +# ---------------------------------------- + +# Runs: command-line args win, else the lollipop/Si-monitor test runs below. +if [ "$#" -gt 0 ]; then + RUNS=("$@") +else + RUNS=( 52 53 54 55 63 65 68 69 72 74 75 76 78 92 168 170 226 262 324 327 337 334 369 370 371 372 373 374) +fi + +mkdir -p "$OUT_DIR" +OUT_ABS="$(cd "$OUT_DIR" && pwd)" +MACRO_ABS="$(pwd)/$MACRO" +export OUT_ABS MACRO_ABS EVT_DIR PREFIX SUFFIX PAD CUT_DIR + +built_path() { + local run; run=$(printf "%0${PAD}d" "$1") + echo "${EVT_DIR}/${PREFIX}${run}${SUFFIX}" +} +export -f built_path + +# One run -> one MiscEdE_.root, runnable under GNU parallel. +misc_ede_one() { + local r="$1" + local infile; infile=$(built_path "$r") + if [ ! -f "$infile" ]; then + echo "SKIP: $infile not found" + return + fi + local infile_abs; infile_abs="$(cd "$(dirname "$infile")" && pwd)/$(basename "$infile")" + echo "=== MISC E-dE: run $r ($infile_abs) ===" + + # EXPORT METADATA FOR C++ SCRIPT + export CURRENT_RUN="$r" + + # Macro names its own output MiscEdE_.root in the CWD, so run inside + # OUT_DIR to keep outputs together. Absolute paths avoid relative-path + # surprises from the cd. + ( cd "$OUT_ABS" && root -l -b -q "${MACRO_ABS}(\"${infile_abs}\")" ) +} +export -f misc_ede_one + +rm -f "$OUT_DIR"/*.root + +JOBS="${JOBS:-7}" # match the -j 7 you use for run_once +if command -v parallel >/dev/null 2>&1; then + parallel --bar -j "$JOBS" misc_ede_one ::: "${RUNS[@]}" +else + echo "GNU parallel not found; running serially." + for r in "${RUNS[@]}"; do misc_ede_one "$r"; done +fi + +if [ "$MERGE" -eq 1 ]; then + shopt -s nullglob + parts=() + for f in "$OUT_DIR"/MiscEdE_*.root; do + [ "$(basename "$f")" = "MiscEdE_ALL.root" ] && continue + parts+=("$f") + done + if [ "${#parts[@]}" -gt 1 ]; then + echo "=== merging ${#parts[@]} outputs -> ${OUT_DIR}/MiscEdE_ALL.root ===" + hadd -f "${OUT_DIR}/MiscEdE_ALL.root" "${parts[@]}" + fi +fi + +echo "Done. Outputs in ${OUT_DIR}/" \ No newline at end of file diff --git a/run_tr.sh b/run_tr.sh index 3941d03..2969941 100644 --- a/run_tr.sh +++ b/run_tr.sh @@ -55,7 +55,7 @@ if [[ 1 -eq 0 ]]; then fi # --- Block 3: 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" @@ -70,6 +70,7 @@ if [[ 1 -eq 0 ]]; then unset CATHODE_GAIN unset timecut_low unset pressure_in_torr + exit fi # --- Block 4: 17F Alpha+Gas Runs (18-21) --- diff --git a/scratch/make_prettyplots.C b/scratch/make_prettyplots.C index 338a118..51c7a10 100644 --- a/scratch/make_prettyplots.C +++ b/scratch/make_prettyplots.C @@ -153,6 +153,9 @@ void make_prettyplots(const char *rootFile, h->GetYaxis()->SetTitleOffset(1.4); h->GetXaxis()->CenterTitle(true); h->GetYaxis()->CenterTitle(true); + + // gPad->SetLogz(); + h->Draw("colz"); } else