Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 10:43:49 

0001 /**_________________________________________________________________
0002 class:   RawPCCProducer.cc
0003 
0004 description: Creates a LumiInfo object that will contain the luminosity per bunch crossing,
0005 along with the total luminosity and the statistical error.
0006 
0007 authors:Sam Higginbotham (shigginb@cern.ch), Chris Palmer (capalmer@cern.ch), Jose Benitez (jose.benitez@cern.ch)
0008 
0009 ________________________________________________________________**/
0010 #include <string>
0011 #include <iostream>
0012 #include <fstream>
0013 #include <vector>
0014 #include <mutex>
0015 #include <cmath>
0016 #include "DataFormats/Luminosity/interface/PixelClusterCounts.h"
0017 #include "DataFormats/Luminosity/interface/LumiInfo.h"
0018 #include "DataFormats/Luminosity/interface/LumiConstants.h"
0019 #include "CondFormats/Luminosity/interface/LumiCorrections.h"
0020 #include "CondFormats/DataRecord/interface/LumiCorrectionsRcd.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "FWCore/Framework/interface/MakerMacros.h"
0023 #include "FWCore/Framework/interface/ConsumesCollector.h"
0024 #include "FWCore/Framework/interface/ESHandle.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/global/EDProducer.h"
0028 #include "FWCore/Framework/interface/Event.h"
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/Utilities/interface/EDGetToken.h"
0031 #include "FWCore/ServiceRegistry/interface/Service.h"
0032 #include "FWCore/Framework/interface/LuminosityBlock.h"
0033 
0034 class RawPCCProducer : public edm::global::EDProducer<edm::EndLuminosityBlockProducer> {
0035 public:
0036   explicit RawPCCProducer(const edm::ParameterSet&);
0037   ~RawPCCProducer() override;
0038 
0039 private:
0040   void globalEndLuminosityBlockProduce(edm::LuminosityBlock& lumiSeg, const edm::EventSetup& iSetup) const final;
0041   void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const final;
0042 
0043   //input object labels
0044   edm::EDGetTokenT<reco::PixelClusterCounts> pccToken_;
0045   //background corrections from DB
0046   const edm::ESGetToken<LumiCorrections, LumiCorrectionsRcd> lumiCorrectionsToken_;
0047   //The list of modules to skip in the lumi calc.
0048   const std::vector<int> modVeto_;
0049   //background corrections
0050   const bool applyCorr_;
0051   //Output average values
0052   const std::string takeAverageValue_;
0053 
0054   //output object labels
0055   const edm::EDPutTokenT<LumiInfo> putToken_;
0056 
0057   //produce csv lumi file
0058   const bool saveCSVFile_;
0059   const std::string csvOutLabel_;
0060   mutable std::mutex fileLock_;
0061 };
0062 
0063 //--------------------------------------------------------------------------------------------------
0064 RawPCCProducer::RawPCCProducer(const edm::ParameterSet& iConfig)
0065     : pccToken_(consumes<reco::PixelClusterCounts, edm::InLumi>(edm::InputTag(
0066           iConfig.getParameter<edm::ParameterSet>("RawPCCProducerParameters").getParameter<std::string>("inputPccLabel"),
0067           iConfig.getParameter<edm::ParameterSet>("RawPCCProducerParameters").getParameter<std::string>("ProdInst")))),
0068       lumiCorrectionsToken_(esConsumes<edm::Transition::EndLuminosityBlock>()),
0069       modVeto_(
0070           iConfig.getParameter<edm::ParameterSet>("RawPCCProducerParameters").getParameter<std::vector<int>>("modVeto")),
0071       applyCorr_(iConfig.getParameter<edm::ParameterSet>("RawPCCProducerParameters")
0072                      .getUntrackedParameter<bool>("ApplyCorrections", false)),
0073       takeAverageValue_(iConfig.getParameter<edm::ParameterSet>("RawPCCProducerParameters")
0074                             .getUntrackedParameter<std::string>("OutputValue", std::string("Average"))),
0075       putToken_(produces<LumiInfo, edm::Transition::EndLuminosityBlock>(
0076           iConfig.getParameter<edm::ParameterSet>("RawPCCProducerParameters")
0077               .getUntrackedParameter<std::string>("outputProductName", "alcaLumi"))),
0078       saveCSVFile_(iConfig.getParameter<edm::ParameterSet>("RawPCCProducerParameters")
0079                        .getUntrackedParameter<bool>("saveCSVFile", false)),
0080       csvOutLabel_(iConfig.getParameter<edm::ParameterSet>("RawPCCProducerParameters")
0081                        .getUntrackedParameter<std::string>("label", std::string("rawPCC.csv"))) {}
0082 
0083 //--------------------------------------------------------------------------------------------------
0084 RawPCCProducer::~RawPCCProducer() {}
0085 
0086 //--------------------------------------------------------------------------------------------------
0087 void RawPCCProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {}
0088 
0089 //--------------------------------------------------------------------------------------------------
0090 void RawPCCProducer::globalEndLuminosityBlockProduce(edm::LuminosityBlock& lumiSeg,
0091                                                      const edm::EventSetup& iSetup) const {
0092   //The total raw luminosity from the pixel clusters - not scaled
0093   float totalLumi = 0.0;
0094   //the statistical error on the lumi - large num ie sqrt(N)
0095   float statErrOnLumi = 0.0;
0096 
0097   //new vector containing clusters per bxid
0098   std::vector<int> clustersPerBXOutput(LumiConstants::numBX, 0);
0099   //new vector containing clusters per bxid with afterglow corrections
0100   std::vector<float> corrClustersPerBXOutput(LumiConstants::numBX, 0);
0101 
0102   //////////////////////////////////
0103   /// read input , clusters per module per bx
0104   /////////////////////////////////
0105   const edm::Handle<reco::PixelClusterCounts> pccHandle = lumiSeg.getHandle(pccToken_);
0106   const reco::PixelClusterCounts& inputPcc = *(pccHandle.product());
0107   //vector with Module IDs 1-1 map to bunch x-ing in clusers
0108   auto modID = inputPcc.readModID();
0109   //vector with total events at each bxid.
0110   auto events = inputPcc.readEvents();
0111   //cluster counts per module per bx
0112   auto clustersPerBXInput = inputPcc.readCounts();
0113 
0114   ////////////////////////////
0115   ///Apply the module veto
0116   ///////////////////////////
0117   std::vector<int> goodMods;
0118   for (unsigned int i = 0; i < modID.size(); i++) {
0119     if (std::find(modVeto_.begin(), modVeto_.end(), modID.at(i)) == modVeto_.end()) {
0120       goodMods.push_back(i);
0121     }
0122   }
0123   for (int bx = 0; bx < int(LumiConstants::numBX); bx++) {
0124     for (unsigned int i = 0; i < goodMods.size(); i++) {
0125       clustersPerBXOutput.at(bx) += clustersPerBXInput.at(goodMods.at(i) * int(LumiConstants::numBX) + bx);
0126     }
0127   }
0128 
0129   //////////////////////////////
0130   //// Apply afterglow corrections
0131   //////////////////////////////
0132   std::vector<float> correctionScaleFactors;
0133   if (applyCorr_) {
0134     const auto pccCorrections = &iSetup.getData(lumiCorrectionsToken_);
0135     correctionScaleFactors = pccCorrections->getCorrectionsBX();
0136   } else {
0137     correctionScaleFactors.resize(LumiConstants::numBX, 1.0);
0138   }
0139 
0140   for (unsigned int i = 0; i < clustersPerBXOutput.size(); i++) {
0141     if (events.at(i) != 0) {
0142       corrClustersPerBXOutput[i] = clustersPerBXOutput[i] * correctionScaleFactors[i];
0143     } else {
0144       corrClustersPerBXOutput[i] = 0.0;
0145     }
0146     totalLumi += corrClustersPerBXOutput[i];
0147     statErrOnLumi += float(events[i]);
0148   }
0149 
0150   std::vector<float> errorPerBX;  //Stat error (or number of events)
0151   errorPerBX.assign(events.begin(), events.end());
0152 
0153   //////////////////////////////////
0154   /// Compute average number of clusters per event
0155   ////////////////////////////////
0156   if (takeAverageValue_ == "Average") {
0157     unsigned int NActiveBX = 0;
0158     for (int bx = 0; bx < int(LumiConstants::numBX); bx++) {
0159       if (events[bx] > 0) {
0160         NActiveBX++;
0161         corrClustersPerBXOutput[bx] /= float(events[bx]);
0162         errorPerBX[bx] = 1 / sqrt(float(events[bx]));
0163       }
0164     }
0165     if (statErrOnLumi != 0) {
0166       totalLumi = totalLumi / statErrOnLumi * float(NActiveBX);
0167       statErrOnLumi = 1 / sqrt(statErrOnLumi) * totalLumi;
0168     }
0169   }
0170 
0171   ///////////////////////////////////////////////////////
0172   ///Lumi saved in the LuminosityBlocks
0173   LumiInfo outputLumiInfo;
0174   outputLumiInfo.setTotalInstLumi(totalLumi);
0175   outputLumiInfo.setTotalInstStatError(statErrOnLumi);
0176   outputLumiInfo.setErrorLumiAllBX(errorPerBX);
0177   outputLumiInfo.setInstLumiAllBX(corrClustersPerBXOutput);
0178   lumiSeg.emplace(putToken_, std::move(outputLumiInfo));
0179 
0180   //Lumi saved in the csv file
0181   if (saveCSVFile_) {
0182     std::lock_guard<std::mutex> lock(fileLock_);
0183     std::ofstream csfile(csvOutLabel_, std::ios_base::app);
0184     csfile << std::to_string(lumiSeg.run()) << ",";
0185     csfile << std::to_string(lumiSeg.luminosityBlock()) << ",";
0186     csfile << std::to_string(totalLumi);
0187 
0188     for (unsigned int bx = 0; bx < LumiConstants::numBX; bx++)
0189       csfile << "," << std::to_string(corrClustersPerBXOutput[bx]);
0190     csfile << std::endl;
0191 
0192     csfile.close();
0193   }
0194 }
0195 
0196 DEFINE_FWK_MODULE(RawPCCProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!