Luminosity producer using files from brilcalc

RecoLuminosity/LumiProducer/README.md

@@ -0,0 +1,19 @@
+# RecoLuminosity/LumiProducer
+
+This package contains utilities for producing luminosity information in CMSSW. Most of the plugins here are obsolete and only for Run 1; they are kept here for backwards compatibility but are not actively maintained.
+
+For Run 2, the plugin `LumiProducerFromBrilcalc` is available. This allows you to add luminosity information to your CMSSW job by reading in a CSV output file produced by the `brilcalc` utility. In order to use it, add the following to your `cfg` file:
+
+```
+process.LumiInfo = cms.EDProducer('LumiProducerFromBrilcalc',
+                                  lumiFile = cms.string("./myLumiFile.csv"),
+                                  throwIfNotFound = cms.bool(False),
+                                  doBunchByBunch = cms.bool(False))
+```
+where `lumiFile` is the output file created by `brilcalc`, `throwIfNotFound` will determine the behavior if the csv file does not contain information for an event in your input file (if `True`, an exception will be thrown; if `False`, the luminosity will just be taken to be zero for that event), and `doBunchByBunch` should remain `False`, since bunch-by-bunch luminosity is not currently supported.
+
+Add `LumiInfo` to your path, and then you can access the `LumiInfo` object produced under the input tag `("LumiInfo", "brilcalc")`.
+
+For more information on the proper way to use `brilcalc` to produce a luminosity csv file, please consult the [LumiPOG twiki](https://twiki.cern.ch/twiki/bin/view/CMS/TWikiLUM#LumiCMSSW).
+
+

RecoLuminosity/LumiProducer/plugins/BuildFile.xml

@@ -106,3 +106,9 @@
   <use name="FWCore/PluginManager"/>
   <flags   EDM_PLUGIN="1"/>
 </library>
+<library file="LumiProducerFromBrilcalc.cc" name="LumiProducerFromBrilcalc">
+  <use name="FWCore/PluginManager"/>
+  <use name="FWCore/ParameterSet"/>
+  <use name="DataFormats/Luminosity"/>
+  <flags EDM_PLUGIN="1"/>
+</library>

RecoLuminosity/LumiProducer/plugins/LumiProducerFromBrilcalc.cc

@@ -0,0 +1,174 @@
+// -*- C++ -*-
+//
+// Package:    RecoLuminosity/LumiProducer
+// Class:      LumiProducerFromBrilcalc
+//
+/**\class LumiProducerFromBrilcalc LumiProducerFromBrilcalc.cc RecoLuminosity/LumiProducer/plugins/LumiProducerFromBrilcalc.cc
+
+   Description: Takes a csv file with luminosity information produced by brilcalc and reads it into the event.
+
+   Implementation:
+   [Notes on implementation]
+*/
+//
+// Original Author:  Paul Lujan
+//         Created:  Fri, 28 Feb 2020 16:32:38 GMT
+//
+//
+
+// system include files
+#include <memory>
+#include <sstream>
+#include <fstream>
+
+// user include files
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/global/EDProducer.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/StreamID.h"
+
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+#include "DataFormats/Luminosity/interface/LumiInfo.h"
+
+//
+// class declaration
+//
+
+class LumiProducerFromBrilcalc : public edm::global::EDProducer<> {
+public:
+  explicit LumiProducerFromBrilcalc(const edm::ParameterSet&);
+  ~LumiProducerFromBrilcalc() = default;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  virtual void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
+
+  // ----------member data ---------------------------
+  const std::string lumiFile_;
+  const bool throwIfNotFound_;
+  const bool doBunchByBunch_;
+  std::map<std::pair<int, int>, std::pair<float, float> > lumiData_;
+};
+
+//
+// constructor
+//
+
+LumiProducerFromBrilcalc::LumiProducerFromBrilcalc(const edm::ParameterSet& iConfig)
+    : lumiFile_(iConfig.getParameter<std::string>("lumiFile")),
+      throwIfNotFound_(iConfig.getParameter<bool>("throwIfNotFound")),
+      doBunchByBunch_(iConfig.getParameter<bool>("doBunchByBunch")) {
+  //register your products
+  produces<LumiInfo>("brilcalc");
+
+  //now do what ever other initialization is needed
+  if (doBunchByBunch_) {
+    throw cms::Exception("LumiProducerFromBrilcalc")
+        << "Sorry, bunch-by-bunch luminosity is not yet supported. Please bug your friendly lumi expert!";
+  }
+
+  // Read luminosity data and store it in lumiData_.
+  edm::LogInfo("LumiProducerFromBrilcalc") << "Reading luminosity data from " << lumiFile_ << "...one moment...";
+  std::ifstream lumiFile(lumiFile_);
+  if (!lumiFile.is_open()) {
+    throw cms::Exception("LumiProducerFromBrilcalc") << "Failed to open input luminosity file " << lumiFile_;
+  }
+
+  int nLS = 0;
+  std::string line;
+  while (1) {
+    std::getline(lumiFile, line);
+    if (lumiFile.eof() || lumiFile.fail())
+      break;
+    if (line.empty())
+      continue;  // skip blank lines
+    if (line.at(0) == '#')
+      continue;  // skip comment lines
+
+    // Break into fields. These should be, in order: run:fill, brills:cmsls, time, beam status, beam energy,
+    // delivered lumi, recorded lumi, average pileup, source.
+    std::stringstream ss(line);
+    std::string field;
+    std::vector<std::string> fields;
+
+    while (std::getline(ss, field, ','))
+      fields.push_back(field);
+
+    if (fields.size() != 9) {
+      edm::LogWarning("LumiProducerFromBrilcalc") << "Malformed line in csv file: " << line;
+      continue;
+    }
+
+    // Extract the run number from fields[0] and the lumisection number from fields[1]. Fortunately since the
+    // thing we want is before the colon, we don't have to split them again.
+    int run, ls;
+    std::stringstream runfill(fields[0]);
+    runfill >> run;
+    std::stringstream lsls(fields[1]);
+    lsls >> ls;
+
+    // Extract the delivered and recorded lumi from fields[5] and fields[6].
+    float lumiDel, lumiRec, dtFrac;
+    std::stringstream lumiDelString(fields[5]);
+    lumiDelString >> lumiDel;
+    std::stringstream lumiRecString(fields[6]);
+    lumiRecString >> lumiRec;
+
+    // Calculate the deadtime fraction
+    dtFrac = 1.0 - lumiRec / lumiDel;
+
+    // Finally, store it all.
+    lumiData_[std::make_pair(run, ls)] = std::make_pair(lumiDel, dtFrac);
+    nLS++;
+  }
+  edm::LogInfo("LumiProducerFromBrilcalc") << "Read " << nLS << " lumisections from " << lumiFile_;
+  lumiFile.close();
+}
+
+//
+// member functions
+//
+
+// ------------ method called to produce the data  ------------
+void LumiProducerFromBrilcalc::produce(edm::StreamID iStreamID,
+                                       edm::Event& iEvent,
+                                       const edm::EventSetup& iSetup) const {
+  std::vector<float> bxlumi(3564, 0);
+  std::pair<int, int> runls = std::make_pair(iEvent.run(), iEvent.luminosityBlock());
+  if (lumiData_.count(runls) == 1) {
+    // if we have data for this run/LS, put it in the event
+    LogDebug("LumiProducerFromBrilcalc") << "Filling for run " << runls.first << " ls " << runls.second
+                                         << " with delivered " << lumiData_.at(runls).first << " dt "
+                                         << lumiData_.at(runls).second;
+    iEvent.put(std::make_unique<LumiInfo>(lumiData_.at(runls).second, bxlumi, lumiData_.at(runls).first), "brilcalc");
+  } else {
+    if (throwIfNotFound_) {
+      throw cms::Exception("LumiProducerFromBrilcalc")
+          << "Failed to find luminosity for run " << runls.first << " LS " << runls.second;
+    } else {
+      // just put in zeroes
+      edm::LogWarning("LumiProducerFromBrilcalc")
+          << "Failed to find luminosity for run " << runls.first << " ls " << runls.second;
+      iEvent.put(std::make_unique<LumiInfo>(0, bxlumi, 0), "brilcalc");
+    }
+  }
+}
+
+// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
+void LumiProducerFromBrilcalc::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  // Allowed parameters
+  edm::ParameterSetDescription desc;
+  desc.add<std::string>("lumiFile");
+  desc.add<bool>("throwIfNotFound", false);
+  desc.add<bool>("doBunchByBunch", false);
+  descriptions.addDefault(desc);
+}
+
+//define this as a plug-in
+DEFINE_FWK_MODULE(LumiProducerFromBrilcalc);

