Expanded analysis capabilities and script

Armory/README_anasenMS.md

@@ -14,6 +14,10 @@
 
 In `Armory` directory:
 
+the directory contains a make file
+
+Run `make AnasenMS` and it will automatically run
+
 ```bash
 g++ -O2 -o anasenMS anasenMS.cpp ClassTransfer.cpp ClassAnasen.cpp ... `root-config --cflags --libs`
 ```
@@ -23,14 +27,21 @@ g++ -O2 -o anasenMS anasenMS.cpp ClassTransfer.cpp ClassAnasen.cpp ... `root-con
 ## Run
 
 ```bash
-./anasenMS [numEvents]
+./anasenMS [numEvents] vis(optional)
 ```
 
 - `numEvents`: optional integer, default `1000000`
-- Outputs: `SimAnasen1.root` containing `tree`
+- Outputs: `SimAnasen1.root` containing `tree1` and `tree2`
+tree1 contains pre-energy loss calculations
+tree2 contains post-energy loss calculations (subject to change)
 
 ## What the code does
 
+- Detector geometry is built within ClassAnasen.h
+
+-To assign dead channels (anode/cathode/SX3), add ID's to IsDead<detector> boolean functions at top of simulation,
+  inside the set
+
 - Initializes reaction: `TransferReaction transfer`
   - `SetA(24,12,0)` target
   - `SetIncidentEnergyAngle(10,0,0)` beam energy and direction
@@ -58,7 +69,6 @@ g++ -O2 -o anasenMS anasenMS.cpp ClassTransfer.cpp ClassAnasen.cpp ... `root-con
 
 ## Notes
 
-- The code now has expanded inline comments explaining each stage.
 - Important methods are from:
   - `ClassTransfer` (`SetA`, `SetIncidentEnergyAngle`, `Seta`, `Setb`, `SetExA`, `SetExB`, `CalReactionConstant`, `Event`)
   - `ClassAnasen` / `SX3` / `PW` (`FindWireID`, `FindSX3Pos`, `GetHitInfo`, `CalTrack`, `CalTrack2`, `GetTrackTheta`, `GetTrackPhi`, `GetZ0`, `GetHitPosWithSigma`, `GetID`, etc.)

Armory/analyze.C

@@ -0,0 +1,112 @@
+#include "TFile.h"
+#include "TTree.h"
+#include "TCanvas.h"
+#include "TH1D.h"
+#include "TLegend.h"
+#include "TString.h"
+#include "TStyle.h"
+#include <iostream>
+
+void analyze(const char* filename = "SimAnasen1.root")
+{
+    gStyle->SetOptStat(1);
+
+    TFile *f = TFile::Open(filename);
+    if (!f || f->IsZombie()) {
+        std::cerr << "ERROR: cannot open file " << filename << std::endl;
+        return;
+    }
+
+    TTree *tree = (TTree*)f->Get("tree");
+    TTree *tree2 = (TTree*)f->Get("tree2");
+
+    if (!tree) {
+        std::cerr << "ERROR: tree not found\n";
+        return;
+    }
+
+    std::cout << "\n===== BASIC TREE INFO =====\n";
+    tree->Print();
+
+    // Event inspection tools
+    std::cout << "\n===== FIRST 10 EVENTS =====\n";
+    tree->Scan("Tb:TB:anodeID[0]:cathodeID[0]:sx3ID", "", "", 10);
+
+    std::cout << "\n===== SINGLE EVENT EXAMPLE (0) =====\n";
+    tree->Show(0);
+
+    // Quick detector gating examples
+    std::cout << "\n===== GATED STATS =====\n";
+    std::cout << "Events with anodeID[0]==5: "
+              << tree->GetEntries("anodeID[0]==5") << std::endl;
+
+    std::cout << "Events with sx3ID>=0: "
+              << tree->GetEntries("sx3ID>=0") << std::endl;
+
+    // Tb vs TB comparison histogram
+
+    TCanvas *c1 = new TCanvas("c1","Tb vs TB",800,600);
+    //set min and max from tree values
+    double min = tree->GetMinimum("Tb");
+    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);
+
+    tree->Draw("Tb>>hTb","","goff");
+    tree->Draw("TB>>hTB","","goff");
+
+    hTb->SetLineColor(kRed);
+    hTB->SetLineColor(kBlue);
+
+    hTb->Draw("HIST");
+    hTB->Draw("HIST SAME");
+
+    TLegend *leg = new TLegend(0.65,0.75,0.88,0.88);
+    leg->AddEntry(hTb,"Tb (light)","l");
+    leg->AddEntry(hTB,"TB (heavy)","l");
+    leg->Draw();
+
+    c1->SaveAs("Tb_TB_compare.png");
+
+    // 4. Detector-gated histogram
+
+    TCanvas *c2 = new TCanvas("c2","Anode gated Tb",800,600);
+    double min2 = tree->GetMinimum("Tb");
+    double max2 = tree->GetMaximum("Tb");
+    min2 = min2 - max2*0.1;
+    max2 = max2 * 1.1;
+    TH1D *hGate = new TH1D("hGate","Tb (anodeID[0]==5);Energy;Counts",200,min2,max2);
+
+    tree->Draw("Tb>>hGate","anodeID[0]==5","goff");
+
+    hGate->SetLineColor(kGreen+2);
+    hGate->Draw("HIST");
+
+    c2->SaveAs("Tb_anode5.png");
+
+    // Tb vs TB correlation (with gate)
+
+    TCanvas *c3 = new TCanvas("c3","Tb vs TB",800,600);
+
+    tree->Draw("TB:Tb>>h2(200,min,max,200,min,max)","sx3ID>=0","COLZ");
+
+    c3->SaveAs("Tb_vs_TB.png");
+
+    // Make gated trees 
+
+    TFile *out = new TFile("gated_output.root","RECREATE");
+
+    TTree *t_anode5 = tree->CopyTree("anodeID[0]==5");
+    t_anode5->Write("tree_anode5");
+
+    TTree *t_sx3valid = tree->CopyTree("sx3ID>=0");
+    t_sx3valid->Write("tree_sx3valid");
+
+    out->Close();
+
+    std::cout << "\n===== DONE =====\n";
+    std::cout << "Saved plots + gated trees in gated_output.root\n";
+}

