TrackResiduals.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_TRACKRESIDUALS_H_
#define ALICEO2_TPC_TRACKRESIDUALS_H_
 
#include <memory>
#include <vector>
#include <array>
#include <bitset>
#include <string>
 
#include "DataFormatsTPC/Defs.h"
#include "SpacePoints/SpacePointsCalibParam.h"
#include "SpacePoints/SpacePointsCalibConfParam.h"
#include "SpacePoints/TrackInterpolation.h"
 
#include "TTree.h"
#include "TFile.h"
 
namespace o2
{
namespace tpc
{
 
class TrackResiduals
{
 public:
  TrackResiduals() = default;
 
  TrackResiduals(const TrackResiduals&) = delete;
  TrackResiduals& operator=(const TrackResiduals&) = delete;
 
  enum { VoxZ,          
         VoxF,          
         VoxX,          
         VoxV,          
         VoxDim = 3,    
         VoxHDim = 4 }; 
 
  enum { ResX,     
         ResY,     
         ResZ,     
         ResD,     
         ResDim }; 
 
  enum { DistDone = 1 << 0,   
         DispDone = 1 << 1,   
         SmoothDone = 1 << 2, 
         Masked = 1 << 7 };   
 
  enum class KernelType { Epanechnikov,
                          Gaussian };
 
  struct VoxRes {
    VoxRes() = default;
    std::array<float, ResDim> D{};            
    std::array<float, ResDim> E{};            
    std::array<float, ResDim> DS{};           
    std::array<float, ResDim> DC{};           
    float EXYCorr{0.f};                       
    float dYSigMAD{0.f};                      
    float dZSigLTM{0.f};                      
    std::array<float, VoxHDim> stat{};        
    std::array<unsigned char, VoxDim> bvox{}; 
    unsigned char bsec{0};                    
    unsigned char flags{0};                   
    ClassDefNV(VoxRes, 1);
  };
 
  struct LocalResid {
    LocalResid() = default;
    LocalResid(short dyIn, short dzIn, short tgSlpIn, std::array<unsigned char, VoxDim> bvoxIn) : dy(dyIn), dz(dzIn), tgSlp(tgSlpIn), bvox(bvoxIn) {}
    short dy{0};                              
    short dz{0};                              
    short tgSlp{0};                           
    std::array<unsigned char, VoxDim> bvox{}; 
    ClassDefNV(LocalResid, 1);
  };
 
  struct VoxStats {
    VoxStats() = default;
    std::array<float, VoxDim> meanPos{};
    float nEntries{0.f};
    ClassDefNV(VoxStats, 1);
  };
 
  // -------------------------------------- initialization --------------------------------------------------
  void init(bool doBinning = true);
  void initBinning();
  void initResultsContainer(int iSec);
  void initVoxelStats();
  void reset();
 
  // -------------------------------------- settings --------------------------------------------------
  void setPrintMemoryUsage() { mPrintMem = true; }
  void setKernelType(KernelType kernel = KernelType::Epanechnikov, float bwX = 2.1f, float bwP = 2.1f, float bwZ = 1.7f, float scX = 1.f, float scP = 1.f, float scZ = 1.f);
 
  void setSmoothPol2(int dim, bool flag) { mSmoothPol2[dim] = flag; }
 
  void setVdriftCorr(float corr) { mEffVdriftCorr = corr; }
 
  void setT0Corr(float corr) { mEffT0Corr = corr; }
 
  // -------------------------------------- I/O --------------------------------------------------
 
  std::vector<LocalResid>& getLocalResVec() { return mLocalResidualsIn; }
  std::vector<VoxStats>** getVoxStatPtr() { return &mVoxStatsInPtr; }
 
  const std::array<std::vector<VoxRes>, SECTORSPERSIDE * SIDES>& getVoxelResults() const { return mVoxelResults; }
 
  // -------------------------------------- steering functions --------------------------------------------------
 
  void processSectorResiduals(Int_t iSec);
 
  void processVoxelResiduals(std::vector<float>& dy, std::vector<float>& dz, std::vector<float>& tg, VoxRes& resVox);
 
  void processVoxelDispersions(std::vector<float>& tg, std::vector<float>& dy, VoxRes& resVox);
 
  int validateVoxels(int iSec);
 
  void smooth(int iSec);
 
  // -------------------------------------- statistics --------------------------------------------------
 
