GPUTRDTracker.cxx

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
 
// #define ENABLE_GPUTRDDEBUG
#define ENABLE_WARNING 0
#define ENABLE_INFO 0
 
#include "GPUTRDTracker.h"
#include "GPUTRDTrackletWord.h"
#include "GPUTRDGeometry.h"
#include "GPUTRDTrackerDebug.h"
#include "GPUCommonMath.h"
#include "GPUCommonAlgorithm.h"
#include "GPUConstantMem.h"
#include "GPUTRDRecoParam.h"
 
using namespace o2::gpu;
 
class GPUTPCGMPolynomialField;
 
#ifndef GPUCA_GPUCODE
#include "GPUMemoryResource.h"
#include "GPUReconstruction.h"
#include <chrono>
#include <vector>
 
#include "GPUChainTracking.h"
 
template <class TRDTRK, class PROP>
void GPUTRDTracker_t<TRDTRK, PROP>::SetMaxData(const GPUTrackingInOutPointers& io)
{
  mNMaxTracks = std::max(std::max(io.nOutputTracksTPCO2, io.nTracksTPCITSO2), std::max(io.nMergedTracks, io.nOutputTracksTPCO2)); // TODO: This is a bit stupid, we should just take the correct number, not the max of all
  mNMaxSpacePoints = io.nTRDTracklets;
  mNMaxCollisions = io.nTRDTriggerRecords;
}
 
template <class TRDTRK, class PROP>
void GPUTRDTracker_t<TRDTRK, PROP>::RegisterMemoryAllocation()
{
  AllocateAndInitializeLate();
  mMemoryPermanent = mRec->RegisterMemoryAllocation(this, &GPUTRDTracker_t<TRDTRK, PROP>::SetPointersBase, GPUMemoryResource::MEMORY_PERMANENT, "TRDInitialize");
  mMemoryTracklets = mRec->RegisterMemoryAllocation(this, &GPUTRDTracker_t<TRDTRK, PROP>::SetPointersTracklets, GPUMemoryResource::MEMORY_INPUT, "TRDTracklets");
  mMemoryTracks = mRec->RegisterMemoryAllocation(this, &GPUTRDTracker_t<TRDTRK, PROP>::SetPointersTracks, GPUMemoryResource::MEMORY_INOUT, "TRDTracks");
}
 
template <class TRDTRK, class PROP>
void* GPUTRDTracker_t<TRDTRK, PROP>::SetPointersBase(void* base)
{
  //--------------------------------------------------------------------
  // Allocate memory for fixed size objects (needs to be done only once)
  //--------------------------------------------------------------------
  mMaxBackendThreads = mRec->GetMaxBackendThreads();
  computePointerWithAlignment(base, mR, kNChambers);
  computePointerWithAlignment(base, mHypothesis, mNCandidates * mMaxBackendThreads);
  computePointerWithAlignment(base, mCandidates, mNCandidates * 2 * mMaxBackendThreads);
  return base;
}
 
template <class TRDTRK, class PROP>
void* GPUTRDTracker_t<TRDTRK, PROP>::SetPointersTracklets(void* base)
{
  //--------------------------------------------------------------------
  // Allocate memory for tracklets and space points
  // (size might change for different events)
  //--------------------------------------------------------------------
  if (mGenerateSpacePoints) {
    computePointerWithAlignment(base, mSpacePoints, mNMaxSpacePoints);
  }
  computePointerWithAlignment(base, mTrackletIndexArray, (kNChambers + 1) * mNMaxCollisions);
  return base;
}
 
template <class TRDTRK, class PROP>
void* GPUTRDTracker_t<TRDTRK, PROP>::SetPointersTracks(void* base)
{
  //--------------------------------------------------------------------
  // Allocate memory for tracks (this is done once per event)
  //--------------------------------------------------------------------
  computePointerWithAlignment(base, mTracks, mNMaxTracks);
  computePointerWithAlignment(base, mTrackAttribs, mNMaxTracks);
  return base;
}
 
template <class TRDTRK, class PROP>
GPUTRDTracker_t<TRDTRK, PROP>::GPUTRDTracker_t() : mR(nullptr), mIsInitialized(false), mGenerateSpacePoints(false), mProcessPerTimeFrame(false), mNAngleHistogramBins(25), mAngleHistogramRange(50), mMemoryPermanent(-1), mMemoryTracklets(-1), mMemoryTracks(-1), mNMaxCollisions(0), mNMaxTracks(0), mNMaxSpacePoints(0), mTracks(nullptr), mTrackAttribs(nullptr), mNCandidates(1), mNTracks(0), mNEvents(0), mMaxBackendThreads(100), mTrackletIndexArray(nullptr), mFT0TriggeredBC(nullptr), mNFT0BC(0), mHypothesis(nullptr), mCandidates(nullptr), mSpacePoints(nullptr), mGeo(nullptr), mRecoParam(nullptr), mDebugOutput(false), mMaxEta(0.84f), mRoadZ(18.f), mTPCVdrift(2.58f), mTPCTDriftOffset(0.f), mDebug(new GPUTRDTrackerDebug<TRDTRK>())
{
  //--------------------------------------------------------------------
  // Default constructor
  //--------------------------------------------------------------------
}
 
template <class TRDTRK, class PROP>
GPUTRDTracker_t<TRDTRK, PROP>::~GPUTRDTracker_t()
{
  //--------------------------------------------------------------------
  // Destructor
  //--------------------------------------------------------------------
  delete mDebug;
}
 
template <class TRDTRK, class PROP>
void GPUTRDTracker_t<TRDTRK, PROP>::InitializeProcessor()
{
  //--------------------------------------------------------------------
  // Initialise tracker
  //--------------------------------------------------------------------
  mRecoParam = GetConstantMem()->calibObjects.trdRecoParam;
  UpdateGeometry();
  mDebug->ExpandVectors();
  mIsInitialized = true;
}
 
template <class TRDTRK, class PROP>
void GPUTRDTracker_t<TRDTRK, PROP>::UpdateGeometry()
{
  //--------------------------------------------------------------------
  // Update Geometry of TRDTracker
  //--------------------------------------------------------------------
  mGeo = (const GPUTRDGeometry*)GetConstantMem()->calibObjects.trdGeometry;
  if (!mGeo) {
    GPUFatal("TRD geometry must be provided externally");
  }
  // obtain average radius of TRD chambers
  float x0[kNLayers] = {300.2f, 312.8f, 325.4f, 338.0f, 350.6f, 363.2f}; // used as default value in case no transformation matrix can be obtained
  auto* matrix = mGeo->GetClusterMatrix(0);
  float loc[3] = {mGeo->AnodePos(), 0.f, 0.f};
  float glb[3] = {0.f, 0.f, 0.f};
  for (int32_t iDet = 0; iDet < kNChambers; ++iDet) {
    matrix = mGeo->GetClusterMatrix(iDet);
    if (!matrix) {
      mR[iDet] = x0[mGeo->GetLayer(iDet)];
      continue;
    }
    matrix->LocalToMaster(loc, glb);
    mR[iDet] = glb[0];
  }
}
 
template <class TRDTRK, class PROP>
void GPUTRDTracker_t<TRDTRK, PROP>::Reset()
{
  //--------------------------------------------------------------------
  // Reset tracker
  //--------------------------------------------------------------------
  mNTracks = 0;
}
 
