CMS 3D CMS Logo

Vx3DHLTAnalyzer.cc
Go to the documentation of this file.
1 /*
2  \File Vx3DHLTAnalyzer.cc
3  \Display Beam-spot monitor entirely based on pixel detector information
4  the monitoring is based on a 3D fit to the vertex cloud
5  \Author Mauro Dinardo
6  \Version $ Revision: 3.5 $
7  \Date $ Date: 2010/23/02 13:15:00 $
8 */
9 
11 
14 
15 #include <Math/Minimizer.h>
16 #include <Math/Factory.h>
17 #include <Math/Functor.h>
18 
19 // ### Calling namespaces ###
20 using namespace std;
21 using namespace edm;
22 using namespace reco;
23 
25  debugMode = true;
26  nLumiFit = 2; // Number of integrated lumis to perform the fit
28  15; // If failing fits, this is the maximum number of integrated lumis after which a reset is issued
29  nLumiXaxisRange = 5000; // Correspond to about 32h of data taking: 32h * 60min * 60s / 23s per lumi-block = 5009
30  dataFromFit = true; // The Beam Spot data can be either taken from the histograms or from the fit results
31  minNentries = 20; // Minimum number of good vertices to perform the fit
32  xRange = 0.8; // [cm]
33  xStep = 0.001; // [cm]
34  yRange = 0.8; // [cm]
35  yStep = 0.001; // [cm]
36  zRange = 30.; // [cm]
37  zStep = 0.04; // [cm]
38  VxErrCorr = 1.3;
39  minVxDoF = 10.; // Good-vertex selection cut
40  // For vertex fitter without track-weight: d.o.f. = 2*NTracks - 3
41  // For vertex fitter with track-weight: d.o.f. = sum_NTracks(2*track_weight) - 3
42  minVxWgt = 0.5; // Good-vertex selection cut
43  fileName = "BeamPixelResults.txt";
44 
45  vertexCollection = consumes<VertexCollection>(
46  iConfig.getUntrackedParameter<InputTag>("vertexCollection", InputTag("pixelVertices")));
47  pixelHitCollection = consumes<SiPixelRecHitCollection>(
48  iConfig.getUntrackedParameter<InputTag>("pixelHitCollection", InputTag("siPixelRecHits")));
49 
50  debugMode = iConfig.getParameter<bool>("debugMode");
51  nLumiFit = iConfig.getParameter<unsigned int>("nLumiFit");
52  maxLumiIntegration = iConfig.getParameter<unsigned int>("maxLumiIntegration");
53  nLumiXaxisRange = iConfig.getParameter<unsigned int>("nLumiXaxisRange");
54  dataFromFit = iConfig.getParameter<bool>("dataFromFit");
55  minNentries = iConfig.getParameter<unsigned int>("minNentries");
56  xRange = iConfig.getParameter<double>("xRange");
57  xStep = iConfig.getParameter<double>("xStep");
58  yRange = iConfig.getParameter<double>("yRange");
59  yStep = iConfig.getParameter<double>("yStep");
60  zRange = iConfig.getParameter<double>("zRange");
61  zStep = iConfig.getParameter<double>("zStep");
62  VxErrCorr = iConfig.getParameter<double>("VxErrCorr");
63  minVxDoF = iConfig.getParameter<double>("minVxDoF");
64  minVxWgt = iConfig.getParameter<double>("minVxWgt");
65  fileName = iConfig.getParameter<string>("fileName");
66 
67  // ### Set internal variables ###
68  nParams = 9; // Number of free parameters in the fit
69  internalDebug = false;
70  considerVxCovariance = true; // Deconvolute vertex covariance matrix
71  pi = 3.141592653589793238;
72  // ##############################
73 }
74 
76 
77 void Vx3DHLTAnalyzer::analyze(const Event& iEvent, const EventSetup& iSetup) {
78  Handle<VertexCollection> Vx3DCollection;
79  iEvent.getByToken(vertexCollection, Vx3DCollection);
80 
81  unsigned int i, j;
82  double det;
83  VertexType MyVertex;
84 
85  if (runNumber != iEvent.id().run()) {
86  reset("scratch");
87  runNumber = iEvent.id().run();
88 
89  if (debugMode == true) {
90  stringstream debugFile;
91  string tmp(fileName);
92 
93  if (outputDebugFile.is_open() == true)
94  outputDebugFile.close();
95  tmp.erase(strlen(fileName.c_str()) - 4, 4);
96  debugFile << tmp.c_str() << "_Run" << iEvent.id().run() << ".txt";
97  outputDebugFile.open(debugFile.str().c_str(), ios::out);
98  outputDebugFile.close();
99  outputDebugFile.open(debugFile.str().c_str(), ios::app);
100  }
101 
102  dqmBeginLuminosityBlock(iEvent.getLuminosityBlock(), iSetup);
103  } else if (beginTimeOfFit != 0) {
104  totalHits += HitCounter(iEvent);
105 
106  if (internalDebug == true) {
107  edm::LogInfo("Vx3DHLTAnalyzer") << "\tI found " << totalHits << " pixel hits until now";
108  edm::LogInfo("Vx3DHLTAnalyzer") << "\tIn this event there are " << Vx3DCollection->size() << " vertex cadidates";
109  }
110 
111  for (vector<Vertex>::const_iterator it3DVx = Vx3DCollection->begin(); it3DVx != Vx3DCollection->end(); it3DVx++) {
112  if (internalDebug == true) {
113  edm::LogInfo("Vx3DHLTAnalyzer") << "\tVertex selections:";
114  edm::LogInfo("Vx3DHLTAnalyzer") << "\tEvent ID = " << iEvent.id();
115  edm::LogInfo("Vx3DHLTAnalyzer") << "\tVertex number = " << it3DVx - Vx3DCollection->begin();
116  edm::LogInfo("Vx3DHLTAnalyzer") << "\tisValid = " << it3DVx->isValid();
117  edm::LogInfo("Vx3DHLTAnalyzer") << "\tisFake = " << it3DVx->isFake();
118  edm::LogInfo("Vx3DHLTAnalyzer") << "\tnodof = " << it3DVx->ndof();
119  edm::LogInfo("Vx3DHLTAnalyzer") << "\ttracksSize = " << it3DVx->tracksSize();
120  }
121 
122  if ((it3DVx->isValid() == true) && (it3DVx->isFake() == false) && (it3DVx->ndof() >= minVxDoF) &&
123  (it3DVx->tracksSize() > 0) && ((it3DVx->ndof() + 3.) / ((double)it3DVx->tracksSize()) >= 2. * minVxWgt)) {
124  for (i = 0; i < DIM; i++) {
125  for (j = 0; j < DIM; j++) {
126  MyVertex.Covariance[i][j] = it3DVx->covariance(i, j);
127  if (isNotFinite(MyVertex.Covariance[i][j]) == true)
128  break;
129  }
130 
131  if (j != DIM)
132  break;
133  }
134 
135  if (i == DIM)
136  det = std::fabs(MyVertex.Covariance[0][0]) *
137  (std::fabs(MyVertex.Covariance[1][1]) * std::fabs(MyVertex.Covariance[2][2]) -
138  MyVertex.Covariance[1][2] * MyVertex.Covariance[1][2]) -
139  MyVertex.Covariance[0][1] * (MyVertex.Covariance[0][1] * std::fabs(MyVertex.Covariance[2][2]) -
140  MyVertex.Covariance[0][2] * MyVertex.Covariance[1][2]) +
141  MyVertex.Covariance[0][2] * (MyVertex.Covariance[0][1] * MyVertex.Covariance[1][2] -
142  MyVertex.Covariance[0][2] * std::fabs(MyVertex.Covariance[1][1]));
143 
144  if ((i == DIM) && (det > 0.)) {
145  if (internalDebug == true)
146  edm::LogInfo("Vx3DHLTAnalyzer") << "\tVertex accepted !";
147 
148  MyVertex.x = it3DVx->x();
149  MyVertex.y = it3DVx->y();
150  MyVertex.z = it3DVx->z();
151  Vertices.push_back(MyVertex);
152 
153  Vx_X->Fill(it3DVx->x());
154  Vx_Y->Fill(it3DVx->y());
155  Vx_Z->Fill(it3DVx->z());
156 
157  Vx_ZX->Fill(it3DVx->z(), it3DVx->x());
158  Vx_ZY->Fill(it3DVx->z(), it3DVx->y());
159  Vx_XY->Fill(it3DVx->x(), it3DVx->y());
160 
161  Vx_X_Cum->Fill(it3DVx->x());
162  Vx_Y_Cum->Fill(it3DVx->y());
163  Vx_Z_Cum->Fill(it3DVx->z());
164 
165  Vx_ZX_Cum->Fill(it3DVx->z(), it3DVx->x());
166  Vx_ZY_Cum->Fill(it3DVx->z(), it3DVx->y());
167  Vx_XY_Cum->Fill(it3DVx->x(), it3DVx->y());
168  } else if (internalDebug == true) {
169  edm::LogInfo("Vx3DHLTAnalyzer") << "\tVertex discarded !";
170 
171  for (i = 0; i < DIM; i++)
172  for (j = 0; j < DIM; j++)
173  edm::LogInfo("Vx3DHLTAnalyzer") << "(i,j) --> " << i << "," << j << " --> " << MyVertex.Covariance[i][j];
174  }
175  } else if (internalDebug == true)
176  edm::LogInfo("Vx3DHLTAnalyzer") << "\tVertex discarded !";
177  }
178  }
179 }
180 
182  Handle<SiPixelRecHitCollection> rechitspixel;
183  iEvent.getByToken(pixelHitCollection, rechitspixel);
184 
185  unsigned int counter = 0;
186 
187  for (SiPixelRecHitCollection::const_iterator j = rechitspixel->begin(); j != rechitspixel->end(); j++)
188  for (edmNew::DetSet<SiPixelRecHit>::const_iterator h = j->begin(); h != j->end(); h++)
189  counter += h->cluster()->size();
190 
191  return counter;
192 }
193 
194 string Vx3DHLTAnalyzer::formatTime(const time_t& t) {
195  char ts[25];
196  strftime(ts, sizeof(ts), "%Y.%m.%d %H:%M:%S %Z", gmtime(&t));
197 
198  string ts_string(ts);
199 
200  return ts_string;
201 }
202 
203 double Vx3DHLTAnalyzer::Gauss3DFunc(const double* par) {
204  double K[DIM][DIM]; // Covariance Matrix
205  double M[DIM][DIM]; // K^-1
206  double det;
207  double sumlog = 0.;
208 
209  // par[0] = K(0,0) --> Var[X]
210  // par[1] = K(1,1) --> Var[Y]
211  // par[2] = K(2,2) --> Var[Z]
212  // par[3] = K(0,1) = K(1,0) --> Cov[X,Y]
213  // par[4] = K(1,2) = K(2,1) --> Cov[Y,Z] --> dy/dz
214  // par[5] = K(0,2) = K(2,0) --> Cov[X,Z] --> dx/dz
215  // par[6] = mean x
216  // par[7] = mean y
217  // par[8] = mean z
218 
219  counterVx = 0;
220  for (unsigned int i = 0; i < Vertices.size(); i++) {
221  if ((std::sqrt((Vertices[i].x - xPos) * (Vertices[i].x - xPos) + (Vertices[i].y - yPos) * (Vertices[i].y - yPos)) <=
222  maxTransRadius) &&
223  (std::fabs(Vertices[i].z - zPos) <= maxLongLength)) {
224  if (considerVxCovariance == true) {
225  K[0][0] = std::fabs(par[0]) + VxErrCorr * VxErrCorr * std::fabs(Vertices[i].Covariance[0][0]);
226  K[1][1] = std::fabs(par[1]) + VxErrCorr * VxErrCorr * std::fabs(Vertices[i].Covariance[1][1]);
227  K[2][2] = std::fabs(par[2]) + VxErrCorr * VxErrCorr * std::fabs(Vertices[i].Covariance[2][2]);
228  K[0][1] = K[1][0] = par[3] + VxErrCorr * VxErrCorr * Vertices[i].Covariance[0][1];
229  K[1][2] = K[2][1] = par[4] * (std::fabs(par[2]) - std::fabs(par[1])) - par[5] * par[3] +
231  K[0][2] = K[2][0] = par[5] * (std::fabs(par[2]) - std::fabs(par[0])) - par[4] * par[3] +
233  } else {
234  K[0][0] = std::fabs(par[0]);
235  K[1][1] = std::fabs(par[1]);
236  K[2][2] = std::fabs(par[2]);
237  K[0][1] = K[1][0] = par[3];
238  K[1][2] = K[2][1] = par[4] * (std::fabs(par[2]) - std::fabs(par[1])) - par[5] * par[3];
239  K[0][2] = K[2][0] = par[5] * (std::fabs(par[2]) - std::fabs(par[0])) - par[4] * par[3];
240  }
241 
242  det = K[0][0] * (K[1][1] * K[2][2] - K[1][2] * K[1][2]) - K[0][1] * (K[0][1] * K[2][2] - K[0][2] * K[1][2]) +
243  K[0][2] * (K[0][1] * K[1][2] - K[0][2] * K[1][1]);
244 
245  M[0][0] = (K[1][1] * K[2][2] - K[1][2] * K[1][2]) / det;
246  M[1][1] = (K[0][0] * K[2][2] - K[0][2] * K[0][2]) / det;
247  M[2][2] = (K[0][0] * K[1][1] - K[0][1] * K[0][1]) / det;
248  M[0][1] = M[1][0] = (K[0][2] * K[1][2] - K[0][1] * K[2][2]) / det;
249  M[1][2] = M[2][1] = (K[0][2] * K[0][1] - K[1][2] * K[0][0]) / det;
250  M[0][2] = M[2][0] = (K[0][1] * K[1][2] - K[0][2] * K[1][1]) / det;
251 
252  sumlog += double(DIM) * std::log(2. * pi) + std::log(std::fabs(det)) +
253  (M[0][0] * (Vertices[i].x - par[6]) * (Vertices[i].x - par[6]) +
254  M[1][1] * (Vertices[i].y - par[7]) * (Vertices[i].y - par[7]) +
255  M[2][2] * (Vertices[i].z - par[8]) * (Vertices[i].z - par[8]) +
256  2. * M[0][1] * (Vertices[i].x - par[6]) * (Vertices[i].y - par[7]) +
257  2. * M[1][2] * (Vertices[i].y - par[7]) * (Vertices[i].z - par[8]) +
258  2. * M[0][2] * (Vertices[i].x - par[6]) * (Vertices[i].z - par[8]));
259 
260  counterVx++;
261  }
262  }
263 
264  return sumlog;
265 }
266 
267 int Vx3DHLTAnalyzer::MyFit(vector<double>* vals) {
268  // ###############################################
269  // # RETURN CODE: #
270  // # >0 == NOT OK - fit status (MINUIT manual) #
271  // # 0 == OK #
272  // # -1 == NOT OK - not finite edm #
273  // # -2 == NOT OK - not enough "minNentries" #
274  // # -3 == NOT OK - not finite errors #
275  // # -4 == NOT OK - negative determinant #
276  // # -5 == NOT OK - maxLumiIntegration reached #
277  // # -6 == NOT OK - FatalRootError: @SUB=Minuit2 #
278  // ###############################################
279 
280  if ((vals != nullptr) && (vals->size() == nParams * 2)) {
281  double nSigmaXY = 10.;
282  double nSigmaZ = 10.;
283  double parDistanceXY = 1e-3; // Unit: [cm]
284  double parDistanceZ = 1e-2; // Unit: [cm]
285  double parDistanceddZ = 1e-3; // Unit: [rad]
286  double parDistanceCxy = 1e-5; // Unit: [cm^2]
287  double bestEdm;
288 
289  const unsigned int trials = 4;
290  double largerDist[trials] = {0.1, 5., 10., 100.};
291 
292  double covxz, covyz, det;
293  double deltaMean;
294  int bestMovementX = 1;
295  int bestMovementY = 1;
296  int bestMovementZ = 1;
297  int goodData;
298 
299  double edm;
300 
301  vector<double>::const_iterator it = vals->begin();
302 
303  ROOT::Math::Minimizer* Gauss3D = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
304  Gauss3D->SetErrorDef(1.0);
305  if (internalDebug == true)
306  Gauss3D->SetPrintLevel(3);
307  else
308  Gauss3D->SetPrintLevel(0);
309 
310  ROOT::Math::Functor _Gauss3DFunc(this, &Vx3DHLTAnalyzer::Gauss3DFunc, nParams);
311  Gauss3D->SetFunction(_Gauss3DFunc);
312 
313  if (internalDebug == true)
314  edm::LogInfo("Vx3DHLTAnalyzer") << "\t@@@ START FITTING @@@";
315 
316  // @@@ Fit at X-deltaMean | X | X+deltaMean @@@
317  bestEdm = 1.;
318  for (int i = 0; i < 3; i++) {
319  deltaMean = (double(i) - 1.) * std::sqrt(*(it + 0));
320  if (internalDebug == true)
321  edm::LogInfo("Vx3DHLTAnalyzer") << "\tdeltaMean --> " << deltaMean;
322 
323  Gauss3D->Clear();
324 
325  Gauss3D->SetVariable(0, "var x ", *(it + 0), parDistanceXY * parDistanceXY);
326  Gauss3D->SetVariable(1, "var y ", *(it + 1), parDistanceXY * parDistanceXY);
327  Gauss3D->SetVariable(2, "var z ", *(it + 2), parDistanceZ * parDistanceZ);
328  Gauss3D->SetVariable(3, "cov xy", *(it + 3), parDistanceCxy);
329  Gauss3D->SetVariable(4, "dydz ", *(it + 4), parDistanceddZ);
330  Gauss3D->SetVariable(5, "dxdz ", *(it + 5), parDistanceddZ);
331  Gauss3D->SetVariable(6, "mean x", *(it + 6) + deltaMean, parDistanceXY);
332  Gauss3D->SetVariable(7, "mean y", *(it + 7), parDistanceXY);
333  Gauss3D->SetVariable(8, "mean z", *(it + 8), parDistanceZ);
334 
335  // Set the central positions of the centroid for vertex rejection
336  xPos = *(it + 6) + deltaMean;
337  yPos = *(it + 7);
338  zPos = *(it + 8);
339 
340  // Set dimensions of the centroid for vertex rejection
341  maxTransRadius = nSigmaXY * std::sqrt(std::fabs(*(it + 0)) + std::fabs(*(it + 1)) / 2.);
342  maxLongLength = nSigmaZ * std::sqrt(std::fabs(*(it + 2)));
343 
344  try {
345  Gauss3D->Minimize();
346  goodData = Gauss3D->Status();
347  edm = Gauss3D->Edm();
348  } catch (cms::Exception& er) {
349  edm::LogError("Vx3DHLTAnalyzer") << "\tCaught Minuit2 exception @ Fit at X: \n" << er.what();
350  goodData = -6;
351  edm = 1.;
352  continue;
353  }
354 
355  if (counterVx < minNentries)
356  goodData = -2;
357  else if (isNotFinite(edm) == true) {
358  goodData = -1;
359  if (internalDebug == true)
360  edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite edm !";
361  } else if (goodData != -6)
362  for (unsigned int j = 0; j < nParams; j++)
363  if (isNotFinite(Gauss3D->Errors()[j]) == true) {
364  goodData = -3;
365  if (internalDebug == true)
366  edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite errors !";
367  break;
368  }
369  if (goodData == 0) {
370  covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
371  Gauss3D->X()[5] * Gauss3D->X()[3];
372  covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
373  Gauss3D->X()[4] * Gauss3D->X()[3];
374 
375  det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
376  Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
377  covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
378  if (det < 0.) {
379  goodData = -4;
380  if (internalDebug == true)
381  edm::LogInfo("Vx3DHLTAnalyzer") << "\tNegative determinant !";
382  }
383  }
384 
385  if ((goodData == 0) && (std::fabs(edm) < bestEdm)) {
386  bestEdm = edm;
387  bestMovementX = i;
388  }
389  }
390  if (internalDebug == true)
391  edm::LogInfo("Vx3DHLTAnalyzer") << "\tFound bestMovementX --> " << bestMovementX;
392 
393  // @@@ Fit at Y-deltaMean | Y | Y+deltaMean @@@
394  bestEdm = 1.;
395  for (int i = 0; i < 3; i++) {
396  deltaMean = (double(i) - 1.) * std::sqrt(*(it + 1));
397  if (internalDebug == true) {
398  edm::LogInfo("Vx3DHLTAnalyzer") << "\tdeltaMean --> " << deltaMean;
399  edm::LogInfo("Vx3DHLTAnalyzer") << "\tdeltaMean X --> " << (double(bestMovementX) - 1.) * std::sqrt(*(it + 0));
400  }
401 
402  Gauss3D->Clear();
403 
404  Gauss3D->SetVariable(0, "var x ", *(it + 0), parDistanceXY * parDistanceXY);
405  Gauss3D->SetVariable(1, "var y ", *(it + 1), parDistanceXY * parDistanceXY);
406  Gauss3D->SetVariable(2, "var z ", *(it + 2), parDistanceZ * parDistanceZ);
407  Gauss3D->SetVariable(3, "cov xy", *(it + 3), parDistanceCxy);
408  Gauss3D->SetVariable(4, "dydz ", *(it + 4), parDistanceddZ);
409  Gauss3D->SetVariable(5, "dxdz ", *(it + 5), parDistanceddZ);
410  Gauss3D->SetVariable(6, "mean x", *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0)), parDistanceXY);
411  Gauss3D->SetVariable(7, "mean y", *(it + 7) + deltaMean, parDistanceXY);
412  Gauss3D->SetVariable(8, "mean z", *(it + 8), parDistanceZ);
413 
414  // Set the central positions of the centroid for vertex rejection
415  xPos = *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0));
416  yPos = *(it + 7) + deltaMean;
417  zPos = *(it + 8);
418 
419  // Set dimensions of the centroid for vertex rejection
420  maxTransRadius = nSigmaXY * std::sqrt(std::fabs(*(it + 0)) + std::fabs(*(it + 1)) / 2.);
421  maxLongLength = nSigmaZ * std::sqrt(std::fabs(*(it + 2)));
422 
423  try {
424  Gauss3D->Minimize();
425  goodData = Gauss3D->Status();
426  edm = Gauss3D->Edm();
427  } catch (cms::Exception& er) {
428  edm::LogError("Vx3DHLTAnalyzer") << "\tCaught Minuit2 exception @ Fit at Y: \n" << er.what();
429  goodData = -6;
430  edm = 1.;
431  continue;
432  }
433 
434  if (counterVx < minNentries)
435  goodData = -2;
436  else if (isNotFinite(edm) == true) {
437  goodData = -1;
438  if (internalDebug == true)
439  edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite edm !";
440  } else if (goodData != -6)
441  for (unsigned int j = 0; j < nParams; j++)
442  if (isNotFinite(Gauss3D->Errors()[j]) == true) {
443  goodData = -3;
444  if (internalDebug == true)
445  edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite errors !";
446  break;
447  }
448  if (goodData == 0) {
449  covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
450  Gauss3D->X()[5] * Gauss3D->X()[3];
451  covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
452  Gauss3D->X()[4] * Gauss3D->X()[3];
453 
454  det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
455  Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
456  covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
457  if (det < 0.) {
458  goodData = -4;
459  if (internalDebug == true)
460  edm::LogInfo("Vx3DHLTAnalyzer") << "\tNegative determinant !";
461  }
462  }
463 
464  if ((goodData == 0) && (std::fabs(edm) < bestEdm)) {
465  bestEdm = edm;
466  bestMovementY = i;
467  }
468  }
469  if (internalDebug == true)
470  edm::LogInfo("Vx3DHLTAnalyzer") << "\tFound bestMovementY --> " << bestMovementY;
471 
472  // @@@ Fit at Z-deltaMean | Z | Z+deltaMean @@@
473  bestEdm = 1.;
474  for (int i = 0; i < 3; i++) {
475  deltaMean = (double(i) - 1.) * std::sqrt(*(it + 2));
476  if (internalDebug == true) {
477  edm::LogInfo("Vx3DHLTAnalyzer") << "\tdeltaMean --> " << deltaMean;
478  edm::LogInfo("Vx3DHLTAnalyzer") << "\tdeltaMean X --> " << (double(bestMovementX) - 1.) * std::sqrt(*(it + 0));
479  edm::LogInfo("Vx3DHLTAnalyzer") << "\tdeltaMean Y --> " << (double(bestMovementY) - 1.) * std::sqrt(*(it + 1));
480  }
481 
482  Gauss3D->Clear();
483 
484  Gauss3D->SetVariable(0, "var x ", *(it + 0), parDistanceXY * parDistanceXY);
485  Gauss3D->SetVariable(1, "var y ", *(it + 1), parDistanceXY * parDistanceXY);
486  Gauss3D->SetVariable(2, "var z ", *(it + 2), parDistanceZ * parDistanceZ);
487  Gauss3D->SetVariable(3, "cov xy", *(it + 3), parDistanceCxy);
488  Gauss3D->SetVariable(4, "dydz ", *(it + 4), parDistanceddZ);
489  Gauss3D->SetVariable(5, "dxdz ", *(it + 5), parDistanceddZ);
490  Gauss3D->SetVariable(6, "mean x", *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0)), parDistanceXY);
491  Gauss3D->SetVariable(7, "mean y", *(it + 7) + (double(bestMovementY) - 1.) * std::sqrt(*(it + 1)), parDistanceXY);
492  Gauss3D->SetVariable(8, "mean z", *(it + 8) + deltaMean, parDistanceZ);
493 
494  // Set the central positions of the centroid for vertex rejection
495  xPos = *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0));
496  yPos = *(it + 7) + (double(bestMovementY) - 1.) * std::sqrt(*(it + 1));
497  zPos = *(it + 8) + deltaMean;
498 
499  // Set dimensions of the centroid for vertex rejection
500  maxTransRadius = nSigmaXY * std::sqrt(std::fabs(*(it + 0)) + std::fabs(*(it + 1)) / 2.);
501  maxLongLength = nSigmaZ * std::sqrt(std::fabs(*(it + 2)));
502 
503  try {
504  Gauss3D->Minimize();
505  goodData = Gauss3D->Status();
506  edm = Gauss3D->Edm();
507  } catch (cms::Exception& er) {
508  edm::LogError("Vx3DHLTAnalyzer") << "\tCaught Minuit2 exception @ Fit at Z: \n" << er.what();
509  goodData = -6;
510  edm = 1.;
511  continue;
512  }
513 
514  if (counterVx < minNentries)
515  goodData = -2;
516  else if (isNotFinite(edm) == true) {
517  goodData = -1;
518  if (internalDebug == true)
519  edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite edm !";
520  } else if (goodData != -6)
521  for (unsigned int j = 0; j < nParams; j++)
522  if (isNotFinite(Gauss3D->Errors()[j]) == true) {
523  goodData = -3;
524  if (internalDebug == true)
525  edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite errors !";
526  break;
527  }
528  if (goodData == 0) {
529  covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
530  Gauss3D->X()[5] * Gauss3D->X()[3];
531  covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
532  Gauss3D->X()[4] * Gauss3D->X()[3];
533 
534  det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
535  Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
536  covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
537  if (det < 0.) {
538  goodData = -4;
539  if (internalDebug == true)
540  edm::LogInfo("Vx3DHLTAnalyzer") << "\tNegative determinant !";
541  }
542  }
543 
544  if ((goodData == 0) && (std::fabs(edm) < bestEdm)) {
545  bestEdm = edm;
546  bestMovementZ = i;
547  }
548  }
549  if (internalDebug == true)
550  edm::LogInfo("Vx3DHLTAnalyzer") << "\tFound bestMovementZ --> " << bestMovementZ;
551 
552  Gauss3D->Clear();
553 
554  // @@@ FINAL FIT @@@
555  Gauss3D->SetVariable(0, "var x ", *(it + 0), parDistanceXY * parDistanceXY);
556  Gauss3D->SetVariable(1, "var y ", *(it + 1), parDistanceXY * parDistanceXY);
557  Gauss3D->SetVariable(2, "var z ", *(it + 2), parDistanceZ * parDistanceZ);
558  Gauss3D->SetVariable(3, "cov xy", *(it + 3), parDistanceCxy);
559  Gauss3D->SetVariable(4, "dydz ", *(it + 4), parDistanceddZ);
560  Gauss3D->SetVariable(5, "dxdz ", *(it + 5), parDistanceddZ);
561  Gauss3D->SetVariable(6, "mean x", *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0)), parDistanceXY);
562  Gauss3D->SetVariable(7, "mean y", *(it + 7) + (double(bestMovementY) - 1.) * std::sqrt(*(it + 1)), parDistanceXY);
563  Gauss3D->SetVariable(8, "mean z", *(it + 8) + (double(bestMovementZ) - 1.) * std::sqrt(*(it + 2)), parDistanceZ);
564 
565  // Set the central positions of the centroid for vertex rejection
566  xPos = *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0));
567  yPos = *(it + 7) + (double(bestMovementY) - 1.) * std::sqrt(*(it + 1));
568  zPos = *(it + 8) + (double(bestMovementZ) - 1.) * std::sqrt(*(it + 2));
569 
570  // Set dimensions of the centroid for vertex rejection
571  maxTransRadius = nSigmaXY * std::sqrt(std::fabs(*(it + 0)) + std::fabs(*(it + 1)) / 2.);
572  maxLongLength = nSigmaZ * std::sqrt(std::fabs(*(it + 2)));
573 
574  try {
575  Gauss3D->Minimize();
576  goodData = Gauss3D->Status();
577  edm = Gauss3D->Edm();
578  } catch (cms::Exception& er) {
579  edm::LogError("Vx3DHLTAnalyzer") << "\tCaught Minuit2 exception @ Final fit: \n" << er.what();
580  goodData = -6;
581  edm = 1.;
582  }
583 
584  if (counterVx < minNentries)
585  goodData = -2;
586  else if (isNotFinite(edm) == true) {
587  goodData = -1;
588  if (internalDebug == true)
589  edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite edm !";
590  } else if (goodData != -6)
591  for (unsigned int j = 0; j < nParams; j++)
592  if (isNotFinite(Gauss3D->Errors()[j]) == true) {
593  goodData = -3;
594  if (internalDebug == true)
595  edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite errors !";
596  break;
597  }
598  if (goodData == 0) {
599  covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
600  Gauss3D->X()[5] * Gauss3D->X()[3];
601  covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
602  Gauss3D->X()[4] * Gauss3D->X()[3];
603 
604  det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
605  Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
606  covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
607  if (det < 0.) {
608  goodData = -4;
609  if (internalDebug == true)
610  edm::LogInfo("Vx3DHLTAnalyzer") << "\tNegative determinant !";
611  }
612  }
613 
614  // @@@ FIT WITH DIFFERENT PARAMETER DISTANCES @@@
615  for (unsigned int i = 0; i < trials; i++) {
616  if ((goodData != 0) && (goodData != -2)) {
617  Gauss3D->Clear();
618 
619  if (internalDebug == true)
620  edm::LogInfo("Vx3DHLTAnalyzer") << "\tFIT WITH DIFFERENT PARAMETER DISTANCES - STEP " << i + 1;
621 
622  Gauss3D->SetVariable(0, "var x ", *(it + 0), parDistanceXY * parDistanceXY * largerDist[i]);
623  Gauss3D->SetVariable(1, "var y ", *(it + 1), parDistanceXY * parDistanceXY * largerDist[i]);
624  Gauss3D->SetVariable(2, "var z ", *(it + 2), parDistanceZ * parDistanceZ * largerDist[i]);
625  Gauss3D->SetVariable(3, "cov xy", *(it + 3), parDistanceCxy * largerDist[i]);
626  Gauss3D->SetVariable(4, "dydz ", *(it + 4), parDistanceddZ * largerDist[i]);
627  Gauss3D->SetVariable(5, "dxdz ", *(it + 5), parDistanceddZ * largerDist[i]);
628  Gauss3D->SetVariable(6,
629  "mean x",
630  *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0)),
631  parDistanceXY * largerDist[i]);
632  Gauss3D->SetVariable(7,
633  "mean y",
634  *(it + 7) + (double(bestMovementY) - 1.) * std::sqrt(*(it + 1)),
635  parDistanceXY * largerDist[i]);
636  Gauss3D->SetVariable(
637  8, "mean z", *(it + 8) + (double(bestMovementZ) - 1.) * std::sqrt(*(it + 2)), parDistanceZ * largerDist[i]);
638 
639  // Set the central positions of the centroid for vertex rejection
640  xPos = *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0));
641  yPos = *(it + 7) + (double(bestMovementY) - 1.) * std::sqrt(*(it + 1));
642  zPos = *(it + 8) + (double(bestMovementZ) - 1.) * std::sqrt(*(it + 2));
643 
644  // Set dimensions of the centroid for vertex rejection
645  maxTransRadius = nSigmaXY * std::sqrt(std::fabs(*(it + 0)) + std::fabs(*(it + 1)) / 2.);
646  maxLongLength = nSigmaZ * std::sqrt(std::fabs(*(it + 2)));
647 
648  try {
649  Gauss3D->Minimize();
650  goodData = Gauss3D->Status();
651  edm = Gauss3D->Edm();
652  } catch (cms::Exception& er) {
653  edm::LogError("Vx3DHLTAnalyzer") << "\tCaught Minuit2 exception @ Fit with different distances: \n"
654  << er.what();
655  goodData = -6;
656  edm = 1.;
657  continue;
658  }
659 
660  if (counterVx < minNentries)
661  goodData = -2;
662  else if (isNotFinite(edm) == true) {
663  goodData = -1;
664  if (internalDebug == true)
665  edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite edm !";
666  } else if (goodData != -6)
667  for (unsigned int j = 0; j < nParams; j++)
668  if (isNotFinite(Gauss3D->Errors()[j]) == true) {
669  goodData = -3;
670  if (internalDebug == true)
671  edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite errors !";
672  break;
673  }
674  if (goodData == 0) {
675  covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
676  Gauss3D->X()[5] * Gauss3D->X()[3];
677  covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
678  Gauss3D->X()[4] * Gauss3D->X()[3];
679 
680  det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
681  Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
682  covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
683  if (det < 0.) {
684  goodData = -4;
685  if (internalDebug == true)
686  edm::LogInfo("Vx3DHLTAnalyzer") << "\tNegative determinant !";
687  }
688  }
689  } else
690  break;
691  }
692 
693  if (goodData == 0)
694  for (unsigned int i = 0; i < nParams; i++) {
695  vals->operator[](i) = Gauss3D->X()[i];
696  vals->operator[](i + nParams) = Gauss3D->Errors()[i];
697  }
698 
699  delete Gauss3D;
700  return goodData;
701  }
702 
703  return -1;
704 }
705 
706 void Vx3DHLTAnalyzer::reset(string ResetType) {
707  if ((debugMode == true) && (outputDebugFile.is_open() == true)) {
708  outputDebugFile << "Runnumber " << runNumber << endl;
709  outputDebugFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
710  outputDebugFile << "BeginLumiRange " << beginLumiOfFit << endl;
711  outputDebugFile << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " " << (endTimeOfFit >> 32) << endl;
712  outputDebugFile << "EndLumiRange " << endLumiOfFit << endl;
713  outputDebugFile << "LumiCounter " << lumiCounter << endl;
714  outputDebugFile << "LastLumiOfFit " << lastLumiOfFit << endl;
715  }
716 
717  if (ResetType == "scratch") {
718  runNumber = 0;
719  numberGoodFits = 0;
720  numberFits = 0;
721  lastLumiOfFit = 0;
722 
723  Vx_X->Reset();
724  Vx_Y->Reset();
725  Vx_Z->Reset();
726 
727  Vx_X_Fit->Reset();
728  Vx_Y_Fit->Reset();
729  Vx_Z_Fit->Reset();
730 
731  Vx_ZX->Reset();
732  Vx_ZY->Reset();
733  Vx_XY->Reset();
734 
735  Vx_X_Cum->Reset();
736  Vx_Y_Cum->Reset();
737  Vx_Z_Cum->Reset();
738 
739  Vx_ZX_Cum->Reset();
740  Vx_ZY_Cum->Reset();
741  Vx_XY_Cum->Reset();
742 
743  mXlumi->Reset();
744  mYlumi->Reset();
745  mZlumi->Reset();
746 
747  sXlumi->Reset();
748  sYlumi->Reset();
749  sZlumi->Reset();
750 
751  dxdzlumi->Reset();
752  dydzlumi->Reset();
753 
754  hitCounter->Reset();
755  goodVxCounter->Reset();
756  statusCounter->Reset();
757  fitResults->Reset();
758 
759  reportSummary->Fill(-1);
760  reportSummaryMap->getTH1()->SetBinContent(1, 1, -1);
761 
762  Vertices.clear();
763 
764  lumiCounter = 0;
765  totalHits = 0;
766  beginTimeOfFit = 0;
767  endTimeOfFit = 0;
768  beginLumiOfFit = 0;
769  endLumiOfFit = 0;
770 
771  if (internalDebug == true)
772  edm::LogInfo("Vx3DHLTAnalyzer") << "\tReset issued: scratch";
773  if ((debugMode == true) && (outputDebugFile.is_open() == true))
774  outputDebugFile << "Reset -scratch- issued\n" << endl;
775  } else if (ResetType == "whole") {
776  Vx_X->Reset();
777  Vx_Y->Reset();
778  Vx_Z->Reset();
779 
780  Vx_ZX->Reset();
781  Vx_ZY->Reset();
782  Vx_XY->Reset();
783 
784  Vertices.clear();
785 
786  lumiCounter = 0;
787  totalHits = 0;
788  beginTimeOfFit = 0;
789  endTimeOfFit = 0;
790  beginLumiOfFit = 0;
791  endLumiOfFit = 0;
792 
793  if (internalDebug == true)
794  edm::LogInfo("Vx3DHLTAnalyzer") << "\tReset issued: whole";
795  if ((debugMode == true) && (outputDebugFile.is_open() == true))
796  outputDebugFile << "Reset -whole- issued\n" << endl;
797  } else if (ResetType == "fit") {
798  Vx_X_Fit->Reset();
799  Vx_Y_Fit->Reset();
800  Vx_Z_Fit->Reset();
801 
802  if (internalDebug == true)
803  edm::LogInfo("Vx3DHLTAnalyzer") << "\tReset issued: fit";
804  if ((debugMode == true) && (outputDebugFile.is_open() == true))
805  outputDebugFile << "Reset -fit- issued\n" << endl;
806  } else if (ResetType == "hitCounter") {
807  totalHits = 0;
808 
809  if (internalDebug == true)
810  edm::LogInfo("Vx3DHLTAnalyzer") << "\tReset issued: hitCounter";
811  if ((debugMode == true) && (outputDebugFile.is_open() == true))
812  outputDebugFile << "Reset -hitCounter- issued\n" << endl;
813  }
814 }
815 
816 void Vx3DHLTAnalyzer::writeToFile(vector<double>* vals,
817  TimeValue_t BeginTimeOfFit,
818  TimeValue_t EndTimeOfFit,
819  unsigned int BeginLumiOfFit,
820  unsigned int EndLumiOfFit,
821  int dataType) {
822  stringstream BufferString;
823  BufferString.precision(5);
824 
825  outputFile.open(fileName.c_str(), ios::out);
826 
827  if ((outputFile.is_open() == true) && (vals != nullptr) && (vals->size() == (nParams - 1) * 2)) {
828  vector<double>::const_iterator it = vals->begin();
829 
830  outputFile << "Runnumber " << runNumber << endl;
831  outputFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
832  outputFile << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " " << (endTimeOfFit >> 32) << endl;
833  outputFile << "LumiRange " << beginLumiOfFit << " - " << endLumiOfFit << endl;
834  outputFile << "Type " << dataType << endl;
835  // 3D Vertexing with Pixel Tracks:
836  // Good data = Type 3
837  // Bad data = Type -1
838 
839  BufferString << *(it + 0);
840  outputFile << "X0 " << BufferString.str().c_str() << endl;
841  BufferString.str("");
842 
843  BufferString << *(it + 1);
844  outputFile << "Y0 " << BufferString.str().c_str() << endl;
845  BufferString.str("");
846 
847  BufferString << *(it + 2);
848  outputFile << "Z0 " << BufferString.str().c_str() << endl;
849  BufferString.str("");
850 
851  BufferString << *(it + 5);
852  outputFile << "sigmaZ0 " << BufferString.str().c_str() << endl;
853  BufferString.str("");
854 
855  BufferString << *(it + 6);
856  outputFile << "dxdz " << BufferString.str().c_str() << endl;
857  BufferString.str("");
858 
859  BufferString << *(it + 7);
860  outputFile << "dydz " << BufferString.str().c_str() << endl;
861  BufferString.str("");
862 
863  BufferString << *(it + 3);
864  outputFile << "BeamWidthX " << BufferString.str().c_str() << endl;
865  BufferString.str("");
866 
867  BufferString << *(it + 4);
868  outputFile << "BeamWidthY " << BufferString.str().c_str() << endl;
869  BufferString.str("");
870 
871  outputFile << "Cov(0,j) " << *(it + 8) << " 0 0 0 0 0 0" << endl;
872  outputFile << "Cov(1,j) 0 " << *(it + 9) << " 0 0 0 0 0" << endl;
873  outputFile << "Cov(2,j) 0 0 " << *(it + 10) << " 0 0 0 0" << endl;
874  outputFile << "Cov(3,j) 0 0 0 " << *(it + 13) << " 0 0 0" << endl;
875  outputFile << "Cov(4,j) 0 0 0 0 " << *(it + 14) << " 0 0" << endl;
876  outputFile << "Cov(5,j) 0 0 0 0 0 " << *(it + 15) << " 0" << endl;
877  outputFile << "Cov(6,j) 0 0 0 0 0 0 "
878  << ((*(it + 11)) + (*(it + 12)) + 2. * std::sqrt((*(it + 11)) * (*(it + 12)))) / 4. << endl;
879 
880  outputFile << "EmittanceX 0" << endl;
881  outputFile << "EmittanceY 0" << endl;
882  outputFile << "BetaStar 0" << endl;
883  }
884  outputFile.close();
885 
886  if ((debugMode == true) && (outputDebugFile.is_open() == true) && (vals != nullptr) &&
887  (vals->size() == (nParams - 1) * 2)) {
888  vector<double>::const_iterator it = vals->begin();
889 
890  outputDebugFile << "Runnumber " << runNumber << endl;
891  outputDebugFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
892  outputDebugFile << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " " << (endTimeOfFit >> 32) << endl;
893  outputDebugFile << "LumiRange " << beginLumiOfFit << " - " << endLumiOfFit << endl;
894  outputDebugFile << "Type " << dataType << endl;
895  // 3D Vertexing with Pixel Tracks:
896  // Good data = Type 3
897  // Bad data = Type -1
898 
899  BufferString << *(it + 0);
900  outputDebugFile << "X0 " << BufferString.str().c_str() << endl;
901  BufferString.str("");
902 
903  BufferString << *(it + 1);
904  outputDebugFile << "Y0 " << BufferString.str().c_str() << endl;
905  BufferString.str("");
906 
907  BufferString << *(it + 2);
908  outputDebugFile << "Z0 " << BufferString.str().c_str() << endl;
909  BufferString.str("");
910 
911  BufferString << *(it + 5);
912  outputDebugFile << "sigmaZ0 " << BufferString.str().c_str() << endl;
913  BufferString.str("");
914 
915  BufferString << *(it + 6);
916  outputDebugFile << "dxdz " << BufferString.str().c_str() << endl;
917  BufferString.str("");
918 
919  BufferString << *(it + 7);
920  outputDebugFile << "dydz " << BufferString.str().c_str() << endl;
921  BufferString.str("");
922 
923  BufferString << *(it + 3);
924  outputDebugFile << "BeamWidthX " << BufferString.str().c_str() << endl;
925  BufferString.str("");
926 
927  BufferString << *(it + 4);
928  outputDebugFile << "BeamWidthY " << BufferString.str().c_str() << endl;
929  BufferString.str("");
930 
931  outputDebugFile << "Cov(0,j) " << *(it + 8) << " 0 0 0 0 0 0" << endl;
932  outputDebugFile << "Cov(1,j) 0 " << *(it + 9) << " 0 0 0 0 0" << endl;
933  outputDebugFile << "Cov(2,j) 0 0 " << *(it + 10) << " 0 0 0 0" << endl;
934  outputDebugFile << "Cov(3,j) 0 0 0 " << *(it + 13) << " 0 0 0" << endl;
935  outputDebugFile << "Cov(4,j) 0 0 0 0 " << *(it + 14) << " 0 0" << endl;
936  outputDebugFile << "Cov(5,j) 0 0 0 0 0 " << *(it + 15) << " 0" << endl;
937  outputDebugFile << "Cov(6,j) 0 0 0 0 0 0 "
938  << ((*(it + 11)) + (*(it + 12)) + 2. * std::sqrt((*(it + 11)) * (*(it + 12)))) / 4. << endl;
939 
940  outputDebugFile << "EmittanceX 0" << endl;
941  outputDebugFile << "EmittanceY 0" << endl;
942  outputDebugFile << "BetaStar 0" << endl;
943 
944  outputDebugFile << "\nUsed vertices: " << counterVx << "\n" << endl;
945  }
946 }
947 
948 void Vx3DHLTAnalyzer::printFitParams(const vector<double>& fitResults) {
949  edm::LogInfo("Vx3DHLTAnalyzer") << "var x --> " << fitResults[0] << " +/- " << fitResults[0 + nParams];
950  edm::LogInfo("Vx3DHLTAnalyzer") << "var y --> " << fitResults[1] << " +/- " << fitResults[1 + nParams];
951  edm::LogInfo("Vx3DHLTAnalyzer") << "var z --> " << fitResults[2] << " +/- " << fitResults[2 + nParams];
952  edm::LogInfo("Vx3DHLTAnalyzer") << "cov xy --> " << fitResults[3] << " +/- " << fitResults[3 + nParams];
953  edm::LogInfo("Vx3DHLTAnalyzer") << "dydz --> " << fitResults[4] << " +/- " << fitResults[4 + nParams];
954  edm::LogInfo("Vx3DHLTAnalyzer") << "dxdz --> " << fitResults[5] << " +/- " << fitResults[5 + nParams];
955  edm::LogInfo("Vx3DHLTAnalyzer") << "mean x --> " << fitResults[6] << " +/- " << fitResults[6 + nParams];
956  edm::LogInfo("Vx3DHLTAnalyzer") << "mean y --> " << fitResults[7] << " +/- " << fitResults[7 + nParams];
957  edm::LogInfo("Vx3DHLTAnalyzer") << "mean z --> " << fitResults[8] << " +/- " << fitResults[8 + nParams];
958 }
959 
961  // @@@ If statement to avoid problems with non-sequential lumisections @@@
962  if ((lumiCounter == 0) && (lumiBlock.luminosityBlock() > lastLumiOfFit)) {
963  beginTimeOfFit = lumiBlock.beginTime().value();
964  beginLumiOfFit = lumiBlock.luminosityBlock();
965  lumiCounter++;
966  } else if ((lumiCounter != 0) && (lumiBlock.luminosityBlock() >= (beginLumiOfFit + lumiCounter)))
967  lumiCounter++;
968  else
969  reset("scratch");
970 }
971 
973  stringstream histTitle;
974  double minXfit, maxXfit;
975  int goodData;
976 
977  if ((nLumiFit != 0) && (lumiCounter % nLumiFit == 0) && (beginTimeOfFit != 0) && (runNumber != 0)) {
978  endTimeOfFit = lumiBlock.endTime().value();
979  endLumiOfFit = lumiBlock.luminosityBlock();
980  lastLumiOfFit = endLumiOfFit;
981  vector<double> vals;
982 
983  hitCounter->getTH1()->SetBinContent(lastLumiOfFit, (double)totalHits);
984  hitCounter->getTH1()->SetBinError(
985  lastLumiOfFit,
986  (totalHits != 0 ? 1.
987  : 0.)); // It's not sqrt(n) because we want to weight all entries in the same way for the fit
988 
989  if (dataFromFit == true) {
990  vector<double> fitResults;
991 
992  fitResults.push_back(Vx_X->getTH1()->GetRMS() * Vx_X->getTH1()->GetRMS());
993  fitResults.push_back(Vx_Y->getTH1()->GetRMS() * Vx_Y->getTH1()->GetRMS());
994  fitResults.push_back(Vx_Z->getTH1()->GetRMS() * Vx_Z->getTH1()->GetRMS());
995  fitResults.push_back(0.0);
996  fitResults.push_back(0.0);
997  fitResults.push_back(0.0);
998  fitResults.push_back(Vx_X->getTH1()->GetMean());
999  fitResults.push_back(Vx_Y->getTH1()->GetMean());
1000  fitResults.push_back(Vx_Z->getTH1()->GetMean());
1001  for (unsigned int i = 0; i < nParams; i++)
1002  fitResults.push_back(0.0);
1003 
1004  if (internalDebug == true) {
1005  edm::LogInfo("Vx3DHLTAnalyzer") << "\t@@@ Beam Spot parameters - prefit @@@";
1006 
1007  printFitParams(fitResults);
1008 
1009  edm::LogInfo("Vx3DHLTAnalyzer") << "Runnumber " << runNumber;
1010  edm::LogInfo("Vx3DHLTAnalyzer") << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " "
1011  << (beginTimeOfFit >> 32);
1012  edm::LogInfo("Vx3DHLTAnalyzer") << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " "
1013  << (endTimeOfFit >> 32);
1014  edm::LogInfo("Vx3DHLTAnalyzer") << "LumiRange " << beginLumiOfFit << " - " << endLumiOfFit;
1015  }
1016 
1017  goodData = MyFit(&fitResults);
1018 
1019  if (internalDebug == true) {
1020  edm::LogInfo("Vx3DHLTAnalyzer") << "\t@@@ Beam Spot parameters - postfit @@@";
1021 
1022  printFitParams(fitResults);
1023 
1024  edm::LogInfo("Vx3DHLTAnalyzer") << "goodData --> " << goodData;
1025  edm::LogInfo("Vx3DHLTAnalyzer") << "Used vertices --> " << counterVx;
1026  }
1027 
1028  if (goodData == 0) {
1029  vals.push_back(fitResults[6]);
1030  vals.push_back(fitResults[7]);
1031  vals.push_back(fitResults[8]);
1032  vals.push_back(std::sqrt(std::fabs(fitResults[0])));
1033  vals.push_back(std::sqrt(std::fabs(fitResults[1])));
1034  vals.push_back(std::sqrt(std::fabs(fitResults[2])));
1035  vals.push_back(fitResults[5]);
1036  vals.push_back(fitResults[4]);
1037 
1038  vals.push_back(std::pow(fitResults[6 + nParams], 2.));
1039  vals.push_back(std::pow(fitResults[7 + nParams], 2.));
1040  vals.push_back(std::pow(fitResults[8 + nParams], 2.));
1041  vals.push_back(std::pow(std::fabs(fitResults[0 + nParams]) / (2. * std::sqrt(std::fabs(fitResults[0]))), 2.));
1042  vals.push_back(std::pow(std::fabs(fitResults[1 + nParams]) / (2. * std::sqrt(std::fabs(fitResults[1]))), 2.));
1043  vals.push_back(std::pow(std::fabs(fitResults[2 + nParams]) / (2. * std::sqrt(std::fabs(fitResults[2]))), 2.));
1044  vals.push_back(std::pow(fitResults[5 + nParams], 2.));
1045  vals.push_back(std::pow(fitResults[4 + nParams], 2.));
1046  } else
1047  for (unsigned int i = 0; i < (nParams - 1) * 2; i++)
1048  vals.push_back(0.0);
1049 
1050  fitResults.clear();
1051  } else {
1052  counterVx = Vx_X->getTH1F()->GetEntries();
1053 
1054  if (Vx_X->getTH1F()->GetEntries() >= minNentries) {
1055  goodData = 0;
1056 
1057  vals.push_back(Vx_X->getTH1F()->GetMean());
1058  vals.push_back(Vx_Y->getTH1F()->GetMean());
1059  vals.push_back(Vx_Z->getTH1F()->GetMean());
1060  vals.push_back(Vx_X->getTH1F()->GetRMS());
1061  vals.push_back(Vx_Y->getTH1F()->GetRMS());
1062  vals.push_back(Vx_Z->getTH1F()->GetRMS());
1063  vals.push_back(0.0);
1064  vals.push_back(0.0);
1065 
1066  vals.push_back(std::pow(Vx_X->getTH1F()->GetMeanError(), 2.));
1067  vals.push_back(std::pow(Vx_Y->getTH1F()->GetMeanError(), 2.));
1068  vals.push_back(std::pow(Vx_Z->getTH1F()->GetMeanError(), 2.));
1069  vals.push_back(std::pow(Vx_X->getTH1F()->GetRMSError(), 2.));
1070  vals.push_back(std::pow(Vx_Y->getTH1F()->GetRMSError(), 2.));
1071  vals.push_back(std::pow(Vx_Z->getTH1F()->GetRMSError(), 2.));
1072  vals.push_back(0.0);
1073  vals.push_back(0.0);
1074  } else {
1075  goodData = -2;
1076  for (unsigned int i = 0; i < (nParams - 1) * 2; i++)
1077  vals.push_back(0.0);
1078  }
1079  }
1080 
1081  // vals[0] = X0
1082  // vals[1] = Y0
1083  // vals[2] = Z0
1084  // vals[3] = sigmaX0
1085  // vals[4] = sigmaY0
1086  // vals[5] = sigmaZ0
1087  // vals[6] = dxdz
1088  // vals[7] = dydz
1089 
1090  // vals[8] = err^2 X0
1091  // vals[9] = err^2 Y0
1092  // vals[10] = err^2 Z0
1093  // vals[11] = err^2 sigmaX0
1094  // vals[12] = err^2 sigmaY0
1095  // vals[13] = err^2 sigmaZ0
1096  // vals[14] = err^2 dxdz
1097  // vals[15] = err^2 dydz
1098 
1099  numberFits++;
1100  writeToFile(&vals, beginTimeOfFit, endTimeOfFit, beginLumiOfFit, endLumiOfFit, 3);
1101  if (internalDebug == true)
1102  edm::LogInfo("Vx3DHLTAnalyzer") << "\tUsed vertices: " << counterVx;
1103 
1104  statusCounter->getTH1()->SetBinContent(lastLumiOfFit, (double)goodData);
1105  statusCounter->getTH1()->SetBinError(lastLumiOfFit, 1e-3);
1106 
1107  // Copy vertex position histograms into to-fit histograms
1108  if (goodData == 0)
1109  reset("fit");
1110  else if (lumiCounter >= maxLumiIntegration) {
1111  reset("fit");
1112  reset("whole");
1113  }
1114 
1115  for (int i = 0; i < Vx_X_Fit->getTH1()->GetNbinsX(); i++) {
1116  Vx_X_Fit->getTH1()->SetBinContent(
1117  i + 1, Vx_X_Fit->getTH1()->GetBinContent(i + 1) + Vx_X->getTH1()->GetBinContent(i + 1));
1118  Vx_X_Fit->getTH1()->SetBinError(i + 1, sqrt(Vx_X_Fit->getTH1()->GetBinContent(i + 1)));
1119  }
1120 
1121  for (int i = 0; i < Vx_Y_Fit->getTH1()->GetNbinsX(); i++) {
1122  Vx_Y_Fit->getTH1()->SetBinContent(
1123  i + 1, Vx_Y_Fit->getTH1()->GetBinContent(i + 1) + Vx_Y->getTH1()->GetBinContent(i + 1));
1124  Vx_Y_Fit->getTH1()->SetBinError(i + 1, sqrt(Vx_Y_Fit->getTH1()->GetBinContent(i + 1)));
1125  }
1126 
1127  for (int i = 0; i < Vx_Z_Fit->getTH1()->GetNbinsX(); i++) {
1128  Vx_Z_Fit->getTH1()->SetBinContent(
1129  i + 1, Vx_Z_Fit->getTH1()->GetBinContent(i + 1) + Vx_Z->getTH1()->GetBinContent(i + 1));
1130  Vx_Z_Fit->getTH1()->SetBinError(i + 1, sqrt(Vx_Z_Fit->getTH1()->GetBinContent(i + 1)));
1131  }
1132 
1133  // Check data quality
1134  if (goodData == 0) {
1135  numberGoodFits++;
1136 
1137  histTitle << "Ongoing: fitted lumis " << beginLumiOfFit << " - " << endLumiOfFit;
1138  reset("whole");
1139  } else {
1140  if (goodData == -2)
1141  histTitle << "Ongoing: not enough evts (" << lumiCounter << " - " << maxLumiIntegration << " lumis)";
1142  else
1143  histTitle << "Ongoing: temporary problems (" << lumiCounter << " - " << maxLumiIntegration << " lumis)";
1144 
1145  if (lumiCounter >= maxLumiIntegration) {
1146  statusCounter->getTH1()->SetBinContent(lastLumiOfFit, -5);
1147  statusCounter->getTH1()->SetBinError(lastLumiOfFit, 1e-3);
1148  } else
1149  reset("hitCounter");
1150  }
1151 
1152  reportSummary->Fill((numberFits != 0 ? ((double)numberGoodFits) / ((double)numberFits) : -1));
1153  reportSummaryMap->getTH1()->SetBinContent(
1154  1, 1, (numberFits != 0 ? ((double)numberGoodFits) / ((double)numberFits) : -1));
1155 
1156  fitResults->setAxisTitle(histTitle.str(), 1);
1157 
1158  fitResults->setBinContent(1, 9, vals[0]);
1159  fitResults->setBinContent(1, 8, vals[1]);
1160  fitResults->setBinContent(1, 7, vals[2]);
1161  fitResults->setBinContent(1, 6, vals[3]);
1162  fitResults->setBinContent(1, 5, vals[4]);
1163  fitResults->setBinContent(1, 4, vals[5]);
1164  fitResults->setBinContent(1, 3, vals[6]);
1165  fitResults->setBinContent(1, 2, vals[7]);
1166  fitResults->setBinContent(1, 1, counterVx);
1167 
1168  fitResults->setBinContent(2, 9, std::sqrt(vals[8]));
1169  fitResults->setBinContent(2, 8, std::sqrt(vals[9]));
1170  fitResults->setBinContent(2, 7, std::sqrt(vals[10]));
1171  fitResults->setBinContent(2, 6, std::sqrt(vals[11]));
1172  fitResults->setBinContent(2, 5, std::sqrt(vals[12]));
1173  fitResults->setBinContent(2, 4, std::sqrt(vals[13]));
1174  fitResults->setBinContent(2, 3, std::sqrt(vals[14]));
1175  fitResults->setBinContent(2, 2, std::sqrt(vals[15]));
1176  fitResults->setBinContent(2, 1, std::sqrt(counterVx));
1177 
1178  // Linear fit to the historical plots
1179  TF1* myLinFit = new TF1(
1180  "myLinFit", "[0] + [1]*x", mXlumi->getTH1()->GetXaxis()->GetXmin(), mXlumi->getTH1()->GetXaxis()->GetXmax());
1181  myLinFit->SetLineColor(2);
1182  myLinFit->SetLineWidth(2);
1183  myLinFit->SetParName(0, "Inter.");
1184  myLinFit->SetParName(1, "Slope");
1185 
1186  mXlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[0]);
1187  mXlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[8]));
1188  myLinFit->SetParameter(0, mXlumi->getTH1()->GetMean(2));
1189  myLinFit->SetParameter(1, 0.0);
1190  mXlumi->getTH1()->Fit(myLinFit, "QR");
1191 
1192  mYlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[1]);
1193  mYlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[9]));
1194  myLinFit->SetParameter(0, mYlumi->getTH1()->GetMean(2));
1195  myLinFit->SetParameter(1, 0.0);
1196  mYlumi->getTH1()->Fit(myLinFit, "QR");
1197 
1198  mZlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[2]);
1199  mZlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[10]));
1200  myLinFit->SetParameter(0, mZlumi->getTH1()->GetMean(2));
1201  myLinFit->SetParameter(1, 0.0);
1202  mZlumi->getTH1()->Fit(myLinFit, "QR");
1203 
1204  sXlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[3]);
1205  sXlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[11]));
1206  myLinFit->SetParameter(0, sXlumi->getTH1()->GetMean(2));
1207  myLinFit->SetParameter(1, 0.0);
1208  sXlumi->getTH1()->Fit(myLinFit, "QR");
1209 
1210  sYlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[4]);
1211  sYlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[12]));
1212  myLinFit->SetParameter(0, sYlumi->getTH1()->GetMean(2));
1213  myLinFit->SetParameter(1, 0.0);
1214  sYlumi->getTH1()->Fit(myLinFit, "QR");
1215 
1216  sZlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[5]);
1217  sZlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[13]));
1218  myLinFit->SetParameter(0, sZlumi->getTH1()->GetMean(2));
1219  myLinFit->SetParameter(1, 0.0);
1220  sZlumi->getTH1()->Fit(myLinFit, "QR");
1221 
1222  dxdzlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[6]);
1223  dxdzlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[14]));
1224  myLinFit->SetParameter(0, dxdzlumi->getTH1()->GetMean(2));
1225  myLinFit->SetParameter(1, 0.0);
1226  dxdzlumi->getTH1()->Fit(myLinFit, "QR");
1227 
1228  dydzlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[7]);
1229  dydzlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[15]));
1230  myLinFit->SetParameter(0, dydzlumi->getTH1()->GetMean(2));
1231  myLinFit->SetParameter(1, 0.0);
1232  dydzlumi->getTH1()->Fit(myLinFit, "QR");
1233 
1234  myLinFit->SetParameter(0, hitCounter->getTH1()->GetMean(2));
1235  myLinFit->SetParameter(1, 0.0);
1236  hitCounter->getTH1()->Fit(myLinFit, "QR");
1237 
1238  goodVxCounter->getTH1()->SetBinContent(lastLumiOfFit, (double)counterVx);
1239  goodVxCounter->getTH1()->SetBinError(
1240  lastLumiOfFit,
1241  (counterVx != 0 ? 1.
1242  : 0.)); // It's not sqrt(n) because we want to weight all entries in the same way for the fit
1243  myLinFit->SetParameter(0, goodVxCounter->getTH1()->GetMean(2));
1244  myLinFit->SetParameter(1, 0.0);
1245  goodVxCounter->getTH1()->Fit(myLinFit, "QR");
1246 
1247  delete myLinFit;
1248  vals.clear();
1249 
1250  // Gaussian fit to 1D vertex coordinate distributions
1251  TF1* myGaussFit = new TF1("myGaussFit",
1252  "[0]*exp(-(x-[1])*(x-[1])/(2*[2]*[2]))",
1253  Vx_Z_Fit->getTH1()->GetXaxis()->GetXmin(),
1254  Vx_Z_Fit->getTH1()->GetXaxis()->GetXmax());
1255  myGaussFit->SetLineColor(2);
1256  myGaussFit->SetLineWidth(2);
1257  myGaussFit->SetParName(0, "Ampl.");
1258  myGaussFit->SetParName(1, "#mu");
1259  myGaussFit->SetParName(2, "#sigma");
1260 
1261  myGaussFit->SetParameter(0, Vx_X_Fit->getTH1()->GetMaximum());
1262  myGaussFit->SetParameter(1, Vx_X_Fit->getTH1()->GetMean());
1263  myGaussFit->SetParameter(2, Vx_X_Fit->getTH1()->GetRMS());
1264  minXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(1);
1265  for (int i = 0; i < Vx_X_Fit->getTH1()->GetNbinsX(); i++) {
1266  if (Vx_X_Fit->getTH1()->GetBinContent(i + 1) > 0) {
1267  minXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(i + 1);
1268  break;
1269  }
1270  }
1271  maxXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(Vx_X_Fit->getTH1()->GetNbinsX());
1272  for (int i = Vx_X_Fit->getTH1()->GetNbinsX(); i > 0; i--) {
1273  if (Vx_X_Fit->getTH1()->GetBinContent(i) > 0) {
1274  maxXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(i);
1275  break;
1276  }
1277  }
1278  myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1279  if (Vx_X_Fit->getTH1()->GetEntries() > 0)
1280  Vx_X_Fit->getTH1()->Fit(myGaussFit, "QRL");
1281 
1282  myGaussFit->SetParameter(0, Vx_Y_Fit->getTH1()->GetMaximum());
1283  myGaussFit->SetParameter(1, Vx_Y_Fit->getTH1()->GetMean());
1284  myGaussFit->SetParameter(2, Vx_Y_Fit->getTH1()->GetRMS());
1285  minXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(1);
1286  for (int i = 0; i < Vx_Y_Fit->getTH1()->GetNbinsX(); i++) {
1287  if (Vx_Y_Fit->getTH1()->GetBinContent(i + 1) > 0) {
1288  minXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(i + 1);
1289  break;
1290  }
1291  }
1292  maxXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(Vx_Y_Fit->getTH1()->GetNbinsX());
1293  for (int i = Vx_Y_Fit->getTH1()->GetNbinsX(); i > 0; i--) {
1294  if (Vx_Y_Fit->getTH1()->GetBinContent(i) > 0) {
1295  maxXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(i);
1296  break;
1297  }
1298  }
1299  myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1300  if (Vx_Y_Fit->getTH1()->GetEntries() > 0)
1301  Vx_Y_Fit->getTH1()->Fit(myGaussFit, "QRL");
1302 
1303  myGaussFit->SetParameter(0, Vx_Z_Fit->getTH1()->GetMaximum());
1304  myGaussFit->SetParameter(1, Vx_Z_Fit->getTH1()->GetMean());
1305  myGaussFit->SetParameter(2, Vx_Z_Fit->getTH1()->GetRMS());
1306  minXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(1);
1307  for (int i = 0; i < Vx_Z_Fit->getTH1()->GetNbinsX(); i++) {
1308  if (Vx_Z_Fit->getTH1()->GetBinContent(i + 1) > 0) {
1309  minXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(i + 1);
1310  break;
1311  }
1312  }
1313  maxXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(Vx_Z_Fit->getTH1()->GetNbinsX());
1314  for (int i = Vx_Z_Fit->getTH1()->GetNbinsX(); i > 0; i--) {
1315  if (Vx_Z_Fit->getTH1()->GetBinContent(i) > 0) {
1316  maxXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(i);
1317  break;
1318  }
1319  }
1320  myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1321  if (Vx_Z_Fit->getTH1()->GetEntries() > 0)
1322  Vx_Z_Fit->getTH1()->Fit(myGaussFit, "QRL");
1323 
1324  delete myGaussFit;
1325  } else if ((nLumiFit != 0) && (lumiCounter % nLumiFit != 0) && (beginTimeOfFit != 0) && (runNumber != 0)) {
1326  histTitle << "Ongoing: accumulating evts (" << lumiCounter % nLumiFit << " - " << nLumiFit << " in " << lumiCounter
1327  << " - " << maxLumiIntegration << " lumis)";
1328  fitResults->setAxisTitle(histTitle.str(), 1);
1329  if ((debugMode == true) && (outputDebugFile.is_open() == true)) {
1330  outputDebugFile << "\n"
1331  << "Runnumber " << runNumber << endl;
1332  outputDebugFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
1333  outputDebugFile << "BeginLumiRange " << beginLumiOfFit << endl;
1334  outputDebugFile << histTitle.str().c_str() << "\n" << endl;
1335  }
1336  } else if ((nLumiFit == 0) || (beginTimeOfFit == 0) || (runNumber == 0)) {
1337  histTitle << "Ongoing: no ongoing fits";
1338  fitResults->setAxisTitle(histTitle.str(), 1);
1339  if ((debugMode == true) && (outputDebugFile.is_open() == true))
1340  outputDebugFile << histTitle.str().c_str() << "\n" << endl;
1341 
1342  endLumiOfFit = lumiBlock.luminosityBlock();
1343 
1344  hitCounter->getTH1()->SetBinContent(endLumiOfFit, (double)totalHits);
1345  hitCounter->getTH1()->SetBinError(endLumiOfFit, std::sqrt((double)totalHits));
1346 
1347  reset("whole");
1348  }
1349 
1350  if (internalDebug == true)
1351  edm::LogInfo("Vx3DHLTAnalyzer") << "::\tHistogram title: " << histTitle.str();
1352 }
1353 
1354 void Vx3DHLTAnalyzer::bookHistograms(DQMStore::IBooker& ibooker, Run const& iRun, EventSetup const& /* iSetup */) {
1355  ibooker.setCurrentFolder("BeamPixel");
1356 
1357  Vx_X = ibooker.book1D(
1358  "F - vertex x", "Primary Vertex X Distribution", int(rint(xRange / xStep)), -xRange / 2., xRange / 2.);
1359  Vx_Y = ibooker.book1D(
1360  "F - vertex y", "Primary Vertex Y Distribution", int(rint(yRange / yStep)), -yRange / 2., yRange / 2.);
1361  Vx_Z = ibooker.book1D(
1362  "F - vertex z", "Primary Vertex Z Distribution", int(rint(zRange / zStep)), -zRange / 2., zRange / 2.);
1363  Vx_X->setAxisTitle("Primary Vertices X [cm]", 1);
1364  Vx_X->setAxisTitle("Entries [#]", 2);
1365  Vx_Y->setAxisTitle("Primary Vertices Y [cm]", 1);
1366  Vx_Y->setAxisTitle("Entries [#]", 2);
1367  Vx_Z->setAxisTitle("Primary Vertices Z [cm]", 1);
1368  Vx_Z->setAxisTitle("Entries [#]", 2);
1369 
1370  Vx_X_Fit = ibooker.book1D("G - vertex x fit",
1371  "Primary Vertex X Distribution (For Fit)",
1372  int(rint(xRange / xStep)),
1373  -xRange / 2.,
1374  xRange / 2.);
1375  Vx_Y_Fit = ibooker.book1D("G - vertex y fit",
1376  "Primary Vertex Y Distribution (For Fit)",
1377  int(rint(yRange / yStep)),
1378  -yRange / 2.,
1379  yRange / 2.);
1380  Vx_Z_Fit = ibooker.book1D("G - vertex z fit",
1381  "Primary Vertex Z Distribution (For Fit)",
1382  int(rint(zRange / zStep)),
1383  -zRange / 2.,
1384  zRange / 2.);
1385  Vx_X_Fit->setAxisTitle("Primary Vertices X [cm]", 1);
1386  Vx_X_Fit->setAxisTitle("Entries [#]", 2);
1387  Vx_Y_Fit->setAxisTitle("Primary Vertices Y [cm]", 1);
1388  Vx_Y_Fit->setAxisTitle("Entries [#]", 2);
1389  Vx_Z_Fit->setAxisTitle("Primary Vertices Z [cm]", 1);
1390  Vx_Z_Fit->setAxisTitle("Entries [#]", 2);
1391 
1392  Vx_X_Cum = ibooker.book1D("I - vertex x cum",
1393  "Primary Vertex X Distribution (Cumulative)",
1394  int(rint(xRange / xStep)),
1395  -xRange / 2.,
1396  xRange / 2.);
1397  Vx_Y_Cum = ibooker.book1D("I - vertex y cum",
1398  "Primary Vertex Y Distribution (Cumulative)",
1399  int(rint(yRange / yStep)),
1400  -yRange / 2.,
1401  yRange / 2.);
1402  Vx_Z_Cum = ibooker.book1D("I - vertex z cum",
1403  "Primary Vertex Z Distribution (Cumulative)",
1404  int(rint(zRange / zStep)),
1405  -zRange / 2.,
1406  zRange / 2.);
1407  Vx_X_Cum->setAxisTitle("Primary Vertices X [cm]", 1);
1408  Vx_X_Cum->setAxisTitle("Entries [#]", 2);
1409  Vx_Y_Cum->setAxisTitle("Primary Vertices Y [cm]", 1);
1410  Vx_Y_Cum->setAxisTitle("Entries [#]", 2);
1411  Vx_Z_Cum->setAxisTitle("Primary Vertices Z [cm]", 1);
1412  Vx_Z_Cum->setAxisTitle("Entries [#]", 2);
1413 
1414  mXlumi = ibooker.book1D(
1415  "B - muX vs lumi", "#mu_{x} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1416  mYlumi = ibooker.book1D(
1417  "B - muY vs lumi", "#mu_{y} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1418  mZlumi = ibooker.book1D(
1419  "B - muZ vs lumi", "#mu_{z} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1420  mXlumi->setAxisTitle("Lumisection [#]", 1);
1421  mXlumi->setAxisTitle("#mu_{x} [cm]", 2);
1422  mXlumi->getTH1()->SetOption("E1");
1423  mYlumi->setAxisTitle("Lumisection [#]", 1);
1424  mYlumi->setAxisTitle("#mu_{y} [cm]", 2);
1425  mYlumi->getTH1()->SetOption("E1");
1426  mZlumi->setAxisTitle("Lumisection [#]", 1);
1427  mZlumi->setAxisTitle("#mu_{z} [cm]", 2);
1428  mZlumi->getTH1()->SetOption("E1");
1429 
1430  sXlumi = ibooker.book1D(
1431  "C - sigmaX vs lumi", "#sigma_{x} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1432  sYlumi = ibooker.book1D(
1433  "C - sigmaY vs lumi", "#sigma_{y} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1434  sZlumi = ibooker.book1D(
1435  "C - sigmaZ vs lumi", "#sigma_{z} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1436  sXlumi->setAxisTitle("Lumisection [#]", 1);
1437  sXlumi->setAxisTitle("#sigma_{x} [cm]", 2);
1438  sXlumi->getTH1()->SetOption("E1");
1439  sYlumi->setAxisTitle("Lumisection [#]", 1);
1440  sYlumi->setAxisTitle("#sigma_{y} [cm]", 2);
1441  sYlumi->getTH1()->SetOption("E1");
1442  sZlumi->setAxisTitle("Lumisection [#]", 1);
1443  sZlumi->setAxisTitle("#sigma_{z} [cm]", 2);
1444  sZlumi->getTH1()->SetOption("E1");
1445 
1446  dxdzlumi = ibooker.book1D(
1447  "D - dxdz vs lumi", "dX/dZ vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1448  dydzlumi = ibooker.book1D(
1449  "D - dydz vs lumi", "dY/dZ vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1450  dxdzlumi->setAxisTitle("Lumisection [#]", 1);
1451  dxdzlumi->setAxisTitle("dX/dZ [rad]", 2);
1452  dxdzlumi->getTH1()->SetOption("E1");
1453  dydzlumi->setAxisTitle("Lumisection [#]", 1);
1454  dydzlumi->setAxisTitle("dY/dZ [rad]", 2);
1455  dydzlumi->getTH1()->SetOption("E1");
1456 
1457  Vx_ZX = ibooker.book2D("E - vertex zx",
1458  "Primary Vertex ZX Distribution",
1459  int(rint(zRange / zStep)),
1460  -zRange / 2.,
1461  zRange / 2.,
1462  int(rint(xRange / xStep)),
1463  -xRange / 2.,
1464  xRange / 2.);
1465  Vx_ZY = ibooker.book2D("E - vertex zy",
1466  "Primary Vertex ZY Distribution",
1467  int(rint(zRange / zStep)),
1468  -zRange / 2.,
1469  zRange / 2.,
1470  int(rint(yRange / yStep)),
1471  -yRange / 2.,
1472  yRange / 2.);
1473  Vx_XY = ibooker.book2D("E - vertex xy",
1474  "Primary Vertex XY Distribution",
1475  int(rint(xRange / xStep)),
1476  -xRange / 2.,
1477  xRange / 2.,
1478  int(rint(yRange / yStep)),
1479  -yRange / 2.,
1480  yRange / 2.);
1481  Vx_ZX->setAxisTitle("Primary Vertices Z [cm]", 1);
1482  Vx_ZX->setAxisTitle("Primary Vertices X [cm]", 2);
1483  Vx_ZX->setAxisTitle("Entries [#]", 3);
1484  Vx_ZY->setAxisTitle("Primary Vertices Z [cm]", 1);
1485  Vx_ZY->setAxisTitle("Primary Vertices Y [cm]", 2);
1486  Vx_ZY->setAxisTitle("Entries [#]", 3);
1487  Vx_XY->setAxisTitle("Primary Vertices X [cm]", 1);
1488  Vx_XY->setAxisTitle("Primary Vertices Y [cm]", 2);
1489  Vx_XY->setAxisTitle("Entries [#]", 3);
1490 
1491  Vx_ZX_Cum = ibooker.book2D("H - vertex zx cum",
1492  "Primary Vertex ZX Distribution (Cumulative)",
1493  int(rint(zRange / zStep)),
1494  -zRange / 2.,
1495  zRange / 2.,
1496  int(rint(xRange / xStep)),
1497  -xRange / 2.,
1498  xRange / 2.);
1499  Vx_ZY_Cum = ibooker.book2D("H - vertex zy cum",
1500  "Primary Vertex ZY Distribution (Cumulative)",
1501  int(rint(zRange / zStep)),
1502  -zRange / 2.,
1503  zRange / 2.,
1504  int(rint(yRange / yStep)),
1505  -yRange / 2.,
1506  yRange / 2.);
1507  Vx_XY_Cum = ibooker.book2D("H - vertex xy cum",
1508  "Primary Vertex XY Distribution (Cumulative)",
1509  int(rint(xRange / xStep)),
1510  -xRange / 2.,
1511  xRange / 2.,
1512  int(rint(yRange / yStep)),
1513  -yRange / 2.,
1514  yRange / 2.);
1515  Vx_ZX_Cum->setAxisTitle("Primary Vertices Z [cm]", 1);
1516  Vx_ZX_Cum->setAxisTitle("Primary Vertices X [cm]", 2);
1517  Vx_ZX_Cum->setAxisTitle("Entries [#]", 3);
1518  Vx_ZY_Cum->setAxisTitle("Primary Vertices Z [cm]", 1);
1519  Vx_ZY_Cum->setAxisTitle("Primary Vertices Y [cm]", 2);
1520  Vx_ZY_Cum->setAxisTitle("Entries [#]", 3);
1521  Vx_XY_Cum->setAxisTitle("Primary Vertices X [cm]", 1);
1522  Vx_XY_Cum->setAxisTitle("Primary Vertices Y [cm]", 2);
1523  Vx_XY_Cum->setAxisTitle("Entries [#]", 3);
1524 
1525  hitCounter = ibooker.book1D(
1526  "J - pixelHits vs lumi", "# Pixel-Hits vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1527  hitCounter->setAxisTitle("Lumisection [#]", 1);
1528  hitCounter->setAxisTitle("Pixel-Hits [#]", 2);
1529  hitCounter->getTH1()->SetOption("E1");
1530 
1531  goodVxCounter = ibooker.book1D("K - good vertices vs lumi",
1532  "# Good vertices vs. Lumisection",
1534  0.5,
1535  ((double)nLumiXaxisRange) + 0.5);
1536  goodVxCounter->setAxisTitle("Lumisection [#]", 1);
1537  goodVxCounter->setAxisTitle("Good vertices [#]", 2);
1538  goodVxCounter->getTH1()->SetOption("E1");
1539 
1540  statusCounter = ibooker.book1D(
1541  "L - status vs lumi", "App. Status vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1542  statusCounter->setAxisTitle("Lumisection [#]", 1);
1543  statusCounter->getTH1()->SetOption("E1");
1544  statusCounter->getTH1()->GetYaxis()->Set(12, -6.5, 5.5);
1545  statusCounter->getTH1()->GetYaxis()->SetBinLabel(1, "FatalRootError");
1546  statusCounter->getTH1()->GetYaxis()->SetBinLabel(2, "Max Lumi.");
1547  statusCounter->getTH1()->GetYaxis()->SetBinLabel(3, "Neg. det.");
1548  statusCounter->getTH1()->GetYaxis()->SetBinLabel(4, "Infinite err.");
1549  statusCounter->getTH1()->GetYaxis()->SetBinLabel(5, "No vtx.");
1550  statusCounter->getTH1()->GetYaxis()->SetBinLabel(6, "Infinite EDM");
1551  statusCounter->getTH1()->GetYaxis()->SetBinLabel(7, "OK");
1552  statusCounter->getTH1()->GetYaxis()->SetBinLabel(8, "MINUIT stat.");
1553  statusCounter->getTH1()->GetYaxis()->SetBinLabel(9, "MINUIT stat.");
1554  statusCounter->getTH1()->GetYaxis()->SetBinLabel(10, "MINUIT stat.");
1555  statusCounter->getTH1()->GetYaxis()->SetBinLabel(11, "MINUIT stat.");
1556  statusCounter->getTH1()->GetYaxis()->SetBinLabel(12, "MINUIT stat.");
1557 
1558  fitResults = ibooker.book2D("A - fit results", "Results of Beam Spot Fit", 2, 0., 2., 9, 0., 9.);
1559  fitResults->setAxisTitle("Ongoing: bootstrapping", 1);
1560  fitResults->setBinLabel(9, "X[cm]", 2);
1561  fitResults->setBinLabel(8, "Y[cm]", 2);
1562  fitResults->setBinLabel(7, "Z[cm]", 2);
1563  fitResults->setBinLabel(6, "#sigma_{X}[cm]", 2);
1564  fitResults->setBinLabel(5, "#sigma_{Y}[cm]", 2);
1565  fitResults->setBinLabel(4, "#sigma_{Z}[cm]", 2);
1566  fitResults->setBinLabel(3, "#frac{dX}{dZ}[rad]", 2);
1567  fitResults->setBinLabel(2, "#frac{dY}{dZ}[rad]", 2);
1568  fitResults->setBinLabel(1, "Vtx[#]", 2);
1569  fitResults->setBinLabel(1, "Value", 1);
1570  fitResults->setBinLabel(2, "Error (stat)", 1);
1571  fitResults->getTH1()->SetOption("text");
1572 
1573  ibooker.setCurrentFolder("BeamPixel/EventInfo");
1574 
1575  reportSummary = ibooker.bookFloat("reportSummary");
1576  reportSummary->Fill(-1);
1577  reportSummaryMap = ibooker.book2D("reportSummaryMap", "Pixel-Vertices Beam Spot: % Good Fits", 1, 0., 1., 1, 0., 1.);
1578  reportSummaryMap->getTH1()->SetBinContent(1, 1, -1);
1579 
1580  ibooker.setCurrentFolder("BeamPixel/EventInfo/reportSummaryContents");
1581 
1582  // Convention for reportSummary and reportSummaryMap:
1583  // - -1% at the moment of creation of the histogram (i.e. white histogram)
1584  // - n% numberGoodFits / numberFits
1585 
1586  reset("scratch"); // Initialize histograms after creation
1587 }
1588 
1589 // Define this as a plug-in
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
MonitorElement * bookFloat(TString const &name, FUNC onbooking=NOOP())
Definition: DQMStore.h:80
virtual void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:36
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
int MyFit(std::vector< double > *vals)
void analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) override
void dqmBeginLuminosityBlock(const edm::LuminosityBlock &lumiBlock, const edm::EventSetup &iSetup) override
#define DIM(a)
void dqmEndLuminosityBlock(const edm::LuminosityBlock &lumiBlock, const edm::EventSetup &iSetup) override
data_type const * const_iterator
Definition: DetSetNew.h:31
Log< level::Error, false > LogError
unsigned int HitCounter(const edm::Event &iEvent)
double Gauss3DFunc(const double *par)
void writeToFile(std::vector< double > *vals, edm::TimeValue_t BeginTimeOfFit, edm::TimeValue_t EndTimeOfFit, unsigned int BeginLumiOfFit, unsigned int EndLumiOfFit, int dataType)
const_iterator end(bool update=false) const
const Double_t pi
T getUntrackedParameter(std::string const &, T const &) const
void Fill(long long x)
double Covariance[3][3]
Vx3DHLTAnalyzer(const edm::ParameterSet &)
int iEvent
Definition: GenABIO.cc:224
std::string formatTime(const time_t &t)
T sqrt(T t)
Definition: SSEVec.h:23
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
void reset(std::string ResetType)
unsigned long long TimeValue_t
Definition: Timestamp.h:21
Log< level::Info, false > LogInfo
boost::transform_iterator< IterHelp, const_IdIter > const_iterator
const_iterator begin(bool update=false) const
TimeValue_t value() const
Definition: Timestamp.h:38
Timestamp const & beginTime() const
Timestamp const & endTime() const
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:221
~Vx3DHLTAnalyzer() override
virtual TH1 * getTH1() const
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
fixed size matrix
HLT enums.
void printFitParams(const std::vector< double > &fitResults)
static std::atomic< unsigned int > counter
float x
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98
char const * what() const noexcept override
Definition: Exception.cc:107
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
tmp
align.sh
Definition: createJobs.py:716
LuminosityBlockNumber_t luminosityBlock() const
void reset(double vett[256])
Definition: TPedValues.cc:11
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
Definition: Run.h:45
virtual void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)