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