Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-24 22:51:39

0001 // -*- C++ -*-
0002 //
0003 // Package:    CalibTracker/SiPixelLorentzAnglePCLHarvester
0004 // Class:      SiPixelLorentzAnglePCLHarvester
0005 //
0006 /**\class SiPixelLorentzAnglePCLHarvester SiPixelLorentzAnglePCLHarvester.cc CalibTracker/SiPixelLorentzAngle/src/SiPixelLorentzAnglePCLHarvester.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 2D histograms of the drift vs depth created by SiPixelLorentzAnglePCLWorker modules and generates 1D profiles which are then fit
0010      with a 5th order polinomial. The extracted value of the tan(theta_L)/B are stored in an output sqlite file which is then uploaded to the conditions database
0011 */
0012 //
0013 // Original Author:  mmusich
0014 //         Created:  Sat, 29 May 2021 14:46:19 GMT
0015 //
0016 //
0017 
0018 // system includes
0019 #include <fmt/format.h>
0020 #include <fmt/printf.h>
0021 #include <fstream>
0022 
0023 // user includes
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  * Auxilliary struct to store fit results
0044  */
0045 namespace SiPixelLAHarvest {
0046   struct fitResults {
0047   public:
0048     fitResults() {
0049       // set all parameters to default
0050       p0 = p1 = p2 = p3 = p4 = p5 = 0.;
0051       e0 = e1 = e2 = e3 = e4 = e5 = 0.;
0052       chi2 = prob = redChi2 = tan_LA = error_LA = -9999.;
0053       ndf = -999;
0054     };
0055 
0056     double p0;
0057     double e0;
0058     double p1;
0059     double e1;
0060     double p2;
0061     double e2;
0062     double p3;
0063     double e3;
0064     double p4;
0065     double e4;
0066     double p5;
0067     double e5;
0068     double chi2;
0069     int ndf;
0070     double prob;
0071     double redChi2;
0072     double tan_LA;
0073     double error_LA;
0074   };
0075 }  // namespace SiPixelLAHarvest
0076 
0077 //------------------------------------------------------------------------------
0078 class SiPixelLorentzAnglePCLHarvester : public DQMEDHarvester {
0079 public:
0080   SiPixelLorentzAnglePCLHarvester(const edm::ParameterSet&);
0081   ~SiPixelLorentzAnglePCLHarvester() override = default;
0082   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0083 
0084   static void fillDescriptions(edm::ConfigurationDescriptions&);
0085 
0086 private:
0087   void dqmEndJob(DQMStore::IBooker&, DQMStore::IGetter&) override;
0088   void endRun(const edm::Run&, const edm::EventSetup&) override;
0089   void findMean(MonitorElement* h_drift_depth_adc_slice_, int i, int i_ring);
0090   SiPixelLAHarvest::fitResults fitAndStore(std::shared_ptr<SiPixelLorentzAngle> theLA, int i_idx, int i_lay, int i_mod);
0091 
0092   // es tokens
0093   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomEsToken_;
0094   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoEsTokenBR_, topoEsTokenER_;
0095   edm::ESGetToken<SiPixelLorentzAngle, SiPixelLorentzAngleRcd> siPixelLAEsToken_;
0096   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
0097 
0098   std::vector<std::string> newmodulelist_;
0099   const std::string dqmDir_;
0100   const double fitChi2Cut_;
0101   const int minHitsCut_;
0102   const std::string recordName_;
0103   std::unique_ptr<TF1> f1;
0104   float width_;
0105 
0106   SiPixelLorentzAngleCalibrationHistograms hists;
0107   const SiPixelLorentzAngle* currentLorentzAngle;
0108   const MagneticField* magField;
0109   std::unique_ptr<TrackerTopology> theTrackerTopology;
0110 };
0111 
0112 //------------------------------------------------------------------------------
0113 SiPixelLorentzAnglePCLHarvester::SiPixelLorentzAnglePCLHarvester(const edm::ParameterSet& iConfig)
0114     : geomEsToken_(esConsumes<edm::Transition::BeginRun>()),
0115       topoEsTokenBR_(esConsumes<edm::Transition::BeginRun>()),
0116       topoEsTokenER_(esConsumes<edm::Transition::EndRun>()),
0117       siPixelLAEsToken_(esConsumes<edm::Transition::BeginRun>()),
0118       magneticFieldToken_(esConsumes<edm::Transition::BeginRun>()),
0119       newmodulelist_(iConfig.getParameter<std::vector<std::string>>("newmodulelist")),
0120       dqmDir_(iConfig.getParameter<std::string>("dqmDir")),
0121       fitChi2Cut_(iConfig.getParameter<double>("fitChi2Cut")),
0122       minHitsCut_(iConfig.getParameter<int>("minHitsCut")),
0123       recordName_(iConfig.getParameter<std::string>("record")) {
0124   // first ensure DB output service is available
0125   edm::Service<cond::service::PoolDBOutputService> poolDbService;
0126   if (!poolDbService.isAvailable())
0127     throw cms::Exception("SiPixelLorentzAnglePCLHarvester") << "PoolDBService required";
0128 }
0129 
0130 //------------------------------------------------------------------------------
0131 void SiPixelLorentzAnglePCLHarvester::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0132   // geometry
0133   const TrackerGeometry* geom = &iSetup.getData(geomEsToken_);
0134   const TrackerTopology* tTopo = &iSetup.getData(topoEsTokenBR_);
0135 
0136   magField = &iSetup.getData(magneticFieldToken_);
0137   currentLorentzAngle = &iSetup.getData(siPixelLAEsToken_);
0138 
0139   PixelTopologyMap map = PixelTopologyMap(geom, tTopo);
0140   hists.nlay = geom->numberOfLayers(PixelSubdetector::PixelBarrel);
0141   hists.nModules_.resize(hists.nlay);
0142   hists.nLadders_.resize(hists.nlay);
0143   for (int i = 0; i < hists.nlay; i++) {
0144     hists.nModules_[i] = map.getPXBModules(i + 1);
0145     hists.nLadders_[i] = map.getPXBLadders(i + 1);
0146   }
0147 
0148   // list of modules already filled, then return (we already entered here)
0149   if (!hists.BPixnewDetIds_.empty() || !hists.FPixnewDetIds_.empty())
0150     return;
0151 
0152   if (!newmodulelist_.empty()) {
0153     for (auto const& modulename : newmodulelist_) {
0154       if (modulename.find("BPix_") != std::string::npos) {
0155         PixelBarrelName bn(modulename, true);
0156         const auto& detId = bn.getDetId(tTopo);
0157         hists.BPixnewmodulename_.push_back(modulename);
0158         hists.BPixnewDetIds_.push_back(detId.rawId());
0159         hists.BPixnewModule_.push_back(bn.moduleName());
0160         hists.BPixnewLayer_.push_back(bn.layerName());
0161       } else if (modulename.find("FPix_") != std::string::npos) {
0162         PixelEndcapName en(modulename, true);
0163         const auto& detId = en.getDetId(tTopo);
0164         hists.FPixnewmodulename_.push_back(modulename);
0165         hists.FPixnewDetIds_.push_back(detId.rawId());
0166         hists.FPixnewDisk_.push_back(en.diskName());
0167         hists.FPixnewBlade_.push_back(en.bladeName());
0168       }
0169     }
0170   }
0171 
0172   uint count = 0;
0173   for (const auto& id : hists.BPixnewDetIds_) {
0174     LogDebug("SiPixelLorentzAnglePCLHarvester") << id;
0175     count++;
0176   }
0177   LogDebug("SiPixelLorentzAnglePCLHarvester") << "Stored a total of " << count << " new detIds.";
0178 
0179   // list of modules already filled, return (we already entered here)
0180   if (!hists.detIdsList.empty())
0181     return;
0182 
0183   std::vector<uint32_t> treatedIndices;
0184 
0185   for (auto det : geom->detsPXB()) {
0186     const PixelGeomDetUnit* pixelDet = dynamic_cast<const PixelGeomDetUnit*>(det);
0187     width_ = pixelDet->surface().bounds().thickness();
0188     const auto& layer = tTopo->pxbLayer(pixelDet->geographicalId());
0189     const auto& module = tTopo->pxbModule(pixelDet->geographicalId());
0190     int i_index = module + (layer - 1) * hists.nModules_[layer - 1];
0191 
0192     uint32_t rawId = pixelDet->geographicalId().rawId();
0193 
0194     // if the detId is already accounted for in the special class, do not attach it
0195     if (std::find(hists.BPixnewDetIds_.begin(), hists.BPixnewDetIds_.end(), rawId) != hists.BPixnewDetIds_.end())
0196       continue;
0197 
0198     if (std::find(treatedIndices.begin(), treatedIndices.end(), i_index) != treatedIndices.end()) {
0199       hists.detIdsList[i_index].push_back(rawId);
0200     } else {
0201       hists.detIdsList.insert(std::pair<uint32_t, std::vector<uint32_t>>(i_index, {rawId}));
0202       treatedIndices.push_back(i_index);
0203     }
0204   }
0205 
0206   count = 0;
0207   for (const auto& i : treatedIndices) {
0208     for (const auto& id : hists.detIdsList[i]) {
0209       LogDebug("SiPixelLorentzAnglePCLHarvester") << id;
0210       count++;
0211     };
0212   }
0213   LogDebug("SiPixelLorentzAnglePCLHarvester") << "Stored a total of " << count << " detIds.";
0214 }
0215 
0216 //------------------------------------------------------------------------------
0217 void SiPixelLorentzAnglePCLHarvester::endRun(edm::Run const& run, edm::EventSetup const& isetup) {
0218   if (!theTrackerTopology) {
0219     theTrackerTopology = std::make_unique<TrackerTopology>(isetup.getData(topoEsTokenER_));
0220   }
0221 }
0222 
0223 //------------------------------------------------------------------------------
0224 void SiPixelLorentzAnglePCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter) {
0225   // go in the right directory
0226   iGetter.cd();
0227   iGetter.setCurrentFolder(dqmDir_);
0228 
0229   // fetch the 2D histograms
0230   for (int i_layer = 1; i_layer <= hists.nlay; i_layer++) {
0231     const auto& prefix_ = fmt::sprintf("%s/BPix/BPixLayer%i", dqmDir_, i_layer);
0232     for (int i_module = 1; i_module <= hists.nModules_[i_layer - 1]; i_module++) {
0233       int i_index = i_module + (i_layer - 1) * hists.nModules_[i_layer - 1];
0234 
0235       hists.h_drift_depth_[i_index] =
0236           iGetter.get(fmt::format("{}/h_drift_depth_layer{}_module{}", prefix_, i_layer, i_module));
0237 
0238       if (hists.h_drift_depth_[i_index] == nullptr) {
0239         edm::LogError("SiPixelLorentzAnglePCLHarvester::dqmEndJob")
0240             << "Failed to retrieve electron drift over depth for layer " << i_layer << ", module " << i_module << ".";
0241         continue;
0242       }
0243 
0244       hists.h_drift_depth_adc_[i_index] =
0245           iGetter.get(fmt::format("{}/h_drift_depth_adc_layer{}_module{}", prefix_, i_layer, i_module));
0246 
0247       hists.h_drift_depth_adc2_[i_index] =
0248           iGetter.get(fmt::format("{}/h_drift_depth_adc2_layer{}_module{}", prefix_, i_layer, i_module));
0249 
0250       hists.h_drift_depth_noadc_[i_index] =
0251           iGetter.get(fmt::format("{}/h_drift_depth_noadc_layer{}_module{}", prefix_, i_layer, i_module));
0252 
0253       hists.h_mean_[i_index] = iGetter.get(fmt::format("{}/h_mean_layer{}_module{}", dqmDir_, i_layer, i_module));
0254 
0255       hists.h_drift_depth_[i_index]->divide(
0256           hists.h_drift_depth_adc_[i_index], hists.h_drift_depth_noadc_[i_index], 1., 1., "");
0257     }
0258   }
0259 
0260   // fetch the new modules 2D histograms
0261   for (int i = 0; i < (int)hists.BPixnewDetIds_.size(); i++) {
0262     int new_index = i + 1 + hists.nModules_[hists.nlay - 1] + (hists.nlay - 1) * hists.nModules_[hists.nlay - 1];
0263 
0264     hists.h_drift_depth_[new_index] = iGetter.get(
0265         fmt::format("{}/h_BPixnew_drift_depth_{}", dqmDir_ + "/BPix/NewModules", hists.BPixnewmodulename_[i]));
0266 
0267     if (hists.h_drift_depth_[new_index] == nullptr) {
0268       edm::LogError("SiPixelLorentzAnglePCLHarvester")
0269           << "Failed to retrieve electron drift over depth for new module " << hists.BPixnewmodulename_[i] << ".";
0270       continue;
0271     }
0272 
0273     hists.h_drift_depth_adc_[new_index] = iGetter.get(
0274         fmt::format("{}/h_BPixnew_drift_depth_adc_{}", dqmDir_ + "/BPix/NewModules", hists.BPixnewmodulename_[i]));
0275 
0276     hists.h_drift_depth_adc2_[new_index] = iGetter.get(
0277         fmt::format("{}/h_BPixnew_drift_depth_adc2_{}", dqmDir_ + "/BPix/NewModules", hists.BPixnewmodulename_[i]));
0278 
0279     hists.h_drift_depth_noadc_[new_index] = iGetter.get(
0280         fmt::format("{}/h_BPixnew_drift_depth_noadc_{}", dqmDir_ + "/BPix/NewModules", hists.BPixnewmodulename_[i]));
0281 
0282     hists.h_mean_[new_index] = iGetter.get(fmt::format("{}/h_BPixnew_mean_{}", dqmDir_, hists.BPixnewmodulename_[i]));
0283 
0284     hists.h_drift_depth_[new_index]->divide(
0285         hists.h_drift_depth_adc_[new_index], hists.h_drift_depth_noadc_[new_index], 1., 1., "");
0286   }
0287 
0288   hists.h_bySectOccupancy_ = iGetter.get(fmt::format("{}/h_bySectorOccupancy", dqmDir_ + "/SectorMonitoring"));
0289   if (hists.h_bySectOccupancy_ == nullptr) {
0290     edm::LogError("SiPixelLorentzAnglePCLHarvester") << "Failed to retrieve the hit on track occupancy.";
0291     return;
0292   }
0293 
0294   int hist_drift_;
0295   int hist_depth_;
0296   double min_drift_;
0297   double max_drift_;
0298 
0299   if (hists.h_drift_depth_adc_[1] != nullptr) {
0300     hist_drift_ = hists.h_drift_depth_adc_[1]->getNbinsX();
0301     hist_depth_ = hists.h_drift_depth_adc_[1]->getNbinsY();
0302     min_drift_ = hists.h_drift_depth_adc_[1]->getAxisMin(1);
0303     max_drift_ = hists.h_drift_depth_adc_[1]->getAxisMax(1);
0304   } else {
0305     hist_drift_ = 100;
0306     hist_depth_ = 50;
0307     min_drift_ = -500.;
0308     max_drift_ = 500.;
0309   }
0310 
0311   iBooker.setCurrentFolder(fmt::format("{}Harvesting", dqmDir_));
0312   MonitorElement* h_drift_depth_adc_slice_ =
0313       iBooker.book1D("h_drift_depth_adc_slice", "slice of adc histogram", hist_drift_, min_drift_, max_drift_);
0314 
0315   // book histogram of differences
0316   MonitorElement* h_diffLA = iBooker.book1D(
0317       "h_diffLA", "difference in #mu_{H}; #Delta #mu_{H}/#mu_{H} (old-new)/old [%];n. modules", 100, -150, 150);
0318 
0319   // retrieve the number of bins from the other monitoring histogram
0320   const auto& maxSect = hists.h_bySectOccupancy_->getNbinsX();
0321   const double lo = -0.5;
0322   const double hi = maxSect - 0.5;
0323 
0324   // this will be booked in the Harvesting folder
0325   iBooker.setCurrentFolder(fmt::format("{}Harvesting/SectorMonitoring", dqmDir_));
0326   std::string repText = "%s tan#theta_{LA}/B by sector;pixel sector;%s tan(#theta_{LA})/B [1/T]";
0327   hists.h_bySectMeasLA_ =
0328       iBooker.book1D("h_LAbySector_Measured", fmt::sprintf(repText, "measured", "measured"), maxSect, lo, hi);
0329   hists.h_bySectSetLA_ =
0330       iBooker.book1D("h_LAbySector_Accepted", fmt::sprintf(repText, "accepted", "accepted"), maxSect, lo, hi);
0331   hists.h_bySectRejectLA_ =
0332       iBooker.book1D("h_LAbySector_Rejected", fmt::sprintf(repText, "rejected", "rejected"), maxSect, lo, hi);
0333   hists.h_bySectLA_ = iBooker.book1D("h_LAbySector", fmt::sprintf(repText, "payload", "payload"), maxSect, lo, hi);
0334   hists.h_bySectDeltaLA_ =
0335       iBooker.book1D("h_deltaLAbySector", fmt::sprintf(repText, "#Delta", "#Delta"), maxSect, lo, hi);
0336   hists.h_bySectChi2_ =
0337       iBooker.book1D("h_bySectorChi2", "Fit #chi^{2}/ndf by sector;pixel sector; fit #chi^{2}/ndf", maxSect, lo, hi);
0338 
0339   // copy the bin labels from the occupancy histogram
0340   for (int bin = 1; bin <= maxSect; bin++) {
0341     const auto& binName = hists.h_bySectOccupancy_->getTH1()->GetXaxis()->GetBinLabel(bin);
0342     hists.h_bySectMeasLA_->setBinLabel(bin, binName);
0343     hists.h_bySectSetLA_->setBinLabel(bin, binName);
0344     hists.h_bySectRejectLA_->setBinLabel(bin, binName);
0345     hists.h_bySectLA_->setBinLabel(bin, binName);
0346     hists.h_bySectDeltaLA_->setBinLabel(bin, binName);
0347     hists.h_bySectChi2_->setBinLabel(bin, binName);
0348   }
0349 
0350   // this will be booked in the Harvesting folder
0351   iBooker.setCurrentFolder(fmt::format("{}Harvesting/LorentzAngleMaps", dqmDir_));
0352   for (int i = 0; i < hists.nlay; i++) {
0353     std::string repName = "h2_byLayerLA_%i";
0354     std::string repText = "BPix Layer %i tan#theta_{LA}/B;module number;ladder number;tan#theta_{LA}/B [1/T]";
0355 
0356     hists.h2_byLayerLA_.emplace_back(iBooker.book2D(fmt::sprintf(repName, i + 1),
0357                                                     fmt::sprintf(repText, i + 1),
0358                                                     hists.nModules_[i],
0359                                                     0.5,
0360                                                     hists.nModules_[i] + 0.5,
0361                                                     hists.nLadders_[i],
0362                                                     0.5,
0363                                                     hists.nLadders_[i] + 0.5));
0364 
0365     repName = "h2_byLayerDiff_%i";
0366     repText = "BPix Layer %i #Delta#mu_{H}/#mu_{H};module number;ladder number;#Delta#mu_{H}/#mu_{H} [%%]";
0367     hists.h2_byLayerDiff_.emplace_back(iBooker.book2D(fmt::sprintf(repName, i + 1),
0368                                                       fmt::sprintf(repText, i + 1),
0369                                                       hists.nModules_[i],
0370                                                       0.5,
0371                                                       hists.nModules_[i] + 0.5,
0372                                                       hists.nLadders_[i],
0373                                                       0.5,
0374                                                       hists.nLadders_[i] + 0.5));
0375   }
0376 
0377   // clang-format off
0378   edm::LogPrint("LorentzAngle") << "module" << "\t" << "layer" << "\t"
0379                                 << "offset" << "\t" << "e0" << "\t"
0380                                 << "slope"  << "\t" << "e1" << "\t"
0381                                 << "rel.err" << "\t" << "pull" << "\t"
0382                                 << "p2" << "\t" << "e2" << "\t"
0383                                 << "p3" << "\t" << "e3" << "\t"
0384                                 << "p4" << "\t" << "e4" << "\t"
0385                                 << "p5" << "\t" << "e5" << "\t"
0386                                 << "chi2" << "\t" << "prob" << "\t"
0387                                 << "newDetId" << "\t" << "tan(LA)" << "\t"
0388                                 << "Error(LA)" ;
0389   // clang-format on
0390 
0391   // payload to be written out
0392   std::shared_ptr<SiPixelLorentzAngle> LorentzAngle = std::make_shared<SiPixelLorentzAngle>();
0393 
0394   // fill the map of simulation values
0395   double p1_simul_newmodule = 0.294044;
0396   double p1_simul[hists.nlay + 1][hists.nModules_[hists.nlay - 1]];
0397   for (int i_layer = 1; i_layer <= hists.nlay; i_layer++) {
0398     for (int i_module = 1; i_module <= hists.nModules_[i_layer - 1]; i_module++) {
0399       if (i_layer == 1)
0400         p1_simul[i_layer - 1][i_module - 1] = 0.436848;
0401       else if (i_layer == 2)
0402         p1_simul[i_layer - 1][i_module - 1] = 0.25802;
0403       else if (i_layer == 3 && i_module <= 4)
0404         p1_simul[i_layer - 1][i_module - 1] = 0.29374;
0405       else if (i_layer == 3 && i_module >= 5)
0406         p1_simul[i_layer - 1][i_module - 1] = 0.31084;
0407       else if (i_layer == 4 && i_module <= 4)
0408         p1_simul[i_layer - 1][i_module - 1] = 0.29944;
0409       else
0410         p1_simul[i_layer - 1][i_module - 1] = 0.31426;
0411     }
0412   }
0413   // fictitious n-th layer to store the values of new modules
0414   for (int i_module = 1; i_module <= hists.nModules_[hists.nlay - 1]; i_module++) {
0415     p1_simul[hists.nlay][i_module - 1] = p1_simul_newmodule;
0416   }
0417 
0418   // loop over "new" BPix modules
0419   for (int j = 0; j < (int)hists.BPixnewDetIds_.size(); j++) {
0420     //uint32_t rawId = hists.BPixnewDetIds_[j];
0421     int new_index = j + 1 + hists.nModules_[hists.nlay - 1] + (hists.nlay - 1) * hists.nModules_[hists.nlay - 1];
0422     if (hists.h_drift_depth_adc_[new_index] == nullptr)
0423       continue;
0424     for (int i = 1; i <= hist_depth_; i++) {
0425       findMean(h_drift_depth_adc_slice_, i, new_index);
0426     }
0427 
0428     // fit the distributions and store the LA in the payload
0429     const auto& res = fitAndStore(LorentzAngle, new_index, hists.BPixnewLayer_[j], hists.BPixnewModule_[j]);
0430 
0431     edm::LogPrint("SiPixelLorentzAngle") << std::setprecision(4) << hists.BPixnewModule_[j] << "\t"
0432                                          << hists.BPixnewLayer_[j] << "\t" << res.p0 << "\t" << res.e0 << "\t" << res.p1
0433                                          << std::setprecision(3) << "\t" << res.e1 << "\t" << res.e1 / res.p1 * 100.
0434                                          << "\t" << (res.p1 - p1_simul[hists.nlay][0]) / res.e1 << "\t" << res.p2
0435                                          << "\t" << res.e2 << "\t" << res.p3 << "\t" << res.e3 << "\t" << res.p4 << "\t"
0436                                          << res.e4 << "\t" << res.p5 << "\t" << res.e5 << "\t" << res.chi2 << "\t"
0437                                          << res.prob << "\t" << hists.BPixnewDetIds_[j] << "\t" << res.tan_LA << "\t"
0438                                          << res.error_LA;
0439   }  // loop on BPix new modules
0440 
0441   //loop over modules and layers to fit the lorentz angle
0442   for (int i_layer = 1; i_layer <= hists.nlay; i_layer++) {
0443     for (int i_module = 1; i_module <= hists.nModules_[i_layer - 1]; i_module++) {
0444       int i_index = i_module + (i_layer - 1) * hists.nModules_[i_layer - 1];
0445       if (hists.h_drift_depth_adc_[i_index] == nullptr)
0446         continue;
0447       //loop over bins in depth (z-local-coordinate) (in order to fit slices)
0448       for (int i = 1; i <= hist_depth_; i++) {
0449         findMean(h_drift_depth_adc_slice_, i, i_index);
0450       }  // end loop over bins in depth
0451 
0452       // fit the distributions and store the LA in the payload
0453       const auto& res = fitAndStore(LorentzAngle, i_index, i_layer, i_module);
0454 
0455       edm::LogPrint("SiPixelLorentzAngle")
0456           << std::setprecision(4) << i_module << "\t" << i_layer << "\t" << res.p0 << "\t" << res.e0 << "\t" << res.p1
0457           << std::setprecision(3) << "\t" << res.e1 << "\t" << res.e1 / res.p1 * 100. << "\t"
0458           << (res.p1 - p1_simul[i_layer - 1][i_module - 1]) / res.e1 << "\t" << res.p2 << "\t" << res.e2 << "\t"
0459           << res.p3 << "\t" << res.e3 << "\t" << res.p4 << "\t" << res.e4 << "\t" << res.p5 << "\t" << res.e5 << "\t"
0460           << res.chi2 << "\t" << res.prob << "\t"
0461           << "null"
0462           << "\t" << res.tan_LA << "\t" << res.error_LA;
0463     }
0464   }  // end loop over modules and layers
0465 
0466   // fill the rest of DetIds not filled above (for the moment FPix)
0467   const auto& currentLAMap = currentLorentzAngle->getLorentzAngles();
0468   const auto& newLAMap = LorentzAngle->getLorentzAngles();
0469   std::vector<unsigned int> currentLADets;
0470   std::vector<unsigned int> newLADets;
0471 
0472   std::transform(currentLAMap.begin(),
0473                  currentLAMap.end(),
0474                  std::back_inserter(currentLADets),
0475                  [](const std::map<unsigned int, float>::value_type& pair) { return pair.first; });
0476 
0477   std::transform(newLAMap.begin(),
0478                  newLAMap.end(),
0479                  std::back_inserter(newLADets),
0480                  [](const std::map<unsigned int, float>::value_type& pair) { return pair.first; });
0481 
0482   std::vector<unsigned int> notCommon;
0483   std::set_symmetric_difference(
0484       currentLADets.begin(), currentLADets.end(), newLADets.begin(), newLADets.end(), std::back_inserter(notCommon));
0485 
0486   for (const auto& id : notCommon) {
0487     float fPixLorentzAnglePerTesla_ = currentLorentzAngle->getLorentzAngle(id);
0488     if (!LorentzAngle->putLorentzAngle(id, fPixLorentzAnglePerTesla_)) {
0489       edm::LogError("SiPixelLorentzAnglePCLHarvester")
0490           << "[SiPixelLorentzAnglePCLHarvester::dqmEndJob] filling rest of payload: detid already exists";
0491     }
0492   }
0493 
0494   for (const auto& id : newLADets) {
0495     float deltaMuHoverMuH = (currentLorentzAngle->getLorentzAngle(id) - LorentzAngle->getLorentzAngle(id)) /
0496                             currentLorentzAngle->getLorentzAngle(id);
0497     h_diffLA->Fill(deltaMuHoverMuH * 100.f);
0498   }
0499 
0500   bool isPayloadChanged{false};
0501   // fill the 2D output Lorentz Angle maps and check if the payload is different from the input one
0502   for (const auto& [id, value] : LorentzAngle->getLorentzAngles()) {
0503     DetId ID = DetId(id);
0504     if (ID.subdetId() == PixelSubdetector::PixelBarrel) {
0505       const auto& layer = theTrackerTopology->pxbLayer(id);
0506       const auto& ladder = theTrackerTopology->pxbLadder(id);
0507       const auto& module = theTrackerTopology->pxbModule(id);
0508       hists.h2_byLayerLA_[layer - 1]->setBinContent(module, ladder, value);
0509 
0510       float deltaMuHoverMuH =
0511           (currentLorentzAngle->getLorentzAngle(id) - value) / currentLorentzAngle->getLorentzAngle(id);
0512       hists.h2_byLayerDiff_[layer - 1]->setBinContent(module, ladder, deltaMuHoverMuH * 100.f);
0513 
0514       if (!isPayloadChanged && (deltaMuHoverMuH != 0.f))
0515         isPayloadChanged = true;
0516     }
0517   }
0518 
0519   if (isPayloadChanged) {
0520     // fill the DB object record
0521     edm::Service<cond::service::PoolDBOutputService> mydbservice;
0522     if (mydbservice.isAvailable()) {
0523       try {
0524         mydbservice->writeOneIOV(*LorentzAngle, mydbservice->currentTime(), recordName_);
0525       } catch (const cond::Exception& er) {
0526         edm::LogError("SiPixelLorentzAngleDB") << er.what();
0527       } catch (const std::exception& er) {
0528         edm::LogError("SiPixelLorentzAngleDB") << "caught std::exception " << er.what();
0529       }
0530     } else {
0531       edm::LogError("SiPixelLorentzAngleDB") << "Service is unavailable";
0532     }
0533   } else {
0534     edm::LogPrint("SiPixelLorentzAngleDB") << __PRETTY_FUNCTION__ << " there is no new valid measurement to append! ";
0535   }
0536 }
0537 
0538 //------------------------------------------------------------------------------
0539 void SiPixelLorentzAnglePCLHarvester::findMean(MonitorElement* h_drift_depth_adc_slice_, int i, int i_ring) {
0540   double nentries = 0;
0541   h_drift_depth_adc_slice_->Reset();
0542   int hist_drift_ = h_drift_depth_adc_slice_->getNbinsX();
0543 
0544   // determine sigma and sigma^2 of the adc counts and average adc counts
0545   //loop over bins in drift width
0546   for (int j = 1; j <= hist_drift_; j++) {
0547     if (hists.h_drift_depth_noadc_[i_ring]->getBinContent(j, i) >= 1) {
0548       double adc_error2 = (hists.h_drift_depth_adc2_[i_ring]->getBinContent(j, i) -
0549                            hists.h_drift_depth_adc_[i_ring]->getBinContent(j, i) *
0550                                hists.h_drift_depth_adc_[i_ring]->getBinContent(j, i) /
0551                                hists.h_drift_depth_noadc_[i_ring]->getBinContent(j, i)) /
0552                           hists.h_drift_depth_noadc_[i_ring]->getBinContent(j, i);
0553 
0554       hists.h_drift_depth_adc_[i_ring]->setBinError(j, i, sqrt(adc_error2));
0555       double error2 = adc_error2 / (hists.h_drift_depth_noadc_[i_ring]->getBinContent(j, i) - 1.);
0556       hists.h_drift_depth_[i_ring]->setBinError(j, i, sqrt(error2));
0557     } else {
0558       hists.h_drift_depth_[i_ring]->setBinError(j, i, 0);
0559       hists.h_drift_depth_adc_[i_ring]->setBinError(j, i, 0);
0560     }
0561     h_drift_depth_adc_slice_->setBinContent(j, hists.h_drift_depth_adc_[i_ring]->getBinContent(j, i));
0562     h_drift_depth_adc_slice_->setBinError(j, hists.h_drift_depth_adc_[i_ring]->getBinError(j, i));
0563     nentries += hists.h_drift_depth_noadc_[i_ring]->getBinContent(j, i);
0564   }  // end loop over bins in drift width
0565 
0566   double mean = h_drift_depth_adc_slice_->getMean(1);
0567   double error = 0;
0568   if (nentries != 0) {
0569     error = h_drift_depth_adc_slice_->getRMS(1) / std::sqrt(nentries);
0570   }
0571   hists.h_mean_[i_ring]->setBinContent(i, mean);
0572   hists.h_mean_[i_ring]->setBinError(i, error);
0573 
0574   h_drift_depth_adc_slice_->Reset();  // clear again after extracting the parameters
0575 }
0576 
0577 //------------------------------------------------------------------------------
0578 SiPixelLAHarvest::fitResults SiPixelLorentzAnglePCLHarvester::fitAndStore(
0579     std::shared_ptr<SiPixelLorentzAngle> theLAPayload, int i_index, int i_layer, int i_module) {
0580   // output results
0581   SiPixelLAHarvest::fitResults res;
0582 
0583   // B-field value
0584   // nominalValue returns the magnetic field value in kgauss (1T = 10 kgauss)
0585   float theMagField = magField->nominalValue() / 10.;
0586 
0587   double half_width = width_ * 10000 / 2;  // pixel half thickness in units of micro meter
0588 
0589   f1 = std::make_unique<TF1>("f1", "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x + [5]*x*x*x*x*x", 5., 280.);
0590   f1->SetParName(0, "offset");
0591   f1->SetParName(1, "tan#theta_{LA}");
0592   f1->SetParName(2, "quad term");
0593   f1->SetParName(3, "cubic term");
0594   f1->SetParName(4, "quartic term");
0595   f1->SetParName(5, "quintic term");
0596 
0597   f1->SetParameter(0, 0);
0598   f1->SetParError(0, 0);
0599   f1->SetParameter(1, 0.4);
0600   f1->SetParError(1, 0);
0601   f1->SetParameter(2, 0.0);
0602   f1->SetParError(2, 0);
0603   f1->SetParameter(3, 0.0);
0604   f1->SetParError(3, 0);
0605   f1->SetParameter(4, 0.0);
0606   f1->SetParError(4, 0);
0607   f1->SetParameter(5, 0.0);
0608   f1->SetParError(5, 0);
0609   f1->SetChisquare(0);
0610 
0611   hists.h_mean_[i_index]->getTH1()->Fit(f1.get(), "ERQ");
0612 
0613   res.p0 = f1->GetParameter(0);
0614   res.e0 = f1->GetParError(0);
0615   res.p1 = f1->GetParameter(1);
0616   res.e1 = f1->GetParError(1);
0617   res.p2 = f1->GetParameter(2);
0618   res.e2 = f1->GetParError(2);
0619   res.p3 = f1->GetParameter(3);
0620   res.e3 = f1->GetParError(3);
0621   res.p4 = f1->GetParameter(4);
0622   res.e4 = f1->GetParError(4);
0623   res.p5 = f1->GetParameter(5);
0624   res.e5 = f1->GetParError(5);
0625   res.chi2 = f1->GetChisquare();
0626   res.ndf = f1->GetNDF();
0627   res.prob = f1->GetProb();
0628   res.redChi2 = res.ndf > 0. ? res.chi2 / res.ndf : 0.;
0629 
0630   double f1_halfwidth = res.p0 + res.p1 * half_width + res.p2 * pow(half_width, 2) + res.p3 * pow(half_width, 3) +
0631                         res.p4 * pow(half_width, 4) + res.p5 * pow(half_width, 5);
0632 
0633   double f1_zerowidth = res.p0;
0634 
0635   // tan_LA = (f1(x = half_width) - f1(x = 0)) / (half_width - 0)
0636   res.tan_LA = (f1_halfwidth - f1_zerowidth) / half_width;
0637   double errsq_LA =
0638       (pow(res.e1, 2) + pow((half_width * res.e2), 2) + pow((half_width * half_width * res.e3), 2) +
0639        pow((half_width * half_width * half_width * res.e4), 2) +
0640        pow((half_width * half_width * half_width * half_width * res.e5), 2));  // Propagation of uncertainty
0641   res.error_LA = sqrt(errsq_LA);
0642 
0643   hists.h_bySectMeasLA_->setBinContent(i_index, (res.tan_LA / theMagField));
0644   hists.h_bySectMeasLA_->setBinError(i_index, (res.error_LA / theMagField));
0645   hists.h_bySectChi2_->setBinContent(i_index, res.redChi2);
0646   hists.h_bySectChi2_->setBinError(i_index, 0.);  // no errors
0647 
0648   int nentries = hists.h_bySectOccupancy_->getBinContent(i_index);  // number of on track hits in that sector
0649 
0650   bool isNew = (i_index > hists.nlay * hists.nModules_[hists.nlay - 1]);
0651   int shiftIdx = i_index - hists.nlay * hists.nModules_[hists.nlay - 1] - 1;
0652 
0653   LogDebug("SiPixelLorentzAnglePCLHarvester")
0654       << " isNew: " << isNew << " i_index: " << i_index << " shift index: " << shiftIdx;
0655 
0656   const auto& detIdsToFill =
0657       isNew ? std::vector<unsigned int>({hists.BPixnewDetIds_[shiftIdx]}) : hists.detIdsList[i_index];
0658 
0659   LogDebug("SiPixelLorentzAnglePCLHarvester")
0660       << "index: " << i_index << " i_module: " << i_module << " i_layer: " << i_layer;
0661   for (const auto& id : detIdsToFill) {
0662     LogDebug("SiPixelLorentzAnglePCLHarvester") << id << ",";
0663   }
0664 
0665   // no errors on the following MEs
0666   hists.h_bySectSetLA_->setBinError(i_index, 0.);
0667   hists.h_bySectRejectLA_->setBinError(i_index, 0.);
0668   hists.h_bySectLA_->setBinError(i_index, 0.);
0669   hists.h_bySectDeltaLA_->setBinError(i_index, 0.);
0670 
0671   float LorentzAnglePerTesla_;
0672   float currentLA = currentLorentzAngle->getLorentzAngle(detIdsToFill.front());
0673   // if the fit quality is OK
0674   if ((res.redChi2 != 0.) && (res.redChi2 < fitChi2Cut_) && (nentries > minHitsCut_)) {
0675     LorentzAnglePerTesla_ = res.tan_LA / theMagField;
0676     // fill the LA actually written to payload
0677     hists.h_bySectSetLA_->setBinContent(i_index, LorentzAnglePerTesla_);
0678     hists.h_bySectRejectLA_->setBinContent(i_index, 0.);
0679     hists.h_bySectLA_->setBinContent(i_index, LorentzAnglePerTesla_);
0680 
0681     const auto& deltaLA = (LorentzAnglePerTesla_ - currentLA);
0682     hists.h_bySectDeltaLA_->setBinContent(i_index, deltaLA);
0683 
0684     for (const auto& id : detIdsToFill) {
0685       if (!theLAPayload->putLorentzAngle(id, LorentzAnglePerTesla_)) {
0686         edm::LogError("SiPixelLorentzAnglePCLHarvester") << "[SiPixelLorentzAnglePCLHarvester::fitAndStore]: detid ("
0687                                                          << i_layer << "," << i_module << ") already exists";
0688       }
0689     }
0690   } else {
0691     // just copy the values from the existing payload
0692     hists.h_bySectSetLA_->setBinContent(i_index, 0.);
0693     hists.h_bySectRejectLA_->setBinContent(i_index, (res.tan_LA / theMagField));
0694     hists.h_bySectLA_->setBinContent(i_index, currentLA);
0695     hists.h_bySectDeltaLA_->setBinContent(i_index, 0.);
0696 
0697     for (const auto& id : detIdsToFill) {
0698       LorentzAnglePerTesla_ = currentLorentzAngle->getLorentzAngle(id);
0699       if (!theLAPayload->putLorentzAngle(id, LorentzAnglePerTesla_)) {
0700         edm::LogError("SiPixelLorentzAnglePCLHarvester") << "[SiPixelLorentzAnglePCLHarvester::fitAndStore]: detid ("
0701                                                          << i_layer << "," << i_module << ") already exists";
0702       }
0703     }
0704   }
0705   // return the struct of fit details
0706   return res;
0707 }
0708 
0709 //------------------------------------------------------------------------------
0710 void SiPixelLorentzAnglePCLHarvester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0711   edm::ParameterSetDescription desc;
0712   desc.setComment("Harvester module of the SiPixel Lorentz Angle PCL monitoring workflow");
0713   desc.add<std::vector<std::string>>("newmodulelist", {})->setComment("the list of DetIds for new sensors");
0714   desc.add<std::string>("dqmDir", "AlCaReco/SiPixelLorentzAngle")->setComment("the directory of PCL Worker output");
0715   desc.add<double>("fitChi2Cut", 20.)->setComment("cut on fit chi2/ndof to accept measurement");
0716   desc.add<int>("minHitsCut", 10000)->setComment("cut on minimum number of on-track hits to accept measurement");
0717   desc.add<std::string>("record", "SiPixelLorentzAngleRcd")->setComment("target DB record");
0718   descriptions.addWithDefaultLabel(desc);
0719 }
0720 
0721 DEFINE_FWK_MODULE(SiPixelLorentzAnglePCLHarvester);