Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:35:12

0001 // -*- C++ -*-
0002 //
0003 // Package:    CalibTracker/SiPixelLorentzAnglePCLHarvesterMCS
0004 // Class:      SiPixelLorentzAnglePCLHarvesterMCS
0005 //
0006 /**\class SiPixelLorentzAnglePCLHarvesterMCS SiPixelLorentzAnglePCLHarvesterMCS.cc CalibTracker/SiPixelLorentzAngle/src/SiPixelLorentzAnglePCLHarvesterMCS.cc
0007  Description: reads the intermediate ALCAPROMPT DQMIO-like dataset and performs the fitting of the SiPixel Lorentz Angle in the Prompt Calibration Loop
0008  Implementation:
0009  Reads the 16 2D histograms of the cluster size x/y vs cot(alpha/beta) and 16*3 histograms of the magnetic field components created by SiPixelLorentzAnglePCLWorker module. The cluster size x/y vs cot(alpha/beta) histograms are used to generate 1D profiles (average cluster size x/y vs cot(alpha/beta)) which are then fit and the values of the cot (alpha/beta) for which the cluster sizes are minimal are determined. The obtained cot(alpha/beta)_min value for z- and z+ side are used to perform fit and the muH for different rings and panels of the Pixel Forward Phase 1 detector using the formulas: 
0010  cot(alpha)_min = vx/vz  = (muHBy + muH^2*Bz*Bx)/(1+muH^2*Bz^2)
0011  cot(beta)_min = vy/vz  = -(muHBx - muH^2*Bz*Bx)/(1+muH^2*Bz^2)
0012   
0013 The extracted value of the muH are stored in an output sqlite file which is then uploaded to the conditions database.
0014 */
0015 //
0016 // Original Author:  tsusa
0017 //         Created:  Sat, 14 Jan 2021 10:12:21 GMT
0018 //
0019 //
0020 
0021 // system includes
0022 
0023 #include <fmt/format.h>
0024 #include <fmt/printf.h>
0025 #include <fstream>
0026 
0027 // user includes
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 // for ROOT fits
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;  // Shift in x direction
0070     double shifty;  // Shift in y direction
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 }  // namespace SiPixelLAHarvestMCS
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   // es tokens
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   // initialize the fit range
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   // first ensure DB output service is available
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   // geometry
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   // list of modules already filled, then return (we already entered here)
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   // list of modules already filled, return (we already entered here)
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     // if the detId is already accounted for in the special class, do not attach it
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   // get mean and size-angle hists, book summary hists
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         // set labels & init summary histos
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   // book hervesting hists
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   // make and fit profile hists, fit muH
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         }  // if fut converged
0505       }    // loop over z sides
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);  // use det-id phaseq
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     }  // if endcap
0546     if (!LorentzAngle->putLorentzAngle(id, muHForDB)) {
0547       edm::LogError("SiPixelLorentzAnglePCLHarvesterMCS")
0548           << "[SiPixelLorentzAnglePCLHarvesterMCS::dqmEndJob]: detid (" << id << ") already exists";
0549     }
0550   }
0551 
0552   if (isPayloadChanged) {
0553     // fill the DB object record
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     //loop over bins in size
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   }  // end loop over bins in depth
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);