diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index 933fdac..9b322a1 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -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³) TGraph* elossLight = LoadELoss("../ELoss/Eloss_HeAlpha"); // for light particle (alpha) 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 g = elossLight; 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) 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("phib", &phib, "phib/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("phiB", &phiB, "phiB/D"); tree->Branch("TB", &TB, "TB/D"); @@ -299,7 +297,6 @@ int main(int argc, char **argv){ rePhi1 = pw->GetTrackPhi() * TMath::RadToDeg(); z0 = pw->GetZ0(); - dEb = 0; // initialize energy loss variable for tree storage, will be updated after loss calculation tree->Fill(); //Energy loss 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); } double tb_temp = Tb; - //double dE_light = dedxLight * dl * density / 1000.0; // adjust for units (example scaling) double dx = 0; double counter = 0; 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 double EkinStep = Tb; // kinetic energy at current step double dedxStep = elossLight->Eval(EkinStep); // dE/dx at current energy @@ -328,7 +324,6 @@ int main(int argc, char **argv){ break; } } - dEb = tb_temp - Tb; // total energy loss for this event // fill tree2 with energy loss adjusted data if (Tb != 0) { tree2->Fill(); @@ -354,7 +349,6 @@ int main(int argc, char **argv){ reTheta1 = TMath::QuietNaN(); rePhi1 = 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 // fill tree with original data (no energy loss for these events) //tree->Fill();