File indexing completed on 2024-04-06 11:59:35
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023 #include <fmt/format.h>
0024 #include <fmt/printf.h>
0025 #include <fstream>
0026
0027
0028 #include "CalibTracker/SiPixelLorentzAngle/interface/SiPixelLorentzAngleCalibrationStruct.h"
0029 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0030 #include "CondFormats/DataRecord/interface/SiPixelLorentzAngleRcd.h"
0031 #include "CondFormats/SiPixelObjects/interface/SiPixelLorentzAngle.h"
0032 #include "DQMServices/Core/interface/DQMEDHarvester.h"
0033 #include "DataFormats/TrackerCommon/interface/PixelBarrelName.h"
0034 #include "DataFormats/TrackerCommon/interface/PixelEndcapName.h"
0035 #include "FWCore/Framework/interface/EventSetup.h"
0036 #include "FWCore/Framework/interface/Frameworkfwd.h"
0037 #include "FWCore/Framework/interface/MakerMacros.h"
0038 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0039 #include "FWCore/ServiceRegistry/interface/Service.h"
0040 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0041 #include "Geometry/TrackerGeometryBuilder/interface/PixelTopologyMap.h"
0042 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0043 #include "MagneticField/Engine/interface/MagneticField.h"
0044 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0045
0046
0047 #include "TFitResult.h"
0048 #include "TVirtualFitter.h"
0049 #include "TH1D.h"
0050
0051 namespace SiPixelLAHarvestMCS {
0052
0053 enum fitStatus { kFitNotPerformed = -2, kNoFitResult = -1, kFitConverged = 0, kFitFailed = 1 };
0054
0055 Double_t MCSFitFunction(Double_t* x, Double_t* par) {
0056 Double_t arg;
0057
0058 if (x[0] < par[3]) {
0059 arg = par[1] * par[1] + par[2] * par[2] * (x[0] - par[3]) * (x[0] - par[3]);
0060 } else {
0061 arg = par[1] * par[1] + par[4] * par[4] * (x[0] - par[3]) * (x[0] - par[3]);
0062 }
0063 Double_t fitval = par[0] + sqrt(arg);
0064 return fitval;
0065 }
0066
0067 struct FPixMuH {
0068 double bfield[3];
0069 double shiftx;
0070 double shifty;
0071 double shiftx_err;
0072 double shifty_err;
0073 };
0074
0075 void fcn_func(Int_t& npar, Double_t* deriv, Double_t& f, Double_t* par, Int_t flag);
0076
0077 class FitFPixMuH : public TObject {
0078 public:
0079 FitFPixMuH() : muH(0), muHErr(0) {}
0080
0081
0082 ~FitFPixMuH() override = default;
0083
0084 friend void fcn_func(Int_t& npar, Double_t* deriv, Double_t& f, Double_t* par, Int_t flag);
0085
0086 void add(const FPixMuH& fmh) { cmvar.push_back(fmh); }
0087 int fit(double initVal) {
0088 minuit = TVirtualFitter::Fitter(this, 1);
0089 minuit->SetFCN(fcn_func);
0090 double tanlorpertesla = initVal;
0091 minuit->SetParameter(0, "muH", tanlorpertesla, 0.1, 0., 0.);
0092
0093 double arglist[100];
0094 arglist[0] = 3.;
0095 minuit->ExecuteCommand("SET PRINT", arglist, 0);
0096 double up = 1.0;
0097 minuit->SetErrorDef(up);
0098 arglist[0] = 100000.;
0099 int status = minuit->ExecuteCommand("MIGRAD", arglist, 0);
0100 muH = minuit->GetParameter(0);
0101 muHErr = minuit->GetParError(0);
0102
0103 return status;
0104 }
0105
0106 double getMuH() { return muH; }
0107 double getMuHErr() { return muHErr; }
0108 unsigned int size() { return cmvar.size(); }
0109
0110 private:
0111 TVirtualFitter* minuit;
0112 std::vector<FPixMuH> cmvar;
0113 double muH;
0114 double muHErr;
0115
0116 double calcChi2(double par_0) {
0117 double tanlorpertesla = par_0;
0118 double tlpt2 = tanlorpertesla * tanlorpertesla;
0119 double v[3], xshift, yshift;
0120
0121 int n = cmvar.size();
0122
0123 double chi2 = 0.0;
0124 for (int i = 0; i < n; i++) {
0125 v[0] = -(tanlorpertesla * cmvar[i].bfield[1] + tlpt2 * cmvar[i].bfield[0] * cmvar[i].bfield[2]);
0126 v[1] = tanlorpertesla * cmvar[i].bfield[0] - tlpt2 * cmvar[i].bfield[1] * cmvar[i].bfield[2];
0127 v[2] = -(1. + tlpt2 * cmvar[i].bfield[2] * cmvar[i].bfield[2]);
0128
0129 xshift = v[0] / v[2];
0130 yshift = v[1] / v[2];
0131
0132 chi2 += (xshift - cmvar[i].shiftx) * (xshift - cmvar[i].shiftx) / cmvar[i].shiftx_err / cmvar[i].shiftx_err +
0133 (yshift - cmvar[i].shifty) * (yshift - cmvar[i].shifty) / cmvar[i].shifty_err / cmvar[i].shifty_err;
0134 }
0135 return chi2;
0136 }
0137 };
0138
0139 void fcn_func(Int_t& npar, Double_t* deriv, Double_t& f, Double_t* par, Int_t flag) {
0140 f = ((dynamic_cast<FitFPixMuH*>((TVirtualFitter::GetFitter())->GetObjectFit()))->calcChi2(par[0]));
0141 }
0142 }
0143
0144
0145 class SiPixelLorentzAnglePCLHarvesterMCS : public DQMEDHarvester {
0146 public:
0147 SiPixelLorentzAnglePCLHarvesterMCS(const edm::ParameterSet&);
0148 ~SiPixelLorentzAnglePCLHarvesterMCS() override = default;
0149 using FPixCotAngleFitResults = std::unordered_map<uint32_t, std::pair<double, double>>;
0150 using FpixMuHResults = std::unordered_map<std::string, std::pair<double, double>>;
0151 using FitParametersInitValuesMuHFitMap = std::unordered_map<std::string, double>;
0152
0153 void beginRun(const edm::Run&, const edm::EventSetup&) override;
0154
0155 static void fillDescriptions(edm::ConfigurationDescriptions&);
0156
0157 private:
0158 void dqmEndJob(DQMStore::IBooker&, DQMStore::IGetter&) override;
0159 void findMean(dqm::reco::MonitorElement* h_2D, dqm::reco::MonitorElement* h_mean, TH1D* h_slice);
0160 int getIndex(bool isBetaAngle, int r, int p, int s);
0161
0162
0163 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomEsToken_;
0164 const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoEsTokenBR_, topoEsTokenER_;
0165 const edm::ESGetToken<SiPixelLorentzAngle, SiPixelLorentzAngleRcd> siPixelLAEsToken_;
0166 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
0167
0168 const std::string dqmDir_;
0169 std::vector<std::string> newmodulelist_;
0170 const std::vector<double> fitRange_;
0171 const std::vector<double> fitParametersInitValues_;
0172 const std::vector<double> fitParametersInitValuesMuHFit_;
0173 FitParametersInitValuesMuHFitMap fitParametersInitValuesMuHFitMap_;
0174
0175 const int minHitsCut_;
0176 const std::string recordName_;
0177 SiPixelLorentzAngleCalibrationHistograms hists_;
0178 SiPixelLAHarvestMCS::fitStatus fitMCSHistogram(dqm::reco::MonitorElement* h_mean);
0179 std::pair<double, double> theFitRange_{-1.5, 1.5};
0180
0181 std::unique_ptr<TF1> f1_;
0182 FpixMuHResults fpixMuHResults;
0183
0184 const SiPixelLorentzAngle* currentLorentzAngle_;
0185 const TrackerTopology* tTopo;
0186 };
0187
0188
0189 SiPixelLorentzAnglePCLHarvesterMCS::SiPixelLorentzAnglePCLHarvesterMCS(const edm::ParameterSet& iConfig)
0190 : geomEsToken_(esConsumes<edm::Transition::BeginRun>()),
0191 topoEsTokenBR_(esConsumes<edm::Transition::BeginRun>()),
0192 topoEsTokenER_(esConsumes<edm::Transition::EndRun>()),
0193 siPixelLAEsToken_(esConsumes<edm::Transition::BeginRun>()),
0194 dqmDir_(iConfig.getParameter<std::string>("dqmDir")),
0195 newmodulelist_(iConfig.getParameter<std::vector<std::string>>("newmodulelist")),
0196 fitRange_(iConfig.getParameter<std::vector<double>>("fitRange")),
0197 fitParametersInitValues_(iConfig.getParameter<std::vector<double>>("fitParameters")),
0198 fitParametersInitValuesMuHFit_(iConfig.getParameter<std::vector<double>>("fitParametersMuHFit")),
0199 minHitsCut_(iConfig.getParameter<int>("minHitsCut")),
0200 recordName_(iConfig.getParameter<std::string>("record")) {
0201
0202
0203 if (fitRange_.size() == 2) {
0204 theFitRange_.first = fitRange_[0];
0205 theFitRange_.second = fitRange_[1];
0206 } else {
0207 throw cms::Exception("SiPixelLorentzAnglePCLHarvesterMCS") << "Wrong number of fit range parameters specified";
0208 }
0209
0210 if (fitParametersInitValues_.size() != 5) {
0211 throw cms::Exception("SiPixelLorentzAnglePCLHarvesterMCS")
0212 << "Wrong number of initial values for fit parameters specified";
0213 }
0214
0215 if (fitParametersInitValuesMuHFit_.size() != 4) {
0216 throw cms::Exception("SiPixelLorentzAnglePCLHarvesterMCS")
0217 << "Wrong number of initial values for fit parameters specified";
0218 }
0219
0220 fitParametersInitValuesMuHFitMap_["R1_P1"] = fitParametersInitValuesMuHFit_[0];
0221 fitParametersInitValuesMuHFitMap_["R1_P2"] = fitParametersInitValuesMuHFit_[1];
0222 fitParametersInitValuesMuHFitMap_["R2_P1"] = fitParametersInitValuesMuHFit_[2];
0223 fitParametersInitValuesMuHFitMap_["R2_P2"] = fitParametersInitValuesMuHFit_[3];
0224
0225
0226 edm::Service<cond::service::PoolDBOutputService> poolDbService;
0227 if (!poolDbService.isAvailable())
0228 throw cms::Exception("SiPixelLorentzAnglePCLHarvesterMCS") << "PoolDBService required";
0229 }
0230
0231
0232 void SiPixelLorentzAnglePCLHarvesterMCS::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0233
0234
0235 const TrackerGeometry* geom = &iSetup.getData(geomEsToken_);
0236 tTopo = &iSetup.getData(topoEsTokenBR_);
0237 currentLorentzAngle_ = &iSetup.getData(siPixelLAEsToken_);
0238 PixelTopologyMap map = PixelTopologyMap(geom, tTopo);
0239 hists_.nlay = geom->numberOfLayers(PixelSubdetector::PixelBarrel);
0240 hists_.nModules_.resize(hists_.nlay);
0241 hists_.nLadders_.resize(hists_.nlay);
0242 for (int i = 0; i < hists_.nlay; i++) {
0243 hists_.nModules_[i] = map.getPXBModules(i + 1);
0244 hists_.nLadders_[i] = map.getPXBLadders(i + 1);
0245 }
0246
0247
0248 if (!hists_.BPixnewDetIds_.empty() || !hists_.FPixnewDetIds_.empty())
0249 return;
0250
0251 if (!newmodulelist_.empty()) {
0252 for (auto const& modulename : newmodulelist_) {
0253 if (modulename.find("BPix_") != std::string::npos) {
0254 PixelBarrelName bn(modulename, true);
0255 const auto& detId = bn.getDetId(tTopo);
0256 hists_.BPixnewmodulename_.push_back(modulename);
0257 hists_.BPixnewDetIds_.push_back(detId.rawId());
0258 hists_.BPixnewModule_.push_back(bn.moduleName());
0259 hists_.BPixnewLayer_.push_back(bn.layerName());
0260 } else if (modulename.find("FPix_") != std::string::npos) {
0261 PixelEndcapName en(modulename, true);
0262 const auto& detId = en.getDetId(tTopo);
0263 hists_.FPixnewmodulename_.push_back(modulename);
0264 hists_.FPixnewDetIds_.push_back(detId.rawId());
0265 hists_.FPixnewDisk_.push_back(en.diskName());
0266 hists_.FPixnewBlade_.push_back(en.bladeName());
0267 }
0268 }
0269 }
0270
0271 uint count = 0;
0272 for (const auto& id : hists_.BPixnewDetIds_) {
0273 LogDebug("SiPixelLorentzAnglePCLHarvesterMCS") << id;
0274 count++;
0275 }
0276 LogDebug("SiPixelLorentzAnglePCLHarvesterMCS") << "Stored a total of " << count << " new detIds.";
0277
0278
0279 if (!hists_.detIdsList.empty())
0280 return;
0281
0282 std::vector<uint32_t> treatedIndices;
0283
0284 for (const auto& det : geom->detsPXB()) {
0285 const PixelGeomDetUnit* pixelDet = dynamic_cast<const PixelGeomDetUnit*>(det);
0286 const auto& layer = tTopo->pxbLayer(pixelDet->geographicalId());
0287 const auto& module = tTopo->pxbModule(pixelDet->geographicalId());
0288 int i_index = module + (layer - 1) * hists_.nModules_[layer - 1];
0289
0290 uint32_t rawId = pixelDet->geographicalId().rawId();
0291
0292
0293 if (std::find(hists_.BPixnewDetIds_.begin(), hists_.BPixnewDetIds_.end(), rawId) != hists_.BPixnewDetIds_.end())
0294 continue;
0295 if (std::find(treatedIndices.begin(), treatedIndices.end(), i_index) != treatedIndices.end()) {
0296 hists_.detIdsList[i_index].push_back(rawId);
0297 } else {
0298 hists_.detIdsList.insert(std::pair<uint32_t, std::vector<uint32_t>>(i_index, {rawId}));
0299 treatedIndices.push_back(i_index);
0300 }
0301 }
0302
0303 count = 0;
0304 for (const auto& i : treatedIndices) {
0305 for (const auto& id : hists_.detIdsList[i]) {
0306 LogDebug("SiPixelLorentzAnglePCLHarvesterMCS") << id;
0307 count++;
0308 };
0309 }
0310 LogDebug("SiPixelLorentzAnglePCLHarvesterMCS") << "Stored a total of " << count << " detIds.";
0311 }
0312
0313 void SiPixelLorentzAnglePCLHarvesterMCS::dqmEndJob(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter) {
0314 iGetter.cd();
0315 iGetter.setCurrentFolder(dqmDir_);
0316
0317
0318 std::string hname;
0319 int nBins = hists_.nAngles_ * hists_.nRings_ * hists_.nPanels_ * hists_.nSides_;
0320 hists_.h_fpixMeanHistoFitStatus_ =
0321 iBooker.book1D("fpixMeanHistoFitStatus",
0322 "fit status by angles/rings/panels/sides; angle/ring/panel/side; fit status",
0323 nBins,
0324 -0.5,
0325 nBins - 0.5);
0326 hists_.h_fpixMinClusterSizeCotAngle_ =
0327 iBooker.book1D("fpixMinClusterSizeCotAngle",
0328 "cot angle of minimal cluster size by angles/rings/panels/sides; angle/ring/panel/side; ",
0329 nBins,
0330 -0.5,
0331 nBins - 0.5);
0332 hists_.h_fpixNhitsClusterSizeCotAngle_ =
0333 iBooker.book1D("fpixNhitsClusterSizeCotAngle",
0334 "number of hits by angles/rings/panels/sides; angle/ring/panel/side; ",
0335 nBins,
0336 -0.5,
0337 nBins - 0.5);
0338
0339 std::string binNameAlpha;
0340 std::string binNameBeta;
0341
0342 const auto& prefix_ = fmt::sprintf("%s/FPix", dqmDir_);
0343 int histsCounter = 0;
0344 int nHistoExpected = 0;
0345 for (int r = 0; r < hists_.nRings_; ++r) {
0346 for (int p = 0; p < hists_.nPanels_; ++p) {
0347 for (int s = 0; s < hists_.nSides_; ++s) {
0348 int idx = getIndex(false, r, p, s);
0349 int idxBeta = getIndex(true, r, p, s);
0350 nHistoExpected++;
0351 hname = fmt::format("{}/R{}_P{}_z{}_alphaMean", dqmDir_, r + 1, p + 1, s + 1);
0352 if ((hists_.h_fpixMean_[idx] = iGetter.get(hname)) == nullptr) {
0353 edm::LogError("SiPixelLorentzAnglePCLHarvesterMCS::dqmEndJob")
0354 << "Failed to retrieve " << hname << " histogram";
0355 } else
0356 histsCounter++;
0357
0358 nHistoExpected++;
0359 hname = fmt::format("{}/R{}_P{}_z{}_betaMean", dqmDir_, r + 1, p + 1, s + 1);
0360 if ((hists_.h_fpixMean_[idxBeta] = iGetter.get(hname)) == nullptr) {
0361 edm::LogError("SiPixelLorentzAnglePCLHarvesterMCS::dqmEndJob")
0362 << "Failed to retrieve " << hname << " histogram";
0363 } else
0364 histsCounter++;
0365
0366 nHistoExpected++;
0367 hname = fmt::format("{}/R{}_P{}_z{}_alpha", prefix_, r + 1, p + 1, s + 1);
0368 if ((hists_.h_fpixAngleSize_[idx] = iGetter.get(hname)) == nullptr) {
0369 edm::LogError("SiPixelLorentzAnglePCLHarvesterMCS::dqmEndJob")
0370 << "Failed to retrieve " << hname << " histogram";
0371 } else
0372 histsCounter++;
0373
0374 nHistoExpected++;
0375 hname = fmt::format("{}/R{}_P{}_z{}_beta", prefix_, r + 1, p + 1, s + 1);
0376 if ((hists_.h_fpixAngleSize_[idxBeta] = iGetter.get(hname)) == nullptr) {
0377 edm::LogError("SiPixelLorentzAnglePCLHarvesterMCS::dqmEndJob")
0378 << "Failed to retrieve " << hname << " histogram";
0379 } else
0380 histsCounter++;
0381
0382 for (int m = 0; m < 3; ++m) {
0383 nHistoExpected++;
0384 hname = fmt::format("{}/R{}_P{}_z{}_B{}", prefix_, r + 1, p + 1, s + 1, m);
0385 if ((hists_.h_fpixMagField_[m][idx] = iGetter.get(hname)) == nullptr) {
0386 edm::LogError("SiPixelLorentzAnglePCLHarvesterMCS::dqmEndJob")
0387 << "Failed to retrieve " << hname << " histogram";
0388 } else
0389 histsCounter++;
0390 }
0391
0392
0393 int binAlpha = idx + 1;
0394 int binBeta = idxBeta + 1;
0395 char sign = s == 0 ? '-' : '+';
0396 binNameAlpha = fmt::sprintf("#alpha: R%d_P%d_z%c", r + 1, p + 1, sign);
0397 binNameBeta = fmt::sprintf("#beta:R%d_P%d_z%c", r + 1, p + 1, sign);
0398 hists_.h_fpixMeanHistoFitStatus_->setBinLabel(binAlpha, binNameAlpha);
0399 hists_.h_fpixMeanHistoFitStatus_->setBinLabel(binBeta, binNameBeta);
0400 hists_.h_fpixMeanHistoFitStatus_->setBinContent(binAlpha, SiPixelLAHarvestMCS::kFitNotPerformed);
0401 hists_.h_fpixMeanHistoFitStatus_->setBinContent(binBeta, SiPixelLAHarvestMCS::kFitNotPerformed);
0402
0403 hists_.h_fpixMinClusterSizeCotAngle_->setBinLabel(binAlpha, binNameAlpha);
0404 hists_.h_fpixMinClusterSizeCotAngle_->setBinLabel(binBeta, binNameBeta);
0405 hists_.h_fpixNhitsClusterSizeCotAngle_->setBinLabel(binAlpha, binNameAlpha);
0406 hists_.h_fpixNhitsClusterSizeCotAngle_->setBinLabel(binBeta, binNameBeta);
0407 }
0408 }
0409 }
0410
0411 if (histsCounter != nHistoExpected) {
0412 edm::LogError("SiPixelLorentzAnglePCLHarvesterMCS::dqmEndJob")
0413 << "Failed to retrieve all histograms, expected 56 got " << histsCounter;
0414 return;
0415 }
0416
0417
0418 iBooker.setCurrentFolder(fmt::format("{}Harvesting", dqmDir_));
0419 int nBinsMuH = hists_.nRings_ * hists_.nPanels_;
0420 hists_.h_fpixFitStatusMuH_ = iBooker.book1D(
0421 "fpixFitStatusMuH", "muH fit status by rings/panels; ring/panel; fitStatus", nBinsMuH, -0.5, nBinsMuH - 0.5);
0422 hists_.h_fpixMuH_ =
0423 iBooker.book1D("fpixMuH", "muH by rings/panels; ring/panel; #muH [1/T]", nBinsMuH, -0.5, nBinsMuH - 0.5);
0424 hists_.h_fpixDeltaMuH_ = iBooker.book1D(
0425 "fpixDeltaMuH", "#Delta muH by rings/panels; ring/panel; #Delta #muH [1/T]", nBinsMuH, -0.5, nBinsMuH - 0.5);
0426 hists_.h_fpixRelDeltaMuH_ = iBooker.book1D("fpixRelDeltaMuH",
0427 "#Delta #muH/#muH by rings/panels; ring/panel; #Delta #muH/#MuH",
0428 nBinsMuH,
0429 -0.5,
0430 nBinsMuH - 0.5);
0431 std::string binName;
0432 for (int r = 0; r < hists_.nRings_; ++r) {
0433 for (int p = 0; p < hists_.nPanels_; ++p) {
0434 int idx = r * hists_.nPanels_ + p + 1;
0435 binName = fmt::sprintf("R%d_P%d", r + 1, p + 1);
0436 hists_.h_fpixFitStatusMuH_->setBinLabel(idx, binName);
0437 hists_.h_fpixFitStatusMuH_->setBinContent(idx, SiPixelLAHarvestMCS::kFitNotPerformed);
0438 hists_.h_fpixMuH_->setBinLabel(idx, binName);
0439 hists_.h_fpixDeltaMuH_->setBinLabel(idx, binName);
0440 hists_.h_fpixRelDeltaMuH_->setBinLabel(idx, binName);
0441 }
0442 }
0443
0444
0445 int nSizeBins = hists_.h_fpixAngleSize_[0]->getNbinsY();
0446 double minSize = hists_.h_fpixAngleSize_[0]->getAxisMin(2);
0447 double maxSize = hists_.h_fpixAngleSize_[0]->getAxisMax(2);
0448 TH1D* h_slice = new TH1D("h_slice", "slice of cot_angle histogram", nSizeBins, minSize, maxSize);
0449 f1_ = std::make_unique<TF1>("f1", SiPixelLAHarvestMCS::MCSFitFunction, -3., 3., 5);
0450 f1_->SetParNames("Offset", "RMS Constant", "SlopeL", "cot(angle)_min", "SlopeR");
0451
0452 for (int r = 0; r < hists_.nRings_; ++r) {
0453 for (int p = 0; p < hists_.nPanels_; ++p) {
0454 SiPixelLAHarvestMCS::FitFPixMuH fitMuH;
0455 SiPixelLAHarvestMCS::FPixMuH fmh;
0456 for (int s = 0; s < hists_.nSides_; ++s) {
0457 int idx = getIndex(false, r, p, s);
0458 int idxBeta = getIndex(true, r, p, s);
0459 int binAlpha = idx + 1;
0460 int binBeta = idxBeta + 1;
0461
0462 int entriesAlpha = hists_.h_fpixAngleSize_[idx]->getEntries();
0463 int entriesBeta = hists_.h_fpixAngleSize_[idxBeta]->getEntries();
0464 hists_.h_fpixNhitsClusterSizeCotAngle_->setBinContent(binAlpha, entriesAlpha);
0465 hists_.h_fpixNhitsClusterSizeCotAngle_->setBinContent(binBeta, entriesBeta);
0466 findMean(hists_.h_fpixAngleSize_[idx], hists_.h_fpixMean_[idx], h_slice);
0467 findMean(hists_.h_fpixAngleSize_[idxBeta], hists_.h_fpixMean_[idxBeta], h_slice);
0468
0469 SiPixelLAHarvestMCS::fitStatus statusAlphaFit = entriesAlpha < minHitsCut_
0470 ? SiPixelLAHarvestMCS::kFitNotPerformed
0471 : fitMCSHistogram(hists_.h_fpixMean_[idx]);
0472 SiPixelLAHarvestMCS::fitStatus statusBetaFit = entriesBeta < minHitsCut_
0473 ? SiPixelLAHarvestMCS::kFitNotPerformed
0474 : fitMCSHistogram(hists_.h_fpixMean_[idxBeta]);
0475
0476 hists_.h_fpixMeanHistoFitStatus_->setBinContent(binAlpha, statusAlphaFit);
0477 hists_.h_fpixMeanHistoFitStatus_->setBinContent(binBeta, statusBetaFit);
0478
0479 if (entriesAlpha < minHitsCut_ || entriesBeta < minHitsCut_)
0480 continue;
0481
0482 assert(strcmp(f1_->GetName(), "f1") == 0);
0483
0484 if (statusAlphaFit == SiPixelLAHarvestMCS::kFitConverged &&
0485 statusBetaFit == SiPixelLAHarvestMCS::kFitConverged) {
0486 double shiftX = hists_.h_fpixMean_[idx]->getTH1()->GetFunction("f1")->GetParameter(3);
0487 double errShiftX = hists_.h_fpixMean_[idx]->getTH1()->GetFunction("f1")->GetParError(3);
0488 double shiftY = hists_.h_fpixMean_[idxBeta]->getTH1()->GetFunction("f1")->GetParameter(3);
0489 double errShiftY = hists_.h_fpixMean_[idxBeta]->getTH1()->GetFunction("f1")->GetParError(3);
0490
0491 hists_.h_fpixMinClusterSizeCotAngle_->setBinContent(binAlpha, shiftX);
0492 hists_.h_fpixMinClusterSizeCotAngle_->setBinError(binAlpha, errShiftX);
0493 hists_.h_fpixMinClusterSizeCotAngle_->setBinContent(binBeta, shiftY);
0494 hists_.h_fpixMinClusterSizeCotAngle_->setBinError(binBeta, errShiftY);
0495
0496 fmh.shiftx = shiftX;
0497 fmh.shiftx_err = errShiftX;
0498 fmh.shifty = shiftY;
0499 fmh.shifty_err = errShiftY;
0500 fmh.bfield[0] = hists_.h_fpixMagField_[0][idx]->getMean();
0501 fmh.bfield[1] = hists_.h_fpixMagField_[1][idx]->getMean();
0502 fmh.bfield[2] = hists_.h_fpixMagField_[2][idx]->getMean();
0503 fitMuH.add(fmh);
0504 }
0505 }
0506
0507 if (fitMuH.size() == hists_.nSides_) {
0508 std::string fpixPartNames = "R" + std::to_string(r + 1) + "_P" + std::to_string(p + 1);
0509 double initMuH = fitParametersInitValuesMuHFitMap_[fpixPartNames];
0510 int status = fitMuH.fit(initMuH);
0511 int idxMuH = r * hists_.nPanels_ + p + 1;
0512 double muH = fitMuH.getMuH();
0513 double muHErr = fitMuH.getMuHErr();
0514 hists_.h_fpixFitStatusMuH_->setBinContent(idxMuH, status);
0515 hists_.h_fpixMuH_->setBinContent(idxMuH, muH);
0516 hists_.h_fpixMuH_->setBinError(idxMuH, muHErr);
0517 fpixMuHResults.insert(std::make_pair(fpixPartNames, std::make_pair(muH, muH)));
0518 }
0519 }
0520 }
0521
0522 std::shared_ptr<SiPixelLorentzAngle> LorentzAngle = std::make_shared<SiPixelLorentzAngle>();
0523
0524 bool isPayloadChanged{false};
0525
0526 for (const auto& [id, value] : currentLorentzAngle_->getLorentzAngles()) {
0527 DetId ID = DetId(id);
0528 float muHForDB = value;
0529 if (ID.subdetId() == PixelSubdetector::PixelEndcap) {
0530 PixelEndcapName pen(ID, tTopo, true);
0531 int panel = pen.pannelName();
0532 int ring = pen.ringName();
0533 std::string fpixPartNames = "R" + std::to_string(ring) + "_P" + std::to_string(panel);
0534 if (fpixMuHResults.find(fpixPartNames) != fpixMuHResults.end()) {
0535 double measuredMuH = fpixMuHResults[fpixPartNames].first;
0536 double deltaMuH = value - measuredMuH;
0537 double deltaMuHoverMuH = deltaMuH / value;
0538 int idxMuH = (ring - 1) * hists_.nPanels_ + panel;
0539 hists_.h_fpixDeltaMuH_->setBinContent(idxMuH, deltaMuH);
0540 hists_.h_fpixRelDeltaMuH_->setBinContent(idxMuH, deltaMuHoverMuH);
0541 muHForDB = measuredMuH;
0542 if (!isPayloadChanged && (deltaMuHoverMuH != 0.f))
0543 isPayloadChanged = true;
0544 }
0545 }
0546 if (!LorentzAngle->putLorentzAngle(id, muHForDB)) {
0547 edm::LogError("SiPixelLorentzAnglePCLHarvesterMCS")
0548 << "[SiPixelLorentzAnglePCLHarvesterMCS::dqmEndJob]: detid (" << id << ") already exists";
0549 }
0550 }
0551
0552 if (isPayloadChanged) {
0553
0554 edm::Service<cond::service::PoolDBOutputService> mydbservice;
0555 if (mydbservice.isAvailable()) {
0556 try {
0557 mydbservice->writeOneIOV(*LorentzAngle, mydbservice->currentTime(), recordName_);
0558 } catch (const cond::Exception& er) {
0559 edm::LogError("SiPixelLorentzAngleDB") << er.what();
0560 } catch (const std::exception& er) {
0561 edm::LogError("SiPixelLorentzAngleDB") << "caught std::exception " << er.what();
0562 }
0563 } else {
0564 edm::LogError("SiPixelLorentzAngleDB") << "Service is unavailable";
0565 }
0566 } else {
0567 edm::LogPrint("SiPixelLorentzAngleDB") << __PRETTY_FUNCTION__ << " there is no new valid measurement to append! ";
0568 }
0569 }
0570
0571
0572 void SiPixelLorentzAnglePCLHarvesterMCS::findMean(dqm::reco::MonitorElement* h_2D,
0573 dqm::reco::MonitorElement* h_mean,
0574 TH1D* h_slice) {
0575 int n_x = h_2D->getNbinsX();
0576 int n_y = h_2D->getNbinsY();
0577
0578 for (int i = 1; i <= n_x; i++) {
0579 h_slice->Reset("ICE");
0580
0581
0582
0583 for (int j = 1; j <= n_y; j++) {
0584 h_slice->SetBinContent(j, h_2D->getBinContent(i, j));
0585 }
0586 double mean = h_slice->GetMean();
0587 double error = h_slice->GetMeanError();
0588 h_mean->setBinContent(i, mean);
0589 h_mean->setBinError(i, error);
0590 }
0591 }
0592
0593
0594 SiPixelLAHarvestMCS::fitStatus SiPixelLorentzAnglePCLHarvesterMCS::fitMCSHistogram(dqm::reco::MonitorElement* h_mean) {
0595 SiPixelLAHarvestMCS::fitStatus retVal;
0596
0597 f1_->SetParameters(fitParametersInitValues_[0],
0598 fitParametersInitValues_[1],
0599 fitParametersInitValues_[2],
0600 fitParametersInitValues_[3],
0601 fitParametersInitValues_[4]);
0602
0603 int nFits = 0;
0604 while (nFits < 5) {
0605 nFits++;
0606 double fitMin = f1_->GetParameter(3) + theFitRange_.first;
0607 double fitMax = f1_->GetParameter(3) + theFitRange_.second;
0608
0609 if (fitMin < -3)
0610 fitMin = -3;
0611 if (fitMax > 3)
0612 fitMax = 3;
0613
0614 TFitResultPtr r = h_mean->getTH1()->Fit(f1_.get(), "ERSM", "", fitMin, fitMax);
0615 retVal = r == -1 ? SiPixelLAHarvestMCS::kNoFitResult
0616 : r->IsValid() ? SiPixelLAHarvestMCS::kFitConverged
0617 : SiPixelLAHarvestMCS::kFitFailed;
0618 }
0619 return retVal;
0620 }
0621
0622
0623 void SiPixelLorentzAnglePCLHarvesterMCS::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0624 edm::ParameterSetDescription desc;
0625 desc.setComment("Harvester module of the SiPixel Lorentz Angle PCL monitoring workflow for MinimalClusterSizeMethod");
0626 desc.add<std::string>("dqmDir", "AlCaReco/SiPixelLorentzAngle")->setComment("the directory of PCL Worker output");
0627 desc.add<std::vector<std::string>>("newmodulelist", {})->setComment("the list of DetIds for new sensors");
0628 desc.add<std::vector<double>>("fitRange", {-1.5, 1.5})->setComment("range to perform the fit");
0629 desc.add<std::vector<double>>("fitParameters", {1., 0.1, -1.6, 0., 1.6})
0630 ->setComment("initial values for fit parameters");
0631 desc.add<std::vector<double>>("fitParametersMuHFit", {0.08, 0.08, 0.08, 0.08})
0632 ->setComment("initial values for fit parameters (muH fit)");
0633 desc.add<int>("minHitsCut", 1000)->setComment("cut on minimum number of on-track hits to accept measurement");
0634 desc.add<std::string>("record", "SiPixelLorentzAngleRcdMCS")->setComment("target DB record");
0635 descriptions.addWithDefaultLabel(desc);
0636 }
0637
0638 int SiPixelLorentzAnglePCLHarvesterMCS::getIndex(bool isBetaAngle, int r, int p, int s) {
0639 int idx = hists_.nSides_ * hists_.nPanels_ * r + hists_.nSides_ * p + s;
0640 return (isBetaAngle ? idx + hists_.betaStartIdx_ : idx);
0641 }
0642 DEFINE_FWK_MODULE(SiPixelLorentzAnglePCLHarvesterMCS);
