Base/src/BMixingModule.cc

0001 // File: BMixingModule.cc
0002 // Description:  see BMixingModule.h
0003 // Author:  Ursula Berthon, LLR Palaiseau, Bill Tanenbaum
0004 //
0005 //--------------------------------------------
0006
0007 #include "Mixing/Base/interface/BMixingModule.h"
0008 #include "FWCore/Version/interface/GetReleaseVersion.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/EventPrincipal.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/Utilities/interface/ExceptionCollector.h"
0015 #include "DataFormats/Common/interface/Handle.h"
0016
0017 #include "TFile.h"
0018 #include "TH1F.h"
0019 #include <iostream>
0020 #include <memory>
0021
0022 const unsigned int edm::BMixingModule::maxNbSources_ = 4;
0023
0024 namespace {
0025   std::shared_ptr<edm::PileUpConfig> maybeConfigPileUp(
0026       edm::ParameterSet const& ps, std::string sourceName, const int minb, const int maxb, const bool playback) {
0027     std::shared_ptr<edm::PileUpConfig> pileupconfig;  // value to be returned
0028     // Make sure we have a parameter named 'sourceName'
0029     if (ps.exists(sourceName)) {
0030       // We have the parameter
0031       // and if we have either averageNumber or cfg by luminosity... make the PileUp
0032       double averageNumber;
0033       std::string histoFileName = " ";
0034       std::string histoName = " ";
0035       std::unique_ptr<TH1F> h(new TH1F("h", "h", 10, 0, 10));
0036       std::vector<int> dataProbFunctionVar;
0037       std::vector<double> dataProb;
0038
0039       const edm::ParameterSet& psin = ps.getParameter<edm::ParameterSet>(sourceName);
0040       std::string type_ = psin.getParameter<std::string>("type");
0041       if (ps.exists("readDB") && ps.getParameter<bool>("readDB")) {
0042         //in case of DB access, do not try to load anything from the PSet, but wait for beginRun.
0043         edm::LogError("BMixingModule") << "Will read from DB: reset to a dummy PileUp object.";
0044         std::unique_ptr<TH1F> h;
0045         pileupconfig = std::make_shared<edm::PileUpConfig>(sourceName, 0.0, h, playback);
0046         return pileupconfig;
0047       }
0048       if (type_ != "none") {
0049         if (psin.exists("nbPileupEvents")) {
0050           edm::ParameterSet psin_average = psin.getParameter<edm::ParameterSet>("nbPileupEvents");
0051           if (psin_average.exists("averageNumber")) {
0052             averageNumber = psin_average.getParameter<double>("averageNumber");
0053             pileupconfig = std::make_shared<edm::PileUpConfig>(sourceName, averageNumber, h, playback);
0054             edm::LogInfo("MixingModule") << " Created source " << sourceName << " with averageNumber " << averageNumber;
0055           } else if (psin_average.exists("fileName") && psin_average.exists("histoName")) {
0056             std::string histoFileName = psin_average.getUntrackedParameter<std::string>("fileName");
0057             std::string histoName = psin_average.getUntrackedParameter<std::string>("histoName");
0058
0059             std::unique_ptr<TFile> infile(new TFile(histoFileName.c_str()));
0060             std::unique_ptr<TH1F> h((TH1F*)infile->Get(histoName.c_str()));
0061
0062             // Check if the histogram exists
0063             if (!h) {
0064               throw cms::Exception("HistogramNotFound") << " Could not find the histogram " << histoName
0065                                                         << "in the file " << histoFileName << "." << std::endl;
0066             } else {
0067               edm::LogInfo("MixingModule")
0068                   << "Open a root file " << histoFileName << " containing the probability distribution histogram "
0069                   << histoName << std::endl;
0070               edm::LogInfo("MixingModule")
0071                   << "The PileUp number to be added will be chosen randomly from this histogram" << std::endl;
0072             }
0073
0074             // Check if the histogram is normalized
0075             if (std::abs(h->Integral() - 1) > 1.0e-02)
0076               throw cms::Exception("BadHistoDistribution") << "The histogram should be normalized!" << std::endl;
0077
0078             // Get the averageNumber from the histo
0079             averageNumber = h->GetMean();
0080
0081             pileupconfig = std::make_shared<edm::PileUpConfig>(sourceName, averageNumber, h, playback);
0082             edm::LogInfo("MixingModule") << " Created source " << sourceName << " with averageNumber " << averageNumber;
0083
0084           } else if (psin_average.exists("probFunctionVariable") && psin_average.exists("probValue") &&
0085                      psin_average.exists("histoFileName")) {
0086             if (type_ != "probFunction") {
0087               edm::LogError("MisConfiguration")
0088                   << "type is set to: " << type_ << " while parameters implies probFunction; changing.";
0089               type_ = "probFunction";
0090             }
0091
0092             dataProbFunctionVar = psin_average.getParameter<std::vector<int> >("probFunctionVariable");
0093             dataProb = psin_average.getParameter<std::vector<double> >("probValue");
0094             histoFileName = psin_average.getUntrackedParameter<std::string>("histoFileName");
0095
0096             int varSize = (int)dataProbFunctionVar.size();
0097             int probSize = (int)dataProb.size();
0098
0099             if ((dataProbFunctionVar[0] != 0) || (dataProbFunctionVar[varSize - 1] != (varSize - 1)))
0100               throw cms::Exception("BadProbFunction")
0101                   << "Please, check the variables of the probability function! The first variable should be 0 and the "
0102                      "difference between two variables should be 1."
0103                   << std::endl;
0104
0105             // Complete the vector containing the probability  function data
0106             // with the values "0"
0107             if (probSize < varSize) {
0108               edm::LogInfo("MixingModule")
0109                   << " The probability function data will be completed with " << (varSize - probSize) << " values 0.";
0110
0111               for (int i = 0; i < (varSize - probSize); i++)
0112                 dataProb.push_back(0);
0113
0114               probSize = dataProb.size();
0115               edm::LogInfo("MixingModule")
0116                   << " The number of the P(x) data set after adding the values 0 is " << probSize;
0117             }
0118
0119             // Create an histogram with the data from the probability function provided by the user
0120             int xmin = (int)dataProbFunctionVar[0];
0121             int xmax = (int)dataProbFunctionVar[varSize - 1] + 1;  // need upper edge to be one beyond last value
0122             int numBins = varSize;
0123
0124             edm::LogInfo("MixingModule") << "An histogram will be created with " << numBins << " bins in the range ("
0125                                          << xmin << "," << xmax << ")." << std::endl;
0126
0127             std::unique_ptr<TH1F> hprob(
0128                 new TH1F("h", "Histo from the user's probability function", numBins, xmin, xmax));
0129
0130             LogDebug("MixingModule") << "Filling histogram with the following data:" << std::endl;
0131
0132             for (int j = 0; j < numBins; j++) {
0133               LogDebug("MixingModule") << " x = " << dataProbFunctionVar[j] << " P(x) = " << dataProb[j];
0134               hprob->Fill(dataProbFunctionVar[j] + 0.5,
0135                           dataProb[j]);  // assuming integer values for the bins, fill bin centers, not edges
0136             }
0137
0138             // Check if the histogram is normalized
0139             if (std::abs(hprob->Integral() - 1) > 1.0e-02) {
0140               throw cms::Exception("BadProbFunction")
0141                   << "The probability function should be normalized!!! " << std::endl;
0142             }
0143
0144             averageNumber = hprob->GetMean();
0145
0146             // Write the created histogram into a root file
0147             edm::LogInfo("MixingModule")
0148                 << " The histogram created from the x, P(x) values will be written into the root file "
0149                 << histoFileName;
0150
0151             TFile* outfile = new TFile(histoFileName.c_str(), "RECREATE");
0152             hprob->Write();
0153             outfile->Write();
0154             outfile->Close();
0155             outfile->Delete();
0156
0157             pileupconfig = std::make_shared<edm::PileUpConfig>(sourceName, averageNumber, hprob, playback);
0158             edm::LogInfo("MixingModule") << " Created source " << sourceName << " with averageNumber " << averageNumber;
0159           }
0160           //special for pileup input
0161           else if (sourceName == "input" && psin_average.exists("Lumi") && psin_average.exists("sigmaInel")) {
0162             averageNumber = psin_average.getParameter<double>("Lumi") * psin_average.getParameter<double>("sigmaInel") *
0163                             ps.getParameter<int>("bunchspace") / 1000 * 3564. / 2808.;  //FIXME
0164             pileupconfig = std::make_shared<edm::PileUpConfig>(sourceName, averageNumber, h, playback);
0165             edm::LogInfo("MixingModule") << " Created source " << sourceName << " with minBunch,maxBunch " << minb
0166                                          << " " << maxb;
0167             edm::LogInfo("MixingModule") << " Luminosity configuration, average number used is " << averageNumber;
0168           }
0169         }
0170       }
0171     }
0172     return pileupconfig;
0173   }
0174 }  // namespace
0175
0176 namespace edm {
0177
0178   // Constructor
0179   BMixingModule::BMixingModule(const edm::ParameterSet& pset, MixingCache::Config const* globalConf)
0180       : bunchSpace_(globalConf->bunchSpace_),
0181         vertexOffset_(0),
0182         minBunch_(globalConf->minBunch_),
0183         maxBunch_(globalConf->maxBunch_),
0184         mixProdStep1_(pset.getParameter<bool>("mixProdStep1")),
0185         mixProdStep2_(pset.getParameter<bool>("mixProdStep2")),
0186         readDB_(globalConf->configFromDB_),
0187         playback_(globalConf->playback_) {
0188     for (size_t makeIdx = 0; makeIdx < maxNbSources_; makeIdx++) {
0189       if (globalConf->inputConfigs_[makeIdx]) {
0190         const edm::ParameterSet& psin =
0191             pset.getParameter<edm::ParameterSet>(globalConf->inputConfigs_[makeIdx]->sourcename_);
0192         inputSources_.push_back(
0193             std::make_shared<PileUp>(psin, globalConf->inputConfigs_[makeIdx], consumesCollector(), readDB_));
0194         inputSources_.back()->input(makeIdx);
0195       } else {
0196         inputSources_.push_back(nullptr);
0197       }
0198     }
0199   }
0200
0201   // Virtual destructor needed.
0202   BMixingModule::~BMixingModule() { ; }
0203
0204   void BMixingModule::registerLateConsumes(eventsetup::ESRecordsToProductResolverIndices const& iES) {
0205     for (size_t endIdx = 0; endIdx < maxNbSources_; ++endIdx) {
0206       if (inputSources_[endIdx])
0207         inputSources_[endIdx]->beginJob(iES);
0208     }
0209   }
0210
0211   namespace MixingCache {
0212     Config::Config(edm::ParameterSet const& pset, unsigned int maxNbSources)
0213         : bunchSpace_(pset.getParameter<int>("bunchspace")),
0214           minBunch_((pset.getParameter<int>("minBunch") * 25) / bunchSpace_),
0215           maxBunch_((pset.getParameter<int>("maxBunch") * 25) / bunchSpace_),
0216           playback_(pset.getUntrackedParameter<bool>("playback", false)) {
0217       if (playback_) {
0218         //this could be explicitly checked
0219         LogInfo("MixingModule") << " ATTENTION:Mixing will be done in playback mode! \n"
0220                                 << " ATTENTION:Mixing Configuration must be the same as for the original mixing!";
0221       }
0222
0223       // Just for debugging print out.
0224       sourceNames_.push_back("input");
0225       sourceNames_.push_back("cosmics");
0226       sourceNames_.push_back("beamhalo_plus");
0227       sourceNames_.push_back("beamhalo_minus");
0228
0229       for (size_t makeIdx = 0; makeIdx < maxNbSources; makeIdx++) {
0230         inputConfigs_.push_back(maybeConfigPileUp(pset, sourceNames_[makeIdx], minBunch_, maxBunch_, playback_));
0231       }
0232
0233       if (pset.exists("readDB"))
0234         configFromDB_ = pset.getParameter<bool>("readDB");
0235     }
0236   }  // namespace MixingCache
0237
0238   std::unique_ptr<MixingCache::Config> BMixingModule::initializeGlobalCache(edm::ParameterSet const& pset) {
0239     return std::make_unique<MixingCache::Config>(pset, maxNbSources_);
0240   }
0241
0242   // update method call at begin run/lumi to reload the mixing configuration
0243   void BMixingModule::beginLuminosityBlock(edm::LuminosityBlock const& lumi, edm::EventSetup const& setup) {
0244     update(setup);
0245     for (size_t endIdx = 0; endIdx < maxNbSources_; ++endIdx) {
0246       if (inputSources_[endIdx])
0247         inputSources_[endIdx]->beginLuminosityBlock(lumi, setup);
0248     }
0249   }
0250
0251   void BMixingModule::beginRun(edm::Run const& run, edm::EventSetup const& setup) {
0252     for (size_t endIdx = 0; endIdx < maxNbSources_; ++endIdx) {
0253       if (inputSources_[endIdx])
0254         inputSources_[endIdx]->beginRun(run, setup);
0255     }
0256   }
0257
0258   void BMixingModule::endLuminosityBlock(edm::LuminosityBlock const& lumi, edm::EventSetup const& setup) {
0259     for (size_t endIdx = 0; endIdx < maxNbSources_; ++endIdx) {
0260       if (inputSources_[endIdx])
0261         inputSources_[endIdx]->endLuminosityBlock(lumi, setup);
0262     }
0263   }
0264
0265   void BMixingModule::endRun(edm::Run const& run, edm::EventSetup const& setup) {
0266     for (size_t endIdx = 0; endIdx < maxNbSources_; ++endIdx) {
0267       if (inputSources_[endIdx])
0268         inputSources_[endIdx]->endRun(run, setup);
0269     }
0270   }
0271
0272   void BMixingModule::update(const edm::EventSetup& setup) {
0273     if (readDB_ && parameterWatcher_.check(setup)) {
0274       for (size_t makeIdx = 0; makeIdx < maxNbSources_; makeIdx++) {
0275         if (inputSources_[makeIdx])
0276           inputSources_[makeIdx]->reload(setup);
0277       }
0278       reload(setup);
0279     }
0280   }
0281
0282   // Functions that get called by framework every event
0283   void BMixingModule::produce(edm::Event& e, const edm::EventSetup& setup) {
0284     // Check if the signal is present in the root file
0285     // for all the objects we want to mix
0286     checkSignal(e);
0287
0288     // Create EDProduct
0289     createnewEDProduct();
0290
0291     initializeEvent(e, setup);
0292
0293     // Add signals
0294     if (!mixProdStep1_) {
0295       addSignals(e, setup);
0296     }
0297
0298     doPileUp(e, setup);
0299
0300     // Includes putting digi products into the edm::Event.
0301     finalizeEvent(e, setup);
0302
0303     // Put output into event (here only playback info)
0304     put(e, setup);
0305   }
0306
0307   void BMixingModule::setupPileUpEvent(const edm::EventSetup& setup) {
0308     for (size_t dropIdx = 0; dropIdx < maxNbSources_; ++dropIdx) {
0309       if (inputSources_[dropIdx])
0310         inputSources_[dropIdx]->setupPileUpEvent(setup);
0311     }
0312   }
0313
0314   void BMixingModule::dropUnwantedBranches(std::vector<std::string> const& wantedBranches) {
0315     for (size_t dropIdx = 0; dropIdx < maxNbSources_; ++dropIdx) {
0316       if (inputSources_[dropIdx])
0317         inputSources_[dropIdx]->dropUnwantedBranches(wantedBranches);
0318     }
0319   }
0320
0321   void BMixingModule::beginStream(edm::StreamID streamID) {
0322     for (size_t endIdx = 0; endIdx < maxNbSources_; ++endIdx) {
0323       if (inputSources_[endIdx])
0324         inputSources_[endIdx]->beginStream(streamID);
0325     }
0326   }
0327
0328   void BMixingModule::endStream() {
0329     ExceptionCollector exceptionCollector(
0330         "Multiple exceptions were thrown while executing endStream and endJob for mixing modules. "
0331         "An exception message follows for each.\n");
0332
0333     for (size_t endIdx = 0; endIdx < maxNbSources_; ++endIdx) {
0334       if (inputSources_[endIdx])
0335         inputSources_[endIdx]->endStream(exceptionCollector);
0336     }
0337     if (exceptionCollector.hasThrown()) {
0338       exceptionCollector.rethrow();
0339     }
0340   }
0341
0342   void BMixingModule::createnewEDProduct() {
0343     edm::LogWarning("MixingModule") << "BMixingModule::createnewEDProduct must be overwritten!";
0344   }
0345
0346   void BMixingModule::checkSignal(const edm::Event& e) {
0347     edm::LogWarning("MixingModule") << "BMixingModule::checkSignal must be overwritten!";
0348   }
0349
0350   void BMixingModule::setBcrOffset() {
0351     edm::LogWarning("MixingModule") << "BMixingModule::setBcrOffset must be overwritten!";
0352   }
0353
0354   void BMixingModule::setSourceOffset(const unsigned int s) {
0355     edm::LogWarning("MixingModule") << "BMixingModule::setSourceOffset must be overwritten!";
0356   }
0357
0358   void BMixingModule::doPileUp(edm::Event& e, const edm::EventSetup& c) {
0359     edm::LogWarning("MixingModule") << "BMixingModule::doPileUp must be overwritten!";
0360   }
0361
0362 }  // namespace edm