From e8a43509b97bac0131facc787bcccbadfb1329d9 Mon Sep 17 00:00:00 2001 From: James Szalkie Date: Mon, 14 Sep 2026 13:21:55 -0400 Subject: [PATCH] increase smudge --- Armory/anasenMS.cpp | 79 ++++++++++++++++------ Armory/anasen_anode_cathode_hyperboloids.h | 2 +- Armory/testMacro.C | 18 ++++- 3 files changed, 76 insertions(+), 23 deletions(-) diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index dca7cc7..bada364 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -70,7 +70,7 @@ bool IsDeadCathode(int id){ bool IsDeadSX3(int id){ static std::set dead = { - 0, 2, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23 + //0, 2, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23 }; // add dead SX3 IDs here, 0-23 1,7,9,3 return dead.count(id); } @@ -91,7 +91,7 @@ bool IsDeadSX3FrontDnChannel(int sx3ID, int chDn){ bool IsDeadSX3BackChannel(int sx3ID, int chBk){ static std::set> dead = { - {1, 10} + //{1, 10} // {sx3ID, back-channel} }; return dead.count({sx3ID, chBk}); @@ -436,6 +436,11 @@ 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_recon, vY_recon, vZ_recon; + 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"); + // reconstructed vertex Z from PW fit double z0; tree1->Branch("z0", &z0, "reconstucted_Z/D"); @@ -483,6 +488,17 @@ int main(int argc, char **argv){ detZ = TMath::QuietNaN(); hitTheta = TMath::QuietNaN(); originalEnergy = TMath::QuietNaN(); + aX = TMath::QuietNaN(); + aY = TMath::QuietNaN(); + aZ = TMath::QuietNaN(); + cX = TMath::QuietNaN(); + cY = TMath::QuietNaN(); + cZ = TMath::QuietNaN(); + theta_recon = TMath::QuietNaN(); + vX_recon = TMath::QuietNaN(); + vY_recon = TMath::QuietNaN(); + vZ_recon = TMath::QuietNaN(); + Ex_recon = TMath::QuietNaN(); sx3ID = -1; sx3Up = -1; sx3Dn = -1; @@ -511,8 +527,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); // mean and standard deviation - vertexY = gRandom->Gaus(0, vertexRangeY); // mean and standard deviation + vertexX = gRandom->Gaus(-5, vertexRangeX); // mean and standard deviation + vertexY = gRandom->Gaus(7, vertexRangeY); // mean and standard deviation //vertexZ = (vertexZRange[1]- vertexZRange[0])*gRandom->Rndm() + vertexZRange[0]; @@ -529,7 +545,8 @@ int main(int argc, char **argv){ //KEA = beamEnergy / beamA; //KEA = gRandom->Uniform(0, beamE); //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 + double beamTheta = gRandom->Gaus(0, 10 * TMath::DegToRad()); + transfer.SetIncidentEnergyAngle(KEA, beamTheta, 0); //arguments are (kinetic energy, polar angle, azimuthal angle) of the incident particle in the lab frame transfer.CalReactionConstant(); // isotropic CM direction @@ -662,14 +679,9 @@ int main(int argc, char **argv){ //TVector3 hitSigma = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L); TVector3 hitSigma = sx3->GetHitPosWithSigma(0, 0); - double hitX = std::lround(hitSigma.X()); - double hitY = std::lround(hitSigma.Y()); - double hitZ = std::lround(hitSigma.Z()); - - originalEnergy = CalculateOriginalEnergy(hitX, hitY, hitZ, - 0, 0, std::lround(vertexZ / 10) * 10, - b, "He", Esx3, - distance_sx3); + double hitX = hitSigma.X(); + double hitY = hitSigma.Y(); + double hitZ = hitSigma.Z(); double distance_A; Eanode = CalculateEnergyLoss(vertexX, vertexY, vertexZ, @@ -698,12 +710,41 @@ 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 - //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); + 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_sigma, sx3Y - aY_sigma, sx3Z - 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 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; + + originalEnergy = CalculateOriginalEnergy(hitX, hitY, hitZ, + r_rhoMin_fix.X(), r_rhoMin_fix.Y(), r_rhoMin_fix.Z(), + //0, 0, r_rhoMin_fix.Z(), + //0, 0, vertexZ, + b, "He", Esx3, + distance_sx3); + + vX_recon = r_rhoMin_fix.X(); + vY_recon = r_rhoMin_fix.Y(); + vZ_recon = r_rhoMin_fix.Z(); + + Ex_recon = apkin_27Al.getExc(originalEnergy, theta_recon); + /*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(); + Ex_recon = TMath::QuietNaN(); + originalEnergy = 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); //EBeam_Kin_3_4 = apkin_27Al.getEbeam_givenQ(Esx3, 3.4, thetab); diff --git a/Armory/anasen_anode_cathode_hyperboloids.h b/Armory/anasen_anode_cathode_hyperboloids.h index d1a88e0..9b80058 100644 --- a/Armory/anasen_anode_cathode_hyperboloids.h +++ b/Armory/anasen_anode_cathode_hyperboloids.h @@ -31,7 +31,7 @@ inline std::tuple find_PC_PathLength(const TVector3& x TVector3 anode_intersect = onesheet_hyperboloid_intersect(32.0429,301.895); TVector3 cathode_intersect = onesheet_hyperboloid_intersect(37.239045,301.895); TVector3 gw_intersect = onesheet_hyperboloid_intersect(27.712,301.895); - if(anode_intersect.Z()!=54321 && cathode_intersect.Z()!=54321) + if(anode_intersect.Z()!=54321 && cathode_intersect.Z()!=54321 && gw_intersect.Z()!=54321) return std::tuple(cathode_intersect,gw_intersect,(cathode_intersect-gw_intersect).Mag()*0.1); else return std::tuple(TVector3(0,0,0), TVector3(0,0,0), 54321); diff --git a/Armory/testMacro.C b/Armory/testMacro.C index 1810b30..ee0e470 100644 --- a/Armory/testMacro.C +++ b/Armory/testMacro.C @@ -8,6 +8,7 @@ TH2F *h2 = new TH2F("EPCvEsx3", "EPC x sin(thetab) vs Esx3;Esx3;EPC * sin(thetab * TMath::DegToRad())", 200, 0, 0, 200, 0, 0); //arguments are (name, title, nbinsX, xlow, xup, nbinsY, ylow, yup) tree1->Draw("EPC*sin(thetab * TMath::DegToRad()):Esx3 >> EPCvEsx3", "Esx3 > 0", "colz"); + /* new TCanvas(); TH2F *h3 = new TH2F( "hExReconVsVZ", @@ -23,7 +24,7 @@ 100, 0, 0, 100, -1, 0 ); - tree1->Draw("Ex:vZ>>hExVsVZ", "sx3ID >= 0", "box same"); + tree1->Draw("Ex:vZ>>hExVsVZ", "sx3ID >= 0", "box same");*/ /* new TCanvas(); tree1->Draw("originalEnergy:Tb", "sx3ID >=0", ""); @@ -31,8 +32,19 @@ f_diagonal->SetLineColor(kRed); 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)", 300, 0, 0, 300, 0, 0); //arguments are (name, title, nbinsX, xlow, xup, nbinsY, ylow, yup) + tree1->Draw("Ex_recon:vZ >> hExReconVsVZ", "sx3ID >= 0", "colz"); + 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", "colz"); + tree1->Draw("Ex:vZ", "sx3ID >=0", "same"); + c1->SaveAs("Ex_recon.png"); + + new TCanvas("c2", "Canvas", 900, 600); + tree1->Draw("theta_recon:thetab", "sx3ID >=0", ""); + c2->SaveAs("theta_recon_vs_thetab.png"); - 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"); }