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

rearranged energy loss settings, added macro script, updated mass plotter

Browse Source
This commit is contained in:
James Szalkie 2026-04-13 15:44:39 -04:00
parent 9ae1c99200
commit 2bd5b18f63
3 changed files with 150 additions and 110 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
Armory/aarootscript.C Normal file
View File

@ -0,0 +1,34 @@
#include "TFile.h"
#include "TTree.h"
#include "TGraph.h"
#include "TLegend.h"
#include "TCanvas.h"
#include "TH1D.h"
#include "TObjArray.h"
#include "TBranch.h"
#include <iostream>
void aarootscript() {
TFile *f = new TFile("SimAnasen1.root");
TTree *tree1 = (TTree*)f->Get("tree1");
TTree *tree2 = (TTree*)f->Get("tree2");
TTreeReader reader("tree");
TTreeReaderValue<double> branchValue(reader, "Tb");
double totalSum = 0;
while (reader.Next()) {
totalSum += *branchValue;
}
std::cout << "Total Sum: " << totalSum << std::endl;
TTreeReader reader2("tree2");
TTreeReaderValue<double> branchValue2(reader2, "Tb");
double totalSum2 = 0;
while (reader2.Next()) {
totalSum2 += *branchValue2;
}
std::cout << "Total Sum: " << totalSum2 << std::endl;
std::cout << "Difference: " << totalSum - totalSum2 << std::endl;
}

40
Armory/anasenMS.cpp
View File

@ -270,22 +270,6 @@ int main(int argc, char **argv){
visTrackHitPos.push_back(hitPos);
visTrackWires.push_back({anodeID[0], cathodeID[0]});
}
/*
// apply energy loss from vertex to SX3 hit position (for light particle)
double dl = (hitPos - vertex).Mag() / 10.0; // path length in cm (positions in mm)
double EkinLight = Pb.E() - Pb.M();
double dedxLight = elossLight->Eval(EkinLight); // interpolate dE/dx
double dE_light = dedxLight * dl * density / 1000.0; // adjust for units (example scaling)
if (dE_light < EkinLight) {
Pb.SetE(Pb.E() - dE_light);
// update momentum to conserve direction
double p_new = TMath::Sqrt(Pb.E()*Pb.E() - Pb.M()*Pb.M());
TVector3 dir_new = Pb.Vect().Unit() * p_new;
Pb.SetPxPyPzE(dir_new.X(), dir_new.Y(), dir_new.Z(), Pb.E());
}
// update kinetic energy after loss
Tb = Pb.E() - Pb.M();*/
// reconstruct track from PW readings + SX3 hit
pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false);
reTheta = pw->GetTrackTheta() * TMath::RadToDeg();
@ -298,22 +282,26 @@ int main(int argc, char **argv){
z0 = pw->GetZ0();
tree->Fill();
//printf("Event %d: Tb before energy loss = %.2f MeV\n", i, Tb);
// apply energy loss from vertex to SX3 hit position (for light particle)
double dl = (hitPos - vertex).Mag() / 10.0; // path length in cm (positions in mm)
double EkinLight = Pb.E() - Pb.M();
double dl = (hitPos - vertex).Mag() / 10; // path length in cm (positions in mm)
//printf("Event %d: dl = %f cm\n", i, dl);
double EkinLight = Pb.E() - Pb.M();// kinetic energy of light particle before loss
double dedxLight = elossLight->Eval(EkinLight); // interpolate dE/dx
//printf("Event %d: dE/dx = %f MeV/(mg/cm^2)\n", i, dedxLight);
double dE_light = dedxLight * dl * density / 1000.0; // adjust for units (example scaling)
if (dE_light < EkinLight) {
Pb.SetE(Pb.E() - dE_light);
// update momentum to conserve direction
double p_new = TMath::Sqrt(Pb.E()*Pb.E() - Pb.M()*Pb.M());
TVector3 dir_new = Pb.Vect().Unit() * p_new;
Pb.SetPxPyPzE(dir_new.X(), dir_new.Y(), dir_new.Z(), Pb.E());
}
//printf("Event %d: dE_light = %f MeV\n", i, dE_light);
//if (dE_light < EkinLight) {
Pb.SetE(Pb.E() - dE_light); // update energy after loss
double p_new = TMath::Sqrt(Pb.E()*Pb.E() - Pb.M()*Pb.M()); // update momentum to conserve direction
TVector3 dir_new = Pb.Vect().Unit() * p_new; // maintain original direction, adjust magnitude
Pb.SetPxPyPzE(dir_new.X(), dir_new.Y(), dir_new.Z(), Pb.E()); // update 4-vector with new energy and momentum
//}
// update kinetic energy after loss
Tb = Pb.E() - Pb.M();
Tb = Tb - dE_light;
// fill tree2 with energy loss adjusted data
tree2->Fill();
//printf("Event %d: Tb after energy loss = %.2f MeV\n", i, Tb);
}else{
// no valid SX3 hit: mark clearly invalid

186
Armory/histcomp.C
View File

@ -1,89 +1,107 @@
void histcomp() {
gROOT->SetBatch(kTRUE);
// Open file
TFile *f = new TFile("SimAnasen1.root");
// 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");
// 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
tree1->Draw(name, "", "goff");
double min = tree1->GetMinimum(name);
double max = tree1->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
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, 800, 600);
h1->SetTitle(name + ";"+name+";Counts");
h1->Draw("HIST");
h2->Draw("HIST SAME");
// Legend
TLegend *leg = new TLegend(0.65,0.75,0.88,0.88);
leg->AddEntry(h1, "tree1", "l");
leg->AddEntry(h2, "tree2", "l");
leg->Draw();
// Save plot (overwrite each run)
TString filename = "plots/" + name + ".png";
c->SaveAs(filename);
// Clean up
delete c;
delete h1;
delete h2;
}
printf("Done! Plots saved in ./plots/\n");
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
tree1->Draw(name, "", "goff");
double min = tree1->GetMinimum(name);
double max = tree1->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
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);
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();
// Save plot (overwrite each run)
TString filename = "plots/" + name + ".png";
c->SaveAs(filename);
// Clean up
delete c;
delete h1;
delete h2;
}
printf("Done! Plots saved in ./plots/\n");
}