diff --git a/Armory/ClassAnasen.h b/Armory/ClassAnasen.h
index db892b7..ff7bd1e 100644
--- a/Armory/ClassAnasen.h
+++ b/Armory/ClassAnasen.h
@@ -108,17 +108,21 @@ inline void ANASEN::Construct3DModel(int anodeID1, int anodeID2, int cathodeID1,
   geom = new TGeoManager("Detector", "ANASEN");
 
   //--- define some materials
-  TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
+  TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0); //name, A, Z, density
+  TGeoMaterial *matHe = new TGeoMaterial("He", 4.0026, 2, 0.000861);
   TGeoMaterial *matAl = new TGeoMaterial("Al", 26.98,13,2.7);
+  TGeoMaterial *matSi = new TGeoMaterial("Si", 28.085,14,2.33);
   //--- define some media
-  TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
+  TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum); //name, number of materials, material
+  TGeoMedium *He = new TGeoMedium("He",2, matHe);
   TGeoMedium *Al = new TGeoMedium("Root Material",2, matAl);
+  TGeoMedium *Si = new TGeoMedium("Si",3, matSi);
 
   //--- make the top container volume 
   Double_t worldx = 200.; //mm
   Double_t worldy = 200.; //mm
   Double_t worldz = 200.; //mm
-  worldBox = geom->MakeBox("ROOT", Vacuum, worldx, worldy, worldz);
+  worldBox = geom->MakeBox("ROOT", He, worldx, worldy, worldz); // name, medium, x half-length, y half-length, z half-length
   geom->SetTopVolume(worldBox);
 
   //--- making axis
@@ -183,7 +187,7 @@ inline void ANASEN::Construct3DModel(int anodeID1, int anodeID2, int cathodeID1,
                                                     new TGeoRotation("rot1", wirePhi , wireTheta, 0.)));
   }
 
