modified: Armory/ClassPW.h

modified:   Armory/PC_StepLadder_Correction.h incorp[orating 4 wire twist]
	modified:   MakeVertex.C

Armory/ClassPW.h

@@ -64,10 +64,9 @@ public:
 
   Coord Crossover[24][24][2];
 
+  inline TVector3 getClosestWirePosAtWirePhi(std::pair<TVector3, TVector3>, double phi);
   inline std::tuple<std::pair<TVector3, TVector3>, double, double, double> GetPseudoWire(const std::vector<std::tuple<int, double, double>> &cluster, std::string type);
 
-  inline TVector3 getClosestWirePosAtWirePhi(std::pair<TVector3, TVector3> awire, double sx3phi_radian);
-
   inline std::tuple<TVector3, double, double, double, double, double, double, double>
   FindCrossoverProperties(const std::vector<std::tuple<int, double, double>> &a_cluster, const std::vector<std::tuple<int, double, double>> &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<TVector3, TVector3> 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<std::vector<std::tuple<int, double, double>>>
 PW::Make_Clusters(std::unordered_map<int, std::tuple<int, double, double>> wireEvents)
 {
@@ -356,40 +398,6 @@ PW::GetPseudoWire(const std::vector<std::tuple<int, double, double>> &cluster, s
   return std::tuple(avgvec, sumEnergy, maxEnergy, tsMaxEnergy);
 }
 
-inline TVector3 PW::getClosestWirePosAtWirePhi(std::pair<TVector3, TVector3> 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<TVector3, double, double, double, double, double, double, double> PW::FindCrossoverProperties(const std::vector<std::tuple<int, double, double>> &a_cluster,
                                                                                                                 const std::vector<std::tuple<int, double, double>> &c_cluster)
 {

Armory/PC_StepLadder_Correction.h

@@ -1,36 +1,83 @@
 #include <TF1.h>
-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 f1a("model_inv",model_a1c1,-200,200,2);
+    eqline.SetNpx(10000);
+    f1a.SetNpx(10000);
+    std::vector<double> 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<double> pars = {0.0,1.};
+    //std::vector<double> 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());
-
 }*/

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);
 
-                // ==========================================================
+                double z_entrance = -274.3;
-                // DYNAMIC BEAM ENERGY CALCULATION (MUST BE HERE!)
-                // ==========================================================
-                double z_entrance = -174.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);
 
-                // ==========================================================
+                double z_entrance = -274.3;
-                // DYNAMIC BEAM ENERGY CALCULATION (MUST BE HERE TOO!)
-                // ==========================================================
-                double z_entrance = -174.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(