OfflineValidation/plugins/TrackerOfflineValidation.cc

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