diff --git a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h index 9a6cd8488f7c8..e5889cc302b38 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h @@ -46,9 +46,10 @@ namespace o2::tpc /// c.setMembers(tpcTrackClIdxVecInput, clusterIndex, tpcTracks); // set the member variables: TrackTPC, TPCClRefElem, o2::tpc::ClusterNativeAccess /// c.setRefit(); // set the refit pointer to perform refitting of tracks, otherwise setPropagateTrack to true /// start looping over the tracks -/// c.calculatedEdx(track, output, 0.01, 0.6, CorrectionFlags::TopologyPol | CorrectionFlags::GainFull | CorrectionFlags::GainResidual | CorrectionFlags::dEdxResidual) // this will fill the dEdxInfo output for given track +/// c.calculatedEdx(track, output, 0.015, 0.60, CorrectionFlags::TopologyPol | CorrectionFlags::dEdxResidual, ClusterFlags::ExcludeEdgeCl) // this will fill the dEdxInfo output for given track enum class CorrectionFlags : unsigned short { + None = 0, TopologySimple = 1 << 0, ///< flag for simple analytical topology correction TopologyPol = 1 << 1, ///< flag for topology correction from polynomials GainFull = 1 << 2, ///< flag for full gain map from calibration container @@ -56,10 +57,23 @@ enum class CorrectionFlags : unsigned short { dEdxResidual = 1 << 4, ///< flag for residual dEdx correction }; +enum class ClusterFlags : unsigned short { + None = 0, + ExcludeSingleCl = 1 << 0, ///< flag to exclude single clusters in dEdx calculation + ExcludeSplitCl = 1 << 1, ///< flag to exclude split clusters in dEdx calculation + ExcludeEdgeCl = 1 << 2, ///< flag to exclude sector edge clusters in dEdx calculation + ExcludeSubthresholdCl = 1 << 3, ///< flag to exclude subthreshold clusters in dEdx calculation + ExcludeSectorBoundaries = 1 << 4, ///< flag to exclude sector boundary clusters in subthreshold cluster treatment +}; + inline CorrectionFlags operator&(CorrectionFlags a, CorrectionFlags b) { return static_cast(static_cast(a) & static_cast(b)); } inline CorrectionFlags operator~(CorrectionFlags a) { return static_cast(~static_cast(a)); } inline CorrectionFlags operator|(CorrectionFlags a, CorrectionFlags b) { return static_cast(static_cast(a) | static_cast(b)); } +inline ClusterFlags operator&(ClusterFlags a, ClusterFlags b) { return static_cast(static_cast(a) & static_cast(b)); } +inline ClusterFlags operator~(ClusterFlags a) { return static_cast(~static_cast(a)); } +inline ClusterFlags operator|(ClusterFlags a, ClusterFlags b) { return static_cast(static_cast(a) | static_cast(b)); } + class CalculatedEdx { public: @@ -74,7 +88,7 @@ class CalculatedEdx void setMembers(std::vector* tpcTrackClIdxVecInput, const o2::tpc::ClusterNativeAccess& clIndex, std::vector* vTPCTracksArrayInp); /// set the refitter - void setRefit(); + void setRefit(const unsigned int nHbfPerTf = 32); /// \param propagate propagate the tracks to extract the track parameters instead of performing a refit void setPropagateTrack(const bool propagate) { mPropagateTrack = propagate; } @@ -88,8 +102,14 @@ class CalculatedEdx /// \param maxMissingCl maximum number of missing clusters for subthreshold check void setMaxMissingCl(int maxMissingCl) { mMaxMissingCl = maxMissingCl; } + /// \param minChargeTotThreshold upper limit for the possible minimum charge tot in subthreshold treatment + void setMinChargeTotThreshold(float minChargeTotThreshold) { mMinChargeTotThreshold = minChargeTotThreshold; } + + /// \param minChargeMaxThreshold upper limit for the possible minimum charge max in subthreshold treatment + void setMinChargeMaxThreshold(float minChargeMaxThreshold) { mMinChargeMaxThreshold = minChargeMaxThreshold; } + /// set the debug streamer - void setStreamer() { mStreamer = std::make_unique("dEdxDebug.root", "recreate"); }; + void setStreamer(const char* debugRootFile) { mStreamer = std::make_unique(debugRootFile, "recreate"); }; /// \return returns magnetic field in kG float getFieldNominalGPUBz() { return mFieldNominalGPUBz; } @@ -97,7 +117,13 @@ class CalculatedEdx /// \return returns maxMissingCl for subthreshold cluster treatment int getMaxMissingCl() { return mMaxMissingCl; } - /// fill missing clusters with minimum charge (method=0) or minimum charge/2 (method=1) + /// \return returns the upper limit for the possible minimum charge tot in subthreshold treatment + float getMinChargeTotThreshold() { return mMinChargeTotThreshold; } + + /// \return returns the upper limit for the possible minimum charge max in subthreshold treatment + float getMinChargeMaxThreshold() { return mMinChargeMaxThreshold; } + + /// fill missing clusters with minimum charge (method=0) or minimum charge/2 (method=1) or Landau (method=2) void fillMissingClusters(int missingClusters[4], float minChargeTot, float minChargeMax, int method); /// get the truncated mean for the input track with the truncation range, charge type, region and corrections @@ -108,11 +134,11 @@ class CalculatedEdx /// \param high higher cluster cut /// \param mask to apply different corrections: TopologySimple = simple analytical topology correction, TopologyPol = topology correction from polynomials, GainFull = full gain map from calibration container, /// GainResidual = residuals gain map from calibration container, dEdxResidual = residual dEdx correction - void calculatedEdx(TrackTPC& track, dEdxInfo& output, float low = 0.05f, float high = 0.6f, CorrectionFlags mask = CorrectionFlags::TopologyPol | CorrectionFlags::GainFull | CorrectionFlags::GainResidual | CorrectionFlags::dEdxResidual); + void calculatedEdx(TrackTPC& track, dEdxInfo& output, float low = 0.015f, float high = 0.6f, CorrectionFlags correctionMask = CorrectionFlags::TopologyPol | CorrectionFlags::dEdxResidual, ClusterFlags clusterMask = ClusterFlags::None, int subthresholdMethod = 0, const char* debugRootFile = "dEdxDebug.root"); /// get the truncated mean for the input charge vector and the truncation range low*nCl& charge, float low, float high) const; @@ -136,14 +162,64 @@ class CalculatedEdx /// \param runNumberOrTimeStamp run number or time stamp void loadCalibsFromCCDB(long runNumberOrTimeStamp); + /// load calibration objects from local CCDB folder + /// \param localCCDBFolder local CCDB folder + void loadCalibsFromLocalCCDBFolder(const char* localCCDBFolder); + + /// load track topology correction from a local file + /// \param folder folder path without a trailing / + /// \param file file path starting with / + /// \param object name of the object to load + void setTrackTopologyCorrectionFromFile(const char* folder, const char* file, const char* object); + + /// load gain map from a local file + /// \param folder folder path without a trailing / + /// \param file file path starting with / + /// \param object name of the object to load + void setGainMapFromFile(const char* folder, const char* file, const char* object); + + /// load gain map residual from a local file + /// \param folder folder path without a trailing / + /// \param file file path starting with / + /// \param object name of the object to load + void setGainMapResidualFromFile(const char* folder, const char* file, const char* object); + + /// load dEdx residual correction from a local file + /// \param folder folder path without a trailing / + /// \param file file path starting with / + /// \param object name of the object to load + void setResidualCorrectionFromFile(const char* folder, const char* file, const char* object); + + /// load zero suppression threshold from a local file + /// \param folder folder path without a trailing / + /// \param file file path starting with / + /// \param object name of the object to load + void setZeroSuppressionThresholdFromFile(const char* folder, const char* file, const char* object); + + /// load magnetic field from a local file + /// \param folder folder path without a trailing / + /// \param file file path starting with / + /// \param object name of the object to load + void setMagneticFieldFromFile(const char* folder, const char* file, const char* object); + + /// load propagator from a local file + /// \param folder folder path without a trailing / + /// \param file file path starting with / + /// \param object name of the object to load + void setPropagatorFromFile(const char* folder, const char* file, const char* object); + private: - std::vector* mTracks{nullptr}; ///< vector containing the tpc tracks which will be processed. + std::vector* mTracks{nullptr}; ///< vector containing the tpc tracks which will be processed std::vector* mTPCTrackClIdxVecInput{nullptr}; ///< input vector with TPC tracks cluster indicies const o2::tpc::ClusterNativeAccess* mClusterIndex{nullptr}; ///< needed to access clusternative with tpctracks - o2::gpu::CorrectionMapsHelper mTPCCorrMapsHelper; ///< cluster corrections map helper + o2::gpu::CorrectionMapsHelper mTPCCorrMapsHelper; ///< cluster correction maps helper + std::vector mTPCRefitterShMap; ///< externally set TPC clusters sharing map + std::vector mTPCRefitterOccMap; ///< externally set TPC clusters occupancy map std::unique_ptr mRefit{nullptr}; ///< TPC refitter used for TPC tracks refit during the reconstruction - int mMaxMissingCl{2}; ///< maximum number of missing clusters for subthreshold check + int mMaxMissingCl{1}; ///< maximum number of missing clusters for subthreshold check + float mMinChargeTotThreshold{50}; ///< upper limit for minimum charge tot value in subthreshold treatment, i.e for a high dEdx track adding a minimum value of 500 to track as a virtual charge doesn't make sense + float mMinChargeMaxThreshold{50}; ///< upper limit for minimum charge max value in subthreshold treatment, i.e for a high dEdx track adding a minimum value of 500 to track as a virtual charge doesn't make sense float mFieldNominalGPUBz{5}; ///< magnetic field in kG, used for track propagation bool mPropagateTrack{false}; ///< propagating the track instead of performing a refit bool mDebug{false}; ///< use the debug streamer @@ -156,4 +232,4 @@ class CalculatedEdx } // namespace o2::tpc -#endif +#endif \ No newline at end of file diff --git a/Detectors/TPC/calibration/src/CalculatedEdx.cxx b/Detectors/TPC/calibration/src/CalculatedEdx.cxx index b70be764078b8..00c8a4c8aa743 100644 --- a/Detectors/TPC/calibration/src/CalculatedEdx.cxx +++ b/Detectors/TPC/calibration/src/CalculatedEdx.cxx @@ -23,7 +23,6 @@ #include "CCDB/BasicCCDBManager.h" #include "TPCBase/CDBInterface.h" #include "TPCReconstruction/TPCFastTransformHelperO2.h" -#include "TPCCalibration/CalibPadGainTracksBase.h" #include "CalibdEdxTrackTopologyPol.h" #include "DataFormatsParameters/GRPMagField.h" #include "GPUO2InterfaceUtils.h" @@ -43,47 +42,44 @@ void CalculatedEdx::setMembers(std::vector* tpcTrackClIdx mClusterIndex = &clIndex; } -void CalculatedEdx::setRefit() +void CalculatedEdx::setRefit(const unsigned int nHbfPerTf) { - mRefit = std::make_unique(mClusterIndex, &mTPCCorrMapsHelper, mFieldNominalGPUBz, mTPCTrackClIdxVecInput->data(), 0, nullptr, nullptr, -1, mTracks); + mTPCRefitterShMap.reserve(mClusterIndex->nClustersTotal); + auto sizeOcc = o2::gpu::GPUO2InterfaceRefit::fillOccupancyMapGetSize(nHbfPerTf, nullptr); + mTPCRefitterOccMap.resize(sizeOcc); + std::fill(mTPCRefitterOccMap.begin(), mTPCRefitterOccMap.end(), 0); + o2::gpu::GPUO2InterfaceRefit::fillSharedClustersAndOccupancyMap(mClusterIndex, *mTracks, mTPCTrackClIdxVecInput->data(), mTPCRefitterShMap.data(), mTPCRefitterOccMap.data(), nHbfPerTf); + mRefit = std::make_unique(mClusterIndex, &mTPCCorrMapsHelper, mFieldNominalGPUBz, mTPCTrackClIdxVecInput->data(), nHbfPerTf, mTPCRefitterShMap.data(), mTPCRefitterOccMap.data(), mTPCRefitterOccMap.size()); } void CalculatedEdx::fillMissingClusters(int missingClusters[4], float minChargeTot, float minChargeMax, int method) { - // adding minimum charge - if (method == 0) { - for (int roc = 0; roc < 4; roc++) { - for (int i = 0; i < missingClusters[roc]; i++) { - mChargeTotROC[roc].emplace_back(minChargeTot); - mChargeTotROC[4].emplace_back(minChargeTot); - - mChargeMaxROC[roc].emplace_back(minChargeMax); - mChargeMaxROC[4].emplace_back(minChargeMax); - } - } + if (method != 0 && method != 1) { + LOGP(info, "Unrecognized subthreshold cluster treatment. Not adding virtual charges to the track!"); + return; } - // adding minimum charge/2 - else if (method == 1) { - for (int roc = 0; roc < 4; roc++) { - for (int i = 0; i < missingClusters[roc]; i++) { - mChargeTotROC[roc].emplace_back(minChargeTot / 2); - mChargeTotROC[4].emplace_back(minChargeTot / 2); + for (int roc = 0; roc < 4; roc++) { + for (int i = 0; i < missingClusters[roc]; i++) { + float chargeTot = (method == 1) ? minChargeTot / 2.f : minChargeTot; + float chargeMax = (method == 1) ? minChargeMax / 2.f : minChargeMax; - mChargeMaxROC[roc].emplace_back(minChargeMax / 2); - mChargeMaxROC[4].emplace_back(minChargeMax / 2); - } + mChargeTotROC[roc].emplace_back(chargeTot); + mChargeTotROC[4].emplace_back(chargeTot); + + mChargeMaxROC[roc].emplace_back(chargeMax); + mChargeMaxROC[4].emplace_back(chargeMax); } } } -void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, float low, float high, CorrectionFlags mask) +void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, float low, float high, CorrectionFlags correctionMask, ClusterFlags clusterMask, int subthresholdMethod, const char* debugRootFile) { // get number of clusters const int nClusters = track.getNClusterReferences(); - int nClsROC[4] = {0}; - int nClsSubThreshROC[4] = {0}; + int nClsROC[4] = {0, 0, 0, 0}; + int nClsSubThreshROC[4] = {0, 0, 0, 0}; mChargeTotROC[0].clear(); mChargeTotROC[1].clear(); @@ -98,11 +94,15 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl mChargeMaxROC[4].clear(); // debug vectors + std::vector excludeClVector; std::vector regionVector; std::vector rowIndexVector; std::vector padVector; - std::vector stackVector; std::vector sectorVector; + std::vector stackVector; + std::vector localXVector; + std::vector localYVector; + std::vector offsPadVector; std::vector topologyCorrVector; std::vector topologyCorrTotVector; @@ -112,19 +112,19 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl std::vector residualCorrTotVector; std::vector residualCorrMaxVector; - std::vector xPositionVector; - std::vector localYVector; - std::vector offsPadVector; - std::vector trackVector; std::vector clVector; if (mDebug) { + excludeClVector.reserve(nClusters); regionVector.reserve(nClusters); rowIndexVector.reserve(nClusters); padVector.reserve(nClusters); stackVector.reserve(nClusters); sectorVector.reserve(nClusters); + localXVector.reserve(nClusters); + localYVector.reserve(nClusters); + offsPadVector.reserve(nClusters); topologyCorrVector.reserve(nClusters); topologyCorrTotVector.reserve(nClusters); topologyCorrMaxVector.reserve(nClusters); @@ -132,9 +132,6 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl gainResidualVector.reserve(nClusters); residualCorrTotVector.reserve(nClusters); residualCorrMaxVector.reserve(nClusters); - xPositionVector.reserve(nClusters); - localYVector.reserve(nClusters); - offsPadVector.reserve(nClusters); trackVector.reserve(nClusters); clVector.reserve(nClusters); } @@ -157,7 +154,33 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl // set sectorIndex, rowIndex, clusterIndexNumb track.getClusterReference(*mTPCTrackClIdxVecInput, iCl, sectorIndex, rowIndex, clusterIndexNumb); - // get x position of the track + // get region, pad, stack and stack ID + const int region = Mapper::REGION[rowIndex]; + const unsigned char pad = std::clamp(static_cast(cl.getPad() + 0.5f), static_cast(0), Mapper::PADSPERROW[region][Mapper::getLocalRowFromGlobalRow(rowIndex)] - 1); // the left side of the pad is defined at e.g. 3.5 and the right side at 4.5 + const CRU cru(Sector(sectorIndex), region); + const auto stack = cru.gemStack(); + StackID stackID{sectorIndex, stack}; + // the stack number for debugging + const int stackNumber = static_cast(stack); + + // get local coordinates, offset and flags + const float localX = o2::tpc::Mapper::instance().getPadCentre(PadPos(rowIndex, pad)).X(); + const float localY = Mapper::instance().getPadCentre(PadPos(rowIndex, pad)).Y(); + const float offsPad = (cl.getPad() - pad) * o2::tpc::Mapper::instance().getPadRegionInfo(Mapper::REGION[rowIndex]).getPadWidth(); + const auto flagsCl = cl.getFlags(); + + int excludeCl = 0; // works as a bit mask + if (((clusterMask & ClusterFlags::ExcludeSingleCl) == ClusterFlags::ExcludeSingleCl) && ((flagsCl & ClusterNative::flagSingle) == ClusterNative::flagSingle)) { + excludeCl += 0b001; // 1 for single cluster + } + if (((clusterMask & ClusterFlags::ExcludeSplitCl) == ClusterFlags::ExcludeSplitCl) && (((flagsCl & ClusterNative::flagSplitPad) == ClusterNative::flagSplitPad) || ((flagsCl & ClusterNative::flagSplitTime) == ClusterNative::flagSplitTime))) { + excludeCl += 0b010; // 2 for split cluster + } + if (((clusterMask & ClusterFlags::ExcludeEdgeCl) == ClusterFlags::ExcludeEdgeCl) && ((flagsCl & ClusterNative::flagEdge) == ClusterNative::flagEdge)) { + excludeCl += 0b100; // 4 for edge cluster + } + + // get the x position of the track const float xPosition = Mapper::instance().getPadCentre(PadPos(rowIndex, 0)).X(); bool check = true; @@ -170,44 +193,82 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl } else { // propagate this track to the plane X=xk (cm) in the field "b" (kG) track.rotate(o2::math_utils::detail::sector2Angle(sectorIndex)); - check = o2::base::Propagator::Instance()->PropagateToXBxByBz(track, xPosition, 0.9f, 2., o2::base::Propagator::MatCorrType::USEMatCorrLUT); + check = o2::base::Propagator::Instance()->PropagateToXBxByBz(track, xPosition, 0.999f, 2., o2::base::Propagator::MatCorrType::USEMatCorrLUT); } if (!check || std::isnan(track.getParam(1))) { + excludeCl += 0b1000; // 8 for failure of track propagation or refit + } + + if (excludeCl != 0) { + // for debugging + if (mDebug) { + excludeClVector.emplace_back(excludeCl); + regionVector.emplace_back(region); + rowIndexVector.emplace_back(rowIndex); + padVector.emplace_back(pad); + sectorVector.emplace_back(sectorIndex); + stackVector.emplace_back(stackNumber); + localXVector.emplace_back(localX); + localYVector.emplace_back(localY); + offsPadVector.emplace_back(offsPad); + trackVector.emplace_back(track); + clVector.emplace_back(cl); + + topologyCorrVector.emplace_back(-999.f); + topologyCorrTotVector.emplace_back(-999.f); + topologyCorrMaxVector.emplace_back(-999.f); + gainVector.emplace_back(-999.f); + gainResidualVector.emplace_back(-999.f); + residualCorrTotVector.emplace_back(-999.f); + residualCorrMaxVector.emplace_back(-999.f); + } + // to avoid counting the skipped cluster as a subthreshold cluster rowIndexOld = rowIndex; sectorIndexOld = sectorIndex; continue; } - // get region and charge value - const int region = Mapper::REGION[rowIndex]; + // get charge values float chargeTot = cl.qTot; float chargeMax = cl.qMax; - // get pad and threshold - const unsigned char pad = std::clamp(static_cast(cl.getPad() + 0.5f), static_cast(0), Mapper::PADSPERROW[region][Mapper::getLocalRowFromGlobalRow(rowIndex)] - 1); // the left side of the pad is defined at e.g. 3.5 and the right side at 4.5 + // get threshold const float threshold = mCalibCont.getZeroSupressionThreshold(sectorIndex, rowIndex, pad); - // get stack and stack ID - const CRU cru(Sector(sectorIndex), region); - const auto stack = cru.gemStack(); - StackID stackID{sectorIndex, stack}; - // find missing clusters int missingClusters = rowIndexOld - rowIndex - 1; - if ((missingClusters > 0) && (missingClusters <= mMaxMissingCl) && (sectorIndexOld == sectorIndex)) { - if (stack == GEMstack::IROCgem) { - nClsSubThreshROC[0] += missingClusters; - nClsROC[0] += missingClusters; - } else if (stack == GEMstack::OROC1gem) { - nClsSubThreshROC[1] += missingClusters; - nClsROC[1] += missingClusters; - } else if (stack == GEMstack::OROC2gem) { - nClsSubThreshROC[2] += missingClusters; - nClsROC[2] += missingClusters; - } else if (stack == GEMstack::OROC3gem) { - nClsSubThreshROC[3] += missingClusters; - nClsROC[3] += missingClusters; + if ((missingClusters > 0) && (missingClusters <= mMaxMissingCl)) { + if ((clusterMask & ClusterFlags::ExcludeSectorBoundaries) == ClusterFlags::ExcludeSectorBoundaries) { + if (sectorIndexOld == sectorIndex) { + if (stack == GEMstack::IROCgem) { + nClsSubThreshROC[0] += missingClusters; + nClsROC[0] += missingClusters; + } else if (stack == GEMstack::OROC1gem) { + nClsSubThreshROC[1] += missingClusters; + nClsROC[1] += missingClusters; + } else if (stack == GEMstack::OROC2gem) { + nClsSubThreshROC[2] += missingClusters; + nClsROC[2] += missingClusters; + } else if (stack == GEMstack::OROC3gem) { + nClsSubThreshROC[3] += missingClusters; + nClsROC[3] += missingClusters; + } + } + } else { + if (stack == GEMstack::IROCgem) { + nClsSubThreshROC[0] += missingClusters; + nClsROC[0] += missingClusters; + } else if (stack == GEMstack::OROC1gem) { + nClsSubThreshROC[1] += missingClusters; + nClsROC[1] += missingClusters; + } else if (stack == GEMstack::OROC2gem) { + nClsSubThreshROC[2] += missingClusters; + nClsROC[2] += missingClusters; + } else if (stack == GEMstack::OROC3gem) { + nClsSubThreshROC[3] += missingClusters; + nClsROC[3] += missingClusters; + } } }; rowIndexOld = rowIndex; @@ -217,12 +278,12 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl float effectiveLength = 1.0f; float effectiveLengthTot = 1.0f; float effectiveLengthMax = 1.0f; - if ((mask & CorrectionFlags::TopologySimple) == CorrectionFlags::TopologySimple) { + if ((correctionMask & CorrectionFlags::TopologySimple) == CorrectionFlags::TopologySimple) { effectiveLength = getTrackTopologyCorrection(track, region); chargeTot /= effectiveLength; chargeMax /= effectiveLength; }; - if ((mask & CorrectionFlags::TopologyPol) == CorrectionFlags::TopologyPol) { + if ((correctionMask & CorrectionFlags::TopologyPol) == CorrectionFlags::TopologyPol) { effectiveLengthTot = getTrackTopologyCorrectionPol(track, cl, region, chargeTot, ChargeType::Tot, threshold); effectiveLengthMax = getTrackTopologyCorrectionPol(track, cl, region, chargeMax, ChargeType::Max, threshold); chargeTot /= effectiveLengthTot; @@ -232,10 +293,10 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl // get gain float gain = 1.0f; float gainResidual = 1.0f; - if ((mask & CorrectionFlags::GainFull) == CorrectionFlags::GainFull) { + if ((correctionMask & CorrectionFlags::GainFull) == CorrectionFlags::GainFull) { gain = mCalibCont.getGain(sectorIndex, rowIndex, pad); }; - if ((mask & CorrectionFlags::GainResidual) == CorrectionFlags::GainResidual) { + if ((correctionMask & CorrectionFlags::GainResidual) == CorrectionFlags::GainResidual) { gainResidual = mCalibCont.getResidualGain(sectorIndex, rowIndex, pad); }; chargeTot /= gain * gainResidual; @@ -244,7 +305,7 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl // get dEdx correction on tgl and sector plane float corrTot = 1.0f; float corrMax = 1.0f; - if ((mask & CorrectionFlags::dEdxResidual) == CorrectionFlags::dEdxResidual) { + if ((correctionMask & CorrectionFlags::dEdxResidual) == CorrectionFlags::dEdxResidual) { corrTot = mCalibCont.getResidualCorrection(stackID, ChargeType::Tot, track.getTgl(), track.getSnp()); corrMax = mCalibCont.getResidualCorrection(stackID, ChargeType::Max, track.getTgl(), track.getSnp()); if (corrTot > 0) { @@ -287,22 +348,17 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl // for debugging if (mDebug) { - // mapping for the stack info - std::map map; - map[GEMstack::IROCgem] = 0; - map[GEMstack::OROC1gem] = 1; - map[GEMstack::OROC2gem] = 2; - map[GEMstack::OROC3gem] = 3; - - const float localY = o2::tpc::Mapper::instance().getPadCentre(o2::tpc::PadPos(rowIndex, pad)).Y(); - const float offsPad = (cl.getPad() - pad) * o2::tpc::Mapper::instance().getPadRegionInfo(Mapper::REGION[rowIndex]).getPadWidth(); - - // filling debug vectors + excludeClVector.emplace_back(0); // cl is successfully processed regionVector.emplace_back(region); rowIndexVector.emplace_back(rowIndex); padVector.emplace_back(pad); - stackVector.emplace_back(map[stack]); sectorVector.emplace_back(sectorIndex); + stackVector.emplace_back(stackNumber); + localXVector.emplace_back(localX); + localYVector.emplace_back(localY); + offsPadVector.emplace_back(offsPad); + trackVector.emplace_back(track); + clVector.emplace_back(cl); topologyCorrVector.emplace_back(effectiveLength); topologyCorrTotVector.emplace_back(effectiveLengthTot); @@ -311,32 +367,32 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl gainResidualVector.emplace_back(gainResidual); residualCorrTotVector.emplace_back(corrTot); residualCorrMaxVector.emplace_back(corrMax); - - xPositionVector.emplace_back(xPosition); - localYVector.emplace_back(localY); - offsPadVector.emplace_back(offsPad); - - trackVector.emplace_back(track); - clVector.emplace_back(cl); }; } // number of clusters - output.NHitsIROC = nClsROC[0] - nClsSubThreshROC[0]; - output.NHitsOROC1 = nClsROC[1] - nClsSubThreshROC[1]; - output.NHitsOROC2 = nClsROC[2] - nClsSubThreshROC[2]; - output.NHitsOROC3 = nClsROC[3] - nClsSubThreshROC[3]; - output.NHitsSubThresholdIROC = nClsROC[0]; output.NHitsSubThresholdOROC1 = nClsROC[1]; output.NHitsSubThresholdOROC2 = nClsROC[2]; output.NHitsSubThresholdOROC3 = nClsROC[3]; - // fill subthreshold clusters - fillMissingClusters(nClsSubThreshROC, minChargeTot, minChargeMax, 0); + // check if the lost clusters are subthreshold clusters based on the charge thresholds + if (minChargeTot <= mMinChargeTotThreshold && minChargeMax <= mMinChargeMaxThreshold) { + output.NHitsIROC = nClsROC[0] - nClsSubThreshROC[0]; + output.NHitsOROC1 = nClsROC[1] - nClsSubThreshROC[1]; + output.NHitsOROC2 = nClsROC[2] - nClsSubThreshROC[2]; + output.NHitsOROC3 = nClsROC[3] - nClsSubThreshROC[3]; - auto chargeTotVector = mChargeTotROC[4]; - auto chargeMaxVector = mChargeMaxROC[4]; + // fill subthreshold clusters if not excluded + if (((clusterMask & ClusterFlags::ExcludeSubthresholdCl) == ClusterFlags::None)) { + fillMissingClusters(nClsSubThreshROC, minChargeTot, minChargeMax, subthresholdMethod); + } + } else { + output.NHitsIROC = nClsROC[0]; + output.NHitsOROC1 = nClsROC[1]; + output.NHitsOROC2 = nClsROC[2]; + output.NHitsOROC3 = nClsROC[3]; + } // calculate dEdx output.dEdxTotIROC = getTruncMean(mChargeTotROC[0], low, high); @@ -354,15 +410,17 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl // for debugging if (mDebug) { if (mStreamer == nullptr) { - setStreamer(); + setStreamer(debugRootFile); } (*mStreamer) << "dEdxDebug" + << "Ncl=" << nClusters + << "excludeClVector=" << excludeClVector << "regionVector=" << regionVector << "rowIndexVector=" << rowIndexVector << "padVector=" << padVector - << "stackVector=" << stackVector << "sectorVector=" << sectorVector + << "stackVector=" << stackVector << "topologyCorrVector=" << topologyCorrVector << "topologyCorrTotVector=" << topologyCorrTotVector << "topologyCorrMaxVector=" << topologyCorrMaxVector @@ -370,13 +428,11 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl << "gainResidualVector=" << gainResidualVector << "residualCorrTotVector=" << residualCorrTotVector << "residualCorrMaxVector=" << residualCorrMaxVector - << "xPositionVector=" << xPositionVector + << "localXVector=" << localXVector << "localYVector=" << localYVector << "offsPadVector=" << offsPadVector << "trackVector=" << trackVector << "clVector=" << clVector - << "chargeTotVector=" << chargeTotVector - << "chargeMaxVector=" << chargeMaxVector << "minChargeTot=" << minChargeTot << "minChargeMax=" << minChargeMax << "output=" << output @@ -487,3 +543,85 @@ void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp) const o2::base::MatLayerCylSet* matLut = o2::base::MatLayerCylSet::rectifyPtrFromFile(cm.get("GLO/Param/MatLUT")); propagator->setMatLUT(matLut); } + +void CalculatedEdx::loadCalibsFromLocalCCDBFolder(const char* localCCDBFolder) +{ + setTrackTopologyCorrectionFromFile(localCCDBFolder, "/TPC/Calib/TopologyGainPiecewise/snapshot.root", "ccdb_object"); + setGainMapFromFile(localCCDBFolder, "/TPC/Calib/PadGainFull/snapshot.root", "ccdb_object"); + setGainMapResidualFromFile(localCCDBFolder, "/TPC/Calib/PadGainResidual/snapshot.root", "ccdb_object"); + setResidualCorrectionFromFile(localCCDBFolder, "/TPC/Calib/TimeGain/snapshot.root", "ccdb_object"); + setZeroSuppressionThresholdFromFile(localCCDBFolder, "/TPC/Config/FEEPad/snapshot.root", "ccdb_object"); + setMagneticFieldFromFile(localCCDBFolder, "/GLO/Config/GRPMagField/snapshot.root", "ccdb_object"); + setPropagatorFromFile(localCCDBFolder, "/GLO/Param/MatLUT/snapshot.root", "ccdb_object"); +} + +void CalculatedEdx::setTrackTopologyCorrectionFromFile(const char* folder, const char* file, const char* object) +{ + o2::tpc::CalibdEdxTrackTopologyPol calibTrackTopology; + calibTrackTopology.loadFromFile(fmt::format("{}{}", folder, file).data(), object); + mCalibCont.setPolTopologyCorrection(calibTrackTopology); +} + +void CalculatedEdx::setGainMapFromFile(const char* folder, const char* file, const char* object) +{ + std::unique_ptr gainMapFile(TFile::Open(fmt::format("{}{}", folder, file).data())); + if (!gainMapFile->IsZombie()) { + LOGP(info, "Using file: {}", gainMapFile->GetName()); + o2::tpc::CalDet* gainMap = (o2::tpc::CalDet*)gainMapFile->Get(object); + mCalibCont.setGainMap(*gainMap, 0., 2.); + } +} + +void CalculatedEdx::setGainMapResidualFromFile(const char* folder, const char* file, const char* object) +{ + std::unique_ptr gainMapResidualFile(TFile::Open(fmt::format("{}{}", folder, file).data())); + if (!gainMapResidualFile->IsZombie()) { + LOGP(info, "Using file: {}", gainMapResidualFile->GetName()); + std::unordered_map>* gainMapResidual = (std::unordered_map>*)gainMapResidualFile->Get(object); + mCalibCont.setGainMapResidual(gainMapResidual->at("GainMap")); + } +} + +void CalculatedEdx::setResidualCorrectionFromFile(const char* folder, const char* file, const char* object) +{ + std::unique_ptr calibdEdxResidualFile(TFile::Open(fmt::format("{}{}", folder, file).data())); + if (!calibdEdxResidualFile->IsZombie()) { + LOGP(info, "Using file: {}", calibdEdxResidualFile->GetName()); + o2::tpc::CalibdEdxCorrection* calibdEdxResidual = (o2::tpc::CalibdEdxCorrection*)calibdEdxResidualFile->Get(object); + mCalibCont.setResidualCorrection(*calibdEdxResidual); + } +} + +void CalculatedEdx::setZeroSuppressionThresholdFromFile(const char* folder, const char* file, const char* object) +{ + std::unique_ptr zeroSuppressionFile(TFile::Open(fmt::format("{}{}", folder, file).data())); + if (!zeroSuppressionFile->IsZombie()) { + LOGP(info, "Using file: {}", zeroSuppressionFile->GetName()); + std::unordered_map>* zeroSupressionThresholdMap = (std::unordered_map>*)zeroSuppressionFile->Get(object); + mCalibCont.setZeroSupresssionThreshold(zeroSupressionThresholdMap->at("ThresholdMap")); + } +} + +void CalculatedEdx::setMagneticFieldFromFile(const char* folder, const char* file, const char* object) +{ + std::unique_ptr magFile(TFile::Open(fmt::format("{}{}", folder, file).data())); + if (!magFile->IsZombie()) { + LOGP(info, "Using file: {}", magFile->GetName()); + o2::parameters::GRPMagField* magField = (o2::parameters::GRPMagField*)magFile->Get(object); + o2::base::Propagator::initFieldFromGRP(magField); + float bz = GPUO2InterfaceUtils::getNominalGPUBz(*magField); + LOGP(info, "Magnetic field: {}", bz); + setFieldNominalGPUBz(bz); + } +} + +void CalculatedEdx::setPropagatorFromFile(const char* folder, const char* file, const char* object) +{ + auto propagator = o2::base::Propagator::Instance(); + std::unique_ptr matLutFile(TFile::Open(fmt::format("{}{}", folder, file).data())); + if (!matLutFile->IsZombie()) { + LOGP(info, "Using file: {}", matLutFile->GetName()); + o2::base::MatLayerCylSet* matLut = o2::base::MatLayerCylSet::rectifyPtrFromFile((o2::base::MatLayerCylSet*)matLutFile->Get(object)); + propagator->setMatLUT(matLut); + } +} \ No newline at end of file