[ALICE3] Add ms & eloss to fasttracker (AliceO2Group#9292)

Author: jesgum

ALICE3/Core/FastTracker.cxx

@@ -26,7 +26,7 @@ namespace fastsim
 FastTracker::FastTracker()
 {
   // base constructor
-  magneticField = 5; // in kiloGauss
+  magneticField = 20; // in kiloGauss
   applyZacceptance = false;
   covMatFactor = 0.99f;
   verboseLevel = 0;
@@ -45,6 +45,20 @@ void FastTracker::AddLayer(TString name, float r, float z, float x0, float xrho,
   layers.push_back(newLayer);
 }
 
+DetLayer FastTracker::GetLayer(int layer, bool ignoreBarrelLayers)
+{
+  int layerIdx = layer;
+  if (ignoreBarrelLayers) {
+    for (int il = 0, trackingLayerIdx = 0; trackingLayerIdx <= layer; il++) {
+      if (layers[il].type == 0)
+        continue;
+      trackingLayerIdx++;
+      layerIdx = il;
+    }
+  }
+  return layers[layerIdx];
+}
+
 void FastTracker::Print()
 {
   // print out layer setup
@@ -170,31 +184,58 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa
   std::array<float, 3> posIni; // provision for != PV
   inputTrack.getXYZGlo(posIni);
   float initialRadius = std::hypot(posIni[0], posIni[1]);
+  float kTrackingMargin = 0.1;
+  int xrhosteps = 100;
+  bool applyAngularCorrection = true;
 
   // +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+
   // 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++) {
     // check if layer is doable
     if (layers[il].r < initialRadius)
       continue; // this layer should not be attempted, but go ahead
-    if (layers[il].type == 0)
-      continue; // inert layer, skip
 
     // check if layer is reached
     float targetX = 1e+3;
+    bool ok = true;
     inputTrack.getXatLabR(layers[il].r, targetX, magneticField);
     if (targetX > 999)
       break; // failed to find intercept
 
-    if (!inputTrack.propagateTo(targetX, magneticField)) {
-      break; // failed to propagate
+    ok = inputTrack.propagateTo(targetX, magneticField);
+    if (ok && applyMSCorrection && layers[il].x0 > 0) {
+      ok = inputTrack.correctForMaterial(layers[il].x0, 0, applyAngularCorrection);
+    }
+    if (ok && applyElossCorrection && layers[il].xrho > 0) { // correct in small steps
+      for (int ise = xrhosteps; ise--;) {
+        ok = inputTrack.correctForMaterial(0, -layers[il].xrho / xrhosteps, applyAngularCorrection);
+        if (!ok)
+          break;
+      }
+    }
+
+    // was there a problem on this layer?
+    if (!ok && il > 0) { // may fail to reach target layer due to the eloss
+      float rad2 = inputTrack.getX() * inputTrack.getX() + inputTrack.getY() * inputTrack.getY();
+      float fMinRadTrack = 132.;
+      float maxR = layers[il - 1].r + kTrackingMargin * 2;
+      float minRad = (fMinRadTrack > 0 && fMinRadTrack < maxR) ? fMinRadTrack : maxR;
+      if (rad2 - minRad * minRad < kTrackingMargin * kTrackingMargin) { // check previously reached layer
+        return -5;                                                      // did not reach min requested layer
+      } else {
+        break;
+      }
     }
     if (std::abs(inputTrack.getZ()) > layers[il].z && applyZacceptance) {
       break; // out of acceptance bounds
     }
 
+    if (layers[il].type == 0)
+      continue; // inert layer, skip
+
     // layer is reached
     if (firstLayerReached < 0)
       firstLayerReached = il;
@@ -204,7 +245,7 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa
 
   // +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+
   // initialize track at outer point
-  new (&outputTrack)(o2::track::TrackParCov)(inputTrack);
+  o2::track::TrackParCov inwardTrack(inputTrack);
 
   // Enlarge covariance matrix
   std::array<float, 5> trPars = {0.};
@@ -241,14 +282,12 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa
   largeCov[kSnp2] = largeCov[kTgl2] = kLargeErr2Dir;
   largeCov[kPtI2] = kLargeErr2PtI * trPars[kPtI] * trPars[kPtI];
 
-  outputTrack.setCov(largeCov);
-  outputTrack.checkCovariance();
+  inwardTrack.setCov(largeCov);
+  inwardTrack.checkCovariance();
 
   // +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+
   // Inward pass to calculate covariances
   for (int il = lastLayerReached; il >= firstLayerReached; il--) {
-    if (layers[il].type == 0)
-      continue; // inert layer, skip
 
     float targetX = 1e+3;
     inputTrack.getXatLabR(layers[il].r, targetX, magneticField);
@@ -258,6 +297,7 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa
     if (!inputTrack.propagateTo(targetX, magneticField)) {
       continue; // failed to propagate
     }
+
     if (std::abs(inputTrack.getZ()) > layers[il].z && applyZacceptance) {
       continue; // out of acceptance bounds but continue inwards
     }
@@ -268,20 +308,42 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa
     std::vector<float> thisHit = {spacePoint[0], spacePoint[1], spacePoint[2]};
 
     // towards adding cluster: move to track alpha
-    double alpha = outputTrack.getAlpha();
+    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],
       spacePoint[2]};
-    if (!(outputTrack.propagateTo(xyz1[0], magneticField)))
+    if (!inwardTrack.propagateTo(xyz1[0], magneticField))
       continue;
 
-    const o2::track::TrackParametrization<float>::dim2_t hitpoint = {
-      static_cast<float>(xyz1[1]),
-      static_cast<float>(xyz1[2])};
-    const o2::track::TrackParametrization<float>::dim3_t hitpointcov = {layers[il].resRPhi * layers[il].resRPhi, 0.f, layers[il].resZ * layers[il].resZ};
-    outputTrack.update(hitpoint, hitpointcov);
-    outputTrack.checkCovariance();
+    if (layers[il].type != 0) { // only update covm for tracker hits
+      const o2::track::TrackParametrization<float>::dim2_t hitpoint = {
+        static_cast<float>(xyz1[1]),
+        static_cast<float>(xyz1[2])};
+      const o2::track::TrackParametrization<float>::dim3_t hitpointcov = {layers[il].resRPhi * layers[il].resRPhi, 0.f, layers[il].resZ * layers[il].resZ};
+
+      inwardTrack.update(hitpoint, hitpointcov);
+      inwardTrack.checkCovariance();
+    }
+
+    if (applyMSCorrection && layers[il].x0 > 0) {
+      if (!inputTrack.correctForMaterial(layers[il].x0, 0, applyAngularCorrection)) {
+        return -6;
+      }
+      if (!inwardTrack.correctForMaterial(layers[il].x0, 0, applyAngularCorrection)) {
+        return -6;
+      }
+    }
+    if (applyElossCorrection && layers[il].xrho > 0) {
+      for (int ise = xrhosteps; ise--;) { // correct in small steps
+        if (!inputTrack.correctForMaterial(0, layers[il].xrho / xrhosteps, applyAngularCorrection)) {
+          return -7;
+        }
+        if (!inwardTrack.correctForMaterial(0, layers[il].xrho / xrhosteps, applyAngularCorrection)) {
+          return -7;
+        }
+      }
+    }
 
     if (layers[il].type == 1)
       nSiliconPoints++; // count silicon hits
@@ -293,18 +355,21 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa
 
   // backpropagate to original radius
   float finalX = 1e+3;
-  outputTrack.getXatLabR(initialRadius, finalX, magneticField);
+  inwardTrack.getXatLabR(initialRadius, finalX, magneticField);
   if (finalX > 999)
     return -3; // failed to find intercept
 
-  if (!outputTrack.propagateTo(finalX, magneticField)) {
+  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;
 
+  outputTrack.setCov(inwardTrack.getCov());
+  outputTrack.checkCovariance();
+
   // Use covariance matrix based smearing
   std::array<double, 15> covMat = {0.};
   for (int ii = 0; ii < 15; ii++)

ALICE3/Core/FastTracker.h

@@ -35,6 +35,7 @@ class FastTracker
   virtual ~FastTracker() {}
 
   void AddLayer(TString name, float r, float z, float x0, float xrho, float resRPhi = 0.0f, float resZ = 0.0f, float eff = 0.0f, int type = 0);
+  DetLayer GetLayer(const int layer, bool ignoreBarrelLayers = true);
 
   void AddSiliconALICE3v4();
   void AddSiliconALICE3v1();
@@ -48,10 +49,12 @@ class FastTracker
   std::vector<std::vector<float>> hits; // bookkeep last added hits
 
   // operational
-  float magneticField;   // in kiloGauss (5 = 0.5T, etc)
-  bool applyZacceptance; // check z acceptance or not
-  float covMatFactor;    // covmat off-diagonal factor to use for covmat fix (negative: no factor)
-  int verboseLevel;      // 0: not verbose, >0 more verbose
+  float magneticField;       // in kiloGauss (5 = 0.5T, etc)
+  bool applyZacceptance;     // check z acceptance or not
+  float covMatFactor;        // covmat off-diagonal factor to use for covmat fix (negative: no factor)
+  int verboseLevel;          // 0: not verbose, >0 more verbose
+  bool applyMSCorrection;    // Apply correction for multiple scattering
+  bool applyElossCorrection; // Apply correction for eloss (requires MS correction)
 
   uint64_t covMatOK;    // cov mat has negative eigenvals
   uint64_t covMatNotOK; // cov mat has negative eigenvals