TRD: updates in vdrift and ExB calibration + possibility to use slope in chi2 matching

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

@@ -33,12 +33,52 @@ class CalGain
 
   void setMPVdEdx(int iDet, float mpv) { mMPVdEdx[iDet] = mpv; }
 
-  float getMPVdEdx(int iDet) const { return mMPVdEdx[iDet]; }
+  float getMPVdEdx(int iDet, bool defaultAvg = true) const
+  {
+    // if defaultAvg = false, we take the value stored whatever it is
+    // if defaultAvg = true and we have default value or bad value stored, we take the average on all chambers instead
+    if (!defaultAvg || isGoodGain(iDet))
+      return mMPVdEdx[iDet];
+    else {
+      if (TMath::Abs(mMeanGain + 999.) < 1e-6)
+        mMeanGain = getAverageGain();
+      return mMeanGain;
+    }
+  }
+
+  float getAverageGain() const
+  {
+    float averageGain = 0.;
+    int ngood = 0;
+
+    for (int iDet = 0; iDet < constants::MAXCHAMBER; iDet++) {
+      if (isGoodGain(iDet)) {
+        // The chamber has correct calibration
+        ngood++;
+        averageGain += mMPVdEdx[iDet];
+      }
+    }
+    if (ngood == 0) {
+      // we should make sure it never happens
+      return constants::MPVDEDXDEFAULT;
+    }
+    averageGain /= ngood;
+    return averageGain;
+  }
+
+  bool isGoodGain(int iDet) const
+  {
+    if (TMath::Abs(mMPVdEdx[iDet] - constants::MPVDEDXDEFAULT) > 1e-6)
+      return true;
+    else
+      return false;
+  }
 
  private:
   std::array<float, constants::MAXCHAMBER> mMPVdEdx{}; ///< Most probable value of dEdx distribution per TRD chamber
+  mutable float mMeanGain{-999.};                      ///! average gain, calculated only once
 
-  ClassDefNV(CalGain, 1);
+  ClassDefNV(CalGain, 2);
 };
 
 } // namespace trd

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

@@ -34,14 +34,102 @@ class CalVdriftExB
   void setVdrift(int iDet, float vd) { mVdrift[iDet] = vd; }
   void setExB(int iDet, float exb) { mExB[iDet] = exb; }
 
