expanded analysis

Armory/ClassPW.h

@@ -7,8 +7,8 @@
 #include <TVector3.h>
 #include <TRandom.h>
 
-std::vector<int> skipAnodes = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23}; 
+std::vector<int> skipAnodes = {}; 
-std::vector<int> skipCathodes = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23};
+std::vector<int> skipCathodes = {};
 
 struct PWHitInfo
 {

Armory/anasenMS.cpp

@@ -46,8 +46,6 @@ bool IsDeadSX3(int id){
 }
 
 // Simulate sequential two-body decay of an unstable parent in its rest frame.
-// The parent is boosted from the lab frame, the daughter (A1,Z1) is returned in lab frame,
-// and the emitted ejectile (A2,Z2) is written to ejectileOut.
 TLorentzVector SimulateSequentialDecay(const TLorentzVector &parent,
                                        int daughterA, int daughterZ,
                                        int ejectA, int ejectZ,
@@ -62,28 +60,28 @@ TLorentzVector SimulateSequentialDecay(const TLorentzVector &parent,
   double sqM = M * M;
   double sum = mD + mE;
   double diff = mD - mE;
-  double p2 = (sqM - sum*sum) * (sqM - diff*diff) / (4.0 * sqM);
+  double p2 = (sqM - sum*sum) * (sqM - diff*diff) / (4.0 * sqM); // two-body decay momentum squared
-  if( p2 < 0 ) p2 = 0;
+  if( p2 < 0 ) p2 = 0; // handle unphysical case where parent mass is less than sum of daughter and ejectile masses
-  double p = TMath::Sqrt(p2);
+  double p = TMath::Sqrt(p2); // two-body decay momentum
 
-  double cosTheta = 2.0 * gRandom->Rndm() - 1.0;
+  double cosTheta = 2.0 * gRandom->Rndm() - 1.0; // isotropic decay in parent rest frame
-  double theta = TMath::ACos(cosTheta);
+  double theta = TMath::ACos(cosTheta); // polar angle of daughter in parent rest frame
-  double phi = gRandom->Rndm() * TMath::TwoPi();
+  double phi = gRandom->Rndm() * TMath::TwoPi(); // azimuthal angle of daughter in parent rest frame
 
-  TVector3 v;
+  TVector3 v; // momentum vector of daughter in parent rest frame
-  v.SetMagThetaPhi(p, theta, phi);
+  v.SetMagThetaPhi(p, theta, phi); // daughter momentum in parent rest frame
 
-  TLorentzVector daughterLab;
+  TLorentzVector daughterLab; // daughter 4-vector in lab frame, initialized with momentum from decay and mass of daughter
-  daughterLab.SetVectM(v, mD);
+  daughterLab.SetVectM(v, mD); // set daughter 4-vector in parent rest frame, then boost to lab frame
 
-  TLorentzVector ejectileLab;
+  TLorentzVector ejectileLab; // ejectile 4-vector in lab frame, initialized with momentum opposite to daughter and mass of ejectile
-  ejectileLab.SetVectM(-v, mE);
+  ejectileLab.SetVectM(-v, mE); //  set ejectile 4-vector in parent rest frame, then boost to lab frame
 
-  TVector3 boost = parent.BoostVector();
+  TVector3 boost = parent.BoostVector(); // boost vector to go from parent rest frame to lab frame
-  daughterLab.Boost(boost);
+  daughterLab.Boost(boost); // boost daughter to lab frame
-  ejectileLab.Boost(boost);
+  ejectileLab.Boost(boost); // boost ejectile to lab frame
 
-  ejectileOut = ejectileLab;
+  ejectileOut = ejectileLab; // return ejectile in lab frame
   return daughterLab;
 }
 
@@ -96,17 +94,15 @@ int main(int argc, char **argv){
   // number of events can be overridden from command line
   int numEvent = 1000000;
   if( argc >= 2 ) numEvent = atoi(argv[1]);
-  // Reaction setup: 18Ne + 4He -> p + 21Na*
-  // then sequential decay of 21Na* -> p + 20Ne
   TransferReaction transfer;
 
-  transfer.SetA(18, 10, 0);            // 18Ne projectile
+  transfer.SetA(14, 7, 0);            // 18Ne projectile
-  transfer.SetIncidentEnergyAngle((4.4), 0, 0); // KEA in MeV/u, theta and phi in degree
+  transfer.SetIncidentEnergyAngle((42.82/14), 0, 0); // KEA in MeV/u, theta and phi in rad
   transfer.Seta(4, 2);                 // 4He target
   transfer.Setb(1, 1);                 // outgoing proton from the primary transfer
-  transfer.SetB(21, 11);               // 21Na* heavy product
+  transfer.SetB(17, 8);               // 21Na* heavy product
 
-  bool enableSequentialDecay = true; // turning to false to disable sequential decay for now, can be set to true to enable
+  bool enableSequentialDecay = false; // turning to false to disable sequential decay for now, can be set to true to enable
   const int decayDaughterA = 20;
   const int decayDaughterZ = 10;
   const int decayEjectA = 1;
@@ -114,7 +110,7 @@ int main(int argc, char **argv){
 
   // Excited state lists (projectile and heavy-product excitation states)
   std::vector<float> ExAList = {0};              // 18Ne projectile excitations in MeV
-  std::vector<float> ExList = {2.5};             // 21Na* excitation in MeV
+  std::vector<float> ExList = {0};             // 21Na* excitation in MeV
 
   // define vertex position uniform distribution ranges (mm)
   double vertexXRange[2] = {   -5,    5}; // mm
@@ -122,8 +118,8 @@ int main(int argc, char **argv){
   double vertexZRange[2] = { -100,  100};
 
   // detector resolution / uncertainty parameters
-  double sigmaSX3_W = -1; // mm, if < 0 use mid-point (no spread in SX3 horizontal dimension)
+  double sigmaSX3_W = 0; // mm, if < 0 use mid-point (no spread in SX3 horizontal dimension)
-  double sigmaSX3_L = 3;  // mm, vertical spread for SX3
+  double sigmaSX3_L = 0;  // mm, vertical spread for SX3
   double sigmaPW_A  = 0;  // normalized anode uncertainty term (0-1)
   double sigmaPW_C  = 0;  // normalized cathode uncertainty term (0-1)
 
@@ -324,7 +320,7 @@ int main(int argc, char **argv){
     transfer.CalReactionConstant();
 
     // isotropic CM direction
-    thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ;
+    thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ; // polar angle in CM frame
     phiCM   = (gRandom->Rndm() - 0.5) * TMath::TwoPi();
 
     //==== Calculate reaction kinematics in lab frame for the primary transfer
@@ -341,7 +337,7 @@ int main(int argc, char **argv){
     T[0] = Tb;
     T[1] = TB;
 
-    // prepare secondary proton from 21Na* sequential decay
+    //secondary decay
     TLorentzVector decayProton;
     TLorentzVector heavy20;
     if(enableSequentialDecay){

ELoss/PCEnergyAnalysis.py

@@ -298,7 +298,7 @@ def calculate_distance_tree2(vz, theta, z_max=34.86):
 
 
 def load_tree_arrays(tree, treename, max_events=None):
-    branches = ["Tb", "thetab", "sx3Z", "vX", "vY", "vZ"]
+    branches = ["Tb", "thetab", "sx3Z", "vX", "vY", "vZ", "sx3ID", "aID", "cID", "aDist", "cDist"]
     if treename == 'tree1':
         branches.extend(["sx3X", "sx3Y"])
     return tree.arrays(branches, library="np", entry_stop=max_events)
@@ -307,6 +307,8 @@ def load_tree_arrays(tree, treename, max_events=None):
 def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
     data = load_tree_arrays(tree, treename, max_events)
 
+    mask = data["thetab"] >= 0
+    data = {key: value[mask] for key, value in data.items()}
     Ei = data["Tb"]
     theta = np.radians(data["thetab"])
     vX = data["vX"]
@@ -339,10 +341,38 @@ def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
     sin_theta = np.sin(theta)
     sin_theta = np.where(sin_theta != 0, sin_theta, 1e-10)
     radii = np.array([3.7, 4.2])
-    dA = radii[0] / sin_theta
+    dA = (radii[0] - np.sqrt((vX/10)**2 + (vY/10)**2))/ sin_theta
-    dC = radii[1] / sin_theta
+    dC = (radii[1] - np.sqrt((vX/10)**2 + (vY/10)**2))/ sin_theta
+
+    # Filter out unphysical distances (negative or unreasonably small)
+    # These typically occur at large angles where trajectory doesn't properly intersect proportional counters
+    sx3ID = data["sx3ID"]
+    aID = data["aID"]
+    cID = data["cID"]
+    
+    # Keep only events with valid PW wire assignments and positive distances
+    # sx3ID >= 0 means SX3 was hit, aID/cID valid means tracks found
+    # Handle both 1D and 2D array structures from ROOT
+    if aID.ndim == 2:
+        aID_valid = aID[:, 0] >= 0
+        cID_valid = cID[:, 0] >= 0
+    else:
+        aID_valid = aID >= 0
+        cID_valid = cID >= 0
+    
+    distance_mask = (dA > 0.1) & (dC > 0.1) & (sx3ID >= 0) & aID_valid & cID_valid
+    
+    Ei = Ei[distance_mask]
+    theta = theta[distance_mask]
+    dA = dA[distance_mask]
+    dC = dC[distance_mask]
+    dsx3 = dsx3[distance_mask]
+    vX = vX[distance_mask]
+    vY = vY[distance_mask]
+    vZ = vZ[distance_mask]
 
     print(f"Computing energies for {particle} ({treename})...")
+    print(f"  Retained {np.sum(distance_mask)} / {len(distance_mask)} events after distance filter")
 
     EA = energy_loss(particle, Ei, dA)
     EC = energy_loss(particle, Ei, dC)
@@ -359,7 +389,10 @@ def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
         "EC": EC,
         "Esx3": Esx3,
         "Elost": Elost,
-        "Eprop": Eprop
+        "Eprop": Eprop,
+        "dA": dA,
+        "dC": dC,
+        "thetab": data["thetab"]
     }
 
 
@@ -377,13 +410,60 @@ def process_file(filename, treename, particle=None, max_events=None):
     print(f"File {filename} particle {particle}, tree {treename}")
     return prepare_tree_data(tree, treename, particle, max_events=max_events)
 
+
+def power_fit_and_plot(x, y, label, color=None):
+
+    x = np.array(x)
+    y = np.array(y)
+
+    mask = (
+        np.isfinite(x) &
+        np.isfinite(y) &
+        (x > 0) &
+        (y > 0)
+    )
+
+    x = x[mask]
+    y = y[mask]
+
+    logx = np.log(x)
+    logy = np.log(y)
+
+    b, loga = np.polyfit(logx, logy, 1)
+
+    a = np.exp(loga)
+
+    y_pred = a * x**b
+
+    ss_res = np.sum((y - y_pred)**2)
+    ss_tot = np.sum((y - np.mean(y))**2)
+
+    r2 = 1 - (ss_res / ss_tot)
+
+    x_fit = np.linspace(np.min(x), np.max(x), 1000)
+    y_fit = a * x_fit**b
+
+    plt.plot(
+        x_fit,
+        y_fit,
+        linewidth=3,
+        color=color,
+        label=f'{label} fit ($R^2$={r2:.4f})'
+    )
+
+    print(f"\n{label} power-law fit:")
+    print(f"y = {a:.6f} * x^{b:.6f}")
+    print(f"R^2 = {r2:.6f}")
+
+    return a, b, r2
+
 class MyInteractiveApp(cmd.Cmd):
     def __init__(self):
         super().__init__()
         # Initial value set when the script starts
         self.buckets = 500
         self.T = 293.15
-        self.P = 400
+        self.P = 379
         self.temp_particle = [0, 0.0, 0.0, ""]
         self.rootFile = "SimAnasen1.root"
         self.file = None
@@ -821,7 +901,7 @@ class MyInteractiveApp(cmd.Cmd):
         global EA
         EA = energy_loss(particle, Ei, dA)
         global EC
-        EC = energy_loss(particle, Ei, dC)
+        EC = energy_loss(particle, Ei, ClassSX3dC)
         global Esx3
         Esx3 = energy_loss(particle, Ei, dsx3)
         global Eprop
@@ -839,7 +919,7 @@ class MyInteractiveApp(cmd.Cmd):
 
         print(f"sx3 average energy: {np.mean(Esx3):.3f} MeV")
 
-        print(f"Average total energy loss to sx3: {np.mean(Elost):.3f} MeV")
+        print(f"Average total energy loClassSX3ss to sx3: {np.mean(Elost):.3f} MeV")
 
         print(f"Maximum total energy loss to sx3: {np.max(Elost):.3f} MeV")
 
@@ -930,13 +1010,17 @@ class MyInteractiveApp(cmd.Cmd):
 
         data = prepare_tree_data(self.tree, treename, particle, max_events=max_events)
 
-        Ei = data["Ei"]
+        mask = data["thetab"] >= 0
-        sx3Z = data["sx3Z"]
+
-        EA = data["EA"]
+        Ei = data["Ei"][mask]
-        EC = data["EC"]
+        sx3Z = data["sx3Z"][mask]
-        Esx3 = data["Esx3"]
+        EA = data["EA"][mask]
-        Elost = data["Elost"]
+        EC = data["EC"][mask]
-        Eprop = data["Eprop"]
+        Esx3 = data["Esx3"][mask]
+        Elost = data["Elost"][mask]
+        Eprop = data["Eprop"][mask]
+        dA = data["dA"][mask]
+        dC = data["dC"][mask]
 
         update_plot_data(f"{particle}_{treename}_Ei", Ei)
         update_plot_data(f"{particle}_{treename}_sx3Z", sx3Z)
@@ -951,48 +1035,6 @@ class MyInteractiveApp(cmd.Cmd):
 
         print(f"Saving plots to folder: {base} ({treename})")
 
-        branch_names = []
-        for key in self.tree.keys():
-            if isinstance(key, bytes):
-                branch_names.append(key.decode("ascii"))
-            elif hasattr(key, "name"):
-                branch_names.append(key.name)
-            else:
-                branch_names.append(str(key))
-        branch_names = list(dict.fromkeys(branch_names))
-
-        if branch_names:
-            print(f"Creating histograms for {len(branch_names)} branches...")
-            all_branches = self.tree.arrays(branch_names, library="np", entry_stop=max_events)
-
-            for branch in branch_names:
-                values = all_branches[branch]
-                try:
-                    values = np.asarray(values, dtype=float)
-                except Exception:
-                    print(f"Skipping non-numeric branch: {branch}")
-                    continue
-                if values.ndim != 1:
-                    print(f"Skipping non-scalar branch: {branch}")
-                    continue
-                values = values[~np.isnan(values)]
-                if values.size == 0:
-                    print(f"Skipping empty branch: {branch}")
-                    continue
-
-                plt.figure(figsize=(7,5))
-                plt.hist(values, bins=100)
-                plt.xlabel(branch)
-                plt.ylabel("Counts")
-                plt.title(f"{particle} ({treename}) {branch} distribution")
-                plt.grid(True)
-                plt.tight_layout()
-                safe_name = re.sub(r"[^0-9A-Za-z_-]", "_", branch)
-                plt.savefig(f"{base}/{safe_name}_hist.png", dpi=300)
-                plt.close()
-        else:
-            print("No branches found to histogram.")
-
         plt.figure(figsize=(7,5))
         plt.hist(Elost, bins=200)
         plt.xlabel("Energy Loss (MeV)")
@@ -1028,6 +1070,18 @@ class MyInteractiveApp(cmd.Cmd):
         except ValueError:
             print("Value error")
             
+        plt.figure(figsize=(7,5))
+        plt.hist(dA, bins=100, label='dA', histtype='step')
+        plt.hist(dC, bins=100, label='dC', histtype='step')
+        plt.xlabel("Distance")
+        plt.ylabel("Counts")
+        plt.title(f"{particle} ({treename}) Anode/ Cathode distance")
+        plt.legend(loc='upper right')
+        plt.grid(True)
+        plt.tight_layout()
+        plt.savefig(f"{base}/d_hist.png", dpi=300)
+        plt.show()
+
         try:
             plt.figure(figsize=(7,6))
             plt.hist2d(sx3Z, Elost, bins=200)
@@ -1075,14 +1129,85 @@ class MyInteractiveApp(cmd.Cmd):
             plt.xlabel("SX3 Energy (MeV)")
             plt.title(f"{particle} ({treename}) Energy Propagation Difference vs SX3 Energy")
             plt.colorbar(label="Counts")
-            plt.xlim(0,30)
+            #plt.xlim(0,30)
-            plt.ylim(0,0.45)
+            #plt.ylim(0,0.45)
             plt.tight_layout()
             plt.savefig(f"{base}/Eprop_vs_Esx3.png", dpi=300)
             plt.show()
         except ValueError:
             print("Value error")
             
+        try:
+            mask1 = ~np.isnan(Ei) & (Esx3 > 0)
+            plt.figure(figsize=(7,6))
+            plt.hist2d(data["thetab"], Ei, bins=200)
+            plt.xlabel("thetab")
+            plt.ylabel("Event Energy (MeV)")
+            plt.title(f"{particle} ({treename}) Energy vs Theta")
+            plt.colorbar(label="Counts")
+            #plt.xlim(0,30)
+            #plt.ylim(0,0.45)
+            plt.tight_layout()
+            plt.savefig(f"{base}/E_vs_theta.png", dpi=300)
+            plt.show()
+        except ValueError:
+            print("Value error")
+            
+            
+        branch_names = []
+        for key in self.tree.keys():
+            if isinstance(key, bytes):
+                branch_names.append(key.decode("ascii"))
+            elif hasattr(key, "name"):
+                branch_names.append(key.name)
+            else:
+                branch_names.append(str(key))
+        branch_names = list(dict.fromkeys(branch_names))
+        
+        if branch_names:
+            print(f"Creating histograms for {len(branch_names)} branches...")
+            all_branches = self.tree.arrays(branch_names, library="np", entry_stop=max_events)
+        
+            # Keep only events with thetab >= 45
+            #thetab_mask = all_branches["thetab"] >= 0
+            thetab_mask = mask
+        
+            for branch in branch_names:
+                values = all_branches[branch]
+        
+                try:
+                    values = np.asarray(values, dtype=float)
+                except Exception:
+                    print(f"Skipping non-numeric branch: {branch}")
+                    continue
+        
+                if values.ndim != 1:
+                    print(f"Skipping non-scalar branch: {branch}")
+                    continue
+        
+                # Apply the thetab cut
+                values = values[thetab_mask]
+        
+                values = values[~np.isnan(values)]
+        
+                if values.size == 0:
+                    print(f"Skipping empty branch: {branch}")
+                    continue
+        
+                plt.figure(figsize=(7,5))
+                plt.hist(values, bins=100)
+                plt.xlabel(branch)
+                plt.ylabel("Counts")
+                plt.title(f"{particle} ({treename}) {branch} distribution")
+                plt.grid(True)
+                plt.tight_layout()
+        
+                safe_name = re.sub(r"[^0-9A-Za-z_-]", "_", branch)
+                plt.savefig(f"{base}/{safe_name}_hist.png", dpi=300)
+                plt.close()
+        else:
+            print("No branches found to histogram.")
+    
         print(f"Plotting complete ({treename}).") 
                 
     def do_hist_comp(self, arg):
@@ -1302,10 +1427,40 @@ class MyInteractiveApp(cmd.Cmd):
             data2["Esx3"][mask2]
         ])
         combined_Eprop = np.concatenate([
-            data1["Eprop"][mask1],
+            data1["Eprop"][mask1] * 3,
-            data2["Eprop"][mask2]
+            data2["Eprop"][mask2] * 2.7
         ])
         
+        combined_Esx3 = (
+            combined_Esx3
+            + np.random.normal(0,0.08,len(combined_Esx3))
+        )
+        
+        if True:
+            sigma = 0.02 * np.sqrt(np.maximum(combined_Eprop, 0))
+            
+            combined_Eprop += np.random.normal(
+                0,
+                sigma,
+                len(combined_Eprop)
+            )
+            
+            combined_Eprop += np.random.normal(
+                0,
+                0.02*np.sqrt(combined_Eprop),
+                len(combined_Eprop)
+            )
+        
+        mask = (
+            np.isfinite(combined_Esx3)
+            & np.isfinite(combined_Eprop)
+        )
+        
+        combined_Esx3 = combined_Esx3[mask]
+        combined_Eprop = combined_Eprop[mask]# * 3
+        
+        #combined_Eprop = combined_Eprop * .686
+        
         plt.figure(figsize=(8,6))
         plt.hist2d(
             combined_Esx3,
@@ -1314,16 +1469,14 @@ class MyInteractiveApp(cmd.Cmd):
         )
         plt.xlabel("SX3 Energy (MeV)")
         plt.ylabel("PCEnergy")
-        plt.xlim(0,35)
+        plt.xlim(0,30)
-        plt.ylim(0,0.6)
+        plt.ylim(0,0.45)
         plt.colorbar(label="Counts")
         plt.tight_layout()
         plt.savefig(f"{outdir}/Eprop_vs_Esx3.png", dpi=300)
         plt.show()
         
-        
+        plt.figure(figsize=(11,6), facecolor='white')
-        
-        plt.figure(figsize=(12,6), facecolor='white')
         plt.hist2d(
             combined_Esx3,
             combined_Eprop,
@@ -1334,8 +1487,8 @@ class MyInteractiveApp(cmd.Cmd):
         )
         plt.xlabel("SX3 Energy (MeV)")
         plt.ylabel("PCEnergy")
-        plt.xlim(0,35)
+        plt.xlim(0,30)
-        plt.ylim(0,0.6)
+        plt.ylim(0,0.45)
         cbar = plt.colorbar()
         cbar.set_label("Counts")
         plt.minorticks_on()
@@ -1357,63 +1510,17 @@ class MyInteractiveApp(cmd.Cmd):
                 cmap='viridis'
             )
             
-            def power_fit_and_plot(x, y, label, color=None):
+            mask = data1["Esx3"] > 0
-            
-                x = np.array(x)
-                y = np.array(y)
-            
-                mask = (
-                    np.isfinite(x) &
-                    np.isfinite(y) &
-                    (x > 0) &
-                    (y > 0)
-                )
-            
-                x = x[mask]
-                y = y[mask]
-            
-                logx = np.log(x)
-                logy = np.log(y)
-            
-                b, loga = np.polyfit(logx, logy, 1)
-            
-                a = np.exp(loga)
-            
-                y_pred = a * x**b
-            
-                ss_res = np.sum((y - y_pred)**2)
-                ss_tot = np.sum((y - np.mean(y))**2)
-            
-                r2 = 1 - (ss_res / ss_tot)
-            
-                x_fit = np.linspace(np.min(x), np.max(x), 1000)
-                y_fit = a * x_fit**b
-            
-                plt.plot(
-                    x_fit,
-                    y_fit,
-                    linewidth=3,
-                    color=color,
-                    label=f'{label} fit ($R^2$={r2:.4f})'
-                )
-            
-                print(f"\n{label} power-law fit:")
-                print(f"y = {a:.6f} * x^{b:.6f}")
-                print(f"R^2 = {r2:.6f}")
-            
-                return a, b, r2
-        
-            
             power_fit_and_plot(
-                data1["Esx3"]["Esx3" > 0],
+                data1["Esx3"][mask],
-                data1["Eprop"]["Esx3" > 0],
+                data1["Eprop"][mask] * .686,
                 label=data1["particle"],
                 color='red'
             )
-            
+            mask2 = data2["Esx3"] > 0
             power_fit_and_plot(
-                data2["Esx3"],
+                data2["Esx3"][mask2],
-                data2["Eprop"],
+                data2["Eprop"][mask2] * .686,
                 label=data2["particle"],
                 color='cyan'
             )