File indexing completed on 2024-09-07 04:35:10
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 <numeric>
0021 #include <sstream>
0022 #include <string>
0023 #include <vector>
0024
0025
0026 #include "CalibTracker/Records/interface/SiStripDependentRecords.h"
0027 #include "CalibTracker/SiStripCommon/interface/TkDetMap.h"
0028 #include "CalibTracker/SiStripLorentzAngle/interface/SiStripLorentzAngleCalibrationHelpers.h"
0029 #include "CalibTracker/SiStripLorentzAngle/interface/SiStripLorentzAngleCalibrationStruct.h"
0030 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0031 #include "CondFormats/SiStripObjects/interface/SiStripLorentzAngle.h"
0032 #include "DQM/SiStripCommon/interface/TkHistoMap.h"
0033 #include "DQMServices/Core/interface/DQMEDHarvester.h"
0034 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0035 #include "FWCore/Framework/interface/EventSetup.h"
0036 #include "FWCore/Framework/interface/Frameworkfwd.h"
0037 #include "FWCore/Framework/interface/MakerMacros.h"
0038 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0039 #include "FWCore/ServiceRegistry/interface/Service.h"
0040 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0041 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0042 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0043 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0044 #include "MagneticField/Engine/interface/MagneticField.h"
0045 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0046
0047
0048 class SiStripLorentzAnglePCLHarvester : public DQMEDHarvester {
0049 public:
0050 SiStripLorentzAnglePCLHarvester(const edm::ParameterSet&);
0051 ~SiStripLorentzAnglePCLHarvester() override = default;
0052 void beginRun(const edm::Run&, const edm::EventSetup&) override;
0053
0054 static void fillDescriptions(edm::ConfigurationDescriptions&);
0055
0056 private:
0057 void dqmEndJob(DQMStore::IBooker&, DQMStore::IGetter&) override;
0058 void endRun(const edm::Run&, const edm::EventSetup&) override;
0059 std::string getStem(const std::string& histoName, bool isFolder);
0060
0061
0062 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomEsToken_;
0063 const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoEsTokenBR_, topoEsTokenER_;
0064 const edm::ESGetToken<TkDetMap, TrackerTopologyRcd> tkDetMapTokenER_;
0065
0066 const edm::ESGetToken<SiStripLatency, SiStripLatencyRcd> latencyToken_;
0067 const edm::ESGetToken<SiStripLorentzAngle, SiStripLorentzAngleDepRcd> siStripLAEsToken_;
0068 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
0069
0070
0071
0072 bool mismatchedBField_;
0073 bool mismatchedLatency_;
0074
0075 const bool debug_;
0076 SiStripLorentzAngleCalibrationHistograms iHists_;
0077
0078 const std::string dqmDir_;
0079 const std::vector<double> fitRange_;
0080 const std::string recordName_;
0081 float theMagField_{0.f};
0082
0083 static constexpr float teslaToInverseGeV_ = 2.99792458e-3f;
0084 std::pair<double, double> theFitRange_{0., 0.};
0085
0086 const SiStripLorentzAngle* currentLorentzAngle_;
0087 std::unique_ptr<TrackerTopology> theTrackerTopology_;
0088 std::unique_ptr<TkDetMap> tkDetMap_;
0089
0090
0091 struct LATkMap {
0092 LATkMap() : hInputLA(nullptr), hOutputLA(nullptr) {}
0093 LATkMap(std::unique_ptr<TkHistoMap>&& input, std::unique_ptr<TkHistoMap>&& output)
0094 : hInputLA(std::move(input)), hOutputLA(std::move(output)) {}
0095
0096 void fill(uint32_t id, float inputLA, float outputLA) {
0097 hInputLA->fill(id, inputLA);
0098 hOutputLA->fill(id, outputLA);
0099 }
0100
0101 void beautify() {
0102 const auto& allInputMaps = hInputLA->getAllMaps();
0103
0104 for (size_t i = 1; i < allInputMaps.size(); i++) {
0105 hInputLA->getMap(i)->setOption("colz");
0106 hOutputLA->getMap(i)->setOption("colz");
0107 }
0108 }
0109
0110 std::unique_ptr<TkHistoMap> hInputLA, hOutputLA;
0111 };
0112
0113 LATkMap h_modulesLA;
0114 };
0115
0116
0117 SiStripLorentzAnglePCLHarvester::SiStripLorentzAnglePCLHarvester(const edm::ParameterSet& iConfig)
0118 : geomEsToken_(esConsumes<edm::Transition::BeginRun>()),
0119 topoEsTokenBR_(esConsumes<edm::Transition::BeginRun>()),
0120 topoEsTokenER_(esConsumes<edm::Transition::EndRun>()),
0121 tkDetMapTokenER_(esConsumes<edm::Transition::EndRun>()),
0122 latencyToken_(esConsumes<edm::Transition::BeginRun>()),
0123 siStripLAEsToken_(esConsumes<edm::Transition::BeginRun>()),
0124 magneticFieldToken_(esConsumes<edm::Transition::BeginRun>()),
0125 mismatchedBField_{false},
0126 mismatchedLatency_{false},
0127 debug_(iConfig.getParameter<bool>("debugMode")),
0128 dqmDir_(iConfig.getParameter<std::string>("dqmDir")),
0129 fitRange_(iConfig.getParameter<std::vector<double>>("fitRange")),
0130 recordName_(iConfig.getParameter<std::string>("record")) {
0131
0132 if (fitRange_.size() == 2) {
0133 theFitRange_.first = fitRange_[0];
0134 theFitRange_.second = fitRange_[1];
0135 } else {
0136 throw cms::Exception("SiStripLorentzAnglePCLHarvester") << "Too many fit range parameters specified";
0137 }
0138
0139
0140 edm::Service<cond::service::PoolDBOutputService> poolDbService;
0141 if (!poolDbService.isAvailable())
0142 throw cms::Exception("SiStripLorentzAnglePCLHarvester") << "PoolDBService required";
0143 }
0144
0145
0146 void SiStripLorentzAnglePCLHarvester::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0147
0148 const TrackerGeometry* geom = &iSetup.getData(geomEsToken_);
0149 const TrackerTopology* tTopo = &iSetup.getData(topoEsTokenBR_);
0150
0151 const MagneticField* magField = &iSetup.getData(magneticFieldToken_);
0152 currentLorentzAngle_ = &iSetup.getData(siStripLAEsToken_);
0153
0154
0155
0156
0157
0158
0159
0160 theMagField_ = 1.f / (magField->inverseBzAtOriginInGeV() * teslaToInverseGeV_);
0161
0162 if (iHists_.bfield_.empty()) {
0163 iHists_.bfield_ = siStripLACalibration::fieldAsString(theMagField_);
0164 } else {
0165 if (iHists_.bfield_ != siStripLACalibration::fieldAsString(theMagField_)) {
0166 mismatchedBField_ = true;
0167 }
0168 }
0169
0170 const SiStripLatency* apvlat = &iSetup.getData(latencyToken_);
0171 if (iHists_.apvmode_.empty()) {
0172 iHists_.apvmode_ = siStripLACalibration::apvModeAsString(apvlat);
0173 } else {
0174 if (iHists_.apvmode_ != siStripLACalibration::apvModeAsString(apvlat)) {
0175 mismatchedLatency_ = true;
0176 }
0177 }
0178
0179 auto dets = geom->detsTIB();
0180 dets.insert(dets.end(), geom->detsTID().begin(), geom->detsTID().end());
0181 dets.insert(dets.end(), geom->detsTOB().begin(), geom->detsTOB().end());
0182 dets.insert(dets.end(), geom->detsTEC().begin(), geom->detsTEC().end());
0183
0184 for (auto det : dets) {
0185 auto detid = det->geographicalId().rawId();
0186 const StripGeomDetUnit* stripDet = dynamic_cast<const StripGeomDetUnit*>(geom->idToDet(det->geographicalId()));
0187 if (stripDet) {
0188 iHists_.la_db_[detid] = currentLorentzAngle_->getLorentzAngle(detid);
0189 iHists_.moduleLocationType_[detid] = siStripLACalibration::moduleLocationType(detid, tTopo);
0190 }
0191 }
0192 }
0193
0194
0195 void SiStripLorentzAnglePCLHarvester::endRun(edm::Run const& run, edm::EventSetup const& isetup) {
0196 if (!theTrackerTopology_) {
0197 theTrackerTopology_ = std::make_unique<TrackerTopology>(isetup.getData(topoEsTokenER_));
0198 }
0199 if (!tkDetMap_) {
0200 tkDetMap_ = std::make_unique<TkDetMap>(isetup.getData(tkDetMapTokenER_));
0201 }
0202 }
0203
0204
0205 void SiStripLorentzAnglePCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter) {
0206 if (mismatchedBField_) {
0207 throw cms::Exception("SiStripLorentzAnglePCLHarvester") << "Trying to harvest runs with different B-field values!";
0208 }
0209
0210 if (mismatchedLatency_) {
0211 throw cms::Exception("SiStripLorentzAnglePCLHarvester") << "Trying to harvest runs with diffent APV modes!";
0212 }
0213
0214
0215 iGetter.cd();
0216 std::string bvalue = (iHists_.bfield_ == "3.8") ? "B-ON" : "B-OFF";
0217 std::string folderToHarvest = fmt::format("{}/{}_{}", dqmDir_, bvalue, iHists_.apvmode_);
0218 edm::LogPrint(moduleDescription().moduleName()) << "Harvesting in " << folderToHarvest;
0219 iGetter.setCurrentFolder(folderToHarvest);
0220
0221
0222 iHists_.nlayers_["TIB"] = 4;
0223 iHists_.nlayers_["TOB"] = 6;
0224 iHists_.modtypes_.push_back("s");
0225 iHists_.modtypes_.push_back("a");
0226
0227 std::vector<std::string> MEtoHarvest = {"tanthcosphtrk_nstrip",
0228 "thetatrk_nstrip",
0229 "tanthcosphtrk_var2",
0230 "tanthcosphtrk_var3",
0231 "thcosphtrk_var2",
0232 "thcosphtrk_var3"};
0233
0234
0235 for (auto& layers : iHists_.nlayers_) {
0236 std::string subdet = layers.first;
0237 for (int l = 1; l <= layers.second; ++l) {
0238 for (auto& t : iHists_.modtypes_) {
0239
0240 if (l > 2 && t == "s")
0241 continue;
0242 std::string locationtype = Form("%s_L%d%s", subdet.c_str(), l, t.c_str());
0243 for (const auto& toHarvest : MEtoHarvest) {
0244 const char* address = Form(
0245 "%s/%s/L%d/%s_%s", folderToHarvest.c_str(), subdet.c_str(), l, locationtype.c_str(), toHarvest.c_str());
0246
0247 LogDebug(moduleDescription().moduleName()) << "harvesting at: " << address << std::endl;
0248
0249 iHists_.h2_[Form("%s_%s", locationtype.c_str(), toHarvest.c_str())] = iGetter.get(address);
0250 if (iHists_.h2_[Form("%s_%s", locationtype.c_str(), toHarvest.c_str())] == nullptr) {
0251 edm::LogError(moduleDescription().moduleName())
0252 << "could not retrieve: " << Form("%s_%s", locationtype.c_str(), toHarvest.c_str());
0253 }
0254 }
0255 }
0256 }
0257 }
0258
0259
0260 iBooker.setCurrentFolder(fmt::format("{}Harvesting/LorentzAngleMaps", dqmDir_));
0261
0262
0263 auto sumValues = [](int accumulator, const std::pair<std::string, int>& element) {
0264 return accumulator + element.second;
0265 };
0266
0267
0268 int totalLayers = std::accumulate(iHists_.nlayers_.begin(), iHists_.nlayers_.end(), 0, sumValues);
0269
0270
0271 auto setHistoLabels = [](TH2F* histogram, const std::map<std::string, int>& nlayers) {
0272
0273 histogram->SetOption("colz1");
0274 histogram->SetStats(false);
0275 histogram->GetYaxis()->SetLabelSize(0.05);
0276 histogram->GetXaxis()->SetLabelSize(0.05);
0277
0278
0279 histogram->GetXaxis()->SetBinLabel(1, "r-#phi");
0280 histogram->GetXaxis()->SetBinLabel(2, "stereo");
0281
0282
0283 int binCounter = 1;
0284 for (const auto& subdet : {"TIB", "TOB"}) {
0285 for (int layer = 1; layer <= nlayers.at(subdet); ++layer) {
0286 std::string label = Form("%s L%d", subdet, layer);
0287 histogram->GetYaxis()->SetBinLabel(binCounter++, label.c_str());
0288 }
0289 }
0290 histogram->GetXaxis()->LabelsOption("h");
0291 };
0292
0293 std::string d_name = "h2_byLayerSiStripLA";
0294 std::string d_text = "SiStrip tan#theta_{LA}/B;module type (r-#phi/stereo);layer number;tan#theta_{LA}/B [1/T]";
0295 iHists_.h2_byLayerLA_ =
0296 iBooker.book2D(d_name.c_str(), d_text.c_str(), 2, -0.5, 1.5, totalLayers, -0.5, totalLayers - 0.5);
0297
0298 setHistoLabels(iHists_.h2_byLayerLA_->getTH2F(), iHists_.nlayers_);
0299
0300 d_name = "h2_byLayerSiStripLADiff";
0301 d_text = "SiStrip #Delta#mu_{H}/#mu_{H};module type (r-#phi/stereo);ladder number;#Delta#mu_{H}/#mu_{H} [%%]";
0302 iHists_.h2_byLayerDiff_ =
0303 iBooker.book2D(d_name.c_str(), d_text.c_str(), 2, -0.5, 1.5, totalLayers, -0.5, totalLayers - 0.5);
0304
0305 setHistoLabels(iHists_.h2_byLayerDiff_->getTH2F(), iHists_.nlayers_);
0306
0307
0308 for (const auto& ME : iHists_.h2_) {
0309 if (!ME.second)
0310 continue;
0311
0312 ME.second->setOption("colz");
0313 TProfile* hp = (TProfile*)ME.second->getTH2F()->ProfileX();
0314 iBooker.setCurrentFolder(folderToHarvest + "/" + getStem(ME.first, true));
0315 iHists_.p_[hp->GetName()] = iBooker.bookProfile(hp->GetName(), hp);
0316 iHists_.p_[hp->GetName()]->setAxisTitle(ME.second->getAxisTitle(2), 2);
0317 delete hp;
0318 }
0319
0320 if (iHists_.p_.empty()) {
0321 edm::LogError(moduleDescription().moduleName()) << "None of the input histograms could be retrieved. Aborting";
0322 return;
0323 }
0324
0325 std::map<std::string, std::pair<double, double>> LAMap_;
0326
0327
0328 for (const auto& prof : iHists_.p_) {
0329
0330 if ((prof.first).find("thetatrk_nstrip") != std::string::npos) {
0331
0332
0333
0334 double low = prof.second->getAxisMin(1);
0335 double high = prof.second->getAxisMax(1);
0336 if (theFitRange_.first > theFitRange_.second) {
0337 low = theFitRange_.first;
0338 high = theFitRange_.second;
0339 }
0340
0341 TF1* fitFunc = new TF1("fitFunc", siStripLACalibration::fitFunction, low, high, 3);
0342
0343
0344 prof.second->getTProfile()->Fit(fitFunc, "F");
0345
0346
0347 Double_t a_fit = fitFunc->GetParameter(0);
0348 Double_t thetaL_fit = fitFunc->GetParameter(1);
0349 Double_t b_fit = fitFunc->GetParameter(2);
0350
0351 Double_t a_fitError = fitFunc->GetParError(0);
0352 Double_t thetaL_fitError = fitFunc->GetParError(1);
0353 Double_t b_fitError = fitFunc->GetParError(2);
0354
0355 if (debug_) {
0356 edm::LogInfo(moduleDescription().moduleName())
0357 << prof.first << " fit result: "
0358 << " a=" << a_fit << " +/ " << a_fitError << " theta_L=" << thetaL_fit << " +/- " << thetaL_fitError
0359 << " b=" << b_fit << " +/- " << b_fitError;
0360 }
0361
0362 LAMap_[getStem(prof.first, false)] = std::make_pair(thetaL_fit, thetaL_fitError);
0363 }
0364 }
0365
0366 if (debug_) {
0367 for (const auto& element : LAMap_) {
0368 edm::LogInfo(moduleDescription().moduleName())
0369 << element.first << " thetaLA = " << element.second.first << "+/-" << element.second.second;
0370 }
0371 }
0372
0373
0374 std::shared_ptr<SiStripLorentzAngle> OutLorentzAngle = std::make_shared<SiStripLorentzAngle>();
0375
0376 bool isPayloadChanged{false};
0377 std::vector<std::pair<int, int>> treatedIndices;
0378 for (const auto& loc : iHists_.moduleLocationType_) {
0379 if (debug_) {
0380 edm::LogInfo(moduleDescription().moduleName()) << "modId: " << loc.first << " " << loc.second;
0381 }
0382
0383 if (!(loc.second).empty() && theMagField_ != 0.f) {
0384 OutLorentzAngle->putLorentzAngle(loc.first, std::abs(LAMap_[loc.second].first / theMagField_));
0385 } else {
0386 OutLorentzAngle->putLorentzAngle(loc.first, iHists_.la_db_[loc.first]);
0387 }
0388
0389
0390 if ((loc.second).empty()) {
0391 continue;
0392 }
0393
0394 const auto& index2D = siStripLACalibration::locationTypeIndex(loc.second);
0395 LogDebug("SiStripLorentzAnglePCLHarvester")
0396 << loc.first << " : " << loc.second << " index: " << index2D.first << "-" << index2D.second << std::endl;
0397
0398
0399 if (index2D != std::make_pair(-1, -1)) {
0400
0401
0402 auto it = std::find(treatedIndices.begin(), treatedIndices.end(), index2D);
0403 if (it == treatedIndices.end()) {
0404
0405 LogTrace("SiStripLorentzAnglePCLHarvester") << "accepted " << loc.first << " : " << loc.second << " bin ("
0406 << index2D.first << "," << index2D.second << ")";
0407
0408 const auto& outputLA = OutLorentzAngle->getLorentzAngle(loc.first);
0409 const auto& inputLA = currentLorentzAngle_->getLorentzAngle(loc.first);
0410
0411 LogTrace("SiStripLorentzAnglePCLHarvester") << "inputLA: " << inputLA << " outputLA: " << outputLA;
0412
0413 iHists_.h2_byLayerLA_->setBinContent(index2D.first, index2D.second, outputLA);
0414
0415 float deltaMuHoverMuH = (inputLA != 0.f) ? (inputLA - outputLA) / inputLA : 0.f;
0416 iHists_.h2_byLayerDiff_->setBinContent(index2D.first, index2D.second, deltaMuHoverMuH * 100.f);
0417 treatedIndices.emplace_back(index2D);
0418
0419
0420
0421
0422
0423 if (deltaMuHoverMuH != 0.f && outputLA != 0.f) {
0424 isPayloadChanged = true;
0425 LogDebug("SiStripLorentzAnglePCLHarvester")
0426 << "accepted " << loc.first << " : " << loc.second << " bin (" << index2D.first << "," << index2D.second
0427 << ") " << deltaMuHoverMuH;
0428 } else {
0429 edm::LogWarning("SiStripLorentzAnglePCLHarvester")
0430 << "Discarded " << loc.first << " : " << loc.second << " bin (" << index2D.first << "," << index2D.second
0431 << ") | delta muH/muH = " << deltaMuHoverMuH << " [1/T], output muH = " << outputLA << " [1/T]";
0432 }
0433 }
0434 } else {
0435 throw cms::Exception("SiStripLorentzAnglePCLHarvester")
0436 << "Trying to fill an inexistent module location from " << loc.second << "!";
0437 }
0438 }
0439
0440
0441 const auto tkDetMapFolderIn = (fmt::format("{}Harvesting/TkDetMaps_LAInput", dqmDir_));
0442 const auto tkDetMapFolderOut = (fmt::format("{}Harvesting/TkDetMaps_LAOutput", dqmDir_));
0443 h_modulesLA = LATkMap(
0444 std::make_unique<TkHistoMap>(tkDetMap_.get(), iBooker, tkDetMapFolderIn, "inputLorentzAngle", 0, false, true),
0445 std::make_unique<TkHistoMap>(tkDetMap_.get(), iBooker, tkDetMapFolderOut, "outputLorentzAngle", 0, false, true));
0446
0447
0448 for (const auto& loc : iHists_.moduleLocationType_) {
0449 const auto& outputLA = OutLorentzAngle->getLorentzAngle(loc.first);
0450 const auto& inputLA = currentLorentzAngle_->getLorentzAngle(loc.first);
0451 h_modulesLA.fill(loc.first, inputLA, outputLA);
0452 }
0453
0454
0455 h_modulesLA.beautify();
0456
0457 if (isPayloadChanged) {
0458
0459 edm::Service<cond::service::PoolDBOutputService> mydbservice;
0460 if (mydbservice.isAvailable()) {
0461 try {
0462 mydbservice->writeOneIOV(*OutLorentzAngle, mydbservice->currentTime(), recordName_);
0463 } catch (const cond::Exception& er) {
0464 edm::LogError("SiStripLorentzAngleDB") << er.what();
0465 } catch (const std::exception& er) {
0466 edm::LogError("SiStripLorentzAngleDB") << "caught std::exception " << er.what();
0467 }
0468 } else {
0469 edm::LogError("SiStripLorentzAngleDB") << "Service is unavailable!";
0470 }
0471 } else {
0472 edm::LogPrint("SiStripLorentzAngleDB")
0473 << "****** WARNING ******\n* " << __PRETTY_FUNCTION__
0474 << "\n* There is no new valid measurement to append!\n* Will NOT update the DB!\n*********************";
0475 }
0476 }
0477
0478
0479 std::string SiStripLorentzAnglePCLHarvester::getStem(const std::string& histoName, bool isFolder) {
0480 std::vector<std::string> tokens;
0481
0482 std::string output{};
0483
0484 std::istringstream iss(histoName);
0485
0486 std::string token;
0487 while (std::getline(iss, token, '_')) {
0488
0489 tokens.push_back(token);
0490 }
0491
0492 if (isFolder) {
0493 output = tokens[0] + "/" + tokens[1];
0494 output.pop_back();
0495 } else {
0496 output = tokens[0] + "_" + tokens[1];
0497 }
0498 return output;
0499 }
0500
0501
0502 void SiStripLorentzAnglePCLHarvester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0503 edm::ParameterSetDescription desc;
0504 desc.setComment("Harvester module of the SiStrip Lorentz Angle PCL monitoring workflow");
0505 desc.add<bool>("debugMode", false)->setComment("determines if it's in debug mode");
0506 desc.add<std::string>("dqmDir", "AlCaReco/SiStripLorentzAngle")->setComment("the directory of PCL Worker output");
0507 desc.add<std::vector<double>>("fitRange", {0., 0.})->setComment("range of depths to perform the LA fit");
0508 desc.add<std::string>("record", "SiStripLorentzAngleRcd")->setComment("target DB record");
0509 descriptions.addWithDefaultLabel(desc);
0510 }
0511
0512 DEFINE_FWK_MODULE(SiStripLorentzAnglePCLHarvester);
