Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-11-02 23:48:42

0001 /*
0002  * \file BeamMonitorBx.cc
0003  * \author Geng-yuan Jeng/UC Riverside
0004  *         Francisco Yumiceva/FNAL
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 // constructors and destructor
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   // create and cd into new folder
0103   dbe_->setCurrentFolder(monitorName_ + "FitBx");
0104   // Results of good fit:
0105   BookTables(1, varMap, "");
0106   //if (resetFitNLumi_ > 0) BookTables(1,varMap,"all");
0107 
0108   // create and cd into new folders
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   // to rebin histograms when number of bx increases
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:  // BS vs time
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:  // PV +/- sigmaPV vs lumi
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:  // PV +/- sigmaPV vs time
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:  // BS vs lumi
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       // check if already exist
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     }  //End of variable loop
0315   }    // End of type loop (lumi, time)
0316 
0317   // num of PVs(#Bx) per LS
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       // Tables
0419       FillTables(bx, nthBin, varMap, bs, "");
0420       // Histograms
0421       FillTrendHistos(bx, nPVs, varMap1, bs, "Trending");
0422     }
0423     // Calculate weighted beam spots
0424     weight(fbspotMap, bsmap);
0425     // Fill the results
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   //   else
0436   //     edm::LogInfo("BeamMonitorBx") << "Bad Fit!!!" << endl;
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);