yaml

Armory/EnergyLoss.h

@@ -20,13 +20,18 @@
 
 class EnergyLossTable {
 public:
+  static std::string& TableDir() { static std::string dir = "../ELoss/"; return dir; }
+  static void SetTableDir(const std::string& dir) {
+    TableDir() = dir.empty() || dir.back() == '/' ? dir : dir + "/";
+  }
+
   static EnergyLossTable& Get(const std::string& particle, const std::string& medium) {
     std::string key = medium + "_" + particle;
     static std::map<std::string, EnergyLossTable*> cache;
     auto it = cache.find(key);
     if (it != cache.end()) return *it->second;
 
-    std::string path = "../ELoss/" + medium + "Loss/E_vs_x_" + particle + ".dat";
+    std::string path = TableDir() + medium + "Loss/E_vs_x_" + particle + ".dat";
     EnergyLossTable* table = new EnergyLossTable(path);
     cache[key] = table;
     return *table;

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 Kinematics.h constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassQQQ.h ClassPW.h ClassAnasen.h EnergyLoss.h AutoHist2D.h anasen_anode_cathode_hyperboloids.h
+AnasenMS: anasenMS.cpp SimConfig.h Kinematics.h constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassQQQ.h ClassPW.h ClassAnasen.h EnergyLoss.h AutoHist2D.h anasen_anode_cathode_hyperboloids.h
 	@echo "--------- making ANASEN Monte Carlo"
 	$(CXX) $(CXXFLAGS) anasenMS.cpp -o AnasenMS $(LDFLAGS)

Armory/SimConfig.h

@@ -0,0 +1,142 @@
+#ifndef SIMCONFIG_H
+#define SIMCONFIG_H
+
+// Minimal YAML-subset reader for simulation configs, plus dated history logging.
+// Supported: nested maps by indentation, "key: value", inline lists "[a, b]", "#" comments.
+
+#include <string>
+#include <vector>
+#include <map>
+#include <fstream>
+#include <sstream>
+#include <iostream>
+#include <cstdlib>
+#include <ctime>
+#include <sys/stat.h>
+
+class SimConfig {
+public:
+  std::string path;     // file the config was loaded from
+  std::string rawText;  // full file contents, used for history/ROOT embedding
+  bool loaded = false;
+
+  bool Load(const std::string& file) {
+    std::ifstream in(file);
+    if (!in) return false;
+    path = file;
+    std::stringstream ss;
+    ss << in.rdbuf();
+    rawText = ss.str();
+    Parse(rawText);
+    loaded = true;
+    return true;
+  }
+
+  bool Has(const std::string& key) const { return values_.count(key) > 0; }
+
+  std::string Str(const std::string& key, const std::string& def) const {
+    auto it = values_.find(key);
+    return it == values_.end() ? def : it->second;
+  }
+  double Double(const std::string& key, double def) const {
+    auto it = values_.find(key);
+    return it == values_.end() ? def : std::atof(it->second.c_str());
+  }
+  int Int(const std::string& key, int def) const {
+    auto it = values_.find(key);
+    return it == values_.end() ? def : std::atoi(it->second.c_str());
+  }
+  bool Bool(const std::string& key, bool def) const {
+    auto it = values_.find(key);
+    if (it == values_.end()) return def;
+    return it->second == "true" || it->second == "yes" || it->second == "1";
+  }
+  std::vector<double> DoubleList(const std::string& key, const std::vector<double>& def) const {
+    auto it = values_.find(key);
+    if (it == values_.end()) return def;
+    std::vector<double> out;
+    for (auto& s : SplitList(it->second)) out.push_back(std::atof(s.c_str()));
+    return out;
+  }
+
+  // Saves a timestamped copy under historyDir and appends a line to historyDir/index.log.
+  // Returns the path of the saved copy ("" on failure).
+  std::string LogToHistory(const std::string& historyDir, int numEvents,
+                           const std::string& outputFile) const {
+    mkdir(historyDir.c_str(), 0755);
+    char stamp[32];
+    std::time_t now = std::time(nullptr);
+    std::strftime(stamp, sizeof(stamp), "%Y-%m-%d_%H%M%S", std::localtime(&now));
+    std::string copy = historyDir + "/" + stamp + ".yaml";
+    std::ofstream out(copy);
+    if (!out) return "";
+    out << "# run: " << stamp << "  events: " << numEvents << "  output: " << outputFile
+        << "  source: " << (loaded ? path : "<built-in defaults>") << "\n"
+        << rawText;
+    if (!rawText.empty() && rawText.back() != '\n') out << "\n";
+    std::ofstream idx(historyDir + "/index.log", std::ios::app);
+    idx << stamp << "\t" << numEvents << "\t" << outputFile << "\t" << copy << "\n";
+    return copy;
+  }
+
+private:
+  std::map<std::string, std::string> values_;  // dotted key -> raw scalar/list text
+
+  static std::string Trim(const std::string& s) {
+    size_t a = s.find_first_not_of(" \t\r");
+    if (a == std::string::npos) return "";
+    size_t b = s.find_last_not_of(" \t\r");
+    return s.substr(a, b - a + 1);
+  }
+
+  static std::string StripQuotes(const std::string& s) {
+    if (s.size() >= 2 && ((s.front() == '"' && s.back() == '"') || (s.front() == '\'' && s.back() == '\'')))
+      return s.substr(1, s.size() - 2);
+    return s;
+  }
+
+  static std::vector<std::string> SplitList(const std::string& s) {
+    std::string t = Trim(s);
+    if (!t.empty() && t.front() == '[') t.erase(0, 1);
+    if (!t.empty() && t.back() == ']') t.pop_back();
+    std::vector<std::string> out;
+    std::stringstream ss(t);
+    std::string item;
+    while (std::getline(ss, item, ',')) {
+      item = Trim(item);
+      if (!item.empty()) out.push_back(StripQuotes(item));
+    }
+    return out;
+  }
+
+  void Parse(const std::string& text) {
+    std::vector<std::pair<int, std::string>> stack;  // (indent, key)
+    std::stringstream ss(text);
+    std::string line;
+    while (std::getline(ss, line)) {
+      // drop comments (a '#' outside quotes)
+      bool inQ = false;
+      char q = 0;
+      for (size_t i = 0; i < line.size(); i++) {
+        char c = line[i];
+        if (inQ) { if (c == q) inQ = false; }
+        else if (c == '"' || c == '\'') { inQ = true; q = c; }
+        else if (c == '#') { line.erase(i); break; }
+      }
+      if (Trim(line).empty()) continue;
+      int indent = line.find_first_not_of(' ');
+      size_t colon = line.find(':', indent);
+      if (colon == std::string::npos) continue;
+      std::string key = Trim(line.substr(indent, colon - indent));
+      std::string val = Trim(line.substr(colon + 1));
+      while (!stack.empty() && stack.back().first >= indent) stack.pop_back();
+      std::string full;
+      for (auto& p : stack) full += p.second + ".";
+      full += key;
+      if (val.empty()) stack.push_back({indent, key});
+      else values_[full] = StripQuotes(val);
+    }
+  }
+};
+
+#endif

Armory/anasenMS.cpp

@@ -29,6 +29,8 @@
 #include "TBranch.h"
 #include <iostream>
 #include <fstream>
+#include "TNamed.h"
+#include "SimConfig.h"   // YAML config + dated history
 bool quit=false;
 void handler(int sig){
   quit=true;
@@ -149,30 +151,57 @@ int main(int argc, char **argv){
   printf("===       ANASEN Monte Carlo          ===\n");
   printf("=========================================\n");
 
-  // number of events can be overridden from command line
+  // optional "--config file.yaml"; stripped from argv so positional args keep working
-  int numEvent = 1000000;
+  std::string configPath = "anasen_config.yaml";
-  if( argc >= 2 ) numEvent = atoi(argv[1]);
+  for( int a = 1; a < argc; a++ ){
+    if( strcmp(argv[a], "--config") == 0 && a + 1 < argc ){
+      configPath = argv[a+1];
+      for( int k = a; k + 2 < argc; k++ ) argv[k] = argv[k+2];
+      argc -= 2;
+      break;
+    }
+  }
+  SimConfig cfg;
+  if( cfg.Load(configPath) ) printf("Loaded configuration: %s\n", configPath.c_str());
+  else printf("Config %s not found, using built-in defaults\n", configPath.c_str());
+
+  // number of events: config value, overridden by first command-line argument
+  int numEvent = cfg.Int("simulation.events", 1000000);
+  if( argc >= 2 && atoi(argv[1]) > 0 ) numEvent = atoi(argv[1]);
+  if( cfg.Int("simulation.seed", -1) >= 0 ) gRandom->SetSeed(cfg.Int("simulation.seed", -1));
   TransferReaction transfer;
 
   // Register signal handler
   std::signal(SIGINT, handler);
 
   //To set beam energy loss, use energy loss app, and create table with target isotope, set Initial beam energy as max energy 
-  transfer.SetA(1, 1, 0);            // 22Mg projectile
+  const auto rxA = cfg.DoubleList("reaction.projectile", {1, 1});
-  transfer.Seta(4, 2);                 // 4He target
+  const auto rxa = cfg.DoubleList("reaction.target", {4, 2});
-  transfer.Setb(4, 2);                 // outgoing proton from the primary transfer
+  const auto rxb = cfg.DoubleList("reaction.light_ejectile", {4, 2});
-  transfer.SetB(1, 1);               // 30Si* heavy product
+  const auto rxB = cfg.DoubleList("reaction.heavy_recoil", {1, 1});
-  double beamE = 7; //56.1
+  transfer.SetA((int)rxA[0], (int)rxA[1], 0);
+  transfer.Seta((int)rxa[0], (int)rxa[1]);
+  transfer.Setb((int)rxb[0], (int)rxb[1]);
+  transfer.SetB((int)rxB[0], (int)rxB[1]);
+  double beamE = cfg.Double("beam.energy_MeV", 7);
+  const double beamEMin = cfg.Double("beam.min_energy_MeV", 6.56);
+  const std::string lightMedium = cfg.Str("eloss.light_medium", "He");
+  EnergyLossTable::SetTableDir(cfg.Str("eloss.table_dir", "../ELoss/"));
   const ReactionConfig reactionConfig = transfer.GetRectionConfig();
   const double beamA = reactionConfig.beamA;                   // mass number of 14N beam
+  const double beamZ = reactionConfig.beamZ;                   // atomic number of 14N beam
+  const double targetA = reactionConfig.targetA;                   // mass number of target
+  const double targetZ = reactionConfig.targetZ;                   // atomic number of target
 
   // Excited state lists (projectile and heavy-product excitation states)
-  std::vector<float> ExAList = {0};              // Beam excited energy
+  const auto exAD = cfg.DoubleList("reaction.projectile_excitation_MeV", {0});
-  std::vector<float> ExList = {0};//{0, 2.235, 3.498, 3.774, 4.831, 4.809, 5.279, 5.614, 5.951, 6.550};  // a/p aluminum
+  const auto exBD = cfg.DoubleList("reaction.heavy_excitation_MeV", {0});
-  //std::vector<float> ExList = {0.84, 1.01, 2.21, 2.73, 2.98, 3}; 
+  std::vector<float> ExAList(exAD.begin(), exAD.end());
+  std::vector<float> ExList(exBD.begin(), exBD.end());
 
-  TGraph* elossBeam = LoadELoss("../ELoss/HeLoss/E_vs_x_proton.dat"); // x = path length (cm), y = beam energy (MeV)
+  const std::string beamTable = cfg.Str("eloss.beam_table", "../ELoss/HeLoss/E_vs_x_proton.dat");
-  TGraph* sigmaXBeam = LoadSigmaXVsEnergy("../ELoss/HeLoss/E_vs_x_proton.dat"); // x = beam energy (MeV), y = distance straggle sigma_x (cm)
+  TGraph* elossBeam = LoadELoss(beamTable); // x = path length (cm), y = beam energy (MeV)
+  TGraph* sigmaXBeam = LoadSigmaXVsEnergy(beamTable); // x = beam energy (MeV), y = distance straggle sigma_x (cm)
   //TGraph* sigmaABeam = LoadSigmaAVsEnergy("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = beam energy (MeV), y = distance straggle sigma_a (cm)
   // Build a temporary inverse (energy -> path) to locate the path at beamE.
   TGraph* elossBeamInverseRaw = new TGraph(elossBeam->GetN());
@@ -210,7 +239,7 @@ int main(int argc, char **argv){
   const int kMHeavy = reactionConfig.recoilHeavyA; // mass number of heavy product
 
   // Nuclear masses are in MeV/c^2; cache them rather than rereading mass20.txt per event.
-  const Isotope al27Iso(27, 13);
+  const Isotope al27Iso(beamZ, beamA);
   const Isotope alphaIso(4, 2);
 
   //const int kZBeam = reactionConfig.beamZ; // atomic number of beam
@@ -236,18 +265,22 @@ int main(int argc, char **argv){
 
   //TGraph* elossLight = LoadELoss("../ELoss/HeLoss/E_vs_x_" + b + ".dat");
   // define vertex position uniform distribution ranges (mm)
-  double vertexXRange[2] = {   -1.75, 1.75}; // mm - 5, 5// 2.5
+  const auto vx = cfg.DoubleList("beam.vertex_x_mm", {-1.75, 1.75});
-  double vertexYRange[2] = {   -1.75, 1.75}; // -5, 5
+  const auto vy = cfg.DoubleList("beam.vertex_y_mm", {-1.75, 1.75});
-  double vertexZRange[2] = { -174.3, 174.3}; // -174.3, 174.3 (full length of gas volume, centered at 0)
+  const auto vz = cfg.DoubleList("beam.vertex_z_mm", {-174.3, 174.3});
-  const double beamEntranceZ = -280 - 174.3; //vertexZRange[0]; // mm, assumed beam entrance into the gas
+  double vertexXRange[2] = { vx[0], vx[1] }; // mm
-  
+  double vertexYRange[2] = { vy[0], vy[1] };
-
+  double vertexZRange[2] = { vz[0], vz[1] }; // full length of gas volume
+  const double beamEntranceZ = cfg.Double("beam.entrance_z_mm", -280 - 174.3); // mm, beam entrance into the gas
 
   // detector resolution / uncertainty parameters
-  double sigmaSX3_W = -1; // mm, if < 0 use mid-point (no spread in SX3 horizontal dimension)
+  double sigmaSX3_W = cfg.Double("resolution.sx3_sigma_w_mm", 2); // mm
-  double sigmaSX3_L = -1;  // mm, vertical spread for SX3
+  double sigmaSX3_L = cfg.Double("resolution.sx3_sigma_l_mm", 2); // mm
-  double sigmaPW_A  = 0;  // normalized anode uncertainty term (0-1)
+  double sigmaPW_A  = cfg.Double("resolution.pw_anode_sigma", 0);  // normalized anode uncertainty term (0-1)
-  double sigmaPW_C  = 0;  // normalized cathode uncertainty term (0-1)
+  double sigmaPW_C  = cfg.Double("resolution.pw_cathode_sigma", 0);  // normalized cathode uncertainty term (0-1)
+  const double anodeSigma = cfg.Double("resolution.anode_position_sigma_mm", 2.3);
+  const double qqqSigmaR = cfg.Double("resolution.qqq_sigma_r_mm", 2);
+  const double qqqSigmaPhi = cfg.Double("resolution.qqq_sigma_phi_mm", 2);
 
   // status printout
   printf("------------ Vertex :\n");
@@ -285,9 +318,12 @@ int main(int argc, char **argv){
   QQQ * qqq = anasen->GetQQQ();         // optional QQQ detector part, not used in this simulation but can be enabled for visualization
 
   // output file + trees
-  TString saveFileName = "SimAnasen1.root";
+  TString saveFileName = cfg.Str("simulation.output_file", "SimAnasen1.root").c_str();
+  const std::string histCopy = cfg.LogToHistory(cfg.Str("simulation.history_dir", "config_history"), numEvent, saveFileName.Data());
+  if( !histCopy.empty() ) printf("Configuration logged to %s\n", histCopy.c_str());
   printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data());
   TFile * saveFile = new TFile(saveFileName, "recreate");
+  TNamed("simConfig", cfg.rawText.c_str()).Write(); // config used for this file
   TTree * tree1 = new TTree("tree1", "tree1");
 
   // beam and CM variables saved in tree
@@ -542,7 +578,7 @@ int main(int argc, char **argv){
     transfer.CalReactionConstant();
 
     // vertex position in target volume
-    beamEnergy = gRandom->Uniform(6.56, beamE); // MeV, sample beam energy at vertex from uniform distribution between 0 and initial beam energy
+    beamEnergy = gRandom->Uniform(beamEMin, beamE); // MeV, uniform beam energy at the vertex
     KEA = beamEnergy / beamA;
     beamPath_cm = elossBeamInverse->Eval(beamEnergy); // beamE maps to x=0 after path shift
     vertexZ = beamEntranceZ + beamPath_cm * 10.0; // cm -> mm
@@ -696,7 +732,7 @@ int main(int argc, char **argv){
       double distance_sx3;
       Esx3 = CalculateEnergyLoss(vertexX, vertexY, vertexZ,
                sx3X, sx3Y, sx3Z,
-               b, "He", Tb,
+               b, lightMedium, Tb,
            distance_sx3);
       dl = distance_sx3;
 
@@ -704,26 +740,26 @@ int main(int argc, char **argv){
       //TVector3 hitSigma = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L);
 
       //smudging
-      TVector3 hitSigma = sx3->GetHitPosWithSigma(2, 2); // args are sigmaSX3_W, sigmaSX3_L, currently set to 0 for no smearing
+      TVector3 hitSigma = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L);
 
       double hitX = hitSigma.X();
       double hitY = hitSigma.Y();
       double hitZ = hitSigma.Z();
 
-      double aX_sigma = gRandom->Gaus(aX, 2.3);
+      double aX_sigma = gRandom->Gaus(aX, anodeSigma);
-      double aY_sigma = gRandom->Gaus(aY, 2.3);
+      double aY_sigma = gRandom->Gaus(aY, anodeSigma);
-      double aZ_sigma = gRandom->Gaus(aZ, 2.3);
+      double aZ_sigma = gRandom->Gaus(aZ, anodeSigma);
                                                                                              
       double distance_A;
       Eanode = CalculateEnergyLoss(vertexX, vertexY, vertexZ,
                aX, aY, aZ,
-               b, "He", Tb,
+               b, lightMedium, Tb,
            distance_A);
 
       double distance_C;
       Ecathode = CalculateEnergyLoss(vertexX, vertexY, vertexZ,
                cX, cY, cZ,
-               b, "He", Tb,
+               b, lightMedium, Tb,
            distance_C);
 
       if (Esx3 <= 0 || Eanode <= 0 || Ecathode <= 0) {
@@ -755,13 +791,13 @@ int main(int argc, char **argv){
         vY_track = r_rhoMin_fix.Y();
         vZ_track = r_rhoMin_fix.Z();
 
-        EBeam_track = elossBeam->Eval((vZ_track + 454.3) / 10);
+        EBeam_track = elossBeam->Eval((vZ_track - beamEntranceZ) / 10);
 
         Tb_recon = CalculateOriginalEnergy(hitX, hitY, hitZ,
                                             vX_track, vY_track, vZ_track,
                                             //0, 0, r_rhoMin_fix.Z(),
                                             //0, 0, vertexZ,
-                                            b, "He", Esx3,
+                                            b, lightMedium, Esx3,
                                             distance_sx3);
 
         const TVector3 projectileDirection = TVector3(hitX - aX_sigma, hitY - aY_sigma, hitZ - aZ_sigma);
@@ -772,7 +808,7 @@ int main(int argc, char **argv){
           
           //kinematic reconstruction
           EBeam_kin = apkin_27Al.getEbeam_givenQ(Tb_recon, Ex, theta_recon); //args are (reaction energy, Ex, thetab) still tracking
-          vZ_kin = elossBeamInverse->Eval(EBeam_kin) * 10 - 454.3;
+          vZ_kin = elossBeamInverse->Eval(EBeam_kin) * 10 + beamEntranceZ;
 
           Ex_kin = apkin_27Al.getExc(Tb_recon, theta_recon);
           /*Checklist: anode smudge, sx3 smudge, beam position off axis, beam angle*/
@@ -803,7 +839,7 @@ int main(int argc, char **argv){
             }
 
             vZ_kin = std::isfinite(EBeam_kin) && EBeam_kin >= 0.0
-                         ? elossBeamInverse->Eval(EBeam_kin) * 10 - 454.3
+                         ? elossBeamInverse->Eval(EBeam_kin) * 10 + beamEntranceZ
                          : TMath::QuietNaN();
 
         }
@@ -856,18 +892,18 @@ int main(int argc, char **argv){
       double distance_qqq;
       Eqqq = CalculateEnergyLoss(vertexX, vertexY, vertexZ,
                            qqqX, qqqY, qqqZ,
-                           b, "He",
+                           b, lightMedium,
                Tb, distance_qqq);
       dl = distance_qqq;
       double distance_A;
       Eanode = CalculateEnergyLoss(vertexX, vertexY, vertexZ,
                            aX, aY, aZ,
-                           b, "He",
+                           b, lightMedium,
                Tb, distance_A);
       double distance_C;
       Ecathode = CalculateEnergyLoss(vertexX, vertexY, vertexZ,
                            cX, cY, cZ,
-                           b, "He",
+                           b, lightMedium,
                Tb, distance_C);
 
       if (Eqqq <= 0 || Eanode <= 0 || Ecathode <= 0) {
@@ -889,16 +925,16 @@ int main(int argc, char **argv){
 
       // QQQ sits at fixed z: smear in-plane (radial, azimuthal), z stays fixed
       const double qHitPhi = hitPos.Phi();
-      const double qSigR = 2.0, qSigPhi = 2.0; // mm
+      const double qSigR = qqqSigmaR, qSigPhi = qqqSigmaPhi;
       const double dR   = gRandom->Gaus(0, qSigR);
       const double dTan = gRandom->Gaus(0, qSigPhi);
       const double qHitX = hitPos.X() + dR * std::cos(qHitPhi) - dTan * std::sin(qHitPhi);
       const double qHitY = hitPos.Y() + dR * std::sin(qHitPhi) + dTan * std::cos(qHitPhi);
       const double qHitZ = hitPos.Z();
 
-      const double qAX = gRandom->Gaus(aX, 2.3);
+      const double qAX = gRandom->Gaus(aX, anodeSigma);
-      const double qAY = gRandom->Gaus(aY, 2.3);
+      const double qAY = gRandom->Gaus(aY, anodeSigma);
-      const double qAZ = gRandom->Gaus(aZ, 2.3);
+      const double qAZ = gRandom->Gaus(aZ, anodeSigma);
 
       if (qqqID >= 0 && !std::isnan(aX) && !std::isnan(aY) && !std::isnan(aZ) && aZ != 54321 && (aX != 0 || aY != 0 || aZ != 0)) {
         // line through QQQ hit and anode point; vertex = closest approach to the beam (z) axis
@@ -914,11 +950,11 @@ int main(int argc, char **argv){
           vY_track = rVtx.Y();
           vZ_track = rVtx.Z();
 
-          EBeam_track = elossBeam->Eval((vZ_track + 454.3) / 10);
+          EBeam_track = elossBeam->Eval((vZ_track - beamEntranceZ) / 10);
 
           Tb_recon = CalculateOriginalEnergy(qHitX, qHitY, qHitZ,
                                               vX_track, vY_track, vZ_track,
-                                              b, "He", Eqqq,
+                                              b, lightMedium, Eqqq,
                                               distance_qqq);
 
           // polar angle of the track relative to the beam (z) axis
@@ -928,7 +964,7 @@ int main(int argc, char **argv){
           if (reactionConfig.recoilLightA == 1 && !TMath::IsNaN(EBeam_track)) {
             Kinematics apkin_27Al(26.981538408,4.00260325413,1.00782503224,29.973770136,EBeam_track/26.981538408);
             EBeam_kin = apkin_27Al.getEbeam_givenQ(Tb_recon, Ex, theta_recon);
-            vZ_kin = elossBeamInverse->Eval(EBeam_kin) * 10 - 454.3;
+            vZ_kin = elossBeamInverse->Eval(EBeam_kin) * 10 + beamEntranceZ;
             Ex_kin = apkin_27Al.getExc(Tb_recon, theta_recon);
           } else if (reactionConfig.recoilLightA == 4) {
             double beamMomentum = std::sqrt((EBeam_track + al27Iso.Mass) * (EBeam_track + al27Iso.Mass) - al27Iso.Mass * al27Iso.Mass);
@@ -951,7 +987,7 @@ int main(int argc, char **argv){
             }
 
             vZ_kin = std::isfinite(EBeam_kin) && EBeam_kin >= 0.0
-                         ? elossBeamInverse->Eval(EBeam_kin) * 10 - 454.3
+                         ? elossBeamInverse->Eval(EBeam_kin) * 10 + beamEntranceZ
                          : TMath::QuietNaN();
           }
         }

Armory/anasen_config.yaml

@@ -0,0 +1,40 @@
+# ANASEN simulation configuration. Run: ./AnasenMS [events] [vis] [--config file.yaml]
+# Each run's config is copied to config_history/<date_time>.yaml and embedded in the ROOT file.
+
+simulation:
+  events: 1000000        # overridden by first command-line argument
+  output_file: SimAnasen1.root
+  seed: -1               # >= 0 sets the ROOT random seed for reproducibility
+  history_dir: config_history
+
+reaction:                # [A, Z]
+  projectile: [1, 1]
+  target: [4, 2]
+  light_ejectile: [4, 2]
+  heavy_recoil: [1, 1]
+  projectile_excitation_MeV: [0]
+  heavy_excitation_MeV: [0]
+
+beam:
+  energy_MeV: 7          # maximum beam energy (table is shifted so this is the entrance energy)
+  min_energy_MeV: 6.56   # lowest energy sampled uniformly at the vertex
+  entrance_z_mm: -454.3  # beam entrance into the gas
+  vertex_x_mm: [-1.75, 1.75]   # beam spot half-width = upper value, added in quadrature to straggle
+  vertex_y_mm: [-1.75, 1.75]
+  vertex_z_mm: [-174.3, 174.3]
+
+eloss:
+  table_dir: ../ELoss/                 # light-particle tables: <table_dir><medium>Loss/E_vs_x_<particle>.dat
+  beam_table: ../ELoss/HeLoss/E_vs_x_proton.dat   # path (cm) vs energy and straggle for the beam
+  light_medium: He
+
+resolution:
+  sx3_sigma_w_mm: 2
+  sx3_sigma_l_mm: 2
+  qqq_sigma_r_mm: 2
+  qqq_sigma_phi_mm: 2
+  pw_anode_sigma: 0      # normalized anode term (0-1) used in PW track fit
+  pw_cathode_sigma: 0    # normalized cathode term (0-1) used in PW track fit
+  anode_position_sigma_mm: 2.3   # smearing of reconstructed anode crossing
+
+dead_channel: []   # list of dead channels, e.g., ["sx3_1", "qqq_2"]

Armory/config_history/2026-10-15_125111.yaml

@@ -0,0 +1,41 @@
+# run: 2026-10-15_125111  events: 200  output: SimAnasen1.root  source: anasen_config.yaml
+# ANASEN simulation configuration. Run: ./AnasenMS [events] [vis] [--config file.yaml]
+# Each run's config is copied to config_history/<date_time>.yaml and embedded in the ROOT file.
+
+simulation:
+  events: 1000000        # overridden by first command-line argument
+  output_file: SimAnasen1.root
+  seed: -1               # >= 0 sets the ROOT random seed for reproducibility
+  history_dir: config_history
+
+reaction:                # [A, Z]
+  projectile: [1, 1]
+  target: [4, 2]
+  light_ejectile: [4, 2]
+  heavy_recoil: [1, 1]
+  projectile_excitation_MeV: [0]
+  heavy_excitation_MeV: [0]
+
+beam:
+  energy_MeV: 7          # maximum beam energy (table is shifted so this is the entrance energy)
+  min_energy_MeV: 6.56   # lowest energy sampled uniformly at the vertex
+  entrance_z_mm: -454.3  # beam entrance into the gas
+  vertex_x_mm: [-1.75, 1.75]   # beam spot half-width = upper value, added in quadrature to straggle
+  vertex_y_mm: [-1.75, 1.75]
+  vertex_z_mm: [-174.3, 174.3]
+
+eloss:
+  table_dir: ../ELoss/                 # light-particle tables: <table_dir><medium>Loss/E_vs_x_<particle>.dat
+  beam_table: ../ELoss/HeLoss/E_vs_x_proton.dat   # path (cm) vs energy and straggle for the beam
+  light_medium: He
+
+resolution:
+  sx3_sigma_w_mm: 2
+  sx3_sigma_l_mm: 2
+  qqq_sigma_r_mm: 2
+  qqq_sigma_phi_mm: 2
+  pw_anode_sigma: 0      # normalized anode term (0-1) used in PW track fit
+  pw_cathode_sigma: 0    # normalized cathode term (0-1) used in PW track fit
+  anode_position_sigma_mm: 2.3   # smearing of reconstructed anode crossing
+
+dead_channel: []   # list of dead channels, e.g., ["sx3_1", "qqq_2"]

Armory/config_history/2026-10-15_125209.yaml

@@ -0,0 +1,41 @@
+# run: 2026-10-15_125209  events: 1000000  output: SimAnasen1.root  source: anasen_config.yaml
+# ANASEN simulation configuration. Run: ./AnasenMS [events] [vis] [--config file.yaml]
+# Each run's config is copied to config_history/<date_time>.yaml and embedded in the ROOT file.
+
+simulation:
+  events: 1000000        # overridden by first command-line argument
+  output_file: SimAnasen1.root
+  seed: -1               # >= 0 sets the ROOT random seed for reproducibility
+  history_dir: config_history
+
+reaction:                # [A, Z]
+  projectile: [1, 1]
+  target: [4, 2]
+  light_ejectile: [4, 2]
+  heavy_recoil: [1, 1]
+  projectile_excitation_MeV: [0]
+  heavy_excitation_MeV: [0]
+
+beam:
+  energy_MeV: 7          # maximum beam energy (table is shifted so this is the entrance energy)
+  min_energy_MeV: 6.56   # lowest energy sampled uniformly at the vertex
+  entrance_z_mm: -454.3  # beam entrance into the gas
+  vertex_x_mm: [-1.75, 1.75]   # beam spot half-width = upper value, added in quadrature to straggle
+  vertex_y_mm: [-1.75, 1.75]
+  vertex_z_mm: [-174.3, 174.3]
+
+eloss:
+  table_dir: ../ELoss/                 # light-particle tables: <table_dir><medium>Loss/E_vs_x_<particle>.dat
+  beam_table: ../ELoss/HeLoss/E_vs_x_proton.dat   # path (cm) vs energy and straggle for the beam
+  light_medium: He
+
+resolution:
+  sx3_sigma_w_mm: 2
+  sx3_sigma_l_mm: 2
+  qqq_sigma_r_mm: 2
+  qqq_sigma_phi_mm: 2
+  pw_anode_sigma: 0      # normalized anode term (0-1) used in PW track fit
+  pw_cathode_sigma: 0    # normalized cathode term (0-1) used in PW track fit
+  anode_position_sigma_mm: 2.3   # smearing of reconstructed anode crossing
+
+dead_channel: []   # list of dead channels, e.g., ["sx3_1", "qqq_2"]

Armory/config_history/index.log

@@ -0,0 +1,2 @@
+2026-10-15_125111	200	SimAnasen1.root	config_history/2026-10-15_125111.yaml
+2026-10-15_125209	1000000	SimAnasen1.root	config_history/2026-10-15_125209.yaml