[EMCAL-566] Added function for fast projection of boost histograms

jokonig · web-flow · commit b1827f026fb3 · 2022-01-27T21:31:27.000+04:00
diff --git a/Common/Utils/include/CommonUtils/BoostHistogramUtils.h b/Common/Utils/include/CommonUtils/BoostHistogramUtils.h
@@ -427,6 +427,31 @@ auto ProjectBoostHistoX(boost::histogram::histogram<axes...>& hist2d, const int
   return histoProj;
 }
 
+/// \brief Function to project 2d boost histogram onto x-axis
+/// \param hist2d 2d boost histogram
+/// \param binLow lower bin in y for projection
+/// \param binHigh lower bin in y for projection
+/// \return result
+///      1d boost histogram from projection of the input 2d boost histogram
+template <typename... axes>
+auto ProjectBoostHistoXFast(boost::histogram::histogram<axes...>& hist2d, const int binLow, const int binHigh)
+{
+  unsigned int nbins = hist2d.axis(0).size();
+  double binStartX = hist2d.axis(0).bin(0).lower();
+  double binEndX = hist2d.axis(0).bin(nbins - 1).upper();
+  auto histoProj = boost::histogram::make_histogram(boost::histogram::axis::regular<>(nbins, binStartX, binEndX));
+
+  // Now rewrite the bin content of the 1d histogram to get the summed bin content in the specified range
+  for (int x = 0; x < nbins; ++x) {
+    histoProj.at(x) = 0;
+    for (int y = binLow; y < binHigh; ++y) {
+      histoProj.at(x) = histoProj.at(x) + hist2d.at(x, y);
+    }
+  }
+
+  return histoProj;
+}
+
 } // end namespace utils
 } // end namespace o2
 
diff --git a/Detectors/EMCAL/calibration/include/EMCALCalibration/EMCALCalibExtractor.h b/Detectors/EMCAL/calibration/include/EMCALCalibration/EMCALCalibExtractor.h
@@ -93,7 +93,10 @@ class EMCALCalibExtractor
     return mOutputBCM;
   }
 
