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