Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-08-21 22:41:46

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