SiStripLorentzAngle/plugins/SiStripLorentzAnglePCLHarvester.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    CalibTracker/SiStripLorentzAnglePCLHarvester
0004 // Class:      SiStripLorentzAnglePCLHarvester
0005 //
0006 /**\class SiStripLorentzAnglePCLHarvester SiStripLorentzAnglePCLHarvester.cc CalibTracker/SiStripLorentzAngle/plugins/SiStripLorentzAnglePCLHarvester.cc
0007  Description: reads the intermediate ALCAPROMPT DQMIO-like dataset and performs the fitting of the SiStrip Lorentz Angle in the Prompt Calibration Loop
0008 */
0009 //
0010 // Original Author:  mmusich
0011 //         Created:  Sat, 29 May 2021 14:46:19 GMT
0012 //
0013 //
0014
0015 // system includes
0016 #include <fmt/format.h>
0017 #include <fmt/printf.h>
0018 #include <fstream>
0019 #include <iostream>
0020 #include <sstream>
0021 #include <string>
0022 #include <vector>
0023
0024 // user includes
0025 #include "CalibTracker/Records/interface/SiStripDependentRecords.h"
0026 #include "CalibTracker/SiStripLorentzAngle/interface/SiStripLorentzAngleCalibrationHelpers.h"
0027 #include "CalibTracker/SiStripLorentzAngle/interface/SiStripLorentzAngleCalibrationStruct.h"
0028 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0029 #include "CondFormats/SiStripObjects/interface/SiStripLorentzAngle.h"
0030 #include "DQMServices/Core/interface/DQMEDHarvester.h"
0031 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0032 #include "FWCore/Framework/interface/EventSetup.h"
0033 #include "FWCore/Framework/interface/Frameworkfwd.h"
0034 #include "FWCore/Framework/interface/MakerMacros.h"
0035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0036 #include "FWCore/ServiceRegistry/interface/Service.h"
0037 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0038 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0039 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0040 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0041 #include "MagneticField/Engine/interface/MagneticField.h"
0042 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0043
0044 //------------------------------------------------------------------------------
0045 class SiStripLorentzAnglePCLHarvester : public DQMEDHarvester {
0046 public:
0047   SiStripLorentzAnglePCLHarvester(const edm::ParameterSet&);
0048   ~SiStripLorentzAnglePCLHarvester() override = default;
0049   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0050
0051   static void fillDescriptions(edm::ConfigurationDescriptions&);
0052
0053 private:
0054   void dqmEndJob(DQMStore::IBooker&, DQMStore::IGetter&) override;
0055   void endRun(const edm::Run&, const edm::EventSetup&) override;
0056   std::string getStem(const std::string& histoName, bool isFolder);
0057
0058   // es tokens
0059   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomEsToken_;
0060   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoEsTokenBR_, topoEsTokenER_;
0061   const edm::ESGetToken<SiStripLatency, SiStripLatencyRcd> latencyToken_;
0062   const edm::ESGetToken<SiStripLorentzAngle, SiStripLorentzAngleDepRcd> siStripLAEsToken_;
0063   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
0064
0065   // member data
0066
0067   bool mismatchedBField_;
0068   bool mismatchedLatency_;
0069
0070   const bool debug_;
0071   SiStripLorentzAngleCalibrationHistograms iHists_;
0072
0073   const std::string dqmDir_;
0074   const std::vector<double> fitRange_;
0075   const std::string recordName_;
0076   float theMagField_{0.f};
0077
0078   static constexpr float teslaToInverseGeV_ = 2.99792458e-3f;
0079   std::pair<double, double> theFitRange_{0., 0.};
0080
0081   SiStripLorentzAngleCalibrationHistograms hists_;
0082   const SiStripLorentzAngle* currentLorentzAngle_;
0083   std::unique_ptr<TrackerTopology> theTrackerTopology_;
0084 };
0085
0086 //------------------------------------------------------------------------------
0087 SiStripLorentzAnglePCLHarvester::SiStripLorentzAnglePCLHarvester(const edm::ParameterSet& iConfig)
0088     : geomEsToken_(esConsumes<edm::Transition::BeginRun>()),
0089       topoEsTokenBR_(esConsumes<edm::Transition::BeginRun>()),
0090       topoEsTokenER_(esConsumes<edm::Transition::EndRun>()),
0091       latencyToken_(esConsumes<edm::Transition::BeginRun>()),
0092       siStripLAEsToken_(esConsumes<edm::Transition::BeginRun>()),
0093       magneticFieldToken_(esConsumes<edm::Transition::BeginRun>()),
0094       mismatchedBField_{false},
0095       mismatchedLatency_{false},
0096       debug_(iConfig.getParameter<bool>("debugMode")),
0097       dqmDir_(iConfig.getParameter<std::string>("dqmDir")),
0098       fitRange_(iConfig.getParameter<std::vector<double>>("fitRange")),
0099       recordName_(iConfig.getParameter<std::string>("record")) {
0100   // initialize the fit range
0101   if (fitRange_.size() == 2) {
0102     theFitRange_.first = fitRange_[0];
0103     theFitRange_.second = fitRange_[1];
0104   } else {
0105     throw cms::Exception("SiStripLorentzAnglePCLHarvester") << "Too many fit range parameters specified";
0106   }
0107
0108   // first ensure DB output service is available
0109   edm::Service<cond::service::PoolDBOutputService> poolDbService;
0110   if (!poolDbService.isAvailable())
0111     throw cms::Exception("SiStripLorentzAnglePCLHarvester") << "PoolDBService required";
0112 }
0113
0114 //------------------------------------------------------------------------------
0115 void SiStripLorentzAnglePCLHarvester::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0116   // geometry
0117   const TrackerGeometry* geom = &iSetup.getData(geomEsToken_);
0118   const TrackerTopology* tTopo = &iSetup.getData(topoEsTokenBR_);
0119
0120   const MagneticField* magField = &iSetup.getData(magneticFieldToken_);
0121   currentLorentzAngle_ = &iSetup.getData(siStripLAEsToken_);
0122
0123   // B-field value
0124   // inverseBzAtOriginInGeV() returns the inverse of field z component for this map in GeV
0125   // for the conversion please consult https://github.com/cms-sw/cmssw/blob/master/MagneticField/Engine/src/MagneticField.cc#L17
0126   // theInverseBzAtOriginInGeV = 1.f / (at0z * 2.99792458e-3f);
0127   // ==> at0z = 1.f / (theInverseBzAtOriginInGeV * 2.99792458e-3f)
0128
0129   theMagField_ = 1.f / (magField->inverseBzAtOriginInGeV() * teslaToInverseGeV_);
0130
0131   if (iHists_.bfield_.empty()) {
0132     iHists_.bfield_ = siStripLACalibration::fieldAsString(theMagField_);
0133   } else {
0134     if (iHists_.bfield_ != siStripLACalibration::fieldAsString(theMagField_)) {
0135       mismatchedBField_ = true;
0136     }
0137   }
0138
0139   const SiStripLatency* apvlat = &iSetup.getData(latencyToken_);
0140   if (iHists_.apvmode_.empty()) {
0141     iHists_.apvmode_ = siStripLACalibration::apvModeAsString(apvlat);
0142   } else {
0143     if (iHists_.apvmode_ != siStripLACalibration::apvModeAsString(apvlat)) {
0144       mismatchedLatency_ = true;
0145     }
0146   }
0147
0148   auto dets = geom->detsTIB();
0149   dets.insert(dets.end(), geom->detsTID().begin(), geom->detsTID().end());
0150   dets.insert(dets.end(), geom->detsTOB().begin(), geom->detsTOB().end());
0151   dets.insert(dets.end(), geom->detsTEC().begin(), geom->detsTEC().end());
0152
0153   for (auto det : dets) {
0154     auto detid = det->geographicalId().rawId();
0155     const StripGeomDetUnit* stripDet = dynamic_cast<const StripGeomDetUnit*>(geom->idToDet(det->geographicalId()));
0156     if (stripDet) {
0157       iHists_.la_db_[detid] = currentLorentzAngle_->getLorentzAngle(detid);
0158       iHists_.moduleLocationType_[detid] = siStripLACalibration::moduleLocationType(detid, tTopo);
0159     }
0160   }
0161 }
0162
0163 //------------------------------------------------------------------------------
0164 void SiStripLorentzAnglePCLHarvester::endRun(edm::Run const& run, edm::EventSetup const& isetup) {
0165   if (!theTrackerTopology_) {
0166     theTrackerTopology_ = std::make_unique<TrackerTopology>(isetup.getData(topoEsTokenER_));
0167   }
0168 }
0169
0170 //------------------------------------------------------------------------------
0171 void SiStripLorentzAnglePCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter) {
0172   if (mismatchedBField_) {
0173     throw cms::Exception("SiStripLorentzAnglePCLHarvester") << "Trying to harvest runs with different B-field values!";
0174   }
0175
0176   if (mismatchedLatency_) {
0177     throw cms::Exception("SiStripLorentzAnglePCLHarvester") << "Trying to harvest runs with diffent APV modes!";
0178   }
0179
0180   // go in the right directory
0181   iGetter.cd();
0182   std::string bvalue = (iHists_.bfield_ == "3.8") ? "B-ON" : "B-OFF";
0183   std::string folderToHarvest = fmt::format("{}/{}_{}", dqmDir_, bvalue, iHists_.apvmode_);
0184   edm::LogPrint(moduleDescription().moduleName()) << "Harvesting in " << folderToHarvest;
0185   iGetter.setCurrentFolder(folderToHarvest);
0186
0187   // fill in the module types
0188   iHists_.nlayers_["TIB"] = 4;
0189   iHists_.nlayers_["TOB"] = 6;
0190   iHists_.modtypes_.push_back("s");
0191   iHists_.modtypes_.push_back("a");
0192
0193   std::vector<std::string> MEtoHarvest = {"tanthcosphtrk_nstrip",
0194                                           "thetatrk_nstrip",
0195                                           "tanthcosphtrk_var2",
0196                                           "tanthcosphtrk_var3",
0197                                           "thcosphtrk_var2",
0198                                           "thcosphtrk_var3"};
0199
0200   // prepare the histograms to be harvested
0201   for (auto& layers : iHists_.nlayers_) {
0202     std::string subdet = layers.first;
0203     for (int l = 1; l <= layers.second; ++l) {
0204       for (auto& t : iHists_.modtypes_) {
0205         // do not fill stereo where there aren't
0206         if (l > 2 && t == "s")
0207           continue;
0208         std::string locationtype = Form("%s_L%d%s", subdet.c_str(), l, t.c_str());
0209         for (const auto& toHarvest : MEtoHarvest) {
0210           const char* address = Form(
0211               "%s/%s/L%d/%s_%s", folderToHarvest.c_str(), subdet.c_str(), l, locationtype.c_str(), toHarvest.c_str());
0212
0213           LogDebug(moduleDescription().moduleName()) << "harvesting at: " << address << std::endl;
0214
0215           iHists_.h2_[Form("%s_%s", locationtype.c_str(), toHarvest.c_str())] = iGetter.get(address);
0216           if (iHists_.h2_[Form("%s_%s", locationtype.c_str(), toHarvest.c_str())] == nullptr) {
0217             edm::LogError(moduleDescription().moduleName())
0218                 << "could not retrieve: " << Form("%s_%s", locationtype.c_str(), toHarvest.c_str());
0219           }
0220         }
0221       }
0222     }
0223   }
0224
0225   // prepare the profiles
0226   for (const auto& ME : iHists_.h2_) {
0227     if (!ME.second)
0228       continue;
0229     // draw colored 2D plots
0230     ME.second->setOption("colz");
0231     TProfile* hp = (TProfile*)ME.second->getTH2F()->ProfileX();
0232     iBooker.setCurrentFolder(folderToHarvest + "/" + getStem(ME.first, /* isFolder = */ true));
0233     iHists_.p_[hp->GetName()] = iBooker.bookProfile(hp->GetName(), hp);
0234     iHists_.p_[hp->GetName()]->setAxisTitle(ME.second->getAxisTitle(2), 2);
0235     delete hp;
0236   }
0237
0238   if (iHists_.p_.empty()) {
0239     edm::LogError(moduleDescription().moduleName()) << "None of the input histograms could be retrieved. Aborting";
0240     return;
0241   }
0242
0243   std::map<std::string, std::pair<double, double>> LAMap_;
0244
0245   // do the fits
0246   for (const auto& prof : iHists_.p_) {
0247     //fit only this type of profile
0248     if ((prof.first).find("thetatrk_nstrip") != std::string::npos) {
0249       // Create the TF1 function
0250
0251       // fitting range (take full axis by default)
0252       double low = prof.second->getAxisMin(1);
0253       double high = prof.second->getAxisMax(1);
0254       if (theFitRange_.first > theFitRange_.second) {
0255         low = theFitRange_.first;
0256         high = theFitRange_.second;
0257       }
0258
0259       TF1* fitFunc = new TF1("fitFunc", siStripLACalibration::fitFunction, low, high, 3);
0260
0261       // Fit the function to the data
0262       prof.second->getTProfile()->Fit(fitFunc, "F");  // "F" option performs a least-squares fit
0263
0264       // Get the fit results
0265       Double_t a_fit = fitFunc->GetParameter(0);
0266       Double_t thetaL_fit = fitFunc->GetParameter(1);
0267       Double_t b_fit = fitFunc->GetParameter(2);
0268
0269       Double_t a_fitError = fitFunc->GetParError(0);
0270       Double_t thetaL_fitError = fitFunc->GetParError(1);
0271       Double_t b_fitError = fitFunc->GetParError(2);
0272
0273       if (debug_) {
0274         edm::LogInfo(moduleDescription().moduleName())
0275             << prof.first << " fit result: "
0276             << " a=" << a_fit << " +/ " << a_fitError << " theta_L=" << thetaL_fit << " +/- " << thetaL_fitError
0277             << " b=" << b_fit << " +/- " << b_fitError;
0278       }
0279
0280       LAMap_[getStem(prof.first, /* isFolder = */ false)] = std::make_pair(thetaL_fit, thetaL_fitError);
0281     }
0282   }
0283
0284   if (debug_) {
0285     for (const auto& element : LAMap_) {
0286       edm::LogInfo(moduleDescription().moduleName())
0287           << element.first << " thetaLA = " << element.second.first << "+/-" << element.second.second;
0288     }
0289   }
0290
0291   // now prepare the output LA
0292   std::shared_ptr<SiStripLorentzAngle> OutLorentzAngle = std::make_shared<SiStripLorentzAngle>();
0293
0294   for (const auto& loc : iHists_.moduleLocationType_) {
0295     if (debug_) {
0296       edm::LogInfo(moduleDescription().moduleName()) << "modId: " << loc.first << " " << loc.second;
0297     }
0298
0299     if (!(loc.second).empty()) {
0300       OutLorentzAngle->putLorentzAngle(loc.first, std::abs(LAMap_[loc.second].first / theMagField_));
0301     } else {
0302       OutLorentzAngle->putLorentzAngle(loc.first, iHists_.la_db_[loc.first]);
0303     }
0304   }
0305
0306   edm::Service<cond::service::PoolDBOutputService> mydbservice;
0307   if (mydbservice.isAvailable()) {
0308     try {
0309       mydbservice->writeOneIOV(*OutLorentzAngle, mydbservice->currentTime(), recordName_);
0310     } catch (const cond::Exception& er) {
0311       edm::LogError("SiStripLorentzAngleDB") << er.what();
0312     } catch (const std::exception& er) {
0313       edm::LogError("SiStripLorentzAngleDB") << "caught std::exception " << er.what();
0314     }
0315   } else {
0316     edm::LogError("SiStripLorentzAngleDB") << "Service is unavailable";
0317   }
0318 }
0319
0320 //------------------------------------------------------------------------------
0321 std::string SiStripLorentzAnglePCLHarvester::getStem(const std::string& histoName, bool isFolder) {
0322   std::vector<std::string> tokens;
0323
0324   std::string output{};
0325   // Create a string stream from the input string
0326   std::istringstream iss(histoName);
0327
0328   std::string token;
0329   while (std::getline(iss, token, '_')) {
0330     // Add each token to the vector
0331     tokens.push_back(token);
0332   }
0333
0334   if (isFolder) {
0335     output = tokens[0] + "/" + tokens[1];
0336     output.pop_back();
0337   } else {
0338     output = tokens[0] + "_" + tokens[1];
0339   }
0340   return output;
0341 }
0342
0343 //------------------------------------------------------------------------------
0344 void SiStripLorentzAnglePCLHarvester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0345   edm::ParameterSetDescription desc;
0346   desc.setComment("Harvester module of the SiStrip Lorentz Angle PCL monitoring workflow");
0347   desc.add<bool>("debugMode", false)->setComment("determines if it's in debug mode");
0348   desc.add<std::string>("dqmDir", "AlCaReco/SiStripLorentzAngle")->setComment("the directory of PCL Worker output");
0349   desc.add<std::vector<double>>("fitRange", {0., 0.})->setComment("range of depths to perform the LA fit");
0350   desc.add<std::string>("record", "SiStripLorentzAngleRcd")->setComment("target DB record");
0351   descriptions.addWithDefaultLabel(desc);
0352 }
0353
0354 DEFINE_FWK_MODULE(SiStripLorentzAnglePCLHarvester);