-  TGeoVolume * sx3Det = geom->MakeBox("box", Al, 0.1, sx3->GetWidth()/2, sx3->GetLength()/2);
+  TGeoVolume * sx3Det = geom->MakeBox("box", Si, 0.1, sx3->GetWidth()/2, sx3->GetLength()/2);
   sx3Det->SetLineColor(kGreen+3);
 
   for( int i = 0; i < sx3->GetNumDet(); i++){
@@ -243,7 +247,7 @@ inline  void ANASEN::DrawTrack(TVector3 pos, TVector3 direction, bool drawEstima
   Track->SetLineColor(kRed);
   worldBox->AddNode(Track, 1, new TGeoCombiTrans( pos.X(), pos.Y(), pos.Z(), new TGeoRotation("rotA", phi +  90, theta, 0.)));
 
-  TGeoVolume * startPos = geom->MakeSphere("startPos", 0, 0, 3);
+  TGeoVolume * startPos = geom->MakeSphere("startPos", 0, 0, 3); 
   startPos->SetLineColor(kBlack);
   worldBox->AddNode(startPos, 3, new TGeoCombiTrans( pos.X(), pos.Y(), pos.Z(), new TGeoRotation("rotA", 0, 0, 0.)));
 
diff --git a/Armory/ClassTransfer.h b/Armory/ClassTransfer.h
index c12c43f..c744431 100644
--- a/Armory/ClassTransfer.h
+++ b/Armory/ClassTransfer.h
@@ -40,7 +40,7 @@ public:
   std::string beamStoppingPowerFile;         ///stopping_power_for_beam
   std::string recoilLightStoppingPowerFile;  ///stopping_power_for_light_recoil
   std::string recoilHeavyStoppingPowerFile;  ///stopping_power_for_heavy_recoil
-  bool isDecay;        ///isDacay
+  bool isDecay;        ///isDecay
   int heavyDecayA;     ///decayNucleus_A
   int heavyDecayZ;     ///decayNucleus_Z
   bool isRedo;         ///isReDo
diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp
index 9e5c91c..d2d1185 100644
--- a/Armory/anasenMS.cpp
+++ b/Armory/anasenMS.cpp
@@ -34,15 +34,15 @@ int main(int argc, char **argv){
   // load energy loss tables (assume units: E in MeV, dE/dx in MeV/(mg/cm²), density in mg/cm³)
   TGraph* elossLight = LoadELoss("../ELoss/Eloss_alpha"); // for light particle (alpha)
   TGraph* elossHeavy = LoadELoss("../ELoss/Eloss_p");     // for heavy particle (proton)
-  double density = 2.7; // example for aluminum target, adjust as needed
+  double density = 0.0000876; // example for aluminum target, adjust as needed
 
   // Reaction setup: projectile + target configuration, energy, and product IDs
   TransferReaction transfer;
 
   transfer.SetA(24,12, 0);           // e.g., 24Mg (Z=12) with 0 excitation
   transfer.SetIncidentEnergyAngle(10, 0, 0); // 10 MeV beam, 0 polar and azimuthal angle
-  transfer.Seta( 4, 2);              // identify reaction product a in internal indexing
-  transfer.Setb( 1, 1);              // identify reaction product b
+  transfer.Seta( 4, 2);              // identify reaction product a in internal indexing e.g., 4He (alpha)
+  transfer.Setb( 1, 1);              // identify reaction product b e.g., 1H (proton)
 
   // TODO add alpha source or alternative reaction channel selection
 
@@ -100,7 +100,7 @@ int main(int argc, char **argv){
   printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data());
   TFile * saveFile = new TFile(saveFileName, "recreate");
   TTree * tree = new TTree("tree", "tree");
-  TTree* tree2 = tree->CloneTree(0); // for 2D histograms or alternative data structure if needed
+   
 
   // beam and CM variables saved in tree
   double KEA;
@@ -175,6 +175,8 @@ int main(int argc, char **argv){
   double z0;
   tree->Branch("z0",     &z0, "reconstucted_Z/D");
 
+  TTree* tree2 = tree->CloneTree(0);
+  tree2->SetName("tree2");
 
   //========timer
   TBenchmark clock;
@@ -268,10 +270,34 @@ int main(int argc, char **argv){
         visTrackHitPos.push_back(hitPos);
         visTrackWires.push_back({anodeID[0], cathodeID[0]});
       }
+      /*
+      // apply energy loss from vertex to SX3 hit position (for light particle)
+      double dl = (hitPos - vertex).Mag() / 10.0; // path length in cm (positions in mm)
+      double EkinLight = Pb.E() - Pb.M();
+      double dedxLight = elossLight->Eval(EkinLight); // interpolate dE/dx
+      double dE_light = dedxLight * dl * density / 1000.0; // adjust for units (example scaling)
+      if (dE_light < EkinLight) {
+        Pb.SetE(Pb.E() - dE_light);
+        // update momentum to conserve direction
+        double p_new = TMath::Sqrt(Pb.E()*Pb.E() - Pb.M()*Pb.M());
+        TVector3 dir_new = Pb.Vect().Unit() * p_new;
+        Pb.SetPxPyPzE(dir_new.X(), dir_new.Y(), dir_new.Z(), Pb.E());
+      }
+      // update kinetic energy after loss
+      Tb = Pb.E() - Pb.M();*/
 
-      // fill tree with original data before energy loss
+      // reconstruct track from PW readings + SX3 hit
+      pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false);
+      reTheta = pw->GetTrackTheta() * TMath::RadToDeg();
+      rePhi   = pw->GetTrackPhi() * TMath::RadToDeg();
+
+      // alternative track algorithm with uncertainty parameters
+      pw->CalTrack2(hitPos, hitInfo, sigmaPW_A, sigmaPW_C, false);
+      reTheta1 = pw->GetTrackTheta() * TMath::RadToDeg();
+      rePhi1   = pw->GetTrackPhi() * TMath::RadToDeg();
+
+      z0 = pw->GetZ0();
       tree->Fill();
-
       // apply energy loss from vertex to SX3 hit position (for light particle)
       double dl = (hitPos - vertex).Mag() / 10.0; // path length in cm (positions in mm)
       double EkinLight = Pb.E() - Pb.M();
@@ -286,19 +312,6 @@ int main(int argc, char **argv){
       }
       // update kinetic energy after loss
       Tb = Pb.E() - Pb.M();
-
-      // reconstruct track from PW readings + SX3 hit
-      pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false);
-      reTheta = pw->GetTrackTheta() * TMath::RadToDeg();
-      rePhi   = pw->GetTrackPhi() * TMath::RadToDeg();
-
-      // alternative track algorithm with uncertainty parameters
-      pw->CalTrack2(hitPos, hitInfo, sigmaPW_A, sigmaPW_C, false);
-      reTheta1 = pw->GetTrackTheta() * TMath::RadToDeg();
-      rePhi1   = pw->GetTrackPhi() * TMath::RadToDeg();
-
-      z0 = pw->GetZ0();
-
       // fill tree2 with energy loss adjusted data
       tree2->Fill();
 
@@ -321,6 +334,7 @@ int main(int argc, char **argv){
 
       // fill tree with original data (no energy loss for these events)
       tree->Fill();
+      tree2->Fill();
     }
 
     //#################################################################### Timer
@@ -343,8 +357,10 @@ int main(int argc, char **argv){
   }
 
   // write results to ROOT file and close
-  tree->Write();
-  tree2->Write();
+  //tree->Write();
+  //tree2->Write();
+  tree->Write("", TObject::kOverwrite);
+  tree2->Write("", TObject::kOverwrite);
   int count = tree->GetEntries();
   int count2 = tree2->GetEntries();
   saveFile->Close();
diff --git a/Armory/histcomp.C b/Armory/histcomp.C
new file mode 100644
index 0000000..bd3393b
--- /dev/null
+++ b/Armory/histcomp.C
@@ -0,0 +1,89 @@
+void histcomp() {
+
+// Open file
+TFile *f = new TFile("SimAnasen1.root");
+
+// Get trees (MAKE SURE names are correct)
+TTree *tree1 = (TTree*)f->Get("tree");
+TTree *tree2 = (TTree*)f->Get("tree2");
+
+if (!tree1 || !tree2) {
+printf("Error: could not find trees. Check names!\n");
+return;
+}
+
+// Create output directory (overwrite-safe)
+gSystem->Exec("mkdir -p plots");
+
+// Get list of branches
+TObjArray *branches = tree1->GetListOfBranches();
+int nBranches = branches->GetEntries();
+//int nBranches = 1;
+
+// Loop over branches
+for (int i = 0; i < nBranches; i++) {
+    TBranch *br = (TBranch*)branches->At(i);
+    TString name = br->GetName();
+
+    printf("Processing branch: %s\n", name.Data());
+
+    // Create histograms (auto-range using Draw first)
+    TString h1name = "h1_" + name;
+    TString h2name = "h2_" + name;
+
+    // Temporary draw to get range
+    tree1->Draw(name, "", "goff");
+    double min = tree1->GetMinimum(name);
+    double max = tree1->GetMaximum(name);
+
+    //if (min == max) continue; // skip constant branches
+
+    // Expand range slightly
+    double margin = 0.1 * (max - min);
+    min -= margin;
+    max += margin;
+
+    TH1D *h1 = new TH1D(h1name, name, 100, min, max);
+    TH1D *h2 = new TH1D(h2name, name, 100, min, max);
+
+    // Fill histograms
+    tree1->Draw(name + ">>" + h1name, "", "goff");
+    tree2->Draw(name + ">>" + h2name, "", "goff");
+
+    // Style
+    h1->SetLineColor(kRed);
+    h1->SetLineWidth(2);
+
+    h2->SetLineColor(kBlue);
+    h2->SetLineWidth(2);
+
+    // Normalize (optional but useful)
+    //if (h1->GetEntries() > 0) h1->Scale(1.0 / h1->GetEntries());
+    //if (h2->GetEntries() > 0) h2->Scale(1.0 / h2->GetEntries());
+
+    // Canvas
+    TCanvas *c = new TCanvas("c", name, 800, 600);
+
+    h1->SetTitle(name + ";"+name+";Counts");
+    h1->Draw("HIST");
+    h2->Draw("HIST SAME");
+
+    // Legend
+    TLegend *leg = new TLegend(0.65,0.75,0.88,0.88);
+    leg->AddEntry(h1, "tree1", "l");
+    leg->AddEntry(h2, "tree2", "l");
+    leg->Draw();
+
+    // Save plot (overwrite each run)
+    TString filename = "plots/" + name + ".png";
+    c->SaveAs(filename);
+
+    // Clean up
+    delete c;
+    delete h1;
+    delete h2;
+
+}
+
+printf("Done! Plots saved in ./plots/\n");
+}