Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 12:45:58 

0001 // -*- C++ -*-
0002 //
0003 // Package:    TrackerOfflineValidation
0004 // Class:      TrackerOfflineValidation
0005 //
0006 /**\class TrackerOfflineValidation TrackerOfflineValidation.cc Alignment/Validator/src/TrackerOfflineValidation.cc
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Erik Butz
0015 //         Created:  Tue Dec 11 14:03:05 CET 2007
0016 // $Id: TrackerOfflineValidation.cc,v 1.55 2012/09/28 20:48:06 flucke Exp $
0017 //
0018 //
0019 
0020 // system include files
0021 #include <memory>
0022 #include <map>
0023 #include <sstream>
0024 #include <cmath>
0025 #include <tuple>
0026 #include <utility>
0027 #include <vector>
0028 #include <iostream>
0029 
0030 // ROOT includes
0031 #include "TH1.h"
0032 #include "TH2.h"
0033 #include "TProfile.h"
0034 #include "TFile.h"
0035 #include "TTree.h"
0036 #include "TF1.h"
0037 #include "TMath.h"
0038 
0039 // user include files
0040 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0041 #include "FWCore/Framework/interface/ConsumesCollector.h"
0042 #include "FWCore/Framework/interface/Frameworkfwd.h"
0043 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0044 #include "FWCore/Framework/interface/EventSetup.h"
0045 #include "FWCore/Framework/interface/Event.h"
0046 #include "FWCore/Framework/interface/MakerMacros.h"
0047 #include "FWCore/ServiceRegistry/interface/Service.h"
0048 
0049 #include "DataFormats/DetId/interface/DetId.h"
0050 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0051 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0052 #include "DataFormats/Math/interface/deltaPhi.h"
0053 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0054 
0055 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0056 #include "CommonTools/Utils/interface/TFileDirectory.h"
0057 
0058 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0059 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0060 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0061 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0062 
0063 #include "DQMServices/Core/interface/DQMStore.h"
0064 
0065 #include "Alignment/CommonAlignment/interface/Alignable.h"
0066 #include "Alignment/CommonAlignment/interface/Utilities.h"
0067 #include "Alignment/CommonAlignment/interface/AlignableObjectId.h"
0068 #include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
0069 #include "Alignment/TrackerAlignment/interface/TrackerAlignableId.h"
0070 
0071 #include "Alignment/OfflineValidation/interface/TrackerValidationVariables.h"
0072 #include "Alignment/OfflineValidation/interface/TkOffTreeVariables.h"
0073 
0074 //
0075 // class declaration
0076 //
0077 class TrackerOfflineValidation : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0078 public:
0079   typedef dqm::legacy::DQMStore DQMStore;
0080   explicit TrackerOfflineValidation(const edm::ParameterSet&);
0081   ~TrackerOfflineValidation() override;
0082 
0083   enum HistogramType {
0084     XResidual,
0085     NormXResidual,
0086     YResidual, /*NormYResidual, */
0087     XprimeResidual,
0088     NormXprimeResidual,
0089     YprimeResidual,
0090     NormYprimeResidual,
0091     XResidualProfile,
0092     YResidualProfile
0093   };
0094 
0095 private:
0096   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
0097   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0098 
0099   struct ModuleHistos {
0100     ModuleHistos()
0101         : ResHisto(),
0102           NormResHisto(),
0103           ResYHisto(), /*NormResYHisto(),*/
0104           ResXprimeHisto(),
0105           NormResXprimeHisto(),
0106           ResYprimeHisto(),
0107           NormResYprimeHisto(),
0108           ResXvsXProfile(),
0109           ResXvsYProfile(),
0110           ResYvsXProfile(),
0111           ResYvsYProfile(),
0112           LocalX(),
0113           LocalY() {}
0114     TH1* ResHisto;
0115     TH1* NormResHisto;
0116     TH1* ResYHisto;
0117     /* TH1* NormResYHisto; */
0118     TH1* ResXprimeHisto;
0119     TH1* NormResXprimeHisto;
0120     TH1* ResYprimeHisto;
0121     TH1* NormResYprimeHisto;
0122 
0123     TProfile* ResXvsXProfile;
0124     TProfile* ResXvsYProfile;
0125     TProfile* ResYvsXProfile;
0126     TProfile* ResYvsYProfile;
0127 
0128     TH1* LocalX;
0129     TH1* LocalY;
0130   };
0131 
0132   // container struct to organize collection of histograms during endJob
0133   struct SummaryContainer {
0134     SummaryContainer()
0135         : sumXResiduals_(),
0136           summaryXResiduals_(),
0137           sumNormXResiduals_(),
0138           summaryNormXResiduals_(),
0139           sumYResiduals_(),
0140           summaryYResiduals_(),
0141           sumNormYResiduals_(),
0142           summaryNormYResiduals_(),
0143           sumResXvsXProfile_(),
0144           sumResXvsYProfile_(),
0145           sumResYvsXProfile_(),
0146           sumResYvsYProfile_() {}
0147 
0148     TH1* sumXResiduals_;
0149     TH1* summaryXResiduals_;
0150     TH1* sumNormXResiduals_;
0151     TH1* summaryNormXResiduals_;
0152     TH1* sumYResiduals_;
0153     TH1* summaryYResiduals_;
0154     TH1* sumNormYResiduals_;
0155     TH1* summaryNormYResiduals_;
0156 
0157     TH1* sumResXvsXProfile_;
0158     TH1* sumResXvsYProfile_;
0159     TH1* sumResYvsXProfile_;
0160     TH1* sumResYvsYProfile_;
0161   };
0162 
0163   struct DirectoryWrapper {
0164     DirectoryWrapper(const DirectoryWrapper& upDir,
0165                      const std::string& newDir,
0166                      const std::string& basedir,
0167                      bool useDqmMode)
0168         : tfd(nullptr), dqmMode(useDqmMode), theDbe(nullptr) {
0169       if (newDir.length() != 0) {
0170         if (upDir.directoryString.length() != 0)
0171           directoryString = upDir.directoryString + "/" + newDir;
0172         else
0173           directoryString = newDir;
0174       } else
0175         directoryString = upDir.directoryString;
0176 
0177       if (!dqmMode) {
0178         if (newDir.length() == 0)
0179           tfd.reset(&(*upDir.tfd));
0180         else
0181           tfd = std::make_unique<TFileDirectory>(upDir.tfd->mkdir(newDir));
0182       } else {
0183         theDbe = edm::Service<DQMStore>().operator->();
0184       }
0185     }
0186 
0187     DirectoryWrapper(const std::string& newDir, const std::string& basedir, bool useDqmMode)
0188         : tfd(nullptr), dqmMode(useDqmMode), theDbe(nullptr) {
0189       if (!dqmMode) {
0190         edm::Service<TFileService> fs;
0191         if (newDir.length() == 0) {
0192           tfd = std::make_unique<TFileDirectory>(fs->tFileDirectory());
0193         } else {
0194           tfd = std::make_unique<TFileDirectory>(fs->mkdir(newDir));
0195           directoryString = newDir;
0196         }
0197       } else {
0198         if (newDir.length() != 0) {
0199           if (basedir.length() != 0)
0200             directoryString = basedir + "/" + newDir;
0201           else
0202             directoryString = newDir;
0203         } else
0204           directoryString = basedir;
0205         theDbe = edm::Service<DQMStore>().operator->();
0206       }
0207     }
0208     // Generalization of Histogram Booking; allows switch between TFileService and DQMStore
0209     template <typename T>
0210     TH1* make(const char* name, const char* title, int nBinX, double minBinX, double maxBinX);
0211     template <typename T>
0212     TH1* make(const char* name, const char* title, int nBinX, double* xBins);  //variable bin size in x for profile histo
0213     template <typename T>
0214     TH1* make(const char* name,
0215               const char* title,
0216               int nBinX,
0217               double minBinX,
0218               double maxBinX,
0219               int nBinY,
0220               double minBinY,
0221               double maxBinY);
0222     template <typename T>
0223     TH1* make(const char* name,
0224               const char* title,
0225               int nBinX,
0226               double minBinX,
0227               double maxBinX,
0228               double minBinY,
0229               double maxBinY);  // at present not used
0230 
0231     std::unique_ptr<TFileDirectory> tfd;
0232     std::string directoryString;
0233     const bool dqmMode;
0234     DQMStore* theDbe;
0235   };
0236 
0237   //
0238   // ------------- private member function -------------
0239   //
0240   void analyze(const edm::Event&, const edm::EventSetup&) override;
0241   void endJob() override;
0242 
0243   virtual void checkBookHists(const edm::EventSetup& setup);
0244 
0245   void bookGlobalHists(DirectoryWrapper& tfd);
0246   void bookDirHists(DirectoryWrapper& tfd, const Alignable& ali, const TrackerTopology* tTopo);
0247   void bookHists(
0248       DirectoryWrapper& tfd, const Alignable& ali, const TrackerTopology* tTopo, align::StructureType type, int i);
0249 
0250   void prepareSummaryHists(DirectoryWrapper& tfd,
0251                            const Alignable& ali,
0252                            std::vector<TrackerOfflineValidation::SummaryContainer>& vLevelProfiles);
0253   void collateSummaryHists();
0254 
0255   void setUpTreeMembers(const std::map<int, TrackerOfflineValidation::ModuleHistos>& moduleHist_,
0256                         const TrackerGeometry& tkgeom,
0257                         const TrackerTopology* tTopo);
0258   void fillTree(TTree& tree,
0259                 TkOffTreeVariables& treeMem,
0260                 const std::map<int, TrackerOfflineValidation::ModuleHistos>& moduleHist_);
0261 
0262   TrackerOfflineValidation::SummaryContainer bookSummaryHists(DirectoryWrapper& tfd,
0263                                                               const Alignable& ali,
0264                                                               align::StructureType type,
0265                                                               int i);
0266 
0267   ModuleHistos& getHistStructFromMap(const DetId& detid);
0268 
0269   bool isBarrel(uint32_t subDetId);
0270   bool isEndCap(uint32_t subDetId);
0271   bool isPixel(uint32_t subDetId);
0272   bool isDetOrDetUnit(align::StructureType type);
0273 
0274   TH1* bookTH1F(bool isTransient,
0275                 DirectoryWrapper& tfd,
0276                 const char* histName,
0277                 const char* histTitle,
0278                 int nBinsX,
0279                 double lowX,
0280                 double highX);
0281 
0282   TProfile* bookTProfile(bool isTransient,
0283                          DirectoryWrapper& tfd,
0284                          const char* histName,
0285                          const char* histTitle,
0286                          int nBinsX,
0287                          double lowX,
0288                          double highX);
0289 
0290   TProfile* bookTProfile(bool isTransient,
0291                          DirectoryWrapper& tfd,
0292                          const char* histName,
0293                          const char* histTitle,
0294                          int nBinsX,
0295                          double lowX,
0296                          double highX,
0297                          double lowY,
0298                          double highY);
0299 
0300   void getBinning(uint32_t subDetId,
0301                   TrackerOfflineValidation::HistogramType residualtype,
0302                   int& nBinsX,
0303                   double& lowerBoundX,
0304                   double& upperBoundX);
0305 
0306   void summarizeBinInContainer(int bin, SummaryContainer& targetContainer, SummaryContainer& sourceContainer);
0307 
0308   void summarizeBinInContainer(int bin,
0309                                uint32_t subDetId,
0310                                SummaryContainer& targetContainer,
0311                                ModuleHistos& sourceContainer);
0312 
0313   void setSummaryBin(int bin, TH1* targetHist, TH1* sourceHist);
0314 
0315   float Fwhm(const TH1* hist) const;
0316   std::pair<float, float> fitResiduals(TH1* hist) const;  //, float meantmp, float rmstmp);
0317   float getMedian(const TH1* hist) const;
0318 
0319   // From MillePedeAlignmentMonitor: Get Index for Arbitary vector<class> by name
0320   template <class OBJECT_TYPE>
0321   int GetIndex(const std::vector<OBJECT_TYPE*>& vec, const TString& name);
0322 
0323   // ---------- member data ---------------------------
0324 
0325   const edm::ParameterSet parSet_;
0326   edm::ESHandle<TrackerGeometry> tkGeom_;
0327   const TrackerGeometry* bareTkGeomPtr_;  // ugly hack to book hists only once, but check
0328   std::unique_ptr<AlignableTracker> alignableTracker_;
0329 
0330   // parameters from cfg to steer
0331   const int compressionSettings_;
0332   const bool lCoorHistOn_;
0333   const bool moduleLevelHistsTransient_;
0334   const bool moduleLevelProfiles_;
0335   const bool stripYResiduals_;
0336   const bool useFwhm_;
0337   const bool useFit_;
0338   const bool useOverflowForRMS_;
0339   const bool dqmMode_;
0340   const std::string moduleDirectory_;
0341 
0342   const int chargeCut_;
0343 
0344   // a vector to keep track which pointers should be deleted at the very end
0345   std::vector<TH1*> vDeleteObjects_;
0346 
0347   std::vector<TH1*> vTrackHistos_;
0348   std::vector<TH1*> vTrackProfiles_;
0349   std::vector<TH1*> vTrack2DHistos_;
0350 
0351   std::map<int, TrackerOfflineValidation::ModuleHistos> mPxbResiduals_;
0352   std::map<int, TrackerOfflineValidation::ModuleHistos> mPxeResiduals_;
0353   std::map<int, TrackerOfflineValidation::ModuleHistos> mTibResiduals_;
0354   std::map<int, TrackerOfflineValidation::ModuleHistos> mTidResiduals_;
0355   std::map<int, TrackerOfflineValidation::ModuleHistos> mTobResiduals_;
0356   std::map<int, TrackerOfflineValidation::ModuleHistos> mTecResiduals_;
0357 
0358   std::map<int, TkOffTreeVariables> mTreeMembers_;
0359 
0360   //There are two types of summary histograms.  The first contains, for each component, a bin per subcomponent.
0361   //These are set up through summarizeBinInContainer().  The second type is just the sum of the lower level histograms.
0362   //Prepare the filling of both types at the beginning, when the tracker topology is available through the eventsetup.
0363   //They are not actually filled until the end, but at that time eventsetup is no longer accessible.
0364 
0365   //To fill the summary hists, store the arguments of setSummaryBin()
0366   //At the end, call the function
0367   std::vector<std::tuple<int, TH1*, TH1*> > summaryBins_;
0368   //To fill the sum hists, just store pairs of TH1*.  At the end, first->Add(second).
0369   std::vector<std::pair<TH1*, TH1*> > sumHistStructure_;
0370   //sum hists are fit at the end using fitResiduals()
0371   std::vector<TH1*> toFit_;
0372 
0373   unsigned long long nTracks_;
0374   const unsigned long long maxTracks_;
0375 
0376   TrackerValidationVariables avalidator_;
0377 };
0378 
0379 //
0380 // constants, enums and typedefs
0381 //
0382 
0383 //
0384 // static data member definitions
0385 //
0386 
0387 template <class OBJECT_TYPE>
0388 int TrackerOfflineValidation::GetIndex(const std::vector<OBJECT_TYPE*>& vec, const TString& name) {
0389   int result = 0;
0390   for (typename std::vector<OBJECT_TYPE*>::const_iterator iter = vec.begin(), iterEnd = vec.end(); iter != iterEnd;
0391        ++iter, ++result) {
0392     if (*iter && (*iter)->GetName() == name)
0393       return result;
0394   }
0395   edm::LogError("Alignment") << "@SUB=TrackerOfflineValidation::GetIndex"
0396                              << " could not find " << name;
0397   return -1;
0398 }
0399 
0400 template <>
0401 TH1* TrackerOfflineValidation::DirectoryWrapper::make<TH1F>(
0402     const char* name, const char* title, int nBinX, double minBinX, double maxBinX) {
0403   if (dqmMode) {
0404     theDbe->setCurrentFolder(directoryString);
0405     return theDbe->book1D(name, title, nBinX, minBinX, maxBinX)->getTH1();
0406   } else {
0407     return tfd->make<TH1F>(name, title, nBinX, minBinX, maxBinX);
0408   }
0409 }
0410 
0411 template <>
0412 TH1* TrackerOfflineValidation::DirectoryWrapper::make<TProfile>(const char* name,
0413                                                                 const char* title,
0414                                                                 int nBinX,
0415                                                                 double* xBins) {
0416   if (dqmMode) {
0417     theDbe->setCurrentFolder(directoryString);
0418     //DQM profile requires y-bins for construction... using TProfile creator by hand...
0419     TProfile* tmpProfile = new TProfile(name, title, nBinX, xBins);
0420     tmpProfile->SetDirectory(nullptr);
0421     return theDbe->bookProfile(name, tmpProfile)->getTH1();
0422   } else {
0423     return tfd->make<TProfile>(name, title, nBinX, xBins);
0424   }
0425 }
0426 
0427 template <>
0428 TH1* TrackerOfflineValidation::DirectoryWrapper::make<TProfile>(
0429     const char* name, const char* title, int nBinX, double minBinX, double maxBinX) {
0430   if (dqmMode) {
0431     theDbe->setCurrentFolder(directoryString);
0432     //DQM profile requires y-bins for construction... using TProfile creator by hand...
0433     TProfile* tmpProfile = new TProfile(name, title, nBinX, minBinX, maxBinX);
0434     tmpProfile->SetDirectory(nullptr);
0435     return theDbe->bookProfile(name, tmpProfile)->getTH1();
0436   } else {
0437     return tfd->make<TProfile>(name, title, nBinX, minBinX, maxBinX);
0438   }
0439 }
0440 
0441 template <>
0442 TH1* TrackerOfflineValidation::DirectoryWrapper::make<TProfile>(
0443     const char* name, const char* title, int nbinX, double minX, double maxX, double minY, double maxY) {
0444   if (dqmMode) {
0445     theDbe->setCurrentFolder(directoryString);
0446     int dummy(0);  // DQMProfile wants Y channels... does not use them!
0447     return (theDbe->bookProfile(name, title, nbinX, minX, maxX, dummy, minY, maxY)->getTH1());
0448   } else {
0449     return tfd->make<TProfile>(name, title, nbinX, minX, maxX, minY, maxY);
0450   }
0451 }
0452 
0453 template <>
0454 TH1* TrackerOfflineValidation::DirectoryWrapper::make<TH2F>(const char* name,
0455                                                             const char* title,
0456                                                             int nBinX,
0457                                                             double minBinX,
0458                                                             double maxBinX,
0459                                                             int nBinY,
0460                                                             double minBinY,
0461                                                             double maxBinY) {
0462   if (dqmMode) {
0463     theDbe->setCurrentFolder(directoryString);
0464     return theDbe->book2D(name, title, nBinX, minBinX, maxBinX, nBinY, minBinY, maxBinY)->getTH1();
0465   } else {
0466     return tfd->make<TH2F>(name, title, nBinX, minBinX, maxBinX, nBinY, minBinY, maxBinY);
0467   }
0468 }
0469 
0470 //
0471 // constructors and destructor
0472 //
0473 TrackerOfflineValidation::TrackerOfflineValidation(const edm::ParameterSet& iConfig)
0474     : geomToken_(esConsumes()),
0475       topoToken_(esConsumes()),
0476       parSet_(iConfig),
0477       bareTkGeomPtr_(nullptr),
0478       compressionSettings_(parSet_.getUntrackedParameter<int>("compressionSettings", -1)),
0479       lCoorHistOn_(parSet_.getParameter<bool>("localCoorHistosOn")),
0480       moduleLevelHistsTransient_(parSet_.getParameter<bool>("moduleLevelHistsTransient")),
0481       moduleLevelProfiles_(parSet_.getParameter<bool>("moduleLevelProfiles")),
0482       stripYResiduals_(parSet_.getParameter<bool>("stripYResiduals")),
0483       useFwhm_(parSet_.getParameter<bool>("useFwhm")),
0484       useFit_(parSet_.getParameter<bool>("useFit")),
0485       useOverflowForRMS_(parSet_.getParameter<bool>("useOverflowForRMS")),
0486       dqmMode_(parSet_.getParameter<bool>("useInDqmMode")),
0487       moduleDirectory_(parSet_.getParameter<std::string>("moduleDirectoryInOutput")),
0488       chargeCut_(parSet_.getParameter<int>("chargeCut")),
0489       nTracks_(0),
0490       maxTracks_(parSet_.getParameter<unsigned long long>("maxTracks")),
0491       avalidator_(iConfig, consumesCollector()) {
0492   usesResource(TFileService::kSharedResource);
0493 }
0494 
0495 TrackerOfflineValidation::~TrackerOfflineValidation() {
0496   // do anything here that needs to be done at desctruction time
0497   // (e.g. close files, deallocate resources etc.)
0498   for (std::vector<TH1*>::const_iterator it = vDeleteObjects_.begin(), itEnd = vDeleteObjects_.end(); it != itEnd; ++it)
0499     delete *it;
0500 }
0501 
0502 //
0503 // member functions
0504 //
0505 
0506 // ------------ method called once each job just before starting event loop  ------------
0507 void TrackerOfflineValidation::checkBookHists(const edm::EventSetup& es) {
0508   tkGeom_ = es.getHandle(geomToken_);
0509   const TrackerGeometry* newBareTkGeomPtr = &(*tkGeom_);
0510 
0511   if (newBareTkGeomPtr == bareTkGeomPtr_)
0512     return;  // already booked hists, nothing changed
0513 
0514   if (!bareTkGeomPtr_) {  // pointer not yet set: called the first time => book hists
0515 
0516     //Retrieve tracker topology from geometry
0517     const TrackerTopology* const tTopo = &es.getData(topoToken_);
0518 
0519     // construct alignable tracker to get access to alignable hierarchy
0520     alignableTracker_ = std::make_unique<AlignableTracker>(&(*tkGeom_), tTopo);
0521 
0522     edm::LogInfo("TrackerOfflineValidation")
0523         << "There are " << newBareTkGeomPtr->detIds().size() << " dets in the Geometry record.\n"
0524         << "Out of these " << newBareTkGeomPtr->detUnitIds().size() << " are detUnits";
0525 
0526     // Book histograms for global track quantities
0527     std::string globDir("GlobalTrackVariables");
0528     DirectoryWrapper trackglobal(globDir, moduleDirectory_, dqmMode_);
0529     this->bookGlobalHists(trackglobal);
0530 
0531     // recursively book histograms on lowest level
0532     DirectoryWrapper tfdw("", moduleDirectory_, dqmMode_);
0533     this->bookDirHists(tfdw, *alignableTracker_, tTopo);
0534     // and prepare the higher level histograms
0535     std::vector<TrackerOfflineValidation::SummaryContainer> vTrackerprofiles;
0536     this->prepareSummaryHists(tfdw, *alignableTracker_, vTrackerprofiles);
0537 
0538     setUpTreeMembers(mPxbResiduals_, *newBareTkGeomPtr, tTopo);
0539     setUpTreeMembers(mPxeResiduals_, *newBareTkGeomPtr, tTopo);
0540     setUpTreeMembers(mTibResiduals_, *newBareTkGeomPtr, tTopo);
0541     setUpTreeMembers(mTidResiduals_, *newBareTkGeomPtr, tTopo);
0542     setUpTreeMembers(mTobResiduals_, *newBareTkGeomPtr, tTopo);
0543     setUpTreeMembers(mTecResiduals_, *newBareTkGeomPtr, tTopo);
0544   } else {  // histograms booked, but changed TrackerGeometry?
0545     edm::LogWarning("GeometryChange") << "@SUB=checkBookHists"
0546                                       << "TrackerGeometry changed, but will not re-book hists!";
0547   }
0548   bareTkGeomPtr_ = newBareTkGeomPtr;
0549 }
0550 
0551 void TrackerOfflineValidation::bookGlobalHists(DirectoryWrapper& tfd) {
0552   vTrackHistos_.push_back(tfd.make<TH1F>("h_tracketa", "Track #eta;#eta_{Track};Number of Tracks", 90, -3., 3.));
0553   vTrackHistos_.push_back(tfd.make<TH1F>("h_trackphi", "Track #phi;#phi_{Track};Number of Tracks", 90, -3.15, 3.15));
0554   vTrackHistos_.push_back(tfd.make<TH1F>(
0555       "h_trackNumberOfValidHits", "Track # of valid hits;# of valid hits _{Track};Number of Tracks", 40, 0., 40.));
0556   vTrackHistos_.push_back(tfd.make<TH1F>(
0557       "h_trackNumberOfLostHits", "Track # of lost hits;# of lost hits _{Track};Number of Tracks", 10, 0., 10.));
0558   vTrackHistos_.push_back(
0559       tfd.make<TH1F>("h_curvature", "Curvature #kappa;#kappa_{Track};Number of Tracks", 100, -.05, .05));
0560   vTrackHistos_.push_back(tfd.make<TH1F>(
0561       "h_curvature_pos", "Curvature |#kappa| Positive Tracks;|#kappa_{pos Track}|;Number of Tracks", 100, .0, .05));
0562   vTrackHistos_.push_back(tfd.make<TH1F>(
0563       "h_curvature_neg", "Curvature |#kappa| Negative Tracks;|#kappa_{neg Track}|;Number of Tracks", 100, .0, .05));
0564   vTrackHistos_.push_back(
0565       tfd.make<TH1F>("h_diff_curvature",
0566                      "Curvature |#kappa| Tracks Difference;|#kappa_{Track}|;# Pos Tracks - # Neg Tracks",
0567                      100,
0568                      .0,
0569                      .05));
0570   vTrackHistos_.push_back(tfd.make<TH1F>("h_chi2", "#chi^{2};#chi^{2}_{Track};Number of Tracks", 500, -0.01, 500.));
0571   vTrackHistos_.push_back(
0572       tfd.make<TH1F>("h_chi2Prob", "#chi^{2} probability;#chi^{2}prob_{Track};Number of Tracks", 100, 0.0, 1.));
0573   vTrackHistos_.push_back(
0574       tfd.make<TH1F>("h_normchi2", "#chi^{2}/ndof;#chi^{2}/ndof;Number of Tracks", 100, -0.01, 10.));
0575   vTrackHistos_.push_back(tfd.make<TH1F>("h_pt", "p_{T}^{track};p_{T}^{track} [GeV];Number of Tracks", 250, 0., 250));
0576   vTrackHistos_.push_back(tfd.make<TH1F>(
0577       "h_ptResolution", "#delta_{p_{T}}/p_{T}^{track};#delta_{p_{T}}/p_{T}^{track};Number of Tracks", 100, 0., 0.5));
0578   vTrackProfiles_.push_back(tfd.make<TProfile>(
0579       "p_d0_vs_phi", "Transverse Impact Parameter vs. #phi;#phi_{Track};#LT d_{0} #GT [cm]", 100, -3.15, 3.15));
0580   vTrackProfiles_.push_back(tfd.make<TProfile>(
0581       "p_dz_vs_phi", "Longitudinal Impact Parameter vs. #phi;#phi_{Track};#LT d_{z} #GT [cm]", 100, -3.15, 3.15));
0582   vTrackProfiles_.push_back(tfd.make<TProfile>(
0583       "p_d0_vs_eta", "Transverse Impact Parameter vs. #eta;#eta_{Track};#LT d_{0} #GT [cm]", 100, -3.15, 3.15));
0584   vTrackProfiles_.push_back(tfd.make<TProfile>(
0585       "p_dz_vs_eta", "Longitudinal Impact Parameter vs. #eta;#eta_{Track};#LT d_{z} #GT [cm]", 100, -3.15, 3.15));
0586   vTrackProfiles_.push_back(
0587       tfd.make<TProfile>("p_chi2_vs_phi", "#chi^{2} vs. #phi;#phi_{Track};#LT #chi^{2} #GT", 100, -3.15, 3.15));
0588   vTrackProfiles_.push_back(tfd.make<TProfile>(
0589       "p_chi2Prob_vs_phi", "#chi^{2} probablility vs. #phi;#phi_{Track};#LT #chi^{2} probability#GT", 100, -3.15, 3.15));
0590   vTrackProfiles_.push_back(tfd.make<TProfile>(
0591       "p_chi2Prob_vs_d0", "#chi^{2} probablility vs. |d_{0}|;|d_{0}|[cm];#LT #chi^{2} probability#GT", 100, 0, 80));
0592   vTrackProfiles_.push_back(tfd.make<TProfile>(
0593       "p_normchi2_vs_phi", "#chi^{2}/ndof vs. #phi;#phi_{Track};#LT #chi^{2}/ndof #GT", 100, -3.15, 3.15));
0594   vTrackProfiles_.push_back(
0595       tfd.make<TProfile>("p_chi2_vs_eta", "#chi^{2} vs. #eta;#eta_{Track};#LT #chi^{2} #GT", 100, -3.15, 3.15));
0596   //variable binning for chi2/ndof vs. pT
0597   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.};
0598   vTrackProfiles_.push_back(tfd.make<TProfile>(
0599       "p_normchi2_vs_pt", "norm #chi^{2} vs. p_{T}_{Track}; p_{T}_{Track};#LT #chi^{2}/ndof #GT", 18, xBins));
0600 
0601   vTrackProfiles_.push_back(
0602       tfd.make<TProfile>("p_normchi2_vs_p", "#chi^{2}/ndof vs. p_{Track};p_{Track};#LT #chi^{2}/ndof #GT", 18, xBins));
0603   vTrackProfiles_.push_back(tfd.make<TProfile>(
0604       "p_chi2Prob_vs_eta", "#chi^{2} probability vs. #eta;#eta_{Track};#LT #chi^{2} probability #GT", 100, -3.15, 3.15));
0605   vTrackProfiles_.push_back(tfd.make<TProfile>(
0606       "p_normchi2_vs_eta", "#chi^{2}/ndof vs. #eta;#eta_{Track};#LT #chi^{2}/ndof #GT", 100, -3.15, 3.15));
0607   vTrackProfiles_.push_back(
0608       tfd.make<TProfile>("p_kappa_vs_phi", "#kappa vs. #phi;#phi_{Track};#kappa", 100, -3.15, 3.15));
0609   vTrackProfiles_.push_back(
0610       tfd.make<TProfile>("p_kappa_vs_eta", "#kappa vs. #eta;#eta_{Track};#kappa", 100, -3.15, 3.15));
0611   vTrackProfiles_.push_back(tfd.make<TProfile>("p_ptResolution_vs_phi",
0612                                                "#delta_{p_{T}}/p_{T}^{track};#phi^{track};#delta_{p_{T}}/p_{T}^{track}",
0613                                                100,
0614                                                -3.15,
0615                                                3.15));
0616   vTrackProfiles_.push_back(tfd.make<TProfile>("p_ptResolution_vs_eta",
0617                                                "#delta_{p_{T}}/p_{T}^{track};#eta^{track};#delta_{p_{T}}/p_{T}^{track}",
0618                                                100,
0619                                                -3.15,
0620                                                3.15));
0621 
0622   vTrack2DHistos_.push_back(tfd.make<TH2F>(
0623       "h2_d0_vs_phi", "Transverse Impact Parameter vs. #phi;#phi_{Track};d_{0} [cm]", 100, -3.15, 3.15, 100, -1., 1.));
0624   vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_dz_vs_phi",
0625                                            "Longitudinal Impact Parameter vs. #phi;#phi_{Track};d_{z} [cm]",
0626                                            100,
0627                                            -3.15,
0628                                            3.15,
0629                                            100,
0630                                            -100.,
0631                                            100.));
0632   vTrack2DHistos_.push_back(tfd.make<TH2F>(
0633       "h2_d0_vs_eta", "Transverse Impact Parameter vs. #eta;#eta_{Track};d_{0} [cm]", 100, -3.15, 3.15, 100, -1., 1.));
0634   vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_dz_vs_eta",
0635                                            "Longitudinal Impact Parameter vs. #eta;#eta_{Track};d_{z} [cm]",
0636                                            100,
0637                                            -3.15,
0638                                            3.15,
0639                                            100,
0640                                            -100.,
0641                                            100.));
0642   vTrack2DHistos_.push_back(
0643       tfd.make<TH2F>("h2_chi2_vs_phi", "#chi^{2} vs. #phi;#phi_{Track};#chi^{2}", 100, -3.15, 3.15, 500, 0., 500.));
0644   vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_chi2Prob_vs_phi",
0645                                            "#chi^{2} probability vs. #phi;#phi_{Track};#chi^{2} probability",
0646                                            100,
0647                                            -3.15,
0648                                            3.15,
0649                                            100,
0650                                            0.,
0651                                            1.));
0652   vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_chi2Prob_vs_d0",
0653                                            "#chi^{2} probability vs. |d_{0}|;|d_{0}| [cm];#chi^{2} probability",
0654                                            100,
0655                                            0,
0656                                            80,
0657                                            100,
0658                                            0.,
0659                                            1.));
0660   vTrack2DHistos_.push_back(tfd.make<TH2F>(
0661       "h2_normchi2_vs_phi", "#chi^{2}/ndof vs. #phi;#phi_{Track};#chi^{2}/ndof", 100, -3.15, 3.15, 100, 0., 10.));
0662   vTrack2DHistos_.push_back(
0663       tfd.make<TH2F>("h2_chi2_vs_eta", "#chi^{2} vs. #eta;#eta_{Track};#chi^{2}", 100, -3.15, 3.15, 500, 0., 500.));
0664   vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_chi2Prob_vs_eta",
0665                                            "#chi^{2} probaility vs. #eta;#eta_{Track};#chi^{2} probability",
0666                                            100,
0667                                            -3.15,
0668                                            3.15,
0669                                            100,
0670                                            0.,
0671                                            1.));
0672   vTrack2DHistos_.push_back(tfd.make<TH2F>(
0673       "h2_normchi2_vs_eta", "#chi^{2}/ndof vs. #eta;#eta_{Track};#chi^{2}/ndof", 100, -3.15, 3.15, 100, 0., 10.));
0674   vTrack2DHistos_.push_back(
0675       tfd.make<TH2F>("h2_kappa_vs_phi", "#kappa vs. #phi;#phi_{Track};#kappa", 100, -3.15, 3.15, 100, .0, .05));
0676   vTrack2DHistos_.push_back(
0677       tfd.make<TH2F>("h2_kappa_vs_eta", "#kappa vs. #eta;#eta_{Track};#kappa", 100, -3.15, 3.15, 100, .0, .05));
0678   vTrack2DHistos_.push_back(tfd.make<TH2F>(
0679       "h2_normchi2_vs_kappa", "#kappa vs. #chi^{2}/ndof;#chi^{2}/ndof;#kappa", 100, 0., 10, 100, -.03, .03));
0680 
0681   /****************** Definition of 2-D Histos of ResX vs momenta ****************************/
0682   vTrack2DHistos_.push_back(tfd.make<TH2F>(
0683       "p_vs_resXprime_pixB", "res_{x'} vs momentum in BPix;p [GeV]; res_{x'}", 15, 0., 15., 200, -0.1, 0.1));
0684   vTrack2DHistos_.push_back(tfd.make<TH2F>(
0685       "p_vs_resXprime_pixE", "res_{x'} vs momentum in FPix;p [GeV]; res_{x'}", 15, 0., 15., 200, -0.1, 0.1));
0686   vTrack2DHistos_.push_back(tfd.make<TH2F>(
0687       "p_vs_resXprime_TIB", "res_{x'} vs momentum in TIB;p [GeV]; res_{x'}", 15, 0., 15., 200, -0.1, 0.1));
0688   vTrack2DHistos_.push_back(tfd.make<TH2F>(
0689       "p_vs_resXprime_TID", "res_{x'} vs momentum in TID;p [GeV]; res_{x'}", 15, 0., 15., 200, -0.1, 0.1));
0690   vTrack2DHistos_.push_back(tfd.make<TH2F>(
0691       "p_vs_resXprime_TOB", "res_{x'} vs momentum in TOB;p [GeV]; res_{x'}", 15, 0., 15., 200, -0.1, 0.1));
0692   vTrack2DHistos_.push_back(tfd.make<TH2F>(
0693       "p_vs_resXprime_TEC", "res_{x'} vs momentum in TEC;p [GeV]; res_{x'}", 15, 0., 15., 200, -0.1, 0.1));
0694 
0695   /****************** Definition of 2-D Histos of ResY vs momenta ****************************/
0696   vTrack2DHistos_.push_back(tfd.make<TH2F>(
0697       "p_vs_resYprime_pixB", "res_{y'} vs momentum in BPix;p [GeV]; res_{y'}", 15, 0., 15., 200, -0.1, 0.1));
0698   vTrack2DHistos_.push_back(tfd.make<TH2F>(
0699       "p_vs_resYprime_pixE", "res_{y'} vs momentum in FPix;p [GeV]; res_{y'}", 15, 0., 15., 200, -0.1, 0.1));
0700 }
0701 
0702 void TrackerOfflineValidation::bookDirHists(DirectoryWrapper& tfd, const Alignable& ali, const TrackerTopology* tTopo) {
0703   const auto& alivec = ali.components();
0704   for (int i = 0, iEnd = ali.components().size(); i < iEnd; ++i) {
0705     std::string structurename = alignableTracker_->objectIdProvider().idToString((alivec)[i]->alignableObjectId());
0706     LogDebug("TrackerOfflineValidation") << "StructureName = " << structurename;
0707     std::stringstream dirname;
0708     dirname << structurename;
0709     // add no suffix counter to Strip and Pixel, just aesthetics
0710     if (structurename != "Strip" && structurename != "Pixel")
0711       dirname << "_" << i + 1;
0712 
0713     if (structurename.find("Endcap", 0) != std::string::npos) {
0714       DirectoryWrapper f(tfd, dirname.str(), moduleDirectory_, dqmMode_);
0715       bookHists(f, *(alivec)[i], tTopo, ali.alignableObjectId(), i);
0716       bookDirHists(f, *(alivec)[i], tTopo);
0717     } else if (!(this->isDetOrDetUnit((alivec)[i]->alignableObjectId())) || alivec[i]->components().size() > 1) {
0718       DirectoryWrapper f(tfd, dirname.str(), moduleDirectory_, dqmMode_);
0719       bookHists(tfd, *(alivec)[i], tTopo, ali.alignableObjectId(), i);
0720       bookDirHists(f, *(alivec)[i], tTopo);
0721     } else {
0722       bookHists(tfd, *(alivec)[i], tTopo, ali.alignableObjectId(), i);
0723     }
0724   }
0725 }
0726 
0727 void TrackerOfflineValidation::bookHists(
0728     DirectoryWrapper& tfd, const Alignable& ali, const TrackerTopology* tTopo, align::StructureType type, int i) {
0729   TrackerAlignableId aliid;
0730   const DetId id = ali.id();
0731 
0732   // comparing subdetandlayer to subdetIds gives a warning at compile time
0733   // -> subdetandlayer could also be pair<uint,uint> but this has to be adapted
0734   // in AlignableObjId
0735   std::pair<int, int> subdetandlayer = aliid.typeAndLayerFromDetId(id, tTopo);
0736 
0737   align::StructureType subtype = align::invalid;
0738 
0739   // are we on or just above det, detunit level respectively?
0740   if (type == align::AlignableDetUnit)
0741     subtype = type;
0742   else if (this->isDetOrDetUnit(ali.alignableObjectId()))
0743     subtype = ali.alignableObjectId();
0744 
0745   // construct histogram title and name
0746   std::stringstream histoname, histotitle, normhistoname, normhistotitle, yhistoname, yhistotitle, xprimehistoname,
0747       xprimehistotitle, normxprimehistoname, normxprimehistotitle, yprimehistoname, yprimehistotitle,
0748       normyprimehistoname, normyprimehistotitle, localxname, localxtitle, localyname, localytitle, resxvsxprofilename,
0749       resxvsxprofiletitle, resyvsxprofilename, resyvsxprofiletitle, resxvsyprofilename, resxvsyprofiletitle,
0750       resyvsyprofilename, resyvsyprofiletitle;
0751 
0752   std::string wheel_or_layer;
0753 
0754   if (this->isEndCap(static_cast<uint32_t>(subdetandlayer.first)))
0755     wheel_or_layer = "_wheel_";
0756   else if (this->isBarrel(static_cast<uint32_t>(subdetandlayer.first)))
0757     wheel_or_layer = "_layer_";
0758   else
0759     edm::LogWarning("TrackerOfflineValidation") << "@SUB=TrackerOfflineValidation::bookHists"
0760                                                 << "Unknown subdetid: " << subdetandlayer.first;
0761 
0762   histoname << "h_residuals_subdet_" << subdetandlayer.first << wheel_or_layer << subdetandlayer.second << "_module_"
0763             << id.rawId();
0764   yhistoname << "h_y_residuals_subdet_" << subdetandlayer.first << wheel_or_layer << subdetandlayer.second << "_module_"
0765              << id.rawId();
0766   xprimehistoname << "h_xprime_residuals_subdet_" << subdetandlayer.first << wheel_or_layer << subdetandlayer.second
0767                   << "_module_" << id.rawId();
0768   yprimehistoname << "h_yprime_residuals_subdet_" << subdetandlayer.first << wheel_or_layer << subdetandlayer.second
0769                   << "_module_" << id.rawId();
0770   normhistoname << "h_normresiduals_subdet_" << subdetandlayer.first << wheel_or_layer << subdetandlayer.second
0771                 << "_module_" << id.rawId();
0772   normxprimehistoname << "h_normxprimeresiduals_subdet_" << subdetandlayer.first << wheel_or_layer
0773                       << subdetandlayer.second << "_module_" << id.rawId();
0774   normyprimehistoname << "h_normyprimeresiduals_subdet_" << subdetandlayer.first << wheel_or_layer
0775                       << subdetandlayer.second << "_module_" << id.rawId();
0776   histotitle << "X Residual for module " << id.rawId() << ";x_{tr} - x_{hit} [cm]";
0777   yhistotitle << "Y Residual for module " << id.rawId() << ";y_{tr} - y_{hit} [cm]";
0778   normhistotitle << "Normalized Residual for module " << id.rawId() << ";x_{tr} - x_{hit}/#sigma";
0779   xprimehistotitle << "X' Residual for module " << id.rawId() << ";(x_{tr} - x_{hit})' [cm]";
0780   normxprimehistotitle << "Normalized X' Residual for module " << id.rawId() << ";(x_{tr} - x_{hit})'/#sigma";
0781   yprimehistotitle << "Y' Residual for module " << id.rawId() << ";(y_{tr} - y_{hit})' [cm]";
0782   normyprimehistotitle << "Normalized Y' Residual for module " << id.rawId() << ";(y_{tr} - y_{hit})'/#sigma";
0783 
0784   if (moduleLevelProfiles_) {
0785     localxname << "h_localx_subdet_" << subdetandlayer.first << wheel_or_layer << subdetandlayer.second << "_module_"
0786                << id.rawId();
0787     localyname << "h_localy_subdet_" << subdetandlayer.first << wheel_or_layer << subdetandlayer.second << "_module_"
0788                << id.rawId();
0789     localxtitle << "u local for module " << id.rawId() << "; u_{tr,r}";
0790     localytitle << "v local for module " << id.rawId() << "; v_{tr,r}";
0791 
0792     resxvsxprofilename << "p_residuals_x_vs_x_subdet_" << subdetandlayer.first << wheel_or_layer
0793                        << subdetandlayer.second << "_module_" << id.rawId();
0794     resyvsxprofilename << "p_residuals_y_vs_x_subdet_" << subdetandlayer.first << wheel_or_layer
0795                        << subdetandlayer.second << "_module_" << id.rawId();
0796     resxvsyprofilename << "p_residuals_x_vs_y_subdet_" << subdetandlayer.first << wheel_or_layer
0797                        << subdetandlayer.second << "_module_" << id.rawId();
0798     resyvsyprofilename << "p_residuals_y_vs_y_subdet_" << subdetandlayer.first << wheel_or_layer
0799                        << subdetandlayer.second << "_module_" << id.rawId();
0800     resxvsxprofiletitle << "U Residual vs u for module " << id.rawId()
0801                         << "; u_{tr,r} ;(u_{tr} - u_{hit})/tan#alpha [cm]";
0802     resyvsxprofiletitle << "V Residual vs u for module " << id.rawId()
0803                         << "; u_{tr,r} ;(v_{tr} - v_{hit})/tan#beta  [cm]";
0804     resxvsyprofiletitle << "U Residual vs v for module " << id.rawId()
0805                         << "; v_{tr,r} ;(u_{tr} - u_{hit})/tan#alpha [cm]";
0806     resyvsyprofiletitle << "V Residual vs v for module " << id.rawId()
0807                         << "; v_{tr,r} ;(v_{tr} - v_{hit})/tan#beta  [cm]";
0808   }
0809 
0810   if (this->isDetOrDetUnit(subtype)) {
0811     ModuleHistos& histStruct = this->getHistStructFromMap(id);
0812     int nbins = 0;
0813     double xmin = 0., xmax = 0.;
0814     double ymin = -0.1, ymax = 0.1;
0815 
0816     // do not allow transient hists in DQM mode
0817     bool moduleLevelHistsTransient(moduleLevelHistsTransient_);
0818     if (dqmMode_)
0819       moduleLevelHistsTransient = false;
0820 
0821     // decide via cfg if hists in local coordinates should be booked
0822     if (lCoorHistOn_) {
0823       this->getBinning(id.subdetId(), XResidual, nbins, xmin, xmax);
0824       histStruct.ResHisto = this->bookTH1F(
0825           moduleLevelHistsTransient, tfd, histoname.str().c_str(), histotitle.str().c_str(), nbins, xmin, xmax);
0826       this->getBinning(id.subdetId(), NormXResidual, nbins, xmin, xmax);
0827       histStruct.NormResHisto = this->bookTH1F(
0828           moduleLevelHistsTransient, tfd, normhistoname.str().c_str(), normhistotitle.str().c_str(), nbins, xmin, xmax);
0829     }
0830     this->getBinning(id.subdetId(), XprimeResidual, nbins, xmin, xmax);
0831     histStruct.ResXprimeHisto = this->bookTH1F(moduleLevelHistsTransient,
0832                                                tfd,
0833                                                xprimehistoname.str().c_str(),
0834                                                xprimehistotitle.str().c_str(),
0835                                                nbins,
0836                                                xmin,
0837                                                xmax);
0838     this->getBinning(id.subdetId(), NormXprimeResidual, nbins, xmin, xmax);
0839     histStruct.NormResXprimeHisto = this->bookTH1F(moduleLevelHistsTransient,
0840                                                    tfd,
0841                                                    normxprimehistoname.str().c_str(),
0842                                                    normxprimehistotitle.str().c_str(),
0843                                                    nbins,
0844                                                    xmin,
0845                                                    xmax);
0846 
0847     if (moduleLevelProfiles_) {
0848       this->getBinning(id.subdetId(), XResidualProfile, nbins, xmin, xmax);
0849 
0850       histStruct.LocalX = this->bookTH1F(
0851           moduleLevelHistsTransient, tfd, localxname.str().c_str(), localxtitle.str().c_str(), nbins, xmin, xmax);
0852       histStruct.LocalY = this->bookTH1F(
0853           moduleLevelHistsTransient, tfd, localyname.str().c_str(), localytitle.str().c_str(), nbins, xmin, xmax);
0854       histStruct.ResXvsXProfile = this->bookTProfile(moduleLevelHistsTransient,
0855                                                      tfd,
0856                                                      resxvsxprofilename.str().c_str(),
0857                                                      resxvsxprofiletitle.str().c_str(),
0858                                                      nbins,
0859                                                      xmin,
0860                                                      xmax,
0861                                                      ymin,
0862                                                      ymax);
0863       histStruct.ResXvsXProfile->Sumw2();  // to be filled with weights, so uncertainties need sum of square of weights
0864       histStruct.ResXvsYProfile = this->bookTProfile(moduleLevelHistsTransient,
0865                                                      tfd,
0866                                                      resxvsyprofilename.str().c_str(),
0867                                                      resxvsyprofiletitle.str().c_str(),
0868                                                      nbins,
0869                                                      xmin,
0870                                                      xmax,
0871                                                      ymin,
0872                                                      ymax);
0873       histStruct.ResXvsYProfile->Sumw2();  // to be filled with weights, so uncertainties need sum of square of weights
0874     }
0875 
0876     if (this->isPixel(subdetandlayer.first) || stripYResiduals_) {
0877       this->getBinning(id.subdetId(), YprimeResidual, nbins, xmin, xmax);
0878       histStruct.ResYprimeHisto = this->bookTH1F(moduleLevelHistsTransient,
0879                                                  tfd,
0880                                                  yprimehistoname.str().c_str(),
0881                                                  yprimehistotitle.str().c_str(),
0882                                                  nbins,
0883                                                  xmin,
0884                                                  xmax);
0885       if (lCoorHistOn_) {  // un-primed y-residual
0886         this->getBinning(id.subdetId(), YResidual, nbins, xmin, xmax);
0887         histStruct.ResYHisto = this->bookTH1F(
0888             moduleLevelHistsTransient, tfd, yhistoname.str().c_str(), yhistotitle.str().c_str(), nbins, xmin, xmax);
0889       }
0890       this->getBinning(id.subdetId(), NormYprimeResidual, nbins, xmin, xmax);
0891       histStruct.NormResYprimeHisto = this->bookTH1F(moduleLevelHistsTransient,
0892                                                      tfd,
0893                                                      normyprimehistoname.str().c_str(),
0894                                                      normyprimehistotitle.str().c_str(),
0895                                                      nbins,
0896                                                      xmin,
0897                                                      xmax);
0898       // Here we could add un-primed normalised y-residuals if(lCoorHistOn_)...
0899       if (moduleLevelProfiles_) {
0900         this->getBinning(id.subdetId(), YResidualProfile, nbins, xmin, xmax);
0901 
0902         histStruct.ResYvsXProfile = this->bookTProfile(moduleLevelHistsTransient,
0903                                                        tfd,
0904                                                        resyvsxprofilename.str().c_str(),
0905                                                        resyvsxprofiletitle.str().c_str(),
0906                                                        nbins,
0907                                                        xmin,
0908                                                        xmax,
0909                                                        ymin,
0910                                                        ymax);
0911         histStruct.ResYvsXProfile->Sumw2();  // to be filled with weights, so uncertainties need sum of square of weights
0912         histStruct.ResYvsYProfile = this->bookTProfile(moduleLevelHistsTransient,
0913                                                        tfd,
0914                                                        resyvsyprofilename.str().c_str(),
0915                                                        resyvsyprofiletitle.str().c_str(),
0916                                                        nbins,
0917                                                        xmin,
0918                                                        xmax,
0919                                                        ymin,
0920                                                        ymax);
0921         histStruct.ResYvsYProfile->Sumw2();  // to be filled with weights, so uncertainties need sum of square of weights
0922       }
0923     }
0924   }
0925 }
0926 
0927 TH1* TrackerOfflineValidation::bookTH1F(bool isTransient,
0928                                         DirectoryWrapper& tfd,
0929                                         const char* histName,
0930                                         const char* histTitle,
0931                                         int nBinsX,
0932                                         double lowX,
0933                                         double highX) {
0934   if (isTransient) {
0935     vDeleteObjects_.push_back(new TH1F(histName, histTitle, nBinsX, lowX, highX));
0936     return vDeleteObjects_.back();  // return last element of vector
0937   } else
0938     return tfd.make<TH1F>(histName, histTitle, nBinsX, lowX, highX);
0939 }
0940 
0941 TProfile* TrackerOfflineValidation::bookTProfile(bool isTransient,
0942                                                  DirectoryWrapper& tfd,
0943                                                  const char* histName,
0944                                                  const char* histTitle,
0945                                                  int nBinsX,
0946                                                  double lowX,
0947                                                  double highX) {
0948   if (isTransient) {
0949     TProfile* profile = new TProfile(histName, histTitle, nBinsX, lowX, highX);
0950     vDeleteObjects_.push_back(profile);
0951     return profile;
0952   } else
0953     return (TProfile*)tfd.make<TProfile>(histName, histTitle, nBinsX, lowX, highX);
0954 }
0955 
0956 TProfile* TrackerOfflineValidation::bookTProfile(bool isTransient,
0957                                                  DirectoryWrapper& tfd,
0958                                                  const char* histName,
0959                                                  const char* histTitle,
0960                                                  int nBinsX,
0961                                                  double lowX,
0962                                                  double highX,
0963                                                  double lowY,
0964                                                  double highY) {
0965   if (isTransient) {
0966     TProfile* profile = new TProfile(histName, histTitle, nBinsX, lowX, highX, lowY, highY);
0967     vDeleteObjects_.push_back(profile);
0968     return profile;
0969   } else
0970     return (TProfile*)tfd.make<TProfile>(histName, histTitle, nBinsX, lowX, highX, lowY, highY);
0971 }
0972 
0973 bool TrackerOfflineValidation::isBarrel(uint32_t subDetId) {
0974   return (subDetId == StripSubdetector::TIB || subDetId == StripSubdetector::TOB ||
0975           subDetId == PixelSubdetector::PixelBarrel);
0976 }
0977 
0978 bool TrackerOfflineValidation::isEndCap(uint32_t subDetId) {
0979   return (subDetId == StripSubdetector::TID || subDetId == StripSubdetector::TEC ||
0980           subDetId == PixelSubdetector::PixelEndcap);
0981 }
0982 
0983 bool TrackerOfflineValidation::isPixel(uint32_t subDetId) {
0984   return (subDetId == PixelSubdetector::PixelBarrel || subDetId == PixelSubdetector::PixelEndcap);
0985 }
0986 
0987 bool TrackerOfflineValidation::isDetOrDetUnit(align::StructureType type) {
0988   return (type == align::AlignableDet || type == align::AlignableDetUnit);
0989 }
0990 
0991 void TrackerOfflineValidation::getBinning(uint32_t subDetId,
0992                                           TrackerOfflineValidation::HistogramType residualType,
0993                                           int& nBinsX,
0994                                           double& lowerBoundX,
0995                                           double& upperBoundX) {
0996   // determine if
0997   const bool isPixel = this->isPixel(subDetId);
0998 
0999   edm::ParameterSet binningPSet;
1000 
1001   switch (residualType) {
1002     case XResidual:
1003       if (isPixel)
1004         binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1XResPixelModules");
1005       else
1006         binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1XResStripModules");
1007       break;
1008     case NormXResidual:
1009       if (isPixel)
1010         binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1NormXResPixelModules");
1011       else
1012         binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1NormXResStripModules");
1013       break;
1014     case XprimeResidual:
1015       if (isPixel)
1016         binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1XprimeResPixelModules");
1017       else
1018         binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1XprimeResStripModules");
1019       break;
1020     case NormXprimeResidual:
1021       if (isPixel)
1022         binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1NormXprimeResPixelModules");
1023       else
1024         binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1NormXprimeResStripModules");
1025       break;
1026     case YResidual:  // borrow y-residual binning from yprime
1027     case YprimeResidual:
1028       if (isPixel)
1029         binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1YResPixelModules");
1030       else
1031         binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1YResStripModules");
1032       break;
1033       /* case NormYResidual :*/
1034     case NormYprimeResidual:
1035       if (isPixel)
1036         binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1NormYResPixelModules");
1037       else
1038         binningPSet = parSet_.getParameter<edm::ParameterSet>("TH1NormYResStripModules");
1039       break;
1040     case XResidualProfile:
1041       if (isPixel)
1042         binningPSet = parSet_.getParameter<edm::ParameterSet>("TProfileXResPixelModules");
1043       else
1044         binningPSet = parSet_.getParameter<edm::ParameterSet>("TProfileXResStripModules");
1045       break;
1046     case YResidualProfile:
1047       if (isPixel)
1048         binningPSet = parSet_.getParameter<edm::ParameterSet>("TProfileYResPixelModules");
1049       else
1050         binningPSet = parSet_.getParameter<edm::ParameterSet>("TProfileYResStripModules");
1051       break;
1052   }
1053   nBinsX = binningPSet.getParameter<int32_t>("Nbinx");
1054   lowerBoundX = binningPSet.getParameter<double>("xmin");
1055   upperBoundX = binningPSet.getParameter<double>("xmax");
1056 }
1057 
1058 void TrackerOfflineValidation::setSummaryBin(int bin, TH1* targetHist, TH1* sourceHist) {
1059   if (targetHist && sourceHist) {
1060     targetHist->SetBinContent(bin, sourceHist->GetMean(1));
1061     if (useFwhm_)
1062       targetHist->SetBinError(bin, Fwhm(sourceHist) / 2.);
1063     else
1064       targetHist->SetBinError(bin, sourceHist->GetRMS(1));
1065   } else
1066     return;
1067 }
1068 
1069 void TrackerOfflineValidation::summarizeBinInContainer(int bin,
1070                                                        SummaryContainer& targetContainer,
1071                                                        SummaryContainer& sourceContainer) {
1072   summaryBins_.emplace_back(bin, targetContainer.summaryXResiduals_, sourceContainer.sumXResiduals_);
1073   summaryBins_.emplace_back(bin, targetContainer.summaryNormXResiduals_, sourceContainer.sumNormXResiduals_);
1074   // If no y-residual hists, just returns:
1075   summaryBins_.emplace_back(bin, targetContainer.summaryYResiduals_, sourceContainer.sumYResiduals_);
1076   summaryBins_.emplace_back(bin, targetContainer.summaryNormYResiduals_, sourceContainer.sumNormYResiduals_);
1077 }
1078 
1079 void TrackerOfflineValidation::summarizeBinInContainer(int bin,
1080                                                        uint32_t subDetId,
1081                                                        SummaryContainer& targetContainer,
1082                                                        ModuleHistos& sourceContainer) {
1083   // takes two summary Containers and sets summaryBins for all histograms
1084   summaryBins_.emplace_back(bin, targetContainer.summaryXResiduals_, sourceContainer.ResXprimeHisto);
1085   summaryBins_.emplace_back(bin, targetContainer.summaryNormXResiduals_, sourceContainer.NormResXprimeHisto);
1086   if (this->isPixel(subDetId) || stripYResiduals_) {
1087     summaryBins_.emplace_back(bin, targetContainer.summaryYResiduals_, sourceContainer.ResYprimeHisto);
1088     summaryBins_.emplace_back(bin, targetContainer.summaryNormYResiduals_, sourceContainer.NormResYprimeHisto);
1089   }
1090 }
1091 
1092 TrackerOfflineValidation::ModuleHistos& TrackerOfflineValidation::getHistStructFromMap(const DetId& detid) {
1093   // get a struct with histograms from the respective map
1094   // if no object exist, the reference is automatically created by the map
1095   // throw exception if non-tracker id is passed
1096   uint subdetid = detid.subdetId();
1097   if (subdetid == PixelSubdetector::PixelBarrel) {
1098     return mPxbResiduals_[detid.rawId()];
1099   } else if (subdetid == PixelSubdetector::PixelEndcap) {
1100     return mPxeResiduals_[detid.rawId()];
1101   } else if (subdetid == StripSubdetector::TIB) {
1102     return mTibResiduals_[detid.rawId()];
1103   } else if (subdetid == StripSubdetector::TID) {
1104     return mTidResiduals_[detid.rawId()];
1105   } else if (subdetid == StripSubdetector::TOB) {
1106     return mTobResiduals_[detid.rawId()];
1107   } else if (subdetid == StripSubdetector::TEC) {
1108     return mTecResiduals_[detid.rawId()];
1109   } else {
1110     throw cms::Exception("Geometry Error")
1111         << "[TrackerOfflineValidation] Error, tried to get reference for non-tracker subdet " << subdetid
1112         << " from detector " << detid.det();
1113     return mPxbResiduals_[0];
1114   }
1115 }
1116 
1117 // ------------ method called to for each event  ------------
1118 void TrackerOfflineValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
1119   if (maxTracks_ > 0 && nTracks_ >= maxTracks_)
1120     return;  //don't do anything after hitting the max number of tracks
1121              //(could just rely on break below, but this way saves fillTrackQuantities)
1122 
1123   if (useOverflowForRMS_)
1124     TH1::StatOverflows(kTRUE);
1125   this->checkBookHists(iSetup);  // check whether hists are booked and do so if not yet done
1126 
1127   std::vector<TrackerValidationVariables::AVTrackStruct> vTrackstruct;
1128   avalidator_.fillTrackQuantities(iEvent, iSetup, vTrackstruct);
1129 
1130   for (std::vector<TrackerValidationVariables::AVTrackStruct>::const_iterator itT = vTrackstruct.begin();
1131        itT != vTrackstruct.end();
1132        ++itT) {
1133     if (itT->charge * chargeCut_ < 0)
1134       continue;
1135 
1136     if (maxTracks_ > 0 && nTracks_++ >= maxTracks_)
1137       break;  //exit the loop after hitting the max number of tracks
1138     // Fill 1D track histos
1139     static const int etaindex = this->GetIndex(vTrackHistos_, "h_tracketa");
1140     vTrackHistos_[etaindex]->Fill(itT->eta);
1141     static const int phiindex = this->GetIndex(vTrackHistos_, "h_trackphi");
1142     vTrackHistos_[phiindex]->Fill(itT->phi);
1143     static const int numOfValidHitsindex = this->GetIndex(vTrackHistos_, "h_trackNumberOfValidHits");
1144     vTrackHistos_[numOfValidHitsindex]->Fill(itT->numberOfValidHits);
1145     static const int numOfLostHitsindex = this->GetIndex(vTrackHistos_, "h_trackNumberOfLostHits");
1146     vTrackHistos_[numOfLostHitsindex]->Fill(itT->numberOfLostHits);
1147     static const int kappaindex = this->GetIndex(vTrackHistos_, "h_curvature");
1148     vTrackHistos_[kappaindex]->Fill(itT->kappa);
1149     static const int kappaposindex = this->GetIndex(vTrackHistos_, "h_curvature_pos");
1150     if (itT->charge > 0)
1151       vTrackHistos_[kappaposindex]->Fill(fabs(itT->kappa));
1152     static const int kappanegindex = this->GetIndex(vTrackHistos_, "h_curvature_neg");
1153     if (itT->charge < 0)
1154       vTrackHistos_[kappanegindex]->Fill(fabs(itT->kappa));
1155     static const int normchi2index = this->GetIndex(vTrackHistos_, "h_normchi2");
1156     vTrackHistos_[normchi2index]->Fill(itT->normchi2);
1157     static const int chi2index = this->GetIndex(vTrackHistos_, "h_chi2");
1158     vTrackHistos_[chi2index]->Fill(itT->chi2);
1159     static const int chi2Probindex = this->GetIndex(vTrackHistos_, "h_chi2Prob");
1160     vTrackHistos_[chi2Probindex]->Fill(itT->chi2Prob);
1161     static const int ptindex = this->GetIndex(vTrackHistos_, "h_pt");
1162     vTrackHistos_[ptindex]->Fill(itT->pt);
1163     if (itT->ptError != 0.) {
1164       static const int ptResolutionindex = this->GetIndex(vTrackHistos_, "h_ptResolution");
1165       vTrackHistos_[ptResolutionindex]->Fill(itT->ptError / itT->pt);
1166     }
1167     // Fill track profiles
1168     static const int d0phiindex = this->GetIndex(vTrackProfiles_, "p_d0_vs_phi");
1169     vTrackProfiles_[d0phiindex]->Fill(itT->phi, itT->d0);
1170     static const int dzphiindex = this->GetIndex(vTrackProfiles_, "p_dz_vs_phi");
1171     vTrackProfiles_[dzphiindex]->Fill(itT->phi, itT->dz);
1172     static const int d0etaindex = this->GetIndex(vTrackProfiles_, "p_d0_vs_eta");
1173     vTrackProfiles_[d0etaindex]->Fill(itT->eta, itT->d0);
1174     static const int dzetaindex = this->GetIndex(vTrackProfiles_, "p_dz_vs_eta");
1175     vTrackProfiles_[dzetaindex]->Fill(itT->eta, itT->dz);
1176     static const int chiphiindex = this->GetIndex(vTrackProfiles_, "p_chi2_vs_phi");
1177     vTrackProfiles_[chiphiindex]->Fill(itT->phi, itT->chi2);
1178     static const int chiProbphiindex = this->GetIndex(vTrackProfiles_, "p_chi2Prob_vs_phi");
1179     vTrackProfiles_[chiProbphiindex]->Fill(itT->phi, itT->chi2Prob);
1180     static const int chiProbabsd0index = this->GetIndex(vTrackProfiles_, "p_chi2Prob_vs_d0");
1181     vTrackProfiles_[chiProbabsd0index]->Fill(fabs(itT->d0), itT->chi2Prob);
1182     static const int normchiphiindex = this->GetIndex(vTrackProfiles_, "p_normchi2_vs_phi");
1183     vTrackProfiles_[normchiphiindex]->Fill(itT->phi, itT->normchi2);
1184     static const int chietaindex = this->GetIndex(vTrackProfiles_, "p_chi2_vs_eta");
1185     vTrackProfiles_[chietaindex]->Fill(itT->eta, itT->chi2);
1186     static const int normchiptindex = this->GetIndex(vTrackProfiles_, "p_normchi2_vs_pt");
1187     vTrackProfiles_[normchiptindex]->Fill(itT->pt, itT->normchi2);
1188     static const int normchipindex = this->GetIndex(vTrackProfiles_, "p_normchi2_vs_p");
1189     vTrackProfiles_[normchipindex]->Fill(itT->p, itT->normchi2);
1190     static const int chiProbetaindex = this->GetIndex(vTrackProfiles_, "p_chi2Prob_vs_eta");
1191     vTrackProfiles_[chiProbetaindex]->Fill(itT->eta, itT->chi2Prob);
1192     static const int normchietaindex = this->GetIndex(vTrackProfiles_, "p_normchi2_vs_eta");
1193     vTrackProfiles_[normchietaindex]->Fill(itT->eta, itT->normchi2);
1194     static const int kappaphiindex = this->GetIndex(vTrackProfiles_, "p_kappa_vs_phi");
1195     vTrackProfiles_[kappaphiindex]->Fill(itT->phi, itT->kappa);
1196     static const int kappaetaindex = this->GetIndex(vTrackProfiles_, "p_kappa_vs_eta");
1197     vTrackProfiles_[kappaetaindex]->Fill(itT->eta, itT->kappa);
1198     static const int ptResphiindex = this->GetIndex(vTrackProfiles_, "p_ptResolution_vs_phi");
1199     vTrackProfiles_[ptResphiindex]->Fill(itT->phi, itT->ptError / itT->pt);
1200     static const int ptResetaindex = this->GetIndex(vTrackProfiles_, "p_ptResolution_vs_eta");
1201     vTrackProfiles_[ptResetaindex]->Fill(itT->eta, itT->ptError / itT->pt);
1202 
1203     // Fill 2D track histos
1204     static const int d0phiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_d0_vs_phi");
1205     vTrack2DHistos_[d0phiindex_2d]->Fill(itT->phi, itT->d0);
1206     static const int dzphiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_dz_vs_phi");
1207     vTrack2DHistos_[dzphiindex_2d]->Fill(itT->phi, itT->dz);
1208     static const int d0etaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_d0_vs_eta");
1209     vTrack2DHistos_[d0etaindex_2d]->Fill(itT->eta, itT->d0);
1210     static const int dzetaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_dz_vs_eta");
1211     vTrack2DHistos_[dzetaindex_2d]->Fill(itT->eta, itT->dz);
1212     static const int chiphiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_chi2_vs_phi");
1213     vTrack2DHistos_[chiphiindex_2d]->Fill(itT->phi, itT->chi2);
1214     static const int chiProbphiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_chi2Prob_vs_phi");
1215     vTrack2DHistos_[chiProbphiindex_2d]->Fill(itT->phi, itT->chi2Prob);
1216     static const int chiProbabsd0index_2d = this->GetIndex(vTrack2DHistos_, "h2_chi2Prob_vs_d0");
1217     vTrack2DHistos_[chiProbabsd0index_2d]->Fill(fabs(itT->d0), itT->chi2Prob);
1218     static const int normchiphiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_normchi2_vs_phi");
1219     vTrack2DHistos_[normchiphiindex_2d]->Fill(itT->phi, itT->normchi2);
1220     static const int chietaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_chi2_vs_eta");
1221     vTrack2DHistos_[chietaindex_2d]->Fill(itT->eta, itT->chi2);
1222     static const int chiProbetaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_chi2Prob_vs_eta");
1223     vTrack2DHistos_[chiProbetaindex_2d]->Fill(itT->eta, itT->chi2Prob);
1224     static const int normchietaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_normchi2_vs_eta");
1225     vTrack2DHistos_[normchietaindex_2d]->Fill(itT->eta, itT->normchi2);
1226     static const int kappaphiindex_2d = this->GetIndex(vTrack2DHistos_, "h2_kappa_vs_phi");
1227     vTrack2DHistos_[kappaphiindex_2d]->Fill(itT->phi, itT->kappa);
1228     static const int kappaetaindex_2d = this->GetIndex(vTrack2DHistos_, "h2_kappa_vs_eta");
1229     vTrack2DHistos_[kappaetaindex_2d]->Fill(itT->eta, itT->kappa);
1230     static const int normchi2kappa_2d = this->GetIndex(vTrack2DHistos_, "h2_normchi2_vs_kappa");
1231     vTrack2DHistos_[normchi2kappa_2d]->Fill(itT->normchi2, itT->kappa);
1232 
1233     // hit quantities: residuals, normalized residuals
1234     for (std::vector<TrackerValidationVariables::AVHitStruct>::const_iterator itH = itT->hits.begin();
1235          itH != itT->hits.end();
1236          ++itH) {
1237       DetId detid(itH->rawDetId);
1238       ModuleHistos& histStruct = this->getHistStructFromMap(detid);
1239 
1240       // fill histos in local coordinates if set in cf
1241       if (lCoorHistOn_) {
1242         histStruct.ResHisto->Fill(itH->resX);
1243         if (itH->resErrX != 0)
1244           histStruct.NormResHisto->Fill(itH->resX / itH->resErrX);
1245         if (this->isPixel(detid.subdetId()) || stripYResiduals_) {
1246           histStruct.ResYHisto->Fill(itH->resY);
1247           // here add un-primed normalised y-residuals if wanted
1248         }
1249       }
1250       if (itH->resXprime != -999.) {
1251         histStruct.ResXprimeHisto->Fill(itH->resXprime);
1252 
1253         /******************************* Fill 2-D histo ResX vs momenta *****************************/
1254         if (detid.subdetId() == PixelSubdetector::PixelBarrel) {
1255           static const int resXvsPindex_2d = this->GetIndex(vTrack2DHistos_, "p_vs_resXprime_pixB");
1256           vTrack2DHistos_[resXvsPindex_2d]->Fill(itT->p, itH->resXprime);
1257         }
1258         if (detid.subdetId() == PixelSubdetector::PixelEndcap) {
1259           static const int resXvsPindex_2d = this->GetIndex(vTrack2DHistos_, "p_vs_resXprime_pixE");
1260           vTrack2DHistos_[resXvsPindex_2d]->Fill(itT->p, itH->resXprime);
1261         }
1262         if (detid.subdetId() == StripSubdetector::TIB) {
1263           static const int resXvsPindex_2d = this->GetIndex(vTrack2DHistos_, "p_vs_resXprime_TIB");
1264           vTrack2DHistos_[resXvsPindex_2d]->Fill(itT->p, itH->resXprime);
1265         }
1266         if (detid.subdetId() == StripSubdetector::TID) {
1267           static const int resXvsPindex_2d = this->GetIndex(vTrack2DHistos_, "p_vs_resXprime_TID");
1268           vTrack2DHistos_[resXvsPindex_2d]->Fill(itT->p, itH->resXprime);
1269         }
1270         if (detid.subdetId() == StripSubdetector::TOB) {
1271           static const int resXvsPindex_2d = this->GetIndex(vTrack2DHistos_, "p_vs_resXprime_TOB");
1272           vTrack2DHistos_[resXvsPindex_2d]->Fill(itT->p, itH->resXprime);
1273         }
1274         if (detid.subdetId() == StripSubdetector::TEC) {
1275           static const int resXvsPindex_2d = this->GetIndex(vTrack2DHistos_, "p_vs_resXprime_TEC");
1276           vTrack2DHistos_[resXvsPindex_2d]->Fill(itT->p, itH->resXprime);
1277         }
1278         /******************************************/
1279 
1280         if (moduleLevelProfiles_ && itH->inside) {
1281           float tgalpha = tan(itH->localAlpha);
1282           if (fabs(tgalpha) != 0) {
1283             histStruct.LocalX->Fill(itH->localXnorm, tgalpha * tgalpha);
1284             histStruct.LocalY->Fill(itH->localYnorm, tgalpha * tgalpha);
1285             /*           if (this->isEndCap(detid.subdetId()) && !this->isPixel(detid.subdetId())) {
1286              if((itH->resX)*(itH->resXprime)>0){
1287             histStruct.ResXvsXProfile->Fill(itH->localXnorm, itH->resXatTrkY/tgalpha, tgalpha*tgalpha);
1288             histStruct.ResXvsYProfile->Fill(itH->localYnorm, itH->resXatTrkY/tgalpha, tgalpha*tgalpha);
1289              } else {
1290             histStruct.ResXvsXProfile->Fill(itH->localXnorm, (-1)*itH->resXatTrkY/tgalpha, tgalpha*tgalpha);
1291             histStruct.ResXvsYProfile->Fill(itH->localYnorm, (-1)*itH->resXatTrkY/tgalpha, tgalpha*tgalpha);
1292                }
1293 
1294           }else {  
1295 */
1296             histStruct.ResXvsXProfile->Fill(itH->localXnorm, itH->resXatTrkY / tgalpha, tgalpha * tgalpha);
1297             histStruct.ResXvsYProfile->Fill(itH->localYnorm, itH->resXatTrkY / tgalpha, tgalpha * tgalpha);
1298 
1299             //          }
1300           }
1301         }
1302 
1303         if (itH->resXprimeErr != 0 && itH->resXprimeErr != -999) {
1304           histStruct.NormResXprimeHisto->Fill(itH->resXprime / itH->resXprimeErr);
1305         }
1306       }
1307 
1308       if (itH->resYprime != -999.) {
1309         if (this->isPixel(detid.subdetId()) || stripYResiduals_) {
1310           histStruct.ResYprimeHisto->Fill(itH->resYprime);
1311 
1312           /******************************* Fill 2-D histo ResY vs momenta *****************************/
1313           if (detid.subdetId() == PixelSubdetector::PixelBarrel) {
1314             static const int resYvsPindex_2d = this->GetIndex(vTrack2DHistos_, "p_vs_resYprime_pixB");
1315             vTrack2DHistos_[resYvsPindex_2d]->Fill(itT->p, itH->resYprime);
1316           }
1317           if (detid.subdetId() == PixelSubdetector::PixelEndcap) {
1318             static const int resYvsPindex_2d = this->GetIndex(vTrack2DHistos_, "p_vs_resYprime_pixE");
1319             vTrack2DHistos_[resYvsPindex_2d]->Fill(itT->p, itH->resYprime);
1320           }
1321           /******************************************/
1322 
1323           if (moduleLevelProfiles_ && itH->inside) {
1324             float tgbeta = tan(itH->localBeta);
1325             if (fabs(tgbeta) != 0) {
1326               /*           if (this->isEndCap(detid.subdetId()) && !this->isPixel(detid.subdetId())) {
1327             
1328             if((itH->resY)*(itH->resYprime)>0){
1329             histStruct.ResYvsXProfile->Fill(itH->localXnorm, itH->resYprime/tgbeta, tgbeta*tgbeta);
1330             histStruct.ResYvsYProfile->Fill(itH->localYnorm, itH->resYprime/tgbeta, tgbeta*tgbeta);
1331              } else {
1332             histStruct.ResYvsXProfile->Fill(itH->localXnorm, (-1)*itH->resYprime/tgbeta, tgbeta*tgbeta);
1333             histStruct.ResYvsYProfile->Fill(itH->localYnorm, (-1)*itH->resYprime/tgbeta, tgbeta*tgbeta);
1334                }
1335 
1336           }else {  
1337 */
1338               histStruct.ResYvsXProfile->Fill(itH->localXnorm, itH->resY / tgbeta, tgbeta * tgbeta);
1339               histStruct.ResYvsYProfile->Fill(itH->localYnorm, itH->resY / tgbeta, tgbeta * tgbeta);
1340               //              }
1341             }
1342           }
1343 
1344           if (itH->resYprimeErr != 0 && itH->resYprimeErr != -999.) {
1345             histStruct.NormResYprimeHisto->Fill(itH->resYprime / itH->resYprimeErr);
1346           }
1347         }
1348       }
1349 
1350     }  // finish loop over hit quantities
1351   }    // finish loop over track quantities
1352 
1353   if (useOverflowForRMS_)
1354     TH1::StatOverflows(kFALSE);
1355 }
1356 
1357 // ------------ method called once each job just after ending the event loop  ------------
1358 void TrackerOfflineValidation::endJob() {
1359   if (!tkGeom_.product())
1360     return;
1361 
1362   static const int kappadiffindex = this->GetIndex(vTrackHistos_, "h_diff_curvature");
1363   vTrackHistos_[kappadiffindex]->Add(vTrackHistos_[this->GetIndex(vTrackHistos_, "h_curvature_neg")],
1364                                      vTrackHistos_[this->GetIndex(vTrackHistos_, "h_curvature_pos")],
1365                                      -1,
1366                                      1);
1367 
1368   // Collate Information for Subdetectors
1369   // create summary histograms recursively
1370   this->collateSummaryHists();
1371 
1372   if (dqmMode_)
1373     return;
1374   // Should be excluded in dqmMode, since TTree is not usable
1375   // In dqmMode tree operations are are sourced out to the additional module TrackerOfflineValidationSummary
1376 
1377   edm::Service<TFileService> fs;
1378   if (compressionSettings_ > 0) {
1379     fs->file().SetCompressionSettings(compressionSettings_);
1380   }
1381 
1382   TTree* tree = fs->make<TTree>("TkOffVal", "TkOffVal");
1383 
1384   TkOffTreeVariables* treeMemPtr = new TkOffTreeVariables;
1385   // We create branches for all members of 'TkOffTreeVariables' (even if not needed).
1386   // This works because we have a dictionary for 'TkOffTreeVariables'
1387   // (see src/classes_def.xml and src/classes.h):
1388   tree->Branch("TkOffTreeVariables", &treeMemPtr);  // address of pointer!
1389 
1390   this->fillTree(*tree, *treeMemPtr, mPxbResiduals_);
1391   this->fillTree(*tree, *treeMemPtr, mPxeResiduals_);
1392   this->fillTree(*tree, *treeMemPtr, mTibResiduals_);
1393   this->fillTree(*tree, *treeMemPtr, mTidResiduals_);
1394   this->fillTree(*tree, *treeMemPtr, mTobResiduals_);
1395   this->fillTree(*tree, *treeMemPtr, mTecResiduals_);
1396 
1397   delete treeMemPtr;
1398   treeMemPtr = nullptr;
1399 }
1400 
1401 void TrackerOfflineValidation::prepareSummaryHists(
1402     DirectoryWrapper& tfd,
1403     const Alignable& ali,
1404     std::vector<TrackerOfflineValidation::SummaryContainer>& vLevelProfiles) {
1405   const auto& alivec = ali.components();
1406   if (this->isDetOrDetUnit((alivec)[0]->alignableObjectId()))
1407     return;
1408 
1409   for (int iComp = 0, iCompEnd = ali.components().size(); iComp < iCompEnd; ++iComp) {
1410     std::vector<TrackerOfflineValidation::SummaryContainer> vProfiles;
1411     std::string structurename = alignableTracker_->objectIdProvider().idToString((alivec)[iComp]->alignableObjectId());
1412 
1413     LogDebug("TrackerOfflineValidation") << "StructureName = " << structurename;
1414     std::stringstream dirname;
1415     dirname << structurename;
1416 
1417     // add no suffix counter to strip and pixel -> just aesthetics
1418     if (structurename != "Strip" && structurename != "Pixel")
1419       dirname << "_" << iComp + 1;
1420 
1421     if (!(this->isDetOrDetUnit((alivec)[iComp]->alignableObjectId())) || (alivec)[0]->components().size() > 1) {
1422       DirectoryWrapper f(tfd, dirname.str(), moduleDirectory_, dqmMode_);
1423       this->prepareSummaryHists(f, *(alivec)[iComp], vProfiles);
1424       vLevelProfiles.push_back(this->bookSummaryHists(tfd, *(alivec[iComp]), ali.alignableObjectId(), iComp + 1));
1425       TH1* hX = vLevelProfiles[iComp].sumXResiduals_;
1426       TH1* hNormX = vLevelProfiles[iComp].sumNormXResiduals_;
1427       TH1* hY = vLevelProfiles[iComp].sumYResiduals_;
1428       TH1* hNormY = vLevelProfiles[iComp].sumNormYResiduals_;
1429       TH1* pXX = vLevelProfiles[iComp].sumResXvsXProfile_;
1430       TH1* pXY = vLevelProfiles[iComp].sumResXvsYProfile_;
1431       TH1* pYX = vLevelProfiles[iComp].sumResYvsXProfile_;
1432       TH1* pYY = vLevelProfiles[iComp].sumResYvsYProfile_;
1433       for (uint n = 0; n < vProfiles.size(); ++n) {
1434         this->summarizeBinInContainer(n + 1, vLevelProfiles[iComp], vProfiles[n]);
1435         sumHistStructure_.emplace_back(hX, vProfiles[n].sumXResiduals_);
1436         sumHistStructure_.emplace_back(hNormX, vProfiles[n].sumNormXResiduals_);
1437         if (hY)
1438           sumHistStructure_.emplace_back(hY, vProfiles[n].sumYResiduals_);  // only if existing
1439         if (hNormY)
1440           sumHistStructure_.emplace_back(hNormY, vProfiles[n].sumNormYResiduals_);  // ditto (pxl, stripYResiduals_)
1441         if (pXX)
1442           sumHistStructure_.emplace_back(pXX, vProfiles[n].sumResXvsXProfile_);
1443         if (pXY)
1444           sumHistStructure_.emplace_back(pXY, vProfiles[n].sumResXvsYProfile_);
1445         if (pYX)
1446           sumHistStructure_.emplace_back(pYX, vProfiles[n].sumResYvsXProfile_);
1447         if (pYY)
1448           sumHistStructure_.emplace_back(pYY, vProfiles[n].sumResYvsYProfile_);
1449       }
1450       if (dqmMode_)
1451         continue;  // No fits in dqmMode
1452       //add fit values to stat box
1453       toFit_.push_back(vLevelProfiles[iComp].sumXResiduals_);
1454       toFit_.push_back(vLevelProfiles[iComp].sumNormXResiduals_);
1455       if (hY)
1456         toFit_.push_back(hY);  // only if existing (pixel or stripYResiduals_)
1457       if (hNormY)
1458         toFit_.push_back(hNormY);  // ditto
1459     } else {
1460       // nothing to be done for det or detunits
1461       continue;
1462     }
1463   }
1464 }
1465 
1466 void TrackerOfflineValidation::collateSummaryHists() {
1467   for (std::vector<std::pair<TH1*, TH1*> >::const_iterator it = sumHistStructure_.begin();
1468        it != sumHistStructure_.end();
1469        ++it)
1470     it->first->Add(it->second);
1471 
1472   for (std::vector<std::tuple<int, TH1*, TH1*> >::const_iterator it = summaryBins_.begin(); it != summaryBins_.end();
1473        ++it)
1474     setSummaryBin(std::get<0>(*it), std::get<1>(*it), std::get<2>(*it));
1475 
1476   for (std::vector<TH1*>::const_iterator it = toFit_.begin(); it != toFit_.end(); ++it)
1477     fitResiduals(*it);
1478 }
1479 
1480 TrackerOfflineValidation::SummaryContainer TrackerOfflineValidation::bookSummaryHists(DirectoryWrapper& tfd,
1481                                                                                       const Alignable& ali,
1482                                                                                       align::StructureType type,
1483                                                                                       int i) {
1484   const uint aliSize = ali.components().size();
1485   const align::StructureType alitype = ali.alignableObjectId();
1486   const align::StructureType subtype = ali.components()[0]->alignableObjectId();
1487   const char* aliTypeName = alignableTracker_->objectIdProvider().idToString(alitype);  // lifetime of char* OK
1488   const char* aliSubtypeName = alignableTracker_->objectIdProvider().idToString(subtype);
1489   const char* typeName = alignableTracker_->objectIdProvider().idToString(type);
1490 
1491   const DetId aliDetId = ali.id();
1492   // y residuals only if pixel or specially requested for strip:
1493   const bool bookResidY = this->isPixel(aliDetId.subdetId()) || stripYResiduals_;
1494 
1495   SummaryContainer sumContainer;
1496 
1497   // Book summary hists with one bin per component,
1498   // but special case for Det with two DetUnit that we want to summarize one level up
1499   // (e.g. in TOBRods with 12 bins for 6 stereo and 6 rphi DetUnit.)
1500   //    component of ali is not Det or Det with just one components
1501   const uint subcompSize = ali.components()[0]->components().size();
1502   if (subtype != align::AlignableDet || subcompSize == 1) {  // Det with 1 comp. should not exist anymore...
1503     const TString title(Form("Summary for substructures in %s %d;%s;", aliTypeName, i, aliSubtypeName));
1504 
1505     sumContainer.summaryXResiduals_ = tfd.make<TH1F>(
1506         Form("h_summaryX%s_%d", aliTypeName, i), title + "#LT #Delta x' #GT", aliSize, 0.5, aliSize + 0.5);
1507     sumContainer.summaryNormXResiduals_ = tfd.make<TH1F>(
1508         Form("h_summaryNormX%s_%d", aliTypeName, i), title + "#LT #Delta x'/#sigma #GT", aliSize, 0.5, aliSize + 0.5);
1509 
1510     if (bookResidY) {
1511       sumContainer.summaryYResiduals_ = tfd.make<TH1F>(
1512           Form("h_summaryY%s_%d", aliTypeName, i), title + "#LT #Delta y' #GT", aliSize, 0.5, aliSize + 0.5);
1513       sumContainer.summaryNormYResiduals_ = tfd.make<TH1F>(
1514           Form("h_summaryNormY%s_%d", aliTypeName, i), title + "#LT #Delta y'/#sigma #GT", aliSize, 0.5, aliSize + 0.5);
1515     }
1516 
1517   } else if (subtype == align::AlignableDet && subcompSize > 1) {  // fixed: was aliSize before
1518     if (subcompSize != 2) {                                        // strange... expect only 2 DetUnits in DS layers
1519       // this 2 is hardcoded factor 2 in binning below and also assumed later on
1520       edm::LogError("Alignment") << "@SUB=bookSummaryHists"
1521                                  << "Det with " << subcompSize << " components";
1522     }
1523     // title contains x-title
1524     const TString title(Form(
1525         "Summary for substructures in %s %d;%s;",
1526         aliTypeName,
1527         i,
1528         alignableTracker_->objectIdProvider().idToString(ali.components()[0]->components()[0]->alignableObjectId())));
1529 
1530     sumContainer.summaryXResiduals_ = tfd.make<TH1F>(
1531         Form("h_summaryX%s_%d", aliTypeName, i), title + "#LT #Delta x' #GT", (2 * aliSize), 0.5, 2 * aliSize + 0.5);
1532     sumContainer.summaryNormXResiduals_ = tfd.make<TH1F>(Form("h_summaryNormX%s_%d", aliTypeName, i),
1533                                                          title + "#LT #Delta x'/#sigma #GT",
1534                                                          (2 * aliSize),
1535                                                          0.5,
1536                                                          2 * aliSize + 0.5);
1537 
1538     if (bookResidY) {
1539       sumContainer.summaryYResiduals_ = tfd.make<TH1F>(
1540           Form("h_summaryY%s_%d", aliTypeName, i), title + "#LT #Delta y' #GT", (2 * aliSize), 0.5, 2 * aliSize + 0.5);
1541       sumContainer.summaryNormYResiduals_ = tfd.make<TH1F>(Form("h_summaryNormY%s_%d", aliTypeName, i),
1542                                                            title + "#LT #Delta y'/#sigma #GT",
1543                                                            (2 * aliSize),
1544                                                            0.5,
1545                                                            2 * aliSize + 0.5);
1546     }
1547 
1548   } else {
1549     edm::LogError("TrackerOfflineValidation")
1550         << "@SUB=TrackerOfflineValidation::bookSummaryHists"
1551         << "No summary histogram for hierarchy level " << aliTypeName << " in subdet " << aliDetId.subdetId();
1552   }
1553 
1554   // Now book hists that just sum up the residual histograms from lower levels.
1555   // Axis title is copied from lowest level module of structure.
1556   // Should be safe that y-hists are only touched if non-null pointers...
1557   int nbins = 0;
1558   double xmin = 0., xmax = 0.;
1559   const TString sumTitle(Form("Residual for %s %d in %s;", aliTypeName, i, typeName));
1560   const ModuleHistos& xTitHists = this->getHistStructFromMap(aliDetId);  // for x-axis titles
1561   this->getBinning(aliDetId.subdetId(), XprimeResidual, nbins, xmin, xmax);
1562 
1563   sumContainer.sumXResiduals_ = tfd.make<TH1F>(Form("h_Xprime_%s_%d", aliTypeName, i),
1564                                                sumTitle + xTitHists.ResXprimeHisto->GetXaxis()->GetTitle(),
1565                                                nbins,
1566                                                xmin,
1567                                                xmax);
1568 
1569   this->getBinning(aliDetId.subdetId(), NormXprimeResidual, nbins, xmin, xmax);
1570   sumContainer.sumNormXResiduals_ = tfd.make<TH1F>(Form("h_NormXprime_%s_%d", aliTypeName, i),
1571                                                    sumTitle + xTitHists.NormResXprimeHisto->GetXaxis()->GetTitle(),
1572                                                    nbins,
1573                                                    xmin,
1574                                                    xmax);
1575 
1576   if (moduleLevelProfiles_) {
1577     this->getBinning(aliDetId.subdetId(), XResidualProfile, nbins, xmin, xmax);
1578     sumContainer.sumResXvsXProfile_ = tfd.make<TProfile>(Form("p_resXX_%s_%d", aliTypeName, i),
1579                                                          sumTitle + xTitHists.ResXvsXProfile->GetXaxis()->GetTitle() +
1580                                                              ";" + xTitHists.ResXvsXProfile->GetYaxis()->GetTitle(),
1581                                                          nbins,
1582                                                          xmin,
1583                                                          xmax);
1584     sumContainer.sumResXvsXProfile_->Sumw2();
1585     sumContainer.sumResXvsYProfile_ = tfd.make<TProfile>(Form("p_resXY_%s_%d", aliTypeName, i),
1586                                                          sumTitle + xTitHists.ResXvsYProfile->GetXaxis()->GetTitle() +
1587                                                              ";" + xTitHists.ResXvsYProfile->GetYaxis()->GetTitle(),
1588                                                          nbins,
1589                                                          xmin,
1590                                                          xmax);
1591     sumContainer.sumResXvsYProfile_->Sumw2();
1592   }
1593 
1594   if (bookResidY) {
1595     this->getBinning(aliDetId.subdetId(), YprimeResidual, nbins, xmin, xmax);
1596     sumContainer.sumYResiduals_ = tfd.make<TH1F>(Form("h_Yprime_%s_%d", aliTypeName, i),
1597                                                  sumTitle + xTitHists.ResYprimeHisto->GetXaxis()->GetTitle(),
1598                                                  nbins,
1599                                                  xmin,
1600                                                  xmax);
1601 
1602     this->getBinning(aliDetId.subdetId(), NormYprimeResidual, nbins, xmin, xmax);
1603     sumContainer.sumNormYResiduals_ = tfd.make<TH1F>(Form("h_NormYprime_%s_%d", aliTypeName, i),
1604                                                      sumTitle + xTitHists.NormResYprimeHisto->GetXaxis()->GetTitle(),
1605                                                      nbins,
1606                                                      xmin,
1607                                                      xmax);
1608 
1609     if (moduleLevelProfiles_) {
1610       this->getBinning(aliDetId.subdetId(), YResidualProfile, nbins, xmin, xmax);
1611       sumContainer.sumResYvsXProfile_ = tfd.make<TProfile>(Form("p_resYX_%s_%d", aliTypeName, i),
1612                                                            sumTitle + xTitHists.ResYvsXProfile->GetXaxis()->GetTitle() +
1613                                                                ";" + xTitHists.ResYvsXProfile->GetYaxis()->GetTitle(),
1614                                                            nbins,
1615                                                            xmin,
1616                                                            xmax);
1617       sumContainer.sumResYvsXProfile_->Sumw2();
1618       sumContainer.sumResYvsYProfile_ = tfd.make<TProfile>(Form("p_resYY_%s_%d", aliTypeName, i),
1619                                                            sumTitle + xTitHists.ResYvsYProfile->GetXaxis()->GetTitle() +
1620                                                                ";" + xTitHists.ResYvsYProfile->GetYaxis()->GetTitle(),
1621                                                            nbins,
1622                                                            xmin,
1623                                                            xmax);
1624       sumContainer.sumResYvsYProfile_->Sumw2();
1625     }
1626   }
1627 
1628   // If we are at the lowest level, we already sum up and fill the summary.
1629 
1630   // special case I: For DetUnits and Dets with only one subcomponent start filling summary histos
1631   if ((subtype == align::AlignableDet && subcompSize == 1) || subtype == align::AlignableDetUnit) {
1632     for (uint k = 0; k < aliSize; ++k) {
1633       DetId detid = ali.components()[k]->id();
1634       ModuleHistos& histStruct = this->getHistStructFromMap(detid);
1635       this->summarizeBinInContainer(k + 1, detid.subdetId(), sumContainer, histStruct);
1636       sumHistStructure_.emplace_back(sumContainer.sumXResiduals_, histStruct.ResXprimeHisto);
1637       sumHistStructure_.emplace_back(sumContainer.sumNormXResiduals_, histStruct.NormResXprimeHisto);
1638       if (moduleLevelProfiles_) {
1639         sumHistStructure_.emplace_back(sumContainer.sumResXvsXProfile_, histStruct.ResXvsXProfile);
1640         sumHistStructure_.emplace_back(sumContainer.sumResXvsYProfile_, histStruct.ResXvsYProfile);
1641       }
1642       if (this->isPixel(detid.subdetId()) || stripYResiduals_) {
1643         sumHistStructure_.emplace_back(sumContainer.sumYResiduals_, histStruct.ResYprimeHisto);
1644         sumHistStructure_.emplace_back(sumContainer.sumNormYResiduals_, histStruct.NormResYprimeHisto);
1645         if (moduleLevelProfiles_) {
1646           sumHistStructure_.emplace_back(sumContainer.sumResYvsXProfile_, histStruct.ResYvsXProfile);
1647           sumHistStructure_.emplace_back(sumContainer.sumResYvsYProfile_, histStruct.ResYvsYProfile);
1648         }
1649       }
1650     }
1651   } else if (subtype == align::AlignableDet && subcompSize > 1) {  // fixed: was aliSize before
1652     // special case II: Fill summary histos for dets with two detunits
1653     for (uint k = 0; k < aliSize; ++k) {
1654       for (uint j = 0; j < subcompSize; ++j) {  // assumes all have same size (as binning does)
1655         DetId detid = ali.components()[k]->components()[j]->id();
1656         ModuleHistos& histStruct = this->getHistStructFromMap(detid);
1657         this->summarizeBinInContainer(2 * k + j + 1, detid.subdetId(), sumContainer, histStruct);
1658         sumHistStructure_.emplace_back(sumContainer.sumXResiduals_, histStruct.ResXprimeHisto);
1659         sumHistStructure_.emplace_back(sumContainer.sumNormXResiduals_, histStruct.NormResXprimeHisto);
1660         if (moduleLevelProfiles_) {
1661           sumHistStructure_.emplace_back(sumContainer.sumResXvsXProfile_, histStruct.ResXvsXProfile);
1662           sumHistStructure_.emplace_back(sumContainer.sumResXvsYProfile_, histStruct.ResXvsYProfile);
1663         }
1664         if (this->isPixel(detid.subdetId()) || stripYResiduals_) {
1665           sumHistStructure_.emplace_back(sumContainer.sumYResiduals_, histStruct.ResYprimeHisto);
1666           sumHistStructure_.emplace_back(sumContainer.sumNormYResiduals_, histStruct.NormResYprimeHisto);
1667           if (moduleLevelProfiles_) {
1668             sumHistStructure_.emplace_back(sumContainer.sumResYvsXProfile_, histStruct.ResYvsXProfile);
1669             sumHistStructure_.emplace_back(sumContainer.sumResYvsYProfile_, histStruct.ResYvsYProfile);
1670           }
1671         }
1672       }
1673     }
1674   }
1675   return sumContainer;
1676 }
1677 
1678 float TrackerOfflineValidation::Fwhm(const TH1* hist) const {
1679   float max = hist->GetMaximum();
1680   int left = -1, right = -1;
1681   for (unsigned int i = 1, iEnd = hist->GetNbinsX(); i <= iEnd; ++i) {
1682     if (hist->GetBinContent(i) < max / 2. && hist->GetBinContent(i + 1) > max / 2. && left == -1) {
1683       if (max / 2. - hist->GetBinContent(i) < hist->GetBinContent(i + 1) - max / 2.) {
1684         left = i;
1685         ++i;
1686       } else {
1687         left = i + 1;
1688         ++i;
1689       }
1690     }
1691     if (left != -1 && right == -1) {
1692       if (hist->GetBinContent(i) > max / 2. && hist->GetBinContent(i + 1) < max / 2.) {
1693         if (hist->GetBinContent(i) - max / 2. < max / 2. - hist->GetBinContent(i + 1)) {
1694           right = i;
1695         } else {
1696           right = i + 1;
1697         }
1698       }
1699     }
1700   }
1701   return hist->GetXaxis()->GetBinCenter(right) - hist->GetXaxis()->GetBinCenter(left);
1702 }
1703 
1704 void TrackerOfflineValidation::setUpTreeMembers(const std::map<int, TrackerOfflineValidation::ModuleHistos>& moduleHist_,
1705                                                 const TrackerGeometry& tkgeom,
1706                                                 const TrackerTopology* tTopo) {
1707   for (std::map<int, TrackerOfflineValidation::ModuleHistos>::const_iterator it = moduleHist_.begin(),
1708                                                                              itEnd = moduleHist_.end();
1709        it != itEnd;
1710        ++it) {
1711     TkOffTreeVariables& treeMem = mTreeMembers_[it->first];
1712 
1713     //variables concerning the tracker components/hierarchy levels
1714     DetId detId_ = it->first;
1715     treeMem.moduleId = detId_;
1716     treeMem.subDetId = detId_.subdetId();
1717     treeMem.isDoubleSide = false;
1718 
1719     if (treeMem.subDetId == PixelSubdetector::PixelBarrel) {
1720       unsigned int whichHalfBarrel(1), rawId(detId_.rawId());  //DetId does not know about halfBarrels is PXB ...
1721       if ((rawId >= 302056964 && rawId < 302059300) || (rawId >= 302123268 && rawId < 302127140) ||
1722           (rawId >= 302189572 && rawId < 302194980))
1723         whichHalfBarrel = 2;
1724       treeMem.layer = tTopo->pxbLayer(detId_);
1725       treeMem.half = whichHalfBarrel;
1726       treeMem.rod = tTopo->pxbLadder(detId_);  // ... so, ladder is not per halfBarrel-Layer, but per barrel-layer!
1727       treeMem.module = tTopo->pxbModule(detId_);
1728     } else if (treeMem.subDetId == PixelSubdetector::PixelEndcap) {
1729       unsigned int whichHalfCylinder(1), rawId(detId_.rawId());  //DetId does not kmow about halfCylinders in PXF
1730       if ((rawId >= 352394500 && rawId < 352406032) || (rawId >= 352460036 && rawId < 352471568) ||
1731           (rawId >= 344005892 && rawId < 344017424) || (rawId >= 344071428 && rawId < 344082960))
1732         whichHalfCylinder = 2;
1733       treeMem.layer = tTopo->pxfDisk(detId_);
1734       treeMem.side = tTopo->pxfSide(detId_);
1735       treeMem.half = whichHalfCylinder;
1736       treeMem.blade = tTopo->pxfBlade(detId_);
1737       treeMem.panel = tTopo->pxfPanel(detId_);
1738       treeMem.module = tTopo->pxfModule(detId_);
1739     } else if (treeMem.subDetId == StripSubdetector::TIB) {
1740       unsigned int whichHalfShell(1), rawId(detId_.rawId());  //DetId does not kmow about halfShells in TIB
1741       if ((rawId >= 369120484 && rawId < 369120688) || (rawId >= 369121540 && rawId < 369121776) ||
1742           (rawId >= 369136932 && rawId < 369137200) || (rawId >= 369137988 && rawId < 369138288) ||
1743           (rawId >= 369153396 && rawId < 369153744) || (rawId >= 369154436 && rawId < 369154800) ||
1744           (rawId >= 369169844 && rawId < 369170256) || (rawId >= 369170900 && rawId < 369171344) ||
1745           (rawId >= 369124580 && rawId < 369124784) || (rawId >= 369125636 && rawId < 369125872) ||
1746           (rawId >= 369141028 && rawId < 369141296) || (rawId >= 369142084 && rawId < 369142384) ||
1747           (rawId >= 369157492 && rawId < 369157840) || (rawId >= 369158532 && rawId < 369158896) ||
1748           (rawId >= 369173940 && rawId < 369174352) || (rawId >= 369174996 && rawId < 369175440))
1749         whichHalfShell = 2;
1750       treeMem.layer = tTopo->tibLayer(detId_);
1751       treeMem.side = tTopo->tibStringInfo(detId_)[0];
1752       treeMem.half = whichHalfShell;
1753       treeMem.rod = tTopo->tibStringInfo(detId_)[2];
1754       treeMem.outerInner = tTopo->tibStringInfo(detId_)[1];
1755       treeMem.module = tTopo->tibModule(detId_);
1756       treeMem.isStereo = tTopo->tibStereo(detId_);
1757       treeMem.isDoubleSide = tTopo->tibIsDoubleSide(detId_);
1758     } else if (treeMem.subDetId == StripSubdetector::TID) {
1759       treeMem.layer = tTopo->tidWheel(detId_);
1760       treeMem.side = tTopo->tidSide(detId_);
1761       treeMem.ring = tTopo->tidRing(detId_);
1762       treeMem.outerInner = tTopo->tidModuleInfo(detId_)[0];
1763       treeMem.module = tTopo->tidModuleInfo(detId_)[1];
1764       treeMem.isStereo = tTopo->tidStereo(detId_);
1765       treeMem.isDoubleSide = tTopo->tidIsDoubleSide(detId_);
1766     } else if (treeMem.subDetId == StripSubdetector::TOB) {
1767       treeMem.layer = tTopo->tobLayer(detId_);
1768       treeMem.side = tTopo->tobRodInfo(detId_)[0];
1769       treeMem.rod = tTopo->tobRodInfo(detId_)[1];
1770       treeMem.module = tTopo->tobModule(detId_);
1771       treeMem.isStereo = tTopo->tobStereo(detId_);
1772       treeMem.isDoubleSide = tTopo->tobIsDoubleSide(detId_);
1773     } else if (treeMem.subDetId == StripSubdetector::TEC) {
1774       treeMem.layer = tTopo->tecWheel(detId_);
1775       treeMem.side = tTopo->tecSide(detId_);
1776       treeMem.ring = tTopo->tecRing(detId_);
1777       treeMem.petal = tTopo->tecPetalInfo(detId_)[1];
1778       treeMem.outerInner = tTopo->tecPetalInfo(detId_)[0];
1779       treeMem.module = tTopo->tecModule(detId_);
1780       treeMem.isStereo = tTopo->tecStereo(detId_);
1781       treeMem.isDoubleSide = tTopo->tecIsDoubleSide(detId_);
1782     }
1783 
1784     //variables concerning the tracker geometry
1785     const Surface::PositionType& gPModule = tkgeom.idToDet(detId_)->position();
1786     treeMem.posPhi = gPModule.phi();
1787     treeMem.posEta = gPModule.eta();
1788     treeMem.posR = gPModule.perp();
1789     treeMem.posX = gPModule.x();
1790     treeMem.posY = gPModule.y();
1791     treeMem.posZ = gPModule.z();
1792 
1793     const Surface& surface = tkgeom.idToDet(detId_)->surface();
1794 
1795     //global Orientation of local coordinate system of dets/detUnits
1796     LocalPoint lUDirection(1., 0., 0.), lVDirection(0., 1., 0.), lWDirection(0., 0., 1.);
1797     GlobalPoint gUDirection = surface.toGlobal(lUDirection), gVDirection = surface.toGlobal(lVDirection),
1798                 gWDirection = surface.toGlobal(lWDirection);
1799     double dR(999.), dPhi(999.), dZ(999.);
1800     if (treeMem.subDetId == PixelSubdetector::PixelBarrel || treeMem.subDetId == StripSubdetector::TIB ||
1801         treeMem.subDetId == StripSubdetector::TOB) {
1802       dR = gWDirection.perp() - gPModule.perp();
1803       dPhi = deltaPhi(gUDirection.barePhi(), gPModule.barePhi());
1804       dZ = gVDirection.z() - gPModule.z();
1805       if (dZ >= 0.)
1806         treeMem.rOrZDirection = 1;
1807       else
1808         treeMem.rOrZDirection = -1;
1809     } else if (treeMem.subDetId == PixelSubdetector::PixelEndcap) {
1810       dR = gUDirection.perp() - gPModule.perp();
1811       dPhi = deltaPhi(gVDirection.barePhi(), gPModule.barePhi());
1812       dZ = gWDirection.z() - gPModule.z();
1813       if (dR >= 0.)
1814         treeMem.rOrZDirection = 1;
1815       else
1816         treeMem.rOrZDirection = -1;
1817     } else if (treeMem.subDetId == StripSubdetector::TID || treeMem.subDetId == StripSubdetector::TEC) {
1818       dR = gVDirection.perp() - gPModule.perp();
1819       dPhi = deltaPhi(gUDirection.barePhi(), gPModule.barePhi());
1820       dZ = gWDirection.z() - gPModule.z();
1821       if (dR >= 0.)
1822         treeMem.rOrZDirection = 1;
1823       else
1824         treeMem.rOrZDirection = -1;
1825     }
1826     if (dR >= 0.)
1827       treeMem.rDirection = 1;
1828     else
1829       treeMem.rDirection = -1;
1830     if (dPhi >= 0.)
1831       treeMem.phiDirection = 1;
1832     else
1833       treeMem.phiDirection = -1;
1834     if (dZ >= 0.)
1835       treeMem.zDirection = 1;
1836     else
1837       treeMem.zDirection = -1;
1838   }
1839 }
1840 
1841 void TrackerOfflineValidation::fillTree(TTree& tree,
1842                                         TkOffTreeVariables& treeMem,
1843                                         const std::map<int, TrackerOfflineValidation::ModuleHistos>& moduleHist_) {
1844   for (std::map<int, TrackerOfflineValidation::ModuleHistos>::const_iterator it = moduleHist_.begin(),
1845                                                                              itEnd = moduleHist_.end();
1846        it != itEnd;
1847        ++it) {
1848     treeMem = mTreeMembers_[it->first];
1849 
1850     //mean and RMS values (extracted from histograms(Xprime on module level)
1851     treeMem.entries = static_cast<UInt_t>(it->second.ResXprimeHisto->GetEntries());
1852     treeMem.meanX = it->second.ResXprimeHisto->GetMean();
1853     treeMem.rmsX = it->second.ResXprimeHisto->GetRMS();
1854     //treeMem.sigmaX = Fwhm(it->second.ResXprimeHisto)/2.355;
1855 
1856     if (useFit_) {
1857       //call fit function which returns mean and sigma from the fit
1858       //for absolute residuals
1859       std::pair<float, float> fitResult1 = this->fitResiduals(it->second.ResXprimeHisto);
1860       treeMem.fitMeanX = fitResult1.first;
1861       treeMem.fitSigmaX = fitResult1.second;
1862       //for normalized residuals
1863       std::pair<float, float> fitResult2 = this->fitResiduals(it->second.NormResXprimeHisto);
1864       treeMem.fitMeanNormX = fitResult2.first;
1865       treeMem.fitSigmaNormX = fitResult2.second;
1866     }
1867 
1868     //get median for absolute residuals
1869     treeMem.medianX = this->getMedian(it->second.ResXprimeHisto);
1870 
1871     int numberOfBins = it->second.ResXprimeHisto->GetNbinsX();
1872     treeMem.numberOfUnderflows = it->second.ResXprimeHisto->GetBinContent(0);
1873     treeMem.numberOfOverflows = it->second.ResXprimeHisto->GetBinContent(numberOfBins + 1);
1874     treeMem.numberOfOutliers =
1875         it->second.ResXprimeHisto->GetBinContent(0) + it->second.ResXprimeHisto->GetBinContent(numberOfBins + 1);
1876 
1877     //mean and RMS values (extracted from histograms(normalized Xprime on module level)
1878     treeMem.meanNormX = it->second.NormResXprimeHisto->GetMean();
1879     treeMem.rmsNormX = it->second.NormResXprimeHisto->GetRMS();
1880 
1881     double stats[20];
1882     it->second.NormResXprimeHisto->GetStats(stats);
1883     // GF  treeMem.chi2PerDofX = stats[3]/(stats[0]-1);
1884     if (stats[0])
1885       treeMem.chi2PerDofX = stats[3] / stats[0];
1886 
1887     treeMem.sigmaNormX = Fwhm(it->second.NormResXprimeHisto) / 2.355;
1888     treeMem.histNameX = it->second.ResXprimeHisto->GetName();
1889     treeMem.histNameNormX = it->second.NormResXprimeHisto->GetName();
1890 
1891     // fill tree variables in local coordinates if set in cfg
1892     if (lCoorHistOn_) {
1893       treeMem.meanLocalX = it->second.ResHisto->GetMean();
1894       treeMem.rmsLocalX = it->second.ResHisto->GetRMS();
1895       treeMem.meanNormLocalX = it->second.NormResHisto->GetMean();
1896       treeMem.rmsNormLocalX = it->second.NormResHisto->GetRMS();
1897 
1898       treeMem.histNameLocalX = it->second.ResHisto->GetName();
1899       treeMem.histNameNormLocalX = it->second.NormResHisto->GetName();
1900       if (it->second.ResYHisto)
1901         treeMem.histNameLocalY = it->second.ResYHisto->GetName();
1902     }
1903 
1904     // mean and RMS values in local y (extracted from histograms(normalized Yprime on module level)
1905     // might exist in pixel only
1906     if (it->second.ResYprimeHisto) {  //(stripYResiduals_)
1907       TH1* h = it->second.ResYprimeHisto;
1908       treeMem.meanY = h->GetMean();
1909       treeMem.rmsY = h->GetRMS();
1910 
1911       if (useFit_) {  // fit function which returns mean and sigma from the fit
1912         std::pair<float, float> fitMeanSigma = this->fitResiduals(h);
1913         treeMem.fitMeanY = fitMeanSigma.first;
1914         treeMem.fitSigmaY = fitMeanSigma.second;
1915       }
1916 
1917       //get median for absolute residuals
1918       treeMem.medianY = this->getMedian(h);
1919 
1920       treeMem.histNameY = h->GetName();
1921     }
1922     if (it->second.NormResYprimeHisto) {
1923       TH1* h = it->second.NormResYprimeHisto;
1924       treeMem.meanNormY = h->GetMean();
1925       treeMem.rmsNormY = h->GetRMS();
1926       h->GetStats(stats);  // stats buffer defined above
1927       if (stats[0])
1928         treeMem.chi2PerDofY = stats[3] / stats[0];
1929 
1930       if (useFit_) {  // fit function which returns mean and sigma from the fit
1931         std::pair<float, float> fitMeanSigma = this->fitResiduals(h);
1932         treeMem.fitMeanNormY = fitMeanSigma.first;
1933         treeMem.fitSigmaNormY = fitMeanSigma.second;
1934       }
1935       treeMem.histNameNormY = h->GetName();
1936     }
1937 
1938     if (moduleLevelProfiles_) {
1939       if (it->second.ResXvsXProfile) {
1940         TH1* h = it->second.ResXvsXProfile;
1941         treeMem.meanResXvsX = h->GetMean();
1942         treeMem.rmsResXvsX = h->GetRMS();
1943         treeMem.profileNameResXvsX = h->GetName();
1944       }
1945       if (it->second.ResXvsYProfile) {
1946         TH1* h = it->second.ResXvsYProfile;
1947         treeMem.meanResXvsY = h->GetMean();
1948         treeMem.rmsResXvsY = h->GetRMS();
1949         treeMem.profileNameResXvsY = h->GetName();
1950       }
1951       if (it->second.ResYvsXProfile) {
1952         TH1* h = it->second.ResYvsXProfile;
1953         treeMem.meanResYvsX = h->GetMean();
1954         treeMem.rmsResYvsX = h->GetRMS();
1955         treeMem.profileNameResYvsX = h->GetName();
1956       }
1957       if (it->second.ResYvsYProfile) {
1958         TH1* h = it->second.ResYvsYProfile;
1959         treeMem.meanResYvsY = h->GetMean();
1960         treeMem.rmsResYvsY = h->GetRMS();
1961         treeMem.profileNameResYvsY = h->GetName();
1962       }
1963     }
1964 
1965     tree.Fill();
1966   }
1967 }
1968 
1969 std::pair<float, float> TrackerOfflineValidation::fitResiduals(TH1* hist) const {
1970   std::pair<float, float> fitResult(9999., 9999.);
1971   if (!hist || hist->GetEntries() < 20)
1972     return fitResult;
1973 
1974   float mean = hist->GetMean();
1975   float sigma = hist->GetRMS();
1976 
1977   try {  // for < CMSSW_2_2_0 since ROOT warnings from fit are converted to exceptions
1978     // Remove the try/catch for more recent CMSSW!
1979     // first fit: two RMS around mean
1980     TF1 func("tmp", "gaus", mean - 2. * sigma, mean + 2. * sigma);
1981     if (0 == hist->Fit(&func, "QNR")) {  // N: do not blow up file by storing fit!
1982       mean = func.GetParameter(1);
1983       sigma = func.GetParameter(2);
1984       // second fit: three sigma of first fit around mean of first fit
1985       func.SetRange(mean - 3. * sigma, mean + 3. * sigma);
1986       // I: integral gives more correct results if binning is too wide
1987       // L: Likelihood can treat empty bins correctly (if hist not weighted...)
1988       if (0 == hist->Fit(&func, "Q0LR")) {
1989         if (hist->GetFunction(func.GetName())) {  // Take care that it is later on drawn:
1990           hist->GetFunction(func.GetName())->ResetBit(TF1::kNotDraw);
1991         }
1992         fitResult.first = func.GetParameter(1);
1993         fitResult.second = func.GetParameter(2);
1994       }
1995     }
1996   } catch (cms::Exception const& e) {
1997     edm::LogWarning("Alignment") << "@SUB=TrackerOfflineValidation::fitResiduals"
1998                                  << "Caught this exception during ROOT fit: " << e.what();
1999   }
2000   return fitResult;
2001 }
2002 
2003 float TrackerOfflineValidation::getMedian(const TH1* histo) const {
2004   float median = 999;
2005   int nbins = histo->GetNbinsX();
2006 
2007   //extract median from histogram
2008   double* x = new double[nbins];
2009   double* y = new double[nbins];
2010   for (int j = 0; j < nbins; j++) {
2011     x[j] = histo->GetBinCenter(j + 1);
2012     y[j] = histo->GetBinContent(j + 1);
2013   }
2014   median = TMath::Median(nbins, x, y);
2015 
2016   delete[] x;
2017   x = nullptr;
2018   delete[] y;
2019   y = nullptr;
2020 
2021   return median;
2022 }
2023 //define this as a plug-in
2024 DEFINE_FWK_MODULE(TrackerOfflineValidation);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!