template <class TRDTRK, class PROP>
void GPUTRDTracker_t<TRDTRK, PROP>::PrepareTracking(GPUChainTracking* chainTracking)
{
  //--------------------------------------------------------------------
  // Prepare tracklet index array and if requested calculate space points
  // in part duplicated from DoTracking() method to allow for calling
  // this function on the host prior to GPU processing
  //--------------------------------------------------------------------
  for (uint32_t iColl = 0; iColl < GetConstantMem()->ioPtrs.nTRDTriggerRecords; ++iColl) {
    if (GetConstantMem()->ioPtrs.trdTrigRecMask && GetConstantMem()->ioPtrs.trdTrigRecMask[iColl] == 0) {
      // this trigger is masked as there is no ITS information available for it
      continue;
    }
    int32_t nTrklts = 0;
    int32_t idxOffset = 0;
    if (mProcessPerTimeFrame) {
      idxOffset = GetConstantMem()->ioPtrs.trdTrackletIdxFirst[iColl];
      nTrklts = (iColl < GetConstantMem()->ioPtrs.nTRDTriggerRecords - 1) ? GetConstantMem()->ioPtrs.trdTrackletIdxFirst[iColl + 1] - GetConstantMem()->ioPtrs.trdTrackletIdxFirst[iColl] : GetConstantMem()->ioPtrs.nTRDTracklets - GetConstantMem()->ioPtrs.trdTrackletIdxFirst[iColl];
    } else {
      nTrklts = GetConstantMem()->ioPtrs.nTRDTracklets;
    }
    const GPUTRDTrackletWord* tracklets = &((GetConstantMem()->ioPtrs.trdTracklets)[idxOffset]);
    int32_t* trkltIndexArray = &mTrackletIndexArray[iColl * (kNChambers + 1) + 1];
    trkltIndexArray[-1] = 0;
    int32_t currDet = 0;
    int32_t nextDet = 0;
    int32_t trkltCounter = 0;
    for (int32_t iTrklt = 0; iTrklt < nTrklts; ++iTrklt) {
      if (tracklets[iTrklt].GetDetector() > currDet) {
        nextDet = tracklets[iTrklt].GetDetector();
        for (int32_t iDet = currDet; iDet < nextDet; ++iDet) {
          trkltIndexArray[iDet] = trkltCounter;
        }
        currDet = nextDet;
      }
      ++trkltCounter;
    }
    for (int32_t iDet = currDet; iDet <= kNChambers; ++iDet) {
      trkltIndexArray[iDet] = trkltCounter;
    }
    if (mGenerateSpacePoints) {
      if (!CalculateSpacePoints(iColl)) {
        GPUError("Space points for at least one chamber could not be calculated (for interaction %i)", iColl);
        break;
      }
    }
  }
  if (mGenerateSpacePoints) {
    chainTracking->mIOPtrs.trdSpacePoints = mSpacePoints;
  }
  mNEvents++;
}
 
template <class TRDTRK, class PROP>
void GPUTRDTracker_t<TRDTRK, PROP>::SetNCandidates(int32_t n)
{
  //--------------------------------------------------------------------
  // set the number of candidates to be used
  //--------------------------------------------------------------------
  if (!mIsInitialized) {
    mNCandidates = n;
  } else {
    GPUError("Cannot change mNCandidates after initialization");
  }
}
 
template <class TRDTRK, class PROP>
void GPUTRDTracker_t<TRDTRK, PROP>::PrintSettings() const
{
  //--------------------------------------------------------------------
  // print current settings to screen
  //--------------------------------------------------------------------
  GPUInfo("##############################################################");
  GPUInfo("Current settings for GPU TRD tracker:");
  GPUInfo(" maxChi2(%.2f), chi2Penalty(%.2f), nCandidates(%i), maxMissingLayers(%i)", Param().rec.trd.maxChi2, Param().rec.trd.penaltyChi2, mNCandidates, Param().rec.trd.stopTrkAfterNMissLy);
  GPUInfo(" ptCut = %.2f GeV, abs(eta) < %.2f", Param().rec.trd.minTrackPt, mMaxEta);
  GPUInfo("##############################################################");
}
 
template <class TRDTRK, class PROP>
void GPUTRDTracker_t<TRDTRK, PROP>::StartDebugging()
{
  mDebug->CreateStreamer();
}
 
#endif 
 
template <>
GPUdi() const GPUTRDPropagatorGPU::propagatorParam* GPUTRDTracker_t<GPUTRDTrackGPU, GPUTRDPropagatorGPU>::getPropagatorParam()
{
  return &Param().polynomialField;
}
 
template <class TRDTRK, class PROP>
GPUdi() const typename PROP::propagatorParam* GPUTRDTracker_t<TRDTRK, PROP>::getPropagatorParam()
{
  return GetConstantMem()->calibObjects.o2Propagator;
}
 
template <class TRDTRK, class PROP>
GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::CheckTrackTRDCandidate(const TRDTRK& trk) const
{
  if (!trk.CheckNumericalQuality()) {
    return false;
  }
  if (CAMath::Abs(trk.getEta()) > mMaxEta) {
    return false;
  }
  if (trk.getPt() < Param().rec.trd.minTrackPt) {
    return false;
  }
  return true;
}
 
template <class TRDTRK, class PROP>
GPUd() int32_t GPUTRDTracker_t<TRDTRK, PROP>::LoadTrack(const TRDTRK& trk, uint32_t tpcTrackId, bool checkTrack, HelperTrackAttributes* attribs)
{
  if (mNTracks >= mNMaxTracks) {
#ifndef GPUCA_GPUCODE
    GPUError("Error: Track dropped (no memory available) -> must not happen");
#endif
    return (1);
  }
  if (checkTrack && !CheckTrackTRDCandidate(trk)) {
    return 2;
  }
  mTracks[mNTracks] = trk;
  mTracks[mNTracks].setRefGlobalTrackIdRaw(tpcTrackId);
  if (attribs) {
    mTrackAttribs[mNTracks] = *attribs;
  }
  mNTracks++;
  return (0);
}
 
template <class TRDTRK, class PROP>
GPUd() void GPUTRDTracker_t<TRDTRK, PROP>::DumpTracks()
{
  //--------------------------------------------------------------------
  // helper function (only for debugging purposes)
  //--------------------------------------------------------------------
  GPUInfo("There are in total %i tracklets loaded", GetConstantMem()->ioPtrs.nTRDTracklets);
  GPUInfo("There are %i tracks loaded. mNMaxTracks(%i)", mNTracks, mNMaxTracks);
  for (int32_t i = 0; i < mNTracks; ++i) {
    auto* trk = &(mTracks[i]);
    GPUInfo("track %i: x=%f, alpha=%f, nTracklets=%i, pt=%f, time=%f", i, trk->getX(), trk->getAlpha(), trk->getNtracklets(), trk->getPt(), mTrackAttribs[i].mTime);
  }
}
 
template <class TRDTRK, class PROP>
GPUd() int32_t GPUTRDTracker_t<TRDTRK, PROP>::GetCollisionIDs(int32_t iTrk, int32_t* collisionIds) const
{
  //--------------------------------------------------------------------
  // Check which TRD trigger times possibly match given input track.
  // If ITS-TPC matches or CE-crossing TPC tracks the time is precisely
  // known and max 1 trigger time can be assigned.
  // For TPC-only tracks the collision IDs are stored in collisionIds array
  // and the number of valid entries in the array is returned
  //--------------------------------------------------------------------
  int32_t nColls = 0;
  for (uint32_t iColl = 0; iColl < GetConstantMem()->ioPtrs.nTRDTriggerRecords; ++iColl) {
    if (GetConstantMem()->ioPtrs.trdTrigRecMask && GetConstantMem()->ioPtrs.trdTrigRecMask[iColl] == 0) {
      continue;
    }
    if (GetConstantMem()->ioPtrs.trdTriggerTimes[iColl] > mTrackAttribs[iTrk].GetTimeMin() && GetConstantMem()->ioPtrs.trdTriggerTimes[iColl] < mTrackAttribs[iTrk].GetTimeMax()) {
      if (nColls == 20) {
        GPUError("Found too many collision candidates for track with tMin(%f) and tMax(%f)", mTrackAttribs[iTrk].GetTimeMin(), mTrackAttribs[iTrk].GetTimeMax());
        return nColls;
      }
      collisionIds[nColls++] = iColl;
    }
  }
  return nColls;
}
 
template <class TRDTRK, class PROP>
GPUd() void GPUTRDTracker_t<TRDTRK, PROP>::DoTrackingThread(int32_t iTrk, int32_t threadId)
{
  //--------------------------------------------------------------------
  // perform the tracking for one track (must be threadsafe)
  //--------------------------------------------------------------------
  int32_t collisionIds[20] = {0}; // due to the dead time there will never exist more possible TRD triggers for a single track
  int32_t nCollisionIds = 1;      // initialize with 1 for AliRoot compatibility
  if (mProcessPerTimeFrame) {
    nCollisionIds = GetCollisionIDs(iTrk, collisionIds);
    if (nCollisionIds == 0) {
      if (ENABLE_INFO) {
        GPUInfo("Did not find TRD data for track %i with t=%f. tMin(%f), tMax(%f)", iTrk, mTrackAttribs[iTrk].mTime, mTrackAttribs[iTrk].GetTimeMin(), mTrackAttribs[iTrk].GetTimeMax());
      }
      // no TRD data available for the bunch crossing this track originates from
      return;
    }
  }
  PROP prop(getPropagatorParam());
  mTracks[iTrk].setChi2(Param().rec.trd.penaltyChi2); // TODO check if this should not be higher
 
  auto trkStart = mTracks[iTrk];
  for (int32_t iColl = 0; iColl < nCollisionIds; ++iColl) {
    // do track following for each collision candidate and keep best track
    auto trkCopy = trkStart;
    prop.setTrack(&trkCopy);
    prop.setFitInProjections(true);
    if (!FollowProlongation(&prop, &trkCopy, iTrk, threadId, collisionIds[iColl])) {
      // track following failed
      continue;
    }
    if (trkCopy.getReducedChi2() < mTracks[iTrk].getReducedChi2()) {
      mTracks[iTrk] = trkCopy; // copy back the resulting track
    }
  }
}
 
