CMS 3D CMS Logo

TrackerOfflineValidationSummary.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: TrackerOfflineValidationSummary
4 // Class: TrackerOfflineValidationSummary
5 //
13 //
14 // Original Author: Johannes Hauk
15 // Created: Sat Aug 22 10:31:34 CEST 2009
16 // $Id: TrackerOfflineValidationSummary.cc,v 1.7 2012/06/30 09:09:47 eulisse Exp $
17 //
18 //
19 
20 // system include files
21 #include <memory>
22 
23 // user include files
26 
29 
31 
36 
37 #include "TTree.h"
38 //#include "Alignment/OfflineValidation/interface/TrackerValidationVariables.h"
40 
48 #include "TH1.h"
49 #include "TMath.h"
50 
53 
54 //
55 // class decleration
56 //
57 
58 class TrackerOfflineValidationSummary : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
59 public:
62  ~TrackerOfflineValidationSummary() override = default;
63 
64  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
65 
66 private:
69 
70  struct ModuleHistos {
73  TH1* ResHisto;
79  };
80 
84  TH1* DmrXprime;
85  TH1* DmrYprime;
86  };
87 
90  HarvestingHierarchy(const std::string& name, const std::string& component, const std::vector<unsigned int>& entries)
91  : hierarchyName(name), componentName(component), treeEntries(entries) {}
92  // Needed for naming of histogram according to residual histograms
94  // "Pixel" or "Strip"
96  // Modules belonging to selected hierarchy
97  std::vector<unsigned int> treeEntries;
99  };
100 
101  void beginRun(const edm::Run&, const edm::EventSetup& iSetup) override{};
102  void analyze(const edm::Event& evt, const edm::EventSetup&) override;
103  void endRun(const edm::Run&, const edm::EventSetup& iSetup) override;
104  void endJob() override;
105 
106  void fillTree(TTree& tree,
107  std::map<int, TrackerOfflineValidationSummary::ModuleHistos>& moduleHist,
108  TkOffTreeVariables& treeMem,
109  const TrackerGeometry& tkgeom,
110  std::map<std::string, std::string>& substructureName,
111  const TrackerTopology* tTopo);
112 
113  std::pair<float, float> fitResiduals(TH1* hist) const;
114  float getMedian(const TH1* hist) const;
115 
118  std::map<std::string, std::string>& substructureName);
119 
120  void collateHarvestingHists(TTree& tree);
121  void applyHarvestingHierarchy(TTree& treeMem);
122  void bookHarvestingHists();
123  void getBinning(const std::string& binningPSetName, int& nBinsX, double& lowerBoundX, double& upperBoundX) const;
124  void fillHarvestingHists(TTree& tree);
125 
126  // ----------member data ---------------------------
127 
130  std::unique_ptr<TrackerTopology> tTopo_;
131 
132  // parameters from cfg to steer
134  const bool useFit_;
135 
137 
139 
140  std::map<int, TrackerOfflineValidationSummary::ModuleHistos> mPxbResiduals_;
141  std::map<int, TrackerOfflineValidationSummary::ModuleHistos> mPxeResiduals_;
142  std::map<int, TrackerOfflineValidationSummary::ModuleHistos> mTibResiduals_;
143  std::map<int, TrackerOfflineValidationSummary::ModuleHistos> mTidResiduals_;
144  std::map<int, TrackerOfflineValidationSummary::ModuleHistos> mTobResiduals_;
145  std::map<int, TrackerOfflineValidationSummary::ModuleHistos> mTecResiduals_;
146 
147  std::vector<HarvestingHierarchy> vHarvestingHierarchy_;
148 };
149 
150 //
151 // constructors and destructor
152 //
154  : geomToken_(esConsumes()),
155  tTopoToken_(esConsumes<edm::Transition::EndRun>()),
156  parSet_(iConfig),
157  tTopo_(nullptr),
158  moduleDirectory_(parSet_.getParameter<std::string>("moduleDirectoryInOutput")),
159  useFit_(parSet_.getParameter<bool>("useFit")),
160  dbe_(nullptr),
161  moduleMapsInitialized(false) {
162  //now do what ever initialization is needed
164 }
165 
168  desc.setComment("Summary of alignment payloads validation by evaluating unbiased track-hit resisuals");
169  desc.add<std::string>("moduleDirectoryInOutput", "Alignment/Tracker");
170  desc.add<bool>("useFit", false);
171  desc.add<bool>("stripYDmrs", false);
172  desc.add<unsigned int>("minEntriesPerModuleForDmr", 100);
173 
174  // fill in the residuals details
175  const std::vector<std::string> listOfDMRPSets = {
176  "TH1DmrXprimeStripModules", "TH1DmrYprimeStripModules", "TH1DmrXprimePixelModules", "TH1DmrYprimePixelModules"};
177 
178  for (const auto& myPSetName : listOfDMRPSets) {
180  myPSet.add<int>("Nbinx", 100);
181  myPSet.add<double>("xmin", -5.f);
182  myPSet.add<double>("xmax", 5.f);
183  desc.add<edm::ParameterSetDescription>(myPSetName, myPSet);
184  }
185  descriptions.addWithDefaultLabel(desc);
186 }
187 
188 //
189 // member functions
190 //
191 
192 // ------------ method called to for each event ------------
194  // Access of EventSetup is needed to get the list of silicon-modules and their IDs
195  // Since they do not change, it is accessed only once
197  return;
198  tkGeom_ = iSetup.getHandle(geomToken_);
199  const TrackerGeometry* bareTkGeomPtr = &(*tkGeom_);
200 
201  const TrackingGeometry::DetIdContainer& detIdContainer = bareTkGeomPtr->detIds();
202  std::vector<DetId>::const_iterator iDet;
203  for (iDet = detIdContainer.begin(); iDet != detIdContainer.end(); ++iDet) {
204  const DetId& detId = *iDet;
205  const uint32_t rawId = detId.rawId();
206  const unsigned int subdetId = detId.subdetId();
207  if (subdetId == PixelSubdetector::PixelBarrel)
208  mPxbResiduals_[rawId];
209  else if (subdetId == PixelSubdetector::PixelEndcap)
210  mPxeResiduals_[rawId];
211  else if (subdetId == StripSubdetector::TIB)
212  mTibResiduals_[rawId];
213  else if (subdetId == StripSubdetector::TID)
214  mTidResiduals_[rawId];
215  else if (subdetId == StripSubdetector::TOB)
216  mTobResiduals_[rawId];
217  else if (subdetId == StripSubdetector::TEC)
218  mTecResiduals_[rawId];
219  else {
220  throw cms::Exception("Geometry Error")
221  << "[TrackerOfflineValidationSummary] Error, tried to get reference for non-tracker subdet " << subdetId
222  << " from detector " << detId.det();
223  }
224  }
225  moduleMapsInitialized = true;
226 }
227 
228 // ------------ method called at each end of Run ------------
230  if (!tTopo_) {
231  tTopo_ = std::make_unique<TrackerTopology>(iSetup.getData(tTopoToken_));
232  }
233 }
234 
235 // ------------ method called once each job just after ending the event loop ------------
237  AlignableTracker aliTracker(&(*tkGeom_), tTopo_.get());
238 
239  TTree* tree = new TTree("TkOffVal", "TkOffVal");
240 
241  TkOffTreeVariables* treeMemPtr = new TkOffTreeVariables;
242  // We create branches for all members of 'TkOffTreeVariables' (even if not needed).
243  // This works because we have a dictionary for 'TkOffTreeVariables'
244  // (see src/classes_def.xml and src/classes.h):
245  tree->Branch("TkOffTreeVariables", &treeMemPtr); // address of pointer!
246  // second branch needed for assigning names and titles to harvesting histograms consistent to others
247  std::map<std::string, std::string>* substructureName = new std::map<std::string, std::string>;
248  tree->Branch("SubstructureName", &substructureName, 32000, 00); // SplitLevel must be set to zero
249 
250  this->fillTree(*tree, mPxbResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo_.get());
251  this->fillTree(*tree, mPxeResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo_.get());
252  this->fillTree(*tree, mTibResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo_.get());
253  this->fillTree(*tree, mTidResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo_.get());
254  this->fillTree(*tree, mTobResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo_.get());
255  this->fillTree(*tree, mTecResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo_.get());
256 
257  //dbe_->save("dqmOut.root");
258 
259  // Method for filling histograms which show summarized values (mean, rms, median ...)
260  // of the module-based histograms from TrackerOfflineValidation
261  this->collateHarvestingHists(*tree);
262 
263  delete tree;
264  tree = nullptr;
265  delete treeMemPtr;
266  treeMemPtr = nullptr;
267  delete substructureName;
268  substructureName = nullptr;
269 }
270 
272  std::map<int, TrackerOfflineValidationSummary::ModuleHistos>& moduleHist,
273  TkOffTreeVariables& treeMem,
274  const TrackerGeometry& tkgeom,
275  std::map<std::string, std::string>& substructureName,
276  const TrackerTopology* tTopo) {
277  for (std::map<int, TrackerOfflineValidationSummary::ModuleHistos>::iterator it = moduleHist.begin(),
278  itEnd = moduleHist.end();
279  it != itEnd;
280  ++it) {
281  treeMem.clear(); // make empty/default
282 
283  //variables concerning the tracker components/hierarchy levels
284  const DetId detId = it->first;
285  treeMem.moduleId = detId;
286  treeMem.subDetId = detId.subdetId();
287 
288  if (treeMem.subDetId == PixelSubdetector::PixelBarrel) {
289  unsigned int whichHalfBarrel(1), rawId(detId.rawId()); //DetId does not know about halfBarrels is PXB ...
290  if ((rawId >= 302056964 && rawId < 302059300) || (rawId >= 302123268 && rawId < 302127140) ||
291  (rawId >= 302189572 && rawId < 302194980))
292  whichHalfBarrel = 2;
293  treeMem.layer = tTopo->pxbLayer(detId);
294  treeMem.half = whichHalfBarrel;
295  treeMem.rod = tTopo->pxbLadder(
296  detId); // ... so, ladder is not per halfBarrel-Layer, but per barrel-layer! Needs complicated calculation in associateModuleHistsWithTree()
297  treeMem.module = tTopo->pxbModule(detId);
298  } else if (treeMem.subDetId == PixelSubdetector::PixelEndcap) {
299  unsigned int whichHalfCylinder(1), rawId(detId.rawId()); //DetId does not kmow about halfCylinders in PXF
300  if ((rawId >= 352394500 && rawId < 352406032) || (rawId >= 352460036 && rawId < 352471568) ||
301  (rawId >= 344005892 && rawId < 344017424) || (rawId >= 344071428 && rawId < 344082960))
302  whichHalfCylinder = 2;
303  treeMem.layer = tTopo->pxfDisk(detId);
304  treeMem.side = tTopo->pxfSide(detId);
305  treeMem.half = whichHalfCylinder;
306  treeMem.blade = tTopo->pxfBlade(detId);
307  treeMem.panel = tTopo->pxfPanel(detId);
308  treeMem.module = tTopo->pxfModule(detId);
309  } else if (treeMem.subDetId == StripSubdetector::TIB) {
310  unsigned int whichHalfShell(1), rawId(detId.rawId()); //DetId does not kmow about halfShells in TIB
311  if ((rawId >= 369120484 && rawId < 369120688) || (rawId >= 369121540 && rawId < 369121776) ||
312  (rawId >= 369136932 && rawId < 369137200) || (rawId >= 369137988 && rawId < 369138288) ||
313  (rawId >= 369153396 && rawId < 369153744) || (rawId >= 369154436 && rawId < 369154800) ||
314  (rawId >= 369169844 && rawId < 369170256) || (rawId >= 369170900 && rawId < 369171344) ||
315  (rawId >= 369124580 && rawId < 369124784) || (rawId >= 369125636 && rawId < 369125872) ||
316  (rawId >= 369141028 && rawId < 369141296) || (rawId >= 369142084 && rawId < 369142384) ||
317  (rawId >= 369157492 && rawId < 369157840) || (rawId >= 369158532 && rawId < 369158896) ||
318  (rawId >= 369173940 && rawId < 369174352) || (rawId >= 369174996 && rawId < 369175440))
319  whichHalfShell = 2;
320  treeMem.layer = tTopo->tibLayer(detId);
321  treeMem.side = tTopo->tibStringInfo(detId)[0];
322  treeMem.half = whichHalfShell;
323  treeMem.rod = tTopo->tibStringInfo(detId)[2];
324  treeMem.outerInner = tTopo->tibStringInfo(detId)[1];
325  treeMem.module = tTopo->tibModule(detId);
326  treeMem.isStereo = tTopo->tibStereo(detId);
327  treeMem.isDoubleSide = tTopo->tibIsDoubleSide(detId);
328  } else if (treeMem.subDetId == StripSubdetector::TID) {
329  treeMem.layer = tTopo->tidWheel(detId);
330  treeMem.side = tTopo->tidSide(detId);
331  treeMem.ring = tTopo->tidRing(detId);
332  treeMem.outerInner = tTopo->tidModuleInfo(detId)[0];
333  treeMem.module = tTopo->tidModuleInfo(detId)[1];
334  treeMem.isStereo = tTopo->tidStereo(detId);
335  treeMem.isDoubleSide = tTopo->tidIsDoubleSide(detId);
336  } else if (treeMem.subDetId == StripSubdetector::TOB) {
337  treeMem.layer = tTopo->tobLayer(detId);
338  treeMem.side = tTopo->tobRodInfo(detId)[0];
339  treeMem.rod = tTopo->tobRodInfo(detId)[1];
340  treeMem.module = tTopo->tobModule(detId);
341  treeMem.isStereo = tTopo->tobStereo(detId);
342  treeMem.isDoubleSide = tTopo->tobIsDoubleSide(detId);
343  } else if (treeMem.subDetId == StripSubdetector::TEC) {
344  treeMem.layer = tTopo->tecWheel(detId);
345  treeMem.side = tTopo->tecSide(detId);
346  treeMem.ring = tTopo->tecRing(detId);
347  treeMem.petal = tTopo->tecPetalInfo(detId)[1];
348  treeMem.outerInner = tTopo->tecPetalInfo(detId)[0];
349  treeMem.module = tTopo->tecModule(detId);
350  treeMem.isStereo = tTopo->tecStereo(detId);
351  treeMem.isDoubleSide = tTopo->tecIsDoubleSide(detId);
352  }
353 
354  //variables concerning the tracker geometry
355 
356  const Surface::PositionType& gPModule = tkgeom.idToDet(detId)->position();
357  treeMem.posPhi = gPModule.phi();
358  treeMem.posEta = gPModule.eta();
359  treeMem.posR = gPModule.perp();
360  treeMem.posX = gPModule.x();
361  treeMem.posY = gPModule.y();
362  treeMem.posZ = gPModule.z();
363 
364  const Surface& surface = tkgeom.idToDet(detId)->surface();
365 
366  //global Orientation of local coordinate system of dets/detUnits
367  LocalPoint lUDirection(1., 0., 0.), lVDirection(0., 1., 0.), lWDirection(0., 0., 1.);
368  GlobalPoint gUDirection = surface.toGlobal(lUDirection), gVDirection = surface.toGlobal(lVDirection),
369  gWDirection = surface.toGlobal(lWDirection);
370  double dR(999.), dPhi(999.), dZ(999.);
372  treeMem.subDetId == StripSubdetector::TOB) {
373  dR = gWDirection.perp() - gPModule.perp();
374  dPhi = deltaPhi(gUDirection.barePhi(), gPModule.barePhi());
375  dZ = gVDirection.z() - gPModule.z();
376  if (dZ >= 0.)
377  treeMem.rOrZDirection = 1;
378  else
379  treeMem.rOrZDirection = -1;
380  } else if (treeMem.subDetId == PixelSubdetector::PixelEndcap) {
381  dR = gUDirection.perp() - gPModule.perp();
382  dPhi = deltaPhi(gVDirection.barePhi(), gPModule.barePhi());
383  dZ = gWDirection.z() - gPModule.z();
384  if (dR >= 0.)
385  treeMem.rOrZDirection = 1;
386  else
387  treeMem.rOrZDirection = -1;
388  } else if (treeMem.subDetId == StripSubdetector::TID || treeMem.subDetId == StripSubdetector::TEC) {
389  dR = gVDirection.perp() - gPModule.perp();
390  dPhi = deltaPhi(gUDirection.barePhi(), gPModule.barePhi());
391  dZ = gWDirection.z() - gPModule.z();
392  if (dR >= 0.)
393  treeMem.rOrZDirection = 1;
394  else
395  treeMem.rOrZDirection = -1;
396  }
397  if (dR >= 0.)
398  treeMem.rDirection = 1;
399  else
400  treeMem.rDirection = -1;
401  if (dPhi >= 0.)
402  treeMem.phiDirection = 1;
403  else
404  treeMem.phiDirection = -1;
405  if (dZ >= 0.)
406  treeMem.zDirection = 1;
407  else
408  treeMem.zDirection = -1;
409 
410  // Assign histos from first step (TrackerOfflineValidation) to the module's entry in the TTree for retrieving mean, rms, median ...
411  const std::string histDir = associateModuleHistsWithTree(treeMem, it->second, substructureName);
412 
413  //mean and RMS values (extracted from histograms Xprime on module level)
414  treeMem.entries = static_cast<UInt_t>(it->second.ResXprimeHisto->GetEntries());
415  treeMem.meanX = it->second.ResXprimeHisto->GetMean();
416  treeMem.rmsX = it->second.ResXprimeHisto->GetRMS();
417  //treeMem.sigmaX = Fwhm(it->second.ResXprimeHisto)/2.355;
418  if (useFit_) {
419  //call fit function which returns mean and sigma from the fit
420  //for absolute residuals
421  std::pair<float, float> fitResult1 = this->fitResiduals(it->second.ResXprimeHisto);
422  treeMem.fitMeanX = fitResult1.first;
423  treeMem.fitSigmaX = fitResult1.second;
424  //for normalized residuals
425  std::pair<float, float> fitResult2 = this->fitResiduals(it->second.NormResXprimeHisto);
426  treeMem.fitMeanNormX = fitResult2.first;
427  treeMem.fitSigmaNormX = fitResult2.second;
428  }
429 
430  //get median for absolute residuals
431  treeMem.medianX = this->getMedian(it->second.ResXprimeHisto);
432 
433  int numberOfBins = it->second.ResXprimeHisto->GetNbinsX();
434  treeMem.numberOfUnderflows = it->second.ResXprimeHisto->GetBinContent(0);
435  treeMem.numberOfOverflows = it->second.ResXprimeHisto->GetBinContent(numberOfBins + 1);
436  treeMem.numberOfOutliers =
437  it->second.ResXprimeHisto->GetBinContent(0) + it->second.ResXprimeHisto->GetBinContent(numberOfBins + 1);
438  //mean and RMS values (extracted from histograms(normalized Xprime on module level)
439  treeMem.meanNormX = it->second.NormResXprimeHisto->GetMean();
440  treeMem.rmsNormX = it->second.NormResXprimeHisto->GetRMS();
441 
442  double stats[20];
443  it->second.NormResXprimeHisto->GetStats(stats);
444  // GF treeMem.chi2PerDofX = stats[3]/(stats[0]-1);
445  if (stats[0])
446  treeMem.chi2PerDofX = stats[3] / stats[0];
447 
448  //treeMem.sigmaNormX = Fwhm(it->second.NormResXprimeHisto)/2.355;
449  treeMem.histNameX = it->second.ResXprimeHisto->GetName();
450  treeMem.histNameNormX = it->second.NormResXprimeHisto->GetName();
451 
452  // fill tree variables in local coordinates if set in cfg of TrackerOfllineValidation
453  if (it->second.ResHisto && it->second.NormResHisto) { // if(lCoorHistOn_) {
454  treeMem.meanLocalX = it->second.ResHisto->GetMean();
455  treeMem.rmsLocalX = it->second.ResHisto->GetRMS();
456  treeMem.meanNormLocalX = it->second.NormResHisto->GetMean();
457  treeMem.rmsNormLocalX = it->second.NormResHisto->GetRMS();
458  treeMem.histNameLocalX = it->second.ResHisto->GetName();
459  treeMem.histNameNormLocalX = it->second.NormResHisto->GetName();
460  }
461 
462  // mean and RMS values in local y (extracted from histograms Yprime on module level)
463  // might exist in pixel only
464  if (it->second.ResYprimeHisto) { //(stripYResiduals_){
465  TH1* h = it->second.ResYprimeHisto;
466  treeMem.meanY = h->GetMean();
467  treeMem.rmsY = h->GetRMS();
468 
469  if (useFit_) { // fit function which returns mean and sigma from the fit
470  std::pair<float, float> fitMeanSigma = this->fitResiduals(h);
471  treeMem.fitMeanY = fitMeanSigma.first;
472  treeMem.fitSigmaY = fitMeanSigma.second;
473  }
474 
475  //get median for absolute residuals
476  treeMem.medianY = this->getMedian(h);
477 
478  treeMem.histNameY = h->GetName();
479  }
480 
481  if (it->second.NormResYprimeHisto) {
482  TH1* h = it->second.NormResYprimeHisto;
483  treeMem.meanNormY = h->GetMean();
484  treeMem.rmsNormY = h->GetRMS();
485  h->GetStats(stats); // stats buffer defined above
486  if (stats[0])
487  treeMem.chi2PerDofY = stats[3] / stats[0];
488 
489  if (useFit_) { // fit function which returns mean and sigma from the fit
490  std::pair<float, float> fitMeanSigma = this->fitResiduals(h);
491  treeMem.fitMeanNormY = fitMeanSigma.first;
492  treeMem.fitSigmaNormY = fitMeanSigma.second;
493  }
494  treeMem.histNameNormY = h->GetName();
495  }
496 
497  tree.Fill();
498  }
499 }
500 
502  const TkOffTreeVariables& treeMem,
504  std::map<std::string, std::string>& substructureName) {
505  std::stringstream histDir, sSubdetName;
506  std::string componentName;
507  if (moduleDirectory_.length() != 0)
508  histDir << moduleDirectory_ << "/";
509  std::string wheelOrLayer("_layer_");
510  if (treeMem.subDetId == PixelSubdetector::PixelBarrel) {
511  unsigned int half(treeMem.half), layer(treeMem.layer), ladder(0);
512  if (layer == 1) {
513  if (half == 2)
514  ladder = treeMem.rod - 5;
515  else if (treeMem.rod > 15)
516  ladder = treeMem.rod - 10;
517  else
518  ladder = treeMem.rod;
519  } else if (layer == 2) {
520  if (half == 2)
521  ladder = treeMem.rod - 8;
522  else if (treeMem.rod > 24)
523  ladder = treeMem.rod - 16;
524  else
525  ladder = treeMem.rod;
526  } else if (layer == 3) {
527  if (half == 2)
528  ladder = treeMem.rod - 11;
529  else if (treeMem.rod > 33)
530  ladder = treeMem.rod - 22;
531  else
532  ladder = treeMem.rod;
533  }
534  componentName = "Pixel";
535  sSubdetName << "TPBBarrel_1";
536  histDir << componentName << "/" << sSubdetName.str() << "/TPBHalfBarrel_" << treeMem.half << "/TPBLayer_"
537  << treeMem.layer << "/TPBLadder_" << ladder;
538  } else if (treeMem.subDetId == PixelSubdetector::PixelEndcap) {
539  unsigned int side(treeMem.side), half(treeMem.half), blade(0);
540  if (side == 1)
541  side = 3;
542  if (half == 2)
543  blade = treeMem.blade - 6;
544  else if (treeMem.blade > 18)
545  blade = treeMem.blade - 12;
546  else
547  blade = treeMem.blade;
548  componentName = "Pixel";
549  sSubdetName << "TPEEndcap_" << side;
550  histDir << componentName << "/" << sSubdetName.str() << "/TPEHalfCylinder_" << treeMem.half << "/TPEHalfDisk_"
551  << treeMem.layer << "/TPEBlade_" << blade << "/TPEPanel_" << treeMem.panel;
552  wheelOrLayer = "_wheel_";
553  } else if (treeMem.subDetId == StripSubdetector::TIB) {
554  unsigned int half(treeMem.half), layer(treeMem.layer), surface(treeMem.outerInner), string(0);
555  if (half == 2) {
556  if (layer == 1) {
557  if (surface == 1)
558  string = treeMem.rod - 13;
559  else if (surface == 2)
560  string = treeMem.rod - 15;
561  }
562  if (layer == 2) {
563  if (surface == 1)
564  string = treeMem.rod - 17;
565  else if (surface == 2)
566  string = treeMem.rod - 19;
567  }
568  if (layer == 3) {
569  if (surface == 1)
570  string = treeMem.rod - 22;
571  else if (surface == 2)
572  string = treeMem.rod - 23;
573  }
574  if (layer == 4) {
575  if (surface == 1)
576  string = treeMem.rod - 26;
577  else if (surface == 2)
578  string = treeMem.rod - 28;
579  }
580  } else
581  string = treeMem.rod;
582  std::stringstream detString;
583  if (treeMem.layer < 3 && !treeMem.isDoubleSide)
584  detString << "/Det_" << treeMem.module;
585  else
586  detString << "";
587  componentName = "Strip";
588  sSubdetName << "TIBBarrel_1";
589  histDir << componentName << "/" << sSubdetName.str() << "/TIBHalfBarrel_" << treeMem.side << "/TIBLayer_"
590  << treeMem.layer << "/TIBHalfShell_" << treeMem.half << "/TIBSurface_" << treeMem.outerInner
591  << "/TIBString_" << string << detString.str();
592  } else if (treeMem.subDetId == StripSubdetector::TID) {
593  unsigned int side(treeMem.side), outerInner(0);
594  if (side == 1)
595  side = 3;
596  if (treeMem.outerInner == 1)
597  outerInner = 2;
598  else if (treeMem.outerInner == 2)
599  outerInner = 1;
600  std::stringstream detString;
601  if (treeMem.ring < 3 && !treeMem.isDoubleSide)
602  detString << "/Det_" << treeMem.module;
603  else
604  detString << "";
605  componentName = "Strip";
606  sSubdetName << "TIDEndcap_" << side;
607  histDir << componentName << "/" << sSubdetName.str() << "/TIDDisk_" << treeMem.layer << "/TIDRing_" << treeMem.ring
608  << "/TIDSide_" << outerInner << detString.str();
609  wheelOrLayer = "_wheel_";
610  } else if (treeMem.subDetId == StripSubdetector::TOB) {
611  std::stringstream detString;
612  if (treeMem.layer < 3 && !treeMem.isDoubleSide)
613  detString << "/Det_" << treeMem.module;
614  else
615  detString << "";
616  componentName = "Strip";
617  sSubdetName << "TOBBarrel_4";
618  histDir << componentName << "/" << sSubdetName.str() << "/TOBHalfBarrel_" << treeMem.side << "/TOBLayer_"
619  << treeMem.layer << "/TOBRod_" << treeMem.rod << detString.str();
620  } else if (treeMem.subDetId == StripSubdetector::TEC) {
621  unsigned int side(0), outerInner(0), ring(0);
622  if (treeMem.side == 1)
623  side = 6;
624  else if (treeMem.side == 2)
625  side = 5;
626  if (treeMem.outerInner == 1)
627  outerInner = 2;
628  else if (treeMem.outerInner == 2)
629  outerInner = 1;
630  if (treeMem.layer > 3 && treeMem.layer < 7)
631  ring = treeMem.ring - 1;
632  else if (treeMem.layer == 7 || treeMem.layer == 8)
633  ring = treeMem.ring - 2;
634  else if (treeMem.layer == 9)
635  ring = treeMem.ring - 3;
636  else
637  ring = treeMem.ring;
638  std::stringstream detString;
639  if ((treeMem.ring < 3 || treeMem.ring == 5) && !treeMem.isDoubleSide)
640  detString << "/Det_" << treeMem.module;
641  else
642  detString << "";
643  componentName = "Strip";
644  sSubdetName << "TECEndcap_" << side;
645  histDir << componentName << "/" << sSubdetName.str() << "/TECDisk_" << treeMem.layer << "/TECSide_" << outerInner
646  << "/TECPetal_" << treeMem.petal << "/TECRing_" << ring << detString.str();
647  wheelOrLayer = "_wheel_";
648  }
649 
650  substructureName["component"] = componentName;
651  substructureName["subdet"] = sSubdetName.str();
652 
653  std::stringstream histName;
654  histName << "residuals_subdet_" << treeMem.subDetId << wheelOrLayer << treeMem.layer << "_module_"
655  << treeMem.moduleId;
656 
658  fullPath = histDir.str() + "/h_xprime_" + histName.str();
659  if (dbe_->get(fullPath))
660  moduleHists.ResXprimeHisto = dbe_->get(fullPath)->getTH1();
661  else {
662  edm::LogError("TrackerOfflineValidationSummary")
663  << "Problem with names in input file produced in TrackerOfflineValidation ...\n"
664  << "This histogram should exist in every configuration, "
665  << "but no histogram with name " << fullPath << " is found!";
666  return "";
667  }
668  fullPath = histDir.str() + "/h_normxprime" + histName.str();
669  if (dbe_->get(fullPath))
670  moduleHists.NormResXprimeHisto = dbe_->get(fullPath)->getTH1();
671  fullPath = histDir.str() + "/h_yprime_" + histName.str();
672  if (dbe_->get(fullPath))
673  moduleHists.ResYprimeHisto = dbe_->get(fullPath)->getTH1();
674  fullPath = histDir.str() + "/h_normyprime" + histName.str();
675  if (dbe_->get(fullPath))
676  moduleHists.NormResYprimeHisto = dbe_->get(fullPath)->getTH1();
677  fullPath = histDir.str() + "/h_" + histName.str();
678  if (dbe_->get(fullPath))
679  moduleHists.ResHisto = dbe_->get(fullPath)->getTH1();
680  fullPath = histDir.str() + "/h_norm" + histName.str();
681  if (dbe_->get(fullPath))
682  moduleHists.NormResHisto = dbe_->get(fullPath)->getTH1();
683 
684  return histDir.str();
685 }
686 
687 std::pair<float, float> TrackerOfflineValidationSummary::fitResiduals(TH1* hist) const {
688  std::pair<float, float> fitResult(9999., 9999.);
689  if (!hist || hist->GetEntries() < 20)
690  return fitResult;
691 
692  float mean = hist->GetMean();
693  float sigma = hist->GetRMS();
694 
695  try { // for < CMSSW_2_2_0 since ROOT warnings from fit are converted to exceptions
696  // Remove the try/catch for more recent CMSSW!
697  // first fit: two RMS around mean
698  TF1 func("tmp", "gaus", mean - 2. * sigma, mean + 2. * sigma);
699  if (0 == hist->Fit(&func, "QNR")) { // N: do not blow up file by storing fit!
700  mean = func.GetParameter(1);
701  sigma = func.GetParameter(2);
702  // second fit: three sigma of first fit around mean of first fit
703  func.SetRange(mean - 3. * sigma, mean + 3. * sigma);
704  // I: integral gives more correct results if binning is too wide
705  // L: Likelihood can treat empty bins correctly (if hist not weighted...)
706  if (0 == hist->Fit(&func, "Q0LR")) {
707  if (hist->GetFunction(func.GetName())) { // Take care that it is later on drawn:
708  hist->GetFunction(func.GetName())->ResetBit(TF1::kNotDraw);
709  }
710  fitResult.first = func.GetParameter(1);
711  fitResult.second = func.GetParameter(2);
712  }
713  }
714  } catch (cms::Exception const& e) {
715  edm::LogWarning("Alignment") << "@SUB=TrackerOfflineValidation::fitResiduals"
716  << "Caught this exception during ROOT fit: " << e.what();
717  }
718  return fitResult;
719 }
720 
722  float median = 999;
723  const int nbins = histo->GetNbinsX();
724 
725  //extract median from histogram
726  double* x = new double[nbins];
727  double* y = new double[nbins];
728  for (int j = 0; j < nbins; j++) {
729  x[j] = histo->GetBinCenter(j + 1);
730  y[j] = histo->GetBinContent(j + 1);
731  }
732  median = TMath::Median(nbins, x, y);
733 
734  delete[] x;
735  x = nullptr;
736  delete[] y;
737  y = nullptr;
738 
739  return median;
740 }
741 
743  this->applyHarvestingHierarchy(tree);
744  this->fillHarvestingHists(tree);
745 }
746 
748  TkOffTreeVariables* treeMemPtr = nullptr;
749  std::map<std::string, std::string>* substructureName = nullptr;
750  tree.SetBranchAddress("TkOffTreeVariables", &treeMemPtr);
751  tree.SetBranchAddress("SubstructureName", &substructureName);
752 
753  // Loop over modules to select accumulation criteria for harvesting plots
754  for (unsigned int iSubdet = 1; iSubdet < 7; ++iSubdet) {
755  std::string hierarchyName("");
756  std::string componentName("");
757  std::vector<unsigned int> treeEntries;
758  for (unsigned int iSide = 1; iSide < 3; ++iSide) {
759  // Set up only one collection for Barrels, not separated for side
760  if (iSide == 1 && (iSubdet == PixelSubdetector::PixelBarrel || iSubdet == StripSubdetector::TIB ||
761  iSubdet == StripSubdetector::TOB))
762  continue;
763  for (int iTree = 0; iTree < tree.GetEntries(); ++iTree) {
764  tree.GetEntry(iTree);
765  // Do not use glued Dets
766  if (treeMemPtr->isDoubleSide)
767  continue;
768  if (treeMemPtr->subDetId == iSubdet) {
769  if (iSide != treeMemPtr->side && (iSubdet == PixelSubdetector::PixelEndcap ||
770  iSubdet == StripSubdetector::TID || iSubdet == StripSubdetector::TEC))
771  continue;
772  treeEntries.push_back(iTree);
773  if (hierarchyName.length() == 0) {
774  hierarchyName = (*substructureName)["subdet"];
775  componentName = (*substructureName)["component"];
776  }
777  }
778  }
779  HarvestingHierarchy harvestingHierarchy(hierarchyName, componentName, treeEntries);
780  vHarvestingHierarchy_.push_back(harvestingHierarchy);
781  hierarchyName = "";
782  componentName = "";
783  treeEntries.clear();
784  }
785  }
786  // Here could be a further separation of the HarvestingHierarchy.
787  // E.g. separate the existing ones by layer and add them to the vector without deleting any element from the vector.
788  // The existing hists will stay and the new ones are added
789 
790  // Now, book the corresponding histos
791  this->bookHarvestingHists();
792 }
793 
795  edm::LogInfo("TrackerOfflineValidationSummary") << "Harvesting histograms will be booked for "
796  << vHarvestingHierarchy_.size() << " different hierarchy selections";
797  for (std::vector<HarvestingHierarchy>::iterator iHier = vHarvestingHierarchy_.begin();
798  iHier != vHarvestingHierarchy_.end();
799  ++iHier) {
800  std::stringstream dmrXprimeHistoName, dmrYprimeHistoName, dmrXprimeHistoTitle, dmrYprimeHistoTitle;
801  dmrXprimeHistoName << "h_DmrXprime_" << iHier->hierarchyName;
802  dmrYprimeHistoName << "h_DmrYprime_" << iHier->hierarchyName;
803  dmrXprimeHistoTitle << "DMR for " << iHier->hierarchyName << ";<#DeltaX> [cm];# modules";
804  dmrYprimeHistoTitle << "DMR for " << iHier->hierarchyName << ";<#DeltaY> [cm];# modules";
805 
806  std::string directoryString(moduleDirectory_);
807  if (directoryString.length() != 0)
808  directoryString += "/";
809  directoryString += iHier->componentName;
810  dbe_->setCurrentFolder(directoryString);
811 
812  int nBinsX(0);
813  double xMin(0.), xMax(0.);
814  if (iHier->componentName == "Pixel") {
815  this->getBinning("TH1DmrXprimePixelModules", nBinsX, xMin, xMax);
816  iHier->harvestingHistos.DmrXprime =
817  dbe_->book1D(dmrXprimeHistoName.str(), dmrXprimeHistoTitle.str(), nBinsX, xMin, xMax)->getTH1();
818  this->getBinning("TH1DmrYprimePixelModules", nBinsX, xMin, xMax);
819  iHier->harvestingHistos.DmrYprime =
820  dbe_->book1D(dmrYprimeHistoName.str(), dmrYprimeHistoTitle.str(), nBinsX, xMin, xMax)->getTH1();
821  } else if (iHier->componentName == "Strip") {
822  this->getBinning("TH1DmrXprimeStripModules", nBinsX, xMin, xMax);
823  iHier->harvestingHistos.DmrXprime =
824  dbe_->book1D(dmrXprimeHistoName.str(), dmrXprimeHistoTitle.str(), nBinsX, xMin, xMax)->getTH1();
825  if (!parSet_.getParameter<bool>("stripYDmrs"))
826  continue;
827  this->getBinning("TH1DmrYprimeStripModules", nBinsX, xMin, xMax);
828  iHier->harvestingHistos.DmrYprime =
829  dbe_->book1D(dmrYprimeHistoName.str(), dmrYprimeHistoTitle.str(), nBinsX, xMin, xMax)->getTH1();
830  }
831  }
832 }
833 
835  int& nBinsX,
836  double& lowerBoundX,
837  double& upperBoundX) const {
838  const edm::ParameterSet& binningPSet = parSet_.getParameter<edm::ParameterSet>(binningPSetName);
839  nBinsX = binningPSet.getParameter<int>("Nbinx");
840  lowerBoundX = binningPSet.getParameter<double>("xmin");
841  upperBoundX = binningPSet.getParameter<double>("xmax");
842 }
843 
845  TkOffTreeVariables* treeMemPtr = nullptr;
846  std::map<std::string, std::string>* substructureName = nullptr;
847  tree.SetBranchAddress("TkOffTreeVariables", &treeMemPtr);
848  tree.SetBranchAddress("SubstructureName", &substructureName);
849 
850  const unsigned int minEntriesPerModule(parSet_.getParameter<unsigned int>("minEntriesPerModuleForDmr"));
851  edm::LogInfo("TrackerOfflineValidationSummary")
852  << "Median of a module is added to DMR plots if it contains at least " << minEntriesPerModule << " hits";
853 
854  for (std::vector<HarvestingHierarchy>::iterator iHier = vHarvestingHierarchy_.begin();
855  iHier != vHarvestingHierarchy_.end();
856  ++iHier) {
857  for (std::vector<unsigned int>::const_iterator iTreeEntries = iHier->treeEntries.begin();
858  iTreeEntries != iHier->treeEntries.end();
859  ++iTreeEntries) {
860  tree.GetEntry(*iTreeEntries);
861  if (treeMemPtr->entries < minEntriesPerModule)
862  continue;
863  iHier->harvestingHistos.DmrXprime->Fill(treeMemPtr->medianX);
864  if (iHier->harvestingHistos.DmrYprime)
865  iHier->harvestingHistos.DmrYprime->Fill(treeMemPtr->medianY);
866  }
867  }
868 }
869 
870 //define this as a plug-in
std::map< int, TrackerOfflineValidationSummary::ModuleHistos > mTobResiduals_
static constexpr auto TEC
bool tibIsDoubleSide(const DetId &id) const
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
void addWithDefaultLabel(ParameterSetDescription const &psetDescription)
bool tecIsDoubleSide(const DetId &id) const
bool tidIsDoubleSide(const DetId &id) const
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
unsigned int tobLayer(const DetId &id) const
unsigned int pxbLayer(const DetId &id) const
T perp() const
Definition: PV3DBase.h:69
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
Definition: EventSetup.h:119
std::vector< unsigned int > tidModuleInfo(const DetId &id) const
unsigned int pxfBlade(const DetId &id) const
T z() const
Definition: PV3DBase.h:61
Put here the histograms created during harvesting.
const DetIdContainer & detIds() const override
Returm a vector of all GeomDet DetIds (including those of GeomDetUnits)
uint32_t tidStereo(const DetId &id) const
unsigned int tibModule(const DetId &id) const
unsigned int tidSide(const DetId &id) const
container to hold data to be written into TTree
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
unsigned int pxfModule(const DetId &id) const
std::vector< unsigned int > tecPetalInfo(const DetId &id) const
unsigned int tidWheel(const DetId &id) const
unsigned int tecWheel(const DetId &id) const
T eta() const
Definition: PV3DBase.h:73
void setCurrentFolder(std::string const &fullpath) override
Definition: DQMStore.h:646
std::vector< HarvestingHierarchy > vHarvestingHierarchy_
void getBinning(const std::string &binningPSetName, int &nBinsX, double &lowerBoundX, double &upperBoundX) const
edm::ESHandle< TrackerGeometry > tkGeom_
unsigned int pxbLadder(const DetId &id) const
Log< level::Error, false > LogError
constexpr Detector det() const
get the detector field from this detid
Definition: DetId.h:46
unsigned int tecRing(const DetId &id) const
ring id
std::map< int, TrackerOfflineValidationSummary::ModuleHistos > mTibResiduals_
bool tobIsDoubleSide(const DetId &id) const
T barePhi() const
Definition: PV3DBase.h:65
numberOfBins
Definition: PV_cfg.py:227
constexpr std::array< uint8_t, layerIndexSize< TrackerTraits > > layer
T x() const
Definition: PV3DBase.h:59
T y() const
Definition: PV3DBase.h:60
unsigned int tecModule(const DetId &id) const
int iEvent
Definition: GenABIO.cc:224
uint32_t tobStereo(const DetId &id) const
unsigned int tecSide(const DetId &id) const
const std::string associateModuleHistsWithTree(const TkOffTreeVariables &treeMem, TrackerOfflineValidationSummary::ModuleHistos &moduleHists, std::map< std::string, std::string > &substructureName)
unsigned int pxfDisk(const DetId &id) const
Transition
Definition: Transition.h:12
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
static constexpr auto TOB
ParameterDescriptionBase * add(U const &iLabel, T const &value)
histDir
Location of output EventInfo/reportSummaryContents.
TrackerOfflineValidationSummary(const edm::ParameterSet &)
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:130
void analyze(const edm::Event &evt, const edm::EventSetup &) override
std::pair< float, float > fitResiduals(TH1 *hist) const
const TrackerGeomDet * idToDet(DetId) const override
std::string histNameNormLocalX
std::map< int, TrackerOfflineValidationSummary::ModuleHistos > mTecResiduals_
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:79
unsigned int pxfPanel(const DetId &id) const
std::map< int, TrackerOfflineValidationSummary::ModuleHistos > mPxbResiduals_
Log< level::Info, false > LogInfo
Definition: DetId.h:17
HarvestingHierarchy(const std::string &name, const std::string &component, const std::vector< unsigned int > &entries)
unsigned int pxfSide(const DetId &id) const
std::vector< unsigned int > tibStringInfo(const DetId &id) const
static constexpr auto TIB
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:37
void beginRun(const edm::Run &, const edm::EventSetup &iSetup) override
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
std::string histNameLocalX
virtual MonitorElement * get(std::string const &fullpath) const
Definition: DQMStore.cc:697
void endRun(const edm::Run &, const edm::EventSetup &iSetup) override
const Surface::PositionType & position() const
The position (origin of the R.F.)
Definition: GeomDet.h:43
std::vector< DetId > DetIdContainer
virtual TH1 * getTH1() const
T median(std::vector< T > values)
Definition: median.h:16
std::vector< unsigned int > tobRodInfo(const DetId &id) const
void clear()
set to empty values
HLT enums.
const edm::ESGetToken< TrackerTopology, TrackerTopologyRcd > tTopoToken_
unsigned int tidRing(const DetId &id) const
std::unique_ptr< TrackerTopology > tTopo_
std::map< int, TrackerOfflineValidationSummary::ModuleHistos > mPxeResiduals_
uint32_t tecStereo(const DetId &id) const
~TrackerOfflineValidationSummary() override=default
const edm::ESGetToken< TrackerGeometry, TrackerDigiGeometryRecord > geomToken_
unsigned int tibLayer(const DetId &id) const
std::map< int, TrackerOfflineValidationSummary::ModuleHistos > mTidResiduals_
unsigned int tobModule(const DetId &id) const
Definition: tree.py:1
unsigned int pxbModule(const DetId &id) const
Log< level::Warning, false > LogWarning
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
void fillTree(TTree &tree, std::map< int, TrackerOfflineValidationSummary::ModuleHistos > &moduleHist, TkOffTreeVariables &treeMem, const TrackerGeometry &tkgeom, std::map< std::string, std::string > &substructureName, const TrackerTopology *tTopo)
static constexpr auto TID
uint32_t tibStereo(const DetId &id) const
Definition: Run.h:45