CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
SiStripHitEffFromCalibTree.cc
Go to the documentation of this file.
1 //Original Author: Christopher Edelmaier
2 // Created: Feb. 11, 2010
3 #include <memory>
4 #include <string>
5 #include <iostream>
6 #include <fstream>
7 #include <sstream>
8 
9 #include "FWCore/Framework/interface/Frameworkfwd.h"
10 #include "FWCore/Framework/interface/Event.h"
11 #include "FWCore/Framework/interface/MakerMacros.h"
12 #include "FWCore/Framework/interface/ESHandle.h"
13 #include "FWCore/Framework/interface/EventSetup.h"
14 #include "FWCore/ParameterSet/interface/ParameterSet.h"
15 #include "FWCore/ParameterSet/interface/FileInPath.h"
16 
17 #include "CalibTracker/SiStripHitEfficiency/interface/HitEff.h"
18 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
19 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
20 #include "DataFormats/Common/interface/Handle.h"
21 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
22 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
23 #include "DataFormats/GeometryVector/interface/LocalVector.h"
24 #include "DataFormats/TrackReco/interface/Track.h"
25 #include "DataFormats/TrackReco/interface/TrackFwd.h"
26 #include "DataFormats/TrackReco/interface/TrackExtra.h"
27 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
28 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
29 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
30 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
31 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
32 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
33 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
34 #include "CalibTracker/SiStripHitEfficiency/interface/TrajectoryAtInvalidHit.h"
35 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
36 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
37 #include "RecoLocalTracker/ClusterParameterEstimator/interface/StripClusterParameterEstimator.h"
38 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
39 #include "DataFormats/TrackReco/interface/DeDxData.h"
40 #include "DataFormats/DetId/interface/DetIdCollection.h"
41 #include "TrackingTools/DetLayers/interface/DetLayer.h"
42 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
43 
44 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
45 #include "RecoLocalTracker/SiStripClusterizer/interface/SiStripClusterInfo.h"
46 #include "CalibTracker/Records/interface/SiStripDetCablingRcd.h"
47 #include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h"
48 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
49 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
50 #include "CalibTracker/SiStripCommon/interface/SiStripDetInfoFileReader.h"
51 #include "DataFormats/Common/interface/DetSetVector.h"
52 #include "DataFormats/Common/interface/DetSetVectorNew.h"
53 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
54 
55 #include "DataFormats/MuonReco/interface/Muon.h"
56 #include "DataFormats/MuonReco/interface/MuonFwd.h"
57 
58 #include "FWCore/ServiceRegistry/interface/Service.h"
59 #include "CommonTools/UtilAlgos/interface/TFileService.h"
60 #include "CommonTools/TrackerMap/interface/TrackerMap.h"
61 #include "CommonTools/ConditionDBWriter/interface/ConditionDBWriter.h"
62 
63 #include "TFile.h"
64 #include "TCanvas.h"
65 #include "TObjString.h"
66 #include "TString.h"
67 #include "TH1F.h"
68 #include "TH2F.h"
69 #include "TProfile.h"
70 #include "TF1.h"
71 #include "TROOT.h"
72 #include "TTree.h"
73 #include "TStyle.h"
74 #include "TLeaf.h"
75 #include "TGaxis.h"
76 #include "TGraphAsymmErrors.h"
77 #include "TLatex.h"
78 #include "TLegend.h"
79 #include "TEfficiency.h"
80 
81 #include <memory>
82 
83 using namespace edm;
84 using namespace reco;
85 using namespace std;
86 
87 struct hit {
88  double x;
89  double y;
90  double z;
91  unsigned int id;
92 };
93 
94 class SiStripHitEffFromCalibTree : public ConditionDBWriter<SiStripBadStrip> {
95 public:
96  explicit SiStripHitEffFromCalibTree(const edm::ParameterSet&);
97  ~SiStripHitEffFromCalibTree() override;
98 
99 private:
100  void algoBeginJob(const edm::EventSetup&) override;
101  void algoEndJob() override;
102  void algoAnalyze(const edm::Event& e, const edm::EventSetup& c) override;
103  void SetBadComponents(int i,
104  int component,
105  SiStripQuality::BadComponent& BC,
106  std::stringstream ssV[4][19],
107  int NBadComponent[4][19][4]);
108  void makeTKMap(bool autoTagging);
109  void makeHotColdMaps();
110  void makeSQLite();
111  void totalStatistics();
112  void makeSummary();
113  void makeSummaryVsBx();
114  void ComputeEff(vector<TH1F*>& vhfound, vector<TH1F*>& vhtotal, string name);
115  void makeSummaryVsLumi();
116  void makeSummaryVsCM();
117  TString GetLayerName(Long_t k);
118  TString GetLayerSideName(Long_t k);
119  float calcPhi(float x, float y);
120 
121  edm::Service<TFileService> fs;
122  SiStripDetInfo _detInfo;
123  edm::FileInPath FileInPath_;
124  SiStripQuality* quality_;
125  std::unique_ptr<SiStripBadStrip> getNewObject() override;
126 
127  TTree* CalibTree;
128  vector<string> CalibTreeFilenames;
129  float threshold;
130  unsigned int nModsMin;
131  unsigned int doSummary;
132  string _badModulesFile;
133  bool _autoIneffModTagging;
134  unsigned int _clusterMatchingMethod;
135  float _ResXSig;
136  float _clusterTrajDist;
137  float _stripsApvEdge;
138  bool _useOnlyHighPurityTracks;
139  unsigned int _bunchx;
140  unsigned int _spaceBetweenTrains;
141  bool _useCM;
142  bool _showEndcapSides;
143  bool _showRings;
144  bool _showTOB6TEC9;
145  bool _showOnlyGoodModules;
146  float _tkMapMin;
147  float _effPlotMin;
148  TString _title;
149 
150  edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> _tkGeomToken;
151  edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> _tTopoToken;
152 
153  unsigned int nTEClayers;
154 
155  TH1F* bxHisto;
156  TH1F* instLumiHisto;
157  TH1F* PUHisto;
158 
159  // for association of informations of the hitEff tree and the event infos tree
160  map<pair<unsigned int, unsigned int>, array<double, 3> > eventInfos;
161 
162  vector<hit> hits[23];
163  vector<TH2F*> HotColdMaps;
164  map<unsigned int, pair<unsigned int, unsigned int> > modCounter[23];
165  TrackerMap* tkmap;
166  TrackerMap* tkmapbad;
167  TrackerMap* tkmapeff;
168  TrackerMap* tkmapnum;
169  TrackerMap* tkmapden;
170  long layerfound[23];
171  long layertotal[23];
172  map<unsigned int, vector<int> > layerfound_perBx;
173  map<unsigned int, vector<int> > layertotal_perBx;
174  vector<TH1F*> layerfound_vsLumi;
175  vector<TH1F*> layertotal_vsLumi;
176  vector<TH1F*> layerfound_vsPU;
177  vector<TH1F*> layertotal_vsPU;
178  vector<TH1F*> layerfound_vsCM;
179  vector<TH1F*> layertotal_vsCM;
180  int goodlayertotal[35];
181  int goodlayerfound[35];
182  int alllayertotal[35];
183  int alllayerfound[35];
184  map<unsigned int, double> BadModules;
185 };
186 
187 SiStripHitEffFromCalibTree::SiStripHitEffFromCalibTree(const edm::ParameterSet& conf)
188  : ConditionDBWriter<SiStripBadStrip>(conf), FileInPath_(SiStripDetInfoFileReader::kDefaultFile) {
189  CalibTreeFilenames = conf.getUntrackedParameter<vector<std::string> >("CalibTreeFilenames");
190  threshold = conf.getParameter<double>("Threshold");
191  nModsMin = conf.getParameter<int>("nModsMin");
192  doSummary = conf.getParameter<int>("doSummary");
193  _badModulesFile = conf.getUntrackedParameter<std::string>("BadModulesFile", "");
194  _autoIneffModTagging = conf.getUntrackedParameter<bool>("AutoIneffModTagging", false);
195  _clusterMatchingMethod = conf.getUntrackedParameter<int>("ClusterMatchingMethod", 0);
196  _ResXSig = conf.getUntrackedParameter<double>("ResXSig", -1);
197  _clusterTrajDist = conf.getUntrackedParameter<double>("ClusterTrajDist", 64.0);
198  _stripsApvEdge = conf.getUntrackedParameter<double>("StripsApvEdge", 10.0);
199  _useOnlyHighPurityTracks = conf.getUntrackedParameter<bool>("UseOnlyHighPurityTracks", true);
200  _bunchx = conf.getUntrackedParameter<int>("BunchCrossing", 0);
201  _spaceBetweenTrains = conf.getUntrackedParameter<int>("SpaceBetweenTrains", 25);
202  _useCM = conf.getUntrackedParameter<bool>("UseCommonMode", false);
203  _showEndcapSides = conf.getUntrackedParameter<bool>("ShowEndcapSides", true);
204  _showRings = conf.getUntrackedParameter<bool>("ShowRings", false);
205  _showTOB6TEC9 = conf.getUntrackedParameter<bool>("ShowTOB6TEC9", false);
206  _showOnlyGoodModules = conf.getUntrackedParameter<bool>("ShowOnlyGoodModules", false);
207  _tkMapMin = conf.getUntrackedParameter<double>("TkMapMin", 0.9);
208  _effPlotMin = conf.getUntrackedParameter<double>("EffPlotMin", 0.9);
209  _title = conf.getParameter<std::string>("Title");
210  _tkGeomToken = esConsumes();
211  _tTopoToken = esConsumes();
212  _detInfo = SiStripDetInfoFileReader::read(FileInPath_.fullPath());
213 
214  nTEClayers = 9; // number of wheels
215  if (_showRings)
216  nTEClayers = 7; // number of rings
217 
218  quality_ = new SiStripQuality(_detInfo);
219 }
220 
221 SiStripHitEffFromCalibTree::~SiStripHitEffFromCalibTree() {}
222 
223 void SiStripHitEffFromCalibTree::algoBeginJob(const edm::EventSetup&) {
224  //I have no idea what goes here
225  //fs->make<TTree>("HitEffHistos","Tree of the inefficient hit histograms");
226 }
227 
228 void SiStripHitEffFromCalibTree::algoEndJob() {
229  //Still have no idea what goes here
230 }
231 
232 void SiStripHitEffFromCalibTree::algoAnalyze(const edm::Event& e, const edm::EventSetup& c) {
233  const auto& tkgeom = c.getData(_tkGeomToken);
234  const auto& tTopo = c.getData(_tTopoToken);
235 
236  // read bad modules to mask
237  ifstream badModules_file;
238  set<uint32_t> badModules_list;
239  if (!_badModulesFile.empty()) {
240  badModules_file.open(_badModulesFile.c_str());
241  uint32_t badmodule_detid;
242  int mods, fiber1, fiber2, fiber3;
243  if (badModules_file.is_open()) {
244  string line;
245  while (getline(badModules_file, line)) {
246  if (badModules_file.eof())
247  continue;
248  stringstream ss(line);
249  ss >> badmodule_detid >> mods >> fiber1 >> fiber2 >> fiber3;
250  if (badmodule_detid != 0 && mods == 1 && (fiber1 == 1 || fiber2 == 1 || fiber3 == 1))
251  badModules_list.insert(badmodule_detid);
252  }
253  badModules_file.close();
254  }
255  }
256  if (!badModules_list.empty())
257  cout << "Remove additionnal bad modules from the analysis: " << endl;
258  set<uint32_t>::iterator itBadMod;
259  for (itBadMod = badModules_list.begin(); itBadMod != badModules_list.end(); ++itBadMod)
260  cout << " " << *itBadMod << endl;
261 
262  // initialze counters and histos
263 
264  bxHisto = fs->make<TH1F>("bx", "bx", 3600, 0, 3600);
265  instLumiHisto = fs->make<TH1F>("instLumi", "inst. lumi.", 250, 0, 25000);
266  PUHisto = fs->make<TH1F>("PU", "PU", 200, 0, 200);
267 
268  for (int l = 0; l < 35; l++) {
269  goodlayertotal[l] = 0;
270  goodlayerfound[l] = 0;
271  alllayertotal[l] = 0;
272  alllayerfound[l] = 0;
273  }
274 
275  TH1F* resolutionPlots[23];
276  for (Long_t ilayer = 0; ilayer < 23; ilayer++) {
277  resolutionPlots[ilayer] =
278  fs->make<TH1F>(Form("resol_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 125, -125, 125);
279  resolutionPlots[ilayer]->GetXaxis()->SetTitle("trajX-clusX [strip unit]");
280 
281  layerfound_vsLumi.push_back(
282  fs->make<TH1F>(Form("layerfound_vsLumi_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 100, 0, 25000));
283  layertotal_vsLumi.push_back(
284  fs->make<TH1F>(Form("layertotal_vsLumi_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 100, 0, 25000));
285  layerfound_vsPU.push_back(
286  fs->make<TH1F>(Form("layerfound_vsPU_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 45, 0, 90));
287  layertotal_vsPU.push_back(
288  fs->make<TH1F>(Form("layertotal_vsPU_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 45, 0, 90));
289 
290  if (_useCM) {
291  layerfound_vsCM.push_back(
292  fs->make<TH1F>(Form("layerfound_vsCM_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 20, 0, 400));
293  layertotal_vsCM.push_back(
294  fs->make<TH1F>(Form("layertotal_vsCM_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 20, 0, 400));
295  }
296  layertotal[ilayer] = 0;
297  layerfound[ilayer] = 0;
298  }
299 
300  if (!_autoIneffModTagging)
301  cout << "A module is bad if efficiency < " << threshold << " and has at least " << nModsMin << " nModsMin." << endl;
302  else
303  cout << "A module is bad if efficiency < the avg in layer - " << threshold << " and has at least " << nModsMin
304  << " nModsMin." << endl;
305 
306  unsigned int run, evt, bx{0};
307  double instLumi, PU;
308 
309  //Open the ROOT Calib Tree
310  for (unsigned int ifile = 0; ifile < CalibTreeFilenames.size(); ifile++) {
311  cout << "Loading file: " << CalibTreeFilenames[ifile] << endl;
312  TFile* CalibTreeFile = TFile::Open(CalibTreeFilenames[ifile].c_str(), "READ");
313 
314  // Get event infos
315  bool foundEventInfos = false;
316  try {
317  CalibTreeFile->cd("eventInfo");
318  } catch (exception& e) {
319  cout << "No event infos tree" << endl;
320  }
321  TTree* EventTree = (TTree*)(gDirectory->Get("tree"));
322 
323  TLeaf* runLf;
324  TLeaf* evtLf;
325  TLeaf* BunchLf;
326  TLeaf* InstLumiLf;
327  TLeaf* PULf;
328  if (EventTree) {
329  cout << "Found event infos tree" << endl;
330 
331  runLf = EventTree->GetLeaf("run");
332  evtLf = EventTree->GetLeaf("event");
333 
334  BunchLf = EventTree->GetLeaf("bx");
335  InstLumiLf = EventTree->GetLeaf("instLumi");
336  PULf = EventTree->GetLeaf("PU");
337 
338  int nevt = EventTree->GetEntries();
339  if (nevt)
340  foundEventInfos = true;
341 
342  for (int j = 0; j < nevt; j++) {
343  EventTree->GetEntry(j);
344  run = runLf->GetValue();
345  evt = evtLf->GetValue();
346  bx = BunchLf->GetValue();
347  instLumi = InstLumiLf->GetValue();
348  PU = PULf->GetValue();
349 
350  bxHisto->Fill(bx);
351  instLumiHisto->Fill(instLumi);
352  PUHisto->Fill(PU);
353 
354  eventInfos[make_pair(run, evt)] = array<double, 3>{{(double)bx, instLumi, PU}};
355  }
356  }
357 
358  // Get hit infos
359  CalibTreeFile->cd("anEff");
360  CalibTree = (TTree*)(gDirectory->Get("traj"));
361 
362  runLf = CalibTree->GetLeaf("run");
363  evtLf = CalibTree->GetLeaf("event");
364  TLeaf* BadLf = CalibTree->GetLeaf("ModIsBad");
365  TLeaf* sistripLf = CalibTree->GetLeaf("SiStripQualBad");
366  TLeaf* idLf = CalibTree->GetLeaf("Id");
367  TLeaf* acceptLf = CalibTree->GetLeaf("withinAcceptance");
368  TLeaf* layerLf = CalibTree->GetLeaf("layer");
369  //TLeaf* nHitsLf = CalibTree->GetLeaf("nHits");
370  TLeaf* highPurityLf = CalibTree->GetLeaf("highPurity");
371  TLeaf* xLf = CalibTree->GetLeaf("TrajGlbX");
372  TLeaf* yLf = CalibTree->GetLeaf("TrajGlbY");
373  TLeaf* zLf = CalibTree->GetLeaf("TrajGlbZ");
374  TLeaf* ResXSigLf = CalibTree->GetLeaf("ResXSig");
375  TLeaf* TrajLocXLf = CalibTree->GetLeaf("TrajLocX");
376  TLeaf* TrajLocYLf = CalibTree->GetLeaf("TrajLocY");
377  TLeaf* ClusterLocXLf = CalibTree->GetLeaf("ClusterLocX");
378  BunchLf = CalibTree->GetLeaf("bunchx");
379  InstLumiLf = CalibTree->GetLeaf("instLumi");
380  PULf = CalibTree->GetLeaf("PU");
381  TLeaf* CMLf = nullptr;
382  if (_useCM)
383  CMLf = CalibTree->GetLeaf("commonMode");
384 
385  int nevents = CalibTree->GetEntries();
386  cout << "Successfully loaded analyze function with " << nevents << " events!\n";
387 
388  map<pair<unsigned int, unsigned int>, array<double, 3> >::iterator itEventInfos;
389 
390  //Loop through all of the events
391  for (int j = 0; j < nevents; j++) {
392  CalibTree->GetEntry(j);
393  run = (unsigned int)runLf->GetValue();
394  evt = (unsigned int)evtLf->GetValue();
395  unsigned int isBad = (unsigned int)BadLf->GetValue();
396  unsigned int quality = (unsigned int)sistripLf->GetValue();
397  unsigned int id = (unsigned int)idLf->GetValue();
398  unsigned int accept = (unsigned int)acceptLf->GetValue();
399  unsigned int layer_wheel = (unsigned int)layerLf->GetValue();
400  unsigned int layer = layer_wheel;
401  if (_showRings && layer > 10) { // use rings instead of wheels
402  if (layer < 14)
403  layer = 10 + ((id >> 9) & 0x3); //TID 3 disks and also 3 rings -> use the same container
404  else
405  layer = 13 + ((id >> 5) & 0x7); //TEC
406  }
407  //unsigned int nHits = (unsigned int)nHitsLf->GetValue();
408  bool highPurity = (bool)highPurityLf->GetValue();
409  double x = xLf->GetValue();
410  double y = yLf->GetValue();
411  double z = zLf->GetValue();
412  double resxsig = ResXSigLf->GetValue();
413  double TrajLocX = TrajLocXLf->GetValue();
414  double TrajLocY = TrajLocYLf->GetValue();
415  double ClusterLocX = ClusterLocXLf->GetValue();
416  double TrajLocXMid;
417  double stripTrajMid;
418  double stripCluster;
419  bool badquality = false;
420 
421  instLumi = 0;
422  PU = 0;
423 
424  // if no special tree with event infos, they may be stored in the hit eff tree
425  if (!foundEventInfos) {
426  bx = (unsigned int)BunchLf->GetValue();
427  if (InstLumiLf != nullptr)
428  instLumi = InstLumiLf->GetValue(); // branch not filled by default
429  if (PULf != nullptr)
430  PU = PULf->GetValue(); // branch not filled by default
431  }
432  int CM = -100;
433  if (_useCM)
434  CM = CMLf->GetValue();
435 
436  // Get infos from eventInfos if they exist
437  if (foundEventInfos) {
438  itEventInfos = eventInfos.find(make_pair(run, evt));
439  if (itEventInfos != eventInfos.end()) {
440  bx = itEventInfos->second[0];
441  instLumi = itEventInfos->second[1];
442  PU = itEventInfos->second[2];
443  }
444  }
445 
446  //We have two things we want to do, both an XY color plot, and the efficiency measurement
447  //First, ignore anything that isn't in acceptance and isn't good quality
448 
449  if (_bunchx > 0 && _bunchx != bx)
450  continue;
451 
452  //if(quality == 1 || accept != 1 || nHits < 8) continue;
453  if (accept != 1)
454  continue;
455  if (_useOnlyHighPurityTracks && !highPurity)
456  continue;
457  if (quality == 1)
458  badquality = true;
459 
460  // don't compute efficiencies in modules from TOB6 and TEC9
461  if (!_showTOB6TEC9 && (layer_wheel == 10 || layer_wheel == 22))
462  continue;
463 
464  // don't use bad modules given in the bad module list
465  itBadMod = badModules_list.find(id);
466  if (itBadMod != badModules_list.end())
467  continue;
468 
469  //Now that we have a good event, we need to look at if we expected it or not, and the location
470  //if we didn't
471  //Fill the missing hit information first
472  bool badflag = false;
473 
474  // By default uses the old matching method
475  if (_ResXSig < 0) {
476  if (isBad == 1)
477  badflag = true; // isBad set to false in the tree when resxsig<999.0
478  } else {
479  if (isBad == 1 || resxsig > _ResXSig)
480  badflag = true;
481  }
482 
483  // Conversion of positions in strip unit
484  int nstrips = -9;
485  float Pitch = -9.0;
486 
487  if (resxsig == 1000.0) { // special treatment, no GeomDetUnit associated in some cases when no cluster found
488  Pitch = 0.0205; // maximum
489  nstrips = 768; // maximum
490  stripTrajMid = TrajLocX / Pitch + nstrips / 2.0;
491  stripCluster = ClusterLocX / Pitch + nstrips / 2.0;
492  } else {
493  DetId ClusterDetId(id);
494  const StripGeomDetUnit* stripdet = (const StripGeomDetUnit*)tkgeom.idToDetUnit(ClusterDetId);
495  const StripTopology& Topo = stripdet->specificTopology();
496  nstrips = Topo.nstrips();
497  Pitch = stripdet->surface().bounds().width() / Topo.nstrips();
498  stripTrajMid = TrajLocX / Pitch + nstrips / 2.0; //layer01->10
499  stripCluster = ClusterLocX / Pitch + nstrips / 2.0;
500 
501  // For trapezoidal modules: extrapolation of x trajectory position to the y middle of the module
502  // for correct comparison with cluster position
503  float hbedge = 0;
504  float htedge = 0;
505  float hapoth = 0;
506  if (layer >= 11) {
507  const BoundPlane& plane = stripdet->surface();
508  const TrapezoidalPlaneBounds* trapezoidalBounds(
509  dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
510  std::array<const float, 4> const& parameters = (*trapezoidalBounds).parameters();
511  hbedge = parameters[0];
512  htedge = parameters[1];
513  hapoth = parameters[3];
514  TrajLocXMid = TrajLocX / (1 + (htedge - hbedge) * TrajLocY / (htedge + hbedge) /
515  hapoth); // radialy extrapolated x loc position at middle
516  stripTrajMid = TrajLocXMid / Pitch + nstrips / 2.0;
517  }
518  }
519 
520  if (!badquality && layer < 23) {
521  if (resxsig != 1000.0)
522  resolutionPlots[layer]->Fill(stripTrajMid - stripCluster);
523  else
524  resolutionPlots[layer]->Fill(1000);
525  }
526 
527  // New matching methods
528  int tapv = -9;
529  int capv = -9;
530  float stripInAPV = 64.;
531 
532  if (_clusterMatchingMethod >= 1) {
533  badflag = false; // reset
534  if (resxsig == 1000.0) { // default value when no cluster found in the module
535  badflag = true; // consider the module inefficient in this case
536  } else {
537  if (_clusterMatchingMethod == 2 ||
538  _clusterMatchingMethod == 4) { // check the distance between cluster and trajectory position
539  if (abs(stripCluster - stripTrajMid) > _clusterTrajDist)
540  badflag = true;
541  }
542  if (_clusterMatchingMethod == 3 ||
543  _clusterMatchingMethod ==
544  4) { // cluster and traj have to be in the same APV (don't take edges into accounts)
545  tapv = (int)stripTrajMid / 128;
546  capv = (int)stripCluster / 128;
547  stripInAPV = stripTrajMid - tapv * 128;
548 
549  if (stripInAPV < _stripsApvEdge || stripInAPV > 128 - _stripsApvEdge)
550  continue;
551  if (tapv != capv)
552  badflag = true;
553  }
554  }
555  }
556 
557  if (badflag && !badquality) {
558  hit temphit;
559  temphit.x = x;
560  temphit.y = y;
561  temphit.z = z;
562  temphit.id = id;
563  hits[layer].push_back(temphit);
564  }
565  pair<unsigned int, unsigned int> newgoodpair(1, 1);
566  pair<unsigned int, unsigned int> newbadpair(1, 0);
567  //First, figure out if the module already exists in the map of maps
568  map<unsigned int, pair<unsigned int, unsigned int> >::iterator it = modCounter[layer].find(id);
569  if (!badquality) {
570  if (it == modCounter[layer].end()) {
571  if (badflag)
572  modCounter[layer][id] = newbadpair;
573  else
574  modCounter[layer][id] = newgoodpair;
575  } else {
576  ((*it).second.first)++;
577  if (!badflag)
578  ((*it).second.second)++;
579  }
580 
581  if (layerfound_perBx.find(bx) == layerfound_perBx.end()) {
582  layerfound_perBx[bx] = vector<int>(23, 0);
583  layertotal_perBx[bx] = vector<int>(23, 0);
584  }
585  if (!badflag)
586  layerfound_perBx[bx][layer]++;
587  layertotal_perBx[bx][layer]++;
588 
589  if (!badflag)
590  layerfound_vsLumi[layer]->Fill(instLumi);
591  layertotal_vsLumi[layer]->Fill(instLumi);
592  if (!badflag)
593  layerfound_vsPU[layer]->Fill(PU);
594  layertotal_vsPU[layer]->Fill(PU);
595 
596  if (_useCM) {
597  if (!badflag)
598  layerfound_vsCM[layer]->Fill(CM);
599  layertotal_vsCM[layer]->Fill(CM);
600  }
601 
602  //Have to do the decoding for which side to go on (ugh)
603  if (layer <= 10) {
604  if (!badflag)
605  goodlayerfound[layer]++;
606  goodlayertotal[layer]++;
607  } else if (layer > 10 && layer < 14) {
608  if (((id >> 13) & 0x3) == 1) {
609  if (!badflag)
610  goodlayerfound[layer]++;
611  goodlayertotal[layer]++;
612  } else if (((id >> 13) & 0x3) == 2) {
613  if (!badflag)
614  goodlayerfound[layer + 3]++;
615  goodlayertotal[layer + 3]++;
616  }
617  } else if (layer > 13 && layer <= 22) {
618  if (((id >> 18) & 0x3) == 1) {
619  if (!badflag)
620  goodlayerfound[layer + 3]++;
621  goodlayertotal[layer + 3]++;
622  } else if (((id >> 18) & 0x3) == 2) {
623  if (!badflag)
624  goodlayerfound[layer + 3 + nTEClayers]++;
625  goodlayertotal[layer + 3 + nTEClayers]++;
626  }
627  }
628  }
629  //Do the one where we don't exclude bad modules!
630  if (layer <= 10) {
631  if (!badflag)
632  alllayerfound[layer]++;
633  alllayertotal[layer]++;
634  } else if (layer > 10 && layer < 14) {
635  if (((id >> 13) & 0x3) == 1) {
636  if (!badflag)
637  alllayerfound[layer]++;
638  alllayertotal[layer]++;
639  } else if (((id >> 13) & 0x3) == 2) {
640  if (!badflag)
641  alllayerfound[layer + 3]++;
642  alllayertotal[layer + 3]++;
643  }
644  } else if (layer > 13 && layer <= 22) {
645  if (((id >> 18) & 0x3) == 1) {
646  if (!badflag)
647  alllayerfound[layer + 3]++;
648  alllayertotal[layer + 3]++;
649  } else if (((id >> 18) & 0x3) == 2) {
650  if (!badflag)
651  alllayerfound[layer + 3 + nTEClayers]++;
652  alllayertotal[layer + 3 + nTEClayers]++;
653  }
654  }
655  //At this point, both of our maps are loaded with the correct information
656  }
657  } // go to next CalibTreeFile
658 
659  makeHotColdMaps();
660  makeTKMap(_autoIneffModTagging);
661  makeSQLite();
662  totalStatistics();
663  makeSummary();
664  makeSummaryVsBx();
665  makeSummaryVsLumi();
666  if (_useCM)
667  makeSummaryVsCM();
668 
670  //try to write out what's in the quality record
672  int NTkBadComponent[4]; //k: 0=BadModule, 1=BadFiber, 2=BadApv, 3=BadStrips
673  int NBadComponent[4][19][4];
674  //legend: NBadComponent[i][j][k]= SubSystem i, layer/disk/wheel j, BadModule/Fiber/Apv k
675  // i: 0=TIB, 1=TID, 2=TOB, 3=TEC
676  // k: 0=BadModule, 1=BadFiber, 2=BadApv, 3=BadStrips
677  std::stringstream ssV[4][19];
678 
679  for (int i = 0; i < 4; ++i) {
680  NTkBadComponent[i] = 0;
681  for (int j = 0; j < 19; ++j) {
682  ssV[i][j].str("");
683  for (int k = 0; k < 4; ++k)
684  NBadComponent[i][j][k] = 0;
685  }
686  }
687 
688  std::vector<SiStripQuality::BadComponent> BC = quality_->getBadComponentList();
689 
690  for (size_t i = 0; i < BC.size(); ++i) {
691  //&&&&&&&&&&&&&
692  //Full Tk
693  //&&&&&&&&&&&&&
694 
695  if (BC[i].BadModule)
696  NTkBadComponent[0]++;
697  if (BC[i].BadFibers)
698  NTkBadComponent[1] += ((BC[i].BadFibers >> 2) & 0x1) + ((BC[i].BadFibers >> 1) & 0x1) + ((BC[i].BadFibers) & 0x1);
699  if (BC[i].BadApvs)
700  NTkBadComponent[2] += ((BC[i].BadApvs >> 5) & 0x1) + ((BC[i].BadApvs >> 4) & 0x1) + ((BC[i].BadApvs >> 3) & 0x1) +
701  ((BC[i].BadApvs >> 2) & 0x1) + ((BC[i].BadApvs >> 1) & 0x1) + ((BC[i].BadApvs) & 0x1);
702 
703  //&&&&&&&&&&&&&&&&&
704  //Single SubSystem
705  //&&&&&&&&&&&&&&&&&
706 
707  int component;
708  SiStripDetId a(BC[i].detid);
709  if (a.subdetId() == SiStripDetId::TIB) {
710  //&&&&&&&&&&&&&&&&&
711  //TIB
712  //&&&&&&&&&&&&&&&&&
713 
714  component = tTopo.tibLayer(BC[i].detid);
715  SetBadComponents(0, component, BC[i], ssV, NBadComponent);
716 
717  } else if (a.subdetId() == SiStripDetId::TID) {
718  //&&&&&&&&&&&&&&&&&
719  //TID
720  //&&&&&&&&&&&&&&&&&
721 
722  component = tTopo.tidSide(BC[i].detid) == 2 ? tTopo.tidWheel(BC[i].detid) : tTopo.tidWheel(BC[i].detid) + 3;
723  SetBadComponents(1, component, BC[i], ssV, NBadComponent);
724 
725  } else if (a.subdetId() == SiStripDetId::TOB) {
726  //&&&&&&&&&&&&&&&&&
727  //TOB
728  //&&&&&&&&&&&&&&&&&
729 
730  component = tTopo.tobLayer(BC[i].detid);
731  SetBadComponents(2, component, BC[i], ssV, NBadComponent);
732 
733  } else if (a.subdetId() == SiStripDetId::TEC) {
734  //&&&&&&&&&&&&&&&&&
735  //TEC
736  //&&&&&&&&&&&&&&&&&
737 
738  component = tTopo.tecSide(BC[i].detid) == 2 ? tTopo.tecWheel(BC[i].detid) : tTopo.tecWheel(BC[i].detid) + 9;
739  SetBadComponents(3, component, BC[i], ssV, NBadComponent);
740  }
741  }
742 
743  //&&&&&&&&&&&&&&&&&&
744  // Single Strip Info
745  //&&&&&&&&&&&&&&&&&&
746  float percentage = 0;
747 
748  SiStripQuality::RegistryIterator rbegin = quality_->getRegistryVectorBegin();
749  SiStripQuality::RegistryIterator rend = quality_->getRegistryVectorEnd();
750 
751  for (SiStripBadStrip::RegistryIterator rp = rbegin; rp != rend; ++rp) {
752  unsigned int detid = rp->detid;
753 
754  int subdet = -999;
755  int component = -999;
756  SiStripDetId a(detid);
757  if (a.subdetId() == 3) {
758  subdet = 0;
759  component = tTopo.tibLayer(detid);
760  } else if (a.subdetId() == 4) {
761  subdet = 1;
762  component = tTopo.tidSide(detid) == 2 ? tTopo.tidWheel(detid) : tTopo.tidWheel(detid) + 3;
763  } else if (a.subdetId() == 5) {
764  subdet = 2;
765  component = tTopo.tobLayer(detid);
766  } else if (a.subdetId() == 6) {
767  subdet = 3;
768  component = tTopo.tecSide(detid) == 2 ? tTopo.tecWheel(detid) : tTopo.tecWheel(detid) + 9;
769  }
770 
771  SiStripQuality::Range sqrange =
772  SiStripQuality::Range(quality_->getDataVectorBegin() + rp->ibegin, quality_->getDataVectorBegin() + rp->iend);
773 
774  percentage = 0;
775  for (int it = 0; it < sqrange.second - sqrange.first; it++) {
776  unsigned int range = quality_->decode(*(sqrange.first + it)).range;
777  NTkBadComponent[3] += range;
778  NBadComponent[subdet][0][3] += range;
779  NBadComponent[subdet][component][3] += range;
780  percentage += range;
781  }
782  if (percentage != 0)
783  percentage /= 128. * _detInfo.getNumberOfApvsAndStripLength(detid).first;
784  if (percentage > 1)
785  edm::LogError("SiStripQualityStatistics") << "PROBLEM detid " << detid << " value " << percentage << std::endl;
786  }
787  //&&&&&&&&&&&&&&&&&&
788  // printout
789  //&&&&&&&&&&&&&&&&&&
790 
791  cout << "\n-----------------\nNew IOV starting from run " << e.id().run() << " event " << e.id().event()
792  << " lumiBlock " << e.luminosityBlock() << " time " << e.time().value() << "\n-----------------\n";
793  cout << "\n-----------------\nGlobal Info\n-----------------";
794  cout << "\nBadComponent \t Modules \tFibers "
795  "\tApvs\tStrips\n----------------------------------------------------------------";
796  cout << "\nTracker:\t\t" << NTkBadComponent[0] << "\t" << NTkBadComponent[1] << "\t" << NTkBadComponent[2] << "\t"
797  << NTkBadComponent[3];
798  cout << endl;
799  cout << "\nTIB:\t\t\t" << NBadComponent[0][0][0] << "\t" << NBadComponent[0][0][1] << "\t" << NBadComponent[0][0][2]
800  << "\t" << NBadComponent[0][0][3];
801  cout << "\nTID:\t\t\t" << NBadComponent[1][0][0] << "\t" << NBadComponent[1][0][1] << "\t" << NBadComponent[1][0][2]
802  << "\t" << NBadComponent[1][0][3];
803  cout << "\nTOB:\t\t\t" << NBadComponent[2][0][0] << "\t" << NBadComponent[2][0][1] << "\t" << NBadComponent[2][0][2]
804  << "\t" << NBadComponent[2][0][3];
805  cout << "\nTEC:\t\t\t" << NBadComponent[3][0][0] << "\t" << NBadComponent[3][0][1] << "\t" << NBadComponent[3][0][2]
806  << "\t" << NBadComponent[3][0][3];
807  cout << "\n";
808 
809  for (int i = 1; i < 5; ++i)
810  cout << "\nTIB Layer " << i << " :\t\t" << NBadComponent[0][i][0] << "\t" << NBadComponent[0][i][1] << "\t"
811  << NBadComponent[0][i][2] << "\t" << NBadComponent[0][i][3];
812  cout << "\n";
813  for (int i = 1; i < 4; ++i)
814  cout << "\nTID+ Disk " << i << " :\t\t" << NBadComponent[1][i][0] << "\t" << NBadComponent[1][i][1] << "\t"
815  << NBadComponent[1][i][2] << "\t" << NBadComponent[1][i][3];
816  for (int i = 4; i < 7; ++i)
817  cout << "\nTID- Disk " << i - 3 << " :\t\t" << NBadComponent[1][i][0] << "\t" << NBadComponent[1][i][1] << "\t"
818  << NBadComponent[1][i][2] << "\t" << NBadComponent[1][i][3];
819  cout << "\n";
820  for (int i = 1; i < 7; ++i)
821  cout << "\nTOB Layer " << i << " :\t\t" << NBadComponent[2][i][0] << "\t" << NBadComponent[2][i][1] << "\t"
822  << NBadComponent[2][i][2] << "\t" << NBadComponent[2][i][3];
823  cout << "\n";
824  for (int i = 1; i < 10; ++i)
825  cout << "\nTEC+ Disk " << i << " :\t\t" << NBadComponent[3][i][0] << "\t" << NBadComponent[3][i][1] << "\t"
826  << NBadComponent[3][i][2] << "\t" << NBadComponent[3][i][3];
827  for (int i = 10; i < 19; ++i)
828  cout << "\nTEC- Disk " << i - 9 << " :\t\t" << NBadComponent[3][i][0] << "\t" << NBadComponent[3][i][1] << "\t"
829  << NBadComponent[3][i][2] << "\t" << NBadComponent[3][i][3];
830  cout << "\n";
831 
832  cout << "\n----------------------------------------------------------------\n\t\t Detid \tModules Fibers "
833  "Apvs\n----------------------------------------------------------------";
834  for (int i = 1; i < 5; ++i)
835  cout << "\nTIB Layer " << i << " :" << ssV[0][i].str();
836  cout << "\n";
837  for (int i = 1; i < 4; ++i)
838  cout << "\nTID+ Disk " << i << " :" << ssV[1][i].str();
839  for (int i = 4; i < 7; ++i)
840  cout << "\nTID- Disk " << i - 3 << " :" << ssV[1][i].str();
841  cout << "\n";
842  for (int i = 1; i < 7; ++i)
843  cout << "\nTOB Layer " << i << " :" << ssV[2][i].str();
844  cout << "\n";
845  for (int i = 1; i < 10; ++i)
846  cout << "\nTEC+ Disk " << i << " :" << ssV[3][i].str();
847  for (int i = 10; i < 19; ++i)
848  cout << "\nTEC- Disk " << i - 9 << " :" << ssV[3][i].str();
849 
850  // store also bad modules in log file
851  ofstream badModules;
852  badModules.open("BadModules.log");
853  badModules << "\n----------------------------------------------------------------\n\t\t Detid \tModules Fibers "
854  "Apvs\n----------------------------------------------------------------";
855  for (int i = 1; i < 5; ++i)
856  badModules << "\nTIB Layer " << i << " :" << ssV[0][i].str();
857  badModules << "\n";
858  for (int i = 1; i < 4; ++i)
859  badModules << "\nTID+ Disk " << i << " :" << ssV[1][i].str();
860  for (int i = 4; i < 7; ++i)
861  badModules << "\nTID- Disk " << i - 3 << " :" << ssV[1][i].str();
862  badModules << "\n";
863  for (int i = 1; i < 7; ++i)
864  badModules << "\nTOB Layer " << i << " :" << ssV[2][i].str();
865  badModules << "\n";
866  for (int i = 1; i < 10; ++i)
867  badModules << "\nTEC+ Disk " << i << " :" << ssV[3][i].str();
868  for (int i = 10; i < 19; ++i)
869  badModules << "\nTEC- Disk " << i - 9 << " :" << ssV[3][i].str();
870  badModules.close();
871 }
872 
873 void SiStripHitEffFromCalibTree::makeHotColdMaps() {
874  cout << "Entering hot cold map generation!\n";
875  TStyle* gStyle = new TStyle("gStyle", "myStyle");
876  gStyle->cd();
877  gStyle->SetPalette(1);
878  gStyle->SetCanvasColor(kWhite);
879  gStyle->SetOptStat(0);
880  //Here we make the hot/cold color maps that we love so very much
881  //Already have access to the data as a private variable
882  //Create all of the histograms in the TFileService
883  TH2F* temph2;
884  for (Long_t maplayer = 1; maplayer <= 22; maplayer++) {
885  //Initialize all of the histograms
886  if (maplayer > 0 && maplayer <= 4) {
887  //We are in the TIB
888  temph2 = fs->make<TH2F>(Form("%s%i", "TIB", (int)(maplayer)), "TIB", 100, -1, 361, 100, -100, 100);
889  temph2->GetXaxis()->SetTitle("Phi");
890  temph2->GetXaxis()->SetBinLabel(1, TString("360"));
891  temph2->GetXaxis()->SetBinLabel(50, TString("180"));
892  temph2->GetXaxis()->SetBinLabel(100, TString("0"));
893  temph2->GetYaxis()->SetTitle("Global Z");
894  temph2->SetOption("colz");
895  HotColdMaps.push_back(temph2);
896  } else if (maplayer > 4 && maplayer <= 10) {
897  //We are in the TOB
898  temph2 = fs->make<TH2F>(Form("%s%i", "TOB", (int)(maplayer - 4)), "TOB", 100, -1, 361, 100, -120, 120);
899  temph2->GetXaxis()->SetTitle("Phi");
900  temph2->GetXaxis()->SetBinLabel(1, TString("360"));
901  temph2->GetXaxis()->SetBinLabel(50, TString("180"));
902  temph2->GetXaxis()->SetBinLabel(100, TString("0"));
903  temph2->GetYaxis()->SetTitle("Global Z");
904  temph2->SetOption("colz");
905  HotColdMaps.push_back(temph2);
906  } else if (maplayer > 10 && maplayer <= 13) {
907  //We are in the TID
908  //Split by +/-
909  temph2 = fs->make<TH2F>(Form("%s%i", "TID-", (int)(maplayer - 10)), "TID-", 100, -100, 100, 100, -100, 100);
910  temph2->GetXaxis()->SetTitle("Global Y");
911  temph2->GetXaxis()->SetBinLabel(1, TString("+Y"));
912  temph2->GetXaxis()->SetBinLabel(50, TString("0"));
913  temph2->GetXaxis()->SetBinLabel(100, TString("-Y"));
914  temph2->GetYaxis()->SetTitle("Global X");
915  temph2->GetYaxis()->SetBinLabel(1, TString("-X"));
916  temph2->GetYaxis()->SetBinLabel(50, TString("0"));
917  temph2->GetYaxis()->SetBinLabel(100, TString("+X"));
918  temph2->SetOption("colz");
919  HotColdMaps.push_back(temph2);
920  temph2 = fs->make<TH2F>(Form("%s%i", "TID+", (int)(maplayer - 10)), "TID+", 100, -100, 100, 100, -100, 100);
921  temph2->GetXaxis()->SetTitle("Global Y");
922  temph2->GetXaxis()->SetBinLabel(1, TString("+Y"));
923  temph2->GetXaxis()->SetBinLabel(50, TString("0"));
924  temph2->GetXaxis()->SetBinLabel(100, TString("-Y"));
925  temph2->GetYaxis()->SetTitle("Global X");
926  temph2->GetYaxis()->SetBinLabel(1, TString("-X"));
927  temph2->GetYaxis()->SetBinLabel(50, TString("0"));
928  temph2->GetYaxis()->SetBinLabel(100, TString("+X"));
929  temph2->SetOption("colz");
930  HotColdMaps.push_back(temph2);
931  } else if (maplayer > 13) {
932  //We are in the TEC
933  //Split by +/-
934  temph2 = fs->make<TH2F>(Form("%s%i", "TEC-", (int)(maplayer - 13)), "TEC-", 100, -120, 120, 100, -120, 120);
935  temph2->GetXaxis()->SetTitle("Global Y");
936  temph2->GetXaxis()->SetBinLabel(1, TString("+Y"));
937  temph2->GetXaxis()->SetBinLabel(50, TString("0"));
938  temph2->GetXaxis()->SetBinLabel(100, TString("-Y"));
939  temph2->GetYaxis()->SetTitle("Global X");
940  temph2->GetYaxis()->SetBinLabel(1, TString("-X"));
941  temph2->GetYaxis()->SetBinLabel(50, TString("0"));
942  temph2->GetYaxis()->SetBinLabel(100, TString("+X"));
943  temph2->SetOption("colz");
944  HotColdMaps.push_back(temph2);
945  temph2 = fs->make<TH2F>(Form("%s%i", "TEC+", (int)(maplayer - 13)), "TEC+", 100, -120, 120, 100, -120, 120);
946  temph2->GetXaxis()->SetTitle("Global Y");
947  temph2->GetXaxis()->SetBinLabel(1, TString("+Y"));
948  temph2->GetXaxis()->SetBinLabel(50, TString("0"));
949  temph2->GetXaxis()->SetBinLabel(100, TString("-Y"));
950  temph2->GetYaxis()->SetTitle("Global X");
951  temph2->GetYaxis()->SetBinLabel(1, TString("-X"));
952  temph2->GetYaxis()->SetBinLabel(50, TString("0"));
953  temph2->GetYaxis()->SetBinLabel(100, TString("+X"));
954  temph2->SetOption("colz");
955  HotColdMaps.push_back(temph2);
956  }
957  }
958  for (Long_t mylayer = 1; mylayer <= 22; mylayer++) {
959  //Determine what kind of plot we want to write out
960  //Loop through the entirety of each layer
961  //Create an array of the histograms
962  vector<hit>::const_iterator iter;
963  for (iter = hits[mylayer].begin(); iter != hits[mylayer].end(); iter++) {
964  //Looping over the particular layer
965  //Fill by 360-x to get the proper location to compare with TKMaps of phi
966  //Also global xy is messed up
967  if (mylayer > 0 && mylayer <= 4) {
968  //We are in the TIB
969  float phi = calcPhi(iter->x, iter->y);
970  HotColdMaps[mylayer - 1]->Fill(360. - phi, iter->z, 1.);
971  } else if (mylayer > 4 && mylayer <= 10) {
972  //We are in the TOB
973  float phi = calcPhi(iter->x, iter->y);
974  HotColdMaps[mylayer - 1]->Fill(360. - phi, iter->z, 1.);
975  } else if (mylayer > 10 && mylayer <= 13) {
976  //We are in the TID
977  //There are 2 different maps here
978  int side = (((iter->id) >> 13) & 0x3);
979  if (side == 1)
980  HotColdMaps[(mylayer - 1) + (mylayer - 11)]->Fill(-iter->y, iter->x, 1.);
981  else if (side == 2)
982  HotColdMaps[(mylayer - 1) + (mylayer - 10)]->Fill(-iter->y, iter->x, 1.);
983  //if(side == 1) HotColdMaps[(mylayer - 1) + (mylayer - 11)]->Fill(iter->x,iter->y,1.);
984  //else if(side == 2) HotColdMaps[(mylayer - 1) + (mylayer - 10)]->Fill(iter->x,iter->y,1.);
985  } else if (mylayer > 13) {
986  //We are in the TEC
987  //There are 2 different maps here
988  int side = (((iter->id) >> 18) & 0x3);
989  if (side == 1)
990  HotColdMaps[(mylayer + 2) + (mylayer - 14)]->Fill(-iter->y, iter->x, 1.);
991  else if (side == 2)
992  HotColdMaps[(mylayer + 2) + (mylayer - 13)]->Fill(-iter->y, iter->x, 1.);
993  //if(side == 1) HotColdMaps[(mylayer + 2) + (mylayer - 14)]->Fill(iter->x,iter->y,1.);
994  //else if(side == 2) HotColdMaps[(mylayer + 2) + (mylayer - 13)]->Fill(iter->x,iter->y,1.);
995  }
996  }
997  }
998  cout << "Finished HotCold Map Generation\n";
999 }
1000 
1001 void SiStripHitEffFromCalibTree::makeTKMap(bool autoTagging = false) {
1002  cout << "Entering TKMap generation!\n";
1003  tkmap = new TrackerMap(" Detector Inefficiency ");
1004  tkmapbad = new TrackerMap(" Inefficient Modules ");
1005  tkmapeff = new TrackerMap(_title.Data());
1006  tkmapnum = new TrackerMap(" Detector numerator ");
1007  tkmapden = new TrackerMap(" Detector denominator ");
1008 
1009  double myeff, mynum, myden;
1010  double eff_limit = 0;
1011 
1012  for (Long_t i = 1; i <= 22; i++) {
1013  //Loop over every layer, extracting the information from
1014  //the map of the efficiencies
1015  layertotal[i] = 0;
1016  layerfound[i] = 0;
1017  TH1F* hEffInLayer =
1018  fs->make<TH1F>(Form("eff_layer%i", int(i)), Form("Module efficiency in layer %i", int(i)), 201, 0, 1.005);
1019 
1020  map<unsigned int, pair<unsigned int, unsigned int> >::const_iterator ih;
1021  for (ih = modCounter[i].begin(); ih != modCounter[i].end(); ih++) {
1022  //We should be in the layer in question, and looping over all of the modules in said layer
1023  //Generate the list for the TKmap, and the bad module list
1024  mynum = (double)(((*ih).second).second);
1025  myden = (double)(((*ih).second).first);
1026  if (myden > 0)
1027  myeff = mynum / myden;
1028  else
1029  myeff = 0;
1030  hEffInLayer->Fill(myeff);
1031 
1032  if (!autoTagging) {
1033  if ((myden >= nModsMin) && (myeff < threshold)) {
1034  //We have a bad module, put it in the list!
1035  BadModules[(*ih).first] = myeff;
1036  tkmapbad->fillc((*ih).first, 255, 0, 0);
1037  cout << "Layer " << i << " (" << GetLayerName(i) << ") module " << (*ih).first << " efficiency: " << myeff
1038  << " , " << mynum << "/" << myden << endl;
1039  } else {
1040  //Fill the bad list with empty results for every module
1041  tkmapbad->fillc((*ih).first, 255, 255, 255);
1042  }
1043  if (myden < nModsMin) {
1044  cout << "Layer " << i << " (" << GetLayerName(i) << ") module " << (*ih).first << " is under occupancy at "
1045  << myden << endl;
1046  }
1047  }
1048 
1049  //Put any module into the TKMap
1050  tkmap->fill((*ih).first, 1. - myeff);
1051  tkmapeff->fill((*ih).first, myeff);
1052  tkmapnum->fill((*ih).first, mynum);
1053  tkmapden->fill((*ih).first, myden);
1054 
1055  //Add the number of hits in the layer
1056  layertotal[i] += long(myden);
1057  layerfound[i] += long(mynum);
1058  }
1059 
1060  if (autoTagging) {
1061  //Compute threshold to use for each layer
1062  hEffInLayer->GetXaxis()->SetRange(3, hEffInLayer->GetNbinsX() + 1); // Remove from the avg modules below 1%
1063  eff_limit = hEffInLayer->GetMean() - threshold;
1064  cout << "Layer " << i << " threshold for bad modules: " << eff_limit << endl;
1065  hEffInLayer->GetXaxis()->SetRange(1, hEffInLayer->GetNbinsX() + 1);
1066 
1067  for (ih = modCounter[i].begin(); ih != modCounter[i].end(); ih++) {
1068  // Second loop over modules to tag inefficient ones
1069  mynum = (double)(((*ih).second).second);
1070  myden = (double)(((*ih).second).first);
1071  if (myden > 0)
1072  myeff = mynum / myden;
1073  else
1074  myeff = 0;
1075  if ((myden >= nModsMin) && (myeff < eff_limit)) {
1076  //We have a bad module, put it in the list!
1077  BadModules[(*ih).first] = myeff;
1078  tkmapbad->fillc((*ih).first, 255, 0, 0);
1079  cout << "Layer " << i << " (" << GetLayerName(i) << ") module " << (*ih).first << " efficiency: " << myeff
1080  << " , " << mynum << "/" << myden << endl;
1081  } else {
1082  //Fill the bad list with empty results for every module
1083  tkmapbad->fillc((*ih).first, 255, 255, 255);
1084  }
1085  if (myden < nModsMin) {
1086  cout << "Layer " << i << " (" << GetLayerName(i) << ") module " << (*ih).first << " layer " << i
1087  << " is under occupancy at " << myden << endl;
1088  }
1089  }
1090  }
1091  }
1092  tkmap->save(true, 0, 0, "SiStripHitEffTKMap.png");
1093  tkmapbad->save(true, 0, 0, "SiStripHitEffTKMapBad.png");
1094  tkmapeff->save(true, _tkMapMin, 1., "SiStripHitEffTKMapEff.png");
1095  tkmapnum->save(true, 0, 0, "SiStripHitEffTKMapNum.png");
1096  tkmapden->save(true, 0, 0, "SiStripHitEffTKMapDen.png");
1097  cout << "Finished TKMap Generation\n";
1098 }
1099 
1100 void SiStripHitEffFromCalibTree::makeSQLite() {
1101  //Generate the SQLite file for use in the Database of the bad modules!
1102  cout << "Entering SQLite file generation!\n";
1103  std::vector<unsigned int> BadStripList;
1104  unsigned short NStrips;
1105  unsigned int id1;
1106  std::unique_ptr<SiStripQuality> pQuality = std::make_unique<SiStripQuality>(_detInfo);
1107  //This is the list of the bad strips, use to mask out entire APVs
1108  //Now simply go through the bad hit list and mask out things that
1109  //are bad!
1110  map<unsigned int, double>::const_iterator it;
1111  for (it = BadModules.begin(); it != BadModules.end(); it++) {
1112  //We need to figure out how many strips are in this particular module
1113  //To Mask correctly!
1114  NStrips = _detInfo.getNumberOfApvsAndStripLength((*it).first).first * 128;
1115  cout << "Number of strips module " << (*it).first << " is " << NStrips << endl;
1116  BadStripList.push_back(pQuality->encode(0, NStrips, 0));
1117  //Now compact into a single bad module
1118  id1 = (unsigned int)(*it).first;
1119  cout << "ID1 shoudl match list of modules above " << id1 << endl;
1120  quality_->compact(id1, BadStripList);
1121  SiStripQuality::Range range(BadStripList.begin(), BadStripList.end());
1122  quality_->put(id1, range);
1123  BadStripList.clear();
1124  }
1125  //Fill all the bad components now
1126  quality_->fillBadComponents();
1127 }
1128 
1129 void SiStripHitEffFromCalibTree::totalStatistics() {
1130  //Calculate the statistics by layer
1131  int totalfound = 0;
1132  int totaltotal = 0;
1133  double layereff;
1134  int subdetfound[5];
1135  int subdettotal[5];
1136 
1137  for (Long_t i = 1; i < 5; i++) {
1138  subdetfound[i] = 0;
1139  subdettotal[i] = 0;
1140  }
1141 
1142  for (Long_t i = 1; i <= 22; i++) {
1143  layereff = double(layerfound[i]) / double(layertotal[i]);
1144  cout << "Layer " << i << " (" << GetLayerName(i) << ") has total efficiency " << layereff << " " << layerfound[i]
1145  << "/" << layertotal[i] << endl;
1146  totalfound += layerfound[i];
1147  totaltotal += layertotal[i];
1148  if (i < 5) {
1149  subdetfound[1] += layerfound[i];
1150  subdettotal[1] += layertotal[i];
1151  }
1152  if (i >= 5 && i < 11) {
1153  subdetfound[2] += layerfound[i];
1154  subdettotal[2] += layertotal[i];
1155  }
1156  if (i >= 11 && i < 14) {
1157  subdetfound[3] += layerfound[i];
1158  subdettotal[3] += layertotal[i];
1159  }
1160  if (i >= 14) {
1161  subdetfound[4] += layerfound[i];
1162  subdettotal[4] += layertotal[i];
1163  }
1164  }
1165 
1166  cout << "The total efficiency is " << double(totalfound) / double(totaltotal) << endl;
1167  cout << " TIB: " << double(subdetfound[1]) / subdettotal[1] << " " << subdetfound[1] << "/" << subdettotal[1]
1168  << endl;
1169  cout << " TOB: " << double(subdetfound[2]) / subdettotal[2] << " " << subdetfound[2] << "/" << subdettotal[2]
1170  << endl;
1171  cout << " TID: " << double(subdetfound[3]) / subdettotal[3] << " " << subdetfound[3] << "/" << subdettotal[3]
1172  << endl;
1173  cout << " TEC: " << double(subdetfound[4]) / subdettotal[4] << " " << subdetfound[4] << "/" << subdettotal[4]
1174  << endl;
1175 }
1176 
1177 void SiStripHitEffFromCalibTree::makeSummary() {
1178  //setTDRStyle();
1179 
1180  int nLayers = 34;
1181  if (_showRings)
1182  nLayers = 30;
1183  if (!_showEndcapSides) {
1184  if (!_showRings)
1185  nLayers = 22;
1186  else
1187  nLayers = 20;
1188  }
1189 
1190  TH1F* found = fs->make<TH1F>("found", "found", nLayers + 1, 0, nLayers + 1);
1191  TH1F* all = fs->make<TH1F>("all", "all", nLayers + 1, 0, nLayers + 1);
1192  TH1F* found2 = fs->make<TH1F>("found2", "found2", nLayers + 1, 0, nLayers + 1);
1193  TH1F* all2 = fs->make<TH1F>("all2", "all2", nLayers + 1, 0, nLayers + 1);
1194  // first bin only to keep real data off the y axis so set to -1
1195  found->SetBinContent(0, -1);
1196  all->SetBinContent(0, 1);
1197 
1198  // new ROOT version: TGraph::Divide don't handle null or negative values
1199  for (Long_t i = 1; i < nLayers + 2; ++i) {
1200  found->SetBinContent(i, 1e-6);
1201  all->SetBinContent(i, 1);
1202  found2->SetBinContent(i, 1e-6);
1203  all2->SetBinContent(i, 1);
1204  }
1205 
1206  TCanvas* c7 = new TCanvas("c7", " test ", 10, 10, 800, 600);
1207  c7->SetFillColor(0);
1208  c7->SetGrid();
1209 
1210  int nLayers_max = nLayers + 1; // barrel+endcap
1211  if (!_showEndcapSides)
1212  nLayers_max = 11; // barrel
1213  for (Long_t i = 1; i < nLayers_max; ++i) {
1214  cout << "Fill only good modules layer " << i << ": S = " << goodlayerfound[i] << " B = " << goodlayertotal[i]
1215  << endl;
1216  if (goodlayertotal[i] > 5) {
1217  found->SetBinContent(i, goodlayerfound[i]);
1218  all->SetBinContent(i, goodlayertotal[i]);
1219  }
1220 
1221  cout << "Filling all modules layer " << i << ": S = " << alllayerfound[i] << " B = " << alllayertotal[i]
1222  << endl;
1223  if (alllayertotal[i] > 5) {
1224  found2->SetBinContent(i, alllayerfound[i]);
1225  all2->SetBinContent(i, alllayertotal[i]);
1226  }
1227  }
1228 
1229  // endcap - merging sides
1230  if (!_showEndcapSides) {
1231  for (Long_t i = 11; i < 14; ++i) { // TID disks
1232  cout << "Fill only good modules layer " << i << ": S = " << goodlayerfound[i] + goodlayerfound[i + 3]
1233  << " B = " << goodlayertotal[i] + goodlayertotal[i + 3] << endl;
1234  if (goodlayertotal[i] + goodlayertotal[i + 3] > 5) {
1235  found->SetBinContent(i, goodlayerfound[i] + goodlayerfound[i + 3]);
1236  all->SetBinContent(i, goodlayertotal[i] + goodlayertotal[i + 3]);
1237  }
1238  cout << "Filling all modules layer " << i << ": S = " << alllayerfound[i] + alllayerfound[i + 3]
1239  << " B = " << alllayertotal[i] + alllayertotal[i + 3] << endl;
1240  if (alllayertotal[i] + alllayertotal[i + 3] > 5) {
1241  found2->SetBinContent(i, alllayerfound[i] + alllayerfound[i + 3]);
1242  all2->SetBinContent(i, alllayertotal[i] + alllayertotal[i + 3]);
1243  }
1244  }
1245  for (Long_t i = 17; i < 17 + nTEClayers; ++i) { // TEC disks
1246  cout << "Fill only good modules layer " << i - 3
1247  << ": S = " << goodlayerfound[i] + goodlayerfound[i + nTEClayers]
1248  << " B = " << goodlayertotal[i] + goodlayertotal[i + nTEClayers] << endl;
1249  if (goodlayertotal[i] + goodlayertotal[i + nTEClayers] > 5) {
1250  found->SetBinContent(i - 3, goodlayerfound[i] + goodlayerfound[i + nTEClayers]);
1251  all->SetBinContent(i - 3, goodlayertotal[i] + goodlayertotal[i + nTEClayers]);
1252  }
1253  cout << "Filling all modules layer " << i - 3 << ": S = " << alllayerfound[i] + alllayerfound[i + nTEClayers]
1254  << " B = " << alllayertotal[i] + alllayertotal[i + nTEClayers] << endl;
1255  if (alllayertotal[i] + alllayertotal[i + nTEClayers] > 5) {
1256  found2->SetBinContent(i - 3, alllayerfound[i] + alllayerfound[i + nTEClayers]);
1257  all2->SetBinContent(i - 3, alllayertotal[i] + alllayertotal[i + nTEClayers]);
1258  }
1259  }
1260  }
1261 
1262  found->Sumw2();
1263  all->Sumw2();
1264 
1265  found2->Sumw2();
1266  all2->Sumw2();
1267 
1268  TGraphAsymmErrors* gr = fs->make<TGraphAsymmErrors>(nLayers + 1);
1269  gr->SetName("eff_good");
1270  gr->BayesDivide(found, all);
1271 
1272  TGraphAsymmErrors* gr2 = fs->make<TGraphAsymmErrors>(nLayers + 1);
1273  gr2->SetName("eff_all");
1274  gr2->BayesDivide(found2, all2);
1275 
1276  for (int j = 0; j < nLayers + 1; j++) {
1277  gr->SetPointError(j, 0., 0., gr->GetErrorYlow(j), gr->GetErrorYhigh(j));
1278  gr2->SetPointError(j, 0., 0., gr2->GetErrorYlow(j), gr2->GetErrorYhigh(j));
1279  }
1280 
1281  gr->GetXaxis()->SetLimits(0, nLayers);
1282  gr->SetMarkerColor(2);
1283  gr->SetMarkerSize(1.2);
1284  gr->SetLineColor(2);
1285  gr->SetLineWidth(4);
1286  gr->SetMarkerStyle(20);
1287  gr->SetMinimum(_effPlotMin);
1288  gr->SetMaximum(1.001);
1289  gr->GetYaxis()->SetTitle("Efficiency");
1290  gStyle->SetTitleFillColor(0);
1291  gStyle->SetTitleBorderSize(0);
1292  gr->SetTitle(_title);
1293 
1294  gr2->GetXaxis()->SetLimits(0, nLayers);
1295  gr2->SetMarkerColor(1);
1296  gr2->SetMarkerSize(1.2);
1297  gr2->SetLineColor(1);
1298  gr2->SetLineWidth(4);
1299  gr2->SetMarkerStyle(21);
1300  gr2->SetMinimum(_effPlotMin);
1301  gr2->SetMaximum(1.001);
1302  gr2->GetYaxis()->SetTitle("Efficiency");
1303  gr2->SetTitle(_title);
1304 
1305  for (Long_t k = 1; k < nLayers + 1; k++) {
1306  TString label;
1307  if (_showEndcapSides)
1308  label = GetLayerSideName(k);
1309  else
1310  label = GetLayerName(k);
1311  if (!_showTOB6TEC9) {
1312  if (k == 10)
1313  label = "";
1314  if (!_showRings && k == nLayers)
1315  label = "";
1316  if (!_showRings && _showEndcapSides && k == 25)
1317  label = "";
1318  }
1319  if (!_showRings) {
1320  if (_showEndcapSides) {
1321  gr->GetXaxis()->SetBinLabel(((k + 1) * 100 + 2) / (nLayers)-4, label);
1322  gr2->GetXaxis()->SetBinLabel(((k + 1) * 100 + 2) / (nLayers)-4, label);
1323  } else {
1324  gr->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-6, label);
1325  gr2->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-6, label);
1326  }
1327  } else {
1328  if (_showEndcapSides) {
1329  gr->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-4, label);
1330  gr2->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-4, label);
1331  } else {
1332  gr->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-7, label);
1333  gr2->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-7, label);
1334  }
1335  }
1336  }
1337 
1338  gr->Draw("AP");
1339  gr->GetXaxis()->SetNdivisions(36);
1340 
1341  c7->cd();
1342  TPad* overlay = new TPad("overlay", "", 0, 0, 1, 1);
1343  overlay->SetFillStyle(4000);
1344  overlay->SetFillColor(0);
1345  overlay->SetFrameFillStyle(4000);
1346  overlay->Draw("same");
1347  overlay->cd();
1348  if (!_showOnlyGoodModules)
1349  gr2->Draw("AP");
1350 
1351  TLegend* leg = new TLegend(0.70, 0.27, 0.88, 0.40);
1352  leg->AddEntry(gr, "Good Modules", "p");
1353  if (!_showOnlyGoodModules)
1354  leg->AddEntry(gr2, "All Modules", "p");
1355  leg->SetTextSize(0.020);
1356  leg->SetFillColor(0);
1357  leg->Draw("same");
1358 
1359  c7->SaveAs("Summary.png");
1360 }
1361 
1362 void SiStripHitEffFromCalibTree::makeSummaryVsBx() {
1363  cout << "Computing efficiency vs bx" << endl;
1364 
1365  unsigned int nLayers = 22;
1366  if (_showRings)
1367  nLayers = 20;
1368 
1369  for (unsigned int ilayer = 1; ilayer < nLayers; ilayer++) {
1370  TH1F* hfound = fs->make<TH1F>(Form("foundVsBx_layer%i", ilayer), Form("layer %i", ilayer), 3565, 0, 3565);
1371  TH1F* htotal = fs->make<TH1F>(Form("totalVsBx_layer%i", ilayer), Form("layer %i", ilayer), 3565, 0, 3565);
1372 
1373  for (unsigned int ibx = 0; ibx < 3566; ibx++) {
1374  hfound->SetBinContent(ibx, 1e-6);
1375  htotal->SetBinContent(ibx, 1);
1376  }
1377  map<unsigned int, vector<int> >::iterator iterMapvsBx;
1378  for (iterMapvsBx = layerfound_perBx.begin(); iterMapvsBx != layerfound_perBx.end(); ++iterMapvsBx)
1379  hfound->SetBinContent(iterMapvsBx->first, iterMapvsBx->second[ilayer]);
1380  for (iterMapvsBx = layertotal_perBx.begin(); iterMapvsBx != layertotal_perBx.end(); ++iterMapvsBx)
1381  if (iterMapvsBx->second[ilayer] > 0)
1382  htotal->SetBinContent(iterMapvsBx->first, iterMapvsBx->second[ilayer]);
1383 
1384  hfound->Sumw2();
1385  htotal->Sumw2();
1386 
1387  TGraphAsymmErrors* geff = fs->make<TGraphAsymmErrors>(3564);
1388  geff->SetName(Form("effVsBx_layer%i", ilayer));
1389  geff->SetTitle("Hit Efficiency vs bx - " + GetLayerName(ilayer));
1390  geff->BayesDivide(hfound, htotal);
1391 
1392  //Average over trains
1393  TGraphAsymmErrors* geff_avg = fs->make<TGraphAsymmErrors>();
1394  geff_avg->SetName(Form("effVsBxAvg_layer%i", ilayer));
1395  geff_avg->SetTitle("Hit Efficiency vs bx - " + GetLayerName(ilayer));
1396  geff_avg->SetMarkerStyle(20);
1397  int ibx = 0;
1398  int previous_bx = -80;
1399  int delta_bx = 0;
1400  int nbx = 0;
1401  int found = 0;
1402  int total = 0;
1403  double sum_bx = 0;
1404  int ipt = 0;
1405  float low, up, eff;
1406  int firstbx = 0;
1407  for (iterMapvsBx = layertotal_perBx.begin(); iterMapvsBx != layertotal_perBx.end(); ++iterMapvsBx) {
1408  ibx = iterMapvsBx->first;
1409  delta_bx = ibx - previous_bx;
1410  // consider a new train
1411  if (delta_bx > (int)_spaceBetweenTrains && nbx > 0 && total > 0) {
1412  eff = found / (float)total;
1413  //cout<<"new train "<<ipt<<" "<<sum_bx/nbx<<" "<<eff<<endl;
1414  geff_avg->SetPoint(ipt, sum_bx / nbx, eff);
1415  low = TEfficiency::Bayesian(total, found, .683, 1, 1, false);
1416  up = TEfficiency::Bayesian(total, found, .683, 1, 1, true);
1417  geff_avg->SetPointError(ipt, sum_bx / nbx - firstbx, previous_bx - sum_bx / nbx, eff - low, up - eff);
1418  ipt++;
1419  sum_bx = 0;
1420  found = 0;
1421  total = 0;
1422  nbx = 0;
1423  firstbx = ibx;
1424  }
1425  sum_bx += ibx;
1426  found += hfound->GetBinContent(ibx);
1427  total += htotal->GetBinContent(ibx);
1428  nbx++;
1429 
1430  previous_bx = ibx;
1431  }
1432  //last train
1433  eff = found / (float)total;
1434  //cout<<"new train "<<ipt<<" "<<sum_bx/nbx<<" "<<eff<<endl;
1435  geff_avg->SetPoint(ipt, sum_bx / nbx, eff);
1436  low = TEfficiency::Bayesian(total, found, .683, 1, 1, false);
1437  up = TEfficiency::Bayesian(total, found, .683, 1, 1, true);
1438  geff_avg->SetPointError(ipt, sum_bx / nbx - firstbx, previous_bx - sum_bx / nbx, eff - low, up - eff);
1439  }
1440 }
1441 
1442 TString SiStripHitEffFromCalibTree::GetLayerName(Long_t k) {
1443  TString layername = "";
1444  TString ringlabel = "D";
1445  if (_showRings)
1446  ringlabel = "R";
1447  if (k > 0 && k < 5) {
1448  layername = TString("TIB L") + k;
1449  } else if (k > 4 && k < 11) {
1450  layername = TString("TOB L") + (k - 4);
1451  } else if (k > 10 && k < 14) {
1452  layername = TString("TID ") + ringlabel + (k - 10);
1453  } else if (k > 13 && k < 14 + nTEClayers) {
1454  layername = TString("TEC ") + ringlabel + (k - 13);
1455  }
1456 
1457  return layername;
1458 }
1459 
1460 void SiStripHitEffFromCalibTree::ComputeEff(vector<TH1F*>& vhfound, vector<TH1F*>& vhtotal, string name) {
1461  unsigned int nLayers = 22;
1462  if (_showRings)
1463  nLayers = 20;
1464 
1465  TH1F* hfound;
1466  TH1F* htotal;
1467 
1468  for (unsigned int ilayer = 1; ilayer < nLayers; ilayer++) {
1469  hfound = vhfound[ilayer];
1470  htotal = vhtotal[ilayer];
1471 
1472  hfound->Sumw2();
1473  htotal->Sumw2();
1474 
1475  // new ROOT version: TGraph::Divide don't handle null or negative values
1476  for (Long_t i = 0; i < hfound->GetNbinsX() + 1; ++i) {
1477  if (hfound->GetBinContent(i) == 0)
1478  hfound->SetBinContent(i, 1e-6);
1479  if (htotal->GetBinContent(i) == 0)
1480  htotal->SetBinContent(i, 1);
1481  }
1482 
1483  TGraphAsymmErrors* geff = fs->make<TGraphAsymmErrors>(hfound->GetNbinsX());
1484  geff->SetName(Form("%s_layer%i", name.c_str(), ilayer));
1485  geff->BayesDivide(hfound, htotal);
1486  if (name == "effVsLumi")
1487  geff->SetTitle("Hit Efficiency vs inst. lumi. - " + GetLayerName(ilayer));
1488  if (name == "effVsPU")
1489  geff->SetTitle("Hit Efficiency vs pileup - " + GetLayerName(ilayer));
1490  if (name == "effVsCM")
1491  geff->SetTitle("Hit Efficiency vs common Mode - " + GetLayerName(ilayer));
1492  geff->SetMarkerStyle(20);
1493  }
1494 }
1495 
1496 void SiStripHitEffFromCalibTree::makeSummaryVsLumi() {
1497  cout << "Computing efficiency vs lumi" << endl;
1498 
1499  if (instLumiHisto->GetEntries()) // from infos per event
1500  cout << "Avg conditions (avg+/-rms): lumi :" << instLumiHisto->GetMean() << "+/-" << instLumiHisto->GetRMS()
1501  << " pu: " << PUHisto->GetMean() << "+/-" << PUHisto->GetRMS() << endl;
1502 
1503  else { // from infos per hit
1504 
1505  unsigned int nLayers = 22;
1506  if (_showRings)
1507  nLayers = 20;
1508  unsigned int nLayersForAvg = 0;
1509  float layerLumi = 0;
1510  float layerPU = 0;
1511  float avgLumi = 0;
1512  float avgPU = 0;
1513 
1514  cout << "Lumi summary: (avg over trajectory measurements)" << endl;
1515  for (unsigned int ilayer = 1; ilayer < nLayers; ilayer++) {
1516  layerLumi = layertotal_vsLumi[ilayer]->GetMean();
1517  layerPU = layertotal_vsPU[ilayer]->GetMean();
1518  //cout<<" layer "<<ilayer<<" lumi: "<<layerLumi<<" pu: "<<layerPU<<endl;
1519  if (layerLumi != 0 && layerPU != 0) {
1520  avgLumi += layerLumi;
1521  avgPU += layerPU;
1522  nLayersForAvg++;
1523  }
1524  }
1525  avgLumi /= nLayersForAvg;
1526  avgPU /= nLayersForAvg;
1527  cout << "Avg conditions: lumi :" << avgLumi << " pu: " << avgPU << endl;
1528  }
1529 
1530  ComputeEff(layerfound_vsLumi, layertotal_vsLumi, "effVsLumi");
1531  ComputeEff(layerfound_vsPU, layertotal_vsPU, "effVsPU");
1532 }
1533 
1534 void SiStripHitEffFromCalibTree::makeSummaryVsCM() {
1535  cout << "Computing efficiency vs CM" << endl;
1536  ComputeEff(layerfound_vsCM, layertotal_vsCM, "effVsCM");
1537 }
1538 
1539 TString SiStripHitEffFromCalibTree::GetLayerSideName(Long_t k) {
1540  TString layername = "";
1541  TString ringlabel = "D";
1542  if (_showRings)
1543  ringlabel = "R";
1544  if (k > 0 && k < 5) {
1545  layername = TString("TIB L") + k;
1546  } else if (k > 4 && k < 11) {
1547  layername = TString("TOB L") + (k - 4);
1548  } else if (k > 10 && k < 14) {
1549  layername = TString("TID- ") + ringlabel + (k - 10);
1550  } else if (k > 13 && k < 17) {
1551  layername = TString("TID+ ") + ringlabel + (k - 13);
1552  } else if (k > 16 && k < 17 + nTEClayers) {
1553  layername = TString("TEC- ") + ringlabel + (k - 16);
1554  } else if (k > 16 + nTEClayers) {
1555  layername = TString("TEC+ ") + ringlabel + (k - 16 - nTEClayers);
1556  }
1557 
1558  return layername;
1559 }
1560 
1561 std::unique_ptr<SiStripBadStrip> SiStripHitEffFromCalibTree::getNewObject() {
1562  //Need this for a Condition DB Writer
1563  //Initialize a return variable
1564  auto obj = std::make_unique<SiStripBadStrip>();
1565 
1566  SiStripBadStrip::RegistryIterator rIter = quality_->getRegistryVectorBegin();
1567  SiStripBadStrip::RegistryIterator rIterEnd = quality_->getRegistryVectorEnd();
1568 
1569  for (; rIter != rIterEnd; ++rIter) {
1570  SiStripBadStrip::Range range(quality_->getDataVectorBegin() + rIter->ibegin,
1571  quality_->getDataVectorBegin() + rIter->iend);
1572  if (!obj->put(rIter->detid, range))
1573  edm::LogError("SiStripHitEffFromCalibTree")
1574  << "[SiStripHitEffFromCalibTree::getNewObject] detid already exists" << std::endl;
1575  }
1576 
1577  return obj;
1578 }
1579 
1580 float SiStripHitEffFromCalibTree::calcPhi(float x, float y) {
1581  float phi = 0;
1582  float Pi = 3.14159;
1583  if ((x >= 0) && (y >= 0))
1584  phi = atan(y / x);
1585  else if ((x >= 0) && (y <= 0))
1586  phi = atan(y / x) + 2 * Pi;
1587  else if ((x <= 0) && (y >= 0))
1588  phi = atan(y / x) + Pi;
1589  else
1590  phi = atan(y / x) + Pi;
1591  phi = phi * 180.0 / Pi;
1592 
1593  return phi;
1594 }
1595 
1596 void SiStripHitEffFromCalibTree::SetBadComponents(
1597  int i, int component, SiStripQuality::BadComponent& BC, std::stringstream ssV[4][19], int NBadComponent[4][19][4]) {
1598  int napv = _detInfo.getNumberOfApvsAndStripLength(BC.detid).first;
1599 
1600  ssV[i][component] << "\n\t\t " << BC.detid << " \t " << BC.BadModule << " \t " << ((BC.BadFibers) & 0x1) << " ";
1601  if (napv == 4)
1602  ssV[i][component] << "x " << ((BC.BadFibers >> 1) & 0x1);
1603 
1604  if (napv == 6)
1605  ssV[i][component] << ((BC.BadFibers >> 1) & 0x1) << " " << ((BC.BadFibers >> 2) & 0x1);
1606  ssV[i][component] << " \t " << ((BC.BadApvs) & 0x1) << " " << ((BC.BadApvs >> 1) & 0x1) << " ";
1607  if (napv == 4)
1608  ssV[i][component] << "x x " << ((BC.BadApvs >> 2) & 0x1) << " " << ((BC.BadApvs >> 3) & 0x1);
1609  if (napv == 6)
1610  ssV[i][component] << ((BC.BadApvs >> 2) & 0x1) << " " << ((BC.BadApvs >> 3) & 0x1) << " "
1611  << ((BC.BadApvs >> 4) & 0x1) << " " << ((BC.BadApvs >> 5) & 0x1) << " ";
1612 
1613  if (BC.BadApvs) {
1614  NBadComponent[i][0][2] += ((BC.BadApvs >> 5) & 0x1) + ((BC.BadApvs >> 4) & 0x1) + ((BC.BadApvs >> 3) & 0x1) +
1615  ((BC.BadApvs >> 2) & 0x1) + ((BC.BadApvs >> 1) & 0x1) + ((BC.BadApvs) & 0x1);
1616  NBadComponent[i][component][2] += ((BC.BadApvs >> 5) & 0x1) + ((BC.BadApvs >> 4) & 0x1) +
1617  ((BC.BadApvs >> 3) & 0x1) + ((BC.BadApvs >> 2) & 0x1) +
1618  ((BC.BadApvs >> 1) & 0x1) + ((BC.BadApvs) & 0x1);
1619  }
1620  if (BC.BadFibers) {
1621  NBadComponent[i][0][1] += ((BC.BadFibers >> 2) & 0x1) + ((BC.BadFibers >> 1) & 0x1) + ((BC.BadFibers) & 0x1);
1622  NBadComponent[i][component][1] +=
1623  ((BC.BadFibers >> 2) & 0x1) + ((BC.BadFibers >> 1) & 0x1) + ((BC.BadFibers) & 0x1);
1624  }
1625  if (BC.BadModule) {
1626  NBadComponent[i][0][0]++;
1627  NBadComponent[i][component][0]++;
1628  }
1629 }
1630 
1631 DEFINE_FWK_MODULE(SiStripHitEffFromCalibTree);
SiStripHitEffFromCalibTree::BadModules
map< unsigned int, double > BadModules
Definition: SiStripHitEffFromCalibTree.cc:184
SiStripBadStrip::data::range
unsigned short range
Definition: SiStripBadStrip.h:34
up
Definition: BitonicSort.h:7
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
compare.overlay
def overlay
Definition: compare.py:122
SiStripHitEffFromCalibTree::layertotal_perBx
map< unsigned int, vector< int > > layertotal_perBx
Definition: SiStripHitEffFromCalibTree.cc:173
StripTopology::nstrips
virtual int nstrips() const =0
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
c
const edm::EventSetup & c
Definition: SiStripLAProfileBooker.cc:66
SiStripHitEffFromCalibTree::getNewObject
std::unique_ptr< SiStripBadStrip > getNewObject() override
Definition: SiStripHitEffFromCalibTree.cc:1561
mps_fire.i
i
Definition: mps_fire.py:428
SiStripQuality::getBadComponentList
const std::vector< BadComponent > & getBadComponentList() const
Definition: SiStripQuality.h:114
gpuClustering::id
uint16_t *__restrict__ id
Definition: gpuClusterChargeCut.h:20
SiStripHitEffFromCalibTree::_tTopoToken
edm::ESGetToken< TrackerTopology, TrackerTopologyRcd > _tTopoToken
Definition: SiStripHitEffFromCalibTree.cc:151
SiStripHitEffFromCalibTree::layerfound_perBx
map< unsigned int, vector< int > > layerfound_perBx
Definition: SiStripHitEffFromCalibTree.cc:172
SiStripDetId::TID
static constexpr auto TID
Definition: SiStripDetId.h:38
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
quality
uint32_t const *__restrict__ Quality * quality
Definition: CAHitNtupletGeneratorKernelsImpl.h:122
Surface::bounds
const Bounds & bounds() const
Definition: Surface.h:87
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
relval_steps.PU
dictionary PU
Definition: relval_steps.py:1605
SiStripHitEffFromCalibTree::SetBadComponents
void SetBadComponents(int i, int component, SiStripQuality::BadComponent &BC, std::stringstream ssV[4][19], int NBadComponent[4][19][4])
Definition: SiStripHitEffFromCalibTree.cc:1596
SiStripBadStrip::RegistryIterator
Registry::const_iterator RegistryIterator
Definition: SiStripBadStrip.h:55
SiStripHitEffFromCalibTree::eventInfos
map< pair< unsigned int, unsigned int >, array< double, 3 > > eventInfos
Definition: SiStripHitEffFromCalibTree.cc:160
edm::LogError
Log< level::Error, false > LogError
Definition: MessageLogger.h:123
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:111
accept
bool accept(const edm::Event &event, const edm::TriggerResults &triggerTable, const std::string &triggerPath)
Definition: TopDQMHelpers.h:31
StripGeomDetUnit::specificTopology
virtual const StripTopology & specificTopology() const
Returns a reference to the strip proxy topology.
Definition: StripGeomDetUnit.cc:17
GeomDet::surface
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:37
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
phase1PixelTopology::layer
constexpr std::array< uint8_t, layerIndexSize > layer
Definition: phase1PixelTopology.h:110
sistrip::SpyUtilities::range
const uint16_t range(const Frame &aFrame)
Definition: SiStripSpyUtilities.cc:76
edm::EventSetup::getData
bool getData(T &iHolder) const
Definition: EventSetup.h:128
AlCaHarvesting_cff.SiStripQuality
tuple SiStripQuality
Definition: AlCaHarvesting_cff.py:254
label
char const * label
Definition: PFTauDecayModeTools.cc:11
SiStripBadStrip::getRegistryVectorEnd
RegistryIterator getRegistryVectorEnd() const
Definition: SiStripBadStrip.h:76
SiStripHitEffFromCalibTree::fs
edm::Service< TFileService > fs
Definition: SiStripHitEffFromCalibTree.cc:121
SiStripHitEffFromCalibTree::ComputeEff
void ComputeEff(vector< TH1F * > &vhfound, vector< TH1F * > &vhtotal, string name)
Definition: SiStripHitEffFromCalibTree.cc:1460
SiStripHitEffFromCalibTree::algoBeginJob
void algoBeginJob(const edm::EventSetup &) override
Definition: SiStripHitEffFromCalibTree.cc:223
TrackerMap::save
void save(bool print_total=true, float minval=0., float maxval=0., std::string s="svgmap.svg", int width=1500, int height=800)
Definition: TrackerMap.cc:811
SiStripQuality::compact
void compact(unsigned int &, std::vector< unsigned int > &)
Definition: SiStripQuality.cc:285
DetId::subdetId
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
compare_using_db.ifile
tuple ifile
Definition: compare_using_db.py:251
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
SiStripHitEffFromCalibTree::_tkGeomToken
edm::ESGetToken< TrackerGeometry, TrackerDigiGeometryRecord > _tkGeomToken
Definition: SiStripHitEffFromCalibTree.cc:150
SiStripDetInfo::getNumberOfApvsAndStripLength
const std::pair< unsigned short, double > getNumberOfApvsAndStripLength(uint32_t detId) const
Definition: SiStripDetInfo.cc:19
SiStripDetInfoFileReader::read
SiStripDetInfo read(std::string filePath)
Definition: SiStripDetInfoFileReader.cc:11
SiStripQuality::fillBadComponents
void fillBadComponents()
Definition: SiStripQuality.cc:490
SiStripDetId::TOB
static constexpr auto TOB
Definition: SiStripDetId.h:39
TrackerMap::fillc
void fillc(int idmod, int RGBcode)
Definition: TrackerMap.h:135
python.cmstools.all
def all
workaround iterator generators for ROOT classes
Definition: cmstools.py:25
hit::id
unsigned int id
Definition: SiStripHitEffFromCalibTree.cc:91
SiStripBadStrip::getDataVectorBegin
ContainerIterator getDataVectorBegin() const
Definition: SiStripBadStrip.h:73
mps_splice.mods
tuple mods
Definition: mps_splice.py:43
SiStripDetId
Detector identifier class for the strip tracker.
Definition: SiStripDetId.h:18
DetId
Definition: DetId.h:17
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
SiStripDetId::TIB
static constexpr auto TIB
Definition: SiStripDetId.h:37
SiStripBadStrip::getRegistryVectorBegin
RegistryIterator getRegistryVectorBegin() const
Definition: SiStripBadStrip.h:75
SiStripHitEffFromCalibTree::algoAnalyze
void algoAnalyze(const edm::Event &e, const edm::EventSetup &c) override
Definition: SiStripHitEffFromCalibTree.cc:232
a
double a
Definition: hdecay.h:119
SiStripBadStrip::Range
std::pair< ContainerIterator, ContainerIterator > Range
Definition: SiStripBadStrip.h:53
dataset.end
string end
Definition: dataset.py:937
edm::FileInPath::fullPath
std::string fullPath() const
Definition: FileInPath.cc:161
SiStripDetInfoFileReader::kDefaultFile
static constexpr char const *const kDefaultFile
Definition: SiStripDetInfoFileReader.h:10
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:144
SiStripBadStrip::put
bool put(const uint32_t &detID, const InputVector &vect)
Definition: SiStripBadStrip.h:65
str
#define str(s)
Definition: TestProcessor.cc:53
Bounds::width
virtual float width() const =0
SiStripHitEffFromCalibTree::SiStripHitEffFromCalibTree
SiStripHitEffFromCalibTree(const edm::ParameterSet &)
Definition: SiStripHitEffFromCalibTree.cc:187
TrackerMap::fill
void fill(int layer, int ring, int nmod, float x)
Definition: TrackerMap.cc:3289
SiStripDetId::TEC
static constexpr auto TEC
Definition: SiStripDetId.h:40
SiStripBadStrip::decode
data decode(const unsigned int &value) const
Definition: SiStripBadStrip.h:78
SiStripHitEffFromCalibTree::modCounter
map< unsigned int, pair< unsigned int, unsigned int > > modCounter[23]
Definition: SiStripHitEffFromCalibTree.cc:164
Generated for CMSSW Reference Manual by doxygen 1.8.5