template <class TRDTRK, class PROP>
GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::CalculateSpacePoints(int32_t iCollision)
{
  //--------------------------------------------------------------------
  // Calculates TRD space points in sector tracking coordinates
  // from online tracklets
  //--------------------------------------------------------------------
 
  bool result = true;
  int32_t idxOffset = iCollision * (kNChambers + 1); // offset for accessing mTrackletIndexArray for collision iCollision
 
  const GPUTRDTrackletWord* tracklets = GetConstantMem()->ioPtrs.trdTracklets;
 
  for (int32_t iDet = 0; iDet < kNChambers; ++iDet) {
    int32_t iFirstTrackletInDet = mTrackletIndexArray[idxOffset + iDet];
    int32_t iFirstTrackletInNextDet = mTrackletIndexArray[idxOffset + iDet + 1];
    int32_t nTrackletsInDet = iFirstTrackletInNextDet - iFirstTrackletInDet;
    if (nTrackletsInDet == 0) {
      continue;
    }
    if (!mGeo->ChamberInGeometry(iDet)) {
      GPUError("Found TRD tracklets in chamber %i which is not included in the geometry", iDet);
      return false;
    }
    auto* matrix = mGeo->GetClusterMatrix(iDet);
    if (!matrix) {
      GPUError("No cluster matrix available for chamber %i. Skipping it...", iDet);
      result = false;
      continue;
    }
    const GPUTRDpadPlane* pp = mGeo->GetPadPlane(iDet);
 
    int32_t trkltIdxOffset = (mProcessPerTimeFrame) ? GetConstantMem()->ioPtrs.trdTrackletIdxFirst[iCollision] : 0; // global index of first tracklet in iCollision
    int32_t trkltIdxStart = trkltIdxOffset + iFirstTrackletInDet;
    for (int32_t trkltIdx = trkltIdxStart; trkltIdx < trkltIdxStart + nTrackletsInDet; ++trkltIdx) {
      int32_t trkltZbin = tracklets[trkltIdx].GetZbin();
      float xTrkltDet[3] = {0.f}; // trklt position in chamber coordinates
      float xTrkltSec[3] = {0.f}; // trklt position in sector coordinates
      xTrkltDet[0] = mGeo->AnodePos() + sRadialOffset;
      xTrkltDet[1] = tracklets[trkltIdx].GetY();
      xTrkltDet[2] = pp->GetRowPos(trkltZbin) - pp->GetRowSize(trkltZbin) / 2.f - pp->GetRowPos(pp->GetNrows() / 2);
      // GPUInfo("Space point local %i: x=%f, y=%f, z=%f", trkltIdx, xTrkltDet[0], xTrkltDet[1], xTrkltDet[2]);
      matrix->LocalToMaster(xTrkltDet, xTrkltSec);
      mSpacePoints[trkltIdx].setX(xTrkltSec[0]);
      mSpacePoints[trkltIdx].setY(xTrkltSec[1]);
      mSpacePoints[trkltIdx].setZ(xTrkltSec[2]);
      mSpacePoints[trkltIdx].setDy(tracklets[trkltIdx].GetdY());
 
      // GPUInfo("Space point global %i: x=%f, y=%f, z=%f", trkltIdx, mSpacePoints[trkltIdx].getX(), mSpacePoints[trkltIdx].getY(), mSpacePoints[trkltIdx].getZ());
    }
  }
  return result;
}
 
