[ALICE3] Update FastTracker

ALICE3/Core/FastTracker.cxx

@@ -9,13 +9,17 @@
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
 
-#include <vector>
-#include <string>
+#include "FastTracker.h"
+
+#include "ReconstructionDataFormats/TrackParametrization.h"
+
 #include "TMath.h"
 #include "TMatrixD.h"
-#include "TRandom.h"
 #include "TMatrixDSymEigen.h"
-#include "FastTracker.h"
+#include "TRandom.h"
+
+#include <string>
+#include <vector>
 
 namespace o2
 {
@@ -24,43 +28,24 @@
 
 // +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+
 
-FastTracker::FastTracker()
-{
-  // base constructor
-  magneticField = 20; // in kiloGauss
-  applyZacceptance = false;
-  applyMSCorrection = true;
-  applyElossCorrection = true;
-  applyEffCorrection = true;
-  covMatFactor = 0.99f;
-  verboseLevel = 0;
-
-  // last fast-tracked track properties
-  covMatOK = 0;
-  covMatNotOK = 0;
-  nIntercepts = 0;
-  nSiliconPoints = 0;
-  nGasPoints = 0;
-
-  maxRadiusSlowDet = 10;
-  integrationTime = 0.02; // ms
-  crossSectionMinB = 8;
-  dNdEtaCent = 2200;
-  dNdEtaMinB = 1;
-  avgRapidity = 0.45;
-  sigmaD = 6.0;
-  luminosity = 1.e27;
-  otherBackground = 0.0; // [0, 1]
-  upcBackgroundMultiplier = 1.0;
-}
-
 void FastTracker::AddLayer(TString name, float r, float z, float x0, float xrho, float resRPhi, float resZ, float eff, int type)
 {
   DetLayer newLayer(name, r, z, x0, xrho, resRPhi, resZ, eff, type);
   // Check that efficient layers are not inert layers
   if (newLayer.getEfficiency() > 0.0f && newLayer.isInert()) {
     LOG(error) << "Layer " << name << " with efficiency > 0.0 should not be inert";
   }
+  // Layers should be ordered by increasing radius, check this
+  if (!layers.empty() && newLayer.getRadius() < layers.back().getRadius()) {
+    LOG(fatal) << "Layer " << newLayer << " is not ordered correctly, it should be after layer " << layers.back();
+  }
+  // Layers should all have different names
+  for (const auto& layer : layers) {
+    if (layer.getName() == newLayer.getName()) {
+      LOG(fatal) << "Layer with name " << newLayer.getName() << " already exists in FastTracker layers";
+    }
+  }
+  // Add the new layer to the layers vector
   layers.push_back(newLayer);
 }
 
@@ -78,7 +63,7 @@
   return layers[layerIdx];
 }
 
