TPC: Add scaling of VDrift with T/P

DataFormats/Detectors/TPC/include/DataFormatsTPC/LtrCalibData.h

@@ -42,6 +42,7 @@ struct LtrCalibData {
   std::vector<uint16_t> matchedLtrIDs; ///< matched laser track IDs
   std::vector<uint16_t> nTrackTF;      ///< number of laser tracks per TF
   std::vector<float> dEdx;             ///< dE/dx of each track
+  float tp{0.f};                       ///< temperature over pressure ratio
 
   bool isValid() const
   {
@@ -138,7 +139,7 @@ struct LtrCalibData {
     dEdx.clear();
   }
 
-  ClassDefNV(LtrCalibData, 4);
+  ClassDefNV(LtrCalibData, 5);
 };
 
 } // namespace o2::tpc

DataFormats/Detectors/TPC/include/DataFormatsTPC/VDriftCorrFact.h

@@ -26,27 +26,35 @@ namespace o2::tpc
 {
 
 struct VDriftCorrFact {
-  long firstTime{};       ///< first time stamp of processed TFs
-  long lastTime{};        ///< last time stamp of processed TFs
-  long creationTime{};    ///< time of creation
-  float corrFact{1.0};    ///< drift velocity correction factor (multiplicative)
-  float corrFactErr{0.0}; ///< stat error of correction factor
-  float refVDrift{0.};    ///< reference vdrift for which factor was extracted
+  long firstTime{};         ///< first time stamp of processed TFs
+  long lastTime{};          ///< last time stamp of processed TFs
+  long creationTime{};      ///< time of creation
+  float corrFact{1.0};      ///< drift velocity correction factor (multiplicative)
+  float corrFactErr{0.0};   ///< stat error of correction factor
+  float refVDrift{0.};      ///< reference vdrift for which factor was extracted
   float refTimeOffset{0.};  ///< additive time offset reference (\mus)
   float timeOffsetCorr{0.}; ///< additive time offset correction (\mus)
+  float refTP{0.};          ///< reference temperature / pressure for which refVDrift was extracted
 
   float getVDrift() const { return refVDrift * corrFact; }
   float getVDriftError() const { return refVDrift * corrFactErr; }
 
   float getTimeOffset() const { return refTimeOffset + timeOffsetCorr; }
 
   // renormalize VDrift reference and correction either to provided new reference (if >0) or to correction 1 wrt current reference
-  void normalize(float newVRef = 0.f)
+  void normalize(float newVRef = 0.f, float tp = 0.f)
   {
+    float normVDrift = newVRef;
     if (newVRef == 0.f) {
-      newVRef = refVDrift * corrFact;
+      normVDrift = refVDrift * corrFact;
+      newVRef = normVDrift;
+      if ((tp > 0) && (refTP > 0)) {
+        // linear scaling based on relative change of T/P
+        normVDrift *= refTP / tp;
+        refTP = tp; // update reference T/P
+      }
     }
-    float fact = refVDrift / newVRef;
+    float fact = refVDrift / normVDrift;
     refVDrift = newVRef;
     corrFactErr *= fact;
     corrFact *= fact;
@@ -66,7 +74,7 @@ struct VDriftCorrFact {
     }
   }
 
-  ClassDefNV(VDriftCorrFact, 2);
+  ClassDefNV(VDriftCorrFact, 3);
 };
 
 } // namespace o2::tpc

Detectors/Calibration/include/DetectorsCalibration/IntegratedClusterCalibrator.h

