File indexing completed on 2021-09-16 23:01:36
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021 #include <memory>
0022
0023
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
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
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
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
0091 std::string hierarchyName;
0092
0093 std::string componentName;
0094
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
0125
0126 const edm::ParameterSet parSet_;
0127 edm::ESHandle<TrackerGeometry> tkGeom_;
0128 std::unique_ptr<TrackerTopology> tTopo_;
0129
0130
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
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
0161 dbe_ = edm::Service<DQMStore>().operator->();
0162 }
0163
0164 TrackerOfflineValidationSummary::~TrackerOfflineValidationSummary() = default;
0165
0166
0167
0168
0169
0170 void TrackerOfflineValidationSummary::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0171
0172
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
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
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
0220
0221
0222 tree->Branch("TkOffTreeVariables", &treeMemPtr);
0223
0224 std::map<std::string, std::string>* substructureName = new std::map<std::string, std::string>;
0225 tree->Branch("SubstructureName", &substructureName, 32000, 00);
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
0235
0236
0237
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();
0259
0260
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());
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);
0274 treeMem.module = tTopo->pxbModule(detId);
0275 } else if (treeMem.subDetId == PixelSubdetector::PixelEndcap) {
0276 unsigned int whichHalfCylinder(1), rawId(detId.rawId());
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());
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
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
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
0388 const std::string histDir = associateModuleHistsWithTree(treeMem, it->second, substructureName);
0389
0390
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
0395 if (useFit_) {
0396
0397
0398 std::pair<float, float> fitResult1 = this->fitResiduals(it->second.ResXprimeHisto);
0399 treeMem.fitMeanX = fitResult1.first;
0400 treeMem.fitSigmaX = fitResult1.second;
0401
0402 std::pair<float, float> fitResult2 = this->fitResiduals(it->second.NormResXprimeHisto);
0403 treeMem.fitMeanNormX = fitResult2.first;
0404 treeMem.fitSigmaNormX = fitResult2.second;
0405 }
0406
0407
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
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
0422 if (stats[0])
0423 treeMem.chi2PerDofX = stats[3] / stats[0];
0424
0425
0426 treeMem.histNameX = it->second.ResXprimeHisto->GetName();
0427 treeMem.histNameNormX = it->second.NormResXprimeHisto->GetName();
0428
0429
0430 if (it->second.ResHisto && it->second.NormResHisto) {
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
0440
0441 if (it->second.ResYprimeHisto) {
0442 TH1* h = it->second.ResYprimeHisto;
0443 treeMem.meanY = h->GetMean();
0444 treeMem.rmsY = h->GetRMS();
0445
0446 if (useFit_) {
0447 std::pair<float, float> fitMeanSigma = this->fitResiduals(h);
0448 treeMem.fitMeanY = fitMeanSigma.first;
0449 treeMem.fitSigmaY = fitMeanSigma.second;
0450 }
0451
0452
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);
0463 if (stats[0])
0464 treeMem.chi2PerDofY = stats[3] / stats[0];
0465
0466 if (useFit_) {
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 {
0673
0674
0675 TF1 func("tmp", "gaus", mean - 2. * sigma, mean + 2. * sigma);
0676 if (0 == hist->Fit(&func, "QNR")) {
0677 mean = func.GetParameter(1);
0678 sigma = func.GetParameter(2);
0679
0680 func.SetRange(mean - 3. * sigma, mean + 3. * sigma);
0681
0682
0683 if (0 == hist->Fit(&func, "Q0LR")) {
0684 if (hist->GetFunction(func.GetName())) {
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
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
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
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
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
0764
0765
0766
0767
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
0848 DEFINE_FWK_MODULE(TrackerOfflineValidationSummary);