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