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

change desnity from pressure

Browse Source
This commit is contained in:
James Szalkie 2026-04-20 17:18:09 -04:00
parent 049365a549
commit 8cf637d0dd
1 changed files with 2 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
Armory/anasenMS.cpp
View File

@ -35,7 +35,7 @@ int main(int argc, char **argv){
// load energy loss tables (assume units: E in MeV, dE/dx in MeV/(mg/cm²), density in mg/cm³) // load energy loss tables (assume units: E in MeV, dE/dx in MeV/(mg/cm²), density in mg/cm³)
TGraph* elossLight = LoadELoss("../ELoss/Eloss_HeAlpha"); // for light particle (alpha) TGraph* elossLight = LoadELoss("../ELoss/Eloss_HeAlpha"); // for light particle (alpha)
TGraph* elossHeavy = LoadELoss("../ELoss/Eloss_p"); // for heavy particle (proton) TGraph* elossHeavy = LoadELoss("../ELoss/Eloss_p"); // for heavy particle (proton)
double density = 0.0000861; // example for aluminum target, adjust as needed double density = (2.1525e-7) * 400; // example for aluminum target, adjust as needed (400 torr is 0.0000861)
auto c1 = new TCanvas("c1", "Graph Example", 800, 600); auto c1 = new TCanvas("c1", "Graph Example", 800, 600);
auto g = elossLight; auto g = elossLight;
g->SetTitle("Energy Loss Table;Kinetic Energy (MeV);dE/dx (MeV/(mg/cm^{2}))"); g->SetTitle("Energy Loss Table;Kinetic Energy (MeV);dE/dx (MeV/(mg/cm^{2}))");
@ -129,11 +129,9 @@ int main(int argc, char **argv){
// outgoing particles in lab frame (light/heavy) // outgoing particles in lab frame (light/heavy)
double thetab, phib, Tb; double thetab, phib, Tb;
double thetaB, phiB, TB; double thetaB, phiB, TB;
double dEb; // energy loss variable for light particle, will be updated after loss calculation
tree->Branch("thetab", &thetab, "thetab/D"); tree->Branch("thetab", &thetab, "thetab/D");
tree->Branch("phib", &phib, "phib/D"); tree->Branch("phib", &phib, "phib/D");
tree->Branch("Tb", &Tb, "Tb/D"); tree->Branch("Tb", &Tb, "Tb/D");
tree->Branch("dEb", &dEb, "dEb/D"); // for energy loss calculation, will be updated with loss
tree->Branch("thetaB", &thetaB, "thetaB/D"); tree->Branch("thetaB", &thetaB, "thetaB/D");
tree->Branch("phiB", &phiB, "phiB/D"); tree->Branch("phiB", &phiB, "phiB/D");
tree->Branch("TB", &TB, "TB/D"); tree->Branch("TB", &TB, "TB/D");
@ -299,7 +297,6 @@ int main(int argc, char **argv){
rePhi1 = pw->GetTrackPhi() * TMath::RadToDeg(); rePhi1 = pw->GetTrackPhi() * TMath::RadToDeg();
z0 = pw->GetZ0(); z0 = pw->GetZ0();
dEb = 0; // initialize energy loss variable for tree storage, will be updated after loss calculation
tree->Fill(); tree->Fill();
//Energy loss //Energy loss
double dl = (hitPos - vertex).Mag() / 10; // path length in cm (positions in mm) double dl = (hitPos - vertex).Mag() / 10; // path length in cm (positions in mm)
@ -307,11 +304,10 @@ int main(int argc, char **argv){
printf("Event %d: Ekin before loss = %f MeV, distance = %f cm\n", i, Tb, dl); printf("Event %d: Ekin before loss = %f MeV, distance = %f cm\n", i, Tb, dl);
} }
double tb_temp = Tb; double tb_temp = Tb;
//double dE_light = dedxLight * dl * density / 1000.0; // adjust for units (example scaling)
double dx = 0; double dx = 0;
double counter = 0; double counter = 0;
while(dx < dl){ while(dx < dl){
double step = 0.01; // cm, step size for energy loss calculation double step = 0.1; // cm, step size for energy loss calculation
if(dx + step > dl) step = dl - dx; // adjust last step to end at hit position if(dx + step > dl) step = dl - dx; // adjust last step to end at hit position
double EkinStep = Tb; // kinetic energy at current step double EkinStep = Tb; // kinetic energy at current step
double dedxStep = elossLight->Eval(EkinStep); // dE/dx at current energy double dedxStep = elossLight->Eval(EkinStep); // dE/dx at current energy
@ -328,7 +324,6 @@ int main(int argc, char **argv){
break; break;
} }
} }
dEb = tb_temp - Tb; // total energy loss for this event
// fill tree2 with energy loss adjusted data // fill tree2 with energy loss adjusted data
if (Tb != 0) { if (Tb != 0) {
tree2->Fill(); tree2->Fill();
@ -354,7 +349,6 @@ int main(int argc, char **argv){
reTheta1 = TMath::QuietNaN(); reTheta1 = TMath::QuietNaN();
rePhi1 = TMath::QuietNaN(); rePhi1 = TMath::QuietNaN();
z0 = TMath::QuietNaN(); z0 = TMath::QuietNaN();
dEb = TMath::QuietNaN(); // initialize energy loss variable for no hit case
//Tb = -12354567; // mark kinetic energy as invalid for no hit case //Tb = -12354567; // mark kinetic energy as invalid for no hit case
// fill tree with original data (no energy loss for these events) // fill tree with original data (no energy loss for these events)
//tree->Fill(); //tree->Fill();