new analysis and comments

Armory/ClassPW.h

@@ -509,10 +509,10 @@ inline void PW::CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA, dou
 inline double PW::GetZ0()
 {
 
-  double x = trackPos.X();
+  [[maybe_unused]]double x = trackPos.X();
-  double y = trackPos.Y();
+  [[maybe_unused]]double y = trackPos.Y();
-  double rho = TMath::Sqrt(x * x + y * y);
+  [[maybe_unused]]double rho = TMath::Sqrt(x * x + y * y);
-  double theta = trackVec.Theta();
+  [[maybe_unused]]double theta = trackVec.Theta();
 
   return trackVec.Z();
 }

Armory/ClassTransfer.h

@@ -16,7 +16,6 @@
 #include "Isotope.h"
 
 class ReactionConfig{
-
 public:
 
   ReactionConfig(){}
@@ -47,9 +46,9 @@ public:
   std::vector<float> beamEx;        ///excitation_energy_of_A[MeV]
 
 
-  void SetReaction(int beamA, int beamZ,
+  void SetReaction(int beamA, int beamZ, // projectile
-                   int targetA, int targetZ,
+                   int targetA, int targetZ, // target
-                   int recoilA, int recoilZ, float beamEnergy_AMeV){
+                   int recoilA, int recoilZ, float beamEnergy_AMeV){ // light recoil, e.g. alpha
     this->beamA = beamA;
     this->beamZ = beamZ;
     this->targetA = targetA;
@@ -178,8 +177,8 @@ public:
   void Setb(int A, int Z);
   void SetB(int A, int Z);
   void SetIncidentEnergyAngle(double KEA, double theta, double phi);
-  void SetExA(double Ex);
+  void SetExA(double Ex); // excitation energy of A in MeV
-  void SetExB(double Ex);
+  void SetExB(double Ex); // excitation energy of B in MeV
   void SetReactionFromFile(string settingFile);
   
   TString GetReactionName();
@@ -247,7 +246,7 @@ TransferReaction::TransferReaction(){
    Seta(4,2);
    Setb(1,1);
    SetB(27,13);
-   TA = 2.5;
+   TA = 2.5; // MeV/u
    T = TA * reaction.beamA;
    
    ExA = 0;
@@ -311,7 +310,7 @@ void TransferReaction::SetB(int A, int Z){
   isBSet = true;
 }
 
-void TransferReaction::SetIncidentEnergyAngle(double KEA, double theta, double phi){
+void TransferReaction::SetIncidentEnergyAngle(double KEA, double theta, double phi){ // KEA in MeV/u, theta and phi in degree
   this->TA = KEA;
   this->T = TA * reaction.beamA;
   this->thetaIN = theta;

Armory/analyze.C

@@ -51,7 +51,6 @@ void analyze(const char* filename = "SimAnasen1.root")
     double max = tree->GetMaximum("TB");
     min = min - max*0.1;
     max = max * 1.1;
-    printf("Tb min: %f, TB max: %f\n", min, max);
     TH1D *hTb = new TH1D("hTb","Tb and TB;Energy (MeV);Counts",200,min,max); //arguments are name, title (with axis labels), number of bins, x-min, x-max
     TH1D *hTB = new TH1D("hTB","",200,min,max);
 
@@ -89,9 +88,9 @@ void analyze(const char* filename = "SimAnasen1.root")
 
     // Tb vs TB correlation (with gate)
 
-    TCanvas *c3 = new TCanvas("c3","Tb vs TB",800,600);
+    TCanvas *c3 = new TCanvas("c3","dEb vs SX3z",800,600);
 
-    tree->Draw("TB:Tb>>h2(200,min,max,200,min,max)","sx3ID>=0","COLZ");
+    tree->Draw("TB:Tb>>h2(200,min,max,200,min,max)","","COLZ"); //arguments are "y:x>>histogram(bins,xmin,xmax,bins,ymin,ymax)", "selection", "options"
 
     c3->SaveAs("Tb_vs_TB.png");
 

Armory/anasenMS.cpp

@@ -41,10 +41,12 @@ bool IsDeadCathode(int id){
 }
 
 bool IsDeadSX3(int id){ 
-    static std::set<int> dead = {}; // add dead SX3 IDs here, 0-23
+    static std::set<int> dead = {}; // add dead SX3 IDs here, 0-23 1,7,9,3
     return dead.count(id);
 }
 
+static std::set<pair<int,int>> ReactionProductb = { {1,1} }; // add reaction product b (light particle) A,Z pairs here, e.g. {1,1} for proton, {4,2} for alpha 
+
 int main(int argc, char **argv){
 
   printf("=========================================\n");
@@ -56,43 +58,43 @@ int main(int argc, char **argv){
   if( argc >= 2 ) numEvent = atoi(argv[1]);
 
   // 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_alpha.dat"); // for light particle (alpha)
+  TGraph* elossAlpha = LoadELoss("../ELoss/E_vs_x_alpha.dat"); // for light particle (alpha)
-  TGraph* elossHeavy = LoadELoss("../ELoss/E_vs_x_proton.dat");     // for heavy particle (proton)
+  TGraph* elossProton = LoadELoss("../ELoss/E_vs_x_proton.dat");     // for heavy particle (proton)
-  //TGraph* elossRecoil = LoadELoss("../ELoss/E_vs_x_recoil.dat");     // for recoil particle (if needed)
+  TGraph *invgAlpha = new TGraph(elossAlpha->GetN(), elossAlpha->GetY(), elossAlpha->GetX());
-  TGraph *invgLight = new TGraph(elossLight->GetN(), elossLight->GetY(), elossLight->GetX());
+  TGraph *invgProton = new TGraph(elossProton->GetN(), elossProton->GetY(), elossProton->GetX());
-  TGraph *invgHeavy = new TGraph(elossHeavy->GetN(), elossHeavy->GetY(), elossHeavy->GetX());
+
   
-  /*
   //Plot energy loss tables (sanity check), vis will not work if this is ran without X11 display (e.g. on cluster), so comment out if running in headless mode
   auto c1 = new TCanvas("c1", "Graph Example", 800, 600);
-  auto g = elossLight;
+  auto g = elossAlpha;
-  g->SetTitle("Energy Loss Table (Light);cm;Kinetic Energy (MeV)");
+  g->SetTitle("Energy Loss Table (Alpha);cm;Kinetic Energy (MeV)");
   g->Draw("ALP");
   g->SetLineColor(kRed);
   //c1->SetLogy();
   //c1->SetLogx();
-  c1->Print("eloss_light.png");
+  c1->Print("eloss_alpha.png");
 
   auto c2 = new TCanvas("c2", "Graph Example", 800, 600);
-  auto g2 = elossHeavy;
+  auto g2 = elossProton;
-  g2->SetTitle("Energy Loss Table (Heavy);cm;Kinetic Energy (MeV)");
+  g2->SetTitle("Energy Loss Table (Proton);cm;Kinetic Energy (MeV)");
   g2->Draw("ALP");
   g2->SetLineColor(kBlue);
-  c2->Print("eloss_heavy.png");*/
+  c2->Print("eloss_proton.png");
 
   // Reaction setup: projectile + target configuration, energy, and product IDs
   TransferReaction transfer;
 
-  transfer.SetA(27,13, 0);           // e.g., 24Mg (Z=12) with 0 excitation
+  transfer.SetA(14, 7, 0);           // e.g., 24Mg (Z=12) with 0 excitation
-  transfer.SetIncidentEnergyAngle(72, 0, 0); // 5.46 MeV beam for alpha source, 0 polar and azimuthal angle, 72 for 27Al 
+  transfer.SetIncidentEnergyAngle((42.82/14.0), 0, 0); // arguments are KEA in MeV/u, theta and phi in degree
   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)
+  transfer.Setb(ReactionProductb.begin()->first, ReactionProductb.begin()->second);              // identify reaction product b e.g., 1H (proton)
+  transfer.SetB(17, 8);              // identify reaction product B e.g., 23Na (Z=11)
 
   // TODO add alpha source or alternative reaction channel selection
 
   // Excited state lists (target and projectile/excited products)
-  std::vector<float> ExAList = {0}; // projectile excitation states in MeV, e.g., 0 for ground state, 1.37 for first excited state of 24Mg, etc.
+  std::vector<float> ExAList = {0}; // projectile excitation states in MeV
-  std::vector<float> ExList = {0, 1, 2}; // target excitation states in MeV, e.g., 0 for ground state, 1.37 for first excited state of 24Mg, etc.
+  std::vector<float> ExList = {0}; // target excitation states in MeV
 
   // define vertex position uniform distribution ranges (mm)
   double vertexXRange[2] = {   -5,    5}; // mm
@@ -263,14 +265,14 @@ int main(int argc, char **argv){
     thetab = Pb.Theta() * TMath::RadToDeg();
     thetaB = PB.Theta() * TMath::RadToDeg();
 
-    Tb = Pb.E() - Pb.M();
+    Tb = (Pb.E() - Pb.M()); // kinetic energy of light particle at vertex (before energy loss) units of MeV
-    TB = PB.E() - PB.M();
+    TB = (PB.E() - PB.M());
     T[0] = Tb;
     T[1] = TB;
-    if (Tb < 1.5) {
+    //if (Tb < 1.5) {
-      //skip event if light particle energy after loss is below detection threshold of 1.5 MeV
+    //  //skip event if light particle energy after loss is below detection threshold of 1.5 MeV
-      continue;
+    //  continue;
-    }
+    //}
 
     phib = Pb.Phi() * TMath::RadToDeg();
     phiB = PB.Phi() * TMath::RadToDeg();
@@ -350,39 +352,38 @@ int main(int argc, char **argv){
 
       //Energy loss
       double dl = (hitPos - vertex).Mag(); // path length in units of cm
-
       if (numEvent <= 100){
         //printf("Event %d: Ekin before loss = %f MeV, distance = %f cm\n", i, Tb, dl);
         //printf("Total T before loss: %f MeV\n", T);
       }
-
       double tb_temp = Tb;
-      double tB_temp = TB;
 
       dEb = tb_temp - Tb; // total energy loss
-      dEB = tB_temp - TB; // total energy loss for heavy particle, currently set equal to light particle loss for simplicity
+      if (ReactionProductb.count({4, 2})){ // if light particle is alpha, use alpha energy loss table
-
+        double x0b = invgAlpha->Eval(Tb);
-      double x0light = invgLight->Eval(Tb);
+        x0b = x0b + dl;
-      double x0heavy = invgHeavy->Eval(TB);
+        Tb = elossAlpha->Eval(x0b);
-
+      } else if (ReactionProductb.count({1, 1})){ // if light particle is proton, use proton energy loss table
-      x0light = x0light + dl;
+        double x0b = invgProton->Eval(Tb);
-      x0heavy = x0heavy + dl;
+        x0b = x0b + dl;
-
+        Tb = elossProton->Eval(x0b);
-      Tb = elossLight->Eval(x0light); // kinetic energy corresponding to range at hit position
+      } else {
-      TB = elossHeavy->Eval(x0heavy); // kinetic energy for heavy particle, currently set equal to light particle loss for simplicity
+        // for other particle types, can add additional energy loss tables or use a generic approximation
-
+        // for now, we will just apply a simple linear energy loss as a placeholder
-      if (Tb < 1.5) {
+        double dE_dx = 5; // MeV/cm, placeholder value for energy loss per unit length
-        //skip event if light particle energy after loss is below detection threshold of 1.5 MeV
+        Tb = Tb - dE_dx * dl;
-        continue;
       }
 
+      //if (Tb < 0) {
+      //  Tb = TMath::QuietNaN();
+      //}
+
       dEb = tb_temp - Tb; // total energy loss
-      dEB = tB_temp - TB; // total energy loss for heavy particle, currently set equal to light particle loss for simplicity
       
       // fill tree2 with energy loss adjusted data
       //Fill T so it can make a histogram of both Tb and TB in root script
       T[0] = Tb;
-      T[1] = TB; 
+      T[1] = 0; 
       //to plot both as one histogram in root, can use tree2->Draw("T(0)"); for light particle and tree2->Draw("T(1)") for heavy particle
 
       tree2->Fill();
@@ -390,7 +391,7 @@ int main(int argc, char **argv){
       if (numEvent <= 10){
         //printf("Event %d: Tb after energy loss = %f MeV, energy loss = %f MeV\n", i, Tb, tb_temp - Tb);
       } //to give in scientific notation, use %e instead of %f in the printf format string. For example: printf("Event %d: Tb after energy loss = %e MeV, energy loss = %e MeV\n", i, Tb, tb_temp - Tb);
-      ELossTotal += (tb_temp - Tb) + (tB_temp - TB);
+      ELossTotal += (tb_temp - Tb);
 
     }else{
       // no valid SX3 hit: mark clearly invalid
@@ -446,8 +447,8 @@ int main(int argc, char **argv){
   saveFile->Close();
 
   printf("=============== done. saved as %s. tree entries: %d, tree2 entries: %d\n", saveFileName.Data(), count, count2);
-  printf("Total energy loss across all events: %e MeV\n", (double)ELossTotal);
+  printf("Total energy loss across all events: %f MeV\n", (double)ELossTotal);
-  printf("Average energy loss across events: %e MeV\n", (double)ELossTotal / count);
+  printf("Average energy loss across events: %f MeV\n", (double)ELossTotal / count);
 
   if(enableVis){ // to enable visualization, run with 3rd argument "vis", e.g. "./anasenMC 1000 vis"
     printf("Displaying geometry with %zu tracks from simulation\n", visTrackVertex.size());
@@ -508,8 +509,8 @@ int main(int argc, char **argv){
   }
 
   delete anasen;
-  delete elossLight;
+  delete elossAlpha;
-  delete elossHeavy;
+  delete elossProton;
 
 
   return 0;

Armory/histcomp.C

@@ -143,5 +143,18 @@ void histcomp() {
 
     }
 
+    // dEb on y, SX3z on x
+    TH2D *h2d = new TH2D("h2d", "dEb vs SX3z;SX3z (cm);dEb (MeV)", 500, -110, 110, 500, 0, 12); //arguments are (name, title, xbins, xlow, xup, ybins, ylow, yup)
+    tree2->Draw("dEb:sx3Z>>h2d", "", "goff"); // arguments are "y:x>>histogram", "selection", "options"
+    TCanvas *c2d = new TCanvas("c2d", "dEb vs SX3z", 900, 600);
+    h2d->Draw("COLZ");
+    c2d->SaveAs("plots/dEb_vs_SX3z.png");
+
+    TH2D *h2z = new TH2D("h2z", "dEb vs z0", 500, -1, 1, 500, 0, 12);
+    tree2->Draw("dEb:z0>>h2z", "", "goff"); // arguments are "y:x>>histogram", "selection", "options"
+    TCanvas *c2z = new TCanvas("c2z", "dEb vs z0", 900, 600);
+    h2z->Draw("COLZ");
+    c2z->SaveAs("plots/dEb_vs_z0.png");
+
     printf("Done! Plots saved in ./plots/\n");
 }