@@ -331,6 +331,9 @@ struct TimeSeriesdEdx {
 };
 
 struct TimeSeriesITSTPC {
+  float mVDrift = 0;                          ///< drift velocity in cm/us
+  float mPressure = 0;                        ///< pressure
+  float mTemperature = 0;                     ///< temperature
   TimeSeries mTSTPC;                          ///< TPC standalone DCAs
   TimeSeries mTSITSTPC;                       ///< ITS-TPC standalone DCAs
   ITSTPC_Matching mITSTPCAll;                 ///< ITS-TPC matching efficiency for ITS standalone + afterburner
@@ -499,7 +502,7 @@ struct TimeSeriesITSTPC {
     nVertexContributors_Quantiles.resize(nTotalQ);
   }
 
-  ClassDefNV(TimeSeriesITSTPC, 5);
+  ClassDefNV(TimeSeriesITSTPC, 6);
 };
 
 } // end namespace tpc

Detectors/GlobalTracking/src/MatchTPCITS.cxx

@@ -101,7 +101,7 @@ void MatchTPCITS::run(const o2::globaltracking::RecoContainer& inp,
       break;
     }
     if (mVDriftCalibOn) { // in the beginning of the output vector we send the full and reference VDrift used for this TF
-      calib.emplace_back(mTPCVDrift, mTPCDrift.refVDrift, -999.);
+      calib.emplace_back(mTPCVDrift, mTPCDrift.refVDrift, mTPCDrift.refTP);
       calib.emplace_back(mTPCDriftTimeOffset, mTPCDrift.refTimeOffset, -999.);
     }
 

Detectors/TPC/calibration/include/TPCCalibration/CalibLaserTracks.h

@@ -25,7 +25,6 @@
 #include <gsl/span>
 #include <string_view>
 
-#include "CommonConstants/MathConstants.h"
 #include "CommonUtils/TreeStreamRedirector.h"
 #include "DataFormatsTPC/TrackTPC.h"
 #include "DataFormatsTPC/LaserTrack.h"
@@ -74,10 +73,12 @@ class CalibLaserTracks
   ~CalibLaserTracks() = default;
 
   /// process all tracks of one TF
-  void fill(const gsl::span<const TrackTPC> tracks);
+  /// \param tp ratio of temperature over pressure
+  void fill(const gsl::span<const TrackTPC> tracks, float tp = 0);
 
   /// process all tracks of one TF
-  void fill(std::vector<TrackTPC> const& tracks);
+  /// \param tp ratio of temperature over pressure
+  void fill(std::vector<TrackTPC> const& tracks, float tp = 0);
 
   /// process single track
   void processTrack(const TrackTPC& track);
@@ -163,6 +164,8 @@ class CalibLaserTracks
   float mDriftV{0};                                            ///< drift velocity used during reconstruction
   float mTOffsetMUS{0};                                        ///< time offset in \mus to impose
   float mZbinWidth{0};                                         ///< width of a bin in us
+  float mAvgTP{0};                                             ///< ratio of average temperature over pressure
+  float mAvgDriftV{0};                                         ///< average drift velocity used for the laser track calibration
   uint64_t mTFstart{0};                                        ///< start time of processed time frames
   uint64_t mTFend{0};                                          ///< end time of processed time frames
   LtrCalibData mCalibDataTF{};                                 ///< calibration data for single TF (debugging)
@@ -184,7 +187,7 @@ class CalibLaserTracks
   /// perform fits on the matched z-position pairs to extract the drift velocity correction factor and trigger offset
   void fillCalibData(LtrCalibData& calibData, const std::vector<TimePair>& pairsA, const std::vector<TimePair>& pairsC);
 
-  ClassDefNV(CalibLaserTracks, 1);
+  ClassDefNV(CalibLaserTracks, 2);
 };
 
 } // namespace o2::tpc

Detectors/TPC/calibration/include/TPCCalibration/PressureTemperatureHelper.h

@@ -63,15 +63,35 @@ class PressureTemperatureHelper
   /// get pressure for given time stamp in ms
   float getPressure(const ULong64_t timestamp) const { return interpolate(mPressure.second, mPressure.first, timestamp); }
 
