CalculatedEdx.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_CalculatedEdx_H
#define AliceO2_TPC_CalculatedEdx_H
 
// o2 includes
#include "DataFormatsTPC/TrackTPC.h"
#include "DataFormatsTPC/ClusterNative.h"
#include "DataFormatsTPC/dEdxInfo.h"
#include "DataFormatsTPC/VDriftCorrFact.h"
#include "TPCBase/Mapper.h"
#include "GPUO2InterfaceRefit.h"
#include "CalibdEdxContainer.h"
#include "CorrectionMapsHelper.h"
#include "CommonUtils/TreeStreamRedirector.h"
#include "TPCCalibration/CorrectdEdxDistortions.h"
#include "TPCFastTransformPOD.h"
#include "GPUCommonRtypes.h"
#include "SimulationDataFormat/MCCompLabel.h"
#include <array>
#include <vector>
#include <map>
#include <set>
#include <unordered_map>
#include <string>
#include <utility>
 
namespace o2::gpu
{
class TPCFastTransform;
}
 
namespace o2::tpc
{
 
struct AverageOccupancy {
  double IROC = 0.;
  double OROC1 = 0.;
  double OROC2 = 0.;
  double OROC3 = 0.;
  ClassDefNV(AverageOccupancy, 1);
};
 
 
enum class CorrectionFlags : unsigned short {
  None = 0,
  TopologySimple = 1 << 0, 
  TopologyPol = 1 << 1,    
  GainFull = 1 << 2,       
  GainResidual = 1 << 3,   
  dEdxResidual = 1 << 4,   
  dEdxSC = 1 << 5,         
};
 
enum class ClusterFlags : unsigned short {
  None = 0,
  ExcludeSingleCl = 1 << 0,         
  ExcludeSplitPadCl = 1 << 1,       
  ExcludeSplitTimeCl = 1 << 2,      
  ExcludeSplitCl = 1 << 3,          
  ExcludeEdgeCl = 1 << 4,           
  ExcludeSubthresholdCl = 1 << 5,   
  ExcludeSectorBoundaries = 1 << 6, 
  ExcludeSharedCl = 1 << 7,         
  ExcludeSamePadRowCl = 1 << 8,     
};
 
inline CorrectionFlags operator&(CorrectionFlags a, CorrectionFlags b) { return static_cast<CorrectionFlags>(static_cast<unsigned short>(a) & static_cast<unsigned short>(b)); }
inline CorrectionFlags operator~(CorrectionFlags a) { return static_cast<CorrectionFlags>(~static_cast<unsigned short>(a)); }
inline CorrectionFlags operator|(CorrectionFlags a, CorrectionFlags b) { return static_cast<CorrectionFlags>(static_cast<unsigned short>(a) | static_cast<unsigned short>(b)); }
 
inline ClusterFlags operator&(ClusterFlags a, ClusterFlags b) { return static_cast<ClusterFlags>(static_cast<unsigned short>(a) & static_cast<unsigned short>(b)); }
inline ClusterFlags operator~(ClusterFlags a) { return static_cast<ClusterFlags>(~static_cast<unsigned short>(a)); }
inline ClusterFlags operator|(ClusterFlags a, ClusterFlags b) { return static_cast<ClusterFlags>(static_cast<unsigned short>(a) | static_cast<unsigned short>(b)); }
 
struct dEdxSettings {
  ClusterFlags clusterMask = ClusterFlags::None;                                                 
  CorrectionFlags correctionMask = CorrectionFlags::TopologyPol | CorrectionFlags::dEdxResidual; 
  unsigned short subthresholdMethod = 0;                                                         
  unsigned short stackBoundaryMethod = 0;                                                        
  unsigned short sameRowClusterMethod = 2;                                                       
  float maxSubthresholdChargeTot = 100000.f;                                                     
  float maxSubthresholdChargeMax = 100000.f;                                                     
  float low = 2.f / 128.f;                                                                       
  float high = 77.f / 128.f;                                                                     
  std::string debugRootFile = "dEdxDebug.root";                                                  
};
 
struct ClInfo {
  unsigned char sectorIndex = 0;
  unsigned char rowIndex = 0;
  bool isShared = false;
};
 
using ClInfoVec = std::vector<ClInfo>;
 
class CalculatedEdx
{
 public:
  CalculatedEdx();
 
  ~CalculatedEdx() = default;
 
