CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
FakeBeamMonitor.cc
Go to the documentation of this file.
1 /*
2  * \file FakeBeamMonitor.cc
3  * \author Geng-yuan Jeng/UC Riverside
4  * Francisco Yumiceva/FNAL
5  */
6 
7 /*
8 The code has been modified for running average
9 mode, and it gives results for the last NLS which is
10 configurable.
11 Sushil S. Chauhan /UCDavis
12 Evan Friis /UCDavis
13 The last tag for working versions without this change is
14 V00-03-25
15 */
16 
32 #include <numeric>
33 #include <cmath>
34 #include <memory>
35 #include <TMath.h>
36 #include <TDatime.h>
37 #include <iostream>
38 #include <TStyle.h>
39 #include <ctime>
40 
41 using namespace std;
42 using namespace edm;
43 
44 void FakeBeamMonitor::formatFitTime(char* ts, const time_t& t) {
45  //constexpr int CET(+1);
46  constexpr int CEST(+2);
47 
48  //tm * ptm;
49  //ptm = gmtime ( &t );
50  //int year = ptm->tm_year;
51 
52  //get correct year from ctime
53  time_t currentTime;
54  struct tm* localTime;
55  time(&currentTime); // Get the current time
56  localTime = localtime(&currentTime); // Convert the current time to the local time
57  int year = localTime->tm_year + 1900;
58 
59  tm* ptm;
60  ptm = gmtime(&t);
61 
62  //check if year is ok
63  if (year <= 37)
64  year += 2000;
65  if (year >= 70 && year <= 137)
66  year += 1900;
67 
68  if (year < 1995) {
69  edm::LogError("BadTimeStamp") << "year reported is " << year << " !!" << std::endl;
70  //year = 2015; //overwritten later by BeamFitter.cc for fits but needed here for TH1
71  //edm::LogError("BadTimeStamp") << "Resetting to " <<year<<std::endl;
72  }
73  sprintf(ts,
74  "%4d-%02d-%02d %02d:%02d:%02d",
75  year,
76  ptm->tm_mon + 1,
77  ptm->tm_mday,
78  (ptm->tm_hour + CEST) % 24,
79  ptm->tm_min,
80  ptm->tm_sec);
81 
82 #ifdef STRIP_TRAILING_BLANKS_IN_TIMEZONE
83  unsigned int b = strlen(ts);
84  while (ts[--b] == ' ') {
85  ts[b] = 0;
86  }
87 #endif
88 }
89 
90 static constexpr int buffTime = 23;
91 
92 std::string FakeBeamMonitor::getGMTstring(const time_t& timeToConvert) {
93  char buff[32];
94  std::strftime(buff, sizeof(buff), "%Y.%m.%d %H:%M:%S GMT", gmtime(&timeToConvert));
95  std::string timeStr(buff);
96  return timeStr;
97 }
98 
99 //
100 // constructors and destructor
101 //
103 
105  : dxBin_(ps.getParameter<int>("dxBin")),
106  dxMin_(ps.getParameter<double>("dxMin")),
107  dxMax_(ps.getParameter<double>("dxMax")),
108 
109  vxBin_(ps.getParameter<int>("vxBin")),
110  vxMin_(ps.getParameter<double>("vxMin")),
111  vxMax_(ps.getParameter<double>("vxMax")),
112 
113  phiBin_(ps.getParameter<int>("phiBin")),
114  phiMin_(ps.getParameter<double>("phiMin")),
115  phiMax_(ps.getParameter<double>("phiMax")),
116 
117  dzBin_(ps.getParameter<int>("dzBin")),
118  dzMin_(ps.getParameter<double>("dzMin")),
119  dzMax_(ps.getParameter<double>("dzMax")),
120 
121  countEvt_(0),
122  countLumi_(0),
123  nthBSTrk_(0),
124  nFitElements_(3),
125  resetHistos_(false),
126  StartAverage_(false),
127  firstAverageFit_(0),
128  countGapLumi_(0) {
129  monitorName_ = ps.getUntrackedParameter<string>("monitorName", "YourSubsystemName");
130  recordName_ = ps.getUntrackedParameter<string>("recordName");
131  intervalInSec_ = ps.getUntrackedParameter<int>("timeInterval", 920); //40 LS X 23"
132  fitNLumi_ = ps.getUntrackedParameter<int>("fitEveryNLumi", -1);
133  resetFitNLumi_ = ps.getUntrackedParameter<int>("resetEveryNLumi", -1);
134  fitPVNLumi_ = ps.getUntrackedParameter<int>("fitPVEveryNLumi", -1);
135  resetPVNLumi_ = ps.getUntrackedParameter<int>("resetPVEveryNLumi", -1);
136  deltaSigCut_ = ps.getUntrackedParameter<double>("deltaSignificanceCut", 15);
137  debug_ = ps.getUntrackedParameter<bool>("Debug");
138  onlineMode_ = ps.getUntrackedParameter<bool>("OnlineMode");
139  min_Ntrks_ = ps.getParameter<ParameterSet>("BeamFitter").getUntrackedParameter<int>("MinimumInputTracks");
140  maxZ_ = ps.getParameter<ParameterSet>("BeamFitter").getUntrackedParameter<double>("MaximumZ");
141  minNrVertices_ = ps.getParameter<ParameterSet>("PVFitter").getUntrackedParameter<unsigned int>("minNrVerticesForFit");
142  minVtxNdf_ = ps.getParameter<ParameterSet>("PVFitter").getUntrackedParameter<double>("minVertexNdf");
143  minVtxWgt_ = ps.getParameter<ParameterSet>("PVFitter").getUntrackedParameter<double>("minVertexMeanWeight");
144  useLockRecords_ = ps.getUntrackedParameter<bool>("useLockRecords");
145 
146  if (!monitorName_.empty())
147  monitorName_ = monitorName_ + "/";
148 
149  if (fitNLumi_ <= 0)
150  fitNLumi_ = 1;
152  refBStime[0] = refBStime[1] = refPVtime[0] = refPVtime[1] = 0;
153  maxZ_ = std::fabs(maxZ_);
154  lastlumi_ = 0;
155  nextlumi_ = 0;
156  processed_ = false;
157 
158  rndm_ = new TRandom3(0);
159 }
160 
161 //--------------------------------------------------------
162 namespace {
163  /*The order of the enums is identical to the order in which
164  MonitorElements are added to hs
165  */
166  enum Hists {
167  k_x0_lumi,
168  k_x0_lumi_all,
169  k_y0_lumi,
170  k_y0_lumi_all,
171  k_z0_lumi,
172  k_z0_lumi_all,
173  k_sigmaX0_lumi,
174  k_sigmaX0_lumi_all,
175  k_sigmaY0_lumi,
176  k_sigmaY0_lumi_all,
177  k_sigmaZ0_lumi,
178  k_sigmaZ0_lumi_all,
179  k_x0_time,
180  k_x0_time_all,
181  k_y0_time,
182  k_y0_time_all,
183  k_z0_time,
184  k_z0_time_all,
185  k_sigmaX0_time,
186  k_sigmaX0_time_all,
187  k_sigmaY0_time,
188  k_sigmaY0_time_all,
189  k_sigmaZ0_time,
190  k_sigmaZ0_time_all,
191  k_PVx_lumi,
192  k_PVx_lumi_all,
193  k_PVy_lumi,
194  k_PVy_lumi_all,
195  k_PVz_lumi,
196  k_PVz_lumi_all,
197  k_PVx_time,
198  k_PVx_time_all,
199  k_PVy_time,
200  k_PVy_time_all,
201  k_PVz_time,
202  k_PVz_time_all,
203  kNumHists
204  };
205 } // namespace
206 
210  }
211 }
212 
213 void FakeBeamMonitor::bookHistograms(DQMStore::IBooker& iBooker, edm::Run const& iRun, edm::EventSetup const& iSetup) {
214  frun = iRun.run();
215  ftimestamp = iRun.beginTime().value();
216  tmpTime = ftimestamp >> 32;
218  char eventTime[64];
219  formatFitTime(eventTime, tmpTime);
220  edm::LogInfo("FakeBeamMonitor") << "TimeOffset = " << eventTime << std::endl;
221  TDatime da(eventTime);
222  if (debug_) {
223  edm::LogInfo("FakeBeamMonitor") << "TimeOffset = ";
224  da.Print();
225  }
226  auto daTime = da.Convert(kTRUE);
227 
228  // book some histograms here
229 
230  // create and cd into new folder
231  iBooker.setCurrentFolder(monitorName_ + "Fit");
232 
233  h_nTrk_lumi = iBooker.book1D("nTrk_lumi", "Num. of selected tracks vs lumi (Fit)", 20, 0.5, 20.5);
234  h_nTrk_lumi->setAxisTitle("Lumisection", 1);
235  h_nTrk_lumi->setAxisTitle("Num of Tracks for Fit", 2);
236 
237  //store vtx vs lumi for monitoring why fits fail
238  h_nVtx_lumi = iBooker.book1D("nVtx_lumi", "Num. of selected Vtx vs lumi (Fit)", 20, 0.5, 20.5);
239  h_nVtx_lumi->setAxisTitle("Lumisection", 1);
240  h_nVtx_lumi->setAxisTitle("Num of Vtx for Fit", 2);
241 
242  h_nVtx_lumi_all = iBooker.book1D("nVtx_lumi_all", "Num. of selected Vtx vs lumi (Fit) all", 20, 0.5, 20.5);
243  h_nVtx_lumi_all->getTH1()->SetCanExtend(TH1::kAllAxes);
244  h_nVtx_lumi_all->setAxisTitle("Lumisection", 1);
245  h_nVtx_lumi_all->setAxisTitle("Num of Vtx for Fit", 2);
246 
247  h_d0_phi0 = iBooker.bookProfile(
248  "d0_phi0", "d_{0} vs. #phi_{0} (Selected Tracks)", phiBin_, phiMin_, phiMax_, dxBin_, dxMin_, dxMax_, "");
249  h_d0_phi0->setAxisTitle("#phi_{0} (rad)", 1);
250  h_d0_phi0->setAxisTitle("d_{0} (cm)", 2);
251 
252  h_vx_vy = iBooker.book2D(
253  "trk_vx_vy", "Vertex (PCA) position of selected tracks", vxBin_, vxMin_, vxMax_, vxBin_, vxMin_, vxMax_);
254  h_vx_vy->setOption("COLZ");
255  // h_vx_vy->getTH1()->SetCanExtend(TH1::kAllAxes);
256  h_vx_vy->setAxisTitle("x coordinate of input track at PCA (cm)", 1);
257  h_vx_vy->setAxisTitle("y coordinate of input track at PCA (cm)", 2);
258 
259  {
260  TDatime* da = new TDatime();
261  gStyle->SetTimeOffset(da->Convert(kTRUE));
262  }
263 
264  const int nvar_ = 6;
265  string coord[nvar_] = {"x", "y", "z", "sigmaX", "sigmaY", "sigmaZ"};
266  string label[nvar_] = {
267  "x_{0} (cm)", "y_{0} (cm)", "z_{0} (cm)", "#sigma_{X_{0}} (cm)", "#sigma_{Y_{0}} (cm)", "#sigma_{Z_{0}} (cm)"};
268 
269  hs.clear();
270  hs.reserve(kNumHists);
271  for (int i = 0; i < 4; i++) {
272  iBooker.setCurrentFolder(monitorName_ + "Fit");
273  for (int ic = 0; ic < nvar_; ++ic) {
274  TString histName(coord[ic]);
275  TString histTitle(coord[ic]);
276  string ytitle(label[ic]);
277  string xtitle("");
278  string options("E1");
279  bool createHisto = true;
280  switch (i) {
281  case 1: // BS vs time
282  histName += "0_time";
283  xtitle = "Time [UTC]";
284  if (ic < 3)
285  histTitle += " coordinate of beam spot vs time (Fit)";
286  else
287  histTitle = histTitle.Insert(5, " ") + " of beam spot vs time (Fit)";
288  break;
289  case 2: // PV vs lumi
290  if (ic < 3) {
291  iBooker.setCurrentFolder(monitorName_ + "PrimaryVertex");
292  histName.Insert(0, "PV");
293  histName += "_lumi";
294  histTitle.Insert(0, "Avg. ");
295  histTitle += " position of primary vtx vs lumi";
296  xtitle = "Lumisection";
297  ytitle.insert(0, "PV");
298  ytitle += " #pm #sigma_{PV";
299  ytitle += coord[ic];
300  ytitle += "} (cm)";
301  } else
302  createHisto = false;
303  break;
304  case 3: // PV vs time
305  if (ic < 3) {
306  iBooker.setCurrentFolder(monitorName_ + "PrimaryVertex");
307  histName.Insert(0, "PV");
308  histName += "_time";
309  histTitle.Insert(0, "Avg. ");
310  histTitle += " position of primary vtx vs time";
311  xtitle = "Time [UTC]";
312  ytitle.insert(0, "PV");
313  ytitle += " #pm #sigma_{PV";
314  ytitle += coord[ic];
315  ytitle += "} (cm)";
316  } else
317  createHisto = false;
318  break;
319  default: // BS vs lumi
320  histName += "0_lumi";
321  xtitle = "Lumisection";
322  if (ic < 3)
323  histTitle += " coordinate of beam spot vs lumi (Fit)";
324  else
325  histTitle = histTitle.Insert(5, " ") + " of beam spot vs lumi (Fit)";
326  break;
327  }
328  if (createHisto) {
329  edm::LogInfo("FakeBeamMonitor") << "hitsName = " << histName << "; histTitle = " << histTitle << std::endl;
330  auto tmpHs = iBooker.book1D(histName, histTitle, 40, 0.5, 40.5);
331  hs.push_back(tmpHs);
332  tmpHs->setAxisTitle(xtitle, 1);
333  tmpHs->setAxisTitle(ytitle, 2);
334  tmpHs->getTH1()->SetOption("E1");
335  if (histName.Contains("time")) {
336  //int nbins = (intervalInSec_/23 > 0 ? intervalInSec_/23 : 40);
337  tmpHs->getTH1()->SetBins(intervalInSec_, 0.5, intervalInSec_ + 0.5);
338  tmpHs->setAxisTimeDisplay(1);
339  tmpHs->setAxisTimeFormat("%H:%M:%S", 1);
340  tmpHs->getTH1()->GetXaxis()->SetTimeOffset(daTime);
341  }
342  histName += "_all";
343  histTitle += " all";
344  tmpHs = iBooker.book1D(histName, histTitle, 40, 0.5, 40.5);
345  hs.push_back(tmpHs);
346  tmpHs->getTH1()->SetCanExtend(TH1::kAllAxes);
347  tmpHs->setAxisTitle(xtitle, 1);
348  tmpHs->setAxisTitle(ytitle, 2);
349  tmpHs->getTH1()->SetOption("E1");
350  if (histName.Contains("time")) {
351  //int nbins = (intervalInSec_/23 > 0 ? intervalInSec_/23 : 40);
352  tmpHs->getTH1()->SetBins(intervalInSec_, 0.5, intervalInSec_ + 0.5);
353  tmpHs->setAxisTimeDisplay(1);
354  tmpHs->setAxisTimeFormat("%H:%M:%S", 1);
355  tmpHs->getTH1()->GetXaxis()->SetTimeOffset(daTime);
356  }
357  }
358  }
359  }
360  assert(hs.size() == kNumHists);
361  assert(0 == strcmp(hs[k_sigmaY0_time]->getTH1()->GetName(), "sigmaY0_time"));
362  assert(0 == strcmp(hs[k_PVz_lumi_all]->getTH1()->GetName(), "PVz_lumi_all"));
363 
364  iBooker.setCurrentFolder(monitorName_ + "Fit");
365 
366  h_trk_z0 = iBooker.book1D("trk_z0", "z_{0} of selected tracks", dzBin_, dzMin_, dzMax_);
367  h_trk_z0->setAxisTitle("z_{0} of selected tracks (cm)", 1);
368 
369  h_vx_dz = iBooker.bookProfile(
370  "vx_dz", "v_{x} vs. dz of selected tracks", dzBin_, dzMin_, dzMax_, dxBin_, dxMin_, dxMax_, "");
371  h_vx_dz->setAxisTitle("dz (cm)", 1);
372  h_vx_dz->setAxisTitle("x coordinate of input track at PCA (cm)", 2);
373 
374  h_vy_dz = iBooker.bookProfile(
375  "vy_dz", "v_{y} vs. dz of selected tracks", dzBin_, dzMin_, dzMax_, dxBin_, dxMin_, dxMax_, "");
376  h_vy_dz->setAxisTitle("dz (cm)", 1);
377  h_vy_dz->setAxisTitle("y coordinate of input track at PCA (cm)", 2);
378 
379  h_x0 = iBooker.book1D("BeamMonitorFeedBack_x0", "x coordinate of beam spot (Fit)", 100, -0.01, 0.01);
380  h_x0->setAxisTitle("x_{0} (cm)", 1);
381  h_x0->getTH1()->SetCanExtend(TH1::kAllAxes);
382 
383  h_y0 = iBooker.book1D("BeamMonitorFeedBack_y0", "y coordinate of beam spot (Fit)", 100, -0.01, 0.01);
384  h_y0->setAxisTitle("y_{0} (cm)", 1);
385  h_y0->getTH1()->SetCanExtend(TH1::kAllAxes);
386 
387  h_z0 = iBooker.book1D("BeamMonitorFeedBack_z0", "z coordinate of beam spot (Fit)", dzBin_, dzMin_, dzMax_);
388  h_z0->setAxisTitle("z_{0} (cm)", 1);
389  h_z0->getTH1()->SetCanExtend(TH1::kAllAxes);
390 
391  h_sigmaX0 = iBooker.book1D("BeamMonitorFeedBack_sigmaX0", "sigma x0 of beam spot (Fit)", 100, 0, 0.05);
392  h_sigmaX0->setAxisTitle("#sigma_{X_{0}} (cm)", 1);
393  h_sigmaX0->getTH1()->SetCanExtend(TH1::kAllAxes);
394 
395  h_sigmaY0 = iBooker.book1D("BeamMonitorFeedBack_sigmaY0", "sigma y0 of beam spot (Fit)", 100, 0, 0.05);
396  h_sigmaY0->setAxisTitle("#sigma_{Y_{0}} (cm)", 1);
397  h_sigmaY0->getTH1()->SetCanExtend(TH1::kAllAxes);
398 
399  h_sigmaZ0 = iBooker.book1D("BeamMonitorFeedBack_sigmaZ0", "sigma z0 of beam spot (Fit)", 100, 0, 10);
400  h_sigmaZ0->setAxisTitle("#sigma_{Z_{0}} (cm)", 1);
401  h_sigmaZ0->getTH1()->SetCanExtend(TH1::kAllAxes);
402 
403  // Histograms of all reco tracks (without cuts):
404  h_trkPt = iBooker.book1D("trkPt", "p_{T} of all reco'd tracks (no selection)", 200, 0., 50.);
405  h_trkPt->setAxisTitle("p_{T} (GeV/c)", 1);
406 
407  h_trkVz = iBooker.book1D("trkVz", "Z coordinate of PCA of all reco'd tracks (no selection)", dzBin_, dzMin_, dzMax_);
408  h_trkVz->setAxisTitle("V_{Z} (cm)", 1);
409 
410  cutFlowTable = iBooker.book1D("cutFlowTable", "Cut flow table of track selection", 9, 0, 9);
411 
412  // Results of previous good fit:
413  fitResults = iBooker.book2D("fitResults", "Results of previous good beam fit", 2, 0, 2, 8, 0, 8);
414  fitResults->setAxisTitle("Fitted Beam Spot (cm)", 1);
415  fitResults->setBinLabel(8, "x_{0}", 2);
416  fitResults->setBinLabel(7, "y_{0}", 2);
417  fitResults->setBinLabel(6, "z_{0}", 2);
418  fitResults->setBinLabel(5, "#sigma_{Z}", 2);
419  fitResults->setBinLabel(4, "#frac{dx}{dz} (rad)", 2);
420  fitResults->setBinLabel(3, "#frac{dy}{dz} (rad)", 2);
421  fitResults->setBinLabel(2, "#sigma_{X}", 2);
422  fitResults->setBinLabel(1, "#sigma_{Y}", 2);
423  fitResults->setBinLabel(1, "Mean", 1);
424  fitResults->setBinLabel(2, "Stat. Error", 1);
425  fitResults->getTH1()->SetOption("text");
426 
427  // Histos of PrimaryVertices:
428  iBooker.setCurrentFolder(monitorName_ + "PrimaryVertex");
429 
430  h_nVtx = iBooker.book1D("vtxNbr", "Reconstructed Vertices(non-fake) in all Event", 60, -0.5, 59.5);
431  h_nVtx->setAxisTitle("Num. of reco. vertices", 1);
432 
433  //For one Trigger only
434  h_nVtx_st = iBooker.book1D("vtxNbr_SelectedTriggers", "Reconstructed Vertices(non-fake) in Events", 60, -0.5, 59.5);
435  //h_nVtx_st->setAxisTitle("Num. of reco. vertices for Un-Prescaled Jet Trigger",1);
436 
437  // Monitor only the PV with highest sum pt of assoc. trks:
438  h_PVx[0] = iBooker.book1D("PVX", "x coordinate of Primary Vtx", 50, -0.01, 0.01);
439  h_PVx[0]->setAxisTitle("PVx (cm)", 1);
440  h_PVx[0]->getTH1()->SetCanExtend(TH1::kAllAxes);
441 
442  h_PVy[0] = iBooker.book1D("PVY", "y coordinate of Primary Vtx", 50, -0.01, 0.01);
443  h_PVy[0]->setAxisTitle("PVy (cm)", 1);
444  h_PVy[0]->getTH1()->SetCanExtend(TH1::kAllAxes);
445 
446  h_PVz[0] = iBooker.book1D("PVZ", "z coordinate of Primary Vtx", dzBin_, dzMin_, dzMax_);
447  h_PVz[0]->setAxisTitle("PVz (cm)", 1);
448 
449  h_PVx[1] = iBooker.book1D("PVXFit", "x coordinate of Primary Vtx (Last Fit)", 50, -0.01, 0.01);
450  h_PVx[1]->setAxisTitle("PVx (cm)", 1);
451  h_PVx[1]->getTH1()->SetCanExtend(TH1::kAllAxes);
452 
453  h_PVy[1] = iBooker.book1D("PVYFit", "y coordinate of Primary Vtx (Last Fit)", 50, -0.01, 0.01);
454  h_PVy[1]->setAxisTitle("PVy (cm)", 1);
455  h_PVy[1]->getTH1()->SetCanExtend(TH1::kAllAxes);
456 
457  h_PVz[1] = iBooker.book1D("PVZFit", "z coordinate of Primary Vtx (Last Fit)", dzBin_, dzMin_, dzMax_);
458  h_PVz[1]->setAxisTitle("PVz (cm)", 1);
459 
460  h_PVxz = iBooker.bookProfile("PVxz", "PVx vs. PVz", dzBin_ / 2, dzMin_, dzMax_, dxBin_ / 2, dxMin_, dxMax_, "");
461  h_PVxz->setAxisTitle("PVz (cm)", 1);
462  h_PVxz->setAxisTitle("PVx (cm)", 2);
463 
464  h_PVyz = iBooker.bookProfile("PVyz", "PVy vs. PVz", dzBin_ / 2, dzMin_, dzMax_, dxBin_ / 2, dxMin_, dxMax_, "");
465  h_PVyz->setAxisTitle("PVz (cm)", 1);
466  h_PVyz->setAxisTitle("PVy (cm)", 2);
467 
468  // Results of previous good fit:
469  pvResults = iBooker.book2D("pvResults", "Results of fitting Primary Vertices", 2, 0, 2, 6, 0, 6);
470  pvResults->setAxisTitle("Fitted Primary Vertex (cm)", 1);
471  pvResults->setBinLabel(6, "PVx", 2);
472  pvResults->setBinLabel(5, "PVy", 2);
473  pvResults->setBinLabel(4, "PVz", 2);
474  pvResults->setBinLabel(3, "#sigma_{X}", 2);
475  pvResults->setBinLabel(2, "#sigma_{Y}", 2);
476  pvResults->setBinLabel(1, "#sigma_{Z}", 2);
477  pvResults->setBinLabel(1, "Mean", 1);
478  pvResults->setBinLabel(2, "Stat. Error", 1);
479  pvResults->getTH1()->SetOption("text");
480 
481  // Summary plots:
482  iBooker.setCurrentFolder(monitorName_ + "EventInfo");
483 
484  reportSummary = iBooker.bookFloat("reportSummary");
485  if (reportSummary)
486  reportSummary->Fill(std::numeric_limits<double>::quiet_NaN());
487 
488  char histo[20];
489  iBooker.setCurrentFolder(monitorName_ + "EventInfo/reportSummaryContents");
490  for (int n = 0; n < nFitElements_; n++) {
491  switch (n) {
492  case 0:
493  sprintf(histo, "x0_status");
494  break;
495  case 1:
496  sprintf(histo, "y0_status");
497  break;
498  case 2:
499  sprintf(histo, "z0_status");
500  break;
501  }
502  reportSummaryContents[n] = iBooker.bookFloat(histo);
503  }
504 
505  for (int i = 0; i < nFitElements_; i++) {
506  summaryContent_[i] = 0.;
507  reportSummaryContents[i]->Fill(std::numeric_limits<double>::quiet_NaN());
508  }
509 
510  iBooker.setCurrentFolder(monitorName_ + "EventInfo");
511 
512  reportSummaryMap = iBooker.book2D("reportSummaryMap", "Beam Spot Summary Map", 1, 0, 1, 3, 0, 3);
514  reportSummaryMap->setAxisTitle("Fitted Beam Spot", 2);
515  reportSummaryMap->setBinLabel(1, " ", 1);
516  reportSummaryMap->setBinLabel(1, "x_{0}", 2);
517  reportSummaryMap->setBinLabel(2, "y_{0}", 2);
518  reportSummaryMap->setBinLabel(3, "z_{0}", 2);
519  for (int i = 0; i < nFitElements_; i++) {
520  reportSummaryMap->setBinContent(1, i + 1, -1.);
521  }
522 }
523 
524 //--------------------------------------------------------
526  // start DB logger
527  DBloggerReturn_ = 0;
530  onlineDbService_->logger().logInfo() << "FakeBeamMonitor::beginLuminosityBlock - LS: " << lumiSeg.luminosityBlock()
531  << " - Run: " << lumiSeg.getRun().run();
532  }
533 
534  int nthlumi = lumiSeg.luminosityBlock();
535  const edm::TimeValue_t fbegintimestamp = lumiSeg.beginTime().value();
536  const std::time_t ftmptime = fbegintimestamp >> 32;
537 
538  if (countLumi_ == 0 && (!processed_)) {
540  refBStime[0] = refPVtime[0] = ftmptime;
541  mapBeginBSLS[countLumi_] = nthlumi;
542  mapBeginPVLS[countLumi_] = nthlumi;
543  mapBeginBSTime[countLumi_] = ftmptime;
544  mapBeginPVTime[countLumi_] = ftmptime;
545  } //for the first record
546 
547  if (nthlumi > nextlumi_) {
548  if (processed_) {
549  countLumi_++;
550  //store here them will need when we remove the first one of Last N LS
551  mapBeginBSLS[countLumi_] = nthlumi;
552  mapBeginPVLS[countLumi_] = nthlumi;
553  mapBeginBSTime[countLumi_] = ftmptime;
554  mapBeginPVTime[countLumi_] = ftmptime;
555  } //processed passed but not the first lumi
556  if ((!processed_) && countLumi_ != 0) {
557  mapBeginBSLS[countLumi_] = nthlumi;
558  mapBeginPVLS[countLumi_] = nthlumi;
559  mapBeginBSTime[countLumi_] = ftmptime;
560  mapBeginPVTime[countLumi_] = ftmptime;
561  } //processed fails for last lumi
562  } //nthLumi > nextlumi
563 
564  if (StartAverage_) {
565  //Just Make sure it get rest
566  refBStime[0] = 0;
567  refPVtime[0] = 0;
568  beginLumiOfPVFit_ = 0;
569  beginLumiOfBSFit_ = 0;
570 
571  if (debug_)
572  edm::LogInfo("FakeBeamMonitor") << " beginLuminosityBlock: Size of mapBeginBSLS before = "
573  << mapBeginBSLS.size() << endl;
574  if (nthlumi >
575  nextlumi_) { //this make sure that it does not take into account this lumi for fitting and only look forward for new lumi
576  //as countLumi also remains the same so map value get overwritten once return to normal running.
577  //even if few LS are misssing and DQM module do not sees them then it catchs up again
578  map<int, int>::iterator itbs = mapBeginBSLS.begin();
579  map<int, int>::iterator itpv = mapBeginPVLS.begin();
580  map<int, std::time_t>::iterator itbstime = mapBeginBSTime.begin();
581  map<int, std::time_t>::iterator itpvtime = mapBeginPVTime.begin();
582 
583  if (processed_) { // otherwise if false then LS range of fit get messed up because we don't remove trk/pvs but we remove LS begin value . This prevent it as it happened if LS is there but no event are processed for some reason
584  mapBeginBSLS.erase(itbs);
585  mapBeginPVLS.erase(itpv);
586  mapBeginBSTime.erase(itbstime);
587  mapBeginPVTime.erase(itpvtime);
588  }
589  /*//not sure if want this or not ??
590  map<int, int>::iterator itgapb=mapBeginBSLS.begin();
591  map<int, int>::iterator itgape=mapBeginBSLS.end(); itgape--;
592  countGapLumi_ = ( (itgape->second) - (itgapb->second) );
593  //if we see Gap more than then 2*resetNFitLumi !!!!!!!
594  //for example if 10-15 is fitted and if 16-25 are missing then we next fit will be for range 11-26 but BS can change in between
595  // so better start as fresh and reset everything like starting in the begining!
596  if(countGapLumi_ >= 2*resetFitNLumi_){RestartFitting(); mapBeginBSLS[countLumi_] = nthlumi;}
597  */
598  }
599 
600  if (debug_)
601  edm::LogInfo("FakeBeamMonitor") << " beginLuminosityBlock:: Size of mapBeginBSLS After = " << mapBeginBSLS.size()
602  << endl;
603 
604  map<int, int>::iterator bbs = mapBeginBSLS.begin();
605  map<int, int>::iterator bpv = mapBeginPVLS.begin();
606  map<int, std::time_t>::iterator bbst = mapBeginBSTime.begin();
607  map<int, std::time_t>::iterator bpvt = mapBeginPVTime.begin();
608 
609  if (beginLumiOfPVFit_ == 0)
610  beginLumiOfPVFit_ = bpv->second; //new begin time after removing the LS
611  if (beginLumiOfBSFit_ == 0)
612  beginLumiOfBSFit_ = bbs->second;
613  if (refBStime[0] == 0)
614  refBStime[0] = bbst->second;
615  if (refPVtime[0] == 0)
616  refPVtime[0] = bpvt->second;
617 
618  } //same logic for average fit as above commented line
619 
620  map<int, std::time_t>::iterator nbbst = mapBeginBSTime.begin();
621  map<int, std::time_t>::iterator nbpvt = mapBeginPVTime.begin();
622 
623  if (onlineMode_ && (nthlumi < nextlumi_))
624  return;
625 
626  if (onlineMode_) {
627  if (nthlumi > nextlumi_) {
628  if (countLumi_ != 0 && processed_)
629  FitAndFill(lumiSeg, lastlumi_, nextlumi_, nthlumi);
630  nextlumi_ = nthlumi;
631  edm::LogInfo("FakeBeamMonitor") << "beginLuminosityBlock:: Next Lumi to Fit: " << nextlumi_ << endl;
632  if ((StartAverage_) && refBStime[0] == 0)
633  refBStime[0] = nbbst->second;
634  if ((StartAverage_) && refPVtime[0] == 0)
635  refPVtime[0] = nbpvt->second;
636  }
637  } else {
638  if (processed_)
639  FitAndFill(lumiSeg, lastlumi_, nextlumi_, nthlumi);
640  nextlumi_ = nthlumi;
641  edm::LogInfo("FakeBeamMonitor") << " beginLuminosityBlock:: Next Lumi to Fit: " << nextlumi_ << endl;
642  if ((StartAverage_) && refBStime[0] == 0)
643  refBStime[0] = nbbst->second;
644  if ((StartAverage_) && refPVtime[0] == 0)
645  refPVtime[0] = nbpvt->second;
646  }
647 
648  //countLumi_++;
649  if (processed_)
650  processed_ = false;
651  edm::LogInfo("FakeBeamMonitor") << " beginLuminosityBlock:: Begin of Lumi: " << nthlumi << endl;
652 }
653 
654 // ----------------------------------------------------------
655 void FakeBeamMonitor::analyze(const Event& iEvent, const EventSetup& iSetup) {
656  const int nthlumi = iEvent.luminosityBlock();
657 
658  if (onlineMode_ && (nthlumi < nextlumi_)) {
659  edm::LogInfo("FakeBeamMonitor") << "analyze:: Spilt event from previous lumi section!" << std::endl;
660  return;
661  }
662  if (onlineMode_ && (nthlumi > nextlumi_)) {
663  edm::LogInfo("FakeBeamMonitor") << "analyze:: Spilt event from next lumi section!!!" << std::endl;
664  return;
665  }
666 
667  countEvt_++;
668  // theBeamFitter->readEvent(
669  // iEvent); //Remember when track fitter read the event in the same place the PVFitter read the events !!!!!!!!!
670 
671  // Handle<reco::BeamSpot> recoBeamSpotHandle;
672  // iEvent.getByToken(bsSrc_, recoBeamSpotHandle);
673  // refBS = *recoBeamSpotHandle;
674 
675  // //------Cut Flow Table filled every event!--------------------------------------
676  // {
677  // // Make a copy of the cut flow table from the beam fitter.
678  // auto tmphisto = static_cast<TH1F*>(theBeamFitter->getCutFlow());
679  // cutFlowTable->getTH1()->SetBins(
680  // tmphisto->GetNbinsX(), tmphisto->GetXaxis()->GetXmin(), tmphisto->GetXaxis()->GetXmax());
681  // // Update the bin labels
682  // if (countEvt_ == 1) // SetLabel just once
683  // for (int n = 0; n < tmphisto->GetNbinsX(); n++)
684  // cutFlowTable->setBinLabel(n + 1, tmphisto->GetXaxis()->GetBinLabel(n + 1), 1);
685  // cutFlowTable->Reset();
686  // cutFlowTable->getTH1()->Add(tmphisto);
687  // }
688 
689  //----Reco tracks -------------------------------------
690  // Handle<reco::TrackCollection> TrackCollection;
691  // iEvent.getByToken(tracksLabel_, TrackCollection);
692  // const reco::TrackCollection* tracks = TrackCollection.product();
693  // for (reco::TrackCollection::const_iterator track = tracks->begin(); track != tracks->end(); ++track) {
694  // h_trkPt->Fill(track->pt()); //no need to change here for average bs
695  // h_trkVz->Fill(track->vz());
696  // }
697 
698  //-------HLT Trigger --------------------------------
699  // edm::Handle<TriggerResults> triggerResults;
700  // bool JetTrigPass = false;
701  // if (iEvent.getByToken(hltSrc_, triggerResults)) {
702  // const edm::TriggerNames& trigNames = iEvent.triggerNames(*triggerResults);
703  // for (unsigned int i = 0; i < triggerResults->size(); i++) {
704  // const std::string& trigName = trigNames.triggerName(i);
705  //
706  // if (JetTrigPass)
707  // continue;
708  //
709  // for (size_t t = 0; t < jetTrigger_.size(); ++t) {
710  // if (JetTrigPass)
711  // continue;
712  //
713  // string string_search(jetTrigger_[t]);
714  // size_t found = trigName.find(string_search);
715  //
716  // if (found != string::npos) {
717  // int thisTrigger_ = trigNames.triggerIndex(trigName);
718  // if (triggerResults->accept(thisTrigger_))
719  // JetTrigPass = true;
720  // } //if trigger found
721  // } //for(t=0;..)
722  // } //for(i=0; ..)
723  // } //if trigger colleciton exist)
724 
725  //------ Primary Vertices-------
726  // edm::Handle<reco::VertexCollection> PVCollection;
727 
728  // if (iEvent.getByToken(pvSrc_, PVCollection)) {
729  int nPVcount = 0;
730  int nPVcount_ST = 0; //For Single Trigger(hence ST)
731 
732  // for (reco::VertexCollection::const_iterator pv = PVCollection->begin(); pv != PVCollection->end(); ++pv) {
733  for (int tmp_idx = 0; tmp_idx < 10; tmp_idx++) {
734  //--- vertex selection
735  // if (pv->isFake() || pv->tracksSize() == 0)
736  // continue;
737  nPVcount++; // count non fake pv:
738 
739  //if (JetTrigPass)
740  nPVcount_ST++; //non-fake pv with a specific trigger
741 
742  // if (pv->ndof() < minVtxNdf_ || (pv->ndof() + 3.) / pv->tracksSize() < 2 * minVtxWgt_)
743  // continue;
744 
745  //Fill this map to store xyx for pv so that later we can remove the first one for run aver
746  mapPVx[countLumi_].push_back(tmp_idx);
747  mapPVy[countLumi_].push_back(tmp_idx);
748  mapPVz[countLumi_].push_back(tmp_idx);
749 
750  // if (!StartAverage_) { //for first N LS
751  // h_PVx[0]->Fill(pv->x());
752  // h_PVy[0]->Fill(pv->y());
753  // h_PVz[0]->Fill(pv->z());
754  // h_PVxz->Fill(pv->z(), pv->x());
755  // h_PVyz->Fill(pv->z(), pv->y());
756  // } //for first N LiS
757  // else {
758  // h_PVxz->Fill(pv->z(), pv->x());
759  // h_PVyz->Fill(pv->z(), pv->y());
760  // }
761 
762  } //loop over pvs
763 
764  // h_nVtx->Fill(nPVcount * 1.); //no need to change it for average BS
765 
766  mapNPV[countLumi_].push_back((nPVcount_ST));
767 
768  // if (!StartAverage_) {
769  // h_nVtx_st->Fill(nPVcount_ST * 1.);
770  // }
771 
772  // } //if pv collection is availaable
773 
774  if (StartAverage_) {
775  map<int, std::vector<float> >::iterator itpvx = mapPVx.begin();
776  map<int, std::vector<float> >::iterator itpvy = mapPVy.begin();
777  map<int, std::vector<float> >::iterator itpvz = mapPVz.begin();
778 
779  map<int, std::vector<int> >::iterator itbspvinfo = mapNPV.begin();
780 
781  if ((int)mapPVx.size() > resetFitNLumi_) { //sometimes the events is not there but LS is there!
782  mapPVx.erase(itpvx);
783  mapPVy.erase(itpvy);
784  mapPVz.erase(itpvz);
785  mapNPV.erase(itbspvinfo);
786  } //loop over Last N lumi collected
787 
788  } //StartAverage==true
789 
790  processed_ = true;
791 }
792 
793 //--------------------------------------------------------
795  int nthlumi = lumiSeg.id().luminosityBlock();
796  edm::LogInfo("FakeBeamMonitor") << "endLuminosityBlock:: Lumi of the last event before endLuminosityBlock: "
797  << nthlumi << endl;
798 
799  if (onlineMode_ && nthlumi < nextlumi_)
800  return;
801  const edm::TimeValue_t fendtimestamp = lumiSeg.endTime().value();
802  const std::time_t fendtime = fendtimestamp >> 32;
803  tmpTime = refBStime[1] = refPVtime[1] = fendtime;
804 
805  // end DB logger
807  onlineDbService_->logger().logInfo() << "FakeBeamMonitor::endLuminosityBlock";
809  }
810 }
811 
812 //--------------------------------------------------------
813 void FakeBeamMonitor::FitAndFill(const LuminosityBlock& lumiSeg, int& lastlumi, int& nextlumi, int& nthlumi) {
814  if (onlineMode_ && (nthlumi <= nextlumi))
815  return;
816 
817  //set the correct run number when no event in the LS for fake output
818  // if ((processed_) && theBeamFitter->getRunNumber() != frun)
819  // theBeamFitter->setRun(frun);
820 
821  int currentlumi = nextlumi;
822  edm::LogInfo("FakeBeamMonitor") << "FitAndFill:: Lumi of the current fit: " << currentlumi << endl;
823  lastlumi = currentlumi;
824  endLumiOfBSFit_ = currentlumi;
825  endLumiOfPVFit_ = currentlumi;
826 
827  //---------Fix for Runninv average-------------
828  mapLSPVStoreSize[countLumi_] = 10; //theBeamFitter->getPVvectorSize();
829 
830  // if (StartAverage_) {
831  // std::map<int, std::size_t>::iterator rmLSPVi = mapLSPVStoreSize.begin();
832  // size_t SizeToRemovePV = rmLSPVi->second;
833  // for (std::map<int, std::size_t>::iterator rmLSPVe = mapLSPVStoreSize.end(); ++rmLSPVi != rmLSPVe;)
834  // rmLSPVi->second -= SizeToRemovePV;
835  //
836  // theBeamFitter->resizePVvector(SizeToRemovePV);
837  //
838  // map<int, std::size_t>::iterator tmpItpv = mapLSPVStoreSize.begin();
839  // mapLSPVStoreSize.erase(tmpItpv);
840  // }
841  // if (debug_)
842  // edm::LogInfo("BeamMonitor") << "FitAndFill:: Size of thePVvector After removing the PVs = "
843  // << theBeamFitter->getPVvectorSize() << endl;
844 
845  //lets filt the PV for GUI here: It was in analyzer in preivous versiton but moved here due to absence of event in some lumis, works OK
846  bool resetHistoFlag_ = false;
847  if ((int)mapPVx.size() >= resetFitNLumi_ && (StartAverage_)) {
848  h_PVx[0]->Reset();
849  h_PVy[0]->Reset();
850  h_PVz[0]->Reset();
851  h_nVtx_st->Reset();
852  resetHistoFlag_ = true;
853  }
854 
855  int MaxPVs = 0;
856  int countEvtLastNLS_ = 0;
857  int countTotPV_ = 0;
858 
859  std::map<int, std::vector<int> >::iterator mnpv = mapNPV.begin();
860  std::map<int, std::vector<float> >::iterator mpv2 = mapPVy.begin();
861  std::map<int, std::vector<float> >::iterator mpv3 = mapPVz.begin();
862 
863  for (std::map<int, std::vector<float> >::iterator mpv1 = mapPVx.begin(); mpv1 != mapPVx.end();
864  ++mpv1, ++mpv2, ++mpv3, ++mnpv) {
865  std::vector<float>::iterator mpvs2 = (mpv2->second).begin();
866  std::vector<float>::iterator mpvs3 = (mpv3->second).begin();
867  for (std::vector<float>::iterator mpvs1 = (mpv1->second).begin(); mpvs1 != (mpv1->second).end();
868  ++mpvs1, ++mpvs2, ++mpvs3) {
869  if (resetHistoFlag_) {
870  h_PVx[0]->Fill(*mpvs1); //these histogram are reset after StartAverage_ flag is ON
871  h_PVy[0]->Fill(*mpvs2);
872  h_PVz[0]->Fill(*mpvs3);
873  }
874  } //loop over second
875 
876  //Do the same here for nPV distr.
877  for (std::vector<int>::iterator mnpvs = (mnpv->second).begin(); mnpvs != (mnpv->second).end(); ++mnpvs) {
878  if ((*mnpvs > 0) && (resetHistoFlag_))
879  h_nVtx_st->Fill((*mnpvs) * (1.0));
880  countEvtLastNLS_++;
881  countTotPV_ += (*mnpvs);
882  if ((*mnpvs) > MaxPVs)
883  MaxPVs = (*mnpvs);
884  } //loop over second of mapNPV
885 
886  } //loop over last N lumis
887 
888  char tmpTitlePV[100];
889  sprintf(tmpTitlePV, "%s %i %s %i", "Num. of reco. vertices for LS: ", beginLumiOfPVFit_, " to ", endLumiOfPVFit_);
890  h_nVtx_st->setAxisTitle(tmpTitlePV, 1);
891 
892  std::vector<float> DipPVInfo_;
893  DipPVInfo_.clear();
894  //
895  // if (countTotPV_ != 0) {
896  // DipPVInfo_.push_back((float)countEvtLastNLS_);
897  // DipPVInfo_.push_back(h_nVtx_st->getMean());
898  // DipPVInfo_.push_back(h_nVtx_st->getMeanError());
899  // DipPVInfo_.push_back(h_nVtx_st->getRMS());
900  // DipPVInfo_.push_back(h_nVtx_st->getRMSError());
901  // DipPVInfo_.push_back((float)MaxPVs);
902  // DipPVInfo_.push_back((float)countTotPV_);
903  // MaxPVs = 0;
904  // } else {
905  // for (size_t i = 0; i < 7; i++) {
906  // if (i > 0) {
907  // DipPVInfo_.push_back(0.);
908  // } else {
909  // DipPVInfo_.push_back((float)countEvtLastNLS_);
910  // }
911  // }
912  // }
913  // theBeamFitter->SetPVInfo(DipPVInfo_);
914  DipPVInfo_.push_back(rndm_->Gaus(1000., 100.)); // Events used
915  DipPVInfo_.push_back(rndm_->Gaus(100., 10.)); // Mean PV
916  DipPVInfo_.push_back(rndm_->Gaus(10., 5.)); // Mean PV err
917  DipPVInfo_.push_back(rndm_->Gaus(10., 5.)); // Rms PV
918  DipPVInfo_.push_back(rndm_->Gaus(5., 3.)); // Rms PV err
919  DipPVInfo_.push_back(rndm_->Gaus(100., 10.)); // Max PVs
920  countEvtLastNLS_ = 0;
921 
922  if (onlineMode_) { // filling LS gap
923  // FIXME: need to add protection for the case if the gap is at the resetting LS!
924  const int countLS_bs = hs[k_x0_lumi]->getTH1()->GetEntries();
925  const int countLS_pv = hs[k_PVx_lumi]->getTH1()->GetEntries();
926  edm::LogInfo("FakeBeamMonitor") << "FitAndFill:: countLS_bs = " << countLS_bs << " ; countLS_pv = " << countLS_pv
927  << std::endl;
928  int LSgap_bs = currentlumi / fitNLumi_ - countLS_bs;
929  int LSgap_pv = currentlumi / fitPVNLumi_ - countLS_pv;
930  if (currentlumi % fitNLumi_ == 0)
931  LSgap_bs--;
932  if (currentlumi % fitPVNLumi_ == 0)
933  LSgap_pv--;
934  edm::LogInfo("FakeBeamMonitor") << "FitAndFill:: LSgap_bs = " << LSgap_bs << " ; LSgap_pv = " << LSgap_pv
935  << std::endl;
936  // filling previous fits if LS gap ever exists
937  for (int ig = 0; ig < LSgap_bs; ig++) {
938  hs[k_x0_lumi]->ShiftFillLast(0., 0., fitNLumi_); //x0 , x0err, fitNLumi_; see DQMCore....
939  hs[k_y0_lumi]->ShiftFillLast(0., 0., fitNLumi_);
940  hs[k_z0_lumi]->ShiftFillLast(0., 0., fitNLumi_);
941  hs[k_sigmaX0_lumi]->ShiftFillLast(0., 0., fitNLumi_);
942  hs[k_sigmaY0_lumi]->ShiftFillLast(0., 0., fitNLumi_);
943  hs[k_sigmaZ0_lumi]->ShiftFillLast(0., 0., fitNLumi_);
945  }
946  for (int ig = 0; ig < LSgap_pv; ig++) {
947  hs[k_PVx_lumi]->ShiftFillLast(0., 0., fitPVNLumi_);
948  hs[k_PVy_lumi]->ShiftFillLast(0., 0., fitPVNLumi_);
949  hs[k_PVz_lumi]->ShiftFillLast(0., 0., fitPVNLumi_);
950  }
951  const int previousLS = h_nTrk_lumi->getTH1()->GetEntries();
952  for (int i = 1; i < (currentlumi - previousLS);
953  i++) //if (current-previoius)= 1 then never go inside the for loop!!!!!!!!!!!
955  }
956 
957  edm::LogInfo("FakeBeamMonitor") << "FitAndFill:: Time lapsed since last scroll = " << tmpTime - refTime << std::endl;
958 
959  if (testScroll(tmpTime, refTime)) {
960  scrollTH1(hs[k_x0_time]->getTH1(), refTime);
961  scrollTH1(hs[k_y0_time]->getTH1(), refTime);
962  scrollTH1(hs[k_z0_time]->getTH1(), refTime);
963  scrollTH1(hs[k_sigmaX0_time]->getTH1(), refTime);
964  scrollTH1(hs[k_sigmaY0_time]->getTH1(), refTime);
965  scrollTH1(hs[k_sigmaZ0_time]->getTH1(), refTime);
966  scrollTH1(hs[k_PVx_time]->getTH1(), refTime);
967  scrollTH1(hs[k_PVy_time]->getTH1(), refTime);
968  scrollTH1(hs[k_PVz_time]->getTH1(), refTime);
969  }
970 
971  // bool doPVFit = false;
972  //
973  // if (fitPVNLumi_ > 0) {
974  // if (onlineMode_) {
975  // if (currentlumi % fitPVNLumi_ == 0)
976  // doPVFit = true;
977  // } else if (countLumi_ % fitPVNLumi_ == 0)
978  // doPVFit = true;
979  // } else
980  // doPVFit = true;
981  //
982  // if (doPVFit) {
983  edm::LogInfo("FakeBeamMonitor") << "FitAndFill:: Do PV Fitting for LS = " << beginLumiOfPVFit_ << " to "
984  << endLumiOfPVFit_ << std::endl;
985  // Primary Vertex Fit:
986  if (h_PVx[0]->getTH1()->GetEntries() > minNrVertices_) {
987  pvResults->Reset();
988  char tmpTitle[50];
989  sprintf(tmpTitle, "%s %i %s %i", "Fitted Primary Vertex (cm) of LS: ", beginLumiOfPVFit_, " to ", endLumiOfPVFit_);
990  pvResults->setAxisTitle(tmpTitle, 1);
991 
992  std::unique_ptr<TF1> fgaus{new TF1("fgaus", "gaus")};
993  double mean, width, meanErr, widthErr;
994  fgaus->SetLineColor(4);
995  h_PVx[0]->getTH1()->Fit(fgaus.get(), "QLM0");
996  mean = fgaus->GetParameter(1);
997  width = fgaus->GetParameter(2);
998  meanErr = fgaus->GetParError(1);
999  widthErr = fgaus->GetParError(2);
1000 
1001  hs[k_PVx_lumi]->ShiftFillLast(mean, width, fitPVNLumi_);
1002  hs[k_PVx_lumi_all]->setBinContent(currentlumi, mean);
1003  hs[k_PVx_lumi_all]->setBinError(currentlumi, width);
1004  int nthBin = tmpTime - refTime;
1005  if (nthBin < 0)
1006  edm::LogInfo("FakeBeamMonitor") << "FitAndFill:: Event time outside current range of time histograms!"
1007  << std::endl;
1008  if (nthBin > 0) {
1009  hs[k_PVx_time]->setBinContent(nthBin, mean);
1010  hs[k_PVx_time]->setBinError(nthBin, width);
1011  }
1012  int jthBin = tmpTime - startTime;
1013  if (jthBin > 0) {
1014  hs[k_PVx_time_all]->setBinContent(jthBin, mean);
1015  hs[k_PVx_time_all]->setBinError(jthBin, width);
1016  }
1017  pvResults->setBinContent(1, 6, mean);
1018  pvResults->setBinContent(1, 3, width);
1019  pvResults->setBinContent(2, 6, meanErr);
1020  pvResults->setBinContent(2, 3, widthErr);
1021 
1022  {
1023  // snap shot of the fit
1024  auto tmphisto = h_PVx[0]->getTH1F();
1025  h_PVx[1]->getTH1()->SetBins(
1026  tmphisto->GetNbinsX(), tmphisto->GetXaxis()->GetXmin(), tmphisto->GetXaxis()->GetXmax());
1027  h_PVx[1]->Reset();
1028  h_PVx[1]->getTH1()->Add(tmphisto);
1029  h_PVx[1]->getTH1()->Fit(fgaus.get(), "QLM");
1030  }
1031 
1032  h_PVy[0]->getTH1()->Fit(fgaus.get(), "QLM0");
1033  mean = fgaus->GetParameter(1);
1034  width = fgaus->GetParameter(2);
1035  meanErr = fgaus->GetParError(1);
1036  widthErr = fgaus->GetParError(2);
1037  hs[k_PVy_lumi]->ShiftFillLast(mean, width, fitPVNLumi_);
1038  hs[k_PVy_lumi_all]->setBinContent(currentlumi, mean);
1039  hs[k_PVy_lumi_all]->setBinError(currentlumi, width);
1040  if (nthBin > 0) {
1041  hs[k_PVy_time]->setBinContent(nthBin, mean);
1042  hs[k_PVy_time]->setBinError(nthBin, width);
1043  }
1044  if (jthBin > 0) {
1045  hs[k_PVy_time_all]->setBinContent(jthBin, mean);
1046  hs[k_PVy_time_all]->setBinError(jthBin, width);
1047  }
1048  pvResults->setBinContent(1, 5, mean);
1049  pvResults->setBinContent(1, 2, width);
1050  pvResults->setBinContent(2, 5, meanErr);
1051  pvResults->setBinContent(2, 2, widthErr);
1052  // snap shot of the fit
1053  {
1054  auto tmphisto = h_PVy[0]->getTH1F();
1055  h_PVy[1]->getTH1()->SetBins(
1056  tmphisto->GetNbinsX(), tmphisto->GetXaxis()->GetXmin(), tmphisto->GetXaxis()->GetXmax());
1057  h_PVy[1]->Reset();
1058  h_PVy[1]->getTH1()->Add(tmphisto);
1059  h_PVy[1]->getTH1()->Fit(fgaus.get(), "QLM");
1060  }
1061 
1062  h_PVz[0]->getTH1()->Fit(fgaus.get(), "QLM0");
1063  mean = fgaus->GetParameter(1);
1064  width = fgaus->GetParameter(2);
1065  meanErr = fgaus->GetParError(1);
1066  widthErr = fgaus->GetParError(2);
1067  hs[k_PVz_lumi]->ShiftFillLast(mean, width, fitPVNLumi_);
1068  hs[k_PVz_lumi_all]->setBinContent(currentlumi, mean);
1069  hs[k_PVz_lumi_all]->setBinError(currentlumi, width);
1070  if (nthBin > 0) {
1071  hs[k_PVz_time]->setBinContent(nthBin, mean);
1072  hs[k_PVz_time]->setBinError(nthBin, width);
1073  }
1074  if (jthBin > 0) {
1075  hs[k_PVz_time_all]->setBinContent(jthBin, mean);
1076  hs[k_PVz_time_all]->setBinError(jthBin, width);
1077  }
1078  pvResults->setBinContent(1, 4, mean);
1079  pvResults->setBinContent(1, 1, width);
1080  pvResults->setBinContent(2, 4, meanErr);
1081  pvResults->setBinContent(2, 1, widthErr);
1082  {
1083  // snap shot of the fit
1084  auto tmphisto = h_PVz[0]->getTH1F();
1085  h_PVz[1]->getTH1()->SetBins(
1086  tmphisto->GetNbinsX(), tmphisto->GetXaxis()->GetXmin(), tmphisto->GetXaxis()->GetXmax());
1087  h_PVz[1]->Reset();
1088  h_PVz[1]->getTH1()->Add(tmphisto);
1089  h_PVz[1]->getTH1()->Fit(fgaus.get(), "QLM");
1090  }
1091  } //check if found min Vertices
1092  // } //do PVfit
1093 
1094  if ((resetPVNLumi_ > 0 && countLumi_ == resetPVNLumi_) || StartAverage_) {
1095  beginLumiOfPVFit_ = 0;
1096  refPVtime[0] = 0;
1097  }
1098 
1099  //---------Readjustment of theBSvector, RefTime, beginLSofFit---------
1100  // vector<BSTrkParameters> theBSvector1 = theBeamFitter->getBSvector();
1101  // mapLSBSTrkSize[countLumi_] = (theBSvector1.size());
1102  size_t PreviousRecords = 0; //needed to fill nth record of tracks in GUI
1103 
1104  // if (StartAverage_) {
1105  // size_t SizeToRemove = 0;
1106  // std::map<int, std::size_t>::iterator rmls = mapLSBSTrkSize.begin();
1107  // SizeToRemove = rmls->second;
1108  // if (debug_)
1109  // edm::LogInfo("BeamMonitor") << " The size to remove is = " << SizeToRemove << endl;
1110  // int changedAfterThis = 0;
1111  // for (std::map<int, std::size_t>::iterator rmLS = mapLSBSTrkSize.begin(); rmLS != mapLSBSTrkSize.end();
1112  // ++rmLS, ++changedAfterThis) {
1113  // if (changedAfterThis > 0) {
1114  // (rmLS->second) = (rmLS->second) - SizeToRemove;
1115  // if ((mapLSBSTrkSize.size() - (size_t)changedAfterThis) == 2)
1116  // PreviousRecords = (rmLS->second);
1117  // }
1118  // }
1119  //
1120  // theBeamFitter->resizeBSvector(SizeToRemove);
1121  //
1122  // map<int, std::size_t>::iterator tmpIt = mapLSBSTrkSize.begin();
1123  // mapLSBSTrkSize.erase(tmpIt);
1124  //
1125  // std::pair<int, int> checkfitLS = theBeamFitter->getFitLSRange();
1126  // std::pair<time_t, time_t> checkfitTime = theBeamFitter->getRefTime();
1127  // theBeamFitter->setFitLSRange(beginLumiOfBSFit_, checkfitLS.second);
1128  // theBeamFitter->setRefTime(refBStime[0], checkfitTime.second);
1129  // }
1130 
1131  //Fill the track for this fit
1132  // vector<BSTrkParameters> theBSvector = theBeamFitter->getBSvector();
1133  // h_nTrk_lumi->ShiftFillLast(theBSvector.size());
1134  //
1135  // if (debug_)
1136  // edm::LogInfo("BeamMonitor") << "FitAndFill:: Size of theBSViector.size() After =" << theBSvector.size() << endl;
1137 
1138  // bool countFitting = false;
1139  // if (theBSvector.size() >= PreviousRecords && theBSvector.size() >= min_Ntrks_) {
1140  // countFitting = true;
1141  // }
1142 
1143  //---Fix for Cut Flow Table for Running average in a same way//the previous code has problem for resetting!!!
1144  // mapLSCF[countLumi_] = *theBeamFitter->getCutFlow();
1145  // if (StartAverage_ && !mapLSCF.empty()) {
1146  // const TH1F& cutFlowToSubtract = mapLSCF.begin()->second;
1147  // // Subtract the last cut flow from all of the others.
1148  // std::map<int, TH1F>::iterator cf = mapLSCF.begin();
1149  // // Start on second entry
1150  // for (; cf != mapLSCF.end(); ++cf) {
1151  // cf->second.Add(&cutFlowToSubtract, -1);
1152  // }
1153  // theBeamFitter->subtractFromCutFlow(&cutFlowToSubtract);
1154  // // Remove the obsolete lumi section
1155  // mapLSCF.erase(mapLSCF.begin());
1156  // }
1157 
1158  if (resetHistos_) {
1159  h_d0_phi0->Reset();
1160  h_vx_vy->Reset();
1161  h_vx_dz->Reset();
1162  h_vy_dz->Reset();
1163  h_trk_z0->Reset();
1164  resetHistos_ = false;
1165  }
1166 
1167  if (StartAverage_)
1168  nthBSTrk_ = PreviousRecords; //after average proccess is ON//for 2-6 LS fit PreviousRecords is size from 2-5 LS
1169 
1170  edm::LogInfo("FakeBeamMonitor") << " The Previous Recored for this fit is =" << nthBSTrk_ << endl;
1171 
1172  // unsigned int itrk = 0;
1173  // for (vector<BSTrkParameters>::const_iterator BSTrk = theBSvector.begin(); BSTrk != theBSvector.end();
1174  // ++BSTrk, ++itrk) {
1175  // if (itrk >= nthBSTrk_) { //fill for this record only !!
1176  // h_d0_phi0->Fill(BSTrk->phi0(), BSTrk->d0());
1177  // double vx = BSTrk->vx();
1178  // double vy = BSTrk->vy();
1179  // double z0 = BSTrk->z0();
1180  // h_vx_vy->Fill(vx, vy);
1181  // h_vx_dz->Fill(z0, vx);
1182  // h_vy_dz->Fill(z0, vy);
1183  // h_trk_z0->Fill(z0);
1184  // }
1185  // }
1186 
1187  // nthBSTrk_ = theBSvector.size(); // keep track of num of tracks filled so far
1188 
1189  edm::LogInfo("FakeBeamMonitor") << " The Current Recored for this fit is =" << nthBSTrk_ << endl;
1190 
1191  // if (countFitting)
1192  // edm::LogInfo("FakeBeamMonitor") << "FitAndFill:: Num of tracks collected = " << nthBSTrk_ << endl;
1193 
1194  if (fitNLumi_ > 0) {
1195  if (onlineMode_) {
1196  if (currentlumi % fitNLumi_ != 0) {
1197  // for (std::map<TString,MonitorElement*>::iterator itAll = hs.begin();
1198  // itAll != hs.end(); ++itAll) {
1199  // if ((*itAll).first.Contains("all")) {
1200  // (*itAll).second->setBinContent(currentlumi,0.);
1201  // (*itAll).second->setBinError(currentlumi,0.);
1202  // }
1203  // }
1204  return;
1205  }
1206  } else if (countLumi_ % fitNLumi_ != 0)
1207  return;
1208  }
1209 
1210  edm::LogInfo("FakeBeamMonitor") << "FitAndFill:: [DebugTime] refBStime[0] = " << refBStime[0]
1211  << "; address = " << &refBStime[0] << std::endl;
1212  edm::LogInfo("FakeBeamMonitor") << "FitAndFill:: [DebugTime] refBStime[1] = " << refBStime[1]
1213  << "; address = " << &refBStime[1] << std::endl;
1214 
1215  //Fill for all LS even if fit fails
1216  // h_nVtx_lumi->ShiftFillLast((theBeamFitter->getPVvectorSize()), 0., fitNLumi_);
1217  // h_nVtx_lumi_all->setBinContent(currentlumi, (theBeamFitter->getPVvectorSize()));
1218 
1219  // if (countFitting) {
1220  nFits_++;
1221  // std::pair<int, int> fitLS = theBeamFitter->getFitLSRange();
1222  std::pair<int, int> fitLS(beginLumiOfBSFit_, endLumiOfBSFit_);
1223  // edm::LogInfo("BeamMonitor") << "FitAndFill:: [BeamFitter] Do BeamSpot Fit for LS = " << fitLS.first << " to "
1224  // << fitLS.second << std::endl;
1225  edm::LogInfo("FakeBeamMonitor") << "FitAndFill:: [FakeBeamMonitor] Do BeamSpot Fit for LS = " << beginLumiOfBSFit_
1226  << " to " << endLumiOfBSFit_ << std::endl;
1227 
1228  //Now Run the PV and Track Fitter over the collected tracks and pvs
1229  // if (theBeamFitter->runPVandTrkFitter()) {
1230  // reco::BeamSpot bs = theBeamFitter->getBeamSpot();
1231 
1232  //Create fake BS here
1234  for (int j = 0; j < 7; ++j) {
1235  for (int k = j; k < 7; ++k) {
1236  matrix(j, k) = 0;
1237  }
1238  }
1239 
1240  // random values for fake BeamSpot
1241  float tmp_BSx = rndm_->Gaus(0.1, 0.1); // [cm]
1242  float tmp_BSy = rndm_->Gaus(0.1, 0.1); // [cm]
1243  float tmp_BSz = rndm_->Gaus(0.1, 0.1); // [cm]
1244  float tmp_BSwidthX = rndm_->Gaus(0.001, 0.0005); // [cm]
1245  float tmp_BSwidthY = rndm_->Gaus(0.001, 0.0005); // [cm]
1246  float tmp_BSwidthZ = rndm_->Gaus(3.5, 0.5); // [cm]
1247 
1248  reco::BeamSpot bs(reco::BeamSpot::Point(tmp_BSx, tmp_BSy, tmp_BSz),
1249  tmp_BSwidthZ,
1250  0,
1251  0,
1252  tmp_BSwidthX,
1253  matrix,
1255  bs.setBeamWidthY(tmp_BSwidthY);
1256 
1257  if (bs.type() > 0) // with good beamwidth fit
1258  preBS = bs; // cache good fit results
1259 
1260  edm::LogInfo("FakeBeamMonitor") << "\n RESULTS OF DEFAULT FIT:" << endl;
1261  edm::LogInfo("FakeBeamMonitor") << bs << endl;
1262  edm::LogInfo("FakeBeamMonitor") << "[BeamFitter] fitting done \n" << endl;
1263 
1264  hs[k_x0_lumi]->ShiftFillLast(bs.x0(), bs.x0Error(), fitNLumi_);
1265  hs[k_y0_lumi]->ShiftFillLast(bs.y0(), bs.y0Error(), fitNLumi_);
1266  hs[k_z0_lumi]->ShiftFillLast(bs.z0(), bs.z0Error(), fitNLumi_);
1267  hs[k_sigmaX0_lumi]->ShiftFillLast(bs.BeamWidthX(), bs.BeamWidthXError(), fitNLumi_);
1268  hs[k_sigmaY0_lumi]->ShiftFillLast(bs.BeamWidthY(), bs.BeamWidthYError(), fitNLumi_);
1269  hs[k_sigmaZ0_lumi]->ShiftFillLast(bs.sigmaZ(), bs.sigmaZ0Error(), fitNLumi_);
1270  hs[k_x0_lumi_all]->setBinContent(currentlumi, bs.x0());
1271  hs[k_x0_lumi_all]->setBinError(currentlumi, bs.x0Error());
1272  hs[k_y0_lumi_all]->setBinContent(currentlumi, bs.y0());
1273  hs[k_y0_lumi_all]->setBinError(currentlumi, bs.y0Error());
1274  hs[k_z0_lumi_all]->setBinContent(currentlumi, bs.z0());
1275  hs[k_z0_lumi_all]->setBinError(currentlumi, bs.z0Error());
1276  hs[k_sigmaX0_lumi_all]->setBinContent(currentlumi, bs.BeamWidthX());
1277  hs[k_sigmaX0_lumi_all]->setBinError(currentlumi, bs.BeamWidthXError());
1278  hs[k_sigmaY0_lumi_all]->setBinContent(currentlumi, bs.BeamWidthY());
1279  hs[k_sigmaY0_lumi_all]->setBinError(currentlumi, bs.BeamWidthYError());
1280  hs[k_sigmaZ0_lumi_all]->setBinContent(currentlumi, bs.sigmaZ());
1281  hs[k_sigmaZ0_lumi_all]->setBinError(currentlumi, bs.sigmaZ0Error());
1282 
1283  int nthBin = tmpTime - refTime;
1284  if (nthBin > 0) {
1285  hs[k_x0_time]->setBinContent(nthBin, bs.x0());
1286  hs[k_y0_time]->setBinContent(nthBin, bs.y0());
1287  hs[k_z0_time]->setBinContent(nthBin, bs.z0());
1288  hs[k_sigmaX0_time]->setBinContent(nthBin, bs.BeamWidthX());
1289  hs[k_sigmaY0_time]->setBinContent(nthBin, bs.BeamWidthY());
1290  hs[k_sigmaZ0_time]->setBinContent(nthBin, bs.sigmaZ());
1291  hs[k_x0_time]->setBinError(nthBin, bs.x0Error());
1292  hs[k_y0_time]->setBinError(nthBin, bs.y0Error());
1293  hs[k_z0_time]->setBinError(nthBin, bs.z0Error());
1294  hs[k_sigmaX0_time]->setBinError(nthBin, bs.BeamWidthXError());
1295  hs[k_sigmaY0_time]->setBinError(nthBin, bs.BeamWidthYError());
1296  hs[k_sigmaZ0_time]->setBinError(nthBin, bs.sigmaZ0Error());
1297  }
1298 
1299  int jthBin = tmpTime - startTime;
1300  if (jthBin > 0) {
1301  hs[k_x0_time_all]->setBinContent(jthBin, bs.x0());
1302  hs[k_y0_time_all]->setBinContent(jthBin, bs.y0());
1303  hs[k_z0_time_all]->setBinContent(jthBin, bs.z0());
1304  hs[k_sigmaX0_time_all]->setBinContent(jthBin, bs.BeamWidthX());
1305  hs[k_sigmaY0_time_all]->setBinContent(jthBin, bs.BeamWidthY());
1306  hs[k_sigmaZ0_time_all]->setBinContent(jthBin, bs.sigmaZ());
1307  hs[k_x0_time_all]->setBinError(jthBin, bs.x0Error());
1308  hs[k_y0_time_all]->setBinError(jthBin, bs.y0Error());
1309  hs[k_z0_time_all]->setBinError(jthBin, bs.z0Error());
1310  hs[k_sigmaX0_time_all]->setBinError(jthBin, bs.BeamWidthXError());
1311  hs[k_sigmaY0_time_all]->setBinError(jthBin, bs.BeamWidthYError());
1312  hs[k_sigmaZ0_time_all]->setBinError(jthBin, bs.sigmaZ0Error());
1313  }
1314 
1315  h_x0->Fill(bs.x0());
1316  h_y0->Fill(bs.y0());
1317  h_z0->Fill(bs.z0());
1318  if (bs.type() > 0) { // with good beamwidth fit
1319  h_sigmaX0->Fill(bs.BeamWidthX());
1320  h_sigmaY0->Fill(bs.BeamWidthY());
1321  }
1322  h_sigmaZ0->Fill(bs.sigmaZ());
1323 
1324  if (nthBSTrk_ >= 2 * min_Ntrks_) {
1325  double amp = std::sqrt(bs.x0() * bs.x0() + bs.y0() * bs.y0());
1326  double alpha = std::atan2(bs.y0(), bs.x0());
1327  std::unique_ptr<TF1> f1{new TF1("f1", "[0]*sin(x-[1])", -3.14, 3.14)};
1328  f1->SetParameters(amp, alpha);
1329  f1->SetParLimits(0, amp - 0.1, amp + 0.1);
1330  f1->SetParLimits(1, alpha - 0.577, alpha + 0.577);
1331  f1->SetLineColor(4);
1332  h_d0_phi0->getTProfile()->Fit(f1.get(), "QR");
1333 
1334  double mean = bs.z0();
1335  double width = bs.sigmaZ();
1336  std::unique_ptr<TF1> fgaus{new TF1("fgaus", "gaus")};
1337  fgaus->SetParameters(mean, width);
1338  fgaus->SetLineColor(4);
1339  h_trk_z0->getTH1()->Fit(fgaus.get(), "QLRM", "", mean - 3 * width, mean + 3 * width);
1340  }
1341 
1342  fitResults->Reset();
1343  std::pair<int, int> LSRange(beginLumiOfBSFit_, endLumiOfBSFit_); //= theBeamFitter->getFitLSRange();
1344  char tmpTitle[50];
1345  sprintf(tmpTitle, "%s %i %s %i", "Fitted Beam Spot (cm) of LS: ", LSRange.first, " to ", LSRange.second);
1346  fitResults->setAxisTitle(tmpTitle, 1);
1347  fitResults->setBinContent(1, 8, bs.x0());
1348  fitResults->setBinContent(1, 7, bs.y0());
1349  fitResults->setBinContent(1, 6, bs.z0());
1350  fitResults->setBinContent(1, 5, bs.sigmaZ());
1351  fitResults->setBinContent(1, 4, bs.dxdz());
1352  fitResults->setBinContent(1, 3, bs.dydz());
1353  if (bs.type() > 0) { // with good beamwidth fit
1354  fitResults->setBinContent(1, 2, bs.BeamWidthX());
1355  fitResults->setBinContent(1, 1, bs.BeamWidthY());
1356  } else { // fill cached widths
1359  }
1360 
1361  fitResults->setBinContent(2, 8, bs.x0Error());
1362  fitResults->setBinContent(2, 7, bs.y0Error());
1363  fitResults->setBinContent(2, 6, bs.z0Error());
1364  fitResults->setBinContent(2, 5, bs.sigmaZ0Error());
1365  fitResults->setBinContent(2, 4, bs.dxdzError());
1366  fitResults->setBinContent(2, 3, bs.dydzError());
1367  if (bs.type() > 0) { // with good beamwidth fit
1370  } else { // fill cached width errors
1373  }
1374 
1375  // count good fit
1376  // if (std::fabs(refBS.x0()-bs.x0())/bs.x0Error() < deltaSigCut_) { // disabled temporarily
1377  summaryContent_[0] += 1.;
1378  // }
1379  // if (std::fabs(refBS.y0()-bs.y0())/bs.y0Error() < deltaSigCut_) { // disabled temporarily
1380  summaryContent_[1] += 1.;
1381  // }
1382  // if (std::fabs(refBS.z0()-bs.z0())/bs.z0Error() < deltaSigCut_) { // disabled temporarily
1383  summaryContent_[2] += 1.;
1384  // }
1385 
1386  // Create the BeamSpotOnlineObjects object
1387  BeamSpotOnlineObjects BSOnline;
1388  BSOnline.setLastAnalyzedLumi(LSRange.second);
1389  BSOnline.setLastAnalyzedRun(frun);
1390  BSOnline.setLastAnalyzedFill(0); // To be updated with correct LHC Fill number
1391  BSOnline.setPosition(bs.x0(), bs.y0(), bs.z0());
1392  BSOnline.setSigmaZ(bs.sigmaZ());
1393  BSOnline.setBeamWidthX(bs.BeamWidthX());
1394  BSOnline.setBeamWidthY(bs.BeamWidthY());
1395  BSOnline.setBeamWidthXError(bs.BeamWidthXError());
1396  BSOnline.setBeamWidthYError(bs.BeamWidthYError());
1397  BSOnline.setdxdz(bs.dxdz());
1398  BSOnline.setdydz(bs.dydz());
1399  BSOnline.setType(bs.type());
1400  BSOnline.setEmittanceX(bs.emittanceX());
1401  BSOnline.setEmittanceY(bs.emittanceY());
1402  BSOnline.setBetaStar(bs.betaStar());
1403  for (int i = 0; i < 7; ++i) {
1404  for (int j = 0; j < 7; ++j) {
1405  BSOnline.setCovariance(i, j, bs.covariance(i, j));
1406  }
1407  }
1408  BSOnline.setNumTracks(50);
1409  BSOnline.setNumPVs(10);
1410  BSOnline.setUsedEvents((int)DipPVInfo_[0]);
1411  BSOnline.setMeanPV(DipPVInfo_[1]);
1412  BSOnline.setMeanErrorPV(DipPVInfo_[2]);
1413  BSOnline.setRmsPV(DipPVInfo_[3]);
1414  BSOnline.setRmsErrorPV(DipPVInfo_[4]);
1415  BSOnline.setMaxPVs((int)DipPVInfo_[5]);
1416  auto creationTime =
1417  std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
1418  BSOnline.setCreationTime(creationTime);
1419 
1420  // use fake timestamps from 1970.01.01 00:00:00 to 1970.01.01 00:00:01 GMT
1421  std::pair<time_t, time_t> timeForDIP = std::make_pair(0, 1);
1422  BSOnline.setStartTimeStamp(timeForDIP.first);
1423  BSOnline.setStartTime(getGMTstring(timeForDIP.first));
1424  BSOnline.setEndTimeStamp(timeForDIP.second);
1425  BSOnline.setEndTime(getGMTstring(timeForDIP.second));
1426 
1427  edm::LogInfo("FakeBeamMonitor") << "FitAndFill::[PayloadCreation] BeamSpotOnline object created: \n" << std::endl;
1428  edm::LogInfo("FakeBeamMonitor") << BSOnline << std::endl;
1429 
1430  // Create the payload for BeamSpotOnlineObjects object
1431  if (onlineDbService_.isAvailable()) {
1432  edm::LogInfo("FakeBeamMonitor") << "FitAndFill::[PayloadCreation] onlineDbService available \n" << std::endl;
1433  onlineDbService_->logger().logInfo() << "FakeBeamMonitor::FitAndFill - Lumi of the current fit: " << currentlumi;
1434  onlineDbService_->logger().logInfo() << "FakeBeamMonitor::FitAndFill - Do PV Fitting for LS = " << beginLumiOfPVFit_
1435  << " to " << endLumiOfPVFit_;
1436  onlineDbService_->logger().logInfo() << "FakeBeamMonitor::FitAndFill - [BeamFitter] Do BeamSpot Fit for LS = "
1437  << LSRange.first << " to " << LSRange.second;
1438  onlineDbService_->logger().logInfo() << "FakeBeamMonitor::FitAndFill - [FakeBeamMonitor] Do BeamSpot Fit for LS = "
1439  << beginLumiOfBSFit_ << " to " << endLumiOfBSFit_;
1440  onlineDbService_->logger().logInfo() << "FakeBeamMonitor - RESULTS OF DEFAULT FIT:";
1441  onlineDbService_->logger().logInfo() << "\n" << bs;
1443  << "FakeBeamMonitor::FitAndFill - [PayloadCreation] BeamSpotOnline object created:";
1444  onlineDbService_->logger().logInfo() << "\n" << BSOnline;
1445  onlineDbService_->logger().logInfo() << "FakeBeamMonitor - Additional parameters for DIP:";
1446  onlineDbService_->logger().logInfo() << "Events used in the fit: " << BSOnline.usedEvents();
1447  onlineDbService_->logger().logInfo() << "Mean PV : " << BSOnline.meanPV();
1448  onlineDbService_->logger().logInfo() << "Mean PV Error : " << BSOnline.meanErrorPV();
1449  onlineDbService_->logger().logInfo() << "Rms PV : " << BSOnline.rmsPV();
1450  onlineDbService_->logger().logInfo() << "Rms PV Error : " << BSOnline.rmsErrorPV();
1451  onlineDbService_->logger().logInfo() << "Max PVs : " << BSOnline.maxPVs();
1452  onlineDbService_->logger().logInfo() << "StartTime : " << BSOnline.startTime();
1453  onlineDbService_->logger().logInfo() << "StartTimeStamp : " << BSOnline.startTimeStamp();
1454  onlineDbService_->logger().logInfo() << "EndTime : " << BSOnline.endTime();
1455  onlineDbService_->logger().logInfo() << "EndTimeStamp : " << BSOnline.endTimeStamp();
1456  onlineDbService_->logger().logInfo() << "FakeBeamMonitor::FitAndFill - [PayloadCreation] onlineDbService available";
1457  onlineDbService_->logger().logInfo() << "FakeBeamMonitor::FitAndFill - [PayloadCreation] SetCreationTime: "
1458  << creationTime << " [epoch in microseconds]";
1459  try {
1462  << "FakeBeamMonitor::FitAndFill - [PayloadCreation] writeIOVForNextLumisection executed correctly";
1463  } catch (const std::exception& e) {
1464  onlineDbService_->logger().logError() << "FakeBeamMonitor - Error writing record: " << recordName_
1465  << " for Run: " << frun << " - Lumi: " << LSRange.second;
1466  onlineDbService_->logger().logError() << "Error is: " << e.what();
1467  onlineDbService_->logger().logError() << "RESULTS OF DEFAULT FIT WAS:";
1468  onlineDbService_->logger().logError() << "\n" << bs;
1469  DBloggerReturn_ = 2;
1470  }
1471  }
1472  edm::LogInfo("FakeBeamMonitor") << "FitAndFill::[PayloadCreation] BeamSpotOnline payload created \n" << std::endl;
1473 
1474  // } //if (theBeamFitter->runPVandTrkFitter())
1475  // else { // beam fit fails
1476  // reco::BeamSpot bs = theBeamFitter->getBeamSpot();
1477  // edm::LogInfo("BeamMonitor") << "FitAndFill:: [BeamMonitor] Beam fit fails!!! \n" << endl;
1478  // edm::LogInfo("BeamMonitor") << "FitAndFill:: [BeamMonitor] Output beam spot for DIP \n" << endl;
1479  // edm::LogInfo("BeamMonitor") << bs << endl;
1480  //
1481  // hs[k_sigmaX0_lumi]->ShiftFillLast(bs.BeamWidthX(), bs.BeamWidthXError(), fitNLumi_);
1482  // hs[k_sigmaY0_lumi]->ShiftFillLast(bs.BeamWidthY(), bs.BeamWidthYError(), fitNLumi_);
1483  // hs[k_sigmaZ0_lumi]->ShiftFillLast(bs.sigmaZ(), bs.sigmaZ0Error(), fitNLumi_);
1484  // hs[k_x0_lumi]->ShiftFillLast(bs.x0(), bs.x0Error(), fitNLumi_);
1485  // hs[k_y0_lumi]->ShiftFillLast(bs.y0(), bs.y0Error(), fitNLumi_);
1486  // hs[k_z0_lumi]->ShiftFillLast(bs.z0(), bs.z0Error(), fitNLumi_);
1487  // } // end of beam fit fails
1488 
1489  // } //-------- end of countFitting------------------------------------------
1490  // else { // no fit
1491  // // Overwrite Fit LS and fit time when no event processed or no track selected
1492  // theBeamFitter->setFitLSRange(beginLumiOfBSFit_, endLumiOfBSFit_);
1493  // theBeamFitter->setRefTime(refBStime[0], refBStime[1]);
1494  // if (theBeamFitter->runPVandTrkFitter()) {
1495  // } // Dump fake beam spot for DIP
1496  // reco::BeamSpot bs = theBeamFitter->getBeamSpot();
1497  // edm::LogInfo("BeamMonitor") << "FitAndFill:: [BeamMonitor] No fitting \n" << endl;
1498  // edm::LogInfo("BeamMonitor") << "FitAndFill:: [BeamMonitor] Output fake beam spot for DIP \n" << endl;
1499  // edm::LogInfo("BeamMonitor") << bs << endl;
1500  //
1501  // hs[k_sigmaX0_lumi]->ShiftFillLast(bs.BeamWidthX(), bs.BeamWidthXError(), fitNLumi_);
1502  // hs[k_sigmaY0_lumi]->ShiftFillLast(bs.BeamWidthY(), bs.BeamWidthYError(), fitNLumi_);
1503  // hs[k_sigmaZ0_lumi]->ShiftFillLast(bs.sigmaZ(), bs.sigmaZ0Error(), fitNLumi_);
1504  // hs[k_x0_lumi]->ShiftFillLast(bs.x0(), bs.x0Error(), fitNLumi_);
1505  // hs[k_y0_lumi]->ShiftFillLast(bs.y0(), bs.y0Error(), fitNLumi_);
1506  // hs[k_z0_lumi]->ShiftFillLast(bs.z0(), bs.z0Error(), fitNLumi_);
1507  // }
1508 
1509  // Fill summary report
1510  // if (countFitting) {
1511  for (int n = 0; n < nFitElements_; n++) {
1513  }
1514 
1515  summarySum_ = 0;
1516  for (int ii = 0; ii < nFitElements_; ii++) {
1518  }
1519  reportSummary_ = summarySum_ / (nFitElements_ * nFits_);
1520  if (reportSummary)
1522 
1523  for (int bi = 0; bi < nFitElements_; bi++) {
1524  reportSummaryMap->setBinContent(1, bi + 1, summaryContent_[bi] / (float)nFits_);
1525  }
1526  // }
1527 
1528  if ((resetFitNLumi_ > 0 &&
1529  ((onlineMode_ &&
1530  countLumi_ == resetFitNLumi_) || //OR it should be currentLumi_ (if in sequence then does not mattar)
1531  (!onlineMode_ && countLumi_ == resetFitNLumi_))) ||
1532  (StartAverage_)) {
1533  edm::LogInfo("FakeBeamMonitor") << "FitAndFill:: The flag is ON for running average Beam Spot fit" << endl;
1534  StartAverage_ = true;
1535  firstAverageFit_++;
1536  resetHistos_ = true;
1537  nthBSTrk_ = 0;
1538  beginLumiOfBSFit_ = 0;
1539  refBStime[0] = 0;
1540  }
1541 }
1542 
1543 //--------------------------------------------------------
1545  if (debug_)
1546  edm::LogInfo("FakeBeamMonitor")
1547  << " RestartingFitting:: Restart Beami everything to a fresh start !!! because Gap is > 10 LS" << endl;
1548  //track based fit reset here
1549  resetHistos_ = true;
1550  nthBSTrk_ = 0;
1551  // theBeamFitter->resetTrkVector();
1552  // theBeamFitter->resetLSRange();
1553  // theBeamFitter->resetRefTime();
1554  // theBeamFitter->resetPVFitter();
1555  // theBeamFitter->resetCutFlow();
1556  beginLumiOfBSFit_ = 0;
1557  refBStime[0] = 0;
1558  //pv based fit iis reset here
1559  h_PVx[0]->Reset();
1560  h_PVy[0]->Reset();
1561  h_PVz[0]->Reset();
1562  beginLumiOfPVFit_ = 0;
1563  refPVtime[0] = 0;
1564  //Clear all the Maps here
1565  mapPVx.clear();
1566  mapPVy.clear();
1567  mapPVz.clear();
1568  mapNPV.clear();
1569  mapBeginBSLS.clear();
1570  mapBeginPVLS.clear();
1571  mapBeginBSTime.clear();
1572  mapBeginPVTime.clear();
1573  mapLSBSTrkSize.clear();
1574  mapLSPVStoreSize.clear();
1575  mapLSCF.clear();
1576  countGapLumi_ = 0;
1577  countLumi_ = 0;
1578  StartAverage_ = false;
1579 }
1580 
1581 //-------------------------------------------------------
1583  if (debug_)
1584  edm::LogInfo("FakeBeamMonitor") << "dqmEndRun:: Clearing all the Maps " << endl;
1585  //Clear all the Maps here
1586  mapPVx.clear();
1587  mapPVy.clear();
1588  mapPVz.clear();
1589  mapNPV.clear();
1590  mapBeginBSLS.clear();
1591  mapBeginPVLS.clear();
1592  mapBeginBSTime.clear();
1593  mapBeginPVTime.clear();
1594  mapLSBSTrkSize.clear();
1595  mapLSPVStoreSize.clear();
1596  mapLSCF.clear();
1599  }
1600 }
1601 
1602 //--------------------------------------------------------
1603 void FakeBeamMonitor::scrollTH1(TH1* h, time_t ref) {
1604  char offsetTime[64];
1605  formatFitTime(offsetTime, ref);
1606  TDatime da(offsetTime);
1607  if (lastNZbin > 0) {
1608  double val = h->GetBinContent(lastNZbin);
1609  double valErr = h->GetBinError(lastNZbin);
1610  h->Reset();
1611  h->GetXaxis()->SetTimeOffset(da.Convert(kTRUE));
1612  int bin = (lastNZbin > buffTime ? buffTime : 1);
1613  h->SetBinContent(bin, val);
1614  h->SetBinError(bin, valErr);
1615  } else {
1616  h->Reset();
1617  h->GetXaxis()->SetTimeOffset(da.Convert(kTRUE));
1618  }
1619 }
1620 
1621 //--------------------------------------------------------
1622 // Method to check whether to chane histogram time offset (forward only)
1623 bool FakeBeamMonitor::testScroll(time_t& tmpTime_, time_t& refTime_) {
1624  bool scroll_ = false;
1625  if (tmpTime_ - refTime_ >= intervalInSec_) {
1626  scroll_ = true;
1627  edm::LogInfo("FakeBeamMonitor") << "testScroll:: Reset Time Offset" << std::endl;
1629  for (int bin = intervalInSec_; bin >= 1; bin--) {
1630  if (hs[k_x0_time]->getBinContent(bin) > 0) {
1631  lastNZbin = bin;
1632  break;
1633  }
1634  }
1635  edm::LogInfo("FakeBeamMonitor") << "testScroll:: Last non zero bin = " << lastNZbin << std::endl;
1636  if (tmpTime_ - refTime_ >= intervalInSec_ + lastNZbin) {
1637  edm::LogInfo("FakeBeamMonitor") << "testScroll:: Time difference too large since last readout" << std::endl;
1638  lastNZbin = 0;
1639  refTime_ = tmpTime_ - buffTime;
1640  } else {
1641  edm::LogInfo("FakeBeamMonitor") << "testScroll:: Offset to last record" << std::endl;
1642  int offset = ((lastNZbin > buffTime) ? (lastNZbin - buffTime) : (lastNZbin - 1));
1643  refTime_ += offset;
1644  }
1645  }
1646  return scroll_;
1647 }
1648 
1650 
1651 // Local Variables:
1652 // show-trailing-whitespace: t
1653 // truncate-lines: t
1654 // End:
MonitorElement * pvResults
LuminosityBlockID id() const
#define CEST
math::Error< dimension >::type CovarianceMatrix
Definition: BeamSpot.h:29
void setBeamWidthYError(double val)
set beam width Y error
double z0() const
z coordinate
Definition: BeamSpot.h:65
T getUntrackedParameter(std::string const &, T const &) const
float alpha
Definition: AMPTWrapper.h:105
void setCovariance(int i, int j, double val)
set i,j element of the full covariance matrix 7x7
double sigmaZ0Error() const
error on sigma z
Definition: BeamSpot.h:92
MonitorElement * bookFloat(TString const &name, FUNC onbooking=NOOP())
Definition: DQMStore.h:80
std::map< int, size_t > mapLSPVStoreSize
void setStartTimeStamp(cond::Time_t val)
MonitorElement * h_z0
std::vector< MonitorElement * > hs
const double phiMax_
void setLastAnalyzedLumi(int val)
Setters Methods.
RunNumber_t run() const
Definition: RunBase.h:40
MonitorElement * h_nTrk_lumi
void setEmittanceY(double val)
set emittance
virtual void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:32
MonitorElement * h_sigmaY0
const double vxMax_
void dqmBeginRun(edm::Run const &, edm::EventSetup const &) override
double dydzError() const
error on dydz
Definition: BeamSpot.h:96
void FitAndFill(const edm::LuminosityBlock &lumiSeg, int &, int &, int &)
MonitorElement * cutFlowTable
virtual void setOption(const char *option)
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
MonitorElement * h_vy_dz
const double phiMin_
std::string getGMTstring(const std::time_t &)
const double dxMin_
void dqmEndRun(const edm::Run &r, const edm::EventSetup &c) override
Hists
Definition: BeamMonitor.cc:176
std::map< int, std::time_t > mapBeginPVTime
void setBetaStar(double val)
set beta star
unsigned int nthBSTrk_
FakeBeamMonitor(const edm::ParameterSet &)
bool testScroll(std::time_t &, std::time_t &)
edm::LuminosityBlockNumber_t luminosityBlock() const
Definition: EventBase.h:61
std::map< int, std::size_t > mapLSBSTrkSize
void setEndTimeStamp(cond::Time_t val)
math::XYZPoint Point
point in the space
Definition: BeamSpot.h:27
MonitorElement * reportSummary
MonitorElement * h_nVtx_lumi_all
Timestamp const & beginTime() const
Float_t summaryContent_[3]
Log< level::Error, false > LogError
std::time_t refBStime[2]
int ii
Definition: cuy.py:589
assert(be >=bs)
std::map< int, std::vector< int > > mapNPV
const double vxMin_
void beginLuminosityBlock(const edm::LuminosityBlock &lumiSeg, const edm::EventSetup &context) override
void setType(int type)
set beam type
MonitorElement * h_PVz[2]
void setEndTime(std::string val)
double emittanceX() const
additional information
Definition: BeamSpot.h:129
MonitorElement * h_trkPt
unsigned int minNrVertices_
LuminosityBlockNumber_t luminosityBlock() const
edm::Service< cond::service::OnlineDBOutputService > onlineDbService_
void Fill(long long x)
virtual void Reset()
Remove all data from the ME, keept the empty histogram with all its settings.
DQM_DEPRECATED void ShiftFillLast(double y, double ye=0., int32_t xscale=1)
double betaStar() const
Definition: BeamSpot.h:131
char const * label
double dydz() const
dydz slope
Definition: BeamSpot.h:80
void setEmittanceX(double val)
set emittance
int iEvent
Definition: GenABIO.cc:224
cond::Time_t startTimeStamp() const
void setBeamWidthY(double v)
Definition: BeamSpot.h:105
double emittanceY() const
Definition: BeamSpot.h:130
EchoedLogStream< edm::LogError > logError()
Definition: Logger.cc:165
double dxdzError() const
error on dxdz
Definition: BeamSpot.h:94
std::time_t refPVtime[2]
MonitorElement * h_nVtx
MonitorElement * h_trkVz
MonitorElement * bookProfile(TString const &name, TString const &title, int nchX, double lowX, double highX, int, double lowY, double highY, char const *option="s", FUNC onbooking=NOOP())
Definition: DQMStore.h:322
cond::persistency::Logger & logger()
T sqrt(T t)
Definition: SSEVec.h:19
void formatFitTime(char *, const std::time_t &)
const double dzMin_
std::string startTime() const
Timestamp const & endTime() const
std::time_t startTime
MonitorElement * h_sigmaZ0
void setdydz(double val)
set dydz slope, crossing angle in XZ
bool isAvailable() const
Definition: Service.h:40
std::string monitorName_
MonitorElement * h_x0
std::map< int, std::vector< float > > mapPVy
std::map< int, std::vector< float > > mapPVz
void bookHistograms(DQMStore::IBooker &i, const edm::Run &r, const edm::EventSetup &c) override
MonitorElement * h_d0_phi0
double BeamWidthX() const
beam width X
Definition: BeamSpot.h:82
cond::Time_t endTimeStamp() const
double BeamWidthYError() const
error on beam width Y, assume error in X = Y
Definition: BeamSpot.h:101
cond::Time_t writeIOVForNextLumisection(const PayloadType &payload, const std::string &recordName)
virtual void setBinLabel(int bin, const std::string &label, int axis=1)
set bin label for x, y or z axis (axis=1, 2, 3 respectively)
double BeamWidthXError() const
error on beam width X, assume error in X = Y
Definition: BeamSpot.h:99
const double dxMax_
double z0Error() const
error on z
Definition: BeamSpot.h:90
unsigned long long TimeValue_t
Definition: Timestamp.h:28
std::map< int, std::vector< float > > mapPVx
virtual TProfile * getTProfile()
void setdxdz(double val)
set dxdz slope, crossing angle
double dxdz() const
dxdz slope
Definition: BeamSpot.h:78
std::map< int, std::time_t > mapBeginBSTime
double x0Error() const
error on x
Definition: BeamSpot.h:86
double y0Error() const
error on y
Definition: BeamSpot.h:88
MonitorElement * h_sigmaX0
edm::TimeValue_t ftimestamp
Log< level::Info, false > LogInfo
std::map< int, int > mapBeginPVLS
void scrollTH1(TH1 *, std::time_t)
~FakeBeamMonitor() override
virtual void setBinContent(int binx, double content)
set content of bin (1-D)
void setBeamWidthXError(double val)
set beam width X error
MonitorElement * h_trk_z0
std::time_t refTime
void setSigmaZ(double val)
set sigma Z, RMS bunch length
double sigmaZ() const
sigma z
Definition: BeamSpot.h:76
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
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
MonitorElement * reportSummaryContents[3]
double BeamWidthY() const
beam width Y
Definition: BeamSpot.h:84
double b
Definition: hdecay.h:118
LuminosityBlockNumber_t luminosityBlock() const
MonitorElement * h_nVtx_lumi
MonitorElement * h_vx_vy
MonitorElement * h_vx_dz
std::map< int, TH1F > mapLSCF
Timestamp const & beginTime() const
Definition: RunBase.h:41
const double dzMax_
double covariance(int i, int j) const
(i,j)-th element of error matrix
Definition: BeamSpot.h:108
void endLuminosityBlock(const edm::LuminosityBlock &lumiSeg, const edm::EventSetup &c) override
MonitorElement * h_PVy[2]
std::string recordName_
MonitorElement * h_y0
MonitorElement * h_PVxz
reco::BeamSpot preBS
MonitorElement * h_nVtx_st
void setStartTime(std::string val)
void setPosition(double x, double y, double z)
set XYZ position
string end
Definition: dataset.py:937
double y0() const
y coordinate
Definition: BeamSpot.h:63
void setBeamWidthX(double val)
set average transverse beam width X
std::time_t tmpTime
void setCreationTime(cond::Time_t val)
std::map< int, int > mapBeginBSLS
void analyze(const edm::Event &e, const edm::EventSetup &c) override
MonitorElement * reportSummaryMap
MonitorElement * h_PVyz
MonitorElement * h_PVx[2]
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98
EchoedLogStream< edm::LogInfo > logInfo()
Definition: Logger.cc:157
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
void end(int retCode)
Definition: Logger.cc:80
unsigned int min_Ntrks_
std::string endTime() const
TimeValue_t value() const
Definition: Timestamp.h:45
static constexpr int buffTime
void setBeamWidthY(double val)
set average transverse beam width Y
MonitorElement * fitResults
Definition: Run.h:45
BeamType type() const
return beam type
Definition: BeamSpot.h:122
double x0() const
x coordinate
Definition: BeamSpot.h:61
virtual void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
Run const & getRun() const