CMS 3D CMS Logo
BeamMonitorBx.cc
Go to the documentation of this file.
1 /*
2  * \file BeamMonitorBx.cc
3  * \author Geng-yuan Jeng/UC Riverside
4  * Francisco Yumiceva/FNAL
5  *
6  */
7 
8 #include "DQM/BeamMonitor/plugins/BeamMonitorBx.h"
9 #include "FWCore/ServiceRegistry/interface/Service.h"
10 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
11 #include "DataFormats/Common/interface/View.h"
12 #include "RecoVertex/BeamSpotProducer/interface/BSFitter.h"
13 #include "FWCore/Framework/interface/MakerMacros.h"
14 #include "FWCore/Framework/interface/Run.h"
15 #include "FWCore/Framework/interface/LuminosityBlock.h"
16 #include "FWCore/Framework/interface/ConsumesCollector.h"
17 #include "FWCore/MessageLogger/interface/MessageLogger.h"
18 #include <numeric>
19 #include <cmath>
20 #include <TMath.h>
21 #include <iostream>
22 #include <TStyle.h>
23 
24 using namespace std;
25 using namespace edm;
26 using namespace reco;
27 
28 void BeamMonitorBx::formatFitTime(char* ts, const time_t& t) {
29 #define CET (+1)
30 #define CEST (+2)
31 
32  tm* ptm;
33  ptm = gmtime(&t);
34  sprintf(ts,
35  "%4d-%02d-%02d %02d:%02d:%02d",
36  ptm->tm_year,
37  ptm->tm_mon + 1,
38  ptm->tm_mday,
39  (ptm->tm_hour + CEST) % 24,
40  ptm->tm_min,
41  ptm->tm_sec);
42 
43 #ifdef STRIP_TRAILING_BLANKS_IN_TIMEZONE
44  unsigned int b = strlen(ts);
45  while (ts[--b] == ' ') {
46  ts[b] = 0;
47  }
48 #endif
49 }
50 
51 //
52 // constructors and destructor
53 //
54 BeamMonitorBx::BeamMonitorBx(const ParameterSet& ps) : countBx_(0), countEvt_(0), countLumi_(0), resetHistos_(false) {
55  parameters_ = ps;
56  monitorName_ = parameters_.getUntrackedParameter<string>("monitorName", "YourSubsystemName");
57  bsSrc_ = parameters_.getUntrackedParameter<InputTag>("beamSpot");
58  fitNLumi_ = parameters_.getUntrackedParameter<int>("fitEveryNLumi", -1);
59  resetFitNLumi_ = parameters_.getUntrackedParameter<int>("resetEveryNLumi", -1);
60 
61  dbe_ = Service<DQMStore>().operator->();
62 
63  if (!monitorName_.empty())
64  monitorName_ = monitorName_ + "/";
65 
66  theBeamFitter = new BeamFitter(parameters_, consumesCollector());
67  theBeamFitter->resetTrkVector();
68  theBeamFitter->resetLSRange();
69  theBeamFitter->resetRefTime();
70  theBeamFitter->resetPVFitter();
71 
72  if (fitNLumi_ <= 0)
73  fitNLumi_ = 1;
74  beginLumiOfBSFit_ = endLumiOfBSFit_ = 0;
75  refBStime[0] = refBStime[1] = 0;
76  lastlumi_ = 0;
77  nextlumi_ = 0;
78  firstlumi_ = 0;
79  processed_ = false;
80  countGoodFit_ = 0;
81 }
82 
83 BeamMonitorBx::~BeamMonitorBx() { delete theBeamFitter; }
84 
85 //--------------------------------------------------------
86 void BeamMonitorBx::beginJob() {
87  varMap["x0_bx"] = "X_{0}";
88  varMap["y0_bx"] = "Y_{0}";
89  varMap["z0_bx"] = "Z_{0}";
90  varMap["sigmaX_bx"] = "#sigma_{X}";
91  varMap["sigmaY_bx"] = "#sigma_{Y}";
92  varMap["sigmaZ_bx"] = "#sigma_{Z}";
93 
94  varMap1["x"] = "X_{0} [cm]";
95  varMap1["y"] = "Y_{0} [cm]";
96  varMap1["z"] = "Z_{0} [cm]";
97  varMap1["sigmaX"] = "#sigma_{X} [cm]";
98  varMap1["sigmaY"] = "#sigma_{Y} [cm]";
99  varMap1["sigmaZ"] = "#sigma_{Z} [cm]";
100 
101  // create and cd into new folder
102  dbe_->setCurrentFolder(monitorName_ + "FitBx");
103  // Results of good fit:
104  BookTables(1, varMap, "");
105  //if (resetFitNLumi_ > 0) BookTables(1,varMap,"all");
106 
107  // create and cd into new folders
108  for (std::map<std::string, std::string>::const_iterator varName = varMap1.begin(); varName != varMap1.end();
109  ++varName) {
110  string subDir_ = "FitBx";
111  subDir_ += "/";
112  subDir_ += "All_";
113  subDir_ += varName->first;
114  dbe_->setCurrentFolder(monitorName_ + subDir_);
115  }
116  dbe_->setCurrentFolder(monitorName_ + "FitBx/All_nPVs");
117 }
118 
119 //--------------------------------------------------------
120 void BeamMonitorBx::beginRun(const edm::Run& r, const EventSetup& context) {
121  ftimestamp = r.beginTime().value();
122  tmpTime = ftimestamp >> 32;
123  startTime = refTime = tmpTime;
124 }
125 
126 //--------------------------------------------------------
127 void BeamMonitorBx::beginLuminosityBlock(const LuminosityBlock& lumiSeg, const EventSetup& context) {
128  int nthlumi = lumiSeg.luminosityBlock();
129  const edm::TimeValue_t fbegintimestamp = lumiSeg.beginTime().value();
130  const std::time_t ftmptime = fbegintimestamp >> 32;
131 
132  if (countLumi_ == 0) {
133  beginLumiOfBSFit_ = nthlumi;
134  refBStime[0] = ftmptime;
135  }
136  if (beginLumiOfBSFit_ == 0)
137  beginLumiOfBSFit_ = nextlumi_;
138 
139  if (nthlumi < nextlumi_)
140  return;
141 
142  if (nthlumi > nextlumi_) {
143  if (countLumi_ != 0 && processed_) {
144  FitAndFill(lumiSeg, lastlumi_, nextlumi_, nthlumi);
145  }
146  nextlumi_ = nthlumi;
147  edm::LogInfo("LS|BX|BeamMonitorBx") << "Next Lumi to Fit: " << nextlumi_ << endl;
148  if (refBStime[0] == 0)
149  refBStime[0] = ftmptime;
150  }
151  countLumi_++;
152  if (processed_)
153  processed_ = false;
154  edm::LogInfo("LS|BX|BeamMonitorBx") << "Begin of Lumi: " << nthlumi << endl;
155 }
156 
157 // ----------------------------------------------------------
158 void BeamMonitorBx::analyze(const Event& iEvent, const EventSetup& iSetup) {
159  bool readEvent_ = true;
160  const int nthlumi = iEvent.luminosityBlock();
161  if (nthlumi < nextlumi_) {
162  readEvent_ = false;
163  edm::LogWarning("BX|BeamMonitorBx") << "Spilt event from previous lumi section!" << endl;
164  }
165  if (nthlumi > nextlumi_) {
166  readEvent_ = false;
167  edm::LogWarning("BX|BeamMonitorBx") << "Spilt event from next lumi section!!!" << endl;
168  }
169 
170  if (readEvent_) {
171  countEvt_++;
172  theBeamFitter->readEvent(iEvent);
173  processed_ = true;
174  }
175 }
176 
177 //--------------------------------------------------------
178 void BeamMonitorBx::endLuminosityBlock(const LuminosityBlock& lumiSeg, const EventSetup& iSetup) {
179  int nthlumi = lumiSeg.id().luminosityBlock();
180  edm::LogInfo("LS|BX|BeamMonitorBx") << "Lumi of the last event before endLuminosityBlock: " << nthlumi << endl;
181 
182  if (nthlumi < nextlumi_)
183  return;
184  const edm::TimeValue_t fendtimestamp = lumiSeg.endTime().value();
185  const std::time_t fendtime = fendtimestamp >> 32;
186  tmpTime = refBStime[1] = fendtime;
187 }
188 
189 //--------------------------------------------------------
190 void BeamMonitorBx::BookTables(int nBx, map<string, string>& vMap, string suffix_) {
191  // to rebin histograms when number of bx increases
192  dbe_->cd(monitorName_ + "FitBx");
193 
194  for (std::map<std::string, std::string>::const_iterator varName = vMap.begin(); varName != vMap.end(); ++varName) {
195  string tmpName = varName->first;
196  if (!suffix_.empty()) {
197  tmpName += "_";
198  tmpName += suffix_;
199  }
200 
201  hs[tmpName] = dbe_->book2D(tmpName, varName->second, 3, 0, 3, nBx, 0, nBx);
202  hs[tmpName]->setBinLabel(1, "bx", 1);
203  hs[tmpName]->setBinLabel(2, varName->second, 1);
204  hs[tmpName]->setBinLabel(3, "#Delta " + varName->second, 1);
205  for (int i = 0; i < nBx; i++) {
206  hs[tmpName]->setBinLabel(i + 1, " ", 2);
207  }
208  hs[tmpName]->getTH1()->SetOption("text");
209  hs[tmpName]->Reset();
210  }
211 }
212 
213 //--------------------------------------------------------
214 void BeamMonitorBx::BookTrendHistos(
215  bool plotPV, int nBx, map<string, string>& vMap, string subDir_, const TString& prefix_, const TString& suffix_) {
216  int nType_ = 2;
217  if (plotPV)
218  nType_ = 4;
219  std::ostringstream ss;
220  std::ostringstream ss1;
221  ss << setfill('0') << setw(5) << nBx;
222  ss1 << nBx;
223 
224  for (int i = 0; i < nType_; i++) {
225  for (std::map<std::string, std::string>::const_iterator varName = vMap.begin(); varName != vMap.end(); ++varName) {
226  string tmpDir_ = subDir_ + "/All_" + varName->first;
227  dbe_->cd(monitorName_ + tmpDir_);
228  TString histTitle(varName->first);
229  TString tmpName;
230  if (prefix_ != "")
231  tmpName = prefix_ + "_" + varName->first;
232  if (suffix_ != "")
233  tmpName = tmpName + "_" + suffix_;
234  tmpName = tmpName + "_" + ss.str();
235 
236  TString histName(tmpName);
237  string ytitle(varName->second);
238  string xtitle("");
239  string options("E1");
240  bool createHisto = true;
241  switch (i) {
242  case 1: // BS vs time
243  histName.Insert(histName.Index("_bx_", 4), "_time");
244  xtitle = "Time [UTC] [Bx# " + ss1.str() + "]";
245  if (ytitle.find("sigma") == string::npos)
246  histTitle += " coordinate of beam spot vs time (Fit)";
247  else
248  histTitle = histTitle.Insert(5, " ") + " of beam spot vs time (Fit)";
249  break;
250  case 2: // PV +/- sigmaPV vs lumi
251  if (ytitle.find("sigma") == string::npos) {
252  histName.Insert(0, "PV");
253  histName.Insert(histName.Index("_bx_", 4), "_lumi");
254  histTitle.Insert(0, "Avg. ");
255  histTitle += " position of primary vtx vs lumi";
256  xtitle = "Lumisection [Bx# " + ss1.str() + "]";
257  ytitle.insert(0, "PV");
258  ytitle += " #pm #sigma_{PV";
259  ytitle += varName->first;
260  ytitle += "} (cm)";
261  } else
262  createHisto = false;
263  break;
264  case 3: // PV +/- sigmaPV vs time
265  if (ytitle.find("sigma") == string::npos) {
266  histName.Insert(0, "PV");
267  histName.Insert(histName.Index("_bx_", 4), "_time");
268  histTitle.Insert(0, "Avg. ");
269  histTitle += " position of primary vtx vs time";
270  xtitle = "Time [UTC] [Bx# " + ss1.str() + "]";
271  ytitle.insert(0, "PV");
272  ytitle += " #pm #sigma_{PV";
273  ytitle += varName->first;
274  ytitle += "} (cm)";
275  } else
276  createHisto = false;
277  break;
278  default: // BS vs lumi
279  histName.Insert(histName.Index("_bx_", 4), "_lumi");
280  xtitle = "Lumisection [Bx# " + ss1.str() + "]";
281  if (ytitle.find("sigma") == string::npos)
282  histTitle += " coordinate of beam spot vs lumi (Fit)";
283  else
284  histTitle = histTitle.Insert(5, " ") + " of beam spot vs lumi (Fit)";
285  break;
286  }
287  // check if already exist
288  if (dbe_->get(monitorName_ + tmpDir_ + "/" + string(histName.Data())))
289  continue;
290 
291  if (createHisto) {
292  edm::LogInfo("BX|BeamMonitorBx") << "histName = " << histName << "; histTitle = " << histTitle << std::endl;
293  hst[histName] = dbe_->book1D(histName, histTitle, 40, 0.5, 40.5);
294  hst[histName]->getTH1()->SetCanExtend(TH1::kAllAxes);
295  hst[histName]->setAxisTitle(xtitle, 1);
296  hst[histName]->setAxisTitle(ytitle, 2);
297  hst[histName]->getTH1()->SetOption("E1");
298  if (histName.Contains("time")) {
299  hst[histName]->getTH1()->SetBins(3600, 0.5, 3600 + 0.5);
300  hst[histName]->setAxisTimeDisplay(1);
301  hst[histName]->setAxisTimeFormat("%H:%M:%S", 1);
302 
303  char eventTime[64];
304  formatFitTime(eventTime, startTime);
305  TDatime da(eventTime);
306  if (debug_) {
307  edm::LogInfo("BX|BeamMonitorBx") << "TimeOffset = ";
308  da.Print();
309  }
310  hst[histName]->getTH1()->GetXaxis()->SetTimeOffset(da.Convert(kTRUE));
311  }
312  }
313  } //End of variable loop
314  } // End of type loop (lumi, time)
315 
316  // num of PVs(#Bx) per LS
317  dbe_->cd(monitorName_ + subDir_ + "/All_nPVs");
318  TString histName = "Trending_nPVs_lumi_bx_" + ss.str();
319  string xtitle = "Lumisection [Bx# " + ss1.str() + "]";
320 
321  hst[histName] = dbe_->book1D(histName, "Number of Good Reconstructed Vertices", 40, 0.5, 40.5);
322  hst[histName]->setAxisTitle(xtitle, 1);
323  hst[histName]->getTH1()->SetCanExtend(TH1::kAllAxes);
324  hst[histName]->getTH1()->SetOption("E1");
325 }
326 
327 //--------------------------------------------------------
328 void BeamMonitorBx::FitAndFill(const LuminosityBlock& lumiSeg, int& lastlumi, int& nextlumi, int& nthlumi) {
329  if (nthlumi <= nextlumi)
330  return;
331 
332  int currentlumi = nextlumi;
333  edm::LogInfo("LS|BX|BeamMonitorBx") << "Lumi of the current fit: " << currentlumi << endl;
334  lastlumi = currentlumi;
335  endLumiOfBSFit_ = currentlumi;
336 
337  edm::LogInfo("BX|BeamMonitorBx") << "Time lapsed = " << tmpTime - refTime << std::endl;
338 
339  if (resetHistos_) {
340  edm::LogInfo("BX|BeamMonitorBx") << "Resetting Histograms" << endl;
341  theBeamFitter->resetCutFlow();
342  resetHistos_ = false;
343  }
344 
345  if (fitNLumi_ > 0)
346  if (currentlumi % fitNLumi_ != 0)
347  return;
348 
349  std::pair<int, int> fitLS = theBeamFitter->getFitLSRange();
350  edm::LogInfo("LS|BX|BeamMonitorBx") << " [Fitter] Do BeamSpot Fit for LS = " << fitLS.first << " to " << fitLS.second
351  << endl;
352  edm::LogInfo("LS|BX|BeamMonitorBx") << " [BX] Do BeamSpot Fit for LS = " << beginLumiOfBSFit_ << " to "
353  << endLumiOfBSFit_ << endl;
354 
355  edm::LogInfo("BX|BeamMonitorBx") << " [BxDebugTime] refBStime[0] = " << refBStime[0]
356  << "; address = " << &refBStime[0] << std::endl;
357  edm::LogInfo("BX|BeamMonitorBx") << " [BxDebugTime] refBStime[1] = " << refBStime[1]
358  << "; address = " << &refBStime[1] << std::endl;
359 
360  if (theBeamFitter->runPVandTrkFitter()) {
361  countGoodFit_++;
362  edm::LogInfo("BX|BeamMonitorBx") << "Number of good fit = " << countGoodFit_ << endl;
363  BeamSpotMapBx bsmap = theBeamFitter->getBeamSpotMap();
364  std::map<int, int> npvsmap = theBeamFitter->getNPVsperBX();
365  edm::LogInfo("BX|BeamMonitorBx") << "Number of bx = " << bsmap.size() << endl;
366  if (bsmap.empty())
367  return;
368  if (countBx_ < bsmap.size()) {
369  countBx_ = bsmap.size();
370  BookTables(countBx_, varMap, "");
371  BookTables(countBx_, varMap, "all");
372  for (BeamSpotMapBx::const_iterator abspot = bsmap.begin(); abspot != bsmap.end(); ++abspot) {
373  int bx = abspot->first;
374  BookTrendHistos(false, bx, varMap1, "FitBx", "Trending", "bx");
375  }
376  }
377 
378  std::pair<int, int> LSRange = theBeamFitter->getFitLSRange();
379  char tmpTitle[50];
380  sprintf(tmpTitle, "%s %i %s %i %s", " [cm] (LS: ", LSRange.first, " to ", LSRange.second, ")");
381  for (std::map<std::string, std::string>::const_iterator varName = varMap.begin(); varName != varMap.end();
382  ++varName) {
383  hs[varName->first]->setTitle(varName->second + " " + tmpTitle);
384  hs[varName->first]->Reset();
385  }
386 
387  if (countGoodFit_ == 1)
388  firstlumi_ = LSRange.first;
389 
390  if (resetFitNLumi_ > 0) {
391  char tmpTitle1[60];
392  if (countGoodFit_ > 1)
393  snprintf(tmpTitle1,
394  sizeof(tmpTitle1),
395  "%s %i %s %i %s",
396  " [cm] (LS: ",
397  firstlumi_,
398  " to ",
399  LSRange.second,
400  ") [weighted average]");
401  else
402  snprintf(tmpTitle1, sizeof(tmpTitle1), "%s", "Need at least two fits to calculate weighted average");
403  for (std::map<std::string, std::string>::const_iterator varName = varMap.begin(); varName != varMap.end();
404  ++varName) {
405  TString tmpName = varName->first + "_all";
406  hs[tmpName]->setTitle(varName->second + " " + tmpTitle1);
407  hs[tmpName]->Reset();
408  }
409  }
410 
411  int nthBin = countBx_;
412  for (BeamSpotMapBx::const_iterator abspot = bsmap.begin(); abspot != bsmap.end(); ++abspot, nthBin--) {
413  reco::BeamSpot bs = abspot->second;
414  int bx = abspot->first;
415  int nPVs = npvsmap.find(bx)->second;
416  edm::LogInfo("BeamMonitorBx") << "Number of PVs of bx#[" << bx << "] = " << nPVs << endl;
417  // Tables
418  FillTables(bx, nthBin, varMap, bs, "");
419  // Histograms
420  FillTrendHistos(bx, nPVs, varMap1, bs, "Trending");
421  }
422  // Calculate weighted beam spots
423  weight(fbspotMap, bsmap);
424  // Fill the results
425  nthBin = countBx_;
426  if (resetFitNLumi_ > 0 && countGoodFit_ > 1) {
427  for (BeamSpotMapBx::const_iterator abspot = fbspotMap.begin(); abspot != fbspotMap.end(); ++abspot, nthBin--) {
428  reco::BeamSpot bs = abspot->second;
429  int bx = abspot->first;
430  FillTables(bx, nthBin, varMap, bs, "all");
431  }
432  }
433  }
434  // else
435  // edm::LogInfo("BeamMonitorBx") << "Bad Fit!!!" << endl;
436 
437  if (resetFitNLumi_ > 0 && currentlumi % resetFitNLumi_ == 0) {
438  edm::LogInfo("LS|BX|BeamMonitorBx") << "Reset track collection for beam fit!!!" << endl;
439  resetHistos_ = true;
440  theBeamFitter->resetTrkVector();
441  theBeamFitter->resetLSRange();
442  theBeamFitter->resetRefTime();
443  theBeamFitter->resetPVFitter();
444  beginLumiOfBSFit_ = 0;
445  refBStime[0] = 0;
446  }
447 }
448 
449 //--------------------------------------------------------
450 void BeamMonitorBx::FillTables(int nthbx, int nthbin_, map<string, string>& vMap, reco::BeamSpot& bs_, string suffix_) {
451  map<string, pair<double, double> > val_;
452  val_["x0_bx"] = pair<double, double>(bs_.x0(), bs_.x0Error());
453  val_["y0_bx"] = pair<double, double>(bs_.y0(), bs_.y0Error());
454  val_["z0_bx"] = pair<double, double>(bs_.z0(), bs_.z0Error());
455  val_["sigmaX_bx"] = pair<double, double>(bs_.BeamWidthX(), bs_.BeamWidthXError());
456  val_["sigmaY_bx"] = pair<double, double>(bs_.BeamWidthY(), bs_.BeamWidthYError());
457  val_["sigmaZ_bx"] = pair<double, double>(bs_.sigmaZ(), bs_.sigmaZ0Error());
458 
459  for (std::map<std::string, std::string>::const_iterator varName = vMap.begin(); varName != vMap.end(); ++varName) {
460  TString tmpName = varName->first;
461  if (!suffix_.empty())
462  tmpName += TString("_" + suffix_);
463  hs[tmpName]->setBinContent(1, nthbin_, nthbx);
464  hs[tmpName]->setBinContent(2, nthbin_, val_[varName->first].first);
465  hs[tmpName]->setBinContent(3, nthbin_, val_[varName->first].second);
466  }
467 }
468 
469 //--------------------------------------------------------
470 void BeamMonitorBx::FillTrendHistos(
471  int nthBx, int nPVs_, map<string, string>& vMap, reco::BeamSpot& bs_, const TString& prefix_) {
472  map<TString, pair<double, double> > val_;
473  val_[TString(prefix_ + "_x")] = pair<double, double>(bs_.x0(), bs_.x0Error());
474  val_[TString(prefix_ + "_y")] = pair<double, double>(bs_.y0(), bs_.y0Error());
475  val_[TString(prefix_ + "_z")] = pair<double, double>(bs_.z0(), bs_.z0Error());
476  val_[TString(prefix_ + "_sigmaX")] = pair<double, double>(bs_.BeamWidthX(), bs_.BeamWidthXError());
477  val_[TString(prefix_ + "_sigmaY")] = pair<double, double>(bs_.BeamWidthY(), bs_.BeamWidthYError());
478  val_[TString(prefix_ + "_sigmaZ")] = pair<double, double>(bs_.sigmaZ(), bs_.sigmaZ0Error());
479 
480  std::ostringstream ss;
481  ss << setfill('0') << setw(5) << nthBx;
482  int ntbin_ = tmpTime - startTime;
483  for (map<TString, MonitorElement*>::iterator itHst = hst.begin(); itHst != hst.end(); ++itHst) {
484  if (!(itHst->first.Contains(ss.str())))
485  continue;
486  if (itHst->first.Contains("nPVs"))
487  continue;
488  edm::LogInfo("BX|BeamMonitorBx") << "Filling histogram: " << itHst->first << endl;
489  if (itHst->first.Contains("time")) {
490  int idx = itHst->first.Index("_time", 5);
491  itHst->second->setBinContent(ntbin_, val_[itHst->first(0, idx)].first);
492  itHst->second->setBinError(ntbin_, val_[itHst->first(0, idx)].second);
493  }
494  if (itHst->first.Contains("lumi")) {
495  int idx = itHst->first.Index("_lumi", 5);
496  itHst->second->setBinContent(endLumiOfBSFit_, val_[itHst->first(0, idx)].first);
497  itHst->second->setBinError(endLumiOfBSFit_, val_[itHst->first(0, idx)].second);
498  }
499  }
500  TString histName = "Trending_nPVs_lumi_bx_" + ss.str();
501  if (hst[histName])
502  hst[histName]->setBinContent(endLumiOfBSFit_, nPVs_);
503 }
504 
505 //--------------------------------------------------------------------------------------------------
506 void BeamMonitorBx::weight(BeamSpotMapBx& weightedMap_, const BeamSpotMapBx& newMap_) {
507  for (BeamSpotMapBx::const_iterator it = newMap_.begin(); it != newMap_.end(); ++it) {
508  if (weightedMap_.find(it->first) == weightedMap_.end() || (it->second.type() != 2)) {
509  weightedMap_[it->first] = it->second;
510  continue;
511  }
512 
513  BeamSpot obs = weightedMap_[it->first];
514  double val_[8] = {
515  obs.x0(), obs.y0(), obs.z0(), obs.sigmaZ(), obs.dxdz(), obs.dydz(), obs.BeamWidthX(), obs.BeamWidthY()};
516  double valErr_[8] = {obs.x0Error(),
517  obs.y0Error(),
518  obs.z0Error(),
519  obs.sigmaZ0Error(),
520  obs.dxdzError(),
521  obs.dydzError(),
522  obs.BeamWidthXError(),
523  obs.BeamWidthYError()};
524 
525  reco::BeamSpot::BeamType type = reco::BeamSpot::Unknown;
526  weight(val_[0], valErr_[0], it->second.x0(), it->second.x0Error());
527  weight(val_[1], valErr_[1], it->second.y0(), it->second.y0Error());
528  weight(val_[2], valErr_[2], it->second.z0(), it->second.z0Error());
529  weight(val_[3], valErr_[3], it->second.sigmaZ(), it->second.sigmaZ0Error());
530  weight(val_[4], valErr_[4], it->second.dxdz(), it->second.dxdzError());
531  weight(val_[5], valErr_[5], it->second.dydz(), it->second.dydzError());
532  weight(val_[6], valErr_[6], it->second.BeamWidthX(), it->second.BeamWidthXError());
533  weight(val_[7], valErr_[7], it->second.BeamWidthY(), it->second.BeamWidthYError());
534  if (it->second.type() == reco::BeamSpot::Tracker) {
535  type = reco::BeamSpot::Tracker;
536  }
537 
538  reco::BeamSpot::Point bsPosition(val_[0], val_[1], val_[2]);
539  reco::BeamSpot::CovarianceMatrix error;
540  for (int i = 0; i < 7; ++i) {
541  error(i, i) = valErr_[i] * valErr_[i];
542  }
543  reco::BeamSpot weightedBeamSpot(bsPosition, val_[3], val_[4], val_[5], val_[6], error, type);
544  weightedBeamSpot.setBeamWidthY(val_[7]);
545  LogInfo("BX|BeamMonitorBx") << weightedBeamSpot << endl;
546  weightedMap_.erase(it->first);
547  weightedMap_[it->first] = weightedBeamSpot;
548  }
549 }
550 
551 //--------------------------------------------------------------------------------------------------
552 void BeamMonitorBx::weight(double& mean, double& meanError, const double& val, const double& valError) {
553  double tmpError = 0;
554  if (meanError < 1e-8) {
555  tmpError = 1 / (valError * valError);
556  mean = val * tmpError;
557  } else {
558  tmpError = 1 / (meanError * meanError) + 1 / (valError * valError);
559  mean = mean / (meanError * meanError) + val / (valError * valError);
560  }
561  mean = mean / tmpError;
562  meanError = std::sqrt(1 / tmpError);
563 }
564 
565 //--------------------------------------------------------
566 void BeamMonitorBx::endRun(const Run& r, const EventSetup& context) {}
567 
568 DEFINE_FWK_MODULE(BeamMonitorBx);
BeamMonitorBx::endLuminosityBlock
void endLuminosityBlock(const edm::LuminosityBlock &lumiSeg, const edm::EventSetup &c) override
Definition: BeamMonitorBx.cc:178
BeamMonitorBx::processed_
bool processed_
Definition: BeamMonitorBx.h:84
BeamSpotPI::nPVs
Definition: BeamSpotPayloadInspectorHelper.h:43
BeamFitter::runPVandTrkFitter
bool runPVandTrkFitter()
Definition: BeamFitter.cc:449
BeamFitter::resetRefTime
void resetRefTime()
Definition: BeamFitter.h:61
mps_fire.i
i
Definition: mps_fire.py:428
BeamMonitorBx::BeamMonitorBx
BeamMonitorBx(const edm::ParameterSet &)
Definition: BeamMonitorBx.cc:54
BeamMonitorBx::FitAndFill
void FitAndFill(const edm::LuminosityBlock &lumiSeg, int &, int &, int &)
Definition: BeamMonitorBx.cc:328
BeamFitter::resetPVFitter
void resetPVFitter()
Definition: BeamFitter.h:72
BeamMonitorBx::countGoodFit_
int countGoodFit_
Definition: BeamMonitorBx.h:80
reco::BeamSpot::setBeamWidthY
void setBeamWidthY(double v)
Definition: BeamSpot.h:105
BeamFitter::resetTrkVector
void resetTrkVector()
Definition: BeamFitter.h:58
MessageLogger.h
SiStripPI::mean
Definition: SiStripPayloadInspectorHelper.h:169
reco::BeamSpot::BeamWidthXError
double BeamWidthXError() const
error on beam width X, assume error in X = Y
Definition: BeamSpot.h:99
funct::false
false
Definition: Factorize.h:29
AlcaBeamMonitor_cfi.BeamFitter
BeamFitter
Definition: AlcaBeamMonitor_cfi.py:10
BeamMonitorBx::~BeamMonitorBx
~BeamMonitorBx() override
Definition: BeamMonitorBx.cc:83
BeamMonitorBx::startTime
std::time_t startTime
Definition: BeamMonitorBx.h:99
edm::LuminosityBlock
Definition: LuminosityBlock.h:50
edm::Run
Definition: Run.h:45
reco::BeamSpot::z0
double z0() const
z coordinate
Definition: BeamSpot.h:65
LuminosityBlock.h
edm
HLT enums.
Definition: AlignableModifier.h:19
BeamMonitorBx::lastlumi_
int lastlumi_
Definition: BeamMonitorBx.h:77
edm::LuminosityBlockBase::beginTime
Timestamp const & beginTime() const
Definition: LuminosityBlockBase.h:46
BeamMonitorBx::fbspotMap
BeamSpotMapBx fbspotMap
Definition: BeamMonitorBx.h:86
CEST
#define CEST
l1GtPatternGenerator_cfi.bx
bx
Definition: l1GtPatternGenerator_cfi.py:18
BeamMonitorBx::endRun
void endRun(const edm::Run &r, const edm::EventSetup &c) override
Definition: BeamMonitorBx.cc:566
Calorimetry_cff.da
da
Definition: Calorimetry_cff.py:155
BeamFitter::readEvent
void readEvent(const edm::Event &iEvent)
Definition: BeamFitter.cc:261
edm::EDConsumerBase::consumesCollector
ConsumesCollector consumesCollector()
Use a ConsumesCollector to gather consumes information from helper functions.
Definition: EDConsumerBase.cc:47
BeamMonitorBx::beginLuminosityBlock
void beginLuminosityBlock(const edm::LuminosityBlock &lumiSeg, const edm::EventSetup &context) override
Definition: BeamMonitorBx.cc:127
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
edm::Timestamp::value
TimeValue_t value() const
Definition: Timestamp.h:45
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:125
reco::BeamSpot::dxdz
double dxdz() const
dxdz slope
Definition: BeamSpot.h:78
reco::BeamSpot::dxdzError
double dxdzError() const
error on dxdz
Definition: BeamSpot.h:94
relativeConstraints.error
error
Definition: relativeConstraints.py:53
BeamMonitorBx::bsSrc_
edm::InputTag bsSrc_
Definition: BeamMonitorBx.h:63
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
BeamMonitorBx::beginRun
void beginRun(const edm::Run &r, const edm::EventSetup &c) override
Definition: BeamMonitorBx.cc:120
heavyIonCSV_trainingSettings.idx
idx
Definition: heavyIonCSV_trainingSettings.py:5
reco::BeamSpot::y0Error
double y0Error() const
error on y
Definition: BeamSpot.h:88
BeamMonitorBx::monitorName_
std::string monitorName_
Definition: BeamMonitorBx.h:62
contentValuesCheck.ss
ss
Definition: contentValuesCheck.py:33
BeamMonitorBx::ftimestamp
edm::TimeValue_t ftimestamp
Definition: BeamMonitorBx.h:100
reco::BeamSpot::sigmaZ
double sigmaZ() const
sigma z
Definition: BeamSpot.h:76
BeamMonitorBx::FillTrendHistos
void FillTrendHistos(int, int, std::map< std::string, std::string > &, reco::BeamSpot &, const TString &)
Definition: BeamMonitorBx.cc:470
MakerMacros.h
cms::cuda::bs
bs
Definition: HistoContainer.h:127
BeamMonitorBx
Definition: BeamMonitorBx.h:25
edm::LuminosityBlockBase::endTime
Timestamp const & endTime() const
Definition: LuminosityBlockBase.h:47
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
BeamSpot.h
reco::BeamSpot::dydzError
double dydzError() const
error on dydz
Definition: BeamSpot.h:96
Service.h
BeamMonitorBx::hs
std::map< TString, MonitorElement * > hs
Definition: BeamMonitorBx.h:90
edm::LuminosityBlockBase::luminosityBlock
LuminosityBlockNumber_t luminosityBlock() const
Definition: LuminosityBlockBase.h:40
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
reco::BeamSpot
Definition: BeamSpot.h:21
reco::BeamSpot::Tracker
Definition: BeamSpot.h:24
Run.h
BeamMonitorBx::refBStime
std::time_t refBStime[2]
Definition: BeamMonitorBx.h:81
hgcalPlots.xtitle
xtitle
Definition: hgcalPlots.py:94
BeamMonitorBx::countBx_
unsigned int countBx_
Definition: BeamMonitorBx.h:72
BeamMonitorBx::theBeamFitter
BeamFitter * theBeamFitter
Definition: BeamMonitorBx.h:70
dqm::implementation::DQMStore::setCurrentFolder
void setCurrentFolder(std::string const &fullpath) override
Definition: DQMStore.h:569
b
double b
Definition: hdecay.h:118
BeamMonitorBx::endLumiOfBSFit_
int endLumiOfBSFit_
Definition: BeamMonitorBx.h:76
reco::BeamSpot::Unknown
Definition: BeamSpot.h:24
BeamFitter::getBeamSpotMap
std::map< int, reco::BeamSpot > getBeamSpotMap()
Definition: BeamFitter.h:94
BeamMonitorBx::dbe_
DQMStore * dbe_
Definition: BeamMonitorBx.h:69
BeamMonitorBx::firstlumi_
int firstlumi_
Definition: BeamMonitorBx.h:79
edm::ParameterSet
Definition: ParameterSet.h:47
BeamMonitorBx::BeamSpotMapBx
std::map< BxNum, reco::BeamSpot > BeamSpotMapBx
Definition: BeamMonitorBx.h:34
B2GTnPMonitor_cfi.histTitle
histTitle
Definition: B2GTnPMonitor_cfi.py:112
BeamMonitorBx::beginLumiOfBSFit_
int beginLumiOfBSFit_
Definition: BeamMonitorBx.h:75
type
type
Definition: SiPixelVCal_PayloadInspector.cc:37
BeamMonitorBx::countEvt_
int countEvt_
Definition: BeamMonitorBx.h:73
BeamMonitorBx::analyze
void analyze(const edm::Event &e, const edm::EventSetup &c) override
Definition: BeamMonitorBx.cc:158
BeamMonitorBx::refTime
std::time_t refTime
Definition: BeamMonitorBx.h:98
BeamMonitorBx::countLumi_
int countLumi_
Definition: BeamMonitorBx.h:74
edm::Service
Definition: Service.h:30
iEvent
int iEvent
Definition: GenABIO.cc:224
edm::LuminosityBlockID::luminosityBlock
LuminosityBlockNumber_t luminosityBlock() const
Definition: LuminosityBlockID.h:42
edm::LuminosityBlockBase::id
LuminosityBlockID id() const
Definition: LuminosityBlockBase.h:44
BeamFitter::resetLSRange
void resetLSRange()
Definition: BeamFitter.h:60
BeamMonitorBx::resetFitNLumi_
int resetFitNLumi_
Definition: BeamMonitorBx.h:66
plotBeamSpotDB.ptm
ptm
Definition: plotBeamSpotDB.py:396
reco::BeamSpot::z0Error
double z0Error() const
error on z
Definition: BeamSpot.h:90
reco::BeamSpot::BeamWidthYError
double BeamWidthYError() const
error on beam width Y, assume error in X = Y
Definition: BeamSpot.h:101
edm::EventSetup
Definition: EventSetup.h:58
reco::BeamSpot::CovarianceMatrix
math::Error< dimension >::type CovarianceMatrix
Definition: BeamSpot.h:29
reco::BeamSpot::BeamWidthX
double BeamWidthX() const
beam width X
Definition: BeamSpot.h:82
reco::BeamSpot::BeamType
BeamType
beam spot flags
Definition: BeamSpot.h:24
BeamMonitorBx::parameters_
edm::ParameterSet parameters_
Definition: BeamMonitorBx.h:61
BeamFitter::resetCutFlow
void resetCutFlow()
Definition: BeamFitter.h:104
BeamMonitorBx::beginJob
void beginJob() override
Definition: BeamMonitorBx.cc:86
alignCSCRings.r
r
Definition: alignCSCRings.py:93
reco::BeamSpot::dydz
double dydz() const
dydz slope
Definition: BeamSpot.h:80
BeamMonitorBx::tmpTime
std::time_t tmpTime
Definition: BeamMonitorBx.h:97
BeamMonitorBx::debug_
bool debug_
Definition: BeamMonitorBx.h:67
reco::BeamSpot::x0
double x0() const
x coordinate
Definition: BeamSpot.h:61
BeamMonitorBx::BookTables
void BookTables(int, std::map< std::string, std::string > &, std::string)
Definition: BeamMonitorBx.cc:190
BeamMonitorBx::nextlumi_
int nextlumi_
Definition: BeamMonitorBx.h:78
heppy_batch.val
val
Definition: heppy_batch.py:351
std
Definition: JetResolutionObject.h:76
BeamMonitorBx::formatFitTime
void formatFitTime(char *, const std::time_t &)
Definition: BeamMonitorBx.cc:28
BeamMonitorBx::fitNLumi_
int fitNLumi_
Definition: BeamMonitorBx.h:65
BeamFitter::getFitLSRange
std::pair< int, int > getFitLSRange()
Definition: BeamFitter.h:114
reco::BeamSpot::x0Error
double x0Error() const
error on x
Definition: BeamSpot.h:86
dqm::implementation::IBooker::book2D
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, FUNC onbooking=NOOP())
Definition: DQMStore.h:177
reco::BeamSpot::sigmaZ0Error
double sigmaZ0Error() const
error on sigma z
Definition: BeamSpot.h:92
BeamMonitorBx::varMap
std::map< std::string, std::string > varMap
Definition: BeamMonitorBx.h:87
BeamMonitorBx::weight
void weight(BeamSpotMapBx &, const BeamSpotMapBx &)
Definition: BeamMonitorBx.cc:506
BeamMonitorBx::BookTrendHistos
void BookTrendHistos(bool, int, std::map< std::string, std::string > &, std::string, const TString &, const TString &)
Definition: BeamMonitorBx.cc:214
dqm::implementation::IGetter::get
virtual MonitorElement * get(std::string const &fullpath) const
Definition: DQMStore.cc:651
AlcaSiPixelAliHarvester0T_cff.options
options
Definition: AlcaSiPixelAliHarvester0T_cff.py:42
BeamMonitorBx::varMap1
std::map< std::string, std::string > varMap1
Definition: BeamMonitorBx.h:88
ConsumesCollector.h
reco::BeamSpot::Point
math::XYZPoint Point
point in the space
Definition: BeamSpot.h:27
View.h
BeamMonitorBx.h
BSFitter.h
reco::BeamSpot::y0
double y0() const
y coordinate
Definition: BeamSpot.h:63
BeamMonitorBx::FillTables
void FillTables(int, int, std::map< std::string, std::string > &, reco::BeamSpot &, std::string)
Definition: BeamMonitorBx.cc:450
BeamFitter::getNPVsperBX
const std::map< int, int > & getNPVsperBX()
Definition: BeamFitter.h:123
edm::Event
Definition: Event.h:73
submitPVValidationJobs.t
string t
Definition: submitPVValidationJobs.py:644
BeamMonitorBx::hst
std::map< TString, MonitorElement * > hst
Definition: BeamMonitorBx.h:91
BeamMonitorBx::resetHistos_
bool resetHistos_
Definition: BeamMonitorBx.h:83
edm::InputTag
Definition: InputTag.h:15
reco::BeamSpot::BeamWidthY
double BeamWidthY() const
beam width Y
Definition: BeamSpot.h:84
edm::TimeValue_t
unsigned long long TimeValue_t
Definition: Timestamp.h:28
dqm::implementation::DQMStore::cd
void cd() override
Definition: DQMStore.h:564
dqm::implementation::IBooker::book1D
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98
hgcalPlots.ytitle
ytitle
Definition: hgcalPlots.py:1110
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16