TRD: improve anchoring using CCDB objects (#15890)

Author: glegras

Detectors/TRD/simulation/include/TRDSimulation/Digitizer.h

@@ -14,14 +14,17 @@
 
 #include "TRDSimulation/Detector.h"
 
-#include "TRDBase/Calibrations.h"
+#include "TRDBase/Geometry.h"
 #include "TRDBase/PadResponse.h"
 #include "TRDSimulation/DiffAndTimeStructEstimator.h"
 #include "TRDSimulation/PileupTool.h"
 
 #include "DataFormatsTRD/Digit.h"
 #include "DataFormatsTRD/SignalArray.h"
 #include "DataFormatsTRD/Constants.h"
+#include "DataFormatsTRD/CalVdriftExB.h"
+#include "DataFormatsTRD/CalGain.h"
+#include "TRDBase/PadCalibrationsAliases.h"
 
 #include "MathUtils/RandomRing.h"
 #include "SimulationDataFormat/MCTruthContainer.h"
@@ -61,19 +64,40 @@ class Digitizer
   void setTriggerTime(double t) { mCurrentTriggerTime = t; }
   void setEventID(int entryID) { mEventID = entryID; }
   void setSrcID(int sourceID) { mSrcID = sourceID; }
-  void setCalibrations(Calibrations* calibrations) { mCalib = calibrations; }
+  void setCalVdriftExB(const CalVdriftExB* calVdriftExB) { mCalVdriftExB = calVdriftExB; }
+  void setFedChamberStatus(int det, int value)
+  {
+    if (det >= constants::MAXCHAMBER) {
+      return;
+    }
+    mFedChamberStatus[det] = value;
+  }
+  void setCalGain(const CalGain* calGain) { mCalGain = calGain; }
+  void setLocalGainFactor(const LocalGainFactor* localGain) { mLocalGainFactor = localGain; }
+  void setPadStatus(const PadStatus* padStatus) { mPadStatus = padStatus; }
   void setCreateSharedDigits(bool flag) { mCreateSharedDigits = flag; }
   int getEventTime() const { return mTime; }
   int getEventID() const { return mEventID; }
   int getSrcID() const { return mSrcID; }
   bool getCreateSharedDigits() const { return mCreateSharedDigits; }
+  int getFedChamberStatus(int det)
+  {
+    if (det >= constants::MAXCHAMBER) {
+      return -1.;
+    }
+    return mFedChamberStatus[det];
+  }
   std::string dumpFlaggedChambers() const;
 
  private:
-  Geometry* mGeo = nullptr;               // access to Geometry
-  PadResponse mPRF{};                     // access to PadResponse
-  SimParam mSimParam{};                   // simulation parameters
-  Calibrations* mCalib = nullptr;         // access to Calibrations in CCDB
+  Geometry* mGeo = nullptr;                                   // access to Geometry
+  PadResponse mPRF{};                                         // access to PadResponse
+  SimParam mSimParam{};                                       // simulation parameters
+  const CalVdriftExB* mCalVdriftExB = nullptr;                // access to time-dependent calibrations CalVdriftExB in CCDB
+  const CalGain* mCalGain = nullptr;                          // access to time-dependent calibrations CalGain in CCDB
+  const PadStatus* mPadStatus = nullptr;                      // access to time-dependent calibrations PadStatus in CCDB
+  const LocalGainFactor* mLocalGainFactor = nullptr;          // access to time-dependent calibrations LocalGainFactor in CCDB
+  std::array<int, constants::MAXCHAMBER> mFedChamberStatus{}; // access to time-dependent FED chamber status in CCDB
   PileupTool pileupTool;
 
   // number of digitizer threads

Detectors/TRD/simulation/include/TRDSimulation/SimParam.h

@@ -66,6 +66,8 @@ class SimParam
   }
 
   // Getters
+  float getT0() const { return mT0; }
+  float getEffVdriftDefault() const { return mEffVdriftDefault; }
   float getGasGain() const { return mGasGain; }
   float getNoise() const { return mNoise; }
   float getChipGain() const { return mChipGain; }
@@ -135,6 +137,9 @@ class SimParam
   float mTRFhi{3.58f};                                               ///<  Higher boundary of the TRF and x-talk
   float mInvTRFwid{static_cast<float>(mTRFbin) / (mTRFhi - mTRFlo)}; ///<  Inverse of the bin width of the integrated TRF and x-talk
   float mGasGain{4000.f};                                            ///< Gas gain
