Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-07-10 22:32:53 

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