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

square PC

Browse Source
This commit is contained in:
James Szalkie 2026-09-15 16:18:27 -04:00
parent e8a43509b9
commit d4e7c18040
3 changed files with 114 additions and 51 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
Armory/anasenMS.cpp
View File

@ -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,17 +723,17 @@ 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*/

129
Armory/anasen_anode_cathode_hyperboloids.h
View File

@ -1,40 +1,107 @@
inline std::tuple<TVector3,TVector3,double> 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
/*
// 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);
}
};
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);
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);
}

10
Armory/testMacro.C
View File

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