new eloss function

Armory/aarootscript.C

@@ -46,10 +46,10 @@ void aarootscript(int argument = 0) {
 
     std::cout << "Zero Count: " << zeroCount << std::endl;
 
-    //std::cout << "Making histograms..." << std::endl;
+    std::cout << "Making histograms..." << std::endl;
 
     gErrorIgnoreLevel = 2001;
-    //gROOT->ProcessLine(".x histcomp.C");
+    gROOT->ProcessLine(".x histcomp.C");
 
     std::cout << "=========================================\n";
 

Armory/anasenMS.cpp

@@ -12,6 +12,8 @@
 #include "ClassTransfer.h" // Reaction kinematics and MC event generation
 #include "ClassAnasen.h"   // ANASEN detector model classes (SX3, PW, etc.)
 #include "TCanvas.h"      // ROOT canvas for drawing
+#include <stdio.h>
+#include <stdlib.h>
 
 //======== Generate light particle based on reaction
 // calculate real and reconstructed tracks and Q-value uncertainty
@@ -31,18 +33,22 @@ int main(int argc, char **argv){
   // number of events can be overridden from command line
   int numEvent = 1000000;
   if( argc >= 2 ) numEvent = atoi(argv[1]);
-
+  //double density = (2.1525e-7) * 1000; // example for aluminum target, adjust as needed (400 torr is 0.0000861)
-  // load energy loss tables (assume units: E in MeV, dE/dx in MeV/(mg/cm²), density in mg/cm³)
+  //char command[256];
-  TGraph* elossLight = LoadELoss("../ELoss/Eloss_HeAlpha"); // for light particle (alpha)
+  //snprintf(command, sizeof(command), "python3 ../ELoss/EvXconverter.py %f", density);
+  //printf("Command: %s\n", command);
+  //system(command); // run the conversion script to generate energy loss tables with correct density
+  // load energy loss tables (assume units: E in MeV, dE/dx in MeV/(mg/cm^2), density in mg/cm^3)
+  TGraph* elossLight = LoadELoss("../ELoss/E_vs_x_light"); // for light particle (alpha)
   TGraph* elossHeavy = LoadELoss("../ELoss/Eloss_p");     // for heavy particle (proton)
-  double density = (2.1525e-7) * 400; // example for aluminum target, adjust as needed (400 torr is 0.0000861)
+
   auto c1 = new TCanvas("c1", "Graph Example", 800, 600);
   auto g = elossLight;
-  g->SetTitle("Energy Loss Table;Kinetic Energy (MeV);dE/dx (MeV/(mg/cm^{2}))");
+  g->SetTitle("Energy Loss Table;cm;Kinetic Energy (MeV)");
   g->Draw("ALP");
   g->SetLineColor(kRed);
-  c1->SetLogy();
+  //c1->SetLogy();
-  c1->SetLogx();
+  //c1->SetLogx();
   c1->Print("eloss_light.png");
   int ELossTotal = 0;
 
@@ -129,12 +135,14 @@ int main(int argc, char **argv){
   // outgoing particles in lab frame (light/heavy)
   double thetab, phib, Tb;
   double thetaB, phiB, TB;
+  double dEb;
   tree->Branch("thetab", &thetab, "thetab/D");
   tree->Branch("phib",     &phib, "phib/D");
   tree->Branch("Tb",         &Tb, "Tb/D");
   tree->Branch("thetaB", &thetaB, "thetaB/D");
   tree->Branch("phiB",     &phiB, "phiB/D");
   tree->Branch("TB",         &TB, "TB/D");
+  tree->Branch("dEb",        &dEb, "dEb/D");
 
   // excitation state identifiers
   int ExAID;
@@ -297,6 +305,7 @@ int main(int argc, char **argv){
       rePhi1   = pw->GetTrackPhi() * TMath::RadToDeg();
 
       z0 = pw->GetZ0();
+      dEb = 0;
       tree->Fill();
       //Energy loss
       double dl = (hitPos - vertex).Mag() / 10; // path length in cm (positions in mm)
@@ -304,9 +313,10 @@ int main(int argc, char **argv){
         printf("Event %d: Ekin before loss = %f MeV, distance = %f cm\n", i, Tb, dl);
       }
       double tb_temp = Tb;
-      double dx = 0;
+      //double dx = 0;
-      double counter = 0;
+      //double counter = 0;
-      while(dx < dl){
+      //energy loss loop
+      /*while(dx < dl){
         double step = 0.1; // cm, step size for energy loss calculation
         if(dx + step > dl) step = dl - dx; // adjust last step to end at hit position
         double EkinStep = Tb; // kinetic energy at current step
@@ -324,6 +334,13 @@ int main(int argc, char **argv){
           break;
         }
       }
+      dEb = tb_temp - Tb; // total energy loss*/
+      TGraph *invg = new TGraph(elossLight->GetN(), elossLight->GetY(), elossLight->GetX());
+      double x0 = invg->Eval(Tb);
+      //double x0 = elossLight->GetX(Tb); // range corresponding to final kinetic energy
+      x0 = x0 + dl;
+      Tb = elossLight->Eval(x0); // kinetic energy corresponding to range at hit position
+      dEb = tb_temp - Tb; // total energy loss
       // fill tree2 with energy loss adjusted data
       if (Tb != 0) {
         tree2->Fill();
@@ -349,6 +366,7 @@ int main(int argc, char **argv){
       reTheta1 = TMath::QuietNaN();
       rePhi1   = TMath::QuietNaN();
       z0       = TMath::QuietNaN();
+      dEb     = 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();