From f2c53940f67e540675459100e5468eb28c63466a Mon Sep 17 00:00:00 2001 From: James Szalkie Date: Thu, 24 Sep 2026 13:46:04 -0400 Subject: [PATCH] proper kinematic technique --- Armory/anasenMS.cpp | 39 +++++++++++++++++++++++---------------- Armory/testMacro.C | 8 ++++---- 2 files changed, 27 insertions(+), 20 deletions(-) diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index cbf630a..98625e5 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -415,12 +415,10 @@ int main(int argc, char **argv){ tree1->Branch("qqqUp", &qqqUp, "qqqUp/I"); tree1->Branch("qqqBk", &qqqBk, "qqqBk/I"); - double EBeam_Kin_gs=NAN, EBeam_Kin_2_2=NAN, EBeam_Kin_3_4=NAN, Ex_recon=NAN, Ebeam_track=NAN; + double EBeam_Kin_gs=NAN, EBeam_Kin_2_2=NAN, EBeam_Kin_3_4=NAN; tree1->Branch("EBeam_Kin", &EBeam_Kin_gs, "EBeam_Kin/D"); tree1->Branch("EBeam_Kin_2.2", &EBeam_Kin_2_2, "EBeam_Kin_2.2/D"); tree1->Branch("EBeam_Kin_3.4", &EBeam_Kin_3_4, "EBeam_Kin_3.4/D"); - tree1->Branch("Ex_recon", &Ex_recon, "Ex_recon/D"); - tree1->Branch("Ebeam_track", &Ebeam_track, "Ebeam_track/D"); // reconstructed angles from PW track fit, method 1 and 2 double reTheta, rePhi; @@ -437,12 +435,19 @@ int main(int argc, char **argv){ tree1->Branch("sigma_a", &sigma_a, "sigma_a/D"); tree1->Branch("theta_recon", &theta_recon, "theta_recon/D"); - double vX_kin, vY_kin, vZ_kin, vZ_track, EBeam_kin; + double vX_track, vY_track, vZ_track, EBeam_track, Ex_track; + tree1->Branch("vX_track", &vX_track, "vX_track/D"); + tree1->Branch("vY_track", &vY_track, "vY_track/D"); + tree1->Branch("vZ_track", &vZ_track, "vZ_track/D"); + tree1->Branch("EBeam_track", &EBeam_track, "EBeam_track/D"); + tree1->Branch("Ex_track", &Ex_track, "Ex_track/D"); + + double vX_kin, vY_kin, vZ_kin, EBeam_kin, Ex_kin; tree1->Branch("vX_kin", &vX_kin, "vX_kin/D"); tree1->Branch("vY_kin", &vY_kin, "vY_kin/D"); tree1->Branch("vZ_kin", &vZ_kin, "vZ_kin/D"); - tree1->Branch("vZ_track", &vZ_track, "vZ_track/D"); tree1->Branch("EBeam_kin", &EBeam_kin, "EBeam_kin/D"); + tree1->Branch("Ex_kin", &Ex_kin, "Ex_kin/D"); // reconstructed vertex Z from PW fit double z0; @@ -501,8 +506,8 @@ int main(int argc, char **argv){ vX_kin = TMath::QuietNaN(); vY_kin = TMath::QuietNaN(); vZ_kin = TMath::QuietNaN(); - Ex_recon = TMath::QuietNaN(); - Ebeam_track = TMath::QuietNaN(); + Ex_kin = TMath::QuietNaN(); + EBeam_track = TMath::QuietNaN(); EBeam_kin = TMath::QuietNaN(); sx3ID = -1; sx3Up = -1; @@ -684,15 +689,15 @@ int main(int argc, char **argv){ //TVector3 hitSigma = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L); //smudging - TVector3 hitSigma = sx3->GetHitPosWithSigma(15,15); // args are sigmaSX3_W, sigmaSX3_L, currently set to 0 for no smearing + TVector3 hitSigma = sx3->GetHitPosWithSigma(10, 10); // args are sigmaSX3_W, sigmaSX3_L, currently set to 0 for no smearing double hitX = hitSigma.X(); double hitY = hitSigma.Y(); double hitZ = hitSigma.Z(); - double aX_sigma = gRandom->Gaus(aX, 15); - double aY_sigma = gRandom->Gaus(aY, 15); - double aZ_sigma = gRandom->Gaus(aZ, 15); + double aX_sigma = gRandom->Gaus(aX, 10); + double aY_sigma = gRandom->Gaus(aY, 10); + double aZ_sigma = gRandom->Gaus(aZ, 10); double distance_A; Eanode = CalculateEnergyLoss(vertexX, vertexY, vertexZ, @@ -723,7 +728,6 @@ int main(int argc, char **argv){ Kinematics apkin_27Al(26.981538408,4.00260325413,1.00782503224,29.973770136,beamEnergy/26.981538408); //m3 is proton if (!std::isnan(aX) && !std::isnan(aY) && !std::isnan(aZ) && aZ != 54321 && (aX != 0 || aY != 0 || aZ != 0)) { //reconstruct theta using anode and sx3 positions, apply a sigma to the anode locations - theta_recon = TVector3(hitX - aX_sigma, hitY - aY_sigma, hitZ - aZ_sigma).Theta() * 180.0 / TMath::Pi(); //calulate vertex like in experiment TVector3 x2f(aX_sigma, aY_sigma, aZ_sigma); TVector3 x1(hitX, hitY, hitZ); @@ -735,6 +739,8 @@ int main(int argc, char **argv){ vY_track = r_rhoMin_fix.Y(); vZ_track = r_rhoMin_fix.Z(); + EBeam_track = elossBeam->Eval((vZ_track + 450) / 10); + Tb_recon = CalculateOriginalEnergy(hitX, hitY, hitZ, vX_track, vY_track, vZ_track, //0, 0, r_rhoMin_fix.Z(), @@ -743,8 +749,9 @@ int main(int argc, char **argv){ distance_sx3); //kinematic reconstruction - Ex_recon = apkin_27Al.getExc(Tb_recon, theta_recon); - EBeam_track = apkin_27Al.getEbeam_givenQ(Tb_recon, Ex_recon, theta_recon); //args are (reaction energy, Ex, thetab) still tracking + theta_recon = TVector3(hitX - aX_sigma, hitY - aY_sigma, hitZ - aZ_sigma).Theta() * 180.0 / TMath::Pi(); + Ex_kin = apkin_27Al.getExc(Tb_recon, theta_recon); + EBeam_kin = apkin_27Al.getEbeam_givenQ(Tb_recon, Ex, theta_recon); //args are (reaction energy, Ex, thetab) still tracking vZ_kin = elossBeamInverse->Eval(EBeam_kin) * 10 - 450; /*Checklist: anode smudge, sx3 smudge, beam position off axis, beam angle*/ @@ -755,8 +762,8 @@ int main(int argc, char **argv){ vX_kin = TMath::QuietNaN(); vY_kin = TMath::QuietNaN(); vZ_kin = TMath::QuietNaN(); - Ex_recon = TMath::QuietNaN(); - Ebeam_track = TMath::QuietNaN(); + Ex_kin = TMath::QuietNaN(); + EBeam_track = TMath::QuietNaN(); Tb_recon = TMath::QuietNaN(); vZ_track = TMath::QuietNaN(); EBeam_kin = TMath::QuietNaN(); diff --git a/Armory/testMacro.C b/Armory/testMacro.C index 598b8a9..1b4e454 100644 --- a/Armory/testMacro.C +++ b/Armory/testMacro.C @@ -34,9 +34,9 @@ tree1->Draw("Ex:vZ >> h5", "sx3ID >= 0 && EPC >= 0.03", "same"); */ new TCanvas("c1", "Canvas", 900, 600); - tree1->Draw("Ex_recon:vZ_kin>>h2(800,-500,300,110,-1,10)", "sx3ID >=0 && !TMath::IsNaN(vZ_kin)", "colz"); - tree1->Draw("Ex:vZ_kin", "sx3ID >=0 && !TMath::IsNaN(vZ_kin)", "same"); - c1->SaveAs("Ex_recon.png"); + tree1->Draw("Ex_kin:vZ_kin>>h2(800,-500,300,110,-1,10)", "sx3ID >=0 && !TMath::IsNaN(vZ_kin)", "colz"); + tree1->Draw("Ex:vZ", "sx3ID >=0 && !TMath::IsNaN(vZ)", "same"); + c1->SaveAs("Ex_kin.png"); /* new TCanvas("c2", "Canvas", 900, 600); tree1->Draw("theta_recon:thetab", "sx3ID >=0 && Esx3 > 0", ""); @@ -55,7 +55,7 @@ c4->SaveAs("vZ_kin_vs_vZ.png"); new TCanvas("c5", "Canvas", 900, 600); - tree1->Draw("Ebeam_track:beamEnergy", "sx3ID >=0 && !TMath::IsNaN(Ebeam_track)", ""); + tree1->Draw("EBeam_track:beamEnergy", "sx3ID >=0 && !TMath::IsNaN(EBeam_track)", ""); TF1 *line_kin_track = new TF1("line_beam_track", "x", gPad->GetUxmin(), gPad->GetUxmax()); line_kin_track->SetLineColor(kRed); line_kin_track->Draw("SAME"); c5->SaveAs("Ebeam_track_vs_beamEnergy.png");