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