rtang/ANASEN_analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

plotting updates

Browse Source
This commit is contained in:
James Szalkie 2026-08-24 10:42:51 -04:00
parent 9e9222283e
commit a0c4e3a7ac
13 changed files with 46 additions and 228 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
.DS_Store vendored
View File

Binary file not shown.

BIN
Armory/.DS_Store vendored
View File

Binary file not shown.

15
Armory/AutoHist2D.h
View File

@ -1,6 +1,7 @@
#ifndef AUTOHIST2D_H #ifndef AUTOHIST2D_H
#define AUTOHIST2D_H #define AUTOHIST2D_H
#include <TCanvas.h>
#include <TH2F.h> #include <TH2F.h>
#include <algorithm> #include <algorithm>
#include <iostream> #include <iostream>
@ -8,6 +9,7 @@
#include <string> #include <string>
#include <unordered_map> #include <unordered_map>
#include <vector> #include <vector>
#include <filesystem>
// Buffers (name -> x,y) pairs during the event loop and, on WriteAll(), // Buffers (name -> x,y) pairs during the event loop and, on WriteAll(),
// books a TH2F per name with a fixed number of bins per axis and a range // books a TH2F per name with a fixed number of bins per axis and a range
@ -49,6 +51,7 @@ public:
TH2F h(name.c_str(), title.c_str(), nbins, xlo, xhi, nbins, ylo, yhi); TH2F h(name.c_str(), title.c_str(), nbins, xlo, xhi, nbins, ylo, yhi);
for (size_t i = 0; i < buf.x.size(); ++i) h.Fill(buf.x[i], buf.y[i]); for (size_t i = 0; i < buf.x.size(); ++i) h.Fill(buf.x[i], buf.y[i]);
h.Write(); h.Write();
SavePNG(h, name);
std::cout << "AutoHist2D: wrote \"" << name << "\" (" << buf.x.size() std::cout << "AutoHist2D: wrote \"" << name << "\" (" << buf.x.size()
<< " entries, " << nbins << "x" << nbins << " bins, range [" << " entries, " << nbins << "x" << nbins << " bins, range ["
@ -58,11 +61,23 @@ public:
} }
private: private:
static constexpr const char* kPlotDir = "Plots";
struct Buffer { struct Buffer {
std::vector<double> x, y; std::vector<double> x, y;
std::string xTitle, yTitle; std::string xTitle, yTitle;
}; };
// Draws 'h' on a throwaway canvas and saves it under kPlotDir/<name>.png,
// creating the directory first if it doesn't already exist; overwrites
// any existing file of that name.
static void SavePNG(TH2F& h, const std::string& name) {
std::filesystem::create_directories(kPlotDir);
TCanvas c(("c_" + name).c_str(), "", 800, 600);
h.Draw("COLZ");
c.SaveAs((std::string(kPlotDir) + "/" + name + ".png").c_str());
}
// Registered histograms, keyed by name. Function-local static avoids // Registered histograms, keyed by name. Function-local static avoids
// needing a separate translation unit for a header-only class. // needing a separate translation unit for a header-only class.
static std::unordered_map<std::string, Buffer>& Registry() { static std::unordered_map<std::string, Buffer>& Registry() {

2
Armory/ClassQQQ.h
View File

@ -274,7 +274,7 @@ inline void QQQ::CalQQQPos(unsigned short ID,
id = ID; id = ID;
chBk = chBack; chBk = chBack;
chDn = chDown; chDn = chDown;
chUp = chUp; //chUp = chUp;
} }
#endif #endif

4
Armory/Isotope.h
View File

@ -109,8 +109,8 @@ private:
lineMass200 = 2774; lineMass200 = 2774;
} }
char * heliosPath; //char * heliosPath;
bool isFindOnce; //bool isFindOnce;
}; };

2
Armory/Makefile
View File

