Skip to content

Commit

Permalink
TPC: adding dEdx space-charge correction
Browse files Browse the repository at this point in the history
- adding check for shared clusters
- adding cluster occupancy
  • Loading branch information
matthias-kleiner committed Oct 2, 2024
1 parent 253c99b commit 535d626
Show file tree
Hide file tree
Showing 5 changed files with 357 additions and 4 deletions.
4 changes: 3 additions & 1 deletion Detectors/TPC/calibration/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@ o2_add_library(TPCCalibration
src/CorrMapParam.cxx
src/TPCMShapeCorrection.cxx
src/DigitAdd.cxx
src/CorrectdEdxDistortions.cxx
PUBLIC_LINK_LIBRARIES O2::DataFormatsTPC O2::TPCBase
O2::TPCReconstruction ROOT::Minuit
Microsoft.GSL::GSL
Expand Down Expand Up @@ -111,7 +112,8 @@ o2_target_root_dictionary(TPCCalibration
include/TPCCalibration/TPCScaler.h
include/TPCCalibration/CorrMapParam.h
include/TPCCalibration/TPCMShapeCorrection.h
include/TPCCalibration/DigitAdd.h)
include/TPCCalibration/DigitAdd.h
include/TPCCalibration/CorrectdEdxDistortions.h)

o2_add_test_root_macro(macro/comparePedestalsAndNoise.C
PUBLIC_LINK_LIBRARIES O2::TPCBase
Expand Down
20 changes: 19 additions & 1 deletion Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@
#include "CalibdEdxContainer.h"
#include "CorrectionMapsHelper.h"
#include "CommonUtils/TreeStreamRedirector.h"
#include "TPCCalibration/CorrectdEdxDistortions.h"

#include <vector>

Expand Down Expand Up @@ -55,6 +56,7 @@ enum class CorrectionFlags : unsigned short {
GainFull = 1 << 2, ///< flag for full gain map from calibration container
GainResidual = 1 << 3, ///< flag for residuals gain map from calibration container
dEdxResidual = 1 << 4, ///< flag for residual dEdx correction
dEdxSC = 1 << 5, ///< flag for space-charge dEdx correction
};

enum class ClusterFlags : unsigned short {
Expand All @@ -64,6 +66,7 @@ enum class ClusterFlags : unsigned short {
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
ExcludeSharedCl = 1 << 5, ///< flag to exclude clusters shared between tracks
};

inline CorrectionFlags operator&(CorrectionFlags a, CorrectionFlags b) { return static_cast<CorrectionFlags>(static_cast<unsigned short>(a) & static_cast<unsigned short>(b)); }
Expand Down Expand Up @@ -111,6 +114,9 @@ class CalculatedEdx
/// set the debug streamer
void setStreamer(const char* debugRootFile) { mStreamer = std::make_unique<o2::utils::TreeStreamRedirector>(debugRootFile, "recreate"); };

/// set the debug streamer of the space-charge dedx correction
void setSCStreamer(const char* debugRootFile = "debug_sc_corrections.root") { mSCdEdxCorrection.setStreamer(debugRootFile); }

/// \return returns magnetic field in kG
float getFieldNominalGPUBz() { return mFieldNominalGPUBz; }

Expand Down Expand Up @@ -160,7 +166,8 @@ class CalculatedEdx

/// load calibration objects from CCDB
/// \param runNumberOrTimeStamp run number or time stamp
void loadCalibsFromCCDB(long runNumberOrTimeStamp);
/// \param isMC set if dEdx space-charge corrections will be loaded for MC or real data
void loadCalibsFromCCDB(long runNumberOrTimeStamp, const bool isMC = false);

/// load calibration objects from local CCDB folder
/// \param localCCDBFolder local CCDB folder
Expand Down Expand Up @@ -208,6 +215,15 @@ class CalculatedEdx
/// \param object name of the object to load
void setPropagatorFromFile(const char* folder, const char* file, const char* object);

/// \param lumi set luminosity for space-charge correction map scaling
void setLumi(const float lumi) { mSCdEdxCorrection.setLumi(lumi); }

/// \return returns space-charge dedx correctin
auto& getSCCorrection() { return mSCdEdxCorrection; }

/// \return returns cluster occupancy for given cluster
unsigned int getOccupancy(const o2::tpc::ClusterNative& cl) const;

private:
std::vector<TrackTPC>* mTracks{nullptr}; ///< vector containing the tpc tracks which will be processed
std::vector<TPCClRefElem>* mTPCTrackClIdxVecInput{nullptr}; ///< input vector with TPC tracks cluster indicies
Expand All @@ -228,6 +244,8 @@ class CalculatedEdx

std::array<std::vector<float>, 5> mChargeTotROC;
std::array<std::vector<float>, 5> mChargeMaxROC;

CorrectdEdxDistortions mSCdEdxCorrection; ///< for space-charge correction of dE/dx
};

} // namespace o2::tpc
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,101 @@
// 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.

