Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:47:53

0001 // Original Author:  Jan Kašpar
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  * \brief Provides LHCInfo data necessary for CTPPS reconstruction (and direct simulation).
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   // get input beta* vs. xangle histogram
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   // parse histogram
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);