TrackInterpolation.h

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.
 
 
#ifndef ALICEO2_TPC_TRACKINTERPOLATION_H_
#define ALICEO2_TPC_TRACKINTERPOLATION_H_
 
#include <gsl/span>
#include "CommonDataFormat/EvIndex.h"
#include "CommonDataFormat/RangeReference.h"
#include "ReconstructionDataFormats/Track.h"
#include "ReconstructionDataFormats/TrackTPCITS.h"
#include "ReconstructionDataFormats/MatchInfoTOF.h"
#include "ReconstructionDataFormats/GlobalTrackID.h"
#include "DataFormatsITSMFT/Cluster.h"
#include "DataFormatsITSMFT/TrkClusRef.h"
#include "DataFormatsITSMFT/TopologyDictionary.h"
#include "DataFormatsITS/TrackITS.h"
#include "DataFormatsTPC/TrackTPC.h"
#include "DataFormatsTPC/Constants.h"
#include "DataFormatsTOF/Cluster.h"
#include "DataFormatsTRD/TrackTRD.h"
#include "DataFormatsTRD/Tracklet64.h"
#include "DataFormatsTRD/CalibratedTracklet.h"
#include "SpacePoints/SpacePointsCalibParam.h"
#include "SpacePoints/SpacePointsCalibConfParam.h"
#include "TPCReconstruction/TPCFastTransformHelperO2.h"
#include "DetectorsBase/Propagator.h"
#include "DataFormatsGlobalTracking/RecoContainer.h"
#include "GPUTRDRecoParam.h"
#include "TRDBase/Geometry.h"
#include "CommonUtils/TreeStreamRedirector.h"
 
class TTree;
 
namespace o2
{
 
namespace gpu
{
class GPUParam;
}
 
namespace tpc
{
class VDriftCorrFact;
 
struct TPCClusterResiduals {
  float dy{};           
  float dz{};           
  float y{};            
  float z{};            
  float snp{};          
  unsigned char sec{};  
  unsigned char dRow{}; 
  unsigned char flags{}; 
  ClassDefNV(TPCClusterResiduals, 6);
};
 
struct UnbinnedResid {
  UnbinnedResid() = default;
  UnbinnedResid(float dyIn, float dzIn, float tgSlpIn, float yIn, float zIn,
                unsigned char rowIn, unsigned char secIn, short chanIn = -1, bool rejFlag = false) : dy(static_cast<short>(dyIn * 0x7fff / param::MaxResid)),
                                                                                                     dz(static_cast<short>(dzIn * 0x7fff / param::MaxResid)),
                                                                                                     tgSlp(static_cast<short>(tgSlpIn * 0x7fff / param::MaxTgSlp)),
                                                                                                     y(static_cast<short>(yIn * 0x7fff / param::MaxY)),
                                                                                                     z(static_cast<short>(zIn * 0x7fff / param::MaxZ)),
                                                                                                     row(rowIn),
                                                                                                     sec(secIn),
                                                                                                     channel(chanIn),
                                                                                                     rejected(rejFlag) {}
  short dy{0};          
  short dz{0};          
  short tgSlp{0};       
  short y{0};           
  short z{0};           
  unsigned char row{0}; 
  unsigned char sec{0}; 
  short channel{-1};    
  bool rejected{false}; 
 
  bool isTPC() const { return row < constants::MAXGLOBALPADROW; }
  bool isTRD() const { return row >= 160 && row < 166; }
  bool isTOF() const { return row == 170; }
  bool isITS() const { return row >= 180; }
  int getDetID() const { return isTPC() ? 1 : (isITS() ? 0 : (isTRD() ? 2 : (isTOF() ? 3 : -1))); }
  int getITSLayer() const { return row - 180; }
  int getTRDLayer() const { return row - 160; }
  float getAlpha() const;
  float getX() const;
 
  static void init(long timestamp = -1);
  static void checkInitDone();
  static bool gInitDone;
 