template <class TRDTRK, class PROP>
GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK* t, int32_t iTrk, int32_t threadId, int32_t collisionId)
{
  //--------------------------------------------------------------------
  // Propagate TPC track layerwise through TRD and pick up closest
  // tracklet(s) on the way
  // -> returns false if prolongation could not be executed fully
  //    or track does not fullfill threshold conditions
  //--------------------------------------------------------------------
 
  if (ENABLE_INFO) {
    GPUInfo("Start track following for track %i at x=%f with pt=%f", t->getRefGlobalTrackIdRaw(), t->getX(), t->getPt());
  }
  mDebug->Reset();
  t->setChi2(0.f);
 
  // Find compatible BC ids
  int32_t nIdxBCMin = -1;
  int32_t nIdxBCMax = -1;
 
  for (int32_t iBC = 0; iBC < mNFT0BC; iBC++) {
    int32_t deltaBC = CAMath::Round(mFT0TriggeredBC[iBC] - GetConstantMem()->ioPtrs.trdTriggerTimes[collisionId] / o2::constants::lhc::LHCBunchSpacingMUS);
    if (nIdxBCMin == -1 && deltaBC > mRecoParam->getPileUpRangeBefore()) {
      nIdxBCMin = iBC;
    }
    if (deltaBC >= mRecoParam->getPileUpRangeAfter()) {
      nIdxBCMax = iBC;
      break;
    }
    if (iBC == mNFT0BC - 1) {
      nIdxBCMax = iBC + 1;
      if (nIdxBCMin == -1) {
        // we did not find the correct BC, so we don't do any pile-up correction
        nIdxBCMin = nIdxBCMax;
      }
    }
  }
 
  float zShiftTrk = 0.f;
  if (mProcessPerTimeFrame) {
    zShiftTrk = (mTrackAttribs[iTrk].mTime - GetConstantMem()->ioPtrs.trdTriggerTimes[collisionId]) * mTPCVdrift * mTrackAttribs[iTrk].mSide;
    // float addZerr = (mTrackAttribs[iTrk].mTimeAddMax + mTrackAttribs[iTrk].mTimeSubMax) * .5f * mTPCVdrift;
    //  increase Z error based on time window
    //  -> this is here since it was done before, but the efficiency seems to be better if the covariance is not updated (more tracklets are attached)
    // t->updateCovZ2(addZerr * addZerr); // TODO check again once detailed performance study tools are available, maybe this can be tuned
  }
  const GPUTRDpadPlane* pad = nullptr;
  const GPUTRDTrackletWord* tracklets = GetConstantMem()->ioPtrs.trdTracklets;
  const GPUTRDSpacePoint* spacePoints = GetConstantMem()->ioPtrs.trdSpacePoints;
 
#ifdef ENABLE_GPUTRDDEBUG
  TRDTRK trackNoUp(*t);
#endif
 
  int32_t candidateIdxOffset = threadId * 2 * mNCandidates;
  int32_t hypothesisIdxOffset = threadId * mNCandidates;
  int32_t trkltIdxOffset = collisionId * (kNChambers + 1);                                                            // offset for accessing mTrackletIndexArray for given collision
  int32_t glbTrkltIdxOffset = (mProcessPerTimeFrame) ? GetConstantMem()->ioPtrs.trdTrackletIdxFirst[collisionId] : 0; // offset of first tracklet in given collision in global tracklet array
 
  auto trkWork = t;
  if (mNCandidates > 1) {
    // copy input track to first candidate
    mCandidates[candidateIdxOffset] = *t;
  }
 
  int32_t nCandidates = 1;     // we always start with one candidate
  int32_t nCurrHypothesis = 0; // the number of track hypothesis in given iLayer
 
  // search window
  float roadY = 0.f;
  float roadZ = 0.f;
  const int32_t nMaxChambersToSearch = 4;
 
  mDebug->SetGeneralInfo(mNEvents, mNTracks, iTrk, t->getPt());
 
  for (int32_t iLayer = 0; iLayer < kNLayers; ++iLayer) {
    nCurrHypothesis = 0;
    bool isOK = false; // if at least one candidate could be propagated or the track was stopped this becomes true
    int32_t currIdx = candidateIdxOffset + iLayer % 2;
    int32_t nextIdx = candidateIdxOffset + (iLayer + 1) % 2;
    pad = mGeo->GetPadPlane(iLayer, 0);
    float tilt = CAMath::Tan(CAMath::Pi() / 180.f * pad->GetTiltingAngle()); // tilt is signed!
    const float zMaxTRD = pad->GetRow0();
 
    // --------------------------------------------------------------------------------
    //
    // for each candidate, propagate to layer radius and look for close tracklets
    //
    // --------------------------------------------------------------------------------
    for (int32_t iCandidate = 0; iCandidate < nCandidates; iCandidate++) {
 
      int32_t det[nMaxChambersToSearch] = {-1, -1, -1, -1}; // TRD chambers to be searched for tracklets
 
      if (mNCandidates > 1) {
        trkWork = &mCandidates[2 * iCandidate + currIdx];
        prop->setTrack(trkWork);
      }
 
      if (trkWork->getIsStopped()) {
        Hypothesis hypo(trkWork->getNlayersFindable(), iCandidate, -1, trkWork->getChi2());
        InsertHypothesis(hypo, nCurrHypothesis, hypothesisIdxOffset);
        isOK = true;
        continue;
      }
 
      // propagate track to average radius of TRD layer iLayer (sector 0, stack 2 is chosen as a reference)
      if (!prop->propagateToX(mR[2 * kNLayers + iLayer], .8f, 2.f)) {
        if (ENABLE_INFO) {
          GPUInfo("Track propagation failed for track %i candidate %i in layer %i (pt=%f, x=%f, mR[layer]=%f)", iTrk, iCandidate, iLayer, trkWork->getPt(), trkWork->getX(), mR[2 * kNLayers + iLayer]);
        }
        continue;
      }
 
      // rotate track in new sector in case of sector crossing
      if (!AdjustSector(prop, trkWork)) {
        if (ENABLE_INFO) {
          GPUInfo("Adjusting sector failed for track %i candidate %i in layer %i", iTrk, iCandidate, iLayer);
        }
        continue;
      }
 
      // check if track is findable
      if (IsGeoFindable(trkWork, iLayer, prop->getAlpha(), zShiftTrk)) {
        trkWork->setIsFindable(iLayer);
      }
 
      // define search window
      roadY = 7.f * CAMath::Sqrt(trkWork->getSigmaY2() + 0.1f * 0.1f) + Param().rec.trd.extraRoadY; // add constant to the road to account for uncertainty due to radial deviations (few mm)
      // roadZ = 7.f * CAMath::Sqrt(trkWork->getSigmaZ2() + 9.f * 9.f / 12.f); // take longest pad length
      roadZ = mRoadZ + Param().rec.trd.extraRoadZ; // simply twice the longest pad length -> efficiency 99.996%
      //
      if (CAMath::Abs(trkWork->getZ() + zShiftTrk) - roadZ >= zMaxTRD) {
        if (ENABLE_INFO) {
          GPUInfo("Track out of TRD acceptance with z=%f in layer %i (eta=%f)", trkWork->getZ() + zShiftTrk, iLayer, trkWork->getEta());
        }
        continue;
      }
 
      // determine chamber(s) to be searched for tracklets
      FindChambersInRoad(trkWork, roadY, roadZ, iLayer, det, zMaxTRD, prop->getAlpha(), zShiftTrk);
 
      // track debug information to be stored in case no matching tracklet can be found
      mDebug->SetTrackParameter(*trkWork, iLayer);
 
      // look for tracklets in chamber(s)
      for (int32_t iDet = 0; iDet < nMaxChambersToSearch; iDet++) {
        int32_t currDet = det[iDet];
        if (currDet == -1) {
          continue;
        }
        pad = mGeo->GetPadPlane(currDet);
        int32_t currSec = mGeo->GetSector(currDet);
        if (currSec != GetSector(prop->getAlpha())) {
          if (!prop->rotate(GetAlphaOfSector(currSec))) {
            if (ENABLE_WARNING) {
              GPUWarning("Track could not be rotated in tracklet coordinate system");
            }
            break;
          }
        }
        if (currSec != GetSector(prop->getAlpha())) {
          GPUError("Track is in sector %i and sector %i is searched for tracklets", GetSector(prop->getAlpha()), currSec);
          continue;
        }
        // propagate track to radius of chamber
        if (!prop->propagateToX(mR[currDet], .8f, .2f)) {
          if (ENABLE_WARNING) {
            GPUWarning("Track parameter for track %i, x=%f at chamber %i x=%f in layer %i cannot be retrieved", iTrk, trkWork->getX(), currDet, mR[currDet], iLayer);
          }
        }
 
        // first propagate track to x of tracklet
        for (int32_t trkltIdx = glbTrkltIdxOffset + mTrackletIndexArray[trkltIdxOffset + currDet]; trkltIdx < glbTrkltIdxOffset + mTrackletIndexArray[trkltIdxOffset + currDet + 1]; ++trkltIdx) {
          if (CAMath::Abs(trkWork->getY() - spacePoints[trkltIdx].getY()) > roadY || CAMath::Abs(trkWork->getZ() + zShiftTrk - spacePoints[trkltIdx].getZ()) > roadZ) {
            // skip tracklets which are too far away
            // although the radii of space points and tracks may differ by ~ few mm the roads are large enough to allow no efficiency loss by this cut
            continue;
          }
          float projY, projZ;
          prop->getPropagatedYZ(spacePoints[trkltIdx].getX(), projY, projZ);
          // correction for tilted pads (only applied if deltaZ < lPad && track z err << lPad)
          float tiltCorr = tilt * (spacePoints[trkltIdx].getZ() - projZ);
          float dyTiltCorr = tilt * trkWork->getTgl() * mGeo->GetCdrHght();
          float lPad = pad->GetRowSize(tracklets[trkltIdx].GetZbin());
          if (!((CAMath::Abs(spacePoints[trkltIdx].getZ() - projZ) < lPad) && (trkWork->getSigmaZ2() < (lPad * lPad / 12.f)))) {
            tiltCorr = 0.f; // will be zero also for TPC tracks which are shifted in z
            dyTiltCorr = 0.f;
          }
 
          // Correction for pile-up: if the track comes from a pile-up event, it should not be extrapolated to the anode plane at t=0,
          // contrary to the assumption when tracklet is reconstructed, so we cancel this extrapolation.
          // The correction is extracted from the most probable trigger. There is also an additional error depending on all the compatible triggers
          float yCorrPileUp = 0.f;
          float yAddErrPileUp2 = 0.f;
          // additional zshiftTrk for TPC only tracks due to pile-up
          float zShiftTrkPileUp = 0.f;
          if (nIdxBCMax - nIdxBCMin >= 2) {
            float maxProb = 0.f;
            // The uncertainty is the RMS wrt the default correction of all possible corrections weighted by their probability
            float sumCorr = 0.f;
            float sumCorr2 = 0.f;
            float sumProb = 0.f;
            // conversion from slope in pad per time bin to slope in cm per BC = tracklets[trkltIdx].getSlopeFloat() * padWidth / BCperTimeBin
            float slopeFactor = tracklets[trkltIdx].GetSlopeFloat() * mGeo->GetPadPlaneWidthIPad(tracklets[trkltIdx].GetDetector()) / 4.f;
            for (int32_t iBC = nIdxBCMin; iBC < nIdxBCMax; iBC++) {
              int32_t deltaBC = CAMath::Round(mFT0TriggeredBC[iBC] - GetConstantMem()->ioPtrs.trdTriggerTimes[collisionId] / o2::constants::lhc::LHCBunchSpacingMUS);
              float probBC = mRecoParam->getPileUpProbTracklet(deltaBC, true, (tracklets[trkltIdx].GetQ0() != 0), (tracklets[trkltIdx].GetQ1() != 0));
              sumCorr += probBC * slopeFactor * deltaBC;
              sumCorr2 += probBC * slopeFactor * deltaBC * slopeFactor * deltaBC;
              sumProb += probBC;
              if (probBC > maxProb) {
                maxProb = probBC;
                yCorrPileUp = -slopeFactor * deltaBC;
                zShiftTrkPileUp = -deltaBC * o2::constants::lhc::LHCBunchSpacingMUS * mTPCVdrift * mTrackAttribs[iTrk].mSide;
              }
            }
            if (sumProb > 1e-6f) {
              yAddErrPileUp2 = sumCorr2 / sumProb - 2 * yCorrPileUp * sumCorr / sumProb + yCorrPileUp * yCorrPileUp;
            }
          }
 
          // number of tracklets within the chamber is the current TRD occupancy estimator
          int nTrackletsChamber = mTrackletIndexArray[trkltIdxOffset + currDet + 1] - mTrackletIndexArray[trkltIdxOffset + currDet];
          float angularPull = GetAngularPull(spacePoints[trkltIdx].getDy() + dyTiltCorr, trkWork->getSnp(), nTrackletsChamber);
 
          // correction for mean z position of tracklet (is not the center of the pad if track eta != 0)
          float zPosCorr = spacePoints[trkltIdx].getZ() + mRecoParam->getZCorrCoeffNRC() * trkWork->getTgl();
          float yPosCorr = spacePoints[trkltIdx].getY() - tiltCorr + yCorrPileUp;
          zPosCorr -= zShiftTrk + zShiftTrkPileUp; // shift tracklet instead of track in order to avoid having to do a re-fit for each collision
 
          // Correction of y position based on angular pull
          if (Param().rec.trd.useAngularPull == 3 || Param().rec.trd.useAngularPull == 4) {
            float corrPull = -angularPull * mRecoParam->getCorrYDy(trkWork->getSnp());
            yPosCorr += corrPull;
          }
 
          float deltaY = yPosCorr - projY;
          float deltaZ = zPosCorr - projZ;
 
          float trkltPosTmpYZ[2] = {yPosCorr, zPosCorr};
          float trkltCovTmp[3] = {0.f};
          if ((CAMath::Abs(deltaY) < roadY) && (CAMath::Abs(deltaZ) < roadZ)) { // TODO: check if this is still necessary after the cut before propagation of track
            // tracklet is in window: get predicted chi2 for update and store tracklet index if best guess
            RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(tracklets[trkltIdx].GetZbin()), (Param().rec.trd.useAngularPull == 2 || Param().rec.trd.useAngularPull == 4 ? angularPull : 0.f), nTrackletsChamber, trkltCovTmp);
            trkltCovTmp[0] += yAddErrPileUp2;
            float chi2 = prop->getPredictedChi2(trkltPosTmpYZ, trkltCovTmp);
 
            if (Param().rec.trd.addDeflectionInChi2 >= 1 && (trkWork->getSnp() < 1.f - 1e-6f) && (trkWork->getSnp() > -1.f + 1e-6f)) {
              // we add the slope in the chi2 calculation
              float trkltCovTmpWithDy[6] = {trkltCovTmp[0], trkltCovTmp[1], trkltCovTmp[2], 0.f, 0.f, 0.f};
              RecalcTrkltCovDy(tilt, trkWork->getSnp(), (Param().rec.trd.useAngularPull == 2 || Param().rec.trd.useAngularPull == 4 ? angularPull : 0.f), nTrackletsChamber, trkltCovTmpWithDy);
              trkltCovTmpWithDy[0] += trkWork->getSigmaY2();
              trkltCovTmpWithDy[1] += trkWork->getSigmaZY();
              trkltCovTmpWithDy[2] += trkWork->getSigmaZ2();
 
              if (Param().rec.trd.useAngularPull == 3 || Param().rec.trd.useAngularPull == 4) {
                // In this case the correlation between tracklet y and tracklet dy is already corrected for so we can cancel it
                trkltCovTmpWithDy[3] = 0.;
                trkltCovTmpWithDy[4] = 0.;
              }
 
              // We add the correlation between track y and track dy
              trkltCovTmpWithDy[3] += trkWork->getSigmaSnpY() * mGeo->GetCdrHght();
              trkltCovTmpWithDy[4] += trkWork->getSigmaSnpZ() * mGeo->GetCdrHght();
 
              // keep these parameters before matrix inversion
              float sigmaZ2 = trkltCovTmpWithDy[2];
              float sigmaDy2 = trkltCovTmpWithDy[5];
 
              // For now, dy uncertainty parametrization also includes track uncertainty, so no need to add additional uncertainty
              if (InvertCov(trkltCovTmpWithDy)) {
                float deltaDy = spacePoints[trkltIdx].getDy() + dyTiltCorr - mRecoParam->convertAngleToDy(trkWork->getSnp());
                if (Param().rec.trd.addDeflectionInChi2 == 2 || Param().rec.trd.addDeflectionInChi2 == 3) {
                  // In this case we take into account the full likelihood, so we replace (deltaDy/sigmaDy)^2 by -2*ln(likelihood), which is the same in the default Gaussian case
                  double likelihood = mRecoParam->getDyLikelihood(trkWork->getSnp(), spacePoints[trkltIdx].getDy() + dyTiltCorr, nTrackletsChamber);
                  if (likelihood < 1e-6f) {
                    continue; // likelihood of 1e-6 is equivalent to 5 sigma deviation, so we can safely cut it to avoid numerical instability in log calculation
                  }
                  deltaDy = CAMath::Sqrt(-2.f * CAMath::Log(likelihood) * sigmaDy2) * (deltaDy > 0.f ? 1.f : -1.f);
                }
                if (Param().rec.trd.addDeflectionInChi2 == 3) {
                  // We do the same for deltaZ
                  double likelihood = mRecoParam->getZLikelihood(deltaZ, pad->GetRowSize(tracklets[trkltIdx].GetZbin()), CAMath::Sqrt(trkWork->getSigmaZ2()));
                  if (likelihood < 1e-6f) {
                    continue;
                  }
                  deltaZ = CAMath::Sqrt(-2.f * CAMath::Log(likelihood) * sigmaZ2) * (deltaZ > 0.f ? 1.f : -1.f);
                }
                chi2 = deltaY * trkltCovTmpWithDy[0] * deltaY + 2 * deltaY * trkltCovTmpWithDy[1] * deltaZ + 2 * deltaY * trkltCovTmpWithDy[3] * deltaDy + deltaZ * trkltCovTmpWithDy[2] * deltaZ + 2 * deltaZ * trkltCovTmpWithDy[4] * deltaDy + deltaDy * trkltCovTmpWithDy[5] * deltaDy;
              }
            }
            // TODO cut on angular pull should be made stricter when proper v-drift calibration for the TRD tracklets is implemented
            if ((chi2 > Param().rec.trd.maxChi2) || (Param().rec.trd.applyDeflectionCut && CAMath::Abs(angularPull) > 4)) {
              continue;
            }
            Hypothesis hypo(trkWork->getNlayersFindable(), iCandidate, trkltIdx, trkWork->getChi2() + chi2);
            InsertHypothesis(hypo, nCurrHypothesis, hypothesisIdxOffset);
          } // end tracklet in window
        } // tracklet loop
      } // chamber loop
 
      // add no update to hypothesis list
      Hypothesis hypoNoUpdate(trkWork->getNlayersFindable(), iCandidate, -1, trkWork->getChi2() + Param().rec.trd.penaltyChi2);
      InsertHypothesis(hypoNoUpdate, nCurrHypothesis, hypothesisIdxOffset);
      isOK = true;
    } // end candidate loop
 
    mDebug->SetChi2Update(mHypothesis[0 + hypothesisIdxOffset].mChi2 - t->getChi2(), iLayer); // only meaningful for ONE candidate!!!
    mDebug->SetRoad(roadY, roadZ, iLayer);                                                    // only meaningful for ONE candidate
    bool wasTrackStored = false;
    // --------------------------------------------------------------------------------
    //
    // loop over the best N_candidates hypothesis
    //
    // --------------------------------------------------------------------------------
    // GPUInfo("nCurrHypothesis=%i, nCandidates=%i", nCurrHypothesis, nCandidates);
    // for (int32_t idx=0; idx<10; ++idx) { GPUInfo("mHypothesis[%i]: candidateId=%i, nLayers=%i, trackletId=%i, chi2=%f", idx, mHypothesis[idx].mCandidateId,  mHypothesis[idx].mLayers, mHypothesis[idx].mTrackletId, mHypothesis[idx].mChi2); }
    for (int32_t iUpdate = 0; iUpdate < nCurrHypothesis && iUpdate < mNCandidates; iUpdate++) {
      if (mHypothesis[iUpdate + hypothesisIdxOffset].mCandidateId == -1) {
        // no more candidates
        if (iUpdate == 0) {
          return false; // no valid candidates for this track (probably propagation failed)
        }
        break; // go to next layer
      }
      nCandidates = iUpdate + 1;
      if (mNCandidates > 1) {
        mCandidates[2 * iUpdate + nextIdx] = mCandidates[2 * mHypothesis[iUpdate + hypothesisIdxOffset].mCandidateId + currIdx];
        trkWork = &mCandidates[2 * iUpdate + nextIdx];
      }
      if (mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId == -1) {
        // no matching tracklet found
        if (trkWork->getIsFindable(iLayer)) {
          if (trkWork->getNmissingConsecLayers(iLayer) > Param().rec.trd.stopTrkAfterNMissLy) {
            trkWork->setIsStopped();
          }
          trkWork->setChi2(trkWork->getChi2() + Param().rec.trd.penaltyChi2);
        }
        if (iUpdate == 0 && mNCandidates > 1) { // TODO: is thie really necessary????? CHECK!
          *t = mCandidates[2 * iUpdate + nextIdx];
        }
        continue;
      }
      // matching tracklet found
      if (mNCandidates > 1) {
        prop->setTrack(trkWork);
      }
      int32_t trkltSec = mGeo->GetSector(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector());
      if (trkltSec != GetSector(prop->getAlpha())) {
        // if after a matching tracklet was found another sector was searched for tracklets the track needs to be rotated back
        prop->rotate(GetAlphaOfSector(trkltSec));
      }
      if (!prop->propagateToX(spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getX(), .8f, 2.f)) {
        if (ENABLE_WARNING) {
          GPUWarning("Final track propagation for track %i update %i in layer %i failed", iTrk, iUpdate, iLayer);
        }
        trkWork->setChi2(trkWork->getChi2() + Param().rec.trd.penaltyChi2);
        if (trkWork->getIsFindable(iLayer)) {
          if (trkWork->getNmissingConsecLayers(iLayer) >= Param().rec.trd.stopTrkAfterNMissLy) {
            trkWork->setIsStopped();
          }
        }
        if (iUpdate == 0 && mNCandidates > 1) {
          *t = mCandidates[2 * iUpdate + nextIdx];
        }
        continue;
      }
 
      pad = mGeo->GetPadPlane(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector());
      float tiltCorrUp = tilt * (spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getZ() - trkWork->getZ());
      float dyTiltCorr = tilt * trkWork->getTgl() * mGeo->GetCdrHght();
 
      float yPosCorrUp = spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getY() - tiltCorrUp;
 
      float zPosCorrUp = spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getZ() + mRecoParam->getZCorrCoeffNRC() * trkWork->getTgl();
      zPosCorrUp -= zShiftTrk;
      float padLength = pad->GetRowSize(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetZbin());
      if (!((trkWork->getSigmaZ2() < (padLength * padLength / 12.f)) && (CAMath::Abs(spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getZ() - trkWork->getZ()) < padLength))) {
        tiltCorrUp = 0.f;
        dyTiltCorr = 0.f;
      }
 
      // Correction for pile-up: if the track comes from a pile-up event, it should not be extrapolated to the anode plane at t=0,
      // contrary to the assumption when tracklet is reconstructed, so we cancel this extrapolation.
      // The correction is extracted from the most probable trigger. There is also an additional error depending on all the compatible triggers
      float yCorrPileUp = 0.f;
      float yAddErrPileUp2 = 0.f;
      float zShiftTrkPileUp = 0.f;
      if (nIdxBCMax - nIdxBCMin >= 2) {
        float maxProb = 0.f;
        // The uncertainty is the RMS wrt the default correction of all possible corrections weighted by their probability
        float sumCorr = 0.f;
        float sumCorr2 = 0.f;
        float sumProb = 0.f;
        // conversion from slope in pad per time bin to slope in cm per BC = tracklets[trkltIdx].getSlopeFloat() * padWidth / BCperTimeBin
        float slopeFactor = tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetSlopeFloat() * mGeo->GetPadPlaneWidthIPad(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector()) / 4.f;
        for (int32_t iBC = nIdxBCMin; iBC < nIdxBCMax; iBC++) {
          int32_t deltaBC = CAMath::Round(mFT0TriggeredBC[iBC] - GetConstantMem()->ioPtrs.trdTriggerTimes[collisionId] / o2::constants::lhc::LHCBunchSpacingMUS);
          float probBC = mRecoParam->getPileUpProbTracklet(deltaBC, true, (tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetQ0() != 0), (tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetQ1() != 0));
          sumCorr += probBC * slopeFactor * deltaBC;
          sumCorr2 += probBC * slopeFactor * deltaBC * slopeFactor * deltaBC;
          sumProb += probBC;
          if (probBC > maxProb) {
            maxProb = probBC;
            yCorrPileUp = -slopeFactor * deltaBC;
            zShiftTrkPileUp = -deltaBC * o2::constants::lhc::LHCBunchSpacingMUS * mTPCVdrift * mTrackAttribs[iTrk].mSide;
          }
        }
        if (sumProb > 1e-6f) {
          yAddErrPileUp2 = sumCorr2 / sumProb - 2 * yCorrPileUp * sumCorr / sumProb + yCorrPileUp * yCorrPileUp;
        }
      }
 
      zPosCorrUp -= zShiftTrkPileUp;
      yPosCorrUp += yCorrPileUp;
 
      const auto currDet = tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector();
      int nTrackletsChamber = mTrackletIndexArray[trkltIdxOffset + currDet + 1] - mTrackletIndexArray[trkltIdxOffset + currDet];
      float angularPull = GetAngularPull(spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getDy() + dyTiltCorr, trkWork->getSnp(), nTrackletsChamber);
 
      // Correction of y position based on angular pull
      if (Param().rec.trd.useAngularPull == 3 || Param().rec.trd.useAngularPull == 4) {
        float corrPull = -angularPull * mRecoParam->getCorrYDy(trkWork->getSnp());
        yPosCorrUp += corrPull;
      }
 
      float trkltPosUp[2] = {yPosCorrUp, zPosCorrUp};
      float trkltCovUp[3] = {0.f};
      RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetZbin()), ((Param().rec.trd.useAngularPull != 0) ? angularPull : 0.f), nTrackletsChamber, trkltCovUp);
      trkltCovUp[0] += yAddErrPileUp2;
 