///
/// @file CorrectdEdxDistortions.h
/// @author Matthias Kleiner, mkleiner@ikf.uni-frankfurt.de
///

#ifndef ALICEO2_TPC_CORRECTDEDXDISTORTIONS_H
#define ALICEO2_TPC_CORRECTDEDXDISTORTIONS_H

#include "GPUTPCGeometry.h"
#include <memory>

// forward declare classes
namespace o2::gpu
{
class TPCFastTransform;
}

namespace o2::utils
{
class TreeStreamRedirector;
}

namespace o2::tpc
{

class CorrectdEdxDistortions
{
public:
/// Default constructor
CorrectdEdxDistortions();

/// Default destructor
~CorrectdEdxDistortions();

/// setting the space-charge corrections from local file
/// \param scAvgFile average space-charge correction map
/// \param scDerFile derivative space-charge correction map
/// \param lumi luminosity which is used for scaling the correction map
void setSCCorrFromFile(const char* scAvgFile, const char* scDerFile, const float lumi = -1);

/// setting the space-charge corrections
/// \param avg average space-charge correction map
/// \param der derivative space-charge correction map
/// \param lumi luminosity which is used for scaling the correction map
void setCorrectionMaps(o2::gpu::TPCFastTransform* avg, o2::gpu::TPCFastTransform* der, const float lumi);

/// setting the space-charge corrections
/// \param avg average space-charge correction map
/// \param der derivative space-charge correction map
void setCorrectionMaps(o2::gpu::TPCFastTransform* avg, o2::gpu::TPCFastTransform* der);

/// \param lumi luminosity which is used for scaling the correction map
void setLumi(float lumi);

/// enable the debug streamer
void setStreamer(const char* debugRootFile = "debug_sc_corrections.root");

/// \return getting the correction value for the cluster charge
/// \param time true time information of the cluster
/// \param sector TPC sector of the cluster
/// \param padrow global padrow
/// \param pad pad number in the padrow
float getCorrection(const float time, unsigned char sector, unsigned char padrow, int pad) const;

/// set minimum allowed lx correction at lower pad
void setMinLX0(const float minLX) { mLX0Min = minLX; }

/// set minimum allowed lx correction upper lower pad
void setMinLX1(const float minLX) { mLX1Min = minLX; }

/// set minimum allowed correction value
void setMinCorr(const float minCorr) { mScCorrMin = minCorr; }

/// set maximum allowed correction value
void setMaxCorr(const float maxCorr) { mScCorrMax = maxCorr; }

private:
float mScaleDer = 0;
std::unique_ptr<o2::gpu::TPCFastTransform> mCorrAvg{nullptr}; ///< correction map for average distortions
std::unique_ptr<o2::gpu::TPCFastTransform> mCorrDer{nullptr}; ///< correction map for distortions
o2::gpu::GPUTPCGeometry mTPCGeometry; ///< geometry information of TPC
std::unique_ptr<o2::utils::TreeStreamRedirector> mStreamer{nullptr}; ///< debug streamer
float mLX0Min{82}; ///< minimum allowed lx position after correction at position of the bottom pad
float mLX1Min{82}; ///< minimum allowed lx position after correction at position of the top pad
float mScCorrMin{0.7}; ///< minimum allowed correction values
float mScCorrMax{1.6}; ///< maximum allowed correction values
};
} // namespace o2::tpc