  ClassDefNV(UnbinnedResid, 3);
};
 
struct DetInfoResid { // detector info associated with residual
  uint32_t word = 0;  // container interpreted in a different way depending on the detector type
  //
  // TPC view: qTot and qMax of the cluster
  uint16_t qTotTPC() const { return static_cast<uint16_t>(word & 0xFFFFu); }
  uint16_t qMaxTPC() const { return static_cast<uint16_t>((word >> 16) & 0xFFFFu); }
  void setTPC(uint16_t qTot, uint16_t qMax) { word = (static_cast<uint32_t>(qMax) << 16) | static_cast<uint32_t>(qTot); }
  //
  // TRD view: q0, q1, q2 + calibrated slope (truncated to in +-3.5 range)
  static constexpr uint32_t TRDQ0NB = 7, TRDQ1NB = 7, TRDQ2NB = 6, TRDSlpNB = 12;
  static constexpr uint32_t TRDQ0Msk = (1 << TRDQ0NB) - 1, TRDQ1Msk = (1 << TRDQ1NB) - 1, TRDQ2Msk = ((1 << TRDQ2NB) - 1), TRDSlpMsk = (1 << TRDSlpNB) - 1;
  static constexpr float TRDMaxSlope = 3.5, TRDSlope2Int = ((1 << TRDSlpNB) - 1) / (2 * TRDMaxSlope), TRDInt2Slope = 1.f / TRDSlope2Int;
  uint16_t q0TRD() const { return static_cast<uint16_t>(word & TRDQ0Msk); }
  uint16_t q1TRD() const { return static_cast<uint16_t>((word >> TRDQ0NB) & TRDQ1Msk); }
  uint16_t q2TRD() const { return static_cast<uint16_t>((word >> (TRDQ0NB + TRDQ1NB)) & TRDQ2Msk); }
  float slopeTRD() const { return ((word >> (TRDQ0NB + TRDQ1NB + TRDQ2NB)) & TRDSlpMsk) * TRDInt2Slope - TRDMaxSlope; }
  void setTRD(uint8_t q0, uint8_t q1, uint8_t q2, float slope)
  {
    float rslope = (slope + TRDMaxSlope) * TRDSlope2Int;
    if (rslope < 0.f) {
      rslope = 0;
    } else if (rslope > TRDSlpMsk) {
      rslope = TRDSlpMsk;
    }
    uint32_t slpI = std::round(rslope);
    word = (static_cast<uint32_t>(slpI << (TRDQ0NB + TRDQ1NB + TRDQ2NB)) |
            static_cast<uint32_t>((q2 & TRDQ2Msk) << (TRDQ0NB + TRDQ1NB)) |
            static_cast<uint32_t>((q1 & TRDQ1Msk) << TRDQ0NB) |
            static_cast<uint32_t>(q0 & TRDQ0Msk));
  }
  //
  // TOF view (time difference in \mus wrt seeding ITS-TPC track)
  float timeTOF() const { return std::bit_cast<float>(word); }
  void setTOF(float t) { word = std::bit_cast<uint32_t>(t); }
  //
  // No info for ITS is stored
  //
  // PV view (time difference in \mus wrt contributing ITS-TPC track)
  float timePV() const { return std::bit_cast<float>(word); }
  void setPV(float t) { word = std::bit_cast<uint32_t>(t); }
 
  ClassDefNV(DetInfoResid, 1);
};
 
struct TrackDataCompact {
  TrackDataCompact() = default;
  TrackDataCompact(uint32_t idx, std::array<uint8_t, 4> mlt, uint8_t nRes, uint8_t source, uint8_t nextraRes = 0, int8_t filt = -1) : idxFirstResidual(idx), multStack{mlt}, nResiduals(nRes), sourceId(source), filterFlag(filt), nExtDetResid(nextraRes) {}
  uint32_t idxFirstResidual; 
  std::array<uint8_t, 4> multStack{}; // multiplicity in the stack packed as asinh(x*0.05)/0.05
  uint8_t nResiduals;        
  uint8_t nExtDetResid = 0;  
  uint8_t sourceId;          
  int8_t filterFlag = -1;    
 
  void setMultStack(float v, int stack)
  {
    uint32_t mltPacked = std::round(std::asinh(v * 0.05) / 0.05);
    multStack[stack] = mltPacked < 0xff ? mltPacked : 0xff;
  }
  float getMultStack(int stack) const
  {
    return std::sinh(multStack[stack] * 0.05) / 0.05;
  }
  float getMultStackPacked(int stack) const { return multStack[stack]; }
 
