From 8d08b0d3555a58ea26ea84093a43f53b8aec7a30 Mon Sep 17 00:00:00 2001 From: James Szalkie Date: Thu, 14 May 2026 14:52:22 -0400 Subject: [PATCH] Expanded analysis capabilities and script --- Armory/README_anasenMS.md | 16 +++++- Armory/analyze.C | 112 ++++++++++++++++++++++++++++++++++++++ Armory/anasenMS.cpp | 29 +++++++--- 3 files changed, 146 insertions(+), 11 deletions(-) create mode 100644 Armory/analyze.C diff --git a/Armory/README_anasenMS.md b/Armory/README_anasenMS.md index afce23d..9b0c535 100644 --- a/Armory/README_anasenMS.md +++ b/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 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.) diff --git a/Armory/analyze.C b/Armory/analyze.C new file mode 100644 index 0000000..4f82470 --- /dev/null +++ b/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 + +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"; +} \ No newline at end of file diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index 0623a8f..097b016 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -40,8 +40,8 @@ bool IsDeadCathode(int id){ return dead.count(id); } -bool IsDeadSX3(int id){ - static std::set dead = {}; // add dead SX3 IDs here, 0-8 for upstream, 9-17 for downstream, 18-23 for whole detector +bool IsDeadSX3(int id){ + static std::set 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.SetIncidentEnergyAngle(72, 0, 0); // 5.46 MeV beam for alpha source, 0 polar and azimuthal angle + 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, 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 ExAList = {0}; - std::vector ExList = {0, 1, 2}; + std::vector ExAList = {0}; // projectile excitation states in MeV, e.g., 0 for ground state, 1.37 for first excited state of 24Mg, etc. + std::vector 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; -} +} \ No newline at end of file