diff --git a/Armory/AutoHist2D.h b/Armory/AutoHist2D.h index 5e928a2..459b4ee 100644 --- a/Armory/AutoHist2D.h +++ b/Armory/AutoHist2D.h @@ -2,7 +2,11 @@ #define AUTOHIST2D_H #include +#include #include +#include +#include +#include #include #include #include @@ -60,6 +64,55 @@ public: Registry().clear(); } + // Writes the standard tree-derived plots from testMacro.C into the + // current ROOT directory, along with the requested vertex-Z residual. + static void WriteTreePlotHistograms(TTree* tree) { + if (!tree) return; + + const auto drawAndWrite = [tree](const char* drawExpression, + const char* selection, + const char* histogramName) { + tree->Draw(drawExpression, selection, "goff"); + TH1* histogram = dynamic_cast(gDirectory->Get(histogramName)); + if (histogram) { + histogram->Write("", TObject::kOverwrite); + } else { + std::cerr << "Could not create ROOT histogram: " << histogramName << std::endl; + } + }; + + drawAndWrite("Ex_kin:vZ_track>>hExKinVsVZTrack(800,-500,300,110,-1,10)", + "sx3ID >= 0 && !TMath::IsNaN(vZ_track) && EPC > 0.02", + "hExKinVsVZTrack"); + drawAndWrite("Ex:vZ>>hExVsVZ", + "sx3ID >= 0 && !TMath::IsNaN(vZ) && EPC > 0.02", + "hExVsVZ"); + drawAndWrite("theta_recon:thetab>>hThetaReconVsThetaB", + "sx3ID >= 0 && Esx3 > 0", + "hThetaReconVsThetaB"); + drawAndWrite("vZ_track:vZ>>hVZTrackVsVZ", + "sx3ID >= 0 && !TMath::IsNaN(vZ_track)", + "hVZTrackVsVZ"); + drawAndWrite("vZ_kin:vZ>>hVZKinVsVZ", + "sx3ID >= 0 && !TMath::IsNaN(vZ_kin)", + "hVZKinVsVZ"); + drawAndWrite("EBeam_track:beamEnergy>>hEBeamTrackVsBeamEnergy", + "sx3ID >= 0 && !TMath::IsNaN(EBeam_track)", + "hEBeamTrackVsBeamEnergy"); + drawAndWrite("EBeam_kin:beamEnergy>>hEBeamKinVsBeamEnergy", + "sx3ID >= 0 && !TMath::IsNaN(EBeam_kin)", + "hEBeamKinVsBeamEnergy"); + drawAndWrite("Ex_kin>>hExKin", + "sx3ID >= 0 && Ex_kin > 0", + "hExKin"); + drawAndWrite("Ex>>hEx", + "sx3ID >= 0 && Ex > 0", + "hEx"); + drawAndWrite("vZ_track-vZ>>hVZTrackMinusVZ", + "sx3ID >= 0 && EPC > 0.03", + "hVZTrackMinusVZ"); + } + private: static constexpr const char* kPlotDir = "Plots"; diff --git a/Armory/Kinematics.h b/Armory/Kinematics.h index e2cb1c0..f5d0db2 100644 --- a/Armory/Kinematics.h +++ b/Armory/Kinematics.h @@ -4,6 +4,9 @@ #include #include #include +#include +#include +#include #include #include "EnergyLoss.h" const double u_MeV = 931.49410372; //u in MeV @@ -226,5 +229,111 @@ double Kinematics::getTheta4(double t3, double angle3) { return theta4; } +double InvertToBeamEnergy(double Ex, + double Tb_recon, + double theta_recon, + const Isotope& al27Iso, + const Isotope& alphaIso) { + + if (!std::isfinite(Ex) || !std::isfinite(Tb_recon) || + !std::isfinite(theta_recon) || !std::isfinite(al27Iso.Mass) || + !std::isfinite(alphaIso.Mass) || Tb_recon < 0.0 || + theta_recon < 0.0 || theta_recon > 180.0 || + al27Iso.Mass <= 0.0 || alphaIso.Mass <= 0.0) { + std::cerr << "[InvertKin] Invalid input: Ex=" << Ex + << ", Tb_recon=" << Tb_recon + << ", theta_recon=" << theta_recon + << ", m(Al27)=" << al27Iso.Mass + << ", m(alpha)=" << alphaIso.Mass << std::endl; + return -1.0; + } + + // Known target state + double Mstar = al27Iso.Mass + Ex; + + // Detected alpha (known per event) + double E3 = Tb_recon + alphaIso.Mass; + double p3sq = E3 * E3 - alphaIso.Mass * alphaIso.Mass; + //double p3 = std::sqrt(p3sq); + double cosT = std::cos(theta_recon * TMath::Pi() / 180.0); + + // Quadratic coefficients (E1 = EBeam_track + al27Iso.Mass) + // C = m2 - E3, K = M*² - m1² + p3² - C² + double m1 = al27Iso.Mass; + double m2 = alphaIso.Mass; // target at rest + + double C = m2 - E3; + double K = Mstar*Mstar - m1*m1 + p3sq - C*C; + + double a = C*C - p3sq * cosT * cosT; + double b = -K * C; + double c = 0.25 * (K*K + 4.0 * p3sq * cosT*cosT * m1*m1); + + double disc = b*b - 4.0*a*c; + if (!std::isfinite(disc) || disc < 0.0) { + std::cerr << "[InvertKin] No real solution: discriminant=" << disc + << ", Ex=" << Ex << ", Talpha=" << Tb_recon + << ", theta=" << theta_recon << std::endl; + return -1.0; // sentinel for bad event + } + + double sqrtDisc = std::sqrt(disc); + + // Solve for total beam energy E1. If the quadratic coefficient is + // effectively zero, use the remaining linear equation instead. + std::vector candidates; + const double aScale = std::max({1.0, std::abs(C * C), + std::abs(p3sq * cosT * cosT)}); + if (std::abs(a) <= 1e-14 * aScale) { + const double bScale = std::max(1.0, std::abs(K * C)); + if (std::abs(b) <= 1e-14 * bScale) { + std::cerr << "[InvertKin] Degenerate equation: a=" << a + << ", b=" << b << ", c=" << c << std::endl; + return -1.0; + } + candidates.push_back(-c / b); + } else { + candidates.push_back((-b + sqrtDisc) / (2.0 * a)); + candidates.push_back((-b - sqrtDisc) / (2.0 * a)); + } + + // Squaring loses the sign of the momentum term, so verify each candidate + // against the unsquared equation K - 2*C*E1 = 2*p1*p3*cos(theta). + auto isPhysical = [&](double E1) -> bool { + if (!std::isfinite(E1) || E1 <= m1) return false; + const double p1sq = E1 * E1 - m1 * m1; + if (p1sq < 0.0) return false; + const double lhs = K - 2.0 * C * E1; + const double rhs = 2.0 * std::sqrt(p1sq * p3sq) * cosT; + const double residualScale = std::max({1.0, std::abs(lhs), std::abs(rhs), std::abs(K)}); + return std::abs(lhs - rhs) <= 1e-8 * residualScale; + }; + + std::vector physicalRoots; + for (const double E1 : candidates) { + if (isPhysical(E1)) physicalRoots.push_back(E1); + } + + if (physicalRoots.size() == 1) return physicalRoots.front() - m1; + if (physicalRoots.size() > 1) { + std::cerr << "[InvertKin] Ambiguous physical roots; selecting higher beam energy." + << std::endl; + return *std::max_element(physicalRoots.begin(), physicalRoots.end()) - m1; + } + + std::cerr << "[InvertKin] No physical root: Ex=" << Ex + << ", Talpha=" << Tb_recon << ", theta=" << theta_recon + << " deg, a=" << a << ", b=" << b << ", c=" << c + << ", disc=" << disc; + for (const double E1 : candidates) { + std::cerr << ", candidate Tbeam=" << E1 - m1 + << " (unsquared lhs=" << K - 2.0 * C * E1 + << ", rhs=" << (E1 > m1 ? 2.0 * std::sqrt(E1 * E1 - m1 * m1) * + std::sqrt(p3sq) * cosT : 0.0) << ")"; + } + std::cerr << std::endl; + return -1.0; +} + #endif diff --git a/Armory/Makefile b/Armory/Makefile index 7f875fe..ebabe3f 100644 --- a/Armory/Makefile +++ b/Armory/Makefile @@ -69,6 +69,6 @@ EventBuilder: EventBuilder.cpp ClassData.h fsuReader.h Hit.h @echo "--------- making EventBuilder" $(CXX) $(CXXFLAGS) EventBuilder.cpp -o EventBuilder $(LDFLAGS) -AnasenMS: anasenMS.cpp constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassQQQ.h ClassPW.h ClassAnasen.h EnergyLoss.h AutoHist2D.h anasen_anode_cathode_hyperboloids.h +AnasenMS: anasenMS.cpp Kinematics.h constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassQQQ.h ClassPW.h ClassAnasen.h EnergyLoss.h AutoHist2D.h anasen_anode_cathode_hyperboloids.h @echo "--------- making ANASEN Monte Carlo" $(CXX) $(CXXFLAGS) anasenMS.cpp -o AnasenMS $(LDFLAGS) \ No newline at end of file diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index acca090..3d0add5 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}); @@ -160,16 +160,16 @@ int main(int argc, char **argv){ //To set beam energy loss, use energy loss app, and create table with target isotope, set Initial beam energy as max energy transfer.SetA(27, 13, 0); // 22Mg projectile transfer.Seta(4, 2); // 4He target - transfer.Setb(4, 2); // outgoing proton from the primary transfer - transfer.SetB(27, 13); // 30Si* heavy product + transfer.Setb(1, 1); // outgoing proton from the primary transfer + transfer.SetB(30, 14); // 30Si* heavy product double beamE = 56.1; //56.1 const ReactionConfig reactionConfig = transfer.GetRectionConfig(); const double beamA = reactionConfig.beamA; // mass number of 14N beam // Excited state lists (projectile and heavy-product excitation states) std::vector ExAList = {0}; // Beam excited energy - //std::vector ExList = {0, 2.235, 3.498, 3.774, 4.831, 4.809, 5.279, 5.614, 5.951, 6.550}; // a/p aluminum - std::vector ExList = {0.84, 1.01, 2.21, 2.73, 2.98, 3}; + std::vector ExList = {0, 2.235, 3.498, 3.774, 4.831, 4.809, 5.279, 5.614, 5.951, 6.550}; // a/p aluminum + //std::vector ExList = {0.84, 1.01, 2.21, 2.73, 2.98, 3}; const int kMBeam = reactionConfig.beamA; // mass number of beam const int kMTarget = reactionConfig.targetA; // mass number of target @@ -203,8 +203,8 @@ int main(int argc, char **argv){ //TGraph* elossLight = LoadELoss("../ELoss/HeLoss/E_vs_x_" + b + ".dat"); // define vertex position uniform distribution ranges (mm) - double vertexXRange[2] = { -5, 5}; // mm - 5, 5 - double vertexYRange[2] = { -5, 5}; // -5, 5 + double vertexXRange[2] = { -3.5 / 2, 3.5 / 2}; // mm - 5, 5// 2.5 + double vertexYRange[2] = { -3.5 / 2, 3.5 / 2}; // -5, 5 double vertexZRange[2] = { -174.3, 174.3}; // -174.3, 174.3 (full length of gas volume, centered at 0) const double beamEntranceZ = -280 - 174.3; //vertexZRange[0]; // mm, assumed beam entrance into the gas @@ -545,8 +545,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, 0);//vertexRangeX); // mean and standard deviation - vertexY = gRandom->Gaus(0, 0);//vertexRangeY); // mean and standard deviation + vertexX = gRandom->Gaus(0, vertexRangeX); // and standard deviation + vertexY = gRandom->Gaus(0, vertexRangeY); // mean and standard deviation //vertexZ = (vertexZRange[1]- vertexZRange[0])*gRandom->Rndm() + vertexZRange[0]; @@ -697,15 +697,15 @@ int main(int argc, char **argv){ //TVector3 hitSigma = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L); //smudging - TVector3 hitSigma = sx3->GetHitPosWithSigma(0, 0); // args are sigmaSX3_W, sigmaSX3_L, currently set to 0 for no smearing + TVector3 hitSigma = sx3->GetHitPosWithSigma(2, 2); // 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, 0); - double aY_sigma = gRandom->Gaus(aY, 0); - double aZ_sigma = gRandom->Gaus(aZ, 0); + double aX_sigma = gRandom->Gaus(aX, 2.3); + double aY_sigma = gRandom->Gaus(aY, 2.3); + double aZ_sigma = gRandom->Gaus(aZ, 2.3); double distance_A; Eanode = CalculateEnergyLoss(vertexX, vertexY, vertexZ, @@ -723,6 +723,7 @@ int main(int argc, char **argv){ Esx3 = NAN; beamEnergy = NAN; Ex = NAN; + sx3ID = -1; //continue; } sx3Tb = Tb; // for simplicity, using the same kinetic energy for SX3 hit events, can be modified to simulate energy loss if desired @@ -759,9 +760,8 @@ int main(int argc, char **argv){ const TVector3 projectileDirection = TVector3(hitX - aX_sigma, hitY - aY_sigma, hitZ - aZ_sigma); theta_recon = projectileDirection.Theta() * 180.0 / TMath::Pi(); - if (reactionConfig.recoilLightA == 1) { - //Kinematics apkin_27Al(26.981538408,4.00260325413,1.00782503224,29.973770136,EBeam_track/26.981538408); //m3 is proton - Kinematics apkin_27Al(26.981538408,4.00260325413,4.00260325413,26.981538408,EBeam_track/26.981538408); //args are (m1, m2, m3, m4, Ebeam/m1) + if (reactionConfig.recoilLightA == 1 && !TMath::IsNaN(EBeam_track)) { + Kinematics apkin_27Al(26.981538408,4.00260325413,1.00782503224,29.973770136,EBeam_track/26.981538408); //m3 is proton //kinematic reconstruction EBeam_kin = apkin_27Al.getEbeam_givenQ(Tb_recon, Ex, theta_recon); //args are (reaction energy, Ex, thetab) still tracking @@ -783,11 +783,21 @@ int main(int argc, char **argv){ alphaMomentum * std::cos(theta_recon * TMath::Pi() / 180.0), alphaEnergy); TLorentzVector recoil4 = beam4 + helium4 - alpha4; - double recoilMass = recoil4.M(); Ex_kin = recoilMass - al27Iso.Mass; - vZ_kin = vZ_track; + + if (std::isfinite(Ex) && std::isfinite(Tb_recon) && + std::isfinite(theta_recon)) { + EBeam_kin = InvertToBeamEnergy(Ex, Tb_recon, theta_recon, + al27Iso, alphaIso); + } else { + EBeam_kin = TMath::QuietNaN(); + } + + vZ_kin = std::isfinite(EBeam_kin) && EBeam_kin >= 0.0 + ? elossBeamInverse->Eval(EBeam_kin) * 10 - 454.3 + : TMath::QuietNaN(); } @@ -917,6 +927,8 @@ int main(int argc, char **argv){ // write results to ROOT file and close AutoHist2D::WriteAll(); // books, fills and writes all registered 2D histograms into saveFile + saveFile->cd(); + AutoHist2D::WriteTreePlotHistograms(tree1); tree1->Write("", TObject::kOverwrite); int count1 = tree1->GetEntries(); //int count3 = tree3->GetEntries(); diff --git a/Armory/anasen_anode_cathode_hyperboloids.h b/Armory/anasen_anode_cathode_hyperboloids.h index 08843d6..2c094d6 100644 --- a/Armory/anasen_anode_cathode_hyperboloids.h +++ b/Armory/anasen_anode_cathode_hyperboloids.h @@ -28,8 +28,8 @@ inline std::tuple find_PC_PathLength(const TVector3& x //TODO: Magic numbers here describing waist 'a', and flare 'c' will need better treatment. //Currently, these are derived by fitting the crossover points to R^2/a^2 - z^2/c^2 = 1 for anodes // Cathode a, c values are found by scaling up the anode waist by 43/37, the ratio of the outermost radii - TVector3 anode_intersect = onesheet_hyperboloid_intersect(32.0429,301.895); - TVector3 cathode_intersect = onesheet_hyperboloid_intersect(37.239045,301.895); + TVector3 anode_intersect = onesheet_hyperboloid_intersect(32.0429,301.895); // anode hyperboloid parameters: a=32.0429, c=301.895 + TVector3 cathode_intersect = onesheet_hyperboloid_intersect(37.239045,301.895); // cathode hyperboloid parameters: a=37.239045, c=301.895 if(anode_intersect.Z()!=54321 && cathode_intersect.Z()!=54321) return std::tuple(cathode_intersect,anode_intersect,(cathode_intersect-anode_intersect).Mag()*0.1); else diff --git a/Armory/efficiency study 8-27.zip b/Armory/efficiency study 8-27.zip deleted file mode 100644 index 5951916..0000000 Binary files a/Armory/efficiency study 8-27.zip and /dev/null differ diff --git a/Armory/testMacro.C b/Armory/testMacro.C index b306304..3aca5e0 100644 --- a/Armory/testMacro.C +++ b/Armory/testMacro.C @@ -4,18 +4,45 @@ TTree *tree1 = (TTree*)(file0->Get("tree1")); new TBrowser(); + /* + new TCanvas(); + TH2F *h3 = new TH2F( + "hExReconVsVZ", + "Ex_recon vs vZ;vZ (mm);Excitation energy (MeV)", + 100, 0, 0, + 100, -1, 0 + ); + tree1->Draw("Ex_recon:vZ>>hExReconVsVZ", "sx3ID >= 0", "colz"); + TH2F *h4 = new TH2F( "hExVsVZ", "", 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", ""); + TF1 *f_diagonal = new TF1("f_diagonal", "x", 0, 200); // assuming the range for Tb is 0 to 200 + 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_kin:vZ_kin>>h2(800,-500,300,110,-1,10)", "sx3ID >=0 && !TMath::IsNaN(vZ_kin) && EPC > .02", "colz"); - tree1->Draw("Ex:vZ", "sx3ID >=0 && !TMath::IsNaN(vZ) && EPC > .02", "same"); + tree1->Draw("Ex_kin:vZ_track>>h2(800,-500,300,110,-1,10)", "sx3ID >=0 && EPC > .02", "colz"); + //tree1->Draw("Ex:vZ", "sx3ID >=0 && EPC > .02", "same"); c1->SaveAs("Ex_kin_vz.png"); + //{0, 2.235, 3.498, 3.774, 4.831, 4.809, 5.279, 5.614, 5.951, 6.550}; +/* + new TCanvas("c2", "Canvas", 900, 600); + tree1->Draw("theta_recon:thetab", "sx3ID >=0 && Esx3 > 0", ""); + c2->SaveAs("theta_recon_vs_thetab.png"); new TCanvas("c3", "Canvas", 900, 600); tree1->Draw("vZ_track:vZ", "sx3ID >=0 && !TMath::IsNaN(vZ_track)", ""); @@ -39,15 +66,16 @@ 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"); + c6->SaveAs("Ebeam_kin_vs_beamEnergy.png");*/ -/* new TCanvas("c7", "Canvas", 900, 600); //c7->SetLogy(); - tree1->Draw("Ex_kin", "sx3ID >= 0"); - tree1->Draw("Ex", "sx3ID >= 0", "same"); + //tree1->Draw("Ex_kin", "sx3ID >= 0 && EPC > .04"); + tree1->Draw("Ex", "sx3ID >= 0 && EPC > .02"); tree1->GetHistogram()->SetLineColor(kRed); + tree1->GetHistogram()->SetLineStyle(2); + tree1->Draw("Ex_kin", "sx3ID >= 0 && EPC > .02", "same"); c7->SaveAs("Ex_kin.png"); -*/ + }