ELoss into MC

Armory/ANASEN_ML.py

@@ -23,7 +23,7 @@ INPUT_SCALER = "input_scaler.pkl"
 OUTPUT_SCALER = "output_scaler.pkl"
 
 # Candidate excitation energies (MeV) for nearest-state snapping
-EX_CANDIDATES = np.array([0, 0.3, 1.7, 2.4, 2.8], dtype=np.float32)
+#EX_CANDIDATES = np.array([0, 0.3, 1.7, 2.4, 2.8], dtype=np.float32)
 
 INPUT_BRANCHES = [
     "Tb",
@@ -72,8 +72,8 @@ beam = pred[:,0]
 Ex = pred[:,1]
 
 # Snap each predicted excitation to the nearest candidate energy.
-nearest_idx = np.argmin(np.abs(Ex[:, None] - EX_CANDIDATES[None, :]), axis=1)
+#nearest_idx = np.argmin(np.abs(Ex[:, None] - EX_CANDIDATES[None, :]), axis=1)
-Ex_snapped = EX_CANDIDATES[nearest_idx]
+#Ex_snapped = EX_CANDIDATES[nearest_idx]
 
 # See whether truth branches exist
 
@@ -106,15 +106,15 @@ if truth_Ex is not None:
     Ex_mae = mean_absolute_error(truth_Ex, Ex)
     Ex_rmse = np.sqrt(np.mean((truth_Ex - Ex)**2))
 
-    Ex_snap_mae = mean_absolute_error(truth_Ex, Ex_snapped)
+    #Ex_snap_mae = mean_absolute_error(truth_Ex, Ex_snapped)
-    Ex_snap_rmse = np.sqrt(np.mean((truth_Ex - Ex_snapped)**2))
+    #Ex_snap_rmse = np.sqrt(np.mean((truth_Ex - Ex_snapped)**2))
 
     print(f"Excitation MAE  : {Ex_mae:.4f} MeV")
     print(f"Excitation RMSE : {Ex_rmse:.4f} MeV")
-    print(f"Ex Snapped MAE  : {Ex_snap_mae:.4f} MeV")
+    #print(f"Ex Snapped MAE  : {Ex_snap_mae:.4f} MeV")
-    print(f"Ex Snapped RMSE : {Ex_snap_rmse:.4f} MeV")
+    #print(f"Ex Snapped RMSE : {Ex_snap_rmse:.4f} MeV")
 
-print("\nEx candidate states (MeV):", ", ".join(f"{x:.3f}" for x in EX_CANDIDATES))
+#print("\nEx candidate states (MeV):", ", ".join(f"{x:.3f}" for x in EX_CANDIDATES))
 
 # Plot Beam Energy
 plt.figure(figsize=(8,6))
@@ -153,14 +153,14 @@ plt.hist(
     linewidth=2,
     label="Predicted (raw)",
 )
-
+"""
 plt.hist(
     Ex_snapped,
     bins=250,
     histtype="step",
     linewidth=2,
     label="Predicted (snapped)",
-)
+)"""
 
 if truth_Ex is not None:
     plt.hist(
@@ -203,6 +203,7 @@ plt.xlabel("True Excitation Energy (MeV)")
 plt.ylabel("Predicted Excitation Energy (MeV)")
 plt.title("Excitation Energy Reconstruction")
 
+"""
 if truth_Ex is not None:
     plt.figure(figsize=(6,6))
 
@@ -215,4 +216,4 @@ if truth_Ex is not None:
 
     plt.xlabel("True Excitation Energy (MeV)")
     plt.ylabel("Snapped Excitation Energy (MeV)")
-    plt.title("Excitation Energy Reconstruction (Snapped)")
+    plt.title("Excitation Energy Reconstruction (Snapped)")"""

Armory/ClassPW.h

@@ -117,7 +117,7 @@ private:
   TVector3 trackVec;
 
   const int nWire = 24;
-  const int wireShift = 3;
+  const int wireShift = 4;
   //const float zLen = 380; // mm
   const float zLen = 348.6; // mm
   const float radiusA = 37;

Armory/Makefile

@@ -69,6 +69,6 @@ EventBuilder: EventBuilder.cpp ClassData.h fsuReader.h Hit.h
 	@echo "--------- making EventBuilder"
 	$(CXX) $(CXXFLAGS) EventBuilder.cpp -o EventBuilder $(LDFLAGS)
 
-AnasenMS: anasenMS.cpp constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassPW.h ClassAnasen.h
+AnasenMS: anasenMS.cpp constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassPW.h ClassAnasen.h EnergyLoss.h
 	@echo "--------- making ANASEN Monte Carlo"
 	$(CXX) $(CXXFLAGS) anasenMS.cpp -o AnasenMS $(LDFLAGS)

Armory/anasenMS.cpp

@@ -12,6 +12,10 @@
 #include "ClassTransfer.h" // Reaction kinematics and MC event generation
 #include "ClassAnasen.h"   // ANASEN detector model classes (SX3, PW, etc.)
 #include "ClassQQQ.h"      // QQQ detector model class
+#include "anasen_anode_cathode_hyperboloids.h"
+#include "EnergyLoss.h"    // energy loss lookup between two positions in a medium
+#include "HistPlotter.h"
+#include "AutoHist2D.h"    // auto-ranged, auto-binned 2D histograms written alongside tree1
 #include <stdio.h>
 #include <stdlib.h>
 #include <set>
@@ -27,8 +31,8 @@
 // calculate real and reconstructed tracks and Q-value uncertainty
 
 // Function to load energy loss table from file
-TGraph* LoadELoss(const char* filename) {
+TGraph* LoadELoss(const std::string& filename) {
-  TGraph* g = new TGraph(filename, "%lg %lg");
+  TGraph* g = new TGraph(filename.c_str(), "%lg %lg");
   return g;
 }
 
@@ -43,7 +47,7 @@ bool IsDeadCathode(int id){
 }
 
 bool IsDeadSX3(int id){ 
-    static std::set<int> dead = {}; //{0, 2, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23}; // add dead SX3 IDs here, 0-23 1,7,9,3
+    static std::set<int> dead = {};//{0, 2, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23}; // add dead SX3 IDs here, 0-23 1,7,9,3
     return dead.count(id);
 }
 
@@ -63,7 +67,7 @@ bool IsDeadSX3FrontDnChannel(int sx3ID, int chDn){
 
 bool IsDeadSX3BackChannel(int sx3ID, int chBk){
   static std::set<std::pair<int, int>> dead = {
-    //{9, 10}
+    //{1, 10}
     // {sx3ID, back-channel} 
   };
   return dead.count({sx3ID, chBk});
@@ -127,18 +131,18 @@ int main(int argc, char **argv){
   TransferReaction transfer;
 
   //To set beam energy loss, use energy loss app, and create table with target isotope, set Initial beam energy as max energy 
-  transfer.SetA(18, 10 0);            // 22Mg projectile
+  transfer.SetA(27, 13, 0);            // 22Mg projectile
-  TGraph* elossBeam = LoadELoss("../ELoss/HeLoss/E_vs_x_Ne-18.dat");
+  TGraph* elossBeam = LoadELoss("../ELoss/HeLoss/E_vs_x_Al-27.dat");
   transfer.Seta(4, 2);                 // 4He target
   transfer.Setb(1, 1);                 // outgoing proton from the primary transfer
-  transfer.SetB(21, 11);               // 30Si* heavy product
+  transfer.SetB(30, 14);               // 30Si* heavy product
   const ReactionConfig reactionConfig = transfer.GetRectionConfig();
   const double beamA = reactionConfig.beamA;                   // mass number of 14N beam
-  //const double beamE = 72 / beamA;                // beam energy in MeV
+  const double beamE = 72 / beamA;                // beam energy in MeV
 
   // Excited state lists (projectile and heavy-product excitation states)
   std::vector<float> ExAList = {0};              // Beam excited energy
-  std::vector<float> ExList = {0, 0.3, 1.7, 2.4, 2.8};             // Heavy product excited energy
+  std::vector<float> ExList = {0};             // Heavy product excited energy
 
   const int kMBeam = reactionConfig.beamA; // mass number of beam
   const int kMTarget = reactionConfig.targetA; // mass number of target
@@ -154,11 +158,23 @@ int main(int argc, char **argv){
   const int decayEjectA = 1;
   const int decayEjectZ = 1;
 
+  std::string b;
+  if (reactionConfig.recoilLightA == 1) {
+    b = "proton";
+  } else if (reactionConfig.recoilLightA == 2) {
+    b = "deuteron";
+  } else if (reactionConfig.recoilLightA == 3) {
+    b = "triton";
+  } else if (reactionConfig.recoilLightA == 4) {
+    b = "alpha";
+  }
+
+  TGraph* elossLight = LoadELoss("../ELoss/HeLoss/E_vs_x_" + b + ".dat");
   // define vertex position uniform distribution ranges (mm)
-  double vertexXRange[2] = {   -5, 5}; // mm - 5, 5
+  double vertexXRange[2] = {   0,0}; // mm - 5, 5
-  double vertexYRange[2] = {   -5, 5}; // -5, 5
+  double vertexYRange[2] = {   0,0}; // -5, 5
   double vertexZRange[2] = { -174.3, 174.3}; // -174.3, 174.3 (full length of gas volume, centered at 0)
-  const double beamEntranceZ = -280; //vertexZRange[0]; // mm, assumed beam entrance into the gas
+  const double beamEntranceZ = -280 - 174.3; //vertexZRange[0]; // mm, assumed beam entrance into the gas
 
 
   // detector resolution / uncertainty parameters
@@ -206,17 +222,11 @@ int main(int argc, char **argv){
   TString saveFileName = "SimAnasen1.root";
   printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data());
   TFile * saveFile = new TFile(saveFileName, "recreate");
-  //TFile * saveFile2 = new TFile("SimAnasen2.root", "recreate");
   TTree * tree1 = new TTree("tree1", "tree1");
-  TTree * tree2 = new TTree("tree2", "tree2");
-  //TTree * tree3 = new TTree("tree3", "tree3"); // only includes Tb and thetab
 
   // beam and CM variables saved in tree
   double KEA;
-  double KEA2;
   double beamPath_cm;
-  double beamEnergy;
-  double beamEnergyLoss;
   int MBeamOut;
   int MTargetOut;
   int MLightOut;
@@ -225,14 +235,10 @@ int main(int argc, char **argv){
   int ZTargetOut;
   int ZLightOut;
   int ZHeavyOut;
+  double beamEnergy;
   tree1->Branch("beamKEA",      &KEA, "beamKEA/D");
-  tree2->Branch("beamKEA",      &KEA2, "beamKEA/D");
   tree1->Branch("beamPath_cm",  &beamPath_cm, "beamPath_cm/D");
-  tree2->Branch("beamPath_cm",  &beamPath_cm, "beamPath_cm/D");
   tree1->Branch("beamEnergy",   &beamEnergy, "beamEnergy/D");
-  tree2->Branch("beamEnergy",   &beamEnergy, "beamEnergy/D");
-  tree1->Branch("beamEnergyLoss", &beamEnergyLoss, "beamEnergyLoss/D");
-  tree2->Branch("beamEnergyLoss", &beamEnergyLoss, "beamEnergyLoss/D");
   tree1->Branch("MBeam", &MBeamOut, "MBeam/I");
   tree1->Branch("MTarget", &MTargetOut, "MTarget/I");
   tree1->Branch("MLight", &MLightOut, "MLight/I");
@@ -241,10 +247,6 @@ int main(int argc, char **argv){
   tree1->Branch("ZTarget", &ZTargetOut, "ZTarget/I");
   tree1->Branch("ZLight", &ZLightOut, "ZLight/I");
   tree1->Branch("ZHeavy", &ZHeavyOut, "ZHeavy/I");
-  tree2->Branch("MBeam", &MBeamOut, "MBeam/I");
-  tree2->Branch("MTarget", &MTargetOut, "MTarget/I");
-  tree2->Branch("MLight", &MLightOut, "MLight/I");
-  tree2->Branch("MHeavy", &MHeavyOut, "MHeavy/I");
 
   // constant reaction mass numbers stored in every event entry
   MBeamOut = kMBeam;
@@ -253,11 +255,8 @@ int main(int argc, char **argv){
   MHeavyOut = kMHeavy;
 
   double thetaCM, phiCM;
-  double thetaCM2, phiCM2;
   tree1->Branch("thetaCM", &thetaCM, "thetaCM/D");
   tree1->Branch("phiCM",     &phiCM, "phiCM/D");
-  tree2->Branch("thetaCM", &thetaCM2, "thetaCM/D");
-  tree2->Branch("phiCM",     &phiCM2, "phiCM/D");
 
   // outgoing particles in lab frame (light/heavy)
   double thetab, phib, Tb, qqqTb;
@@ -273,21 +272,15 @@ int main(int argc, char **argv){
   tree1->Branch("qqqTb", &qqqTb, "qqqTb/D"); // kinetic energy of light particle at vertex (before energy loss) for events where the light particle hits the QQQ, currently set to 0 for simplicity
   tree1->Branch("qqqTB", &qqqTB, "qqqTB/D"); // kinetic energy of heavy particle at vertex (before energy loss) for events where the light
 
-  double thetab2, phib2, Tb2, qqqTb2;
+  double Esx3, Eqqq, Edet;
-  double thetaB2, phiB2, TB2, qqqTB2;
+  tree1->Branch("Esx3", &Esx3, "Esx3/D");
-  std::array<double, 2> T2;
+  tree1->Branch("Eqqq", &Eqqq, "Eqqq/D");
-  tree2->Branch("thetab", &thetab2, "thetab/D");
+  tree1->Branch("Edet", &Edet, "Edet/D");
-  tree2->Branch("phib",     &phib2, "phib/D");
-  tree2->Branch("Tb",         &Tb2, "Tb/D");
-  tree2->Branch("thetaB", &thetaB2, "thetaB/D");
-  tree2->Branch("phiB",     &phiB2, "phiB/D");
-  tree2->Branch("TB",         &TB2, "TB/D");
-  tree2->Branch("T",          &T2, "T/D");
-  tree2->Branch("qqqTb", &qqqTb2, "qqqTb/D");
-  tree2->Branch("qqqTB", &qqqTB2, "qqqTB/D");
 
-  //tree3->Branch("Tb", &Tb, "Tb/D");
+  double Eanode, Ecathode, EPC;
-  //tree3->Branch("thetab", &thetab, "thetab/D");
+  tree1->Branch("Eanode", &Eanode, "Eanode/D");
+  tree1->Branch("Ecathode", &Ecathode, "Ecathode/D");
+  tree1->Branch("EPC", &EPC, "EPC/D");
 
   // excitation state identifiers
   int ExAID;
@@ -295,107 +288,79 @@ int main(int argc, char **argv){
   tree1->Branch("ExAID", &ExAID, "ExAID/I"); // projectile excitation state ID
   tree1->Branch("ExA",     &ExA, "ExA/D"); // projectile excitation energy in MeV
 
-  int ExAID2;
-  double ExA2;
-  tree2->Branch("ExAID", &ExAID2, "ExAID/I");
-  tree2->Branch("ExA",     &ExA2, "ExA/D");
-
   int ExID;
   double Ex;
   tree1->Branch("ExID", &ExID, "ExID/I"); // target excitation state ID
   tree1->Branch("Ex",     &Ex, "Ex/D"); // target excitation energy in MeV
 
-  int ExID2;
-  double Ex2;
-  tree2->Branch("ExID", &ExID2, "ExID/I");
-  tree2->Branch("Ex",     &Ex2, "Ex/D");
-
   // true vertex position in target volume
-  double vertexX, vertexY, vertexZ;
+  double vertexX, vertexY, vertexZ, beamDistance;
+  tree1->Branch("beamDistance", &beamDistance, "beamDistance/D"); // distance of the beam in the target volume in mm
   tree1->Branch("vX",     &vertexX, "VertexX/D"); // true vertex X position in mm
   tree1->Branch("vY",     &vertexY, "VertexY/D"); // true vertex Y position in mm
   tree1->Branch("vZ",     &vertexZ, "VertexZ/D"); // true vertex Z position in mm
 
-  double vertexX2, vertexY2, vertexZ2;
-  tree2->Branch("vX",     &vertexX2, "VertexX/D");
-  tree2->Branch("vY",     &vertexY2, "VertexY/D");
-  tree2->Branch("vZ",     &vertexZ2, "VertexZ/D");
-
   // reconstructed SX3 hit position
   double sx3X, sx3Y, sx3Z;
   tree1->Branch("sx3X",     &sx3X, "sx3X/D"); // reconstructed X position from SX3 (with optional smearing) in mm
   tree1->Branch("sx3Y",     &sx3Y, "sx3Y/D"); // reconstructed Y position from SX3 (with optional smearing)
   tree1->Branch("sx3Z",     &sx3Z, "sx3Z/D"); // reconstructed Z position from SX3 (with optional smearing)
 
-  double sx3X2, sx3Y2, sx3Z2;
-  tree2->Branch("sx3X",     &sx3X2, "sx3X/D");
-  tree2->Branch("sx3Y",     &sx3Y2, "sx3Y/D");
-  tree2->Branch("sx3Z",     &sx3Z2, "sx3Z/D");
-
   double qqqX, qqqY, qqqZ;
   tree1->Branch("qqqX",     &qqqX, "qqqX/D"); // reconstructed X position from QQQ (with optional smearing) in mm
   tree1->Branch("qqqY",     &qqqY, "qqqY/D"); // reconstructed Y position from QQQ (with optional smearing)
   tree1->Branch("qqqZ",     &qqqZ, "qqqZ/D"); // reconstructed Z position from QQQ (with optional smearing)
 
-  double qqqX2, qqqY2, qqqZ2;
+  double detX, detY, detZ;
-  tree2->Branch("qqqX",     &qqqX2, "qqqX/D");
+  tree1->Branch("detX",     &detX, "detX/D");
-  tree2->Branch("qqqY",     &qqqY2, "qqqY/D");
+  tree1->Branch("detY",     &detY, "detY/D");
-  tree2->Branch("qqqZ",     &qqqZ2, "qqqZ/D");
+  tree1->Branch("detZ",     &detZ, "detZ/D");
+
+  double aX, aY, aZ;
+  tree1->Branch("aX",     &aX, "aX/D");
+  tree1->Branch("aY",     &aY, "aY/D");
+  tree1->Branch("aZ",     &aZ, "aZ/D");
+
+  double cX, cY, cZ;
+  tree1->Branch("cX",     &cX, "cX/D");
+  tree1->Branch("cY",     &cY, "cY/D");
+  tree1->Branch("cZ",     &cZ, "cZ/D");
+
+  double dl;
+  tree1->Branch("dl", &dl, "dl/D");
   
   // PW nearest and next nearest wires
   int anodeID[2], cathodeID[2];
-  int anodeID2[2], cathodeID2[2];
   tree1->Branch("aID", anodeID, "anodeID/I"); // anodeID[0] is nearest anode wire, anodeID[1] is next nearest anode wire
   tree1->Branch("cID", cathodeID, "cathodeID/I"); // cathodeID[0] is nearest cathode wire, cathodeID[1] is next nearest cathode wire
-  tree2->Branch("aID", anodeID2, "anodeID/I");
-  tree2->Branch("cID", cathodeID2, "cathodeID/I");
 
   // distances to nearest wires
   double anodeDist[2], cathodeDist[2];
-  double anodeDist2[2], cathodeDist2[2];
   tree1->Branch("aDist",   anodeDist, "anodeDist/D");
   tree1->Branch("cDist",  cathodeDist, "cathodeDist/D");
-  tree2->Branch("aDist",   anodeDist2, "anodeDist/D");
-  tree2->Branch("cDist",  cathodeDist2, "cathodeDist/D");
 
   // SX3 channel assignment and Z fraction (depth) information
   int sx3ID, sx3Up, sx3Dn, sx3Bk, qqqID;
   double sx3ZFrac;
-  int sx3ID2, sx3Up2, sx3Dn2, sx3Bk2, qqqID2;
-  double sx3ZFrac2;
   tree1->Branch("sx3ID",   &sx3ID, "sx3ID/I");
   tree1->Branch("qqqID",   &qqqID, "qqqID/I");
   tree1->Branch("sx3Up",   &sx3Up, "sx3Up/I");
   tree1->Branch("sx3Dn",   &sx3Dn, "sx3Dn/I");
   tree1->Branch("sx3Bk",   &sx3Bk, "sx3Bk/I");
   tree1->Branch("sx3ZFrac", &sx3ZFrac, "sx3ZFrac/D");
-  tree2->Branch("sx3ID",   &sx3ID2, "sx3ID/I");
-  tree2->Branch("qqqID",   &qqqID2, "qqqID/I");
-  tree2->Branch("sx3Up",   &sx3Up2, "sx3Up/I");
-  tree2->Branch("sx3Dn",   &sx3Dn2, "sx3Dn/I");
-  tree2->Branch("sx3Bk",   &sx3Bk2, "sx3Bk/I");
-  tree2->Branch("sx3ZFrac", &sx3ZFrac2, "sx3ZFrac/D");
 
   // reconstructed angles from PW track fit, method 1 and 2
   double reTheta, rePhi;
-  double reTheta2, rePhi2;
   tree1->Branch("reTheta", &reTheta, "reconstucted_theta/D");
   tree1->Branch("rePhi",     &rePhi, "reconstucted_phi/D");
-  tree2->Branch("reTheta", &reTheta2, "reconstucted_theta/D");
-  tree2->Branch("rePhi",     &rePhi2, "reconstucted_phi/D");
 
   double reTheta1, rePhi1;
-  double reTheta12, rePhi12;
   tree1->Branch("reTheta1", &reTheta1, "reconstucted_theta1/D");
   tree1->Branch("rePhi1",     &rePhi1, "reconstucted_phi1/D");
-  tree2->Branch("reTheta1", &reTheta12, "reconstucted_theta1/D");
-  tree2->Branch("rePhi1",     &rePhi12, "reconstucted_phi1/D");
 
   // reconstructed vertex Z from PW fit
   double z0;
-  double z02;
   tree1->Branch("z0",     &z0, "reconstucted_Z/D");
-  tree2->Branch("z0",     &z02, "reconstucted_Z/D");
 
   //========timer
   TBenchmark clock;
@@ -428,9 +393,10 @@ int main(int argc, char **argv){
 
     // compute beam energy at the event vertex from the gas path length
     beamPath_cm = TVector3(vertexZ - beamEntranceZ, vertexX, vertexY).Mag() * 0.1;
+    beamDistance = vertexZ - beamEntranceZ;
     if( beamPath_cm < 0 ) beamPath_cm = 0;
     beamEnergy = elossBeam->Eval(beamPath_cm); // MeV
-    beamEnergyLoss = elossBeam->Eval(0.0) - beamEnergy;
+    double beamEnergyLoss = elossBeam->Eval(0.0) - beamEnergy;
     KEA = beamEnergy / beamA;
     //KEA = gRandom->Uniform(0, beamE);
     transfer.SetIncidentEnergyAngle(KEA, 0, 0);
@@ -454,31 +420,6 @@ int main(int argc, char **argv){
     T[0] = Tb;
     T[1] = TB;
 
-    //secondary decay
-    TLorentzVector decayProton;
-    TLorentzVector heavy20;
-    if(enableSequentialDecay){
-      heavy20 = SimulateSequentialDecay(PB, decayDaughterA, decayDaughterZ,
-                                       decayEjectA, decayEjectZ, decayProton);
-      thetab2 = decayProton.Theta() * TMath::RadToDeg();
-      phib2 = decayProton.Phi() * TMath::RadToDeg();
-      Tb2 = decayProton.E() - decayProton.M();
-      thetaB2 = heavy20.Theta() * TMath::RadToDeg();
-      phiB2 = heavy20.Phi() * TMath::RadToDeg();
-      TB2 = heavy20.E() - heavy20.M();
-      T2[0] = Tb2;
-      T2[1] = TB2;
-    } else {
-      thetab2 = TMath::QuietNaN();
-      phib2  = TMath::QuietNaN();
-      Tb2    = TMath::QuietNaN();
-      thetaB2 = TMath::QuietNaN();
-      phiB2 = TMath::QuietNaN();
-      TB2 = TMath::QuietNaN();
-      T2[0] = TMath::QuietNaN();
-      T2[1] = TMath::QuietNaN();
-    }
-
     delete [] output;
 
     // set direction vector from lab angle
@@ -547,101 +488,63 @@ int main(int argc, char **argv){
       rePhi1   = pw->GetTrackPhi() * TMath::RadToDeg();
 
       z0 = pw->GetZ0();
+
+      pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false);
+
+      const TVector3 trackDir = pw->GetTrackVec();
+      const double transverseDirection2 =
+          trackDir.X() * trackDir.X() + trackDir.Y() * trackDir.Y();
+
+      if (transverseDirection2 > 0.0) {
+        const double pathToRhoMin =
+            -(hitPos.X() * trackDir.X() + hitPos.Y() * trackDir.Y())
+            / transverseDirection2;
+
+        //const TVector3 rhoMin = hitPos + pathToRhoMin * trackDir;
+        const TVector3 rhoMin = vertex;
+
+        // Original two-wire-layer model, kept for reference:
+        auto [cathodeIntersection, anodeIntersection, pcPathLengthCm] =
+            find_PC_PathLength(rhoMin, hitPos);
+        
+        // The helper returns cathode first, then anode.
+        aX = anodeIntersection.X();
+        aY = anodeIntersection.Y();
+        aZ = anodeIntersection.Z();
+        
+        cX = cathodeIntersection.X();
+        cY = cathodeIntersection.Y();
+        cZ = cathodeIntersection.Z();
+      }
+
+      //Energy loss calculations
+
+      double distance_sx3;
+      Esx3 = CalculateEnergyLoss(vertexX, vertexY, vertexZ,
+               sx3X, sx3Y, sx3Z,
+               b, "He", Tb,
+           distance_sx3);
+      dl = distance_sx3;
+
+      double distance_A;
+      Eanode = CalculateEnergyLoss(vertexX, vertexY, vertexZ,
+               aX, aY, aZ,
+               b, "He", Tb,
+           distance_A);
+
+      double distance_C;
+      Ecathode = CalculateEnergyLoss(vertexX, vertexY, vertexZ,
+               cX, cY, cZ,
+               b, "He", Tb,
+           distance_C);
+
+      EPC = Eanode - Ecathode;
+
+      if (Esx3 <= 0 || Eanode <= 0 || Ecathode <= 0) {
+        Esx3 = NAN;
+      }
+      Edet = Esx3;
       tree1->Fill();
-      //tree3->Fill();
-
-      // fill tree2 using the secondary proton from 21Na* decay
-      TVector3 dir2(1, 0, 0);
-      dir2.SetTheta(thetab2 * TMath::DegToRad());
-      dir2.SetPhi(phib2 * TMath::DegToRad());
-
-      pw->FindWireID(vertex, dir2, false);
-      sx3->FindSX3Pos(vertex, dir2, false);
-      PWHitInfo hitInfo2 = pw->GetHitInfo();
-
-      anodeID2[0]   = hitInfo2.nearestWire.first;
-      cathodeID2[0] = hitInfo2.nearestWire.second;
-      anodeID2[1]   = hitInfo2.nextNearestWire.first;
-      cathodeID2[1] = hitInfo2.nextNearestWire.second;
-
-      anodeDist2[1]   = hitInfo2.nextNearestDist.first;
-      cathodeDist2[1] = hitInfo2.nextNearestDist.second;
-
-      if(IsDeadAnode(anodeID2[0]) || IsDeadCathode(cathodeID2[0])){
-        sx3ID2 = -1;
-      } else {
-        sx3ID2 = sx3->GetID();
-      }
-
-      if(sx3ID2 < 0 || IsDeadSX3(sx3ID2)){
-        sx3ID2 = -1;
-        sx3Up2 = -1;
-        sx3Dn2 = -1;
-        sx3Bk2 = -1;
-        sx3ZFrac2 = TMath::QuietNaN();
-        sx3X2 = TMath::QuietNaN();
-        sx3Y2 = TMath::QuietNaN();
-        sx3Z2 = TMath::QuietNaN();
-        anodeDist2[0] = TMath::QuietNaN();
-        cathodeDist2[0] = TMath::QuietNaN();
-        reTheta2 = TMath::QuietNaN();
-        rePhi2 = TMath::QuietNaN();
-        reTheta12 = TMath::QuietNaN();
-        rePhi12 = TMath::QuietNaN();
-        z02 = TMath::QuietNaN();
-      } else {
-        anodeDist2[0]   = hitInfo2.nearestDist.first;
-        cathodeDist2[0] = hitInfo2.nearestDist.second;
-        sx3Up2    = sx3->GetChUp();
-        sx3Dn2    = sx3->GetChDn();
-        sx3Bk2    = sx3->GetChBk();
-        if(IsDeadSX3ChannelCombo(sx3ID2, sx3Up2, sx3Dn2, sx3Bk2)){
-          sx3ID2 = -1;
-          sx3Up2 = -1;
-          sx3Dn2 = -1;
-          sx3Bk2 = -1;
-          sx3ZFrac2 = TMath::QuietNaN();
-          sx3X2 = TMath::QuietNaN();
-          sx3Y2 = TMath::QuietNaN();
-          sx3Z2 = TMath::QuietNaN();
-          anodeDist2[0] = TMath::QuietNaN();
-          cathodeDist2[0] = TMath::QuietNaN();
-          reTheta2 = TMath::QuietNaN();
-          rePhi2 = TMath::QuietNaN();
-          reTheta12 = TMath::QuietNaN();
-          rePhi12 = TMath::QuietNaN();
-          z02 = TMath::QuietNaN();
-        } else {
-        sx3ZFrac2 = sx3->GetZFrac();
-        TVector3 hitPos2 = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L);
-        sx3X2 = hitPos2.X();
-        sx3Y2 = hitPos2.Y();
-        sx3Z2 = hitPos2.Z();
-        pw->CalTrack(hitPos2, anodeID2[0], cathodeID2[0], false);
-        reTheta2 = pw->GetTrackTheta() * TMath::RadToDeg();
-        rePhi2   = pw->GetTrackPhi() * TMath::RadToDeg();
-        pw->CalTrack2(hitPos2, hitInfo2, sigmaPW_A, sigmaPW_C, false);
-        reTheta12 = pw->GetTrackTheta() * TMath::RadToDeg();
-        rePhi12   = pw->GetTrackPhi() * TMath::RadToDeg();
-        z02       = pw->GetZ0();
-        }
-      }
-
-      KEA2 = KEA;
-      thetaCM2 = thetaCM;
-      phiCM2 = phiCM;
-      ExAID2 = ExAID;
-      ExA2 = ExA;
-      ExID2 = ExID;
-      Ex2 = Ex;
-      vertexX2 = vertexX;
-      vertexY2 = vertexY;
-      vertexZ2 = vertexZ;
-      qqqX = TMath::QuietNaN();
-      qqqY = TMath::QuietNaN();
-      qqqZ = TMath::QuietNaN();
-
-      tree2->Fill();
 
     }else if (qqqID >= 0){
       // handle QQQ hit case
@@ -672,6 +575,15 @@ int main(int argc, char **argv){
       qqqY = hitPos.Y();
       qqqZ = hitPos.Z();
 
+        auto [cathodeIntersection, anodeIntersection, pcPathLengthCm] =
+          find_PC_PathLength(vertex, hitPos);
+        aX = anodeIntersection.X();
+        aY = anodeIntersection.Y();
+        aZ = anodeIntersection.Z();
+        cX = cathodeIntersection.X();
+        cY = cathodeIntersection.Y();
+        cZ = cathodeIntersection.Z();
+
       if(enableVis){
         visTrackVertex.push_back(vertex);
         visTrackDir.push_back(dir);
@@ -679,47 +591,34 @@ int main(int argc, char **argv){
         //visTrackWires.push_back({anodeID[0], cathodeID[0]});
       }
 
-      tree1->Fill();
+      //Energy loss calculations
-      //tree3->Fill();
 
-      TVector3 dir2(1, 0, 0);
+      double distance_qqq;
-      dir2.SetTheta(thetab2 * TMath::DegToRad());
+      Eqqq = CalculateEnergyLoss(vertexX, vertexY, vertexZ,
-      dir2.SetPhi(phib2 * TMath::DegToRad());
+                           qqqX, qqqY, qqqZ,
+                           b, "He",
+               qqqTb, distance_qqq);
+      dl = distance_qqq;
+      double distance_A;
+      Eanode = CalculateEnergyLoss(vertexX, vertexY, vertexZ,
+                           aX, aY, aZ,
+                           b, "He",
+               qqqTb, distance_A);
+      double distance_C;
+      Ecathode = CalculateEnergyLoss(vertexX, vertexY, vertexZ,
+                           cX, cY, cZ,
+                           b, "He",
+               qqqTb, distance_C);
 
-      qqq->FindQQQPos(vertex, dir2, false);
+      if (Eqqq <= 0 || Eanode <= 0 || Ecathode <= 0) {
-
+          Eqqq = NAN;
-      if(qqqID2 < 0){
-        qqqID2 = -1;
-        qqqX2 = TMath::QuietNaN();
-        qqqY2 = TMath::QuietNaN();
-        qqqZ2 = TMath::QuietNaN();
-        anodeDist2[0] = TMath::QuietNaN();
-        cathodeDist2[0] = TMath::QuietNaN();
-        reTheta2 = TMath::QuietNaN();
-        rePhi2 = TMath::QuietNaN();
-        reTheta12 = TMath::QuietNaN();
-        rePhi12 = TMath::QuietNaN();
-        z02 = TMath::QuietNaN();
-        anodeID2[0]   = TMath::QuietNaN(); // no valid anode wire for QQQ hit case
-        cathodeID2[0] = TMath::QuietNaN(); // no valid cathode wire for QQQ hit case
-        anodeID2[1]   = TMath::QuietNaN(); // no valid next nearest anode wire for QQQ hit case
-        cathodeID2[1] = TMath::QuietNaN(); // no valid next nearest cathode wire for QQQ hit case
-        anodeDist2[1]   = TMath::QuietNaN();
-        cathodeDist2[1] = TMath::QuietNaN();
       }
+      Edet = Eqqq;
+      EPC = Eanode - Ecathode;
 
-      KEA2 = KEA;
+      tree1->Fill();
-      thetaCM2 = thetaCM;
+      AutoHist2D::Fill("beamEnergy_vs_vZ", beamEnergy, vertexZ, "beamEnergy (MeV)", "vZ (mm)");
-      phiCM2 = phiCM;
-      ExAID2 = ExAID;
-      ExA2 = ExA;
-      ExID2 = ExID;
-      Ex2 = Ex;
-      vertexX2 = vertexX;
-      vertexY2 = vertexY;
-      vertexZ2 = vertexZ;
     
-      tree2->Fill();
     }else{
       // no valid SX3 hit: mark clearly invalid
       sx3Up = -1;
@@ -763,15 +662,13 @@ int main(int argc, char **argv){
   }
 
   // write results to ROOT file and close
+  AutoHist2D::WriteAll(); // books, fills and writes all registered 2D histograms into saveFile
   tree1->Write("", TObject::kOverwrite);
-  tree2->Write("", TObject::kOverwrite);
-  //tree3->Write("", TObject::kOverwrite);
   int count1 = tree1->GetEntries();
-  int count2 = tree2->GetEntries();
   //int count3 = tree3->GetEntries();
   saveFile->Close();
 
-  printf("=============== done. saved as %s. tree1 entries: %d, tree2 entries: %d\n", saveFileName.Data(), count1, count2);
+  printf("=============== done. saved as %s. tree1 entries: %d\n", saveFileName.Data(), count1);
 
   if(enableVis){ // to enable visualization, run with 3rd argument "vis", e.g. "./anasenMC 1000 vis"
     printf("Displaying geometry with %zu tracks from simulation\n", visTrackVertex.size());

ELoss/EXtable.py

@@ -35,6 +35,35 @@ material_def = [
 gas_mix = catima.Material(material_def)
 gas_mix.density(rho_g_cm3)
 
+# MATERIAL BANK - Additional Materials for Energy Loss Calculations
+
+# Kapton (C22H10N2O5)
+# Density: 1.42 g/cm3
+# Molecular weight: 22*12 + 10*1 + 2*14 + 5*16 = 264 + 10 + 28 + 80 = 382 g/mol
+m_h  = 1.0078
+m_n  = 14.0067
+kapton_molar_mass = 22*m_c + 10*m_h + 2*m_n + 5*m_o  # ~382 g/mol
+kapton_material_def = [
+    (m_c,  6, 22/382 * kapton_molar_mass / m_c),
+    (m_h,  1, 10/382 * kapton_molar_mass / m_h),
+    (m_n,  7, 2/382 * kapton_molar_mass / m_n),
+    (m_o,  8, 5/382 * kapton_molar_mass / m_o)
+]
+kapton = catima.Material(kapton_material_def)
+kapton.density(1.42)  # g/cm3
+
+# Mylar (C10H8O4, polyethylene terephthalate)
+# Density: 1.39 g/cm3
+# Molecular weight: 10*12 + 8*1 + 4*16 = 120 + 8 + 64 = 192 g/mol
+mylar_molar_mass = 10*m_c + 8*m_h + 4*m_o  # ~192 g/mol
+mylar_material_def = [
+    (m_c,  6, 10/192 * mylar_molar_mass / m_c),
+    (m_h,  1, 8/192 * mylar_molar_mass / m_h),
+    (m_o,  8, 4/192 * mylar_molar_mass / m_o)
+]
+mylar = catima.Material(mylar_material_def)
+mylar.density(1.39)  # g/cm3
+
 # FUNCTION
 
 def make_E_vs_x(
@@ -42,9 +71,13 @@ def make_E_vs_x(
     mass_u,
     emax_mev,
     label,
-    npoints=500
+    npoints=500,
+    material=None
 ):
 
+    if material is None:
+        material = gas_mix
+
     projectile = catima.Projectile(mass_u, z)
 
     # Energy grid
@@ -58,7 +91,7 @@ def make_E_vs_x(
         projectile.T(energy / mass_u)
 
         # MeV / (g/cm^2)
-        S_mass[i] = catima.dedx(projectile, gas_mix)
+        S_mass[i] = catima.dedx(projectile, material)
 
     # Convert to MeV/cm
     S_linear = S_mass * rho_g_cm3
@@ -115,6 +148,42 @@ x, E = make_E_vs_x(
     label="alpha"
 )
 
+# Generate tables for kapton
+print("\n=== Generating Kapton Energy Loss Tables ===")
+x, E = make_E_vs_x(
+    z=1,
+    mass_u=1.0078,
+    emax_mev=20,
+    label="proton_kapton",
+    material=kapton
+)
+
+x, E = make_E_vs_x(
+    z=2,
+    mass_u=4.0026,
+    emax_mev=40,
+    label="alpha_kapton",
+    material=kapton
+)
+
+# Generate tables for mylar
+print("\n=== Generating Mylar Energy Loss Tables ===")
+x, E = make_E_vs_x(
+    z=1,
+    mass_u=1.0078,
+    emax_mev=20,
+    label="proton_mylar",
+    material=mylar
+)
+
+x, E = make_E_vs_x(
+    z=2,
+    mass_u=4.0026,
+    emax_mev=40,
+    label="alpha_mylar",
+    material=mylar
+)
+
 # PLOT
 
 plt.figure(figsize=(8,6))

ELoss/PCEnergyAnalysis.py

@@ -23,6 +23,7 @@ import re
 from matplotlib.colors import LinearSegmentedColormap
 from mpl_toolkits.mplot3d import Axes3D
 import nbformat as nbf
+import shutil
 
 # ROOT-like styling
 plt.rcParams.update({
@@ -54,7 +55,6 @@ plt.rcParams.update({
 
 base_cmap = plt.cm.jet
 
-# Stop before the red region
 colors = base_cmap(np.linspace(0.2, 0.65, 256))
 
 yellow_jet = LinearSegmentedColormap.from_list(
@@ -78,8 +78,6 @@ SX3_SI_THICKNESS_CM = 0.1
 QQQ_SI_THICKNESS_CM = 0.1
 SX3_THETA_MIN_DEG = 0.0
 
-# Optional dead-channel masks for SX3.
-# Format for channel masks is: (sx3ID, channel)
 DEAD_SX3_IDS = set() #Ex {9}
 DEAD_SX3_FRONT_UP = set() # Ex {(9, 3)}
 DEAD_SX3_FRONT_DN = set() # Ex {(9, 3)}
@@ -176,7 +174,7 @@ def make_E_vs_x(
         material_def = [(m_he, 2, 1.0)]
 
         #m_mix_avg = 0.96 * m_he + 0.04 * (m_c + 2 * m_o)
-        m_mix_avg = 1.0 * m_he
+        m_mix_avg = 0.97 * m_he + 0.03 * (m_c + 2 * m_o)
         rho_g_cm3 = (molar_density * m_mix_avg) / 1e6
         
         gas = catima.Material(material_def)
@@ -190,6 +188,44 @@ def make_E_vs_x(
         gas = catima.Material([(m_si, 14, 1.0)])
         gas.density(rho_g_cm3)
 
+    elif medium == "kapton":
+        # Kapton (C22H10N2O5)
+        # Density: 1.42 g/cm3
+        m_h = 1.0078
+        m_n = 14.0067
+        m_c = 12.0000
+        m_o = 15.9949
+        kapton_molar_mass = 22*m_c + 10*m_h + 2*m_n + 5*m_o
+        
+        material_def = [
+            (m_c, 6, 22/kapton_molar_mass * kapton_molar_mass / m_c),
+            (m_h, 1, 10/kapton_molar_mass * kapton_molar_mass / m_h),
+            (m_n, 7, 2/kapton_molar_mass * kapton_molar_mass / m_n),
+            (m_o, 8, 5/kapton_molar_mass * kapton_molar_mass / m_o)
+        ]
+        
+        rho_g_cm3 = 1.42
+        gas = catima.Material(material_def)
+        gas.density(rho_g_cm3)
+
+    elif medium == "mylar":
+        # Mylar (C10H8O4, polyethylene terephthalate)
+        # Density: 1.39 g/cm3
+        m_h = 1.0078
+        m_c = 12.0000
+        m_o = 15.9949
+        mylar_molar_mass = 10*m_c + 8*m_h + 4*m_o
+        
+        material_def = [
+            (m_c, 6, 10/mylar_molar_mass * mylar_molar_mass / m_c),
+            (m_h, 1, 8/mylar_molar_mass * mylar_molar_mass / m_h),
+            (m_o, 8, 4/mylar_molar_mass * mylar_molar_mass / m_o)
+        ]
+        
+        rho_g_cm3 = 1.39
+        gas = catima.Material(material_def)
+        gas.density(rho_g_cm3)
+
     else:
         raise ValueError("Unsupported medium")
 
@@ -199,7 +235,6 @@ def make_E_vs_x(
     E = np.linspace(0.1, emax_mev, npoints)
     S_mass = np.zeros_like(E)
     for i, energy in enumerate(E):
-        # CATIMA expects projectile energy in MeV/u, while the table axis stays in total MeV.
         projectile.T(energy / mass_u)
         S_mass[i] = catima.dedx(projectile, gas)  # MeV/(g/cm^2)
 
@@ -645,6 +680,8 @@ def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
     EfinalQ = energy_loss(particle, "Si", Eqqq, qqq_silicon_path_cm)
     Edet = (Esx3 - Efinal)
     EdetQ = Eqqq - EfinalQ
+    sx3_reached_mask = np.isfinite(Esx3) & (Esx3 > 0)
+    qqq_reached_mask = np.isfinite(Eqqq) & (Eqqq > 0)
 
     #Edet = np.where(Prange > lsx3, Efinal + Esx3, Esx3)
 
@@ -681,10 +718,13 @@ def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
         "lsx3": sx3_silicon_path_cm,
         "lqqq": qqq_silicon_path_cm,
         "EdetQ": EdetQ,
-        "beamEnergy": beamEnergy,
+        "beamEnergy": beamEnergy[sx3_event_mask],
-        "vZ": vertex_z,
+        "vZ": vertex_z[sx3_event_mask],
-        "beamEnergyLoss": data["beamEnergyLoss"],
+        "Ex": data["Ex"][sx3_event_mask],
-        "Ex": data["Ex"],
+        "sx3Reached": sx3_reached_mask,
+        "beamEnergyQqq": beamEnergy[qqq_event_mask],
+        "vZQqq": vertex_z[qqq_event_mask],
+        "qqqReached": qqq_reached_mask,
         "vZsx3": sx3_vertex_z
     }
 
@@ -1125,197 +1165,6 @@ class MyInteractiveApp(cmd.Cmd):
         look up 'uproot' for more details"""
         self.run_command_line()
         
-    def do_energy_analysis(self, arg): #geometry needs update
-    
-        args = shlex.split(arg)
-    
-        try:
-            particle = args[0]
-        except IndexError:
-            print("Please indicate reactant for analysis")
-            return
-    
-        try:
-            max_events = int(args[1])
-        except IndexError:
-            max_events = None
-    
-        if self.tree is None:
-            self.do_set_tree("")
-            print(f"Using TTree: {self.tree}")
-            if self.tree is None:
-                print("No tree is loaded. Use 'uproot file <filename>' then 'set tree tree1' or 'set tree tree2'.")
-                return
-    
-        treename = self.tree.name if hasattr(self.tree, 'name') else "tree1"
-
-        branches = [
-            "Tb",
-            "thetab",
-            "sx3Z",
-            "vX",
-            "vY",
-            "vZ",
-            "sx3X",
-            "sx3Y",
-            "sx3ID",
-            "sx3Up",
-            "sx3Dn",
-            "sx3Bk",
-            "sx3XExit",
-            "sx3YExit",
-            "sx3ZExit",
-        ]
-        initial_energy_key = "Tb"
-        polar_angle_key = "thetab"
-
-    
-        if max_events:
-            n_events = max_events
-        else:
-            n_events = self.tree.num_entries
-    
-        print(f"Loading {n_events} events from {treename}...")
-    
-        data = self.tree.arrays(
-            branches,
-            library="np",
-            entry_stop=max_events
-        )
-        initial_energy = data[initial_energy_key]
-        global sx3Z
-        sx3Z = data["sx3Z"]
-        sx3XExit = data["sx3XExit"] if "sx3XExit" in data else None
-        sx3YExit = data["sx3YExit"] if "sx3YExit" in data else None
-        sx3ZExit = data["sx3ZExit"] if "sx3ZExit" in data else None
-    
-        polar_angle_rad = np.radians(data[polar_angle_key])
-        vertex_x = data["vX"]
-        vertex_y = data["vY"]
-        vertex_z = data["vZ"]
-        
-        if treename == 'tree1':
-            sx3_hit_x = data["sx3X"]
-            sx3_hit_y = data["sx3Y"]
-            sx3_live_mask = build_sx3_live_mask(
-                data["sx3ID"],
-                data["sx3Up"],
-                data["sx3Dn"],
-                data["sx3Bk"],
-            )
-            event_mask = ~np.isnan(sx3_hit_x) & ~np.isnan(sx3_hit_y) & ~np.isnan(sx3Z) & sx3_live_mask
-        else:
-            event_mask = ~np.isnan(initial_energy) & ~np.isnan(polar_angle_rad)
-    
-        initial_energy = initial_energy[event_mask]
-        polar_angle_rad = polar_angle_rad[event_mask]
-        sx3Z = sx3Z[event_mask]
-        vertex_z_masked = vertex_z[event_mask]
-        
-        sin_theta = np.sin(polar_angle_rad)
-        sin_theta = np.where(sin_theta != 0, sin_theta, 1e-10)
-    
-        if treename == 'tree1':
-            vertex_x_masked = vertex_x[event_mask]
-            vertex_y_masked = vertex_y[event_mask]
-            sx3_hit_x_masked = sx3_hit_x[event_mask]
-            sx3_hit_y_masked = sx3_hit_y[event_mask]
-            sx3_track_distance_cm = calculate_distance_tree1(
-                vertex_x_masked,
-                vertex_y_masked,
-                vertex_z_masked,
-                sx3_hit_x_masked,
-                sx3_hit_y_masked,
-                sx3Z,
-            )
-            sx3_silicon_path_cm = calculate_sx3_thickness_path(
-                    vertex_x_masked,
-                    vertex_y_masked,
-                    vertex_z_masked,
-                    sx3_hit_x_masked,
-                    sx3_hit_y_masked,
-                    sx3Z,
-                    SX3_SI_THICKNESS_CM,)
-        else:
-            sx3_track_distance_cm = calculate_distance_tree2(vertex_z_masked, polar_angle_rad, z_max=34.86)
-            sx3_silicon_path_cm = np.full_like(sx3_track_distance_cm, np.nan, dtype=float)
-        
-        radii = np.array([3.7, 4.2])
-        anode_path_cm = radii[0] / sin_theta
-        cathode_path_cm = radii[1] / sin_theta
-    
-        print("Calculating energy losses...")
-        global EA
-        EA = energy_loss(particle, "He", initial_energy, anode_path_cm)
-        global EC
-        EC = energy_loss(particle, "He", initial_energy, cathode_path_cm)
-        global Esx3
-        Esx3 = energy_loss(particle, "He", initial_energy, sx3_track_distance_cm)
-        global Eprop
-        Eprop = EA - EC
-        global Elost
-        Elost = initial_energy - Esx3
-
-        print("Analysis complete")
-        print(f"Processed events: {len(initial_energy)}")
-        print(f"Anode average energy: {np.mean(EA):.3f} MeV")
-        print(f"Cathode average energy: {np.mean(EC):.3f} MeV")
-        print(f"sx3 average energy: {np.mean(Esx3):.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")
-        print(f"Minimum total energy loss to sx3: {np.min(Elost):.3f} MeV")
-        print(f"Proportion counter average energy difference: {np.mean(Eprop):.3f} MeV")
-        print(f"Maximum proportion counter energy difference: {np.max(Eprop):.3f} MeV")
-        print(f"Minimum proportion counter energy difference: {np.min(Eprop):.3f} MeV")
-
-
-        output_filename = "energy_analysis.root"
-
-        print(f"Writing new tree to {output_filename}")
-
-        # Load ALL original branches
-        all_data = self.tree.arrays(library="np",entry_stop=max_events)
-
-        # Create full-length arrays initialized to NaN
-        n_total = len(data["Tb"])
-
-        EA_full = np.full(n_total, np.nan)
-        EC_full = np.full(n_total, np.nan)
-        Esx3_full = np.full(n_total, np.nan)
-        Eprop_full = np.full(n_total, np.nan)
-        Elost_full = np.full(n_total, np.nan)
-
-        # Put values back into valid entries
-        EA_full[event_mask] = EA
-        EC_full[event_mask] = EC
-        Esx3_full[event_mask] = Esx3
-        Eprop_full[event_mask] = Eprop
-        Elost_full[event_mask] = Elost
-
-        # Add new branches
-        all_data["EA"] = EA_full
-        all_data["EC"] = EC_full
-        all_data["Esx3"] = Esx3_full
-        all_data["Eprop"] = Eprop_full
-        all_data["Elost"] = Elost_full
-        if treename == 'tree1':
-            lsx3_full = np.full(n_total, np.nan)
-            lsx3_full[event_mask] = sx3_silicon_path_cm
-            all_data["lsx3"] = lsx3_full
-        # Write new ROOT file as a classic TTree
-        with uproot.recreate(output_filename) as fout:
-
-            branch_types = {
-                name: array.dtype
-                for name, array in all_data.items()
-            }
-
-            fout.mktree("tree", branch_types)
-
-            fout["tree"].extend(all_data)
-
-        print("Finished writing augmented ROOT file")   
-        
     def do_make_plots(self, arg):
 
         import os
@@ -1362,9 +1211,13 @@ class MyInteractiveApp(cmd.Cmd):
         EdetQ = data["EdetQ"]
         beamEnergy = data["beamEnergy"]
         vZ = data["vZ"]
-        beamEnergyLoss = data["beamEnergyLoss"]
+        #beamEnergyLoss = data["beamEnergyLoss"]
         Ex = data["Ex"]
         vZsx3 = data["vZsx3"]
+        sx3_reached_mask = data["sx3Reached"]
+        beamEnergyQqq = data["beamEnergyQqq"]
+        vZQqq = data["vZQqq"]
+        qqq_reached_mask = data["qqqReached"]
 
         update_plot_data(f"{particle}_{treename}_Ei", Ei)
         update_plot_data(f"{particle}_{treename}_sx3Z", sx3Z)
@@ -1376,10 +1229,57 @@ class MyInteractiveApp(cmd.Cmd):
         update_plot_data(f"{particle}_treename_Edert", Edet)
 
         base = f"{particle}_{treename}_plots"
+        
+        if os.path.exists(base):
+            shutil.rmtree(base)
+        
+        # Create a fresh, empty plot folder
         os.makedirs(base, exist_ok=True)
 
         print(f"Saving plots to folder: {base} ({treename})")
 
+        # Keep only events with positive residual energy at the SX3 entrance.
+        # Events stopping in the gas never reach the detector and are excluded.
+        finite_sx3_observables = (
+            np.isfinite(beamEnergy)
+            & np.isfinite(vZ)
+            & np.isfinite(thetab)
+        )
+        finite_qqq_observables = (
+            np.isfinite(beamEnergyQqq)
+            & np.isfinite(vZQqq)
+            & np.isfinite(thetabqqq)
+        )
+        sx3_reached_hist_mask = sx3_reached_mask & finite_sx3_observables
+        qqq_reached_hist_mask = qqq_reached_mask & finite_qqq_observables
+        print(
+            f"Detector-reachable events: SX3 {np.sum(sx3_reached_hist_mask)} / "
+            f"{len(sx3_reached_hist_mask)}, QQQ {np.sum(qqq_reached_hist_mask)} / "
+            f"{len(qqq_reached_hist_mask)}"
+        )
+
+        filtered_histograms = (
+            ("SX3", "beamEnergy", beamEnergy, sx3_reached_hist_mask, "Beam Energy (MeV)"),
+            ("SX3", "vZ", vZ, sx3_reached_hist_mask, "Vertex Z (mm)"),
+            ("SX3", "thetab", thetab, sx3_reached_hist_mask, "thetab (deg)"),
+            ("QQQ", "beamEnergy", beamEnergyQqq, qqq_reached_hist_mask, "Beam Energy (MeV)"),
+            ("QQQ", "vZ", vZQqq, qqq_reached_hist_mask, "Vertex Z (mm)"),
+            ("QQQ", "thetab", thetabqqq, qqq_reached_hist_mask, "thetab (deg)"),
+        )
+        for detector_name, name, values, mask, xlabel in filtered_histograms:
+            plt.figure(figsize=(7, 5))
+            plt.hist(values[mask], bins=100, histtype="stepfilled")
+            plt.xlabel(xlabel)
+            plt.ylabel("Counts")
+            plt.title(
+                f"{particle} ({treename}) {detector_name} {name} "
+                "for detector-reachable events"
+            )
+            plt.grid(True)
+            plt.tight_layout()
+            plt.savefig(f"{base}/{detector_name}_{name}_reachable_hist.png", dpi=300)
+            plt.close()
+        
         # --- clean data ---
         x = np.asarray(sx3X, dtype=float)
         y = np.asarray(sx3Y, dtype=float)
@@ -1446,7 +1346,7 @@ class MyInteractiveApp(cmd.Cmd):
             
         mask1 = ~np.isnan(Ei) & ~np.isnan(thetab) & sx3_theta_plot_mask
         plt.figure(figsize=(7,6))
-        plt.hist2d(np.radians(thetab[mask1]), Ei[mask1], bins=200)
+        plt.hist2d(thetab[mask1], Ei[mask1], bins=200)
         plt.xlabel("thetab")
         plt.ylabel("Event Energy (MeV)")
         plt.title(f"{particle} ({treename}) Energy vs Theta")
@@ -1788,17 +1688,6 @@ class MyInteractiveApp(cmd.Cmd):
         plt.savefig(f"{base}/EBeam_vs_Z", dpi=300)
         plt.show()
         
-        mask1 = (beamEnergyLoss > 0) & ~np.isnan(vZ) & ~np.isnan(beamEnergy)
-        plt.figure(figsize=(7,6))
-        plt.hist2d(vZ[mask1], beamEnergyLoss[mask1], bins=200)
-        plt.ylabel("Beam Energy Loss")
-        plt.xlabel("Z")
-        plt.title(f"{particle} ({treename}) Beam Energy vs. Z")
-        plt.colorbar(label="Counts")
-        plt.tight_layout()
-        plt.savefig(f"{base}/EBeamLoss_vs_Z", dpi=300)
-        plt.show()
-        
         plt.figure(figsize=(7,6))
         plt.hist2d(vZ, Ex, bins=200)
         plt.ylabel("Excitation Energy")
@@ -1834,45 +1723,64 @@ class MyInteractiveApp(cmd.Cmd):
         
         if branch_names and True:
             print(f"Creating histograms for {len(branch_names)} branches...")
-            all_branches = self.tree.arrays(branch_names, library="np", entry_stop=max_events)
+            #all_branches = self.tree.arrays(branch_names, library="np", entry_stop=max_events)
 
-            for branch in branch_names:
+            for branch in self.tree.keys():
-                values = all_branches[branch]
+            
+                branch = branch.decode() if isinstance(branch, bytes) else str(branch)
             
                 try:
+                    values = self.tree[branch].array(
+                        library="np",
+                        entry_stop=max_events
+                    )
+            
                     values = np.asarray(values, dtype=float)
-                except Exception:
+            
-                    print(f"Skipping non-numeric branch: {branch}")
+                except Exception as e:
+                    print(f"Skipping branch {branch}: {e}")
                     continue
             
                 if values.ndim != 1:
                     print(f"Skipping non-scalar branch: {branch}")
                     continue
+                if branch == "qqqZ":
+                    continue
             
-                # Apply the thetab cut
+                values = values[np.isfinite(values)]
-                #values = values[mask]
-        
-                values = values[~np.isnan(values)]
             
                 if values.size == 0:
                     print(f"Skipping empty branch: {branch}")
                     continue
             
-                plt.figure(figsize=(7,5))
+                plt.figure(figsize=(7, 5))
-                try:
-                    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.hist(
-                    plt.savefig(f"{base}/{safe_name}_hist.png", dpi=300)
+                    values,
-                    plt.close()
+                    bins=100,
-                except:
+                    histtype="stepfilled"
-                    print(f"Can not print branch {branch}")
+                )
-                    continue
+            
+                plt.xlabel(branch)
+                plt.ylabel("Counts")
+                plt.title(f"{particle} {branch} distribution")
+                plt.grid(True)
+                plt.tight_layout()
+            
+                safe_name = re.sub(r"[^0-9A-Za-z_-]", "_", branch)
+                
+                filename = f"{base}/{safe_name}_hist.png"
+                
+                counter = 1
+                while os.path.exists(filename):
+                    filename = f"{base}/{safe_name}_{counter}_hist.png"
+                    counter += 1
+                
+                plt.savefig(filename, dpi=300)
+                plt.close()
+                #print(f"Saved: {filename}")
+            
+                #plt.show()
 
         else:
             print("No branches found to histogram.")
@@ -1930,7 +1838,8 @@ class MyInteractiveApp(cmd.Cmd):
                 plt.tight_layout()
                 safe_name = re.sub(r"[^0-9A-Za-z_-]", "_", branch)
                 plt.savefig(f"{base}/{safe_name}_hist.png", dpi=300)
-                plt.show()
+                #plt.show()
+                plt.close()
         else:
             print("No branches found to histogram.")