  void setMembers(std::vector<o2::tpc::TPCClRefElem>* tpcTrackClIdxVecInput, const o2::tpc::ClusterNativeAccess& clIndex, std::vector<o2::tpc::TrackTPC>* vTPCTracksArrayInp);
 
  void setRefit(const unsigned int nHbfPerTf = 32);
 
  void setTPCCorrMap(const o2::gpu::TPCFastTransform& corrMap, const o2::gpu::TPCFastTransform* corrMapRef = nullptr, float lumiScale = -1.f);
 
  void setTPCVDrift(const o2::tpc::VDriftCorrFact& v);
 
  void setPropagateTrack(const bool propagate) { mPropagateTrack = propagate; }
 
  void setPropagateParams(const bool propagate) { mPropagateParams = propagate; }
 
  void setPropagateToClusterX(const bool propagate) { mPropagateToClusterX = propagate; }
 
  void setDebug(const bool debug) { mDebug = debug; }
 
  void setFieldNominalGPUBz(const float field) { mFieldNominalGPUBz = field; }
 
  void setMaxMissingCl(int maxMissingCl) { mMaxMissingCl = maxMissingCl; }
 
  //           set to 0 to fill every 1-row gap regardless of its position on the track
  void setSubThreshEdgeRows(int n) { mSubThreshEdgeRows = n; }
 
  void setSameRowMaxPadDiff(float d) { mSameRowMaxPadDiff = d; }
 
  void setSameRowMaxTimeDiff(float d) { mSameRowMaxTimeDiff = d; }
 
  void setStreamer(const char* debugRootFile)
  {
    auto& streamer = mStreamers[debugRootFile];
    if (!streamer) {
      streamer = std::make_unique<o2::utils::TreeStreamRedirector>(debugRootFile, "recreate");
    }
  };
 
  void setSCStreamer(const char* debugRootFile = "debug_sc_corrections.root") { mSCdEdxCorrection.setStreamer(debugRootFile); }
 
  void setLumi(const float lumi) { mSCdEdxCorrection.setLumi(lumi); }
 
  float getFieldNominalGPUBz() { return mFieldNominalGPUBz; }
 
  int getMaxMissingCl() { return mMaxMissingCl; }
 
  int getSubThreshEdgeRows() const { return mSubThreshEdgeRows; }
 
  float getSameRowMaxPadDiff() const { return mSameRowMaxPadDiff; }
 
  float getSameRowMaxTimeDiff() const { return mSameRowMaxTimeDiff; }
 
  long getNPropagationFailed() const { return mNPropagationFailed; }
 
  long getNRefitFallback() const { return mNRefitFallback; }
 
  long getNRowsProcessed() const { return mNRowsProcessed; }
 
  const std::vector<long>& getNSubThresholdFilledPerSettings() const { return mNSubThresholdFilledPerSettings; }
 
  void resetDebugCounters()
  {
    mNPropagationFailed = 0;
    mNRefitFallback = 0;
    mNRowsProcessed = 0;
    mNSubThresholdFilledPerSettings.clear();
  }
 
  void fillMissingClusters(int missingClusters[4], const float minChargeTot[4], const float minChargeMax[4], int method, std::array<std::vector<float>, 5>& chargeTotROC, std::array<std::vector<float>, 5>& chargeMaxROC);
 
  void handleSameRowClusters(o2::tpc::TrackTPC& track, std::vector<std::pair<unsigned char, unsigned char>>& rowOrder, std::map<std::pair<unsigned char, unsigned char>, std::vector<int>>& clustersByRow, std::set<std::pair<unsigned char, unsigned char>>& mergeableRows, std::map<int, std::tuple<unsigned char, unsigned char, unsigned int>>& clusterReferencesByIndex);
 
  void handleSameRowClusters(const std::vector<o2::tpc::ClusterNative>& clusters, const ClInfoVec& clusterInfos, std::vector<std::pair<unsigned char, unsigned char>>& rowOrder, std::map<std::pair<unsigned char, unsigned char>, std::vector<int>>& clustersByRow, std::set<std::pair<unsigned char, unsigned char>>& mergeableRows);
 
  void calculatedEdx(TrackTPC& track, dEdxInfo& output, AverageOccupancy& averageOcc, float low = 2.f / 128.f, float high = 77.f / 128.f, CorrectionFlags correctionMask = CorrectionFlags::TopologyPol | CorrectionFlags::dEdxResidual, ClusterFlags clusterMask = ClusterFlags::None, int subthresholdMethod = 0, int stackBoundaryMethod = 0, const char* debugRootFile = "dEdxDebug.root", float maxSubthresholdChargeTot = 100000.f, float maxSubthresholdChargeMax = 100000.f, int sameRowClusterMethod = 2);
 
