ESProducers/plugins/CTPPSInterpolatedOpticalFunctionsESSource.cc

0001 // authors: Jan Kaspar (jan.kaspar@gmail.com)
0002
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/ModuleFactory.h"
0007 #include "FWCore/Framework/interface/ESProducer.h"
0008
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include "CondFormats/DataRecord/interface/LHCInfoRcd.h"
0012 #include "CondFormats/DataRecord/interface/CTPPSOpticsRcd.h"
0013 #include "CondFormats/DataRecord/interface/CTPPSInterpolatedOpticsRcd.h"
0014
0015 #include "CondFormats/RunInfo/interface/LHCInfo.h"
0016 #include "CondFormats/PPSObjects/interface/LHCOpticalFunctionsSetCollection.h"
0017 #include "CondFormats/PPSObjects/interface/LHCInterpolatedOpticalFunctionsSetCollection.h"
0018
0019 class CTPPSInterpolatedOpticalFunctionsESSource : public edm::ESProducer {
0020 public:
0021   CTPPSInterpolatedOpticalFunctionsESSource(const edm::ParameterSet &);
0022   ~CTPPSInterpolatedOpticalFunctionsESSource() override {}
0023
0024   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0025
0026   std::shared_ptr<LHCInterpolatedOpticalFunctionsSetCollection> produce(const CTPPSInterpolatedOpticsRcd &);
0027
0028 private:
0029   edm::ESGetToken<LHCOpticalFunctionsSetCollection, CTPPSOpticsRcd> opticsToken_;
0030   edm::ESGetToken<LHCInfo, LHCInfoRcd> lhcInfoToken_;
0031   std::shared_ptr<LHCInterpolatedOpticalFunctionsSetCollection> currentData_;
0032   float currentCrossingAngle_;
0033   bool currentDataValid_;
0034 };
0035
0036 //----------------------------------------------------------------------------------------------------
0037 //----------------------------------------------------------------------------------------------------
0038
0039 CTPPSInterpolatedOpticalFunctionsESSource::CTPPSInterpolatedOpticalFunctionsESSource(const edm::ParameterSet &iConfig)
0040     : currentCrossingAngle_(-1.), currentDataValid_(false) {
0041   auto cc = setWhatProduced(this, iConfig.getParameter<std::string>("opticsLabel"));
0042   opticsToken_ = cc.consumes(edm::ESInputTag("", iConfig.getParameter<std::string>("opticsLabel")));
0043   lhcInfoToken_ = cc.consumes(edm::ESInputTag("", iConfig.getParameter<std::string>("lhcInfoLabel")));
0044 }
0045
0046 //----------------------------------------------------------------------------------------------------
0047
0048 void CTPPSInterpolatedOpticalFunctionsESSource::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0049   edm::ParameterSetDescription desc;
0050
0051   desc.add<std::string>("lhcInfoLabel", "")->setComment("label of the LHCInfo record");
0052   desc.add<std::string>("opticsLabel", "")->setComment("label of the optics records");
0053
0054   descriptions.add("ctppsInterpolatedOpticalFunctionsESSource", desc);
0055 }
0056
0057 //----------------------------------------------------------------------------------------------------
0058
0059 std::shared_ptr<LHCInterpolatedOpticalFunctionsSetCollection> CTPPSInterpolatedOpticalFunctionsESSource::produce(
0060     const CTPPSInterpolatedOpticsRcd &iRecord) {
0061   // get the input data
0062   LHCOpticalFunctionsSetCollection const &ofColl = iRecord.get(opticsToken_);
0063
0064   LHCInfo const &lhcInfo = iRecord.get(lhcInfoToken_);
0065
0066   // is there anything to do?
0067   if (currentDataValid_ && lhcInfo.crossingAngle() == currentCrossingAngle_)
0068     return currentData_;
0069
0070   // is crossing angle reasonable (LHCInfo is correctly filled in DB)?
0071   if (lhcInfo.crossingAngle() == 0.) {
0072     edm::LogInfo("CTPPSInterpolatedOpticalFunctionsESSource")
0073         << "Invalid crossing angle, no optical functions produced.";
0074
0075     currentDataValid_ = false;
0076     currentCrossingAngle_ = -1;
0077     currentData_ = std::make_shared<LHCInterpolatedOpticalFunctionsSetCollection>();
0078
0079     return currentData_;
0080   }
0081
0082   // set new crossing angle
0083   currentCrossingAngle_ = lhcInfo.crossingAngle();
0084   edm::LogInfo("CTPPSInterpolatedOpticalFunctionsESSource")
0085       << "Crossing angle has changed to " << currentCrossingAngle_ << ".";
0086
0087   // is input optics available ?
0088   if (ofColl.empty()) {
0089     edm::LogInfo("CTPPSInterpolatedOpticalFunctionsESSource")
0090         << "No input optics available, no optical functions produced.";
0091
0092     currentDataValid_ = false;
0093     currentData_ = std::make_shared<LHCInterpolatedOpticalFunctionsSetCollection>();
0094
0095     return currentData_;
0096   }
0097
0098   // regular case with single-xangle input
0099   if (ofColl.size() == 1) {
0100     const auto &it = ofColl.begin();
0101
0102     // does the input xangle correspond to the actual one?
0103     if (fabs(currentCrossingAngle_ - it->first) > 1e-6)
0104       throw cms::Exception("CTPPSInterpolatedOpticalFunctionsESSource")
0105           << "Cannot interpolate: input given only for xangle " << it->first << " while interpolation requested for "
0106           << currentCrossingAngle_ << ".";
0107
0108     currentData_ = std::make_shared<LHCInterpolatedOpticalFunctionsSetCollection>();
0109     for (const auto &rp_p : it->second) {
0110       const auto rpId = rp_p.first;
0111       LHCInterpolatedOpticalFunctionsSet iof(rp_p.second);
0112       iof.initializeSplines();
0113       currentData_->emplace(rpId, std::move(iof));
0114     }
0115
0116     currentDataValid_ = true;
0117   }
0118
0119   // regular case with multi-xangle input
0120   if (ofColl.size() > 1) {
0121     // find the closest xangle points for interpolation
0122     auto it1 = ofColl.begin();
0123     auto it2 = std::next(it1);
0124
0125     if (currentCrossingAngle_ > it1->first) {
0126       for (; it1 != ofColl.end(); ++it1) {
0127         it2 = std::next(it1);
0128
0129         if (it2 == ofColl.end()) {
0130           it2 = it1;
0131           it1 = std::prev(it1);
0132           break;
0133         }
0134
0135         if (it1->first <= currentCrossingAngle_ && currentCrossingAngle_ < it2->first)
0136           break;
0137       }
0138     }
0139
0140     const auto &xangle1 = it1->first;
0141     const auto &xangle2 = it2->first;
0142
0143     const auto &ofs1 = it1->second;
0144     const auto &ofs2 = it2->second;
0145
0146     // do the interpoaltion RP by RP
0147     currentData_ = std::make_shared<LHCInterpolatedOpticalFunctionsSetCollection>();
0148     for (const auto &rp_p : ofs1) {
0149       const auto rpId = rp_p.first;
0150       const auto &rp_it2 = ofs2.find(rpId);
0151       if (rp_it2 == ofs2.end())
0152         throw cms::Exception("CTPPSInterpolatedOpticalFunctionsESSource") << "RP mismatch between ofs1 and ofs2.";
0153
0154       const auto &of1 = rp_p.second;
0155       const auto &of2 = rp_it2->second;
0156
0157       const size_t num_xi_vals1 = of1.getXiValues().size();
0158       const size_t num_xi_vals2 = of2.getXiValues().size();
0159
0160       if (num_xi_vals1 != num_xi_vals2)
0161         throw cms::Exception("CTPPSInterpolatedOpticalFunctionsESSource") << "Size mismatch between ofs1 and ofs2.";
0162
0163       const size_t num_xi_vals = num_xi_vals1;
0164
0165       LHCInterpolatedOpticalFunctionsSet iof;
0166       iof.m_z = of1.getScoringPlaneZ();
0167       iof.m_fcn_values.resize(LHCInterpolatedOpticalFunctionsSet::nFunctions);
0168       iof.m_xi_values.resize(num_xi_vals);
0169
0170       for (size_t fi = 0; fi < of1.getFcnValues().size(); ++fi) {
0171         iof.m_fcn_values[fi].resize(num_xi_vals);
0172
0173         for (size_t pi = 0; pi < num_xi_vals; ++pi) {
0174           double xi = of1.getXiValues()[pi];
0175           double xi_control = of2.getXiValues()[pi];
0176
0177           if (fabs(xi - xi_control) > 1e-6)
0178             throw cms::Exception("CTPPSInterpolatedOpticalFunctionsESSource") << "Xi mismatch between ofs1 and ofs2.";
0179
0180           iof.m_xi_values[pi] = xi;
0181
0182           double v1 = of1.getFcnValues()[fi][pi];
0183           double v2 = of2.getFcnValues()[fi][pi];
0184           iof.m_fcn_values[fi][pi] = v1 + (v2 - v1) / (xangle2 - xangle1) * (currentCrossingAngle_ - xangle1);
0185         }
0186       }
0187
0188       iof.initializeSplines();
0189
0190       currentDataValid_ = true;
0191       currentData_->emplace(rpId, std::move(iof));
0192     }
0193   }
0194
0195   return currentData_;
0196 }
0197
0198 //----------------------------------------------------------------------------------------------------
0199
0200 DEFINE_FWK_EVENTSETUP_MODULE(CTPPSInterpolatedOpticalFunctionsESSource);