  float fitPoly1Robust(std::vector<float>& x, std::vector<float>& y, std::array<float, 2>& res, std::array<float, 3>& err, float cutLTM) const;
 
  float getMAD2Sigma(std::vector<float> data) const;
 
  void medFit(int nPoints, int offset, const std::vector<float>& x, const std::vector<float>& y, float& a, float& b, std::array<float, 3>& err) const;
 
  float roFunc(int nPoints, int offset, const std::vector<float>& x, const std::vector<float>& y, float b, float& aa) const;
 
  float selectKthMin(const int k, std::vector<float>& data) const;
 
  bool getSmoothEstimate(int iSec, float x, float p, float z, std::array<float, ResDim>& res, int whichDim = 0);
 
  double getKernelWeight(std::array<double, 3> u2vec) const;
 
  static void fitCircle(int nCl, std::array<float, param::NPadRows>& x, std::array<float, param::NPadRows>& y, float& xc, float& yc, float& r, std::array<float, param::NPadRows>& residHelixY);
 
  static bool fitPoly1(int nCl, std::array<float, param::NPadRows>& x, std::array<float, param::NPadRows>& y, std::array<float, 2>& res);
 
  // -------------------------------------- binning / geometry --------------------------------------------------
 
  void setNXBins(int nBins) { mNXBins = nBins; }
  int getNXBins() const { return mNXBins; }
 
  void setNY2XBins(int nBins) { mNY2XBins = nBins; }
  int getNY2XBins() const { return mNY2XBins; }
 
  void setNZ2XBins(int nBins) { mNZ2XBins = nBins; }
  int getNZ2XBins() const { return mNZ2XBins; }
 
  int getNVoxelsPerSector() const { return mNVoxPerSector; }
 
  void setY2XBinning(const std::vector<float>& binning);
 
  void setZ2XBinning(const std::vector<float>& binning);
 
  size_t getGlbVoxBin(int ix, int ip, int iz) const;
 
  size_t getGlbVoxBin(const std::array<unsigned char, VoxDim>& bvox) const;
 
  void getVoxelCoordinates(int isec, int ix, int ip, int iz, float& x, float& p, float& z) const;
 
  float getX(int ix) const;
 
  float getMaxY2X(int ix) const;
 
  float getY2X(int ix, int ip) const;
 
  float getZ2X(int iz) const;
 
  bool getXBinIgnored(int iSec, int bin) const { return mXBinsIgnore[iSec].test(bin); }
 
  void findVoxel(float x, float y2x, float z2x, int& ix, int& ip, int& iz) const;
 
  bool findVoxelBin(int secID, float x, float y, float z, std::array<unsigned char, VoxDim>& bvox) const;
 
  int getXBin(float x) const;
 
  int getXBinExact(float x) const;
 
  int getRowID(float x) const;
 
  int getY2XBinExact(float y2x, int ix) const;
 
  int getY2XBin(float y2x, int ix) const;
 
  int getZ2XBinExact(float z2x) const;
 
  int getZ2XBin(float z2x) const;
 
  float getDXI(int ix) const;
 
  float getDY2XI(int ix, int iy = 0) const;
 
  float getDZ2X(int iz = 0) const;
 
  float getDZ2XI(int iz = 0) const;
 
  // -------------------------------------- debugging --------------------------------------------------
 
  void printMem() const;
 
  void dumpResults(int iSec);
 
  void createOutputFile(const char* filename = "debugVoxRes.root");
 
  void closeOutputFile();
 
  TFile* getOutputFilePtr() { return mFileOut.get(); }
 
  void setStats(const std::vector<TrackResiduals::VoxStats>& statsIn, int iSec);
 
  void fillStats(int iSec);
 
  void clear();
 
  TTree* getOutputTree() { return mTreeOut.get(); }
 
  void setAdhocScalingFactorX(const std::array<float, 2>& scaling) { mAdhocScalingX = scaling; }
 
  void doAdhocCorrectionZ2X(bool corr) { mDoAdhocCorrectionZ2X = corr; }
 
 private:
  std::bitset<SECTORSPERSIDE * SIDES> mInitResultsContainer{};
 
  // some constants
  static constexpr float sFloatEps{1.e-7f}; 
  static constexpr float sDeadZone{1.5f};   
  static constexpr int sSmtLinDim{4};       
  static constexpr int sMaxSmtDim{7};       
 
  // settings
  const SpacePointsCalibConfParam* mParams = nullptr;
 
