From 8b67fc83a778f8bc47673964043427745e9a7495 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner <48915672+matthias-kleiner@users.noreply.github.com> Date: Fri, 1 Nov 2024 08:54:34 +0100 Subject: [PATCH] TPC: Adding option to make timegain calib with gausian fits (#13641) * TPC: Adding option to make timegain calib with gausian fits - setting tighter default momentum range for filtering MIPs to reject kaons and electrons - adding option to write CalibdEdx to local file by converting the boost histogram to a ROOT histogram - lowering minimum dEdx to 10 for better stabillity in A16 other changes: - TPC: make calculation of dedx parallelisable * TPC: Fixing variable names --- .../DataFormatsTPC/CalibdEdxCorrection.h | 4 +- .../Detectors/TPC/src/CalibdEdxCorrection.cxx | 8 +- .../include/TPCCalibration/CalculatedEdx.h | 5 +- .../include/TPCCalibration/CalibdEdx.h | 35 +- .../include/TPCCalibration/CalibratordEdx.h | 4 +- .../TPC/calibration/src/CalculatedEdx.cxx | 76 ++-- Detectors/TPC/calibration/src/CalibdEdx.cxx | 340 +++++++++++++++++- .../TPC/calibration/src/CalibratordEdx.cxx | 2 +- Detectors/TPC/workflow/src/CalibdEdxSpec.cxx | 12 +- .../TPC/workflow/src/CalibratordEdxSpec.cxx | 7 +- .../TPC/workflow/src/MIPTrackFilterSpec.cxx | 6 +- 11 files changed, 424 insertions(+), 75 deletions(-) diff --git a/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxCorrection.h b/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxCorrection.h index 482c577ec9083..22ee80992f432 100644 --- a/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxCorrection.h +++ b/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxCorrection.h @@ -91,8 +91,8 @@ class CalibdEdxCorrection void clear(); - void writeToFile(std::string_view fileName) const; - void loadFromFile(std::string_view fileName); + void writeToFile(std::string_view fileName, std::string_view objName = "CalibdEdxCorrection") const; + void loadFromFile(std::string_view fileName, std::string_view objName = "CalibdEdxCorrection"); /// \param outFileName name of the output file void dumpToTree(const char* outFileName = "calib_dedx.root") const; diff --git a/DataFormats/Detectors/TPC/src/CalibdEdxCorrection.cxx b/DataFormats/Detectors/TPC/src/CalibdEdxCorrection.cxx index 1a66eb962cdf0..c8224aca5b344 100644 --- a/DataFormats/Detectors/TPC/src/CalibdEdxCorrection.cxx +++ b/DataFormats/Detectors/TPC/src/CalibdEdxCorrection.cxx @@ -36,16 +36,16 @@ void CalibdEdxCorrection::clear() mDims = -1; } -void CalibdEdxCorrection::writeToFile(std::string_view fileName) const +void CalibdEdxCorrection::writeToFile(std::string_view fileName, std::string_view objName) const { std::unique_ptr file(TFile::Open(fileName.data(), "recreate")); - file->WriteObject(this, "CalibdEdxCorrection"); + file->WriteObject(this, objName.data()); } -void CalibdEdxCorrection::loadFromFile(std::string_view fileName) +void CalibdEdxCorrection::loadFromFile(std::string_view fileName, std::string_view objName) { std::unique_ptr file(TFile::Open(fileName.data())); - auto tmp = file->Get("CalibdEdxCorrection"); + auto tmp = file->Get(objName.data()); if (tmp != nullptr) { *this = *tmp; } diff --git a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h index 701c840e060eb..3a744d2b1cfb4 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h @@ -130,7 +130,7 @@ class CalculatedEdx 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); + void fillMissingClusters(int missingClusters[4], float minChargeTot, float minChargeMax, int method, std::array, 5>& chargeTotROC, std::array, 5>& chargeMaxROC); /// get the truncated mean for the input track with the truncation range, charge type, region and corrections /// the cluster charge is normalized by effective length*gain, you can turn off the normalization by setting all corrections to false @@ -242,9 +242,6 @@ class CalculatedEdx CalibdEdxContainer mCalibCont; ///< calibration container std::unique_ptr mStreamer{nullptr}; ///< debug streamer - std::array, 5> mChargeTotROC; - std::array, 5> mChargeMaxROC; - CorrectdEdxDistortions mSCdEdxCorrection; ///< for space-charge correction of dE/dx }; diff --git a/Detectors/TPC/calibration/include/TPCCalibration/CalibdEdx.h b/Detectors/TPC/calibration/include/TPCCalibration/CalibdEdx.h index 470e2c84b5b60..b40daa7b6e61f 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/CalibdEdx.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/CalibdEdx.h @@ -33,6 +33,8 @@ #include #include "THn.h" +class TLinearFitter; + namespace o2::tpc { @@ -74,7 +76,7 @@ class CalibdEdx /// \param angularBins number of bins for Tgl and Snp /// \param fitSnp enable Snp correction - CalibdEdx(int dEdxBins = 60, float mindEdx = 20, float maxdEdx = 90, int angularBins = 36, bool fitSnp = false); + CalibdEdx(int dEdxBins = 70, float mindEdx = 10, float maxdEdx = 90, int angularBins = 36, bool fitSnp = false); void setCuts(const TrackCuts& cuts) { mCuts = cuts; } void setApplyCuts(bool apply) { mApplyCuts = apply; } @@ -124,7 +126,8 @@ class CalibdEdx /// Compute MIP position from dEdx histograms and save result in the correction container. /// To retrieve the correction call `CalibdEdx::getCalib()` - void finalize(); + /// \param useGausFits make gaussian fits of dEdx vs tgl instead of fitting the mean dEdx + void finalize(const bool useGausFits = true); /// Return calib data histogram const Hist& getHist() const { return mHist; } @@ -170,6 +173,17 @@ class CalibdEdx constexpr static float scaleTgl(float tgl, GEMstack rocType) { return tgl / conf_dedx_corr::TglScale[rocType]; } constexpr static float recoverTgl(float scaledTgl, GEMstack rocType) { return scaledTgl * conf_dedx_corr::TglScale[rocType]; } + /// dump this object to a file - the boost histogram is converted to a ROOT histogram - + void dumpToFile(const char* outFile, const char* outName) const; + + /// read the object from a file + static CalibdEdx readFromFile(const char* inFile, const char* inName); + + /// set lower and upper range in units of sigma which are used for the gaussian fits + /// \param lowerSigma low sigma range + /// \param upperSigma upper sigma range + void setSigmaFitRange(const float lowerSigma, const float upperSigma); + private: // ATTENTION: Adjust copy constructor bool mFitSnp{}; @@ -182,15 +196,26 @@ class CalibdEdx float mFitLowCutFactor = 1.5; ///< dEdx cut multiplier for the lower dE/dx range int mFitPasses = 3; ///< number of fit passes used to remove electron tracks TFIDInfo mTFID{}; ///< current TFID + float mSigmaUpper = 1.; ///< mSigma*sigma_gaus used for cutting electrons in case gaussian fits are performed + float mSigmaLower = 1.5; ///< mSigma*sigma_gaus used for cutting electrons in case gaussian fits are performed - Hist mHist; ///< dEdx multidimensional histogram + Hist mHist; /// mDebugOutputStreamer; ///< Debug output streamer - ClassDefNV(CalibdEdx, 4); + std::unique_ptr mDebugOutputStreamer; /// values) { mElectronCut = values; } void setMaterialType(o2::base::Propagator::MatCorrType materialType) { mMatType = materialType; } + void setMakeGaussianFits(const bool makeGaussianFits) { mMakeGaussianFits = makeGaussianFits; } /// \brief Check if there are enough data to compute the calibration. /// \return false if any of the histograms has less entries than mMinEntries @@ -117,6 +118,7 @@ class CalibratordEdx final : public o2::calibration::TimeSlotCalibration mElectronCut{}; ///< Values passed to CalibdEdx::setElectronCut TrackCuts mCuts; ///< Cut object o2::base::Propagator::MatCorrType mMatType{}; ///< material type for track propagation + bool mMakeGaussianFits{}; ///< fit mean of gaussian fits instead of mean dedx TFinterval mTFIntervals; ///< start and end time frame IDs of each calibration time slots TimeInterval mTimeIntervals; ///< start and end times of each calibration time slots @@ -124,7 +126,7 @@ class CalibratordEdx final : public o2::calibration::TimeSlotCalibration mDebugOutputStreamer; ///< Debug output streamer - ClassDefOverride(CalibratordEdx, 2); + ClassDefOverride(CalibratordEdx, 3); }; } // namespace o2::tpc diff --git a/Detectors/TPC/calibration/src/CalculatedEdx.cxx b/Detectors/TPC/calibration/src/CalculatedEdx.cxx index 9809cc454e94f..2ac3b44938bce 100644 --- a/Detectors/TPC/calibration/src/CalculatedEdx.cxx +++ b/Detectors/TPC/calibration/src/CalculatedEdx.cxx @@ -53,7 +53,7 @@ void CalculatedEdx::setRefit(const unsigned int 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) +void CalculatedEdx::fillMissingClusters(int missingClusters[4], float minChargeTot, float minChargeMax, int method, std::array, 5>& chargeTotROC, std::array, 5>& chargeMaxROC) { if (method != 0 && method != 1) { LOGP(info, "Unrecognized subthreshold cluster treatment. Not adding virtual charges to the track!"); @@ -65,11 +65,11 @@ void CalculatedEdx::fillMissingClusters(int missingClusters[4], float minChargeT float chargeTot = (method == 1) ? minChargeTot / 2.f : minChargeTot; float chargeMax = (method == 1) ? minChargeMax / 2.f : minChargeMax; - mChargeTotROC[roc].emplace_back(chargeTot); - mChargeTotROC[4].emplace_back(chargeTot); + chargeTotROC[roc].emplace_back(chargeTot); + chargeTotROC[4].emplace_back(chargeTot); - mChargeMaxROC[roc].emplace_back(chargeMax); - mChargeMaxROC[4].emplace_back(chargeMax); + chargeMaxROC[roc].emplace_back(chargeMax); + chargeMaxROC[4].emplace_back(chargeMax); } } } @@ -82,17 +82,13 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl int nClsROC[4] = {0, 0, 0, 0}; int nClsSubThreshROC[4] = {0, 0, 0, 0}; - mChargeTotROC[0].clear(); - mChargeTotROC[1].clear(); - mChargeTotROC[2].clear(); - mChargeTotROC[3].clear(); - mChargeTotROC[4].clear(); - - mChargeMaxROC[0].clear(); - mChargeMaxROC[1].clear(); - mChargeMaxROC[2].clear(); - mChargeMaxROC[3].clear(); - mChargeMaxROC[4].clear(); + const int nType = 5; + std::array, nType> chargeTotROC; + std::array, nType> chargeMaxROC; + for (int i = 0; i < nType; ++i) { + chargeTotROC[i].reserve(Mapper::PADROWS); + chargeMaxROC[i].reserve(Mapper::PADROWS); + } // debug vectors std::vector excludeClVector; @@ -356,25 +352,25 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl } if (stack == GEMstack::IROCgem) { - mChargeTotROC[0].emplace_back(chargeTot); - mChargeMaxROC[0].emplace_back(chargeMax); + chargeTotROC[0].emplace_back(chargeTot); + chargeMaxROC[0].emplace_back(chargeMax); nClsROC[0]++; } else if (stack == GEMstack::OROC1gem) { - mChargeTotROC[1].emplace_back(chargeTot); - mChargeMaxROC[1].emplace_back(chargeMax); + chargeTotROC[1].emplace_back(chargeTot); + chargeMaxROC[1].emplace_back(chargeMax); nClsROC[1]++; } else if (stack == GEMstack::OROC2gem) { - mChargeTotROC[2].emplace_back(chargeTot); - mChargeMaxROC[2].emplace_back(chargeMax); + chargeTotROC[2].emplace_back(chargeTot); + chargeMaxROC[2].emplace_back(chargeMax); nClsROC[2]++; } else if (stack == GEMstack::OROC3gem) { - mChargeTotROC[3].emplace_back(chargeTot); - mChargeMaxROC[3].emplace_back(chargeMax); + chargeTotROC[3].emplace_back(chargeTot); + chargeMaxROC[3].emplace_back(chargeMax); nClsROC[3]++; }; - mChargeTotROC[4].emplace_back(chargeTot); - mChargeMaxROC[4].emplace_back(chargeMax); + chargeTotROC[4].emplace_back(chargeTot); + chargeMaxROC[4].emplace_back(chargeMax); // for debugging if (mDebug) { @@ -418,7 +414,7 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl // fill subthreshold clusters if not excluded if (((clusterMask & ClusterFlags::ExcludeSubthresholdCl) == ClusterFlags::None)) { - fillMissingClusters(nClsSubThreshROC, minChargeTot, minChargeMax, subthresholdMethod); + fillMissingClusters(nClsSubThreshROC, minChargeTot, minChargeMax, subthresholdMethod, chargeTotROC, chargeMaxROC); } } else { output.NHitsIROC = nClsROC[0]; @@ -428,21 +424,21 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl } // copy corrected cluster charges - auto chargeTotVector = mChargeTotROC[4]; - auto chargeMaxVector = mChargeMaxROC[4]; + auto chargeTotVector = mDebug ? chargeTotROC[4] : std::vector(); + auto chargeMaxVector = mDebug ? chargeMaxROC[4] : std::vector(); // calculate dEdx - output.dEdxTotIROC = getTruncMean(mChargeTotROC[0], low, high); - output.dEdxTotOROC1 = getTruncMean(mChargeTotROC[1], low, high); - output.dEdxTotOROC2 = getTruncMean(mChargeTotROC[2], low, high); - output.dEdxTotOROC3 = getTruncMean(mChargeTotROC[3], low, high); - output.dEdxTotTPC = getTruncMean(mChargeTotROC[4], low, high); - - output.dEdxMaxIROC = getTruncMean(mChargeMaxROC[0], low, high); - output.dEdxMaxOROC1 = getTruncMean(mChargeMaxROC[1], low, high); - output.dEdxMaxOROC2 = getTruncMean(mChargeMaxROC[2], low, high); - output.dEdxMaxOROC3 = getTruncMean(mChargeMaxROC[3], low, high); - output.dEdxMaxTPC = getTruncMean(mChargeMaxROC[4], low, high); + output.dEdxTotIROC = getTruncMean(chargeTotROC[0], low, high); + output.dEdxTotOROC1 = getTruncMean(chargeTotROC[1], low, high); + output.dEdxTotOROC2 = getTruncMean(chargeTotROC[2], low, high); + output.dEdxTotOROC3 = getTruncMean(chargeTotROC[3], low, high); + output.dEdxTotTPC = getTruncMean(chargeTotROC[4], low, high); + + output.dEdxMaxIROC = getTruncMean(chargeMaxROC[0], low, high); + output.dEdxMaxOROC1 = getTruncMean(chargeMaxROC[1], low, high); + output.dEdxMaxOROC2 = getTruncMean(chargeMaxROC[2], low, high); + output.dEdxMaxOROC3 = getTruncMean(chargeMaxROC[3], low, high); + output.dEdxMaxTPC = getTruncMean(chargeMaxROC[4], low, high); // for debugging if (mDebug) { diff --git a/Detectors/TPC/calibration/src/CalibdEdx.cxx b/Detectors/TPC/calibration/src/CalibdEdx.cxx index 114081f57c2f0..bf86db78664d7 100644 --- a/Detectors/TPC/calibration/src/CalibdEdx.cxx +++ b/Detectors/TPC/calibration/src/CalibdEdx.cxx @@ -21,6 +21,7 @@ #include #include #include +#include // o2 includes #include "CommonConstants/PhysicsConstants.h" @@ -31,6 +32,7 @@ #include "Framework/Logger.h" #include "TPCBase/ParameterGas.h" #include "TPCBase/Utils.h" +#include "CommonUtils/BoostHistogramUtils.h" // root includes #include "TFile.h" @@ -139,7 +141,7 @@ void CalibdEdx::fill(const TrackTPC& track) static bool reported = false; if (!reported && mCalibIn.getDims() >= 0) { const auto meanParamTot = mCalibIn.getMeanParams(ChargeType::Tot); - LOGP(info, "Undoing previously apllied corrections with mean qTot Params {}", utils::elementsToString(meanParamTot)); + LOGP(info, "Undoing previously applied corrections with mean qTot Params {}", utils::elementsToString(meanParamTot)); reported = true; } @@ -179,9 +181,11 @@ void CalibdEdx::merge(const CalibdEdx* other) template void fitHist(const Hist& hist, CalibdEdxCorrection& corr, TLinearFitter& fitter, - const float dEdxCut, const float dEdxLowCutFactor, const int passes, const CalibdEdxCorrection* stackMean = nullptr) + const float dEdxCut, const float dEdxLowCutFactor, const int passes, const CalibdEdxCorrection* stackMean = nullptr, o2::utils::TreeStreamRedirector* streamer = nullptr) { + using timer = std::chrono::high_resolution_clock; using ax = CalibdEdx::Axis; + auto start = timer::now(); // number of bins per stack int stackBins = 1; @@ -242,6 +246,28 @@ void fitHist(const Hist& hist, CalibdEdxCorrection& corr, TLinearFitter& fitter, } const double error = 1. / sqrt(counts); fitter.AddPoint(inputs, dEdx, error); + + if (streamer) { + float oldCorr = corr.getCorrection(id, charge, inputs[0], inputs[1]); + float lowerCut = (1.f - dEdxLowCutFactor * dEdxCut) * oldCorr; + float upperCut = (1.f + dEdxCut) * oldCorr; + + (*streamer) << "fit_standard" + << "dedx=" << dEdx + << "itgl=" << hist.axis(ax::Tgl).index(entry->bin(ax::Tgl).center()) + << "snp=" << inputs[1] + << "iter=" << fitPass + << "ifit=" << fit + << "bin=" << bin + << "isector=" << int(id.sector) + << "istack=" << int(id.type) + << "icharge=" << int(charge) + << "counts=" << counts + << "oldCorr=" << oldCorr + << "lowerCut=" << lowerCut + << "upperCut=" << upperCut + << "\n"; + } } } fitter.Eval(); @@ -276,9 +302,201 @@ void fitHist(const Hist& hist, CalibdEdxCorrection& corr, TLinearFitter& fitter, id.sector, int(id.type), int(charge), fitPass, (float)outliers / (float)entries * 100, entries, fitter.GetNpoints(), params[0]); } } + auto stop = timer::now(); + std::chrono::duration time = stop - start; + LOGP(info, "Calibration fits took: {}", time.count()); } -void CalibdEdx::finalize() +template +auto ProjectBoostHistoXFast(const Hist& hist, std::vector& bin_indices, int axis) +{ + const unsigned int nbins = hist.axis(axis).size(); + const double binStartX = hist.axis(axis).bin(0).lower(); + const double binEndX = hist.axis(axis).bin(nbins - 1).upper(); + auto histoProj = boost::histogram::make_histogram(CalibdEdx::FloatAxis(nbins, binStartX, binEndX)); + + // loop over all bins of the requested axis + for (int i = 0; i < nbins; ++i) { + // setting bin of the requested axis + bin_indices[axis] = i; + + // access the bin content specified by bin_indices + histoProj.at(i) = hist.at(bin_indices); + } + + return histoProj; +} + +void CalibdEdx::fitHistGaus(TLinearFitter& fitter, CalibdEdxCorrection& corr, const CalibdEdxCorrection* stackMean) +{ + using timer = std::chrono::high_resolution_clock; + using ax = CalibdEdx::Axis; + auto start = timer::now(); + const bool projSectors = stackMean != nullptr; + constexpr int sectors = SECTORSPERSIDE * SIDES; + for (int iSnp = 0; iSnp < mHist.axis(ax::Snp).size(); ++iSnp) { + for (int iSec = 0; iSec < mHist.axis(ax::Sector).size(); ++iSec) { + for (int iStack = 0; iStack < mHist.axis(ax::Stack).size(); ++iStack) { + for (int iCharge = 0; iCharge < mHist.axis(ax::Charge).size(); ++iCharge) { + + fitter.ClearPoints(); + StackID id{}; + id.type = static_cast(mHist.axis(ax::Stack).bin(iStack).center()); + id.sector = static_cast(mHist.axis(ax::Sector).bin(iSec).center()); + const auto charge = static_cast(mHist.axis(ax::Charge).bin(iCharge).center()); + int entries = 0; + + for (int iTgl = 0; iTgl < mHist.axis(ax::Tgl).size(); ++iTgl) { + // calculate sigma vs tgl in first iteration + // apply cut in n sigma in second iteration + float sigma_vs_tgl = 0; + float mean_vs_tgl = 0; + std::vector bin_indices(ax::Size); + bin_indices[ax::Tgl] = iTgl; + bin_indices[ax::Snp] = iSnp; + bin_indices[ax::Sector] = iSec; + bin_indices[ax::Stack] = iStack; + bin_indices[ax::Charge] = iCharge; + + for (int iter = 0; iter < 2; ++iter) { + auto boostHist1d = ProjectBoostHistoXFast(mHist, bin_indices, ax::dEdx); + + float lowerCut = 0; + float upperCut = 0; + + // make gaussian fit + if (iter == 0) { + int maxElementIndex = std::max_element(boostHist1d.begin(), boostHist1d.end()) - boostHist1d.begin() - 1; + if (maxElementIndex < 0) { + maxElementIndex = 0; + } + float maxElementCenter = 0.5 * (boostHist1d.axis(0).bin(maxElementIndex).upper() + boostHist1d.axis(0).bin(maxElementIndex).lower()); + + lowerCut = (1.f - mFitLowCutFactor * mFitCut) * maxElementCenter; + upperCut = (1.f + mFitCut) * maxElementCenter; + } else { + lowerCut = mean_vs_tgl - sigma_vs_tgl * mSigmaLower; + upperCut = mean_vs_tgl + sigma_vs_tgl * mSigmaUpper; + } + + // Restrict fit range to maximum +- restrictFitRangeToMax + double max_range = mHist.axis(ax::dEdx).bin(mHist.axis(ax::dEdx).size() - 1).lower(); + double min_range = mHist.axis(ax::dEdx).bin(0).lower(); + if ((upperCut <= lowerCut) || (lowerCut > max_range) || (upperCut < min_range)) { + break; + } + + // remove up and low bins + boostHist1d = boost::histogram::algorithm::reduce(boostHist1d, boost::histogram::algorithm::shrink(lowerCut, upperCut)); + + try { + auto fitValues = o2::utils::fitGaus(boostHist1d.begin(), boostHist1d.end(), o2::utils::BinCenterView(boostHist1d.axis(0).begin()), false); + + // add the mean from gaus fit to the fitter + double dEdx = fitValues[1]; + double inputs[] = { + CalibdEdx::recoverTgl(mHist.axis(ax::Tgl).bin(iTgl).center(), id.type), + mHist.axis(ax::Snp).bin(iSnp).center()}; + + // scale fitted dEdx using the stacks mean + if (stackMean != nullptr) { + dEdx /= stackMean->getCorrection(id, charge); + } + + const auto fitNPoints = fitValues[3]; + const float sigma = fitValues[2]; + const double fitMeanErr = (fitNPoints > 0) ? (sigma / std::sqrt(fitNPoints)) : -1; + if (iter == 0) { + sigma_vs_tgl = sigma; + mean_vs_tgl = dEdx; + } else { + entries += fitNPoints; + if (fitMeanErr > 0) { + fitter.AddPoint(inputs, dEdx, fitMeanErr); + } + } + + if (mDebugOutputStreamer) { + const int nbinsx = boostHist1d.axis(0).size(); + std::vector binCenter(nbinsx); + std::vector dedx(nbinsx); + for (int ix = 0; ix < nbinsx; ix++) { + binCenter[ix] = boostHist1d.axis(0).bin(ix).center(); + dedx[ix] = boostHist1d.at(ix); + } + + (*mDebugOutputStreamer) << "fit_gaus" + << "fitConstant=" << fitValues[0] + << "fitMean=" << dEdx + << "fitMeanErr=" << fitMeanErr + << "fitSigma=" << sigma_vs_tgl + << "fitSum=" << fitNPoints + << "dedx=" << binCenter + << "counts=" << dedx + << "itgl=" << bin_indices[1] + << "isnp=" << bin_indices[2] + << "isector=" << bin_indices[3] + << "istack=" << bin_indices[4] + << "icharge=" << bin_indices[5] + << "upperCut=" << upperCut + << "lowerCut=" << lowerCut + << "mFitCut=" << mFitCut + << "mFitLowCutFactor=" << mFitLowCutFactor + << "iter=" << iter + << "mSigmaLower=" << mSigmaLower + << "mSigmaUpper=" << mSigmaUpper + << "\n"; + } + } catch (o2::utils::FitGausError_t err) { + LOGP(warning, "Skipping bin: iTgl {} iSnp {} iSec {} iStack {} iCharge {} iter {}", iTgl, iSnp, iSec, iStack, iCharge, iter); + LOG(warning) << createErrorMessageFitGaus(err); + break; + } + } + } + + const int fitStatus = fitter.Eval(); + if (fitStatus > 0) { + LOGP(warning, "Fit failed"); + continue; + } + + const auto paramSize = CalibdEdxCorrection::ParamSize; + float params[paramSize] = {0}; + for (int param = 0; param < fitter.GetNumberFreeParameters(); ++param) { + params[param] = fitter.GetParameter(param); + } + + // with a projected hist, copy the fit to every sector + if (projSectors) { + for (int i = 0; i < sectors; ++i) { + id.sector = i; + const float mean = stackMean->getCorrection(id, charge); + + // rescale the params to get the true correction + float scaledParams[paramSize]; + for (int i = 0; i < paramSize; ++i) { + scaledParams[i] = params[i] * mean; + } + corr.setParams(id, charge, scaledParams); + corr.setChi2(id, charge, fitter.GetChisquare()); + corr.setEntries(id, charge, entries); + } + } else { + corr.setParams(id, charge, params); + corr.setChi2(id, charge, fitter.GetChisquare()); + corr.setEntries(id, charge, entries); + } + } + } + } + } + auto stop = timer::now(); + std::chrono::duration time = stop - start; + LOGP(info, "Calibration fits took: {}", time.count()); +} + +void CalibdEdx::finalize(const bool useGausFits) { const float entries = minStackEntries(); mCalib.clear(); @@ -295,11 +513,15 @@ void CalibdEdx::finalize() fitter.SetFormula("1"); mCalib.setDims(0); } - LOGP(info, "Fitting {}D dE/dx correction for GEM stacks", mCalib.getDims()); + LOGP(info, "Fitting {}D dE/dx correction for GEM stacks with gaussian fits {}", mCalib.getDims(), useGausFits); // if entries below minimum sector threshold, integrate all sectors if (mCalib.getDims() == 0 || entries >= mSectorThreshold) { - fitHist(mHist, mCalib, fitter, mFitCut, mFitLowCutFactor, mFitPasses); + if (!useGausFits) { + fitHist(mHist, mCalib, fitter, mFitCut, mFitLowCutFactor, mFitPasses, nullptr, mDebugOutputStreamer.get()); + } else { + fitHistGaus(fitter, mCalib, nullptr); + } } else { LOGP(info, "Integrating GEM stacks sectors in dE/dx correction due to low statistics"); @@ -308,10 +530,17 @@ void CalibdEdx::finalize() meanCorr.setDims(0); TLinearFitter meanFitter(0); meanFitter.SetFormula("1"); - fitHist(mHist, meanCorr, meanFitter, mFitCut, mFitLowCutFactor, mFitPasses); + if (!useGausFits) { + fitHist(mHist, meanCorr, meanFitter, mFitCut, mFitLowCutFactor, mFitPasses); - // get higher dimension corrections with projected sectors - fitHist(mHist, mCalib, fitter, mFitCut, mFitLowCutFactor, mFitPasses, &meanCorr); + // get higher dimension corrections with projected sectors + fitHist(mHist, mCalib, fitter, mFitCut, mFitLowCutFactor, mFitPasses, &meanCorr); + } else { + fitHistGaus(meanFitter, meanCorr, nullptr); + + // get higher dimension corrections with projected sectors + fitHistGaus(fitter, mCalib, &meanCorr); + } } } @@ -330,7 +559,7 @@ bool CalibdEdx::hasEnoughData(float minEntries) const return minStackEntries() >= minEntries; } -THnF* CalibdEdx::getRootHist() const +THnF* CalibdEdx::getTHnF() const { std::vector bins{}; std::vector axisMin{}; @@ -346,6 +575,13 @@ THnF* CalibdEdx::getRootHist() const } auto hn = new THnF("hdEdxMIP", "MIP dEdx per GEM stack", histRank, bins.data(), axisMin.data(), axisMax.data()); + return hn; +} + +THnF* CalibdEdx::getRootHist() const +{ + auto hn = getTHnF(); + const size_t histRank = mHist.rank(); std::vector xs(histRank); for (auto&& entry : bh::indexed(mHist)) { if (*entry == 0) { @@ -360,6 +596,60 @@ THnF* CalibdEdx::getRootHist() const return hn; } +void CalibdEdx::setFromRootHist(const THnF* hist) +{ + // Get the number of dimensions + int n_dim = hist->GetNdimensions(); + + // Vectors to store axis ranges and bin counts + std::vector> axis_ranges(n_dim); // Min and max of each axis + std::vector bin_counts(n_dim); // Number of bins in each dimension + + // Loop over each dimension to extract the bin edges and ranges + for (int dim = 0; dim < n_dim; ++dim) { + TAxis* axis = hist->GetAxis(dim); + int bins = axis->GetNbins(); + double min = axis->GetXmin(); + double max = axis->GetXmax(); + bin_counts[dim] = bins; + axis_ranges[dim] = {min, max}; // Store the min and max range for the axis + } + + // Define a Boost histogram using the bin edges + mHist = bh::make_histogram( + FloatAxis(bin_counts[0], axis_ranges[0].first, axis_ranges[0].second, "dEdx"), // dE/dx + FloatAxis(bin_counts[1], axis_ranges[1].first, axis_ranges[1].second, "Tgl"), // Tgl + FloatAxis(bin_counts[2], axis_ranges[2].first, axis_ranges[2].second, "Snp"), // snp + IntAxis(0, bin_counts[3], "sector"), // sector + IntAxis(0, bin_counts[4], "stackType"), // stack type + IntAxis(0, bin_counts[5], "charge") // charge type + ); + + // Fill the Boost histogram with the bin contents from the ROOT histogram + int total_bins = hist->GetNbins(); + for (int bin = 0; bin < total_bins; ++bin) { + std::vector bin_indices(n_dim); + double content = hist->GetBinContent(bin, bin_indices.data()); // Get bin coordinates + + // Check if any coordinate is in the underflow (0) or overflow (nbins+1) bins + bool is_underflow_or_overflow = false; + for (int dim = 0; dim < n_dim; ++dim) { + if ((bin_indices[dim] == 0) || (bin_indices[dim] == bin_counts[dim] + 1)) { + is_underflow_or_overflow = true; + break; + } + } + + // Skip underflow/overflow bins + if (is_underflow_or_overflow) { + continue; + } + + // Assign the content to the corresponding bin in the Boost histogram + mHist.at(bin_indices[0] - 1, bin_indices[1] - 1, bin_indices[2] - 1, bin_indices[3] - 1, bin_indices[4] - 1, bin_indices[5] - 1) = content; + } +} + void CalibdEdx::print() const { const int uniqueEntries = std::accumulate(mHist.begin(), mHist.end(), 0.0) / GEMSTACKSPERSECTOR / 2; @@ -418,3 +708,35 @@ void CalibdEdx::finalizeDebugOutput() const mDebugOutputStreamer->Close(); } } + +void CalibdEdx::dumpToFile(const char* outFile, const char* outName) const +{ + TFile f(outFile, "RECREATE"); + f.WriteObject(this, outName); + const auto* thn = getRootHist(); + f.WriteObject(thn, "histogram_data"); +} + +CalibdEdx CalibdEdx::readFromFile(const char* inFile, const char* inName) +{ + TFile f(inFile, "READ"); + auto* obj = (CalibdEdx*)f.Get(inName); + if (!obj) { + CalibdEdx calTmp; + return calTmp; + } + CalibdEdx cal(*obj); + THnF* hTmp = (THnF*)f.Get("histogram_data"); + if (!hTmp) { + CalibdEdx calTmp; + return calTmp; + } + cal.setFromRootHist(hTmp); + return cal; +} + +void CalibdEdx::setSigmaFitRange(const float lowerSigma, const float upperSigma) +{ + mSigmaUpper = upperSigma; + mSigmaLower = lowerSigma; +} diff --git a/Detectors/TPC/calibration/src/CalibratordEdx.cxx b/Detectors/TPC/calibration/src/CalibratordEdx.cxx index 15a7335ef90d3..7599e2f5d4472 100644 --- a/Detectors/TPC/calibration/src/CalibratordEdx.cxx +++ b/Detectors/TPC/calibration/src/CalibratordEdx.cxx @@ -40,7 +40,7 @@ void CalibratordEdx::finalizeSlot(Slot& slot) // compute calibration values from histograms CalibdEdx* container = slot.getContainer(); - container->finalize(); + container->finalize(mMakeGaussianFits); container->finalizeDebugOutput(); mCalibs.push_back(container->getCalib()); diff --git a/Detectors/TPC/workflow/src/CalibdEdxSpec.cxx b/Detectors/TPC/workflow/src/CalibdEdxSpec.cxx index 046585e8c96b6..a32a4a1bb3089 100644 --- a/Detectors/TPC/workflow/src/CalibdEdxSpec.cxx +++ b/Detectors/TPC/workflow/src/CalibdEdxSpec.cxx @@ -57,6 +57,7 @@ class CalibdEdxDevice : public Task const auto maxdEdx = ic.options().get("max-dedx"); const auto angularBins = ic.options().get("angularbins"); const auto fitSnp = ic.options().get("fit-snp"); + mMakeGaussianFits = !ic.options().get("disable-gaussian-fits"); mDumpToFile = ic.options().get("file-dump"); @@ -103,12 +104,13 @@ class CalibdEdxDevice : public Task void endOfStream(EndOfStreamContext& eos) final { LOGP(info, "Finalizing calibration"); - mCalib->finalize(); + mCalib->finalize(mMakeGaussianFits); mCalib->print(); sendOutput(eos.outputs()); if (mDumpToFile) { - mCalib->getCalib().writeToFile("calibdEdx.root"); + mCalib->dumpToFile("calibdEdx_Obj.root", "calib"); + mCalib->getCalib().writeToFile("calibdEdx.root", "ccdb_object"); if (mDumpToFile > 1) { mCalib->writeTTree("calibdEdx.histo.tree.root"); } @@ -141,6 +143,7 @@ class CalibdEdxDevice : public Task uint64_t mRunNumber{0}; ///< processed run number uint64_t mTimeStampStart{0}; ///< time stamp for first TF for CCDB output std::unique_ptr mCalib; + bool mMakeGaussianFits{true}; ///< make gaussian fits or take the mean }; DataProcessorSpec getCalibdEdxSpec(const o2::base::Propagator::MatCorrType matType) @@ -175,11 +178,12 @@ DataProcessorSpec getCalibdEdxSpec(const o2::base::Propagator::MatCorrType matTy {"fit-threshold", VariantType::Float, 0.2f, {"dEdx width around the MIP peak used in the fit"}}, {"fit-threshold-low-factor", VariantType::Float, 1.5f, {"factor for low dEdx width around the MIP peak used in the fit"}}, - {"dedxbins", VariantType::Int, 60, {"number of dEdx bins"}}, - {"min-dedx", VariantType::Float, 20.0f, {"minimum value for dEdx histograms"}}, + {"dedxbins", VariantType::Int, 70, {"number of dEdx bins"}}, + {"min-dedx", VariantType::Float, 10.0f, {"minimum value for dEdx histograms"}}, {"max-dedx", VariantType::Float, 90.0f, {"maximum value for dEdx histograms"}}, {"angularbins", VariantType::Int, 36, {"number of angular bins: Tgl and Snp"}}, {"fit-snp", VariantType::Bool, false, {"enable Snp correction"}}, + {"disable-gaussian-fits", VariantType::Bool, false, {"disable calibration with gaussian fits and use mean instead"}}, {"file-dump", VariantType::Int, 0, {"directly dump calibration to file"}}}}; } diff --git a/Detectors/TPC/workflow/src/CalibratordEdxSpec.cxx b/Detectors/TPC/workflow/src/CalibratordEdxSpec.cxx index 70b27443018bb..6e477084d992c 100644 --- a/Detectors/TPC/workflow/src/CalibratordEdxSpec.cxx +++ b/Detectors/TPC/workflow/src/CalibratordEdxSpec.cxx @@ -69,6 +69,7 @@ class CalibratordEdxDevice : public Task const auto dumpData = ic.options().get("file-dump"); const auto dumpHistograms = ic.options().get("dump-histograms"); const auto trackDebug = ic.options().get("track-debug"); + const bool makeGaussianFits = !ic.options().get("disable-gaussian-fits"); mCalibrator = std::make_unique(); mCalibrator->setHistParams(dEdxBins, mindEdx, maxdEdx, angularBins, fitSnp); @@ -82,6 +83,7 @@ class CalibratordEdxDevice : public Task mCalibrator->setMaterialType(mMatType); mCalibrator->setDumpHistograms(dumpHistograms); mCalibrator->setTrackDebug(trackDebug); + mCalibrator->setMakeGaussianFits(makeGaussianFits); if (dumpData) { const auto dumpDataName = ic.options().get("file-dump-name"); @@ -218,8 +220,8 @@ DataProcessorSpec getCalibratordEdxSpec(const o2::base::Propagator::MatCorrType {"fit-threshold", VariantType::Float, 0.2f, {"dEdx width around the MIP peak used in the fit"}}, {"fit-threshold-low-factor", VariantType::Float, 1.5f, {"factor for low dEdx width around the MIP peak used in the fit"}}, - {"dedxbins", VariantType::Int, 60, {"number of dEdx bins"}}, - {"min-dedx", VariantType::Float, 20.0f, {"minimum value for dEdx histograms"}}, + {"dedxbins", VariantType::Int, 70, {"number of dEdx bins"}}, + {"min-dedx", VariantType::Float, 10.0f, {"minimum value for dEdx histograms"}}, {"max-dedx", VariantType::Float, 90.0f, {"maximum value for dEdx histograms"}}, {"angularbins", VariantType::Int, 36, {"number of angular bins: Tgl and Snp"}}, {"fit-snp", VariantType::Bool, false, {"enable Snp correction"}}, @@ -228,6 +230,7 @@ DataProcessorSpec getCalibratordEdxSpec(const o2::base::Propagator::MatCorrType {"file-dump", VariantType::Bool, false, {"directly dump calibration to file"}}, {"file-dump-name", VariantType::String, "calibratordEdx.root", {"name of the file dump output file"}}, {"track-debug", VariantType::Bool, false, {"track dEdx debugging"}}, + {"disable-gaussian-fits", VariantType::Bool, false, {"disable calibration with gaussian fits and use mean instead"}}, }}; } diff --git a/Detectors/TPC/workflow/src/MIPTrackFilterSpec.cxx b/Detectors/TPC/workflow/src/MIPTrackFilterSpec.cxx index 7d8d2439e7295..e3970012d1373 100644 --- a/Detectors/TPC/workflow/src/MIPTrackFilterSpec.cxx +++ b/Detectors/TPC/workflow/src/MIPTrackFilterSpec.cxx @@ -163,9 +163,9 @@ DataProcessorSpec getMIPTrackFilterSpec() outputs, adaptFromTask(), Options{ - {"min-momentum", VariantType::Double, 0.3, {"minimum momentum cut"}}, - {"max-momentum", VariantType::Double, 0.7, {"maximum momentum cut"}}, - {"min-dedx", VariantType::Double, 20., {"minimum dEdx cut"}}, + {"min-momentum", VariantType::Double, 0.35, {"minimum momentum cut"}}, + {"max-momentum", VariantType::Double, 0.55, {"maximum momentum cut"}}, + {"min-dedx", VariantType::Double, 10., {"minimum dEdx cut"}}, {"max-dedx", VariantType::Double, 200., {"maximum dEdx cut"}}, {"min-clusters", VariantType::Int, 60, {"minimum number of clusters in a track"}}, {"processEveryNthTF", VariantType::Int, 1, {"Using only a fraction of the data: 1: Use every TF, 10: Process only every tenth TF."}},