-  float getVdrift(int iDet) const { return mVdrift[iDet]; }
-  float getExB(int iDet) const { return mExB[iDet]; }
+  float getVdrift(int iDet, bool defaultAvg = true) const
+  {
+    // if defaultAvg = false, we take the value stored whatever it is
+    // if defaultAvg = true and we have default value or bad value stored, we take the average on all chambers instead
+    if (!defaultAvg || (isGoodExB(iDet) && isGoodVdrift(iDet)))
+      return mVdrift[iDet];
+    else {
+      if (TMath::Abs(mMeanVdrift + 999.) < 1e-6)
+        mMeanVdrift = getAverageVdrift();
+      return mMeanVdrift;
+    }
+  }
+  float getExB(int iDet, bool defaultAvg = true) const
+  {
+    if (!defaultAvg || (isGoodExB(iDet) && isGoodVdrift(iDet)))
+      return mExB[iDet];
+    else {
+      if (TMath::Abs(mMeanExB + 999.) < 1e-6)
+        mMeanExB = getAverageExB();
+      return mMeanExB;
+    }
+  }
+
+  float getAverageVdrift() const
+  {
+    float averageVdrift = 0.;
+    int ngood = 0;
+
+    for (int iDet = 0; iDet < constants::MAXCHAMBER; iDet++) {
+      if (isGoodExB(iDet) && isGoodVdrift(iDet)) {
+        // Both values need to be correct to declare a chamber as well calibrated
+        ngood++;
+        averageVdrift += mVdrift[iDet];
+      }
+    }
+    if (ngood == 0) {
+      // we should make sure it never happens
+      return constants::VDRIFTDEFAULT;
+    }
+    averageVdrift /= ngood;
+    return averageVdrift;
+  }
+
+  float getAverageExB() const
+  {
+    float averageExB = 0.;
+    int ngood = 0;
+
+    for (int iDet = 0; iDet < constants::MAXCHAMBER; iDet++) {
+      if (isGoodExB(iDet) && isGoodVdrift(iDet)) {
+        // Both values need to be correct to declare a chamber as well calibrated
+        ngood++;
+        averageExB += mExB[iDet];
+      }
+    }
+    if (ngood == 0) {
+      // we should make sure it never happens
+      return constants::EXBDEFAULT;
+    }
+    averageExB /= ngood;
+    return averageExB;
+  }
+
+  bool isGoodExB(int iDet) const
+  {
+    // check if value is well calibrated or not
+    // default calibration if not enough entries
+    // close to boundaries indicate a failed fit
+    if (TMath::Abs(mExB[iDet] - constants::EXBDEFAULT) > 1e-6 &&
+        TMath::Abs(mExB[iDet] - constants::EXBMIN) > 0.01 &&
+        TMath::Abs(mExB[iDet] - constants::EXBMAX) > 0.01)
+      return true;
+    else
+      return false;
+  }
+
+  bool isGoodVdrift(int iDet) const
+  {
+    // check if value is well calibrated or not
+    // default calibration if not enough entries
+    // close to boundaries indicate a failed fit
+    if (TMath::Abs(mVdrift[iDet] - constants::VDRIFTDEFAULT) > 1e-6 &&
+        TMath::Abs(mVdrift[iDet] - constants::VDRIFTMIN) > 0.1 &&
+        TMath::Abs(mVdrift[iDet] - constants::VDRIFTMAX) > 0.1)
+      return true;
+    else
+      return false;
+  }
 
  private:
   std::array<float, constants::MAXCHAMBER> mVdrift{}; ///< calibrated drift velocity per TRD chamber
   std::array<float, constants::MAXCHAMBER> mExB{};    ///< calibrated Lorentz angle per TRD chamber
+  mutable float mMeanVdrift{-999.};                   ///! average drift velocity, calculated only once
+  mutable float mMeanExB{-999.};                      ///! average lorentz angle, calculated only once
 
-  ClassDefNV(CalVdriftExB, 1);
+  ClassDefNV(CalVdriftExB, 2);
 };
 
 } // namespace trd

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

@@ -75,7 +75,11 @@ constexpr int TIMEBINS = 30;            ///< the number of time bins
 constexpr float MAXIMPACTANGLE = 25.f;  ///< the maximum impact angle for tracks relative to the TRD detector plane to be considered for vDrift and ExB calibration
 constexpr int NBINSANGLEDIFF = 25;      ///< the number of bins for the track angle used for the vDrift and ExB calibration based on the tracking
 constexpr double VDRIFTDEFAULT = 1.546; ///< default value for vDrift
+constexpr double VDRIFTMIN = 0.4;       ///< min value for vDrift
+constexpr double VDRIFTMAX = 2.0;       ///< max value for vDrift
 constexpr double EXBDEFAULT = 0.0;      ///< default value for LorentzAngle
+constexpr double EXBMIN = -0.4;         ///< min value for LorentzAngle
+constexpr double EXBMAX = 0.4;          ///< max 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

Detectors/TRD/base/src/TrackletTransformer.cxx