  // input data
  std::vector<LocalResid> mLocalResidualsIn;                        
  std::vector<VoxStats> mVoxStatsIn, *mVoxStatsInPtr{&mVoxStatsIn}; 
  // output data
  std::unique_ptr<TFile> mFileOut; 
  std::unique_ptr<TTree> mTreeOut; 
  // status flags
  bool mIsInitialized{}; 
  bool mPrintMem{};      
  // binning
  int mNXBins{param::NPadRows};                            
  int mNY2XBins{param::NY2XBins};                          
  int mNZ2XBins{param::NZ2XBins};                          
  int mNVoxPerSector{};                                    
  float mDX{};                                             
  float mDXI{};                                            
  std::vector<float> mMaxY2X{};                            
  std::vector<float> mDY2X{};                              
  std::vector<float> mDY2XI{};                             
  std::vector<float> mY2XBinsDH{};                         
  std::vector<float> mY2XBinsDI{};                         
  std::vector<float> mY2XBinsCenter{};                     
  float mDZ2X{};                                           
  float mDZ2XI{};                                          
  std::vector<float> mZ2XBinsDH{};                         
  std::vector<float> mZ2XBinsDI{};                         
  std::vector<float> mZ2XBinsCenter{};                     
  float mMaxZ2X{1.f};                                      
  std::array<bool, VoxDim> mUniformBins{true, true, true}; 
  // smoothing
  KernelType mKernelType{KernelType::Epanechnikov};                
  bool mUseErrInSmoothing{true};                                   
  std::array<bool, VoxDim> mSmoothPol2{};                          
  std::array<int, SECTORSPERSIDE * SIDES> mNSmoothingFailedBins{}; 
  std::array<int, VoxDim> mStepKern{};                             
  std::array<float, VoxDim> mKernelScaleEdge{};                    
  std::array<float, VoxDim> mKernelWInv{};                         
  std::array<double, ResDim * sMaxSmtDim> mLastSmoothingRes{};     
  // calibrated parameters
  float mEffVdriftCorr{0.f}; 
  float mEffT0Corr{0.f};     
  // (intermediate) results
  std::array<std::bitset<param::NPadRows>, SECTORSPERSIDE * SIDES> mXBinsIgnore{};          
  std::array<std::array<float, param::NPadRows>, SECTORSPERSIDE * SIDES> mValidFracXBins{}; 
  std::array<std::vector<VoxRes>, SECTORSPERSIDE * SIDES> mVoxelResults{};                  
  VoxRes mVoxelResultsOut{};                                                                
  VoxRes* mVoxelResultsOutPtr{&mVoxelResultsOut};                                           
  std::array<float, 2> mAdhocScalingX{0, 0};                                                
  bool mDoAdhocCorrectionZ2X{false};                                                        
 
