Fully removed polynomial magnetic field

martenole · davidrohr · commit fe2c48f131a6 · 2020-11-04T21:37:41.000+01:00
- assuming constant Bz: entropy constant, 7% speedup
- removed function for propagation with bxbybz
diff --git a/GPU/GPUTracking/DataCompression/GPUTPCCompressionTrackModel.cxx b/GPU/GPUTracking/DataCompression/GPUTPCCompressionTrackModel.cxx
@@ -126,7 +126,7 @@ GPUd() void GPUTPCCompressionTrackModel::Init(float x, float y, float z, float a
   mP[4] = (qPt - 127.f) * (20.f / 127.f);
   resetCovariance();
   mNDF = -5;
-  mField = &param.polynomialField;
+  mBz = param.ConstBz;
   float pti = CAMath::Abs(mP[4]);
   if (pti < 1.e-4f) {
     pti = 1.e-4f; // set 10.000 GeV momentum for straight track
@@ -163,22 +163,20 @@ GPUd() int GPUTPCCompressionTrackModel::Propagate(float x, float alpha)
     return 0;
   }
 
-  float Bz = getBz(mTrk.x, mTrk.y, mTrk.z);
-
   // propagate mTrk to t0e
   PhysicalTrackModel t0e(mTrk);
   float dLp = 0;
   if (CAMath::Abs(x - t0e.x) < 1.e-8f) {
     return 0;
   }
-  if (propagateToXBzLightNoUpdate(t0e, x, Bz, dLp)) {
+  if (propagateToXBzLightNoUpdate(t0e, x, mBz, dLp)) {
     return 1;
   }
   updatePhysicalTrackValues(t0e);
   if (CAMath::Abs(t0e.sinphi) >= MaxSinPhi) {
     return -3;
   }
-  return followLinearization(t0e, Bz, dLp);
+  return followLinearization(t0e, mBz, dLp);
 }
 
 GPUd() int GPUTPCCompressionTrackModel::Filter(float y, float z, int iRow)
@@ -294,11 +292,7 @@ GPUd() int GPUTPCCompressionTrackModel::Filter(float y, float z, int iRow)
 
 GPUd() int GPUTPCCompressionTrackModel::Mirror()
 {
-  float Bz = getBz(mTrk.x, mTrk.y, mTrk.z);
-  if (CAMath::Abs(Bz) < 1.e-8f) {
-    Bz = 1.e-8f;
-  }
-  float dy = -2.f * mTrk.q * mTrk.px / Bz;
+  float dy = -2.f * mTrk.q * mTrk.px / mBz;
   float dS; // path in XY
   {
     float chord = dy;        // chord to the extrapolated point == |dy|*sign(x direction)
@@ -361,25 +355,6 @@ GPUd() int GPUTPCCompressionTrackModel::Mirror()
   return 0;
 }
 
-GPUd() void GPUTPCCompressionTrackModel::getBxByBz(float cosAlpha, float sinAlpha, float x, float y, float z, float b[3]) const
-{
-  float xGlb = x * cosAlpha - y * sinAlpha;
-  float yGlb = x * sinAlpha + y * cosAlpha;
-  float bb[3];
-  mField->GetField(xGlb, yGlb, z, bb);
-  // rotate field to local coordinates
-  b[0] = bb[0] * cosAlpha + bb[1] * sinAlpha;
-  b[1] = -bb[0] * sinAlpha + bb[1] * cosAlpha;
-  b[2] = bb[2];
-}
-
-GPUd() float GPUTPCCompressionTrackModel::getBz(float x, float y, float z) const
-{
-  float xGlb = x * mCosAlpha - y * mSinAlpha;
-  float yGlb = x * mSinAlpha + y * mCosAlpha;
-  return mField->GetFieldBz(xGlb, yGlb, z);
-}
-
 GPUd() void GPUTPCCompressionTrackModel::updatePhysicalTrackValues(PhysicalTrackModel& trk)
 {
   float px = trk.px;
@@ -460,12 +435,11 @@ GPUd() int GPUTPCCompressionTrackModel::rotateToAlpha(float newAlpha)
   float trackX = x0 * cc + ss * mP[0];
 
   // transport t0 to trackX
-  float B[3];
-  getBxByBz(CAMath::Cos(newAlpha), CAMath::Sin(newAlpha), t0.x, t0.y, t0.z, B);
   float dLp = 0;
-  if (propagateToXBxByBz(t0, trackX, B[0], B[1], B[2], dLp)) {
+  if (propagateToXBzLightNoUpdate(t0, trackX, mBz, dLp)) {
     return -1;
   }
+  updatePhysicalTrackValues(t0);
 
   if (CAMath::Abs(t0.sinphi) >= MaxSinPhi) {
     return -1;
@@ -531,7 +505,7 @@ GPUd() int GPUTPCCompressionTrackModel::rotateToAlpha(float newAlpha)
   // only covariance changes. Use rotated and transported t0 for linearisation
   float j3 = -t0.py / t0.px;
   float j4 = -t0.pz / t0.px;
-  float j5 = t0.qpt * B[2];
+  float j5 = t0.qpt * mBz;
 
   //                    Y  Z Sin DzDs q/p  X
   // Jacobian J1 = { {  1, 0, 0,  0,  0,  j3 }, // Y
@@ -566,98 +540,6 @@ GPUd() int GPUTPCCompressionTrackModel::rotateToAlpha(float newAlpha)
   return 0;
 }
 
-GPUd() int GPUTPCCompressionTrackModel::propagateToXBxByBz(PhysicalTrackModel& t, float x, float Bx, float By, float Bz, float& dLp)
-{
-  //
-  // transport the track to X=x in magnetic field B = ( Bx, By, Bz )[kG*0.000299792458]
-  // xyzPxPyPz as well as all the additional values will change. No need to call UpdateValues() afterwards.
-  // the method returns error code (0 == no error)
-  //
-  dLp = 0.f;
-
-  // Rotate to the system where Bx=By=0.
-
-  float bt = CAMath::Sqrt(Bz * Bz + By * By);
-  float bb = CAMath::Sqrt(Bx * Bx + By * By + Bz * Bz);
-
-  float c1 = 1.f, s1 = 0.f;
-  float c2 = 1.f, s2 = 0.f;
-
-  if (bt > 1.e-4f) {
-    c1 = Bz / bt;
-    s1 = By / bt;
-    c2 = bt / bb;
-    s2 = -Bx / bb;
-  }
-
-  // rotation matrix: first around x, then around y'
-  // after the first rotation: Bx'==Bx, By'==0, Bz'==Bt, X'==X
-  // after the second rotation: Bx''==0, By''==0, Bz''==B, X'' axis is as close as possible to the original X
-
-  //
-  //     ( c2 0 s2 )   ( 1  0   0 )
-  // R = (  0 1 0  ) X ( 0 c1 -s1 )
-  //     (-s2 0 c2 )   ( 0 s1  c1 )
-  //
-
-  float R0[3] = {c2, s1 * s2, c1 * s2};
-  float R1[3] = {0, c1, -s1};
-  float R2[3] = {-s2, s1 * c2, c1 * c2};
-
-  // parameters and the extrapolation point in the rotated coordinate system
-  {
-    float lx = t.x, ly = t.y, lz = t.z, lpx = t.px, lpy = t.py, lpz = t.pz;
-
-    t.x = R0[0] * lx + R0[1] * ly + R0[2] * lz;
-    t.y = R1[0] * lx + R1[1] * ly + R1[2] * lz;
-    t.z = R2[0] * lx + R2[1] * ly + R2[2] * lz;
-
-    t.px = R0[0] * lpx + R0[1] * lpy + R0[2] * lpz;
-    t.py = R1[0] * lpx + R1[1] * lpy + R1[2] * lpz;
-    t.pz = R2[0] * lpx + R2[1] * lpy + R2[2] * lpz;
-  }
-
-  float dx = x - mX;
-  float xe = t.x + dx; // propagate on same dx in rotated system
-
-  // transport in rotated coordinate system to X''=xe:
-
-  if (t.px < (1.f - MaxSinPhi)) {
-    t.px = 1.f - MaxSinPhi;
-  }
-  if (propagateToXBzLightNoUpdate(t, xe, bb, dLp) != 0) {
-    return -1;
-  }
-
-  // rotate coordinate system back to the original R{-1}==R{T}
-  {
-    float lx = t.x, ly = t.y, lz = t.z, lpx = t.px, lpy = t.py, lpz = t.pz;
-
-    t.x = R0[0] * lx + R1[0] * ly + R2[0] * lz;
-    t.y = R0[1] * lx + R1[1] * ly + R2[1] * lz;
-    t.z = R0[2] * lx + R1[2] * ly + R2[2] * lz;
-
-    t.px = R0[0] * lpx + R1[0] * lpy + R2[0] * lpz;
-    t.py = R0[1] * lpx + R1[1] * lpy + R2[1] * lpz;
-    t.pz = R0[2] * lpx + R1[2] * lpy + R2[2] * lpz;
-  }
-
-  // a small (hopefully) additional step to X=x. Perhaps it may be replaced by linear extrapolation.
-
-  float ddLp = 0;
-  if (t.px < (1.f - MaxSinPhi)) {
-    t.px = 1.f - MaxSinPhi;
-  }
-  if (propagateToXBzLightNoUpdate(t, x, Bz, ddLp) != 0) {
-    return -1;
-  }
-
-  dLp += ddLp;
-
-  updatePhysicalTrackValues(t);
-  return 0;
-}
-
 GPUd() int GPUTPCCompressionTrackModel::propagateToXBzLightNoUpdate(PhysicalTrackModel& t, float x, float Bz, float& dLp)
 {
   //
diff --git a/GPU/GPUTracking/DataCompression/GPUTPCCompressionTrackModel.h b/GPU/GPUTracking/DataCompression/GPUTPCCompressionTrackModel.h
@@ -28,7 +28,6 @@
 #include "GPUTPCTrackParam.h"
 
 #else // Default internal track model for compression
-#include "GPUTPCGMPolynomialField.h"
 #endif
 
 namespace GPUCA_NAMESPACE::gpu
@@ -89,12 +88,9 @@ class GPUTPCCompressionTrackModel
   GPUd() float Z() const { return mP[1]; }
 
   // helper functions for standalone propagation and update methods
-  GPUd() void getBxByBz(float cosAlpha, float sinAlpha, float x, float y, float z, float b[3]) const;
-  GPUd() float getBz(float x, float y, float z) const;
   GPUd() void updatePhysicalTrackValues(PhysicalTrackModel& trk);
   GPUd() void changeDirection();
   GPUd() int rotateToAlpha(float newAlpha);
-  GPUd() int propagateToXBxByBz(PhysicalTrackModel& t, float x, float Bx, float By, float Bz, float& dLp);
   GPUd() int propagateToXBzLightNoUpdate(PhysicalTrackModel& t, float x, float Bz, float& dLp);
   GPUd() bool setDirectionAlongX(PhysicalTrackModel& t);
   GPUd() int followLinearization(const PhysicalTrackModel& t0e, float Bz, float dLp);
@@ -150,7 +146,7 @@ class GPUTPCCompressionTrackModel
   float mSinAlpha;
 
   // propagation parameters
-  const GPUTPCGMPolynomialField* mField = nullptr;
+  float mBz;
   MaterialCorrection mMaterial;
 
   PhysicalTrackModel mTrk;
