Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-26 05:05:42 

0001 //Original Author:  Christopher Edelmaier
0002 //        Created:  Feb. 11, 2010
0003 
0004 // system includes
0005 #include <memory>
0006 #include <string>
0007 #include <iostream>
0008 #include <fstream>
0009 #include <sstream>
0010 
0011 // user includes
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 // ROOT includes
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 // custom made printout
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   // overridden from ConditionDBWriter
0101   void algoAnalyze(const edm::Event& e, const edm::EventSetup& c) override;
0102   std::unique_ptr<SiStripBadStrip> getNewObject() override;
0103 
0104   // native methods
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   // to be used everywhere
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   // for association of informations of the hitEff tree and the event infos tree
0167   map<pair<unsigned int, unsigned int>, array<double, 3> > eventInfos;
0168   // for using events after number of tracks cut
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;  // number of wheels
0228   if (showRings_)
0229     nTEClayers = 7;  // number of rings
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   // read bad modules to mask
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   // initialze counters and histos
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   //Open the ROOT Calib Tree
0323   for (const auto& calibTreeFileName : calibTreeFileNames_) {
0324     LOGPRINT << "Loading file: " << calibTreeFileName;
0325     TFile* CalibTreeFile = TFile::Open(calibTreeFileName.c_str(), "READ");
0326 
0327     // Get event infos
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     // Get hit infos
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     //TLeaf* nHitsLf = calibTree_->GetLeaf("nHits");
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     //Loop through all of the events
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) {  // use rings instead of wheels
0417         if (layer < 14)
0418           layer = 10 + ((id >> 9) & 0x3);  //TID   3 disks and also 3 rings -> use the same container
0419         else
0420           layer = 13 + ((id >> 5) & 0x7);  //TEC
0421       }
0422       //unsigned int nHits = (unsigned int)nHitsLf->GetValue();
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       // if no special tree with event infos, they may be stored in the hit eff tree
0440       if (!foundEventInfos) {
0441         bx = (unsigned int)BunchLf->GetValue();
0442         if (InstLumiLf != nullptr)
0443           instLumi = InstLumiLf->GetValue();  // branch not filled by default
0444         if (PULf != nullptr)
0445           PU = PULf->GetValue();  // branch not filled by default
0446         // Mark new event
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       // Get infos from eventInfos if they exist
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       // Fill event info for events from the anEff tree
0466       // They can differ from the eventInfo tree due to an optional cut on the #tracks when filling the anEff tree
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       //We have two things we want to do, both an XY color plot, and the efficiency measurement
0478       //First, ignore anything that isn't in acceptance and isn't good quality
0479 
0480       if (bunchX_ > 0 && bunchX_ != bx)
0481         continue;
0482 
0483       //if(quality == 1 || accept != 1 || nHits < 8) continue;
0484       if (accept != 1)
0485         continue;
0486       if (useOnlyHighPurityTracks_ && !highPurity)
0487         continue;
0488       if (quality == 1)
0489         badquality = true;
0490 
0491       // don't compute efficiencies in modules from TOB6 and TEC9
0492       if (!showTOB6TEC9_ && (layer_wheel == 10 || layer_wheel == siStripLayers_))
0493         continue;
0494 
0495       // don't use bad modules given in the bad module list
0496       itBadMod = badModules_list.find(id);
0497       if (itBadMod != badModules_list.end())
0498         continue;
0499 
0500       //Now that we have a good event, we need to look at if we expected it or not, and the location
0501       //if we didn't
0502       //Fill the missing hit information first
0503       bool badflag = false;
0504 
0505       // By default uses the old matching method
0506       if (resXSig_ < 0) {
0507         if (isBad == 1)
0508           badflag = true;  // isBad set to false in the tree when resxsig<999.0
0509       } else {
0510         if (isBad == 1 || resxsig > resXSig_)
0511           badflag = true;
0512       }
0513 
0514       // Conversion of positions in strip unit
0515       int nstrips = -9;
0516       float Pitch = -9.0;
0517 
0518       if (resxsig == 1000.0) {  // special treatment, no GeomDetUnit associated in some cases when no cluster found
0519         Pitch = 0.0205;         // maximum
0520         nstrips = 768;          // maximum
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;  //layer01->10
0530         stripCluster = ClusterLocX / Pitch + nstrips / 2.0;
0531 
0532         // For trapezoidal modules: extrapolation of x trajectory position to the y middle of the module
0533         //  for correct comparison with cluster position
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);  // radialy extrapolated x loc position at middle
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       // New matching methods
0559       int tapv = -9;
0560       int capv = -9;
0561       float stripInAPV = 64.;
0562 
0563       if (clusterMatchingMethod_ >= 1) {
0564         badflag = false;          // reset
0565         if (resxsig == 1000.0) {  // default value when no cluster found in the module
0566           badflag = true;         // consider the module inefficient in this case
0567         } else {
0568           if (clusterMatchingMethod_ == 2 ||
0569               clusterMatchingMethod_ == 4) {  // check the distance between cluster and trajectory position
0570             if (abs(stripCluster - stripTrajMid) > clusterTrajDist_)
0571               badflag = true;
0572           }
0573           if (clusterMatchingMethod_ == 3 ||
0574               clusterMatchingMethod_ ==
0575                   4) {  // cluster and traj have to be in the same APV (don't take edges into accounts)
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       //First, figure out if the module already exists in the map of maps
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         //Have to do the decoding for which side to go on (ugh)
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       //Do the one where we don't exclude bad modules!
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       //At this point, both of our maps are loaded with the correct information
0687     }
0688   }  // go to next CalibTreeFile
0689 
0690   makeHotColdMaps();
0691   makeTKMap(autoIneffModTagging_);
0692   makeSQLite();
0693   totalStatistics();
0694   makeSummary();
0695   makeSummaryVsBx();
0696   makeSummaryVsLumi();
0697   if (useCM_)
0698     makeSummaryVsCM();
0699 
0700   ////////////////////////////////////////////////////////////////////////
0701   //try to write out what's in the quality record
0702   /////////////////////////////////////////////////////////////////////////////
0703   int NTkBadComponent[4];  //k: 0=BadModule, 1=BadFiber, 2=BadApv, 3=BadStrips
0704   int NBadComponent[4][19][4];
0705   //legend: NBadComponent[i][j][k]= SubSystem i, layer/disk/wheel j, BadModule/Fiber/Apv k
0706   //     i: 0=TIB, 1=TID, 2=TOB, 3=TEC
0707   //     k: 0=BadModule, 1=BadFiber, 2=BadApv, 3=BadStrips
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     //Full Tk
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     //Single SubSystem
0736     //&&&&&&&&&&&&&&&&&
0737 
0738     int component;
0739     SiStripDetId a(BC[i].detid);
0740     if (a.subdetId() == SiStripDetId::TIB) {
0741       //&&&&&&&&&&&&&&&&&
0742       //TIB
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       //TID
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       //TOB
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       //TEC
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   // Single Strip Info
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   // printout
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   // store also bad modules in log file
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   //Here we make the hot/cold color maps that we love so very much
0914   //Already have access to the data as a private variable
0915   //Create all of the histograms in the TFileService
0916   TH2F* temph2;
0917   for (Long_t maplayer = 1; maplayer <= siStripLayers_; maplayer++) {
0918     //Initialize all of the histograms
0919     if (maplayer > 0 && maplayer <= 4) {
0920       //We are in the TIB
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       //We are in the TOB
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       //We are in the TID
0941       //Split by +/-
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       //We are in the TEC
0966       //Split by +/-
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     //Determine what kind of plot we want to write out
0993     //Loop through the entirety of each layer
0994     //Create an array of the histograms
0995     vector<hit>::const_iterator iter;
0996     for (iter = hits[mylayer].begin(); iter != hits[mylayer].end(); iter++) {
0997       //Looping over the particular layer
0998       //Fill by 360-x to get the proper location to compare with TKMaps of phi
0999       //Also global xy is messed up
1000       if (mylayer > 0 && mylayer <= 4) {
1001         //We are in the TIB
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         //We are in the TOB
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         //We are in the TID
1010         //There are 2 different maps here
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         //if(side == 1) HotColdMaps[(mylayer - 1) + (mylayer - 11)]->Fill(iter->x,iter->y,1.);
1017         //else if(side == 2) HotColdMaps[(mylayer - 1) + (mylayer - 10)]->Fill(iter->x,iter->y,1.);
1018       } else if (mylayer > 13) {
1019         //We are in the TEC
1020         //There are 2 different maps here
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         //if(side == 1) HotColdMaps[(mylayer + 2) + (mylayer - 14)]->Fill(iter->x,iter->y,1.);
1027         //else if(side == 2) HotColdMaps[(mylayer + 2) + (mylayer - 13)]->Fill(iter->x,iter->y,1.);
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     //Loop over every layer, extracting the information from
1063     //the map of the efficiencies
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       //We should be in the layer in question, and looping over all of the modules in said layer
1071       //Generate the list for the TKmap, and the bad module list
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           //We have a bad module, put it in the list!
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           //Fill the bad list with empty results for every module
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       //Put any module into the TKMap
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       //Add the number of hits in the layer
1117       layertotal[i] += long(myden);
1118       layerfound[i] += long(mynum);
1119     }
1120 
1121     if (autoTagging) {
1122       //Compute threshold to use for each layer
1123       hEffInLayer->GetXaxis()->SetRange(3, hEffInLayer->GetNbinsX() + 1);  // Remove from the avg modules below 1%
1124       layer_min_eff =
1125           hEffInLayer->GetMean() - 2.5 * hEffInLayer->GetRMS();  // uses RMS in case the distribution is wide
1126       if (threshold_ > 2.5 * hEffInLayer->GetRMS())
1127         layer_min_eff = hEffInLayer->GetMean() - threshold_;  // otherwise uses the parameter '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         // Second loop over modules to tag inefficient ones
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         // upper limit on the efficiency
1142         myeff_up = TEfficiency::Bayesian(myden, mynum, .99, 1, 1, true);
1143         if ((myden >= nModsMin_) && (myeff_up < layer_min_eff)) {
1144           //We have a bad module, put it in the list!
1145           BadModules[ih.first] = myeff;
1146           tkmapbad->fillc(ih.first, 255, 0, 0);
1147           t_isTaggedIneff = true;
1148         } else {
1149           //Fill the bad list with empty results for every module
1150           tkmapbad->fillc(ih.first, 255, 255, 255);
1151           t_isTaggedIneff = false;
1152         }
1153         if (myeff_up < layer_min_eff + 0.08)  // printing message also for modules slighly above (8%) the limit
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   //Generate the SQLite file for use in the Database of the bad modules!
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   //This is the list of the bad strips, use to mask out entire APVs
1190   //Now simply go through the bad hit list and mask out things that
1191   //are bad!
1192   for (const auto& it : BadModules) {
1193     //We need to figure out how many strips are in this particular module
1194     //To Mask correctly!
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     //Now compact into a single bad module
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   //Fill all the bad components now
1207   quality_->fillBadComponents();
1208 }
1209 
1210 void SiStripHitEffFromCalibTree::totalStatistics() {
1211   //Calculate the statistics by layer
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   //setTDRStyle();
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   // first bin only to keep real data off the y axis so set to -1
1276   found->SetBinContent(0, -1);
1277   all->SetBinContent(0, 1);
1278 
1279   // new ROOT version: TGraph::Divide don't handle null or negative values
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;  // barrel+endcap
1292   if (!showEndcapSides_)
1293     nLayers_max = 11;  // barrel
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   // endcap - merging sides
1310   if (!showEndcapSides_) {
1311     for (Long_t i = 11; i < 14; ++i) {  // TID disks
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) {  // TEC disks
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     //Average over trains
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       // consider a new train
1493       if (delta_bx > (int)spaceBetweenTrains_ && nbx > 0 && total > 0) {
1494         eff = found / (float)total;
1495         //LOGPRINT<<"new train "<<ipt<<" "<<sum_bx/nbx<<" "<<eff<<endl;
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     //last train
1515     eff = found / (float)total;
1516     //LOGPRINT<<"new train "<<ipt<<" "<<sum_bx/nbx<<" "<<eff<<endl;
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     // new ROOT version: TGraph::Divide don't handle null or negative values
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())  // from infos per event
1566     LOGPRINT << "Avg conditions (avg+/-rms):   lumi :" << instLumiHisto->GetMean() << "+/-" << instLumiHisto->GetRMS()
1567              << "   pu: " << PUHisto->GetMean() << "+/-" << PUHisto->GetRMS();
1568 
1569   else {  // from infos per hit
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       //LOGPRINT<<" layer "<<ilayer<<"  lumi: "<<layerLumi<<"  pu: "<<layerPU<<endl;
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   //Need this for a Condition DB Writer
1629   //Initialize a return variable
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!