RecoLuminosity/LumiProducer/test/BuildFile.xml

@@ -105,3 +105,6 @@
   <flags EDM_PLUGIN="1"/>
 </library>
 
+<library file="TestLumiProducerFromBrilcalc.cc" name="TestLumiProducerFromBrilcalc">
+  <flags EDM_PLUGIN="1"/>
+</library>

RecoLuminosity/LumiProducer/test/TestLumiProducerFromBrilcalc.cc

@@ -0,0 +1,109 @@
+// -*- C++ -*-
+//
+// Package:    RecoLuminosity/TestLumiProducerFromBrilcalc
+// Class:      TestLumiProducerFromBrilcalc
+//
+/**\class TestLumiProducerFromBrilcalc TestLumiProducerFromBrilcalc.cc RecoLuminosity/LumiProducer/test/TestLumiProducerFromBrilcalc.cc
+
+   Description: A simple analyzer class to test the functionality of TestLumiProducerFromBrilcalc.
+
+   Implementation:
+   Get the luminosity and prints it out.
+*/
+//
+// Original Author:  Paul Lujan
+//         Created:  Fri, 20 Mar 2020 09:32:27 GMT
+//
+//
+
+// system include files
+#include <memory>
+#include <iostream>
+
+// user include files
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/one/EDAnalyzer.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+
+#include "DataFormats/Luminosity/interface/LumiInfo.h"
+
+//
+// class declaration
+//
+
+// If the analyzer does not use TFileService, please remove
+// the template argument to the base class so the class inherits
+// from edm::one::EDAnalyzer<>
+// This will improve performance in multithreaded jobs.
+
+class TestLumiProducerFromBrilcalc : public edm::one::EDAnalyzer<> {
+public:
+  explicit TestLumiProducerFromBrilcalc(const edm::ParameterSet&);
+  ~TestLumiProducerFromBrilcalc();
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  virtual void beginJob() override;
+  virtual void analyze(const edm::Event&, const edm::EventSetup&) override;
+  virtual void endJob() override;
+
+  // ----------member data ---------------------------
+  edm::InputTag inputTag_;
+  edm::EDGetTokenT<LumiInfo> lumiToken_;
+};
+
+//
+// constants, enums and typedefs
+//
+
+//
+// static data member definitions
+//
+
+//
+// constructors and destructor
+//
+TestLumiProducerFromBrilcalc::TestLumiProducerFromBrilcalc(const edm::ParameterSet& iConfig)
+    : inputTag_(iConfig.getUntrackedParameter<edm::InputTag>("inputTag")), lumiToken_(consumes<LumiInfo>(inputTag_)) {}
+
+TestLumiProducerFromBrilcalc::~TestLumiProducerFromBrilcalc() {
+  // do anything here that needs to be done at destruction time
+  // (e.g. close files, deallocate resources etc.)
+}
+
+//
+// member functions
+//
+
+// ------------ method called for each event  ------------
+void TestLumiProducerFromBrilcalc::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  using namespace edm;
+
+  const LumiInfo& lumi = iEvent.get(lumiToken_);
+
+  std::cout << "Luminosity for " << iEvent.run() << " LS " << iEvent.luminosityBlock() << " is "
+            << lumi.getTotalInstLumi() << std::endl;
+}
+
+// ------------ method called once each job just before starting event loop  ------------
+void TestLumiProducerFromBrilcalc::beginJob() {}
+
+// ------------ method called once each job just after ending the event loop  ------------
+void TestLumiProducerFromBrilcalc::endJob() {}
+
+// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
+void TestLumiProducerFromBrilcalc::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  // Allowed parameters
+  edm::ParameterSetDescription desc;
+  desc.addUntracked<edm::InputTag>("inputTag");
+  descriptions.addDefault(desc);
+}
+
+//define this as a plug-in
+DEFINE_FWK_MODULE(TestLumiProducerFromBrilcalc);