-  /// \brief For now a dummy function to calibrate the time.
+  /// \brief Calibrate time for all cells
+  /// \param hist -- 2d boost histogram: cell-time vs. cell-ID
+  /// \param minTime -- min. time considered for fit
+  /// \param maxTime -- max. time considered for fit
   template <typename... axes>
   o2::emcal::TimeCalibrationParams calibrateTime(boost::histogram::histogram<axes...>& hist, double minTime = 0, double maxTime = 1000)
   {
@@ -102,6 +105,8 @@ class EMCALCalibExtractor
 
     o2::emcal::TimeCalibrationParams TCP;
 
+    double mean = 0;
+
 #if (defined(WITH_OPENMP) && !defined(__CLING__))
     if (mNThreads < 1) {
       mNThreads = std::min(omp_get_max_threads(), NCELLS);
@@ -114,18 +119,18 @@ class EMCALCalibExtractor
 #endif
 
     for (unsigned int i = 0; i < NCELLS; ++i) {
+
       // project boost histogram to 1d just for 1 cell
-      auto boostHist1d = o2::utils::ProjectBoostHistoX(histReduced, i, i + 1);
-      // maybe cut histo to max +- 25 ns
+      auto boostHist1d = o2::utils::ProjectBoostHistoXFast(histReduced, i, i + 1);
+      // ToDo: maybe cut histo to max +- 25 ns
 
-      // fit with gaussian to extract mean
       try {
         auto fitValues = o2::utils::fitBoostHistoWithGaus<double>(boostHist1d);
-        double mean = fitValues.at(1);
+        mean = fitValues.at(1);
         // add mean to time calib params
         TCP.addTimeCalibParam(i, mean, 0);
       } catch (o2::utils::FitGausError_t) {
-        TCP.addTimeCalibParam(i, 400, 0); // 400 ns shift default value
+        TCP.addTimeCalibParam(i, mean, 0); // take calib value of last cell; or 400 ns shift default value
       }
     }
     return TCP;
diff --git a/Detectors/EMCAL/calibration/run/runCalibOffline.cxx b/Detectors/EMCAL/calibration/run/runCalibOffline.cxx
@@ -44,15 +44,19 @@ int main(int argc, char** argv)
   std::string ccdbServerPath;
   bool doBadChannelCalib;
   std::string nameCalibInputHist; // hCellIdVsTimeAbove300 for time, hCellIdVsEnergy for bad channel
+  std::string namePathStoreLocal; // name for path + histogram to store the calibration locally in root TH1 format
 
   unsigned int nthreads; // number of threads used by openMP
 
   unsigned long rangestart; //30/10/2021, 01:02:32 for run 505566 -> 1635548552000
   unsigned long rangeend;   // 30/10/2021, 02:31:10 for run 505566 -> 1635553870000
 
+  double timeRangeLow;
+  double timeRangeHigh;
+
   try {
     bpo::options_description desc("Allowed options");
-    desc.add_options()("help", "Print this help message")("CalibInputPath", bpo::value<std::string>()->required(), "Set root input histogram")("ccdbServerPath", bpo::value<std::string>()->default_value(o2::base::NameConf::getCCDBServer()), "Set path to ccdb server")("mode", bpo::value<std::string>()->required(), "Set if time or bad channel calib")("nameInputHisto", bpo::value<std::string>()->default_value("hCellIdVsTimeAbove300"), "Set name of input histogram")("nthreads", bpo::value<unsigned int>()->default_value(1), "Set number of threads for OpenMP")("timestampStart", bpo::value<unsigned long>()->default_value(1635548552000), "Set timestamp from start of run")("timestampEnd", bpo::value<unsigned long>()->default_value(1635553870000), "Set timestamp from end of run");
+    desc.add_options()("help", "Print this help message")("CalibInputPath", bpo::value<std::string>()->required(), "Set root input histogram")("ccdbServerPath", bpo::value<std::string>()->default_value(o2::base::NameConf::getCCDBServer()), "Set path to ccdb server")("mode", bpo::value<std::string>()->required(), "Set if time or bad channel calib")("nameInputHisto", bpo::value<std::string>()->default_value("hCellIdVsTimeAbove300"), "Set name of input histogram")("nthreads", bpo::value<unsigned int>()->default_value(1), "Set number of threads for OpenMP")("timestampStart", bpo::value<unsigned long>()->default_value(1635548552000), "Set timestamp from start of run")("timestampEnd", bpo::value<unsigned long>()->default_value(1635553870000), "Set timestamp from end of run")("namePathStoreLocal", bpo::value<std::string>()->default_value(""), "Set path to store histo of time calib locally")("timeRangeLow", bpo::value<double>()->default_value(1), "Set lower boundary of fit interval for time calibration (in ns)")("timeRangeHigh", bpo::value<double>()->default_value(1000), "Set upper boundary of fit interval for time calibration (in ns)");
 
     bpo::store(bpo::parse_command_line(argc, argv, desc), vm);
 
@@ -116,6 +120,23 @@ int main(int argc, char** argv)
       rangeend = vm["timestampEnd"].as<unsigned long>();
     }
 
+    if (vm.count("namePathStoreLocal")) {
+      std::cout << "namePathStoreLocal was set to "
+                << vm["namePathStoreLocal"].as<std::string>() << ".\n";
+      namePathStoreLocal = vm["namePathStoreLocal"].as<std::string>();
+    }
+
+    if (vm.count("timeRangeLow")) {
+      std::cout << "timeRangeLow was set to "
+                << vm["timeRangeLow"].as<double>() << ".\n";
+      timeRangeLow = vm["timeRangeLow"].as<double>();
+    }
+    if (vm.count("timeRangeHigh")) {
+      std::cout << "timeRangeHigh was set to "
+                << vm["timeRangeHigh"].as<double>() << ".\n";
+      timeRangeHigh = vm["timeRangeHigh"].as<double>();
+    }
+
   } catch (bpo::error& e) {
     std::cerr << "ERROR: " << e.what() << std::endl
               << std::endl;
@@ -159,10 +180,19 @@ int main(int argc, char** argv)
 
     // calibrate the time
     o2::emcal::TimeCalibrationParams TCparams;
-    TCparams = CalibExtractor.calibrateTime(hCalibInputHist);
+    TCparams = CalibExtractor.calibrateTime(hCalibInputHist, timeRangeLow, timeRangeHigh);
 
     // store parameters in ccdb via emcal calibdb
     std::map<std::string, std::string> metadata;
     calibdb.storeTimeCalibParam(&TCparams, metadata, rangestart, rangeend);
+
+    if (namePathStoreLocal.find(".root") != std::string::npos) {
+      TFile fLocalStorage(namePathStoreLocal.c_str(), "update");
+      fLocalStorage.cd();
+      TH1F* histTCparams = (TH1F*)TCparams.getHistogramRepresentation(false);
+      std::string nameTCHist = "TCParams_" + std::to_string(rangestart) + "_" + std::to_string(rangeend);
+      histTCparams->Write(nameTCHist.c_str(), TObject::kOverwrite);
+      fLocalStorage.Close();
+    }
   }
 }
