Energy deposition

Armory/README_anasenMS.md

@@ -39,6 +39,7 @@ g++ -O2 -o anasenMS anasenMS.cpp ClassTransfer.cpp ClassAnasen.cpp ... `root-con
 - Vertex and resolution settings:
   - `vertexX/Y/Z` ranges
   - `sigmaSX3_W`, `sigmaSX3_L`, `sigmaPW_A`, `sigmaPW_C`
+- Loads energy loss tables from `../ELoss/` using `TGraph` for interpolation
 - Prepares ROOT output tree and branches for truth/reconstructed
 - Loop over events:
   - Sample excitation and CM direction
@@ -50,6 +51,7 @@ g++ -O2 -o anasenMS anasenMS.cpp ClassTransfer.cpp ClassAnasen.cpp ... `root-con
     - `sx3->FindSX3Pos(...)`
   - Read out wire hits and SX3 channel + depth
   - Apply position smearing for SX3
+  - **Apply energy loss** to light particle using interpolated dE/dx from table, based on path length from vertex to hit
   - Reconstruct track via `pw->CalTrack` and `pw->CalTrack2`
   - Fill ROOT tree
 - At end: write tree, close file, clean up

Armory/anasenMS.cpp

@@ -5,14 +5,19 @@
 #include "TH2.h"          // 2D histograms
 #include "TStyle.h"       // ROOT plotting style controls
 #include "TCanvas.h"      // ROOT canvas drawing
-#include "TBenchmark.h"   // timing measurement
-
+#include "TBenchmark.h"   // timing measurement#include "TGraph.h"       // for energy loss interpolation
 #include "ClassTransfer.h" // Reaction kinematics and MC event generation
 #include "ClassAnasen.h"   // ANASEN detector model classes (SX3, PW, etc.)
 
 //======== Generate light particle based on reaction
 // calculate real and reconstructed tracks and Q-value uncertainty
 
+// Function to load energy loss table from file
+TGraph* LoadELoss(const char* filename) {
+  TGraph* g = new TGraph(filename, "%lg %lg");
+  return g;
+}
+
 int main(int argc, char **argv){
 
   printf("=========================================\n");
@@ -23,6 +28,11 @@ int main(int argc, char **argv){
   int numEvent = 1000000;
   if( argc >= 2 ) numEvent = atoi(argv[1]);
 
+  // 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
+
   // Reaction setup: projectile + target configuration, energy, and product IDs
   TransferReaction transfer;
 
@@ -236,6 +246,21 @@ int main(int argc, char **argv){
       sx3Y = hitPos.Y();
       sx3Z = hitPos.Z();
 
+      // 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();
+
       // reconstruct track from PW readings + SX3 hit
       pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false);
       reTheta = pw->GetTrackTheta() * TMath::RadToDeg();
@@ -295,6 +320,8 @@ int main(int argc, char **argv){
   printf("=============== done. saved as %s. count(hit==1) : %d\n", saveFileName.Data(), count);
 
   delete anasen;
+  delete elossLight;
+  delete elossHeavy;
 
   return 0;