CMS 3D CMS Logo

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