Merge pull request #28160 from CMSTrackerDPG/apv_dynamical_gain_simul…

…ation_10_6_X

[10.6.X]  Dynamical APV gain simulation (backport)

CalibTracker/SiStripESProducers/plugins/real/SiStripApvSimulationParametersESProducer.cc

@@ -0,0 +1,186 @@
+#include "FWCore/Framework/interface/ESProducer.h"
+#include "FWCore/Framework/interface/EventSetupRecordIntervalFinder.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "CondFormats/DataRecord/interface/SiStripCondDataRecords.h"
+#include "CondFormats/SiStripObjects/interface/SiStripApvSimulationParameters.h"
+#include <fstream>
+#include <boost/range/adaptor/indexed.hpp>
+
+class SiStripApvSimulationParametersESSource : public edm::ESProducer, public edm::EventSetupRecordIntervalFinder {
+public:
+  explicit SiStripApvSimulationParametersESSource(const edm::ParameterSet& conf);
+  ~SiStripApvSimulationParametersESSource() override {}
+
+  void setIntervalFor(const edm::eventsetup::EventSetupRecordKey&,
+                      const edm::IOVSyncValue& iov,
+                      edm::ValidityInterval& iValidity) override;
+
+  std::unique_ptr<SiStripApvSimulationParameters> produce(const SiStripApvSimulationParametersRcd& record);
+
+private:
+  std::vector<edm::FileInPath> baselineFiles_TOB_;
+  std::vector<edm::FileInPath> baselineFiles_TIB_;
+  std::vector<edm::FileInPath> baselineFiles_TID_;
+  std::vector<edm::FileInPath> baselineFiles_TEC_;
+  unsigned int baseline_nBins_;
+  float baseline_min_;
+  float baseline_max_;
+  std::vector<float> puBinEdges_;
+  std::vector<float> zBinEdges_;
+  std::vector<float> rBinEdgesTID_;
+  std::vector<float> rBinEdgesTEC_;
+
+  SiStripApvSimulationParameters::LayerParameters makeLayerParameters(const std::string& apvBaselinesFileName,
+                                                                      const std::vector<float>& rzBinEdges) const;
+};
+
+SiStripApvSimulationParametersESSource::SiStripApvSimulationParametersESSource(const edm::ParameterSet& conf)
+    : baseline_nBins_(conf.getUntrackedParameter<unsigned int>("apvBaselines_nBinsPerBaseline")),
+      baseline_min_(conf.getUntrackedParameter<double>("apvBaselines_minBaseline")),
+      baseline_max_(conf.getUntrackedParameter<double>("apvBaselines_maxBaseline")) {
+  setWhatProduced(this);
+  findingRecord<SiStripApvSimulationParametersRcd>();
+  for (const auto x : conf.getUntrackedParameter<std::vector<double>>("apvBaselines_puBinEdges")) {
+    puBinEdges_.push_back(x);
+  }
+  for (const auto x : conf.getUntrackedParameter<std::vector<double>>("apvBaselines_zBinEdges")) {
+    zBinEdges_.push_back(x);
+  }
+  for (const auto x : conf.getUntrackedParameter<std::vector<double>>("apvBaselines_rBinEdges_TID")) {
+    rBinEdgesTID_.push_back(x);
+  }
+  for (const auto x : conf.getUntrackedParameter<std::vector<double>>("apvBaselines_rBinEdges_TEC")) {
+    rBinEdgesTEC_.push_back(x);
+  }
+  baselineFiles_TIB_ = {conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tib1"),
+                        conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tib2"),
+                        conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tib3"),
+                        conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tib4")};
+  baselineFiles_TOB_ = {conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tob1"),
+                        conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tob2"),
+                        conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tob3"),
+                        conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tob4"),
+                        conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tob5"),
+                        conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tob6")};
+
+  if(rBinEdgesTID_.size()>0){
+    baselineFiles_TID_ = {conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tid1"),
+			  conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tid2"),
+			  conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tid3")};
+  }
+
+  if( rBinEdgesTEC_.size()>0){
+    baselineFiles_TEC_ = {conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tec1"),
+			  conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tec2"),
+			  conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tec3"),
+			  conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tec4"),
+			  conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tec5"),
+			  conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tec6"),
+			  conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tec7"),
+			  conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tec8"),
+			  conf.getUntrackedParameter<edm::FileInPath>("apvBaselinesFile_tec9")};
+  }
+}
+
+void SiStripApvSimulationParametersESSource::setIntervalFor(const edm::eventsetup::EventSetupRecordKey&,
+    const edm::IOVSyncValue& iov,
+    edm::ValidityInterval& iValidity) {
+  iValidity = edm::ValidityInterval{iov.beginOfTime(), iov.endOfTime()};
+}
+
+SiStripApvSimulationParameters::LayerParameters SiStripApvSimulationParametersESSource::makeLayerParameters(
+    const std::string& apvBaselinesFileName, const std::vector<float>& rzBinEdges) const {
+  // Prepare histograms
+  unsigned int nZBins = rzBinEdges.size();
+  unsigned int nPUBins = puBinEdges_.size();
+
+  if (nPUBins == 0 || nZBins == 0 || baseline_nBins_ == 0) {
+    throw cms::Exception("MissingInput") << "The parameters for the APV simulation are not correctly configured\n";
+  }
+  std::vector<float> baselineBinEdges{};
+  const auto baseline_binWidth = (baseline_max_ - baseline_min_) / baseline_nBins_;
+  for (unsigned i{0}; i != baseline_nBins_; ++i) {
+    baselineBinEdges.push_back(baseline_min_ + i * baseline_binWidth);
+  }
+  baselineBinEdges.push_back(baseline_max_);
+
+  SiStripApvSimulationParameters::LayerParameters layerParams{baselineBinEdges, puBinEdges_, rzBinEdges};
+
+  // Read apv baselines from text files
+  std::vector<double> theAPVBaselines;
+  std::ifstream apvBaselineFile(apvBaselinesFileName.c_str());
+  if (!apvBaselineFile.good()) {
+    throw cms::Exception("FileError") << "Problem opening APV baselines file: " << apvBaselinesFileName;
+  }
+  std::string line;
+  while (std::getline(apvBaselineFile, line)) {
+    if (!line.empty()) {
+      std::istringstream lStr{line};
+      double value;
+      while (lStr >> value) {
+        theAPVBaselines.push_back(value);
+      }
+    }
+  }
+  if (theAPVBaselines.empty()) {
+    throw cms::Exception("WrongAPVBaselines")
+        << "Problem reading from APV baselines file " << apvBaselinesFileName << ": no values read in";
+  }
+
+  if (theAPVBaselines.size() != nZBins * nPUBins * baseline_nBins_) {
+    throw cms::Exception("WrongAPVBaselines") << "Problem reading from APV baselines file " << apvBaselinesFileName
+                                              << ": number of baselines read different to that expected i.e. nZBins * "
+                                                 "nPUBins * apvBaselines_nBinsPerBaseline_";
+  }
+
+  // Put baselines into histograms
+  for (auto const& apvBaseline : theAPVBaselines | boost::adaptors::indexed(0)) {
+    unsigned int binInCurrentHistogram = apvBaseline.index() % baseline_nBins_ + 1;
+    unsigned int binInZ = int(apvBaseline.index()) / (nPUBins * baseline_nBins_);
+    unsigned int binInPU = int(apvBaseline.index() - binInZ * (nPUBins)*baseline_nBins_) / baseline_nBins_;
+
+    layerParams.setBinContent(binInCurrentHistogram, binInPU+1, binInZ+1, apvBaseline.value());
+  }
+
+  return layerParams;
+}
+
+std::unique_ptr<SiStripApvSimulationParameters> SiStripApvSimulationParametersESSource::produce(
+    const SiStripApvSimulationParametersRcd& record) {
+  auto apvSimParams = std::make_unique<SiStripApvSimulationParameters>(
+      baselineFiles_TIB_.size(), baselineFiles_TOB_.size(), baselineFiles_TID_.size(), baselineFiles_TEC_.size());
+  for (unsigned int i{0}; i != baselineFiles_TIB_.size(); ++i) {
+
+    if (!apvSimParams->putTIB(i + 1, makeLayerParameters(baselineFiles_TIB_[i].fullPath(), zBinEdges_))) {
+      throw cms::Exception("SiStripApvSimulationParameters") << "Could not add parameters for TIB layer " << (i + 1);
+    } else {
+      LogDebug("SiStripApvSimulationParameters") << "Added parameters for TIB layer " << (i+1);
+    }
+  }
+  for (unsigned int i{0}; i != baselineFiles_TOB_.size(); ++i) {
+    if (!apvSimParams->putTOB(i + 1, makeLayerParameters(baselineFiles_TOB_[i].fullPath(), zBinEdges_))) {
+      throw cms::Exception("SiStripApvSimulationParameters") << "Could not add parameters for TOB layer " << (i + 1);
+    } else {
+      LogDebug("SiStripApvSimulationParameters") << "Added parameters for TOB layer " << (i+1);
+    }
+  }
+  for (unsigned int i{0}; i != baselineFiles_TID_.size(); ++i) {
+    if (!apvSimParams->putTID(i + 1, makeLayerParameters(baselineFiles_TID_[i].fullPath(), rBinEdgesTID_))) {
+      throw cms::Exception("SiStripApvSimulationParameters") << "Could not add parameters for TID wheel " << (i + 1);
+    } else {
+      LogDebug("SiStripApvSimulationParameters") << "Added parameters for TID wheel " << (i + 1);
+    }
+  }
+  for (unsigned int i{0}; i != baselineFiles_TEC_.size(); ++i) {
+    if (!apvSimParams->putTEC(i + 1, makeLayerParameters(baselineFiles_TEC_[i].fullPath(), rBinEdgesTEC_))) {
+      throw cms::Exception("SiStripApvSimulationParameters") << "Could not add parameters for TEC wheel " << (i + 1);
+    } else {
+      LogDebug("SiStripApvSimulationParameters") << "Added parameters for TEC wheel " << (i + 1);
+    }
+  }
+  return apvSimParams;
+}
+
+#include "FWCore/Framework/interface/SourceFactory.h"
+DEFINE_FWK_EVENTSETUP_SOURCE(SiStripApvSimulationParametersESSource);

CondCore/SiStripPlugins/plugins/plugin.cc

@@ -65,6 +65,16 @@ REGISTER_PLUGIN(SiStripSummaryRcd,SiStripSummary);
 #include "CondFormats/SiStripObjects/interface/SiStripConfObject.h"
 REGISTER_PLUGIN(SiStripConfObjectRcd, SiStripConfObject);
 
+#include "CondFormats/SiStripObjects/interface/SiStripApvSimulationParameters.h"
+namespace {
+  struct initializeApvSimulationParameters {
+    void operator()(SiStripApvSimulationParameters& param) { param.calculateIntegrals(); }
+  };
+}  // namespace
+REGISTER_PLUGIN_INIT(SiStripApvSimulationParametersRcd,
+                     SiStripApvSimulationParameters,
+                     initializeApvSimulationParameters);
+
 #include "CondFormats/SiStripObjects/interface/Phase2TrackerCabling.h"
 namespace {
   struct initializeCabling {void operator()(Phase2TrackerCabling& c){ c.initializeCabling();}};

CondCore/Utilities/plugins/Module_2XML.cc

@@ -250,6 +250,7 @@ PAYLOAD_2XML_MODULE( pluginUtilities_payload2xml ){
   PAYLOAD_2XML_CLASS( SiPixelQualityProbabilities );
   PAYLOAD_2XML_CLASS( SiPixelTemplateDBObject );
   PAYLOAD_2XML_CLASS( SiStripApvGain );
+  PAYLOAD_2XML_CLASS( SiStripApvSimulationParameters );
   PAYLOAD_2XML_CLASS( SiStripBackPlaneCorrection );
   PAYLOAD_2XML_CLASS( SiStripBadStrip );
   PAYLOAD_2XML_CLASS( SiStripConfObject );

CondCore/Utilities/src/CondDBFetch.cc

@@ -283,6 +283,7 @@ namespace cond {
       FETCH_PAYLOAD_CASE( SiPixelTemplateDBObject )
       FETCH_PAYLOAD_CASE( SiPixel2DTemplateDBObject )
       FETCH_PAYLOAD_CASE( SiStripApvGain )
+      FETCH_PAYLOAD_CASE( SiStripApvSimulationParameters )
       FETCH_PAYLOAD_CASE( SiStripBackPlaneCorrection )
       FETCH_PAYLOAD_CASE( SiStripBadStrip )
       FETCH_PAYLOAD_CASE( SiStripConfObject )

CondCore/Utilities/src/CondDBImport.cc

@@ -304,6 +304,7 @@ namespace cond {
       IMPORT_PAYLOAD_CASE( SiPixelTemplateDBObject )
       IMPORT_PAYLOAD_CASE( SiPixel2DTemplateDBObject )
       IMPORT_PAYLOAD_CASE( SiStripApvGain )
+      IMPORT_PAYLOAD_CASE( SiStripApvSimulationParameters )
       IMPORT_PAYLOAD_CASE( SiStripBadStrip )
       IMPORT_PAYLOAD_CASE( SiStripBackPlaneCorrection )
       IMPORT_PAYLOAD_CASE( SiStripConfObject )

CondCore/Utilities/src/CondFormats.h

@@ -151,6 +151,7 @@
 #include "CondFormats/SiStripObjects/interface/SiStripFedCabling.h"
 #include "CondFormats/SiStripObjects/interface/SiStripThreshold.h"
 #include "CondFormats/SiStripObjects/interface/SiStripBackPlaneCorrection.h"
+#include "CondFormats/SiStripObjects/interface/SiStripApvSimulationParameters.h"
 #include "CondFormats/CSCObjects/interface/CSCDBCrosstalk.h"
 #include "CondFormats/CSCObjects/interface/CSCDBL1TPParameters.h"
 #include "CondFormats/CastorObjects/interface/CastorChannelQuality.h"

CondFormats/DataRecord/interface/SiStripCondDataRecords.h

@@ -46,6 +46,10 @@ class SiStripClusterThresholdRcd : public edm::eventsetup::EventSetupRecordImple
 /*Record for the configuration object*/
 class SiStripConfObjectRcd : public edm::eventsetup::EventSetupRecordImplementation<SiStripConfObjectRcd> {};
 
+/*Record for the APV simulation parameters*/
+class SiStripApvSimulationParametersRcd
+    : public edm::eventsetup::EventSetupRecordImplementation<SiStripApvSimulationParametersRcd> {};
+
 /*Records for upgrade */
 class Phase2TrackerCablingRcd : public edm::eventsetup::EventSetupRecordImplementation<Phase2TrackerCablingRcd> {};
 

CondFormats/DataRecord/src/SiStripCondDataRecords.cc

@@ -39,5 +39,7 @@ EVENTSETUP_RECORD_REG(SiStripClusterThresholdRcd);
 
 EVENTSETUP_RECORD_REG(SiStripConfObjectRcd);
 
+EVENTSETUP_RECORD_REG(SiStripApvSimulationParametersRcd);
+
 EVENTSETUP_RECORD_REG(Phase2TrackerCablingRcd);
 

CondFormats/PhysicsToolsObjects/interface/Histogram3D.icc

@@ -59,7 +59,7 @@ Histogram3D<Value_t, AxisX_t, AxisY_t, AxisZ_t>::Histogram3D(
 	strideX(binULimitsX.size() + 1),
         strideY(binULimitsY.size() + 1),
 	binULimitsX(binULimitsX), binULimitsY(binULimitsY), binULimitsZ(binULimitsZ),
-	binValues((binULimitsY.size() + 1) * strideX * strideY),
+	binValues((binULimitsZ.size() + 1) * strideX * strideY),
 	limitsX(AxisX_t(), AxisX_t()), limitsY(AxisY_t(), AxisY_t()), limitsZ(AxisZ_t(), AxisZ_t()),
 	total(Value_t()),
 	sliceTotal(binULimitsZ.size() + 1), rowTotal(binULimitsY.size() + 1), columnTotal(binULimitsX.size() + 1)

CondFormats/SiStripObjects/BuildFile.xml

@@ -1,5 +1,6 @@
 <use   name="CondFormats/Common"/>
 <use   name="CondFormats/Serialization"/>
+<use   name="CondFormats/PhysicsToolsObjects"/>
 <use   name="DataFormats/TrackerCommon"/>
 <use   name="DataFormats/SiStripCommon"/>
 <use   name="DataFormats/SiStripDetId"/>

CondFormats/SiStripObjects/interface/SiStripApvSimulationParameters.h

@@ -0,0 +1,77 @@
+#ifndef SiStripApvSimulationParameters_h
+#define SiStripApvSimulationParameters_h
+
+#include "CondFormats/Serialization/interface/Serializable.h"
+
+#include <vector>
+#include "CondFormats/PhysicsToolsObjects/interface/Histogram2D.h"
+#include "CondFormats/PhysicsToolsObjects/interface/Histogram3D.h"
+
+namespace CLHEP {
+  class HepRandomEngine;
+}
+
+/**
+ * Stores a histogram binned in PU, z and baseline voltage, for every barrel layer of the strip tracker
+ */
+class SiStripApvSimulationParameters {
+public:
+  using layerid = unsigned int;
+  using LayerParameters = PhysicsTools::Calibration::HistogramF3D;
+
+  SiStripApvSimulationParameters(layerid nTIB, layerid nTOB, layerid nTID, layerid nTEC)
+      : m_nTIB(nTIB), m_nTOB(nTOB), m_nTID(nTID), m_nTEC(nTEC) {
+    m_barrelParam.resize(m_nTIB + m_nTOB);
+    m_barrelParam_xInt.resize(m_nTIB + m_nTOB);
+    m_endcapParam.resize(m_nTID + m_nTEC);
+    m_endcapParam_xInt.resize(m_nTID + m_nTEC);
+  }
+  SiStripApvSimulationParameters() {}
+  ~SiStripApvSimulationParameters() {}
+
+  void calculateIntegrals();  // make sure integrals have been calculated
+
+  bool putTIB(layerid layer, const LayerParameters& params) { return putTIB(layer, LayerParameters(params)); }
+  bool putTIB(layerid layer, LayerParameters&& params);
+
+  bool putTOB(layerid layer, const LayerParameters& params) { return putTOB(layer, LayerParameters(params)); }
+  bool putTOB(layerid layer, LayerParameters&& params);
+
+  bool putTID(layerid wheel, const LayerParameters& params) { return putTID(wheel, LayerParameters(params)); }
+  bool putTID(layerid wheel, LayerParameters&& params);
+
+  bool putTEC(layerid wheel, const LayerParameters& params) { return putTEC(wheel, LayerParameters(params)); }
+  bool putTEC(layerid wheel, LayerParameters&& params);
+
+  const LayerParameters& getTIB(layerid layer) const { return m_barrelParam[layer - 1]; }
+  const LayerParameters& getTOB(layerid layer) const { return m_barrelParam[m_nTIB + layer - 1]; }
+  const LayerParameters& getTID(layerid wheel) const { return m_endcapParam[wheel - 1]; }
+  const LayerParameters& getTEC(layerid wheel) const { return m_endcapParam[m_nTID + wheel - 1]; }
+
+  float sampleTIB(layerid layer, float z, float pu, CLHEP::HepRandomEngine* engine) const {
+    return sampleBarrel(layer - 1, z, pu, engine);
+  }
+  float sampleTOB(layerid layer, float z, float pu, CLHEP::HepRandomEngine* engine) const {
+    return sampleBarrel(m_nTIB + layer - 1, z, pu, engine);
+  };
+  float sampleTID(layerid wheel, float r, float pu, CLHEP::HepRandomEngine* engine) const {
+    return sampleEndcap(wheel - 1, r, pu, engine);
+  }
+  float sampleTEC(layerid wheel, float r, float pu, CLHEP::HepRandomEngine* engine) const {
+    return sampleEndcap(m_nTID + wheel - 1, r, pu, engine);
+  }
+
+private:
+  layerid m_nTIB, m_nTOB, m_nTID, m_nTEC;
+  std::vector<PhysicsTools::Calibration::HistogramF3D> m_barrelParam;
+  std::vector<PhysicsTools::Calibration::HistogramF2D> m_barrelParam_xInt;
+  std::vector<PhysicsTools::Calibration::HistogramF3D> m_endcapParam;
+  std::vector<PhysicsTools::Calibration::HistogramF2D> m_endcapParam_xInt;
+
+  float sampleBarrel(layerid layerIdx, float z, float pu, CLHEP::HepRandomEngine* engine) const;
+  float sampleEndcap(layerid wheelIdx, float r, float pu, CLHEP::HepRandomEngine* engine) const;
+
+  COND_SERIALIZABLE;
+};
+
+#endif