From fe112a60047e37b62a0d9f06f014057a16c90367 Mon Sep 17 00:00:00 2001
From: James Szalkie <jamesszalkie122@gmail.com>
Date: Mon, 20 Apr 2026 13:37:31 -0400
Subject: [PATCH] updated root script

---
 Armory/aarootscript.C | 46 ++++++++++++++++++++++++++++++++++++++++---
 Armory/anasenMS.cpp   | 15 ++++++++++----
 Armory/histcomp.C     | 14 ++++++++++---
 3 files changed, 65 insertions(+), 10 deletions(-)

diff --git a/Armory/aarootscript.C b/Armory/aarootscript.C
index 2c82588..4ff7989 100644
--- a/Armory/aarootscript.C
+++ b/Armory/aarootscript.C
@@ -8,9 +8,18 @@
 #include "TBranch.h"
 #include <iostream>
 
-void aarootscript() {
+void aarootscript(int argument = 0) {
+    std::cout << "\n\n\n";
+    std::cout << "=========================================\n";
+    std::cout << "=========   ANASEN Root Script  =========\n";
+    std::cout << "=========================================\n";
+
     TFile *f = new TFile("SimAnasen1.root");
-    TTree *tree1 = (TTree*)f->Get("tree1");
+    TTree *tree1 = (TTree*)f->Get("tree");
+    if (!tree1) {
+        std::cerr << "Error: tree1 not found in the file!" << std::endl;
+        return;
+    }
     TTree *tree2 = (TTree*)f->Get("tree2");
 
     TTreeReader reader("tree");
@@ -24,11 +33,42 @@ void aarootscript() {
     TTreeReader reader2("tree2");
     TTreeReaderValue<double> branchValue2(reader2, "Tb");
     double totalSum2 = 0;
+    int zeroCount = 0;
     while (reader2.Next()) {
         totalSum2 += *branchValue2;
+        if (*branchValue2 == 0) {
+            zeroCount++;
+        }
     }
     std::cout << "Total Sum: " << totalSum2 << std::endl;
 
     std::cout << "Difference: " << totalSum - totalSum2 << std::endl;
 
-}
\ No newline at end of file
+    std::cout << "Zero Count: " << zeroCount << std::endl;
+
+    std::cout << "Making histograms..." << std::endl;
+
+    gErrorIgnoreLevel = 2001;
+    gROOT->ProcessLine(".x histcomp.C");
+
+    std::cout << "=========================================\n";
+
+    if (argument == 1) {
+        if (tree1) {
+            gROOT->ProcessLine("tree->Print();");
+        } else {
+            std::cout << "Tree1 not found!" << std::endl;
+        }
+    }
+
+    if (argument == 2) {
+        if (tree2) {
+            gROOT->ProcessLine("tree2->Print();");
+        } else {
+            std::cout << "Tree2 not found!" << std::endl;
+        }
+    }
+    std::cout << "\n\n\n";
+
+}
+
diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp
index d6a0b7c..e9dafa5 100644
--- a/Armory/anasenMS.cpp
+++ b/Armory/anasenMS.cpp
@@ -44,6 +44,7 @@ int main(int argc, char **argv){
   c1->SetLogy();
   c1->SetLogx();
   c1->Print("eloss_light.png");
+  int ELossTotal = 0;
 
   auto g2 = elossHeavy;
   g2->SetTitle("Energy Loss Table;Kinetic Energy (MeV);dE/dx (MeV/(mg/cm^{2}))");
@@ -299,8 +300,6 @@ int main(int argc, char **argv){
       tree->Fill();
       //Energy loss
       double dl = (hitPos - vertex).Mag() / 10; // path length in cm (positions in mm)
-      //double EkinLight = Pb.E() - Pb.M();// kinetic energy of light particle before loss
-      //double dedxLight = elossLight->Eval(EkinLight); // interpolate dE/dx
       if (numEvent <= 10){
         printf("Event %d: Ekin before loss = %f MeV, distance = %f cm\n", i, Tb, dl);
       }
@@ -321,12 +320,19 @@ int main(int argc, char **argv){
         dx += step;
         counter++;
         Tb = Tb - dE_step; // update kinetic energy for tree storage
+        if (Tb < 0) {
+          Tb = 0; // prevent negative kinetic energy
+          break;
+        }
       }
       // fill tree2 with energy loss adjusted data
-      tree2->Fill();
+      if (Tb != 0) {
+        tree2->Fill();
+      }
       if (numEvent <= 10){
         printf("Event %d: Tb after energy loss = %f MeV, energy loss = %f MeV\n", i, Tb, tb_temp - Tb);
       }
+      ELossTotal += (tb_temp - Tb);
 
     }else{
       // no valid SX3 hit: mark clearly invalid
@@ -344,7 +350,7 @@ int main(int argc, char **argv){
       reTheta1 = TMath::QuietNaN();
       rePhi1   = TMath::QuietNaN();
       z0       = TMath::QuietNaN();
-
+      //Tb = -12354567; // mark kinetic energy as invalid for no hit case
       // fill tree with original data (no energy loss for these events)
       //tree->Fill();
       //tree2->Fill();
@@ -379,6 +385,7 @@ int main(int argc, char **argv){
   saveFile->Close();
 
   printf("=============== done. saved as %s. tree entries: %d, tree2 entries: %d\n", saveFileName.Data(), count, count2);
+  printf("Average energy loss for light particle: %f MeV\n", (double)ELossTotal / count);
 
   if(enableVis){
     printf("Displaying geometry with %zu tracks from simulation\n", visTrackVertex.size());
diff --git a/Armory/histcomp.C b/Armory/histcomp.C
index 66214fc..b4f813c 100644
--- a/Armory/histcomp.C
+++ b/Armory/histcomp.C
@@ -26,7 +26,7 @@ void histcomp() {
         TBranch *br = (TBranch*)branches->At(i);
         TString name = br->GetName();
 
-        printf("Processing branch: %s\n", name.Data());
+        //printf("Processing branch: %s\n", name.Data());
 
         // Create histograms (auto-range using Draw first)
         TString h1name = "h1_" + name;
@@ -34,8 +34,8 @@ void histcomp() {
 
         // Temporary draw to get range
         tree1->Draw(name, "", "goff");
-        double min = tree1->GetMinimum(name);
-        double max = tree1->GetMaximum(name);
+        double min = fmax(tree1->GetMinimum(name), tree2->GetMinimum(name));
+        double max = fmin(tree1->GetMaximum(name), tree2->GetMaximum(name));
 
         //if (min == max) continue; // skip constant branches
 
@@ -95,6 +95,14 @@ void histcomp() {
         // Save plot (overwrite each run)
         TString filename = "plots/" + name + ".png";
         c->SaveAs(filename);
+        //c->SetLogx(0);
+        //c->SetLogy(0);
+        // Optional: save log plots as well
+        c->SetLogx(1);
+        c->SetLogy(1);
+        h1->SetTitle(name + " (log);"+name+";Counts");
+        c->SaveAs("plots/" + name + "_log.png");
+
 
         // Clean up
         delete c;