  ClassDefNV(TrackDataCompact, 4);
};
 
// TODO add to UnbinnedResid::sec flag if cluster was used or not
// TODO add to TrackDataCompact::sourceId flag if track was used or not (if possible)
 
struct TrackDataExtended {
  o2::dataformats::GlobalTrackID gid{};              
  o2::its::TrackITS trkITS{};                        
  o2::tpc::TrackTPC trkTPC{};                        
  o2::trd::TrackTRD trkTRD{};                        
  o2::track::TrackPar trkOuter{};                    
  o2::dataformats::MatchInfoTOF matchTOF{};          
  std::vector<o2::BaseCluster<float>> clsITS{};      
  std::vector<o2::trd::Tracklet64> trkltTRD{};       
  std::vector<o2::trd::CalibratedTracklet> clsTRD{}; 
  o2::tof::Cluster clsTOF{};                         
  o2::dataformats::RangeReference<> clIdx{};         
  uint8_t nExtDetResid = 0;                          
  int8_t filterFlag = -1;                            
  ClassDefNV(TrackDataExtended, 4);
};
 
struct TrackData {
  o2::dataformats::GlobalTrackID gid{};      
  o2::track::TrackParCov par{};              
  float dEdxTPC{};                           
  float chi2TPC{};                           
  float chi2ITS{};                           
  float chi2TRD{};                           
  float deltaTOF{};                          
 
  unsigned short nClsTPC{};                  
  unsigned short nClsITS{};                  
  unsigned short nTrkltsTRD{};               
  unsigned short clAvailTOF{};               
  short TRDTrkltSlope[6] = {};               
  uint8_t nExtDetResid = 0;                  
  int8_t filterFlag = -1;                    
  o2::dataformats::RangeReference<> clIdx{}; 
  std::array<uint8_t, 4> multStack{};        // multiplicity in the stack packed as asinh(x*0.05)/0.05
  float getT0Error() const { return float(clAvailTOF); }
  bool isTOFAvail() const { return clAvailTOF != 0; }
 
  void setMultStack(float v, int stack)
  {
    uint32_t mltPacked = std::round(std::asinh(v * 0.05) / 0.05);
    multStack[stack] = mltPacked < 0xff ? mltPacked : 0xff;
  }
  float getMultStack(int stack) const
  {
    return std::sinh(multStack[stack] * 0.05) / 0.05;
  }
  float getMultStackPacked(int stack) const { return multStack[stack]; }
 
  ClassDefNV(TrackData, 12);
};
 
class TrackInterpolation
{
 public:
  using MatCorrType = o2::base::Propagator::MatCorrType;
 
  TrackInterpolation() = default;
 
  // since this class has pointer members, we should explicitly delete copy and assignment operators
  TrackInterpolation(const TrackInterpolation&) = delete;
  TrackInterpolation& operator=(const TrackInterpolation&) = delete;
 
  ~TrackInterpolation();
 
  enum {
    ExtOut = 0, 
    ExtIn,      
    Int,        
    NIndices    
  };
 
  struct CacheStruct {
    std::array<float, NIndices> y{};
    std::array<float, NIndices> z{};
    std::array<float, NIndices> sy2{};
    std::array<float, NIndices> szy{};
    std::array<float, NIndices> sz2{};
    std::array<float, NIndices> snp{};
    float clY{0.f};
    float clZ{0.f};
    float clAngle{0.f};
    unsigned short clAvailable{0};
    unsigned char clSec{0};
    unsigned char clFlags{0};
  };
 
  struct ValidationPoint {
    float xTrk{0.f};
    float yTrk{0.f};
    float zTrk{0.f};
    float xLab{0.f};
    float yLab{0.f};
    float sPath{0.f};
    float dy{0.f};
    float dz{0.f};
    float tglArr{0.f};
    float residHelixY{0.f};
    float residHelixZ{0.f};
    float diffYSmooth{0.f};
    float diffZSmooth{0.f};
    int8_t sec{0};
    bool flagRej{false};
    ClassDefNV(ValidationPoint, 1);
  };
 
