From 574506acb72c6045635a1e88dbe632c754f8284c Mon Sep 17 00:00:00 2001 From: vsitaraman Date: Fri, 15 May 2026 11:14:52 -0400 Subject: [PATCH] modified: Armory/ClassPW.h modified: Armory/PC_StepLadder_Correction.h incorp[orating 4 wire twist] modified: MakeVertex.C --- Armory/ClassPW.h | 84 +++++++++++++++++-------------- Armory/PC_StepLadder_Correction.h | 62 ++++++++++++++++++++--- MakeVertex.C | 11 +--- 3 files changed, 102 insertions(+), 55 deletions(-) diff --git a/Armory/ClassPW.h b/Armory/ClassPW.h index de9fde9..1db9d69 100755 --- a/Armory/ClassPW.h +++ b/Armory/ClassPW.h @@ -64,10 +64,9 @@ public: Coord Crossover[24][24][2]; + inline TVector3 getClosestWirePosAtWirePhi(std::pair, double phi); inline std::tuple, double, double, double> GetPseudoWire(const std::vector> &cluster, std::string type); - inline TVector3 getClosestWirePosAtWirePhi(std::pair awire, double sx3phi_radian); - inline std::tuple FindCrossoverProperties(const std::vector> &a_cluster, const std::vector> &c_cluster); @@ -117,7 +116,7 @@ private: TVector3 trackVec; const int nWire = 24; - const int wireShift = 3; + const int wireShift = 4; // const float zLen = 380; // mm // const float zLen = 348.6; // mm const float zLen = 174.3 * 2; // mm @@ -157,6 +156,8 @@ inline void PW::ConstructGeo() double k = TMath::TwoPi() / 24.; // 48 solder thru holes, wires in every other one double offset_a1 = -6 * k - 3 * k; double offset_c1 = -4 * k - 2 * k - TMath::TwoPi() / 48; // correct for a half-turn + // std::cerr << "Here!" << std::endl; + // #include "../scratch/testing.h" double offset_a2 = offset_a1 + wireShift * k; double offset_c2 = offset_c1 - wireShift * k; @@ -208,9 +209,14 @@ inline void PW::ConstructGeo() Crossover[i][j][0].y = An[i].second.Y() + alpha * a.Y(); Crossover[i][j][0].z = An[i].second.Z() + alpha * a.Z(); - if (Crossover[i][j][0].z < -190 || Crossover[i][j][0].z > 190 || (i + j) % 24 == 12) + if (Crossover[i][j][0].z < -190 || Crossover[i][j][0].z > 190) + { + // std::cout << "Weird crossover but ok" << std::endl; + } + if ((i + j) % 24 == 12 || Crossover[i][j][0].z < -190 || Crossover[i][j][0].z > 190) { Crossover[i][j][0].z = 9999999; + // std::cout << "Weird crossover" << std::endl; } Crossover[i][j][1].x = alpha; @@ -223,6 +229,42 @@ inline void PW::ConstructGeo() cathodeLength = TMath::Sqrt(zLen * zLen + TMath::Power(2 * radiusC * TMath::Sin(dAngle / 2), 2)); // chord length subtending an angle alpha is 2rsin(alpha/2) } +inline TVector3 PW::getClosestWirePosAtWirePhi(std::pair awire, double sx3phi_radian) +{ + // 1. Get wire geometry + TVector3 a1 = awire.first; // Top of the wire + TVector3 a2 = awire.second; // Bottom of the wire + TVector3 wireVec = a2 - a1; // Vector pointing down the wire + + // Variables to track our minimums during the scan + double min_delta_phi = 9999.0; + double best_t = -1.0; + TVector3 best_pcz_intersect; + + // 2. THE SCAN: Walk down the wire in 1000 tiny steps + // (For a 380mm wire, this is checking every 0.38 mm) + int num_steps = 1000; + for (int i = 0; i <= num_steps; ++i) + { + double t_test = (double)i / num_steps; // Ranges from 0.0 to 1.0 + TVector3 test_pt = a1 + t_test * wireVec; // The 3D point at this step + + // Calculate absolute Delta Phi between Si hit and this specific point on the wire + if (TMath::IsNaN(sx3phi_radian - test_pt.Phi())) + continue; + double dPhi = TMath::Abs(TVector2::Phi_mpi_pi(sx3phi_radian - test_pt.Phi())); // Phi_mpi_pi just puts the angle in the range -180 to 180 + + // If this is the smallest Delta Phi we've seen so far, save it! + if (dPhi < min_delta_phi) + { + min_delta_phi = dPhi; + best_t = t_test; + best_pcz_intersect = test_pt; + } + } + return best_pcz_intersect; +} + inline std::vector>> PW::Make_Clusters(std::unordered_map> wireEvents) { @@ -356,40 +398,6 @@ PW::GetPseudoWire(const std::vector> &cluster, s return std::tuple(avgvec, sumEnergy, maxEnergy, tsMaxEnergy); } -inline TVector3 PW::getClosestWirePosAtWirePhi(std::pair awire, double sx3phi_radian) -{ - // 1. Get wire geometry - TVector3 a1 = awire.first; // Top of the wire - TVector3 a2 = awire.second; // Bottom of the wire - TVector3 wireVec = a2 - a1; // Vector pointing down the wire - - // Variables to track our minimums during the scan - double min_delta_phi = 9999.0; - double best_t = -1.0; - TVector3 best_pcz_intersect; - - // 2. THE SCAN: Walk down the wire in 1000 tiny steps - // (For a 380mm wire, this is checking every 0.38 mm) - int num_steps = 1000; - for (int i = 0; i <= num_steps; ++i) - { - double t_test = (double)i / num_steps; // Ranges from 0.0 to 1.0 - TVector3 test_pt = a1 + t_test * wireVec; // The 3D point at this step - - // Calculate absolute Delta Phi between Si hit and this specific point on the wire - double dPhi = TMath::Abs(TVector2::Phi_mpi_pi(sx3phi_radian - test_pt.Phi())); // Phi_mpi_pi just puts the angle in the range -180 to 180 - - // If this is the smallest Delta Phi we've seen so far, save it! - if (dPhi < min_delta_phi) - { - min_delta_phi = dPhi; - best_t = t_test; - best_pcz_intersect = test_pt; - } - } - return best_pcz_intersect; -} - inline std::tuple PW::FindCrossoverProperties(const std::vector> &a_cluster, const std::vector> &c_cluster) { diff --git a/Armory/PC_StepLadder_Correction.h b/Armory/PC_StepLadder_Correction.h index 1740743..6618596 100644 --- a/Armory/PC_StepLadder_Correction.h +++ b/Armory/PC_StepLadder_Correction.h @@ -1,36 +1,83 @@ #include -double model(double* x, double* p) { +/*double model(double* x, double* p) { double result = x[0]; double factor = 29.0; double slope = 0.7; if(TMath::Abs(x[0]) < 16.2) result=x[0]*slope; else if(TMath::Abs(x[0]) < 49.8 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor; else if(TMath::Abs(x[0]) < 85.2 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*2; - else result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*2; + else result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*3; return result; } double model_invert(double *y, double *q) { double result=y[0]; double slope = 0.7; - double factor = 40.0; + double factor = 0.0; if(TMath::Abs(y[0]) < 16.2/slope) result = y[0]/slope; else if(TMath::Abs(y[0]) < 49.8/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor; else if(TMath::Abs(y[0]) < 85.2/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor*2; - else result=y[0]/slope-TMath::Sign(1.0,y[0])*factor*2; + else result=y[0]/slope-TMath::Sign(1.0,y[0])*factor*3; + return result; +}*/ + +double model_invert(double* y, double* p) { + double result = y[0]; + double slope = 0.6; + double z_grid[8] = {147.998,101.946,59.7634,19.6965,-19.6965,-59.7634,-101.946,-147.998}; + for(int i=0;i<7;i++) { + if(y[0] <= z_grid[i] && y[0] > z_grid[i+1]) { + double zavg = (z_grid[i] + z_grid[i+1])*0.5; //midpoint about which we pivot + result = (y[0]-zavg)/slope + zavg; + break; + } + } + return result+80; +} + +double model_a1c1(double* x, double* p) { + double result = x[0]; + double factor = 29.0; + double slope = 0.0; + if(TMath::Abs(x[0]) < 16.2) result=x[0]*slope; + else if(TMath::Abs(x[0]) < 49.8 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor; + else if(TMath::Abs(x[0]) < 85.2 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*2; + else result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*3; + return result; +} + +double model_invert_a1c1(double *y, double *q) { + double result=y[0]; +/* double slope = 1.0; + double factor = 5.0; + if(TMath::Abs(y[0]) < 16.2/slope) result = y[0]/slope; + else if(TMath::Abs(y[0]) < 49.8/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor; + else if(TMath::Abs(y[0]) < 85.2/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor*2; + else result=y[0]/slope-TMath::Sign(1.0,y[0])*factor**/; return result+40; } + /*void testmodel() { TF1 eqline("x","x",-200,200); eqline.Draw(""); eqline.SetLineStyle(kDashed); - //TF1 f1("model",model,-200,200,2); - TF1 f1("model_inv",model_invert,-200,200,2); + //TF1 f1("model",model,-200,200,2); + TF1 f1a("model_inv",model_a1c1,-200,200,2); + eqline.SetNpx(10000); + f1a.SetNpx(10000); + std::vector pars = {0.0,1.}; + f1a.SetParameters(pars.data()); + f1a.SetLineColor(kGreen+2); + f1a.SetLineStyle(kLine); + f1a.Draw("L SAME"); + + TF1 f1("model",model,-200,200,2); + //TF1 f1("model_inv",model_invert,-200,200,2); eqline.SetNpx(10000); f1.SetNpx(10000); - std::vector pars = {0.0,1.}; + //std::vector pars = {0.0,1.}; f1.SetParameters(pars.data()); f1.SetLineColor(kGreen+2); f1.SetLineStyle(kLine); @@ -38,5 +85,4 @@ double model_invert(double *y, double *q) { gPad->Modified(); gPad->Update(); while(gPad->WaitPrimitive()); - }*/ diff --git a/MakeVertex.C b/MakeVertex.C index 6c1e048..7a0d97a 100755 --- a/MakeVertex.C +++ b/MakeVertex.C @@ -1296,13 +1296,9 @@ Bool_t MakeVertex::Process(Long64_t entry) // 3. Reconstruct Vertex Z double deltaRho = sx3event.pos.Perp() - pcz_intersect.Perp(); double deltaZ = sx3event.pos.Z() - pcz_intersect.Z(); - double vertex_recon = sx3event.pos.Z() - sx3event.pos.Perp() * (deltaZ / deltaRho); - // ========================================================== - // DYNAMIC BEAM ENERGY CALCULATION (MUST BE HERE!) - // ========================================================== - double z_entrance = -174.3; + double z_entrance = -274.3; double beam_path_length = TMath::Abs(vertex_recon - z_entrance) * 0.1; double initial_beam_energy = 72.0; double beam_energy_at_vertex = cm_to_MeV_beam->Eval( @@ -1401,10 +1397,7 @@ Bool_t MakeVertex::Process(Long64_t entry) double vertex_recon = qqqevent.pos.Z() - qqqevent.pos.Perp() * (deltaZ / deltaRho); - // ========================================================== - // DYNAMIC BEAM ENERGY CALCULATION (MUST BE HERE TOO!) - // ========================================================== - double z_entrance = -174.3; + double z_entrance = -274.3; double beam_path_length = TMath::Abs(vertex_recon - z_entrance) * 0.1; double initial_beam_energy = 72.0; double beam_energy_at_vertex = cm_to_MeV_beam->Eval(