  void calculatedEdxMultipleSettings(TrackTPC& track, std::vector<dEdxInfo>& outputs, AverageOccupancy& averageOcc, const std::vector<dEdxSettings>& settingsList, const MCCompLabel* mcLabel = nullptr);
 
  void calculatedEdx(TrackTPC& track, const std::vector<o2::tpc::ClusterNative>& clusters, const ClInfoVec& clusterInfos, dEdxInfo& output, AverageOccupancy& averageOcc, float low = 2.f / 128.f, float high = 77.f / 128.f, CorrectionFlags correctionMask = CorrectionFlags::TopologyPol | CorrectionFlags::dEdxResidual, ClusterFlags clusterMask = ClusterFlags::None, int subthresholdMethod = 0, int stackBoundaryMethod = 0, const char* debugRootFile = "dEdxDebug.root", float maxSubthresholdChargeTot = 100000.f, float maxSubthresholdChargeMax = 100000.f, int sameRowClusterMethod = 2);
 
  void calculatedEdxMultipleSettings(TrackTPC& track, const std::vector<o2::tpc::ClusterNative>& clusters, const ClInfoVec& clusterInfos, std::vector<dEdxInfo>& outputs, AverageOccupancy& averageOcc, const std::vector<dEdxSettings>& settingsList, const MCCompLabel* mcLabel = nullptr);
 
  float getTruncMean(std::vector<float>& charge, float low, float high) const;
 
  float getTrackTopologyCorrection(const o2::tpc::TrackTPC& track, const unsigned int region) const;
 
  float getTrackTopologyCorrectionPol(const o2::tpc::TrackTPC& track, const o2::tpc::ClusterNative& cl, const unsigned int region, const float charge, ChargeType chargeType, const float threshold) const;
 
  auto& getSCCorrection() { return mSCdEdxCorrection; }
 
  unsigned int getOccupancy(float clTime) const;
 
  bool isInStackBoundaries(int stackNumber, unsigned char rowIndex, int stackBoundaryMethod);
 
  void loadCalibsFromCCDB(long runNumberOrTimeStamp, const bool isMC = false, const bool loadSCCorrMap = true, const bool useSCCorrMap = true, const bool useVDrift = true, const float corrMapLumiScale = -1.f);
 
  float getTPCScalerLumiScale(long timeStamp, float meanIDC) const;
 
  void loadCalibsFromLocalCCDBFolder(const char* localCCDBFolder, const bool useSCCorrMap = true, const bool useVDrift = true);
 
  void setTrackTopologyCorrectionFromFile(const char* folder, const char* file, const char* object);
 
  void setGainMapFromFile(const char* folder, const char* file, const char* object);
 
  void setGainMapResidualFromFile(const char* folder, const char* file, const char* object);
 
  void setResidualCorrectionFromFile(const char* folder, const char* file, const char* object);
 
  void setZeroSuppressionThresholdFromFile(const char* folder, const char* file, const char* object);
 
  void setMagneticFieldFromFile(const char* folder, const char* file, const char* object);
 
  void setPropagatorFromFile(const char* folder, const char* file, const char* object);
 
  void setTPCCorrMapFromFile(const char* folder, const char* file, const char* object);
 
  void setVDriftFromFile(const char* folder, const char* file, const char* object);
 
 private:
  struct RowFragment {
    o2::tpc::ClusterNative cl;
    unsigned char pad; 
    float threshold;
    float gain;
    float gainResidual;
    unsigned int occupancy;
    bool isShared;
  };
 
  static const RowFragment& pickDominantFragment(const std::vector<RowFragment>& fragments);
 
  struct RowClusterData {
    std::vector<RowFragment> fragments; 
    bool mergeable;                     
    RowFragment mergedFragment{};       
    o2::tpc::TrackTPC trackSnapshot;    
    unsigned char sectorIndex;
    unsigned char rowIndex;
    unsigned int region;
    GEMstack stack;
    int stackNumber;
    StackID stackID;
    bool propagationFailed;                    
    bool refitFellBack;                        
    int missingClusters;                       
    std::array<int, 4> missingClustersStack{}; 
    bool sameSectorAsPrevRow;                  
    bool missingClusterGapDeadOrEdge;          
    std::vector<int> inputClusterIndices;      
  };
 