#ifdef ENABLE_GPUTRDDEBUG
      prop->setTrack(&trackNoUp);
      prop->rotate(GetAlphaOfSector(trkltSec));
      // prop->propagateToX(spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getX(), .8f, 2.f);
      prop->propagateToX(mR[tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector()], .8f, 2.f);
      prop->setTrack(trkWork);
#endif
 
      if (!wasTrackStored) {
#ifdef ENABLE_GPUTRDDEBUG
        mDebug->SetTrackParameterNoUp(trackNoUp, iLayer);
#endif
        mDebug->SetTrackParameter(*trkWork, iLayer);
        mDebug->SetRawTrackletPosition(spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getX(), spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getY(), spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getZ(), iLayer);
        mDebug->SetCorrectedTrackletPosition(trkltPosUp, iLayer);
        mDebug->SetTrackletCovariance(trkltCovUp, iLayer);
        mDebug->SetTrackletProperties(spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getDy(), tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector(), iLayer);
        wasTrackStored = true;
      }
 
      if (!prop->update(trkltPosUp, trkltCovUp)) {
        if (ENABLE_WARNING) {
          GPUWarning("Failed to update track %i with space point in layer %i", iTrk, iLayer);
        }
        trkWork->setChi2(trkWork->getChi2() + Param().rec.trd.penaltyChi2);
        if (trkWork->getIsFindable(iLayer)) {
          if (trkWork->getNmissingConsecLayers(iLayer) >= Param().rec.trd.stopTrkAfterNMissLy) {
            trkWork->setIsStopped();
          }
        }
        if (iUpdate == 0 && mNCandidates > 1) {
          *t = mCandidates[2 * iUpdate + nextIdx];
        }
        continue;
      }
      if (!trkWork->CheckNumericalQuality()) {
        if (ENABLE_INFO) {
          GPUInfo("Track %i has invalid covariance matrix. Aborting track following\n", iTrk);
        }
        return false;
      }
      trkWork->addTracklet(iLayer, mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId);
      trkWork->setChi2(mHypothesis[iUpdate + hypothesisIdxOffset].mChi2);
      trkWork->setIsFindable(iLayer);
      trkWork->setCollisionId(collisionId);
      // check if track crosses padrows
      float projZEntry, projYEntry;
      // Get Z for Entry of Track
      prop->getPropagatedYZ(trkWork->getX() - mGeo->GetCdrHght(), projYEntry, projZEntry);
      // Get Padrow number for Exit&Entry and compare. If is not equal mark
      // as padrow crossing. While simple, this comes with the cost of not
      // properly handling Tracklets that hit a different Pad but still get
      // attached to the Track. It is estimated that the probability for
      // this is low.
      const auto padrowEntry = pad->GetPadRowNumber(projZEntry);
      const auto padrowExit = pad->GetPadRowNumber(trkWork->getZ());
      if (padrowEntry != padrowExit) {
        trkWork->setIsCrossingNeighbor(iLayer);
        trkWork->setHasPadrowCrossing();
      }
 
      // Mark tracklets as Padrow crossing if they have a neighboring tracklet.
      for (int32_t trkltIdx = glbTrkltIdxOffset + mTrackletIndexArray[trkltIdxOffset + currDet]; trkltIdx < glbTrkltIdxOffset + mTrackletIndexArray[trkltIdxOffset + currDet + 1]; ++trkltIdx) {
        // skip orig tracklet
        if (mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId == trkltIdx) {
          continue;
        }
        if (CAMath::Abs(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetZbin() - tracklets[trkltIdx].GetZbin()) == 1 &&
            CAMath::Abs(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetY() - tracklets[trkltIdx].GetY()) < 1) {
          trkWork->setIsCrossingNeighbor(iLayer);
          trkWork->setHasNeighbor();
          break;
        }
      }
      if (iUpdate == 0 && mNCandidates > 1) {
        *t = mCandidates[2 * iUpdate + nextIdx];
      }
 
    } // end update loop
 
    if (!isOK) {
      if (ENABLE_INFO) {
        GPUInfo("Track %i cannot be followed. Stopped in layer %i", iTrk, iLayer);
      }
      return false;
    }
  } // end layer loop
 
  // --------------------------------------------------------------------------------
  // add some debug information (compare labels of attached tracklets to track label)
  // and store full track information
  // --------------------------------------------------------------------------------
  if (mDebugOutput) {
    mDebug->SetTrack(*t);
    mDebug->Output();
  }
  if (ENABLE_INFO) {
    GPUInfo("Ended track following for track %i at x=%f with pt=%f. Attached %i tracklets", t->getRefGlobalTrackIdRaw(), t->getX(), t->getPt(), t->getNtracklets());
  }
  if (nCurrHypothesis > 1) {
    if (CAMath::Abs(mHypothesis[hypothesisIdxOffset + 1].GetReducedChi2() - mHypothesis[hypothesisIdxOffset].GetReducedChi2()) < Param().rec.trd.chi2SeparationCut) {
      t->setIsAmbiguous();
    }
  }
  return true;
}
 
