From c6c8be8cdee52eef7e5965668f845b6cbabe5c6b Mon Sep 17 00:00:00 2001 From: James Szalkie Date: Wed, 9 Sep 2026 13:24:21 -0400 Subject: [PATCH] anode location blur --- Armory/anasenMS.cpp | 33 ++++++++++++++++++++------------- Armory/testMacro.C | 2 ++ 2 files changed, 22 insertions(+), 13 deletions(-) diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index 8f5358d..dca7cc7 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -168,7 +168,7 @@ int main(int argc, char **argv){ // Excited state lists (projectile and heavy-product excitation states) std::vector ExAList = {0}; // Beam excited energy - std::vector ExList = {0, 2.2, 3.4, 6.0}; // Heavy product excited energy + std::vector ExList = {0, 2.235, 3.498, 4.809, 5.614, 6.550}; // Heavy product excited energy const int kMBeam = reactionConfig.beamA; // mass number of beam const int kMTarget = reactionConfig.targetA; // mass number of target @@ -205,7 +205,7 @@ int main(int argc, char **argv){ TGraph* elossBeam = LoadELoss("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = path length (cm), y = beam energy (MeV) TGraph* sigmaXBeam = LoadSigmaXVsEnergy("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = beam energy (MeV), y = distance straggle sigma_x (cm) - TGraph* sigmaABeam = LoadSigmaAVsEnergy("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = beam energy (MeV), y = distance straggle sigma_a (cm) + //TGraph* sigmaABeam = LoadSigmaAVsEnergy("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = beam energy (MeV), y = distance straggle sigma_a (cm) // Build a temporary inverse (energy -> path) to locate the path at beamE. TGraph* elossBeamInverseRaw = new TGraph(elossBeam->GetN()); for( int p = 0; p < elossBeam->GetN(); p++ ){ @@ -430,10 +430,11 @@ int main(int argc, char **argv){ tree1->Branch("reTheta1", &reTheta1, "reconstucted_theta1/D"); tree1->Branch("rePhi1", &rePhi1, "reconstucted_phi1/D"); - double hitTheta, originalEnergy, sigma_a; + double hitTheta, originalEnergy, sigma_a, theta_recon; tree1->Branch("hitTheta", &hitTheta, "hitTheta/D"); tree1->Branch("originalEnergy", &originalEnergy, "originalEnergy/D"); tree1->Branch("sigma_a", &sigma_a, "sigma_a/D"); + tree1->Branch("theta_recon", &theta_recon, "theta_recon/D"); // reconstructed vertex Z from PW fit double z0; @@ -510,8 +511,8 @@ int main(int argc, char **argv){ double vertexRangeX = std::sqrt((vertexXRange[1] * vertexXRange[1]) + (sigmaX_mm * sigmaX_mm)); double vertexRangeY = std::sqrt((vertexYRange[1] * vertexYRange[1]) + (sigmaX_mm * sigmaX_mm)); //std::cout << "vertexRangeX: " << vertexRangeX << ", vertexRangeY: " << vertexRangeY << std::endl; - vertexX = gRandom->Gaus(0, vertexRangeX) + 5; // mean and standard deviation - vertexY = gRandom->Gaus(0, vertexRangeY) - 5; // mean and standard deviation + vertexX = gRandom->Gaus(0, vertexRangeX); // mean and standard deviation + vertexY = gRandom->Gaus(0, vertexRangeY); // mean and standard deviation //vertexZ = (vertexZRange[1]- vertexZRange[0])*gRandom->Rndm() + vertexZRange[0]; @@ -527,8 +528,8 @@ int main(int argc, char **argv){ //double beamEnergyLoss = elossBeam->Eval(0.0) - beamEnergy; //KEA = beamEnergy / beamA; //KEA = gRandom->Uniform(0, beamE); - sigma_a = std::max(0.0, sigmaABeam->Eval(beamEnergy)); - transfer.SetIncidentEnergyAngle(KEA, sigma_a, 0); //arguments are (kinetic energy, polar angle, azimuthal angle) of the incident particle in the lab frame + //sigma_a = std::max(0.0, sigmaABeam->Eval(beamEnergy)); + transfer.SetIncidentEnergyAngle(KEA, 0, 0); //arguments are (kinetic energy, polar angle, azimuthal angle) of the incident particle in the lab frame transfer.CalReactionConstant(); // isotropic CM direction @@ -658,13 +659,14 @@ int main(int argc, char **argv){ dl = distance_sx3; //std::cout << std::lround(sx3X / 10) * 10 << " " << std::lround(sx3Y / 10) * 10 << " " << std::lround(sx3Z / 10) * 10 << std::endl; - TVector3 hitSigma = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L); + //TVector3 hitSigma = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L); + TVector3 hitSigma = sx3->GetHitPosWithSigma(0, 0); - double hitX = hitSigma.X(); - double hitY = hitSigma.Y(); - double hitZ = hitSigma.Z(); + double hitX = std::lround(hitSigma.X()); + double hitY = std::lround(hitSigma.Y()); + double hitZ = std::lround(hitSigma.Z()); - originalEnergy = CalculateOriginalEnergy(std::lround(hitX), std::lround(hitY), std::lround(hitZ), + originalEnergy = CalculateOriginalEnergy(hitX, hitY, hitZ, 0, 0, std::lround(vertexZ / 10) * 10, b, "He", Esx3, distance_sx3); @@ -696,7 +698,12 @@ int main(int argc, char **argv){ //tree1->Fill(); Kinematics apkin_27Al(26.981538408,4.00260325413,1.00782503224,29.973770136,beamEnergy/26.981538408); //m3 is proton - Ex_recon = apkin_27Al.getExc(originalEnergy, std::lround(thetab)); + //reconstruct theta using anode and sx3 positions, apply a sigma to the anode locations + double aX_sigma = gRandom->Gaus(aX, 3); + double aY_sigma = gRandom->Gaus(aY, 3); + double aZ_sigma = gRandom->Gaus(aZ, 3); + theta_recon = std::asin(TVector3(sx3X - aX_sigma, sx3Y - aY_sigma, 0).Mag() / TVector3(sx3X - aX_sigma, sx3Y - aY_sigma, sx3Z - aZ_sigma).Mag()) * 180.0 / TMath::Pi(); + Ex_recon = apkin_27Al.getExc(originalEnergy, theta_recon); //EBeam_Kin_gs = apkin_27Al.getEbeam_givenQ(Esx3, 0.0, thetab); //EBeam_Kin_2_2 = apkin_27Al.getEbeam_givenQ(Esx3, 2.2, thetab); //EBeam_Kin_3_4 = apkin_27Al.getEbeam_givenQ(Esx3, 3.4, thetab); diff --git a/Armory/testMacro.C b/Armory/testMacro.C index c7a01f9..1810b30 100644 --- a/Armory/testMacro.C +++ b/Armory/testMacro.C @@ -32,5 +32,7 @@ f_diagonal->SetLineStyle(2); f_diagonal->Draw("same");*/ + TH2F *h5 = new TH2F("Excitation Energy vs vZ", "Excitation Energy vs vZ;vZ (mm);Excitation energy (MeV)", 200, 0, 0, 200, 0, 0); //arguments are (name, title, nbinsX, xlow, xup, nbinsY, ylow, yup) + tree1->Draw("Ex:vZ >> Excitation Energy vs vZ", "sx3ID >= 0", "colz"); }