  void gatherRowClusterData(o2::tpc::TrackTPC& track, std::vector<RowClusterData>& rowData, AverageOccupancy& averageOcc);
 
  void gatherRowClusterData(o2::tpc::TrackTPC& track, const std::vector<o2::tpc::ClusterNative>& clusters, const ClInfoVec& clusterInfos, std::vector<RowClusterData>& rowData, AverageOccupancy& averageOcc);
 
  void gatherRowClusterDataForRow(o2::tpc::TrackTPC& track, const std::vector<o2::tpc::ClusterNative>& fragmentClusters, const std::vector<bool>& fragmentIsShared, unsigned char sectorIndex, unsigned char rowIndex, bool mergeable, unsigned char rowIndexOld, unsigned char sectorIndexOld, std::array<std::vector<unsigned int>, 4>& occupancyROC, std::vector<RowClusterData>& rowData, bool refitFailed);
 
  bool refitTrack(o2::tpc::TrackTPC& track);
 
  bool refitTrack(o2::tpc::TrackTPC& track, const std::vector<o2::tpc::ClusterNative>& clusters, const ClInfoVec& clusterInfos);
 
  bool propagateTrackToX(o2::track::TrackParCov& track, float xPosition, unsigned char sectorIndex) const;
 
  void rebuildTPCCorrMapPOD();
 
  void calculatedEdxFromRowData(const std::vector<RowClusterData>& rowData, const dEdxSettings& settings, size_t settingsIndex, float trackTime0, const o2::tpc::TrackTPC& trackOrig, const AverageOccupancy& averageOcc, dEdxInfo& output, const MCCompLabel* mcLabel = nullptr);
 
  std::vector<TrackTPC>* mTracks{nullptr};                    
  std::vector<TPCClRefElem>* mTPCTrackClIdxVecInput{nullptr}; 
  const o2::tpc::ClusterNativeAccess* mClusterIndex{nullptr}; 
  const o2::gpu::TPCFastTransformPOD* mTPCCorrMap{nullptr};   
  o2::gpu::aligned_unique_buffer_ptr<o2::gpu::TPCFastTransformPOD> mTPCCorrMapBuffer;
  std::unique_ptr<o2::gpu::TPCFastTransform> mTPCCorrMapFull;               
  std::vector<unsigned char> mTPCRefitterShMap;                             
  std::vector<unsigned int> mTPCRefitterOccMap;                             
  std::unique_ptr<o2::gpu::GPUO2InterfaceRefit> mRefit{nullptr};            
  std::unique_ptr<o2::gpu::GPUO2InterfaceRefit> mRefitExt{nullptr};         
  std::unique_ptr<o2::tpc::ClusterNativeAccess> mExtClusterAccess{nullptr}; 
  std::vector<o2::tpc::ClusterNative> mExtClusters;                         
  std::vector<o2::tpc::TPCClRefElem> mExtClusterRefs;                       
  std::vector<unsigned char> mExtShMap;                                     
  bool mUseRefit{false};                                                    
  unsigned int mNHbfPerTf{32};                                              
 
  int mMaxMissingCl{1};                                                                         
  int mSubThreshEdgeRows{30};                                                                   
  float mSameRowMaxPadDiff{3.f};                                                                
  float mSameRowMaxTimeDiff{4.f};                                                               
  float mFieldNominalGPUBz{5};                                                                  
  bool mPropagateTrack{false};                                                                  
  bool mPropagateParams{false};                                                                 
  bool mPropagateToClusterX{true};                                                              
  bool mDebug{false};                                                                           
  CalibdEdxContainer mCalibCont;                                                                
  std::unordered_map<std::string, std::unique_ptr<o2::utils::TreeStreamRedirector>> mStreamers; 
  long mDebugTrackIndex{-1};                                                                    
  bool mWarnedMissingRefit{false};                                                              
  long mNPropagationFailed{0};                                                                  
  long mNRefitFallback{0};                                                                      
  long mNRowsProcessed{0};                                                                      
  std::vector<long> mNSubThresholdFilledPerSettings;                                            
 
  CorrectdEdxDistortions mSCdEdxCorrection; 
 
  std::array<std::vector<unsigned char>, 4> mStackBoundaries = {{{0, 62}, {63, 96}, {97, 126}, {127, 151}}}; // for excluding stack boundaries in dEdx calculation
};
 
} // namespace o2::tpc
 
#endif