template <class TRDTRK, class PROP>
GPUd() void GPUTRDTracker_t<TRDTRK, PROP>::InsertHypothesis(Hypothesis hypo, int32_t& nCurrHypothesis, int32_t idxOffset)
{
  // Insert hypothesis into the array. If the array is full and the reduced chi2 is worse
  // than the worst hypothesis in the array it is dropped.
  // The hypothesis array is always sorted.
 
  if (nCurrHypothesis == 0) {
    // this is the first hypothesis in the array
    mHypothesis[idxOffset] = hypo;
    ++nCurrHypothesis;
  } else if (nCurrHypothesis > 0 && nCurrHypothesis < mNCandidates) {
    // insert the hypothesis into the right position and shift all worse hypothesis to the right
    for (int32_t i = idxOffset; i < nCurrHypothesis + idxOffset; ++i) {
      if (hypo.GetReducedChi2() < mHypothesis[i].GetReducedChi2()) {
        for (int32_t k = nCurrHypothesis + idxOffset; k > i; --k) {
          mHypothesis[k] = mHypothesis[k - 1];
        }
        mHypothesis[i] = hypo;
        ++nCurrHypothesis;
        return;
      }
    }
    mHypothesis[nCurrHypothesis + idxOffset] = hypo;
    ++nCurrHypothesis;
    return;
  } else {
    // array is already full, check if new hypothesis should be inserted
    int32_t i = nCurrHypothesis + idxOffset - 1;
    for (; i >= idxOffset; --i) {
      if (mHypothesis[i].GetReducedChi2() < hypo.GetReducedChi2()) {
        break;
      }
    }
    if (i < (nCurrHypothesis + idxOffset - 1)) {
      // new hypothesis should be inserted into the array
      for (int32_t k = nCurrHypothesis + idxOffset - 1; k > i + 1; --k) {
        mHypothesis[k] = mHypothesis[k - 1];
      }
      mHypothesis[i + 1] = hypo;
    }
  }
}
 