+  float mT0{-0.130f};                                                ///< time delay before amplification peak
+
+  float mEffVdriftDefault{1.22f}; ///< The CCDB contains an effective drift velocity, and not the true drift velocity, so we anchor only the variations. This default parameter is the one reproducing data the best for run 568191.
 
   ClassDefNV(SimParam, 2); // The TRD simulation parameters
 };

Detectors/TRD/simulation/include/TRDSimulation/TRDSimParams.h

@@ -41,6 +41,7 @@ struct TRDSimParams : public o2::conf::ConfigurableParamHelper<TRDSimParams> {
   int trf = 0;                  // Sampled TRF function. 0: default TRF, 1: TRF described in TRF TDR, 2: No TRF, 3: Landau dist as TRF (parameters specified below)
   float trf_landau_mu = 0.;     // Mu of the Landau distribution used to describe TRF
   float trf_landau_sigma = .03; // Sigma of the Landau distribution used to describe TRF
+  float t0 = -0.130;            // Time delay before start of amplification peak, in us
   O2ParamDef(TRDSimParams, "TRDSimParams");
 };
 

Detectors/TRD/simulation/src/Digitizer.cxx

@@ -139,9 +139,6 @@ void Digitizer::clearContainers()
 
 void Digitizer::process(std::vector<Hit> const& hits)
 {
-  if (!mCalib) {
-    LOG(fatal) << "TRD Calibration database not available";
-  }
 
   // Get the a hit container for all the hits in a given detector then call convertHits for a given detector (0 - 539)
   std::array<std::vector<Hit>, MAXCHAMBER> hitsPerDetector;
@@ -160,7 +157,8 @@ void Digitizer::process(std::vector<Hit> const& hits)
     auto& signalsMap = mSignalsMapCollection[det];
     // Jump to the next detector if the detector is
     // switched off, not installed, etc
-    if (mCalib->getChamberStatus()->isNoData(det)) {
+    // In the chamber Fed status, 3 corresponds to good chamber
+    if (getFedChamberStatus(det) != 3) {
       continue;
     }
     if (!mGeo->chamberInGeometry(det)) {
@@ -199,7 +197,7 @@ bool Digitizer::convertHits(const int det, const std::vector<Hit>& hits, SignalC
 
   double padSignal[mNpad];
 
-  const double calExBDetValue = mCalib->getExB(det); // T * V/cm (check units)
+  const double calExBDetValue = mCalVdriftExB->getExB(det, true); // T * V/cm (check units)
   const PadPlane* padPlane = mGeo->getPadPlane(det);
   const int layer = mGeo->getLayer(det);
   const float rowEndROC = padPlane->getRowEndROC();
@@ -253,8 +251,11 @@ bool Digitizer::convertHits(const int det, const std::vector<Hit>& hits, SignalC
       absDriftLength /= std::sqrt(1 / (1 + calExBDetValue * calExBDetValue));
     }
 
-    float driftVelocity = mCalib->getVDrift(det, colE, rowE); // The drift velocity
-    float t0 = mCalib->getT0(det, colE, rowE);                // The T0 velocity
+    float driftVelocity = constants::VDRIFTDEFAULT; // The default drift velocity
+    if (TMath::Abs(mCalVdriftExB->getVdrift(det, true) - constants::VDRIFTDEFAULT) > 1e-6) {
+      driftVelocity = mCalVdriftExB->getVdrift(det, true) * constants::VDRIFTDEFAULT / mSimParam.getEffVdriftDefault(); // If they are available in the CCDB, we anchor the vdrift variations
+    }
+    float t0 = mSimParam.getT0(); // The T0 velocity
 
     // Loop over all created electrons
     const int nElectrons = std::fabs(qTotal);
@@ -277,7 +278,8 @@ bool Digitizer::convertHits(const int det, const std::vector<Hit>& hits, SignalC
 
       // Apply E x B effects
       if (mSimParam.isExBOn()) {
-        locCd = locCd + calExBDetValue * driftLength;
+        // minus sign is necessary to be compatible with how Vdrift and ExB are calibrated
+        locCd = locCd - TMath::Tan(calExBDetValue) * driftLength;
       }
       // The electron position after diffusion and ExB in pad coordinates.
       rowE = padPlane->getPadRowNumberROC(locRd);
@@ -292,8 +294,6 @@ bool Digitizer::convertHits(const int det, const std::vector<Hit>& hits, SignalC
         continue;
       }
       const double colOffset = padPlane->getPadColOffset(colE, locCd + offsetTilt);
-      driftVelocity = mCalib->getVDrift(det, colE, rowE); // The drift velocity for the updated col and row
-      t0 = mCalib->getT0(det, colE, rowE);                // The T0 velocity for the updated col and row
       // Convert the position to drift time [mus], using either constant drift velocity or
       // time structure of drift cells (non-isochronity, GARFIELD calculation).
       // Also add absolute time of hits to take pile-up events into account properly
@@ -313,7 +313,8 @@ bool Digitizer::convertHits(const int det, const std::vector<Hit>& hits, SignalC
       }
 
       // Apply the gas gain including fluctuations
-      const double signal = -(mSimParam.getGasGain()) * mLogRandomRings[thread].getNextValue();
+      double timeDepGainFactor = mCalGain->getMPVdEdx(det, true) / constants::MPVDEDXDEFAULT;
+      const double signal = -(mSimParam.getGasGain()) * timeDepGainFactor * mLogRandomRings[thread].getNextValue();
 
       // Apply the pad response
       if (mSimParam.prfOn()) {
@@ -332,7 +333,7 @@ bool Digitizer::convertHits(const int det, const std::vector<Hit>& hits, SignalC
       }
 
       // The time bin (always positive), with t0 distortion
-      double timeBinIdeal = driftTime * mSamplingRate + t0;
+      double timeBinIdeal = (driftTime + t0) * mSamplingRate;
       // Protection
       if (std::fabs(timeBinIdeal) > (2 * mMaxTimeBins)) { // OS: why 2*mMaxTimeBins?
         timeBinIdeal = 2 * mMaxTimeBins;
@@ -423,19 +424,13 @@ bool Digitizer::convertSignalsToADC(SignalContainer& signalMapCont, DigitContain
     int mcm = (int)(col / 18);               // current group of 18 col pads
     int halfchamberside = (mcm > 3) ? 1 : 0; // 0=Aside, 1=Bside
 
-    // Halfchambers that are switched off, masked by mCalib
-    if ((halfchamberside == 0 && mCalib->getChamberStatus()->isNoDataSideA(det)) ||
-        (halfchamberside == 1 && mCalib->getChamberStatus()->isNoDataSideB(det))) {
-      continue;
-    }
-
     // Check whether pad is masked
     // Bridged pads are not considered yet!!!
-    if (mCalib->getPadStatus()->isMasked(det, col, row) || mCalib->getPadStatus()->isNotConnected(det, col, row)) {
+    if (mPadStatus->isMasked(det, col, row) || mPadStatus->isNotConnected(det, col, row)) {
       continue;
     }
 
-    float padgain = mCalib->getPadGainFactor(det, row, col); // The gain factor
+    float padgain = std::fabs(mLocalGainFactor->getValue(det, col, row)); // The gain factor
     if (padgain <= 0) {
       LOG(fatal) << "Not a valid gain " << padgain << ", " << det << ", " << col << ", " << row;
     }

Detectors/TRD/simulation/src/SimParam.cxx

@@ -29,6 +29,7 @@ SimParam::SimParam()
   mTRF = TRDSimParams::Instance().trf;
   mMu = TRDSimParams::Instance().trf_landau_mu;
   mSigma = TRDSimParams::Instance().trf_landau_sigma;
+  mT0 = TRDSimParams::Instance().t0;
   init();
 }
 

Detectors/TRD/workflow/src/TRDDigitizerSpec.cxx

@@ -15,6 +15,7 @@
 #include "Framework/DataProcessorSpec.h"
 #include "Framework/DataRefUtils.h"
 #include "Framework/Lifetime.h"
+#include "Framework/CCDBParamSpec.h"
 #include "Headers/DataHeader.h"
 #include "TStopwatch.h"
 #include "TChain.h"
@@ -64,11 +65,11 @@ class TRDDPLDigitizerTask : public o2::base::BaseDPLDigitizer
 
     bool mctruth = pc.outputs().isAllowed({"TRD", "LABELS", 0});
 
-    Calibrations simcal;
-    // the timestamp can be extracted from the DPL header (it is set in SimReader)
-    auto creationTime = pc.services().get<o2::framework::TimingInfo>().creation;
-    simcal.getCCDBObjects(creationTime);
-    mDigitizer.setCalibrations(&simcal);
+    pc.inputs().get<o2::trd::CalVdriftExB*>("calvdexb");
+    pc.inputs().get<o2::trd::CalGain*>("calgain");
+    pc.inputs().get<std::array<int, constants::MAXCHAMBER>*>("fedchamberstatus");
+    pc.inputs().get<o2::trd::PadStatus*>("padstatus");
+    pc.inputs().get<o2::trd::LocalGainFactor*>("localgainfactor");
 
     // read collision context from input
     auto context = pc.inputs().get<o2::steer::DigitizationContext*>("collisioncontext");
@@ -207,6 +208,38 @@ class TRDDPLDigitizerTask : public o2::base::BaseDPLDigitizer
     finished = true;
   }
 
+  void finaliseCCDB(ConcreteDataMatcher& matcher, void* obj)
+  {
+    if (matcher == ConcreteDataMatcher("TRD", "CALVDRIFTEXB", 0)) {
+      LOG(info) << "CalVdriftExB object has been updated";
+      mDigitizer.setCalVdriftExB((const o2::trd::CalVdriftExB*)obj);
+      return;
+    }
+    if (matcher == ConcreteDataMatcher("TRD", "CALGAIN", 0)) {
+      LOG(info) << "CalGain object has been updated";
+      mDigitizer.setCalGain((const o2::trd::CalGain*)obj);
+      return;
+    }
+    if (matcher == ConcreteDataMatcher("TRD", "FEDCHAMBERSTATUS", 0)) {
+      LOG(info) << "Fed chamber status has been updated";
+      std::array<int, constants::MAXCHAMBER>* fedStatus = (std::array<int, constants::MAXCHAMBER>*)obj;
+      for (int det = 0; det < constants::MAXCHAMBER; det++) {
+        mDigitizer.setFedChamberStatus(det, (*fedStatus)[det]);
+      }
+      return;
+    }
+    if (matcher == ConcreteDataMatcher("TRD", "PADSTATUS", 0)) {
+      LOG(info) << "PadStatus object has been updated";
+      mDigitizer.setPadStatus((const o2::trd::PadStatus*)obj);
+      return;
+    }
+    if (matcher == ConcreteDataMatcher("TRD", "LOCALGAINFACTOR", 0)) {
+      LOG(info) << "LocalGainFactor object has been updated";
+      mDigitizer.setLocalGainFactor((const o2::trd::LocalGainFactor*)obj);
+      return;
+    }
+  }
+
  private:
   Digitizer mDigitizer;
   std::vector<TChain*> mSimChains;
@@ -230,10 +263,18 @@ o2::framework::DataProcessorSpec getTRDDigitizerSpec(int channel, bool mctruth)
   }
   outputs.emplace_back("TRD", "ROMode", 0, Lifetime::Timeframe);
 
+  std::vector<InputSpec> inputs;
+  inputs.emplace_back("collisioncontext", "SIM", "COLLISIONCONTEXT", static_cast<SubSpecificationType>(channel), Lifetime::Timeframe);
+  inputs.emplace_back("calvdexb", "TRD", "CALVDRIFTEXB", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/CalVdriftExB"));
+  inputs.emplace_back("calgain", "TRD", "CALGAIN", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/CalGain"));
+  inputs.emplace_back("fedchamberstatus", "TRD", "FEDCHAMBERSTATUS", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/DCSDPsFedChamberStatus"));
+  inputs.emplace_back("padstatus", "TRD", "PADSTATUS", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/PadStatus"));
+  inputs.emplace_back("localgainfactor", "TRD", "LOCALGAINFACTOR", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/LocalGainFactor"));
+
   return DataProcessorSpec{
     "TRDDigitizer",
-    Inputs{InputSpec{"collisioncontext", "SIM", "COLLISIONCONTEXT", static_cast<SubSpecificationType>(channel), Lifetime::Timeframe}},
-
+    // Inputs{InputSpec{"collisioncontext", "SIM", "COLLISIONCONTEXT", static_cast<SubSpecificationType>(channel), Lifetime::Timeframe}},
+    inputs,
     outputs,
 
     AlgorithmSpec{adaptFromTask<TRDDPLDigitizerTask>()},