Armory/anasenMS.cpp

@@ -40,8 +40,8 @@ bool IsDeadCathode(int id){
     return dead.count(id);
 }
 
-bool IsDeadSX3(int id){
+bool IsDeadSX3(int id){ 
-    static std::set<int> dead = {}; // add dead SX3 IDs here, 0-8 for upstream, 9-17 for downstream, 18-23 for whole detector
+    static std::set<int> dead = {}; // add dead SX3 IDs here, 0-23
     return dead.count(id);
 }
 
@@ -83,16 +83,16 @@ int main(int argc, char **argv){
   // Reaction setup: projectile + target configuration, energy, and product IDs
   TransferReaction transfer;
 
-  transfer.SetA(17,13, 0);           // e.g., 24Mg (Z=12) with 0 excitation
+  transfer.SetA(27,13, 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
+  transfer.SetIncidentEnergyAngle(72, 0, 0); // 5.46 MeV beam for alpha source, 0 polar and azimuthal angle, 72 for 27Al 
   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)
 
   // TODO add alpha source or alternative reaction channel selection
 
   // Excited state lists (target and projectile/excited products)
-  std::vector<float> ExAList = {0};
+  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> ExList = {0, 1, 2};
+  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.
 
   // define vertex position uniform distribution ranges (mm)
   double vertexXRange[2] = {   -5,    5}; // mm
@@ -267,6 +267,10 @@ int main(int argc, char **argv){
     TB = PB.E() - PB.M();
     T[0] = Tb;
     T[1] = TB;
+    if (Tb < 1.5) {
+      //skip event if light particle energy after loss is below detection threshold of 1.5 MeV
+      continue;
+    }
 
     phib = Pb.Phi() * TMath::RadToDeg();
     phiB = PB.Phi() * TMath::RadToDeg();
@@ -367,12 +371,20 @@ int main(int argc, char **argv){
       Tb = elossLight->Eval(x0light); // kinetic energy corresponding to range at hit position
       TB = elossHeavy->Eval(x0heavy); // kinetic energy for heavy particle, currently set equal to light particle loss for simplicity
 
+      if (Tb < 1.5) {
+        //skip event if light particle energy after loss is below detection threshold of 1.5 MeV
+        continue;
+      }
+
       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] = TB; 
+      //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();
       
       if (numEvent <= 10){
@@ -400,6 +412,7 @@ int main(int argc, char **argv){
       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)
+      //comment out tree fill for no hit case  
       //tree->Fill();
       //tree2->Fill();
     }
@@ -501,4 +514,4 @@ int main(int argc, char **argv){
 
   return 0;
 
-}
+}