Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-01 03:53:37

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 etotalcounts = 0.0;
0415   double totalevents = 0.0;
0416   double lumi_factor_per_bx = 0.0;
0417   if (useInnerBarrelLayer)
0418     lumi_factor_per_bx = FREQ_ORBIT * SECONDS_PER_LS * fResetIntervalInLumiSections / XSEC_PIXEL_CLUSTER;
0419   else
0420     lumi_factor_per_bx = FREQ_ORBIT * SECONDS_PER_LS * fResetIntervalInLumiSections / rXSEC_PIXEL_CLUSTER;
0421 
0422   for (std::map<int, PixelClusterCount>::iterator it = fNumPixelClusters.begin(); it != fNumPixelClusters.end(); it++) {
0423     // Sum all clusters for this BX.
0424     unsigned int total = (*it).second.numB.at(1) + (*it).second.numB.at(2) + (*it).second.numFP.at(0) +
0425                          (*it).second.numFP.at(1) + (*it).second.numFM.at(0) + (*it).second.numFM.at(1);
0426     if (useInnerBarrelLayer)
0427       total += (*it).second.numB.at(0);
0428     totalcounts += total;
0429     double etotal = (*it).second.dnumB.at(1) + (*it).second.dnumB.at(2) + (*it).second.dnumFP.at(0) +
0430                     (*it).second.dnumFP.at(1) + (*it).second.dnumFM.at(0) + (*it).second.dnumFM.at(1);
0431     if (useInnerBarrelLayer)
0432       etotal = (*it).second.dnumB.at(0);
0433     etotalcounts += etotal;
0434     etotal = sqrt(etotal);
0435 
0436     fHistClusterCountByBxLastLumi->setBinContent((*it).first, total);
0437     fHistClusterCountByBxLastLumi->setBinError((*it).first, etotal);
0438     fHistClusterCountByBxCumulative->setBinContent((*it).first,
0439                                                    fHistClusterCountByBxCumulative->getBinContent((*it).first) + total);
0440 
0441     unsigned int events_per_bx = fHistBunchCrossingsLastLumi->getBinContent((*it).first);
0442     totalevents += events_per_bx;
0443     double average_cluster_count = events_per_bx != 0 ? double(total) / events_per_bx : 0.;
0444     double average_cluster_count_unc = events_per_bx != 0 ? etotal / events_per_bx : 0.;
0445     double pixel_bx_lumi_per_ls = lumi_factor_per_bx * average_cluster_count / CM2_TO_NANOBARN;
0446     double pixel_bx_lumi_per_ls_unc = 0.0;
0447     if (useInnerBarrelLayer)
0448       pixel_bx_lumi_per_ls_unc = sqrt(lumi_factor_per_bx * lumi_factor_per_bx *
0449                                       (average_cluster_count_unc * average_cluster_count_unc +
0450                                        (average_cluster_count * XSEC_PIXEL_CLUSTER_UNC / XSEC_PIXEL_CLUSTER) *
0451                                            (average_cluster_count * XSEC_PIXEL_CLUSTER / XSEC_PIXEL_CLUSTER))) /
0452                                  CM2_TO_NANOBARN;
0453     else
0454       pixel_bx_lumi_per_ls_unc = sqrt(lumi_factor_per_bx * lumi_factor_per_bx *
0455                                       (average_cluster_count_unc * average_cluster_count_unc +
0456                                        (average_cluster_count * rXSEC_PIXEL_CLUSTER_UNC / rXSEC_PIXEL_CLUSTER) *
0457                                            (average_cluster_count * rXSEC_PIXEL_CLUSTER / rXSEC_PIXEL_CLUSTER))) /
0458                                  CM2_TO_NANOBARN;
0459 
0460     fHistRecordedByBxLastLumi->setBinContent((*it).first, pixel_bx_lumi_per_ls);
0461     fHistRecordedByBxLastLumi->setBinError((*it).first, pixel_bx_lumi_per_ls_unc);
0462 
0463     fHistRecordedByBxCumulative->setBinContent(
0464         (*it).first, fHistRecordedByBxCumulative->getBinContent((*it).first) + pixel_bx_lumi_per_ls);
0465 
0466     /*
0467       if(fHistRecordedByBxLastLumi->getBinContent((*it).first)!=0.)
0468       fHistRecordedByBxLastLumi->getBinContent((*it).first));
0469       if(fHistRecordedByBxCumulative->getBinContent((*it).first)!=0.)
0470       fHistRecordedByBxCumulative->getBinContent((*it).first));
0471     */
0472 
0473     nBClus[0] += (*it).second.numB.at(0);
0474     nBClus[1] += (*it).second.numB.at(1);
0475     nBClus[2] += (*it).second.numB.at(2);
0476     nFPClus[0] += (*it).second.numFP.at(0);
0477     nFPClus[1] += (*it).second.numFP.at(1);
0478     nFMClus[0] += (*it).second.numFM.at(0);
0479     nFMClus[1] += (*it).second.numFM.at(1);
0480 
0481     // Reset counters
0482     (*it).second.Reset();
0483 
0484     // edm::LogWarning("pixelLumi") << "bx="<< (*it).first << " clusters=" <<
0485     // (*it).second.numB.at(0));
0486   }
0487 
0488   if ((filledAndUnmaskedBunches = calculateBunchMask(fHistClusterCountByBxCumulative, bunchTriggerMask)) != 0) {
0489     for (unsigned int i = 0; i <= lastBunchCrossing; i++) {
0490       if (bunchTriggerMask[i]) {
0491         double err = fHistRecordedByBxLastLumi->getBinError(i);
0492         total_recorded += fHistRecordedByBxLastLumi->getBinContent(i);
0493         total_recorded_unc_square += err * err;
0494       }
0495     }
0496 
0497     // Replace the total obtained by summing over BXs with the average per BX
0498     // from the total cluster count and rescale
0499 
0500     if (totalevents > 10) {
0501       total_recorded = lumi_factor_per_bx * totalcounts / totalevents / CM2_TO_NANOBARN;
0502     } else
0503       total_recorded = 0.0;
0504 
0505     edm::LogWarning("pixelLumi") << " Total recorded " << total_recorded;
0506     fHistTotalRecordedLumiByLS->setBinContent(ls, total_recorded);
0507     fHistTotalRecordedLumiByLS->setBinError(ls, sqrt(total_recorded_unc_square));
0508   }
0509   // fill cluster counts by detector regions for sanity checks
0510   unsigned int all_detectors_counts = 0;
0511   for (unsigned int i = 0; i < 3; i++) {
0512     all_detectors_counts += nBClus[i];
0513   }
0514   for (unsigned int i = 0; i < 2; i++) {
0515     all_detectors_counts += nFPClus[i];
0516   }
0517   for (unsigned int i = 0; i < 2; i++) {
0518     all_detectors_counts += nFMClus[i];
0519   }
0520 
0521   fHistClusByLS->setBinContent(ls, all_detectors_counts);
0522 
0523   for (unsigned int i = 0; i < 3; i++) {
0524     fHistnBClusVsLS[i]->setBinContent(ls, float(nBClus[i]) / float(all_detectors_counts));
0525   }
0526   for (unsigned int i = 0; i < 2; i++) {
0527     fHistnFPClusVsLS[i]->setBinContent(ls, float(nFPClus[i]) / float(all_detectors_counts));
0528   }
0529   for (unsigned int i = 0; i < 2; i++) {
0530     fHistnFMClusVsLS[i]->setBinContent(ls, float(nFMClus[i]) / float(all_detectors_counts));
0531   }
0532 
0533   logFile_.open(fLogFileName_.c_str(), std::ios_base::trunc);
0534 
0535   timeval tv;
0536   gettimeofday(&tv, nullptr);
0537   tm *ts = gmtime(&tv.tv_sec);
0538   char datestring[256];
0539   strftime(datestring, sizeof(datestring), "%Y.%m.%d %T GMT %s", ts);
0540   logFile_ << "RunNumber " << fRunNo << std::endl;
0541   logFile_ << "EndTimeOfFit " << datestring << std::endl;
0542   logFile_ << "LumiRange " << ls << "-" << ls << std::endl;
0543   logFile_ << "Fill " << -99 << std::endl;
0544   logFile_ << "ActiveBunchCrossings " << filledAndUnmaskedBunches << std::endl;
0545   logFile_ << "PixelLumi " << fHistTotalRecordedLumiByLS->getBinContent(ls) * 0.98 << std::endl;
0546   logFile_ << "HFLumi " << -99 << std::endl;
0547   logFile_ << "Ratio " << -99 << std::endl;
0548   logFile_.close();
0549 }
0550 
0551 unsigned int PixelLumiDQM::calculateBunchMask(MonitorElement *e, std::vector<bool> &mask) {
0552   unsigned int nbins = e->getNbinsX();
0553   std::vector<float> ar(nbins + 1, 0.);
0554   for (unsigned int i = 1; i <= nbins; i++) {
0555     ar[i] = e->getBinContent(i);
0556   }
0557   return calculateBunchMask(ar, nbins, mask);
0558 }
0559 unsigned int PixelLumiDQM::calculateBunchMask(std::vector<float> &e, unsigned int nbins, std::vector<bool> &mask) {
0560   // Take the cumulative cluster count histogram and find max and average of
0561   // non-empty bins.
0562   unsigned int active_count = 0;
0563   double maxc = 0.0;
0564   double ave = 0.0;  // Average of non-empty bins
0565   unsigned int non_empty_bins = 0;
0566 
0567   for (unsigned int i = 1; i <= nbins; i++) {
0568     double bin = e[i];
0569     if (bin != 0.0) {
0570       if (maxc < bin)
0571         maxc = bin;
0572       ave += bin;
0573       non_empty_bins++;
0574     }
0575   }
0576 
0577   ave /= non_empty_bins;
0578   edm::LogWarning("pixelLumi") << "Bunch mask finder - non empty bins " << non_empty_bins
0579                                << " average of non empty bins " << ave << " max content of one bin " << maxc;
0580   double mean = 0.;
0581   double sigma = 0.;
0582   if (non_empty_bins < 50) {
0583     mean = maxc;
0584     sigma = (maxc) / 20;
0585   } else {
0586     TH1F dist("dist", "dist", 500, 0., maxc + (maxc / 500.) * 20.);
0587     for (unsigned int i = 1; i <= nbins; i++) {
0588       double bin = e[i];
0589       dist.Fill(bin);
0590     }
0591     TF1 fit("plgaus", "gaus");
0592     dist.Fit(&fit, "", "", fmax(0., ave - (maxc - ave) / 5.), maxc);
0593     mean = fit.GetParameter("Mean");
0594     sigma = fit.GetParameter("Sigma");
0595   }
0596   edm::LogWarning("pixelLumi") << "Bunch mask will use mean" << mean << " sigma " << sigma;
0597   // Active BX defined as those which have nclus within fixed standard
0598   // deviations of peak.
0599   for (unsigned int i = 1; i <= nbins; i++) {
0600     double bin = e[i];
0601     if (bin > 0. && std::abs(bin - mean) < 5. * (sigma)) {
0602       mask[i] = true;
0603       active_count++;
0604     }
0605   }
0606   edm::LogWarning("pixelLumi") << "Bunch mask finds " << active_count << " active bunch crossings ";
0607   //   edm::LogWarning("pixelLumi") << "this is the full bx mask " ;
0608   //   for(unsigned int i = 1; i<= nbins; i++)
0609   //     edm::LogWarning("pixelLumi") << ((mask[i]) ? 1:0);
0610   return active_count;
0611 }
0612 // Define this as a CMSSW plug-in.
0613 DEFINE_FWK_MODULE(PixelLumiDQM);