CMS 3D CMS Logo

TrackerOfflineValidation.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: TrackerOfflineValidation
4 // Class: TrackerOfflineValidation
5 //
13 //
14 // Original Author: Erik Butz
15 // Created: Tue Dec 11 14:03:05 CET 2007
16 // $Id: TrackerOfflineValidation.cc,v 1.55 2012/09/28 20:48:06 flucke Exp $
17 //
18 //
19 
20 // system include files
21 #include <memory>
22 #include <map>
23 #include <sstream>
24 #include <cmath>
25 #include <tuple>
26 #include <utility>
27 #include <vector>
28 #include <iostream>
29 
30 // ROOT includes
31 #include "TH1.h"
32 #include "TH2.h"
33 #include "TProfile.h"
34 #include "TFile.h"
35 #include "TTree.h"
36 #include "TF1.h"
37 #include "TMath.h"
38 
39 // user include files
48 
54 
57 
62 
64 
70 
73 
74 //
75 // class declaration
76 //
77 class TrackerOfflineValidation : public edm::one::EDAnalyzer<edm::one::SharedResources> {
78 public:
81  ~TrackerOfflineValidation() override;
82  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
83 
87  YResidual, /*NormYResidual, */
94  };
95 
96 private:
99 
100  struct ModuleHistos {
102  : ResHisto(),
103  NormResHisto(),
104  ResYHisto(), /*NormResYHisto(),*/
105  ResXprimeHisto(),
107  ResYprimeHisto(),
109  ResXvsXProfile(),
110  ResXvsYProfile(),
111  ResYvsXProfile(),
112  ResYvsYProfile(),
113  LocalX(),
114  LocalY() {}
115  TH1* ResHisto;
117  TH1* ResYHisto;
118  /* TH1* NormResYHisto; */
123 
124  TProfile* ResXvsXProfile;
125  TProfile* ResXvsYProfile;
126  TProfile* ResYvsXProfile;
127  TProfile* ResYvsYProfile;
128 
129  TH1* LocalX;
130  TH1* LocalY;
131 
132  unsigned int EntriesInt;
133  };
134 
135  // container struct to organize collection of histograms during endJob
138  : sumXResiduals_(),
142  sumYResiduals_(),
149  sumResYvsYProfile_() {}
150 
159 
164  };
165 
168  const std::string& newDir,
169  const std::string& basedir,
170  bool useDqmMode)
171  : tfd(nullptr), dqmMode(useDqmMode), theDbe(nullptr) {
172  if (newDir.length() != 0) {
173  if (upDir.directoryString.length() != 0)
174  directoryString = upDir.directoryString + "/" + newDir;
175  else
176  directoryString = newDir;
177  } else
179 
180  if (!dqmMode) {
181  if (newDir.length() == 0)
182  tfd.reset(&(*upDir.tfd));
183  else
184  tfd = std::make_unique<TFileDirectory>(upDir.tfd->mkdir(newDir));
185  } else {
187  }
188  }
189 
190  DirectoryWrapper(const std::string& newDir, const std::string& basedir, bool useDqmMode)
191  : tfd(nullptr), dqmMode(useDqmMode), theDbe(nullptr) {
192  if (!dqmMode) {
194  if (newDir.length() == 0) {
195  tfd = std::make_unique<TFileDirectory>(fs->tFileDirectory());
196  } else {
197  tfd = std::make_unique<TFileDirectory>(fs->mkdir(newDir));
198  directoryString = newDir;
199  }
200  } else {
201  if (newDir.length() != 0) {
202  if (basedir.length() != 0)
203  directoryString = basedir + "/" + newDir;
204  else
205  directoryString = newDir;
206  } else
209  }
210  }
211  // Generalization of Histogram Booking; allows switch between TFileService and DQMStore
212  template <typename T>
213  TH1* make(const char* name, const char* title, int nBinX, double minBinX, double maxBinX);
214  template <typename T>
215  TH1* make(const char* name, const char* title, int nBinX, double* xBins); //variable bin size in x for profile histo
216  template <typename T>
217  TH1* make(const char* name,
218  const char* title,
219  int nBinX,
220  double minBinX,
221  double maxBinX,
222  int nBinY,
223  double minBinY,
224  double maxBinY);
225  template <typename T>
226  TH1* make(const char* name,
227  const char* title,
228  int nBinX,
229  double minBinX,
230  double maxBinX,
231  double minBinY,
232  double maxBinY); // at present not used
233 
234  std::unique_ptr<TFileDirectory> tfd;
236  const bool dqmMode;
238  };
239 
240  //
241  // ------------- private member function -------------
242  //
243  void analyze(const edm::Event&, const edm::EventSetup&) override;
244  void endJob() override;
245 
246  virtual void checkBookHists(const edm::EventSetup& setup);
247 
249  void bookDirHists(DirectoryWrapper& tfd, const Alignable& ali, const TrackerTopology* tTopo);
250  void bookHists(
251  DirectoryWrapper& tfd, const Alignable& ali, const TrackerTopology* tTopo, align::StructureType type, int i);
252 
254  const Alignable& ali,
255  std::vector<TrackerOfflineValidation::SummaryContainer>& vLevelProfiles);
256  void collateSummaryHists();
257 
258  void setUpTreeMembers(const std::map<int, TrackerOfflineValidation::ModuleHistos>& moduleHist_,
259  const TrackerGeometry& tkgeom,
260  const TrackerTopology* tTopo);
261  void fillTree(TTree& tree,
262  TkOffTreeVariables& treeMem,
263  const std::map<int, TrackerOfflineValidation::ModuleHistos>& moduleHist_);
264 
266  const Alignable& ali,
268  int i);
269 
271 
272  bool isBarrel(uint32_t subDetId);
273  bool isEndCap(uint32_t subDetId);
274  bool isPixel(uint32_t subDetId);
276 
277  TH1* bookTH1F(bool isTransient,
278  DirectoryWrapper& tfd,
279  const char* histName,
280  const char* histTitle,
281  int nBinsX,
282  double lowX,
283  double highX);
284 
285  TProfile* bookTProfile(bool isTransient,
286  DirectoryWrapper& tfd,
287  const char* histName,
288  const char* histTitle,
289  int nBinsX,
290  double lowX,
291  double highX);
292 
293  TProfile* bookTProfile(bool isTransient,
294  DirectoryWrapper& tfd,
295  const char* histName,
296  const char* histTitle,
297  int nBinsX,
298  double lowX,
299  double highX,
300  double lowY,
301  double highY);
302 
303  void getBinning(uint32_t subDetId,
305  int& nBinsX,
306  double& lowerBoundX,
307  double& upperBoundX);
308 
309  void summarizeBinInContainer(int bin, SummaryContainer& targetContainer, SummaryContainer& sourceContainer);
310 
311  void summarizeBinInContainer(int bin,
312  uint32_t subDetId,
313  SummaryContainer& targetContainer,
314  ModuleHistos& sourceContainer);
315 
316  void setSummaryBin(int bin, TH1* targetHist, TH1* sourceHist);
317 
318  float Fwhm(const TH1* hist) const;
319  std::pair<float, float> fitResiduals(TH1* hist) const; //, float meantmp, float rmstmp);
320  float getMedian(const TH1* hist) const;
321 
322  // From MillePedeAlignmentMonitor: Get Index for Arbitary vector<class> by name
323  template <class OBJECT_TYPE>
324  int GetIndex(const std::vector<OBJECT_TYPE*>& vec, const TString& name);
325 
326  // ---------- member data ---------------------------
327 
330  const TrackerGeometry* bareTkGeomPtr_; // ugly hack to book hists only once, but check
331  std::unique_ptr<AlignableTracker> alignableTracker_;
332 
333  // parameters from cfg to steer
335  const bool lCoorHistOn_;
338  const bool stripYResiduals_;
339  const bool useFwhm_;
340  const bool useFit_;
341  const bool useOverflowForRMS_;
342  const bool dqmMode_;
344 
345  const int chargeCut_;
346 
347  // a vector to keep track which pointers should be deleted at the very end
348  std::vector<TH1*> vDeleteObjects_;
349 
350  std::vector<TH1*> vTrackHistos_;
351  std::vector<TH1*> vTrackProfiles_;
352  std::vector<TH1*> vTrack2DHistos_;
353 
354  std::map<int, TrackerOfflineValidation::ModuleHistos> mPxbResiduals_;
355  std::map<int, TrackerOfflineValidation::ModuleHistos> mPxeResiduals_;
356  std::map<int, TrackerOfflineValidation::ModuleHistos> mTibResiduals_;
357  std::map<int, TrackerOfflineValidation::ModuleHistos> mTidResiduals_;
358  std::map<int, TrackerOfflineValidation::ModuleHistos> mTobResiduals_;
359  std::map<int, TrackerOfflineValidation::ModuleHistos> mTecResiduals_;
360 
361  std::map<int, TkOffTreeVariables> mTreeMembers_;
362 
363  //There are two types of summary histograms. The first contains, for each component, a bin per subcomponent.
364  //These are set up through summarizeBinInContainer(). The second type is just the sum of the lower level histograms.
365  //Prepare the filling of both types at the beginning, when the tracker topology is available through the eventsetup.
366  //They are not actually filled until the end, but at that time eventsetup is no longer accessible.
367 
368  //To fill the summary hists, store the arguments of setSummaryBin()
369  //At the end, call the function
370  std::vector<std::tuple<int, TH1*, TH1*> > summaryBins_;
371  //To fill the sum hists, just store pairs of TH1*. At the end, first->Add(second).
372  std::vector<std::pair<TH1*, TH1*> > sumHistStructure_;
373  //sum hists are fit at the end using fitResiduals()
374  std::vector<TH1*> toFit_;
375 
376  unsigned long long nTracks_;
377  const unsigned long long maxTracks_;
378  const unsigned int maxEntriesPerModuleForDmr_;
380 };
381 
382 //
383 // constants, enums and typedefs
384 //
385 
386 //
387 // static data member definitions
388 //
389 
390 template <class OBJECT_TYPE>
391 int TrackerOfflineValidation::GetIndex(const std::vector<OBJECT_TYPE*>& vec, const TString& name) {
392  int result = 0;
393  for (typename std::vector<OBJECT_TYPE*>::const_iterator iter = vec.begin(), iterEnd = vec.end(); iter != iterEnd;
394  ++iter, ++result) {
395  if (*iter && (*iter)->GetName() == name)
396  return result;
397  }
398  edm::LogError("Alignment") << "@SUB=TrackerOfflineValidation::GetIndex"
399  << " could not find " << name;
400  return -1;
401 }
402 
403 template <>
404 TH1* TrackerOfflineValidation::DirectoryWrapper::make<TH1F>(
405  const char* name, const char* title, int nBinX, double minBinX, double maxBinX) {
406  if (dqmMode) {
407  theDbe->setCurrentFolder(directoryString);
408  return theDbe->book1D(name, title, nBinX, minBinX, maxBinX)->getTH1();
409  } else {
410  return tfd->make<TH1F>(name, title, nBinX, minBinX, maxBinX);
411  }
412 }
413 
414 template <>
415 TH1* TrackerOfflineValidation::DirectoryWrapper::make<TProfile>(const char* name,
416  const char* title,
417  int nBinX,
418  double* xBins) {
419  if (dqmMode) {
420  theDbe->setCurrentFolder(directoryString);
421  //DQM profile requires y-bins for construction... using TProfile creator by hand...
422  TProfile* tmpProfile = new TProfile(name, title, nBinX, xBins);
423  tmpProfile->SetDirectory(nullptr);
424  return theDbe->bookProfile(name, tmpProfile)->getTH1();
425  } else {
426  return tfd->make<TProfile>(name, title, nBinX, xBins);
427  }
428 }
429 
430 template <>
431 TH1* TrackerOfflineValidation::DirectoryWrapper::make<TProfile>(
432  const char* name, const char* title, int nBinX, double minBinX, double maxBinX) {
433  if (dqmMode) {
434  theDbe->setCurrentFolder(directoryString);
435  //DQM profile requires y-bins for construction... using TProfile creator by hand...
436  TProfile* tmpProfile = new TProfile(name, title, nBinX, minBinX, maxBinX);
437  tmpProfile->SetDirectory(nullptr);
438  return theDbe->bookProfile(name, tmpProfile)->getTH1();
439  } else {
440  return tfd->make<TProfile>(name, title, nBinX, minBinX, maxBinX);
441  }
442 }
443 
444 template <>
445 TH1* TrackerOfflineValidation::DirectoryWrapper::make<TProfile>(
446  const char* name, const char* title, int nbinX, double minX, double maxX, double minY, double maxY) {
447  if (dqmMode) {
448  theDbe->setCurrentFolder(directoryString);
449  int dummy(0); // DQMProfile wants Y channels... does not use them!
450  return (theDbe->bookProfile(name, title, nbinX, minX, maxX, dummy, minY, maxY)->getTH1());
451  } else {
452  return tfd->make<TProfile>(name, title, nbinX, minX, maxX, minY, maxY);
453  }
454 }
455 
456 template <>
457 TH1* TrackerOfflineValidation::DirectoryWrapper::make<TH2F>(const char* name,
458  const char* title,
459  int nBinX,
460  double minBinX,
461  double maxBinX,
462  int nBinY,
463  double minBinY,
464  double maxBinY) {
465  if (dqmMode) {
466  theDbe->setCurrentFolder(directoryString);
467  return theDbe->book2D(name, title, nBinX, minBinX, maxBinX, nBinY, minBinY, maxBinY)->getTH1();
468  } else {
469  return tfd->make<TH2F>(name, title, nBinX, minBinX, maxBinX, nBinY, minBinY, maxBinY);
470  }
471 }
472 
473 //
474 // constructors and destructor
475 //
477  : geomToken_(esConsumes()),
478  topoToken_(esConsumes()),
479  parSet_(iConfig),
480  bareTkGeomPtr_(nullptr),
481  compressionSettings_(parSet_.getUntrackedParameter<int>("compressionSettings", -1)),
482  lCoorHistOn_(parSet_.getParameter<bool>("localCoorHistosOn")),
483  moduleLevelHistsTransient_(parSet_.getParameter<bool>("moduleLevelHistsTransient")),
484  moduleLevelProfiles_(parSet_.getParameter<bool>("moduleLevelProfiles")),
485  stripYResiduals_(parSet_.getParameter<bool>("stripYResiduals")),
486  useFwhm_(parSet_.getParameter<bool>("useFwhm")),
487  useFit_(parSet_.getParameter<bool>("useFit")),
488  useOverflowForRMS_(parSet_.getParameter<bool>("useOverflowForRMS")),
489  dqmMode_(parSet_.getParameter<bool>("useInDqmMode")),
490  moduleDirectory_(parSet_.getParameter<std::string>("moduleDirectoryInOutput")),
491  chargeCut_(parSet_.getParameter<int>("chargeCut")),
492  nTracks_(0),
493  maxTracks_(parSet_.getParameter<unsigned long long>("maxTracks")),
494  maxEntriesPerModuleForDmr_(parSet_.getParameter<unsigned int>("maxEntriesPerModuleForDmr")),
495  avalidator_(iConfig, consumesCollector()) {
496  usesResource(TFileService::kSharedResource);
497 }
498 
501  desc.setComment("Validates alignment payloads by evaluating unbiased track-hit resisuals");
503  desc.addUntracked<int>("compressionSettings", -1);
504  desc.add<bool>("localCoorHistosOn", false);
505  desc.add<bool>("moduleLevelHistsTransient", false);
506  desc.add<bool>("moduleLevelProfiles", false);
507  desc.add<bool>("stripYResiduals", false);
508  desc.add<bool>("useFwhm", true);
509  desc.add<bool>("useFit", false);
510  desc.add<bool>("useOverflowForRMS", false);
511  desc.add<bool>("useInDqmMode", false);
512  desc.add<std::string>("moduleDirectoryInOutput", {});
513  desc.add<int>("chargeCut", 0);
514  desc.add<unsigned long long>("maxTracks", 0);
515  desc.add<unsigned int>("maxEntriesPerModuleForDmr", 0);
516 
517  // fill in the residuals details
518  std::vector<std::string> listOfResidualsPSets = {"TH1XResPixelModules",
519  "TH1XResStripModules",
520  "TH1NormXResPixelModules",
521  "TH1NormXResStripModules",
522  "TH1XprimeResPixelModules",
523  "TH1XprimeResStripModules",
524  "TH1NormXprimeResPixelModules",
525  "TH1NormXprimeResStripModules",
526  "TH1YResPixelModules",
527  "TH1YResStripModules",
528  "TH1NormYResPixelModules",
529  "TH1NormYResStripModules",
530  "TProfileXResPixelModules",
531  "TProfileXResStripModules",
532  "TProfileYResPixelModules",
533  "TProfileYResStripModules"};
534 
535  for (const auto& myPSetName : listOfResidualsPSets) {
537  myPSet.add<int>("Nbinx", 100);
538  myPSet.add<double>("xmin", -5.f);
539  myPSet.add<double>("xmax", 5.f);
540  desc.add<edm::ParameterSetDescription>(myPSetName, myPSet);
541  }
542  descriptions.addWithDefaultLabel(desc);
543 }
544 
546  // do anything here that needs to be done at desctruction time
547  // (e.g. close files, deallocate resources etc.)
548  for (std::vector<TH1*>::const_iterator it = vDeleteObjects_.begin(), itEnd = vDeleteObjects_.end(); it != itEnd; ++it)
549  delete *it;
550 }
551 
552 //
553 // member functions
554 //
555 
556 // ------------ method called once each job just before starting event loop ------------
559  const TrackerGeometry* newBareTkGeomPtr = &(*tkGeom_);
560 
561  if (newBareTkGeomPtr == bareTkGeomPtr_)
562  return; // already booked hists, nothing changed
563 
564  if (!bareTkGeomPtr_) { // pointer not yet set: called the first time => book hists
565 
566  //Retrieve tracker topology from geometry
567  const TrackerTopology* const tTopo = &es.getData(topoToken_);
568 
569  // construct alignable tracker to get access to alignable hierarchy
570  alignableTracker_ = std::make_unique<AlignableTracker>(&(*tkGeom_), tTopo);
571 
572  edm::LogInfo("TrackerOfflineValidation")
573  << "There are " << newBareTkGeomPtr->detIds().size() << " dets in the Geometry record.\n"
574  << "Out of these " << newBareTkGeomPtr->detUnitIds().size() << " are detUnits";
575 
576  // Book histograms for global track quantities
577  std::string globDir("GlobalTrackVariables");
578  DirectoryWrapper trackglobal(globDir, moduleDirectory_, dqmMode_);
579  this->bookGlobalHists(trackglobal);
580 
581  // recursively book histograms on lowest level
583  this->bookDirHists(tfdw, *alignableTracker_, tTopo);
584  // and prepare the higher level histograms
585  std::vector<TrackerOfflineValidation::SummaryContainer> vTrackerprofiles;
586  this->prepareSummaryHists(tfdw, *alignableTracker_, vTrackerprofiles);
587 
588  setUpTreeMembers(mPxbResiduals_, *newBareTkGeomPtr, tTopo);
589  setUpTreeMembers(mPxeResiduals_, *newBareTkGeomPtr, tTopo);
590  setUpTreeMembers(mTibResiduals_, *newBareTkGeomPtr, tTopo);
591  setUpTreeMembers(mTidResiduals_, *newBareTkGeomPtr, tTopo);
592  setUpTreeMembers(mTobResiduals_, *newBareTkGeomPtr, tTopo);
593  setUpTreeMembers(mTecResiduals_, *newBareTkGeomPtr, tTopo);
594  } else { // histograms booked, but changed TrackerGeometry?
595  edm::LogWarning("GeometryChange") << "@SUB=checkBookHists"
596  << "TrackerGeometry changed, but will not re-book hists!";
597  }
598  bareTkGeomPtr_ = newBareTkGeomPtr;
599 }
600 
602  vTrackHistos_.push_back(tfd.make<TH1F>("h_tracketa", "Track #eta;#eta_{Track};Number of Tracks", 90, -3., 3.));
603  vTrackHistos_.push_back(tfd.make<TH1F>("h_trackphi", "Track #phi;#phi_{Track};Number of Tracks", 90, -3.15, 3.15));
604  vTrackHistos_.push_back(tfd.make<TH1F>(
605  "h_trackNumberOfValidHits", "Track # of valid hits;# of valid hits _{Track};Number of Tracks", 40, 0., 40.));
606  vTrackHistos_.push_back(tfd.make<TH1F>(
607  "h_trackNumberOfLostHits", "Track # of lost hits;# of lost hits _{Track};Number of Tracks", 10, 0., 10.));
608  vTrackHistos_.push_back(
609  tfd.make<TH1F>("h_curvature", "Curvature #kappa;#kappa_{Track};Number of Tracks", 100, -.05, .05));
610  vTrackHistos_.push_back(tfd.make<TH1F>(
611  "h_curvature_pos", "Curvature |#kappa| Positive Tracks;|#kappa_{pos Track}|;Number of Tracks", 100, .0, .05));
612  vTrackHistos_.push_back(tfd.make<TH1F>(
613  "h_curvature_neg", "Curvature |#kappa| Negative Tracks;|#kappa_{neg Track}|;Number of Tracks", 100, .0, .05));
614  vTrackHistos_.push_back(
615  tfd.make<TH1F>("h_diff_curvature",
616  "Curvature |#kappa| Tracks Difference;|#kappa_{Track}|;# Pos Tracks - # Neg Tracks",
617  100,
618  .0,
619  .05));
620  vTrackHistos_.push_back(tfd.make<TH1F>("h_chi2", "#chi^{2};#chi^{2}_{Track};Number of Tracks", 500, -0.01, 500.));
621  vTrackHistos_.push_back(
622  tfd.make<TH1F>("h_chi2Prob", "#chi^{2} probability;#chi^{2}prob_{Track};Number of Tracks", 100, 0.0, 1.));
623  vTrackHistos_.push_back(
624  tfd.make<TH1F>("h_normchi2", "#chi^{2}/ndof;#chi^{2}/ndof;Number of Tracks", 100, -0.01, 10.));
625  vTrackHistos_.push_back(tfd.make<TH1F>("h_pt", "p_{T}^{track};p_{T}^{track} [GeV];Number of Tracks", 250, 0., 250));
626  vTrackHistos_.push_back(tfd.make<TH1F>(
627  "h_ptResolution", "#delta_{p_{T}}/p_{T}^{track};#delta_{p_{T}}/p_{T}^{track};Number of Tracks", 100, 0., 0.5));
628  vTrackProfiles_.push_back(tfd.make<TProfile>(
629  "p_d0_vs_phi", "Transverse Impact Parameter vs. #phi;#phi_{Track};#LT d_{0} #GT [cm]", 100, -3.15, 3.15));
630  vTrackProfiles_.push_back(tfd.make<TProfile>(
631  "p_dz_vs_phi", "Longitudinal Impact Parameter vs. #phi;#phi_{Track};#LT d_{z} #GT [cm]", 100, -3.15, 3.15));
632  vTrackProfiles_.push_back(tfd.make<TProfile>(
633  "p_d0_vs_eta", "Transverse Impact Parameter vs. #eta;#eta_{Track};#LT d_{0} #GT [cm]", 100, -3.15, 3.15));
634  vTrackProfiles_.push_back(tfd.make<TProfile>(
635  "p_dz_vs_eta", "Longitudinal Impact Parameter vs. #eta;#eta_{Track};#LT d_{z} #GT [cm]", 100, -3.15, 3.15));
636  vTrackProfiles_.push_back(
637  tfd.make<TProfile>("p_chi2_vs_phi", "#chi^{2} vs. #phi;#phi_{Track};#LT #chi^{2} #GT", 100, -3.15, 3.15));
638  vTrackProfiles_.push_back(tfd.make<TProfile>(
639  "p_chi2Prob_vs_phi", "#chi^{2} probablility vs. #phi;#phi_{Track};#LT #chi^{2} probability#GT", 100, -3.15, 3.15));
640  vTrackProfiles_.push_back(tfd.make<TProfile>(
641  "p_chi2Prob_vs_d0", "#chi^{2} probablility vs. |d_{0}|;|d_{0}|[cm];#LT #chi^{2} probability#GT", 100, 0, 80));
642  vTrackProfiles_.push_back(tfd.make<TProfile>(
643  "p_normchi2_vs_phi", "#chi^{2}/ndof vs. #phi;#phi_{Track};#LT #chi^{2}/ndof #GT", 100, -3.15, 3.15));
644  vTrackProfiles_.push_back(
645  tfd.make<TProfile>("p_chi2_vs_eta", "#chi^{2} vs. #eta;#eta_{Track};#LT #chi^{2} #GT", 100, -3.15, 3.15));
646  //variable binning for chi2/ndof vs. pT
647  double xBins[19] = {0., 0.15, 0.5, 1., 1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5., 7., 10., 15., 25., 40., 100., 200.};
648  vTrackProfiles_.push_back(tfd.make<TProfile>(
649  "p_normchi2_vs_pt", "norm #chi^{2} vs. p_{T}_{Track}; p_{T}_{Track};#LT #chi^{2}/ndof #GT", 18, xBins));
650 
651  vTrackProfiles_.push_back(
652  tfd.make<TProfile>("p_normchi2_vs_p", "#chi^{2}/ndof vs. p_{Track};p_{Track};#LT #chi^{2}/ndof #GT", 18, xBins));
653  vTrackProfiles_.push_back(tfd.make<TProfile>(
654  "p_chi2Prob_vs_eta", "#chi^{2} probability vs. #eta;#eta_{Track};#LT #chi^{2} probability #GT", 100, -3.15, 3.15));
655  vTrackProfiles_.push_back(tfd.make<TProfile>(
656  "p_normchi2_vs_eta", "#chi^{2}/ndof vs. #eta;#eta_{Track};#LT #chi^{2}/ndof #GT", 100, -3.15, 3.15));
657  vTrackProfiles_.push_back(
658  tfd.make<TProfile>("p_kappa_vs_phi", "#kappa vs. #phi;#phi_{Track};#kappa", 100, -3.15, 3.15));
659  vTrackProfiles_.push_back(
660  tfd.make<TProfile>("p_kappa_vs_eta", "#kappa vs. #eta;#eta_{Track};#kappa", 100, -3.15, 3.15));
661  vTrackProfiles_.push_back(tfd.make<TProfile>("p_ptResolution_vs_phi",
662  "#delta_{p_{T}}/p_{T}^{track};#phi^{track};#delta_{p_{T}}/p_{T}^{track}",
663  100,
664  -3.15,
665  3.15));
666  vTrackProfiles_.push_back(tfd.make<TProfile>("p_ptResolution_vs_eta",
667  "#delta_{p_{T}}/p_{T}^{track};#eta^{track};#delta_{p_{T}}/p_{T}^{track}",
668  100,
669  -3.15,
670  3.15));
671 
672  vTrack2DHistos_.push_back(tfd.make<TH2F>(
673  "h2_d0_vs_phi", "Transverse Impact Parameter vs. #phi;#phi_{Track};d_{0} [cm]", 100, -3.15, 3.15, 100, -1., 1.));
674  vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_dz_vs_phi",
675  "Longitudinal Impact Parameter vs. #phi;#phi_{Track};d_{z} [cm]",
676  100,
677  -3.15,
678  3.15,
679  100,
680  -100.,
681  100.));
682  vTrack2DHistos_.push_back(tfd.make<TH2F>(
683  "h2_d0_vs_eta", "Transverse Impact Parameter vs. #eta;#eta_{Track};d_{0} [cm]", 100, -3.15, 3.15, 100, -1., 1.));
684  vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_dz_vs_eta",
685  "Longitudinal Impact Parameter vs. #eta;#eta_{Track};d_{z} [cm]",
686  100,
687  -3.15,
688  3.15,
689  100,
690  -100.,
691  100.));
692  vTrack2DHistos_.push_back(
693  tfd.make<TH2F>("h2_chi2_vs_phi", "#chi^{2} vs. #phi;#phi_{Track};#chi^{2}", 100, -3.15, 3.15, 500, 0., 500.));
694  vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_chi2Prob_vs_phi",
695  "#chi^{2} probability vs. #phi;#phi_{Track};#chi^{2} probability",
696  100,
697  -3.15,
698  3.15,
699  100,
700  0.,
701  1.));
702  vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_chi2Prob_vs_d0",
703  "#chi^{2} probability vs. |d_{0}|;|d_{0}| [cm];#chi^{2} probability",
704  100,
705  0,
706  80,
707  100,
708  0.,
709  1.));
710  vTrack2DHistos_.push_back(tfd.make<TH2F>(
711  "h2_normchi2_vs_phi", "#chi^{2}/ndof vs. #phi;#phi_{Track};#chi^{2}/ndof", 100, -3.15, 3.15, 100, 0., 10.));
712  vTrack2DHistos_.push_back(
713  tfd.make<TH2F>("h2_chi2_vs_eta", "#chi^{2} vs. #eta;#eta_{Track};#chi^{2}", 100, -3.15, 3.15, 500, 0., 500.));
714  vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_chi2Prob_vs_eta",
715  "#chi^{2} probaility vs. #eta;#eta_{Track};#chi^{2} probability",
716  100,
717  -3.15,
718  3.15,
719  100,
720  0.,
721  1.));
722  vTrack2DHistos_.push_back(tfd.make<TH2F>(
723  "h2_normchi2_vs_eta", "#chi^{2}/ndof vs. #eta;#eta_{Track};#chi^{2}/ndof", 100, -3.15, 3.15, 100, 0., 10.));
724  vTrack2DHistos_.push_back(
725  tfd.make<TH2F>("h2_kappa_vs_phi", "#kappa vs. #phi;#phi_{Track};#kappa", 100, -3.15, 3.15, 100, .0, .05));
726  vTrack2DHistos_.push_back(
727  tfd.make<TH2F>("h2_kappa_vs_eta", "#kappa vs. #eta;#eta_{Track};#kappa", 100, -3.15, 3.15, 100, .0, .05));
728  vTrack2DHistos_.push_back(tfd.make<TH2F>(
729  "h2_normchi2_vs_kappa", "#kappa vs. #chi^{2}/ndof;#chi^{2}/ndof;#kappa", 100, 0., 10, 100, -.03, .03));
730 
731  /****************** Definition of 2-D Histos of ResX vs momenta ****************************/
732  vTrack2DHistos_.push_back(tfd.make<TH2F>(
733  "p_vs_resXprime_pixB", "res_{x'} vs momentum in BPix;p [GeV]; res_{x'}", 15, 0., 15., 200, -0.1, 0.1));
734  vTrack2DHistos_.push_back(tfd.make<TH2F>(
735  "p_vs_resXprime_pixE", "res_{x'} vs momentum in FPix;p [GeV]; res_{x'}", 15, 0., 15., 200, -0.1, 0.1));
736  vTrack2DHistos_.push_back(tfd.make<TH2F>(
737  "p_vs_resXprime_TIB", "res_{x'} vs momentum in TIB;p [GeV]; res_{x'}", 15, 0., 15., 200, -0.1, 0.1));
738  vTrack2DHistos_.push_back(tfd.make<TH2F>(
739  "p_vs_resXprime_TID", "res_{x'} vs momentum in TID;p [GeV]; res_{x'}", 15, 0., 15., 200, -0.1, 0.1));
740  vTrack2DHistos_.push_back(tfd.make<TH2F>(
741  "p_vs_resXprime_TOB", "res_{x'} vs momentum in TOB;p [GeV]; res_{x'}", 15, 0., 15., 200, -0.1, 0.1));
742  vTrack2DHistos_.push_back(tfd.make<TH2F>(
743  "p_vs_resXprime_TEC", "res_{x'} vs momentum in TEC;p [GeV]; res_{x'}", 15, 0., 15., 200, -0.1, 0.1));
744 
745  /****************** Definition of 2-D Histos of ResY vs momenta ****************************/
746  vTrack2DHistos_.push_back(tfd.make<TH2F>(
747  "p_vs_resYprime_pixB", "res_{y'} vs momentum in BPix;p [GeV]; res_{y'}", 15, 0., 15., 200, -0.1, 0.1));
748  vTrack2DHistos_.push_back(tfd.make<TH2F>(
749  "p_vs_resYprime_pixE", "res_{y'} vs momentum in FPix;p [GeV]; res_{y'}", 15, 0., 15., 200, -0.1, 0.1));
750 }
751 
753  const auto& alivec = ali.components();
754  for (int i = 0, iEnd = ali.components().size(); i < iEnd; ++i) {
755  std::string structurename = alignableTracker_->objectIdProvider().idToString((alivec)[i]->alignableObjectId());
756  LogDebug("TrackerOfflineValidation") << "StructureName = " << structurename;
757  std::stringstream dirname;
758  dirname << structurename;
759  // add no suffix counter to Strip and Pixel, just aesthetics
760  if (structurename != "Strip" && structurename != "Pixel")
761  dirname << "_" << i + 1;
762 
763  if (structurename.find("Endcap", 0) != std::string::npos) {
765  bookHists(f, *(alivec)[i], tTopo, ali.alignableObjectId(), i);
766  bookDirHists(f, *(alivec)[i], tTopo);
767  } else if (!(this->isDetOrDetUnit((alivec)[i]->alignableObjectId())) || alivec[i]->components().size() > 1) {
769  bookHists(tfd, *(alivec)[i], tTopo, ali.alignableObjectId(), i);
770  bookDirHists(f, *(alivec)[i], tTopo);
771  } else {
772  bookHists(tfd, *(alivec)[i], tTopo, ali.alignableObjectId(), i);
773  }
774  }
775 }
776 
778  DirectoryWrapper& tfd, const Alignable& ali, const TrackerTopology* tTopo, align::StructureType type, int i) {
779  TrackerAlignableId aliid;
780  const DetId id = ali.id();
781 
782  // comparing subdetandlayer to subdetIds gives a warning at compile time
783  // -> subdetandlayer could also be pair<uint,uint> but this has to be adapted
784  // in AlignableObjId
785  std::pair<int, int> subdetandlayer = aliid.typeAndLayerFromDetId(id, tTopo);
786 
788 
789  // are we on or just above det, detunit level respectively?
791  subtype = type;
792  else if (this->isDetOrDetUnit(ali.alignableObjectId()))
793  subtype = ali.alignableObjectId();
794 
795  // construct histogram title and name
796  std::stringstream histoname, histotitle, normhistoname, normhistotitle, yhistoname, yhistotitle, xprimehistoname,
797  xprimehistotitle, normxprimehistoname, normxprimehistotitle, yprimehistoname, yprimehistotitle,
798  normyprimehistoname, normyprimehistotitle, localxname, localxtitle, localyname, localytitle, resxvsxprofilename,
799  resxvsxprofiletitle, resyvsxprofilename, resyvsxprofiletitle, resxvsyprofilename, resxvsyprofiletitle,
800  resyvsyprofilename, resyvsyprofiletitle;
801 
802  std::string wheel_or_layer;
803 
804  if (this->isEndCap(static_cast<uint32_t>(subdetandlayer.first)))
805  wheel_or_layer = "_wheel_";
806  else if (this->isBarrel(static_cast<uint32_t>(subdetandlayer.first)))
807  wheel_or_layer = "_layer_";
808  else
809  edm::LogWarning("TrackerOfflineValidation") << "@SUB=TrackerOfflineValidation::bookHists"
810  << "Unknown subdetid: " << subdetandlayer.first;
811 
812  histoname << "h_residuals_subdet_" << subdetandlayer.first << wheel_or_layer << subdetandlayer.second << "_module_"
813  << id.rawId();
814  yhistoname << "h_y_residuals_subdet_" << subdetandlayer.first << wheel_or_layer << subdetandlayer.second << "_module_"
815  << id.rawId();
816  xprimehistoname << "h_xprime_residuals_subdet_" << subdetandlayer.first << wheel_or_layer << subdetandlayer.second
817  << "_module_" << id.rawId();
818  yprimehistoname << "h_yprime_residuals_subdet_" << subdetandlayer.first << wheel_or_layer << subdetandlayer.second
819  << "_module_" << id.rawId();
820  normhistoname << "h_normresiduals_subdet_" << subdetandlayer.first << wheel_or_layer << subdetandlayer.second
821  << "_module_" << id.rawId();
822  normxprimehistoname << "h_normxprimeresiduals_subdet_" << subdetandlayer.first << wheel_or_layer
823  << subdetandlayer.second << "_module_" << id.rawId();
824  normyprimehistoname << "h_normyprimeresiduals_subdet_" << subdetandlayer.first << wheel_or_layer
825  << subdetandlayer.second << "_module_" << id.rawId();
826  histotitle << "X Residual for module " << id.rawId() << ";x_{tr} - x_{hit} [cm]";
827  yhistotitle << "Y Residual for module " << id.rawId() << ";y_{tr} - y_{hit} [cm]";
828  normhistotitle << "Normalized Residual for module " << id.rawId() << ";x_{tr} - x_{hit}/#sigma";
829  xprimehistotitle << "X' Residual for module " << id.rawId() << ";(x_{tr} - x_{hit})' [cm]";
830  normxprimehistotitle << "Normalized X' Residual for module " << id.rawId() << ";(x_{tr} - x_{hit})'/#sigma";
831  yprimehistotitle << "Y' Residual for module " << id.rawId() << ";(y_{tr} - y_{hit})' [cm]";
832  normyprimehistotitle << "Normalized Y' Residual for module " << id.rawId() << ";(y_{tr} - y_{hit})'/#sigma";
833 
834  if (moduleLevelProfiles_) {
835  localxname << "h_localx_subdet_" << subdetandlayer.first << wheel_or_layer << subdetandlayer.second << "_module_"
836  << id.rawId();
837  localyname << "h_localy_subdet_" << subdetandlayer.first << wheel_or_layer << subdetandlayer.second << "_module_"
838  << id.rawId();
839  localxtitle << "u local for module " << id.rawId() << "; u_{tr,r}";
840  localytitle << "v local for module " << id.rawId() << "; v_{tr,r}";
841 
842  resxvsxprofilename << "p_residuals_x_vs_x_subdet_" << subdetandlayer.first << wheel_or_layer
843  << subdetandlayer.second << "_module_" << id.rawId();
844  resyvsxprofilename << "p_residuals_y_vs_x_subdet_" << subdetandlayer.first << wheel_or_layer
845  << subdetandlayer.second << "_module_" << id.rawId();
846  resxvsyprofilename << "p_residuals_x_vs_y_subdet_" << subdetandlayer.first << wheel_or_layer
847  << subdetandlayer.second << "_module_" << id.rawId();
848  resyvsyprofilename << "p_residuals_y_vs_y_subdet_" << subdetandlayer.first << wheel_or_layer
849  << subdetandlayer.second << "_module_" << id.rawId();
850  resxvsxprofiletitle << "U Residual vs u for module " << id.rawId()
851  << "; u_{tr,r} ;(u_{tr} - u_{hit})/tan#alpha [cm]";
852  resyvsxprofiletitle << "V Residual vs u for module " << id.rawId()
853  << "; u_{tr,r} ;(v_{tr} - v_{hit})/tan#beta [cm]";
854  resxvsyprofiletitle << "U Residual vs v for module " << id.rawId()
855  << "; v_{tr,r} ;(u_{tr} - u_{hit})/tan#alpha [cm]";
856  resyvsyprofiletitle << "V Residual vs v for module " << id.rawId()
857  << "; v_{tr,r} ;(v_{tr} - v_{hit})/tan#beta [cm]";
858  }
859 
860  if (this->isDetOrDetUnit(subtype)) {
861  ModuleHistos& histStruct = this->getHistStructFromMap(id);
862  int nbins = 0;
863  double xmin = 0., xmax = 0.;
864  double ymin = -0.1, ymax = 0.1;
865 
866  // do not allow transient hists in DQM mode
868  if (dqmMode_)
869  moduleLevelHistsTransient = false;
870 
871  // decide via cfg if hists in local coordinates should be booked
872  if (lCoorHistOn_) {
873  this->getBinning(id.subdetId(), XResidual, nbins, xmin, xmax);
874  histStruct.ResHisto = this->bookTH1F(
875  moduleLevelHistsTransient, tfd, histoname.str().c_str(), histotitle.str().c_str(), nbins, xmin, xmax);
876  this->getBinning(id.subdetId(), NormXResidual, nbins, xmin, xmax);
877  histStruct.NormResHisto = this->bookTH1F(
878  moduleLevelHistsTransient, tfd, normhistoname.str().c_str(), normhistotitle.str().c_str(), nbins, xmin, xmax);
879  }
880  this->getBinning(id.subdetId(), XprimeResidual, nbins, xmin, xmax);
881  histStruct.ResXprimeHisto = this->bookTH1F(moduleLevelHistsTransient,
882  tfd,
883  xprimehistoname.str().c_str(),
884  xprimehistotitle.str().c_str(),
885  nbins,
886  xmin,
887  xmax);
888  this->getBinning(id.subdetId(), NormXprimeResidual, nbins, xmin, xmax);
889  histStruct.NormResXprimeHisto = this->bookTH1F(moduleLevelHistsTransient,
890  tfd,
891  normxprimehistoname.str().c_str(),
892  normxprimehistotitle.str().c_str(),
893  nbins,
894  xmin,
895  xmax);
896 
897  if (moduleLevelProfiles_) {
898  this->getBinning(id.subdetId(), XResidualProfile, nbins, xmin, xmax);
899 
900  histStruct.LocalX = this->bookTH1F(
901  moduleLevelHistsTransient, tfd, localxname.str().c_str(), localxtitle.str().c_str(), nbins, xmin, xmax);
902  histStruct.LocalY = this->bookTH1F(
903  moduleLevelHistsTransient, tfd, localyname.str().c_str(), localytitle.str().c_str(), nbins, xmin, xmax);
904  histStruct.ResXvsXProfile = this->bookTProfile(moduleLevelHistsTransient,
905  tfd,
906  resxvsxprofilename.str().c_str(),
907  resxvsxprofiletitle.str().c_str(),
908  nbins,
909  xmin,
910  xmax,
911  ymin,
912  ymax);
913  histStruct.ResXvsXProfile->Sumw2(); // to be filled with weights, so uncertainties need sum of square of weights
914  histStruct.ResXvsYProfile = this->bookTProfile(moduleLevelHistsTransient,
915  tfd,
916  resxvsyprofilename.str().c_str(),
917  resxvsyprofiletitle.str().c_str(),
918  nbins,
919  xmin,
920  xmax,
921  ymin,
922  ymax);
923  histStruct.ResXvsYProfile->Sumw2(); // to be filled with weights, so uncertainties need sum of square of weights
924  }
925 
926  if (this->isPixel(subdetandlayer.first) || stripYResiduals_) {
927  this->getBinning(id.subdetId(), YprimeResidual, nbins, xmin, xmax);
928  histStruct.ResYprimeHisto = this->bookTH1F(moduleLevelHistsTransient,
929  tfd,
930  yprimehistoname.str().c_str(),
931  yprimehistotitle.str().c_str(),
932  nbins,
933  xmin,
934  xmax);
935  if (lCoorHistOn_) { // un-primed y-residual
936  this->getBinning(id.subdetId(), YResidual, nbins, xmin, xmax);
937  histStruct.ResYHisto = this->bookTH1F(
938  moduleLevelHistsTransient, tfd, yhistoname.str().c_str(), yhistotitle.str().c_str(), nbins, xmin, xmax);
939  }
940  this->getBinning(id.subdetId(), NormYprimeResidual, nbins, xmin, xmax);
941  histStruct.NormResYprimeHisto = this->bookTH1F(moduleLevelHistsTransient,
942  tfd,
943  normyprimehistoname.str().c_str(),
944  normyprimehistotitle.str().c_str(),
945  nbins,
946  xmin,
947  xmax);
948  // Here we could add un-primed normalised y-residuals if(lCoorHistOn_)...
949  if (moduleLevelProfiles_) {
950  this->getBinning(id.subdetId(), YResidualProfile, nbins, xmin, xmax);
951 
952  histStruct.ResYvsXProfile = this->bookTProfile(moduleLevelHistsTransient,
953  tfd,
954  resyvsxprofilename.str().c_str(),
955  resyvsxprofiletitle.str().c_str(),
956  nbins,
957  xmin,
958  xmax,
959  ymin,
960  ymax);
961  histStruct.ResYvsXProfile->Sumw2(); // to be filled with weights, so uncertainties need sum of square of weights
962  histStruct.ResYvsYProfile = this->bookTProfile(moduleLevelHistsTransient,
963  tfd,
964  resyvsyprofilename.str().c_str(),
965  resyvsyprofiletitle.str().c_str(),
966  nbins,
967  xmin,
968  xmax,
969  ymin,
970  ymax);
971  histStruct.ResYvsYProfile->Sumw2(); // to be filled with weights, so uncertainties need sum of square of weights
972  }
973  }
974  }
975 }
976 
978  DirectoryWrapper& tfd,
979  const char* histName,
980  const char* histTitle,
981  int nBinsX,
982  double lowX,
983  double highX) {
984  if (isTransient) {
985  vDeleteObjects_.push_back(new TH1F(histName, histTitle, nBinsX, lowX, highX));
986  return vDeleteObjects_.back(); // return last element of vector
987  } else
988  return tfd.make<TH1F>(histName, histTitle, nBinsX, lowX, highX);
989 }
990 
991 TProfile* TrackerOfflineValidation::bookTProfile(bool isTransient,
992  DirectoryWrapper& tfd,
993  const char* histName,
994  const char* histTitle,
995  int nBinsX,
996  double lowX,
997  double highX) {
998  if (isTransient) {
999  TProfile* profile = new TProfile(histName, histTitle, nBinsX, lowX, highX);
1000  vDeleteObjects_.push_back(profile);
1001  return profile;
1002  } else
1003  return (TProfile*)tfd.make<TProfile>(histName, histTitle, nBinsX, lowX, highX);
1004 }
1005 
1006 TProfile* TrackerOfflineValidation::bookTProfile(bool isTransient,
1007  DirectoryWrapper& tfd,
1008  const char* histName,
1009  const char* histTitle,
1010  int nBinsX,
1011  double lowX,
1012  double highX,
1013  double lowY,
1014  double highY) {
1015  if (isTransient) {
1016  TProfile* profile = new TProfile(histName, histTitle, nBinsX, lowX, highX, lowY, highY);
1017  vDeleteObjects_.push_back(profile);
1018  return profile;
1019  } else
1020  return (TProfile*)tfd.make<TProfile>(histName, histTitle, nBinsX, lowX, highX, lowY, highY);
1021 }
1022 
1026 }
1027 
1031 }
1032 
1035 }
1036 
1039 }
1040 
1043  int& nBinsX,
1044  double& lowerBoundX,
1045  double& upperBoundX) {
1046  // determine if
1047  const bool isPixel = this->isPixel(subDetId);
1048 
1049  edm::ParameterSet binningPSet;
1050 
1051  switch (residualType) {
1052  case XResidual:
1053  if (isPixel)
1054  binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1XResPixelModules");
1055  else
1056  binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1XResStripModules");
1057  break;
1058  case NormXResidual:
1059  if (isPixel)
1060  binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1NormXResPixelModules");
1061  else
1062  binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1NormXResStripModules");
1063  break;
1064  case XprimeResidual:
1065  if (isPixel)
1066  binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1XprimeResPixelModules");
1067  else
1068  binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1XprimeResStripModules");
1069  break;
1070  case NormXprimeResidual:
1071  if (isPixel)
1072  binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1NormXprimeResPixelModules");
1073  else
1074  binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1NormXprimeResStripModules");
1075  break;
1076  case YResidual: // borrow y-residual binning from yprime
1077  case YprimeResidual:
1078  if (isPixel)
1079  binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1YResPixelModules");
1080  else
1081  binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1YResStripModules");
1082  break;
1083  /* case NormYResidual :*/
1084  case NormYprimeResidual:
1085  if (isPixel)
1086  binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1NormYResPixelModules");
1087  else
1088  binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1NormYResStripModules");
1089  break;
1090  case XResidualProfile:
1091  if (isPixel)
1092  binningPSet = parSet_.getParameter<edm::ParameterSet>("TProfileXResPixelModules");
1093  else
1094  binningPSet = parSet_.getParameter<edm::ParameterSet>("TProfileXResStripModules");
1095  break;
1096  case YResidualProfile:
1097  if (isPixel)
1098  binningPSet = parSet_.getParameter<edm::ParameterSet>("TProfileYResPixelModules");
1099  else
1100  binningPSet = parSet_.getParameter<edm::ParameterSet>("TProfileYResStripModules");
1101  break;
1102  }
1103  nBinsX = binningPSet.getParameter<int32_t>("Nbinx");
1104  lowerBoundX = binningPSet.getParameter<double>("xmin");
1105  upperBoundX = binningPSet.getParameter<double>("xmax");
1106 }
1107 
1108 void TrackerOfflineValidation::setSummaryBin(int bin, TH1* targetHist, TH1* sourceHist) {
1109  if (targetHist && sourceHist) {
1110  targetHist->SetBinContent(bin, sourceHist->GetMean(1));
1111  if (useFwhm_)
1112  targetHist->SetBinError(bin, Fwhm(sourceHist) / 2.);
1113  else
1114  targetHist->SetBinError(bin, sourceHist->GetRMS(1));
1115  } else
1116  return;
1117 }
1118 
1120  SummaryContainer& targetContainer,
1121  SummaryContainer& sourceContainer) {
1122  summaryBins_.emplace_back(bin, targetContainer.summaryXResiduals_, sourceContainer.sumXResiduals_);
1123  summaryBins_.emplace_back(bin, targetContainer.summaryNormXResiduals_, sourceContainer.sumNormXResiduals_);
1124  // If no y-residual hists, just returns:
1125  summaryBins_.emplace_back(bin, targetContainer.summaryYResiduals_, sourceContainer.sumYResiduals_);
1126  summaryBins_.emplace_back(bin, targetContainer.summaryNormYResiduals_, sourceContainer.sumNormYResiduals_);
1127 }
1128 
1130  uint32_t subDetId,
1131  SummaryContainer& targetContainer,
1132  ModuleHistos& sourceContainer) {
1133  // takes two summary Containers and sets summaryBins for all histograms
1134  summaryBins_.emplace_back(bin, targetContainer.summaryXResiduals_, sourceContainer.ResXprimeHisto);
1135  summaryBins_.emplace_back(bin, targetContainer.summaryNormXResiduals_, sourceContainer.NormResXprimeHisto);
1136  if (this->isPixel(subDetId) || stripYResiduals_) {
1137  summaryBins_.emplace_back(bin, targetContainer.summaryYResiduals_, sourceContainer.ResYprimeHisto);
1138  summaryBins_.emplace_back(bin, targetContainer.summaryNormYResiduals_, sourceContainer.NormResYprimeHisto);
1139  }
1140 }
1141 
1143  // get a struct with histograms from the respective map
1144  // if no object exist, the reference is automatically created by the map
1145  // throw exception if non-tracker id is passed
1146  uint subdetid = detid.subdetId();
1147  if (subdetid == PixelSubdetector::PixelBarrel) {
1148  return mPxbResiduals_[detid.rawId()];
1149  } else if (subdetid == PixelSubdetector::PixelEndcap) {
1150  return mPxeResiduals_[detid.rawId()];
1151  } else if (subdetid == StripSubdetector::TIB) {
1152  return mTibResiduals_[detid.rawId()];
1153  } else if (subdetid == StripSubdetector::TID) {
1154  return mTidResiduals_[detid.rawId()];
1155  } else if (subdetid == StripSubdetector::TOB) {
1156  return mTobResiduals_[detid.rawId()];
1157  } else if (subdetid == StripSubdetector::TEC) {
1158  return mTecResiduals_[detid.rawId()];
1159  } else {
1160  throw cms::Exception("Geometry Error")
1161  << "[TrackerOfflineValidation] Error, tried to get reference for non-tracker subdet " << subdetid
1162  << " from detector " << detid.det();
1163  return mPxbResiduals_[0];
1164  }
1165 }
1166 
1167 // ------------ method called to for each event ------------
1169  if (maxTracks_ > 0 && nTracks_ >= maxTracks_)
1170  return; //don't do anything after hitting the max number of tracks
1171  //(could just rely on break below, but this way saves fillTrackQuantities)
1172 
1173  if (useOverflowForRMS_)
1174  TH1::StatOverflows(kTRUE);
1175  this->checkBookHists(iSetup); // check whether hists are booked and do so if not yet done
1176 
1177  std::vector<TrackerValidationVariables::AVTrackStruct> vTrackstruct;
1178  avalidator_.fillTrackQuantities(iEvent, iSetup, vTrackstruct);
1179 
1180  for (std::vector<TrackerValidationVariables::AVTrackStruct>::const_iterator itT = vTrackstruct.begin();
1181  itT != vTrackstruct.end();
1182  ++itT) {
1183  if (itT->charge * chargeCut_ < 0)
1184  continue;
1185 
1186  if (maxTracks_ > 0 && nTracks_++ >= maxTracks_)
1187  break; //exit the loop after hitting the max number of tracks
1188  // Fill 1D track histos
1189  static const int etaindex = this->GetIndex(vTrackHistos_, "h_tracketa");
1190  vTrackHistos_[etaindex]->Fill(itT->eta);
1191  static const int phiindex = this->GetIndex(vTrackHistos_, "h_trackphi");
1192  vTrackHistos_[phiindex]->Fill(itT->phi);
1193  static const int numOfValidHitsindex = this->GetIndex(vTrackHistos_, "h_trackNumberOfValidHits");
1194  vTrackHistos_[numOfValidHitsindex]->Fill(itT->numberOfValidHits);
1195  static const int numOfLostHitsindex = this->GetIndex(vTrackHistos_, "h_trackNumberOfLostHits");
1196  vTrackHistos_[numOfLostHitsindex]->Fill(itT->numberOfLostHits);
1197  static const int kappaindex = this->GetIndex(vTrackHistos_, "h_curvature");
1198  vTrackHistos_[kappaindex]->Fill(itT->kappa);
1199  static const int kappaposindex = this->GetIndex(vTrackHistos_, "h_curvature_pos");
1200  if (itT->charge > 0)
1201  vTrackHistos_[kappaposindex]->Fill(fabs(itT->kappa));
1202  static const int kappanegindex = this->GetIndex(vTrackHistos_, "h_curvature_neg");
1203  if (itT->charge < 0)
1204  vTrackHistos_[kappanegindex]->Fill(fabs(itT->kappa));
1205  static const int normchi2index = this->GetIndex(vTrackHistos_, "h_normchi2");
1206  vTrackHistos_[normchi2index]->Fill(itT->normchi2);
1207  static const int chi2index = this->GetIndex(vTrackHistos_, "h_chi2");
1208  vTrackHistos_[chi2index]->Fill(itT->chi2);
1209  static const int chi2Probindex = this->GetIndex(vTrackHistos_, "h_chi2Prob");
1210  vTrackHistos_[chi2Probindex]->Fill(itT->chi2Prob);
1211  static const int ptindex = this->GetIndex(vTrackHistos_, "h_pt");
1212  vTrackHistos_[ptindex]->Fill(itT->pt);
1213  if (itT->ptError != 0.) {
1214  static const int ptResolutionindex = this->GetIndex(vTrackHistos_, "h_ptResolution");
1215  vTrackHistos_[ptResolutionindex]->Fill(itT->ptError / itT->pt);
1216  }
1217  // Fill track profiles
1218  static const int d0phiindex = this->GetIndex(vTrackProfiles_, "p_d0_vs_phi");
1219  vTrackProfiles_[d0phiindex]->Fill(itT->phi, itT->d0);
1220  static const int dzphiindex = this->GetIndex(vTrackProfiles_, "p_dz_vs_phi");
1221  vTrackProfiles_[dzphiindex]->Fill(itT->phi, itT->dz);
1222  static const int d0etaindex = this->GetIndex(vTrackProfiles_, "p_d0_vs_eta");
1223  vTrackProfiles_[d0etaindex]->Fill(itT->eta, itT->d0);
1224  static const int dzetaindex = this->GetIndex(vTrackProfiles_, "p_dz_vs_eta");
1225  vTrackProfiles_[dzetaindex]->Fill(itT->eta, itT->dz);
1226  static const int chiphiindex = this->GetIndex(vTrackProfiles_, "p_chi2_vs_phi");
1227  vTrackProfiles_[chiphiindex]->Fill(itT->phi, itT->chi2);
1228  static const int chiProbphiindex = this->GetIndex(vTrackProfiles_, "p_chi2Prob_vs_phi");
1229  vTrackProfiles_[chiProbphiindex]->Fill(itT->phi, itT->chi2Prob);
1230  static const int chiProbabsd0index = this->GetIndex(vTrackProfiles_, "p_chi2Prob_vs_d0");
1231  vTrackProfiles_[chiProbabsd0index]->Fill(fabs(itT->d0), itT->chi2Prob);
1232  static const int normchiphiindex = this->GetIndex(vTrackProfiles_, "p_normchi2_vs_phi");
1233  vTrackProfiles_[normchiphiindex]->Fill(itT->phi, itT->normchi2);
1234  static const int chietaindex = this->GetIndex(vTrackProfiles_, "p_chi2_vs_eta");
1235  vTrackProfiles_[chietaindex]->Fill(itT->eta, itT->chi2);
1236  static const int normchiptindex = this->GetIndex(vTrackProfiles_, "p_normchi2_vs_pt");
1237  vTrackProfiles_[normchiptindex]->Fill(itT->pt, itT->normchi2);
1238  static const int normchipindex = this->GetIndex(vTrackProfiles_, "p_normchi2_vs_p");
1239  vTrackProfiles_[normchipindex]->Fill(itT->p, itT->normchi2);
1240  static const int chiProbetaindex = this->GetIndex(vTrackProfiles_, "p_chi2Prob_vs_eta");
1241  vTrackProfiles_[chiProbetaindex]->Fill(itT->eta, itT->chi2Prob);
1242  static const int normchietaindex = this->GetIndex(vTrackProfiles_, "p_normchi2_vs_eta");
1243  vTrackProfiles_[normchietaindex]->Fill(itT->eta, itT->normchi2);
1244  static const int kappaphiindex = this->GetIndex(vTrackProfiles_, "p_kappa_vs_phi");
1245  vTrackProfiles_[kappaphiindex]->Fill(itT->phi, itT->kappa);
1246  static const int kappaetaindex = this->GetIndex(vTrackProfiles_, "p_kappa_vs_eta");
1247  vTrackProfiles_[kappaetaindex]->Fill(itT->eta, itT->kappa);
1248  static const int ptResphiindex = this->GetIndex(vTrackProfiles_, "p_ptResolution_vs_phi");
1249  vTrackProfiles_[ptResphiindex]->Fill(itT->phi, itT->ptError / itT->pt);
1250  static const int ptResetaindex = this->GetIndex(vTrackProfiles_, "p_ptResolution_vs_eta");
1251  vTrackProfiles_[ptResetaindex]->Fill(itT->eta, itT->ptError / itT->pt);
1252 
1253  // Fill 2D track histos
1254  static const int d0phiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_d0_vs_phi");
1255  vTrack2DHistos_[d0phiindex_2d]->Fill(itT->phi, itT->d0);
1256  static const int dzphiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_dz_vs_phi");
1257  vTrack2DHistos_[dzphiindex_2d]->Fill(itT->phi, itT->dz);
1258  static const int d0etaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_d0_vs_eta");
1259  vTrack2DHistos_[d0etaindex_2d]->Fill(itT->eta, itT->d0);
1260  static const int dzetaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_dz_vs_eta");
1261  vTrack2DHistos_[dzetaindex_2d]->Fill(itT->eta, itT->dz);
1262  static const int chiphiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_chi2_vs_phi");
1263  vTrack2DHistos_[chiphiindex_2d]->Fill(itT->phi, itT->chi2);
1264  static const int chiProbphiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_chi2Prob_vs_phi");
1265  vTrack2DHistos_[chiProbphiindex_2d]->Fill(itT->phi, itT->chi2Prob);
1266  static const int chiProbabsd0index_2d = this->GetIndex(vTrack2DHistos_, "h2_chi2Prob_vs_d0");
1267  vTrack2DHistos_[chiProbabsd0index_2d]->Fill(fabs(itT->d0), itT->chi2Prob);
1268  static const int normchiphiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_normchi2_vs_phi");
1269  vTrack2DHistos_[normchiphiindex_2d]->Fill(itT->phi, itT->normchi2);
1270  static const int chietaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_chi2_vs_eta");
1271  vTrack2DHistos_[chietaindex_2d]->Fill(itT->eta, itT->chi2);
1272  static const int chiProbetaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_chi2Prob_vs_eta");
1273  vTrack2DHistos_[chiProbetaindex_2d]->Fill(itT->eta, itT->chi2Prob);
1274  static const int normchietaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_normchi2_vs_eta");
1275  vTrack2DHistos_[normchietaindex_2d]->Fill(itT->eta, itT->normchi2);
1276  static const int kappaphiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_kappa_vs_phi");
1277  vTrack2DHistos_[kappaphiindex_2d]->Fill(itT->phi, itT->kappa);
1278  static const int kappaetaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_kappa_vs_eta");
1279  vTrack2DHistos_[kappaetaindex_2d]->Fill(itT->eta, itT->kappa);
1280  static const int normchi2kappa_2d = this->GetIndex(vTrack2DHistos_, "h2_normchi2_vs_kappa");
1281  vTrack2DHistos_[normchi2kappa_2d]->Fill(itT->normchi2, itT->kappa);
1282 
1283  // hit quantities: residuals, normalized residuals
1284  for (std::vector<TrackerValidationVariables::AVHitStruct>::const_iterator itH = itT->hits.begin();
1285  itH != itT->hits.end();
1286  ++itH) {
1287  DetId detid(itH->rawDetId);
1288  ModuleHistos& histStruct = this->getHistStructFromMap(detid);
1289 
1290  // fill histos in local coordinates if set in cf
1291  if (lCoorHistOn_) {
1292  histStruct.ResHisto->Fill(itH->resX);
1293  if (itH->resErrX != 0)
1294  histStruct.NormResHisto->Fill(itH->resX / itH->resErrX);
1295  if (this->isPixel(detid.subdetId()) || stripYResiduals_) {
1296  histStruct.ResYHisto->Fill(itH->resY);
1297  // here add un-primed normalised y-residuals if wanted
1298  }
1299  }
1300  if (itH->resXprime != -999.) {
1301  histStruct.ResXprimeHisto->Fill(itH->resXprime);
1302 
1303  /******************************* Fill 2-D histo ResX vs momenta *****************************/
1304  if (detid.subdetId() == PixelSubdetector::PixelBarrel) {
1305  static const int resXvsPindex_2d = this->GetIndex(vTrack2DHistos_, "p_vs_resXprime_pixB");
1306  vTrack2DHistos_[resXvsPindex_2d]->Fill(itT->p, itH->resXprime);
1307  }
1308  if (detid.subdetId() == PixelSubdetector::PixelEndcap) {
1309  static const int resXvsPindex_2d = this->GetIndex(vTrack2DHistos_, "p_vs_resXprime_pixE");
1310  vTrack2DHistos_[resXvsPindex_2d]->Fill(itT->p, itH->resXprime);
1311  }
1312  if (detid.subdetId() == StripSubdetector::TIB) {
1313  static const int resXvsPindex_2d = this->GetIndex(vTrack2DHistos_, "p_vs_resXprime_TIB");
1314  vTrack2DHistos_[resXvsPindex_2d]->Fill(itT->p, itH->resXprime);
1315  }
1316  if (detid.subdetId() == StripSubdetector::TID) {
1317  static const int resXvsPindex_2d = this->GetIndex(vTrack2DHistos_, "p_vs_resXprime_TID");
1318  vTrack2DHistos_[resXvsPindex_2d]->Fill(itT->p, itH->resXprime);
1319  }
1320  if (detid.subdetId() == StripSubdetector::TOB) {
1321  static const int resXvsPindex_2d = this->GetIndex(vTrack2DHistos_, "p_vs_resXprime_TOB");
1322  vTrack2DHistos_[resXvsPindex_2d]->Fill(itT->p, itH->resXprime);
1323  }
1324  if (detid.subdetId() == StripSubdetector::TEC) {
1325  static const int resXvsPindex_2d = this->GetIndex(vTrack2DHistos_, "p_vs_resXprime_TEC");
1326  vTrack2DHistos_[resXvsPindex_2d]->Fill(itT->p, itH->resXprime);
1327  }
1328  /******************************************/
1329 
1330  if (moduleLevelProfiles_ && itH->inside) {
1331  float tgalpha = tan(itH->localAlpha);
1332  histStruct.EntriesInt = histStruct.LocalX->GetEntries();
1334  continue;
1335  if (fabs(tgalpha) != 0) {
1336  histStruct.LocalX->Fill(itH->localXnorm, tgalpha * tgalpha);
1337  histStruct.LocalY->Fill(itH->localYnorm, tgalpha * tgalpha);
1338  /* if (this->isEndCap(detid.subdetId()) && !this->isPixel(detid.subdetId())) {
1339  if((itH->resX)*(itH->resXprime)>0){
1340  histStruct.ResXvsXProfile->Fill(itH->localXnorm, itH->resXatTrkY/tgalpha, tgalpha*tgalpha);
1341  histStruct.ResXvsYProfile->Fill(itH->localYnorm, itH->resXatTrkY/tgalpha, tgalpha*tgalpha);
1342  } else {
1343  histStruct.ResXvsXProfile->Fill(itH->localXnorm, (-1)*itH->resXatTrkY/tgalpha, tgalpha*tgalpha);
1344  histStruct.ResXvsYProfile->Fill(itH->localYnorm, (-1)*itH->resXatTrkY/tgalpha, tgalpha*tgalpha);
1345  }
1346 
1347  }else {
1348 */
1349  histStruct.ResXvsXProfile->Fill(itH->localXnorm, itH->resXatTrkY / tgalpha, tgalpha * tgalpha);
1350  histStruct.ResXvsYProfile->Fill(itH->localYnorm, itH->resXatTrkY / tgalpha, tgalpha * tgalpha);
1351 
1352  // }
1353  }
1354  }
1355 
1356  if (itH->resXprimeErr != 0 && itH->resXprimeErr != -999) {
1357  histStruct.NormResXprimeHisto->Fill(itH->resXprime / itH->resXprimeErr);
1358  }
1359  }
1360 
1361  if (itH->resYprime != -999.) {
1362  if (this->isPixel(detid.subdetId()) || stripYResiduals_) {
1363  histStruct.ResYprimeHisto->Fill(itH->resYprime);
1364 
1365  /******************************* Fill 2-D histo ResY vs momenta *****************************/
1366  if (detid.subdetId() == PixelSubdetector::PixelBarrel) {
1367  static const int resYvsPindex_2d = this->GetIndex(vTrack2DHistos_, "p_vs_resYprime_pixB");
1368  vTrack2DHistos_[resYvsPindex_2d]->Fill(itT->p, itH->resYprime);
1369  }
1370  if (detid.subdetId() == PixelSubdetector::PixelEndcap) {
1371  static const int resYvsPindex_2d = this->GetIndex(vTrack2DHistos_, "p_vs_resYprime_pixE");
1372  vTrack2DHistos_[resYvsPindex_2d]->Fill(itT->p, itH->resYprime);
1373  }
1374  /******************************************/
1375 
1376  if (moduleLevelProfiles_ && itH->inside) {
1377  float tgbeta = tan(itH->localBeta);
1378  if (fabs(tgbeta) != 0) {
1379  /* if (this->isEndCap(detid.subdetId()) && !this->isPixel(detid.subdetId())) {
1380 
1381  if((itH->resY)*(itH->resYprime)>0){
1382  histStruct.ResYvsXProfile->Fill(itH->localXnorm, itH->resYprime/tgbeta, tgbeta*tgbeta);
1383  histStruct.ResYvsYProfile->Fill(itH->localYnorm, itH->resYprime/tgbeta, tgbeta*tgbeta);
1384  } else {
1385  histStruct.ResYvsXProfile->Fill(itH->localXnorm, (-1)*itH->resYprime/tgbeta, tgbeta*tgbeta);
1386  histStruct.ResYvsYProfile->Fill(itH->localYnorm, (-1)*itH->resYprime/tgbeta, tgbeta*tgbeta);
1387  }
1388 
1389  }else {
1390 */
1391  histStruct.ResYvsXProfile->Fill(itH->localXnorm, itH->resY / tgbeta, tgbeta * tgbeta);
1392  histStruct.ResYvsYProfile->Fill(itH->localYnorm, itH->resY / tgbeta, tgbeta * tgbeta);
1393  // }
1394  }
1395  }
1396 
1397  if (itH->resYprimeErr != 0 && itH->resYprimeErr != -999.) {
1398  histStruct.NormResYprimeHisto->Fill(itH->resYprime / itH->resYprimeErr);
1399  }
1400  }
1401  }
1402 
1403  } // finish loop over hit quantities
1404  } // finish loop over track quantities
1405 
1406  if (useOverflowForRMS_)
1407  TH1::StatOverflows(kFALSE);
1408 }
1409 
1410 // ------------ method called once each job just after ending the event loop ------------
1412  if (!tkGeom_.product())
1413  return;
1414 
1415  static const int kappadiffindex = this->GetIndex(vTrackHistos_, "h_diff_curvature");
1416  vTrackHistos_[kappadiffindex]->Add(vTrackHistos_[this->GetIndex(vTrackHistos_, "h_curvature_neg")],
1417  vTrackHistos_[this->GetIndex(vTrackHistos_, "h_curvature_pos")],
1418  -1,
1419  1);
1420 
1421  // Collate Information for Subdetectors
1422  // create summary histograms recursively
1423  this->collateSummaryHists();
1424 
1425  if (dqmMode_)
1426  return;
1427  // Should be excluded in dqmMode, since TTree is not usable
1428  // In dqmMode tree operations are are sourced out to the additional module TrackerOfflineValidationSummary
1429 
1431  if (compressionSettings_ > 0) {
1432  fs->file().SetCompressionSettings(compressionSettings_);
1433  }
1434 
1435  TTree* tree = fs->make<TTree>("TkOffVal", "TkOffVal");
1436 
1437  TkOffTreeVariables* treeMemPtr = new TkOffTreeVariables;
1438  // We create branches for all members of 'TkOffTreeVariables' (even if not needed).
1439  // This works because we have a dictionary for 'TkOffTreeVariables'
1440  // (see src/classes_def.xml and src/classes.h):
1441  tree->Branch("TkOffTreeVariables", &treeMemPtr); // address of pointer!
1442 
1443  this->fillTree(*tree, *treeMemPtr, mPxbResiduals_);
1444  this->fillTree(*tree, *treeMemPtr, mPxeResiduals_);
1445  this->fillTree(*tree, *treeMemPtr, mTibResiduals_);
1446  this->fillTree(*tree, *treeMemPtr, mTidResiduals_);
1447  this->fillTree(*tree, *treeMemPtr, mTobResiduals_);
1448  this->fillTree(*tree, *treeMemPtr, mTecResiduals_);
1449 
1450  delete treeMemPtr;
1451  treeMemPtr = nullptr;
1452 }
1453 
1455  DirectoryWrapper& tfd,
1456  const Alignable& ali,
1457  std::vector<TrackerOfflineValidation::SummaryContainer>& vLevelProfiles) {
1458  const auto& alivec = ali.components();
1459  if (this->isDetOrDetUnit((alivec)[0]->alignableObjectId()))
1460  return;
1461 
1462  for (int iComp = 0, iCompEnd = ali.components().size(); iComp < iCompEnd; ++iComp) {
1463  std::vector<TrackerOfflineValidation::SummaryContainer> vProfiles;
1464  std::string structurename = alignableTracker_->objectIdProvider().idToString((alivec)[iComp]->alignableObjectId());
1465 
1466  LogDebug("TrackerOfflineValidation") << "StructureName = " << structurename;
1467  std::stringstream dirname;
1468  dirname << structurename;
1469 
1470  // add no suffix counter to strip and pixel -> just aesthetics
1471  if (structurename != "Strip" && structurename != "Pixel")
1472  dirname << "_" << iComp + 1;
1473 
1474  if (!(this->isDetOrDetUnit((alivec)[iComp]->alignableObjectId())) || (alivec)[0]->components().size() > 1) {
1476  this->prepareSummaryHists(f, *(alivec)[iComp], vProfiles);
1477  vLevelProfiles.push_back(this->bookSummaryHists(tfd, *(alivec[iComp]), ali.alignableObjectId(), iComp + 1));
1478  TH1* hX = vLevelProfiles[iComp].sumXResiduals_;
1479  TH1* hNormX = vLevelProfiles[iComp].sumNormXResiduals_;
1480  TH1* hY = vLevelProfiles[iComp].sumYResiduals_;
1481  TH1* hNormY = vLevelProfiles[iComp].sumNormYResiduals_;
1482  TH1* pXX = vLevelProfiles[iComp].sumResXvsXProfile_;
1483  TH1* pXY = vLevelProfiles[iComp].sumResXvsYProfile_;
1484  TH1* pYX = vLevelProfiles[iComp].sumResYvsXProfile_;
1485  TH1* pYY = vLevelProfiles[iComp].sumResYvsYProfile_;
1486  for (uint n = 0; n < vProfiles.size(); ++n) {
1487  this->summarizeBinInContainer(n + 1, vLevelProfiles[iComp], vProfiles[n]);
1488  sumHistStructure_.emplace_back(hX, vProfiles[n].sumXResiduals_);
1489  sumHistStructure_.emplace_back(hNormX, vProfiles[n].sumNormXResiduals_);
1490  if (hY)
1491  sumHistStructure_.emplace_back(hY, vProfiles[n].sumYResiduals_); // only if existing
1492  if (hNormY)
1493  sumHistStructure_.emplace_back(hNormY, vProfiles[n].sumNormYResiduals_); // ditto (pxl, stripYResiduals_)
1494  if (pXX)
1495  sumHistStructure_.emplace_back(pXX, vProfiles[n].sumResXvsXProfile_);
1496  if (pXY)
1497  sumHistStructure_.emplace_back(pXY, vProfiles[n].sumResXvsYProfile_);
1498  if (pYX)
1499  sumHistStructure_.emplace_back(pYX, vProfiles[n].sumResYvsXProfile_);
1500  if (pYY)
1501  sumHistStructure_.emplace_back(pYY, vProfiles[n].sumResYvsYProfile_);
1502  }
1503  if (dqmMode_)
1504  continue; // No fits in dqmMode
1505  //add fit values to stat box
1506  toFit_.push_back(vLevelProfiles[iComp].sumXResiduals_);
1507  toFit_.push_back(vLevelProfiles[iComp].sumNormXResiduals_);
1508  if (hY)
1509  toFit_.push_back(hY); // only if existing (pixel or stripYResiduals_)
1510  if (hNormY)
1511  toFit_.push_back(hNormY); // ditto
1512  } else {
1513  // nothing to be done for det or detunits
1514  continue;
1515  }
1516  }
1517 }
1518 
1520  for (std::vector<std::pair<TH1*, TH1*> >::const_iterator it = sumHistStructure_.begin();
1521  it != sumHistStructure_.end();
1522  ++it)
1523  it->first->Add(it->second);
1524 
1525  for (std::vector<std::tuple<int, TH1*, TH1*> >::const_iterator it = summaryBins_.begin(); it != summaryBins_.end();
1526  ++it)
1527  setSummaryBin(std::get<0>(*it), std::get<1>(*it), std::get<2>(*it));
1528 
1529  for (std::vector<TH1*>::const_iterator it = toFit_.begin(); it != toFit_.end(); ++it)
1530  fitResiduals(*it);
1531 }
1532 
1534  const Alignable& ali,
1536  int i) {
1537  const uint aliSize = ali.components().size();
1538  const align::StructureType alitype = ali.alignableObjectId();
1539  const align::StructureType subtype = ali.components()[0]->alignableObjectId();
1540  const char* aliTypeName = alignableTracker_->objectIdProvider().idToString(alitype); // lifetime of char* OK
1541  const char* aliSubtypeName = alignableTracker_->objectIdProvider().idToString(subtype);
1542  const char* typeName = alignableTracker_->objectIdProvider().idToString(type);
1543 
1544  const DetId aliDetId = ali.id();
1545  // y residuals only if pixel or specially requested for strip:
1546  const bool bookResidY = this->isPixel(aliDetId.subdetId()) || stripYResiduals_;
1547 
1548  SummaryContainer sumContainer;
1549 
1550  // Book summary hists with one bin per component,
1551  // but special case for Det with two DetUnit that we want to summarize one level up
1552  // (e.g. in TOBRods with 12 bins for 6 stereo and 6 rphi DetUnit.)
1553  // component of ali is not Det or Det with just one components
1554  const uint subcompSize = ali.components()[0]->components().size();
1555  if (subtype != align::AlignableDet || subcompSize == 1) { // Det with 1 comp. should not exist anymore...
1556  const TString title(Form("Summary for substructures in %s %d;%s;", aliTypeName, i, aliSubtypeName));
1557 
1558  sumContainer.summaryXResiduals_ = tfd.make<TH1F>(
1559  Form("h_summaryX%s_%d", aliTypeName, i), title + "#LT #Delta x' #GT", aliSize, 0.5, aliSize + 0.5);
1560  sumContainer.summaryNormXResiduals_ = tfd.make<TH1F>(
1561  Form("h_summaryNormX%s_%d", aliTypeName, i), title + "#LT #Delta x'/#sigma #GT", aliSize, 0.5, aliSize + 0.5);
1562 
1563  if (bookResidY) {
1564  sumContainer.summaryYResiduals_ = tfd.make<TH1F>(
1565  Form("h_summaryY%s_%d", aliTypeName, i), title + "#LT #Delta y' #GT", aliSize, 0.5, aliSize + 0.5);
1566  sumContainer.summaryNormYResiduals_ = tfd.make<TH1F>(
1567  Form("h_summaryNormY%s_%d", aliTypeName, i), title + "#LT #Delta y'/#sigma #GT", aliSize, 0.5, aliSize + 0.5);
1568  }
1569 
1570  } else if (subtype == align::AlignableDet && subcompSize > 1) { // fixed: was aliSize before
1571  if (subcompSize != 2) { // strange... expect only 2 DetUnits in DS layers
1572  // this 2 is hardcoded factor 2 in binning below and also assumed later on
1573  edm::LogError("Alignment") << "@SUB=bookSummaryHists"
1574  << "Det with " << subcompSize << " components";
1575  }
1576  // title contains x-title
1577  const TString title(Form(
1578  "Summary for substructures in %s %d;%s;",
1579  aliTypeName,
1580  i,
1581  alignableTracker_->objectIdProvider().idToString(ali.components()[0]->components()[0]->alignableObjectId())));
1582 
1583  sumContainer.summaryXResiduals_ = tfd.make<TH1F>(
1584  Form("h_summaryX%s_%d", aliTypeName, i), title + "#LT #Delta x' #GT", (2 * aliSize), 0.5, 2 * aliSize + 0.5);
1585  sumContainer.summaryNormXResiduals_ = tfd.make<TH1F>(Form("h_summaryNormX%s_%d", aliTypeName, i),
1586  title + "#LT #Delta x'/#sigma #GT",
1587  (2 * aliSize),
1588  0.5,
1589  2 * aliSize + 0.5);
1590 
1591  if (bookResidY) {
1592  sumContainer.summaryYResiduals_ = tfd.make<TH1F>(
1593  Form("h_summaryY%s_%d", aliTypeName, i), title + "#LT #Delta y' #GT", (2 * aliSize), 0.5, 2 * aliSize + 0.5);
1594  sumContainer.summaryNormYResiduals_ = tfd.make<TH1F>(Form("h_summaryNormY%s_%d", aliTypeName, i),
1595  title + "#LT #Delta y'/#sigma #GT",
1596  (2 * aliSize),
1597  0.5,
1598  2 * aliSize + 0.5);
1599  }
1600 
1601  } else {
1602  edm::LogError("TrackerOfflineValidation")
1603  << "@SUB=TrackerOfflineValidation::bookSummaryHists"
1604  << "No summary histogram for hierarchy level " << aliTypeName << " in subdet " << aliDetId.subdetId();
1605  }
1606 
1607  // Now book hists that just sum up the residual histograms from lower levels.
1608  // Axis title is copied from lowest level module of structure.
1609  // Should be safe that y-hists are only touched if non-null pointers...
1610  int nbins = 0;
1611  double xmin = 0., xmax = 0.;
1612  const TString sumTitle(Form("Residual for %s %d in %s;", aliTypeName, i, typeName));
1613  const ModuleHistos& xTitHists = this->getHistStructFromMap(aliDetId); // for x-axis titles
1614  this->getBinning(aliDetId.subdetId(), XprimeResidual, nbins, xmin, xmax);
1615 
1616  sumContainer.sumXResiduals_ = tfd.make<TH1F>(Form("h_Xprime_%s_%d", aliTypeName, i),
1617  sumTitle + xTitHists.ResXprimeHisto->GetXaxis()->GetTitle(),
1618  nbins,
1619  xmin,
1620  xmax);
1621 
1622  this->getBinning(aliDetId.subdetId(), NormXprimeResidual, nbins, xmin, xmax);
1623  sumContainer.sumNormXResiduals_ = tfd.make<TH1F>(Form("h_NormXprime_%s_%d", aliTypeName, i),
1624  sumTitle + xTitHists.NormResXprimeHisto->GetXaxis()->GetTitle(),
1625  nbins,
1626  xmin,
1627  xmax);
1628 
1629  if (moduleLevelProfiles_) {
1630  this->getBinning(aliDetId.subdetId(), XResidualProfile, nbins, xmin, xmax);
1631  sumContainer.sumResXvsXProfile_ = tfd.make<TProfile>(Form("p_resXX_%s_%d", aliTypeName, i),
1632  sumTitle + xTitHists.ResXvsXProfile->GetXaxis()->GetTitle() +
1633  ";" + xTitHists.ResXvsXProfile->GetYaxis()->GetTitle(),
1634  nbins,
1635  xmin,
1636  xmax);
1637  sumContainer.sumResXvsXProfile_->Sumw2();
1638  sumContainer.sumResXvsYProfile_ = tfd.make<TProfile>(Form("p_resXY_%s_%d", aliTypeName, i),
1639  sumTitle + xTitHists.ResXvsYProfile->GetXaxis()->GetTitle() +
1640  ";" + xTitHists.ResXvsYProfile->GetYaxis()->GetTitle(),
1641  nbins,
1642  xmin,
1643  xmax);
1644  sumContainer.sumResXvsYProfile_->Sumw2();
1645  }
1646 
1647  if (bookResidY) {
1648  this->getBinning(aliDetId.subdetId(), YprimeResidual, nbins, xmin, xmax);
1649  sumContainer.sumYResiduals_ = tfd.make<TH1F>(Form("h_Yprime_%s_%d", aliTypeName, i),
1650  sumTitle + xTitHists.ResYprimeHisto->GetXaxis()->GetTitle(),
1651  nbins,
1652  xmin,
1653  xmax);
1654 
1655  this->getBinning(aliDetId.subdetId(), NormYprimeResidual, nbins, xmin, xmax);
1656  sumContainer.sumNormYResiduals_ = tfd.make<TH1F>(Form("h_NormYprime_%s_%d", aliTypeName, i),
1657  sumTitle + xTitHists.NormResYprimeHisto->GetXaxis()->GetTitle(),
1658  nbins,
1659  xmin,
1660  xmax);
1661 
1662  if (moduleLevelProfiles_) {
1663  this->getBinning(aliDetId.subdetId(), YResidualProfile, nbins, xmin, xmax);
1664  sumContainer.sumResYvsXProfile_ = tfd.make<TProfile>(Form("p_resYX_%s_%d", aliTypeName, i),
1665  sumTitle + xTitHists.ResYvsXProfile->GetXaxis()->GetTitle() +
1666  ";" + xTitHists.ResYvsXProfile->GetYaxis()->GetTitle(),
1667  nbins,
1668  xmin,
1669  xmax);
1670  sumContainer.sumResYvsXProfile_->Sumw2();
1671  sumContainer.sumResYvsYProfile_ = tfd.make<TProfile>(Form("p_resYY_%s_%d", aliTypeName, i),
1672  sumTitle + xTitHists.ResYvsYProfile->GetXaxis()->GetTitle() +
1673  ";" + xTitHists.ResYvsYProfile->GetYaxis()->GetTitle(),
1674  nbins,
1675  xmin,
1676  xmax);
1677  sumContainer.sumResYvsYProfile_->Sumw2();
1678  }
1679  }
1680 
1681  // If we are at the lowest level, we already sum up and fill the summary.
1682 
1683  // special case I: For DetUnits and Dets with only one subcomponent start filling summary histos
1684  if ((subtype == align::AlignableDet && subcompSize == 1) || subtype == align::AlignableDetUnit) {
1685  for (uint k = 0; k < aliSize; ++k) {
1686  DetId detid = ali.components()[k]->id();
1687  ModuleHistos& histStruct = this->getHistStructFromMap(detid);
1688  this->summarizeBinInContainer(k + 1, detid.subdetId(), sumContainer, histStruct);
1689  sumHistStructure_.emplace_back(sumContainer.sumXResiduals_, histStruct.ResXprimeHisto);
1690  sumHistStructure_.emplace_back(sumContainer.sumNormXResiduals_, histStruct.NormResXprimeHisto);
1691  if (moduleLevelProfiles_) {
1692  sumHistStructure_.emplace_back(sumContainer.sumResXvsXProfile_, histStruct.ResXvsXProfile);
1693  sumHistStructure_.emplace_back(sumContainer.sumResXvsYProfile_, histStruct.ResXvsYProfile);
1694  }
1695  if (this->isPixel(detid.subdetId()) || stripYResiduals_) {
1696  sumHistStructure_.emplace_back(sumContainer.sumYResiduals_, histStruct.ResYprimeHisto);
1697  sumHistStructure_.emplace_back(sumContainer.sumNormYResiduals_, histStruct.NormResYprimeHisto);
1698  if (moduleLevelProfiles_) {
1699  sumHistStructure_.emplace_back(sumContainer.sumResYvsXProfile_, histStruct.ResYvsXProfile);
1700  sumHistStructure_.emplace_back(sumContainer.sumResYvsYProfile_, histStruct.ResYvsYProfile);
1701  }
1702  }
1703  }
1704  } else if (subtype == align::AlignableDet && subcompSize > 1) { // fixed: was aliSize before
1705  // special case II: Fill summary histos for dets with two detunits
1706  for (uint k = 0; k < aliSize; ++k) {
1707  for (uint j = 0; j < subcompSize; ++j) { // assumes all have same size (as binning does)
1708  DetId detid = ali.components()[k]->components()[j]->id();
1709  ModuleHistos& histStruct = this->getHistStructFromMap(detid);
1710  this->summarizeBinInContainer(2 * k + j + 1, detid.subdetId(), sumContainer, histStruct);
1711  sumHistStructure_.emplace_back(sumContainer.sumXResiduals_, histStruct.ResXprimeHisto);
1712  sumHistStructure_.emplace_back(sumContainer.sumNormXResiduals_, histStruct.NormResXprimeHisto);
1713  if (moduleLevelProfiles_) {
1714  sumHistStructure_.emplace_back(sumContainer.sumResXvsXProfile_, histStruct.ResXvsXProfile);
1715  sumHistStructure_.emplace_back(sumContainer.sumResXvsYProfile_, histStruct.ResXvsYProfile);
1716  }
1717  if (this->isPixel(detid.subdetId()) || stripYResiduals_) {
1718  sumHistStructure_.emplace_back(sumContainer.sumYResiduals_, histStruct.ResYprimeHisto);
1719  sumHistStructure_.emplace_back(sumContainer.sumNormYResiduals_, histStruct.NormResYprimeHisto);
1720  if (moduleLevelProfiles_) {
1721  sumHistStructure_.emplace_back(sumContainer.sumResYvsXProfile_, histStruct.ResYvsXProfile);
1722  sumHistStructure_.emplace_back(sumContainer.sumResYvsYProfile_, histStruct.ResYvsYProfile);
1723  }
1724  }
1725  }
1726  }
1727  }
1728  return sumContainer;
1729 }
1730 
1731 float TrackerOfflineValidation::Fwhm(const TH1* hist) const {
1732  float max = hist->GetMaximum();
1733  int left = -1, right = -1;
1734  for (unsigned int i = 1, iEnd = hist->GetNbinsX(); i <= iEnd; ++i) {
1735  if (hist->GetBinContent(i) < max / 2. && hist->GetBinContent(i + 1) > max / 2. && left == -1) {
1736  if (max / 2. - hist->GetBinContent(i) < hist->GetBinContent(i + 1) - max / 2.) {
1737  left = i;
1738  ++i;
1739  } else {
1740  left = i + 1;
1741  ++i;
1742  }
1743  }
1744  if (left != -1 && right == -1) {
1745  if (hist->GetBinContent(i) > max / 2. && hist->GetBinContent(i + 1) < max / 2.) {
1746  if (hist->GetBinContent(i) - max / 2. < max / 2. - hist->GetBinContent(i + 1)) {
1747  right = i;
1748  } else {
1749  right = i + 1;
1750  }
1751  }
1752  }
1753  }
1754  return hist->GetXaxis()->GetBinCenter(right) - hist->GetXaxis()->GetBinCenter(left);
1755 }
1756 
1757 void TrackerOfflineValidation::setUpTreeMembers(const std::map<int, TrackerOfflineValidation::ModuleHistos>& moduleHist_,
1758  const TrackerGeometry& tkgeom,
1759  const TrackerTopology* tTopo) {
1760  for (std::map<int, TrackerOfflineValidation::ModuleHistos>::const_iterator it = moduleHist_.begin(),
1761  itEnd = moduleHist_.end();
1762  it != itEnd;
1763  ++it) {
1764  TkOffTreeVariables& treeMem = mTreeMembers_[it->first];
1765 
1766  //variables concerning the tracker components/hierarchy levels
1767  DetId detId_ = it->first;
1768  treeMem.moduleId = detId_;
1769  treeMem.subDetId = detId_.subdetId();
1770  treeMem.isDoubleSide = false;
1771 
1772  if (treeMem.subDetId == PixelSubdetector::PixelBarrel) {
1773  unsigned int whichHalfBarrel(1), rawId(detId_.rawId()); //DetId does not know about halfBarrels is PXB ...
1774  if ((rawId >= 302056964 && rawId < 302059300) || (rawId >= 302123268 && rawId < 302127140) ||
1775  (rawId >= 302189572 && rawId < 302194980))
1776  whichHalfBarrel = 2;
1777  treeMem.layer = tTopo->pxbLayer(detId_);
1778  treeMem.half = whichHalfBarrel;
1779  treeMem.rod = tTopo->pxbLadder(detId_); // ... so, ladder is not per halfBarrel-Layer, but per barrel-layer!
1780  treeMem.module = tTopo->pxbModule(detId_);
1781  } else if (treeMem.subDetId == PixelSubdetector::PixelEndcap) {
1782  unsigned int whichHalfCylinder(1), rawId(detId_.rawId()); //DetId does not kmow about halfCylinders in PXF
1783  if ((rawId >= 352394500 && rawId < 352406032) || (rawId >= 352460036 && rawId < 352471568) ||
1784  (rawId >= 344005892 && rawId < 344017424) || (rawId >= 344071428 && rawId < 344082960))
1785  whichHalfCylinder = 2;
1786  treeMem.layer = tTopo->pxfDisk(detId_);
1787  treeMem.side = tTopo->pxfSide(detId_);
1788  treeMem.half = whichHalfCylinder;
1789  treeMem.blade = tTopo->pxfBlade(detId_);
1790  treeMem.panel = tTopo->pxfPanel(detId_);
1791  treeMem.module = tTopo->pxfModule(detId_);
1792  } else if (treeMem.subDetId == StripSubdetector::TIB) {
1793  unsigned int whichHalfShell(1), rawId(detId_.rawId()); //DetId does not kmow about halfShells in TIB
1794  if ((rawId >= 369120484 && rawId < 369120688) || (rawId >= 369121540 && rawId < 369121776) ||
1795  (rawId >= 369136932 && rawId < 369137200) || (rawId >= 369137988 && rawId < 369138288) ||
1796  (rawId >= 369153396 && rawId < 369153744) || (rawId >= 369154436 && rawId < 369154800) ||
1797  (rawId >= 369169844 && rawId < 369170256) || (rawId >= 369170900 && rawId < 369171344) ||
1798  (rawId >= 369124580 && rawId < 369124784) || (rawId >= 369125636 && rawId < 369125872) ||
1799  (rawId >= 369141028 && rawId < 369141296) || (rawId >= 369142084 && rawId < 369142384) ||
1800  (rawId >= 369157492 && rawId < 369157840) || (rawId >= 369158532 && rawId < 369158896) ||
1801  (rawId >= 369173940 && rawId < 369174352) || (rawId >= 369174996 && rawId < 369175440))
1802  whichHalfShell = 2;
1803  treeMem.layer = tTopo->tibLayer(detId_);
1804  treeMem.side = tTopo->tibStringInfo(detId_)[0];
1805  treeMem.half = whichHalfShell;
1806  treeMem.rod = tTopo->tibStringInfo(detId_)[2];
1807  treeMem.outerInner = tTopo->tibStringInfo(detId_)[1];
1808  treeMem.module = tTopo->tibModule(detId_);
1809  treeMem.isStereo = tTopo->tibStereo(detId_);
1810  treeMem.isDoubleSide = tTopo->tibIsDoubleSide(detId_);
1811  } else if (treeMem.subDetId == StripSubdetector::TID) {
1812  treeMem.layer = tTopo->tidWheel(detId_);
1813  treeMem.side = tTopo->tidSide(detId_);
1814  treeMem.ring = tTopo->tidRing(detId_);
1815  treeMem.outerInner = tTopo->tidModuleInfo(detId_)[0];
1816  treeMem.module = tTopo->tidModuleInfo(detId_)[1];
1817  treeMem.isStereo = tTopo->tidStereo(detId_);
1818  treeMem.isDoubleSide = tTopo->tidIsDoubleSide(detId_);
1819  } else if (treeMem.subDetId == StripSubdetector::TOB) {
1820  treeMem.layer = tTopo->tobLayer(detId_);
1821  treeMem.side = tTopo->tobRodInfo(detId_)[0];
1822  treeMem.rod = tTopo->tobRodInfo(detId_)[1];
1823  treeMem.module = tTopo->tobModule(detId_);
1824  treeMem.isStereo = tTopo->tobStereo(detId_);
1825  treeMem.isDoubleSide = tTopo->tobIsDoubleSide(detId_);
1826  } else if (treeMem.subDetId == StripSubdetector::TEC) {
1827  treeMem.layer = tTopo->tecWheel(detId_);
1828  treeMem.side = tTopo->tecSide(detId_);
1829  treeMem.ring = tTopo->tecRing(detId_);
1830  treeMem.petal = tTopo->tecPetalInfo(detId_)[1];
1831  treeMem.outerInner = tTopo->tecPetalInfo(detId_)[0];
1832  treeMem.module = tTopo->tecModule(detId_);
1833  treeMem.isStereo = tTopo->tecStereo(detId_);
1834  treeMem.isDoubleSide = tTopo->tecIsDoubleSide(detId_);
1835  }
1836 
1837  //variables concerning the tracker geometry
1838  const Surface::PositionType& gPModule = tkgeom.idToDet(detId_)->position();
1839  treeMem.posPhi = gPModule.phi();
1840  treeMem.posEta = gPModule.eta();
1841  treeMem.posR = gPModule.perp();
1842  treeMem.posX = gPModule.x();
1843  treeMem.posY = gPModule.y();
1844  treeMem.posZ = gPModule.z();
1845 
1846  const Surface& surface = tkgeom.idToDet(detId_)->surface();
1847 
1848  //global Orientation of local coordinate system of dets/detUnits
1849  LocalPoint lUDirection(1., 0., 0.), lVDirection(0., 1., 0.), lWDirection(0., 0., 1.);
1850  GlobalPoint gUDirection = surface.toGlobal(lUDirection), gVDirection = surface.toGlobal(lVDirection),
1851  gWDirection = surface.toGlobal(lWDirection);
1852  double dR(999.), dPhi(999.), dZ(999.);
1854  treeMem.subDetId == StripSubdetector::TOB) {
1855  dR = gWDirection.perp() - gPModule.perp();
1856  dPhi = deltaPhi(gUDirection.barePhi(), gPModule.barePhi());
1857  dZ = gVDirection.z() - gPModule.z();
1858  if (dZ >= 0.)
1859  treeMem.rOrZDirection = 1;
1860  else
1861  treeMem.rOrZDirection = -1;
1862  } else if (treeMem.subDetId == PixelSubdetector::PixelEndcap) {
1863  dR = gUDirection.perp() - gPModule.perp();
1864  dPhi = deltaPhi(gVDirection.barePhi(), gPModule.barePhi());
1865  dZ = gWDirection.z() - gPModule.z();
1866  if (dR >= 0.)
1867  treeMem.rOrZDirection = 1;
1868  else
1869  treeMem.rOrZDirection = -1;
1870  } else if (treeMem.subDetId == StripSubdetector::TID || treeMem.subDetId == StripSubdetector::TEC) {
1871  dR = gVDirection.perp() - gPModule.perp();
1872  dPhi = deltaPhi(gUDirection.barePhi(), gPModule.barePhi());
1873  dZ = gWDirection.z() - gPModule.z();
1874  if (dR >= 0.)
1875  treeMem.rOrZDirection = 1;
1876  else
1877  treeMem.rOrZDirection = -1;
1878  }
1879  if (dR >= 0.)
1880  treeMem.rDirection = 1;
1881  else
1882  treeMem.rDirection = -1;
1883  if (dPhi >= 0.)
1884  treeMem.phiDirection = 1;
1885  else
1886  treeMem.phiDirection = -1;
1887  if (dZ >= 0.)
1888  treeMem.zDirection = 1;
1889  else
1890  treeMem.zDirection = -1;
1891  }
1892 }
1893 
1895  TkOffTreeVariables& treeMem,
1896  const std::map<int, TrackerOfflineValidation::ModuleHistos>& moduleHist_) {
1897  for (std::map<int, TrackerOfflineValidation::ModuleHistos>::const_iterator it = moduleHist_.begin(),
1898  itEnd = moduleHist_.end();
1899  it != itEnd;
1900  ++it) {
1901  treeMem = mTreeMembers_[it->first];
1902 
1903  //mean and RMS values (extracted from histograms(Xprime on module level)
1904  treeMem.entries = static_cast<UInt_t>(it->second.ResXprimeHisto->GetEntries());
1905  treeMem.meanX = it->second.ResXprimeHisto->GetMean();
1906  treeMem.rmsX = it->second.ResXprimeHisto->GetRMS();
1907  //treeMem.sigmaX = Fwhm(it->second.ResXprimeHisto)/2.355;
1908 
1909  if (useFit_) {
1910  //call fit function which returns mean and sigma from the fit
1911  //for absolute residuals
1912  std::pair<float, float> fitResult1 = this->fitResiduals(it->second.ResXprimeHisto);
1913  treeMem.fitMeanX = fitResult1.first;
1914  treeMem.fitSigmaX = fitResult1.second;
1915  //for normalized residuals
1916  std::pair<float, float> fitResult2 = this->fitResiduals(it->second.NormResXprimeHisto);
1917  treeMem.fitMeanNormX = fitResult2.first;
1918  treeMem.fitSigmaNormX = fitResult2.second;
1919  }
1920 
1921  //get median for absolute residuals
1922  treeMem.medianX = this->getMedian(it->second.ResXprimeHisto);
1923 
1924  int numberOfBins = it->second.ResXprimeHisto->GetNbinsX();
1925  treeMem.numberOfUnderflows = it->second.ResXprimeHisto->GetBinContent(0);
1926  treeMem.numberOfOverflows = it->second.ResXprimeHisto->GetBinContent(numberOfBins + 1);
1927  treeMem.numberOfOutliers =
1928  it->second.ResXprimeHisto->GetBinContent(0) + it->second.ResXprimeHisto->GetBinContent(numberOfBins + 1);
1929 
1930  //mean and RMS values (extracted from histograms(normalized Xprime on module level)
1931  treeMem.meanNormX = it->second.NormResXprimeHisto->GetMean();
1932  treeMem.rmsNormX = it->second.NormResXprimeHisto->GetRMS();
1933 
1934  double stats[20];
1935  it->second.NormResXprimeHisto->GetStats(stats);
1936  // GF treeMem.chi2PerDofX = stats[3]/(stats[0]-1);
1937  if (stats[0])
1938  treeMem.chi2PerDofX = stats[3] / stats[0];
1939 
1940  treeMem.sigmaNormX = Fwhm(it->second.NormResXprimeHisto) / 2.355;
1941  treeMem.histNameX = it->second.ResXprimeHisto->GetName();
1942  treeMem.histNameNormX = it->second.NormResXprimeHisto->GetName();
1943 
1944  // fill tree variables in local coordinates if set in cfg
1945  if (lCoorHistOn_) {
1946  treeMem.meanLocalX = it->second.ResHisto->GetMean();
1947  treeMem.rmsLocalX = it->second.ResHisto->GetRMS();
1948  treeMem.meanNormLocalX = it->second.NormResHisto->GetMean();
1949  treeMem.rmsNormLocalX = it->second.NormResHisto->GetRMS();
1950 
1951  treeMem.histNameLocalX = it->second.ResHisto->GetName();
1952  treeMem.histNameNormLocalX = it->second.NormResHisto->GetName();
1953  if (it->second.ResYHisto)
1954  treeMem.histNameLocalY = it->second.ResYHisto->GetName();
1955  }
1956 
1957  // mean and RMS values in local y (extracted from histograms(normalized Yprime on module level)
1958  // might exist in pixel only
1959  if (it->second.ResYprimeHisto) { //(stripYResiduals_)
1960  TH1* h = it->second.ResYprimeHisto;
1961  treeMem.meanY = h->GetMean();
1962  treeMem.rmsY = h->GetRMS();
1963 
1964  if (useFit_) { // fit function which returns mean and sigma from the fit
1965  std::pair<float, float> fitMeanSigma = this->fitResiduals(h);
1966  treeMem.fitMeanY = fitMeanSigma.first;
1967  treeMem.fitSigmaY = fitMeanSigma.second;
1968  }
1969 
1970  //get median for absolute residuals
1971  treeMem.medianY = this->getMedian(h);
1972 
1973  treeMem.histNameY = h->GetName();
1974  }
1975  if (it->second.NormResYprimeHisto) {
1976  TH1* h = it->second.NormResYprimeHisto;
1977  treeMem.meanNormY = h->GetMean();
1978  treeMem.rmsNormY = h->GetRMS();
1979  h->GetStats(stats); // stats buffer defined above
1980  if (stats[0])
1981  treeMem.chi2PerDofY = stats[3] / stats[0];
1982 
1983  if (useFit_) { // fit function which returns mean and sigma from the fit
1984  std::pair<float, float> fitMeanSigma = this->fitResiduals(h);
1985  treeMem.fitMeanNormY = fitMeanSigma.first;
1986  treeMem.fitSigmaNormY = fitMeanSigma.second;
1987  }
1988  treeMem.histNameNormY = h->GetName();
1989  }
1990 
1991  if (moduleLevelProfiles_) {
1992  if (it->second.ResXvsXProfile) {
1993  TH1* h = it->second.ResXvsXProfile;
1994  treeMem.meanResXvsX = h->GetMean();
1995  treeMem.rmsResXvsX = h->GetRMS();
1996  treeMem.profileNameResXvsX = h->GetName();
1997  }
1998  if (it->second.ResXvsYProfile) {
1999  TH1* h = it->second.ResXvsYProfile;
2000  treeMem.meanResXvsY = h->GetMean();
2001  treeMem.rmsResXvsY = h->GetRMS();
2002  treeMem.profileNameResXvsY = h->GetName();
2003  }
2004  if (it->second.ResYvsXProfile) {
2005  TH1* h = it->second.ResYvsXProfile;
2006  treeMem.meanResYvsX = h->GetMean();
2007  treeMem.rmsResYvsX = h->GetRMS();
2008  treeMem.profileNameResYvsX = h->GetName();
2009  }
2010  if (it->second.ResYvsYProfile) {
2011  TH1* h = it->second.ResYvsYProfile;
2012  treeMem.meanResYvsY = h->GetMean();
2013  treeMem.rmsResYvsY = h->GetRMS();
2014  treeMem.profileNameResYvsY = h->GetName();
2015  }
2016  }
2017 
2018  tree.Fill();
2019  }
2020 }
2021 
2022 std::pair<float, float> TrackerOfflineValidation::fitResiduals(TH1* hist) const {
2023  std::pair<float, float> fitResult(9999., 9999.);
2024  if (!hist || hist->GetEntries() < 20)
2025  return fitResult;
2026 
2027  float mean = hist->GetMean();
2028  float sigma = hist->GetRMS();
2029 
2030  try { // for < CMSSW_2_2_0 since ROOT warnings from fit are converted to exceptions
2031  // Remove the try/catch for more recent CMSSW!
2032  // first fit: two RMS around mean
2033  TF1 func("tmp", "gaus", mean - 2. * sigma, mean + 2. * sigma);
2034  if (0 == hist->Fit(&func, "QNR")) { // N: do not blow up file by storing fit!
2035  mean = func.GetParameter(1);
2036  sigma = func.GetParameter(2);
2037  // second fit: three sigma of first fit around mean of first fit
2038  func.SetRange(mean - 3. * sigma, mean + 3. * sigma);
2039  // I: integral gives more correct results if binning is too wide
2040  // L: Likelihood can treat empty bins correctly (if hist not weighted...)
2041  if (0 == hist->Fit(&func, "Q0LR")) {
2042  if (hist->GetFunction(func.GetName())) { // Take care that it is later on drawn:
2043  hist->GetFunction(func.GetName())->ResetBit(TF1::kNotDraw);
2044  }
2045  fitResult.first = func.GetParameter(1);
2046  fitResult.second = func.GetParameter(2);
2047  }
2048  }
2049  } catch (cms::Exception const& e) {
2050  edm::LogWarning("Alignment") << "@SUB=TrackerOfflineValidation::fitResiduals"
2051  << "Caught this exception during ROOT fit: " << e.what();
2052  }
2053  return fitResult;
2054 }
2055 
2057  float median = 999;
2058  int nbins = histo->GetNbinsX();
2059 
2060  //extract median from histogram
2061  double* x = new double[nbins];
2062  double* y = new double[nbins];
2063  for (int j = 0; j < nbins; j++) {
2064  x[j] = histo->GetBinCenter(j + 1);
2065  y[j] = histo->GetBinContent(j + 1);
2066  }
2067  median = TMath::Median(nbins, x, y);
2068 
2069  delete[] x;
2070  x = nullptr;
2071  delete[] y;
2072  y = nullptr;
2073 
2074  return median;
2075 }
2076 //define this as a plug-in
const TrackerGeometry * bareTkGeomPtr_
size
Write out results.
static const std::string kSharedResource
Definition: TFileService.h:76
static constexpr auto TEC
bool tibIsDoubleSide(const DetId &id) const
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
std::map< int, TrackerOfflineValidation::ModuleHistos > mTidResiduals_
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
unsigned int tobLayer(const DetId &id) const
unsigned int pxbLayer(const DetId &id) const
std::vector< TH1 * > vDeleteObjects_
T perp() const
Definition: PV3DBase.h:69
moduleLevelHistsTransient
Definition: DMR_cfg.py:152
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
Definition: EventSetup.h:119
const unsigned long long maxTracks_
std::vector< unsigned int > tidModuleInfo(const DetId &id) const
TH1 * make(const char *name, const char *title, int nBinX, double minBinX, double maxBinX)
unsigned int pxfBlade(const DetId &id) const
T z() const
Definition: PV3DBase.h:61
const DetIdContainer & detIds() const override
Returm a vector of all GeomDet DetIds (including those of GeomDetUnits)
std::vector< TH1 * > vTrack2DHistos_
TH1 * bookTH1F(bool isTransient, DirectoryWrapper &tfd, const char *histName, const char *histTitle, int nBinsX, double lowX, double highX)
std::unique_ptr< AlignableTracker > alignableTracker_
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
std::string profileNameResXvsX
void fillTree(TTree &tree, TkOffTreeVariables &treeMem, const std::map< int, TrackerOfflineValidation::ModuleHistos > &moduleHist_)
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
constexpr unsigned int subDetId[21]
T eta() const
Definition: PV3DBase.h:73
edm::ESHandle< TrackerGeometry > tkGeom_
unsigned int pxbLadder(const DetId &id) const
Log< level::Error, false > LogError
ModuleHistos & getHistStructFromMap(const DetId &detid)
const edm::ESGetToken< TrackerGeometry, TrackerDigiGeometryRecord > geomToken_
void setSummaryBin(int bin, TH1 *targetHist, TH1 *sourceHist)
TrackerOfflineValidation::SummaryContainer bookSummaryHists(DirectoryWrapper &tfd, const Alignable &ali, align::StructureType type, int i)
unsigned int tecRing(const DetId &id) const
ring id
bool tobIsDoubleSide(const DetId &id) const
const edm::ESGetToken< TrackerTopology, TrackerTopologyRcd > topoToken_
std::map< int, TrackerOfflineValidation::ModuleHistos > mTobResiduals_
T barePhi() const
Definition: PV3DBase.h:65
numberOfBins
Definition: PV_cfg.py:227
std::string profileNameResYvsX
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
T const * product() const
Definition: ESHandle.h:86
std::map< int, TrackerOfflineValidation::ModuleHistos > mPxeResiduals_
uint32_t tobStereo(const DetId &id) const
std::map< int, TkOffTreeVariables > mTreeMembers_
void setUpTreeMembers(const std::map< int, TrackerOfflineValidation::ModuleHistos > &moduleHist_, const TrackerGeometry &tkgeom, const TrackerTopology *tTopo)
unsigned int tecSide(const DetId &id) const
float Fwhm(const TH1 *hist) const
float getMedian(const TH1 *hist) const
std::vector< std::pair< TH1 *, TH1 * > > sumHistStructure_
unsigned int pxfDisk(const DetId &id) const
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
virtual StructureType alignableObjectId() const =0
Return the alignable type identifier.
void getBinning(uint32_t subDetId, TrackerOfflineValidation::HistogramType residualtype, int &nBinsX, double &lowerBoundX, double &upperBoundX)
virtual void checkBookHists(const edm::EventSetup &setup)
TProfile * bookTProfile(bool isTransient, DirectoryWrapper &tfd, const char *histName, const char *histTitle, int nBinsX, double lowX, double highX)
double f[11][100]
static void fillPSetDescription(edm::ParameterSetDescription &descriptions)
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
const unsigned int maxEntriesPerModuleForDmr_
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
void prepareSummaryHists(DirectoryWrapper &tfd, const Alignable &ali, std::vector< TrackerOfflineValidation::SummaryContainer > &vLevelProfiles)
virtual const Alignables & components() const =0
Return vector of all direct components.
ParameterDescriptionBase * add(U const &iLabel, T const &value)
std::string histNameLocalY
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:130
const TrackerGeomDet * idToDet(DetId) const override
std::string histNameNormLocalX
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
bool isDetOrDetUnit(align::StructureType type)
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:79
unsigned int pxfPanel(const DetId &id) const
Log< level::Info, false > LogInfo
DirectoryWrapper(const std::string &newDir, const std::string &basedir, bool useDqmMode)
Definition: DetId.h:17
std::map< int, TrackerOfflineValidation::ModuleHistos > mTecResiduals_
align::ID id() const
Return the ID of Alignable, i.e. DetId of &#39;first&#39; component GeomDet(Unit).
Definition: Alignable.h:180
unsigned int pxfSide(const DetId &id) const
std::vector< std::tuple< int, TH1 *, TH1 * > > summaryBins_
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 fillTrackQuantities(const edm::Event &, const edm::EventSetup &, std::vector< AVTrackStruct > &v_avtrackout)
std::map< int, TrackerOfflineValidation::ModuleHistos > mTibResiduals_
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
int GetIndex(const std::vector< OBJECT_TYPE *> &vec, const TString &name)
std::string histNameLocalX
void summarizeBinInContainer(int bin, SummaryContainer &targetContainer, SummaryContainer &sourceContainer)
bool isEndCap(uint32_t subDetId)
const Surface::PositionType & position() const
The position (origin of the R.F.)
Definition: GeomDet.h:43
std::string profileNameResXvsY
T median(std::vector< T > values)
Definition: median.h:16
std::vector< unsigned int > tobRodInfo(const DetId &id) const
void bookGlobalHists(DirectoryWrapper &tfd)
DirectoryWrapper(const DirectoryWrapper &upDir, const std::string &newDir, const std::string &basedir, bool useDqmMode)
TrackerValidationVariables avalidator_
std::pair< int, int > typeAndLayerFromDetId(const DetId &detId, const TrackerTopology *tTopo) const
std::vector< TH1 * > vTrackProfiles_
unsigned int tidRing(const DetId &id) const
void bookDirHists(DirectoryWrapper &tfd, const Alignable &ali, const TrackerTopology *tTopo)
uint32_t tecStereo(const DetId &id) const
unsigned int tibLayer(const DetId &id) const
unsigned int tobModule(const DetId &id) const
Definition: tree.py:1
unsigned int pxbModule(const DetId &id) const
Log< level::Warning, false > LogWarning
std::pair< float, float > fitResiduals(TH1 *hist) const
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
const edm::ParameterSet parSet_
TrackerOfflineValidation(const edm::ParameterSet &)
static constexpr auto TID
const DetIdContainer & detUnitIds() const override
Returm a vector of all GeomDetUnit DetIds.
std::map< int, TrackerOfflineValidation::ModuleHistos > mPxbResiduals_
uint32_t tibStereo(const DetId &id) const
void analyze(const edm::Event &, const edm::EventSetup &) override
std::string profileNameResYvsY
void bookHists(DirectoryWrapper &tfd, const Alignable &ali, const TrackerTopology *tTopo, align::StructureType type, int i)
#define LogDebug(id)
bool isBarrel(uint32_t subDetId)