CTPPSLHCInfoRandomXangleESSource.cc

Go to the documentation of this file.

// Original Author:  Jan Kašpar
 
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/SourceFactory.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/ESProducer.h"
#include "FWCore/Framework/interface/EventSetupRecordIntervalFinder.h"
#include "FWCore/Framework/interface/ESProducts.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
 
#include "CondFormats/RunInfo/interface/LHCInfo.h"
#include "CondFormats/DataRecord/interface/LHCInfoRcd.h"
 
#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/JamesRandom.h"
 
#include "TFile.h"
#include "TH1D.h"
 
//----------------------------------------------------------------------------------------------------
 
class CTPPSLHCInfoRandomXangleESSource : public edm::ESProducer, public edm::EventSetupRecordIntervalFinder {
public:
  CTPPSLHCInfoRandomXangleESSource(const edm::ParameterSet &);
  edm::ESProducts<std::unique_ptr<LHCInfo>> produce(const LHCInfoRcd &);
  static void fillDescriptions(edm::ConfigurationDescriptions &);
 
private:
  void setIntervalFor(const edm::eventsetup::EventSetupRecordKey &,
                      const edm::IOVSyncValue &,
                      edm::ValidityInterval &) override;
 
  std::string m_label;
 
  unsigned int m_generateEveryNEvents;
 
  double m_beamEnergy;
  double m_betaStar;
 
  std::unique_ptr<CLHEP::HepRandomEngine> m_engine;
 
  struct BinData {
    double min, max, xangle;
  };
 
  std::vector<BinData> binData;
};
 
//----------------------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------------------
 
CTPPSLHCInfoRandomXangleESSource::CTPPSLHCInfoRandomXangleESSource(const edm::ParameterSet &conf)
    : m_label(conf.getParameter<std::string>("label")),
 
      m_generateEveryNEvents(conf.getParameter<unsigned int>("generateEveryNEvents")),
 
      m_beamEnergy(conf.getParameter<double>("beamEnergy")),
      m_betaStar(conf.getParameter<double>("betaStar")),
 
      m_engine(new CLHEP::HepJamesRandom(conf.getParameter<unsigned int>("seed"))) {
  const auto &xangleHistogramFile = conf.getParameter<std::string>("xangleHistogramFile");
  const auto &xangleHistogramObject = conf.getParameter<std::string>("xangleHistogramObject");
 
  TFile *f_in = TFile::Open(xangleHistogramFile.c_str());
  TH1D *h_xangle = (TH1D *)f_in->Get(xangleHistogramObject.c_str());
 
  double s = 0.;
  for (int bi = 1; bi <= h_xangle->GetNbinsX(); ++bi)
    s += h_xangle->GetBinContent(bi);
 
  double cw = 0.;
  for (int bi = 1; bi <= h_xangle->GetNbinsX(); ++bi) {
    double xangle = h_xangle->GetBinCenter(bi);
    double w = h_xangle->GetBinContent(bi) / s;
 
    binData.push_back({cw, cw + w, xangle});
 
    cw += w;
  }
 
  delete f_in;
 
  setWhatProduced(this, m_label);
  findingRecord<LHCInfoRcd>();
}
 
//----------------------------------------------------------------------------------------------------
 
void CTPPSLHCInfoRandomXangleESSource::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
  edm::ParameterSetDescription desc;
 
  desc.add<std::string>("label", "")->setComment("label of the LHCInfo record");
 
  desc.add<unsigned int>("seed", 1)->setComment("random seed");
 
  desc.add<unsigned int>("generateEveryNEvents", 1)->setComment("how often to generate new xangle");
 
  desc.add<std::string>("xangleHistogramFile", "")->setComment("ROOT file with xangle distribution");
  desc.add<std::string>("xangleHistogramObject", "")->setComment("xangle distribution object in the ROOT file");
 
  desc.add<double>("beamEnergy", 0.)->setComment("beam energy");
  desc.add<double>("betaStar", 0.)->setComment("beta*");
 
  descriptions.add("ctppsLHCInfoRandomXangleESSource", desc);
}
 
//----------------------------------------------------------------------------------------------------
 
void CTPPSLHCInfoRandomXangleESSource::setIntervalFor(const edm::eventsetup::EventSetupRecordKey &key,
                                                      const edm::IOVSyncValue &iosv,
                                                      edm::ValidityInterval &oValidity) {
  edm::EventID beginEvent = iosv.eventID();
  edm::EventID endEvent(beginEvent.run(), beginEvent.luminosityBlock(), beginEvent.event() + m_generateEveryNEvents);
  oValidity = edm::ValidityInterval(edm::IOVSyncValue(beginEvent), edm::IOVSyncValue(endEvent));
}
 
//----------------------------------------------------------------------------------------------------
 
edm::ESProducts<std::unique_ptr<LHCInfo>> CTPPSLHCInfoRandomXangleESSource::produce(const LHCInfoRcd &) {
  auto output = std::make_unique<LHCInfo>();
 
  const double u = CLHEP::RandFlat::shoot(m_engine.get(), 0., 1.);
 
  double xangle = 0.;
  for (const auto &d : binData) {
    if (d.min <= u && u <= d.max) {
      xangle = d.xangle;
      break;
    }
  }
 
  output->setEnergy(m_beamEnergy);
  output->setCrossingAngle(xangle);
  output->setBetaStar(m_betaStar);
 
  return edm::es::products(std::move(output));
}
 
//----------------------------------------------------------------------------------------------------
 
DEFINE_FWK_EVENTSETUP_SOURCE(CTPPSLHCInfoRandomXangleESSource);