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