CMS 3D CMS Logo

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