Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-01-04 00:29:06 

0001 // -*- C++ -*-
0002 //
0003 // Package:    TrackerOfflineValidationSummary
0004 // Class:      TrackerOfflineValidationSummary
0005 //
0006 /**\class TrackerOfflineValidationSummary TrackerOfflineValidationSummary.cc Alignment/TrackerOfflineValidationSummary/src/TrackerOfflineValidationSummary.cc
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Johannes Hauk
0015 //         Created:  Sat Aug 22 10:31:34 CEST 2009
0016 // $Id: TrackerOfflineValidationSummary.cc,v 1.7 2012/06/30 09:09:47 eulisse Exp $
0017 //
0018 //
0019 
0020 // system include files
0021 #include <memory>
0022 
0023 // user include files
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0026 
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029 
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 
0032 #include "FWCore/Framework/interface/ESHandle.h"
0033 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0034 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0035 #include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
0036 
0037 #include "TTree.h"
0038 //#include "Alignment/OfflineValidation/interface/TrackerValidationVariables.h"
0039 #include "Alignment/OfflineValidation/interface/TkOffTreeVariables.h"
0040 
0041 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0042 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0043 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0044 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0045 #include "DataFormats/GeometrySurface/interface/Surface.h"
0046 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0047 #include "DataFormats/Math/interface/deltaPhi.h"
0048 #include "TH1.h"
0049 #include "TMath.h"
0050 
0051 #include "DQMServices/Core/interface/DQMStore.h"
0052 #include "FWCore/ServiceRegistry/interface/Service.h"
0053 
0054 //
0055 // class decleration
0056 //
0057 
0058 class TrackerOfflineValidationSummary : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
0059 public:
0060   typedef dqm::legacy::DQMStore DQMStore;
0061   explicit TrackerOfflineValidationSummary(const edm::ParameterSet&);
0062   ~TrackerOfflineValidationSummary() override = default;
0063 
0064   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0065 
0066 private:
0067   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
0068   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0069 
0070   struct ModuleHistos {
0071     ModuleHistos()
0072         : ResHisto(), NormResHisto(), ResXprimeHisto(), NormResXprimeHisto(), ResYprimeHisto(), NormResYprimeHisto() {}
0073     TH1* ResHisto;
0074     TH1* NormResHisto;
0075     TH1* ResXprimeHisto;
0076     TH1* NormResXprimeHisto;
0077     TH1* ResYprimeHisto;
0078     TH1* NormResYprimeHisto;
0079   };
0080 
0081   /// Put here the histograms created during harvesting
0082   struct HarvestingHistos {
0083     HarvestingHistos() : DmrXprime(), DmrYprime() {}
0084     TH1* DmrXprime;
0085     TH1* DmrYprime;
0086   };
0087 
0088   struct HarvestingHierarchy {
0089     HarvestingHierarchy() {}
0090     HarvestingHierarchy(const std::string& name, const std::string& component, const std::vector<unsigned int>& entries)
0091         : hierarchyName(name), componentName(component), treeEntries(entries) {}
0092     // Needed for naming of histogram according to residual histograms
0093     std::string hierarchyName;
0094     // "Pixel" or "Strip"
0095     std::string componentName;
0096     // Modules belonging to selected hierarchy
0097     std::vector<unsigned int> treeEntries;
0098     HarvestingHistos harvestingHistos;
0099   };
0100 
0101   void beginRun(const edm::Run&, const edm::EventSetup& iSetup) override {}
0102   void analyze(const edm::Event& evt, const edm::EventSetup&) override;
0103   void endRun(const edm::Run&, const edm::EventSetup& iSetup) override;
0104   void endJob() override;
0105 
0106   void fillTree(TTree& tree,
0107                 std::map<int, TrackerOfflineValidationSummary::ModuleHistos>& moduleHist,
0108                 TkOffTreeVariables& treeMem,
0109                 const TrackerGeometry& tkgeom,
0110                 std::map<std::string, std::string>& substructureName,
0111                 const TrackerTopology* tTopo);
0112 
0113   std::pair<float, float> fitResiduals(TH1* hist) const;
0114   float getMedian(const TH1* hist) const;
0115 
0116   const std::string associateModuleHistsWithTree(const TkOffTreeVariables& treeMem,
0117                                                  TrackerOfflineValidationSummary::ModuleHistos& moduleHists,
0118                                                  std::map<std::string, std::string>& substructureName);
0119 
0120   void collateHarvestingHists(TTree& tree);
0121   void applyHarvestingHierarchy(TTree& treeMem);
0122   void bookHarvestingHists();
0123   void getBinning(const std::string& binningPSetName, int& nBinsX, double& lowerBoundX, double& upperBoundX) const;
0124   void fillHarvestingHists(TTree& tree);
0125 
0126   // ----------member data ---------------------------
0127 
0128   const edm::ParameterSet parSet_;
0129   edm::ESHandle<TrackerGeometry> tkGeom_;
0130   std::unique_ptr<TrackerTopology> tTopo_;
0131 
0132   // parameters from cfg to steer
0133   const std::string moduleDirectory_;
0134   const bool useFit_;
0135 
0136   DQMStore* dbe_;
0137 
0138   bool moduleMapsInitialized;
0139 
0140   std::map<int, TrackerOfflineValidationSummary::ModuleHistos> mPxbResiduals_;
0141   std::map<int, TrackerOfflineValidationSummary::ModuleHistos> mPxeResiduals_;
0142   std::map<int, TrackerOfflineValidationSummary::ModuleHistos> mTibResiduals_;
0143   std::map<int, TrackerOfflineValidationSummary::ModuleHistos> mTidResiduals_;
0144   std::map<int, TrackerOfflineValidationSummary::ModuleHistos> mTobResiduals_;
0145   std::map<int, TrackerOfflineValidationSummary::ModuleHistos> mTecResiduals_;
0146 
0147   std::vector<HarvestingHierarchy> vHarvestingHierarchy_;
0148 };
0149 
0150 //
0151 // constructors and destructor
0152 //
0153 TrackerOfflineValidationSummary::TrackerOfflineValidationSummary(const edm::ParameterSet& iConfig)
0154     : geomToken_(esConsumes()),
0155       tTopoToken_(esConsumes<edm::Transition::EndRun>()),
0156       parSet_(iConfig),
0157       tTopo_(nullptr),
0158       moduleDirectory_(parSet_.getParameter<std::string>("moduleDirectoryInOutput")),
0159       useFit_(parSet_.getParameter<bool>("useFit")),
0160       dbe_(nullptr),
0161       moduleMapsInitialized(false) {
0162   //now do what ever initialization is needed
0163   dbe_ = edm::Service<DQMStore>().operator->();
0164 }
0165 
0166 void TrackerOfflineValidationSummary::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0167   edm::ParameterSetDescription desc;
0168   desc.setComment("Summary of alignment payloads validation by evaluating unbiased track-hit resisuals");
0169   desc.add<std::string>("moduleDirectoryInOutput", "Alignment/Tracker");
0170   desc.add<bool>("useFit", false);
0171   desc.add<bool>("stripYDmrs", false);
0172   desc.add<unsigned int>("minEntriesPerModuleForDmr", 100);
0173 
0174   // fill in the residuals details
0175   const std::vector<std::string> listOfDMRPSets = {
0176       "TH1DmrXprimeStripModules", "TH1DmrYprimeStripModules", "TH1DmrXprimePixelModules", "TH1DmrYprimePixelModules"};
0177 
0178   for (const auto& myPSetName : listOfDMRPSets) {
0179     edm::ParameterSetDescription myPSet;
0180     myPSet.add<int>("Nbinx", 100);
0181     myPSet.add<double>("xmin", -5.f);
0182     myPSet.add<double>("xmax", 5.f);
0183     desc.add<edm::ParameterSetDescription>(myPSetName, myPSet);
0184   }
0185   descriptions.addWithDefaultLabel(desc);
0186 }
0187 
0188 //
0189 // member functions
0190 //
0191 
0192 // ------------ method called to for each event  ------------
0193 void TrackerOfflineValidationSummary::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0194   // Access of EventSetup is needed to get the list of silicon-modules and their IDs
0195   // Since they do not change, it is accessed only once
0196   if (moduleMapsInitialized)
0197     return;
0198   tkGeom_ = iSetup.getHandle(geomToken_);
0199   const TrackerGeometry* bareTkGeomPtr = &(*tkGeom_);
0200 
0201   const TrackingGeometry::DetIdContainer& detIdContainer = bareTkGeomPtr->detIds();
0202   std::vector<DetId>::const_iterator iDet;
0203   for (iDet = detIdContainer.begin(); iDet != detIdContainer.end(); ++iDet) {
0204     const DetId& detId = *iDet;
0205     const uint32_t rawId = detId.rawId();
0206     const unsigned int subdetId = detId.subdetId();
0207     if (subdetId == PixelSubdetector::PixelBarrel)
0208       mPxbResiduals_[rawId];
0209     else if (subdetId == PixelSubdetector::PixelEndcap)
0210       mPxeResiduals_[rawId];
0211     else if (subdetId == StripSubdetector::TIB)
0212       mTibResiduals_[rawId];
0213     else if (subdetId == StripSubdetector::TID)
0214       mTidResiduals_[rawId];
0215     else if (subdetId == StripSubdetector::TOB)
0216       mTobResiduals_[rawId];
0217     else if (subdetId == StripSubdetector::TEC)
0218       mTecResiduals_[rawId];
0219     else {
0220       throw cms::Exception("Geometry Error")
0221           << "[TrackerOfflineValidationSummary] Error, tried to get reference for non-tracker subdet " << subdetId
0222           << " from detector " << detId.det();
0223     }
0224   }
0225   moduleMapsInitialized = true;
0226 }
0227 
0228 // ------------ method called at each end of Run  ------------
0229 void TrackerOfflineValidationSummary::endRun(const edm::Run&, const edm::EventSetup& iSetup) {
0230   if (!tTopo_) {
0231     tTopo_ = std::make_unique<TrackerTopology>(iSetup.getData(tTopoToken_));
0232   }
0233 }
0234 
0235 // ------------ method called once each job just after ending the event loop  ------------
0236 void TrackerOfflineValidationSummary::endJob() {
0237   AlignableTracker aliTracker(&(*tkGeom_), tTopo_.get());
0238 
0239   TTree* tree = new TTree("TkOffVal", "TkOffVal");
0240 
0241   TkOffTreeVariables* treeMemPtr = new TkOffTreeVariables;
0242   // We create branches for all members of 'TkOffTreeVariables' (even if not needed).
0243   // This works because we have a dictionary for 'TkOffTreeVariables'
0244   // (see src/classes_def.xml and src/classes.h):
0245   tree->Branch("TkOffTreeVariables", &treeMemPtr);  // address of pointer!
0246   // second branch needed for assigning names and titles to harvesting histograms consistent to others
0247   std::map<std::string, std::string>* substructureName = new std::map<std::string, std::string>;
0248   tree->Branch("SubstructureName", &substructureName, 32000, 00);  // SplitLevel must be set to zero
0249 
0250   this->fillTree(*tree, mPxbResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo_.get());
0251   this->fillTree(*tree, mPxeResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo_.get());
0252   this->fillTree(*tree, mTibResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo_.get());
0253   this->fillTree(*tree, mTidResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo_.get());
0254   this->fillTree(*tree, mTobResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo_.get());
0255   this->fillTree(*tree, mTecResiduals_, *treeMemPtr, *tkGeom_, *substructureName, tTopo_.get());
0256 
0257   //dbe_->save("dqmOut.root");
0258 
0259   // Method for filling histograms which show summarized values (mean, rms, median ...)
0260   // of the module-based histograms from TrackerOfflineValidation
0261   this->collateHarvestingHists(*tree);
0262 
0263   delete tree;
0264   tree = nullptr;
0265   delete treeMemPtr;
0266   treeMemPtr = nullptr;
0267   delete substructureName;
0268   substructureName = nullptr;
0269 }
0270 
0271 void TrackerOfflineValidationSummary::fillTree(TTree& tree,
0272                                                std::map<int, TrackerOfflineValidationSummary::ModuleHistos>& moduleHist,
0273                                                TkOffTreeVariables& treeMem,
0274                                                const TrackerGeometry& tkgeom,
0275                                                std::map<std::string, std::string>& substructureName,
0276                                                const TrackerTopology* tTopo) {
0277   for (std::map<int, TrackerOfflineValidationSummary::ModuleHistos>::iterator it = moduleHist.begin(),
0278                                                                               itEnd = moduleHist.end();
0279        it != itEnd;
0280        ++it) {
0281     treeMem.clear();  // make empty/default
0282 
0283     //variables concerning the tracker components/hierarchy levels
0284     const DetId detId = it->first;
0285     treeMem.moduleId = detId;
0286     treeMem.subDetId = detId.subdetId();
0287 
0288     if (treeMem.subDetId == PixelSubdetector::PixelBarrel) {
0289       unsigned int whichHalfBarrel(1), rawId(detId.rawId());  //DetId does not know about halfBarrels is PXB ...
0290       if ((rawId >= 302056964 && rawId < 302059300) || (rawId >= 302123268 && rawId < 302127140) ||
0291           (rawId >= 302189572 && rawId < 302194980))
0292         whichHalfBarrel = 2;
0293       treeMem.layer = tTopo->pxbLayer(detId);
0294       treeMem.half = whichHalfBarrel;
0295       treeMem.rod = tTopo->pxbLadder(
0296           detId);  // ... so, ladder is not per halfBarrel-Layer, but per barrel-layer! Needs complicated calculation in associateModuleHistsWithTree()
0297       treeMem.module = tTopo->pxbModule(detId);
0298     } else if (treeMem.subDetId == PixelSubdetector::PixelEndcap) {
0299       unsigned int whichHalfCylinder(1), rawId(detId.rawId());  //DetId does not kmow about halfCylinders in PXF
0300       if ((rawId >= 352394500 && rawId < 352406032) || (rawId >= 352460036 && rawId < 352471568) ||
0301           (rawId >= 344005892 && rawId < 344017424) || (rawId >= 344071428 && rawId < 344082960))
0302         whichHalfCylinder = 2;
0303       treeMem.layer = tTopo->pxfDisk(detId);
0304       treeMem.side = tTopo->pxfSide(detId);
0305       treeMem.half = whichHalfCylinder;
0306       treeMem.blade = tTopo->pxfBlade(detId);
0307       treeMem.panel = tTopo->pxfPanel(detId);
0308       treeMem.module = tTopo->pxfModule(detId);
0309     } else if (treeMem.subDetId == StripSubdetector::TIB) {
0310       unsigned int whichHalfShell(1), rawId(detId.rawId());  //DetId does not kmow about halfShells in TIB
0311       if ((rawId >= 369120484 && rawId < 369120688) || (rawId >= 369121540 && rawId < 369121776) ||
0312           (rawId >= 369136932 && rawId < 369137200) || (rawId >= 369137988 && rawId < 369138288) ||
0313           (rawId >= 369153396 && rawId < 369153744) || (rawId >= 369154436 && rawId < 369154800) ||
0314           (rawId >= 369169844 && rawId < 369170256) || (rawId >= 369170900 && rawId < 369171344) ||
0315           (rawId >= 369124580 && rawId < 369124784) || (rawId >= 369125636 && rawId < 369125872) ||
0316           (rawId >= 369141028 && rawId < 369141296) || (rawId >= 369142084 && rawId < 369142384) ||
0317           (rawId >= 369157492 && rawId < 369157840) || (rawId >= 369158532 && rawId < 369158896) ||
0318           (rawId >= 369173940 && rawId < 369174352) || (rawId >= 369174996 && rawId < 369175440))
0319         whichHalfShell = 2;
0320       treeMem.layer = tTopo->tibLayer(detId);
0321       treeMem.side = tTopo->tibStringInfo(detId)[0];
0322       treeMem.half = whichHalfShell;
0323       treeMem.rod = tTopo->tibStringInfo(detId)[2];
0324       treeMem.outerInner = tTopo->tibStringInfo(detId)[1];
0325       treeMem.module = tTopo->tibModule(detId);
0326       treeMem.isStereo = tTopo->tibStereo(detId);
0327       treeMem.isDoubleSide = tTopo->tibIsDoubleSide(detId);
0328     } else if (treeMem.subDetId == StripSubdetector::TID) {
0329       treeMem.layer = tTopo->tidWheel(detId);
0330       treeMem.side = tTopo->tidSide(detId);
0331       treeMem.ring = tTopo->tidRing(detId);
0332       treeMem.outerInner = tTopo->tidModuleInfo(detId)[0];
0333       treeMem.module = tTopo->tidModuleInfo(detId)[1];
0334       treeMem.isStereo = tTopo->tidStereo(detId);
0335       treeMem.isDoubleSide = tTopo->tidIsDoubleSide(detId);
0336     } else if (treeMem.subDetId == StripSubdetector::TOB) {
0337       treeMem.layer = tTopo->tobLayer(detId);
0338       treeMem.side = tTopo->tobRodInfo(detId)[0];
0339       treeMem.rod = tTopo->tobRodInfo(detId)[1];
0340       treeMem.module = tTopo->tobModule(detId);
0341       treeMem.isStereo = tTopo->tobStereo(detId);
0342       treeMem.isDoubleSide = tTopo->tobIsDoubleSide(detId);
0343     } else if (treeMem.subDetId == StripSubdetector::TEC) {
0344       treeMem.layer = tTopo->tecWheel(detId);
0345       treeMem.side = tTopo->tecSide(detId);
0346       treeMem.ring = tTopo->tecRing(detId);
0347       treeMem.petal = tTopo->tecPetalInfo(detId)[1];
0348       treeMem.outerInner = tTopo->tecPetalInfo(detId)[0];
0349       treeMem.module = tTopo->tecModule(detId);
0350       treeMem.isStereo = tTopo->tecStereo(detId);
0351       treeMem.isDoubleSide = tTopo->tecIsDoubleSide(detId);
0352     }
0353 
0354     //variables concerning the tracker geometry
0355 
0356     const Surface::PositionType& gPModule = tkgeom.idToDet(detId)->position();
0357     treeMem.posPhi = gPModule.phi();
0358     treeMem.posEta = gPModule.eta();
0359     treeMem.posR = gPModule.perp();
0360     treeMem.posX = gPModule.x();
0361     treeMem.posY = gPModule.y();
0362     treeMem.posZ = gPModule.z();
0363 
0364     const Surface& surface = tkgeom.idToDet(detId)->surface();
0365 
0366     //global Orientation of local coordinate system of dets/detUnits
0367     LocalPoint lUDirection(1., 0., 0.), lVDirection(0., 1., 0.), lWDirection(0., 0., 1.);
0368     GlobalPoint gUDirection = surface.toGlobal(lUDirection), gVDirection = surface.toGlobal(lVDirection),
0369                 gWDirection = surface.toGlobal(lWDirection);
0370     double dR(999.), dPhi(999.), dZ(999.);
0371     if (treeMem.subDetId == PixelSubdetector::PixelBarrel || treeMem.subDetId == StripSubdetector::TIB ||
0372         treeMem.subDetId == StripSubdetector::TOB) {
0373       dR = gWDirection.perp() - gPModule.perp();
0374       dPhi = deltaPhi(gUDirection.barePhi(), gPModule.barePhi());
0375       dZ = gVDirection.z() - gPModule.z();
0376       if (dZ >= 0.)
0377         treeMem.rOrZDirection = 1;
0378       else
0379         treeMem.rOrZDirection = -1;
0380     } else if (treeMem.subDetId == PixelSubdetector::PixelEndcap) {
0381       dR = gUDirection.perp() - gPModule.perp();
0382       dPhi = deltaPhi(gVDirection.barePhi(), gPModule.barePhi());
0383       dZ = gWDirection.z() - gPModule.z();
0384       if (dR >= 0.)
0385         treeMem.rOrZDirection = 1;
0386       else
0387         treeMem.rOrZDirection = -1;
0388     } else if (treeMem.subDetId == StripSubdetector::TID || treeMem.subDetId == StripSubdetector::TEC) {
0389       dR = gVDirection.perp() - gPModule.perp();
0390       dPhi = deltaPhi(gUDirection.barePhi(), gPModule.barePhi());
0391       dZ = gWDirection.z() - gPModule.z();
0392       if (dR >= 0.)
0393         treeMem.rOrZDirection = 1;
0394       else
0395         treeMem.rOrZDirection = -1;
0396     }
0397     if (dR >= 0.)
0398       treeMem.rDirection = 1;
0399     else
0400       treeMem.rDirection = -1;
0401     if (dPhi >= 0.)
0402       treeMem.phiDirection = 1;
0403     else
0404       treeMem.phiDirection = -1;
0405     if (dZ >= 0.)
0406       treeMem.zDirection = 1;
0407     else
0408       treeMem.zDirection = -1;
0409 
0410     // Assign histos from first step (TrackerOfflineValidation) to the module's entry in the TTree for retrieving mean, rms, median ...
0411     const std::string histDir = associateModuleHistsWithTree(treeMem, it->second, substructureName);
0412 
0413     //mean and RMS values (extracted from histograms Xprime on module level)
0414     treeMem.entries = static_cast<UInt_t>(it->second.ResXprimeHisto->GetEntries());
0415     treeMem.meanX = it->second.ResXprimeHisto->GetMean();
0416     treeMem.rmsX = it->second.ResXprimeHisto->GetRMS();
0417     //treeMem.sigmaX = Fwhm(it->second.ResXprimeHisto)/2.355;
0418     if (useFit_) {
0419       //call fit function which returns mean and sigma from the fit
0420       //for absolute residuals
0421       std::pair<float, float> fitResult1 = this->fitResiduals(it->second.ResXprimeHisto);
0422       treeMem.fitMeanX = fitResult1.first;
0423       treeMem.fitSigmaX = fitResult1.second;
0424       //for normalized residuals
0425       std::pair<float, float> fitResult2 = this->fitResiduals(it->second.NormResXprimeHisto);
0426       treeMem.fitMeanNormX = fitResult2.first;
0427       treeMem.fitSigmaNormX = fitResult2.second;
0428     }
0429 
0430     //get median for absolute residuals
0431     treeMem.medianX = this->getMedian(it->second.ResXprimeHisto);
0432 
0433     int numberOfBins = it->second.ResXprimeHisto->GetNbinsX();
0434     treeMem.numberOfUnderflows = it->second.ResXprimeHisto->GetBinContent(0);
0435     treeMem.numberOfOverflows = it->second.ResXprimeHisto->GetBinContent(numberOfBins + 1);
0436     treeMem.numberOfOutliers =
0437         it->second.ResXprimeHisto->GetBinContent(0) + it->second.ResXprimeHisto->GetBinContent(numberOfBins + 1);
0438     //mean and RMS values (extracted from histograms(normalized Xprime on module level)
0439     treeMem.meanNormX = it->second.NormResXprimeHisto->GetMean();
0440     treeMem.rmsNormX = it->second.NormResXprimeHisto->GetRMS();
0441 
0442     double stats[20];
0443     it->second.NormResXprimeHisto->GetStats(stats);
0444     // GF  treeMem.chi2PerDofX = stats[3]/(stats[0]-1);
0445     if (stats[0])
0446       treeMem.chi2PerDofX = stats[3] / stats[0];
0447 
0448     //treeMem.sigmaNormX = Fwhm(it->second.NormResXprimeHisto)/2.355;
0449     treeMem.histNameX = it->second.ResXprimeHisto->GetName();
0450     treeMem.histNameNormX = it->second.NormResXprimeHisto->GetName();
0451 
0452     // fill tree variables in local coordinates if set in cfg of TrackerOfllineValidation
0453     if (it->second.ResHisto && it->second.NormResHisto) {  // if(lCoorHistOn_) {
0454       treeMem.meanLocalX = it->second.ResHisto->GetMean();
0455       treeMem.rmsLocalX = it->second.ResHisto->GetRMS();
0456       treeMem.meanNormLocalX = it->second.NormResHisto->GetMean();
0457       treeMem.rmsNormLocalX = it->second.NormResHisto->GetRMS();
0458       treeMem.histNameLocalX = it->second.ResHisto->GetName();
0459       treeMem.histNameNormLocalX = it->second.NormResHisto->GetName();
0460     }
0461 
0462     // mean and RMS values in local y (extracted from histograms Yprime on module level)
0463     // might exist in pixel only
0464     if (it->second.ResYprimeHisto) {  //(stripYResiduals_){
0465       TH1* h = it->second.ResYprimeHisto;
0466       treeMem.meanY = h->GetMean();
0467       treeMem.rmsY = h->GetRMS();
0468 
0469       if (useFit_) {  // fit function which returns mean and sigma from the fit
0470         std::pair<float, float> fitMeanSigma = this->fitResiduals(h);
0471         treeMem.fitMeanY = fitMeanSigma.first;
0472         treeMem.fitSigmaY = fitMeanSigma.second;
0473       }
0474 
0475       //get median for absolute residuals
0476       treeMem.medianY = this->getMedian(h);
0477 
0478       treeMem.histNameY = h->GetName();
0479     }
0480 
0481     if (it->second.NormResYprimeHisto) {
0482       TH1* h = it->second.NormResYprimeHisto;
0483       treeMem.meanNormY = h->GetMean();
0484       treeMem.rmsNormY = h->GetRMS();
0485       h->GetStats(stats);  // stats buffer defined above
0486       if (stats[0])
0487         treeMem.chi2PerDofY = stats[3] / stats[0];
0488 
0489       if (useFit_) {  // fit function which returns mean and sigma from the fit
0490         std::pair<float, float> fitMeanSigma = this->fitResiduals(h);
0491         treeMem.fitMeanNormY = fitMeanSigma.first;
0492         treeMem.fitSigmaNormY = fitMeanSigma.second;
0493       }
0494       treeMem.histNameNormY = h->GetName();
0495     }
0496 
0497     tree.Fill();
0498   }
0499 }
0500 
0501 const std::string TrackerOfflineValidationSummary::associateModuleHistsWithTree(
0502     const TkOffTreeVariables& treeMem,
0503     TrackerOfflineValidationSummary::ModuleHistos& moduleHists,
0504     std::map<std::string, std::string>& substructureName) {
0505   std::stringstream histDir, sSubdetName;
0506   std::string componentName;
0507   if (!moduleDirectory_.empty())
0508     histDir << moduleDirectory_ << "/";
0509   std::string wheelOrLayer("_layer_");
0510   if (treeMem.subDetId == PixelSubdetector::PixelBarrel) {
0511     unsigned int half(treeMem.half), layer(treeMem.layer), ladder(0);
0512     if (layer == 1) {
0513       if (half == 2)
0514         ladder = treeMem.rod - 5;
0515       else if (treeMem.rod > 15)
0516         ladder = treeMem.rod - 10;
0517       else
0518         ladder = treeMem.rod;
0519     } else if (layer == 2) {
0520       if (half == 2)
0521         ladder = treeMem.rod - 8;
0522       else if (treeMem.rod > 24)
0523         ladder = treeMem.rod - 16;
0524       else
0525         ladder = treeMem.rod;
0526     } else if (layer == 3) {
0527       if (half == 2)
0528         ladder = treeMem.rod - 11;
0529       else if (treeMem.rod > 33)
0530         ladder = treeMem.rod - 22;
0531       else
0532         ladder = treeMem.rod;
0533     }
0534     componentName = "Pixel";
0535     sSubdetName << "TPBBarrel_1";
0536     histDir << componentName << "/" << sSubdetName.str() << "/TPBHalfBarrel_" << treeMem.half << "/TPBLayer_"
0537             << treeMem.layer << "/TPBLadder_" << ladder;
0538   } else if (treeMem.subDetId == PixelSubdetector::PixelEndcap) {
0539     unsigned int side(treeMem.side), half(treeMem.half), blade(0);
0540     if (side == 1)
0541       side = 3;
0542     if (half == 2)
0543       blade = treeMem.blade - 6;
0544     else if (treeMem.blade > 18)
0545       blade = treeMem.blade - 12;
0546     else
0547       blade = treeMem.blade;
0548     componentName = "Pixel";
0549     sSubdetName << "TPEEndcap_" << side;
0550     histDir << componentName << "/" << sSubdetName.str() << "/TPEHalfCylinder_" << treeMem.half << "/TPEHalfDisk_"
0551             << treeMem.layer << "/TPEBlade_" << blade << "/TPEPanel_" << treeMem.panel;
0552     wheelOrLayer = "_wheel_";
0553   } else if (treeMem.subDetId == StripSubdetector::TIB) {
0554     unsigned int half(treeMem.half), layer(treeMem.layer), surface(treeMem.outerInner), string(0);
0555     if (half == 2) {
0556       if (layer == 1) {
0557         if (surface == 1)
0558           string = treeMem.rod - 13;
0559         else if (surface == 2)
0560           string = treeMem.rod - 15;
0561       }
0562       if (layer == 2) {
0563         if (surface == 1)
0564           string = treeMem.rod - 17;
0565         else if (surface == 2)
0566           string = treeMem.rod - 19;
0567       }
0568       if (layer == 3) {
0569         if (surface == 1)
0570           string = treeMem.rod - 22;
0571         else if (surface == 2)
0572           string = treeMem.rod - 23;
0573       }
0574       if (layer == 4) {
0575         if (surface == 1)
0576           string = treeMem.rod - 26;
0577         else if (surface == 2)
0578           string = treeMem.rod - 28;
0579       }
0580     } else
0581       string = treeMem.rod;
0582     std::stringstream detString;
0583     if (treeMem.layer < 3 && !treeMem.isDoubleSide)
0584       detString << "/Det_" << treeMem.module;
0585     else
0586       detString << "";
0587     componentName = "Strip";
0588     sSubdetName << "TIBBarrel_1";
0589     histDir << componentName << "/" << sSubdetName.str() << "/TIBHalfBarrel_" << treeMem.side << "/TIBLayer_"
0590             << treeMem.layer << "/TIBHalfShell_" << treeMem.half << "/TIBSurface_" << treeMem.outerInner
0591             << "/TIBString_" << string << detString.str();
0592   } else if (treeMem.subDetId == StripSubdetector::TID) {
0593     unsigned int side(treeMem.side), outerInner(0);
0594     if (side == 1)
0595       side = 3;
0596     if (treeMem.outerInner == 1)
0597       outerInner = 2;
0598     else if (treeMem.outerInner == 2)
0599       outerInner = 1;
0600     std::stringstream detString;
0601     if (treeMem.ring < 3 && !treeMem.isDoubleSide)
0602       detString << "/Det_" << treeMem.module;
0603     else
0604       detString << "";
0605     componentName = "Strip";
0606     sSubdetName << "TIDEndcap_" << side;
0607     histDir << componentName << "/" << sSubdetName.str() << "/TIDDisk_" << treeMem.layer << "/TIDRing_" << treeMem.ring
0608             << "/TIDSide_" << outerInner << detString.str();
0609     wheelOrLayer = "_wheel_";
0610   } else if (treeMem.subDetId == StripSubdetector::TOB) {
0611     std::stringstream detString;
0612     if (treeMem.layer < 3 && !treeMem.isDoubleSide)
0613       detString << "/Det_" << treeMem.module;
0614     else
0615       detString << "";
0616     componentName = "Strip";
0617     sSubdetName << "TOBBarrel_4";
0618     histDir << componentName << "/" << sSubdetName.str() << "/TOBHalfBarrel_" << treeMem.side << "/TOBLayer_"
0619             << treeMem.layer << "/TOBRod_" << treeMem.rod << detString.str();
0620   } else if (treeMem.subDetId == StripSubdetector::TEC) {
0621     unsigned int side(0), outerInner(0), ring(0);
0622     if (treeMem.side == 1)
0623       side = 6;
0624     else if (treeMem.side == 2)
0625       side = 5;
0626     if (treeMem.outerInner == 1)
0627       outerInner = 2;
0628     else if (treeMem.outerInner == 2)
0629       outerInner = 1;
0630     if (treeMem.layer > 3 && treeMem.layer < 7)
0631       ring = treeMem.ring - 1;
0632     else if (treeMem.layer == 7 || treeMem.layer == 8)
0633       ring = treeMem.ring - 2;
0634     else if (treeMem.layer == 9)
0635       ring = treeMem.ring - 3;
0636     else
0637       ring = treeMem.ring;
0638     std::stringstream detString;
0639     if ((treeMem.ring < 3 || treeMem.ring == 5) && !treeMem.isDoubleSide)
0640       detString << "/Det_" << treeMem.module;
0641     else
0642       detString << "";
0643     componentName = "Strip";
0644     sSubdetName << "TECEndcap_" << side;
0645     histDir << componentName << "/" << sSubdetName.str() << "/TECDisk_" << treeMem.layer << "/TECSide_" << outerInner
0646             << "/TECPetal_" << treeMem.petal << "/TECRing_" << ring << detString.str();
0647     wheelOrLayer = "_wheel_";
0648   }
0649 
0650   substructureName["component"] = componentName;
0651   substructureName["subdet"] = sSubdetName.str();
0652 
0653   std::stringstream histName;
0654   histName << "residuals_subdet_" << treeMem.subDetId << wheelOrLayer << treeMem.layer << "_module_"
0655            << treeMem.moduleId;
0656 
0657   std::string fullPath;
0658   fullPath = histDir.str() + "/h_xprime_" + histName.str();
0659   if (dbe_->get(fullPath))
0660     moduleHists.ResXprimeHisto = dbe_->get(fullPath)->getTH1();
0661   else {
0662     edm::LogError("TrackerOfflineValidationSummary")
0663         << "Problem with names in input file produced in TrackerOfflineValidation ...\n"
0664         << "This histogram should exist in every configuration, "
0665         << "but no histogram with name " << fullPath << " is found!";
0666     return "";
0667   }
0668   fullPath = histDir.str() + "/h_normxprime" + histName.str();
0669   if (dbe_->get(fullPath))
0670     moduleHists.NormResXprimeHisto = dbe_->get(fullPath)->getTH1();
0671   fullPath = histDir.str() + "/h_yprime_" + histName.str();
0672   if (dbe_->get(fullPath))
0673     moduleHists.ResYprimeHisto = dbe_->get(fullPath)->getTH1();
0674   fullPath = histDir.str() + "/h_normyprime" + histName.str();
0675   if (dbe_->get(fullPath))
0676     moduleHists.NormResYprimeHisto = dbe_->get(fullPath)->getTH1();
0677   fullPath = histDir.str() + "/h_" + histName.str();
0678   if (dbe_->get(fullPath))
0679     moduleHists.ResHisto = dbe_->get(fullPath)->getTH1();
0680   fullPath = histDir.str() + "/h_norm" + histName.str();
0681   if (dbe_->get(fullPath))
0682     moduleHists.NormResHisto = dbe_->get(fullPath)->getTH1();
0683 
0684   return histDir.str();
0685 }
0686 
0687 std::pair<float, float> TrackerOfflineValidationSummary::fitResiduals(TH1* hist) const {
0688   std::pair<float, float> fitResult(9999., 9999.);
0689   if (!hist || hist->GetEntries() < 20)
0690     return fitResult;
0691 
0692   float mean = hist->GetMean();
0693   float sigma = hist->GetRMS();
0694 
0695   try {  // for < CMSSW_2_2_0 since ROOT warnings from fit are converted to exceptions
0696     // Remove the try/catch for more recent CMSSW!
0697     // first fit: two RMS around mean
0698     TF1 func("tmp", "gaus", mean - 2. * sigma, mean + 2. * sigma);
0699     if (0 == hist->Fit(&func, "QNR")) {  // N: do not blow up file by storing fit!
0700       mean = func.GetParameter(1);
0701       sigma = func.GetParameter(2);
0702       // second fit: three sigma of first fit around mean of first fit
0703       func.SetRange(mean - 3. * sigma, mean + 3. * sigma);
0704       // I: integral gives more correct results if binning is too wide
0705       // L: Likelihood can treat empty bins correctly (if hist not weighted...)
0706       if (0 == hist->Fit(&func, "Q0LR")) {
0707         if (hist->GetFunction(func.GetName())) {  // Take care that it is later on drawn:
0708           hist->GetFunction(func.GetName())->ResetBit(TF1::kNotDraw);
0709         }
0710         fitResult.first = func.GetParameter(1);
0711         fitResult.second = func.GetParameter(2);
0712       }
0713     }
0714   } catch (cms::Exception const& e) {
0715     edm::LogWarning("Alignment") << "@SUB=TrackerOfflineValidation::fitResiduals"
0716                                  << "Caught this exception during ROOT fit: " << e.what();
0717   }
0718   return fitResult;
0719 }
0720 
0721 float TrackerOfflineValidationSummary::getMedian(const TH1* histo) const {
0722   float median = 999;
0723   const int nbins = histo->GetNbinsX();
0724 
0725   //extract median from histogram
0726   double* x = new double[nbins];
0727   double* y = new double[nbins];
0728   for (int j = 0; j < nbins; j++) {
0729     x[j] = histo->GetBinCenter(j + 1);
0730     y[j] = histo->GetBinContent(j + 1);
0731   }
0732   median = TMath::Median(nbins, x, y);
0733 
0734   delete[] x;
0735   x = nullptr;
0736   delete[] y;
0737   y = nullptr;
0738 
0739   return median;
0740 }
0741 
0742 void TrackerOfflineValidationSummary::collateHarvestingHists(TTree& tree) {
0743   this->applyHarvestingHierarchy(tree);
0744   this->fillHarvestingHists(tree);
0745 }
0746 
0747 void TrackerOfflineValidationSummary::applyHarvestingHierarchy(TTree& tree) {
0748   TkOffTreeVariables* treeMemPtr = nullptr;
0749   std::map<std::string, std::string>* substructureName = nullptr;
0750   tree.SetBranchAddress("TkOffTreeVariables", &treeMemPtr);
0751   tree.SetBranchAddress("SubstructureName", &substructureName);
0752 
0753   // Loop over modules to select accumulation criteria for harvesting plots
0754   for (unsigned int iSubdet = 1; iSubdet < 7; ++iSubdet) {
0755     std::string hierarchyName("");
0756     std::string componentName("");
0757     std::vector<unsigned int> treeEntries;
0758     for (unsigned int iSide = 1; iSide < 3; ++iSide) {
0759       // Set up only one collection for Barrels, not separated for side
0760       if (iSide == 1 && (iSubdet == PixelSubdetector::PixelBarrel || iSubdet == StripSubdetector::TIB ||
0761                          iSubdet == StripSubdetector::TOB))
0762         continue;
0763       for (int iTree = 0; iTree < tree.GetEntries(); ++iTree) {
0764         tree.GetEntry(iTree);
0765         // Do not use glued Dets
0766         if (treeMemPtr->isDoubleSide)
0767           continue;
0768         if (treeMemPtr->subDetId == iSubdet) {
0769           if (iSide != treeMemPtr->side && (iSubdet == PixelSubdetector::PixelEndcap ||
0770                                             iSubdet == StripSubdetector::TID || iSubdet == StripSubdetector::TEC))
0771             continue;
0772           treeEntries.push_back(iTree);
0773           if (hierarchyName.empty()) {
0774             hierarchyName = (*substructureName)["subdet"];
0775             componentName = (*substructureName)["component"];
0776           }
0777         }
0778       }
0779       HarvestingHierarchy harvestingHierarchy(hierarchyName, componentName, treeEntries);
0780       vHarvestingHierarchy_.push_back(harvestingHierarchy);
0781       hierarchyName = "";
0782       componentName = "";
0783       treeEntries.clear();
0784     }
0785   }
0786   // Here could be a further separation of the HarvestingHierarchy.
0787   // E.g. separate the existing ones by layer and add them to the vector without deleting any element from the vector.
0788   // The existing hists will stay and the new ones are added
0789 
0790   // Now, book the corresponding histos
0791   this->bookHarvestingHists();
0792 }
0793 
0794 void TrackerOfflineValidationSummary::bookHarvestingHists() {
0795   edm::LogInfo("TrackerOfflineValidationSummary") << "Harvesting histograms will be booked for "
0796                                                   << vHarvestingHierarchy_.size() << " different hierarchy selections";
0797   for (std::vector<HarvestingHierarchy>::iterator iHier = vHarvestingHierarchy_.begin();
0798        iHier != vHarvestingHierarchy_.end();
0799        ++iHier) {
0800     std::stringstream dmrXprimeHistoName, dmrYprimeHistoName, dmrXprimeHistoTitle, dmrYprimeHistoTitle;
0801     dmrXprimeHistoName << "h_DmrXprime_" << iHier->hierarchyName;
0802     dmrYprimeHistoName << "h_DmrYprime_" << iHier->hierarchyName;
0803     dmrXprimeHistoTitle << "DMR for " << iHier->hierarchyName << ";<#DeltaX> [cm];# modules";
0804     dmrYprimeHistoTitle << "DMR for " << iHier->hierarchyName << ";<#DeltaY> [cm];# modules";
0805 
0806     std::string directoryString(moduleDirectory_);
0807     if (!directoryString.empty())
0808       directoryString += "/";
0809     directoryString += iHier->componentName;
0810     dbe_->setCurrentFolder(directoryString);
0811 
0812     int nBinsX(0);
0813     double xMin(0.), xMax(0.);
0814     if (iHier->componentName == "Pixel") {
0815       this->getBinning("TH1DmrXprimePixelModules", nBinsX, xMin, xMax);
0816       iHier->harvestingHistos.DmrXprime =
0817           dbe_->book1D(dmrXprimeHistoName.str(), dmrXprimeHistoTitle.str(), nBinsX, xMin, xMax)->getTH1();
0818       this->getBinning("TH1DmrYprimePixelModules", nBinsX, xMin, xMax);
0819       iHier->harvestingHistos.DmrYprime =
0820           dbe_->book1D(dmrYprimeHistoName.str(), dmrYprimeHistoTitle.str(), nBinsX, xMin, xMax)->getTH1();
0821     } else if (iHier->componentName == "Strip") {
0822       this->getBinning("TH1DmrXprimeStripModules", nBinsX, xMin, xMax);
0823       iHier->harvestingHistos.DmrXprime =
0824           dbe_->book1D(dmrXprimeHistoName.str(), dmrXprimeHistoTitle.str(), nBinsX, xMin, xMax)->getTH1();
0825       if (!parSet_.getParameter<bool>("stripYDmrs"))
0826         continue;
0827       this->getBinning("TH1DmrYprimeStripModules", nBinsX, xMin, xMax);
0828       iHier->harvestingHistos.DmrYprime =
0829           dbe_->book1D(dmrYprimeHistoName.str(), dmrYprimeHistoTitle.str(), nBinsX, xMin, xMax)->getTH1();
0830     }
0831   }
0832 }
0833 
0834 void TrackerOfflineValidationSummary::getBinning(const std::string& binningPSetName,
0835                                                  int& nBinsX,
0836                                                  double& lowerBoundX,
0837                                                  double& upperBoundX) const {
0838   const edm::ParameterSet& binningPSet = parSet_.getParameter<edm::ParameterSet>(binningPSetName);
0839   nBinsX = binningPSet.getParameter<int>("Nbinx");
0840   lowerBoundX = binningPSet.getParameter<double>("xmin");
0841   upperBoundX = binningPSet.getParameter<double>("xmax");
0842 }
0843 
0844 void TrackerOfflineValidationSummary::fillHarvestingHists(TTree& tree) {
0845   TkOffTreeVariables* treeMemPtr = nullptr;
0846   std::map<std::string, std::string>* substructureName = nullptr;
0847   tree.SetBranchAddress("TkOffTreeVariables", &treeMemPtr);
0848   tree.SetBranchAddress("SubstructureName", &substructureName);
0849 
0850   const unsigned int minEntriesPerModule(parSet_.getParameter<unsigned int>("minEntriesPerModuleForDmr"));
0851   edm::LogInfo("TrackerOfflineValidationSummary")
0852       << "Median of a module is added to DMR plots if it contains at least " << minEntriesPerModule << " hits";
0853 
0854   for (std::vector<HarvestingHierarchy>::iterator iHier = vHarvestingHierarchy_.begin();
0855        iHier != vHarvestingHierarchy_.end();
0856        ++iHier) {
0857     for (std::vector<unsigned int>::const_iterator iTreeEntries = iHier->treeEntries.begin();
0858          iTreeEntries != iHier->treeEntries.end();
0859          ++iTreeEntries) {
0860       tree.GetEntry(*iTreeEntries);
0861       if (treeMemPtr->entries < minEntriesPerModule)
0862         continue;
0863       iHier->harvestingHistos.DmrXprime->Fill(treeMemPtr->medianX);
0864       if (iHier->harvestingHistos.DmrYprime)
0865         iHier->harvestingHistos.DmrYprime->Fill(treeMemPtr->medianY);
0866     }
0867   }
0868 }
0869 
0870 //define this as a plug-in
0871 DEFINE_FWK_MODULE(TrackerOfflineValidationSummary);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!