#endif
63 changes: 61 additions & 2 deletions Detectors/TPC/calibration/src/CalculatedEdx.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
#include "CalibdEdxTrackTopologyPol.h"
#include "DataFormatsParameters/GRPMagField.h"
#include "GPUO2InterfaceUtils.h"
#include "GPUTPCGMMergedTrackHit.h"

using namespace o2::tpc;

Expand Down Expand Up @@ -111,9 +112,12 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
std::vector<float> gainResidualVector;
std::vector<float> residualCorrTotVector;
std::vector<float> residualCorrMaxVector;
std::vector<float> scCorrVector;

std::vector<o2::tpc::TrackTPC> trackVector;
std::vector<o2::tpc::ClusterNative> clVector;
std::vector<unsigned int> occupancyVector;
std::vector<bool> isClusterShared;

if (mDebug) {
excludeClVector.reserve(nClusters);
Expand All @@ -134,6 +138,9 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
residualCorrMaxVector.reserve(nClusters);
trackVector.reserve(nClusters);
clVector.reserve(nClusters);
scCorrVector.reserve(nClusters);
occupancyVector.reserve(nClusters);
isClusterShared.reserve(nClusters);
}

// for missing clusters
Expand All @@ -154,6 +161,10 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
// set sectorIndex, rowIndex, clusterIndexNumb
track.getClusterReference(*mTPCTrackClIdxVecInput, iCl, sectorIndex, rowIndex, clusterIndexNumb);

// check if the cluster is shared
const unsigned int absoluteIndex = mClusterIndex->clusterOffset[sectorIndex][rowIndex] + clusterIndexNumb;
const bool isShared = mRefit ? (mTPCRefitterShMap[absoluteIndex] & GPUCA_NAMESPACE::gpu::GPUTPCGMMergedTrackHit::flagShared) : 0;

// get region, pad, stack and stack ID
const int region = Mapper::REGION[rowIndex];
const unsigned char pad = std::clamp(static_cast<unsigned int>(cl.getPad() + 0.5f), static_cast<unsigned int>(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
Expand All @@ -179,6 +190,9 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
if (((clusterMask & ClusterFlags::ExcludeEdgeCl) == ClusterFlags::ExcludeEdgeCl) && ((flagsCl & ClusterNative::flagEdge) == ClusterNative::flagEdge)) {
excludeCl += 0b100; // 4 for edge cluster
}
if (((clusterMask & ClusterFlags::ExcludeSharedCl) == ClusterFlags::ExcludeSharedCl) && isShared) {
excludeCl += 0b10000; // for shared cluster
}

// get the x position of the track
const float xPosition = Mapper::instance().getPadCentre(PadPos(rowIndex, 0)).X();
Expand Down Expand Up @@ -214,6 +228,8 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
offsPadVector.emplace_back(offsPad);
trackVector.emplace_back(track);
clVector.emplace_back(cl);
occupancyVector.emplace_back(getOccupancy(cl));
isClusterShared.emplace_back(isShared);

topologyCorrVector.emplace_back(-999.f);
topologyCorrTotVector.emplace_back(-999.f);
Expand All @@ -222,6 +238,7 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
gainResidualVector.emplace_back(-999.f);
residualCorrTotVector.emplace_back(-999.f);
residualCorrMaxVector.emplace_back(-999.f);
scCorrVector.emplace_back(-999.f);
}
// to avoid counting the skipped cluster as a subthreshold cluster
rowIndexOld = rowIndex;
Expand Down Expand Up @@ -325,6 +342,19 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
minChargeMax = chargeMax;
};

// space-charge dEdx corrections
const float time = cl.getTime() - track.getTime0(); // ToDo: get correct time from ITS-TPC track if possible
float scCorr = 1.0f;
if ((correctionMask & CorrectionFlags::dEdxSC) == CorrectionFlags::dEdxSC) {
scCorr = mSCdEdxCorrection.getCorrection(time, sectorIndex, rowIndex, pad);
if (scCorr > 0) {
chargeTot /= scCorr;
};
if (corrMax > 0) {
chargeMax /= scCorr;
};
}