  struct TrackValidationData {
    float qpt{0.f};
    float tgl{0.f};
    float xcLab{0.f};
    float ycLab{0.f};
    float r{0.f};
    float zOffs{0.f};
    uint8_t nRej = 0;
    std::vector<ValidationPoint> points;
    void clear()
    {
      points.clear();
      nRej = 0;
    }
    ClassDefNV(TrackValidationData, 1);
  };
 
  // -------------------------------------- processing functions --------------------------------------------------
 
  void init(o2::dataformats::GlobalTrackID::mask_t src, o2::dataformats::GlobalTrackID::mask_t srcMap);
 
  bool isInputTrackAccepted(const o2::dataformats::GlobalTrackID& gid, const o2::globaltracking::RecoContainer::GlobalIDSet& gidTable, const o2::dataformats::PrimaryVertex& pv) const;
 
  o2::dataformats::GlobalTrackID::Source findValidSource(const o2::dataformats::GlobalTrackID::mask_t mask, const o2::dataformats::GlobalTrackID::Source src) const;
 
  void prepareInputTrackSample(const o2::globaltracking::RecoContainer& inp);
 
  bool isTrackSelected(const o2::track::TrackParCov& trk) const;
 
  void process();
 
  void extrapolateTrack(int iSeed);
 
  void interpolateTrack(int iSeed);
 
  void reset();
 
  // refit ITS track taking PID (unless already refitted) from the seed and reassign to the seed
  bool refITSTrack(o2::dataformats::GlobalTrackID, int iSeed);
  // -------------------------------------- outlier rejection --------------------------------------------------
 
  int8_t validateTrack(const TrackData& trk, TrackValidationData& params, const std::vector<TPCClusterResiduals>& clsRes, bool interpol);
 
  bool outlierFiltering(const TrackData& trk, TrackValidationData& params, const std::vector<TPCClusterResiduals>& clsRes);
 
  float checkResiduals(const TrackData& trk, TrackValidationData& params, const std::vector<TPCClusterResiduals>& clsRes);
 
  bool compareToHelix(const TrackData& trk, TrackValidationData& params, const std::vector<TPCClusterResiduals>& clsRes);
 
  void diffToLocLine(TrackValidationData& params, int start, int np);
 
  void diffToMA(const int np, const std::array<float, param::NPadRows>& y, std::array<float, param::NPadRows>& diffMA);
 
  // -------------------------------------- settings --------------------------------------------------
  void setNHBPerTF(int n) { mNHBPerTF = n; }
 
  void setTPCVDrift(const o2::tpc::VDriftCorrFact& v);
 
  void setMatCorr(MatCorrType matCorr) { mMatCorr = matCorr; }
 
  void setMaxTracksPerTF(int n) { mMaxTracksPerTF = n; }
 
  void setAddTracksForMapPerTF(int n) { mAddTracksForMapPerTF = n; }
 
  void setDumpTrackPoints() { mDumpTrackPoints = true; }
 
  void setITSClusterDictionary(const o2::itsmft::TopologyDictionary* dict) { mITSDict = dict; }
 
  void setProcessSeeds() { mProcessSeeds = true; }
 
  void setProcessITSTPConly() { mProcessITSTPConly = true; }
 
  void setSqrtS(float s) { mSqrtS = s; }
 
  void setExtDetResid(bool v) { mExtDetResid = v; }
 
  int processTRDLayer(const o2::trd::TrackTRD& trkTRD, int iLayer, o2::track::TrackParCov& trkWork, std::array<float, 2>* trkltTRDYZ = nullptr,
                      std::array<float, 3>* trkltTRDCov = nullptr, TrackData* trkData = nullptr,
                      o2::trd::Tracklet64* trk64 = nullptr, o2::trd::CalibratedTracklet* trkCalib = nullptr);
 
  // --------------------------------- output ---------------------------------------------
  std::vector<UnbinnedResid>& getClusterResiduals() { return mClRes; }
  std::vector<DetInfoResid>& getClusterResidualsDetInfo() { return mDetInfoRes; }
  std::vector<TrackDataCompact>& getTrackDataCompact() { return mTrackDataCompact; }
  std::vector<TrackDataExtended>& getTrackDataExtended() { return mTrackDataExtended; }
  std::vector<TrackData>& getReferenceTracks() { return mTrackData; }
 
