Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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