CMS 3D CMS Logo

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