+  /// manually set the pressure
+  void setPressure(const std::pair<std::vector<float>, std::vector<ULong64_t>>& pressure) { mPressure = pressure; }
+
+  /// manually set the temperature
+  void setTemperature(const std::pair<std::vector<float>, std::vector<ULong64_t>>& temperatureA, const std::pair<std::vector<float>, std::vector<ULong64_t>>& temperatureC)
+  {
+    mTemperatureA = temperatureA;
+    mTemperatureC = temperatureC;
+  }
+
   /// get temperature for given time stamp in ms
   dataformats::Pair<float, float> getTemperature(const ULong64_t timestamp) const { return dataformats::Pair<float, float>{interpolate(mTemperatureA.second, mTemperatureA.first, timestamp), interpolate(mTemperatureC.second, mTemperatureC.first, timestamp)}; }
 
+  /// get mean temperature over A and C side
+  float getMeanTemperature(const ULong64_t timestamp) const;
+
+  // get ratio of temperature over pressure for given time stamp
+  float getTP(int64_t ts) const;
+
   static constexpr o2::header::DataDescription getDataDescriptionPressure() { return o2::header::DataDescription{"pressure"}; }
   static constexpr o2::header::DataDescription getDataDescriptionTemperature() { return o2::header::DataDescription{"temperature"}; }
 
+  /// get minimum and maximum time stamps of the pressure and temperature data
+  std::pair<ULong64_t, ULong64_t> getMinMaxTime() const;
+
  protected:
   static void addInput(std::vector<o2::framework::InputSpec>& inputs, o2::framework::InputSpec&& isp);
   static void addOutput(std::vector<o2::framework::OutputSpec>& outputs, o2::framework::OutputSpec&& osp);
+  static constexpr float toKelvin(float celsius) { return celsius + 273.15f; } // convert Celsius to Kelvin
 
   std::pair<std::vector<float>, std::vector<ULong64_t>> mPressure;     ///< pressure values for both measurements
   std::pair<std::vector<float>, std::vector<ULong64_t>> mTemperatureA; ///< temperature values A-side

Detectors/TPC/calibration/include/TPCCalibration/TPCVDriftTglCalibration.h

