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

Cleaned up old Eloss functions in favor of new program

Browse Source
This commit is contained in:
James Szalkie 2026-05-27 13:18:51 -04:00
parent b4eb81a0ec
commit d2953531ef
1 changed files with 2 additions and 92 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

94
Armory/anasenMS.cpp
View File

@ -45,8 +45,6 @@ bool IsDeadSX3(int id){
return dead.count(id);
}
static std::set<pair<int,int>> 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<double, 2> 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;