RecoLuminosity/LumiProducer/test/TestLumiProducerFromBrilcalc_cfg.py

@@ -0,0 +1,30 @@
+import FWCore.ParameterSet.Config as cms
+
+process = cms.Process("TestLumiProducerFromBrilcalc")
+
+process.load("FWCore.MessageLogger.MessageLogger_cfi")
+process.MessageLogger = cms.Service("MessageLogger",
+                                    destinations   = cms.untracked.vstring("cout"),
+                                    categories     = cms.untracked.vstring("LumiProducerFromBrilcalc"),
+                                    debugModules   = cms.untracked.vstring("LumiInfo"),
+                                    cout           = cms.untracked.PSet(
+                                        threshold  = cms.untracked.string('DEBUG')
+                                    )
+)
+
+# just use a random relval which has meaningless run/LS numbers, and then a corresponding test file
+process.source = cms.Source("PoolSource",
+                            fileNames = cms.untracked.vstring("/store/relval/CMSSW_10_2_7/RelValMinBias_13/GEN-SIM/102X_mc2017_realistic_v4_AC_v01-v1/20000/BB654FAE-5375-164F-BBFE-B330713759C6.root")
+)
+process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(10))
+
+process.LumiInfo = cms.EDProducer('LumiProducerFromBrilcalc',
+                                  lumiFile = cms.string("./testLumiFile.csv"),
+                                  throwIfNotFound = cms.bool(False),
+                                  doBunchByBunch = cms.bool(False))
+
+process.test = cms.EDAnalyzer('TestLumiProducerFromBrilcalc',
+                              inputTag = cms.untracked.InputTag("LumiInfo", "brilcalc"))
+
+process.p = cms.Path(process.LumiInfo*process.test)
+

RecoLuminosity/LumiProducer/test/testLumiFile.csv

@@ -0,0 +1,4 @@
+# This is a dummy file for TestLumiProducerFromBrilcalc.
+#run:fill,ls,time,beamstatus,E(GeV),delivered(hz/ub),recorded(hz/ub),avgpu,source
+1:1,1:1,04/17/18 10:54:40,STABLE BEAMS,6500,16.507442576,16.080153561,58.7,HFOC
+1:1,2:2,04/17/18 10:55:03,STABLE BEAMS,6500,16.483134426,1.668723696,58.6,HFOC