Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-04-22 22:55:18 

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