  void setLane(int lID, int nL)
  {
    mLaneID = lID;
    mNLanes = nL;
  }
 
  void finalize();
 
 private:
  static constexpr float sFloatEps{1.e-7f}; 
  static constexpr int NSTACKS = 4;
  static constexpr std::array<int, NSTACKS + 1> STACKROWS{0, 63, 97, 127, 152};
  // parameters + settings
  const SpacePointsCalibConfParam* mParams = nullptr;
  std::shared_ptr<o2::gpu::GPUParam> mTPCParam = nullptr;
  int mLaneID = 0;
  int mNLanes = 1;
  int mNHBPerTF = 32;
  int mNTPCOccBinLength = 16;                                   
  float mNTPCOccBinLengthInv = 1.f / 16;                        
  float mTPCTimeBinMUS{.2f};                                    
  float mTPCVDriftRef = -1.;                                    
  float mTPCDriftTimeOffsetRef = 0.;                            
  float mSqrtS{13600.f};                                        
  MatCorrType mMatCorr{MatCorrType::USEMatCorrNONE};            
  int mMaxTracksPerTF{-1};                                      
  int mAddTracksForMapPerTF{0};                                 
  bool mDumpTrackPoints{false};                                 
  bool mExtDetResid{true};                                      
  bool mProcessSeeds{false};                                    
  bool mProcessITSTPConly{false};                               
  o2::dataformats::GlobalTrackID::mask_t mSourcesConfigured;    
  o2::dataformats::GlobalTrackID::mask_t mSourcesConfiguredMap; 
  bool mSingleSourcesConfigured{true};                          
 
  // input
  const o2::globaltracking::RecoContainer* mRecoCont = nullptr;             
  std::vector<o2::dataformats::GlobalTrackID> mGIDs{};                      
  std::vector<o2::globaltracking::RecoContainer::GlobalIDSet> mGIDtables{}; 
  std::vector<float> mTrackTimes{};                                         
  std::vector<int> mTrackPVID{};                                            
  std::vector<o2::track::TrackParCov> mSeeds{};                             
  std::vector<int> mParentID{};                                             
  std::map<int, int> mTrackTypes;                                           
  std::array<std::vector<uint32_t>, 4> mTrackIndices;                       
  gsl::span<const TPCClRefElem> mTPCTrackClusIdx;                           
  gsl::span<const unsigned char> mTPCShClassMap;                            
  const ClusterNativeAccess* mTPCClusterIdxStruct = nullptr;                
 
  // ITS specific input only needed for debugging
  gsl::span<const int> mITSTrackClusIdx;                    
  std::vector<o2::BaseCluster<float>> mITSClustersArray;    
  std::vector<int> mITSRefitSeedID;                         
  const o2::itsmft::TopologyDictionary* mITSDict = nullptr; 
 
  // output
  std::vector<TrackData> mTrackData{};                 
  std::vector<TrackDataCompact> mTrackDataCompact{};   
  std::vector<TrackDataExtended> mTrackDataExtended{}; 
  std::vector<UnbinnedResid> mClRes{};                 
  std::vector<DetInfoResid> mDetInfoRes{};             
  std::unique_ptr<o2::utils::TreeStreamRedirector> mDBGOut;
 
  // cache
  std::array<CacheStruct, constants::MAXGLOBALPADROW> mCache{{}}; 
  std::vector<o2::dataformats::GlobalTrackID> mGIDsSuccess;       
 
  TrackValidationData mTrackValidation;
 
  // helpers
  o2::gpu::GPUTRDRecoParam mRecoParam;                
  o2::trd::Geometry* mGeoTRD;                         
  std::unique_ptr<TPCFastTransform> mFastTransform{}; 
  float mBz;                                          
  bool mInitDone{false};                              
  size_t mRejectedResiduals{};                        
  size_t mNRejRefit = 0;
  size_t mNRejProp = 0;
  size_t mNRejLoop = 0;
 
  ClassDefNV(TrackInterpolation, 1);
};
 
} // namespace tpc
 
} // namespace o2
 
#endif