ZDC - Fixes to intercalibration, waveform extraction. New workflow to…

… parse and analyze CTF data (#13623)

* Fixes + orthogonal regression

* Various fixes

* Small updates

* Digits parser workflow

* Digits parser workflow

* Digits parser workflow

* WIP

* WIP

* Debugging line shapes and hits directly from CTF

* Fixes

* Fix bug in hit checking

* Please consider the following formatting changes (#89)

* Fix missing case

* Please consider the following formatting changes (#90)

* Small updates

* Comments

* Please consider the following formatting changes (#91)

* Fix compilation error

* Just to force new tests

---------

Co-authored-by: ALICE Builder <alibuild@users.noreply.github.com>

DataFormats/Detectors/ZDC/include/DataFormatsZDC/BCData.h

@@ -55,8 +55,9 @@ struct BCData {
   o2::dataformats::RangeRefComp<6> ref;
   o2::InteractionRecord ir;
   std::array<uint16_t, NModules> moduleTriggers{};
+  // N.B. channels and triggers have geographical addressing (0x1 << (NChPerModule * im + ic)
   uint32_t channels = 0;    // pattern of channels it refers to
-  uint32_t triggers = 0;    // pattern of triggered channels (not necessarily stored) in this BC
+  uint32_t triggers = 0;    // pattern of triggered channels (not necessarily stored) in this BC (i.e. with Hit bit on)
   uint8_t ext_triggers = 0; // pattern of ALICE triggers
 
   BCData() = default;

Detectors/ZDC/calib/include/ZDCCalib/CalibParamZDC.h

@@ -27,6 +27,7 @@ namespace o2
 namespace zdc
 {
 struct CalibParamZDC : public o2::conf::ConfigurableParamHelper<CalibParamZDC> {
+  bool dumpCalib = false;                            // Dump partial calibration object
   bool debugOutput = false;                          // Debug output
   bool rootOutput = true;                            // Output histograms to EOS
   std::string outputDir = "./";                      // ROOT files output directory

Detectors/ZDC/calib/include/ZDCCalib/InterCalib.h

@@ -80,6 +80,8 @@ class InterCalib
   void setInterCalibConfig(const InterCalibConfig* param) { mInterCalibConfig = param; };
   const InterCalibConfig* getInterCalibConfig() const { return mInterCalibConfig; };
 
+  InterCalibData& getData() { return mData; };
+
   void setVerbosity(int v) { mVerbosity = v; }
   int getVerbosity() const { return mVerbosity; }
 

Detectors/ZDC/calib/include/ZDCCalib/InterCalibConfig.h

@@ -36,9 +36,13 @@ struct InterCalibConfig {
   // Meaningful values are in the range of tower x centers i.e. from
   // 2.8 to 19.6 If one puts less than 2.8 then the computation will be
   // the same as for ZPA/ZPC with no cuts
-  double xcut_ZPA = 6;
-  double xcut_ZPC = 6;
-  double tower_cut_ZP = 0;
+  double xcut_ZPA = 0;
+  double xcut_ZPC = 0;
+  double rms_cut_ZP = 0;                                                                                                                                                                         // RMS of ZP centroid can go from 0 to 8.4 cm
+  double towerCutLow_ZPA[4] = {0, 0, 0, 0};                                                                                                                                                      // Applied to all ZP fits except ZPI
+  double towerCutLow_ZPC[4] = {0, 0, 0, 0};                                                                                                                                                      // Applied to all ZP fits except ZPI
+  double towerCutHigh_ZPA[4] = {std::numeric_limits<float>::infinity(), std::numeric_limits<float>::infinity(), std::numeric_limits<float>::infinity(), std::numeric_limits<float>::infinity()}; // Applied to all ZP fits except ZPI
+  double towerCutHigh_ZPC[4] = {std::numeric_limits<float>::infinity(), std::numeric_limits<float>::infinity(), std::numeric_limits<float>::infinity(), std::numeric_limits<float>::infinity()}; // Applied to all ZP fits except ZPI
   bool cross_check = false;
 
   int nb1[NH] = {0};      /// 1D histogram: number of bins
@@ -87,7 +91,7 @@ struct InterCalibConfig {
     enabled[7] = c7;
     enabled[8] = c8;
   }
-  ClassDefNV(InterCalibConfig, 4);
+  ClassDefNV(InterCalibConfig, 5);
 };
 } // namespace zdc
 } // namespace o2

Detectors/ZDC/calib/include/ZDCCalib/WaveformCalibData.h

@@ -87,6 +87,7 @@ struct WaveformCalibData {
   void setCreationTime(uint64_t ctime);
   void setN(int n);
   int saveDebugHistos(const std::string fn);
+  int dumpCalib(const std::string fn);
   ClassDefNV(WaveformCalibData, 1);
 };
 

Detectors/ZDC/calib/include/ZDCCalib/WaveformCalibEPN.h

@@ -37,9 +37,11 @@ class WaveformCalibEPN
               const gsl::span<const o2::zdc::ZDCWaveform>& wave);
   int endOfRun();
   int saveDebugHistos(const std::string fn = "ZDCWaveformCalibEPN.root");
+  int dumpCalib(const std::string fn = "ZDCWaveformCalibEPNDump.root");
   void setConfig(const WaveformCalibConfig* param) { mConfig = param; };
   const WaveformCalibConfig* getConfig() const { return mConfig; };
   void setSaveDebugHistos() { mSaveDebugHistos = true; }
+  void setDumpCalib() { mDumpCalib = true; }
   void setDontSaveDebugHistos() { mSaveDebugHistos = false; }
   void setVerbosity(int val) { mVerbosity = val; }
   WaveformCalibData mData;
@@ -48,6 +50,7 @@ class WaveformCalibEPN
  private:
   bool mInitDone = false;
   bool mSaveDebugHistos = false;
+  bool mDumpCalib = false;
   int32_t mNBin = 0;
   int32_t mVerbosity = DbgMinimal;
   const WaveformCalibConfig* mConfig = nullptr; /// Configuration of intercalibration

Detectors/ZDC/calib/src/InterCalib.cxx

@@ -210,12 +210,12 @@ void InterCalib::assign(int ih, bool ismod)
   } else if (ih == 5) {
     nid = 1;
     id = id_5;
-    LOG(warn) << "InterCalib::assign unimplemented coefficient ih = " << ih;
+    LOG(warn) << "InterCalib::assign is not implemented for coefficient ih = " << ih;
     return;
   } else if (ih == 6) {
     nid = 1;
     id = id_6;
-    LOG(warn) << "InterCalib::assign unimplemented coefficient ih = " << ih;
+    LOG(warn) << "InterCalib::assign is not implemented for coefficient ih = " << ih;
     return;
   } else if (ih == 7 || ih == 8) {
     nid = 4;
@@ -246,15 +246,23 @@ void InterCalib::assign(int ih, bool ismod)
       if (oldval > 0) {
         val = val * mPar[ih][iid + 1];
       }
-      if (mVerbosity > DbgZero) {
+      if (mTowerParamUpd.modified[ich]) {
+        LOGF(warn, "%s OVERWRITING MODIFIED PARAMETER %8.6f", ChannelNames[ich].data(), mTowerParamUpd.getTowerCalib(ich));
+        LOGF(info, "%s updated %8.6f -> %8.6f", ChannelNames[ich].data(), oldval, val);
+      } else if (mVerbosity > DbgZero) {
         LOGF(info, "%s updated %8.6f -> %8.6f", ChannelNames[ich].data(), oldval, val);
       }
       mTowerParamUpd.setTowerCalib(ich, val, true);
     } else {
-      if (mVerbosity > DbgZero) {
-        LOGF(info, "%s NOT CHANGED %8.6f", ChannelNames[ich].data(), oldval);
+      // Check if another fit has already modified the parameters
+      if (mTowerParamUpd.modified[ich]) {
+        LOGF(warn, "%s NOT OVERWRITING MODIFIED PARAMETER %8.6f", ChannelNames[ich].data(), mTowerParamUpd.getTowerCalib(ich));
+      } else {
+        if (mVerbosity > DbgZero) {
+          LOGF(info, "%s NOT CHANGED %8.6f", ChannelNames[ich].data(), oldval);
+        }
+        mTowerParamUpd.setTowerCalib(ich, oldval, false);
       }
-      mTowerParamUpd.setTowerCalib(ich, oldval, false);
     }
   }
 }
@@ -294,6 +302,10 @@ int InterCalib::process(const char* hname, int ic)
     ih = HidZNI;
   } else if (hn.EqualTo("hZPI")) {
     ih = HidZPI;
+  } else if (hn.EqualTo("hZPAX")) {
+    ih = HidZPAX;
+  } else if (hn.EqualTo("hZPCX")) {
+    ih = HidZPCX;
   } else {
     LOGF(error, "Not recognized histogram name: %s\n", hname);
     return -1;
@@ -434,18 +446,32 @@ void InterCalib::add(int ih, o2::dataformats::FlatHisto2D<float>& h2)
 
 void InterCalib::cumulate(int ih, double tc, double t1, double t2, double t3, double t4, double w = 1)
 {
+  constexpr double minfty = -std::numeric_limits<double>::infinity();
   if (tc < mInterCalibConfig->cutLow[ih] || tc > mInterCalibConfig->cutHigh[ih]) {
     return;
   }
-  double val[NPAR] = {0, 0, 0, 0, 0, 1};
-  val[0] = tc;
-  val[1] = t1;
-  val[2] = t2;
-  val[3] = t3;
-  val[4] = t4;
-  for (int32_t i = 0; i < NPAR; i++) {
-    for (int32_t j = i; j < NPAR; j++) {
-      mData.mSum[ih][i][j] += val[i] * val[j] * w;
+  if ((ih == HidZPA || ih == HidZPAX)) {
+    if (t1 < mInterCalibConfig->towerCutLow_ZPA[0] || t2 < mInterCalibConfig->towerCutLow_ZPA[1] || t3 < mInterCalibConfig->towerCutLow_ZPA[2] || t4 < mInterCalibConfig->towerCutLow_ZPA[3]) {
+      return;
+    }
+    if (t1 > mInterCalibConfig->towerCutHigh_ZPA[0] || t2 > mInterCalibConfig->towerCutHigh_ZPA[1] || t3 > mInterCalibConfig->towerCutHigh_ZPA[2] || t4 > mInterCalibConfig->towerCutHigh_ZPA[3]) {
+      return;
+    }
+  }
+  if (ih == HidZPC || ih == HidZPCX) {
+    if (t1 < mInterCalibConfig->towerCutLow_ZPC[0] || t2 < mInterCalibConfig->towerCutLow_ZPC[1] || t3 < mInterCalibConfig->towerCutLow_ZPC[2] || t4 < mInterCalibConfig->towerCutLow_ZPC[3]) {
+      return;
+    }
+    if (t1 > mInterCalibConfig->towerCutHigh_ZPC[0] || t2 > mInterCalibConfig->towerCutHigh_ZPC[1] || t3 > mInterCalibConfig->towerCutHigh_ZPC[2] || t4 > mInterCalibConfig->towerCutHigh_ZPC[3]) {
+      return;
+    }
+  }
+  double val[NPAR] = {tc, t1, t2, t3, t4, 1};
+  if (tc > minfty && t1 > minfty && t2 > minfty && t3 > minfty && t4 > minfty) {
+    for (int32_t i = 0; i < NPAR; i++) {
+      for (int32_t j = i; j < NPAR; j++) {
+        mData.mSum[ih][i][j] += val[i] * val[j] * w;
+      }
     }
   }
   // mData.mSum[ih][5][5] contains the number of analyzed events
@@ -470,6 +496,9 @@ void InterCalib::fcn(int& npar, double* gin, double& chi, double* par, int iflag
       chi += (i == 0 ? par[i] : -par[i]) * (j == 0 ? par[j] : -par[j]) * mAdd[i][j];
     }
   }
+  // Following line modifies the chisquare computation (sum of squares of residuals)
+  // to perform orthogonal least squares instead of ordinary least squares minimization
+  chi = chi / (1 + par[1] * par[1] + par[2] * par[2] + par[3] * par[3] + par[4] * par[4]);
 }
 
 int InterCalib::mini(int ih)
@@ -498,15 +527,11 @@ int InterCalib::mini(int ih)
   // Calibration cvoefficient is forced to and step is forced to zero
   mMn[ih]->mnparm(0, "c0", 1., 0., 1., 1., ierflg);
 
-  // Special fit for proton calorimeters: fit least exposed towers with using previous
-  // fit of all towers
+  // Special fit for proton calorimeters: fit least exposed towers
+  // starting from parameters of previous fit to all towers
 
   // Tower 1
-  if (ih == HidZPCX) {
-    mMn[ih]->mnparm(1, "c1", mPar[HidZPC][1], 0, l_bnd, u_bnd, ierflg);
-  } else {
-    mMn[ih]->mnparm(1, "c1", start, step, l_bnd, u_bnd, ierflg);
-  }
+  mMn[ih]->mnparm(1, "c1", start, step, l_bnd, u_bnd, ierflg);
 
   // Tower 2
   // Only two ZEM calorimeters: equalize response
@@ -518,20 +543,11 @@ int InterCalib::mini(int ih)
     step = 0;
   }
 
-  if (ih == HidZPCX) {
-    mMn[ih]->mnparm(2, "c2", mPar[HidZPC][2], 0, l_bnd, u_bnd, ierflg);
-  } else {
-    mMn[ih]->mnparm(2, "c2", start, step, l_bnd, u_bnd, ierflg);
-  }
+  mMn[ih]->mnparm(2, "c2", start, step, l_bnd, u_bnd, ierflg);
 
   // Towers 3 and 4
-  if (ih == HidZPAX) {
-    mMn[ih]->mnparm(3, "c3", mPar[HidZPA][3], 0, l_bnd, u_bnd, ierflg);
-    mMn[ih]->mnparm(4, "c4", mPar[HidZPA][4], 0, l_bnd, u_bnd, ierflg);
-  } else {
-    mMn[ih]->mnparm(3, "c3", start, step, l_bnd, u_bnd, ierflg);
-    mMn[ih]->mnparm(4, "c4", start, step, l_bnd, u_bnd, ierflg);
-  }
+  mMn[ih]->mnparm(3, "c3", start, step, l_bnd, u_bnd, ierflg);
+  mMn[ih]->mnparm(4, "c4", start, step, l_bnd, u_bnd, ierflg);
 
   // Offset
   l_bnd = mInterCalibConfig->l_bnd_o[ih];
@@ -551,13 +567,15 @@ int InterCalib::mini(int ih)
       l_bnd = mInterCalibConfig->l_bnd[ih];
       u_bnd = mInterCalibConfig->u_bnd[ih];
       for (int i = 1; i <= 4; i++) {
-        if (TMath::Abs(mPar[ih][i] - l_bnd) < 1e-3 || TMath::Abs(mPar[ih][i] - u_bnd) < 1e-3) {
+        if (TMath::Abs(mPar[ih][i] - l_bnd) < 1e-2 || TMath::Abs(mPar[ih][i] - u_bnd) < 1e-2) {
           retry = true;
           LOG(warn) << "ih=" << ih << " par " << i << " too close to boundaries";
           if (ih == 1 || ih == 7) {
-            mMn[ih]->mnparm(i, parn[i], mTowerParam->tower_calib[IdZPAC + i], 0, l_bnd, u_bnd, ierflg);
+            // mMn[ih]->mnparm(i, parn[i], mTowerParam->tower_calib[IdZPAC + i], 0, l_bnd, u_bnd, ierflg);
+            mMn[ih]->mnparm(i, parn[i], mInterCalibConfig->start[ih], 0, l_bnd, u_bnd, ierflg);
           } else if (ih == 3 || ih == 8) {
-            mMn[ih]->mnparm(i, parn[i], mTowerParam->tower_calib[IdZPCC + i], 0, l_bnd, u_bnd, ierflg);
+            // mMn[ih]->mnparm(i, parn[i], mTowerParam->tower_calib[IdZPCC + i], 0, l_bnd, u_bnd, ierflg);
+            mMn[ih]->mnparm(i, parn[i], mInterCalibConfig->start[ih], 0, l_bnd, u_bnd, ierflg);
           } else {
             LOG(fatal) << "ERROR on InterCalib minimization ih=" << ih;
           }

Detectors/ZDC/calib/src/InterCalibConfig.cxx

@@ -28,7 +28,11 @@ void InterCalibConfig::print() const
   }
   LOG(info) << "xcut_ZPA = " << xcut_ZPA;
   LOG(info) << "xcut_ZPC = " << xcut_ZPC;
-  LOG(info) << "tower_cut_ZP = " << tower_cut_ZP;
+  LOG(info) << "towerCutLow_ZPA = {" << towerCutLow_ZPA[0] << ", " << towerCutLow_ZPA[1] << ", " << towerCutLow_ZPA[2] << ", " << towerCutLow_ZPA[3] << "};";
+  LOG(info) << "towerCutHigh_ZPA = {" << towerCutHigh_ZPA[0] << ", " << towerCutHigh_ZPA[1] << ", " << towerCutHigh_ZPA[2] << ", " << towerCutHigh_ZPA[3] << "};";
+  LOG(info) << "towerCutLow_ZPC = {" << towerCutLow_ZPC[0] << ", " << towerCutLow_ZPC[1] << ", " << towerCutLow_ZPC[2] << ", " << towerCutLow_ZPC[3] << "};";
+  LOG(info) << "towerCutHigh_ZPC = {" << towerCutHigh_ZPC[0] << ", " << towerCutHigh_ZPC[1] << ", " << towerCutHigh_ZPC[2] << ", " << towerCutHigh_ZPC[3] << "};";
+  LOG(info) << "rms_cut_ZP = " << rms_cut_ZP;
   if (cross_check) {
     LOG(warn) << "THIS IS A CROSS CHECK CONFIGURATION (vs SUM)";
   }

Detectors/ZDC/calib/src/InterCalibEPN.cxx

@@ -106,7 +106,7 @@ int InterCalibEPN::process(const gsl::span<const o2::zdc::BCRecData>& RecBC,
         float x, rms;
         ev.centroidZPA(x, rms);
         cumulate(HidZPA, ev.EZDC(IdZPAC), ev.EZDC(IdZPA1), ev.EZDC(IdZPA2), ev.EZDC(IdZPA3), ev.EZDC(IdZPA4), 1.);
-        if (x < -(mInterCalibConfig->xcut_ZPA)) {
+        if (x < -(mInterCalibConfig->xcut_ZPA) && rms >= mInterCalibConfig->rms_cut_ZP) {
           cumulate(HidZPAX, ev.EZDC(IdZPAC), ev.EZDC(IdZPA1), ev.EZDC(IdZPA2), ev.EZDC(IdZPA3), ev.EZDC(IdZPA4), 1.);
         }
       }
@@ -117,7 +117,7 @@ int InterCalibEPN::process(const gsl::span<const o2::zdc::BCRecData>& RecBC,
         float x, rms;
         ev.centroidZPC(x, rms);
         cumulate(HidZPC, ev.EZDC(IdZPCC), ev.EZDC(IdZPC1), ev.EZDC(IdZPC2), ev.EZDC(IdZPC3), ev.EZDC(IdZPC4), 1.);
-        if (x > (mInterCalibConfig->xcut_ZPC)) {
+        if (x > (mInterCalibConfig->xcut_ZPC) && rms >= mInterCalibConfig->rms_cut_ZP) {
           cumulate(HidZPCX, ev.EZDC(IdZPCC), ev.EZDC(IdZPC1), ev.EZDC(IdZPC2), ev.EZDC(IdZPC3), ev.EZDC(IdZPC4), 1.);
         }
       }
@@ -266,22 +266,33 @@ void InterCalibEPN::clear(int ih)
 
 void InterCalibEPN::cumulate(int ih, double tc, double t1, double t2, double t3, double t4, double w = 1)
 {
+  constexpr double minfty = -std::numeric_limits<double>::infinity();
   // printf("%s: ih=%d tc=%g t1=%g t2=%g t3=%g t4=%g w=%g\n",__func__,ih, tc, t1, t2, t3, t4, w); fflush(stdout);
   if (tc < mInterCalibConfig->cutLow[ih] || tc > mInterCalibConfig->cutHigh[ih]) {
     return;
   }
-  if ((ih == 7 || ih == 8) && (t1 < mInterCalibConfig->tower_cut_ZP || t2 < mInterCalibConfig->tower_cut_ZP || t3 < mInterCalibConfig->tower_cut_ZP || t4 < mInterCalibConfig->tower_cut_ZP)) {
-    return;
+  if ((ih == HidZPA || ih == HidZPAX)) {
+    if (t1 < mInterCalibConfig->towerCutLow_ZPA[0] || t2 < mInterCalibConfig->towerCutLow_ZPA[1] || t3 < mInterCalibConfig->towerCutLow_ZPA[2] || t4 < mInterCalibConfig->towerCutLow_ZPA[3]) {
+      return;
+    }
+    if (t1 > mInterCalibConfig->towerCutHigh_ZPA[0] || t2 > mInterCalibConfig->towerCutHigh_ZPA[1] || t3 > mInterCalibConfig->towerCutHigh_ZPA[2] || t4 > mInterCalibConfig->towerCutHigh_ZPA[3]) {
+      return;
+    }
   }
-  double val[NPAR] = {0, 0, 0, 0, 0, 1};
-  val[0] = tc;
-  val[1] = t1;
-  val[2] = t2;
-  val[3] = t3;
-  val[4] = t4;
-  for (int32_t i = 0; i < NPAR; i++) {
-    for (int32_t j = i; j < NPAR; j++) {
-      mData.mSum[ih][i][j] += val[i] * val[j] * w;
+  if (ih == HidZPC || ih == HidZPCX) {
+    if (t1 < mInterCalibConfig->towerCutLow_ZPC[0] || t2 < mInterCalibConfig->towerCutLow_ZPC[1] || t3 < mInterCalibConfig->towerCutLow_ZPC[2] || t4 < mInterCalibConfig->towerCutLow_ZPC[3]) {
+      return;
+    }
+    if (t1 > mInterCalibConfig->towerCutHigh_ZPC[0] || t2 > mInterCalibConfig->towerCutHigh_ZPC[1] || t3 > mInterCalibConfig->towerCutHigh_ZPC[2] || t4 > mInterCalibConfig->towerCutHigh_ZPC[3]) {
+      return;
+    }
+  }
+  double val[NPAR] = {tc, t1, t2, t3, t4, 1};
+  if (tc > minfty && t1 > minfty && t2 > minfty && t3 > minfty && t4 > minfty) {
+    for (int32_t i = 0; i < NPAR; i++) {
+      for (int32_t j = i; j < NPAR; j++) {
+        mData.mSum[ih][i][j] += val[i] * val[j] * w;
+      }
     }
   }
   // mData.mSum[ih][5][5] contains the number of analyzed events

Detectors/ZDC/calib/src/WaveformCalibConfig.cxx

@@ -20,9 +20,10 @@ WaveformCalibConfig::WaveformCalibConfig()
     cutLow[isig] = -std::numeric_limits<float>::infinity();
     cutHigh[isig] = std::numeric_limits<float>::infinity();
   }
+  // Firmware aligns signals within one sample
   for (int itdc = 0; itdc < NTDCChannels; itdc++) {
-    cutTimeLow[itdc] = -1.25;
-    cutTimeHigh[itdc] = 1.25;
+    cutTimeHigh[itdc] = o2::constants::lhc::LHCBunchSpacingNS / NTimeBinsPerBC;
+    cutTimeLow[itdc] = -cutTimeHigh[itdc];
   }
 }
 

Detectors/ZDC/calib/src/WaveformCalibData.cxx

@@ -187,6 +187,21 @@ int WaveformCalibData::saveDebugHistos(const std::string fn)
   return 0;
 }
 
+//______________________________________________________________________________
+int WaveformCalibData::dumpCalib(const std::string fn)
+{
+  TDirectory* cwd = gDirectory;
+  TFile* f = new TFile(fn.data(), "recreate");
+  if (f->IsZombie()) {
+    LOG(error) << "Cannot create file: " << fn;
+    return 1;
+  }
+  f->WriteObjectAny((void*)this, o2::zdc::WaveformCalibData::Class(), "WaveformCalibData");
+  f->Close();
+  cwd->cd();
+  return 0;
+}
+
 //______________________________________________________________________________
 void WaveformCalibData::clear()
 {