@@ -32,6 +32,7 @@ struct TPCVDTglContainer {
   double driftVFullMean = 0.;
   static float tOffsetRef;
   static float driftVRef;
+  float tp = 0;
 
   TPCVDTglContainer(int ntgl, float tglMax, int ndtgl, float dtglMax)
   {
@@ -42,9 +43,12 @@ struct TPCVDTglContainer {
   {
     histo = std::make_unique<o2::dataformats::FlatHisto2D_f>(*(src.histo.get()));
     entries = src.entries;
+    tp = src.tp;
+    driftVFullMean = src.driftVFullMean;
   }
 
-  void fill(const gsl::span<const o2::dataformats::Triplet<float, float, float>> data)
+  /// \param tp ratio of temperature over pressure
+  void fill(const gsl::span<const o2::dataformats::Triplet<float, float, float>> data, float currentTemperaturePressure = 0)
   {
     if (data.size() < 3) { // first 2 entres always contains the {full and reference VDrift} and {full and reference DriftTimeOffset} used for the TF
       return;
@@ -59,12 +63,14 @@ struct TPCVDTglContainer {
     }
     //
     float vfull = data[0].first, vref = data[0].second;
+    const float temperaturePressure = (data[0].third == 0) ? currentTemperaturePressure : data[0].third;
     if (driftVRef == 0.f) {
       driftVRef = vref;
     } else if (driftVRef != vref) {
       LOGP(warn, "data with VDriftRef={} were received while initially was set to {}, keep old one", vref, driftVRef);
     }
     driftVFullMean = (driftVFullMean * nTFProc + vfull) / (nTFProc + 1);
+    tp = (tp * nTFProc + temperaturePressure) / (nTFProc + 1);
     if (tOffsetRef == 0.f) {
       tOffsetRef = data[1].first; // assign 1st full toffset as a reference
     }
@@ -73,6 +79,11 @@ struct TPCVDTglContainer {
 
   void merge(const TPCVDTglContainer* other)
   {
+    const int norm = nTFProc + other->nTFProc;
+    if (norm > 0) {
+      tp = (tp * nTFProc + other->tp * other->nTFProc) / norm;
+      driftVFullMean = (driftVFullMean * nTFProc + other->driftVFullMean * other->nTFProc) / norm;
+    }
     entries += other->entries;
     histo->add(*(other->histo));
     LOGP(debug, "Old entries:{} New entries:{} oldSum: {} newSum: {}", other->entries, entries, other->histo->getSum(), histo->getSum());
@@ -82,7 +93,7 @@ struct TPCVDTglContainer {
   {
     LOG(info) << "Nentries = " << entries;
   }
-  ClassDefNV(TPCVDTglContainer, 1);
+  ClassDefNV(TPCVDTglContainer, 2);
 };
 
 class TPCVDriftTglCalibration : public o2::calibration::TimeSlotCalibration<TPCVDTglContainer>

Detectors/TPC/calibration/include/TPCCalibration/VDriftHelper.h

@@ -18,6 +18,7 @@
 
 #include "GPUCommonRtypes.h"
 #include "DataFormatsTPC/VDriftCorrFact.h"
+#include "TPCCalibration/PressureTemperatureHelper.h"
 #include <array>
 #include <vector>
 #include <string_view>
@@ -56,23 +57,28 @@ class VDriftHelper
   Source getSource() const { return mSource; }
   static std::string_view getSourceName(Source s) { return SourceNames[s]; }
   std::string_view getSourceName() const { return SourceNames[mSource]; }
+  const auto& getPTHelper() const { return mPTHelper; }
 
   bool accountCCDBInputs(const o2::framework::ConcreteDataMatcher& matcher, void* obj);
   void extractCCDBInputs(o2::framework::ProcessingContext& pc, bool laser = true, bool itstpcTgl = true);
   static void requestCCDBInputs(std::vector<o2::framework::InputSpec>& inputs, bool laser = true, bool itstpcTgl = true);
 
  protected:
   static void addInput(std::vector<o2::framework::InputSpec>& inputs, o2::framework::InputSpec&& isp);
+  bool extractTPForVDrift(VDriftCorrFact& vdrift, int64_t tsStepMS = 100 * 1000);
   VDriftCorrFact mVDLaser{};
   VDriftCorrFact mVDTPCITSTgl{};
   VDriftCorrFact mVD{};
-  Source mSource{Source::Param}; // update source
-  bool mUpdated = false;  // signal update, must be reset once new value is fetched
+  Source mSource{Source::Param};       // update source
+  bool mUpdated = false;               // signal update, must be reset once new value is fetched
+  bool mIsTPScalingPossible = false;   // if T/P scaling is possible always perform the updating
   bool mForceParamDrift = false;       // enforce vdrift from gasParam
   bool mForceParamOffset = false;      // enforce offset from DetectorParam
+  bool mForceTPScaling = false;        // enforce T/P scaling from gasParam (scaling disabled by negative T or P)
   uint32_t mMayRenormSrc = 0xffffffff; // if starting VDrift correction != 1, we will renorm reference in such a way that initial correction is 1.0, flag per source
+  PressureTemperatureHelper mPTHelper; // helper to extract pressure and temperature from CCDB
 
-  ClassDefNV(VDriftHelper, 1);
+  ClassDefNV(VDriftHelper, 2);
 };
 } // namespace o2::tpc
 #endif

Detectors/TPC/calibration/src/CalibLaserTracks.cxx

@@ -24,13 +24,13 @@
 #include <string_view>
 
 using namespace o2::tpc;
-void CalibLaserTracks::fill(std::vector<TrackTPC> const& tracks)
+void CalibLaserTracks::fill(std::vector<TrackTPC> const& tracks, float tp)
 {
-  fill(gsl::span(tracks.data(), tracks.size()));
+  fill(gsl::span(tracks.data(), tracks.size()), tp);
 }
 
 //______________________________________________________________________________
-void CalibLaserTracks::fill(const gsl::span<const TrackTPC> tracks)
+void CalibLaserTracks::fill(const gsl::span<const TrackTPC> tracks, float tp)
 {
   // ===| clean up TF data |===
   mZmatchPairsTFA.clear();
@@ -63,6 +63,9 @@ void CalibLaserTracks::fill(const gsl::span<const TrackTPC> tracks)
   mCalibDataTF.firstTime = mTFstart;
   mCalibDataTF.lastTime = tfEnd;
 
+  mAvgTP = (mAvgTP * mCalibData.processedTFs + tp) / (mCalibData.processedTFs + 1);
+  mAvgDriftV = (mAvgDriftV * mCalibData.processedTFs + mDriftV) / (mCalibData.processedTFs + 1);
+
   // ===| TF counters |===
   ++mCalibData.processedTFs;
   ++mCalibDataTF.processedTFs;
@@ -147,6 +150,8 @@ void CalibLaserTracks::processTrack(const TrackTPC& track)
                   << "ltr=" << ltr              // matched ideal laser track
                   << "trOutLtr=" << parOutAtLtr // track rotated and propagated to ideal track position
                   << "TPCTracks=" << writeTrack // original TPC track
+                  << "mDriftV=" << mDriftV
+                  << "laserTrackID=" << laserTrackID
                   << "\n";
   }
 }
@@ -277,6 +282,18 @@ void CalibLaserTracks::merge(const CalibLaserTracks* other)
   mCalibData.firstTime = std::min(mCalibData.firstTime, other->mCalibData.firstTime);
   mCalibData.lastTime = std::max(mCalibData.lastTime, other->mCalibData.lastTime);
 
+  if ((mAvgTP > 0) && (other->mAvgTP > 0)) {
+    mAvgTP = (mAvgTP + other->mAvgTP) / 2.0;
+  } else if (other->mAvgTP > 0) {
+    mAvgTP = other->mAvgTP;
+  }
+
+  if ((mAvgDriftV > 0) && (other->mAvgDriftV > 0)) {
+    mAvgDriftV = (mAvgDriftV + other->mAvgDriftV) / 2.0;
+  } else if (other->mAvgDriftV > 0) {
+    mAvgDriftV = other->mAvgDriftV;
+  }
+
   sort(mZmatchPairsA);
   sort(mZmatchPairsC);
 
@@ -296,6 +313,7 @@ void CalibLaserTracks::endTF()
                     << "zPairsA=" << mZmatchPairsTFA
                     << "zPairsC=" << mZmatchPairsTFC
                     << "calibData=" << mCalibDataTF
+                    << "mDriftV=" << mDriftV
                     << "\n";
   }
 }
@@ -330,7 +348,7 @@ void CalibLaserTracks::fillCalibData(LtrCalibData& calibData, const std::vector<
   auto dvA = fit(pairsA, "A-Side");
   auto dvC = fit(pairsC, "C-Side");
   calibData.creationTime = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
-  calibData.refVDrift = mDriftV;
+  calibData.refVDrift = mAvgDriftV;
   calibData.dvOffsetA = dvA.x1;
   calibData.dvCorrectionA = dvA.x2;
   calibData.nTracksA = uint16_t(pairsA.size());
@@ -340,6 +358,7 @@ void CalibLaserTracks::fillCalibData(LtrCalibData& calibData, const std::vector<
   calibData.nTracksC = uint16_t(pairsC.size());
 
   calibData.refTimeOffset = mTOffsetMUS;
+  calibData.tp = mAvgTP;
 }
 
 //______________________________________________________________________________