From 17bbf30c4209f400e24bce7784c3aef83bbd7185 Mon Sep 17 00:00:00 2001 From: James Szalkie Date: Thu, 24 Sep 2026 10:28:13 -0400 Subject: [PATCH] organized labels --- Armory/anasenMS.cpp | 91 ++++++++++++++++++++++++++++----------------- Armory/testMacro.C | 33 ++++++++++++---- 2 files changed, 82 insertions(+), 42 deletions(-) diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index c305af0..049d7b2 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -415,12 +415,12 @@ 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, EbeamRecon=NAN; + double EBeam_Kin_gs=NAN, EBeam_Kin_2_2=NAN, EBeam_Kin_3_4=NAN, Ex_recon=NAN, Ebeam_track=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("EbeamRecon", &EbeamRecon, "EbeamRecon/D"); + tree1->Branch("Ebeam_track", &Ebeam_track, "Ebeam_track/D"); // reconstructed angles from PW track fit, method 1 and 2 double reTheta, rePhi; @@ -431,17 +431,18 @@ 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, theta_recon; + double hitTheta, Tb_recon, sigma_a, theta_recon; tree1->Branch("hitTheta", &hitTheta, "hitTheta/D"); - tree1->Branch("originalEnergy", &originalEnergy, "originalEnergy/D"); + tree1->Branch("Tb_recon", &Tb_recon, "Tb_recon/D"); tree1->Branch("sigma_a", &sigma_a, "sigma_a/D"); tree1->Branch("theta_recon", &theta_recon, "theta_recon/D"); - double vX_recon, vY_recon, vZ_recon, vZ_from_beam; - tree1->Branch("vX_recon", &vX_recon, "vX_recon/D"); - tree1->Branch("vY_recon", &vY_recon, "vY_recon/D"); - tree1->Branch("vZ_recon", &vZ_recon, "vZ_recon/D"); - tree1->Branch("vZ_from_beam", &vZ_from_beam, "vZ_from_beam/D"); + double vX_kin, vY_kin, vZ_kin, vZ_track, EBeam_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"); // reconstructed vertex Z from PW fit double z0; @@ -489,7 +490,7 @@ int main(int argc, char **argv){ detY = TMath::QuietNaN(); detZ = TMath::QuietNaN(); hitTheta = TMath::QuietNaN(); - originalEnergy = TMath::QuietNaN(); + Tb_recon = TMath::QuietNaN(); aX = TMath::QuietNaN(); aY = TMath::QuietNaN(); aZ = TMath::QuietNaN(); @@ -497,11 +498,12 @@ int main(int argc, char **argv){ cY = TMath::QuietNaN(); cZ = TMath::QuietNaN(); theta_recon = TMath::QuietNaN(); - vX_recon = TMath::QuietNaN(); - vY_recon = TMath::QuietNaN(); - vZ_recon = TMath::QuietNaN(); + vX_kin = TMath::QuietNaN(); + vY_kin = TMath::QuietNaN(); + vZ_kin = TMath::QuietNaN(); Ex_recon = TMath::QuietNaN(); - EbeamRecon = TMath::QuietNaN(); + Ebeam_track = TMath::QuietNaN(); + EBeam_kin = TMath::QuietNaN(); sx3ID = -1; sx3Up = -1; sx3Dn = -1; @@ -680,11 +682,17 @@ int main(int argc, char **argv){ //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(0, 0); + + //smudging + TVector3 hitSigma = sx3->GetHitPosWithSigma(15,15); // 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 distance_A; Eanode = CalculateEnergyLoss(vertexX, vertexY, vertexZ, @@ -715,43 +723,58 @@ 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 - double aX_sigma = gRandom->Gaus(aX, 0); - double aY_sigma = gRandom->Gaus(aY, 0); - double aZ_sigma = gRandom->Gaus(aZ, 0); - theta_recon = TVector3(sx3X - aX, sx3Y - aY, sx3Z - aZ).Theta() * 180.0 / TMath::Pi(); + 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(sx3X, sx3Y, sx3Z); + TVector3 x1(hitX, hitY, hitZ); TVector3 v = x2f - x1; double t_minimum = -1.0*(x1.X()*v.X()+x1.Y()*v.Y())/(v.X()*v.X()+v.Y()*v.Y()); TVector3 r_rhoMin_fix = x1 + t_minimum*v; - vX_recon = r_rhoMin_fix.X(); - vY_recon = r_rhoMin_fix.Y(); - vZ_recon = r_rhoMin_fix.Z(); + vX_kin = r_rhoMin_fix.X(); + vY_kin = r_rhoMin_fix.Y(); + vZ_kin = r_rhoMin_fix.Z(); - originalEnergy = CalculateOriginalEnergy(hitX, hitY, hitZ, - vX_recon, vY_recon, vZ_recon, + Tb_recon = CalculateOriginalEnergy(hitX, hitY, hitZ, + vX_kin, vY_kin, vZ_kin, //0, 0, r_rhoMin_fix.Z(), //0, 0, vertexZ, b, "He", Esx3, distance_sx3); - Ex_recon = apkin_27Al.getExc(originalEnergy, theta_recon); - EbeamRecon = apkin_27Al.getEbeam_givenQ(originalEnergy, Ex, theta_recon); - vZ_from_beam = elossBeamInverse->Eval(EbeamRecon) - 450; + //kinematic reconstruction + Ex_recon = apkin_27Al.getExc(Tb_recon, theta_recon); + EBeam_kin = apkin_27Al.getEbeam_givenQ(Tb_recon, Ex_recon, theta_recon); //args are (reaction energy, Ex, thetab) + + /* + vZ_kin = elossBeamInverse->Eval(EBeam_kin) * 10 - 450; + + Tb_recon = CalculateOriginalEnergy(hitX, hitY, hitZ, + vX_kin, vY_kin, vZ_kin, + //0, 0, r_rhoMin_fix.Z(), + //0, 0, vertexZ, + b, "He", Esx3, + distance_sx3); + + theta_recon = TVector3(hitX - vX_kin, hitY - vY_kin, hitZ - vZ_kin).Theta() * 180.0 / TMath::Pi();*/ + + //tracking reconstruction + Ebeam_track = apkin_27Al.getEbeam_givenQ(Tb_recon, Ex, theta_recon); + vZ_track = elossBeamInverse->Eval(Ebeam_track) * 10 - 450; + /*Checklist: anode smudge, sx3 smudge, beam position off axis, beam angle*/ } else { theta_recon = TMath::QuietNaN(); - vX_recon = TMath::QuietNaN(); - vY_recon = TMath::QuietNaN(); - vZ_recon = TMath::QuietNaN(); + vX_kin = TMath::QuietNaN(); + vY_kin = TMath::QuietNaN(); + vZ_kin = TMath::QuietNaN(); Ex_recon = TMath::QuietNaN(); - EbeamRecon = TMath::QuietNaN(); - originalEnergy = TMath::QuietNaN(); - vZ_from_beam = TMath::QuietNaN(); + Ebeam_track = TMath::QuietNaN(); + Tb_recon = TMath::QuietNaN(); + vZ_track = TMath::QuietNaN(); + EBeam_kin = TMath::QuietNaN(); } //EBeam_Kin_gs = apkin_27Al.getEbeam_givenQ(Esx3, 0.0, thetab); //EBeam_Kin_2_2 = apkin_27Al.getEbeam_givenQ(Esx3, 2.2, thetab); diff --git a/Armory/testMacro.C b/Armory/testMacro.C index b7edee9..598b8a9 100644 --- a/Armory/testMacro.C +++ b/Armory/testMacro.C @@ -34,20 +34,37 @@ tree1->Draw("Ex:vZ >> h5", "sx3ID >= 0 && EPC >= 0.03", "same"); */ new TCanvas("c1", "Canvas", 900, 600); - tree1->Draw("Ex_recon:vZ_recon>>h2(800,-500,300,110,-1,10)", "sx3ID >=0 && !TMath::IsNaN(vZ_recon)", "colz"); - tree1->Draw("Ex:vZ", "sx3ID >=0 && !TMath::IsNaN(vZ_recon)", "same"); + 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"); - + /* new TCanvas("c2", "Canvas", 900, 600); tree1->Draw("theta_recon:thetab", "sx3ID >=0 && Esx3 > 0", ""); - c2->SaveAs("theta_recon_vs_thetab.png"); + c2->SaveAs("theta_recon_vs_thetab.png");*/ new TCanvas("c3", "Canvas", 900, 600); - tree1->Draw("vZ_recon:vZ", "sx3ID >=0 && !TMath::IsNaN(vZ_recon)", ""); - TF1 *line = new TF1("line", "x", gPad->GetUxmin(), gPad->GetUxmax()); - line->SetLineColor(kRed); line->Draw("SAME"); + tree1->Draw("vZ_track:vZ", "sx3ID >=0 && !TMath::IsNaN(vZ_track)", ""); + TF1 *line_track = new TF1("line_track", "x", gPad->GetUxmin(), gPad->GetUxmax()); + line_track->SetLineColor(kRed); line_track->Draw("SAME"); + c3->SaveAs("vZ_track_vs_vZ.png"); - c3->SaveAs("vZ_recon_vs_vZ.png"); + new TCanvas("c4", "Canvas", 900, 600); + tree1->Draw("vZ_kin:vZ", "sx3ID >=0 && !TMath::IsNaN(vZ_kin)", ""); + TF1 *line_kin = new TF1("line_kin", "x", gPad->GetUxmin(), gPad->GetUxmax()); + line_kin->SetLineColor(kRed); line_kin->Draw("SAME"); + c4->SaveAs("vZ_kin_vs_vZ.png"); + + new TCanvas("c5", "Canvas", 900, 600); + 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"); + + new TCanvas("c6", "Canvas", 900, 600); + tree1->Draw("EBeam_kin:beamEnergy", "sx3ID >=0 && !TMath::IsNaN(EBeam_kin)", ""); + TF1 *line_kin_kin = new TF1("line_beam_kin", "x", gPad->GetUxmin(), gPad->GetUxmax()); + line_kin_kin->SetLineColor(kRed); line_kin_kin->Draw("SAME"); + c6->SaveAs("Ebeam_kin_vs_beamEnergy.png"); }