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