expanded dual-plotter into multi-plotter

Armory/anasenMS.cpp

@@ -97,13 +97,13 @@ int main(int argc, char **argv){
   if( argc >= 2 ) numEvent = atoi(argv[1]);
   TransferReaction transfer;
 
-  transfer.SetA(1, 1, 0);            // 18Ne projectile
+  transfer.SetA(27, 13, 0);            // 18Ne projectile
   //transfer.SetIncidentEnergyAngle(0, 0, 0); // KEA in MeV/u, theta and phi in rad
-  TGraph* elossBeam = LoadELoss("../ELoss/E_vs_x_proton.dat"); 
+  TGraph* elossBeam = LoadELoss("../ELoss/E_vs_x_Al-27.dat"); 
   transfer.Seta(4, 2);                 // 4He target
-  transfer.Setb(1, 1);                 // outgoing proton from the primary transfer
+  transfer.Setb(2, 1);                 // outgoing proton from the primary transfer
-  transfer.SetB(4, 2);               // 21Na* heavy product
+  transfer.SetB(29, 14);               // 21Na* heavy product
-  const double beamA = 1.0;                   // mass number of 27Al beam
+  const double beamA = 27;                   // mass number of 27Al beam
 
   bool enableSequentialDecay = false; // turning to false to disable sequential decay for now, can be set to true to enable
   const int decayDaughterA = 20;
@@ -119,7 +119,7 @@ int main(int argc, char **argv){
   double vertexXRange[2] = {   -5, 5}; // mm - 5, 5
   double vertexYRange[2] = {   -5, 5}; // -5, 5
   double vertexZRange[2] = { -174.3, 174.3}; // -174.3, 174.3 (full length of gas volume, centered at 0)
-  const double beamEntranceZ = -174.3; //vertexZRange[0]; // mm, assumed beam entrance into the gas
+  const double beamEntranceZ = -280; //vertexZRange[0]; // mm, assumed beam entrance into the gas
 
 
   // detector resolution / uncertainty parameters
@@ -553,7 +553,7 @@ int main(int argc, char **argv){
 
       tree2->Fill();
 
-    }else if (qqqID >= 0){
+    }else if (false){//(qqqID >= 0){
       // handle QQQ hit case
       sx3Up = -1;
       sx3Dn = -1;

ELoss/Eloss.cpp

@@ -0,0 +1,248 @@
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <string>
+#include <map>
+#include <cmath>
+#include <memory>
+#include <sstream>
+#include <algorithm>
+#include <filesystem>
+
+// ROOT
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1D.h"
+#include "TH2D.h"
+#include "TCanvas.h"
+
+// Catima
+#include <catima/catima.h>
+
+using namespace std;
+
+/* =========================
+   GLOBAL STRUCTURES
+========================= */
+
+struct Particle {
+    int z;
+    double mass_u;
+    double emax;
+    string name;
+};
+
+map<string, Particle> particles = {
+    {"alpha",    {2, 4.0026, 40, "alpha"}},
+    {"proton",   {1, 1.0078, 20, "proton"}},
+    {"deuteron", {1, 2.0141, 30, "deuteron"}}
+};
+
+/* =========================
+   INTERPOLATION STORAGE
+========================= */
+
+struct EnergyTable {
+    vector<double> x;
+    vector<double> E;
+};
+
+map<string, EnergyTable> table_cache;
+
+/* Linear interpolation */
+double interp(const vector<double>& x, const vector<double>& y, double xi) {
+    if (xi <= x.front()) return y.front();
+    if (xi >= x.back()) return y.back();
+
+    auto it = upper_bound(x.begin(), x.end(), xi);
+    int i = distance(x.begin(), it) - 1;
+
+    double x0 = x[i], x1 = x[i+1];
+    double y0 = y[i], y1 = y[i+1];
+
+    return y0 + (y1 - y0) * (xi - x0) / (x1 - x0);
+}
+
+/* =========================
+   ENERGY TABLE GENERATION
+========================= */
+
+void make_E_vs_x(int z, double mass_u, double emax, string label, int npoints, double P_torr, double T) {
+
+    double R = 8.3144;
+
+    double p_pa = P_torr * 133.322;
+    double molar_density = p_pa / (R * T);
+
+    double m_he = 4.0026;
+    double m_c = 12.0;
+    double m_o = 15.9949;
+
+    double m_mix = 0.96*m_he + 0.04*(m_c + 2*m_o);
+
+    double rho = (molar_density * m_mix) / 1e6;
+
+    catima::Material gas({
+        {m_he, 2, 0.96},
+        {m_c,  6, 0.04},
+        {m_o,  8, 0.08}
+    });
+    gas.density(rho);
+
+    catima::Projectile proj(mass_u, z);
+
+    vector<double> E(npoints);
+    vector<double> S(npoints);
+
+    for (int i = 0; i < npoints; i++) {
+        E[i] = 0.01 + i * (emax / npoints);
+
+        proj.T(E[i] / mass_u);
+        S[i] = catima::dedx(proj, gas) * rho;
+    }
+
+    reverse(E.begin(), E.end());
+    reverse(S.begin(), S.end());
+
+    vector<double> x(npoints, 0.0);
+
+    for (int i = 1; i < npoints; i++) {
+        double invS = 1.0 / S[i];
+        x[i] = x[i-1] + 0.5 * (invS + 1.0/S[i-1]) * (E[i] - E[i-1]);
+    }
+
+    ofstream out("E_vs_x_" + label + ".dat");
+    for (int i = 0; i < npoints; i++) {
+        out << x[i] << "\t" << E[i] << "\n";
+    }
+
+    cout << "Saved E_vs_x_" << label << ".dat\n";
+}
+
+/* =========================
+   LOAD TABLE
+========================= */
+
+EnergyTable load_table(string fname) {
+    EnergyTable t;
+    ifstream in(fname);
+
+    double x, E;
+    while (in >> x >> E) {
+        t.x.push_back(x);
+        t.E.push_back(E);
+    }
+    return t;
+}
+
+/* =========================
+   GET TABLE
+========================= */
+
+EnergyTable& get_table(string particle) {
+    if (!table_cache.count(particle)) {
+        table_cache[particle] = load_table("E_vs_x_" + particle + ".dat");
+    }
+    return table_cache[particle];
+}
+
+/* =========================
+   ENERGY FUNCTIONS
+========================= */
+
+double energy_loss(string p, double Ei, double dx) {
+    auto &t = get_table(p);
+
+    double xi = interp(t.E, t.x, Ei);
+    double xf = xi + dx;
+
+    double Ef = interp(t.x, t.E, xf);
+    return max(Ef, 0.0);
+}
+
+double energy_reconstruction(string p, double Ef, double dx) {
+    auto &t = get_table(p);
+
+    double xf = interp(t.E, t.x, Ef);
+    double xi = xf - dx;
+
+    double Ei = interp(t.x, t.E, xi);
+    return max(Ei, 0.0);
+}
+
+double energy_distance(string p, double Ei, double Ef) {
+    auto &t = get_table(p);
+
+    double xi = interp(t.E, t.x, Ei);
+    double xf = interp(t.E, t.x, Ef);
+
+    return fabs(xf - xi);
+}
+
+/* =========================
+   ROOT FILE ANALYSIS
+========================= */
+
+void process_file(string filename) {
+
+    TFile f(filename.c_str());
+    TTree *tree = (TTree*)f.Get("tree");
+
+    double Tb, thetab, sx3Z;
+
+    tree->SetBranchAddress("Tb", &Tb);
+    tree->SetBranchAddress("thetab", &thetab);
+    tree->SetBranchAddress("sx3Z", &sx3Z);
+
+    vector<double> Ei, theta, sx3;
+
+    Long64_t n = tree->GetEntries();
+
+    for (Long64_t i = 0; i < n; i++) {
+        tree->GetEntry(i);
+
+        double th = thetab * M_PI / 180.0;
+        if (sin(th) == 0) continue;
+
+        Ei.push_back(Tb);
+        theta.push_back(th);
+        sx3.push_back(sx3Z);
+    }
+
+    cout << "Processed " << Ei.size() << " events\n";
+}
+
+/* =========================
+   MAIN CLI (simplified)
+========================= */
+
+int main() {
+
+    string cmd;
+
+    cout << "C++ PCEnergy Shell\n";
+
+    while (true) {
+        cout << ">> ";
+        getline(cin, cmd);
+
+        if (cmd == "exit") break;
+
+        if (cmd.rfind("make table", 0) == 0) {
+            string name = cmd.substr(11);
+            auto p = particles[name];
+
+            make_E_vs_x(p.z, p.mass_u, p.emax, p.name, 500, 400, 293.15);
+        }
+
+        else if (cmd.rfind("energy loss", 0) == 0) {
+            cout << "Use API call version in compiled mode\n";
+        }
+
+        else {
+            cout << "Unknown command\n";
+        }
+    }
+
+    return 0;
+}

ELoss/PCEnergyAnalysis.py

@@ -1345,46 +1345,89 @@ class MyInteractiveApp(cmd.Cmd):
         
         
     def do_dual_plotter(self, arg):
+    
         args = shlex.split(arg)
     
-        if len(args) < 2:
+        # Default files if none provided
-            try:
+        if len(args) == 0:
-                args = ['SimAnasenProton.root', 'SimAnasenAlpha.root']
+            files = [
-            except:
+                "SimAnasenProton.root",
-                print("Usage: make dual plots proton_data.root alpha_data.root")
+                "SimAnasenAlpha.root"
-                return
+            ]
+        else:
+            files = args
     
-        file1 = args[0]
-        file2 = args[1]
-        #tree1_name = args[2] if len(args) > 2 else 'tree2'
-        #tree2_name = args[3] if len(args) > 3 else 'tree1'
         outdir = "dual_plots"
-        
-        # load both trees for file1 and combine their arrays into a single dataset
-        data1_tree1 = process_file(os.path.join("..", "Armory", file1), "tree1")
-        data1_tree2 = process_file(os.path.join("..", "Armory", file1), "tree1")
-        # concatenate matching array entries
-        data1 = {"particle": f"{data1_tree1['particle']}_combined"}
-        for key in data1_tree1:
-            if key == "particle":
-                continue
-            try:
-                data1[key] = np.concatenate([data1_tree1[key], data1_tree2[key]])
-            except Exception:
-                # fallback: prefer tree1 value if concat fails
-                data1[key] = data1_tree1[key]
-
-        data2 = process_file(os.path.join("..", "Armory", file2), "tree1")
-
-        #print(f"File one {file1} ({tree1_name}) length {len(data1['Ei'])} \nFile two {file2} ({tree2_name}) length {len(data2['Ei'])}")
         os.makedirs(outdir, exist_ok=True)
     
         print(f"Saving plots to: {outdir}")
     
-        #Overlay histogram: Elost
+        datasets = []
-        plt.figure(figsize=(8,6))
+    
-        plt.hist(data1["Elost"],bins=200,histtype='step',linewidth=2,density=True,label=data1["particle"])
+        for file in files:
-        plt.hist(data2["Elost"],bins=200,histtype='step',linewidth=2,density=True,label=data2["particle"])
+    
+            try:
+                # If you want to combine tree1 + tree2:
+                tree1 = process_file(
+                    os.path.join("..", "Armory", file),
+                    "tree1"
+                )
+    
+                try:
+                    tree2 = process_file(
+                        os.path.join("..", "Armory", file),
+                        "tree2"
+                    )
+    
+                    data = {
+                        "particle":
+                        f"{tree1['particle']}_combined"
+                    }
+    
+                    for key in tree1:
+    
+                        if key == "particle":
+                            continue
+    
+                        try:
+                            data[key] = np.concatenate(
+                                [tree1[key], tree2[key]]
+                            )
+                        except Exception:
+                            data[key] = tree1[key]
+    
+                except Exception:
+                    data = tree1
+    
+                datasets.append(data)
+    
+                print(
+                    f"Loaded {file} "
+                    f"({len(data['Ei'])} events)"
+                )
+    
+            except Exception as e:
+                print(f"Failed to load {file}: {e}")
+    
+        if len(datasets) == 0:
+            print("No valid datasets loaded.")
+            return
+    
+        # --------------------------------------------------
+        # Elost overlay
+        # --------------------------------------------------
+    
+        plt.figure(figsize=(8, 6))
+    
+        for data in datasets:
+            plt.hist(
+                data["Elost"],
+                bins=200,
+                histtype="step",
+                linewidth=2,
+                density=True,
+                label=data["particle"]
+            )
     
         plt.xlabel("Energy Loss (MeV)")
         plt.ylabel("Normalized Counts")
@@ -1393,194 +1436,218 @@ class MyInteractiveApp(cmd.Cmd):
         plt.grid(True)
         plt.tight_layout()
     
-        plt.savefig(f"{outdir}/Elost_overlay.png", dpi=300)
+        plt.savefig(
+            f"{outdir}/Elost_overlay.png",
+            dpi=300
+        )
+    
         plt.show()
     
-        #Overlay histogram: sx3Z
+        # --------------------------------------------------
-        plt.figure(figsize=(8,6))
+        # SX3 position overlay
-        plt.hist(data1["sx3Z"],bins=150,histtype='step',linewidth=2,density=True,label=data1["particle"])
+        # --------------------------------------------------
-        plt.hist(data2["sx3Z"],bins=150,histtype='step',linewidth=2,density=True,label=data2["particle"])
+    
+        plt.figure(figsize=(8, 6))
+    
+        for data in datasets:
+            plt.hist(
+                data["sx3Z"],
+                bins=150,
+                histtype="step",
+                linewidth=2,
+                density=True,
+                label=data["particle"]
+            )
+    
         plt.xlabel("SX3 Z")
         plt.ylabel("Normalized Counts")
         plt.title("SX3 Position Comparison")
         plt.legend()
         plt.grid(True)
         plt.tight_layout()
-        plt.savefig(f"{outdir}/sx3Z_overlay.png", dpi=300)
+    
+        plt.savefig(
+            f"{outdir}/sx3Z_overlay.png",
+            dpi=300
+        )
+    
         plt.show()
-        try:
-            #Side-by-side 2D plots
-            fig, axes = plt.subplots(1, 2, figsize=(14,6))
     
-            h1 = axes[0].hist2d(data1["sx3Z"],data1["Elost"],bins=200)
+        # --------------------------------------------------
+        # Individual Elost vs SX3 plots
+        # --------------------------------------------------
     
-            axes[0].set_title(f'{data1["particle"]} Elost vs SX3')
+        n = len(datasets)
-            axes[0].set_xlabel("SX3 Z")
-            axes[0].set_ylabel("Energy Loss (MeV)")
     
-            h2 = axes[1].hist2d(data2["sx3Z"],data2["Elost"],bins=200)
+        fig, axes = plt.subplots(
+            1,
+            n,
+            figsize=(7 * n, 6)
+        )
     
-            axes[1].set_title(f'{data2["particle"]} Elost vs SX3')
+        if n == 1:
-            axes[1].set_xlabel("SX3 Z")
+            axes = [axes]
-            axes[1].set_ylabel("Energy Loss (MeV)")
     
-            fig.colorbar(h1[3], ax=axes[0], label="Counts")
+        for ax, data in zip(axes, datasets):
-            fig.colorbar(h2[3], ax=axes[1], label="Counts")
     
-            plt.tight_layout()
+            h = ax.hist2d(
+                data["sx3Z"],
+                data["Elost"],
+                bins=200
+            )
     
-            plt.savefig(f"{outdir}/Elost_vs_sx3_comparison.png", dpi=300)
+            ax.set_title(
-            plt.show()
+                f'{data["particle"]}\nElost vs SX3'
-        except:
+            )
-            print("Error with side-by-side plots")
+    
+            ax.set_xlabel("SX3 Z")
+            ax.set_ylabel("Energy Loss (MeV)")
+    
+            fig.colorbar(
+                h[3],
+                ax=ax,
+                label="Counts"
+            )
+    
+        plt.tight_layout()
+    
+        plt.savefig(
+            f"{outdir}/Elost_vs_sx3_comparison.png",
+            dpi=300
+        )
+    
+        plt.show()
+    
+        # --------------------------------------------------
+        # EA vs Esx3 overlay scatter
+        # --------------------------------------------------
+    
+        plt.figure(figsize=(8, 6))
+    
+        for data in datasets:
+    
+            plt.scatter(
+                data["EA"],
+                data["Esx3"],
+                s=1,
+                alpha=0.3,
+                label=data["particle"]
+            )
     
-        #EA vs Esx3 overlay scatter
-        plt.figure(figsize=(8,6))
-        plt.scatter(data1["EA"],data1["Esx3"],s=1,alpha=0.3,label=data1["particle"])
-        plt.scatter(data2["EA"],data2["Esx3"],s=1,alpha=0.3,label=data2["particle"])
         plt.xlabel("EA (MeV)")
         plt.ylabel("Esx3 (MeV)")
         plt.title("Anode vs SX3 Energy")
         plt.legend()
         plt.grid(True)
         plt.tight_layout()
-        plt.savefig(f"{outdir}/EA_vs_Esx3_overlay.png", dpi=300)
-        plt.show()
     
-        #PCE vs SiE
+        plt.savefig(
-        
+            f"{outdir}/EA_vs_Esx3_overlay.png",
-        mask1 = data1["Esx3"] > 1
+            dpi=300
-        mask2 = data2["Esx3"] > 1
-        
-        thetab1 = np.deg2rad(data1["thetab"][mask1])
-        thetab2 = np.deg2rad(data2["thetab"][mask2])
-        
-        #theta smear (spatial uncertainty)
-        if False:
-            sigma = 0.2 * np.sqrt(np.maximum(thetab1, 0))
-            
-            thetab1 += np.random.normal(0,sigma,len(thetab1))
-            thetab1 += np.random.normal(0,0.02*np.sqrt(thetab1),len(thetab1))
-            
-            sigma = 0.2 * np.sqrt(np.maximum(thetab2, 0))
-            thetab2 += np.random.normal(0,sigma,len(thetab2))
-            thetab2 += np.random.normal(0,0.02*np.sqrt(thetab2),len(thetab2))
-        
-        combined_Esx3 = np.concatenate([
-            (data1["Esx3"][mask1]),
-            data2["Esx3"][mask2]
-        ])
-        combined_Eprop = np.concatenate([
-            data1["Eprop"][mask1] * np.sin(thetab1) * 3,
-            data2["Eprop"][mask2] * np.sin(thetab2) * 3
-        ])
-        
-        combined_Esx3 = (
-            combined_Esx3
-            + np.random.normal(0,0.08,len(combined_Esx3))
         )
     
-        #Artifical smear
+        plt.show()
-        if False:
-            sigma = 0.02 * np.sqrt(np.maximum(combined_Eprop, 0))
     
-            combined_Eprop += np.random.normal(
+        # --------------------------------------------------
-                0,
+        # Combined PCE vs SX3
-                sigma,
+        # --------------------------------------------------
-                len(combined_Eprop)
+    
+        all_Esx3 = []
+        all_Eprop = []
+    
+        for data in datasets:
+    
+            mask = data["Esx3"] > 1
+    
+            thetab = np.deg2rad(
+                data["thetab"][mask]
             )
     
-            combined_Eprop += np.random.normal(
+            all_Esx3.append(
-                0,
+                data["Esx3"][mask]
-                0.02*np.sqrt(combined_Eprop),
-                len(combined_Eprop)
             )
     
-        mask = (np.isfinite(combined_Esx3)& np.isfinite(combined_Eprop))
+            all_Eprop.append(
+                data["Eprop"][mask]
+                * np.sin(thetab)
+                * 3
+            )
+    
+        combined_Esx3 = np.concatenate(all_Esx3)
+        combined_Eprop = np.concatenate(all_Eprop)
+    
+        combined_Esx3 += np.random.normal(
+            0,
+            0.08,
+            len(combined_Esx3)
+        )
+    
+        mask = (
+            np.isfinite(combined_Esx3)
+            &
+            np.isfinite(combined_Eprop)
+        )
+    
         combined_Esx3 = combined_Esx3[mask]
         combined_Eprop = combined_Eprop[mask]
     
-        #combined_Eprop = combined_Eprop * .686
+        plt.figure(figsize=(8, 6))
     
-        plt.figure(figsize=(8,6))
         plt.hist2d(
             combined_Esx3,
             combined_Eprop,
             bins=200
         )
+    
         plt.xlabel("SX3 Energy (MeV)")
         plt.ylabel("PCEnergy x Sin(theta)")
-        #plt.xlim(0,30)
-        #plt.ylim(0,0.45)
         plt.colorbar(label="Counts")
+    
         plt.tight_layout()
-        plt.savefig(f"{outdir}/Eprop_vs_Esx3.png", dpi=300)
+    
+        plt.savefig(
+            f"{outdir}/Eprop_vs_Esx3.png",
+            dpi=300
+        )
+    
         plt.show()
     
-        plt.figure(figsize=(14,6), facecolor='white')
+        plt.figure(
+            figsize=(14, 6),
+            facecolor="white"
+        )
+    
         plt.hist2d(
             combined_Esx3,
             combined_Eprop,
-            bins=[500,500],
+            bins=[500, 500],
-            #range=[[0,35],[0,0.6]],
             cmap=yellow_jet,
-            cmin=1)
+            cmin=1
-        plt.margins(0)
+        )
+    
         plt.xlabel("SX3 Energy (MeV)")
         plt.ylabel("PCEnergy x Sin(theta)")
-        #plt.xlim(0,35)
+    
-        #plt.ylim(0,0.6)
         cbar = plt.colorbar()
         cbar.set_label("Counts")
+    
         plt.minorticks_on()
         plt.grid(False)
+    
         plt.tight_layout()
+    
         plt.savefig(
             f"{outdir}/ROOT_style_plot.png",
             dpi=300,
-            facecolor='white'
+            facecolor="white"
         )
+    
         plt.show()
-        if False:
-            plt.figure(figsize=(7,6))
-            plt.hist2d(
-                combined_Esx3,
-                combined_Eprop,
-                bins=100,
-                range=[[0,30],[0,0.45]],
-                cmap='viridis'
-            )
     
-            mask = data1["Esx3"] > 0
+        print(
-            power_fit_and_plot(
+            f"Completed plotting "
-                data1["Esx3"][mask],
+            f"{len(datasets)} datasets."
-                data1["Eprop"][mask],
+        )
-                label=data1["particle"],
-                color='red'
-            )
-            mask2 = data2["Esx3"] > 0
-            power_fit_and_plot(
-                data2["Esx3"][mask2],
-                data2["Eprop"][mask2],
-                label=data2["particle"],
-                color='cyan'
-            )
-            
-            plt.ylabel("PCEnergy")
-            plt.xlabel("SX3 Energy (MeV)")
-            plt.title(
-                f'{data1["particle"]} + {data2["particle"]} '
-                'Energy Propagation Difference vs SX3 Energy'
-            )
-            plt.colorbar(label="Counts")
-            #plt.xlim(0,30)
-            #plt.ylim(0,0.45)
-            plt.legend()
-            plt.tight_layout()
-            plt.savefig(f"{outdir}/Combined_Eprop_vs_Esx3_fit.png",dpi=300)
-            plt.show()
-        
-        print("Dual plotting complete.")
             
 #exec(open("PCEnergyAnalysis.py").read())