diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index bada364..e13526a 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -162,7 +162,7 @@ int main(int argc, char **argv){ transfer.Seta(4, 2); // 4He target transfer.Setb(1, 1); // outgoing proton from the primary transfer transfer.SetB(30, 14); // 30Si* heavy product - double beamE = 56.1; + double beamE = 100; //56.1 const ReactionConfig reactionConfig = transfer.GetRectionConfig(); const double beamA = reactionConfig.beamA; // mass number of 14N beam @@ -527,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(-5, vertexRangeX); // mean and standard deviation - vertexY = gRandom->Gaus(7, vertexRangeY); // mean and standard deviation + vertexX = gRandom->Gaus(0, vertexRangeX); // mean and standard deviation + vertexY = gRandom->Gaus(0, vertexRangeY); // mean and standard deviation //vertexZ = (vertexZRange[1]- vertexZRange[0])*gRandom->Rndm() + vertexZRange[0]; @@ -545,8 +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)); - 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 + //double beamTheta = gRandom->Gaus(0, 10 * TMath::DegToRad()); + transfer.SetIncidentEnergyAngle(KEA, 0, 0); //arguments are (kinetic energy, polar angle, azimuthal angle) of the incident particle in the lab frame transfer.CalReactionConstant(); // isotropic CM direction @@ -715,7 +715,7 @@ int main(int argc, char **argv){ 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(); + theta_recon = TVector3(sx3X - aX, sx3Y - aY, sx3Z - aZ).Theta() * 180.0 / TMath::Pi(); //calulate vertex like in experiment TVector3 x2f(aX_sigma, aY_sigma, aZ_sigma); TVector3 x1(sx3X, sx3Y, sx3Z); @@ -723,16 +723,16 @@ int main(int argc, char **argv){ 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(); + originalEnergy = CalculateOriginalEnergy(hitX, hitY, hitZ, - r_rhoMin_fix.X(), r_rhoMin_fix.Y(), r_rhoMin_fix.Z(), + vX_recon, vY_recon, vZ_recon, //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*/ diff --git a/Armory/anasen_anode_cathode_hyperboloids.h b/Armory/anasen_anode_cathode_hyperboloids.h index 9b80058..96c1c60 100644 --- a/Armory/anasen_anode_cathode_hyperboloids.h +++ b/Armory/anasen_anode_cathode_hyperboloids.h @@ -1,40 +1,107 @@ inline std::tuple find_PC_PathLength(const TVector3& x1, const TVector3& x2) { /* - Function that finds the path length between anode and cathode surfaces, both one-sheet hyperboloids of form (x*x+y*y)/(a*a) - (z*z)/(c*c) = 1 - * path length found for a given particle moving along a certain direction from x1 to x2 - * Typical arguments here will be x1=r_rhoMin, x2=qqqevent.pos or sx3event.pos + Square-aperture detector geometry replacing the original circular one-sheet hyperboloid. + The active surfaces are centered on the z axis and form squares rather than circles: + - anodes: 40 mm x 40 mm + - cathodes: 42 mm x 42 mm + The beam is kept on-axis, so the square openings enclose the beam while still defining + the anode/cathode boundaries in the x-y plane. */ - TVector3 dx = x2-x1;// direction vector - double t2 = 1.0; //The value of 't' at the destination point, by definition: t=(z(t)-z0)/dz - auto onesheet_hyperboloid_intersect = [&](double a, double c) { - auto A = pow(dx.Perp(),2)/(a*a) - pow(dx.Z(),2)/(c*c); - auto B = 2*(dx.X()*x1.X()+dx.Y()*x1.Y())/(a*a) - 2*(dx.Z()*x1.Z())/(c*c); - auto C = pow(x1.Perp(),2)/(a*a) - pow(x1.Z(),2)/(c*c) - 1.0; - double disc = B*B - 4*A*C; - if(disc<0) - return TVector3(0,0,54321); - else { - double tsol1 = (-B + TMath::Sqrt(disc))/(2*A); - double tsol2 = (-B - TMath::Sqrt(disc))/(2*A); - if(tsol1 >= 0 && tsol1 <= t2) - return x1+tsol1*dx; - else if(tsol2>=0 && tsol2 <= t2) - return x1+tsol2*dx; - else - return TVector3(0,0,54321); - } - }; - //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 gw_intersect = onesheet_hyperboloid_intersect(27.712,301.895); - 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); + /* + // Original circular hyperboloid implementation retained for reference. + // This section is intentionally commented out to avoid using the circular geometry. + TVector3 dx = x2-x1; // direction vector + double t2 = 1.0; //The value of 't' at the destination point, by definition: t=(z(t)-z0)/dz + auto onesheet_hyperboloid_intersect = [&](double a, double c) { + auto A = pow(dx.Perp(),2)/(a*a) - pow(dx.Z(),2)/(c*c); + auto B = 2*(dx.X()*x1.X()+dx.Y()*x1.Y())/(a*a) - 2*(dx.Z()*x1.Z())/(c*c); + auto C = pow(x1.Perp(),2)/(a*a) - pow(x1.Z(),2)/(c*c) - 1.0; + double disc = B*B - 4*A*C; + if(disc<0) + return TVector3(0,0,54321); + else { + double tsol1 = (-B + TMath::Sqrt(disc))/(2*A); + double tsol2 = (-B - TMath::Sqrt(disc))/(2*A); + if(tsol1 >= 0 && tsol1 <= t2) + return x1+tsol1*dx; + else if(tsol2>=0 && tsol2 <= t2) + return x1+tsol2*dx; + else + return TVector3(0,0,54321); + } + }; + + 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 && 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); + */ + + const double anode_half_width = 20.0; // 40 mm square side length + const double cathode_half_width = 21.0; // 42 mm square side length + const double guard_half_width = 18.0; // square guard opening centered on beam axis, still encloses beam + const double z_start = -100.0; // square detector opening begins here along the beam axis + const double invalid_z = 54321.0; + const double eps = 1.0e-12; + + TVector3 dx = x2 - x1; + + auto square_intersection = [&](const TVector3& start, const TVector3& direction, double half_width) { + double best_t = 1.0; + bool found = false; + + auto consider_t = [&](double t, bool x_hit, bool y_hit) { + if (t < -eps || t > 1.0 + eps) return; + + TVector3 point = start + t * direction; + double abs_x = std::abs(point.X()); + double abs_y = std::abs(point.Y()); + + bool inside_square = (abs_x <= half_width + eps) && (abs_y <= half_width + eps); + bool on_boundary = (x_hit && std::abs(abs_x - half_width) <= eps) || (y_hit && std::abs(abs_y - half_width) <= eps); + bool beyond_start = point.Z() >= z_start - eps; + + if (inside_square && on_boundary && beyond_start && t >= 0.0 && t <= 1.0) { + if (!found || t < best_t) { + best_t = t; + found = true; + } + } + }; + + if (std::abs(direction.X()) > eps) { + double t_plus_x = (half_width - start.X()) / direction.X(); + double t_minus_x = (-half_width - start.X()) / direction.X(); + consider_t(t_plus_x, true, false); + consider_t(t_minus_x, true, false); + } + + if (std::abs(direction.Y()) > eps) { + double t_plus_y = (half_width - start.Y()) / direction.Y(); + double t_minus_y = (-half_width - start.Y()) / direction.Y(); + consider_t(t_plus_y, false, true); + consider_t(t_minus_y, false, true); + } + + if (!found) { + return TVector3(0, 0, invalid_z); + } + + return start + best_t * direction; + }; + + TVector3 anode_intersect = square_intersection(x1, dx, anode_half_width); + TVector3 cathode_intersect = square_intersection(x1, dx, cathode_half_width); + TVector3 gw_intersect = square_intersection(x1, dx, guard_half_width); + + if (anode_intersect.Z() != invalid_z && cathode_intersect.Z() != invalid_z && gw_intersect.Z() != invalid_z) + 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), invalid_z); } diff --git a/Armory/testMacro.C b/Armory/testMacro.C index ee0e470..211b073 100644 --- a/Armory/testMacro.C +++ b/Armory/testMacro.C @@ -3,10 +3,6 @@ file0->ls(); TTree *tree1 = (TTree*)(file0->Get("tree1")); new TBrowser(); - tree1->Draw("Tb:thetab","","col"); - tree1->Draw("Tb:thetab","vZ>-140 && vZ<-130", "box same"); - 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(); @@ -38,12 +34,12 @@ 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"); + 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"); c1->SaveAs("Ex_recon.png"); new TCanvas("c2", "Canvas", 900, 600); - tree1->Draw("theta_recon:thetab", "sx3ID >=0", ""); + tree1->Draw("theta_recon:thetab", "sx3ID >=0 && Esx3 > 0", ""); c2->SaveAs("theta_recon_vs_thetab.png");