File indexing completed on 2024-04-06 11:57:10
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 = 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
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
0093 std::string hierarchyName;
0094
0095 std::string componentName;
0096
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
0127
0128 const edm::ParameterSet parSet_;
0129 edm::ESHandle<TrackerGeometry> tkGeom_;
0130 std::unique_ptr<TrackerTopology> tTopo_;
0131
0132
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
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
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
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
0190
0191
0192
0193 void TrackerOfflineValidationSummary::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0194
0195
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
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
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
0243
0244
0245 tree->Branch("TkOffTreeVariables", &treeMemPtr);
0246
0247 std::map<std::string, std::string>* substructureName = new std::map<std::string, std::string>;
0248 tree->Branch("SubstructureName", &substructureName, 32000, 00);
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
0258
0259
0260
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();
0282
0283
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());
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);
0297 treeMem.module = tTopo->pxbModule(detId);
0298 } else if (treeMem.subDetId == PixelSubdetector::PixelEndcap) {
0299 unsigned int whichHalfCylinder(1), rawId(detId.rawId());
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());
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
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
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
0411 const std::string histDir = associateModuleHistsWithTree(treeMem, it->second, substructureName);
0412
0413
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
0418 if (useFit_) {
0419
0420
0421 std::pair<float, float> fitResult1 = this->fitResiduals(it->second.ResXprimeHisto);
0422 treeMem.fitMeanX = fitResult1.first;
0423 treeMem.fitSigmaX = fitResult1.second;
0424
0425 std::pair<float, float> fitResult2 = this->fitResiduals(it->second.NormResXprimeHisto);
0426 treeMem.fitMeanNormX = fitResult2.first;
0427 treeMem.fitSigmaNormX = fitResult2.second;
0428 }
0429
0430
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
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
0445 if (stats[0])
0446 treeMem.chi2PerDofX = stats[3] / stats[0];
0447
0448
0449 treeMem.histNameX = it->second.ResXprimeHisto->GetName();
0450 treeMem.histNameNormX = it->second.NormResXprimeHisto->GetName();
0451
0452
0453 if (it->second.ResHisto && it->second.NormResHisto) {
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
0463
0464 if (it->second.ResYprimeHisto) {
0465 TH1* h = it->second.ResYprimeHisto;
0466 treeMem.meanY = h->GetMean();
0467 treeMem.rmsY = h->GetRMS();
0468
0469 if (useFit_) {
0470 std::pair<float, float> fitMeanSigma = this->fitResiduals(h);
0471 treeMem.fitMeanY = fitMeanSigma.first;
0472 treeMem.fitSigmaY = fitMeanSigma.second;
0473 }
0474
0475
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);
0486 if (stats[0])
0487 treeMem.chi2PerDofY = stats[3] / stats[0];
0488
0489 if (useFit_) {
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_.length() != 0)
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 {
0696
0697
0698 TF1 func("tmp", "gaus", mean - 2. * sigma, mean + 2. * sigma);
0699 if (0 == hist->Fit(&func, "QNR")) {
0700 mean = func.GetParameter(1);
0701 sigma = func.GetParameter(2);
0702
0703 func.SetRange(mean - 3. * sigma, mean + 3. * sigma);
0704
0705
0706 if (0 == hist->Fit(&func, "Q0LR")) {
0707 if (hist->GetFunction(func.GetName())) {
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
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
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
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
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.length() == 0) {
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
0787
0788
0789
0790
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.length() != 0)
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
0871 DEFINE_FWK_MODULE(TrackerOfflineValidationSummary);