File indexing completed on 2024-04-06 11:58:35
0001
0002
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004 #include "FWCore/Framework/interface/SourceFactory.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/Framework/interface/ESProducer.h"
0008 #include "FWCore/Framework/interface/EventSetupRecordIntervalFinder.h"
0009 #include "FWCore/Framework/interface/ESProducts.h"
0010 #include "FWCore/ServiceRegistry/interface/Service.h"
0011
0012 #include "CondFormats/RunInfo/interface/LHCInfo.h"
0013 #include "CondFormats/DataRecord/interface/LHCInfoRcd.h"
0014
0015 #include "CLHEP/Random/RandFlat.h"
0016 #include "CLHEP/Random/JamesRandom.h"
0017
0018 #include "TFile.h"
0019 #include "TH2D.h"
0020
0021
0022
0023
0024
0025
0026 class CTPPSLHCInfoRandomXangleESSource : public edm::ESProducer, public edm::EventSetupRecordIntervalFinder {
0027 public:
0028 CTPPSLHCInfoRandomXangleESSource(const edm::ParameterSet &);
0029 edm::ESProducts<std::unique_ptr<LHCInfo>> produce(const LHCInfoRcd &);
0030 static void fillDescriptions(edm::ConfigurationDescriptions &);
0031
0032 private:
0033 void setIntervalFor(const edm::eventsetup::EventSetupRecordKey &,
0034 const edm::IOVSyncValue &,
0035 edm::ValidityInterval &) override;
0036
0037 std::string m_label;
0038
0039 unsigned int m_generateEveryNEvents;
0040
0041 double m_beamEnergy;
0042
0043 std::unique_ptr<CLHEP::HepRandomEngine> m_engine;
0044
0045 struct BinData {
0046 double min, max;
0047 double xangle, betaStar;
0048 };
0049
0050 std::vector<BinData> xangleBetaStarBins;
0051 };
0052
0053
0054
0055
0056 CTPPSLHCInfoRandomXangleESSource::CTPPSLHCInfoRandomXangleESSource(const edm::ParameterSet &conf)
0057 : m_label(conf.getParameter<std::string>("label")),
0058
0059 m_generateEveryNEvents(conf.getParameter<unsigned int>("generateEveryNEvents")),
0060
0061 m_beamEnergy(conf.getParameter<double>("beamEnergy")),
0062
0063 m_engine(new CLHEP::HepJamesRandom(conf.getParameter<unsigned int>("seed"))) {
0064
0065 const auto &xangleBetaStarHistogramFile = conf.getParameter<std::string>("xangleBetaStarHistogramFile");
0066 const auto &xangleBetaStarHistogramObject = conf.getParameter<std::string>("xangleBetaStarHistogramObject");
0067
0068 edm::FileInPath fip(xangleBetaStarHistogramFile);
0069 std::unique_ptr<TFile> f_in(TFile::Open(fip.fullPath().c_str()));
0070 if (!f_in)
0071 throw cms::Exception("PPS") << "Cannot open input file '" << xangleBetaStarHistogramFile << "'.";
0072
0073 TH2D *h_xangle_beta_star = (TH2D *)f_in->Get(xangleBetaStarHistogramObject.c_str());
0074 if (!h_xangle_beta_star)
0075 throw cms::Exception("PPS") << "Cannot load input object '" << xangleBetaStarHistogramObject << "'.";
0076
0077
0078 double sum = 0.;
0079 for (int x = 1; x <= h_xangle_beta_star->GetNbinsX(); ++x) {
0080 for (int y = 1; y <= h_xangle_beta_star->GetNbinsY(); ++y)
0081 sum += h_xangle_beta_star->GetBinContent(x, y);
0082 }
0083
0084 double cw = 0;
0085 for (int x = 1; x <= h_xangle_beta_star->GetNbinsX(); ++x) {
0086 for (int y = 1; y <= h_xangle_beta_star->GetNbinsY(); ++y) {
0087 const double c = h_xangle_beta_star->GetBinContent(x, y);
0088 const double xangle = h_xangle_beta_star->GetXaxis()->GetBinCenter(x);
0089 const double betaStar = h_xangle_beta_star->GetYaxis()->GetBinCenter(y);
0090
0091 if (c > 0.) {
0092 const double rc = c / sum;
0093 xangleBetaStarBins.push_back({cw, cw + rc, xangle, betaStar});
0094 cw += rc;
0095 }
0096 }
0097 }
0098
0099 setWhatProduced(this, m_label);
0100 findingRecord<LHCInfoRcd>();
0101 }
0102
0103
0104
0105 void CTPPSLHCInfoRandomXangleESSource::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0106 edm::ParameterSetDescription desc;
0107
0108 desc.add<std::string>("label", "")->setComment("label of the LHCInfo record");
0109
0110 desc.add<unsigned int>("seed", 1)->setComment("random seed");
0111
0112 desc.add<unsigned int>("generateEveryNEvents", 1)->setComment("how often to generate new xangle");
0113
0114 desc.add<std::string>("xangleBetaStarHistogramFile", "")->setComment("ROOT file with xangle distribution");
0115 desc.add<std::string>("xangleBetaStarHistogramObject", "")->setComment("xangle distribution object in the ROOT file");
0116
0117 desc.add<double>("beamEnergy", 0.)->setComment("beam energy");
0118
0119 descriptions.add("ctppsLHCInfoRandomXangleESSource", desc);
0120 }
0121
0122
0123
0124 void CTPPSLHCInfoRandomXangleESSource::setIntervalFor(const edm::eventsetup::EventSetupRecordKey &key,
0125 const edm::IOVSyncValue &iosv,
0126 edm::ValidityInterval &oValidity) {
0127 edm::EventID beginEvent = iosv.eventID();
0128 edm::EventID endEvent(beginEvent.run(), beginEvent.luminosityBlock(), beginEvent.event() + m_generateEveryNEvents);
0129 oValidity = edm::ValidityInterval(edm::IOVSyncValue(beginEvent), edm::IOVSyncValue(endEvent));
0130 }
0131
0132
0133
0134 edm::ESProducts<std::unique_ptr<LHCInfo>> CTPPSLHCInfoRandomXangleESSource::produce(const LHCInfoRcd &) {
0135 auto output = std::make_unique<LHCInfo>();
0136
0137 double xangle = 0., betaStar = 0.;
0138 const double u = CLHEP::RandFlat::shoot(m_engine.get(), 0., 1.);
0139 for (const auto &d : xangleBetaStarBins) {
0140 if (d.min <= u && u <= d.max) {
0141 xangle = d.xangle;
0142 betaStar = d.betaStar;
0143 break;
0144 }
0145 }
0146
0147 output->setEnergy(m_beamEnergy);
0148 output->setCrossingAngle(xangle);
0149 output->setBetaStar(betaStar);
0150
0151 return edm::es::products(std::move(output));
0152 }
0153
0154
0155
0156 DEFINE_FWK_EVENTSETUP_SOURCE(CTPPSLHCInfoRandomXangleESSource);
