From d2953531efc5430ae0d11cc6e88f7ac41986ff34 Mon Sep 17 00:00:00 2001 From: James Szalkie Date: Wed, 27 May 2026 13:18:51 -0400 Subject: [PATCH] Cleaned up old Eloss functions in favor of new program --- Armory/anasenMS.cpp | 94 +-------------------------------------------- 1 file changed, 2 insertions(+), 92 deletions(-) diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index 564c4dc..9254824 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -45,8 +45,6 @@ bool IsDeadSX3(int id){ return dead.count(id); } -static std::set> ReactionProductb = { {4,2} }; // 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,38 +54,13 @@ int main(int argc, char **argv){ // number of events can be overridden from command line int numEvent = 1000000; 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* elossAlpha = LoadELoss("../ELoss/E_vs_x_alpha.dat"); // for light particle (alpha) - TGraph* elossProton = LoadELoss("../ELoss/E_vs_x_proton.dat"); // for heavy particle (proton) - TGraph *invgAlpha = new TGraph(elossAlpha->GetN(), elossAlpha->GetY(), elossAlpha->GetX()); - TGraph *invgProton = new TGraph(elossProton->GetN(), elossProton->GetY(), elossProton->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 = elossAlpha; - g->SetTitle("Energy Loss Table (Alpha);cm;Kinetic Energy (MeV)"); - g->Draw("ALP"); - g->SetLineColor(kRed); - //c1->SetLogy(); - //c1->SetLogx(); - c1->Print("eloss_alpha.png"); - - auto c2 = new TCanvas("c2", "Graph Example", 800, 600); - auto g2 = elossProton; - g2->SetTitle("Energy Loss Table (Proton);cm;Kinetic Energy (MeV)"); - g2->Draw("ALP"); - g2->SetLineColor(kBlue); - c2->Print("eloss_proton.png"); - // Reaction setup: projectile + target configuration, energy, and product IDs TransferReaction transfer; transfer.SetA(14, 7, 0); // e.g., 24Mg (Z=12) with 0 excitation 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(ReactionProductb.begin()->first, ReactionProductb.begin()->second); // identify reaction product b e.g., 1H (proton) + transfer.Setb(1, 1); // identify reaction product b e.g., 1H (proton) transfer.SetB(14, 7); // identify reaction product B e.g., 23Na (Z=11) // TODO add alpha source or alternative reaction channel selection @@ -159,8 +132,6 @@ int main(int argc, char **argv){ // outgoing particles in lab frame (light/heavy) double thetab, phib, Tb; double thetaB, phiB, TB; - double dEb; - double dEB; std::array T; tree->Branch("thetab", &thetab, "thetab/D"); // polar angle of light particle in lab frame tree->Branch("phib", &phib, "phib/D"); // azimuthal angle of light particle in lab frame @@ -168,8 +139,6 @@ int main(int argc, char **argv){ tree->Branch("thetaB", &thetaB, "thetaB/D"); tree->Branch("phiB", &phiB, "phiB/D"); tree->Branch("TB", &TB, "TB/D"); // kinetic energy of heavy particle at vertex (before energy loss) - tree->Branch("dEb", &dEb, "dEb/D"); - tree->Branch("dEB", &dEB, "dEB/D"); // placeholder for heavy particle energy loss, currently set equal to light particle loss for simplicity tree->Branch("T", &T, "T/D"); // placeholder for true Q-value, currently set to 0 for simplicity // excitation state identifiers @@ -227,16 +196,12 @@ int main(int argc, char **argv){ double z0; tree->Branch("z0", &z0, "reconstucted_Z/D"); - TTree* tree2 = tree->CloneTree(0); - tree2->SetName("tree2"); - //========timer TBenchmark clock; bool shown ; clock.Reset(); clock.Start("timer"); shown = false; - int ELossTotal = 0; //================================= Calculate event loop for( int i = 0; i < numEvent ; i++){ @@ -346,53 +311,8 @@ int main(int argc, char **argv){ rePhi1 = pw->GetTrackPhi() * TMath::RadToDeg(); z0 = pw->GetZ0(); - dEb = 0; - dEB = 0; tree->Fill(); - //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; - - dEb = tb_temp - Tb; // total energy loss - if (ReactionProductb.count({4, 2})){ // if light particle is alpha, use alpha energy loss table - double x0b = invgAlpha->Eval(Tb); - x0b = x0b + dl; - Tb = elossAlpha->Eval(x0b); - } else if (ReactionProductb.count({1, 1})){ // if light particle is proton, use proton energy loss table - double x0b = invgProton->Eval(Tb); - x0b = x0b + dl; - Tb = elossProton->Eval(x0b); - } else { - // 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 - double dE_dx = 5; // MeV/cm, placeholder value for energy loss per unit length - Tb = Tb - dE_dx * dl; - } - - //if (Tb < 0) { - // Tb = TMath::QuietNaN(); - //} - - dEb = tb_temp - Tb; // total energy loss - - // 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] = 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(); - - 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); - }else{ // no valid SX3 hit: mark clearly invalid sx3Up = -1; @@ -409,13 +329,10 @@ int main(int argc, char **argv){ reTheta1 = TMath::QuietNaN(); rePhi1 = TMath::QuietNaN(); z0 = TMath::QuietNaN(); - dEb = TMath::QuietNaN(); - 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(); } //#################################################################### Timer @@ -439,16 +356,11 @@ int main(int argc, char **argv){ // write results to ROOT file and close //tree->Write(); - //tree2->Write(); tree->Write("", TObject::kOverwrite); - tree2->Write("", TObject::kOverwrite); int count = tree->GetEntries(); - int count2 = tree2->GetEntries(); 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: %f MeV\n", (double)ELossTotal); - printf("Average energy loss across events: %f MeV\n", (double)ELossTotal / count); + printf("=============== done. saved as %s. tree entries: %d\n", saveFileName.Data(), 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()); @@ -509,8 +421,6 @@ int main(int argc, char **argv){ } delete anasen; - delete elossAlpha; - delete elossProton; return 0;