  ClassDefNV(TrackResiduals, 3);
};
 
//_____________________________________________________
inline size_t TrackResiduals::getGlbVoxBin(const std::array<unsigned char, VoxDim>& bvox) const
{
  return bvox[VoxX] + (bvox[VoxF] + bvox[VoxZ] * mNY2XBins) * mNXBins;
}
 
//_____________________________________________________
inline size_t TrackResiduals::getGlbVoxBin(int ix, int ip, int iz) const
{
  return ix + (ip + iz * mNY2XBins) * mNXBins;
}
 
//_____________________________________________________
inline void TrackResiduals::getVoxelCoordinates(int isec, int ix, int ip, int iz, float& x, float& p, float& z) const
{
  x = getX(ix);
  p = getY2X(ix, ip);
  z = getZ2X(iz);
  if (isec >= SECTORSPERSIDE) {
    z = -z;
  }
}
 
//_____________________________________________________
inline float TrackResiduals::getDXI(int ix) const
{
  if (mUniformBins[VoxX]) {
    return mDXI;
  } else {
    if (ix < param::NRowsPerROC[0]) {
      // we are in the IROC
      return 1.f / param::RowDX[0];
    } else if (ix > param::NRowsAccumulated[param::NROCTypes - 2]) {
      // we are in the last OROC
      return 1.f / param::RowDX[param::NROCTypes - 1];
    } else if (ix < param::NRowsAccumulated[1]) {
      // OROC1
      return 1.f / param::RowDX[1];
    } else {
      // OROC2
      return 1.f / param::RowDX[2];
    }
  }
}
 
//_____________________________________________________
inline float TrackResiduals::getX(int ix) const
{
  return mUniformBins[VoxX] ? param::MinX + (ix + 0.5) * mDX : param::RowX[ix];
}
 
//_____________________________________________________
inline float TrackResiduals::getMaxY2X(int ix) const
{
  return mMaxY2X[ix];
}
 
//_____________________________________________________
inline float TrackResiduals::getY2X(int ix, int ip) const
{
  if (mUniformBins[VoxF]) {
    return (0.5f + ip) * mDY2X[ix] - mMaxY2X[ix];
  }
  return mMaxY2X[ix] * mY2XBinsCenter[ip];
}
 
//_____________________________________________________
inline float TrackResiduals::getZ2X(int iz) const
{
  if (mUniformBins[VoxZ]) {
    return (0.5f + iz) * getDZ2X();
  }
  return mZ2XBinsCenter[iz];
}
 
//_____________________________________________________
inline float TrackResiduals::getDY2XI(int ix, int iy) const
{
  if (mUniformBins[VoxF]) {
    return mDY2XI[ix];
  }
  return mY2XBinsDI[iy] / mMaxY2X[ix];
}
 
//_____________________________________________________
inline float TrackResiduals::getDZ2X(int iz) const
{
  if (mUniformBins[VoxZ]) {
    return mDZ2X;
  }
  return 2.f * mZ2XBinsDH[iz];
}
 
//_____________________________________________________
inline float TrackResiduals::getDZ2XI(int iz) const
{
  if (mUniformBins[VoxZ]) {
    return mDZ2XI;
  }
  return mZ2XBinsDI[iz];
}
 
//_____________________________________________________
inline void TrackResiduals::findVoxel(float x, float y2x, float z2x, int& ix, int& ip, int& iz) const
{
  ix = getXBin(x);
  ip = getY2XBin(y2x, ix);
  iz = getZ2XBin(z2x);
}
 
//_____________________________________________________
inline int TrackResiduals::getXBinExact(float x) const
{
  if (mUniformBins[VoxX]) {
    int ix = (x - param::MinX) * mDXI;
    return (ix < 0 || ix >= mNXBins) ? -1 : ix;
  }
  return getRowID(x);
}
 
//_____________________________________________________
inline int TrackResiduals::getXBin(float x) const
{
  int bx = getXBinExact(x);
  return bx > -1 ? (bx < mNXBins ? bx : mNXBins - 1) : 0;
}
 
//_____________________________________________________
inline int TrackResiduals::getY2XBinExact(float y2x, int ix) const
{
  if (y2x < -mMaxY2X[ix]) {
    return -1;
  }
  if (y2x > mMaxY2X[ix]) {
    return mNY2XBins;
  }
  if (mUniformBins[VoxF]) {
    return static_cast<int>((y2x + mMaxY2X[ix]) * getDY2XI(ix));
  }
  y2x /= mMaxY2X[ix];
  for (int iBin = 0; iBin < mNY2XBins; ++iBin) {
    if (y2x < mY2XBinsCenter[iBin] + mY2XBinsDH[iBin]) {
      return iBin;
    }
  }
  return mNY2XBins;
}
 
//_____________________________________________________
inline int TrackResiduals::getY2XBin(float y2x, int ix) const
{
  int bp = getY2XBinExact(y2x, ix);
  return bp > -1 ? (bp < mNY2XBins ? bp : mNY2XBins - 1) : 0;
}
 
//_____________________________________________________
inline int TrackResiduals::getZ2XBinExact(float z2x) const
{
  if (mUniformBins[VoxZ]) {
    float bz = z2x * getDZ2XI();
    if (bz >= mNZ2XBins) {
      return -1;
    }
    if (bz < 0) {
      // accounting for clusters which were moved to the wrong side
      // TODO: how can this happen?
      bz = 0;
    }
    return static_cast<int>(bz);
  }
  for (int iBin = 0; iBin < mNZ2XBins; ++iBin) {
    if (z2x < mZ2XBinsCenter[iBin] + mZ2XBinsDH[iBin]) {
      return iBin;
    }
  }
  return -1;
}
 
//_____________________________________________________
inline int TrackResiduals::getZ2XBin(float z2x) const
{
  int iz = getZ2XBinExact(z2x);
  return iz < 0 ? mNZ2XBins - 1 : iz;
}
 
} // namespace tpc
 
} // namespace o2
 
#endif