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

organized labels

Browse Source
This commit is contained in:
James Szalkie 2026-09-24 10:28:13 -04:00
parent 34b9fd5f3d
commit 17bbf30c42
2 changed files with 82 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

91
Armory/anasenMS.cpp
View File

@ -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,12 +682,18 @@ 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,
aX, aY, aZ,
@ -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);

33
Armory/testMacro.C
View File

@ -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");
}