File indexing completed on 2024-09-26 05:05:42
0001
0002
0003
0004
0005 #include <memory>
0006 #include <string>
0007 #include <iostream>
0008 #include <fstream>
0009 #include <sstream>
0010
0011
0012 #include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h"
0013 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
0014 #include "CalibTracker/Records/interface/SiStripDetCablingRcd.h"
0015 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
0016 #include "CalibTracker/SiStripCommon/interface/SiStripDetInfoFileReader.h"
0017 #include "CalibTracker/SiStripHitEfficiency/interface/SiStripHitEfficiencyHelpers.h"
0018 #include "CalibTracker/SiStripHitEfficiency/interface/TrajectoryAtInvalidHit.h"
0019 #include "CalibTracker/SiStripHitEfficiency/plugins/HitEff.h"
0020 #include "CommonTools/ConditionDBWriter/interface/ConditionDBWriter.h"
0021 #include "CommonTools/TrackerMap/interface/TrackerMap.h"
0022 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0023 #include "DataFormats/Common/interface/DetSetVector.h"
0024 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0025 #include "DataFormats/Common/interface/Handle.h"
0026 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0027 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0028 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0029 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0030 #include "DataFormats/MuonReco/interface/Muon.h"
0031 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0032 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0033 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0034 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0035 #include "DataFormats/TrackReco/interface/DeDxData.h"
0036 #include "DataFormats/TrackReco/interface/Track.h"
0037 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0038 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0039 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0040 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0041 #include "FWCore/Framework/interface/Event.h"
0042 #include "FWCore/Framework/interface/EventSetup.h"
0043 #include "FWCore/Framework/interface/Frameworkfwd.h"
0044 #include "FWCore/Framework/interface/MakerMacros.h"
0045 #include "FWCore/ParameterSet/interface/FileInPath.h"
0046 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0047 #include "FWCore/ServiceRegistry/interface/Service.h"
0048 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0049 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0050 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0051 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0052 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0053 #include "RecoLocalTracker/ClusterParameterEstimator/interface/StripClusterParameterEstimator.h"
0054 #include "RecoLocalTracker/SiStripClusterizer/interface/SiStripClusterInfo.h"
0055 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0056 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0057 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0058 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0059 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0060
0061
0062 #include "TCanvas.h"
0063 #include "TEfficiency.h"
0064 #include "TF1.h"
0065 #include "TFile.h"
0066 #include "TGaxis.h"
0067 #include "TGraphAsymmErrors.h"
0068 #include "TH1F.h"
0069 #include "TH2F.h"
0070 #include "TLatex.h"
0071 #include "TLeaf.h"
0072 #include "TLegend.h"
0073 #include "TObjString.h"
0074 #include "TProfile.h"
0075 #include "TROOT.h"
0076 #include "TString.h"
0077 #include "TStyle.h"
0078 #include "TTree.h"
0079
0080
0081 #define LOGPRINT edm::LogPrint("SiStripHitEffFromCalibTree")
0082
0083 using namespace edm;
0084 using namespace reco;
0085 using namespace std;
0086
0087 struct hit {
0088 double x;
0089 double y;
0090 double z;
0091 unsigned int id;
0092 };
0093
0094 class SiStripHitEffFromCalibTree : public ConditionDBWriter<SiStripBadStrip> {
0095 public:
0096 explicit SiStripHitEffFromCalibTree(const edm::ParameterSet&);
0097 ~SiStripHitEffFromCalibTree() override = default;
0098
0099 private:
0100
0101 void algoAnalyze(const edm::Event& e, const edm::EventSetup& c) override;
0102 std::unique_ptr<SiStripBadStrip> getNewObject() override;
0103
0104
0105 void setBadComponents(int i,
0106 int component,
0107 SiStripQuality::BadComponent& BC,
0108 std::stringstream ssV[4][19],
0109 int NBadComponent[4][19][4]);
0110 void makeTKMap(bool autoTagging);
0111 void makeHotColdMaps();
0112 void makeSQLite();
0113 void totalStatistics();
0114 void makeSummary();
0115 void makeSummaryVsBx();
0116 void computeEff(vector<TH1F*>& vhfound, vector<TH1F*>& vhtotal, string name);
0117 void makeSummaryVsLumi();
0118 void makeSummaryVsCM();
0119 TString getLayerSideName(Long_t k);
0120
0121
0122 static constexpr int siStripLayers_ = 22;
0123 static constexpr double nBxInAnOrbit_ = 3565;
0124
0125 edm::Service<TFileService> fs;
0126 SiStripDetInfo detInfo_;
0127 edm::FileInPath fileInPath_;
0128 SiStripQuality* quality_;
0129
0130 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0131 const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0132
0133 TTree* calibTree_;
0134 vector<string> calibTreeFileNames_;
0135 float threshold_;
0136 unsigned int nModsMin_;
0137 string badModulesFile_;
0138 bool autoIneffModTagging_;
0139 unsigned int clusterMatchingMethod_;
0140 float resXSig_;
0141 float clusterTrajDist_;
0142 float stripsApvEdge_;
0143 bool useOnlyHighPurityTracks_;
0144 unsigned int bunchX_;
0145 unsigned int spaceBetweenTrains_;
0146 bool useCM_;
0147 bool showEndcapSides_;
0148 bool showRings_;
0149 bool showTOB6TEC9_;
0150 bool showOnlyGoodModules_;
0151 float tkMapMin_;
0152 float effPlotMin_;
0153 TString title_;
0154 bool storeModEff_;
0155
0156 unsigned int nTEClayers;
0157
0158 TH1F* bxHisto;
0159 TH1F* instLumiHisto;
0160 TH1F* PUHisto;
0161
0162 TH1F* bxHisto_cutOnTracks;
0163 TH1F* instLumiHisto_cutOnTracks;
0164 TH1F* PUHisto_cutOnTracks;
0165
0166
0167 map<pair<unsigned int, unsigned int>, array<double, 3> > eventInfos;
0168
0169 map<pair<unsigned int, unsigned int>, bool> event_used;
0170
0171 vector<hit> hits[23];
0172 vector<TH2F*> HotColdMaps;
0173 map<unsigned int, pair<unsigned int, unsigned int> > modCounter[23];
0174 TrackerMap* tkmap;
0175 TrackerMap* tkmapbad;
0176 TrackerMap* tkmapeff;
0177 TrackerMap* tkmapnum;
0178 TrackerMap* tkmapden;
0179 long layerfound[23];
0180 long layertotal[23];
0181 map<unsigned int, vector<int> > layerfound_perBx;
0182 map<unsigned int, vector<int> > layertotal_perBx;
0183 vector<TH1F*> layerfound_vsLumi;
0184 vector<TH1F*> layertotal_vsLumi;
0185 vector<TH1F*> layerfound_vsPU;
0186 vector<TH1F*> layertotal_vsPU;
0187 vector<TH1F*> layerfound_vsCM;
0188 vector<TH1F*> layertotal_vsCM;
0189 vector<TH1F*> layerfound_vsBX;
0190 vector<TH1F*> layertotal_vsBX;
0191 int goodlayertotal[35];
0192 int goodlayerfound[35];
0193 int alllayertotal[35];
0194 int alllayerfound[35];
0195 map<unsigned int, double> BadModules;
0196 };
0197
0198 SiStripHitEffFromCalibTree::SiStripHitEffFromCalibTree(const edm::ParameterSet& conf)
0199 : ConditionDBWriter<SiStripBadStrip>(conf),
0200 fileInPath_(SiStripDetInfoFileReader::kDefaultFile),
0201 tkGeomToken_(esConsumes()),
0202 tTopoToken_(esConsumes()) {
0203 usesResource(TFileService::kSharedResource);
0204 calibTreeFileNames_ = conf.getUntrackedParameter<vector<std::string> >("CalibTreeFilenames");
0205 threshold_ = conf.getParameter<double>("Threshold");
0206 nModsMin_ = conf.getParameter<int>("nModsMin");
0207 badModulesFile_ = conf.getUntrackedParameter<std::string>("BadModulesFile", "");
0208 autoIneffModTagging_ = conf.getUntrackedParameter<bool>("AutoIneffModTagging", false);
0209 clusterMatchingMethod_ = conf.getUntrackedParameter<int>("ClusterMatchingMethod", 0);
0210 resXSig_ = conf.getUntrackedParameter<double>("ResXSig", -1);
0211 clusterTrajDist_ = conf.getUntrackedParameter<double>("ClusterTrajDist", 64.0);
0212 stripsApvEdge_ = conf.getUntrackedParameter<double>("StripsApvEdge", 10.0);
0213 useOnlyHighPurityTracks_ = conf.getUntrackedParameter<bool>("UseOnlyHighPurityTracks", true);
0214 bunchX_ = conf.getUntrackedParameter<int>("BunchCrossing", 0);
0215 spaceBetweenTrains_ = conf.getUntrackedParameter<int>("SpaceBetweenTrains", 25);
0216 useCM_ = conf.getUntrackedParameter<bool>("UseCommonMode", false);
0217 showEndcapSides_ = conf.getUntrackedParameter<bool>("ShowEndcapSides", true);
0218 showRings_ = conf.getUntrackedParameter<bool>("ShowRings", false);
0219 showTOB6TEC9_ = conf.getUntrackedParameter<bool>("ShowTOB6TEC9", false);
0220 showOnlyGoodModules_ = conf.getUntrackedParameter<bool>("ShowOnlyGoodModules", false);
0221 tkMapMin_ = conf.getUntrackedParameter<double>("TkMapMin", 0.9);
0222 effPlotMin_ = conf.getUntrackedParameter<double>("EffPlotMin", 0.9);
0223 title_ = conf.getParameter<std::string>("Title");
0224 storeModEff_ = conf.getUntrackedParameter<bool>("StoreModuleEff", false);
0225 detInfo_ = SiStripDetInfoFileReader::read(fileInPath_.fullPath());
0226
0227 nTEClayers = 9;
0228 if (showRings_)
0229 nTEClayers = 7;
0230
0231 quality_ = new SiStripQuality(detInfo_);
0232 }
0233
0234 void SiStripHitEffFromCalibTree::algoAnalyze(const edm::Event& e, const edm::EventSetup& c) {
0235 const auto& tkgeom = c.getData(tkGeomToken_);
0236 const auto& tTopo = c.getData(tTopoToken_);
0237
0238
0239 ifstream badModules_file;
0240 set<uint32_t> badModules_list;
0241 if (!badModulesFile_.empty()) {
0242 badModules_file.open(badModulesFile_.c_str());
0243 uint32_t badmodule_detid;
0244 int mods, fiber1, fiber2, fiber3;
0245 if (badModules_file.is_open()) {
0246 string line;
0247 while (getline(badModules_file, line)) {
0248 if (badModules_file.eof())
0249 continue;
0250 stringstream ss(line);
0251 ss >> badmodule_detid >> mods >> fiber1 >> fiber2 >> fiber3;
0252 if (badmodule_detid != 0 && mods == 1 && (fiber1 == 1 || fiber2 == 1 || fiber3 == 1))
0253 badModules_list.insert(badmodule_detid);
0254 }
0255 badModules_file.close();
0256 }
0257 }
0258 if (!badModules_list.empty())
0259 LOGPRINT << "Remove additionnal bad modules from the analysis: ";
0260 set<uint32_t>::iterator itBadMod;
0261 for (const auto& badMod : badModules_list) {
0262 LOGPRINT << " " << badMod;
0263 }
0264
0265
0266
0267 bxHisto = fs->make<TH1F>("bx", "bx", 3600, 0, 3600);
0268 instLumiHisto = fs->make<TH1F>("instLumi", "inst. lumi.", 250, 0, 25000);
0269 PUHisto = fs->make<TH1F>("PU", "PU", 300, 0, 300);
0270
0271 bxHisto_cutOnTracks = fs->make<TH1F>("bx_cutOnTracks", "bx", 3600, 0, 3600);
0272 instLumiHisto_cutOnTracks = fs->make<TH1F>("instLumi_cutOnTracks", "inst. lumi.", 250, 0, 25000);
0273 PUHisto_cutOnTracks = fs->make<TH1F>("PU_cutOnTracks", "PU", 300, 0, 300);
0274
0275 for (int l = 0; l < 35; l++) {
0276 goodlayertotal[l] = 0;
0277 goodlayerfound[l] = 0;
0278 alllayertotal[l] = 0;
0279 alllayerfound[l] = 0;
0280 }
0281
0282 TH1F* resolutionPlots[23];
0283 for (Long_t ilayer = 0; ilayer < 23; ilayer++) {
0284 std::string lyrName = ::layerName(ilayer, showRings_, nTEClayers);
0285
0286 resolutionPlots[ilayer] = fs->make<TH1F>(Form("resol_layer_%i", (int)(ilayer)), lyrName.c_str(), 125, -125, 125);
0287 resolutionPlots[ilayer]->GetXaxis()->SetTitle("trajX-clusX [strip unit]");
0288
0289 layerfound_vsLumi.push_back(
0290 fs->make<TH1F>(Form("layerfound_vsLumi_layer_%i", (int)(ilayer)), lyrName.c_str(), 100, 0, 25000));
0291 layertotal_vsLumi.push_back(
0292 fs->make<TH1F>(Form("layertotal_vsLumi_layer_%i", (int)(ilayer)), lyrName.c_str(), 100, 0, 25000));
0293 layerfound_vsPU.push_back(
0294 fs->make<TH1F>(Form("layerfound_vsPU_layer_%i", (int)(ilayer)), lyrName.c_str(), 45, 0, 90));
0295 layertotal_vsPU.push_back(
0296 fs->make<TH1F>(Form("layertotal_vsPU_layer_%i", (int)(ilayer)), lyrName.c_str(), 45, 0, 90));
0297
0298 layerfound_vsBX.push_back(fs->make<TH1F>(
0299 Form("foundVsBx_layer%i", (int)ilayer), Form("layer %i", (int)ilayer), nBxInAnOrbit_, 0, nBxInAnOrbit_));
0300 layertotal_vsBX.push_back(fs->make<TH1F>(
0301 Form("totalVsBx_layer%i", (int)ilayer), Form("layer %i", (int)ilayer), nBxInAnOrbit_, 0, nBxInAnOrbit_));
0302
0303 if (useCM_) {
0304 layerfound_vsCM.push_back(
0305 fs->make<TH1F>(Form("layerfound_vsCM_layer_%i", (int)(ilayer)), lyrName.c_str(), 20, 0, 400));
0306 layertotal_vsCM.push_back(
0307 fs->make<TH1F>(Form("layertotal_vsCM_layer_%i", (int)(ilayer)), lyrName.c_str(), 20, 0, 400));
0308 }
0309 layertotal[ilayer] = 0;
0310 layerfound[ilayer] = 0;
0311 }
0312
0313 if (!autoIneffModTagging_)
0314 LOGPRINT << "A module is bad if efficiency < " << threshold_ << " and has at least " << nModsMin_ << " nModsMin.";
0315 else
0316 LOGPRINT << "A module is bad if the upper limit on the efficiency is < to the avg in the layer - " << threshold_
0317 << " and has at least " << nModsMin_ << " nModsMin.";
0318
0319 unsigned int run, evt, bx{0};
0320 double instLumi, PU;
0321
0322
0323 for (const auto& calibTreeFileName : calibTreeFileNames_) {
0324 LOGPRINT << "Loading file: " << calibTreeFileName;
0325 TFile* CalibTreeFile = TFile::Open(calibTreeFileName.c_str(), "READ");
0326
0327
0328 bool foundEventInfos = false;
0329 try {
0330 CalibTreeFile->cd("eventInfo");
0331 } catch (exception& e) {
0332 LOGPRINT << "No event infos tree";
0333 }
0334 TTree* EventTree = (TTree*)(gDirectory->Get("tree"));
0335
0336 TLeaf* runLf;
0337 TLeaf* evtLf;
0338 TLeaf* BunchLf;
0339 TLeaf* InstLumiLf;
0340 TLeaf* PULf;
0341 if (EventTree) {
0342 LOGPRINT << "Found event infos tree";
0343
0344 runLf = EventTree->GetLeaf("run");
0345 evtLf = EventTree->GetLeaf("event");
0346
0347 BunchLf = EventTree->GetLeaf("bx");
0348 InstLumiLf = EventTree->GetLeaf("instLumi");
0349 PULf = EventTree->GetLeaf("PU");
0350
0351 int nevt = EventTree->GetEntries();
0352 if (nevt)
0353 foundEventInfos = true;
0354
0355 for (int j = 0; j < nevt; j++) {
0356 EventTree->GetEntry(j);
0357 run = runLf->GetValue();
0358 evt = evtLf->GetValue();
0359 bx = BunchLf->GetValue();
0360 instLumi = InstLumiLf->GetValue();
0361 PU = PULf->GetValue();
0362
0363 bxHisto->Fill(bx);
0364 instLumiHisto->Fill(instLumi);
0365 PUHisto->Fill(PU);
0366
0367 eventInfos[make_pair(run, evt)] = array<double, 3>{{(double)bx, instLumi, PU}};
0368 event_used[make_pair(run, evt)] = false;
0369 }
0370 }
0371
0372
0373 CalibTreeFile->cd("anEff");
0374 calibTree_ = (TTree*)(gDirectory->Get("traj"));
0375
0376 runLf = calibTree_->GetLeaf("run");
0377 evtLf = calibTree_->GetLeaf("event");
0378 TLeaf* BadLf = calibTree_->GetLeaf("ModIsBad");
0379 TLeaf* sistripLf = calibTree_->GetLeaf("SiStripQualBad");
0380 TLeaf* idLf = calibTree_->GetLeaf("Id");
0381 TLeaf* acceptLf = calibTree_->GetLeaf("withinAcceptance");
0382 TLeaf* layerLf = calibTree_->GetLeaf("layer");
0383
0384 TLeaf* highPurityLf = calibTree_->GetLeaf("highPurity");
0385 TLeaf* xLf = calibTree_->GetLeaf("TrajGlbX");
0386 TLeaf* yLf = calibTree_->GetLeaf("TrajGlbY");
0387 TLeaf* zLf = calibTree_->GetLeaf("TrajGlbZ");
0388 TLeaf* ResXSigLf = calibTree_->GetLeaf("ResXSig");
0389 TLeaf* TrajLocXLf = calibTree_->GetLeaf("TrajLocX");
0390 TLeaf* TrajLocYLf = calibTree_->GetLeaf("TrajLocY");
0391 TLeaf* ClusterLocXLf = calibTree_->GetLeaf("ClusterLocX");
0392 BunchLf = calibTree_->GetLeaf("bunchx");
0393 InstLumiLf = calibTree_->GetLeaf("instLumi");
0394 PULf = calibTree_->GetLeaf("PU");
0395 TLeaf* CMLf = nullptr;
0396 if (useCM_)
0397 CMLf = calibTree_->GetLeaf("commonMode");
0398
0399 int nevents = calibTree_->GetEntries();
0400 LOGPRINT << "Successfully loaded analyze function with " << nevents << " events!\n";
0401
0402 map<pair<unsigned int, unsigned int>, array<double, 3> >::iterator itEventInfos;
0403 map<pair<unsigned int, unsigned int>, bool>::iterator itEventUsed;
0404
0405
0406 for (int j = 0; j < nevents; j++) {
0407 calibTree_->GetEntry(j);
0408 run = (unsigned int)runLf->GetValue();
0409 evt = (unsigned int)evtLf->GetValue();
0410 unsigned int isBad = (unsigned int)BadLf->GetValue();
0411 unsigned int quality = (unsigned int)sistripLf->GetValue();
0412 unsigned int id = (unsigned int)idLf->GetValue();
0413 unsigned int accept = (unsigned int)acceptLf->GetValue();
0414 unsigned int layer_wheel = (unsigned int)layerLf->GetValue();
0415 unsigned int layer = layer_wheel;
0416 if (showRings_ && layer > 10) {
0417 if (layer < 14)
0418 layer = 10 + ((id >> 9) & 0x3);
0419 else
0420 layer = 13 + ((id >> 5) & 0x7);
0421 }
0422
0423 bool highPurity = (bool)highPurityLf->GetValue();
0424 double x = xLf->GetValue();
0425 double y = yLf->GetValue();
0426 double z = zLf->GetValue();
0427 double resxsig = ResXSigLf->GetValue();
0428 double TrajLocX = TrajLocXLf->GetValue();
0429 double TrajLocY = TrajLocYLf->GetValue();
0430 double ClusterLocX = ClusterLocXLf->GetValue();
0431 double TrajLocXMid;
0432 double stripTrajMid;
0433 double stripCluster;
0434 bool badquality = false;
0435
0436 instLumi = 0;
0437 PU = 0;
0438
0439
0440 if (!foundEventInfos) {
0441 bx = (unsigned int)BunchLf->GetValue();
0442 if (InstLumiLf != nullptr)
0443 instLumi = InstLumiLf->GetValue();
0444 if (PULf != nullptr)
0445 PU = PULf->GetValue();
0446
0447 itEventUsed = event_used.find(make_pair(run, evt));
0448 if (itEventUsed == event_used.end())
0449 event_used[make_pair(run, evt)] = false;
0450 }
0451 int CM = -100;
0452 if (useCM_)
0453 CM = CMLf->GetValue();
0454
0455
0456 if (foundEventInfos) {
0457 itEventInfos = eventInfos.find(make_pair(run, evt));
0458 if (itEventInfos != eventInfos.end()) {
0459 bx = itEventInfos->second[0];
0460 instLumi = itEventInfos->second[1];
0461 PU = itEventInfos->second[2];
0462 }
0463 }
0464
0465
0466
0467 itEventUsed = event_used.find(make_pair(run, evt));
0468 if (itEventUsed != event_used.end()) {
0469 if (itEventUsed->second == false) {
0470 bxHisto_cutOnTracks->Fill(bx);
0471 instLumiHisto_cutOnTracks->Fill(instLumi);
0472 PUHisto_cutOnTracks->Fill(PU);
0473 itEventUsed->second = true;
0474 }
0475 }
0476
0477
0478
0479
0480 if (bunchX_ > 0 && bunchX_ != bx)
0481 continue;
0482
0483
0484 if (accept != 1)
0485 continue;
0486 if (useOnlyHighPurityTracks_ && !highPurity)
0487 continue;
0488 if (quality == 1)
0489 badquality = true;
0490
0491
0492 if (!showTOB6TEC9_ && (layer_wheel == 10 || layer_wheel == siStripLayers_))
0493 continue;
0494
0495
0496 itBadMod = badModules_list.find(id);
0497 if (itBadMod != badModules_list.end())
0498 continue;
0499
0500
0501
0502
0503 bool badflag = false;
0504
0505
0506 if (resXSig_ < 0) {
0507 if (isBad == 1)
0508 badflag = true;
0509 } else {
0510 if (isBad == 1 || resxsig > resXSig_)
0511 badflag = true;
0512 }
0513
0514
0515 int nstrips = -9;
0516 float Pitch = -9.0;
0517
0518 if (resxsig == 1000.0) {
0519 Pitch = 0.0205;
0520 nstrips = 768;
0521 stripTrajMid = TrajLocX / Pitch + nstrips / 2.0;
0522 stripCluster = ClusterLocX / Pitch + nstrips / 2.0;
0523 } else {
0524 DetId ClusterDetId(id);
0525 const StripGeomDetUnit* stripdet = (const StripGeomDetUnit*)tkgeom.idToDetUnit(ClusterDetId);
0526 const StripTopology& Topo = stripdet->specificTopology();
0527 nstrips = Topo.nstrips();
0528 Pitch = stripdet->surface().bounds().width() / Topo.nstrips();
0529 stripTrajMid = TrajLocX / Pitch + nstrips / 2.0;
0530 stripCluster = ClusterLocX / Pitch + nstrips / 2.0;
0531
0532
0533
0534 float hbedge = 0;
0535 float htedge = 0;
0536 float hapoth = 0;
0537 if (layer >= 11) {
0538 const BoundPlane& plane = stripdet->surface();
0539 const TrapezoidalPlaneBounds* trapezoidalBounds(
0540 dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
0541 std::array<const float, 4> const& parameters = (*trapezoidalBounds).parameters();
0542 hbedge = parameters[0];
0543 htedge = parameters[1];
0544 hapoth = parameters[3];
0545 TrajLocXMid = TrajLocX / (1 + (htedge - hbedge) * TrajLocY / (htedge + hbedge) /
0546 hapoth);
0547 stripTrajMid = TrajLocXMid / Pitch + nstrips / 2.0;
0548 }
0549 }
0550
0551 if (!badquality && layer < 23) {
0552 if (resxsig != 1000.0)
0553 resolutionPlots[layer]->Fill(stripTrajMid - stripCluster);
0554 else
0555 resolutionPlots[layer]->Fill(1000);
0556 }
0557
0558
0559 int tapv = -9;
0560 int capv = -9;
0561 float stripInAPV = 64.;
0562
0563 if (clusterMatchingMethod_ >= 1) {
0564 badflag = false;
0565 if (resxsig == 1000.0) {
0566 badflag = true;
0567 } else {
0568 if (clusterMatchingMethod_ == 2 ||
0569 clusterMatchingMethod_ == 4) {
0570 if (abs(stripCluster - stripTrajMid) > clusterTrajDist_)
0571 badflag = true;
0572 }
0573 if (clusterMatchingMethod_ == 3 ||
0574 clusterMatchingMethod_ ==
0575 4) {
0576 tapv = (int)stripTrajMid / 128;
0577 capv = (int)stripCluster / 128;
0578 stripInAPV = stripTrajMid - tapv * 128;
0579
0580 if (stripInAPV < stripsApvEdge_ || stripInAPV > 128 - stripsApvEdge_)
0581 continue;
0582 if (tapv != capv)
0583 badflag = true;
0584 }
0585 }
0586 }
0587
0588 if (badflag && !badquality) {
0589 hit temphit;
0590 temphit.x = x;
0591 temphit.y = y;
0592 temphit.z = z;
0593 temphit.id = id;
0594 hits[layer].push_back(temphit);
0595 }
0596 pair<unsigned int, unsigned int> newgoodpair(1, 1);
0597 pair<unsigned int, unsigned int> newbadpair(1, 0);
0598
0599 map<unsigned int, pair<unsigned int, unsigned int> >::iterator it = modCounter[layer].find(id);
0600 if (!badquality) {
0601 if (it == modCounter[layer].end()) {
0602 if (badflag)
0603 modCounter[layer][id] = newbadpair;
0604 else
0605 modCounter[layer][id] = newgoodpair;
0606 } else {
0607 ((*it).second.first)++;
0608 if (!badflag)
0609 ((*it).second.second)++;
0610 }
0611
0612 if (layerfound_perBx.find(bx) == layerfound_perBx.end()) {
0613 layerfound_perBx[bx] = vector<int>(23, 0);
0614 layertotal_perBx[bx] = vector<int>(23, 0);
0615 }
0616 if (!badflag)
0617 layerfound_perBx[bx][layer]++;
0618 layertotal_perBx[bx][layer]++;
0619
0620 if (!badflag)
0621 layerfound_vsLumi[layer]->Fill(instLumi);
0622 layertotal_vsLumi[layer]->Fill(instLumi);
0623 if (!badflag)
0624 layerfound_vsPU[layer]->Fill(PU);
0625 layertotal_vsPU[layer]->Fill(PU);
0626
0627 if (useCM_) {
0628 if (!badflag)
0629 layerfound_vsCM[layer]->Fill(CM);
0630 layertotal_vsCM[layer]->Fill(CM);
0631 }
0632
0633
0634 if (layer <= 10) {
0635 if (!badflag)
0636 goodlayerfound[layer]++;
0637 goodlayertotal[layer]++;
0638 } else if (layer > 10 && layer < 14) {
0639 if (((id >> 13) & 0x3) == 1) {
0640 if (!badflag)
0641 goodlayerfound[layer]++;
0642 goodlayertotal[layer]++;
0643 } else if (((id >> 13) & 0x3) == 2) {
0644 if (!badflag)
0645 goodlayerfound[layer + 3]++;
0646 goodlayertotal[layer + 3]++;
0647 }
0648 } else if (layer > 13 && layer <= siStripLayers_) {
0649 if (((id >> 18) & 0x3) == 1) {
0650 if (!badflag)
0651 goodlayerfound[layer + 3]++;
0652 goodlayertotal[layer + 3]++;
0653 } else if (((id >> 18) & 0x3) == 2) {
0654 if (!badflag)
0655 goodlayerfound[layer + 3 + nTEClayers]++;
0656 goodlayertotal[layer + 3 + nTEClayers]++;
0657 }
0658 }
0659 }
0660
0661 if (layer <= 10) {
0662 if (!badflag)
0663 alllayerfound[layer]++;
0664 alllayertotal[layer]++;
0665 } else if (layer > 10 && layer < 14) {
0666 if (((id >> 13) & 0x3) == 1) {
0667 if (!badflag)
0668 alllayerfound[layer]++;
0669 alllayertotal[layer]++;
0670 } else if (((id >> 13) & 0x3) == 2) {
0671 if (!badflag)
0672 alllayerfound[layer + 3]++;
0673 alllayertotal[layer + 3]++;
0674 }
0675 } else if (layer > 13 && layer <= siStripLayers_) {
0676 if (((id >> 18) & 0x3) == 1) {
0677 if (!badflag)
0678 alllayerfound[layer + 3]++;
0679 alllayertotal[layer + 3]++;
0680 } else if (((id >> 18) & 0x3) == 2) {
0681 if (!badflag)
0682 alllayerfound[layer + 3 + nTEClayers]++;
0683 alllayertotal[layer + 3 + nTEClayers]++;
0684 }
0685 }
0686
0687 }
0688 }
0689
0690 makeHotColdMaps();
0691 makeTKMap(autoIneffModTagging_);
0692 makeSQLite();
0693 totalStatistics();
0694 makeSummary();
0695 makeSummaryVsBx();
0696 makeSummaryVsLumi();
0697 if (useCM_)
0698 makeSummaryVsCM();
0699
0700
0701
0702
0703 int NTkBadComponent[4];
0704 int NBadComponent[4][19][4];
0705
0706
0707
0708 std::stringstream ssV[4][19];
0709
0710 for (int i = 0; i < 4; ++i) {
0711 NTkBadComponent[i] = 0;
0712 for (int j = 0; j < 19; ++j) {
0713 ssV[i][j].str("");
0714 for (int k = 0; k < 4; ++k)
0715 NBadComponent[i][j][k] = 0;
0716 }
0717 }
0718
0719 std::vector<SiStripQuality::BadComponent> BC = quality_->getBadComponentList();
0720
0721 for (size_t i = 0; i < BC.size(); ++i) {
0722
0723
0724
0725
0726 if (BC[i].BadModule)
0727 NTkBadComponent[0]++;
0728 if (BC[i].BadFibers)
0729 NTkBadComponent[1] += ((BC[i].BadFibers >> 2) & 0x1) + ((BC[i].BadFibers >> 1) & 0x1) + ((BC[i].BadFibers) & 0x1);
0730 if (BC[i].BadApvs)
0731 NTkBadComponent[2] += ((BC[i].BadApvs >> 5) & 0x1) + ((BC[i].BadApvs >> 4) & 0x1) + ((BC[i].BadApvs >> 3) & 0x1) +
0732 ((BC[i].BadApvs >> 2) & 0x1) + ((BC[i].BadApvs >> 1) & 0x1) + ((BC[i].BadApvs) & 0x1);
0733
0734
0735
0736
0737
0738 int component;
0739 SiStripDetId a(BC[i].detid);
0740 if (a.subdetId() == SiStripDetId::TIB) {
0741
0742
0743
0744
0745 component = tTopo.tibLayer(BC[i].detid);
0746 setBadComponents(0, component, BC[i], ssV, NBadComponent);
0747
0748 } else if (a.subdetId() == SiStripDetId::TID) {
0749
0750
0751
0752
0753 component = tTopo.tidSide(BC[i].detid) == 2 ? tTopo.tidWheel(BC[i].detid) : tTopo.tidWheel(BC[i].detid) + 3;
0754 setBadComponents(1, component, BC[i], ssV, NBadComponent);
0755
0756 } else if (a.subdetId() == SiStripDetId::TOB) {
0757
0758
0759
0760
0761 component = tTopo.tobLayer(BC[i].detid);
0762 setBadComponents(2, component, BC[i], ssV, NBadComponent);
0763
0764 } else if (a.subdetId() == SiStripDetId::TEC) {
0765
0766
0767
0768
0769 component = tTopo.tecSide(BC[i].detid) == 2 ? tTopo.tecWheel(BC[i].detid) : tTopo.tecWheel(BC[i].detid) + 9;
0770 setBadComponents(3, component, BC[i], ssV, NBadComponent);
0771 }
0772 }
0773
0774
0775
0776
0777 float percentage = 0;
0778
0779 SiStripQuality::RegistryIterator rbegin = quality_->getRegistryVectorBegin();
0780 SiStripQuality::RegistryIterator rend = quality_->getRegistryVectorEnd();
0781
0782 for (SiStripBadStrip::RegistryIterator rp = rbegin; rp != rend; ++rp) {
0783 unsigned int detid = rp->detid;
0784
0785 int subdet = -999;
0786 int component = -999;
0787 SiStripDetId a(detid);
0788 if (a.subdetId() == 3) {
0789 subdet = 0;
0790 component = tTopo.tibLayer(detid);
0791 } else if (a.subdetId() == 4) {
0792 subdet = 1;
0793 component = tTopo.tidSide(detid) == 2 ? tTopo.tidWheel(detid) : tTopo.tidWheel(detid) + 3;
0794 } else if (a.subdetId() == 5) {
0795 subdet = 2;
0796 component = tTopo.tobLayer(detid);
0797 } else if (a.subdetId() == 6) {
0798 subdet = 3;
0799 component = tTopo.tecSide(detid) == 2 ? tTopo.tecWheel(detid) : tTopo.tecWheel(detid) + 9;
0800 }
0801
0802 SiStripQuality::Range sqrange =
0803 SiStripQuality::Range(quality_->getDataVectorBegin() + rp->ibegin, quality_->getDataVectorBegin() + rp->iend);
0804
0805 percentage = 0;
0806 for (int it = 0; it < sqrange.second - sqrange.first; it++) {
0807 unsigned int range = quality_->decode(*(sqrange.first + it)).range;
0808 NTkBadComponent[3] += range;
0809 NBadComponent[subdet][0][3] += range;
0810 NBadComponent[subdet][component][3] += range;
0811 percentage += range;
0812 }
0813 if (percentage != 0)
0814 percentage /= 128. * detInfo_.getNumberOfApvsAndStripLength(detid).first;
0815 if (percentage > 1)
0816 edm::LogError("SiStripQualityStatistics") << "PROBLEM detid " << detid << " value " << percentage << std::endl;
0817 }
0818
0819
0820
0821 std::ostringstream ss;
0822
0823 ss << "\n-----------------\nNew IOV starting from run " << e.id().run() << " event " << e.id().event()
0824 << " lumiBlock " << e.luminosityBlock() << " time " << e.time().value() << "\n-----------------\n";
0825 ss << "\n-----------------\nGlobal Info\n-----------------";
0826 ss << "\nBadComponent \t Modules \tFibers "
0827 "\tApvs\tStrips\n----------------------------------------------------------------";
0828 ss << "\nTracker:\t\t" << NTkBadComponent[0] << "\t" << NTkBadComponent[1] << "\t" << NTkBadComponent[2] << "\t"
0829 << NTkBadComponent[3];
0830 ss << "\nTIB:\t\t\t" << NBadComponent[0][0][0] << "\t" << NBadComponent[0][0][1] << "\t" << NBadComponent[0][0][2]
0831 << "\t" << NBadComponent[0][0][3];
0832 ss << "\nTID:\t\t\t" << NBadComponent[1][0][0] << "\t" << NBadComponent[1][0][1] << "\t" << NBadComponent[1][0][2]
0833 << "\t" << NBadComponent[1][0][3];
0834 ss << "\nTOB:\t\t\t" << NBadComponent[2][0][0] << "\t" << NBadComponent[2][0][1] << "\t" << NBadComponent[2][0][2]
0835 << "\t" << NBadComponent[2][0][3];
0836 ss << "\nTEC:\t\t\t" << NBadComponent[3][0][0] << "\t" << NBadComponent[3][0][1] << "\t" << NBadComponent[3][0][2]
0837 << "\t" << NBadComponent[3][0][3];
0838 ss << "\n";
0839
0840 for (int i = 1; i < 5; ++i)
0841 ss << "\nTIB Layer " << i << " :\t\t" << NBadComponent[0][i][0] << "\t" << NBadComponent[0][i][1] << "\t"
0842 << NBadComponent[0][i][2] << "\t" << NBadComponent[0][i][3];
0843 ss << "\n";
0844 for (int i = 1; i < 4; ++i)
0845 ss << "\nTID+ Disk " << i << " :\t\t" << NBadComponent[1][i][0] << "\t" << NBadComponent[1][i][1] << "\t"
0846 << NBadComponent[1][i][2] << "\t" << NBadComponent[1][i][3];
0847 for (int i = 4; i < 7; ++i)
0848 ss << "\nTID- Disk " << i - 3 << " :\t\t" << NBadComponent[1][i][0] << "\t" << NBadComponent[1][i][1] << "\t"
0849 << NBadComponent[1][i][2] << "\t" << NBadComponent[1][i][3];
0850 ss << "\n";
0851 for (int i = 1; i < 7; ++i)
0852 ss << "\nTOB Layer " << i << " :\t\t" << NBadComponent[2][i][0] << "\t" << NBadComponent[2][i][1] << "\t"
0853 << NBadComponent[2][i][2] << "\t" << NBadComponent[2][i][3];
0854 ss << "\n";
0855 for (int i = 1; i < 10; ++i)
0856 ss << "\nTEC+ Disk " << i << " :\t\t" << NBadComponent[3][i][0] << "\t" << NBadComponent[3][i][1] << "\t"
0857 << NBadComponent[3][i][2] << "\t" << NBadComponent[3][i][3];
0858 for (int i = 10; i < 19; ++i)
0859 ss << "\nTEC- Disk " << i - 9 << " :\t\t" << NBadComponent[3][i][0] << "\t" << NBadComponent[3][i][1] << "\t"
0860 << NBadComponent[3][i][2] << "\t" << NBadComponent[3][i][3];
0861 ss << "\n";
0862
0863 ss << "\n----------------------------------------------------------------\n\t\t Detid \tModules Fibers "
0864 "Apvs\n----------------------------------------------------------------";
0865 for (int i = 1; i < 5; ++i)
0866 ss << "\nTIB Layer " << i << " :" << ssV[0][i].str();
0867 ss << "\n";
0868 for (int i = 1; i < 4; ++i)
0869 ss << "\nTID+ Disk " << i << " :" << ssV[1][i].str();
0870 for (int i = 4; i < 7; ++i)
0871 ss << "\nTID- Disk " << i - 3 << " :" << ssV[1][i].str();
0872 ss << "\n";
0873 for (int i = 1; i < 7; ++i)
0874 ss << "\nTOB Layer " << i << " :" << ssV[2][i].str();
0875 ss << "\n";
0876 for (int i = 1; i < 10; ++i)
0877 ss << "\nTEC+ Disk " << i << " :" << ssV[3][i].str();
0878 for (int i = 10; i < 19; ++i)
0879 ss << "\nTEC- Disk " << i - 9 << " :" << ssV[3][i].str();
0880
0881 LOGPRINT << ss.str();
0882
0883
0884 ofstream badModules;
0885 badModules.open("BadModules.log");
0886 badModules << "\n----------------------------------------------------------------\n\t\t Detid \tModules Fibers "
0887 "Apvs\n----------------------------------------------------------------";
0888 for (int i = 1; i < 5; ++i)
0889 badModules << "\nTIB Layer " << i << " :" << ssV[0][i].str();
0890 badModules << "\n";
0891 for (int i = 1; i < 4; ++i)
0892 badModules << "\nTID+ Disk " << i << " :" << ssV[1][i].str();
0893 for (int i = 4; i < 7; ++i)
0894 badModules << "\nTID- Disk " << i - 3 << " :" << ssV[1][i].str();
0895 badModules << "\n";
0896 for (int i = 1; i < 7; ++i)
0897 badModules << "\nTOB Layer " << i << " :" << ssV[2][i].str();
0898 badModules << "\n";
0899 for (int i = 1; i < 10; ++i)
0900 badModules << "\nTEC+ Disk " << i << " :" << ssV[3][i].str();
0901 for (int i = 10; i < 19; ++i)
0902 badModules << "\nTEC- Disk " << i - 9 << " :" << ssV[3][i].str();
0903 badModules.close();
0904 }
0905
0906 void SiStripHitEffFromCalibTree::makeHotColdMaps() {
0907 LOGPRINT << "Entering hot cold map generation!\n";
0908 TStyle* gStyle = new TStyle("gStyle", "myStyle");
0909 gStyle->cd();
0910 gStyle->SetPalette(1);
0911 gStyle->SetCanvasColor(kWhite);
0912 gStyle->SetOptStat(0);
0913
0914
0915
0916 TH2F* temph2;
0917 for (Long_t maplayer = 1; maplayer <= siStripLayers_; maplayer++) {
0918
0919 if (maplayer > 0 && maplayer <= 4) {
0920
0921 temph2 = fs->make<TH2F>(Form("%s%i", "TIB", (int)(maplayer)), "TIB", 100, -1, 361, 100, -100, 100);
0922 temph2->GetXaxis()->SetTitle("Phi");
0923 temph2->GetXaxis()->SetBinLabel(1, TString("360"));
0924 temph2->GetXaxis()->SetBinLabel(50, TString("180"));
0925 temph2->GetXaxis()->SetBinLabel(100, TString("0"));
0926 temph2->GetYaxis()->SetTitle("Global Z");
0927 temph2->SetOption("colz");
0928 HotColdMaps.push_back(temph2);
0929 } else if (maplayer > 4 && maplayer <= 10) {
0930
0931 temph2 = fs->make<TH2F>(Form("%s%i", "TOB", (int)(maplayer - 4)), "TOB", 100, -1, 361, 100, -120, 120);
0932 temph2->GetXaxis()->SetTitle("Phi");
0933 temph2->GetXaxis()->SetBinLabel(1, TString("360"));
0934 temph2->GetXaxis()->SetBinLabel(50, TString("180"));
0935 temph2->GetXaxis()->SetBinLabel(100, TString("0"));
0936 temph2->GetYaxis()->SetTitle("Global Z");
0937 temph2->SetOption("colz");
0938 HotColdMaps.push_back(temph2);
0939 } else if (maplayer > 10 && maplayer <= 13) {
0940
0941
0942 temph2 = fs->make<TH2F>(Form("%s%i", "TID-", (int)(maplayer - 10)), "TID-", 100, -100, 100, 100, -100, 100);
0943 temph2->GetXaxis()->SetTitle("Global Y");
0944 temph2->GetXaxis()->SetBinLabel(1, TString("+Y"));
0945 temph2->GetXaxis()->SetBinLabel(50, TString("0"));
0946 temph2->GetXaxis()->SetBinLabel(100, TString("-Y"));
0947 temph2->GetYaxis()->SetTitle("Global X");
0948 temph2->GetYaxis()->SetBinLabel(1, TString("-X"));
0949 temph2->GetYaxis()->SetBinLabel(50, TString("0"));
0950 temph2->GetYaxis()->SetBinLabel(100, TString("+X"));
0951 temph2->SetOption("colz");
0952 HotColdMaps.push_back(temph2);
0953 temph2 = fs->make<TH2F>(Form("%s%i", "TID+", (int)(maplayer - 10)), "TID+", 100, -100, 100, 100, -100, 100);
0954 temph2->GetXaxis()->SetTitle("Global Y");
0955 temph2->GetXaxis()->SetBinLabel(1, TString("+Y"));
0956 temph2->GetXaxis()->SetBinLabel(50, TString("0"));
0957 temph2->GetXaxis()->SetBinLabel(100, TString("-Y"));
0958 temph2->GetYaxis()->SetTitle("Global X");
0959 temph2->GetYaxis()->SetBinLabel(1, TString("-X"));
0960 temph2->GetYaxis()->SetBinLabel(50, TString("0"));
0961 temph2->GetYaxis()->SetBinLabel(100, TString("+X"));
0962 temph2->SetOption("colz");
0963 HotColdMaps.push_back(temph2);
0964 } else if (maplayer > 13) {
0965
0966
0967 temph2 = fs->make<TH2F>(Form("%s%i", "TEC-", (int)(maplayer - 13)), "TEC-", 100, -120, 120, 100, -120, 120);
0968 temph2->GetXaxis()->SetTitle("Global Y");
0969 temph2->GetXaxis()->SetBinLabel(1, TString("+Y"));
0970 temph2->GetXaxis()->SetBinLabel(50, TString("0"));
0971 temph2->GetXaxis()->SetBinLabel(100, TString("-Y"));
0972 temph2->GetYaxis()->SetTitle("Global X");
0973 temph2->GetYaxis()->SetBinLabel(1, TString("-X"));
0974 temph2->GetYaxis()->SetBinLabel(50, TString("0"));
0975 temph2->GetYaxis()->SetBinLabel(100, TString("+X"));
0976 temph2->SetOption("colz");
0977 HotColdMaps.push_back(temph2);
0978 temph2 = fs->make<TH2F>(Form("%s%i", "TEC+", (int)(maplayer - 13)), "TEC+", 100, -120, 120, 100, -120, 120);
0979 temph2->GetXaxis()->SetTitle("Global Y");
0980 temph2->GetXaxis()->SetBinLabel(1, TString("+Y"));
0981 temph2->GetXaxis()->SetBinLabel(50, TString("0"));
0982 temph2->GetXaxis()->SetBinLabel(100, TString("-Y"));
0983 temph2->GetYaxis()->SetTitle("Global X");
0984 temph2->GetYaxis()->SetBinLabel(1, TString("-X"));
0985 temph2->GetYaxis()->SetBinLabel(50, TString("0"));
0986 temph2->GetYaxis()->SetBinLabel(100, TString("+X"));
0987 temph2->SetOption("colz");
0988 HotColdMaps.push_back(temph2);
0989 }
0990 }
0991 for (Long_t mylayer = 1; mylayer <= siStripLayers_; mylayer++) {
0992
0993
0994
0995 vector<hit>::const_iterator iter;
0996 for (iter = hits[mylayer].begin(); iter != hits[mylayer].end(); iter++) {
0997
0998
0999
1000 if (mylayer > 0 && mylayer <= 4) {
1001
1002 float phi = ::calcPhi(iter->x, iter->y);
1003 HotColdMaps[mylayer - 1]->Fill(360. - phi, iter->z, 1.);
1004 } else if (mylayer > 4 && mylayer <= 10) {
1005
1006 float phi = ::calcPhi(iter->x, iter->y);
1007 HotColdMaps[mylayer - 1]->Fill(360. - phi, iter->z, 1.);
1008 } else if (mylayer > 10 && mylayer <= 13) {
1009
1010
1011 int side = (((iter->id) >> 13) & 0x3);
1012 if (side == 1)
1013 HotColdMaps[(mylayer - 1) + (mylayer - 11)]->Fill(-iter->y, iter->x, 1.);
1014 else if (side == 2)
1015 HotColdMaps[(mylayer - 1) + (mylayer - 10)]->Fill(-iter->y, iter->x, 1.);
1016
1017
1018 } else if (mylayer > 13) {
1019
1020
1021 int side = (((iter->id) >> 18) & 0x3);
1022 if (side == 1)
1023 HotColdMaps[(mylayer + 2) + (mylayer - 14)]->Fill(-iter->y, iter->x, 1.);
1024 else if (side == 2)
1025 HotColdMaps[(mylayer + 2) + (mylayer - 13)]->Fill(-iter->y, iter->x, 1.);
1026
1027
1028 }
1029 }
1030 }
1031 LOGPRINT << "Finished HotCold Map Generation\n";
1032 }
1033
1034 void SiStripHitEffFromCalibTree::makeTKMap(bool autoTagging = false) {
1035 LOGPRINT << "Entering TKMap generation!\n";
1036
1037 TTree* tree = nullptr;
1038 unsigned int t_DetId, t_found, t_total;
1039 unsigned char t_layer;
1040 bool t_isTaggedIneff;
1041 float t_threshold;
1042 if (storeModEff_) {
1043 tree = fs->make<TTree>("ModEff", "ModEff");
1044 tree->Branch("DetId", &t_DetId, "DetId/i");
1045 tree->Branch("Layer", &t_layer, "Layer/b");
1046 tree->Branch("FoundHits", &t_found, "FoundHits/i");
1047 tree->Branch("AllHits", &t_total, "AllHits/i");
1048 tree->Branch("IsTaggedIneff", &t_isTaggedIneff, "IsTaggedIneff/O");
1049 tree->Branch("TagThreshold", &t_threshold, "TagThreshold/F");
1050 }
1051
1052 tkmap = new TrackerMap(" Detector Inefficiency ");
1053 tkmapbad = new TrackerMap(" Inefficient Modules ");
1054 tkmapeff = new TrackerMap(title_.Data());
1055 tkmapnum = new TrackerMap(" Detector numerator ");
1056 tkmapden = new TrackerMap(" Detector denominator ");
1057
1058 double myeff, mynum, myden, myeff_up;
1059 double layer_min_eff = 0;
1060
1061 for (Long_t i = 1; i <= siStripLayers_; i++) {
1062
1063
1064 layertotal[i] = 0;
1065 layerfound[i] = 0;
1066 TH1F* hEffInLayer =
1067 fs->make<TH1F>(Form("eff_layer%i", int(i)), Form("Module efficiency in layer %i", int(i)), 201, 0, 1.005);
1068
1069 for (const auto& ih : modCounter[i]) {
1070
1071
1072 mynum = (double)((ih.second).second);
1073 myden = (double)((ih.second).first);
1074 if (myden > 0)
1075 myeff = mynum / myden;
1076 else
1077 myeff = 0;
1078 hEffInLayer->Fill(myeff);
1079
1080 if (!autoTagging) {
1081 if ((myden >= nModsMin_) && (myeff < threshold_)) {
1082
1083 BadModules[ih.first] = myeff;
1084 tkmapbad->fillc(ih.first, 255, 0, 0);
1085 LOGPRINT << "Layer " << i << " (" << ::layerName(i, showRings_, nTEClayers) << ") module " << ih.first
1086 << " efficiency: " << myeff << " , " << mynum << "/" << myden;
1087 } else {
1088
1089 tkmapbad->fillc(ih.first, 255, 255, 255);
1090 }
1091 if (myeff < threshold_)
1092 LOGPRINT << "Layer " << i << " (" << ::layerName(i, showRings_, nTEClayers) << ") module " << ih.first
1093 << " efficiency: " << myeff << " , " << mynum << "/" << myden;
1094 if (myden < nModsMin_) {
1095 LOGPRINT << "Layer " << i << " (" << ::layerName(i, showRings_, nTEClayers) << ") module " << ih.first
1096 << " is under occupancy at " << myden;
1097 }
1098 }
1099
1100 if (storeModEff_) {
1101 t_DetId = ih.first;
1102 t_layer = i;
1103 t_found = mynum;
1104 t_total = myden;
1105 t_isTaggedIneff = false;
1106 t_threshold = 0;
1107 tree->Fill();
1108 }
1109
1110
1111 tkmap->fill(ih.first, 1. - myeff);
1112 tkmapeff->fill(ih.first, myeff);
1113 tkmapnum->fill(ih.first, mynum);
1114 tkmapden->fill(ih.first, myden);
1115
1116
1117 layertotal[i] += long(myden);
1118 layerfound[i] += long(mynum);
1119 }
1120
1121 if (autoTagging) {
1122
1123 hEffInLayer->GetXaxis()->SetRange(3, hEffInLayer->GetNbinsX() + 1);
1124 layer_min_eff =
1125 hEffInLayer->GetMean() - 2.5 * hEffInLayer->GetRMS();
1126 if (threshold_ > 2.5 * hEffInLayer->GetRMS())
1127 layer_min_eff = hEffInLayer->GetMean() - threshold_;
1128 LOGPRINT << "Layer " << i << " threshold for bad modules: <" << layer_min_eff
1129 << " (layer mean: " << hEffInLayer->GetMean() << " rms: " << hEffInLayer->GetRMS() << ")";
1130
1131 hEffInLayer->GetXaxis()->SetRange(1, hEffInLayer->GetNbinsX() + 1);
1132
1133 for (const auto& ih : modCounter[i]) {
1134
1135 mynum = (double)((ih.second).second);
1136 myden = (double)((ih.second).first);
1137 if (myden > 0)
1138 myeff = mynum / myden;
1139 else
1140 myeff = 0;
1141
1142 myeff_up = TEfficiency::Bayesian(myden, mynum, .99, 1, 1, true);
1143 if ((myden >= nModsMin_) && (myeff_up < layer_min_eff)) {
1144
1145 BadModules[ih.first] = myeff;
1146 tkmapbad->fillc(ih.first, 255, 0, 0);
1147 t_isTaggedIneff = true;
1148 } else {
1149
1150 tkmapbad->fillc(ih.first, 255, 255, 255);
1151 t_isTaggedIneff = false;
1152 }
1153 if (myeff_up < layer_min_eff + 0.08)
1154 LOGPRINT << "Layer " << i << " (" << ::layerName(i, showRings_, nTEClayers) << ") module " << ih.first
1155 << " efficiency: " << myeff << " , " << mynum << "/" << myden << " , upper limit: " << myeff_up;
1156 if (myden < nModsMin_) {
1157 LOGPRINT << "Layer " << i << " (" << ::layerName(i, showRings_, nTEClayers) << ") module " << ih.first
1158 << " layer " << i << " is under occupancy at " << myden;
1159 }
1160
1161 if (storeModEff_) {
1162 t_DetId = ih.first;
1163 t_layer = i;
1164 t_found = mynum;
1165 t_total = myden;
1166 t_threshold = layer_min_eff;
1167 tree->Fill();
1168 }
1169 }
1170 }
1171 }
1172 tkmap->save(true, 0, 0, "SiStripHitEffTKMap.png");
1173 tkmapbad->save(true, 0, 0, "SiStripHitEffTKMapBad.png");
1174 tkmapeff->save(true, tkMapMin_, 1., "SiStripHitEffTKMapEff.png");
1175 tkmapnum->save(true, 0, 0, "SiStripHitEffTKMapNum.png");
1176 tkmapden->save(true, 0, 0, "SiStripHitEffTKMapDen.png");
1177 LOGPRINT << "Finished TKMap Generation\n";
1178 if (storeModEff_)
1179 LOGPRINT << "Modules efficiencies stored in tree\n";
1180 }
1181
1182 void SiStripHitEffFromCalibTree::makeSQLite() {
1183
1184 LOGPRINT << "Entering SQLite file generation!\n";
1185 std::vector<unsigned int> BadStripList;
1186 unsigned short NStrips;
1187 unsigned int id1;
1188 std::unique_ptr<SiStripQuality> pQuality = std::make_unique<SiStripQuality>(detInfo_);
1189
1190
1191
1192 for (const auto& it : BadModules) {
1193
1194
1195 NStrips = detInfo_.getNumberOfApvsAndStripLength(it.first).first * 128;
1196 LOGPRINT << "Number of strips module " << it.first << " is " << NStrips;
1197 BadStripList.push_back(pQuality->encode(0, NStrips, 0));
1198
1199 id1 = (unsigned int)it.first;
1200 LOGPRINT << "ID1 shoudl match list of modules above " << id1;
1201 quality_->compact(id1, BadStripList);
1202 SiStripQuality::Range range(BadStripList.begin(), BadStripList.end());
1203 quality_->put(id1, range);
1204 BadStripList.clear();
1205 }
1206
1207 quality_->fillBadComponents();
1208 }
1209
1210 void SiStripHitEffFromCalibTree::totalStatistics() {
1211
1212 int totalfound = 0;
1213 int totaltotal = 0;
1214 double layereff;
1215 int subdetfound[5];
1216 int subdettotal[5];
1217
1218 for (Long_t i = 1; i < 5; i++) {
1219 subdetfound[i] = 0;
1220 subdettotal[i] = 0;
1221 }
1222
1223 for (Long_t i = 1; i <= siStripLayers_; i++) {
1224 layereff = double(layerfound[i]) / double(layertotal[i]);
1225 LOGPRINT << "Layer " << i << " (" << ::layerName(i, showRings_, nTEClayers) << ") has total efficiency " << layereff
1226 << " " << layerfound[i] << "/" << layertotal[i];
1227 totalfound += layerfound[i];
1228 totaltotal += layertotal[i];
1229 if (i < 5) {
1230 subdetfound[1] += layerfound[i];
1231 subdettotal[1] += layertotal[i];
1232 }
1233 if (i >= 5 && i < 11) {
1234 subdetfound[2] += layerfound[i];
1235 subdettotal[2] += layertotal[i];
1236 }
1237 if (i >= 11 && i < 14) {
1238 subdetfound[3] += layerfound[i];
1239 subdettotal[3] += layertotal[i];
1240 }
1241 if (i >= 14) {
1242 subdetfound[4] += layerfound[i];
1243 subdettotal[4] += layertotal[i];
1244 }
1245 }
1246
1247 LOGPRINT << "The total efficiency is " << double(totalfound) / double(totaltotal);
1248 LOGPRINT << " TIB: " << double(subdetfound[1]) / subdettotal[1] << " " << subdetfound[1] << "/"
1249 << subdettotal[1];
1250 LOGPRINT << " TOB: " << double(subdetfound[2]) / subdettotal[2] << " " << subdetfound[2] << "/"
1251 << subdettotal[2];
1252 LOGPRINT << " TID: " << double(subdetfound[3]) / subdettotal[3] << " " << subdetfound[3] << "/"
1253 << subdettotal[3];
1254 LOGPRINT << " TEC: " << double(subdetfound[4]) / subdettotal[4] << " " << subdetfound[4] << "/"
1255 << subdettotal[4];
1256 }
1257
1258 void SiStripHitEffFromCalibTree::makeSummary() {
1259
1260
1261 int nLayers = 34;
1262 if (showRings_)
1263 nLayers = 30;
1264 if (!showEndcapSides_) {
1265 if (!showRings_)
1266 nLayers = siStripLayers_;
1267 else
1268 nLayers = 20;
1269 }
1270
1271 TH1F* found = fs->make<TH1F>("found", "found", nLayers + 1, 0, nLayers + 1);
1272 TH1F* all = fs->make<TH1F>("all", "all", nLayers + 1, 0, nLayers + 1);
1273 TH1F* found2 = fs->make<TH1F>("found2", "found2", nLayers + 1, 0, nLayers + 1);
1274 TH1F* all2 = fs->make<TH1F>("all2", "all2", nLayers + 1, 0, nLayers + 1);
1275
1276 found->SetBinContent(0, -1);
1277 all->SetBinContent(0, 1);
1278
1279
1280 for (Long_t i = 1; i < nLayers + 2; ++i) {
1281 found->SetBinContent(i, 1e-6);
1282 all->SetBinContent(i, 1);
1283 found2->SetBinContent(i, 1e-6);
1284 all2->SetBinContent(i, 1);
1285 }
1286
1287 TCanvas* c7 = new TCanvas("c7", " test ", 10, 10, 800, 600);
1288 c7->SetFillColor(0);
1289 c7->SetGrid();
1290
1291 int nLayers_max = nLayers + 1;
1292 if (!showEndcapSides_)
1293 nLayers_max = 11;
1294 for (Long_t i = 1; i < nLayers_max; ++i) {
1295 LOGPRINT << "Fill only good modules layer " << i << ": S = " << goodlayerfound[i]
1296 << " B = " << goodlayertotal[i];
1297 if (goodlayertotal[i] > 5) {
1298 found->SetBinContent(i, goodlayerfound[i]);
1299 all->SetBinContent(i, goodlayertotal[i]);
1300 }
1301
1302 LOGPRINT << "Filling all modules layer " << i << ": S = " << alllayerfound[i] << " B = " << alllayertotal[i];
1303 if (alllayertotal[i] > 5) {
1304 found2->SetBinContent(i, alllayerfound[i]);
1305 all2->SetBinContent(i, alllayertotal[i]);
1306 }
1307 }
1308
1309
1310 if (!showEndcapSides_) {
1311 for (Long_t i = 11; i < 14; ++i) {
1312 LOGPRINT << "Fill only good modules layer " << i << ": S = " << goodlayerfound[i] + goodlayerfound[i + 3]
1313 << " B = " << goodlayertotal[i] + goodlayertotal[i + 3];
1314 if (goodlayertotal[i] + goodlayertotal[i + 3] > 5) {
1315 found->SetBinContent(i, goodlayerfound[i] + goodlayerfound[i + 3]);
1316 all->SetBinContent(i, goodlayertotal[i] + goodlayertotal[i + 3]);
1317 }
1318 LOGPRINT << "Filling all modules layer " << i << ": S = " << alllayerfound[i] + alllayerfound[i + 3]
1319 << " B = " << alllayertotal[i] + alllayertotal[i + 3];
1320 if (alllayertotal[i] + alllayertotal[i + 3] > 5) {
1321 found2->SetBinContent(i, alllayerfound[i] + alllayerfound[i + 3]);
1322 all2->SetBinContent(i, alllayertotal[i] + alllayertotal[i + 3]);
1323 }
1324 }
1325 for (Long_t i = 17; i < 17 + nTEClayers; ++i) {
1326 LOGPRINT << "Fill only good modules layer " << i - 3
1327 << ": S = " << goodlayerfound[i] + goodlayerfound[i + nTEClayers]
1328 << " B = " << goodlayertotal[i] + goodlayertotal[i + nTEClayers];
1329 if (goodlayertotal[i] + goodlayertotal[i + nTEClayers] > 5) {
1330 found->SetBinContent(i - 3, goodlayerfound[i] + goodlayerfound[i + nTEClayers]);
1331 all->SetBinContent(i - 3, goodlayertotal[i] + goodlayertotal[i + nTEClayers]);
1332 }
1333 LOGPRINT << "Filling all modules layer " << i - 3 << ": S = " << alllayerfound[i] + alllayerfound[i + nTEClayers]
1334 << " B = " << alllayertotal[i] + alllayertotal[i + nTEClayers];
1335 if (alllayertotal[i] + alllayertotal[i + nTEClayers] > 5) {
1336 found2->SetBinContent(i - 3, alllayerfound[i] + alllayerfound[i + nTEClayers]);
1337 all2->SetBinContent(i - 3, alllayertotal[i] + alllayertotal[i + nTEClayers]);
1338 }
1339 }
1340 }
1341
1342 found->Sumw2();
1343 all->Sumw2();
1344
1345 found2->Sumw2();
1346 all2->Sumw2();
1347
1348 TGraphAsymmErrors* gr = fs->make<TGraphAsymmErrors>(nLayers + 1);
1349 gr->SetName("eff_good");
1350 gr->BayesDivide(found, all);
1351
1352 TGraphAsymmErrors* gr2 = fs->make<TGraphAsymmErrors>(nLayers + 1);
1353 gr2->SetName("eff_all");
1354 gr2->BayesDivide(found2, all2);
1355
1356 for (int j = 0; j < nLayers + 1; j++) {
1357 gr->SetPointError(j, 0., 0., gr->GetErrorYlow(j), gr->GetErrorYhigh(j));
1358 gr2->SetPointError(j, 0., 0., gr2->GetErrorYlow(j), gr2->GetErrorYhigh(j));
1359 }
1360
1361 gr->GetXaxis()->SetLimits(0, nLayers);
1362 gr->SetMarkerColor(2);
1363 gr->SetMarkerSize(1.2);
1364 gr->SetLineColor(2);
1365 gr->SetLineWidth(4);
1366 gr->SetMarkerStyle(20);
1367 gr->SetMinimum(effPlotMin_);
1368 gr->SetMaximum(1.001);
1369 gr->GetYaxis()->SetTitle("Efficiency");
1370 gStyle->SetTitleFillColor(0);
1371 gStyle->SetTitleBorderSize(0);
1372 gr->SetTitle(title_);
1373
1374 gr2->GetXaxis()->SetLimits(0, nLayers);
1375 gr2->SetMarkerColor(1);
1376 gr2->SetMarkerSize(1.2);
1377 gr2->SetLineColor(1);
1378 gr2->SetLineWidth(4);
1379 gr2->SetMarkerStyle(21);
1380 gr2->SetMinimum(effPlotMin_);
1381 gr2->SetMaximum(1.001);
1382 gr2->GetYaxis()->SetTitle("Efficiency");
1383 gr2->SetTitle(title_);
1384
1385 for (Long_t k = 1; k < nLayers + 1; k++) {
1386 TString label;
1387 if (showEndcapSides_)
1388 label = getLayerSideName(k);
1389 else
1390 label = ::layerName(k, showRings_, nTEClayers);
1391 if (!showTOB6TEC9_) {
1392 if (k == 10)
1393 label = "";
1394 if (!showRings_ && k == nLayers)
1395 label = "";
1396 if (!showRings_ && showEndcapSides_ && k == 25)
1397 label = "";
1398 }
1399 if (!showRings_) {
1400 if (showEndcapSides_) {
1401 gr->GetXaxis()->SetBinLabel(((k + 1) * 100 + 2) / (nLayers)-4, label);
1402 gr2->GetXaxis()->SetBinLabel(((k + 1) * 100 + 2) / (nLayers)-4, label);
1403 } else {
1404 gr->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-6, label);
1405 gr2->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-6, label);
1406 }
1407 } else {
1408 if (showEndcapSides_) {
1409 gr->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-4, label);
1410 gr2->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-4, label);
1411 } else {
1412 gr->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-7, label);
1413 gr2->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-7, label);
1414 }
1415 }
1416 }
1417
1418 gr->Draw("AP");
1419 gr->GetXaxis()->SetNdivisions(36);
1420
1421 c7->cd();
1422 TPad* overlay = new TPad("overlay", "", 0, 0, 1, 1);
1423 overlay->SetFillStyle(4000);
1424 overlay->SetFillColor(0);
1425 overlay->SetFrameFillStyle(4000);
1426 overlay->Draw("same");
1427 overlay->cd();
1428 if (!showOnlyGoodModules_)
1429 gr2->Draw("AP");
1430
1431 TLegend* leg = new TLegend(0.70, 0.27, 0.88, 0.40);
1432 leg->AddEntry(gr, "Good Modules", "p");
1433 if (!showOnlyGoodModules_)
1434 leg->AddEntry(gr2, "All Modules", "p");
1435 leg->SetTextSize(0.020);
1436 leg->SetFillColor(0);
1437 leg->Draw("same");
1438
1439 c7->SaveAs("Summary.png");
1440 c7->SaveAs("Summary.C");
1441 }
1442
1443 void SiStripHitEffFromCalibTree::makeSummaryVsBx() {
1444 LOGPRINT << "Computing efficiency vs bx";
1445
1446 unsigned int nLayers = siStripLayers_;
1447 if (showRings_)
1448 nLayers = 20;
1449
1450 for (unsigned int ilayer = 1; ilayer < nLayers; ilayer++) {
1451 for (unsigned int ibx = 0; ibx <= nBxInAnOrbit_; ibx++) {
1452 layerfound_vsBX[ilayer]->SetBinContent(ibx, 1e-6);
1453 layertotal_vsBX[ilayer]->SetBinContent(ibx, 1);
1454 }
1455
1456 for (const auto& iterMapvsBx : layerfound_perBx) {
1457 layerfound_vsBX[ilayer]->SetBinContent(iterMapvsBx.first, iterMapvsBx.second[ilayer]);
1458 }
1459 for (const auto& iterMapvsBx : layertotal_perBx) {
1460 if (iterMapvsBx.second[ilayer] > 0) {
1461 layertotal_vsBX[ilayer]->SetBinContent(iterMapvsBx.first, iterMapvsBx.second[ilayer]);
1462 }
1463 }
1464
1465 layerfound_vsBX[ilayer]->Sumw2();
1466 layertotal_vsBX[ilayer]->Sumw2();
1467
1468 TGraphAsymmErrors* geff = fs->make<TGraphAsymmErrors>(nBxInAnOrbit_ - 1);
1469 geff->SetName(Form("effVsBx_layer%i", ilayer));
1470
1471 geff->SetTitle(fmt::format("Hit Efficiency vs bx - {}", ::layerName(ilayer, showRings_, nTEClayers)).c_str());
1472 geff->BayesDivide(layerfound_vsBX[ilayer], layertotal_vsBX[ilayer]);
1473
1474
1475 TGraphAsymmErrors* geff_avg = fs->make<TGraphAsymmErrors>();
1476 geff_avg->SetName(Form("effVsBxAvg_test_layer%i", ilayer));
1477 geff_avg->SetTitle(fmt::format("Hit Efficiency vs bx - {}", ::layerName(ilayer, showRings_, nTEClayers)).c_str());
1478 geff_avg->SetMarkerStyle(20);
1479 int ibx = 0;
1480 int previous_bx = -80;
1481 int delta_bx = 0;
1482 int nbx = 0;
1483 int found = 0;
1484 int total = 0;
1485 double sum_bx = 0;
1486 int ipt = 0;
1487 float low, up, eff;
1488 int firstbx = 0;
1489 for (const auto& iterMapvsBx : layertotal_perBx) {
1490 ibx = iterMapvsBx.first;
1491 delta_bx = ibx - previous_bx;
1492
1493 if (delta_bx > (int)spaceBetweenTrains_ && nbx > 0 && total > 0) {
1494 eff = found / (float)total;
1495
1496 geff_avg->SetPoint(ipt, sum_bx / nbx, eff);
1497 low = TEfficiency::Bayesian(total, found, .683, 1, 1, false);
1498 up = TEfficiency::Bayesian(total, found, .683, 1, 1, true);
1499 geff_avg->SetPointError(ipt, sum_bx / nbx - firstbx, previous_bx - sum_bx / nbx, eff - low, up - eff);
1500 ipt++;
1501 sum_bx = 0;
1502 found = 0;
1503 total = 0;
1504 nbx = 0;
1505 firstbx = ibx;
1506 }
1507 sum_bx += ibx;
1508 found += layerfound_vsBX[ilayer]->GetBinContent(ibx);
1509 total += layertotal_vsBX[ilayer]->GetBinContent(ibx);
1510 nbx++;
1511
1512 previous_bx = ibx;
1513 }
1514
1515 eff = found / (float)total;
1516
1517 geff_avg->SetPoint(ipt, sum_bx / nbx, eff);
1518 low = TEfficiency::Bayesian(total, found, .683, 1, 1, false);
1519 up = TEfficiency::Bayesian(total, found, .683, 1, 1, true);
1520 geff_avg->SetPointError(ipt, sum_bx / nbx - firstbx, previous_bx - sum_bx / nbx, eff - low, up - eff);
1521 }
1522 }
1523
1524 void SiStripHitEffFromCalibTree::computeEff(vector<TH1F*>& vhfound, vector<TH1F*>& vhtotal, string name) {
1525 unsigned int nLayers = siStripLayers_;
1526 if (showRings_)
1527 nLayers = 20;
1528
1529 TH1F* hfound;
1530 TH1F* htotal;
1531
1532 for (unsigned int ilayer = 1; ilayer < nLayers; ilayer++) {
1533 hfound = vhfound[ilayer];
1534 htotal = vhtotal[ilayer];
1535
1536 hfound->Sumw2();
1537 htotal->Sumw2();
1538
1539
1540 for (Long_t i = 0; i < hfound->GetNbinsX() + 1; ++i) {
1541 if (hfound->GetBinContent(i) == 0)
1542 hfound->SetBinContent(i, 1e-6);
1543 if (htotal->GetBinContent(i) == 0)
1544 htotal->SetBinContent(i, 1);
1545 }
1546
1547 TGraphAsymmErrors* geff = fs->make<TGraphAsymmErrors>(hfound->GetNbinsX());
1548 geff->SetName(Form("%s_test_layer%i", name.c_str(), ilayer));
1549 geff->BayesDivide(hfound, htotal);
1550 if (name == "effVsLumi")
1551 geff->SetTitle(
1552 fmt::format("Hit Efficiency vs inst. lumi. - {}", ::layerName(ilayer, showRings_, nTEClayers)).c_str());
1553 if (name == "effVsPU")
1554 geff->SetTitle(fmt::format("Hit Efficiency vs pileup - {}", ::layerName(ilayer, showRings_, nTEClayers)).c_str());
1555 if (name == "effVsCM")
1556 geff->SetTitle(
1557 fmt::format("Hit Efficiency vs common Mode - {}", ::layerName(ilayer, showRings_, nTEClayers)).c_str());
1558 geff->SetMarkerStyle(20);
1559 }
1560 }
1561
1562 void SiStripHitEffFromCalibTree::makeSummaryVsLumi() {
1563 LOGPRINT << "Computing efficiency vs lumi";
1564
1565 if (instLumiHisto->GetEntries())
1566 LOGPRINT << "Avg conditions (avg+/-rms): lumi :" << instLumiHisto->GetMean() << "+/-" << instLumiHisto->GetRMS()
1567 << " pu: " << PUHisto->GetMean() << "+/-" << PUHisto->GetRMS();
1568
1569 else {
1570
1571 unsigned int nLayers = siStripLayers_;
1572 if (showRings_)
1573 nLayers = 20;
1574 unsigned int nLayersForAvg = 0;
1575 float layerLumi = 0;
1576 float layerPU = 0;
1577 float avgLumi = 0;
1578 float avgPU = 0;
1579
1580 LOGPRINT << "Lumi summary: (avg over trajectory measurements)";
1581 for (unsigned int ilayer = 1; ilayer < nLayers; ilayer++) {
1582 layerLumi = layertotal_vsLumi[ilayer]->GetMean();
1583 layerPU = layertotal_vsPU[ilayer]->GetMean();
1584
1585 if (layerLumi != 0 && layerPU != 0) {
1586 avgLumi += layerLumi;
1587 avgPU += layerPU;
1588 nLayersForAvg++;
1589 }
1590 }
1591 avgLumi /= nLayersForAvg;
1592 avgPU /= nLayersForAvg;
1593 LOGPRINT << "Avg conditions: lumi :" << avgLumi << " pu: " << avgPU;
1594 }
1595
1596 computeEff(layerfound_vsLumi, layertotal_vsLumi, "effVsLumi");
1597 computeEff(layerfound_vsPU, layertotal_vsPU, "effVsPU");
1598 }
1599
1600 void SiStripHitEffFromCalibTree::makeSummaryVsCM() {
1601 LOGPRINT << "Computing efficiency vs CM";
1602 computeEff(layerfound_vsCM, layertotal_vsCM, "effVsCM");
1603 }
1604
1605 TString SiStripHitEffFromCalibTree::getLayerSideName(Long_t k) {
1606 TString layername = "";
1607 TString ringlabel = "D";
1608 if (showRings_)
1609 ringlabel = "R";
1610 if (k > 0 && k < 5) {
1611 layername = TString("TIB L") + k;
1612 } else if (k > 4 && k < 11) {
1613 layername = TString("TOB L") + (k - 4);
1614 } else if (k > 10 && k < 14) {
1615 layername = TString("TID- ") + ringlabel + (k - 10);
1616 } else if (k > 13 && k < 17) {
1617 layername = TString("TID+ ") + ringlabel + (k - 13);
1618 } else if (k > 16 && k < 17 + nTEClayers) {
1619 layername = TString("TEC- ") + ringlabel + (k - 16);
1620 } else if (k > 16 + nTEClayers) {
1621 layername = TString("TEC+ ") + ringlabel + (k - 16 - nTEClayers);
1622 }
1623
1624 return layername;
1625 }
1626
1627 std::unique_ptr<SiStripBadStrip> SiStripHitEffFromCalibTree::getNewObject() {
1628
1629
1630 auto obj = std::make_unique<SiStripBadStrip>();
1631
1632 SiStripBadStrip::RegistryIterator rIter = quality_->getRegistryVectorBegin();
1633 SiStripBadStrip::RegistryIterator rIterEnd = quality_->getRegistryVectorEnd();
1634
1635 for (; rIter != rIterEnd; ++rIter) {
1636 SiStripBadStrip::Range range(quality_->getDataVectorBegin() + rIter->ibegin,
1637 quality_->getDataVectorBegin() + rIter->iend);
1638 if (!obj->put(rIter->detid, range))
1639 edm::LogError("SiStripHitEffFromCalibTree")
1640 << "[SiStripHitEffFromCalibTree::getNewObject] detid already exists" << std::endl;
1641 }
1642
1643 return obj;
1644 }
1645
1646 void SiStripHitEffFromCalibTree::setBadComponents(
1647 int i, int component, SiStripQuality::BadComponent& BC, std::stringstream ssV[4][19], int NBadComponent[4][19][4]) {
1648 int napv = detInfo_.getNumberOfApvsAndStripLength(BC.detid).first;
1649
1650 ssV[i][component] << "\n\t\t " << BC.detid << " \t " << BC.BadModule << " \t " << ((BC.BadFibers) & 0x1) << " ";
1651 if (napv == 4)
1652 ssV[i][component] << "x " << ((BC.BadFibers >> 1) & 0x1);
1653
1654 if (napv == 6)
1655 ssV[i][component] << ((BC.BadFibers >> 1) & 0x1) << " " << ((BC.BadFibers >> 2) & 0x1);
1656 ssV[i][component] << " \t " << ((BC.BadApvs) & 0x1) << " " << ((BC.BadApvs >> 1) & 0x1) << " ";
1657 if (napv == 4)
1658 ssV[i][component] << "x x " << ((BC.BadApvs >> 2) & 0x1) << " " << ((BC.BadApvs >> 3) & 0x1);
1659 if (napv == 6)
1660 ssV[i][component] << ((BC.BadApvs >> 2) & 0x1) << " " << ((BC.BadApvs >> 3) & 0x1) << " "
1661 << ((BC.BadApvs >> 4) & 0x1) << " " << ((BC.BadApvs >> 5) & 0x1) << " ";
1662
1663 if (BC.BadApvs) {
1664 NBadComponent[i][0][2] += ((BC.BadApvs >> 5) & 0x1) + ((BC.BadApvs >> 4) & 0x1) + ((BC.BadApvs >> 3) & 0x1) +
1665 ((BC.BadApvs >> 2) & 0x1) + ((BC.BadApvs >> 1) & 0x1) + ((BC.BadApvs) & 0x1);
1666 NBadComponent[i][component][2] += ((BC.BadApvs >> 5) & 0x1) + ((BC.BadApvs >> 4) & 0x1) +
1667 ((BC.BadApvs >> 3) & 0x1) + ((BC.BadApvs >> 2) & 0x1) +
1668 ((BC.BadApvs >> 1) & 0x1) + ((BC.BadApvs) & 0x1);
1669 }
1670 if (BC.BadFibers) {
1671 NBadComponent[i][0][1] += ((BC.BadFibers >> 2) & 0x1) + ((BC.BadFibers >> 1) & 0x1) + ((BC.BadFibers) & 0x1);
1672 NBadComponent[i][component][1] +=
1673 ((BC.BadFibers >> 2) & 0x1) + ((BC.BadFibers >> 1) & 0x1) + ((BC.BadFibers) & 0x1);
1674 }
1675 if (BC.BadModule) {
1676 NBadComponent[i][0][0]++;
1677 NBadComponent[i][component][0]++;
1678 }
1679 }
1680
1681 DEFINE_FWK_MODULE(SiStripHitEffFromCalibTree);