File indexing completed on 2024-09-07 04:35:06
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 #include <fmt/format.h>
0020 #include <fmt/printf.h>
0021 #include <fstream>
0022
0023
0024 #include "CalibTracker/SiPixelLorentzAngle/interface/SiPixelLorentzAngleCalibrationStruct.h"
0025 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0026 #include "CondFormats/DataRecord/interface/SiPixelLorentzAngleRcd.h"
0027 #include "CondFormats/SiPixelObjects/interface/SiPixelLorentzAngle.h"
0028 #include "DQMServices/Core/interface/DQMEDHarvester.h"
0029 #include "DataFormats/TrackerCommon/interface/PixelBarrelName.h"
0030 #include "DataFormats/TrackerCommon/interface/PixelEndcapName.h"
0031 #include "FWCore/Framework/interface/EventSetup.h"
0032 #include "FWCore/Framework/interface/Frameworkfwd.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035 #include "FWCore/ServiceRegistry/interface/Service.h"
0036 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0037 #include "Geometry/TrackerGeometryBuilder/interface/PixelTopologyMap.h"
0038 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0039 #include "MagneticField/Engine/interface/MagneticField.h"
0040 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0041
0042
0043 #include "TFitResult.h"
0044
0045
0046
0047
0048 namespace SiPixelLAHarvest {
0049
0050 enum covStatus { kUndefined = -1, kNotCalculated = 0, kApproximated = 1, kMadePosDef = 2, kAccurate = 3 };
0051 enum cuts { kZeroChi2 = 0, kChi2Cut = 1, kCovStatus = 2, kNentries = 3 };
0052
0053 struct fitResults {
0054 public:
0055 fitResults() {
0056
0057 p0 = p1 = p2 = p3 = p4 = p5 = 0.;
0058 e0 = e1 = e2 = e3 = e4 = e5 = 0.;
0059 chi2 = prob = dSq = redChi2 = -9999.;
0060 tan_LA = error_LA = -9999.;
0061 fitStatus = covMatrixStatus = ndf = nentries = -999;
0062 quality = {0b0000};
0063 };
0064
0065 void setQualityCutBit(const SiPixelLAHarvest::cuts& cut) {
0066 switch (cut) {
0067 case kZeroChi2:
0068 quality.set(0);
0069 break;
0070 case kChi2Cut:
0071 quality.set(1);
0072 break;
0073 case kCovStatus:
0074 quality.set(2);
0075 break;
0076 case kNentries:
0077 quality.set(3);
0078 break;
0079 default:
0080 throw cms::Exception("Inconsistent Data") << "Passed inexistent cut value";
0081 }
0082 }
0083
0084 double p0, e0;
0085 double p1, e1;
0086 double p2, e2;
0087 double p3, e3;
0088 double p4, e4;
0089 double p5, e5;
0090 double chi2;
0091 int ndf;
0092 int nentries;
0093 double prob;
0094 int fitStatus, covMatrixStatus;
0095 double dSq;
0096 double redChi2;
0097 double tan_LA;
0098 double error_LA;
0099 std::bitset<4> quality;
0100 };
0101 }
0102
0103
0104 class SiPixelLorentzAnglePCLHarvester : public DQMEDHarvester {
0105 public:
0106 SiPixelLorentzAnglePCLHarvester(const edm::ParameterSet&);
0107 ~SiPixelLorentzAnglePCLHarvester() override = default;
0108 void beginRun(const edm::Run&, const edm::EventSetup&) override;
0109
0110 static void fillDescriptions(edm::ConfigurationDescriptions&);
0111
0112 private:
0113 void dqmEndJob(DQMStore::IBooker&, DQMStore::IGetter&) override;
0114 void endRun(const edm::Run&, const edm::EventSetup&) override;
0115 void findMean(MonitorElement* h_drift_depth_adc_slice_, int i, int i_ring);
0116 SiPixelLAHarvest::fitResults fitAndStore(std::shared_ptr<SiPixelLorentzAngle> theLA, int i_idx, int i_lay, int i_mod);
0117 bool isFitGood(SiPixelLAHarvest::fitResults& res);
0118
0119
0120 edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomEsToken_;
0121 edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoEsTokenBR_, topoEsTokenER_;
0122 edm::ESGetToken<SiPixelLorentzAngle, SiPixelLorentzAngleRcd> siPixelLAEsToken_;
0123 edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
0124
0125 std::vector<std::string> newmodulelist_;
0126 const std::string dqmDir_;
0127 const bool doChebyshevFit_;
0128 const int order_;
0129 const double fitChi2Cut_;
0130 const std::vector<double> fitRange_;
0131 const int minHitsCut_;
0132 const std::string recordName_;
0133 std::unique_ptr<TF1> f1_;
0134 float width_;
0135 float theMagField_{0.f};
0136
0137 static constexpr float teslaToInverseGeV_ = 2.99792458e-3f;
0138 std::pair<double, double> theFitRange_{5., 280.};
0139
0140 SiPixelLorentzAngleCalibrationHistograms hists_;
0141 const SiPixelLorentzAngle* currentLorentzAngle_;
0142 std::unique_ptr<TrackerTopology> theTrackerTopology_;
0143 };
0144
0145
0146 SiPixelLorentzAnglePCLHarvester::SiPixelLorentzAnglePCLHarvester(const edm::ParameterSet& iConfig)
0147 : geomEsToken_(esConsumes<edm::Transition::BeginRun>()),
0148 topoEsTokenBR_(esConsumes<edm::Transition::BeginRun>()),
0149 topoEsTokenER_(esConsumes<edm::Transition::EndRun>()),
0150 siPixelLAEsToken_(esConsumes<edm::Transition::BeginRun>()),
0151 magneticFieldToken_(esConsumes<edm::Transition::BeginRun>()),
0152 newmodulelist_(iConfig.getParameter<std::vector<std::string>>("newmodulelist")),
0153 dqmDir_(iConfig.getParameter<std::string>("dqmDir")),
0154 doChebyshevFit_(iConfig.getParameter<bool>("doChebyshevFit")),
0155 order_(iConfig.getParameter<int>("order")),
0156 fitChi2Cut_(iConfig.getParameter<double>("fitChi2Cut")),
0157 fitRange_(iConfig.getParameter<std::vector<double>>("fitRange")),
0158 minHitsCut_(iConfig.getParameter<int>("minHitsCut")),
0159 recordName_(iConfig.getParameter<std::string>("record")) {
0160
0161 if (fitRange_.size() == 2) {
0162 theFitRange_.first = fitRange_[0];
0163 theFitRange_.second = fitRange_[1];
0164 } else {
0165 throw cms::Exception("SiPixelLorentzAnglePCLHarvester") << "Too many fit range parameters specified";
0166 }
0167
0168
0169 edm::Service<cond::service::PoolDBOutputService> poolDbService;
0170 if (!poolDbService.isAvailable())
0171 throw cms::Exception("SiPixelLorentzAnglePCLHarvester") << "PoolDBService required";
0172 }
0173
0174
0175 void SiPixelLorentzAnglePCLHarvester::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0176
0177 const TrackerGeometry* geom = &iSetup.getData(geomEsToken_);
0178 const TrackerTopology* tTopo = &iSetup.getData(topoEsTokenBR_);
0179
0180 const MagneticField* magField = &iSetup.getData(magneticFieldToken_);
0181 currentLorentzAngle_ = &iSetup.getData(siPixelLAEsToken_);
0182
0183
0184
0185
0186
0187
0188
0189 theMagField_ = 1.f / (magField->inverseBzAtOriginInGeV() * teslaToInverseGeV_);
0190
0191 PixelTopologyMap map = PixelTopologyMap(geom, tTopo);
0192 hists_.nlay = geom->numberOfLayers(PixelSubdetector::PixelBarrel);
0193 hists_.nModules_.resize(hists_.nlay);
0194 hists_.nLadders_.resize(hists_.nlay);
0195 for (int i = 0; i < hists_.nlay; i++) {
0196 hists_.nModules_[i] = map.getPXBModules(i + 1);
0197 hists_.nLadders_[i] = map.getPXBLadders(i + 1);
0198 }
0199
0200
0201 if (!hists_.BPixnewDetIds_.empty() || !hists_.FPixnewDetIds_.empty())
0202 return;
0203
0204 if (!newmodulelist_.empty()) {
0205 for (auto const& modulename : newmodulelist_) {
0206 if (modulename.find("BPix_") != std::string::npos) {
0207 PixelBarrelName bn(modulename, true);
0208 const auto& detId = bn.getDetId(tTopo);
0209 hists_.BPixnewmodulename_.push_back(modulename);
0210 hists_.BPixnewDetIds_.push_back(detId.rawId());
0211 hists_.BPixnewModule_.push_back(bn.moduleName());
0212 hists_.BPixnewLayer_.push_back(bn.layerName());
0213 } else if (modulename.find("FPix_") != std::string::npos) {
0214 PixelEndcapName en(modulename, true);
0215 const auto& detId = en.getDetId(tTopo);
0216 hists_.FPixnewmodulename_.push_back(modulename);
0217 hists_.FPixnewDetIds_.push_back(detId.rawId());
0218 hists_.FPixnewDisk_.push_back(en.diskName());
0219 hists_.FPixnewBlade_.push_back(en.bladeName());
0220 }
0221 }
0222 }
0223
0224 uint count = 0;
0225 for (const auto& id : hists_.BPixnewDetIds_) {
0226 LogDebug("SiPixelLorentzAnglePCLHarvester") << id;
0227 count++;
0228 }
0229 LogDebug("SiPixelLorentzAnglePCLHarvester") << "Stored a total of " << count << " new detIds.";
0230
0231
0232 if (!hists_.detIdsList.empty())
0233 return;
0234
0235 std::vector<uint32_t> treatedIndices;
0236
0237 for (const auto& det : geom->detsPXB()) {
0238 const PixelGeomDetUnit* pixelDet = dynamic_cast<const PixelGeomDetUnit*>(det);
0239 width_ = pixelDet->surface().bounds().thickness();
0240 const auto& layer = tTopo->pxbLayer(pixelDet->geographicalId());
0241 const auto& module = tTopo->pxbModule(pixelDet->geographicalId());
0242 int i_index = module + (layer - 1) * hists_.nModules_[layer - 1];
0243
0244 uint32_t rawId = pixelDet->geographicalId().rawId();
0245
0246
0247 if (std::find(hists_.BPixnewDetIds_.begin(), hists_.BPixnewDetIds_.end(), rawId) != hists_.BPixnewDetIds_.end())
0248 continue;
0249
0250 if (std::find(treatedIndices.begin(), treatedIndices.end(), i_index) != treatedIndices.end()) {
0251 hists_.detIdsList[i_index].push_back(rawId);
0252 } else {
0253 hists_.detIdsList.insert(std::pair<uint32_t, std::vector<uint32_t>>(i_index, {rawId}));
0254 treatedIndices.push_back(i_index);
0255 }
0256 }
0257
0258 count = 0;
0259 for (const auto& i : treatedIndices) {
0260 for (const auto& id : hists_.detIdsList[i]) {
0261 LogDebug("SiPixelLorentzAnglePCLHarvester") << id;
0262 count++;
0263 };
0264 }
0265 LogDebug("SiPixelLorentzAnglePCLHarvester") << "Stored a total of " << count << " detIds.";
0266 }
0267
0268
0269 void SiPixelLorentzAnglePCLHarvester::endRun(edm::Run const& run, edm::EventSetup const& isetup) {
0270 if (!theTrackerTopology_) {
0271 theTrackerTopology_ = std::make_unique<TrackerTopology>(isetup.getData(topoEsTokenER_));
0272 }
0273 }
0274
0275
0276 void SiPixelLorentzAnglePCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter) {
0277
0278 iGetter.cd();
0279 iGetter.setCurrentFolder(dqmDir_);
0280
0281
0282 for (int i_layer = 1; i_layer <= hists_.nlay; i_layer++) {
0283 const auto& prefix_ = fmt::sprintf("%s/BPix/BPixLayer%i", dqmDir_, i_layer);
0284 for (int i_module = 1; i_module <= hists_.nModules_[i_layer - 1]; i_module++) {
0285 int i_index = i_module + (i_layer - 1) * hists_.nModules_[i_layer - 1];
0286
0287 hists_.h_drift_depth_[i_index] =
0288 iGetter.get(fmt::format("{}/h_drift_depth_layer{}_module{}", prefix_, i_layer, i_module));
0289
0290 if (hists_.h_drift_depth_[i_index] == nullptr) {
0291 edm::LogError("SiPixelLorentzAnglePCLHarvester::dqmEndJob")
0292 << "Failed to retrieve electron drift over depth for layer " << i_layer << ", module " << i_module << ".";
0293 continue;
0294 }
0295
0296 hists_.h_drift_depth_adc_[i_index] =
0297 iGetter.get(fmt::format("{}/h_drift_depth_adc_layer{}_module{}", prefix_, i_layer, i_module));
0298
0299 hists_.h_drift_depth_adc2_[i_index] =
0300 iGetter.get(fmt::format("{}/h_drift_depth_adc2_layer{}_module{}", prefix_, i_layer, i_module));
0301
0302 hists_.h_drift_depth_noadc_[i_index] =
0303 iGetter.get(fmt::format("{}/h_drift_depth_noadc_layer{}_module{}", prefix_, i_layer, i_module));
0304
0305 hists_.h_mean_[i_index] = iGetter.get(fmt::format("{}/h_mean_layer{}_module{}", dqmDir_, i_layer, i_module));
0306
0307 hists_.h_drift_depth_[i_index]->divide(
0308 hists_.h_drift_depth_adc_[i_index], hists_.h_drift_depth_noadc_[i_index], 1., 1., "");
0309 }
0310 }
0311
0312
0313 for (int i = 0; i < (int)hists_.BPixnewDetIds_.size(); i++) {
0314 int new_index = i + 1 + hists_.nModules_[hists_.nlay - 1] + (hists_.nlay - 1) * hists_.nModules_[hists_.nlay - 1];
0315
0316 hists_.h_drift_depth_[new_index] = iGetter.get(
0317 fmt::format("{}/h_BPixnew_drift_depth_{}", dqmDir_ + "/BPix/NewModules", hists_.BPixnewmodulename_[i]));
0318
0319 if (hists_.h_drift_depth_[new_index] == nullptr) {
0320 edm::LogError("SiPixelLorentzAnglePCLHarvester")
0321 << "Failed to retrieve electron drift over depth for new module " << hists_.BPixnewmodulename_[i] << ".";
0322 continue;
0323 }
0324
0325 hists_.h_drift_depth_adc_[new_index] = iGetter.get(
0326 fmt::format("{}/h_BPixnew_drift_depth_adc_{}", dqmDir_ + "/BPix/NewModules", hists_.BPixnewmodulename_[i]));
0327
0328 hists_.h_drift_depth_adc2_[new_index] = iGetter.get(
0329 fmt::format("{}/h_BPixnew_drift_depth_adc2_{}", dqmDir_ + "/BPix/NewModules", hists_.BPixnewmodulename_[i]));
0330
0331 hists_.h_drift_depth_noadc_[new_index] = iGetter.get(
0332 fmt::format("{}/h_BPixnew_drift_depth_noadc_{}", dqmDir_ + "/BPix/NewModules", hists_.BPixnewmodulename_[i]));
0333
0334 hists_.h_mean_[new_index] = iGetter.get(fmt::format("{}/h_BPixnew_mean_{}", dqmDir_, hists_.BPixnewmodulename_[i]));
0335
0336 hists_.h_drift_depth_[new_index]->divide(
0337 hists_.h_drift_depth_adc_[new_index], hists_.h_drift_depth_noadc_[new_index], 1., 1., "");
0338 }
0339
0340 hists_.h_bySectOccupancy_ = iGetter.get(fmt::format("{}/h_bySectorOccupancy", dqmDir_ + "/SectorMonitoring"));
0341 if (hists_.h_bySectOccupancy_ == nullptr) {
0342 edm::LogError("SiPixelLorentzAnglePCLHarvester") << "Failed to retrieve the hit on track occupancy.";
0343 return;
0344 }
0345
0346 int hist_drift_;
0347 int hist_depth_;
0348 double min_drift_;
0349 double max_drift_;
0350
0351 if (hists_.h_drift_depth_adc_[1] != nullptr) {
0352 hist_drift_ = hists_.h_drift_depth_adc_[1]->getNbinsX();
0353 hist_depth_ = hists_.h_drift_depth_adc_[1]->getNbinsY();
0354 min_drift_ = hists_.h_drift_depth_adc_[1]->getAxisMin(1);
0355 max_drift_ = hists_.h_drift_depth_adc_[1]->getAxisMax(1);
0356 } else {
0357 hist_drift_ = 100;
0358 hist_depth_ = 50;
0359 min_drift_ = -500.;
0360 max_drift_ = 500.;
0361 }
0362
0363 iBooker.setCurrentFolder(fmt::format("{}Harvesting", dqmDir_));
0364 MonitorElement* h_drift_depth_adc_slice_ =
0365 iBooker.book1D("h_drift_depth_adc_slice", "slice of adc histogram", hist_drift_, min_drift_, max_drift_);
0366
0367
0368 MonitorElement* h_diffLA = iBooker.book1D(
0369 "h_diffLA", "difference in #mu_{H}; #Delta #mu_{H}/#mu_{H} (old-new)/old [%];n. modules", 100, -150, 150);
0370
0371
0372 const auto& maxSect = hists_.h_bySectOccupancy_->getNbinsX();
0373 const double lo = -0.5;
0374 const double hi = maxSect - 0.5;
0375
0376
0377 iBooker.setCurrentFolder(fmt::format("{}Harvesting/SectorMonitoring", dqmDir_));
0378 std::string repText = "%s tan#theta_{LA}/B by sector;pixel sector;%s tan(#theta_{LA})/B [1/T]";
0379 hists_.h_bySectMeasLA_ =
0380 iBooker.book1D("h_LAbySector_Measured", fmt::sprintf(repText, "measured", "measured"), maxSect, lo, hi);
0381 hists_.h_bySectSetLA_ =
0382 iBooker.book1D("h_LAbySector_Accepted", fmt::sprintf(repText, "accepted", "accepted"), maxSect, lo, hi);
0383 hists_.h_bySectRejectLA_ =
0384 iBooker.book1D("h_LAbySector_Rejected", fmt::sprintf(repText, "rejected", "rejected"), maxSect, lo, hi);
0385 hists_.h_bySectLA_ = iBooker.book1D("h_LAbySector", fmt::sprintf(repText, "payload", "payload"), maxSect, lo, hi);
0386 hists_.h_bySectDeltaLA_ =
0387 iBooker.book1D("h_deltaLAbySector", fmt::sprintf(repText, "#Delta", "#Delta"), maxSect, lo, hi);
0388 hists_.h_bySectChi2_ =
0389 iBooker.book1D("h_bySectorChi2", "Fit #chi^{2}/ndf by sector;pixel sector; fit #chi^{2}/ndf", maxSect, lo, hi);
0390
0391 hists_.h_bySectFitStatus_ =
0392 iBooker.book1D("h_bySectFitStatus", "Fit Status by sector;pixel sector; fit status", maxSect, lo, hi);
0393
0394 hists_.h_bySectCovMatrixStatus_ = iBooker.book1D(
0395 "h_bySectorCovMatrixStatus", "Fit Covariance Matrix Status by sector;pixel sector; fit status", maxSect, lo, hi);
0396 hists_.h_bySectDriftError_ =
0397 iBooker.book1D("h_bySectorDriftError",
0398 "square error on the measured drift at half-width by sector;pixel sector;#Delta d^{2}(t/2)",
0399 maxSect,
0400 lo,
0401 hi);
0402
0403
0404 std::vector<std::string> labels = {"#chi^{2}!=0", "#chi^{2} cut", "covStat>=2", "n. entries cut"};
0405 hists_.h_bySectFitQuality_ = iBooker.book2D("h_bySectFitQuality",
0406 "Fit quality by sector;pixel sector;quality cut",
0407 maxSect,
0408 lo,
0409 hi,
0410 labels.size(),
0411 -0.5,
0412 labels.size() - 0.5);
0413
0414
0415 for (int bin = 1; bin <= maxSect; bin++) {
0416 const auto& binName = hists_.h_bySectOccupancy_->getTH1()->GetXaxis()->GetBinLabel(bin);
0417 hists_.h_bySectMeasLA_->setBinLabel(bin, binName);
0418 hists_.h_bySectSetLA_->setBinLabel(bin, binName);
0419 hists_.h_bySectRejectLA_->setBinLabel(bin, binName);
0420 hists_.h_bySectLA_->setBinLabel(bin, binName);
0421 hists_.h_bySectDeltaLA_->setBinLabel(bin, binName);
0422 hists_.h_bySectChi2_->setBinLabel(bin, binName);
0423 hists_.h_bySectFitStatus_->setBinLabel(bin, binName);
0424 hists_.h_bySectCovMatrixStatus_->setBinLabel(bin, binName);
0425 hists_.h_bySectDriftError_->setBinLabel(bin, binName);
0426 hists_.h_bySectFitQuality_->setBinLabel(bin, binName, 1);
0427 }
0428
0429 for (unsigned int binY = 1; binY <= labels.size(); binY++) {
0430 hists_.h_bySectFitQuality_->setBinLabel(binY, labels[binY - 1], 2);
0431 }
0432
0433
0434 iBooker.setCurrentFolder(fmt::format("{}Harvesting/LorentzAngleMaps", dqmDir_));
0435 for (int i = 0; i < hists_.nlay; i++) {
0436 std::string repName = "h2_byLayerLA_%i";
0437 std::string repText = "BPix Layer %i tan#theta_{LA}/B;module number;ladder number;tan#theta_{LA}/B [1/T]";
0438
0439 hists_.h2_byLayerLA_.emplace_back(iBooker.book2D(fmt::sprintf(repName, i + 1),
0440 fmt::sprintf(repText, i + 1),
0441 hists_.nModules_[i],
0442 0.5,
0443 hists_.nModules_[i] + 0.5,
0444 hists_.nLadders_[i],
0445 0.5,
0446 hists_.nLadders_[i] + 0.5));
0447
0448 repName = "h2_byLayerDiff_%i";
0449 repText = "BPix Layer %i #Delta#mu_{H}/#mu_{H};module number;ladder number;#Delta#mu_{H}/#mu_{H} [%%]";
0450 hists_.h2_byLayerDiff_.emplace_back(iBooker.book2D(fmt::sprintf(repName, i + 1),
0451 fmt::sprintf(repText, i + 1),
0452 hists_.nModules_[i],
0453 0.5,
0454 hists_.nModules_[i] + 0.5,
0455 hists_.nLadders_[i],
0456 0.5,
0457 hists_.nLadders_[i] + 0.5));
0458 }
0459
0460
0461 edm::LogPrint("LorentzAngle") << "module" << "\t" << "layer" << "\t"
0462 << "offset" << "\t" << "e0" << "\t"
0463 << "slope" << "\t" << "e1" << "\t"
0464 << "rel.err" << "\t" << "pull" << "\t"
0465 << "p2" << "\t" << "e2" << "\t"
0466 << "p3" << "\t" << "e3" << "\t"
0467 << "p4" << "\t" << "e4" << "\t"
0468 << "p5" << "\t" << "e5" << "\t"
0469 << "chi2" << "\t" << "prob" << "\t"
0470 << "newDetId" << "\t" << "tan(LA)" << "\t"
0471 << "Error(LA)" ;
0472
0473
0474
0475 std::shared_ptr<SiPixelLorentzAngle> LorentzAngle = std::make_shared<SiPixelLorentzAngle>();
0476
0477
0478 double p1_simul_newmodule = 0.294044;
0479 double p1_simul[hists_.nlay + 1][hists_.nModules_[hists_.nlay - 1]];
0480 for (int i_layer = 1; i_layer <= hists_.nlay; i_layer++) {
0481 for (int i_module = 1; i_module <= hists_.nModules_[i_layer - 1]; i_module++) {
0482 if (i_layer == 1)
0483 p1_simul[i_layer - 1][i_module - 1] = 0.436848;
0484 else if (i_layer == 2)
0485 p1_simul[i_layer - 1][i_module - 1] = 0.25802;
0486 else if (i_layer == 3 && i_module <= 4)
0487 p1_simul[i_layer - 1][i_module - 1] = 0.29374;
0488 else if (i_layer == 3 && i_module >= 5)
0489 p1_simul[i_layer - 1][i_module - 1] = 0.31084;
0490 else if (i_layer == 4 && i_module <= 4)
0491 p1_simul[i_layer - 1][i_module - 1] = 0.29944;
0492 else
0493 p1_simul[i_layer - 1][i_module - 1] = 0.31426;
0494 }
0495 }
0496
0497 for (int i_module = 1; i_module <= hists_.nModules_[hists_.nlay - 1]; i_module++) {
0498 p1_simul[hists_.nlay][i_module - 1] = p1_simul_newmodule;
0499 }
0500
0501
0502 for (int j = 0; j < (int)hists_.BPixnewDetIds_.size(); j++) {
0503
0504 int new_index = j + 1 + hists_.nModules_[hists_.nlay - 1] + (hists_.nlay - 1) * hists_.nModules_[hists_.nlay - 1];
0505 if (hists_.h_drift_depth_adc_[new_index] == nullptr)
0506 continue;
0507 for (int i = 1; i <= hist_depth_; i++) {
0508 findMean(h_drift_depth_adc_slice_, i, new_index);
0509 }
0510
0511
0512 const auto& res = fitAndStore(LorentzAngle, new_index, hists_.BPixnewLayer_[j], hists_.BPixnewModule_[j]);
0513
0514 edm::LogPrint("SiPixelLorentzAngle") << std::setprecision(4) << hists_.BPixnewModule_[j] << "\t"
0515 << hists_.BPixnewLayer_[j] << "\t" << res.p0 << "\t" << res.e0 << "\t"
0516 << res.p1 << std::setprecision(3) << "\t" << res.e1 << "\t"
0517 << res.e1 / res.p1 * 100. << "\t"
0518 << (res.p1 - p1_simul[hists_.nlay][0]) / res.e1 << "\t" << res.p2 << "\t"
0519 << res.e2 << "\t" << res.p3 << "\t" << res.e3 << "\t" << res.p4 << "\t"
0520 << res.e4 << "\t" << res.p5 << "\t" << res.e5 << "\t" << res.chi2 << "\t"
0521 << res.prob << "\t" << hists_.BPixnewDetIds_[j] << "\t" << res.tan_LA << "\t"
0522 << res.error_LA;
0523 }
0524
0525
0526 for (int i_layer = 1; i_layer <= hists_.nlay; i_layer++) {
0527 for (int i_module = 1; i_module <= hists_.nModules_[i_layer - 1]; i_module++) {
0528 int i_index = i_module + (i_layer - 1) * hists_.nModules_[i_layer - 1];
0529 if (hists_.h_drift_depth_adc_[i_index] == nullptr)
0530 continue;
0531
0532 for (int i = 1; i <= hist_depth_; i++) {
0533 findMean(h_drift_depth_adc_slice_, i, i_index);
0534 }
0535
0536
0537 const auto& res = fitAndStore(LorentzAngle, i_index, i_layer, i_module);
0538
0539 edm::LogPrint("SiPixelLorentzAngle")
0540 << std::setprecision(4) << i_module << "\t" << i_layer << "\t" << res.p0 << "\t" << res.e0 << "\t" << res.p1
0541 << std::setprecision(3) << "\t" << res.e1 << "\t" << res.e1 / res.p1 * 100. << "\t"
0542 << (res.p1 - p1_simul[i_layer - 1][i_module - 1]) / res.e1 << "\t" << res.p2 << "\t" << res.e2 << "\t"
0543 << res.p3 << "\t" << res.e3 << "\t" << res.p4 << "\t" << res.e4 << "\t" << res.p5 << "\t" << res.e5 << "\t"
0544 << res.chi2 << "\t" << res.prob << "\t"
0545 << "null"
0546 << "\t" << res.tan_LA << "\t" << res.error_LA;
0547 }
0548 }
0549
0550
0551 const auto& currentLAMap = currentLorentzAngle_->getLorentzAngles();
0552 const auto& newLAMap = LorentzAngle->getLorentzAngles();
0553 std::vector<unsigned int> currentLADets;
0554 std::vector<unsigned int> newLADets;
0555
0556 std::transform(currentLAMap.begin(),
0557 currentLAMap.end(),
0558 std::back_inserter(currentLADets),
0559 [](const std::map<unsigned int, float>::value_type& pair) { return pair.first; });
0560
0561 std::transform(newLAMap.begin(),
0562 newLAMap.end(),
0563 std::back_inserter(newLADets),
0564 [](const std::map<unsigned int, float>::value_type& pair) { return pair.first; });
0565
0566 std::vector<unsigned int> notCommon;
0567 std::set_symmetric_difference(
0568 currentLADets.begin(), currentLADets.end(), newLADets.begin(), newLADets.end(), std::back_inserter(notCommon));
0569
0570 for (const auto& id : notCommon) {
0571 float fPixLorentzAnglePerTesla_ = currentLorentzAngle_->getLorentzAngle(id);
0572 if (!LorentzAngle->putLorentzAngle(id, fPixLorentzAnglePerTesla_)) {
0573 edm::LogError("SiPixelLorentzAnglePCLHarvester")
0574 << "[SiPixelLorentzAnglePCLHarvester::dqmEndJob] filling rest of payload: detid already exists";
0575 }
0576 }
0577
0578 for (const auto& id : newLADets) {
0579 float deltaMuHoverMuH = (currentLorentzAngle_->getLorentzAngle(id) - LorentzAngle->getLorentzAngle(id)) /
0580 currentLorentzAngle_->getLorentzAngle(id);
0581 h_diffLA->Fill(deltaMuHoverMuH * 100.f);
0582 }
0583
0584 bool isPayloadChanged{false};
0585
0586 for (const auto& [id, value] : LorentzAngle->getLorentzAngles()) {
0587 DetId ID = DetId(id);
0588 if (ID.subdetId() == PixelSubdetector::PixelBarrel) {
0589 const auto& layer = theTrackerTopology_->pxbLayer(id);
0590 const auto& ladder = theTrackerTopology_->pxbLadder(id);
0591 const auto& module = theTrackerTopology_->pxbModule(id);
0592 hists_.h2_byLayerLA_[layer - 1]->setBinContent(module, ladder, value);
0593
0594 float deltaMuHoverMuH =
0595 (currentLorentzAngle_->getLorentzAngle(id) - value) / currentLorentzAngle_->getLorentzAngle(id);
0596 hists_.h2_byLayerDiff_[layer - 1]->setBinContent(module, ladder, deltaMuHoverMuH * 100.f);
0597
0598 if (!isPayloadChanged && (deltaMuHoverMuH != 0.f))
0599 isPayloadChanged = true;
0600 }
0601 }
0602
0603 if (isPayloadChanged) {
0604
0605 edm::Service<cond::service::PoolDBOutputService> mydbservice;
0606 if (mydbservice.isAvailable()) {
0607 try {
0608 mydbservice->writeOneIOV(*LorentzAngle, mydbservice->currentTime(), recordName_);
0609 } catch (const cond::Exception& er) {
0610 edm::LogError("SiPixelLorentzAngleDB") << er.what();
0611 } catch (const std::exception& er) {
0612 edm::LogError("SiPixelLorentzAngleDB") << "caught std::exception " << er.what();
0613 }
0614 } else {
0615 edm::LogError("SiPixelLorentzAngleDB") << "Service is unavailable";
0616 }
0617 } else {
0618 edm::LogPrint("SiPixelLorentzAngleDB") << __PRETTY_FUNCTION__ << " there is no new valid measurement to append! ";
0619 }
0620 }
0621
0622
0623 void SiPixelLorentzAnglePCLHarvester::findMean(MonitorElement* h_drift_depth_adc_slice_, int i, int i_ring) {
0624 double nentries = 0;
0625 h_drift_depth_adc_slice_->Reset();
0626 int hist_drift_ = h_drift_depth_adc_slice_->getNbinsX();
0627
0628
0629
0630 for (int j = 1; j <= hist_drift_; j++) {
0631 if (hists_.h_drift_depth_noadc_[i_ring]->getBinContent(j, i) >= 1) {
0632 double adc_error2 = (hists_.h_drift_depth_adc2_[i_ring]->getBinContent(j, i) -
0633 hists_.h_drift_depth_adc_[i_ring]->getBinContent(j, i) *
0634 hists_.h_drift_depth_adc_[i_ring]->getBinContent(j, i) /
0635 hists_.h_drift_depth_noadc_[i_ring]->getBinContent(j, i)) /
0636 hists_.h_drift_depth_noadc_[i_ring]->getBinContent(j, i);
0637
0638 hists_.h_drift_depth_adc_[i_ring]->setBinError(j, i, sqrt(adc_error2));
0639 double error2 = adc_error2 / (hists_.h_drift_depth_noadc_[i_ring]->getBinContent(j, i) - 1.);
0640 hists_.h_drift_depth_[i_ring]->setBinError(j, i, sqrt(error2));
0641 } else {
0642 hists_.h_drift_depth_[i_ring]->setBinError(j, i, 0);
0643 hists_.h_drift_depth_adc_[i_ring]->setBinError(j, i, 0);
0644 }
0645 h_drift_depth_adc_slice_->setBinContent(j, hists_.h_drift_depth_adc_[i_ring]->getBinContent(j, i));
0646 h_drift_depth_adc_slice_->setBinError(j, hists_.h_drift_depth_adc_[i_ring]->getBinError(j, i));
0647 nentries += hists_.h_drift_depth_noadc_[i_ring]->getBinContent(j, i);
0648 }
0649
0650 double mean = h_drift_depth_adc_slice_->getMean(1);
0651 double error = 0;
0652 if (nentries != 0) {
0653 error = h_drift_depth_adc_slice_->getRMS(1) / std::sqrt(nentries);
0654 }
0655 hists_.h_mean_[i_ring]->setBinContent(i, mean);
0656 hists_.h_mean_[i_ring]->setBinError(i, error);
0657
0658 h_drift_depth_adc_slice_->Reset();
0659 }
0660
0661
0662 SiPixelLAHarvest::fitResults SiPixelLorentzAnglePCLHarvester::fitAndStore(
0663 std::shared_ptr<SiPixelLorentzAngle> theLAPayload, int i_index, int i_layer, int i_module) {
0664
0665 SiPixelLAHarvest::fitResults res;
0666
0667 double half_width = width_ * siPixelLACalibration::cmToum / 2.f;
0668
0669 if (doChebyshevFit_) {
0670 const int npar = order_ + 1;
0671 auto cheb = std::make_unique<siPixelLACalibration::Chebyshev>(order_, theFitRange_.first, theFitRange_.second);
0672 f1_ = std::make_unique<TF1>("f1", cheb.release(), theFitRange_.first, theFitRange_.second, npar, "Chebyshev");
0673 } else {
0674 f1_ = std::make_unique<TF1>("f1",
0675 "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x + [5]*x*x*x*x*x",
0676 theFitRange_.first,
0677 theFitRange_.second);
0678 }
0679
0680
0681
0682
0683
0684 assert(strcmp(f1_->GetName(), "f1") == 0);
0685
0686 f1_->SetParName(0, "offset");
0687 f1_->SetParName(1, "tan#theta_{LA}");
0688 f1_->SetParName(2, "quad term");
0689 f1_->SetParName(3, "cubic term");
0690 f1_->SetParName(4, "quartic term");
0691 f1_->SetParName(5, "quintic term");
0692
0693 f1_->SetParameter(0, 0);
0694 f1_->SetParError(0, 0);
0695 f1_->SetParameter(1, 0.4);
0696 f1_->SetParError(1, 0);
0697 f1_->SetParameter(2, 0.0);
0698 f1_->SetParError(2, 0);
0699 f1_->SetParameter(3, 0.0);
0700 f1_->SetParError(3, 0);
0701 f1_->SetParameter(4, 0.0);
0702 f1_->SetParError(4, 0);
0703 f1_->SetParameter(5, 0.0);
0704 f1_->SetParError(5, 0);
0705 f1_->SetChisquare(0);
0706
0707 TFitResultPtr result = hists_.h_mean_[i_index]->getTH1()->Fit(f1_.get(), "ERQS");
0708 if (result.Get()) {
0709 res.fitStatus = result->Status();
0710 edm::LogInfo("SiPixelLorentzAnglePCLHarvester") << "Fit status: " << res.fitStatus;
0711 } else {
0712 edm::LogError("SiPixelLorentzAnglePCLHarvester") << "Could not retrieve the fit status! Setting it to 0 by default";
0713 res.fitStatus = -1;
0714 }
0715
0716 if (res.fitStatus >= 0) {
0717 res.covMatrixStatus = result->CovMatrixStatus();
0718 } else {
0719 res.covMatrixStatus = SiPixelLAHarvest::kUndefined;
0720 }
0721
0722
0723 float dSq{0.};
0724 float cov00{0.};
0725 if (!doChebyshevFit_) {
0726 if (res.covMatrixStatus > SiPixelLAHarvest::kNotCalculated) {
0727 for (int k = 0; k < order_; k++) {
0728 for (int l = 0; l < order_; l++) {
0729 dSq += (std::pow(half_width, k) * std::pow(half_width, l) * result->CovMatrix(k, l));
0730 }
0731 }
0732 cov00 = result->CovMatrix(0, 0);
0733 }
0734 }
0735
0736 res.dSq = dSq;
0737
0738 res.p0 = f1_->GetParameter(0);
0739 res.e0 = f1_->GetParError(0);
0740 res.p1 = f1_->GetParameter(1);
0741 res.e1 = f1_->GetParError(1);
0742 res.p2 = f1_->GetParameter(2);
0743 res.e2 = f1_->GetParError(2);
0744 res.p3 = f1_->GetParameter(3);
0745 res.e3 = f1_->GetParError(3);
0746 res.p4 = f1_->GetParameter(4);
0747 res.e4 = f1_->GetParError(4);
0748 res.p5 = f1_->GetParameter(5);
0749 res.e5 = f1_->GetParError(5);
0750 res.chi2 = f1_->GetChisquare();
0751 res.ndf = f1_->GetNDF();
0752 res.prob = f1_->GetProb();
0753 res.redChi2 = res.ndf > 0. ? res.chi2 / res.ndf : 0.;
0754
0755 double f1_halfwidth = res.p0 + res.p1 * half_width + res.p2 * pow(half_width, 2) + res.p3 * pow(half_width, 3) +
0756 res.p4 * pow(half_width, 4) + res.p5 * pow(half_width, 5);
0757
0758 double f1_zerowidth = res.p0;
0759
0760
0761 res.tan_LA = (f1_halfwidth - f1_zerowidth) / half_width;
0762 double errsq_LA =
0763 (pow(res.e1, 2) + pow((half_width * res.e2), 2) + pow((half_width * half_width * res.e3), 2) +
0764 pow((half_width * half_width * half_width * res.e4), 2) +
0765 pow((half_width * half_width * half_width * half_width * res.e5), 2));
0766
0767 res.error_LA = doChebyshevFit_ ? sqrt(errsq_LA) : sqrt(res.dSq + cov00) / half_width;
0768
0769 hists_.h_bySectMeasLA_->setBinContent(i_index, (res.tan_LA / theMagField_));
0770 hists_.h_bySectMeasLA_->setBinError(i_index, (res.error_LA / theMagField_));
0771 hists_.h_bySectChi2_->setBinContent(i_index, res.redChi2);
0772 hists_.h_bySectChi2_->setBinError(i_index, 0.);
0773
0774 hists_.h_bySectFitStatus_->setBinContent(i_index, res.fitStatus);
0775 hists_.h_bySectFitStatus_->setBinError(i_index, 0.);
0776 hists_.h_bySectCovMatrixStatus_->setBinContent(i_index, res.covMatrixStatus);
0777 hists_.h_bySectCovMatrixStatus_->setBinError(i_index, 0.);
0778 hists_.h_bySectDriftError_->setBinContent(i_index, res.dSq);
0779 hists_.h_bySectDriftError_->setBinError(i_index, 0.);
0780
0781 res.nentries = hists_.h_bySectOccupancy_->getBinContent(i_index);
0782
0783 bool isNew = (i_index > hists_.nlay * hists_.nModules_[hists_.nlay - 1]);
0784 int shiftIdx = i_index - hists_.nlay * hists_.nModules_[hists_.nlay - 1] - 1;
0785
0786 LogDebug("SiPixelLorentzAnglePCLHarvester")
0787 << " isNew: " << isNew << " i_index: " << i_index << " shift index: " << shiftIdx;
0788
0789 const auto& detIdsToFill =
0790 isNew ? std::vector<unsigned int>({hists_.BPixnewDetIds_[shiftIdx]}) : hists_.detIdsList[i_index];
0791
0792 LogDebug("SiPixelLorentzAnglePCLHarvester")
0793 << "index: " << i_index << " i_module: " << i_module << " i_layer: " << i_layer;
0794 for (const auto& id : detIdsToFill) {
0795 LogDebug("SiPixelLorentzAnglePCLHarvester") << id << ",";
0796 }
0797
0798
0799 hists_.h_bySectSetLA_->setBinError(i_index, 0.);
0800 hists_.h_bySectRejectLA_->setBinError(i_index, 0.);
0801 hists_.h_bySectLA_->setBinError(i_index, 0.);
0802 hists_.h_bySectDeltaLA_->setBinError(i_index, 0.);
0803
0804 float LorentzAnglePerTesla_;
0805 float currentLA = currentLorentzAngle_->getLorentzAngle(detIdsToFill.front());
0806
0807
0808 bool goodFit = isFitGood(res);
0809 for (unsigned int i_bin = 0; i_bin < res.quality.size(); i_bin++) {
0810 if (res.quality[i_bin]) {
0811 hists_.h_bySectFitQuality_->setBinContent(i_index, i_bin + 1, 1.);
0812 } else {
0813 hists_.h_bySectFitQuality_->setBinContent(i_index, i_bin + 1, 0.);
0814 }
0815 }
0816
0817
0818 if (goodFit) {
0819 LorentzAnglePerTesla_ = res.tan_LA / theMagField_;
0820
0821 hists_.h_bySectSetLA_->setBinContent(i_index, LorentzAnglePerTesla_);
0822 hists_.h_bySectRejectLA_->setBinContent(i_index, 0.);
0823 hists_.h_bySectLA_->setBinContent(i_index, LorentzAnglePerTesla_);
0824
0825 const auto& deltaLA = (LorentzAnglePerTesla_ - currentLA);
0826 hists_.h_bySectDeltaLA_->setBinContent(i_index, deltaLA);
0827
0828 for (const auto& id : detIdsToFill) {
0829 if (!theLAPayload->putLorentzAngle(id, LorentzAnglePerTesla_)) {
0830 edm::LogError("SiPixelLorentzAnglePCLHarvester") << "[SiPixelLorentzAnglePCLHarvester::fitAndStore]: detid ("
0831 << i_layer << "," << i_module << ") already exists";
0832 }
0833 }
0834 } else {
0835
0836 hists_.h_bySectSetLA_->setBinContent(i_index, 0.);
0837 hists_.h_bySectRejectLA_->setBinContent(i_index, (res.tan_LA / theMagField_));
0838 hists_.h_bySectLA_->setBinContent(i_index, currentLA);
0839 hists_.h_bySectDeltaLA_->setBinContent(i_index, 0.);
0840
0841 for (const auto& id : detIdsToFill) {
0842 LorentzAnglePerTesla_ = currentLorentzAngle_->getLorentzAngle(id);
0843 if (!theLAPayload->putLorentzAngle(id, LorentzAnglePerTesla_)) {
0844 edm::LogError("SiPixelLorentzAnglePCLHarvester") << "[SiPixelLorentzAnglePCLHarvester::fitAndStore]: detid ("
0845 << i_layer << "," << i_module << ") already exists";
0846 }
0847 }
0848 }
0849
0850 return res;
0851 }
0852
0853
0854 bool SiPixelLorentzAnglePCLHarvester::isFitGood(SiPixelLAHarvest::fitResults& res) {
0855
0856 const bool notZeroCut = (res.redChi2 != 0.);
0857
0858 const bool chi2Cut = doChebyshevFit_ ? (res.redChi2 < fitChi2Cut_) : true;
0859
0860 const bool covMatrixStatusCut = (res.covMatrixStatus >= SiPixelLAHarvest::kMadePosDef);
0861
0862 const bool nentriesCut = (res.nentries > minHitsCut_);
0863
0864 if (notZeroCut) {
0865 res.setQualityCutBit(SiPixelLAHarvest::kZeroChi2);
0866 }
0867 if (chi2Cut) {
0868 res.setQualityCutBit(SiPixelLAHarvest::kChi2Cut);
0869 }
0870 if (covMatrixStatusCut) {
0871 res.setQualityCutBit(SiPixelLAHarvest::kCovStatus);
0872 }
0873 if (nentriesCut) {
0874 res.setQualityCutBit(SiPixelLAHarvest::kNentries);
0875 }
0876
0877
0878 return (res.quality.all());
0879 }
0880
0881
0882 void SiPixelLorentzAnglePCLHarvester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0883 edm::ParameterSetDescription desc;
0884 desc.setComment("Harvester module of the SiPixel Lorentz Angle PCL monitoring workflow");
0885 desc.add<std::vector<std::string>>("newmodulelist", {})->setComment("the list of DetIds for new sensors");
0886 desc.add<std::string>("dqmDir", "AlCaReco/SiPixelLorentzAngle")->setComment("the directory of PCL Worker output");
0887 desc.add<bool>("doChebyshevFit", false)->setComment("use Chebyshev polynomials for the dript vs depth fit");
0888 desc.add<int>("order", 5)->setComment("order of the Chebyshev polynomial used for the fit");
0889 desc.add<double>("fitChi2Cut", 20.)->setComment("cut on fit chi2/ndof to accept measurement");
0890 desc.add<std::vector<double>>("fitRange", {5., 280.})->setComment("range of depths to perform the LA fit");
0891 desc.add<int>("minHitsCut", 10000)->setComment("cut on minimum number of on-track hits to accept measurement");
0892 desc.add<std::string>("record", "SiPixelLorentzAngleRcd")->setComment("target DB record");
0893 descriptions.addWithDefaultLabel(desc);
0894 }
0895
0896 DEFINE_FWK_MODULE(SiPixelLorentzAnglePCLHarvester);
