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