Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:09:21

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