if (stack == GEMstack::IROCgem) {
mChargeTotROC[0].emplace_back(chargeTot);
mChargeMaxROC[0].emplace_back(chargeMax);
Expand Down Expand Up @@ -359,6 +389,8 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
offsPadVector.emplace_back(offsPad);
trackVector.emplace_back(track);
clVector.emplace_back(cl);
occupancyVector.emplace_back(getOccupancy(cl));
isClusterShared.emplace_back(isShared);

topologyCorrVector.emplace_back(effectiveLength);
topologyCorrTotVector.emplace_back(effectiveLengthTot);
Expand All @@ -367,6 +399,7 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
gainResidualVector.emplace_back(gainResidual);
residualCorrTotVector.emplace_back(corrTot);
residualCorrMaxVector.emplace_back(corrMax);
scCorrVector.emplace_back(scCorr);
};
}

Expand Down Expand Up @@ -394,6 +427,10 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
output.NHitsOROC3 = nClsROC[3];
}

// copy corrected cluster charges
auto chargeTotVector = mChargeTotROC[4];
auto chargeMaxVector = mChargeMaxROC[4];

// calculate dEdx
output.dEdxTotIROC = getTruncMean(mChargeTotROC[0], low, high);
output.dEdxTotOROC1 = getTruncMean(mChargeTotROC[1], low, high);
Expand Down Expand Up @@ -428,6 +465,7 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
<< "gainResidualVector=" << gainResidualVector
<< "residualCorrTotVector=" << residualCorrTotVector
<< "residualCorrMaxVector=" << residualCorrMaxVector
<< "scCorrVector=" << scCorrVector
<< "localXVector=" << localXVector
<< "localYVector=" << localYVector
<< "offsPadVector=" << offsPadVector
Expand All @@ -436,6 +474,10 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl
<< "minChargeTot=" << minChargeTot
<< "minChargeMax=" << minChargeMax
<< "output=" << output
<< "occupancy=" << occupancyVector
<< "chargeTotVector=" << chargeTotVector
<< "chargeMaxVector=" << chargeMaxVector
<< "isClusterShared=" << isClusterShared
<< "\n";
}
}
Expand Down Expand Up @@ -492,7 +534,7 @@ float CalculatedEdx::getTrackTopologyCorrectionPol(const o2::tpc::TrackTPC& trac
return effectiveLength;
}

void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp)
void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp, const bool isMC)
{
// setup CCDB manager
auto& cm = o2::ccdb::BasicCCDBManager::instance();
Expand Down Expand Up @@ -542,6 +584,15 @@ void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp)
auto propagator = o2::base::Propagator::Instance();
const o2::base::MatLayerCylSet* matLut = o2::base::MatLayerCylSet::rectifyPtrFromFile(cm.get<o2::base::MatLayerCylSet>("GLO/Param/MatLUT"));
propagator->setMatLUT(matLut);

// load sc correction maps
auto avgMap = isMC ? cm.getForTimeStamp<o2::gpu::TPCFastTransform>(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrMapMC), tRun) : cm.getForTimeStamp<o2::gpu::TPCFastTransform>(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrMap), tRun);
avgMap->rectifyAfterReadingFromFile();

auto derMap = isMC ? cm.getForTimeStamp<o2::gpu::TPCFastTransform>(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrDerivMapMC), tRun) : cm.getForTimeStamp<o2::gpu::TPCFastTransform>(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrDerivMap), tRun);
derMap->rectifyAfterReadingFromFile();

mSCdEdxCorrection.setCorrectionMaps(avgMap, derMap);
}

void CalculatedEdx::loadCalibsFromLocalCCDBFolder(const char* localCCDBFolder)
Expand Down Expand Up @@ -624,4 +675,12 @@ void CalculatedEdx::setPropagatorFromFile(const char* folder, const char* file,
o2::base::MatLayerCylSet* matLut = o2::base::MatLayerCylSet::rectifyPtrFromFile((o2::base::MatLayerCylSet*)matLutFile->Get(object));
propagator->setMatLUT(matLut);
}
}
}

unsigned int CalculatedEdx::getOccupancy(const o2::tpc::ClusterNative& cl) const
{
const int nTimeBinsPerOccupBin = 16;
const int iBinOcc = cl.getTime() / nTimeBinsPerOccupBin + 2;
const unsigned int occupancy = mTPCRefitterOccMap.empty() ? -1 : mTPCRefitterOccMap[iBinOcc];
return occupancy;
}
Loading

0 comments on commit 535d626

Please sign in to comment.