Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-09-16 23:01:36 

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