@ -69,6 +69,6 @@ EventBuilder: EventBuilder.cpp ClassData.h fsuReader.h Hit.h
@echo "--------- making EventBuilder" @echo "--------- making EventBuilder"
$(CXX) $(CXXFLAGS) EventBuilder.cpp -o EventBuilder $(LDFLAGS) $(CXX) $(CXXFLAGS) EventBuilder.cpp -o EventBuilder $(LDFLAGS)
AnasenMS: anasenMS.cpp constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassPW.h ClassAnasen.h EnergyLoss.h AnasenMS: anasenMS.cpp constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassPW.h ClassAnasen.h EnergyLoss.h AutoHist2D.h
@echo "--------- making ANASEN Monte Carlo" @echo "--------- making ANASEN Monte Carlo"
$(CXX) $(CXXFLAGS) anasenMS.cpp -o AnasenMS $(LDFLAGS) $(CXX) $(CXXFLAGS) anasenMS.cpp -o AnasenMS $(LDFLAGS)

BIN
Armory/Training_Data/.DS_Store → Armory/Plots/.DS_Store vendored
View File

Binary file not shown.

45
Armory/anasenMS.cpp
View File

@ -53,7 +53,7 @@ bool IsDeadCathode(int id){
} }
bool IsDeadSX3(int id){ bool IsDeadSX3(int id){
static std::set<int> dead = {};//{0, 2, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23}; // add dead SX3 IDs here, 0-23 1,7,9,3 static std::set<int> dead = {0, 2, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23}; // add dead SX3 IDs here, 0-23 1,7,9,3
return dead.count(id); return dead.count(id);
} }
@ -73,7 +73,7 @@ bool IsDeadSX3FrontDnChannel(int sx3ID, int chDn){
bool IsDeadSX3BackChannel(int sx3ID, int chBk){ bool IsDeadSX3BackChannel(int sx3ID, int chBk){
static std::set<std::pair<int, int>> dead = { static std::set<std::pair<int, int>> dead = {
//{1, 10} {1, 10}
// {sx3ID, back-channel} // {sx3ID, back-channel}
}; };
return dead.count({sx3ID, chBk}); return dead.count({sx3ID, chBk});
@ -155,25 +155,26 @@ int main(int argc, char **argv){
transfer.SetB(30, 14); // 30Si* heavy product transfer.SetB(30, 14); // 30Si* heavy product
const ReactionConfig reactionConfig = transfer.GetRectionConfig(); const ReactionConfig reactionConfig = transfer.GetRectionConfig();
const double beamA = reactionConfig.beamA; // mass number of 14N beam const double beamA = reactionConfig.beamA; // mass number of 14N beam
const double beamE = 72 / beamA; // beam energy in MeV //const double beamE = 72 / beamA; // beam energy in MeV
// Excited state lists (projectile and heavy-product excitation states) // Excited state lists (projectile and heavy-product excitation states)
std::vector<float> ExAList = {0}; // Beam excited energy std::vector<float> ExAList = {0}; // Beam excited energy
std::vector<float> ExList = {3.4}; // Heavy product excited energy std::vector<float> ExList = {0.0, 2.2, 3.4, 6.0}; // Heavy product excited energy
const int kMBeam = reactionConfig.beamA; // mass number of beam const int kMBeam = reactionConfig.beamA; // mass number of beam
const int kMTarget = reactionConfig.targetA; // mass number of target const int kMTarget = reactionConfig.targetA; // mass number of target
const int kMLight = reactionConfig.recoilLightA; // mass number of light ejectile const int kMLight = reactionConfig.recoilLightA; // mass number of light ejectile
const int kMHeavy = reactionConfig.recoilHeavyA; // mass number of heavy product const int kMHeavy = reactionConfig.recoilHeavyA; // mass number of heavy product
const int kZBeam = reactionConfig.beamZ; // atomic number of beam
const int kZTarget = reactionConfig.targetZ; // atomic number of target //const int kZBeam = reactionConfig.beamZ; // atomic number of beam
const int kZLight = reactionConfig.recoilLightZ; // atomic number of light ejectile //const int kZTarget = reactionConfig.targetZ; // atomic number of target
const int kZHeavy = reactionConfig.recoilHeavyZ; // atomic number of heavy product //const int kZLight = reactionConfig.recoilLightZ; // atomic number of light ejectile
bool enableSequentialDecay = false; // turning to false to disable sequential decay for now, can be set to true to enable //const int kZHeavy = reactionConfig.recoilHeavyZ; // atomic number of heavy product
const int decayDaughterA = 20; //bool enableSequentialDecay = false; // turning to false to disable sequential decay for now, can be set to true to enable
const int decayDaughterZ = 10; //const int decayDaughterA = 20;
const int decayEjectA = 1; //const int decayDaughterZ = 10;
const int decayEjectZ = 1; //const int decayEjectA = 1;
//const int decayEjectZ = 1;
std::string b; std::string b;
if (reactionConfig.recoilLightA == 1) { if (reactionConfig.recoilLightA == 1) {
@ -186,7 +187,7 @@ int main(int argc, char **argv){
b = "alpha"; b = "alpha";
} }
TGraph* elossLight = LoadELoss("../ELoss/HeLoss/E_vs_x_" + b + ".dat"); //TGraph* elossLight = LoadELoss("../ELoss/HeLoss/E_vs_x_" + b + ".dat");
// define vertex position uniform distribution ranges (mm) // define vertex position uniform distribution ranges (mm)
double vertexXRange[2] = { -5, 5}; // mm - 5, 5 double vertexXRange[2] = { -5, 5}; // mm - 5, 5
double vertexYRange[2] = { -5, 5}; // -5, 5 double vertexYRange[2] = { -5, 5}; // -5, 5
@ -521,9 +522,9 @@ int main(int argc, char **argv){
trackDir.X() * trackDir.X() + trackDir.Y() * trackDir.Y(); trackDir.X() * trackDir.X() + trackDir.Y() * trackDir.Y();
if (transverseDirection2 > 0.0) { if (transverseDirection2 > 0.0) {
const double pathToRhoMin = /*const double pathToRhoMin =
-(hitPos.X() * trackDir.X() + hitPos.Y() * trackDir.Y()) -(hitPos.X() * trackDir.X() + hitPos.Y() * trackDir.Y())
/ transverseDirection2; / transverseDirection2;*/
//const TVector3 rhoMin = hitPos + pathToRhoMin * trackDir; //const TVector3 rhoMin = hitPos + pathToRhoMin * trackDir;
const TVector3 rhoMin = vertex; const TVector3 rhoMin = vertex;
@ -566,14 +567,18 @@ int main(int argc, char **argv){
if (Esx3 <= 0 || Eanode <= 0 || Ecathode <= 0) { if (Esx3 <= 0 || Eanode <= 0 || Ecathode <= 0) {
Esx3 = NAN; Esx3 = NAN;
beamEnergy = NAN;
Ex = NAN;
continue;
} }
Edet = Esx3; Edet = Esx3;
Eqqq = TMath::QuietNaN(); // mark QQQ energy as invalid for SX3 hit case
AutoHist2D::Fill("beamEnergy_vs_vZ", vertexZ / 10, beamEnergy, "vZ (cm)", "beamEnergy (MeV)"); AutoHist2D::Fill("beamEnergy_vs_vZ", vertexZ / 10, beamEnergy, "vZ (cm)", "beamEnergy (MeV)");
AutoHist2D::Fill("EPC x sin(theta) vs Esx3", Esx3, EPC * sin(reTheta * TMath::DegToRad()), "Esx3 (MeV)", "EPC x sin(theta) (MeV)"); AutoHist2D::Fill("EPC x sin(theta) vs Esx3", Esx3, EPC * sin(thetab * TMath::DegToRad()), "Esx3 (MeV)", "EPC x sin(theta) (MeV)");
tree1->Fill(); tree1->Fill();
}else if (false){//(qqqID >= 0){ }else if (qqqID >= 0){
// handle QQQ hit case // handle QQQ hit case
sx3Up = -1; sx3Up = -1;
sx3Dn = -1; sx3Dn = -1;
@ -639,10 +644,14 @@ int main(int argc, char **argv){
if (Eqqq <= 0 || Eanode <= 0 || Ecathode <= 0) { if (Eqqq <= 0 || Eanode <= 0 || Ecathode <= 0) {
Eqqq = NAN; Eqqq = NAN;
Ex = NAN;
continue;
} }
Esx3 = TMath::QuietNaN(); // mark SX3 energy as invalid for QQQ hit case
Edet = Eqqq; Edet = Eqqq;
EPC = Eanode - Ecathode; EPC = Eanode - Ecathode;
AutoHist2D::Fill("beamEnergy_vs_vZ", vertexZ / 10, beamEnergy, "vZ (cm)", "beamEnergy (MeV)"); AutoHist2D::Fill("beamEnergy_vs_vZ", vertexZ / 10, beamEnergy, "vZ (cm)", "beamEnergy (MeV)");
AutoHist2D::Fill("EPC x sin(theta) vs Eqqq", Eqqq, EPC * sin(thetab * TMath::DegToRad()), "Eqqq (MeV)", "EPC x sin(theta) (MeV)");
beamEnergy = TMath::QuietNaN(); // mark beam energy as invalid for QQQ hit case beamEnergy = TMath::QuietNaN(); // mark beam energy as invalid for QQQ hit case
tree1->Fill(); tree1->Fill();

178
Armory/histcomp.C
View File

@ -1,178 +0,0 @@
void histcomp() {
gROOT->SetBatch(kTRUE);
// Open file
TFile *f = new TFile("SimAnasen1.root");
// Get trees (MAKE SURE names are correct)
TTree *tree1 = (TTree*)f->Get("tree");
TTree *tree2 = (TTree*)f->Get("tree2");
if (!tree1 || !tree2) {
printf("Error: could not find trees. Check names!\n");
return;
}
// Create output directory (overwrite-safe)
gSystem->Exec("mkdir -p plots");
// Get list of branches
TObjArray *branches = tree1->GetListOfBranches();
int nBranches = branches->GetEntries();
//int nBranches = 1;
// Loop over branches
for (int i = 0; i < nBranches; i++) {
TBranch *br = (TBranch*)branches->At(i);
TString name = br->GetName();
//printf("Processing branch: %s\n", name.Data());
// Create histograms (auto-range using Draw first)
TString h1name = "h1_" + name;
TString h2name = "h2_" + name;
// Temporary draw to get range
double min, max;
if(name == "T"){
//Get minimum value of T[0] and use as min
min = tree2->GetMinimum("Tb");
max = tree1->GetMaximum("TB");
}else{
tree1->Draw(name, "", "goff");
min = fmin(tree1->GetMinimum(name),
tree2->GetMinimum(name));
max = fmax(tree1->GetMaximum(name),
tree2->GetMaximum(name));
}
//if (min == max) continue; // skip constant branches
// Expand range slightly
double margin = 0.1 * (max - min);
min -= margin;
max += margin;
TH1D *h1 = new TH1D(h1name, name, 100, min, max);
TH1D *h2 = new TH1D(h2name, name, 100, min, max);
// Fill histograms
if(name == "T"){
// Fill both array elements into same histogram
tree1->Draw("Tb>>+" + h1name, "", "goff");
tree1->Draw("TB>>+" + h1name, "", "goff");
tree2->Draw("Tb>>+" + h2name, "", "goff");
tree2->Draw("TB>>+" + h2name, "", "goff");
}else{
tree1->Draw(name + ">>" + h1name, "", "goff");
tree2->Draw(name + ">>" + h2name, "", "goff");
}
// Style
h1->SetLineColor(kRed);
h1->SetLineWidth(2);
h2->SetLineColor(kBlue);
h2->SetLineWidth(2);
// Normalize (optional but useful)
//if (h1->GetEntries() > 0) h1->Scale(1.0 / h1->GetEntries());
//if (h2->GetEntries() > 0) h2->Scale(1.0 / h2->GetEntries());
// Canvas
TCanvas *c = new TCanvas("c", name, 900, 600); //arguments are (name, title, width, height)
c->SetRightMargin(0.18);
c->Modified();
c->Update();
h1->SetTitle(name + ";"+name+";Counts");
h1->Draw("HIST");
h2->Draw("HIST SAME");
gPad->Update();
TPaveStats *st = (TPaveStats*)h1->FindObject("stats");
st->SetX1NDC(0.85); // New X start (left)
st->SetY1NDC(0.5); // New Y start (bottom)
st->SetX2NDC(0.98); // New X end (right)
st->SetY2NDC(0.8); // New Y end (top)
st->Draw();
gPad->Modified();
gPad->Update();
// Legend
TLegend *leg = new TLegend(0.65 + .2,0.75 + .1,0.88 + .1,0.88 + .1);
leg->AddEntry(h1, "tree1", "l");
leg->AddEntry(h2, "tree2", "l");
leg->Draw();
//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
// Save plot (overwrite each run)
TString filename = "plots/" + name + ".png";
c->SaveAs(filename);
// Optional: save log plots as well
if (false) { // set to True to also save log plots
c->SetLogy(1);
h1->SetTitle(name + " (log);"+name+";Counts");
c->SaveAs("plots/" + name + "_logy.png");
c->SetLogy(0);
c->SetLogx(1);
h1->SetTitle(name + " (log);"+name+";Counts");
c->SaveAs("plots/" + name + "_logx.png");
// Clean up
delete c;
delete h1;
delete h2;
}
}
// dEb on y, SX3z on x
TH2D *h2d = new TH2D("h2d", "dEb vs SX3z;SX3z (cm);dEb (MeV)", 500, tree2->GetMinimum("sx3Z"), tree2->GetMaximum("sx3Z"), 500, tree2->GetMinimum("dEb"), tree2->GetMaximum("dEb")); //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, tree2->GetMinimum("z0"), tree2->GetMaximum("z0"), 500, tree2->GetMinimum("dEb"), tree2->GetMaximum("dEb"));
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");
TH2D *h2theta = new TH2D("h2theta", "dEb vs reTheta", 500, tree2->GetMinimum("reTheta"), tree2->GetMaximum("reTheta"), 500, tree2->GetMinimum("dEb"), tree2->GetMaximum("dEb"));
tree2->Draw("dEb:reTheta>>h2theta", "", "goff"); // arguments are "y:x>>histogram", "selection", "options"
TCanvas *c2theta = new TCanvas("c2theta", "dEb vs reTheta", 900, 600);
h2theta->Draw("COLZ");
c2theta->SaveAs("plots/dEb_vs_reTheta.png");
TH2D *h2phi = new TH2D("h2phi", "dEb vs rePhi", 500, tree2->GetMinimum("rePhi"), tree2->GetMaximum("rePhi"), 500, tree2->GetMinimum("dEb"), tree2->GetMaximum("dEb"));
tree2->Draw("dEb:rePhi>>h2phi", "", "goff"); // arguments are "y:x>>histogram", "selection", "options"
TCanvas *c2phi = new TCanvas("c2phi", "dEb vs rePhi", 900, 600);
h2phi->Draw("COLZ");
c2phi->SaveAs("plots/dEb_vs_rePhi.png");
TH2D *h2dE = new TH2D("h2dE", "dEb vs Tb", 500, tree2->GetMinimum("Tb"), tree2->GetMaximum("Tb"), 500, tree2->GetMinimum("dEb"), tree2->GetMaximum("dEb"));
tree2->Draw("dEb:Tb>>h2dE", "", "goff"); // arguments are "y:x>>histogram", "selection", "options"
TCanvas *c2dE = new TCanvas("c2dE", "dEb vs Tb", 900, 600);
h2dE->Draw("COLZ");
c2dE->SaveAs("plots/dEb_vs_Tb.png");
printf("Done! Plots saved in ./plots/\n");
}

4
Armory/run_script.C
View File

@ -1,4 +0,0 @@
.L ANASEN_model.C
.L anasenMS_root.cpp+
ANASEN_model();
Run(10);

BIN
ELoss/.DS_Store vendored
View File

Binary file not shown.

BIN
ELoss/AlSi.zip
View File

Binary file not shown.

24
ELoss/EnergyLoss.py
View File

@ -565,30 +565,6 @@ class MyInteractiveApp(cmd.Cmd):
print("Please input particle, final energy from detector, and distance travelled") print("Please input particle, final energy from detector, and distance travelled")
def do_uproot_file(self, arg):
"""Open a specific root file for inspection"""
args = shlex.split(arg)
if len(args) > 0:
filename = args[0]
else:
filename = self.rootFile
try:
print(f"Opening {filename}")
# Try Armory path first
try:
self.file = uproot.open(f"../Armory/{filename}")
except FileNotFoundError:
self.file = uproot.open(filename)
print("File loaded successfully.")
print("Keys:", self.file.keys())
except Exception as e:
print("Error opening file:", e)
if __name__ == "__main__": if __name__ == "__main__":