Adding TRD t0 fitting procedure

DataFormats/Detectors/TRD/CMakeLists.txt

@@ -14,6 +14,7 @@ o2_add_library(DataFormatsTRD
                        src/LinkRecord.cxx
                        src/AngularResidHistos.cxx
                        src/GainCalibHistos.cxx
+                       src/T0FitHistos.cxx
                        src/DcsCcdbObjects.cxx
                        src/Tracklet64.cxx
                        src/RawData.cxx
@@ -37,6 +38,7 @@ o2_target_root_dictionary(DataFormatsTRD
                        include/DataFormatsTRD/DcsCcdbObjects.h
                        include/DataFormatsTRD/AngularResidHistos.h
                        include/DataFormatsTRD/GainCalibHistos.h
+                       include/DataFormatsTRD/T0FitHistos.h
                        include/DataFormatsTRD/HelperMethods.h
                        include/DataFormatsTRD/Hit.h
                        include/DataFormatsTRD/Digit.h

DataFormats/Detectors/TRD/include/DataFormatsTRD/CalT0.h

@@ -32,11 +32,34 @@ class CalT0
   ~CalT0() = default;
 
   void setT0(int iDet, float t0) { mT0[iDet] = t0; }
+  void setT0av(float t0) { mT0av = t0; }
 
   float getT0(int iDet) const { return mT0[iDet]; }
+  float getT0av() const { return mT0av; }
+  float calcT0av() const
+  {
+    if (mT0.size() == 0) {
+      return -1;
+    }
+    float sum = 0;
+    int counts = 0;
+    for (int iDet = 0; iDet < constants::MAXCHAMBER; ++iDet) {
+      if (mT0[iDet] > -5) {
+        sum += mT0[iDet];
+        ++counts;
+      }
+    }
+
+    if (counts > 0) {
+      return sum / counts;
+    } else {
+      return -2;
+    }
+  }
 
  private:
   std::array<float, constants::MAXCHAMBER> mT0{}; ///< calibrated T0 per TRD chamber
+  float mT0av{-1};                                ///< calibrated average T0
 
   ClassDefNV(CalT0, 1);
 };

DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h

@@ -77,6 +77,7 @@ constexpr double VDRIFTDEFAULT = 1.546; ///< default value for vDrift
 constexpr double EXBDEFAULT = 0.0;      ///< default value for LorentzAngle
 constexpr int NBINSGAINCALIB = 320;     ///< number of bins in the charge (Q0+Q1+Q2) histogram for gain calibration
 constexpr float MPVDEDXDEFAULT = 42.;   ///< default Most Probable Value of TRD dEdx
+constexpr float T0DEFAULT = 1.2;        ///< default value for t0
 
 // array size to store incoming half cru payload.
 constexpr int HBFBUFFERMAX = 1048576;                 ///< max buffer size for data read from a half cru, (all events)

DataFormats/Detectors/TRD/include/DataFormatsTRD/PHData.h

@@ -13,6 +13,7 @@
 #define ALICEO2_TRD_PHDATA_H_
 
 #include <cstdint>
+#include "Rtypes.h"
 
 namespace o2::trd
 {

DataFormats/Detectors/TRD/include/DataFormatsTRD/T0FitHistos.h

@@ -0,0 +1,61 @@
+// 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 T0FitHistos.h
+/// \brief Class to store the TRD PH values for each TRD chamber
+
+#ifndef ALICEO2_T0FITHISTOS_H
+#define ALICEO2_T0FITHISTOS_H
+
+#include "DataFormatsTRD/Constants.h"
+#include "DataFormatsTRD/PHData.h"
+#include "Framework/InputRecord.h"
+#include "Rtypes.h"
+#include <vector>
+#include <memory>
+#include <gsl/span>
+
+namespace o2
+{
+namespace trd
+{
+
+class T0FitHistos
+{
+ public:
+  T0FitHistos() = default;
+  T0FitHistos(const T0FitHistos&) = default;
+  ~T0FitHistos() = default;
+  void init();
+  void reset();
+  auto getDetector(int index) const { return mDet[index]; }
+  auto getTimeBin(int index) const { return mTB[index]; }
+  auto getADC(int index) const { return mADC[index]; }
+  auto getNEntries() const { return mNEntriesTot; }
+
+  void fill(const std::vector<o2::trd::PHData> data);
+  void merge(const T0FitHistos* prev);
+  void print();
+
+ private:
+  std::vector<int> mDet{};
+  std::vector<int> mTB{};
+  std::vector<int> mADC{};
+  size_t mNEntriesTot{0};
+  bool mInitialized{false};
+
+  ClassDefNV(T0FitHistos, 1);
+};
+
+} // namespace trd
+} // namespace o2
+
+#endif // ALICEO2_T0FITHISTOS_H

DataFormats/Detectors/TRD/src/DataFormatsTRDLinkDef.h

@@ -33,6 +33,7 @@
 #pragma link C++ class o2::trd::NoiseStatusMCM + ;
 #pragma link C++ class o2::trd::AngularResidHistos + ;
 #pragma link C++ class o2::trd::GainCalibHistos + ;
+#pragma link C++ class o2::trd::T0FitHistos + ;
 #pragma link C++ class o2::trd::CalVdriftExB + ;
 #pragma link C++ class o2::trd::CalGain + ;
 #pragma link C++ class o2::trd::CalT0 + ;
@@ -51,6 +52,7 @@
 #pragma link C++ class std::vector < o2::trd::Digit> + ;
 #pragma link C++ class std::vector < o2::trd::AngularResidHistos> + ;
 #pragma link C++ class std::vector < o2::trd::GainCalibHistos> + ;
+#pragma link C++ class std::vector < o2::trd::T0FitHistos> + ;
 #pragma link C++ class std::vector < o2::trd::PHData> + ;
 #pragma link C++ class std::vector < o2::trd::KrCluster> + ;
 #pragma link C++ class std::vector < o2::trd::KrClusterTriggerRecord> + ;

DataFormats/Detectors/TRD/src/T0FitHistos.cxx

@@ -0,0 +1,75 @@
+// 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 T0FitHistos.cxx
+/// \brief Class to store the TRD PH values for each chamber
+
+#include "DataFormatsTRD/T0FitHistos.h"
+#include <fairlogger/Logger.h>
+#include <algorithm>
+
+using namespace o2::trd;
+using namespace o2::trd::constants;
+
+void T0FitHistos::init()
+{
+  mDet.resize(0);
+  mTB.resize(0);
+  mADC.resize(0);
+  mInitialized = true;
+}
+
+void T0FitHistos::reset()
+{
+  mDet.resize(0);
+  mTB.resize(0);
+  mADC.resize(0);
+  mNEntriesTot = 0;
+}
+
+void T0FitHistos::fill(const std::vector<o2::trd::PHData> data)
+{
+
+  if (!mInitialized) {
+    init();
+  }
+
+  for (auto ph : data) {
+    int det = ph.getDetector();
+    int tb = ph.getTimebin();
+    int adc = ph.getADC();
+
+    if (ph.getNneighbours() != 2) {
+      continue;
+    }
+
+    mDet.push_back(det);
+    mTB.push_back(tb);
+    mADC.push_back(adc);
+    ++mNEntriesTot;
+  }
+}
+
+void T0FitHistos::merge(const T0FitHistos* prev)
+{
+  auto sizePrev = (int)prev->getNEntries();
+
+  for (int i = 0; i < sizePrev; ++i) {
+    mDet.push_back(prev->getDetector(i));
+    mTB.push_back(prev->getTimeBin(i));
+    mADC.push_back(prev->getADC(i));
+  }
+}
+
+void T0FitHistos::print()
+{
+  LOG(info) << "There are " << mNEntriesTot << " entries in the container";
+}

Detectors/TRD/calibration/CMakeLists.txt

@@ -16,6 +16,7 @@ o2_add_library(TRDCalibration
                        src/CalibratorVdExB.cxx
                        src/CalibratorGain.cxx
                        src/CalibratorNoise.cxx
+                       src/T0Fit.cxx
                        src/PadCalibCCDBBuilder.cxx
                        src/KrClusterFinder.cxx
                        src/DCSProcessor.cxx
@@ -35,6 +36,7 @@ o2_add_library(TRDCalibration
                                    include/TRDCalibration/CalibratorVdExB.h
                                    include/TRDCalibration/CalibratorGain.h
                                    include/TRDCalibration/CalibratorNoise.h
+                                   include/TRDCalibration/T0Fit.h
                                    include/TRDCalibration/CalibrationParams.h
                                    include/TRDCalibration/PadCalibCCDBBuilder.h
                                    include/TRDCalibration/KrClusterFinder.h

Detectors/TRD/calibration/include/TRDCalibration/CalibrationParams.h

@@ -39,6 +39,10 @@ struct TRDCalibParams : public o2::conf::ConfigurableParamHelper<TRDCalibParams>
   float dEdxTPCMin = 30.;
   float dEdxTPCMax = 70.;
 
+  // For t0 fits
+  size_t minEntriesTotalT0Fit = 1'000'000; ///< minimum number of entries in inclusive distribution for (meaningful) t0 fit
+  size_t minEntriesChamberT0Fit = 30'000;  ///< minimum number of entries in one chamber for (meaningful) t0 fit
+
   // Creation of PH plots
   float pileupCut = 0.7;
 

Detectors/TRD/calibration/include/TRDCalibration/T0Fit.h

@@ -0,0 +1,111 @@
+// 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   T0Fit.h
+/// \brief  Fits the TRD PH distributions to extract the t0 value
+/// \author Luisa Bergmann
+
+#ifndef O2_TRD_T0FIT_H
+#define O2_TRD_T0FIT_H
+
+#include "DataFormatsTRD/T0FitHistos.h"
+#include "DetectorsCalibration/TimeSlotCalibration.h"
+#include "DetectorsCalibration/TimeSlot.h"
+#include "DataFormatsTRD/Constants.h"
+#include "CCDB/CcdbObjectInfo.h"
+#include "DataFormatsTRD/CalT0.h"
+#include "TRDCalibration/CalibrationParams.h"
+
+#include "Rtypes.h"
+#include "TH2.h"
+#include "TH1.h"
+#include "TProfile.h"
+#include "TF1.h"
+#include "Fit/Fitter.h"
+#include "TFile.h"
+#include "TTree.h"
+
+#include <array>
+#include <cstdlib>
+#include <memory>
+
+namespace o2
+{
+namespace trd
+{
+//______________________________________________________________________________________________
+struct ErfLandauChi2Functor {
+  double operator()(const double* par) const;
+  std::vector<float> x; ///< x-value (time-bin) of adc profile
+  std::vector<float> y; ///< y-value (av. adc) for corresp. time-bin
+  float lowerBoundFit;  ///< lower bound for fit
+  float upperBoundFit;  ///< upper bound for fit
+};
+
+//______________________________________________________________________________________________
+class T0Fit final : public o2::calibration::TimeSlotCalibration<o2::trd::T0FitHistos>
+{
+  using Slot = o2::calibration::TimeSlot<o2::trd::T0FitHistos>;
+
+ public:
+  T0Fit() = default;
+  ~T0Fit() final = default;
+
+  bool hasEnoughData(const Slot& slot) const final { return slot.getContainer()->getNEntries() >= mMinEntriesTotal; }
+  void initOutput() final;
+  void finalizeSlot(Slot& slot) final;
+  Slot& emplaceNewSlot(bool front, TFType tStart, TFType tEnd) final;
+
+  /// (Re-)Creates a file "trd_t0fit.root". This lets continually fill
+  /// a tree with the fit results.
+  void createOutputFile();
+
+  /// Close the output file. E.g. First write the tree to the file and let the
+  /// smart pointers take care of closing the file.
+  void closeOutputFile();
+
+  const std::vector<o2::trd::CalT0>& getCcdbObjectVector() const { return mObjectVector; }
+  std::vector<o2::ccdb::CcdbObjectInfo>& getCcdbObjectInfoVector() { return mInfoVector; }
+
+  void initProcessing();
+
+  /// Initialize the fit values once with the previous valid ones if they are
+  /// available.
+  void retrievePrev(o2::framework::ProcessingContext& pc);
+
+ private:
+  bool mInitDone{false};                                         ///< flag to avoid creating output etc multiple times
+  const TRDCalibParams& mParams{TRDCalibParams::Instance()};     ///< reference to calibration parameters
+  size_t mMinEntriesTotal{mParams.minEntriesTotalT0Fit};         ///< minimum total number of angular deviations
+  size_t mMinEntriesChamber{mParams.minEntriesChamberT0Fit};     ///< minimum number of angular deviations per chamber for accepting refitted value
+  bool mEnableOutput{false};                                     ///< enable output in a root file instead of the ccdb
+  std::unique_ptr<TFile> mOutFile{nullptr};                      ///< output file
+  std::unique_ptr<TTree> mOutTree{nullptr};                      ///< output tree
+  ErfLandauChi2Functor mFitFunctor;                              ///< used for minimization process, provides chi2 estimate
+  ROOT::Fit::Fitter mFitter;                                     ///< instance of the ROOT fitter
+  std::array<double, 4> mParamsStart;                            ///< Starting parameters for fit
+  std::unique_ptr<TF1> mFuncErfLandau;                           ///< helper function to calculate the t0 value after the fitting procedure
+  std::array<float, o2::trd::constants::MAXCHAMBER> t0_chambers; ///< t0 values of the individual chambers
+  float t0_average{-5};                                          ///< average t0 value across all chambers
+
+  std::vector<o2::ccdb::CcdbObjectInfo> mInfoVector; ///< vector of CCDB infos; each element is filled with CCDB description of accompanying CCDB calibration object
+  std::vector<o2::trd::CalT0> mObjectVector;         ///< vector of CCDB calibration objects; the extracted t0 per chamber and average for given slot
+
+  std::unique_ptr<TProfile> adcProfIncl;                                            ///< profile that holds inclusive PH spectrum
+  std::array<std::unique_ptr<TProfile>, o2::trd::constants::MAXCHAMBER> adcProfDet; ///< array of profiles for PH spectrum of each chamber
+
+  ClassDefNV(T0Fit, 1);
+};
+
+} // namespace trd
+} // namespace o2
+
+#endif // O2_TRD_T0FIT_H