template <class TRDTRK, class PROP>
GPUd() int32_t GPUTRDTracker_t<TRDTRK, PROP>::GetDetectorNumber(const float zPos, const float alpha, const int32_t layer) const
{
  //--------------------------------------------------------------------
  // if track position is within chamber return the chamber number
  // otherwise return -1
  //--------------------------------------------------------------------
  int32_t stack = mGeo->GetStack(zPos, layer);
  if (stack < 0) {
    return -1;
  }
  int32_t sector = GetSector(alpha);
 
  return mGeo->GetDetector(layer, stack, sector);
}
 
template <class TRDTRK, class PROP>
GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::AdjustSector(PROP* prop, TRDTRK* t) const
{
  //--------------------------------------------------------------------
  // rotate track in new sector if necessary and
  // propagate to previous x afterwards
  // cancel if track crosses two sector boundaries
  //--------------------------------------------------------------------
  float alpha = mGeo->GetAlpha();
  float xTmp = t->getX();
  float y = t->getY();
  float yMax = t->getX() * CAMath::Tan(0.5f * alpha);
  float alphaCurr = t->getAlpha();
 
  if (CAMath::Abs(y) > 2.f * yMax) {
    if (ENABLE_INFO) {
      GPUInfo("AdjustSector: Track %i with pT = %f crossing two sector boundaries at x = %f", t->getRefGlobalTrackIdRaw(), t->getPt(), t->getX());
    }
    return false;
  }
 
  int32_t nTries = 0;
  while (CAMath::Abs(y) > yMax) {
    if (nTries >= 2) {
      return false;
    }
    int32_t sign = (y > 0) ? 1 : -1;
    float alphaNew = alphaCurr + alpha * sign;
    if (alphaNew > CAMath::Pi()) {
      alphaNew -= 2 * CAMath::Pi();
    } else if (alphaNew < -CAMath::Pi()) {
      alphaNew += 2 * CAMath::Pi();
    }
    if (!prop->rotate(alphaNew)) {
      return false;
    }
    if (!prop->propagateToX(xTmp, .8f, 2.f)) {
      return false;
    }
    y = t->getY();
    ++nTries;
  }
  return true;
}
 
template <class TRDTRK, class PROP>
GPUd() int32_t GPUTRDTracker_t<TRDTRK, PROP>::GetSector(float alpha) const
{
  //--------------------------------------------------------------------
  // TRD sector number for reference system alpha
  //--------------------------------------------------------------------
  if (alpha < 0) {
    alpha += 2.f * CAMath::Pi();
  } else if (alpha >= 2.f * CAMath::Pi()) {
    alpha -= 2.f * CAMath::Pi();
  }
  return (int32_t)(alpha * (float)kNSectors / (2.f * CAMath::Pi()));
}
 
template <class TRDTRK, class PROP>
GPUd() float GPUTRDTracker_t<TRDTRK, PROP>::GetAlphaOfSector(const int32_t sec) const
{
  //--------------------------------------------------------------------
  // rotation angle for TRD sector sec
  //--------------------------------------------------------------------
  float alpha = 2.0f * CAMath::Pi() / (float)kNSectors * ((float)sec + 0.5f);
  if (alpha > CAMath::Pi()) {
    alpha -= 2 * CAMath::Pi();
  }
  return alpha;
}
 
template <class TRDTRK, class PROP>
GPUd() void GPUTRDTracker_t<TRDTRK, PROP>::RecalcTrkltCov(const float tilt, const float snp, const float rowSize, const float pull, const int occupancy, float (&cov)[3])
{
  //--------------------------------------------------------------------
  // recalculate tracklet covariance taking track phi angle into account
  // store the new values in cov
  //--------------------------------------------------------------------
  float t2 = tilt * tilt;      // tan^2 (tilt)
  float c2 = 1.f / (1.f + t2); // cos^2 (tilt)
  float sy2 = mRecoParam->getRPhiRes(snp, CAMath::Abs(pull), occupancy);
  float sz2 = rowSize * rowSize / 12.f;
  cov[0] = c2 * (sy2 + t2 * sz2);
  cov[1] = c2 * tilt * (sz2 - sy2);
  cov[2] = c2 * (t2 * sy2 + sz2);
}
 
template <class TRDTRK, class PROP>
GPUd() void GPUTRDTracker_t<TRDTRK, PROP>::RecalcTrkltCovDy(const float tilt, const float snp, const float pull, const int occupancy, float (&cov)[6])
{
  float t2 = tilt * tilt;      // tan^2 (tilt)
  float c2 = 1.f / (1.f + t2); // cos^2 (tilt)
  // float sy2 = mRecoParam->getRPhiRes(snp, CAMath::Abs(pull), occupancy);
  float sdy2 = mRecoParam->getDyRes(snp, occupancy);
  cov[3] = mRecoParam->getCorrYDy(snp) * CAMath::Sqrt(sdy2 * c2);
  cov[4] = -tilt * mRecoParam->getCorrYDy(snp) * CAMath::Sqrt(sdy2 * c2);
  cov[5] = sdy2;
}
 
