File indexing completed on 2024-09-08 23:51:38
0001
0002
0003
0004
0005
0006
0007
0008 #include "DQM/BeamMonitor/plugins/BeamMonitorBx.h"
0009 #include "FWCore/ServiceRegistry/interface/Service.h"
0010 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0011 #include "DataFormats/Common/interface/View.h"
0012 #include "RecoVertex/BeamSpotProducer/interface/BSFitter.h"
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014 #include "FWCore/Framework/interface/Run.h"
0015 #include "FWCore/Framework/interface/LuminosityBlock.h"
0016 #include "FWCore/Framework/interface/ConsumesCollector.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include <numeric>
0019 #include <cmath>
0020 #include <TMath.h>
0021 #include <iostream>
0022 #include <TStyle.h>
0023
0024 using namespace std;
0025 using namespace edm;
0026 using namespace reco;
0027
0028 void BeamMonitorBx::formatFitTime(char* ts, const time_t& t) {
0029 #define CET (+1)
0030 #define CEST (+2)
0031
0032 tm* ptm;
0033 ptm = gmtime(&t);
0034 sprintf(ts,
0035 "%4d-%02d-%02d %02d:%02d:%02d",
0036 ptm->tm_year,
0037 ptm->tm_mon + 1,
0038 ptm->tm_mday,
0039 (ptm->tm_hour + CEST) % 24,
0040 ptm->tm_min,
0041 ptm->tm_sec);
0042
0043 #ifdef STRIP_TRAILING_BLANKS_IN_TIMEZONE
0044 unsigned int b = strlen(ts);
0045 while (ts[--b] == ' ') {
0046 ts[b] = 0;
0047 }
0048 #endif
0049 }
0050
0051
0052
0053
0054 BeamMonitorBx::BeamMonitorBx(const ParameterSet& ps) : countBx_(0), countEvt_(0), countLumi_(0), resetHistos_(false) {
0055 parameters_ = ps;
0056 monitorName_ = parameters_.getUntrackedParameter<string>("monitorName", "YourSubsystemName");
0057 bsSrc_ = parameters_.getUntrackedParameter<InputTag>("beamSpot");
0058 fitNLumi_ = parameters_.getUntrackedParameter<int>("fitEveryNLumi", -1);
0059 resetFitNLumi_ = parameters_.getUntrackedParameter<int>("resetEveryNLumi", -1);
0060
0061 usesResource("DQMStore");
0062 dbe_ = Service<DQMStore>().operator->();
0063
0064 if (!monitorName_.empty())
0065 monitorName_ = monitorName_ + "/";
0066
0067 theBeamFitter = new BeamFitter(parameters_, consumesCollector());
0068 theBeamFitter->resetTrkVector();
0069 theBeamFitter->resetLSRange();
0070 theBeamFitter->resetRefTime();
0071 theBeamFitter->resetPVFitter();
0072
0073 if (fitNLumi_ <= 0)
0074 fitNLumi_ = 1;
0075 beginLumiOfBSFit_ = endLumiOfBSFit_ = 0;
0076 refBStime[0] = refBStime[1] = 0;
0077 lastlumi_ = 0;
0078 nextlumi_ = 0;
0079 firstlumi_ = 0;
0080 processed_ = false;
0081 countGoodFit_ = 0;
0082 }
0083
0084 BeamMonitorBx::~BeamMonitorBx() { delete theBeamFitter; }
0085
0086
0087 void BeamMonitorBx::beginJob() {
0088 varMap["x0_bx"] = "X_{0}";
0089 varMap["y0_bx"] = "Y_{0}";
0090 varMap["z0_bx"] = "Z_{0}";
0091 varMap["sigmaX_bx"] = "#sigma_{X}";
0092 varMap["sigmaY_bx"] = "#sigma_{Y}";
0093 varMap["sigmaZ_bx"] = "#sigma_{Z}";
0094
0095 varMap1["x"] = "X_{0} [cm]";
0096 varMap1["y"] = "Y_{0} [cm]";
0097 varMap1["z"] = "Z_{0} [cm]";
0098 varMap1["sigmaX"] = "#sigma_{X} [cm]";
0099 varMap1["sigmaY"] = "#sigma_{Y} [cm]";
0100 varMap1["sigmaZ"] = "#sigma_{Z} [cm]";
0101
0102
0103 dbe_->setCurrentFolder(monitorName_ + "FitBx");
0104
0105 BookTables(1, varMap, "");
0106
0107
0108
0109 for (std::map<std::string, std::string>::const_iterator varName = varMap1.begin(); varName != varMap1.end();
0110 ++varName) {
0111 string subDir_ = "FitBx";
0112 subDir_ += "/";
0113 subDir_ += "All_";
0114 subDir_ += varName->first;
0115 dbe_->setCurrentFolder(monitorName_ + subDir_);
0116 }
0117 dbe_->setCurrentFolder(monitorName_ + "FitBx/All_nPVs");
0118 }
0119
0120
0121 void BeamMonitorBx::beginRun(const edm::Run& r, const EventSetup& context) {
0122 ftimestamp = r.beginTime().value();
0123 tmpTime = ftimestamp >> 32;
0124 startTime = refTime = tmpTime;
0125 }
0126
0127
0128 void BeamMonitorBx::beginLuminosityBlock(const LuminosityBlock& lumiSeg, const EventSetup& context) {
0129 int nthlumi = lumiSeg.luminosityBlock();
0130 const edm::TimeValue_t fbegintimestamp = lumiSeg.beginTime().value();
0131 const std::time_t ftmptime = fbegintimestamp >> 32;
0132
0133 if (countLumi_ == 0) {
0134 beginLumiOfBSFit_ = nthlumi;
0135 refBStime[0] = ftmptime;
0136 }
0137 if (beginLumiOfBSFit_ == 0)
0138 beginLumiOfBSFit_ = nextlumi_;
0139
0140 if (nthlumi < nextlumi_)
0141 return;
0142
0143 if (nthlumi > nextlumi_) {
0144 if (countLumi_ != 0 && processed_) {
0145 FitAndFill(lumiSeg, lastlumi_, nextlumi_, nthlumi);
0146 }
0147 nextlumi_ = nthlumi;
0148 edm::LogInfo("LS|BX|BeamMonitorBx") << "Next Lumi to Fit: " << nextlumi_ << endl;
0149 if (refBStime[0] == 0)
0150 refBStime[0] = ftmptime;
0151 }
0152 countLumi_++;
0153 if (processed_)
0154 processed_ = false;
0155 edm::LogInfo("LS|BX|BeamMonitorBx") << "Begin of Lumi: " << nthlumi << endl;
0156 }
0157
0158
0159 void BeamMonitorBx::analyze(const Event& iEvent, const EventSetup& iSetup) {
0160 bool readEvent_ = true;
0161 const int nthlumi = iEvent.luminosityBlock();
0162 if (nthlumi < nextlumi_) {
0163 readEvent_ = false;
0164 edm::LogWarning("BX|BeamMonitorBx") << "Spilt event from previous lumi section!" << endl;
0165 }
0166 if (nthlumi > nextlumi_) {
0167 readEvent_ = false;
0168 edm::LogWarning("BX|BeamMonitorBx") << "Spilt event from next lumi section!!!" << endl;
0169 }
0170
0171 if (readEvent_) {
0172 countEvt_++;
0173 theBeamFitter->readEvent(iEvent);
0174 processed_ = true;
0175 }
0176 }
0177
0178
0179 void BeamMonitorBx::endLuminosityBlock(const LuminosityBlock& lumiSeg, const EventSetup& iSetup) {
0180 int nthlumi = lumiSeg.id().luminosityBlock();
0181 edm::LogInfo("LS|BX|BeamMonitorBx") << "Lumi of the last event before endLuminosityBlock: " << nthlumi << endl;
0182
0183 if (nthlumi < nextlumi_)
0184 return;
0185 const edm::TimeValue_t fendtimestamp = lumiSeg.endTime().value();
0186 const std::time_t fendtime = fendtimestamp >> 32;
0187 tmpTime = refBStime[1] = fendtime;
0188 }
0189
0190
0191 void BeamMonitorBx::BookTables(int nBx, map<string, string>& vMap, string suffix_) {
0192
0193 dbe_->cd(monitorName_ + "FitBx");
0194
0195 for (std::map<std::string, std::string>::const_iterator varName = vMap.begin(); varName != vMap.end(); ++varName) {
0196 string tmpName = varName->first;
0197 if (!suffix_.empty()) {
0198 tmpName += "_";
0199 tmpName += suffix_;
0200 }
0201
0202 hs[tmpName] = dbe_->book2D(tmpName, varName->second, 3, 0, 3, nBx, 0, nBx);
0203 hs[tmpName]->setBinLabel(1, "bx", 1);
0204 hs[tmpName]->setBinLabel(2, varName->second, 1);
0205 hs[tmpName]->setBinLabel(3, "#Delta " + varName->second, 1);
0206 for (int i = 0; i < nBx; i++) {
0207 hs[tmpName]->setBinLabel(i + 1, " ", 2);
0208 }
0209 hs[tmpName]->getTH1()->SetOption("text");
0210 hs[tmpName]->Reset();
0211 }
0212 }
0213
0214
0215 void BeamMonitorBx::BookTrendHistos(
0216 bool plotPV, int nBx, map<string, string>& vMap, string subDir_, const TString& prefix_, const TString& suffix_) {
0217 int nType_ = 2;
0218 if (plotPV)
0219 nType_ = 4;
0220 std::ostringstream ss;
0221 std::ostringstream ss1;
0222 ss << setfill('0') << setw(5) << nBx;
0223 ss1 << nBx;
0224
0225 for (int i = 0; i < nType_; i++) {
0226 for (std::map<std::string, std::string>::const_iterator varName = vMap.begin(); varName != vMap.end(); ++varName) {
0227 string tmpDir_ = subDir_ + "/All_" + varName->first;
0228 dbe_->cd(monitorName_ + tmpDir_);
0229 TString histTitle(varName->first);
0230 TString tmpName;
0231 if (prefix_ != "")
0232 tmpName = prefix_ + "_" + varName->first;
0233 if (suffix_ != "")
0234 tmpName = tmpName + "_" + suffix_;
0235 tmpName = tmpName + "_" + ss.str();
0236
0237 TString histName(tmpName);
0238 string ytitle(varName->second);
0239 string xtitle("");
0240 string options("E1");
0241 bool createHisto = true;
0242 switch (i) {
0243 case 1:
0244 histName.Insert(histName.Index("_bx_", 4), "_time");
0245 xtitle = "Time [UTC] [Bx# " + ss1.str() + "]";
0246 if (ytitle.find("sigma") == string::npos)
0247 histTitle += " coordinate of beam spot vs time (Fit)";
0248 else
0249 histTitle = histTitle.Insert(5, " ") + " of beam spot vs time (Fit)";
0250 break;
0251 case 2:
0252 if (ytitle.find("sigma") == string::npos) {
0253 histName.Insert(0, "PV");
0254 histName.Insert(histName.Index("_bx_", 4), "_lumi");
0255 histTitle.Insert(0, "Avg. ");
0256 histTitle += " position of primary vtx vs lumi";
0257 xtitle = "Lumisection [Bx# " + ss1.str() + "]";
0258 ytitle.insert(0, "PV");
0259 ytitle += " #pm #sigma_{PV";
0260 ytitle += varName->first;
0261 ytitle += "} (cm)";
0262 } else
0263 createHisto = false;
0264 break;
0265 case 3:
0266 if (ytitle.find("sigma") == string::npos) {
0267 histName.Insert(0, "PV");
0268 histName.Insert(histName.Index("_bx_", 4), "_time");
0269 histTitle.Insert(0, "Avg. ");
0270 histTitle += " position of primary vtx vs time";
0271 xtitle = "Time [UTC] [Bx# " + ss1.str() + "]";
0272 ytitle.insert(0, "PV");
0273 ytitle += " #pm #sigma_{PV";
0274 ytitle += varName->first;
0275 ytitle += "} (cm)";
0276 } else
0277 createHisto = false;
0278 break;
0279 default:
0280 histName.Insert(histName.Index("_bx_", 4), "_lumi");
0281 xtitle = "Lumisection [Bx# " + ss1.str() + "]";
0282 if (ytitle.find("sigma") == string::npos)
0283 histTitle += " coordinate of beam spot vs lumi (Fit)";
0284 else
0285 histTitle = histTitle.Insert(5, " ") + " of beam spot vs lumi (Fit)";
0286 break;
0287 }
0288
0289 if (dbe_->get(monitorName_ + tmpDir_ + "/" + string(histName.Data())))
0290 continue;
0291
0292 if (createHisto) {
0293 edm::LogInfo("BX|BeamMonitorBx") << "histName = " << histName << "; histTitle = " << histTitle << std::endl;
0294 hst[histName] = dbe_->book1D(histName, histTitle, 40, 0.5, 40.5);
0295 hst[histName]->getTH1()->SetCanExtend(TH1::kAllAxes);
0296 hst[histName]->setAxisTitle(xtitle, 1);
0297 hst[histName]->setAxisTitle(ytitle, 2);
0298 hst[histName]->getTH1()->SetOption("E1");
0299 if (histName.Contains("time")) {
0300 hst[histName]->getTH1()->SetBins(3600, 0.5, 3600 + 0.5);
0301 hst[histName]->setAxisTimeDisplay(1);
0302 hst[histName]->setAxisTimeFormat("%H:%M:%S", 1);
0303
0304 char eventTime[64];
0305 formatFitTime(eventTime, startTime);
0306 TDatime da(eventTime);
0307 if (debug_) {
0308 edm::LogInfo("BX|BeamMonitorBx") << "TimeOffset = ";
0309 da.Print();
0310 }
0311 hst[histName]->getTH1()->GetXaxis()->SetTimeOffset(da.Convert(kTRUE));
0312 }
0313 }
0314 }
0315 }
0316
0317
0318 dbe_->cd(monitorName_ + subDir_ + "/All_nPVs");
0319 TString histName = "Trending_nPVs_lumi_bx_" + ss.str();
0320 string xtitle = "Lumisection [Bx# " + ss1.str() + "]";
0321
0322 hst[histName] = dbe_->book1D(histName, "Number of Good Reconstructed Vertices", 40, 0.5, 40.5);
0323 hst[histName]->setAxisTitle(xtitle, 1);
0324 hst[histName]->getTH1()->SetCanExtend(TH1::kAllAxes);
0325 hst[histName]->getTH1()->SetOption("E1");
0326 }
0327
0328
0329 void BeamMonitorBx::FitAndFill(const LuminosityBlock& lumiSeg, int& lastlumi, int& nextlumi, int& nthlumi) {
0330 if (nthlumi <= nextlumi)
0331 return;
0332
0333 int currentlumi = nextlumi;
0334 edm::LogInfo("LS|BX|BeamMonitorBx") << "Lumi of the current fit: " << currentlumi << endl;
0335 lastlumi = currentlumi;
0336 endLumiOfBSFit_ = currentlumi;
0337
0338 edm::LogInfo("BX|BeamMonitorBx") << "Time lapsed = " << tmpTime - refTime << std::endl;
0339
0340 if (resetHistos_) {
0341 edm::LogInfo("BX|BeamMonitorBx") << "Resetting Histograms" << endl;
0342 theBeamFitter->resetCutFlow();
0343 resetHistos_ = false;
0344 }
0345
0346 if (fitNLumi_ > 0)
0347 if (currentlumi % fitNLumi_ != 0)
0348 return;
0349
0350 std::pair<int, int> fitLS = theBeamFitter->getFitLSRange();
0351 edm::LogInfo("LS|BX|BeamMonitorBx") << " [Fitter] Do BeamSpot Fit for LS = " << fitLS.first << " to " << fitLS.second
0352 << endl;
0353 edm::LogInfo("LS|BX|BeamMonitorBx") << " [BX] Do BeamSpot Fit for LS = " << beginLumiOfBSFit_ << " to "
0354 << endLumiOfBSFit_ << endl;
0355
0356 edm::LogInfo("BX|BeamMonitorBx") << " [BxDebugTime] refBStime[0] = " << refBStime[0]
0357 << "; address = " << &refBStime[0] << std::endl;
0358 edm::LogInfo("BX|BeamMonitorBx") << " [BxDebugTime] refBStime[1] = " << refBStime[1]
0359 << "; address = " << &refBStime[1] << std::endl;
0360
0361 if (theBeamFitter->runPVandTrkFitter()) {
0362 countGoodFit_++;
0363 edm::LogInfo("BX|BeamMonitorBx") << "Number of good fit = " << countGoodFit_ << endl;
0364 BeamSpotMapBx bsmap = theBeamFitter->getBeamSpotMap();
0365 std::map<int, int> npvsmap = theBeamFitter->getNPVsperBX();
0366 edm::LogInfo("BX|BeamMonitorBx") << "Number of bx = " << bsmap.size() << endl;
0367 if (bsmap.empty())
0368 return;
0369 if (countBx_ < bsmap.size()) {
0370 countBx_ = bsmap.size();
0371 BookTables(countBx_, varMap, "");
0372 BookTables(countBx_, varMap, "all");
0373 for (BeamSpotMapBx::const_iterator abspot = bsmap.begin(); abspot != bsmap.end(); ++abspot) {
0374 int bx = abspot->first;
0375 BookTrendHistos(false, bx, varMap1, "FitBx", "Trending", "bx");
0376 }
0377 }
0378
0379 std::pair<int, int> LSRange = theBeamFitter->getFitLSRange();
0380 char tmpTitle[50];
0381 sprintf(tmpTitle, "%s %i %s %i %s", " [cm] (LS: ", LSRange.first, " to ", LSRange.second, ")");
0382 for (std::map<std::string, std::string>::const_iterator varName = varMap.begin(); varName != varMap.end();
0383 ++varName) {
0384 hs[varName->first]->setTitle(varName->second + " " + tmpTitle);
0385 hs[varName->first]->Reset();
0386 }
0387
0388 if (countGoodFit_ == 1)
0389 firstlumi_ = LSRange.first;
0390
0391 if (resetFitNLumi_ > 0) {
0392 char tmpTitle1[60];
0393 if (countGoodFit_ > 1)
0394 snprintf(tmpTitle1,
0395 sizeof(tmpTitle1),
0396 "%s %i %s %i %s",
0397 " [cm] (LS: ",
0398 firstlumi_,
0399 " to ",
0400 LSRange.second,
0401 ") [weighted average]");
0402 else
0403 snprintf(tmpTitle1, sizeof(tmpTitle1), "%s", "Need at least two fits to calculate weighted average");
0404 for (std::map<std::string, std::string>::const_iterator varName = varMap.begin(); varName != varMap.end();
0405 ++varName) {
0406 TString tmpName = varName->first + "_all";
0407 hs[tmpName]->setTitle(varName->second + " " + tmpTitle1);
0408 hs[tmpName]->Reset();
0409 }
0410 }
0411
0412 int nthBin = countBx_;
0413 for (BeamSpotMapBx::const_iterator abspot = bsmap.begin(); abspot != bsmap.end(); ++abspot, nthBin--) {
0414 reco::BeamSpot bs = abspot->second;
0415 int bx = abspot->first;
0416 int nPVs = npvsmap.find(bx)->second;
0417 edm::LogInfo("BeamMonitorBx") << "Number of PVs of bx#[" << bx << "] = " << nPVs << endl;
0418
0419 FillTables(bx, nthBin, varMap, bs, "");
0420
0421 FillTrendHistos(bx, nPVs, varMap1, bs, "Trending");
0422 }
0423
0424 weight(fbspotMap, bsmap);
0425
0426 nthBin = countBx_;
0427 if (resetFitNLumi_ > 0 && countGoodFit_ > 1) {
0428 for (BeamSpotMapBx::const_iterator abspot = fbspotMap.begin(); abspot != fbspotMap.end(); ++abspot, nthBin--) {
0429 reco::BeamSpot bs = abspot->second;
0430 int bx = abspot->first;
0431 FillTables(bx, nthBin, varMap, bs, "all");
0432 }
0433 }
0434 }
0435
0436
0437
0438 if (resetFitNLumi_ > 0 && currentlumi % resetFitNLumi_ == 0) {
0439 edm::LogInfo("LS|BX|BeamMonitorBx") << "Reset track collection for beam fit!!!" << endl;
0440 resetHistos_ = true;
0441 theBeamFitter->resetTrkVector();
0442 theBeamFitter->resetLSRange();
0443 theBeamFitter->resetRefTime();
0444 theBeamFitter->resetPVFitter();
0445 beginLumiOfBSFit_ = 0;
0446 refBStime[0] = 0;
0447 }
0448 }
0449
0450
0451 void BeamMonitorBx::FillTables(int nthbx, int nthbin_, map<string, string>& vMap, reco::BeamSpot& bs_, string suffix_) {
0452 map<string, pair<double, double> > val_;
0453 val_["x0_bx"] = pair<double, double>(bs_.x0(), bs_.x0Error());
0454 val_["y0_bx"] = pair<double, double>(bs_.y0(), bs_.y0Error());
0455 val_["z0_bx"] = pair<double, double>(bs_.z0(), bs_.z0Error());
0456 val_["sigmaX_bx"] = pair<double, double>(bs_.BeamWidthX(), bs_.BeamWidthXError());
0457 val_["sigmaY_bx"] = pair<double, double>(bs_.BeamWidthY(), bs_.BeamWidthYError());
0458 val_["sigmaZ_bx"] = pair<double, double>(bs_.sigmaZ(), bs_.sigmaZ0Error());
0459
0460 for (std::map<std::string, std::string>::const_iterator varName = vMap.begin(); varName != vMap.end(); ++varName) {
0461 TString tmpName = varName->first;
0462 if (!suffix_.empty())
0463 tmpName += TString("_" + suffix_);
0464 hs[tmpName]->setBinContent(1, nthbin_, nthbx);
0465 hs[tmpName]->setBinContent(2, nthbin_, val_[varName->first].first);
0466 hs[tmpName]->setBinContent(3, nthbin_, val_[varName->first].second);
0467 }
0468 }
0469
0470
0471 void BeamMonitorBx::FillTrendHistos(
0472 int nthBx, int nPVs_, map<string, string>& vMap, reco::BeamSpot& bs_, const TString& prefix_) {
0473 map<TString, pair<double, double> > val_;
0474 val_[TString(prefix_ + "_x")] = pair<double, double>(bs_.x0(), bs_.x0Error());
0475 val_[TString(prefix_ + "_y")] = pair<double, double>(bs_.y0(), bs_.y0Error());
0476 val_[TString(prefix_ + "_z")] = pair<double, double>(bs_.z0(), bs_.z0Error());
0477 val_[TString(prefix_ + "_sigmaX")] = pair<double, double>(bs_.BeamWidthX(), bs_.BeamWidthXError());
0478 val_[TString(prefix_ + "_sigmaY")] = pair<double, double>(bs_.BeamWidthY(), bs_.BeamWidthYError());
0479 val_[TString(prefix_ + "_sigmaZ")] = pair<double, double>(bs_.sigmaZ(), bs_.sigmaZ0Error());
0480
0481 std::ostringstream ss;
0482 ss << setfill('0') << setw(5) << nthBx;
0483 int ntbin_ = tmpTime - startTime;
0484 for (map<TString, MonitorElement*>::iterator itHst = hst.begin(); itHst != hst.end(); ++itHst) {
0485 if (!(itHst->first.Contains(ss.str())))
0486 continue;
0487 if (itHst->first.Contains("nPVs"))
0488 continue;
0489 edm::LogInfo("BX|BeamMonitorBx") << "Filling histogram: " << itHst->first << endl;
0490 if (itHst->first.Contains("time")) {
0491 int idx = itHst->first.Index("_time", 5);
0492 itHst->second->setBinContent(ntbin_, val_[itHst->first(0, idx)].first);
0493 itHst->second->setBinError(ntbin_, val_[itHst->first(0, idx)].second);
0494 }
0495 if (itHst->first.Contains("lumi")) {
0496 int idx = itHst->first.Index("_lumi", 5);
0497 itHst->second->setBinContent(endLumiOfBSFit_, val_[itHst->first(0, idx)].first);
0498 itHst->second->setBinError(endLumiOfBSFit_, val_[itHst->first(0, idx)].second);
0499 }
0500 }
0501 TString histName = "Trending_nPVs_lumi_bx_" + ss.str();
0502 if (hst[histName])
0503 hst[histName]->setBinContent(endLumiOfBSFit_, nPVs_);
0504 }
0505
0506
0507 void BeamMonitorBx::weight(BeamSpotMapBx& weightedMap_, const BeamSpotMapBx& newMap_) {
0508 for (BeamSpotMapBx::const_iterator it = newMap_.begin(); it != newMap_.end(); ++it) {
0509 if (weightedMap_.find(it->first) == weightedMap_.end() || (it->second.type() != 2)) {
0510 weightedMap_[it->first] = it->second;
0511 continue;
0512 }
0513
0514 BeamSpot obs = weightedMap_[it->first];
0515 double val_[8] = {
0516 obs.x0(), obs.y0(), obs.z0(), obs.sigmaZ(), obs.dxdz(), obs.dydz(), obs.BeamWidthX(), obs.BeamWidthY()};
0517 double valErr_[8] = {obs.x0Error(),
0518 obs.y0Error(),
0519 obs.z0Error(),
0520 obs.sigmaZ0Error(),
0521 obs.dxdzError(),
0522 obs.dydzError(),
0523 obs.BeamWidthXError(),
0524 obs.BeamWidthYError()};
0525
0526 reco::BeamSpot::BeamType type = reco::BeamSpot::Unknown;
0527 weight(val_[0], valErr_[0], it->second.x0(), it->second.x0Error());
0528 weight(val_[1], valErr_[1], it->second.y0(), it->second.y0Error());
0529 weight(val_[2], valErr_[2], it->second.z0(), it->second.z0Error());
0530 weight(val_[3], valErr_[3], it->second.sigmaZ(), it->second.sigmaZ0Error());
0531 weight(val_[4], valErr_[4], it->second.dxdz(), it->second.dxdzError());
0532 weight(val_[5], valErr_[5], it->second.dydz(), it->second.dydzError());
0533 weight(val_[6], valErr_[6], it->second.BeamWidthX(), it->second.BeamWidthXError());
0534 weight(val_[7], valErr_[7], it->second.BeamWidthY(), it->second.BeamWidthYError());
0535 if (it->second.type() == reco::BeamSpot::Tracker) {
0536 type = reco::BeamSpot::Tracker;
0537 }
0538
0539 reco::BeamSpot::Point bsPosition(val_[0], val_[1], val_[2]);
0540 reco::BeamSpot::CovarianceMatrix error;
0541 for (int i = 0; i < 7; ++i) {
0542 error(i, i) = valErr_[i] * valErr_[i];
0543 }
0544 reco::BeamSpot weightedBeamSpot(bsPosition, val_[3], val_[4], val_[5], val_[6], error, type);
0545 weightedBeamSpot.setBeamWidthY(val_[7]);
0546 LogInfo("BX|BeamMonitorBx") << weightedBeamSpot << endl;
0547 weightedMap_.erase(it->first);
0548 weightedMap_[it->first] = weightedBeamSpot;
0549 }
0550 }
0551
0552
0553 void BeamMonitorBx::weight(double& mean, double& meanError, const double& val, const double& valError) {
0554 double tmpError = 0;
0555 if (meanError < 1e-8) {
0556 tmpError = 1 / (valError * valError);
0557 mean = val * tmpError;
0558 } else {
0559 tmpError = 1 / (meanError * meanError) + 1 / (valError * valError);
0560 mean = mean / (meanError * meanError) + val / (valError * valError);
0561 }
0562 mean = mean / tmpError;
0563 meanError = std::sqrt(1 / tmpError);
0564 }
0565
0566
0567 void BeamMonitorBx::endRun(const Run& r, const EventSetup& context) {}
0568
0569 DEFINE_FWK_MODULE(BeamMonitorBx);
