CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TrackerOfflineValidationSummary.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: TrackerOfflineValidationSummary
4 // Class: TrackerOfflineValidationSummary
5 //
13 //
14 // Original Author: Johannes Hauk
15 // Created: Sat Aug 22 10:31:34 CEST 2009
16 // $Id: TrackerOfflineValidationSummary.cc,v 1.7 2012/06/30 09:09:47 eulisse Exp $
17 //
18 //
19 
20 
21 // system include files
22 #include <memory>
23 
24 // user include files
27 
30 
32 
33 
39 
40 
41 #include "TTree.h"
42 //#include "Alignment/OfflineValidation/interface/TrackerValidationVariables.h"
44 
52 #include "TH1.h"
53 #include "TMath.h"
54 
58 
59 //
60 // class decleration
61 //
62 
64  public:
67 
68 
69  private:
70 
71  struct ModuleHistos{
74  TH1* ResHisto;
80  };
81 
85  TH1* DmrXprime;
86  TH1* DmrYprime;
87  };
88 
91  HarvestingHierarchy(const std::string& name, const std::string& component, const std::vector<unsigned int>& entries) : hierarchyName(name), componentName(component), treeEntries(entries) {}
92  // Needed for naming of histogram according to residual histograms
94  // "Pixel" or "Strip"
96  // Modules belonging to selected hierarchy
97  std::vector<unsigned int> treeEntries;
99  };
100 
101  virtual void analyze(const edm::Event& evt, const edm::EventSetup&) override;
102  virtual void endJob() override ;
103 
104  void fillTree(TTree& tree, std::map<int, TrackerOfflineValidationSummary::ModuleHistos>& moduleHist,
105  TkOffTreeVariables& treeMem, const TrackerGeometry& tkgeom,
106  std::map<std::string,std::string>& substructureName, const TrackerTopology* tTopo);
107 
108  std::pair<float,float> fitResiduals(TH1* hist)const;
109  float getMedian(const TH1* hist)const;
110 
111  const std::string associateModuleHistsWithTree(const TkOffTreeVariables& treeMem, TrackerOfflineValidationSummary::ModuleHistos& moduleHists, std::map<std::string,std::string>& substructureName);
112 
113  void collateHarvestingHists(TTree& tree);
114  void applyHarvestingHierarchy(TTree& treeMem);
115  void bookHarvestingHists();
116  void getBinning(const std::string& binningPSetName, int& nBinsX, double& lowerBoundX, double& upperBoundX)const;
117  void fillHarvestingHists(TTree& tree);
118 
119  // ----------member data ---------------------------
120 
123 
124  // parameters from cfg to steer
126  const bool useFit_;
127 
129 
131 
132  std::map<int,TrackerOfflineValidationSummary::ModuleHistos> mPxbResiduals_;
133  std::map<int,TrackerOfflineValidationSummary::ModuleHistos> mPxeResiduals_;
134  std::map<int,TrackerOfflineValidationSummary::ModuleHistos> mTibResiduals_;
135  std::map<int,TrackerOfflineValidationSummary::ModuleHistos> mTidResiduals_;
136  std::map<int,TrackerOfflineValidationSummary::ModuleHistos> mTobResiduals_;
137  std::map<int,TrackerOfflineValidationSummary::ModuleHistos> mTecResiduals_;
138 
139  std::vector<HarvestingHierarchy> vHarvestingHierarchy_;
140 
142 };
143 
144 //
145 // constants, enums and typedefs
146 //
147 
148 //
149 // static data member definitions
150 //
151 
152 //
153 // constructors and destructor
154 //
156  parSet_(iConfig), moduleDirectory_(parSet_.getParameter<std::string>("moduleDirectoryInOutput")),
157  useFit_(parSet_.getParameter<bool>("useFit")), dbe_(0), moduleMapsInitialized(false), lastSetup_(nullptr)
158 {
159  //now do what ever initialization is needed
161 }
162 
163 
165 {
166 
167  // do anything here that needs to be done at desctruction time
168  // (e.g. close files, deallocate resources etc.)
169 
170 }
171 
172 
173 //
174 // member functions
175 //
176 
177 // ------------ method called to for each event ------------
178 void
180 {
181  lastSetup_ = &iSetup;
182 
183  // Access of EventSetup is needed to get the list of silicon-modules and their IDs
184  // Since they do not change, it is accessed only once
185  if(moduleMapsInitialized)return;
186  iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom_ );
187  const TrackerGeometry* bareTkGeomPtr = &(*tkGeom_);
188  const TrackingGeometry::DetIdContainer& detIdContainer = bareTkGeomPtr->detIds();
189  std::vector<DetId>::const_iterator iDet;
190  for(iDet = detIdContainer.begin(); iDet != detIdContainer.end(); ++iDet){
191  const DetId& detId = *iDet;
192  const uint32_t rawId = detId.rawId();
193  const unsigned int subdetId = detId.subdetId();
194  if (subdetId == PixelSubdetector::PixelBarrel) mPxbResiduals_[rawId];
195  else if(subdetId == PixelSubdetector::PixelEndcap) mPxeResiduals_[rawId];
196  else if(subdetId == StripSubdetector::TIB) mTibResiduals_[rawId];
197  else if(subdetId == StripSubdetector::TID) mTidResiduals_[rawId];
198  else if(subdetId == StripSubdetector::TOB) mTobResiduals_[rawId];
199  else if(subdetId == StripSubdetector::TEC) mTecResiduals_[rawId];
200  else {
201  throw cms::Exception("Geometry Error")
202  << "[TrackerOfflineValidationSummary] Error, tried to get reference for non-tracker subdet " << subdetId
203  << " from detector " << detId.det();
204  }
205  }
206  moduleMapsInitialized = true;
207 }
208 
209 
210 // ------------ method called once each job just after ending the event loop ------------
211 void
213 {
214  //Retrieve tracker topology from geometry
215  edm::ESHandle<TrackerTopology> tTopoHandle;
216  lastSetup_->get<IdealGeometryRecord>().get(tTopoHandle);
217  const TrackerTopology* const tTopo = tTopoHandle.product();
218 
219  AlignableTracker aliTracker(&(*tkGeom_), tTopo);
220 
221  TTree* tree = new TTree("TkOffVal","TkOffVal");
222 
223  TkOffTreeVariables *treeMemPtr = new TkOffTreeVariables;
224  // We create branches for all members of 'TkOffTreeVariables' (even if not needed).
225  // This works because we have a dictionary for 'TkOffTreeVariables'
226  // (see src/classes_def.xml and src/classes.h):
227  tree->Branch("TkOffTreeVariables", &treeMemPtr); // address of pointer!
228  // second branch needed for assigning names and titles to harvesting histograms consistent to others
229  std::map<std::string,std::string> *substructureName = new std::map<std::string,std::string>;
230  tree->Branch("SubstructureName", &substructureName, 32000, 00); // SplitLevel must be set to zero
231 
232 
233  this->fillTree(*tree, mPxbResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo);
234  this->fillTree(*tree, mPxeResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo);
235  this->fillTree(*tree, mTibResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo);
236  this->fillTree(*tree, mTidResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo);
237  this->fillTree(*tree, mTobResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo);
238  this->fillTree(*tree, mTecResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo);
239 
240  //dbe_->showDirStructure();
241  //dbe_->save("dqmOut.root");
242 
243  // Method for filling histograms which show summarized values (mean, rms, median ...)
244  // of the module-based histograms from TrackerOfflineValidation
245  this->collateHarvestingHists(*tree);
246 
247  delete tree; tree = 0;
248  delete treeMemPtr; treeMemPtr = 0;
249  delete substructureName; substructureName = 0;
250 }
251 
252 
253 void
254 TrackerOfflineValidationSummary::fillTree(TTree& tree, std::map<int, TrackerOfflineValidationSummary::ModuleHistos>& moduleHist,
255  TkOffTreeVariables& treeMem, const TrackerGeometry& tkgeom,
256  std::map<std::string,std::string>& substructureName,
257  const TrackerTopology* tTopo)
258 {
259  for(std::map<int, TrackerOfflineValidationSummary::ModuleHistos>::iterator it = moduleHist.begin(),
260  itEnd= moduleHist.end(); it != itEnd;++it ) {
261  treeMem.clear(); // make empty/default
262 
263  //variables concerning the tracker components/hierarchy levels
264  const DetId detId = it->first;
265  treeMem.moduleId = detId;
266  treeMem.subDetId = detId.subdetId();
267 
269 
270  unsigned int whichHalfBarrel(1), rawId(detId.rawId()); //DetId does not know about halfBarrels is PXB ...
271  if( (rawId>=302056964 && rawId<302059300) || (rawId>=302123268 && rawId<302127140) || (rawId>=302189572 && rawId<302194980) )whichHalfBarrel=2;
272  treeMem.layer = tTopo->pxbLayer(detId);
273  treeMem.half = whichHalfBarrel;
274  treeMem.rod = tTopo->pxbLadder(detId); // ... so, ladder is not per halfBarrel-Layer, but per barrel-layer! Needs complicated calculation in associateModuleHistsWithTree()
275  treeMem.module = tTopo->pxbModule(detId);
276  } else if(treeMem.subDetId == PixelSubdetector::PixelEndcap){
277 
278  unsigned int whichHalfCylinder(1), rawId(detId.rawId()); //DetId does not kmow about halfCylinders in PXF
279  if( (rawId>=352394500 && rawId<352406032) || (rawId>=352460036 && rawId<352471568) || (rawId>=344005892 && rawId<344017424) || (rawId>=344071428 && rawId<344082960) )whichHalfCylinder=2;
280  treeMem.layer = tTopo->pxfDisk(detId);
281  treeMem.side = tTopo->pxfSide(detId);
282  treeMem.half = whichHalfCylinder;
283  treeMem.blade = tTopo->pxfBlade(detId);
284  treeMem.panel = tTopo->pxfPanel(detId);
285  treeMem.module = tTopo->pxfModule(detId);
286  } else if(treeMem.subDetId == StripSubdetector::TIB){
287 
288  unsigned int whichHalfShell(1), rawId(detId.rawId()); //DetId does not kmow about halfShells in TIB
289  if ( (rawId>=369120484 && rawId<369120688) || (rawId>=369121540 && rawId<369121776) || (rawId>=369136932 && rawId<369137200) || (rawId>=369137988 && rawId<369138288) ||
290  (rawId>=369153396 && rawId<369153744) || (rawId>=369154436 && rawId<369154800) || (rawId>=369169844 && rawId<369170256) || (rawId>=369170900 && rawId<369171344) ||
291  (rawId>=369124580 && rawId<369124784) || (rawId>=369125636 && rawId<369125872) || (rawId>=369141028 && rawId<369141296) || (rawId>=369142084 && rawId<369142384) ||
292  (rawId>=369157492 && rawId<369157840) || (rawId>=369158532 && rawId<369158896) || (rawId>=369173940 && rawId<369174352) || (rawId>=369174996 && rawId<369175440) ) whichHalfShell=2;
293  treeMem.layer = tTopo->tibLayer(detId);
294  treeMem.side = tTopo->tibStringInfo(detId)[0];
295  treeMem.half = whichHalfShell;
296  treeMem.rod = tTopo->tibStringInfo(detId)[2];
297  treeMem.outerInner = tTopo->tibStringInfo(detId)[1];
298  treeMem.module = tTopo->tibModule(detId);
299  treeMem.isStereo = tTopo->tibStereo(detId);
300  treeMem.isDoubleSide = tTopo->tibIsDoubleSide(detId);
301  } else if(treeMem.subDetId == StripSubdetector::TID){
302 
303  treeMem.layer = tTopo->tidWheel(detId);
304  treeMem.side = tTopo->tidSide(detId);
305  treeMem.ring = tTopo->tidRing(detId);
306  treeMem.outerInner = tTopo->tidModuleInfo(detId)[0];
307  treeMem.module = tTopo->tidModuleInfo(detId)[1];
308  treeMem.isStereo = tTopo->tidStereo(detId);
309  treeMem.isDoubleSide = tTopo->tidIsDoubleSide(detId);
310  } else if(treeMem.subDetId == StripSubdetector::TOB){
311 
312  treeMem.layer = tTopo->tobLayer(detId);
313  treeMem.side = tTopo->tobRodInfo(detId)[0];
314  treeMem.rod = tTopo->tobRodInfo(detId)[1];
315  treeMem.module = tTopo->tobModule(detId);
316  treeMem.isStereo = tTopo->tobStereo(detId);
317  treeMem.isDoubleSide = tTopo->tobIsDoubleSide(detId);
318  } else if(treeMem.subDetId == StripSubdetector::TEC) {
319 
320  treeMem.layer = tTopo->tecWheel(detId);
321  treeMem.side = tTopo->tecSide(detId);
322  treeMem.ring = tTopo->tecRing(detId);
323  treeMem.petal = tTopo->tecPetalInfo(detId)[1];
324  treeMem.outerInner = tTopo->tecPetalInfo(detId)[0];
325  treeMem.module = tTopo->tecModule(detId);
326  treeMem.isStereo = tTopo->tecStereo(detId);
327  treeMem.isDoubleSide = tTopo->tecIsDoubleSide(detId);
328  }
329 
330  //variables concerning the tracker geometry
331 
332  const Surface::PositionType &gPModule = tkgeom.idToDet(detId)->position();
333  treeMem.posPhi = gPModule.phi();
334  treeMem.posEta = gPModule.eta();
335  treeMem.posR = gPModule.perp();
336  treeMem.posX = gPModule.x();
337  treeMem.posY = gPModule.y();
338  treeMem.posZ = gPModule.z();
339 
340  const Surface& surface = tkgeom.idToDet(detId)->surface();
341 
342  //global Orientation of local coordinate system of dets/detUnits
343  LocalPoint lUDirection(1.,0.,0.), lVDirection(0.,1.,0.), lWDirection(0.,0.,1.);
344  GlobalPoint gUDirection = surface.toGlobal(lUDirection),
345  gVDirection = surface.toGlobal(lVDirection),
346  gWDirection = surface.toGlobal(lWDirection);
347  double dR(999.), dPhi(999.), dZ(999.);
349  dR = gWDirection.perp() - gPModule.perp();
350  dPhi = deltaPhi(gUDirection.phi(),gPModule.phi());
351  dZ = gVDirection.z() - gPModule.z();
352  if(dZ>=0.)treeMem.rOrZDirection = 1; else treeMem.rOrZDirection = -1;
353  }else if(treeMem.subDetId==PixelSubdetector::PixelEndcap){
354  dR = gUDirection.perp() - gPModule.perp();
355  dPhi = deltaPhi(gVDirection.phi(),gPModule.phi());
356  dZ = gWDirection.z() - gPModule.z();
357  if(dR>=0.)treeMem.rOrZDirection = 1; else treeMem.rOrZDirection = -1;
358  }else if(treeMem.subDetId==StripSubdetector::TID || treeMem.subDetId==StripSubdetector::TEC){
359  dR = gVDirection.perp() - gPModule.perp();
360  dPhi = deltaPhi(gUDirection.phi(),gPModule.phi());
361  dZ = gWDirection.z() - gPModule.z();
362  if(dR>=0.)treeMem.rOrZDirection = 1; else treeMem.rOrZDirection = -1;
363  }
364  if(dR>=0.)treeMem.rDirection = 1; else treeMem.rDirection = -1;
365  if(dPhi>=0.)treeMem.phiDirection = 1; else treeMem.phiDirection = -1;
366  if(dZ>=0.)treeMem.zDirection = 1; else treeMem.zDirection = -1;
367 
368 
369  // Assign histos from first step (TrackerOfflineValidation) to the module's entry in the TTree for retrieving mean, rms, median ...
370  const std::string histDir = associateModuleHistsWithTree(treeMem, it->second, substructureName);
371 
372 
373  //mean and RMS values (extracted from histograms Xprime on module level)
374  treeMem.entries = static_cast<UInt_t>(it->second.ResXprimeHisto->GetEntries());
375  treeMem.meanX = it->second.ResXprimeHisto->GetMean();
376  treeMem.rmsX = it->second.ResXprimeHisto->GetRMS();
377  //treeMem.sigmaX = Fwhm(it->second.ResXprimeHisto)/2.355;
378  if (useFit_) {
379  //call fit function which returns mean and sigma from the fit
380  //for absolute residuals
381  std::pair<float,float> fitResult1 = this->fitResiduals(it->second.ResXprimeHisto);
382  treeMem.fitMeanX = fitResult1.first;
383  treeMem.fitSigmaX = fitResult1.second;
384  //for normalized residuals
385  std::pair<float,float> fitResult2 = this->fitResiduals(it->second.NormResXprimeHisto);
386  treeMem.fitMeanNormX = fitResult2.first;
387  treeMem.fitSigmaNormX = fitResult2.second;
388  }
389 
390  //get median for absolute residuals
391  treeMem.medianX = this->getMedian(it->second.ResXprimeHisto);
392 
393  int numberOfBins=it->second.ResXprimeHisto->GetNbinsX();
394  treeMem.numberOfUnderflows = it->second.ResXprimeHisto->GetBinContent(0);
395  treeMem.numberOfOverflows = it->second.ResXprimeHisto->GetBinContent(numberOfBins+1);
396  treeMem.numberOfOutliers = it->second.ResXprimeHisto->GetBinContent(0)+it->second.ResXprimeHisto->GetBinContent(numberOfBins+1);
397  //mean and RMS values (extracted from histograms(normalized Xprime on module level)
398  treeMem.meanNormX = it->second.NormResXprimeHisto->GetMean();
399  treeMem.rmsNormX = it->second.NormResXprimeHisto->GetRMS();
400 
401  double stats[20];
402  it->second.NormResXprimeHisto->GetStats(stats);
403  // GF treeMem.chi2PerDofX = stats[3]/(stats[0]-1);
404  if (stats[0]) treeMem.chi2PerDofX = stats[3]/stats[0];
405 
406  //treeMem.sigmaNormX = Fwhm(it->second.NormResXprimeHisto)/2.355;
407  treeMem.histNameX = it->second.ResXprimeHisto->GetName();
408  treeMem.histNameNormX = it->second.NormResXprimeHisto->GetName();
409 
410 
411  // fill tree variables in local coordinates if set in cfg of TrackerOfllineValidation
412  if(it->second.ResHisto && it->second.NormResHisto){ // if(lCoorHistOn_) {
413  treeMem.meanLocalX = it->second.ResHisto->GetMean();
414  treeMem.rmsLocalX = it->second.ResHisto->GetRMS();
415  treeMem.meanNormLocalX = it->second.NormResHisto->GetMean();
416  treeMem.rmsNormLocalX = it->second.NormResHisto->GetRMS();
417  treeMem.histNameLocalX = it->second.ResHisto->GetName();
418  treeMem.histNameNormLocalX = it->second.NormResHisto->GetName();
419  }
420 
421  // mean and RMS values in local y (extracted from histograms Yprime on module level)
422  // might exist in pixel only
423  if (it->second.ResYprimeHisto) { //(stripYResiduals_){
424  TH1 *h = it->second.ResYprimeHisto;
425  treeMem.meanY = h->GetMean();
426  treeMem.rmsY = h->GetRMS();
427 
428  if (useFit_) { // fit function which returns mean and sigma from the fit
429  std::pair<float,float> fitMeanSigma = this->fitResiduals(h);
430  treeMem.fitMeanY = fitMeanSigma.first;
431  treeMem.fitSigmaY = fitMeanSigma.second;
432  }
433 
434  //get median for absolute residuals
435  treeMem.medianY = this->getMedian(h);
436 
437  treeMem.histNameY = h->GetName();
438  }
439 
440  if (it->second.NormResYprimeHisto) {
441  TH1 *h = it->second.NormResYprimeHisto;
442  treeMem.meanNormY = h->GetMean();
443  treeMem.rmsNormY = h->GetRMS();
444  h->GetStats(stats); // stats buffer defined above
445  if (stats[0]) treeMem.chi2PerDofY = stats[3]/stats[0];
446 
447  if (useFit_) { // fit function which returns mean and sigma from the fit
448  std::pair<float,float> fitMeanSigma = this->fitResiduals(h);
449  treeMem.fitMeanNormY = fitMeanSigma.first;
450  treeMem.fitSigmaNormY = fitMeanSigma.second;
451  }
452  treeMem.histNameNormY = h->GetName();
453  }
454 
455  // Delete module level hists if set in cfg
456  const bool removeModuleLevelHists(parSet_.getParameter<bool>("removeModuleLevelHists"));
457  if(removeModuleLevelHists){
458  dbe_->setCurrentFolder("");
459  if(it->second.ResHisto) dbe_->removeElement(histDir + "/" + it->second.ResHisto->GetName());
460  if(it->second.NormResHisto) dbe_->removeElement(histDir + "/" + it->second.NormResHisto->GetName());
461  if(it->second.ResXprimeHisto) dbe_->removeElement(histDir + "/" + it->second.ResXprimeHisto->GetName());
462  if(it->second.NormResXprimeHisto)dbe_->removeElement(histDir + "/" + it->second.NormResXprimeHisto->GetName());
463  if(it->second.ResYprimeHisto) dbe_->removeElement(histDir + "/" + it->second.ResYprimeHisto->GetName());
464  if(it->second.NormResYprimeHisto)dbe_->removeElement(histDir + "/" + it->second.NormResYprimeHisto->GetName());
465  }
466 
467  tree.Fill();
468  }
469 }
470 
471 
472 const std::string
473 TrackerOfflineValidationSummary::associateModuleHistsWithTree(const TkOffTreeVariables& treeMem, TrackerOfflineValidationSummary::ModuleHistos& moduleHists, std::map<std::string,std::string>& substructureName){
474  std::stringstream histDir, sSubdetName;
475  std::string componentName;
476  if(moduleDirectory_.length() != 0)histDir<<moduleDirectory_<<"/";
477  std::string wheelOrLayer("_layer_");
479  unsigned int half(treeMem.half), layer(treeMem.layer), ladder(0);
480  if(layer==1){
481  if(half==2)ladder = treeMem.rod -5;
482  else if(treeMem.rod>15)ladder = treeMem.rod -10;
483  else ladder = treeMem.rod;
484  }else if(layer==2){
485  if(half==2)ladder = treeMem.rod -8;
486  else if(treeMem.rod>24)ladder = treeMem.rod -16;
487  else ladder = treeMem.rod;
488  }else if(layer==3){
489  if(half==2)ladder = treeMem.rod -11;
490  else if(treeMem.rod>33)ladder = treeMem.rod -22;
491  else ladder = treeMem.rod;
492  }
493  componentName = "Pixel";
494  sSubdetName<<"TPBBarrel_1";
495  histDir<<componentName<<"/"<<sSubdetName.str()<<"/TPBHalfBarrel_"<<treeMem.half<<"/TPBLayer_"<<treeMem.layer<<"/TPBLadder_"<<ladder;
496  }else if(treeMem.subDetId == PixelSubdetector::PixelEndcap){
497  unsigned int side(treeMem.side), half(treeMem.half), blade(0);
498  if(side==1)side=3;
499  if(half==2)blade = treeMem.blade -6;
500  else if(treeMem.blade>18)blade = treeMem.blade -12;
501  else blade = treeMem.blade;
502  componentName = "Pixel";
503  sSubdetName<<"TPEEndcap_"<<side;
504  histDir<<componentName<<"/"<<sSubdetName.str()<<"/TPEHalfCylinder_"<<treeMem.half<<"/TPEHalfDisk_"<<treeMem.layer<<"/TPEBlade_"<<blade<<"/TPEPanel_"<<treeMem.panel;
505  wheelOrLayer = "_wheel_";
506  }else if(treeMem.subDetId == StripSubdetector::TIB){
507  unsigned int half(treeMem.half), layer(treeMem.layer), surface(treeMem.outerInner), string(0);
508  if(half==2){
509  if(layer==1){
510  if(surface==1)string = treeMem.rod -13;
511  else if(surface==2)string = treeMem.rod -15;
512  }
513  if(layer==2){
514  if(surface==1)string = treeMem.rod -17;
515  else if(surface==2)string = treeMem.rod -19;
516  }
517  if(layer==3){
518  if(surface==1)string = treeMem.rod -22;
519  else if(surface==2)string = treeMem.rod -23;
520  }
521  if(layer==4){
522  if(surface==1)string = treeMem.rod -26;
523  else if(surface==2)string = treeMem.rod -28;
524  }
525  }
526  else string = treeMem.rod;
527  std::stringstream detString;
528  if(treeMem.layer<3 && !treeMem.isDoubleSide)detString<<"/Det_"<<treeMem.module;
529  else detString<<"";
530  componentName = "Strip";
531  sSubdetName<<"TIBBarrel_1";
532  histDir<<componentName<<"/"<<sSubdetName.str()<<"/TIBHalfBarrel_"<<treeMem.side<<"/TIBLayer_"<<treeMem.layer<<"/TIBHalfShell_"<<treeMem.half<<"/TIBSurface_"<<treeMem.outerInner<<"/TIBString_"<<string<<detString.str();
533  }else if(treeMem.subDetId == StripSubdetector::TID){
534  unsigned int side(treeMem.side), outerInner(0);
535  if(side==1)side=3;
536  if(treeMem.outerInner==1)outerInner=2;
537  else if(treeMem.outerInner==2)outerInner=1;
538  std::stringstream detString;
539  if(treeMem.ring<3 && !treeMem.isDoubleSide)detString<<"/Det_"<<treeMem.module;
540  else detString<<"";
541  componentName = "Strip";
542  sSubdetName<<"TIDEndcap_"<<side;
543  histDir<<componentName<<"/"<<sSubdetName.str()<<"/TIDDisk_"<<treeMem.layer<<"/TIDRing_"<<treeMem.ring<<"/TIDSide_"<<outerInner<<detString.str();
544  wheelOrLayer = "_wheel_";
545  }else if(treeMem.subDetId == StripSubdetector::TOB){
546  std::stringstream detString;
547  if(treeMem.layer<3 && !treeMem.isDoubleSide)detString<<"/Det_"<<treeMem.module;
548  else detString<<"";
549  componentName = "Strip";
550  sSubdetName<<"TOBBarrel_4";
551  histDir<<componentName<<"/"<<sSubdetName.str()<<"/TOBHalfBarrel_"<<treeMem.side<<"/TOBLayer_"<<treeMem.layer<<"/TOBRod_"<<treeMem.rod<<detString.str();
552  }else if(treeMem.subDetId == StripSubdetector::TEC) {
553  unsigned int side(0), outerInner(0), ring(0);
554  if(treeMem.side==1)side=6;
555  else if(treeMem.side==2)side=5;
556  if(treeMem.outerInner==1)outerInner = 2;
557  else if(treeMem.outerInner==2)outerInner=1;
558  if(treeMem.layer>3 && treeMem.layer<7)ring = treeMem.ring -1;
559  else if(treeMem.layer==7 || treeMem.layer==8)ring = treeMem.ring -2;
560  else if(treeMem.layer==9)ring = treeMem.ring -3;
561  else ring = treeMem.ring;
562  std::stringstream detString;
563  if((treeMem.ring<3 || treeMem.ring==5) && !treeMem.isDoubleSide)detString<<"/Det_"<<treeMem.module;
564  else detString<<"";
565  componentName = "Strip";
566  sSubdetName<<"TECEndcap_"<<side;
567  histDir<<componentName<<"/"<<sSubdetName.str()<<"/TECDisk_"<<treeMem.layer<<"/TECSide_"<<outerInner<<"/TECPetal_"<<treeMem.petal<<"/TECRing_"<<ring<<detString.str();
568  wheelOrLayer = "_wheel_";
569  }
570 
571  substructureName["component"] = componentName;
572  substructureName["subdet"] = sSubdetName.str();
573 
574  std::stringstream histName;
575  histName<<"residuals_subdet_"<<treeMem.subDetId<<wheelOrLayer<<treeMem.layer<<"_module_"<<treeMem.moduleId;
576 
578  fullPath = histDir.str()+"/h_xprime_"+histName.str();
579  if(dbe_->get(fullPath)) moduleHists.ResXprimeHisto = dbe_->get(fullPath)->getTH1();
580  else{edm::LogError("TrackerOfflineValidationSummary")<<"Problem with names in input file produced in TrackerOfflineValidation ...\n"
581  <<"This histogram should exist in every configuration, "
582  <<"but no histogram with name "<<fullPath<<" is found!";
583  return "";
584  }
585  fullPath = histDir.str()+"/h_normxprime"+histName.str();
586  if(dbe_->get(fullPath)) moduleHists.NormResXprimeHisto = dbe_->get(fullPath)->getTH1();
587  fullPath = histDir.str()+"/h_yprime_"+histName.str();
588  if(dbe_->get(fullPath)) moduleHists.ResYprimeHisto = dbe_->get(fullPath)->getTH1();
589  fullPath = histDir.str()+"/h_normyprime"+histName.str();
590  if(dbe_->get(fullPath)) moduleHists.NormResYprimeHisto = dbe_->get(fullPath)->getTH1();
591  fullPath = histDir.str()+"/h_"+histName.str();
592  if(dbe_->get(fullPath)) moduleHists.ResHisto = dbe_->get(fullPath)->getTH1();
593  fullPath = histDir.str()+"/h_norm"+histName.str();
594  if(dbe_->get(fullPath)) moduleHists.NormResHisto = dbe_->get(fullPath)->getTH1();
595 
596  return histDir.str();
597 }
598 
599 
600 std::pair<float,float>
602 {
603  std::pair<float,float> fitResult(9999., 9999.);
604  if (!hist || hist->GetEntries() < 20) return fitResult;
605 
606  float mean = hist->GetMean();
607  float sigma = hist->GetRMS();
608 
609  try { // for < CMSSW_2_2_0 since ROOT warnings from fit are converted to exceptions
610  // Remove the try/catch for more recent CMSSW!
611  // first fit: two RMS around mean
612  TF1 func("tmp", "gaus", mean - 2.*sigma, mean + 2.*sigma);
613  if (0 == hist->Fit(&func,"QNR")) { // N: do not blow up file by storing fit!
614  mean = func.GetParameter(1);
615  sigma = func.GetParameter(2);
616  // second fit: three sigma of first fit around mean of first fit
617  func.SetRange(mean - 3.*sigma, mean + 3.*sigma);
618  // I: integral gives more correct results if binning is too wide
619  // L: Likelihood can treat empty bins correctly (if hist not weighted...)
620  if (0 == hist->Fit(&func, "Q0LR")) {
621  if (hist->GetFunction(func.GetName())) { // Take care that it is later on drawn:
622  hist->GetFunction(func.GetName())->ResetBit(TF1::kNotDraw);
623  }
624  fitResult.first = func.GetParameter(1);
625  fitResult.second = func.GetParameter(2);
626  }
627  }
628  } catch (cms::Exception const & e) {
629  edm::LogWarning("Alignment") << "@SUB=TrackerOfflineValidation::fitResiduals"
630  << "Caught this exception during ROOT fit: "
631  << e.what();
632  }
633  return fitResult;
634 }
635 
636 
637 float
639 {
640  float median = 999;
641  const int nbins = histo->GetNbinsX();
642 
643  //extract median from histogram
644  double *x = new double[nbins];
645  double *y = new double[nbins];
646  for (int j = 0; j < nbins; j++) {
647  x[j] = histo->GetBinCenter(j+1);
648  y[j] = histo->GetBinContent(j+1);
649  }
650  median = TMath::Median(nbins, x, y);
651 
652  delete[] x; x = 0;
653  delete [] y; y = 0;
654 
655  return median;
656 }
657 
658 
659 void
661 {
662  this->applyHarvestingHierarchy(tree);
663  this->fillHarvestingHists(tree);
664 }
665 
666 
667 void
669 {
670  TkOffTreeVariables *treeMemPtr = 0;
671  std::map<std::string,std::string> *substructureName = 0;
672  tree.SetBranchAddress("TkOffTreeVariables", &treeMemPtr);
673  tree.SetBranchAddress("SubstructureName", &substructureName);
674 
675  // Loop over modules to select accumulation criteria for harvesting plots
676  for(unsigned int iSubdet = 1; iSubdet<7; ++iSubdet){
677  std::string hierarchyName("");
678  std::string componentName("");
679  std::vector<unsigned int> treeEntries;
680  for(unsigned int iSide = 1; iSide<3; ++iSide){
681  // Set up only one collection for Barrels, not separated for side
682  if(iSide==1 && (iSubdet==PixelSubdetector::PixelBarrel || iSubdet==StripSubdetector::TIB || iSubdet==StripSubdetector::TOB))continue;
683  for(int iTree=0; iTree<tree.GetEntries(); ++iTree){
684  tree.GetEntry(iTree);
685  // Do not use glued Dets
686  if(treeMemPtr->isDoubleSide)continue;
687  if(treeMemPtr->subDetId == iSubdet){
688  if(iSide!=treeMemPtr->side && (iSubdet==PixelSubdetector::PixelEndcap || iSubdet==StripSubdetector::TID || iSubdet==StripSubdetector::TEC))continue;
689  treeEntries.push_back(iTree);
690  if(hierarchyName.length() == 0){
691  hierarchyName = (*substructureName)["subdet"];
692  componentName = (*substructureName)["component"];
693  }
694  }
695  }
696  HarvestingHierarchy harvestingHierarchy(hierarchyName,componentName,treeEntries);
697  vHarvestingHierarchy_.push_back(harvestingHierarchy);
698  hierarchyName = ""; componentName = ""; treeEntries.clear();
699  }
700  }
701  // Here could be a further separation of the HarvestingHierarchy.
702  // E.g. separate the existing ones by layer and add them to the vector without deleting any element from the vector.
703  // The existing hists will stay and the new ones are added
704 
705  // Now, book the corresponding histos
706  this->bookHarvestingHists();
707 }
708 
709 
710 void
712 {
713  edm::LogInfo("TrackerOfflineValidationSummary")<<"Harvesting histograms will be booked for "<<vHarvestingHierarchy_.size()<<" different hierarchy selections";
714  for(std::vector<HarvestingHierarchy>::iterator iHier = vHarvestingHierarchy_.begin(); iHier != vHarvestingHierarchy_.end(); ++iHier){
715 
716  std::stringstream dmrXprimeHistoName, dmrYprimeHistoName, dmrXprimeHistoTitle, dmrYprimeHistoTitle;
717  dmrXprimeHistoName << "h_DmrXprime_" << iHier->hierarchyName;
718  dmrYprimeHistoName << "h_DmrYprime_" << iHier->hierarchyName;
719  dmrXprimeHistoTitle<< "DMR for " << iHier->hierarchyName <<";<#DeltaX> [cm];# modules";
720  dmrYprimeHistoTitle<< "DMR for " << iHier->hierarchyName <<";<#DeltaY> [cm];# modules";
721 
722  std::string directoryString(moduleDirectory_);
723  if(directoryString.length()!=0)directoryString += "/";
724  directoryString += iHier->componentName;
725  dbe_->setCurrentFolder(directoryString);
726 
727  int nBinsX(0); double xMin(0.), xMax(0.);
728  if(iHier->componentName == "Pixel"){
729  this->getBinning("TH1DmrXprimePixelModules",nBinsX,xMin,xMax);
730  iHier->harvestingHistos.DmrXprime = dbe_->book1D(dmrXprimeHistoName.str(),dmrXprimeHistoTitle.str(),nBinsX,xMin,xMax)->getTH1();
731  this->getBinning("TH1DmrYprimePixelModules",nBinsX,xMin,xMax);
732  iHier->harvestingHistos.DmrYprime = dbe_->book1D(dmrYprimeHistoName.str(),dmrYprimeHistoTitle.str(),nBinsX,xMin,xMax)->getTH1();
733  }
734  else if(iHier->componentName == "Strip"){
735  this->getBinning("TH1DmrXprimeStripModules",nBinsX,xMin,xMax);
736  iHier->harvestingHistos.DmrXprime = dbe_->book1D(dmrXprimeHistoName.str(),dmrXprimeHistoTitle.str(),nBinsX,xMin,xMax)->getTH1();
737  if(!parSet_.getParameter<bool>("stripYDmrs"))continue;
738  this->getBinning("TH1DmrYprimeStripModules",nBinsX,xMin,xMax);
739  iHier->harvestingHistos.DmrYprime = dbe_->book1D(dmrYprimeHistoName.str(),dmrYprimeHistoTitle.str(),nBinsX,xMin,xMax)->getTH1();
740  }
741  }
742 }
743 
744 
745 void
746 TrackerOfflineValidationSummary::getBinning(const std::string& binningPSetName, int& nBinsX, double& lowerBoundX, double& upperBoundX)const
747 {
748  const edm::ParameterSet& binningPSet = parSet_.getParameter<edm::ParameterSet>(binningPSetName);
749  nBinsX = binningPSet.getParameter<int>("Nbinx");
750  lowerBoundX = binningPSet.getParameter<double>("xmin");
751  upperBoundX = binningPSet.getParameter<double>("xmax");
752 }
753 
754 
755 void
757 {
758  TkOffTreeVariables *treeMemPtr = 0;
759  std::map<std::string,std::string> *substructureName = 0;
760  tree.SetBranchAddress("TkOffTreeVariables", &treeMemPtr);
761  tree.SetBranchAddress("SubstructureName", &substructureName);
762 
763  const unsigned int minEntriesPerModule(parSet_.getParameter<unsigned int>("minEntriesPerModuleForDmr"));
764  edm::LogInfo("TrackerOfflineValidationSummary")<<"Median of a module is added to DMR plots if it contains at least "<<minEntriesPerModule<<" hits";
765 
766  for(std::vector<HarvestingHierarchy>::iterator iHier = vHarvestingHierarchy_.begin(); iHier != vHarvestingHierarchy_.end(); ++iHier){
767  for(std::vector<unsigned int>::const_iterator iTreeEntries = iHier->treeEntries.begin(); iTreeEntries != iHier->treeEntries.end(); ++iTreeEntries){
768  tree.GetEntry(*iTreeEntries);
769  if(treeMemPtr->entries < minEntriesPerModule)continue;
770  iHier->harvestingHistos.DmrXprime->Fill(treeMemPtr->medianX);
771  if(iHier->harvestingHistos.DmrYprime)iHier->harvestingHistos.DmrYprime->Fill(treeMemPtr->medianY);
772  }
773  }
774 }
775 
776 
777 //define this as a plug-in
virtual char const * what() const
Definition: Exception.cc:141
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:114
T getParameter(std::string const &) const
bool tecIsDoubleSide(const DetId &id) const
bool tobIsDoubleSide(const DetId &id) const
bool tibIsDoubleSide(const DetId &id) const
T perp() const
Definition: PV3DBase.h:72
unsigned int tibLayer(const DetId &id) const
unsigned int tidRing(const DetId &id) const
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:971
Put here the histograms created during harvesting.
std::vector< unsigned int > tidModuleInfo(const DetId &id) const
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
unsigned int pxfDisk(const DetId &id) const
container to hold data to be written into TTree
unsigned int tecRing(const DetId &id) const
ring id
std::map< int, TrackerOfflineValidationSummary::ModuleHistos > mPxbResiduals_
uint32_t tobStereo(const DetId &id) const
std::vector< HarvestingHierarchy > vHarvestingHierarchy_
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
unsigned int pxbLadder(const DetId &id) const
T y() const
Definition: PV3DBase.h:63
unsigned int tidWheel(const DetId &id) const
unsigned int pxbModule(const DetId &id) const
edm::ESHandle< TrackerGeometry > tkGeom_
#define nullptr
std::vector< unsigned int > tibStringInfo(const DetId &id) const
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:40
void getBinning(const std::string &binningPSetName, int &nBinsX, double &lowerBoundX, double &upperBoundX) const
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
virtual const DetIdContainer & detIds() const
Returm a vector of all GeomDet DetIds (including those of GeomDetUnits)
std::map< int, TrackerOfflineValidationSummary::ModuleHistos > mTibResiduals_
int iEvent
Definition: GenABIO.cc:230
std::vector< unsigned int > tecPetalInfo(const DetId &id) const
std::map< int, TrackerOfflineValidationSummary::ModuleHistos > mTobResiduals_
const Surface::PositionType & position() const
The position (origin of the R.F.)
Definition: GeomDet.h:46
unsigned int tidSide(const DetId &id) const
std::map< int, TrackerOfflineValidationSummary::ModuleHistos > mPxeResiduals_
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
uint32_t tidStereo(const DetId &id) const
std::vector< unsigned int > tobRodInfo(const DetId &id) const
const std::string associateModuleHistsWithTree(const TkOffTreeVariables &treeMem, TrackerOfflineValidationSummary::ModuleHistos &moduleHists, std::map< std::string, std::string > &substructureName)
T z() const
Definition: PV3DBase.h:64
void removeElement(const std::string &name)
Definition: DQMStore.cc:3166
int j
Definition: DBlmapReader.cc:9
TH1 * getTH1(void) const
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
MonitorElement * get(const std::string &path) const
get ME from full pathname (e.g. &quot;my/long/dir/my_histo&quot;)
Definition: DQMStore.cc:1725
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
std::pair< float, float > fitResiduals(TH1 *hist) const
TrackerOfflineValidationSummary(const edm::ParameterSet &)
virtual void analyze(const edm::Event &evt, const edm::EventSetup &) override
DQMStore * dbe_
std::string histNameNormLocalX
unsigned int tibModule(const DetId &id) const
unsigned int pxfModule(const DetId &id) const
std::map< int, TrackerOfflineValidationSummary::ModuleHistos > mTidResiduals_
unsigned int pxbLayer(const DetId &id) const
unsigned int tecModule(const DetId &id) const
Definition: DetId.h:18
HarvestingHierarchy(const std::string &name, const std::string &component, const std::vector< unsigned int > &entries)
std::string histNameLocalX
const T & get() const
Definition: EventSetup.h:55
bool tidIsDoubleSide(const DetId &id) const
T const * product() const
Definition: ESHandle.h:86
std::map< int, TrackerOfflineValidationSummary::ModuleHistos > mTecResiduals_
unsigned int tobModule(const DetId &id) const
T eta() const
Definition: PV3DBase.h:76
void clear()
set to empty values
uint32_t tecStereo(const DetId &id) const
unsigned int pxfSide(const DetId &id) const
volatile std::atomic< bool > shutdown_flag false
Definition: DDAxes.h:10
uint32_t tibStereo(const DetId &id) const
Detector det() const
get the detector field from this detid
Definition: DetId.h:35
T x() const
Definition: PV3DBase.h:62
void fillTree(TTree &tree, std::map< int, TrackerOfflineValidationSummary::ModuleHistos > &moduleHist, TkOffTreeVariables &treeMem, const TrackerGeometry &tkgeom, std::map< std::string, std::string > &substructureName, const TrackerTopology *tTopo)
unsigned int tecWheel(const DetId &id) const
unsigned int pxfPanel(const DetId &id) const
unsigned int pxfBlade(const DetId &id) const
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:684
unsigned int tobLayer(const DetId &id) const
unsigned int tecSide(const DetId &id) const
virtual const TrackerGeomDet * idToDet(DetId) const
std::vector< DetId > DetIdContainer