File indexing completed on 2023-10-25 09:35:36
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
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
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
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
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
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
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
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
0124
0125
0126
0127
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
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
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
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
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
0226 for (const auto& ME : iHists_.h2_) {
0227 if (!ME.second)
0228 continue;
0229
0230 ME.second->setOption("colz");
0231 TProfile* hp = (TProfile*)ME.second->getTH2F()->ProfileX();
0232 iBooker.setCurrentFolder(folderToHarvest + "/" + getStem(ME.first, 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
0246 for (const auto& prof : iHists_.p_) {
0247
0248 if ((prof.first).find("thetatrk_nstrip") != std::string::npos) {
0249
0250
0251
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
0262 prof.second->getTProfile()->Fit(fitFunc, "F");
0263
0264
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, 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
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
0326 std::istringstream iss(histoName);
0327
0328 std::string token;
0329 while (std::getline(iss, token, '_')) {
0330
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);