File indexing completed on 2023-03-17 10:55:53
0001
0002
0003
0004
0005
0006
0007
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
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
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
0099 std::map<int, PixelClusterCount>::iterator it = fNumPixelClusters.find(fBXNo);
0100 if (it == fNumPixelClusters.end())
0101 fNumPixelClusters[fBXNo] = PixelClusterCount();
0102
0103 const TrackerGeometry *trackerGeo = &iSetup.getData(tkGeomToken_);
0104 if (fIncludePixelClusterInfo) {
0105
0106 edm::Handle<edmNew::DetSetVector<SiPixelCluster>> pixelClusters;
0107 iEvent.getByToken(fPixelClusterLabel, pixelClusters);
0108
0109
0110 for (TrackerGeometry::DetContainer::const_iterator i = trackerGeo->dets().begin(); i != trackerGeo->dets().end();
0111 ++i) {
0112
0113
0114 if (GeomDetEnumerators::isTrackerPixel((*i)->subDetector())) {
0115 DetId detId = (*i)->geographicalId();
0116
0117 edmNew::DetSetVector<SiPixelCluster>::const_iterator iSearch = pixelClusters->find(detId);
0118 if (iSearch != pixelClusters->end()) {
0119
0120
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
0129 assert(numClusters <= iSearch->size());
0130
0131
0132
0133
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
0141 assert(detId.subdetId() == PixelSubdetector::PixelEndcap);
0142
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
0170 if (fIncludePixelQualCheckHistos) {
0171
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
0180 for (TrackerGeometry::DetContainer::const_iterator i = trackerGeo->dets().begin(); i != trackerGeo->dets().end();
0181 ++i) {
0182
0183 if (GeomDetEnumerators::isTrackerPixel((*i)->subDetector())) {
0184 DetId detId = (*i)->geographicalId();
0185
0186
0187 if (filterDeadModules && find(deadModulesBegin, deadModulesEnd, detId()) != deadModulesEnd) {
0188 continue;
0189 }
0190
0191
0192 edmNew::DetSetVector<SiPixelCluster>::const_iterator iSearch = pixelClusters->find(detId);
0193
0194
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
0223 assert(detId.subdetId() == PixelSubdetector::PixelEndcap);
0224
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 & ) {
0252 edm::LogInfo("Status") << "Starting processing of run #" << run.id().run();
0253
0254
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
0291 subfolder = folder + "qualityChecks/";
0292 ibooker.setCurrentFolder(subfolder);
0293
0294 if (fIncludePixelQualCheckHistos) {
0295
0296 edm::LogInfo("Status") << "Creating histograms for run #" << run.id().run();
0297
0298
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
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
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
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
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
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
0373
0374 void PixelLumiDQM::beginLuminosityBlock(edm::LuminosityBlock const &lumiBlock, edm::EventSetup const &) {
0375
0376
0377
0378
0379
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
0391
0392 void PixelLumiDQM::endLuminosityBlock(edm::LuminosityBlock const &lumiBlock, edm::EventSetup const &es) {
0393 unsigned int ls = lumiBlock.luminosityBlockAuxiliary().luminosityBlock();
0394
0395
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
0412
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
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
0466
0467
0468
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
0480 (*it).second.Reset();
0481
0482
0483
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
0496
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
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
0559
0560 unsigned int active_count = 0;
0561 double maxc = 0.0;
0562 double ave = 0.0;
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
0596
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
0606
0607
0608 return active_count;
0609 }
0610
0611 DEFINE_FWK_MODULE(PixelLumiDQM);