-int FastTracker::GetLayerIndex(std::string name) const
+int FastTracker::GetLayerIndex(const std::string& name) const
 {
   int i = 0;
   for (const auto& layer : layers) {
@@ -191,8 +176,8 @@
   for (int i = 0; i < kNPassiveBound; i++) {
     AddLayer(Form("tpc_boundary%d", i), rBoundary[i], zLength, radLBoundary[i], xrhoBoundary[i], 0); // dummy errors
   }
   for (Int_t k = 0; k < tpcRows; k++) {
     Float_t rowRadius = 0;
     if (k < innerRows)
       rowRadius = rowOneRadius + k * tpcInnerRadialPitch;
     else if (k >= innerRows && k < (innerRows + middleRows))
@@ -211,15 +196,15 @@
   // https://github.com/AliceO2Group/DelphesO2/blob/master/src/DetectorK/DetectorK.cxx#L743
   int index = 1;
   int nSteps = 301;
-  double dist = 0.0;
-  double dz0 = (4 * sigmaD - (-4) * sigmaD / (nSteps = 1));
-  double z0 = 0.0;
+  float dist = 0.0;
+  float dz0 = (4 * sigmaD - (-4) * sigmaD / (nSteps = 1));
+  float z0 = 0.0;
   for (int i = 0; i < nSteps; i++) {
     if (i == nSteps - 1)
       index = 1;
     z0 = -4 * sigmaD + i * dz0;
     dist += index * (dz0 / 3.) * (1 / o2::math_utils::sqrt(o2::constants::math::TwoPI) / sigmaD) * std::exp(-z0 * z0 / 2. / sigmaD / sigmaD) * (1 / o2::math_utils::sqrt((z - z0) * (z - z0) + r * r));
     if (index != 4)
       index = 4;
     else
       index = 2;
@@ -243,7 +228,7 @@
   // porting of DetektorK::IntegratedHitDensity
   // see here:
   // https://github.com/AliceO2Group/DelphesO2/blob/master/src/DetectorK/DetectorK.cxx#L712
-  float zdcHz = luminosity * 1.e24 * crossSectionMinB;
+  float zdcHz = luminosity * 1.e24 * mCrossSectionMinB;
   float den = zdcHz * integrationTime / 1000. * multiplicity * Dist(0., radius) / (o2::constants::math::TwoPI * radius);
   if (den < OneEventHitDensity(multiplicity, radius))
     den = OneEventHitDensity(multiplicity, radius);
@@ -301,6 +286,7 @@
 // returns number of intercepts (generic for now)
 int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackParCov& outputTrack, const float nch)
 {
+  dNdEtaCent = nch; // set the number of charged particles per unit rapidity
   hits.clear();
   nIntercepts = 0;
   nSiliconPoints = 0;
@@ -310,38 +296,54 @@
   const float initialRadius = std::hypot(posIni[0], posIni[1]);
   const float kTrackingMargin = 0.1;
   const int kMaxNumberOfDetectors = 20;
+  if (kMaxNumberOfDetectors < layers.size()) {
+    LOG(fatal) << "Too many layers in FastTracker, increase kMaxNumberOfDetectors";
+    return -1; // too many layers
+  }
+  int firstActiveLayer = -1; // first layer that is not inert
+  for (size_t i = 0; i < layers.size(); ++i) {
+    if (!layers[i].isInert()) {
+      firstActiveLayer = i;
+      break;
+    }
+  }
+  if (firstActiveLayer <= 0) {
+    LOG(fatal) << "No active layers found in FastTracker, check layer setup";
+    return -2; // no active layers
+  }
   const int xrhosteps = 100;
   const bool applyAngularCorrection = true;
 
   goodHitProbability.clear();
-  for (int i = 0; i < kMaxNumberOfDetectors; ++i)
+  for (int i = 0; i < kMaxNumberOfDetectors; ++i) {
     goodHitProbability.push_back(-1.);
-  goodHitProbability[0] = 1.;
+  }
+  goodHitProbability[0] = 1.; // we use layer zero to accumulate
 
   // +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+
   // Outward pass to find intercepts
   int firstLayerReached = -1;
   int lastLayerReached = -1;
   new (&outputTrack)(o2::track::TrackParCov)(inputTrack);
-  for (uint32_t il = 0; il < layers.size(); il++) {
+  for (size_t il = 0; il < layers.size(); il++) {
     // check if layer is doable
-    if (layers[il].getRadius() < initialRadius)
+    if (layers[il].getRadius() < initialRadius) {
       continue; // this layer should not be attempted, but go ahead
+    }
 
     // check if layer is reached
     float targetX = 1e+3;
-    bool ok = true;
     inputTrack.getXatLabR(layers[il].getRadius(), targetX, magneticField);
     if (targetX > 999.f) {
       LOGF(debug, "Failed to find intercept for layer %d at radius %.2f cm", il, layers[il].getRadius());
       break; // failed to find intercept
     }
 
-    ok = inputTrack.propagateTo(targetX, magneticField);
-    if (ok && applyMSCorrection && layers[il].getRadiationLength() > 0) {
+    bool ok = inputTrack.propagateTo(targetX, magneticField);
+    if (ok && mApplyMSCorrection && layers[il].getRadiationLength() > 0) {
       ok = inputTrack.correctForMaterial(layers[il].getRadiationLength(), 0, applyAngularCorrection);
     }
-    if (ok && applyElossCorrection && layers[il].getDensity() > 0) { // correct in small steps
+    if (ok && mApplyElossCorrection && layers[il].getDensity() > 0) { // correct in small steps
       for (int ise = xrhosteps; ise--;) {
         ok = inputTrack.correctForMaterial(0, -layers[il].getDensity() / xrhosteps, applyAngularCorrection);
         if (!ok)
@@ -361,7 +363,7 @@
         break;
       }
     }
-    if (std::abs(inputTrack.getZ()) > layers[il].getZ() && applyZacceptance) {
+    if (std::abs(inputTrack.getZ()) > layers[il].getZ() && mApplyZacceptance) {
       break; // out of acceptance bounds
     }
 
@@ -384,39 +386,19 @@
   o2::track::TrackParCov inwardTrack(inputTrack);
 
   // Enlarge covariance matrix
-  std::array<float, 5> trPars = {0.};
-  for (int ip = 0; ip < 5; ip++) {
+  std::array<float, o2::track::kNParams> trPars = {0.};
+  for (int ip = 0; ip < o2::track::kNParams; ip++) {
     trPars[ip] = outputTrack.getParam(ip);
   }
-  std::array<float, 15> largeCov = {0.};
-  enum { kY,
-         kZ,
-         kSnp,
-         kTgl,
-         kPtI }; // track parameter aliases
-  enum { kY2,
-         kYZ,
-         kZ2,
-         kYSnp,
-         kZSnp,
-         kSnp2,
-         kYTgl,
-         kZTgl,
-         kSnpTgl,
-         kTgl2,
-         kYPtI,
-         kZPtI,
-         kSnpPtI,
-         kTglPtI,
-         kPtI2 }; // cov.matrix aliases
-  const double kLargeErr2Coord = 5 * 5;
-  const double kLargeErr2Dir = 0.7 * 0.7;
-  const double kLargeErr2PtI = 30.5 * 30.5;
-  for (int ic = 15; ic--;)
+  static constexpr float kLargeErr2Coord = 5 * 5;
+  static constexpr float kLargeErr2Dir = 0.7 * 0.7;
+  static constexpr float kLargeErr2PtI = 30.5 * 30.5;
+  std::array<float, o2::track::kCovMatSize> largeCov = {0.};
+  for (int ic = o2::track::kCovMatSize; ic--;)
     largeCov[ic] = 0.;
-  largeCov[kY2] = largeCov[kZ2] = kLargeErr2Coord;
-  largeCov[kSnp2] = largeCov[kTgl2] = kLargeErr2Dir;
-  largeCov[kPtI2] = kLargeErr2PtI * trPars[kPtI] * trPars[kPtI];
+  largeCov[o2::track::CovLabels::kSigY2] = largeCov[o2::track::CovLabels::kSigZ2] = kLargeErr2Coord;
+  largeCov[o2::track::CovLabels::kSigSnp2] = largeCov[o2::track::CovLabels::kSigTgl2] = kLargeErr2Dir;
+  largeCov[o2::track::CovLabels::kSigQ2Pt2] = kLargeErr2PtI * trPars[o2::track::ParLabels::kQ2Pt] * trPars[o2::track::ParLabels::kQ2Pt];
 
   inwardTrack.setCov(largeCov);
   inwardTrack.checkCovariance();
@@ -427,14 +409,14 @@
 
     float targetX = 1e+3;
     inputTrack.getXatLabR(layers[il].getRadius(), targetX, magneticField);
     if (targetX > 999)
       continue; // failed to find intercept
 
     if (!inputTrack.propagateTo(targetX, magneticField)) {
       continue; // failed to propagate
     }
 
-    if (std::abs(inputTrack.getZ()) > layers[il].getZ() && applyZacceptance) {
+    if (std::abs(inputTrack.getZ()) > layers[il].getZ() && mApplyZacceptance) {
       continue; // out of acceptance bounds but continue inwards
     }
 
@@ -444,10 +426,10 @@
     std::vector<float> thisHit = {spacePoint[0], spacePoint[1], spacePoint[2]};
 
     // towards adding cluster: move to track alpha
-    double alpha = inwardTrack.getAlpha();
-    double xyz1[3]{
-      TMath::Cos(alpha) * spacePoint[0] + TMath::Sin(alpha) * spacePoint[1],
-      -TMath::Sin(alpha) * spacePoint[0] + TMath::Cos(alpha) * spacePoint[1],
+    float alpha = inwardTrack.getAlpha();
+    float xyz1[3]{
+      std::cos(alpha) * spacePoint[0] + std::sin(alpha) * spacePoint[1],
+      -std::sin(alpha) * spacePoint[0] + std::cos(alpha) * spacePoint[1],
       spacePoint[2]};
     if (!inwardTrack.propagateTo(xyz1[0], magneticField))
       continue;
@@ -462,15 +444,15 @@
       inwardTrack.checkCovariance();
     }
 
-    if (applyMSCorrection && layers[il].getRadiationLength() > 0) {
+    if (mApplyMSCorrection && layers[il].getRadiationLength() > 0) {
       if (!inputTrack.correctForMaterial(layers[il].getRadiationLength(), 0, applyAngularCorrection)) {
         return -6;
       }
       if (!inwardTrack.correctForMaterial(layers[il].getRadiationLength(), 0, applyAngularCorrection)) {
         return -6;
       }
     }
-    if (applyElossCorrection && layers[il].getDensity() > 0) {
+    if (mApplyElossCorrection && layers[il].getDensity() > 0) {
       for (int ise = xrhosteps; ise--;) { // correct in small steps
         if (!inputTrack.correctForMaterial(0, layers[il].getDensity() / xrhosteps, applyAngularCorrection)) {
           return -7;
@@ -488,40 +470,40 @@
 
     hits.push_back(thisHit);
 
-    if (applyEffCorrection && !layers[il].isInert()) { // good hit probability calculation
-      double sigYCmb = o2::math_utils::sqrt(inwardTrack.getSigmaY2() + layers[il].getResolutionRPhi() * layers[il].getResolutionRPhi());
-      double sigZCmb = o2::math_utils::sqrt(inwardTrack.getSigmaZ2() + layers[il].getResolutionZ() * layers[il].getResolutionZ());
+    if (!layers[il].isInert()) { // good hit probability calculation
+      float sigYCmb = o2::math_utils::sqrt(inwardTrack.getSigmaY2() + layers[il].getResolutionRPhi() * layers[il].getResolutionRPhi());
+      float sigZCmb = o2::math_utils::sqrt(inwardTrack.getSigmaZ2() + layers[il].getResolutionZ() * layers[il].getResolutionZ());
       goodHitProbability[il] = ProbGoodChiSqHit(layers[il].getRadius() * 100, sigYCmb * 100, sigZCmb * 100);
       goodHitProbability[0] *= goodHitProbability[il];
     }
   }
 
   // backpropagate to original radius
   float finalX = 1e+3;
-  inwardTrack.getXatLabR(initialRadius, finalX, magneticField);
-  if (finalX > 999)
+  bool inPropStatus = inwardTrack.getXatLabR(initialRadius, finalX, magneticField);
+  if (finalX > 999) {
+    LOG(debug) << "Failed to find intercept for initial radius " << initialRadius << " cm, x = " << finalX << " and status " << inPropStatus << " and sn = " << inwardTrack.getSnp() << " r = " << inwardTrack.getY() * inwardTrack.getY();
     return -3; // failed to find intercept
+  }
 
   if (!inwardTrack.propagateTo(finalX, magneticField)) {
     return -4; // failed to propagate
   }
 
   // only attempt to continue if intercepts are at least four
   if (nIntercepts < 4)
     return nIntercepts;
 
   // generate efficiency
-  if (applyEffCorrection) {
-    dNdEtaCent = nch;
-    float eff = 1.;
-    for (int i = 0; i < kMaxNumberOfDetectors; i++) {
-      float iGoodHit = goodHitProbability[i];
-      if (iGoodHit <= 0)
-        continue;
-
-      eff *= iGoodHit;
-    }
+  float eff = 1.;
+  for (int i = 0; i < kMaxNumberOfDetectors; i++) {
+    float iGoodHit = goodHitProbability[i];
+    if (iGoodHit <= 0)
+      continue;
 
+    eff *= iGoodHit;
+  }
+  if (mApplyEffCorrection) {
     if (gRandom->Uniform() > eff)
       return -8;
   }
@@ -530,11 +512,11 @@
   outputTrack.checkCovariance();
 
   // Use covariance matrix based smearing
-  std::array<double, 15> covMat = {0.};
-  for (int ii = 0; ii < 15; ii++)
+  std::array<float, o2::track::kCovMatSize> covMat = {0.};
+  for (int ii = 0; ii < o2::track::kCovMatSize; ii++)
     covMat[ii] = outputTrack.getCov()[ii];
   TMatrixDSym m(5);
-  double fcovm[5][5];
+  float fcovm[5][5];
 
   for (int ii = 0, k = 0; ii < 5; ++ii) {
     for (int j = 0; j < ii + 1; ++j, ++k) {
@@ -544,7 +526,7 @@
   }
 
   // evaluate ruben's conditional, regularise
-  bool makePositiveDefinite = (covMatFactor > -1e-5); // apply fix
+  const bool makePositiveDefinite = (covMatFactor > -1e-5); // apply fix
   bool rubenConditional = false;
   for (int ii = 0; ii < 5; ii++) {
     for (int jj = 0; jj < 5; jj++) {
@@ -553,7 +535,7 @@
       if (fcovm[ii][jj] * fcovm[ii][jj] > std::abs(fcovm[ii][ii] * fcovm[jj][jj])) {
         rubenConditional = true;
         if (makePositiveDefinite) {
           fcovm[ii][jj] = TMath::Sign(1, fcovm[ii][jj]) * covMatFactor * sqrt(std::abs(fcovm[ii][ii] * fcovm[jj][jj]));
         }
       }
     }
@@ -571,7 +553,7 @@
   }
 
   if (negEigVal && rubenConditional && makePositiveDefinite) {
-    if (verboseLevel > 0) {
+    if (mVerboseLevel > 0) {
       LOG(info) << "WARNING: this diagonalization (at pt = " << inputTrack.getPt() << ") has negative eigenvalues despite Ruben's fix! Please be careful!";
       LOG(info) << "Printing info:";
       LOG(info) << "Kalman updates: " << nIntercepts;
@@ -585,26 +567,26 @@
   covMatOK++;
 
   // transform parameter vector and smear
-  double params_[5];
+  float params_[5];
   for (int ii = 0; ii < 5; ++ii) {
-    double val = 0.;
+    float val = 0.;
     for (int j = 0; j < 5; ++j)
       val += eigVec[j][ii] * outputTrack.getParam(j);
     // smear parameters according to eigenvalues
     params_[ii] = gRandom->Gaus(val, sqrt(eigVal[ii]));
   }
 
   // invert eigenvector matrix
   eigVec.Invert();
   // transform back params vector
   for (int ii = 0; ii < 5; ++ii) {
-    double val = 0.;
+    float val = 0.;
     for (int j = 0; j < 5; ++j)
       val += eigVec[j][ii] * params_[j];
     outputTrack.setParam(val, ii);
   }
   // should make a sanity check that par[2] sin(phi) is in [-1, 1]
   if (fabs(outputTrack.getParam(2)) > 1.) {
     LOG(info) << " --- smearTrack failed sin(phi) sanity check: " << outputTrack.getParam(2);
     return -2;
   }