@@ -50,7 +50,8 @@ float TrackletTransformer::calculateDy(int detector, int slope, const PadPlane*
   // NOTE: check what drift height is used in calibration code to ensure consistency
   // NOTE: check sign convention of Lorentz angle
   // NOTE: confirm the direction in which vDrift is measured/determined. Is it in x or in direction of drift?
-  double lorentzCorrection = TMath::Tan(exb) * mXAnode;
+  // The Lorentz correction have to be applied both at the point of entrance and at the end of the drift region
+  double lorentzCorrection = TMath::Tan(exb) * mGeo->cdrHght();
 
   // assuming angle in Bailhache, fig. 4.17 would be positive in our calibration code
   double calibratedDy = rawDy - lorentzCorrection;

Detectors/TRD/calibration/README.md

@@ -27,15 +27,15 @@ For 'o2-calibration-trd-workflow --vDriftAndExB' there are also the following ke
 
 *Hint: You can get information on the meaning of the parameters by running `o2-calibration-trd-workflow --vDriftAndExB -b --help full`*
 
-If you want to run the calibration from a local file with residuals, trdangreshistos.root, you can run:
+If you want to run the calibration from a local file with residuals, trdcaliboutput.root, you can run:
 
-    o2-calibration-trd-workflow --vDriftAndExB -b --enable-root-input --calib-vdexb-calibration '--tf-per-slot 1' --configKeyValues "TRDCalibParams.minEntriesChamber=100;TRDCalibParams.minEntriesTotal=50000"
+    o2-calibration-trd-workflow --vDriftAndExB -b --enable-root-input --calib-vdexb-calibration '--tf-per-slot 1' --configKeyValues "TRDCalibParams.minEntriesChamber=100;TRDCalibParams.minEntriesTotal=50000" --trd-calib-infile trdcaliboutput.root
 
 Additionally it is possible to perform the calibrations fit manually per chamber if you have TPC-TRD or ITS-TPC-TRD tracks, you can run:
 
     o2-trd-global-tracking -b --enable-trackbased-calib
 
-This produces `trdangreshistos.root` which holds the residuals of the angles and differences.
+This produces `trdcaliboutput.root` which holds the residuals of the angles and differences.
 Then run the macro `Detectors/TRD/calibration/macros/manualCalibFit.C`.
 This produces a file of similar name with the fitted data and prints out the fit results.
 This is equivalent to running:

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

@@ -25,12 +25,13 @@ namespace trd
 /// VDrift and ExB calibration parameters.
 struct TRDCalibParams : public o2::conf::ConfigurableParamHelper<TRDCalibParams> {
   unsigned int nTrackletsMin = 5;  ///< minimum amount of tracklets
+  unsigned int nTrackletsMinLoose = 4; ///< minimum amount of tracklets if two layers with a large lever arm both have a hit
   unsigned int chi2RedMax = 6;     ///< maximum reduced chi2 acceptable for track quality
-  size_t minEntriesChamber = 75;   ///< minimum number of entries per chamber to fit single time slot
-  size_t minEntriesTotal = 40'500; ///< minimum total required for meaningful fits
+  size_t minEntriesChamber = 200;  ///< minimum number of entries per chamber to fit single time slot
+  size_t minEntriesTotal = 400'000; ///< minimum total required for meaningful fits
 
   // For gain calibration
-  unsigned int nTrackletsMinGainCalib = 5;
+  unsigned int nTrackletsMinGainCalib = 3;
   size_t minEntriesChamberGainCalib = 500;     ///< minimum number of entries per chamber to fit single time slot
   size_t minEntriesTotalGainCalib = 1'000'000; ///< minimum total required for meaningful fits
   // Cuts for selecting clean pion candidates for gain calibration

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

@@ -88,8 +88,8 @@ class CalibratorVdExB final : public o2::calibration::TimeSlotCalibration<o2::tr
  private:
   bool mInitDone{false};                                     ///< flag to avoid creating the TProfiles multiple times
   const TRDCalibParams& mParams{TRDCalibParams::Instance()}; ///< reference to calibration parameters
-  size_t mMinEntriesTotal{mParams.minEntriesChamber};        ///< minimum total number of angular deviations (on average ~3 entries per bin for each TRD chamber)
-  size_t mMinEntriesChamber{mParams.minEntriesTotal};        ///< minimum number of angular deviations per chamber for accepting refitted value (~3 per bin)
+  size_t mMinEntriesTotal{mParams.minEntriesTotal};          ///< minimum total number of angular deviations (on average ~3 entries per bin for each TRD chamber)
+  size_t mMinEntriesChamber{mParams.minEntriesChamber};      ///< minimum number of angular deviations per chamber for accepting refitted value (~3 per bin)
   bool mEnableOutput{false};                                 ///< enable output of calibration fits and tprofiles in a root file instead of the ccdb
   std::unique_ptr<TFile> mOutFile{nullptr};                  ///< output file
   std::unique_ptr<TTree> mOutTree{nullptr};                  ///< output tree

Detectors/TRD/calibration/macros/manualCalibFit.C

@@ -30,18 +30,19 @@
 
 // O2 header
 #include <TRDCalibration/CalibratorVdExB.h>
+#include "DetectorsBase/Propagator.h"
 
 #endif
 
 // This root macro reads in 'trdangreshistos.root' and
 // performs the calibration fits manually as in CalibratorVdExB.cxx
 // This can be used for checking if the calibration fits make sense.
-void manualCalibFit()
+void manualCalibFit(int runNumber = 563335, bool usePreCorrFromCCDB = false)
 {
   //----------------------------------------------------
   // TTree and File
   //----------------------------------------------------
-  std::unique_ptr<TFile> inFilePtr(TFile::Open("trdangreshistos.root"));
+  std::unique_ptr<TFile> inFilePtr(TFile::Open("trdcaliboutput.root"));
   if (inFilePtr == nullptr) {
     printf("Input File could not be read!\n'");
     return;
@@ -60,18 +61,46 @@ void manualCalibFit()
   tree->SetBranchAddress("mHistogramEntries[13500]", &mHistogramEntries);
   tree->SetBranchAddress("mNEntriesPerBin[13500]", &mNEntriesPerBin);
 
+  // use precorr values from ccdb
+  // necessary when the angular residuals were calculated already using ccdb calibration (e.g. in a local run)
+
+  o2::trd::CalVdriftExB* calObject;
+  if (usePreCorrFromCCDB) {
+    auto& ccdbmgr = o2::ccdb::BasicCCDBManager::instance();
+
+    o2::ccdb::CcdbApi ccdb;
+    ccdb.init("http://alice-ccdb.cern.ch");
+    auto runDuration = ccdbmgr.getRunDuration(runNumber);
+
+    std::map<std::string, std::string> metadata;
+    std::map<std::string, std::string> headers;
+
+    calObject = ccdb.retrieveFromTFileAny<o2::trd::CalVdriftExB>("TRD/Calib/CalVdriftExB", metadata, runDuration.first + 60000, &headers, "", "", "1689478811721");
+  }
+
   //----------------------------------------------------
   // Configure Fitter
   //----------------------------------------------------
   o2::trd::FitFunctor mFitFunctor;
   std::array<std::unique_ptr<TProfile>, 540> profiles; ///< profile histograms for each TRD chamber
+  int counter = 0;
   for (int iDet = 0; iDet < 540; ++iDet) {
     mFitFunctor.profiles[iDet] = std::make_unique<TProfile>(Form("profAngleDiff_%i", iDet), Form("profAngleDiff_%i", iDet), 25, -25.f, 25.f);
+    if (usePreCorrFromCCDB) {
+      if (calObject->isGoodExB(iDet))
+        counter++;
+      mFitFunctor.vdPreCorr[iDet] = calObject->getVdrift(iDet, true);
+      mFitFunctor.laPreCorr[iDet] = calObject->getExB(iDet, true);
+    }
   }
+  std::cout << counter << " good entries in the CCDB " << std::endl;
+  mFitFunctor.mAnodePlane = 3.35; // don't really care as long as it's not zero, this parameter could  be removed
   mFitFunctor.lowerBoundAngleFit = 80 * TMath::DegToRad();
   mFitFunctor.upperBoundAngleFit = 100 * TMath::DegToRad();
-  mFitFunctor.vdPreCorr.fill(1.546);
-  mFitFunctor.laPreCorr.fill(0.0);
+  if (!usePreCorrFromCCDB) {
+    mFitFunctor.vdPreCorr.fill(1.546);
+    mFitFunctor.laPreCorr.fill(0.0);
+  }
 
   //----------------------------------------------------
   // Loop
@@ -88,15 +117,18 @@ void manualCalibFit()
   //----------------------------------------------------
   // Fill profiles
   //----------------------------------------------------
+  int nEntriesDetTotal[540] = {};
   for (int iDet = 0; iDet < 540; ++iDet) {
     for (int iBin = 0; iBin < 25; ++iBin) {
       auto angleDiffSum = mHistogramEntriesSum[iDet * 25 + iBin];
       auto nEntries = mNEntriesPerBinSum[iDet * 25 + iBin];
+      nEntriesDetTotal[iDet] += nEntries;
       if (nEntries > 0) { // skip entries which have no entries; ?
         // add to the respective profile for fitting later on
         mFitFunctor.profiles[iDet]->Fill(2 * iBin - 25.f, angleDiffSum / nEntries, nEntries);
       }
     }
+    printf("Det %d: nEntries=%d \n", iDet, nEntriesDetTotal[iDet]);
   }
 
   //----------------------------------------------------
@@ -105,16 +137,23 @@ void manualCalibFit()
   printf("-------- Started fits\n");
   std::array<float, 540> laFitResults{};
   std::array<float, 540> vdFitResults{};
+
+  TH1F* hVd = new TH1F("hVd", "v drift", 150, 0.5, 2.);
+  TH1F* hLa = new TH1F("hLa", "lorentz angle", 200, -25., 25.);
+  o2::trd::CalVdriftExB* calObjectOut = new o2::trd::CalVdriftExB();
+
   for (int iDet = 0; iDet < 540; ++iDet) {
+    if (nEntriesDetTotal[iDet] < 75)
+      continue;
     mFitFunctor.currDet = iDet;
     ROOT::Fit::Fitter fitter;
     double paramsStart[2];
-    paramsStart[0] = 0. * TMath::DegToRad();
+    paramsStart[0] = 0.;
     paramsStart[1] = 1.;
     fitter.SetFCN<o2::trd::FitFunctor>(2, mFitFunctor, paramsStart);
     fitter.Config().ParSettings(0).SetLimits(-0.7, 0.7);
     fitter.Config().ParSettings(0).SetStepSize(.01);
-    fitter.Config().ParSettings(1).SetLimits(0., 3.);
+    fitter.Config().ParSettings(1).SetLimits(0.01, 3.);
     fitter.Config().ParSettings(1).SetStepSize(.01);
     ROOT::Math::MinimizerOptions opt;
     opt.SetMinimizerType("Minuit2");
@@ -127,14 +166,28 @@ void manualCalibFit()
     auto fitResult = fitter.Result();
     laFitResults[iDet] = fitResult.Parameter(0);
     vdFitResults[iDet] = fitResult.Parameter(1);
-    printf("Det %d: la=%f\tvd=%f\n", iDet, laFitResults[iDet] * TMath::RadToDeg(), vdFitResults[iDet]);
+    if (fitResult.MinFcnValue() > 0.03)
+      continue;
+    printf("Det %d: la=%.3f \tvd=%.3f \t100*minValue=%f \tentries=%d\n", iDet, laFitResults[iDet] * TMath::RadToDeg(), vdFitResults[iDet], 100 * fitResult.MinFcnValue(), nEntriesDetTotal[iDet]);
+    hVd->Fill(vdFitResults[iDet]);
+    hLa->Fill(laFitResults[iDet] * TMath::RadToDeg());
+    calObjectOut->setVdrift(iDet, vdFitResults[iDet]);
+    calObjectOut->setExB(iDet, laFitResults[iDet]);
   }
   printf("-------- Finished fits\n");
 
+  std::cout << "number of chambers with enough entries: " << hVd->GetEntries() << std::endl;
+  ;
+  std::cout << "vdrift mean: " << hVd->GetMean() << " sigma: " << hVd->GetStdDev() << std::endl;
+  std::cout << "lorentz angle mean: " << hLa->GetMean() << " sigma: " << hLa->GetStdDev() << std::endl;
+
   //----------------------------------------------------
   // Write
   //----------------------------------------------------
   std::unique_ptr<TFile> outFilePtr(TFile::Open("manualCalibFit.root", "RECREATE"));
+  hVd->Write();
+  hLa->Write();
+  outFilePtr->WriteObjectAny(calObjectOut, "o2::trd::CalVdriftExB", "calObject");
   for (auto& p : mFitFunctor.profiles)
     p->Write();
 }