CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
BeamMonitor.cc
Go to the documentation of this file.
1 /*
2  * \file BeamMonitor.cc
3  * \author Geng-yuan Jeng/UC Riverside
4  * Francisco Yumiceva/FNAL
5  * $Date: 2010/10/28 08:46:01 $
6  * $Revision: 1.58 $
7  *
8  */
9 
26 #include <numeric>
27 #include <math.h>
28 #include <TMath.h>
29 #include <iostream>
30 #include <TStyle.h>
31 
32 using namespace std;
33 using namespace edm;
34 
35 const char * BeamMonitor::formatFitTime( const time_t & t ) {
36 #define CET (+1)
37 #define CEST (+2)
38 
39  static char ts[] = "yyyy-Mm-dd hh:mm:ss";
40  tm * ptm;
41  ptm = gmtime ( &t );
42  int year = ptm->tm_year;
43  if (year < 1995) {
44  edm::LogError("BadTimeStamp") << "year reported is " << year << "!! resetting to 2010..." << std::endl;
45  year = 2010;
46  }
47  sprintf( ts, "%4d-%02d-%02d %02d:%02d:%02d", 2010,ptm->tm_mon+1,ptm->tm_mday,(ptm->tm_hour+CEST)%24, ptm->tm_min, ptm->tm_sec);
48 
49 #ifdef STRIP_TRAILING_BLANKS_IN_TIMEZONE
50  unsigned int b = strlen(ts);
51  while (ts[--b] == ' ') {ts[b] = 0;}
52 #endif
53  return ts;
54 
55 }
56 
57 #define buffTime (23)
58 
59 //
60 // constructors and destructor
61 //
63  countEvt_(0),countLumi_(0),nthBSTrk_(0),nFitElements_(3),resetHistos_(false) {
64 
65  parameters_ = ps;
66  monitorName_ = parameters_.getUntrackedParameter<string>("monitorName","YourSubsystemName");
69  intervalInSec_ = parameters_.getUntrackedParameter<int>("timeInterval",920);//40 LS X 23"
70  fitNLumi_ = parameters_.getUntrackedParameter<int>("fitEveryNLumi",-1);
71  resetFitNLumi_ = parameters_.getUntrackedParameter<int>("resetEveryNLumi",-1);
72  fitPVNLumi_ = parameters_.getUntrackedParameter<int>("fitPVEveryNLumi",-1);
73  resetPVNLumi_ = parameters_.getUntrackedParameter<int>("resetPVEveryNLumi",-1);
74  deltaSigCut_ = parameters_.getUntrackedParameter<double>("deltaSignificanceCut",15);
75  debug_ = parameters_.getUntrackedParameter<bool>("Debug");
76  onlineMode_ = parameters_.getUntrackedParameter<bool>("OnlineMode");
77  tracksLabel_ = parameters_.getParameter<ParameterSet>("BeamFitter").getUntrackedParameter<InputTag>("TrackCollection");
78  min_Ntrks_ = parameters_.getParameter<ParameterSet>("BeamFitter").getUntrackedParameter<int>("MinimumInputTracks");
79  maxZ_ = parameters_.getParameter<ParameterSet>("BeamFitter").getUntrackedParameter<double>("MaximumZ");
80  minNrVertices_ = parameters_.getParameter<ParameterSet>("PVFitter").getUntrackedParameter<unsigned int>("minNrVerticesForFit");
81  minVtxNdf_ = parameters_.getParameter<ParameterSet>("PVFitter").getUntrackedParameter<double>("minVertexNdf");
82  minVtxWgt_ = parameters_.getParameter<ParameterSet>("PVFitter").getUntrackedParameter<double>("minVertexMeanWeight");
83 
85 
86  if (monitorName_ != "" ) monitorName_ = monitorName_+"/" ;
87 
93 
94  if (fitNLumi_ <= 0) fitNLumi_ = 1;
96  refBStime[0] = refBStime[1] = refPVtime[0] = refPVtime[1] = 0;
97  maxZ_ = std::fabs(maxZ_);
98  lastlumi_ = 0;
99  nextlumi_ = 0;
100  processed_ = false;
101 }
102 
103 
105  delete theBeamFitter;
106 }
107 
108 
109 //--------------------------------------------------------
111 
112  // book some histograms here
113  const int dxBin = parameters_.getParameter<int>("dxBin");
114  const double dxMin = parameters_.getParameter<double>("dxMin");
115  const double dxMax = parameters_.getParameter<double>("dxMax");
116 
117  const int vxBin = parameters_.getParameter<int>("vxBin");
118  const double vxMin = parameters_.getParameter<double>("vxMin");
119  const double vxMax = parameters_.getParameter<double>("vxMax");
120 
121  const int phiBin = parameters_.getParameter<int>("phiBin");
122  const double phiMin = parameters_.getParameter<double>("phiMin");
123  const double phiMax = parameters_.getParameter<double>("phiMax");
124 
125  const int dzBin = parameters_.getParameter<int>("dzBin");
126  const double dzMin = parameters_.getParameter<double>("dzMin");
127  const double dzMax = parameters_.getParameter<double>("dzMax");
128 
129  // create and cd into new folder
131 
132  h_nTrk_lumi=dbe_->book1D("nTrk_lumi","Num. of selected tracks vs lumi",20,0.5,20.5);
133  h_nTrk_lumi->setAxisTitle("Lumisection",1);
134  h_nTrk_lumi->setAxisTitle("Num of Tracks",2);
135 
136  h_d0_phi0 = dbe_->bookProfile("d0_phi0","d_{0} vs. #phi_{0} (Selected Tracks)",phiBin,phiMin,phiMax,dxBin,dxMin,dxMax,"");
137  h_d0_phi0->setAxisTitle("#phi_{0} (rad)",1);
138  h_d0_phi0->setAxisTitle("d_{0} (cm)",2);
139 
140  h_vx_vy = dbe_->book2D("trk_vx_vy","Vertex (PCA) position of selected tracks",vxBin,vxMin,vxMax,vxBin,vxMin,vxMax);
141  h_vx_vy->getTH2F()->SetOption("COLZ");
142  // h_vx_vy->getTH1()->SetBit(TH1::kCanRebin);
143  h_vx_vy->setAxisTitle("x coordinate of input track at PCA (cm)",1);
144  h_vx_vy->setAxisTitle("y coordinate of input track at PCA (cm)",2);
145 
146  TDatime *da = new TDatime();
147  gStyle->SetTimeOffset(da->Convert(kTRUE));
148 
149  const int nvar_ = 6;
150  string coord[nvar_] = {"x","y","z","sigmaX","sigmaY","sigmaZ"};
151  string label[nvar_] = {"x_{0} (cm)","y_{0} (cm)","z_{0} (cm)",
152  "#sigma_{X_{0}} (cm)","#sigma_{Y_{0}} (cm)","#sigma_{Z_{0}} (cm)"};
153  for (int i = 0; i < 4; i++) {
155  for (int ic=0; ic<nvar_; ++ic) {
156  TString histName(coord[ic]);
157  TString histTitle(coord[ic]);
158  string ytitle(label[ic]);
159  string xtitle("");
160  string options("E1");
161  bool createHisto = true;
162  switch(i) {
163  case 1: // BS vs time
164  histName += "0_time";
165  xtitle = "Time [UTC]";
166  if (ic < 3)
167  histTitle += " coordinate of beam spot vs time (Fit)";
168  else
169  histTitle = histTitle.Insert(5," ") + " of beam spot vs time (Fit)";
170  break;
171  case 2: // PV vs lumi
172  if (ic < 3) {
173  dbe_->setCurrentFolder(monitorName_+"PrimaryVertex");
174  histName.Insert(0,"PV");
175  histName += "_lumi";
176  histTitle.Insert(0,"Avg. ");
177  histTitle += " position of primary vtx vs lumi";
178  xtitle = "Lumisection";
179  ytitle.insert(0,"PV");
180  ytitle += " #pm #sigma_{PV";
181  ytitle += coord[ic];
182  ytitle += "} (cm)";
183  }
184  else createHisto = false;
185  break;
186  case 3: // PV vs time
187  if (ic < 3) {
188  dbe_->setCurrentFolder(monitorName_+"PrimaryVertex");
189  histName.Insert(0,"PV");
190  histName += "_time";
191  histTitle.Insert(0,"Avg. ");
192  histTitle += " position of primary vtx vs time";
193  xtitle = "Time [UTC]";
194  ytitle.insert(0,"PV");
195  ytitle += " #pm #sigma_{PV";
196  ytitle += coord[ic];
197  ytitle += "} (cm)";
198  }
199  else createHisto = false;
200  break;
201  default: // BS vs lumi
202  histName += "0_lumi";
203  xtitle = "Lumisection";
204  if (ic < 3)
205  histTitle += " coordinate of beam spot vs lumi (Fit)";
206  else
207  histTitle = histTitle.Insert(5," ") + " of beam spot vs lumi (Fit)";
208  break;
209  }
210  if (createHisto) {
211  edm::LogInfo("BeamMonitor") << "hitsName = " << histName << "; histTitle = " << histTitle << std::endl;
212  hs[histName] = dbe_->book1D(histName,histTitle,40,0.5,40.5);
213  hs[histName]->setAxisTitle(xtitle,1);
214  hs[histName]->setAxisTitle(ytitle,2);
215  hs[histName]->getTH1()->SetOption("E1");
216  if (histName.Contains("time")) {
217  //int nbins = (intervalInSec_/23 > 0 ? intervalInSec_/23 : 40);
218  hs[histName]->getTH1()->SetBins(intervalInSec_,0.5,intervalInSec_+0.5);
219  hs[histName]->setAxisTimeDisplay(1);
220  hs[histName]->setAxisTimeFormat("%H:%M:%S",1);
221  }
222  histName += "_all";
223  histTitle += " all";
224  hs[histName] = dbe_->book1D(histName,histTitle,40,0.5,40.5);
225  hs[histName]->getTH1()->SetBit(TH1::kCanRebin);
226  hs[histName]->setAxisTitle(xtitle,1);
227  hs[histName]->setAxisTitle(ytitle,2);
228  hs[histName]->getTH1()->SetOption("E1");
229  if (histName.Contains("time")) {
230  //int nbins = (intervalInSec_/23 > 0 ? intervalInSec_/23 : 40);
231  hs[histName]->getTH1()->SetBins(intervalInSec_,0.5,intervalInSec_+0.5);
232  hs[histName]->setAxisTimeDisplay(1);
233  hs[histName]->setAxisTimeFormat("%H:%M:%S",1);
234  }
235  }
236  }
237  }
239 
240  h_trk_z0 = dbe_->book1D("trk_z0","z_{0} of selected tracks",dzBin,dzMin,dzMax);
241  h_trk_z0->setAxisTitle("z_{0} of selected tracks (cm)",1);
242 
243  h_vx_dz = dbe_->bookProfile("vx_dz","v_{x} vs. dz of selected tracks",dzBin,dzMin,dzMax,dxBin,dxMin,dxMax,"");
244  h_vx_dz->setAxisTitle("dz (cm)",1);
245  h_vx_dz->setAxisTitle("x coordinate of input track at PCA (cm)",2);
246 
247  h_vy_dz = dbe_->bookProfile("vy_dz","v_{y} vs. dz of selected tracks",dzBin,dzMin,dzMax,dxBin,dxMin,dxMax,"");
248  h_vy_dz->setAxisTitle("dz (cm)",1);
249  h_vy_dz->setAxisTitle("x coordinate of input track at PCA (cm)",2);
250 
251  h_x0 = dbe_->book1D("BeamMonitorFeedBack_x0","x coordinate of beam spot (Fit)",100,-0.01,0.01);
252  h_x0->setAxisTitle("x_{0} (cm)",1);
253  h_x0->getTH1()->SetBit(TH1::kCanRebin);
254 
255  h_y0 = dbe_->book1D("BeamMonitorFeedBack_y0","y coordinate of beam spot (Fit)",100,-0.01,0.01);
256  h_y0->setAxisTitle("y_{0} (cm)",1);
257  h_y0->getTH1()->SetBit(TH1::kCanRebin);
258 
259  h_z0 = dbe_->book1D("BeamMonitorFeedBack_z0","z coordinate of beam spot (Fit)",dzBin,dzMin,dzMax);
260  h_z0->setAxisTitle("z_{0} (cm)",1);
261  h_z0->getTH1()->SetBit(TH1::kCanRebin);
262 
263  h_sigmaX0 = dbe_->book1D("BeamMonitorFeedBack_sigmaX0","sigma x0 of beam spot (Fit)",100,0,0.05);
264  h_sigmaX0->setAxisTitle("#sigma_{X_{0}} (cm)",1);
265  h_sigmaX0->getTH1()->SetBit(TH1::kCanRebin);
266 
267  h_sigmaY0 = dbe_->book1D("BeamMonitorFeedBack_sigmaY0","sigma y0 of beam spot (Fit)",100,0,0.05);
268  h_sigmaY0->setAxisTitle("#sigma_{Y_{0}} (cm)",1);
269  h_sigmaY0->getTH1()->SetBit(TH1::kCanRebin);
270 
271  h_sigmaZ0 = dbe_->book1D("BeamMonitorFeedBack_sigmaZ0","sigma z0 of beam spot (Fit)",100,0,10);
272  h_sigmaZ0->setAxisTitle("#sigma_{Z_{0}} (cm)",1);
273  h_sigmaZ0->getTH1()->SetBit(TH1::kCanRebin);
274 
275  // Histograms of all reco tracks (without cuts):
276  h_trkPt=dbe_->book1D("trkPt","p_{T} of all reco'd tracks (no selection)",200,0.,50.);
277  h_trkPt->setAxisTitle("p_{T} (GeV/c)",1);
278 
279  h_trkVz=dbe_->book1D("trkVz","Z coordinate of PCA of all reco'd tracks (no selection)",dzBin,dzMin,dzMax);
280  h_trkVz->setAxisTitle("V_{Z} (cm)",1);
281 
282  cutFlowTable = dbe_->book1D("cutFlowTable","Cut flow table of track selection", 9, 0, 9 );
283 
284  // Results of previous good fit:
285  fitResults=dbe_->book2D("fitResults","Results of previous good beam fit",2,0,2,8,0,8);
286  fitResults->setAxisTitle("Fitted Beam Spot (cm)",1);
287  fitResults->setBinLabel(8,"x_{0}",2);
288  fitResults->setBinLabel(7,"y_{0}",2);
289  fitResults->setBinLabel(6,"z_{0}",2);
290  fitResults->setBinLabel(5,"#sigma_{Z}",2);
291  fitResults->setBinLabel(4,"#frac{dx}{dz} (rad)",2);
292  fitResults->setBinLabel(3,"#frac{dy}{dz} (rad)",2);
293  fitResults->setBinLabel(2,"#sigma_{X}",2);
294  fitResults->setBinLabel(1,"#sigma_{Y}",2);
295  fitResults->setBinLabel(1,"Mean",1);
296  fitResults->setBinLabel(2,"Stat. Error",1);
297  fitResults->getTH1()->SetOption("text");
298 
299  // Histos of PrimaryVertices:
300  dbe_->setCurrentFolder(monitorName_+"PrimaryVertex");
301 
302  h_nVtx = dbe_->book1D("vtxNbr","Reconstructed Vertices in Event",20,-0.5,19.5);
303  h_nVtx->setAxisTitle("Num. of reco. vertices",1);
304 
305  // Monitor only the PV with highest sum pt of assoc. trks:
306  h_PVx[0] = dbe_->book1D("PVX","x coordinate of Primary Vtx",50,-0.01,0.01);
307  h_PVx[0]->setAxisTitle("PVx (cm)",1);
308  h_PVx[0]->getTH1()->SetBit(TH1::kCanRebin);
309 
310  h_PVy[0] = dbe_->book1D("PVY","y coordinate of Primary Vtx",50,-0.01,0.01);
311  h_PVy[0]->setAxisTitle("PVy (cm)",1);
312  h_PVy[0]->getTH1()->SetBit(TH1::kCanRebin);
313 
314  h_PVz[0] = dbe_->book1D("PVZ","z coordinate of Primary Vtx",dzBin,dzMin,dzMax);
315  h_PVz[0]->setAxisTitle("PVz (cm)",1);
316 
317  h_PVx[1] = dbe_->book1D("PVXFit","x coordinate of Primary Vtx (Last Fit)",50,-0.01,0.01);
318  h_PVx[1]->setAxisTitle("PVx (cm)",1);
319  h_PVx[1]->getTH1()->SetBit(TH1::kCanRebin);
320 
321  h_PVy[1] = dbe_->book1D("PVYFit","y coordinate of Primary Vtx (Last Fit)",50,-0.01,0.01);
322  h_PVy[1]->setAxisTitle("PVy (cm)",1);
323  h_PVy[1]->getTH1()->SetBit(TH1::kCanRebin);
324 
325  h_PVz[1] = dbe_->book1D("PVZFit","z coordinate of Primary Vtx (Last Fit)",dzBin,dzMin,dzMax);
326  h_PVz[1]->setAxisTitle("PVz (cm)",1);
327 
328  h_PVxz = dbe_->bookProfile("PVxz","PVx vs. PVz",dzBin/2,dzMin,dzMax,dxBin/2,dxMin,dxMax,"");
329  h_PVxz->setAxisTitle("PVz (cm)",1);
330  h_PVxz->setAxisTitle("PVx (cm)",2);
331 
332  h_PVyz = dbe_->bookProfile("PVyz","PVy vs. PVz",dzBin/2,dzMin,dzMax,dxBin/2,dxMin,dxMax,"");
333  h_PVyz->setAxisTitle("PVz (cm)",1);
334  h_PVyz->setAxisTitle("PVy (cm)",2);
335 
336  // Results of previous good fit:
337  pvResults=dbe_->book2D("pvResults","Results of fitting Primary Vertices",2,0,2,6,0,6);
338  pvResults->setAxisTitle("Fitted Primary Vertex (cm)",1);
339  pvResults->setBinLabel(6,"PVx",2);
340  pvResults->setBinLabel(5,"PVy",2);
341  pvResults->setBinLabel(4,"PVz",2);
342  pvResults->setBinLabel(3,"#sigma_{X}",2);
343  pvResults->setBinLabel(2,"#sigma_{Y}",2);
344  pvResults->setBinLabel(1,"#sigma_{Z}",2);
345  pvResults->setBinLabel(1,"Mean",1);
346  pvResults->setBinLabel(2,"Stat. Error",1);
347  pvResults->getTH1()->SetOption("text");
348 
349  // Summary plots:
350  dbe_->setCurrentFolder(monitorName_+"EventInfo");
351  reportSummary = dbe_->get(monitorName_+"EventInfo/reportSummary");
353 
354  reportSummary = dbe_->bookFloat("reportSummary");
355  if(reportSummary) reportSummary->Fill(0./0.);
356 
357  char histo[20];
358  dbe_->setCurrentFolder(monitorName_+"EventInfo/reportSummaryContents");
359  for (int n = 0; n < nFitElements_; n++) {
360  switch(n){
361  case 0 : sprintf(histo,"x0_status"); break;
362  case 1 : sprintf(histo,"y0_status"); break;
363  case 2 : sprintf(histo,"z0_status"); break;
364  }
366  }
367 
368  for (int i = 0; i < nFitElements_; i++) {
369  summaryContent_[i] = 0.;
370  reportSummaryContents[i]->Fill(0./0.);
371  }
372 
373  dbe_->setCurrentFolder(monitorName_+"EventInfo");
374 
375  reportSummaryMap = dbe_->get(monitorName_+"EventInfo/reportSummaryMap");
377 
378  reportSummaryMap = dbe_->book2D("reportSummaryMap", "Beam Spot Summary Map", 1, 0, 1, 3, 0, 3);
380  reportSummaryMap->setAxisTitle("Fitted Beam Spot",2);
381  reportSummaryMap->setBinLabel(1," ",1);
382  reportSummaryMap->setBinLabel(1,"x_{0}",2);
383  reportSummaryMap->setBinLabel(2,"y_{0}",2);
384  reportSummaryMap->setBinLabel(3,"z_{0}",2);
385  for (int i = 0; i < nFitElements_; i++) {
386  reportSummaryMap->setBinContent(1,i+1,-1.);
387  }
388 }
389 
390 //--------------------------------------------------------
391 void BeamMonitor::beginRun(const edm::Run& r, const EventSetup& context) {
392 
393  ftimestamp = r.beginTime().value();
394  tmpTime = ftimestamp >> 32;
396  const char* eventTime = formatFitTime(tmpTime);
397  std::cout << "TimeOffset = " << eventTime << std::endl;
398  TDatime da(eventTime);
399  if (debug_) {
400  edm::LogInfo("BeamMonitor") << "TimeOffset = ";
401  da.Print();
402  }
403  for (std::map<TString,MonitorElement*>::iterator it = hs.begin();
404  it != hs.end(); ++it) {
405  if ((*it).first.Contains("time"))
406  (*it).second->getTH1()->GetXaxis()->SetTimeOffset(da.Convert(kTRUE));
407  }
408 }
409 
410 //--------------------------------------------------------
412  const EventSetup& context) {
413  int nthlumi = lumiSeg.luminosityBlock();
414  const edm::TimeValue_t fbegintimestamp = lumiSeg.beginTime().value();
415  const std::time_t ftmptime = fbegintimestamp >> 32;
416 
417  if (countLumi_ == 0) {
419  refBStime[0] = refPVtime[0] = ftmptime;
420  }
423 
424  if (onlineMode_ && (nthlumi < nextlumi_)) return;
425 
426  if (onlineMode_) {
427  if (nthlumi > nextlumi_) {
428  if (countLumi_ != 0 && processed_) FitAndFill(lumiSeg,lastlumi_,nextlumi_,nthlumi);
429  nextlumi_ = nthlumi;
430  edm::LogInfo("BeamMonitor") << "Next Lumi to Fit: " << nextlumi_ << endl;
431  if (refBStime[0] == 0) refBStime[0] = ftmptime;
432  if (refPVtime[0] == 0) refPVtime[0] = ftmptime;
433  }
434  }
435  else{
436  if (processed_) FitAndFill(lumiSeg,lastlumi_,nextlumi_,nthlumi);
437  nextlumi_ = nthlumi;
438  edm::LogInfo("BeamMonitor") << "Next Lumi to Fit: " << nextlumi_ << endl;
439  if (refBStime[0] == 0) refBStime[0] = ftmptime;
440  if (refPVtime[0] == 0) refPVtime[0] = ftmptime;
441  }
442  countLumi_++;
443  if (processed_) processed_ = false;
444  edm::LogInfo("BeamMonitor") << "Begin of Lumi: " << nthlumi << endl;
445 }
446 
447 // ----------------------------------------------------------
449  const EventSetup& iSetup ) {
450  bool readEvent_ = true;
451  const int nthlumi = iEvent.luminosityBlock();
452  if (onlineMode_ && (nthlumi < nextlumi_)) {
453  readEvent_ = false;
454  edm::LogInfo("BeamMonitor") << "Spilt event from previous lumi section!" << std::endl;
455  }
456  if (onlineMode_ && (nthlumi > nextlumi_)) {
457  readEvent_ = false;
458  edm::LogInfo("BeamMonitor") << "Spilt event from next lumi section!!!" << std::endl;
459  }
460 
461  if (readEvent_) {
462  countEvt_++;
463  theBeamFitter->readEvent(iEvent);
464  Handle<reco::BeamSpot> recoBeamSpotHandle;
465  iEvent.getByLabel(bsSrc_,recoBeamSpotHandle);
466  refBS = *recoBeamSpotHandle;
467 
469  const char* tmpfile;
470  TH1F * tmphisto;
471  tmpfile = (cutFlowTable->getName()).c_str();
472  tmphisto = static_cast<TH1F *>((theBeamFitter->getCutFlow())->Clone("tmphisto"));
473  cutFlowTable->getTH1()->SetBins(tmphisto->GetNbinsX(),tmphisto->GetXaxis()->GetXmin(),tmphisto->GetXaxis()->GetXmax());
474  if (countEvt_ == 1) // SetLabel just once
475  for(int n=0; n < tmphisto->GetNbinsX(); n++)
476  cutFlowTable->setBinLabel(n+1,tmphisto->GetXaxis()->GetBinLabel(n+1),1);
477 
478  cutFlowTable = dbe_->book1D(tmpfile,tmphisto);
479 
481  iEvent.getByLabel(tracksLabel_, TrackCollection);
482  const reco::TrackCollection *tracks = TrackCollection.product();
483  for ( reco::TrackCollection::const_iterator track = tracks->begin();
484  track != tracks->end(); ++track ) {
485  h_trkPt->Fill(track->pt());
486  h_trkVz->Fill(track->vz());
487  }
488 
489  //------ Primary Vertices
491 
492  if (iEvent.getByLabel(pvSrc_, PVCollection )) {
493  int nPVcount = 0;
494  for (reco::VertexCollection::const_iterator pv = PVCollection->begin(); pv != PVCollection->end(); ++pv) {
495  //--- vertex selection
496  if (pv->isFake() || pv->tracksSize()==0) continue;
497  nPVcount++; // count non fake pv
498  if (pv->ndof() < minVtxNdf_ || (pv->ndof()+3.)/pv->tracksSize() < 2*minVtxWgt_) continue;
499  h_PVx[0]->Fill(pv->x());
500  h_PVy[0]->Fill(pv->y());
501  h_PVz[0]->Fill(pv->z());
502  h_PVxz->Fill(pv->z(),pv->x());
503  h_PVyz->Fill(pv->z(),pv->y());
504  }
505  if (nPVcount>0) h_nVtx->Fill(nPVcount*1.);
506  }
507  processed_ = true;
508  }//end of read event
509 
510 }
511 
512 
513 //--------------------------------------------------------
515  const EventSetup& iSetup) {
516  int nthlumi = lumiSeg.id().luminosityBlock();
517  edm::LogInfo("BeamMonitor") << "Lumi of the last event before endLuminosityBlock: " << nthlumi << endl;
518 
519  if (onlineMode_ && nthlumi < nextlumi_) return;
520  const edm::TimeValue_t fendtimestamp = lumiSeg.endTime().value();
521  const std::time_t fendtime = fendtimestamp >> 32;
522  tmpTime = refBStime[1] = refPVtime[1] = fendtime;
523 }
524 
525 //--------------------------------------------------------
526 void BeamMonitor::FitAndFill(const LuminosityBlock& lumiSeg,int &lastlumi,int &nextlumi,int &nthlumi){
527  if (onlineMode_ && (nthlumi <= nextlumi)) return;
528 
529  int currentlumi = nextlumi;
530  edm::LogInfo("BeamMonitor") << "Lumi of the current fit: " << currentlumi << endl;
531  lastlumi = currentlumi;
532  endLumiOfBSFit_ = currentlumi;
533  endLumiOfPVFit_ = currentlumi;
534 
535  if (onlineMode_) { // filling LS gap
536  // FIXME: need to add protection for the case if the gap is at the resetting LS!
537  const int countLS_bs = hs["x0_lumi"]->getTH1()->GetEntries();
538  const int countLS_pv = hs["PVx_lumi"]->getTH1()->GetEntries();
539  edm::LogInfo("BeamMonitor") << "countLS_bs = " << countLS_bs << " ; countLS_pv = " << countLS_pv << std::endl;
540  int LSgap_bs = currentlumi/fitNLumi_ - countLS_bs;
541  int LSgap_pv = currentlumi/fitPVNLumi_ - countLS_pv;
542  if (currentlumi%fitNLumi_ == 0)
543  LSgap_bs--;
544  if (currentlumi%fitPVNLumi_ == 0)
545  LSgap_pv--;
546  edm::LogInfo("BeamMonitor") << "LSgap_bs = " << LSgap_bs << " ; LSgap_pv = " << LSgap_pv << std::endl;
547  // filling previous fits if LS gap ever exists
548  for (int ig = 0; ig < LSgap_bs; ig++) {
549  hs["x0_lumi"]->ShiftFillLast( 0., 0., fitNLumi_ );
550  hs["y0_lumi"]->ShiftFillLast( 0., 0., fitNLumi_ );
551  hs["z0_lumi"]->ShiftFillLast( 0., 0., fitNLumi_ );
552  hs["sigmaX0_lumi"]->ShiftFillLast( 0., 0., fitNLumi_ );
553  hs["sigmaY0_lumi"]->ShiftFillLast( 0., 0., fitNLumi_ );
554  hs["sigmaZ0_lumi"]->ShiftFillLast( 0., 0., fitNLumi_ );
555  }
556  for (int ig = 0; ig < LSgap_pv; ig++) {
557  hs["PVx_lumi"]->ShiftFillLast( 0., 0., fitPVNLumi_ );
558  hs["PVy_lumi"]->ShiftFillLast( 0., 0., fitPVNLumi_ );
559  hs["PVz_lumi"]->ShiftFillLast( 0., 0., fitPVNLumi_ );
560  }
561  const int previousLS = h_nTrk_lumi->getTH1()->GetEntries();
562  for (int i=1;i < (currentlumi - previousLS);i++)
564  }
565 
566  edm::LogInfo("BeamMonitor") << "Time lapsed since last scroll = " << tmpTime - refTime << std:: endl;
567 
568  if (testScroll(tmpTime,refTime)) {
569  scrollTH1(hs["x0_time"]->getTH1(),refTime);
570  scrollTH1(hs["y0_time"]->getTH1(),refTime);
571  scrollTH1(hs["z0_time"]->getTH1(),refTime);
572  scrollTH1(hs["sigmaX0_time"]->getTH1(),refTime);
573  scrollTH1(hs["sigmaY0_time"]->getTH1(),refTime);
574  scrollTH1(hs["sigmaZ0_time"]->getTH1(),refTime);
575  scrollTH1(hs["PVx_time"]->getTH1(),refTime);
576  scrollTH1(hs["PVy_time"]->getTH1(),refTime);
577  scrollTH1(hs["PVz_time"]->getTH1(),refTime);
578  }
579 
580  bool doPVFit = false;
581 
582  if (fitPVNLumi_ > 0) {
583  if (onlineMode_) {
584  if (currentlumi%fitPVNLumi_ == 0)
585  doPVFit = true;
586  }
587  else
588  if (countLumi_%fitPVNLumi_ == 0)
589  doPVFit = true;
590  }
591  else
592  doPVFit = true;
593 
594  if (doPVFit) {
595 
596  edm::LogInfo("BeamMonitor") << "Do PV Fitting for LS = " << beginLumiOfPVFit_ << " to " << endLumiOfPVFit_ << std::endl;
597  // Primary Vertex Fit:
598  if (h_PVx[0]->getTH1()->GetEntries() > minNrVertices_) {
599  pvResults->Reset();
600  char tmpTitle[50];
601  sprintf(tmpTitle,"%s %i %s %i","Fitted Primary Vertex (cm) of LS: ",beginLumiOfPVFit_," to ",endLumiOfPVFit_);
602  pvResults->setAxisTitle(tmpTitle,1);
603 
604  TF1 *fgaus = new TF1("fgaus","gaus");
605  double mean,width,meanErr,widthErr;
606  fgaus->SetLineColor(4);
607  h_PVx[0]->getTH1()->Fit("fgaus","QLM0");
608  mean = fgaus->GetParameter(1);
609  width = fgaus->GetParameter(2);
610  meanErr = fgaus->GetParError(1);
611  widthErr = fgaus->GetParError(2);
612  hs["PVx_lumi"]->ShiftFillLast(mean,width,fitPVNLumi_);
613  hs["PVx_lumi_all"]->setBinContent(currentlumi,mean);
614  hs["PVx_lumi_all"]->setBinError(currentlumi,width);
615  int nthBin = tmpTime - refTime;
616  if (nthBin < 0)
617  edm::LogInfo("BeamMonitor") << "Event time outside current range of time histograms!" << std::endl;
618  if (nthBin > 0) {
619  hs["PVx_time"]->setBinContent(nthBin,mean);
620  hs["PVx_time"]->setBinError(nthBin,width);
621  }
622  int jthBin = tmpTime - startTime;
623  if (jthBin > 0) {
624  hs["PVx_time_all"]->setBinContent(jthBin,mean);
625  hs["PVx_time_all"]->setBinError(jthBin,width);
626  }
627  pvResults->setBinContent(1,6,mean);
628  pvResults->setBinContent(1,3,width);
629  pvResults->setBinContent(2,6,meanErr);
630  pvResults->setBinContent(2,3,widthErr);
631 
632  dbe_->setCurrentFolder(monitorName_+"PrimaryVertex/");
633  const char* tmpfile;
634  TH1D * tmphisto;
635  // snap shot of the fit
636  tmpfile= (h_PVx[1]->getName()).c_str();
637  tmphisto = static_cast<TH1D *>((h_PVx[0]->getTH1())->Clone("tmphisto"));
638  h_PVx[1]->getTH1()->SetBins(tmphisto->GetNbinsX(),tmphisto->GetXaxis()->GetXmin(),tmphisto->GetXaxis()->GetXmax());
639  h_PVx[1] = dbe_->book1D(tmpfile,h_PVx[0]->getTH1F());
640  h_PVx[1]->getTH1()->Fit("fgaus","QLM");
641 
642 
643  h_PVy[0]->getTH1()->Fit("fgaus","QLM0");
644  mean = fgaus->GetParameter(1);
645  width = fgaus->GetParameter(2);
646  meanErr = fgaus->GetParError(1);
647  widthErr = fgaus->GetParError(2);
648  hs["PVy_lumi"]->ShiftFillLast(mean,width,fitPVNLumi_);
649  hs["PVy_lumi_all"]->setBinContent(currentlumi,mean);
650  hs["PVy_lumi_all"]->setBinError(currentlumi,width);
651  if (nthBin > 0) {
652  hs["PVy_time"]->setBinContent(nthBin,mean);
653  hs["PVy_time"]->setBinError(nthBin,width);
654  }
655  if (jthBin > 0) {
656  hs["PVy_time_all"]->setBinContent(jthBin,mean);
657  hs["PVy_time_all"]->setBinError(jthBin,width);
658  }
659  pvResults->setBinContent(1,5,mean);
660  pvResults->setBinContent(1,2,width);
661  pvResults->setBinContent(2,5,meanErr);
662  pvResults->setBinContent(2,2,widthErr);
663  // snap shot of the fit
664  tmpfile= (h_PVy[1]->getName()).c_str();
665  tmphisto = static_cast<TH1D *>((h_PVy[0]->getTH1())->Clone("tmphisto"));
666  h_PVy[1]->getTH1()->SetBins(tmphisto->GetNbinsX(),tmphisto->GetXaxis()->GetXmin(),tmphisto->GetXaxis()->GetXmax());
667  h_PVy[1]->update();
668  h_PVy[1] = dbe_->book1D(tmpfile,h_PVy[0]->getTH1F());
669  h_PVy[1]->getTH1()->Fit("fgaus","QLM");
670 
671 
672  h_PVz[0]->getTH1()->Fit("fgaus","QLM0");
673  mean = fgaus->GetParameter(1);
674  width = fgaus->GetParameter(2);
675  meanErr = fgaus->GetParError(1);
676  widthErr = fgaus->GetParError(2);
677  hs["PVz_lumi"]->ShiftFillLast(mean,width,fitPVNLumi_);
678  hs["PVz_lumi_all"]->setBinContent(currentlumi,mean);
679  hs["PVz_lumi_all"]->setBinError(currentlumi,width);
680  if (nthBin > 0) {
681  hs["PVz_time"]->setBinContent(nthBin,mean);
682  hs["PVz_time"]->setBinError(nthBin,width);
683  }
684  if (jthBin > 0) {
685  hs["PVz_time_all"]->setBinContent(jthBin,mean);
686  hs["PVz_time_all"]->setBinError(jthBin,width);
687  }
688  pvResults->setBinContent(1,4,mean);
689  pvResults->setBinContent(1,1,width);
690  pvResults->setBinContent(2,4,meanErr);
691  pvResults->setBinContent(2,1,widthErr);
692  // snap shot of the fit
693  tmpfile= (h_PVz[1]->getName()).c_str();
694  tmphisto = static_cast<TH1D *>((h_PVz[0]->getTH1())->Clone("tmphisto"));
695  h_PVz[1]->getTH1()->SetBins(tmphisto->GetNbinsX(),tmphisto->GetXaxis()->GetXmin(),tmphisto->GetXaxis()->GetXmax());
696  h_PVz[1]->update();
697  h_PVz[1] = dbe_->book1D(tmpfile,h_PVz[0]->getTH1F());
698  h_PVz[1]->getTH1()->Fit("fgaus","QLM");
699 
700  }
701  }
702 
703  if (resetPVNLumi_ > 0 &&
704  ((onlineMode_ && currentlumi%resetPVNLumi_ == 0) ||
705  (!onlineMode_ && countLumi_%resetPVNLumi_ == 0))) {
706  h_PVx[0]->Reset();
707  h_PVy[0]->Reset();
708  h_PVz[0]->Reset();
709  beginLumiOfPVFit_ = 0;
710  refPVtime[0] = 0;
711  }
712 
713  // Beam Spot Fit:
714  vector<BSTrkParameters> theBSvector = theBeamFitter->getBSvector();
715  h_nTrk_lumi->ShiftFillLast( theBSvector.size() );
716 
717  bool countFitting = false;
718  if (theBSvector.size() > nthBSTrk_ && theBSvector.size() >= min_Ntrks_) {
719  countFitting = true;
720  }
721 
722  if (resetHistos_) {
723  edm::LogInfo("BeamMonitor") << "Resetting Histograms" << endl;
724  h_d0_phi0->Reset();
725  h_vx_vy->Reset();
726  h_vx_dz->Reset();
727  h_vy_dz->Reset();
728  h_trk_z0->Reset();
730  resetHistos_ = false;
731  }
732 
733  edm::LogInfo("BeamMonitor") << "Fill histos, start from " << nthBSTrk_ + 1 << "th record of selected tracks" << endl;
734  unsigned int itrk = 0;
735  for (vector<BSTrkParameters>::const_iterator BSTrk = theBSvector.begin();
736  BSTrk != theBSvector.end();
737  ++BSTrk, ++itrk){
738  if (itrk >= nthBSTrk_){
739  h_d0_phi0->Fill( BSTrk->phi0(), BSTrk->d0() );
740  double vx = BSTrk->vx();
741  double vy = BSTrk->vy();
742  double z0 = BSTrk->z0();
743  h_vx_vy->Fill( vx, vy );
744  h_vx_dz->Fill( z0, vx );
745  h_vy_dz->Fill( z0, vy );
746  h_trk_z0->Fill( z0 );
747  }
748  }
749  nthBSTrk_ = theBSvector.size(); // keep track of num of tracks filled so far
750  if (countFitting) edm::LogInfo("BeamMonitor") << "Num of tracks collected = " << nthBSTrk_ << endl;
751 
752  if (fitNLumi_ > 0) {
753  if (onlineMode_){
754  if (currentlumi%fitNLumi_!=0) {
755 // for (std::map<TString,MonitorElement*>::iterator itAll = hs.begin();
756 // itAll != hs.end(); ++itAll) {
757 // if ((*itAll).first.Contains("all")) {
758 // (*itAll).second->setBinContent(currentlumi,0.);
759 // (*itAll).second->setBinError(currentlumi,0.);
760 // }
761 // }
762  return;
763  }
764  }
765  else
766  if (countLumi_%fitNLumi_!=0) return;
767  }
768 
769  edm::LogInfo("BeamMonitor") << " [DebugTime] refBStime[0] = " << refBStime[0]
770  << "; address = " << &refBStime[0] << std::endl;
771  edm::LogInfo("BeamMonitor") << " [DebugTime] refBStime[1] = " << refBStime[1]
772  << "; address = " << &refBStime[1] << std::endl;
773 
774  if (countFitting) {
775  nFits_++;
776  int * fitLS = theBeamFitter->getFitLSRange();
777  edm::LogInfo("BeamMonitor") << "[BeamFitter] Do BeamSpot Fit for LS = " << fitLS[0] << " to " << fitLS[1] << std::endl;
778  edm::LogInfo("BeamMonitor") << "[BeamMonitor] Do BeamSpot Fit for LS = " << beginLumiOfBSFit_ << " to " << endLumiOfBSFit_ << std::endl;
781  if (bs.type() > 0) // with good beamwidth fit
782  preBS = bs; // cache good fit results
783  edm::LogInfo("BeamMonitor") << "\n RESULTS OF DEFAULT FIT:" << endl;
784  edm::LogInfo("BeamMonitor") << bs << endl;
785  edm::LogInfo("BeamMonitor") << "[BeamFitter] fitting done \n" << endl;
786 
787  hs["x0_lumi"]->ShiftFillLast( bs.x0(), bs.x0Error(), fitNLumi_ );
788  hs["y0_lumi"]->ShiftFillLast( bs.y0(), bs.y0Error(), fitNLumi_ );
789  hs["z0_lumi"]->ShiftFillLast( bs.z0(), bs.z0Error(), fitNLumi_ );
790  hs["sigmaX0_lumi"]->ShiftFillLast( bs.BeamWidthX(), bs.BeamWidthXError(), fitNLumi_ );
791  hs["sigmaY0_lumi"]->ShiftFillLast( bs.BeamWidthY(), bs.BeamWidthYError(), fitNLumi_ );
792  hs["sigmaZ0_lumi"]->ShiftFillLast( bs.sigmaZ(), bs.sigmaZ0Error(), fitNLumi_ );
793  hs["x0_lumi_all"]->setBinContent(currentlumi,bs.x0());
794  hs["x0_lumi_all"]->setBinError(currentlumi,bs.x0Error());
795  hs["y0_lumi_all"]->setBinContent(currentlumi,bs.y0());
796  hs["y0_lumi_all"]->setBinError(currentlumi,bs.y0Error());
797  hs["z0_lumi_all"]->setBinContent(currentlumi,bs.z0());
798  hs["z0_lumi_all"]->setBinError(currentlumi,bs.z0Error());
799  hs["sigmaX0_lumi_all"]->setBinContent(currentlumi, bs.BeamWidthX());
800  hs["sigmaX0_lumi_all"]->setBinError(currentlumi, bs.BeamWidthXError());
801  hs["sigmaY0_lumi_all"]->setBinContent(currentlumi, bs.BeamWidthY());
802  hs["sigmaY0_lumi_all"]->setBinError(currentlumi, bs.BeamWidthYError());
803  hs["sigmaZ0_lumi_all"]->setBinContent(currentlumi, bs.sigmaZ());
804  hs["sigmaZ0_lumi_all"]->setBinError(currentlumi, bs.sigmaZ0Error());
805 
806  int nthBin = tmpTime - refTime;
807  if (nthBin > 0) {
808  hs["x0_time"]->setBinContent(nthBin, bs.x0());
809  hs["y0_time"]->setBinContent(nthBin, bs.y0());
810  hs["z0_time"]->setBinContent(nthBin, bs.z0());
811  hs["sigmaX0_time"]->setBinContent(nthBin, bs.BeamWidthX());
812  hs["sigmaY0_time"]->setBinContent(nthBin, bs.BeamWidthY());
813  hs["sigmaZ0_time"]->setBinContent(nthBin, bs.sigmaZ());
814  hs["x0_time"]->setBinError(nthBin, bs.x0Error());
815  hs["y0_time"]->setBinError(nthBin, bs.y0Error());
816  hs["z0_time"]->setBinError(nthBin, bs.z0Error());
817  hs["sigmaX0_time"]->setBinError(nthBin, bs.BeamWidthXError());
818  hs["sigmaY0_time"]->setBinError(nthBin, bs.BeamWidthYError());
819  hs["sigmaZ0_time"]->setBinError(nthBin, bs.sigmaZ0Error());
820  }
821 
822  int jthBin = tmpTime - startTime;
823  if (jthBin > 0) {
824  hs["x0_time_all"]->setBinContent(jthBin, bs.x0());
825  hs["y0_time_all"]->setBinContent(jthBin, bs.y0());
826  hs["z0_time_all"]->setBinContent(jthBin, bs.z0());
827  hs["sigmaX0_time_all"]->setBinContent(jthBin, bs.BeamWidthX());
828  hs["sigmaY0_time_all"]->setBinContent(jthBin, bs.BeamWidthY());
829  hs["sigmaZ0_time_all"]->setBinContent(jthBin, bs.sigmaZ());
830  hs["x0_time_all"]->setBinError(jthBin, bs.x0Error());
831  hs["y0_time_all"]->setBinError(jthBin, bs.y0Error());
832  hs["z0_time_all"]->setBinError(jthBin, bs.z0Error());
833  hs["sigmaX0_time_all"]->setBinError(jthBin, bs.BeamWidthXError());
834  hs["sigmaY0_time_all"]->setBinError(jthBin, bs.BeamWidthYError());
835  hs["sigmaZ0_time_all"]->setBinError(jthBin, bs.sigmaZ0Error());
836  }
837 
838  h_x0->Fill( bs.x0());
839  h_y0->Fill( bs.y0());
840  h_z0->Fill( bs.z0());
841  if (bs.type() > 0) { // with good beamwidth fit
842  h_sigmaX0->Fill( bs.BeamWidthX());
843  h_sigmaY0->Fill( bs.BeamWidthY());
844  }
845  h_sigmaZ0->Fill( bs.sigmaZ());
846 
847  if (nthBSTrk_ >= 2*min_Ntrks_) {
848  double amp = std::sqrt(bs.x0()*bs.x0()+bs.y0()*bs.y0());
849  double alpha = std::atan2(bs.y0(),bs.x0());
850  TF1 *f1 = new TF1("f1","[0]*sin(x-[1])",-3.14,3.14);
851  f1->SetParameters(amp,alpha);
852  f1->SetParLimits(1,amp-0.1,amp+0.1);
853  f1->SetParLimits(2,alpha-0.577,alpha+0.577);
854  f1->SetLineColor(4);
855  h_d0_phi0->getTProfile()->Fit("f1","QR");
856 
857  double mean = bs.z0();
858  double width = bs.sigmaZ();
859  TF1 *fgaus = new TF1("fgaus","gaus");
860  fgaus->SetParameters(mean,width);
861  fgaus->SetLineColor(4);
862  h_trk_z0->getTH1()->Fit("fgaus","QLRM","",mean-3*width,mean+3*width);
863  }
864 
865  fitResults->Reset();
866  int * LSRange = theBeamFitter->getFitLSRange();
867  char tmpTitle[50];
868  sprintf(tmpTitle,"%s %i %s %i","Fitted Beam Spot (cm) of LS: ",LSRange[0]," to ",LSRange[1]);
869  fitResults->setAxisTitle(tmpTitle,1);
870  fitResults->setBinContent(1,8,bs.x0());
871  fitResults->setBinContent(1,7,bs.y0());
872  fitResults->setBinContent(1,6,bs.z0());
873  fitResults->setBinContent(1,5,bs.sigmaZ());
874  fitResults->setBinContent(1,4,bs.dxdz());
875  fitResults->setBinContent(1,3,bs.dydz());
876  if (bs.type() > 0) { // with good beamwidth fit
879  }
880  else { // fill cached widths
883  }
884 
885  fitResults->setBinContent(2,8,bs.x0Error());
886  fitResults->setBinContent(2,7,bs.y0Error());
887  fitResults->setBinContent(2,6,bs.z0Error());
891  if (bs.type() > 0) { // with good beamwidth fit
894  }
895  else { // fill cached width errors
898  }
899 
900  // count good fit
901  // if (std::fabs(refBS.x0()-bs.x0())/bs.x0Error() < deltaSigCut_) { // disabled temporarily
902  summaryContent_[0] += 1.;
903  // }
904  // if (std::fabs(refBS.y0()-bs.y0())/bs.y0Error() < deltaSigCut_) { // disabled temporarily
905  summaryContent_[1] += 1.;
906  // }
907  // if (std::fabs(refBS.z0()-bs.z0())/bs.z0Error() < deltaSigCut_) { // disabled temporarily
908  summaryContent_[2] += 1.;
909  // }
910  } // beam fit is good
911  else { // beam fit fails
913  edm::LogInfo("BeamMonitor") << "[BeamMonitor] Beam fit fails!!! \n" << endl;
914  edm::LogInfo("BeamMonitor") << "[BeamMonitor] Output beam spot for DIP \n" << endl;
915  edm::LogInfo("BeamMonitor") << bs << endl;
916 
917  hs["sigmaX0_lumi"]->ShiftFillLast( bs.BeamWidthX(), bs.BeamWidthXError(), fitNLumi_ );
918  hs["sigmaY0_lumi"]->ShiftFillLast( bs.BeamWidthY(), bs.BeamWidthYError(), fitNLumi_ );
919  hs["sigmaZ0_lumi"]->ShiftFillLast( bs.sigmaZ(), bs.sigmaZ0Error(), fitNLumi_ );
920  hs["x0_lumi"]->ShiftFillLast( bs.x0(), bs.x0Error(), fitNLumi_ );
921  hs["y0_lumi"]->ShiftFillLast( bs.y0(), bs.y0Error(), fitNLumi_ );
922  hs["z0_lumi"]->ShiftFillLast( bs.z0(), bs.z0Error(), fitNLumi_ );
923  } // end of beam fit fails
924  } // end of countFitting
925  else { // no fit
926  // Overwrite Fit LS and fit time when no event processed or no track selected
929  if (theBeamFitter->runPVandTrkFitter()) {} // Dump fake beam spot for DIP
931  edm::LogInfo("BeamMonitor") << "[BeamMonitor] No fitting \n" << endl;
932  edm::LogInfo("BeamMonitor") << "[BeamMonitor] Output fake beam spot for DIP \n" << endl;
933  edm::LogInfo("BeamMonitor") << bs << endl;
934 
935  hs["sigmaX0_lumi"]->ShiftFillLast( bs.BeamWidthX(), bs.BeamWidthXError(), fitNLumi_ );
936  hs["sigmaY0_lumi"]->ShiftFillLast( bs.BeamWidthY(), bs.BeamWidthYError(), fitNLumi_ );
937  hs["sigmaZ0_lumi"]->ShiftFillLast( bs.sigmaZ(), bs.sigmaZ0Error(), fitNLumi_ );
938  hs["x0_lumi"]->ShiftFillLast( bs.x0(), bs.x0Error(), fitNLumi_ );
939  hs["y0_lumi"]->ShiftFillLast( bs.y0(), bs.y0Error(), fitNLumi_ );
940  hs["z0_lumi"]->ShiftFillLast( bs.z0(), bs.z0Error(), fitNLumi_ );
941  }
942 
943  // Fill summary report
944  if (countFitting) {
945  for (int n = 0; n < nFitElements_; n++) {
947  }
948 
949  summarySum_ = 0;
950  for (int ii = 0; ii < nFitElements_; ii++) {
952  }
953  reportSummary_ = summarySum_ / (nFitElements_ * nFits_);
955 
956  for ( int bi = 0; bi < nFitElements_ ; bi++) {
958  }
959  }
960 
961  if (resetFitNLumi_ > 0 &&
962  ((onlineMode_ && currentlumi%resetFitNLumi_ == 0) ||
963  (!onlineMode_ && countLumi_%resetFitNLumi_ == 0))) {
964  edm::LogInfo("BeamMonitor") << "Reset track collection for beam fit!!!" <<endl;
965  resetHistos_ = true;
966  nthBSTrk_ = 0;
971  beginLumiOfBSFit_ = 0;
972  refBStime[0] = 0;
973  }
974 }
975 
976 //--------------------------------------------------------
977 void BeamMonitor::endRun(const Run& r, const EventSetup& context){
978 
979 }
980 
981 //--------------------------------------------------------
983  const EventSetup& iSetup){
984  if (!onlineMode_) endLuminosityBlock(lumiSeg, iSetup);
985 }
986 
987 //--------------------------------------------------------
988 void BeamMonitor::scrollTH1(TH1 * h, time_t ref) {
989  const char* offsetTime = formatFitTime(ref);
990  TDatime da(offsetTime);
991  if (lastNZbin > 0) {
992  double val = h->GetBinContent(lastNZbin);
993  double valErr = h->GetBinError(lastNZbin);
994  h->Reset();
995  h->GetXaxis()->SetTimeOffset(da.Convert(kTRUE));
996  int bin = (lastNZbin > buffTime ? buffTime : 1);
997  h->SetBinContent(bin,val);
998  h->SetBinError(bin,valErr);
999  }
1000  else {
1001  h->Reset();
1002  h->GetXaxis()->SetTimeOffset(da.Convert(kTRUE));
1003  }
1004 }
1005 
1006 //--------------------------------------------------------
1007 // Method to check whether to chane histogram time offset (forward only)
1008 bool BeamMonitor::testScroll(time_t & tmpTime_, time_t & refTime_){
1009  bool scroll_ = false;
1010  if (tmpTime_ - refTime_ >= intervalInSec_) {
1011  scroll_ = true;
1012  edm::LogInfo("BeamMonitor") << "Reset Time Offset" << std::endl;
1014  for (int bin = intervalInSec_; bin >= 1; bin--) {
1015  if (hs["x0_time"]->getBinContent(bin) > 0) {
1016  lastNZbin = bin;
1017  break;
1018  }
1019  }
1020  edm::LogInfo("BeamMonitor") << "Last non zero bin = " << lastNZbin << std::endl;
1021  if (tmpTime_ - refTime_ >= intervalInSec_ + lastNZbin) {
1022  edm::LogInfo("BeamMonitor") << "Time difference too large since last readout" << std::endl;
1023  lastNZbin = 0;
1024  refTime_ = tmpTime_ - buffTime;
1025  }
1026  else{
1027  edm::LogInfo("BeamMonitor") << "Offset to last record" << std::endl;
1028  int offset = ((lastNZbin > buffTime) ? (lastNZbin - buffTime) : (lastNZbin - 1));
1029  refTime_ += offset;
1030  }
1031  }
1032  return scroll_;
1033 }
1034 
void beginLuminosityBlock(const edm::LuminosityBlock &lumiSeg, const edm::EventSetup &context)
Definition: BeamMonitor.cc:411
LuminosityBlockID id() const
const char * formatFitTime(const std::time_t &)
Definition: BeamMonitor.cc:35
T getParameter(std::string const &) const
BeamMonitor(const edm::ParameterSet &)
Definition: BeamMonitor.cc:62
double z0() const
z coordinate
Definition: BeamSpot.h:69
T getUntrackedParameter(std::string const &, T const &) const
const std::string & getName(void) const
get name of ME
int i
Definition: DBlmapReader.cc:9
#define buffTime
Definition: BeamMonitor.cc:57
float alpha
Definition: AMPTWrapper.h:95
unsigned int minNrVertices_
Definition: BeamMonitor.h:94
void setBinContent(int binx, double content)
set content of bin (1-D)
double sigmaZ0Error() const
error on sigma z
Definition: BeamSpot.h:87
int nFitElements_
Definition: BeamMonitor.h:89
std::time_t refBStime[2]
Definition: BeamMonitor.h:86
MonitorElement * reportSummary
Definition: BeamMonitor.h:135
void analyze(const edm::Event &e, const edm::EventSetup &c)
Definition: BeamMonitor.cc:448
int endLumiOfPVFit_
Definition: BeamMonitor.h:83
const std::string & label
Definition: MVAComputer.cc:186
int resetFitNLumi_
Definition: BeamMonitor.h:69
bool runPVandTrkFitter()
Definition: BeamFitter.cc:394
MonitorElement * h_nVtx
Definition: BeamMonitor.h:122
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:519
double minVtxWgt_
Definition: BeamMonitor.h:96
MonitorElement * cutFlowTable
Definition: BeamMonitor.h:138
double dydzError() const
error on dydz
Definition: BeamSpot.h:91
MonitorElement * h_PVz[2]
Definition: BeamMonitor.h:125
#define CEST
bool testScroll(std::time_t &, std::time_t &)
std::vector< BSTrkParameters > getBSvector()
Definition: BeamFitter.h:62
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
MonitorElement * h_PVxz
Definition: BeamMonitor.h:126
double minVtxNdf_
Definition: BeamMonitor.h:95
bool processed_
Definition: BeamMonitor.h:99
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:10
void setRefTime(std::time_t t0, std::time_t t1)
Definition: BeamFitter.h:52
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)
edm::TimeValue_t ftimestamp
Definition: BeamMonitor.h:145
edm::LuminosityBlockNumber_t luminosityBlock() const
Definition: EventBase.h:59
void resetTrkVector()
Definition: BeamFitter.h:48
void update(void)
Mark the object updated.
Timestamp const & beginTime() const
MonitorElement * h_sigmaZ0
Definition: BeamMonitor.h:121
MonitorElement * h_nTrk_lumi
Definition: BeamMonitor.h:107
tuple histo
Definition: trackerHits.py:12
MonitorElement * fitResults
Definition: BeamMonitor.h:115
int fitPVNLumi_
Definition: BeamMonitor.h:68
void readEvent(const edm::Event &iEvent)
Definition: BeamFitter.cc:211
MonitorElement * h_sigmaY0
Definition: BeamMonitor.h:120
MonitorElement * bookFloat(const char *name)
Book float.
Definition: DQMStore.cc:456
void Fill(long long x)
edm::InputTag bsSrc_
Definition: BeamMonitor.h:63
LuminosityBlockNumber_t luminosityBlock() const
double dydz() const
dydz slope
Definition: BeamSpot.h:75
MonitorElement * h_y0
Definition: BeamMonitor.h:117
void ShiftFillLast(double y, double ye=0., int32_t xscale=1)
int intervalInSec_
Definition: BeamMonitor.h:71
void beginJob()
Definition: BeamMonitor.cc:110
int iEvent
Definition: GenABIO.cc:243
double dxdzError() const
error on dxdz
Definition: BeamSpot.h:89
MonitorElement * h_vy_dz
Definition: BeamMonitor.h:112
T sqrt(T t)
Definition: SSEVec.h:28
Float_t reportSummary_
Definition: BeamMonitor.h:132
Timestamp const & endTime() const
void removeElement(const std::string &name)
Definition: DQMStore.cc:2338
MonitorElement * h_PVyz
Definition: BeamMonitor.h:127
void beginRun(const edm::Run &r, const edm::EventSetup &c)
Definition: BeamMonitor.cc:391
std::time_t refTime
Definition: BeamMonitor.h:144
MonitorElement * pvResults
Definition: BeamMonitor.h:128
int endLumiOfBSFit_
Definition: BeamMonitor.h:81
virtual void endJob()
Definition: EDAnalyzer.h:56
std::time_t refPVtime[2]
Definition: BeamMonitor.h:87
reco::BeamSpot refBS
Definition: BeamMonitor.h:103
TH1 * getTH1(void) const
int * getFitLSRange()
Definition: BeamFitter.h:70
double BeamWidthX() const
beam width X
Definition: BeamSpot.h:77
void setFitLSRange(int ls0, int ls1)
Definition: BeamFitter.h:76
TH1F * getTH1F(std::string name, std::string process, std::string rootfolder, DQMStore *dbe_, bool verb, bool clone)
MonitorElement * bookProfile(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, const char *option="s")
Definition: DQMStore.cc:833
void resetLSRange()
Definition: BeamFitter.h:50
double BeamWidthYError() const
error on beam width Y, assume error in X = Y
Definition: BeamSpot.h:96
unsigned int offset(bool)
MonitorElement * get(const std::string &path) const
get ME from full pathname (e.g. &quot;my/long/dir/my_histo&quot;)
Definition: DQMStore.cc:1270
double BeamWidthXError() const
error on beam width X, assume error in X = Y
Definition: BeamSpot.h:94
double z0Error() const
error on z
Definition: BeamSpot.h:85
unsigned long long TimeValue_t
Definition: Timestamp.h:27
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:359
double dxdz() const
dxdz slope
Definition: BeamSpot.h:73
double x0Error() const
error on x
Definition: BeamSpot.h:81
double y0Error() const
error on y
Definition: BeamSpot.h:83
DQMStore * dbe_
Definition: BeamMonitor.h:75
double maxZ_
Definition: BeamMonitor.h:93
std::time_t startTime
Definition: BeamMonitor.h:143
TimeValue_t value() const
Definition: Timestamp.cc:72
edm::InputTag pvSrc_
Definition: BeamMonitor.h:65
void FitAndFill(const edm::LuminosityBlock &lumiSeg, int &, int &, int &)
Definition: BeamMonitor.cc:526
Float_t summaryContent_[3]
Definition: BeamMonitor.h:134
MonitorElement * h_sigmaX0
Definition: BeamMonitor.h:119
reco::BeamSpot preBS
Definition: BeamMonitor.h:104
MonitorElement * h_PVy[2]
Definition: BeamMonitor.h:124
tuple tracks
Definition: testEve_cfg.py:39
int beginLumiOfPVFit_
Definition: BeamMonitor.h:82
MonitorElement * reportSummaryContents[3]
Definition: BeamMonitor.h:136
double sigmaZ() const
sigma z
Definition: BeamSpot.h:71
double BeamWidthY() const
beam width Y
Definition: BeamSpot.h:79
double b
Definition: hdecay.h:120
LuminosityBlockNumber_t luminosityBlock() const
MonitorElement * h_x0
Definition: BeamMonitor.h:116
MonitorElement * h_trk_z0
Definition: BeamMonitor.h:109
T const * product() const
Definition: Handle.h:74
int resetPVNLumi_
Definition: BeamMonitor.h:70
Timestamp const & beginTime() const
Definition: RunBase.h:45
Float_t summarySum_
Definition: BeamMonitor.h:133
TH1F * getCutFlow()
Definition: BeamFitter.h:63
edm::ParameterSet parameters_
Definition: BeamMonitor.h:61
MonitorElement * h_trkPt
Definition: BeamMonitor.h:113
MonitorElement * h_vx_dz
Definition: BeamMonitor.h:111
MonitorElement * h_z0
Definition: BeamMonitor.h:118
void resetPVFitter()
Definition: BeamFitter.h:56
TProfile * getTProfile(void) const
double y0() const
y coordinate
Definition: BeamSpot.h:67
void resetCutFlow()
Definition: BeamFitter.h:64
unsigned int nthBSTrk_
Definition: BeamMonitor.h:88
tuple cout
Definition: gather_cfg.py:41
unsigned int min_Ntrks_
Definition: BeamMonitor.h:92
void endLuminosityBlock(const edm::LuminosityBlock &lumiSeg, const edm::EventSetup &c)
Definition: BeamMonitor.cc:514
bool resetHistos_
Definition: BeamMonitor.h:98
double deltaSigCut_
Definition: BeamMonitor.h:91
TH2F * getTH2F(void) const
MonitorElement * book2D(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Book 2D histogram.
Definition: DQMStore.cc:647
std::time_t tmpTime
Definition: BeamMonitor.h:142
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
BeamFitter * theBeamFitter
Definition: BeamMonitor.h:76
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
MonitorElement * h_PVx[2]
Definition: BeamMonitor.h:123
void Reset(void)
reset ME (ie. contents, errors, etc)
int beginLumiOfBSFit_
Definition: BeamMonitor.h:80
MonitorElement * h_trkVz
Definition: BeamMonitor.h:114
void endRun(const edm::Run &r, const edm::EventSetup &c)
Definition: BeamMonitor.cc:977
std::string monitorName_
Definition: BeamMonitor.h:62
MonitorElement * h_d0_phi0
Definition: BeamMonitor.h:108
edm::InputTag tracksLabel_
Definition: BeamMonitor.h:64
reco::BeamSpot getBeamSpot()
Definition: BeamFitter.h:60
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:237
Definition: Run.h:31
MonitorElement * h_vx_vy
Definition: BeamMonitor.h:110
int countLumi_
Definition: BeamMonitor.h:79
void scrollTH1(TH1 *, std::time_t)
Definition: BeamMonitor.cc:988
BeamType type() const
return beam type
Definition: BeamSpot.h:120
MonitorElement * reportSummaryMap
Definition: BeamMonitor.h:137
bool onlineMode_
Definition: BeamMonitor.h:73
double x0() const
x coordinate
Definition: BeamSpot.h:65
std::map< TString, MonitorElement * > hs
Definition: BeamMonitor.h:129
void resetRefTime()
Definition: BeamFitter.h:51