sanitized

Armory/anasenMS.cpp

@@ -139,27 +139,18 @@ int main(int argc, char **argv){
   // Register signal handler
   std::signal(SIGINT, handler);
 
-  TGraph* elossBeam = LoadELoss("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = path length (cm), y = beam energy (MeV)
-  TGraph* elossBeamInverse = new TGraph(elossBeam->GetN());
-  for( int p = 0; p < elossBeam->GetN(); p++ ){
-    double x, y;
-    elossBeam->GetPoint(p, x, y);
-    elossBeamInverse->SetPoint(p, y, x);
-  }
-  elossBeamInverse->Sort(); // TGraph::Eval requires ascending x
-
   //To set beam energy loss, use energy loss app, and create table with target isotope, set Initial beam energy as max energy 
   transfer.SetA(27, 13, 0);            // 22Mg projectile
   transfer.Seta(4, 2);                 // 4He target
   transfer.Setb(1, 1);                 // outgoing proton from the primary transfer
   transfer.SetB(30, 14);               // 30Si* heavy product
+  double beamE = 56.1;
   const ReactionConfig reactionConfig = transfer.GetRectionConfig();
   const double beamA = reactionConfig.beamA;                   // mass number of 14N beam
-  //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, 2.2, 3.4, 6.0};             // Heavy product excited energy
+  std::vector<float> ExList = {0.0};             // Heavy product excited energy
 
   const int kMBeam = reactionConfig.beamA; // mass number of beam
   const int kMTarget = reactionConfig.targetA; // mass number of target
@@ -194,6 +185,36 @@ int main(int argc, char **argv){
   double vertexZRange[2] = { -174.3, 174.3}; // -174.3, 174.3 (full length of gas volume, centered at 0)
   const double beamEntranceZ = -280 - 174.3; //vertexZRange[0]; // mm, assumed beam entrance into the gas
   
+  TGraph* elossBeam = LoadELoss("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = path length (cm), y = beam energy (MeV)
+  // Build a temporary inverse (energy -> path) to locate the path at beamE.
+  TGraph* elossBeamInverseRaw = new TGraph(elossBeam->GetN());
+  for( int p = 0; p < elossBeam->GetN(); p++ ){
+    double x, y;
+    elossBeam->GetPoint(p, x, y);
+    elossBeamInverseRaw->SetPoint(p, y, x);
+  }
+  elossBeamInverseRaw->Sort(); // TGraph::Eval requires ascending x
+
+  const double pathAtBeamE = elossBeamInverseRaw->Eval(beamE);
+
+  // Keep only energies <= beamE and shift path so beamE corresponds to x = 0.
+  TGraph* elossBeamFiltered = new TGraph();
+  for( int p = 0; p < elossBeam->GetN(); p++ ){
+    double x, y;
+    elossBeam->GetPoint(p, x, y);
+    if( y <= beamE ){
+      const int n = elossBeamFiltered->GetN();
+      elossBeamFiltered->SetPoint(n, x - pathAtBeamE, y);
+    }
+  }
+
+  TGraph* elossBeamInverse = new TGraph(elossBeamFiltered->GetN());
+  for( int p = 0; p < elossBeamFiltered->GetN(); p++ ){
+    double x, y;
+    elossBeamFiltered->GetPoint(p, x, y);
+    elossBeamInverse->SetPoint(p, y, x);
+  }
+  elossBeamInverse->Sort(); // TGraph::Eval requires ascending x
 
   // detector resolution / uncertainty parameters
   double sigmaSX3_W = 0; // mm, if < 0 use mid-point (no spread in SX3 horizontal dimension)
@@ -277,8 +298,8 @@ int main(int argc, char **argv){
   tree1->Branch("phiCM",     &phiCM, "phiCM/D");
 
   // outgoing particles in lab frame (light/heavy)
-  double thetab, phib, Tb, qqqTb;
+  double thetab, phib, Tb, qqqTb, sx3Tb;
-  double thetaB, phiB, TB, qqqTB;
+  double thetaB, phiB, TB, qqqTB, sx3TB;
   std::array<double, 2> T;
   tree1->Branch("thetab", &thetab, "thetab/D"); // polar angle of light particle in lab frame
   tree1->Branch("phib",     &phib, "phib/D"); // azimuthal angle of light particle in lab frame
@@ -288,7 +309,9 @@ int main(int argc, char **argv){
   tree1->Branch("TB",         &TB, "TB/D"); // kinetic energy of heavy particle at vertex (before energy loss)
   tree1->Branch("T",          &T, "T/D");   // placeholder for true Q-value, currently set to 0 for simplicity
   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
+  tree1->Branch("qqqTB", &qqqTB, "qqqTB/D"); // kinetic energy of heavy particle at vertex (before energy loss) for events where the light particle hits the QQQ, currently set to 0 for simplicity
+  tree1->Branch("sx3Tb", &sx3Tb, "sx3Tb/D"); // kinetic energy of light particle at vertex (before energy loss) for events where the light particle hits the SX3, currently set to 0 for simplicity
+  tree1->Branch("sx3TB", &sx3TB, "sx3TB/D"); // kinetic energy of heavy particle at vertex (before energy loss) for events where the light particle hits the SX3, currently set to 0 for simplicity
 
   double Esx3, Eqqq, Edet;
   tree1->Branch("Esx3", &Esx3, "Esx3/D");
@@ -405,9 +428,9 @@ int main(int argc, char **argv){
     // vertex position in target volume
     vertexX = (vertexXRange[1]- vertexXRange[0])*gRandom->Rndm() +  vertexXRange[0];
     vertexY = (vertexYRange[1]- vertexYRange[0])*gRandom->Rndm() +  vertexYRange[0];
-    beamEnergy = gRandom->Uniform(0, 56.1);
+    beamEnergy = gRandom->Uniform(0, beamE); // MeV, sample beam energy at vertex from uniform distribution between 0 and initial beam energy
     KEA = beamEnergy / beamA;
-    beamPath_cm = elossBeamInverse->Eval(beamEnergy); // interpolate path length (cm) that gives this beam energy
+    beamPath_cm = elossBeamInverse->Eval(beamEnergy); // beamE maps to x=0 after path shift
     vertexZ = beamEntranceZ + beamPath_cm * 10.0; // cm -> mm
 
     //vertexZ = (vertexZRange[1]- vertexZRange[0])*gRandom->Rndm() +  vertexZRange[0];
@@ -563,7 +586,7 @@ int main(int argc, char **argv){
                b, "He", Tb,
            distance_C);
 
-      EPC = Eanode - Ecathode;
+
 
       if (Esx3 <= 0 || Eanode <= 0 || Ecathode <= 0) {
         Esx3 = NAN;
@@ -571,9 +594,13 @@ int main(int argc, char **argv){
         Ex = NAN;
         continue;
       }
-      Edet = Esx3;
+      sx3Tb = Tb; // for simplicity, using the same kinetic energy for SX3 hit events, can be modified to simulate energy loss if desired
+      sx3TB = TB;
+      qqqTb = TMath::QuietNaN(); // mark kinetic energy as invalid for SX3 hit case
+      qqqTB = TMath::QuietNaN();
       Eqqq = TMath::QuietNaN(); // mark QQQ energy as invalid for SX3 hit case
-
+      Edet = Esx3;
+      EPC = Eanode - Ecathode;
       AutoHist2D::Fill("beamEnergy_vs_vZ", vertexZ / 10, beamEnergy, "vZ (cm)", "beamEnergy (MeV)");
       AutoHist2D::Fill("EPC x sin(theta) vs Esx3", Esx3, EPC * sin(thetab * TMath::DegToRad()), "Esx3 (MeV)", "EPC x sin(theta) (MeV)");
       tree1->Fill();
@@ -595,12 +622,6 @@ int main(int argc, char **argv){
       rePhi1   = TMath::QuietNaN();
       z0       = TMath::QuietNaN();
 
-      qqqTb = Tb; // for simplicity, using the same kinetic energy for QQQ hit events, can be modified to simulate energy loss if desired
-      qqqTB = TB;
-
-      Tb = TMath::QuietNaN(); // mark kinetic energy as invalid for SX3 hit case
-      TB = TMath::QuietNaN();
-
       TVector3 hitPos = qqq->GetHitPos();
 
       qqqX = hitPos.X();
@@ -629,24 +650,29 @@ int main(int argc, char **argv){
       Eqqq = CalculateEnergyLoss(vertexX, vertexY, vertexZ,
                            qqqX, qqqY, qqqZ,
                            b, "He",
-               qqqTb, distance_qqq);
+               Tb, distance_qqq);
       dl = distance_qqq;
       double distance_A;
       Eanode = CalculateEnergyLoss(vertexX, vertexY, vertexZ,
                            aX, aY, aZ,
                            b, "He",
-               qqqTb, distance_A);
+               Tb, distance_A);
       double distance_C;
       Ecathode = CalculateEnergyLoss(vertexX, vertexY, vertexZ,
                            cX, cY, cZ,
                            b, "He",
-               qqqTb, distance_C);
+               Tb, distance_C);
 
       if (Eqqq <= 0 || Eanode <= 0 || Ecathode <= 0) {
           Eqqq = NAN;
           Ex = NAN;
           continue;
       }
+
+      qqqTb = Tb; // for simplicity, using the same kinetic energy for QQQ hit events, can be modified to simulate energy loss if desired
+      qqqTB = TB;
+      sx3Tb = TMath::QuietNaN(); // mark kinetic energy as invalid for QQQ hit case
+      sx3TB = TMath::QuietNaN();
       Esx3 = TMath::QuietNaN(); // mark SX3 energy as invalid for QQQ hit case
       Edet = Eqqq;
       EPC = Eanode - Ecathode;
@@ -656,7 +682,6 @@ int main(int argc, char **argv){
       beamEnergy = TMath::QuietNaN(); // mark beam energy as invalid for QQQ hit case
       tree1->Fill();
 
-    
     }else{
       // no valid SX3 hit: mark clearly invalid
       sx3Up = -1;
@@ -675,6 +700,10 @@ int main(int argc, char **argv){
       z0       = TMath::QuietNaN();
       Tb       = TMath::QuietNaN(); // mark kinetic energy as invalid for no hit case
       TB       = TMath::QuietNaN();
+      sx3Tb     = TMath::QuietNaN();
+      sx3TB     = TMath::QuietNaN();
+      qqqTb     = TMath::QuietNaN();
+      qqqTB     = TMath::QuietNaN();
       // fill tree with original data (no energy loss for these events)
       //comment out tree fill for no hit case  
       //tree1->Fill();

ELoss/material_definitions.py

@@ -0,0 +1,73 @@
+"""Shared material definitions for ELoss analysis scripts."""
+
+import pycatima as catima
+
+R_GAS_CONSTANT = 8.3144  # J/mol/K
+
+
+def _gas_molar_density(p_torr, temp_k):
+    """Return molar density in mol/m^3 from pressure (Torr) and temperature (K)."""
+    p_pa = p_torr * 133.322
+    return p_pa / (R_GAS_CONSTANT * temp_k)
+
+
+def build_material(medium, p_torr, temp_k):
+    """
+    Build a CATIMA material and return (material, density_g_cm3).
+
+    Supported media names:
+      - He
+      - PureHe
+      - Si
+      - kapton
+    """
+    medium_key = medium.strip().lower()
+
+    if medium_key == "he":
+        m_he = 4.0026
+        m_c = 12.0000
+        m_o = 15.9949
+
+        material_def = [
+            (m_he, 2, 0.96),
+            (m_c, 6, 0.04),
+            (m_o, 8, 0.08),
+        ]
+
+        molar_density = _gas_molar_density(p_torr, temp_k)
+        m_mix_avg = 0.96 * m_he + 0.04 * (m_c + 2 * m_o)
+        rho_g_cm3 = (molar_density * m_mix_avg) / 1e6
+
+    elif medium_key == "purehe":
+        m_he = 4.0026
+        material_def = [(m_he, 2, 1.0)]
+
+        molar_density = _gas_molar_density(p_torr, temp_k)
+        m_mix_avg = m_he
+        rho_g_cm3 = (molar_density * m_mix_avg) / 1e6
+
+    elif medium_key == "si":
+        m_si = 28.084
+        material_def = [(m_si, 14, 1.0)]
+        rho_g_cm3 = 2.33
+
+    elif medium_key == "kapton":
+        m_h = 1.008
+        m_c = 12.011
+        m_n = 14.007
+        m_o = 15.999
+
+        material_def = [
+            (m_h, 1, 0.026),
+            (m_c, 6, 0.691),
+            (m_n, 7, 0.073),
+            (m_o, 8, 0.209),
+        ]
+        rho_g_cm3 = 1.42
+
+    else:
+        raise ValueError(f"Unsupported medium: {medium}")
+
+    material = catima.Material(material_def)
+    material.density(rho_g_cm3)
+    return material, rho_g_cm3