packages/trackreco/GPUTPCTrackParam.cxx

0001 // Copyright CERN and copyright holders of ALICE O2. This software is
0002 // distributed under the terms of the GNU General Public License v3 (GPL
0003 // Version 3), copied verbatim in the file "COPYING".
0004 //
0005 // See http://alice-o2.web.cern.ch/license for full licensing information.
0006 //
0007 // In applying this license CERN does not waive the privileges and immunities
0008 // granted to it by virtue of its status as an Intergovernmental Organization
0009 // or submit itself to any jurisdiction.
0010
0011 /// \file GPUTPCTrackParam.cxx
0012 /// \author Sergey Gorbunov, Ivan Kisel, David Rohr
0013
0014 #include "GPUTPCTrackLinearisation.h"
0015 #include "GPUTPCTrackParam.h"
0016 #include <cmath>
0017 #include <algorithm>
0018
0019 //
0020 // Circle in XY:
0021 //
0022 // kCLight = 0.000299792458;
0023 // Kappa = -Bz*kCLight*QPt;
0024 // R  = 1/fstd::absf(Kappa);
0025 // Xc = X - sin(Phi)/Kappa;
0026 // Yc = Y + cos(Phi)/Kappa;
0027 //
0028
0029 double GPUTPCTrackParam::GetDist2(const GPUTPCTrackParam& t) const
0030 {
0031   // get squared distance between tracks
0032
0033   double dx = GetX() - t.GetX();
0034   double dy = GetY() - t.GetY();
0035   double dz = GetZ() - t.GetZ();
0036   return dx * dx + dy * dy + dz * dz;
0037 }
0038
0039 double GPUTPCTrackParam::GetDistXZ2(const GPUTPCTrackParam& t) const
0040 {
0041   // get squared distance between tracks in X&Z
0042
0043   double dx = GetX() - t.GetX();
0044   double dz = GetZ() - t.GetZ();
0045   return dx * dx + dz * dz;
0046 }
0047
0048 double GPUTPCTrackParam::GetS(double x, double y, double Bz) const
0049 {
0050   //* Get XY path length to the given point
0051
0052   double k = GetKappa(Bz);
0053   double ex = GetCosPhi();
0054   double ey = GetSinPhi();
0055   x -= GetX();
0056   y -= GetY();
0057   double dS = x * ex + y * ey;
0058   if (std::abs(k) > 1.e-4f) {
0059     dS = atan2(k * dS, 1 + k * (x * ey - y * ex)) / k;
0060   }
0061   return dS;
0062 }
0063
0064 void GPUTPCTrackParam::GetDCAPoint(double x, double y, double z, double& xp, double& yp, double& zp, double Bz) const
0065 {
0066   //* Get the track point closest to the (x,y,z)
0067
0068   double x0 = GetX();
0069   double y0 = GetY();
0070   double k = GetKappa(Bz);
0071   double ex = GetCosPhi();
0072   double ey = GetSinPhi();
0073   double dx = x - x0;
0074   double dy = y - y0;
0075   double ax = dx * k + ey;
0076   double ay = dy * k - ex;
0077   double a = std::sqrt(ax * ax + ay * ay);
0078   xp = x0 + (dx - ey * ((dx * dx + dy * dy) * k - 2 * (-dx * ey + dy * ex)) / (a + 1)) / a;
0079   yp = y0 + (dy + ex * ((dx * dx + dy * dy) * k - 2 * (-dx * ey + dy * ex)) / (a + 1)) / a;
0080   double s = GetS(x, y, Bz);
0081   zp = GetZ() + GetDzDs() * s;
0082   if (std::abs(k) > 1.e-2f) {
0083     double dZ = std::abs(GetDzDs() * 2*M_PI / k);
0084     if (dZ > .1f) {
0085       zp += round((z - zp) / dZ) * dZ;
0086     }
0087   }
0088 }
0089
0090 //*
0091 //* Transport routines
0092 //*
0093
0094 bool GPUTPCTrackParam::TransportToX(double x, GPUTPCTrackLinearisation& t0, double Bz, double maxSinPhi, double* DL)
0095 {
0096   //* Transport the track parameters to X=x, using linearization at t0, and the field value Bz
0097   //* maxSinPhi is the max. allowed value for |t0.SinPhi()|
0098   //* linearisation of trajectory t0 is also transported to X=x,
0099   //* returns 1 if OK
0100   //*
0101
0102   double ex = t0.CosPhi();
0103   double ey = t0.SinPhi();
0104   double k = -t0.QPt() * Bz; // = 1/(radius of curvature)
0105
0106   double R = 1./fabs(k);
0107   double phi_center;
0108   if(k>0.) phi_center = atan2(ey,ex) - M_PI/2.;
0109   else phi_center = atan2(ey,ex) + M_PI/2.;
0110   double xc = X() - R*cos(phi_center);
0111   double yc = Y() - R*sin(phi_center);
0112
0113   double dx = x - X();
0114   double xnew = x - xc;
0115   double y0 = Y() - yc;
0116
0117
0118   double discriminant = R*R*y0*y0-xnew*xnew*y0*y0;
0119   if(discriminant<0.) return 0; // no intersection with new X
0120
0121   double ynew = sqrt(R*R-xnew*xnew);
0122   if(y0<0.) ynew *= -1.;
0123
0124   double dy = ynew - y0;
0125   double new_phi_center = atan2(ynew,xnew);
0126   double ex1, ey1;
0127   if(k>0.)
0128   {
0129     ey1 = sin(new_phi_center+M_PI/2.);
0130     ex1 = cos(new_phi_center+M_PI/2.);
0131   }
0132   else
0133   {
0134     ey1 = sin(new_phi_center-M_PI/2.);
0135     ex1 = cos(new_phi_center-M_PI/2.);
0136   }
0137
0138   if(fabs(ey1) > maxSinPhi) return 0;
0139   double dl = sqrt(dx*dx+dy*dy);
0140 /*
0141   double ey1 = k * dx + ey; // dX/R as rotation angle only works if dY is small
0142   double ex1;
0143
0144   // check for intersection with X=x
0145
0146   if (std::abs(ey1) > maxSinPhi) {
0147     return 0;
0148   }
0149
0150   ex1 = std::sqrt(1 - ey1 * ey1);
0151   if (ex < 0) {
0152     ex1 = -ex1;
0153   }
0154 */
0155   double dx2 = dx * dx;
0156   double ss = ey + ey1;
0157   double cc = ex + ex1;
0158 /*
0159   if (std::abs(cc) < 1.e-4f || std::abs(ex) < 1.e-4f || std::abs(ex1) < 1.e-4f) {
0160     return 0;
0161   }
0162
0163   double tg = ss / cc; // tanf((phi1+phi)/2)
0164
0165   double dy = dx * tg;
0166   double dl = dx * std::sqrt(1 + tg * tg);
0167
0168   if (cc < 0) {
0169     dl = -dl;
0170   }
0171 */
0172   double dSin = dl * k / 2;
0173   if (dSin > 1) {
0174     dSin = 1;
0175   }
0176   if (dSin < -1) {
0177     dSin = -1;
0178   }
0179   double dphi_center = new_phi_center-phi_center;
0180   if(dphi_center>M_PI) dphi_center = 2.*M_PI-dphi_center;
0181   if(dphi_center<-M_PI) dphi_center = 2.*M_PI+dphi_center;
0182   double dS = R*fabs(dphi_center);//(std::abs(k) > 1.e-4f) ? (2 * asin(dSin) / k) : dl;
0183   double dz = dS * t0.DzDs();
0184
0185   if (DL) {
0186     *DL = -dS * std::sqrt(1 + t0.DzDs() * t0.DzDs());
0187   }
0188
0189   double cci = 1.f / cc;
0190   double exi = 1.f / ex;
0191   double ex1i = 1.f / ex1;
0192
0193 //  double d[5] = {0, 0, GetPar(2) - t0.SinPhi(), GetPar(3) - t0.DzDs(), GetPar(4) - t0.QPt()};
0194
0195   // double H0[5] = { 1,0, h2,  0, h4 };
0196   // double H1[5] = { 0, 1, 0, dS,  0 };
0197   // double H2[5] = { 0, 0, 1,  0, dxBz };
0198   // double H3[5] = { 0, 0, 0,  1,  0 };
0199   // double H4[5] = { 0, 0, 0,  0,  1 };
0200
0201   double h2 = dx * (1 + ey * ey1 + ex * ex1) * exi * ex1i * cci;
0202   double h4 = dx2 * (cc + ss * ey1 * ex1i) * cci * cci * (-Bz);
0203   double dxBz = dx * (-Bz);
0204 /*
0205   SetX(X() + dx);
0206   SetPar(0, Y() + dy + h2 * d[2] + h4 * d[4]);
0207   SetPar(1, Z() + dz + dS * d[3]);
0208   SetPar(2, t0.SinPhi() + d[2] + dxBz * d[4]);
0209 */
0210
0211   SetX(X() + dx);
0212   SetPar(0, Y() + dy);
0213   SetPar(1, Z() + dz);
0214   SetPar(2, ey1);
0215
0216   t0.SetCosPhi(ex1);
0217   t0.SetSinPhi(ey1);
0218
0219   double c00 = mC[0];
0220   double c10 = mC[1];
0221   double c11 = mC[2];
0222   double c20 = mC[3];
0223   double c21 = mC[4];
0224   double c22 = mC[5];
0225   double c30 = mC[6];
0226   double c31 = mC[7];
0227   double c32 = mC[8];
0228   double c33 = mC[9];
0229   double c40 = mC[10];
0230   double c41 = mC[11];
0231   double c42 = mC[12];
0232   double c43 = mC[13];
0233   double c44 = mC[14];
0234
0235   mC[0] = c00 + h2 * h2 * c22 + h4 * h4 * c44 + 2 * (h2 * c20 + h4 * c40 + h2 * h4 * c42);
0236
0237   mC[1] = c10 + h2 * c21 + h4 * c41 + dS * (c30 + h2 * c32 + h4 * c43);
0238   mC[2] = c11 + 2 * dS * c31 + dS * dS * c33;
0239
0240   mC[3] = c20 + h2 * c22 + h4 * c42 + dxBz * (c40 + h2 * c42 + h4 * c44);
0241   mC[4] = c21 + dS * c32 + dxBz * (c41 + dS * c43);
0242   mC[5] = c22 + 2 * dxBz * c42 + dxBz * dxBz * c44;
0243
0244   mC[6] = c30 + h2 * c32 + h4 * c43;
0245   mC[7] = c31 + dS * c33;
0246   mC[8] = c32 + dxBz * c43;
0247   mC[9] = c33;
0248
0249   mC[10] = c40 + h2 * c42 + h4 * c44;
0250   mC[11] = c41 + dS * c43;
0251   mC[12] = c42 + dxBz * c44;
0252   mC[13] = c43;
0253   mC[14] = c44;
0254
0255   return 1;
0256 }
0257
0258 bool GPUTPCTrackParam::TransportToX(double x, double sinPhi0, double cosPhi0, double Bz, double maxSinPhi)
0259 {
0260   //* Transport the track parameters to X=x, using linearization at phi0 with 0 curvature,
0261   //* and the field value Bz
0262   //* maxSinPhi is the max. allowed value for |t0.SinPhi()|
0263   //* linearisation of trajectory t0 is also transported to X=x,
0264   //* returns 1 if OK
0265   //*
0266
0267   double ex = cosPhi0;
0268   double ey = sinPhi0;
0269   double dx = x - X();
0270
0271   if (std::abs(ex) < 1.e-4f) {
0272     return 0;
0273   }
0274   double exi = 1.f / ex;
0275
0276   double dxBz = dx * (-Bz);
0277   double dS = dx * exi;
0278   double h2 = dS * exi * exi;
0279   double h4 = .5f * h2 * dxBz;
0280
0281   // double H0[5] = { 1,0, h2,  0, h4 };
0282   // double H1[5] = { 0, 1, 0, dS,  0 };
0283   // double H2[5] = { 0, 0, 1,  0, dxBz };
0284   // double H3[5] = { 0, 0, 0,  1,  0 };
0285   // double H4[5] = { 0, 0, 0,  0,  1 };
0286
0287   double sinPhi = SinPhi() + dxBz * QPt();
0288   if (maxSinPhi > 0 && std::abs(sinPhi) > maxSinPhi) {
0289     return 0;
0290   }
0291
0292   SetX(X() + dx);
0293   SetPar(0, GetPar(0) + dS * ey + h2 * (SinPhi() - ey) + h4 * QPt());
0294   SetPar(1, GetPar(1) + dS * DzDs());
0295   SetPar(2, sinPhi);
0296
0297   double c00 = mC[0];
0298   double c10 = mC[1];
0299   double c11 = mC[2];
0300   double c20 = mC[3];
0301   double c21 = mC[4];
0302   double c22 = mC[5];
0303   double c30 = mC[6];
0304   double c31 = mC[7];
0305   double c32 = mC[8];
0306   double c33 = mC[9];
0307   double c40 = mC[10];
0308   double c41 = mC[11];
0309   double c42 = mC[12];
0310   double c43 = mC[13];
0311   double c44 = mC[14];
0312
0313   // This is not the correct propagation!!! The max ensures the positional error does not decrease during track finding.
0314   // A different propagation method is used for the fit.
0315   mC[0] = std::max(c00, c00 + h2 * h2 * c22 + h4 * h4 * c44 + 2 * (h2 * c20 + h4 * c40 + h2 * h4 * c42));
0316
0317   mC[1] = c10 + h2 * c21 + h4 * c41 + dS * (c30 + h2 * c32 + h4 * c43);
0318   mC[2] = std::max(c11, c11 + 2 * dS * c31 + dS * dS * c33);
0319
0320   mC[3] = c20 + h2 * c22 + h4 * c42 + dxBz * (c40 + h2 * c42 + h4 * c44);
0321   mC[4] = c21 + dS * c32 + dxBz * (c41 + dS * c43);
0322   mC[5] = c22 + 2 * dxBz * c42 + dxBz * dxBz * c44;
0323
0324   mC[6] = c30 + h2 * c32 + h4 * c43;
0325   mC[7] = c31 + dS * c33;
0326   mC[8] = c32 + dxBz * c43;
0327   mC[9] = c33;
0328
0329   mC[10] = c40 + h2 * c42 + h4 * c44;
0330   mC[11] = c41 + dS * c43;
0331   mC[12] = c42 + dxBz * c44;
0332   mC[13] = c43;
0333   mC[14] = c44;
0334
0335   return 1;
0336 }
0337
0338 bool GPUTPCTrackParam::TransportToX(double x, double Bz, double maxSinPhi)
0339 {
0340   //* Transport the track parameters to X=x
0341   GPUTPCTrackLinearisation t0(*this);
0342   return TransportToX(x, t0, Bz, maxSinPhi);
0343 }
0344
0345 bool GPUTPCTrackParam::TransportToXWithMaterial(double x, GPUTPCTrackLinearisation& t0, GPUTPCTrackFitParam& par, double Bz, double maxSinPhi)
0346 {
0347   //* Transport the track parameters to X=x  taking into account material budget
0348
0349 //  const double kRho = 1.025e-3f;  // 0.9e-3;
0350 //  const double kRadLen = 29.532f; // 28.94;
0351   double Ne_nEff = Ne_frac * Ne_Rho / Ne_mA;
0352   double Ar_nEff = Ar_frac * Ar_Rho / Ar_mA; //Scaled effective number of particles in mix (Avogadro's number and density scale cancel in the ratio with total particles)
0353   double CF4_nEff = CF4_frac * CF4_Rho / CF4_mA;
0354   double N2_nEff = N2_frac * N2_Rho / N2_mA;
0355   double isobutane_nEff = isobutane_frac * isobutane_Rho / isobutane_mA;
0356   double nEff = Ar_nEff + CF4_nEff + N2_nEff + isobutane_nEff;
0357
0358   // these are for the sPHENIX TPC gas mixture
0359   const double kRho = Ne_frac * Ne_Rho
0360                     + Ar_frac * Ar_Rho
0361                     + CF4_frac * CF4_Rho
0362                     + N2_frac * N2_Rho
0363                     + isobutane_frac * isobutane_Rho;
0364
0365   const double kRadLen = (1/nEff) * ((Ne_nEff * Ne_RadLen)
0366                                   +  (Ar_nEff * Ar_RadLen)
0367                                   +  (CF4_nEff * CF4_RadLen)
0368                                   +  (N2_nEff * N2_RadLen)
0369                                   +  (isobutane_nEff * isobutane_RadLen));
0370   const double kRhoOverRadLen = kRho / kRadLen;
0371   double dl;
0372
0373   if (!TransportToX(x, t0, Bz, maxSinPhi, &dl)) {
0374     return 0;
0375   }
0376
0377   CorrectForMeanMaterial(dl * kRhoOverRadLen, dl * kRho, par);
0378   return 1;
0379 }
0380
0381 bool GPUTPCTrackParam::TransportToXWithMaterial(double x, GPUTPCTrackFitParam& par, double Bz, double maxSinPhi)
0382 {
0383   //* Transport the track parameters to X=x  taking into account material budget
0384
0385   GPUTPCTrackLinearisation t0(*this);
0386   return TransportToXWithMaterial(x, t0, par, Bz, maxSinPhi);
0387 }
0388
0389 bool GPUTPCTrackParam::TransportToXWithMaterial(double x, double Bz, double maxSinPhi)
0390 {
0391   //* Transport the track parameters to X=x taking into account material budget
0392
0393   GPUTPCTrackFitParam par;
0394   CalculateFitParameters(par);
0395   return TransportToXWithMaterial(x, par, Bz, maxSinPhi);
0396 }
0397
0398 //*
0399 //*  Multiple scattering and energy losses
0400 //*
0401 double GPUTPCTrackParam::BetheBlochGeant(double bg2, double kp0, double kp1, double kp2, double kp3, double kp4)
0402 {
0403   //
0404   // This is the parameterization of the Bethe-Bloch formula inspired by Geant.
0405   //
0406   // bg2  - (beta*gamma)^2
0407   // kp0 - density [g/cm^3]
0408   // kp1 - density effect first junction point
0409   // kp2 - density effect second junction point
0410   // kp3 - mean excitation energy [GeV]
0411   // kp4 - mean Z/A
0412   //
0413   // The default values for the kp* parameters are for silicon.
0414   // The returned value is in [GeV/(g/cm^2)].
0415   //
0416
0417   const double mK = 0.307075e-3f; // [GeV*cm^2/g]
0418   const double me = 0.511e-3f;    // [GeV/c^2]
0419   const double rho = kp0;
0420   const double x0 = kp1 * 2.303f;
0421   const double x1 = kp2 * 2.303f;
0422   const double mI = kp3;
0423   const double mZA = kp4;
0424   const double maxT = 2 * me * bg2; // neglecting the electron mass
0425
0426   //*** Density effect
0427   double d2 = 0.f;
0428   const double x = 0.5f * log(bg2);
0429   const double lhwI = log(28.816f * 1e-9f * std::sqrt(rho * mZA) / mI);
0430   if (x > x1) {
0431     d2 = lhwI + x - 0.5f;
0432   } else if (x > x0) {
0433     const double r = (x1 - x) / (x1 - x0);
0434     d2 = lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r;
0435   }
0436
0437   return mK * mZA * (1 + bg2) / bg2 * (0.5f * log(2 * me * bg2 * maxT / (mI * mI)) - bg2 / (1 + bg2) - d2);
0438 }
0439
0440 double GPUTPCTrackParam::BetheBlochSolid(double bg)
0441 {
0442   //------------------------------------------------------------------
0443   // This is an approximation of the Bethe-Bloch formula,
0444   // reasonable for solid materials.
0445   // All the parameters are, in fact, for Si.
0446   // The returned value is in [GeV]
0447   //------------------------------------------------------------------
0448
0449   return BetheBlochGeant(bg);
0450 }
0451
0452 double GPUTPCTrackParam::BetheBlochGas(double bg)
0453 {
0454   //------------------------------------------------------------------
0455   // This is an approximation of the Bethe-Bloch formula,
0456   // reasonable for gas materials.
0457   // All the parameters are, in fact, for Ne.
0458   // The returned value is in [GeV]
0459   //------------------------------------------------------------------
0460
0461 //  const double rho = 0.9e-3f;
0462 //  const double x0 = 2.f;
0463 //  const double x1 = 4.f;
0464 //  const double mI = 140.e-9f;
0465 //  const double mZA = 0.49555f;
0466
0467   double Ne_nEff = Ne_frac * Ne_Rho / Ne_mA; //scaled effective number of particles in mix (avogadro's number and density scale cancel in the ratio with total particles)
0468   double Ar_nEff = Ar_frac * Ar_Rho / Ar_mA; //scaled effective number of particles in mix (avogadro's number and density scale cancel in the ratio with total particles)
0469   double CF4_nEff = CF4_frac * CF4_Rho / CF4_mA;
0470   double N2_nEff = N2_frac * N2_Rho / N2_mA;
0471   double isobutane_nEff = isobutane_frac * isobutane_Rho / isobutane_mA;
0472   double nEff = Ar_nEff + CF4_nEff + N2_nEff + isobutane_nEff;
0473
0474   // these are for the sPHENIX TPC gas mixture
0475 //  const double rho = 2.485e-3f;
0476 //  const double x0 = 2.f;
0477 //  const double x1 = 4.f;
0478 //  const double mI = 11.6e-9f;
0479 //  const double mZA = 0.46158f;
0480
0481   const double rho = Ne_frac * Ne_Rho
0482                    + Ar_frac * Ar_Rho
0483                    + CF4_frac * CF4_Rho
0484                    + N2_frac * N2_Rho
0485                    + isobutane_frac * isobutane_Rho;
0486
0487   const double x0 = (1/nEff) * ((Ne_nEff * Ne_x0)
0488                              +  (Ar_nEff * Ar_x0)
0489                              +  (CF4_nEff * CF4_x0)
0490                              +  (N2_nEff * N2_x0)
0491                              +  (isobutane_nEff * isobutane_x0));
0492
0493   const double x1 = (1/nEff) * ((Ne_nEff * Ne_x1)
0494                              +  (Ar_nEff * Ar_x1)
0495                              +  (CF4_nEff * CF4_x1)
0496                              +  (N2_nEff * N2_x1)
0497                              +  (isobutane_nEff * isobutane_x1));
0498
0499   const double mI = (1/nEff) * ((Ne_nEff * Ne_mI)
0500                              +  (Ar_nEff * Ar_mI)
0501                              +  (CF4_nEff * CF4_mI)
0502                              +  (N2_nEff * N2_mI)
0503                              +  (isobutane_nEff * isobutane_mI));
0504
0505   const double mZA = (1/nEff) * ((Ne_nEff * Ne_mZA)
0506                               +  (Ar_nEff * Ar_mZA)
0507                               +  (CF4_nEff * CF4_mZA)
0508                               +  (N2_nEff * N2_mZA)
0509                               +  (isobutane_nEff * isobutane_mZA));
0510   return BetheBlochGeant(bg, rho, x0, x1, mI, mZA);
0511 }
0512
0513 double GPUTPCTrackParam::ApproximateBetheBloch(double beta2)
0514 {
0515   //------------------------------------------------------------------
0516   // This is an approximation of the Bethe-Bloch formula with
0517   // the density effect taken into account at beta*gamma > 3.5
0518   // (the approximation is reasonable only for solid materials)
0519   //------------------------------------------------------------------
0520   if (beta2 >= 1) {
0521     return 0;
0522   }
0523
0524   if (beta2 / (1 - beta2) > 3.5f * 3.5f) {
0525     return 0.153e-3f / beta2 * (log(3.5f * 5940) + 0.5f * log(beta2 / (1 - beta2)) - beta2);
0526   }
0527   return 0.153e-3f / beta2 * (log(5940 * beta2 / (1 - beta2)) - beta2);
0528 }
0529
0530 void GPUTPCTrackParam::CalculateFitParameters(GPUTPCTrackFitParam& par, double mass)
0531 {
0532   //*!
0533
0534   double qpt = GetPar(4);
0535   //if (mC[14] >= 1.f) {
0536   //  qpt = 1.f / 0.35f;
0537   //}
0538
0539   double p2 = (1.f + GetPar(3) * GetPar(3));
0540   double k2 = qpt * qpt;
0541   double mass2 = mass * mass;
0542   double beta2 = p2 / (p2 + mass2 * k2);
0543
0544   double pp2 = (k2 > 1.e-8f) ? p2 / k2 : 10000; // impuls 2
0545
0546   par.bethe = BetheBlochGas( beta2/(1+beta2) );
0547   //par.bethe = ApproximateBetheBloch(pp2 / mass2);
0548   par.e = std::sqrt(pp2 + mass2);
0549   par.theta2 = 14.1f * 14.1f / (beta2 * pp2 * 1e6f);
0550   par.EP2 = par.e / pp2;
0551
0552   // Approximate energy loss fluctuation (M.Ivanov)
0553
0554   const double knst = 0.0007f; // To be tuned.
0555   par.sigmadE2 = knst * par.EP2 * qpt;
0556   par.sigmadE2 = par.sigmadE2 * par.sigmadE2;
0557
0558   par.k22 = (1.f + GetPar(3) * GetPar(3));
0559   par.k33 = par.k22 * par.k22;
0560   par.k43 = 0;
0561   par.k44 = GetPar(3) * GetPar(3) * k2;
0562 }
0563
0564 bool GPUTPCTrackParam::CorrectForMeanMaterial(double xOverX0, double xTimesRho, const GPUTPCTrackFitParam& par)
0565 {
0566   //------------------------------------------------------------------
0567   // This function corrects the track parameters for the crossed material.
0568   // "xOverX0"   - X/X0, the thickness in units of the radiation length.
0569   // "xTimesRho" - is the product length*density (g/cm^2).
0570   //------------------------------------------------------------------
0571
0572   double& mC22 = mC[5];
0573   double& mC33 = mC[9];
0574   double& mC40 = mC[10];
0575   double& mC41 = mC[11];
0576   double& mC42 = mC[12];
0577   double& mC43 = mC[13];
0578   double& mC44 = mC[14];
0579
0580   // Energy losses************************
0581
0582   double dE = par.bethe * xTimesRho;
0583   if (std::abs(dE) > 0.3f * par.e) {
0584     return 0; // 30% energy loss is too much!
0585   }
0586   double corr = (1.f - par.EP2 * dE);
0587   if (corr < 0.3f || corr > 1.3f) {
0588     return 0;
0589   }
0590
0591   SetPar(4, GetPar(4) * corr);
0592   mC40 *= corr;
0593   mC41 *= corr;
0594   mC42 *= corr;
0595   mC43 *= corr;
0596   mC44 *= corr * corr;
0597   mC44 += par.sigmadE2 * std::abs(dE);
0598
0599   // Multiple scattering******************
0600
0601   double theta2 = par.theta2 * std::abs(xOverX0);
0602   mC22 += theta2 * par.k22 * (1.f - GetPar(2)) * (1.f + GetPar(2));
0603   mC33 += theta2 * par.k33;
0604   mC43 += theta2 * par.k43;
0605   mC44 += theta2 * par.k44;
0606
0607   return 1;
0608 }
0609
0610 //*
0611 //* Rotation
0612 //*
0613 bool GPUTPCTrackParam::Rotate(double alpha, double maxSinPhi)
0614 {
0615   //* Rotate the coordinate system in XY on the angle alpha
0616
0617   double cA = cos(alpha);
0618   double sA = sin(alpha);
0619   double x = X(), y = Y(), sP = SinPhi(), cP = GetCosPhi();
0620   double cosPhi = cP * cA + sP * sA;
0621   double sinPhi = -cP * sA + sP * cA;
0622
0623   if (std::abs(sinPhi) > maxSinPhi) {// || std::abs(cosPhi) < 1.e-2f || std::abs(cP) < 1.e-2f) {
0624     return 0;
0625   }
0626
0627   double j0 = cP / cosPhi;
0628   double j2 = cosPhi / cP;
0629
0630   SetX(x * cA + y * sA);
0631   SetY(-x * sA + y * cA);
0632   SetSignCosPhi(cosPhi);
0633   SetSinPhi(sinPhi);
0634
0635   // double J[5][5] = { { j0, 0, 0,  0,  0 }, // Y
0636   //                      {  0, 1, 0,  0,  0 }, // Z
0637   //                      {  0, 0, j2, 0,  0 }, // SinPhi
0638   //                    {  0, 0, 0,  1,  0 }, // DzDs
0639   //                    {  0, 0, 0,  0,  1 } }; // Kappa
0640   // cout<<"alpha="<<alpha<<" "<<x<<" "<<y<<" "<<sP<<" "<<cP<<" "<<j0<<" "<<j2<<endl;
0641   // cout<<"      "<<mC[0]<<" "<<mC[1]<<" "<<mC[6]<<" "<<mC[10]<<" "<<mC[4]<<" "<<mC[5]<<" "<<mC[8]<<" "<<mC[12]<<endl;
0642   mC[0] *= j0 * j0;
0643   mC[1] *= j0;
0644   mC[3] *= j0;
0645   mC[6] *= j0;
0646   mC[10] *= j0;
0647
0648   mC[3] *= j2;
0649   mC[4] *= j2;
0650   mC[5] *= j2 * j2;
0651   mC[8] *= j2;
0652   mC[12] *= j2;
0653 /*
0654 //  if (cosPhi < 0) {
0655   if((cP>0 && cosPhi<0) || (cP<0 && cosPhi>0))
0656   {
0657     SetSinPhi(-SinPhi());
0658     SetDzDs(-DzDs());
0659     SetQPt(-QPt());
0660     mC[3] = -mC[3];
0661     mC[4] = -mC[4];
0662     mC[6] = -mC[6];
0663     mC[7] = -mC[7];
0664     mC[10] = -mC[10];
0665     mC[11] = -mC[11];
0666   }
0667 */
0668   // cout<<"      "<<mC[0]<<" "<<mC[1]<<" "<<mC[6]<<" "<<mC[10]<<" "<<mC[4]<<" "<<mC[5]<<" "<<mC[8]<<" "<<mC[12]<<endl;
0669   return 1;
0670 }
0671
0672 bool GPUTPCTrackParam::Rotate(double alpha, GPUTPCTrackLinearisation& t0, double maxSinPhi)
0673 {
0674   //* Rotate the coordinate system in XY on the angle alpha
0675
0676   double cA = cos(alpha);
0677   double sA = sin(alpha);
0678   double x0 = X(), y0 = Y(), sP = t0.SinPhi(), cP = t0.CosPhi();
0679   double cosPhi = cP * cA + sP * sA;
0680   double sinPhi = -cP * sA + sP * cA;
0681
0682   if (std::abs(sinPhi) > maxSinPhi || std::abs(cosPhi) < 1.e-2f || std::abs(cP) < 1.e-2f) {
0683     return 0;
0684   }
0685
0686   // double J[5][5] = { { j0, 0, 0,  0,  0 }, // Y
0687   //                    {  0, 1, 0,  0,  0 }, // Z
0688   //                    {  0, 0, j2, 0,  0 }, // SinPhi
0689   //                  {  0, 0, 0,  1,  0 }, // DzDs
0690   //                  {  0, 0, 0,  0,  1 } }; // Kappa
0691
0692   double j0 = cP / cosPhi;
0693   double j2 = cosPhi / cP;
0694   double d[2] = {Y() - y0, SinPhi() - sP};
0695
0696   SetX(x0 * cA + y0 * sA);
0697   SetY(-x0 * sA + y0 * cA + j0 * d[0]);
0698   t0.SetCosPhi(cosPhi);
0699   t0.SetSinPhi(sinPhi);
0700
0701   SetSinPhi(sinPhi + j2 * d[1]);
0702
0703   mC[0] *= j0 * j0;
0704   mC[1] *= j0;
0705   mC[3] *= j0;
0706   mC[6] *= j0;
0707   mC[10] *= j0;
0708
0709   mC[3] *= j2;
0710   mC[4] *= j2;
0711   mC[5] *= j2 * j2;
0712   mC[8] *= j2;
0713   mC[12] *= j2;
0714
0715   return 1;
0716 }
0717
0718 bool GPUTPCTrackParam::Filter(double y, double z, double err2Y, double err2Z, double maxSinPhi, bool paramOnly)
0719 {
0720   //* Add the y,z measurement with the Kalman filter
0721
0722   double c00 = mC[0], c11 = mC[2], c20 = mC[3], c31 = mC[7], c40 = mC[10];
0723
0724   err2Y += c00;
0725   err2Z += c11;
0726
0727   double z0 = y - GetPar(0), z1 = z - GetPar(1);
0728
0729   if (err2Y < 1.e-8f || err2Z < 1.e-8f) {
0730     return 0;
0731   }
0732
0733   double mS0 = 1.f / err2Y;
0734   double mS2 = 1.f / err2Z;
0735
0736   // K = CHtS
0737
0738   double k00, k11, k20, k31, k40;
0739
0740   k00 = c00 * mS0;
0741   k20 = c20 * mS0;
0742   k40 = c40 * mS0;
0743
0744   k11 = c11 * mS2;
0745   k31 = c31 * mS2;
0746
0747   double sinPhi = GetPar(2) + k20 * z0;
0748
0749   if (maxSinPhi > 0 && std::abs(sinPhi) >= maxSinPhi) {
0750     return 0;
0751   }
0752
0753   SetPar(0, GetPar(0) + k00 * z0);
0754   SetPar(1, GetPar(1) + k11 * z1);
0755   SetPar(2, sinPhi);
0756   SetPar(3, GetPar(3) + k31 * z1);
0757   SetPar(4, GetPar(4) + k40 * z0);
0758   if (paramOnly) {
0759     return true;
0760   }
0761
0762   mNDF += 2;
0763   mChi2 += mS0 * z0 * z0 + mS2 * z1 * z1;
0764
0765   mC[0] -= k00 * c00;
0766   mC[3] -= k20 * c00;
0767   mC[5] -= k20 * c20;
0768   mC[10] -= k40 * c00;
0769   mC[12] -= k40 * c20;
0770   mC[14] -= k40 * c40;
0771
0772   mC[2] -= k11 * c11;
0773   mC[7] -= k31 * c11;
0774   mC[9] -= k31 * c31;
0775
0776   return 1;
0777 }
0778
0779 bool GPUTPCTrackParam::CheckNumericalQuality() const
0780 {
0781   //* Check that the track parameters and covariance matrix are reasonable
0782
0783   bool ok = std::isfinite(GetX()) && std::isfinite(mSignCosPhi) && std::isfinite(mChi2);
0784
0785   const double* c = Cov();
0786   for (int i = 0; i < 15; i++) {
0787     ok = ok && std::isfinite(c[i]);
0788   }
0789   for (int i = 0; i < 5; i++) {
0790     ok = ok && std::isfinite(Par()[i]);
0791   }
0792
0793   if (c[0] <= 0 || c[2] <= 0 || c[5] <= 0 || c[9] <= 0 || c[14] <= 0) {
0794     ok = 0;
0795   }
0796   if (c[0] > 5.f || c[2] > 5.f || c[5] > 2.f || c[9] > 2
0797       //|| ( std::abs( QPt() ) > 1.e-2 && c[14] > 2. )
0798   ) {
0799     ok = 0;
0800   }
0801
0802   if (std::abs(SinPhi()) > GPUCA_MAX_SIN_PHI) {
0803     ok = 0;
0804   }
0805   if (std::abs(QPt()) > 1.f / 0.05f) {
0806     ok = 0;
0807   }
0808   if (ok) {
0809     ok = ok && (c[1] * c[1] <= c[2] * c[0]) && (c[3] * c[3] <= c[5] * c[0]) && (c[4] * c[4] <= c[5] * c[2]) && (c[6] * c[6] <= c[9] * c[0]) && (c[7] * c[7] <= c[9] * c[2]) && (c[8] * c[8] <= c[9] * c[5]) && (c[10] * c[10] <= c[14] * c[0]) && (c[11] * c[11] <= c[14] * c[2]) &&
0810          (c[12] * c[12] <= c[14] * c[5]) && (c[13] * c[13] <= c[14] * c[9]);
0811   }
0812   return ok;
0813 }
0814
0815 #if !defined(GPUCA_GPUCODE)
0816 #include <iostream>
0817 #endif
0818
0819 void GPUTPCTrackParam::Print() const
0820 {
0821   //* print parameters
0822
0823 #if !defined(GPUCA_GPUCODE)
0824   std::cout << "track: x=" << GetX() << " c=" << GetSignCosPhi() << ", P= " << GetY() << " " << GetZ() << " " << GetSinPhi() << " " << GetDzDs() << " " << GetQPt() << std::endl;
0825   std::cout << "errs2: " << GetErr2Y() << " " << GetErr2Z() << " " << GetErr2SinPhi() << " " << GetErr2DzDs() << " " << GetErr2QPt() << std::endl;
0826 #endif
0827 }