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