PixelLumi/plugins/PixelLumiDQM.cc

0001 // -*- C++ -*-
0002
0003 // Package:    PixelLumiDQM
0004 // Class:      PixelLumiDQM
0005
0006 // Author: Amita Raval
0007 // Based on Jeroen Hegeman's code for Pixel Cluster Count Luminosity
0008
0009 #include "PixelLumiDQM.h"
0010
0011 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0012 #include "DQMServices/Core/interface/DQMStore.h"
0013 #include "DataFormats/Common/interface/Handle.h"
0014 #include "DataFormats/DetId/interface/DetId.h"
0015 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0016 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0017 #include "DataFormats/GeometryVector/interface/Pi.h"
0018 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0019 #include "DataFormats/SiPixelDetId/interface/PixelBarrelNameUpgrade.h"
0020 #include "DataFormats/SiPixelDetId/interface/PixelEndcapNameUpgrade.h"
0021 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0022 #include "FWCore/Framework/interface/ESHandle.h"
0023 #include "FWCore/Framework/interface/Event.h"
0024 #include "FWCore/Framework/interface/EventSetup.h"
0025 #include "FWCore/Framework/interface/LuminosityBlock.h"
0026 #include "FWCore/Framework/interface/MakerMacros.h"
0027 #include "FWCore/Framework/interface/Run.h"
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/ServiceRegistry/interface/Service.h"
0031 #include "FWCore/Utilities/interface/EDMException.h"
0032 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0033 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0034 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0035 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0036
0037 #include <ctime>
0038 #include <fstream>
0039 #include <map>
0040 #include <string>
0041 #include <sys/time.h>
0042 #include <vector>
0043
0044 const unsigned int PixelLumiDQM::lastBunchCrossing;
0045
0046 // Constructors and destructor.
0047 PixelLumiDQM::PixelLumiDQM(const edm::ParameterSet &iConfig)
0048     : fPixelClusterLabel(consumes<edmNew::DetSetVector<SiPixelCluster>>(
0049           iConfig.getUntrackedParameter<edm::InputTag>("pixelClusterLabel", edm::InputTag("siPixelClusters")))),
0050       tkGeomToken_(esConsumes()),
0051       fIncludePixelClusterInfo(iConfig.getUntrackedParameter<bool>("includePixelClusterInfo", true)),
0052       fIncludePixelQualCheckHistos(iConfig.getUntrackedParameter<bool>("includePixelQualCheckHistos", true)),
0053       fResetIntervalInLumiSections(iConfig.getUntrackedParameter<int>("resetEveryNLumiSections", 1)),
0054       fDeadModules(iConfig.getUntrackedParameter<std::vector<uint32_t>>("deadModules", std::vector<uint32_t>())),
0055       fMinPixelsPerCluster(iConfig.getUntrackedParameter<int>("minNumPixelsPerCluster", 0)),
0056       fMinClusterCharge(iConfig.getUntrackedParameter<double>("minChargePerCluster", 0)),
0057       bunchTriggerMask(lastBunchCrossing + 1, false),
0058       filledAndUnmaskedBunches(0),
0059       useInnerBarrelLayer(iConfig.getUntrackedParameter<bool>("useInnerBarrelLayer", false)),
0060       fLogFileName_(iConfig.getUntrackedParameter<std::string>("logFileName", "/tmp/pixel_lumi.txt")) {
0061   edm::LogInfo("Configuration") << "PixelLumiDQM looking for pixel clusters in '"
0062                                 << iConfig.getUntrackedParameter<edm::InputTag>("pixelClusterLabel",
0063                                                                                 edm::InputTag("siPixelClusters"))
0064                                 << "'";
0065   edm::LogInfo("Configuration") << "PixelLumiDQM storing pixel cluster info? " << fIncludePixelClusterInfo;
0066   edm::LogInfo("Configuration") << "PixelLumiDQM storing pixel cluster quality check histograms? "
0067                                 << fIncludePixelQualCheckHistos;
0068
0069   if (fDeadModules.empty()) {
0070     edm::LogInfo("Configuration") << "No pixel modules specified to be ignored";
0071   } else {
0072     edm::LogInfo("Configuration") << fDeadModules.size() << " pixel modules specified to be ignored:";
0073     for (std::vector<uint32_t>::const_iterator it = fDeadModules.begin(); it != fDeadModules.end(); ++it) {
0074       edm::LogInfo("Configuration") << "  " << *it;
0075     }
0076   }
0077   edm::LogInfo("Configuration") << "Ignoring pixel clusters with less than " << fMinPixelsPerCluster << " pixels";
0078   edm::LogInfo("Configuration") << "Ignoring pixel clusters with charge less than " << fMinClusterCharge;
0079 }
0080
0081 PixelLumiDQM::~PixelLumiDQM() {}
0082
0083 void PixelLumiDQM::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0084   edm::ParameterSetDescription desc;
0085   desc.setUnknown();
0086   descriptions.addDefault(desc);
0087 }
0088
0089 void PixelLumiDQM::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0090   // Collect all bookkeeping information.
0091   fRunNo = iEvent.id().run();
0092   fEvtNo = iEvent.id().event();
0093   fLSNo = iEvent.getLuminosityBlock().luminosityBlock();
0094   fBXNo = iEvent.bunchCrossing();
0095   fTimestamp = iEvent.time().unixTime();
0096   fHistBunchCrossings->Fill(float(fBXNo));
0097   fHistBunchCrossingsLastLumi->Fill(float(fBXNo));
0098   // This serves as event counter to compute luminosity from cluster counts.
0099   std::map<int, PixelClusterCount>::iterator it = fNumPixelClusters.find(fBXNo);
0100   if (it == fNumPixelClusters.end())
0101     fNumPixelClusters[fBXNo] = PixelClusterCount();
0102   // Find tracker geometry.
0103   const TrackerGeometry *trackerGeo = &iSetup.getData(tkGeomToken_);
0104   if (fIncludePixelClusterInfo) {
0105     // Find pixel clusters.
0106     edm::Handle<edmNew::DetSetVector<SiPixelCluster>> pixelClusters;
0107     iEvent.getByToken(fPixelClusterLabel, pixelClusters);
0108
0109     // Loop over entire tracker geometry.
0110     for (TrackerGeometry::DetContainer::const_iterator i = trackerGeo->dets().begin(); i != trackerGeo->dets().end();
0111          ++i) {
0112       // See if this is a pixel unit(?).
0113
0114       if (GeomDetEnumerators::isTrackerPixel((*i)->subDetector())) {
0115         DetId detId = (*i)->geographicalId();
0116         // Find all clusters on this detector module.
0117         edmNew::DetSetVector<SiPixelCluster>::const_iterator iSearch = pixelClusters->find(detId);
0118         if (iSearch != pixelClusters->end()) {
0119           // Count the number of clusters with at least a minimum
0120           // number of pixels per cluster and at least a minimum charge.
0121           size_t numClusters = 0;
0122           for (edmNew::DetSet<SiPixelCluster>::const_iterator itClus = iSearch->begin(); itClus != iSearch->end();
0123                ++itClus) {
0124             if ((itClus->size() >= fMinPixelsPerCluster) && (itClus->charge() >= fMinClusterCharge)) {
0125               ++numClusters;
0126             }
0127           }
0128           // DEBUG DEBUG DEBUG
0129           assert(numClusters <= iSearch->size());
0130           // DEBUG DEBUG DEBUG end
0131
0132           // Add up the cluster count based on the position of this detector
0133           // element.
0134           if (detId.subdetId() == PixelSubdetector::PixelBarrel) {
0135             PixelBarrelNameUpgrade detName = PixelBarrelNameUpgrade(detId);
0136             int layer = detName.layerName() - kOffsetLayers;
0137             fNumPixelClusters[fBXNo].numB.at(layer) += numClusters;
0138             fNumPixelClusters[fBXNo].dnumB.at(layer) += sqrt(numClusters);
0139           } else {
0140             // DEBUG DEBUG DEBUG
0141             assert(detId.subdetId() == PixelSubdetector::PixelEndcap);
0142             // DEBUG DEBUG DEBUG end
0143
0144             PixelEndcapNameUpgrade detName = PixelEndcapNameUpgrade(detId);
0145             PixelEndcapNameUpgrade::HalfCylinder halfCylinder = detName.halfCylinder();
0146             int disk = detName.diskName() - kOffsetDisks;
0147             switch (halfCylinder) {
0148               case PixelEndcapNameUpgrade::mO:
0149               case PixelEndcapNameUpgrade::mI:
0150                 fNumPixelClusters[fBXNo].numFM.at(disk) += numClusters;
0151                 fNumPixelClusters[fBXNo].dnumFM.at(disk) += sqrt(numClusters);
0152                 break;
0153               case PixelEndcapNameUpgrade::pO:
0154               case PixelEndcapNameUpgrade::pI:
0155                 fNumPixelClusters[fBXNo].numFP.at(disk) += numClusters;
0156                 fNumPixelClusters[fBXNo].dnumFP.at(disk) += sqrt(numClusters);
0157                 break;
0158               default:
0159                 assert(false);
0160                 break;
0161             }
0162           }
0163         }
0164       }
0165     }
0166   }
0167   // ----------
0168
0169   // Fill some pixel cluster quality check histograms if requested.
0170   if (fIncludePixelQualCheckHistos) {
0171     // Find pixel clusters.
0172     edm::Handle<edmNew::DetSetVector<SiPixelCluster>> pixelClusters;
0173     iEvent.getByToken(fPixelClusterLabel, pixelClusters);
0174
0175     bool filterDeadModules = (!fDeadModules.empty());
0176     std::vector<uint32_t>::const_iterator deadModulesBegin = fDeadModules.begin();
0177     std::vector<uint32_t>::const_iterator deadModulesEnd = fDeadModules.end();
0178
0179     // Loop over entire tracker geometry.
0180     for (TrackerGeometry::DetContainer::const_iterator i = trackerGeo->dets().begin(); i != trackerGeo->dets().end();
0181          ++i) {
0182       // See if this is a pixel module.
0183       if (GeomDetEnumerators::isTrackerPixel((*i)->subDetector())) {
0184         DetId detId = (*i)->geographicalId();
0185
0186         // Skip this module if it's on the list of modules to be ignored.
0187         if (filterDeadModules && find(deadModulesBegin, deadModulesEnd, detId()) != deadModulesEnd) {
0188           continue;
0189         }
0190
0191         // Find all clusters in this module.
0192         edmNew::DetSetVector<SiPixelCluster>::const_iterator iSearch = pixelClusters->find(detId);
0193
0194         // Loop over all clusters in this module.
0195         if (iSearch != pixelClusters->end()) {
0196           for (edmNew::DetSet<SiPixelCluster>::const_iterator clus = iSearch->begin(); clus != iSearch->end(); ++clus) {
0197             if ((clus->size() >= fMinPixelsPerCluster) && (clus->charge() >= fMinClusterCharge)) {
0198               PixelGeomDetUnit const *theGeomDet = dynamic_cast<PixelGeomDetUnit const *>(trackerGeo->idToDet(detId));
0199               PixelTopology const *topol = &(theGeomDet->specificTopology());
0200               double x = clus->x();
0201               double y = clus->y();
0202               LocalPoint clustLP = topol->localPosition(MeasurementPoint(x, y));
0203               GlobalPoint clustGP = theGeomDet->surface().toGlobal(clustLP);
0204               double charge = clus->charge() / 1.e3;
0205               int size = clus->size();
0206
0207               if (detId.subdetId() == PixelSubdetector::PixelBarrel) {
0208                 PixelBarrelNameUpgrade detName = PixelBarrelNameUpgrade(detId);
0209                 unsigned int layer = detName.layerName() - kOffsetLayers;
0210                 if (layer < kNumLayers) {
0211                   std::string histName;
0212                   histName = "clusPosBarrel" + std::to_string(layer);
0213                   fHistContainerThisRun[histName]->Fill(clustGP.z(), clustGP.phi());
0214                   histName = "clusChargeBarrel" + std::to_string(layer);
0215                   fHistContainerThisRun[histName]->Fill(iEvent.bunchCrossing(), charge);
0216                   histName = "clusSizeBarrel" + std::to_string(layer);
0217                   fHistContainerThisRun[histName]->Fill(iEvent.bunchCrossing(), size);
0218                 } else {
0219                   edm::LogWarning("pixelLumi") << "higher layer number, " << layer << ", than layers";
0220                 }
0221               } else {
0222                 // DEBUG DEBUG DEBUG
0223                 assert(detId.subdetId() == PixelSubdetector::PixelEndcap);
0224                 // DEBUG DEBUG DEBUG end
0225
0226                 PixelEndcapNameUpgrade detName = PixelEndcapNameUpgrade(detId);
0227                 unsigned int disk = detName.diskName() - kOffsetDisks;
0228                 if (disk < kNumDisks) {
0229                   std::string histName;
0230                   histName = "clusPosEndCap" + std::to_string(disk);
0231                   fHistContainerThisRun[histName]->Fill(clustGP.x(), clustGP.y());
0232                   histName = "clusChargeEndCap" + std::to_string(disk);
0233                   fHistContainerThisRun[histName]->Fill(iEvent.bunchCrossing(), charge);
0234                   histName = "clusSizeEndCap" + std::to_string(disk);
0235                   fHistContainerThisRun[histName]->Fill(iEvent.bunchCrossing(), size);
0236                 } else {
0237                   edm::LogWarning("pixelLumi")
0238                       << "higher disk number, " << disk << ", than disks," << kNumDisks << std::endl;
0239                 }
0240               }
0241             }
0242           }
0243         }
0244       }
0245     }
0246   }
0247 }
0248
0249 void PixelLumiDQM::bookHistograms(DQMStore::IBooker &ibooker,
0250                                   edm::Run const &run,
0251                                   edm::EventSetup const & /* iSetup */) {
0252   edm::LogInfo("Status") << "Starting processing of run #" << run.id().run();
0253
0254   // Top folder containing high-level information about pixel and HF lumi.
0255   std::string folder = "PixelLumi/";
0256   ibooker.setCurrentFolder(folder);
0257
0258   fHistTotalRecordedLumiByLS = ibooker.book1D("totalPixelLumiByLS", "Pixel Lumi in nb vs LS", 8000, 0.5, 8000.5);
0259   fHistRecordedByBxCumulative = ibooker.book1D("PXLumiByBXsum",
0260                                                "Pixel Lumi in nb by BX Cumulative vs LS",
0261                                                lastBunchCrossing,
0262                                                0.5,
0263                                                float(lastBunchCrossing) + 0.5);
0264
0265   std::string subfolder = folder + "lastLS/";
0266   ibooker.setCurrentFolder(subfolder);
0267   fHistRecordedByBxLastLumi = ibooker.book1D(
0268       "PXByBXLastLumi", "Pixel By BX Last Lumi", lastBunchCrossing + 1, -0.5, float(lastBunchCrossing) + 0.5);
0269
0270   subfolder = folder + "ClusterCountingDetails/";
0271   ibooker.setCurrentFolder(subfolder);
0272
0273   fHistnBClusVsLS[0] = ibooker.book1D("nBClusVsLS_0", "Fraction of Clusters vs LS Barrel layer 0", 8000, 0.5, 8000.5);
0274   fHistnBClusVsLS[1] = ibooker.book1D("nBClusVsLS_1", "Fraction of Clusters vs LS Barrel layer 1", 8000, 0.5, 8000.5);
0275   fHistnBClusVsLS[2] = ibooker.book1D("nBClusVsLS_2", "Fraction of Clusters vs LS Barrel layer 2", 8000, 0.5, 8000.5);
0276   fHistnFPClusVsLS[0] = ibooker.book1D("nFPClusVsLS_0", "Fraction of Clusters vs LS Barrel layer 0", 8000, 0.5, 8000.5);
0277   fHistnFPClusVsLS[1] = ibooker.book1D("nFPClusVsLS_1", "Fraction of Clusters vs LS Barrel layer 1", 8000, 0.5, 8000.5);
0278   fHistnFMClusVsLS[0] = ibooker.book1D("nFMClusVsLS_0", "Fraction of Clusters vs LS Barrel layer 0", 8000, 0.5, 8000.5);
0279   fHistnFMClusVsLS[1] = ibooker.book1D("nFMClusVsLS_1", "Fraction of Clusters vs LS Barrel layer 1", 8000, 0.5, 8000.5);
0280   fHistBunchCrossings = ibooker.book1D(
0281       "BunchCrossings", "Cumulative Bunch Crossings", lastBunchCrossing, 0.5, float(lastBunchCrossing) + 0.5);
0282   fHistBunchCrossingsLastLumi = ibooker.book1D(
0283       "BunchCrossingsLL", "Bunch Crossings Last Lumi", lastBunchCrossing, 0.5, float(lastBunchCrossing) + 0.5);
0284   fHistClusterCountByBxLastLumi = ibooker.book1D(
0285       "ClusterCountByBxLL", "Cluster Count by BX Last Lumi", lastBunchCrossing, 0.5, float(lastBunchCrossing) + 0.5);
0286   fHistClusterCountByBxCumulative = ibooker.book1D(
0287       "ClusterCountByBxSum", "Cluster Count by BX Cumulative", lastBunchCrossing, 0.5, float(lastBunchCrossing) + 0.5);
0288   fHistClusByLS = ibooker.book1D("totalClusByLS", "Number of Clusters all dets vs LS", 8000, 0.5, 8000.5);
0289
0290   // Add some pixel cluster quality check histograms (in a subfolder).
0291   subfolder = folder + "qualityChecks/";
0292   ibooker.setCurrentFolder(subfolder);
0293
0294   if (fIncludePixelQualCheckHistos) {
0295     // Create histograms for this run if not already present in our list.
0296     edm::LogInfo("Status") << "Creating histograms for run #" << run.id().run();
0297
0298     // Pixel cluster positions in the barrel - (z, phi).
0299     for (size_t i = 0; i <= kNumLayers; ++i) {
0300       std::stringstream key;
0301       key << "clusPosBarrel" << i;
0302       std::stringstream name;
0303       name << key.str() << "_" << run.run();
0304       std::stringstream title;
0305       title << "Pixel cluster position - barrel layer " << i;
0306       fHistContainerThisRun[key.str()] =
0307           ibooker.book2D(name.str().c_str(), title.str().c_str(), 100, -30., 30., 64, -Geom::pi(), Geom::pi());
0308     }
0309
0310     // Pixel cluster positions in the endcaps (x, y).
0311     for (size_t i = 0; i <= kNumDisks; ++i) {
0312       std::stringstream key;
0313       key << "clusPosEndCap" << i;
0314       std::stringstream name;
0315       name << key.str() << "_" << run.run();
0316       std::stringstream title;
0317       title << "Pixel cluster position - endcap disk " << i;
0318       fHistContainerThisRun[key.str()] =
0319           ibooker.book2D(name.str().c_str(), title.str().c_str(), 100, -20., 20., 100, -20., 20.);
0320     }
0321
0322     // Pixel cluster charge in the barrel, per bx.
0323     for (size_t i = 0; i <= kNumLayers; ++i) {
0324       std::stringstream key;
0325       key << "clusChargeBarrel" << i;
0326       std::stringstream name;
0327       name << key.str() << "_" << run.run();
0328       std::stringstream title;
0329       title << "Pixel cluster charge - barrel layer " << i;
0330       fHistContainerThisRun[key.str()] =
0331           ibooker.book2D(name.str().c_str(), title.str().c_str(), 3564, .5, 3564.5, 100, 0., 100.);
0332     }
0333
0334     // Pixel cluster charge in the endcaps, per bx.
0335     for (size_t i = 0; i <= kNumDisks; ++i) {
0336       std::stringstream key;
0337       key << "clusChargeEndCap" << i;
0338       std::stringstream name;
0339       name << key.str() << "_" << run.run();
0340       std::stringstream title;
0341       title << "Pixel cluster charge - endcap disk " << i;
0342       fHistContainerThisRun[key.str()] =
0343           ibooker.book2D(name.str().c_str(), title.str().c_str(), 3564, .5, 3564.5, 100, 0., 100.);
0344     }
0345
0346     // Pixel cluster size in the barrel, per bx.
0347     for (size_t i = 0; i <= kNumLayers; ++i) {
0348       std::stringstream key;
0349       key << "clusSizeBarrel" << i;
0350       std::stringstream name;
0351       name << key.str() << "_" << run.run();
0352       std::stringstream title;
0353       title << "Pixel cluster size - barrel layer " << i;
0354       fHistContainerThisRun[key.str()] =
0355           ibooker.book2D(name.str().c_str(), title.str().c_str(), 3564, .5, 3564.5, 100, 0., 100.);
0356     }
0357
0358     // Pixel cluster size in the endcaps, per bx.
0359     for (size_t i = 0; i <= kNumDisks; ++i) {
0360       std::stringstream key;
0361       key << "clusSizeEndCap" << i;
0362       std::stringstream name;
0363       name << key.str() << "_" << run.run();
0364       std::stringstream title;
0365       title << "Pixel cluster size - endcap disk " << i;
0366       fHistContainerThisRun[key.str()] =
0367           ibooker.book2D(name.str().c_str(), title.str().c_str(), 3564, .5, 3564.5, 100, 0., 100.);
0368     }
0369   }
0370 }
0371
0372 // ------------ Method called when starting to process a luminosity block.
0373 // ------------
0374 void PixelLumiDQM::beginLuminosityBlock(edm::LuminosityBlock const &lumiBlock, edm::EventSetup const &) {
0375   // Only reset and fill every fResetIntervalInLumiSections (default is 1 LS)
0376   // Return unless the PREVIOUS LS was at the right modulo value
0377   // (e.g. is resetinterval = 5 the rest will only be executed at LS=6
0378   // NB: reset is done here so the histograms by LS are sent before resetting.
0379   // NB: not being used for now since default is 1 LS. There is a bug here.
0380
0381   unsigned int ls = lumiBlock.luminosityBlockAuxiliary().luminosityBlock();
0382
0383   if ((ls - 1) % fResetIntervalInLumiSections == 0) {
0384     fHistBunchCrossingsLastLumi->Reset();
0385     fHistClusterCountByBxLastLumi->Reset();
0386     fHistRecordedByBxLastLumi->Reset();
0387   }
0388 }
0389
0390 // ------------ Method called when ending the processing of a luminosity block.
0391 // ------------
0392 void PixelLumiDQM::endLuminosityBlock(edm::LuminosityBlock const &lumiBlock, edm::EventSetup const &es) {
0393   unsigned int ls = lumiBlock.luminosityBlockAuxiliary().luminosityBlock();
0394
0395   // Only fill every fResetIntervalInLumiSections (default is 1 LS)
0396   if (ls % fResetIntervalInLumiSections != 0)
0397     return;
0398
0399   printf("Lumi Block = %d\n", ls);
0400
0401   if ((ls - 1) % fResetIntervalInLumiSections == 0) {
0402   }
0403
0404   unsigned int nBClus[3] = {0, 0, 0};
0405   unsigned int nFPClus[2] = {0, 0};
0406   unsigned int nFMClus[2] = {0, 0};
0407
0408   double total_recorded = 0.;
0409   double total_recorded_unc_square = 0.;
0410
0411   // Obtain bunch-by-bunch cluster counts and compute totals for lumi
0412   // calculation.
0413   double totalcounts = 0.0;
0414   double totalevents = 0.0;
0415   double lumi_factor_per_bx = 0.0;
0416   if (useInnerBarrelLayer)
0417     lumi_factor_per_bx = FREQ_ORBIT * SECONDS_PER_LS * fResetIntervalInLumiSections / XSEC_PIXEL_CLUSTER;
0418   else
0419     lumi_factor_per_bx = FREQ_ORBIT * SECONDS_PER_LS * fResetIntervalInLumiSections / rXSEC_PIXEL_CLUSTER;
0420
0421   for (std::map<int, PixelClusterCount>::iterator it = fNumPixelClusters.begin(); it != fNumPixelClusters.end(); it++) {
0422     // Sum all clusters for this BX.
0423     unsigned int total = (*it).second.numB.at(1) + (*it).second.numB.at(2) + (*it).second.numFP.at(0) +
0424                          (*it).second.numFP.at(1) + (*it).second.numFM.at(0) + (*it).second.numFM.at(1);
0425     if (useInnerBarrelLayer)
0426       total += (*it).second.numB.at(0);
0427     totalcounts += total;
0428     double etotal = useInnerBarrelLayer
0429                         ? (*it).second.dnumB.at(0)
0430                         : (*it).second.dnumB.at(1) + (*it).second.dnumB.at(2) + (*it).second.dnumFP.at(0) +
0431                               (*it).second.dnumFP.at(1) + (*it).second.dnumFM.at(0) + (*it).second.dnumFM.at(1);
0432     etotal = sqrt(etotal);
0433
0434     fHistClusterCountByBxLastLumi->setBinContent((*it).first, total);
0435     fHistClusterCountByBxLastLumi->setBinError((*it).first, etotal);
0436     fHistClusterCountByBxCumulative->setBinContent((*it).first,
0437                                                    fHistClusterCountByBxCumulative->getBinContent((*it).first) + total);
0438
0439     unsigned int events_per_bx = fHistBunchCrossingsLastLumi->getBinContent((*it).first);
0440     totalevents += events_per_bx;
0441     double average_cluster_count = events_per_bx != 0 ? double(total) / events_per_bx : 0.;
0442     double average_cluster_count_unc = events_per_bx != 0 ? etotal / events_per_bx : 0.;
0443     double pixel_bx_lumi_per_ls = lumi_factor_per_bx * average_cluster_count / CM2_TO_NANOBARN;
0444     double pixel_bx_lumi_per_ls_unc = 0.0;
0445     if (useInnerBarrelLayer)
0446       pixel_bx_lumi_per_ls_unc = sqrt(lumi_factor_per_bx * lumi_factor_per_bx *
0447                                       (average_cluster_count_unc * average_cluster_count_unc +
0448                                        (average_cluster_count * XSEC_PIXEL_CLUSTER_UNC / XSEC_PIXEL_CLUSTER) *
0449                                            (average_cluster_count * XSEC_PIXEL_CLUSTER / XSEC_PIXEL_CLUSTER))) /
0450                                  CM2_TO_NANOBARN;
0451     else
0452       pixel_bx_lumi_per_ls_unc = sqrt(lumi_factor_per_bx * lumi_factor_per_bx *
0453                                       (average_cluster_count_unc * average_cluster_count_unc +
0454                                        (average_cluster_count * rXSEC_PIXEL_CLUSTER_UNC / rXSEC_PIXEL_CLUSTER) *
0455                                            (average_cluster_count * rXSEC_PIXEL_CLUSTER / rXSEC_PIXEL_CLUSTER))) /
0456                                  CM2_TO_NANOBARN;
0457
0458     fHistRecordedByBxLastLumi->setBinContent((*it).first, pixel_bx_lumi_per_ls);
0459     fHistRecordedByBxLastLumi->setBinError((*it).first, pixel_bx_lumi_per_ls_unc);
0460
0461     fHistRecordedByBxCumulative->setBinContent(
0462         (*it).first, fHistRecordedByBxCumulative->getBinContent((*it).first) + pixel_bx_lumi_per_ls);
0463
0464     /*
0465       if(fHistRecordedByBxLastLumi->getBinContent((*it).first)!=0.)
0466       fHistRecordedByBxLastLumi->getBinContent((*it).first));
0467       if(fHistRecordedByBxCumulative->getBinContent((*it).first)!=0.)
0468       fHistRecordedByBxCumulative->getBinContent((*it).first));
0469     */
0470
0471     nBClus[0] += (*it).second.numB.at(0);
0472     nBClus[1] += (*it).second.numB.at(1);
0473     nBClus[2] += (*it).second.numB.at(2);
0474     nFPClus[0] += (*it).second.numFP.at(0);
0475     nFPClus[1] += (*it).second.numFP.at(1);
0476     nFMClus[0] += (*it).second.numFM.at(0);
0477     nFMClus[1] += (*it).second.numFM.at(1);
0478
0479     // Reset counters
0480     (*it).second.Reset();
0481
0482     // edm::LogWarning("pixelLumi") << "bx="<< (*it).first << " clusters=" <<
0483     // (*it).second.numB.at(0));
0484   }
0485
0486   if ((filledAndUnmaskedBunches = calculateBunchMask(fHistClusterCountByBxCumulative, bunchTriggerMask)) != 0) {
0487     for (unsigned int i = 0; i <= lastBunchCrossing; i++) {
0488       if (bunchTriggerMask[i]) {
0489         double err = fHistRecordedByBxLastLumi->getBinError(i);
0490         total_recorded += fHistRecordedByBxLastLumi->getBinContent(i);
0491         total_recorded_unc_square += err * err;
0492       }
0493     }
0494
0495     // Replace the total obtained by summing over BXs with the average per BX
0496     // from the total cluster count and rescale
0497
0498     if (totalevents > 10) {
0499       total_recorded = lumi_factor_per_bx * totalcounts / totalevents / CM2_TO_NANOBARN;
0500     } else
0501       total_recorded = 0.0;
0502
0503     edm::LogWarning("pixelLumi") << " Total recorded " << total_recorded;
0504     fHistTotalRecordedLumiByLS->setBinContent(ls, total_recorded);
0505     fHistTotalRecordedLumiByLS->setBinError(ls, sqrt(total_recorded_unc_square));
0506   }
0507   // fill cluster counts by detector regions for sanity checks
0508   unsigned int all_detectors_counts = 0;
0509   for (unsigned int i = 0; i < 3; i++) {
0510     all_detectors_counts += nBClus[i];
0511   }
0512   for (unsigned int i = 0; i < 2; i++) {
0513     all_detectors_counts += nFPClus[i];
0514   }
0515   for (unsigned int i = 0; i < 2; i++) {
0516     all_detectors_counts += nFMClus[i];
0517   }
0518
0519   fHistClusByLS->setBinContent(ls, all_detectors_counts);
0520
0521   for (unsigned int i = 0; i < 3; i++) {
0522     fHistnBClusVsLS[i]->setBinContent(ls, float(nBClus[i]) / float(all_detectors_counts));
0523   }
0524   for (unsigned int i = 0; i < 2; i++) {
0525     fHistnFPClusVsLS[i]->setBinContent(ls, float(nFPClus[i]) / float(all_detectors_counts));
0526   }
0527   for (unsigned int i = 0; i < 2; i++) {
0528     fHistnFMClusVsLS[i]->setBinContent(ls, float(nFMClus[i]) / float(all_detectors_counts));
0529   }
0530
0531   logFile_.open(fLogFileName_.c_str(), std::ios_base::trunc);
0532
0533   timeval tv;
0534   gettimeofday(&tv, nullptr);
0535   tm *ts = gmtime(&tv.tv_sec);
0536   char datestring[256];
0537   strftime(datestring, sizeof(datestring), "%Y.%m.%d %T GMT %s", ts);
0538   logFile_ << "RunNumber " << fRunNo << std::endl;
0539   logFile_ << "EndTimeOfFit " << datestring << std::endl;
0540   logFile_ << "LumiRange " << ls << "-" << ls << std::endl;
0541   logFile_ << "Fill " << -99 << std::endl;
0542   logFile_ << "ActiveBunchCrossings " << filledAndUnmaskedBunches << std::endl;
0543   logFile_ << "PixelLumi " << fHistTotalRecordedLumiByLS->getBinContent(ls) * 0.98 << std::endl;
0544   logFile_ << "HFLumi " << -99 << std::endl;
0545   logFile_ << "Ratio " << -99 << std::endl;
0546   logFile_.close();
0547 }
0548
0549 unsigned int PixelLumiDQM::calculateBunchMask(MonitorElement *e, std::vector<bool> &mask) {
0550   unsigned int nbins = e->getNbinsX();
0551   std::vector<float> ar(nbins + 1, 0.);
0552   for (unsigned int i = 1; i <= nbins; i++) {
0553     ar[i] = e->getBinContent(i);
0554   }
0555   return calculateBunchMask(ar, nbins, mask);
0556 }
0557 unsigned int PixelLumiDQM::calculateBunchMask(std::vector<float> &e, unsigned int nbins, std::vector<bool> &mask) {
0558   // Take the cumulative cluster count histogram and find max and average of
0559   // non-empty bins.
0560   unsigned int active_count = 0;
0561   double maxc = 0.0;
0562   double ave = 0.0;  // Average of non-empty bins
0563   unsigned int non_empty_bins = 0;
0564
0565   for (unsigned int i = 1; i <= nbins; i++) {
0566     double bin = e[i];
0567     if (bin != 0.0) {
0568       if (maxc < bin)
0569         maxc = bin;
0570       ave += bin;
0571       non_empty_bins++;
0572     }
0573   }
0574
0575   ave /= non_empty_bins;
0576   edm::LogWarning("pixelLumi") << "Bunch mask finder - non empty bins " << non_empty_bins
0577                                << " average of non empty bins " << ave << " max content of one bin " << maxc;
0578   double mean = 0.;
0579   double sigma = 0.;
0580   if (non_empty_bins < 50) {
0581     mean = maxc;
0582     sigma = (maxc) / 20;
0583   } else {
0584     TH1F dist("dist", "dist", 500, 0., maxc + (maxc / 500.) * 20.);
0585     for (unsigned int i = 1; i <= nbins; i++) {
0586       double bin = e[i];
0587       dist.Fill(bin);
0588     }
0589     TF1 fit("plgaus", "gaus");
0590     dist.Fit(&fit, "", "", fmax(0., ave - (maxc - ave) / 5.), maxc);
0591     mean = fit.GetParameter("Mean");
0592     sigma = fit.GetParameter("Sigma");
0593   }
0594   edm::LogWarning("pixelLumi") << "Bunch mask will use mean" << mean << " sigma " << sigma;
0595   // Active BX defined as those which have nclus within fixed standard
0596   // deviations of peak.
0597   for (unsigned int i = 1; i <= nbins; i++) {
0598     double bin = e[i];
0599     if (bin > 0. && std::abs(bin - mean) < 5. * (sigma)) {
0600       mask[i] = true;
0601       active_count++;
0602     }
0603   }
0604   edm::LogWarning("pixelLumi") << "Bunch mask finds " << active_count << " active bunch crossings ";
0605   //   edm::LogWarning("pixelLumi") << "this is the full bx mask " ;
0606   //   for(unsigned int i = 1; i<= nbins; i++)
0607   //     edm::LogWarning("pixelLumi") << ((mask[i]) ? 1:0);
0608   return active_count;
0609 }
0610 // Define this as a CMSSW plug-in.
0611 DEFINE_FWK_MODULE(PixelLumiDQM);