template <class TRDTRK, class PROP>
GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::InvertCov(float (&cov)[6])
{
  // invert a 3*3 symmetric matrix. Adapted from https://root.cern.ch/doc/master/TMatrixTSymCramerInv_8cxx_source.html
 
  float c00 = cov[2] * cov[5] - cov[4] * cov[4];
  float c01 = cov[4] * cov[3] - cov[1] * cov[5];
  float c02 = cov[1] * cov[4] - cov[2] * cov[3];
  float c11 = cov[5] * cov[0] - cov[3] * cov[3];
  float c12 = cov[3] * cov[1] - cov[4] * cov[0];
  float c22 = cov[0] * cov[2] - cov[1] * cov[1];
 
  float t0 = CAMath::Abs(cov[0]);
  float t1 = CAMath::Abs(cov[1]);
  float t2 = CAMath::Abs(cov[3]);
 
  float det;
  float tmp;
 
  if (t0 >= t1) {
    if (t2 >= t0) {
      tmp = cov[3];
      det = c12 * c01 - c11 * c02;
    } else {
      tmp = cov[0];
      det = c11 * c22 - c12 * c12;
    }
  } else if (t2 >= t1) {
    tmp = cov[3];
    det = c12 * c01 - c11 * c02;
  } else {
    tmp = cov[1];
    det = c02 * c12 - c01 * c22;
  }
 
  if (det == 0 || tmp == 0) {
    return false;
  }
 
  float s = tmp / det;
 
  cov[0] = s * c00;
  cov[1] = s * c01;
  cov[3] = s * c02;
  cov[2] = s * c11;
  cov[4] = s * c12;
  cov[5] = s * c22;
 
  return true;
}
 
template <class TRDTRK, class PROP>
GPUd() float GPUTRDTracker_t<TRDTRK, PROP>::GetAngularPull(float dYtracklet, float snp, int occupancy) const
{
  float dYtrack = mRecoParam->convertAngleToDy(snp);
  float dYresolution = mRecoParam->getDyRes(snp, occupancy);
  if (dYresolution < 1e-6f) {
    return 999.f;
  }
  return (dYtracklet - dYtrack) / CAMath::Sqrt(dYresolution);
}
 
template <class TRDTRK, class PROP>
GPUd() int GPUTRDTracker_t<TRDTRK, PROP>::GetNtrackletsChamber(int collisionId, int detector) const
{
  // get the number of tracklets for a given chamber and trigger
  int32_t trkltIdxOffset = collisionId * (kNChambers + 1);
  int nTrackletsChamber = mTrackletIndexArray[trkltIdxOffset + detector + 1] - mTrackletIndexArray[trkltIdxOffset + detector];
  return nTrackletsChamber;
}
 
template <class TRDTRK, class PROP>
GPUd() void GPUTRDTracker_t<TRDTRK, PROP>::FindChambersInRoad(const TRDTRK* t, const float roadY, const float roadZ, const int32_t iLayer, int32_t* det, const float zMax, const float alpha, const float zShiftTrk) const
{
  //--------------------------------------------------------------------
  // determine initial chamber where the track ends up
  // add more chambers of the same sector or (and) neighbouring
  // stack if track is close the edge(s) of the chamber
  //--------------------------------------------------------------------
 
  const float yMax = CAMath::Abs(mGeo->GetCol0(iLayer));
  float zTrk = t->getZ() + zShiftTrk;
 
  int32_t currStack = mGeo->GetStack(zTrk, iLayer);
  int32_t currSec = GetSector(alpha);
  int32_t currDet;
 
  int32_t nDets = 0;
 
  if (currStack > -1) {
    // chamber unambiguous
    currDet = mGeo->GetDetector(iLayer, currStack, currSec);
    det[nDets++] = currDet;
    const GPUTRDpadPlane* pp = mGeo->GetPadPlane(iLayer, currStack);
    int32_t lastPadRow = mGeo->GetRowMax(iLayer, currStack, 0);
    float zCenter = pp->GetRowPos(lastPadRow / 2);
    if ((zTrk + roadZ) > pp->GetRow0() || (zTrk - roadZ) < pp->GetRowEnd()) {
      int32_t addStack = zTrk > zCenter ? currStack - 1 : currStack + 1;
      if (addStack < kNStacks && addStack > -1) {
        det[nDets++] = mGeo->GetDetector(iLayer, addStack, currSec);
      }
    }
  } else {
    if (CAMath::Abs(zTrk) > zMax) {
      // shift track in z so it is in the TRD acceptance
      if (zTrk > 0) {
        currDet = mGeo->GetDetector(iLayer, 0, currSec);
      } else {
        currDet = mGeo->GetDetector(iLayer, kNStacks - 1, currSec);
      }
      det[nDets++] = currDet;
      currStack = mGeo->GetStack(currDet);
    } else {
      // track in between two stacks, add both surrounding chambers
      // gap between two stacks is 4 cm wide
      currDet = GetDetectorNumber(zTrk + 4.0f, alpha, iLayer);
      if (currDet != -1) {
        det[nDets++] = currDet;
      }
      currDet = GetDetectorNumber(zTrk - 4.0f, alpha, iLayer);
      if (currDet != -1) {
        det[nDets++] = currDet;
      }
    }
  }
  // add chamber(s) from neighbouring sector in case the track is close to the boundary
  if ((CAMath::Abs(t->getY()) + roadY) > yMax) {
    const int32_t nStacksToSearch = nDets;
    int32_t newSec;
    if (t->getY() > 0) {
      newSec = (currSec + 1) % kNSectors;
    } else {
      newSec = (currSec > 0) ? currSec - 1 : kNSectors - 1;
    }
    for (int32_t idx = 0; idx < nStacksToSearch; ++idx) {
      currStack = mGeo->GetStack(det[idx]);
      det[nDets++] = mGeo->GetDetector(iLayer, currStack, newSec);
    }
  }
  // skip PHOS hole and non-existing chamber 17_4_4
  for (int32_t iDet = 0; iDet < nDets; iDet++) {
    if (!mGeo->ChamberInGeometry(det[iDet])) {
      det[iDet] = -1;
    }
  }
}
 
template <class TRDTRK, class PROP>
GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::IsGeoFindable(const TRDTRK* t, const int32_t layer, const float alpha, const float zShiftTrk) const
{
  //--------------------------------------------------------------------
  // returns true if track position inside active area of the TRD
  // and not too close to the boundaries
  //--------------------------------------------------------------------
 
  float zTrk = t->getZ() + zShiftTrk;
 
  int32_t det = GetDetectorNumber(zTrk, alpha, layer);
 
  // reject tracks between stacks
  if (det < 0) {
    return false;
  }
 
  // reject tracks in PHOS hole and for non existent chamber 17_4_4
  if (!mGeo->ChamberInGeometry(det)) {
    return false;
  }
 
  // reject tracks in bad chambers
  if (mChamberStatus[det]) {
    return false;
  }
 
  const GPUTRDpadPlane* pp = mGeo->GetPadPlane(det);
  float yMax = pp->GetColEnd();
  float zMax = pp->GetRow0();
  float zMin = pp->GetRowEnd();
 
  float epsY = 5.f;
  float epsZ = 5.f;
 
  // reject tracks closer than epsY cm to pad plane boundary
  if (yMax - CAMath::Abs(t->getY()) < epsY) {
    return false;
  }
  // reject tracks closer than epsZ cm to stack boundary
  if (!((zTrk > zMin + epsZ) && (zTrk < zMax - epsZ))) {
    return false;
  }
 
  // reject tracks if the corresponding pad and neighboring pads (to take into account track uncertainty and charge sharing with neighbours) are masked
  /*int32_t padrow = pp->GetPadRowNumber(zTrk);
  int32_t padcol = pp->GetPadColNumber(t->getY());
  int32_t idxPad = det * kNPadColumns * kNPadRows + padcol * kNPadRows + padrow;
  if (padrow != -1 && padcol != -1 && mPadStatus[idxPad] && (padrow == 0 || mPadStatus[idxPad - kNPadRows]) && (padrow == kNPadRows - 1 || mPadStatus[idxPad + kNPadRows])) {
    return false;
  }*/
 
  return true;
}
 
#ifndef GPUCA_GPUCODE
namespace o2::gpu
{
template class GPUTRDTracker_t<GPUTRDTrack, GPUTRDPropagator>;
template class GPUTRDTracker_t<GPUTRDTrackGPU, GPUTRDPropagatorGPU>;
} // namespace o2::gpu
#endif