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 //
3 // Package: Vx3DHLTAnalyzer
4 // Class: Vx3DHLTAnalyzer
5 //
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 
17 
21 
24 
25 #include <Math/Minimizer.h>
26 #include <Math/Factory.h>
27 #include <Math/Functor.h>
28 
29 
30 using namespace std;
31 using namespace reco;
32 using namespace edm;
33 
34 
36 {
37  debugMode = true;
38  nLumiReset = 1; // Number of integrated lumis to perform the fit
39  dataFromFit = true; // The Beam Spot data can be either taken from the histograms or from the fit results
40  minNentries = 35; // Minimum number of good vertices to perform the fit
41  xRange = 2.; // [cm]
42  xStep = 0.001; // [cm]
43  yRange = 2.; // [cm]
44  yStep = 0.001; // [cm]
45  zRange = 30.; // [cm]
46  zStep = 0.05; // [cm]
47  VxErrCorr = 1.5;
48  fileName = "BeamPixelResults.txt";
49 
50 
51  vertexCollection = consumes<reco::VertexCollection>(iConfig.getUntrackedParameter<edm::InputTag>("vertexCollection", edm::InputTag("pixelVertices")));
52  pixelHitCollection = consumes<SiPixelRecHitCollection>(iConfig.getUntrackedParameter<edm::InputTag>("pixelHitCollection", edm::InputTag("siPixelRecHits")));
53 
54  debugMode = iConfig.getParameter<bool>("debugMode");
55  nLumiReset = iConfig.getParameter<unsigned int>("nLumiReset");
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  fileName = iConfig.getParameter<string>("fileName");
66 }
67 
68 
70 {
71 }
72 
73 
74 void Vx3DHLTAnalyzer::analyze(const Event& iEvent, const EventSetup& iSetup)
75 {
77  iEvent.getByToken(vertexCollection, Vx3DCollection);
78 
79  unsigned int i,j;
80  double det;
81  VertexType MyVertex;
82 
83  if (runNumber != iEvent.id().run())
84  {
85  reset("scratch");
86  runNumber = iEvent.id().run();
87 
88  if (debugMode == true)
89  {
90  stringstream debugFile;
91  string tmp(fileName);
92 
93  if (outputDebugFile.is_open() == true) outputDebugFile.close();
94  tmp.erase(strlen(fileName.c_str())-4,4);
95  debugFile << tmp.c_str() << "_Run" << iEvent.id().run() << ".txt";
96  outputDebugFile.open(debugFile.str().c_str(), ios::out);
97  outputDebugFile.close();
98  outputDebugFile.open(debugFile.str().c_str(), ios::app);
99  }
100 
101  beginLuminosityBlock(iEvent.getLuminosityBlock(),iSetup);
102  }
103  else if (beginTimeOfFit != 0)
104  {
105  totalHits += HitCounter(iEvent);
106 
107  for (vector<Vertex>::const_iterator it3DVx = Vx3DCollection->begin(); it3DVx != Vx3DCollection->end(); it3DVx++)
108  {
109  if ((it3DVx->isValid() == true) &&
110  (it3DVx->isFake() == false) &&
111  (it3DVx->ndof() >= minVxDoF))
112  {
113  for (i = 0; i < DIM; i++)
114  {
115  for (j = 0; j < DIM; j++)
116  {
117  MyVertex.Covariance[i][j] = it3DVx->covariance(i,j);
118  if (edm::isNotFinite(MyVertex.Covariance[i][j]) == true) break;
119  }
120  if (j != DIM) break;
121  }
122  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]) -
123  MyVertex.Covariance[0][1]*(MyVertex.Covariance[0][1]*std::fabs(MyVertex.Covariance[2][2]) - MyVertex.Covariance[0][2]*MyVertex.Covariance[1][2]) +
124  MyVertex.Covariance[0][2]*(MyVertex.Covariance[0][1]*MyVertex.Covariance[1][2] - MyVertex.Covariance[0][2]*std::fabs(MyVertex.Covariance[1][1]));
125  if ((i == DIM) && (det > 0.))
126  {
127  MyVertex.x = it3DVx->x();
128  MyVertex.y = it3DVx->y();
129  MyVertex.z = it3DVx->z();
130  Vertices.push_back(MyVertex);
131  }
132  else if (internalDebug == true)
133  {
134  cout << "Vertex discarded !" << endl;
135  for (i = 0; i < DIM; i++)
136  for (j = 0; j < DIM; j++)
137  cout << "(i,j) --> " << i << "," << j << " --> " << MyVertex.Covariance[i][j] << endl;
138  }
139 
140  Vx_X->Fill(it3DVx->x());
141  Vx_Y->Fill(it3DVx->y());
142  Vx_Z->Fill(it3DVx->z());
143 
144  Vx_ZX->Fill(it3DVx->z(), it3DVx->x());
145  Vx_ZY->Fill(it3DVx->z(), it3DVx->y());
146  Vx_XY->Fill(it3DVx->x(), it3DVx->y());
147  }
148  }
149  }
150 }
151 
152 
154 {
156  iEvent.getByToken(pixelHitCollection, rechitspixel);
157 
158  unsigned int counter = 0;
159 
160  for (SiPixelRecHitCollection::const_iterator j = rechitspixel->begin(); j != rechitspixel->end(); j++)
161  for (edmNew::DetSet<SiPixelRecHit>::const_iterator h = j->begin(); h != j->end(); h++) counter += h->cluster()->size();
162 
163  return counter;
164 }
165 
166 
168 {
169  char ts[25];
170  strftime(ts, sizeof(ts), "%Y.%m.%d %H:%M:%S %Z", gmtime(&t));
171 
172  std::string ts_string(ts);
173 
174  return ts_string;
175 }
176 
177 
178 double Gauss3DFunc(const double* par)
179 {
180  double K[DIM][DIM]; // Covariance Matrix
181  double M[DIM][DIM]; // K^-1
182  double det;
183  double sumlog = 0.;
184 
185 // par[0] = K(0,0) --> Var[X]
186 // par[1] = K(1,1) --> Var[Y]
187 // par[2] = K(2,2) --> Var[Z]
188 // par[3] = K(0,1) = K(1,0) --> Cov[X,Y]
189 // par[4] = K(1,2) = K(2,1) --> Cov[Y,Z] --> dy/dz
190 // par[5] = K(0,2) = K(2,0) --> Cov[X,Z] --> dx/dz
191 // par[6] = mean x
192 // par[7] = mean y
193 // par[8] = mean z
194 
195  counterVx = 0;
196  for (unsigned int i = 0; i < Vertices.size(); i++)
197  {
199  (std::fabs(Vertices[i].z-zPos) <= maxLongLength))
200  {
201  if (considerVxCovariance == true)
202  {
203  K[0][0] = std::fabs(par[0]) + VxErrCorr*VxErrCorr * std::fabs(Vertices[i].Covariance[0][0]);
204  K[1][1] = std::fabs(par[1]) + VxErrCorr*VxErrCorr * std::fabs(Vertices[i].Covariance[1][1]);
205  K[2][2] = std::fabs(par[2]) + VxErrCorr*VxErrCorr * std::fabs(Vertices[i].Covariance[2][2]);
206  K[0][1] = K[1][0] = par[3] + VxErrCorr*VxErrCorr * Vertices[i].Covariance[0][1];
207  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];
208  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];
209  }
210  else
211  {
212  K[0][0] = std::fabs(par[0]);
213  K[1][1] = std::fabs(par[1]);
214  K[2][2] = std::fabs(par[2]);
215  K[0][1] = K[1][0] = par[3];
216  K[1][2] = K[2][1] = par[4]*(std::fabs(par[2])-std::fabs(par[1])) - par[5]*par[3];
217  K[0][2] = K[2][0] = par[5]*(std::fabs(par[2])-std::fabs(par[0])) - par[4]*par[3];
218  }
219 
220  det = K[0][0]*(K[1][1]*K[2][2] - K[1][2]*K[1][2]) -
221  K[0][1]*(K[0][1]*K[2][2] - K[0][2]*K[1][2]) +
222  K[0][2]*(K[0][1]*K[1][2] - K[0][2]*K[1][1]);
223 
224  M[0][0] = (K[1][1]*K[2][2] - K[1][2]*K[1][2]) / det;
225  M[1][1] = (K[0][0]*K[2][2] - K[0][2]*K[0][2]) / det;
226  M[2][2] = (K[0][0]*K[1][1] - K[0][1]*K[0][1]) / det;
227  M[0][1] = M[1][0] = (K[0][2]*K[1][2] - K[0][1]*K[2][2]) / det;
228  M[1][2] = M[2][1] = (K[0][2]*K[0][1] - K[1][2]*K[0][0]) / det;
229  M[0][2] = M[2][0] = (K[0][1]*K[1][2] - K[0][2]*K[1][1]) / det;
230 
231  sumlog += double(DIM)*std::log(2.*pi) + std::log(std::fabs(det)) +
232  (M[0][0]*(Vertices[i].x-par[6])*(Vertices[i].x-par[6]) +
233  M[1][1]*(Vertices[i].y-par[7])*(Vertices[i].y-par[7]) +
234  M[2][2]*(Vertices[i].z-par[8])*(Vertices[i].z-par[8]) +
235  2.*M[0][1]*(Vertices[i].x-par[6])*(Vertices[i].y-par[7]) +
236  2.*M[1][2]*(Vertices[i].y-par[7])*(Vertices[i].z-par[8]) +
237  2.*M[0][2]*(Vertices[i].x-par[6])*(Vertices[i].z-par[8]));
238 
239  counterVx++;
240  }
241  }
242 
243  return sumlog;
244 }
245 
246 
247 int Vx3DHLTAnalyzer::MyFit(vector<double>* vals)
248 {
249  // RETURN CODE:
250  // 0 == OK
251  // -2 == NO OK - not enough "minNentries"
252  // Any other number == NO OK
253  unsigned int nParams = 9;
254 
255  if ((vals != NULL) && (vals->size() == nParams*2))
256  {
257  double nSigmaXY = 100.;
258  double nSigmaZ = 100.;
259  double varFactor = 4./25.; // Take into account the difference between the RMS and sigma (RMS usually greater than sigma)
260  double parDistanceXY = 0.005; // Unit: [cm]
261  double parDistanceZ = 0.5; // Unit: [cm]
262  double parDistanceddZ = 1e-3; // Unit: [rad]
263  double parDistanceCxy = 1e-5; // Unit: [cm^2]
264  double bestEdm = 1e-1;
265 
266  const unsigned int trials = 4;
267  double largerDist[trials] = {0.1, 5., 10., 100.};
268 
269  double covxz,covyz,det;
270  double deltaMean;
271  int bestMovementX = 1;
272  int bestMovementY = 1;
273  int bestMovementZ = 1;
274  int goodData;
275 
276  double edm;
277 
278  vector<double>::const_iterator it = vals->begin();
279 
280  ROOT::Math::Minimizer* Gauss3D = ROOT::Math::Factory::CreateMinimizer("Minuit2","Migrad");
281  Gauss3D->SetMaxFunctionCalls(1e4);
282  Gauss3D->SetTolerance(1e-9); // Tolerance on likelihood
283  if (internalDebug == true) Gauss3D->SetPrintLevel(3);
284  else Gauss3D->SetPrintLevel(0);
285 
286  ROOT::Math::Functor _Gauss3DFunc(&Gauss3DFunc,nParams);
287  Gauss3D->SetFunction(_Gauss3DFunc);
288 
289  if (internalDebug == true) cout << "\n@@@ START FITTING @@@" << endl;
290 
291  // @@@ Fit at X-deltaMean | X | X+deltaMean @@@
292  bestEdm = 1.;
293  for (int i = 0; i < 3; i++)
294  {
295  deltaMean = (double(i)-1.)*std::sqrt((*(it+0))*varFactor);
296  if (internalDebug == true) cout << "deltaMean --> " << deltaMean << endl;
297 
298  Gauss3D->Clear();
299 
300  Gauss3D->SetVariable(0,"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY);
301  Gauss3D->SetVariable(1,"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY);
302  Gauss3D->SetVariable(2,"var z ", *(it+2), parDistanceZ*parDistanceZ);
303  Gauss3D->SetVariable(3,"cov xy", *(it+3), parDistanceCxy);
304  Gauss3D->SetVariable(4,"dydz ", *(it+4), parDistanceddZ);
305  Gauss3D->SetVariable(5,"dxdz ", *(it+5), parDistanceddZ);
306  Gauss3D->SetVariable(6,"mean x", *(it+6)+deltaMean, parDistanceXY);
307  Gauss3D->SetVariable(7,"mean y", *(it+7), parDistanceXY);
308  Gauss3D->SetVariable(8,"mean z", *(it+8), parDistanceZ);
309 
310  // Set the central positions of the centroid for vertex rejection
311  xPos = Gauss3D->X()[6];
312  yPos = Gauss3D->X()[7];
313  zPos = Gauss3D->X()[8];
314 
315  // Set dimensions of the centroid for vertex rejection
316  maxTransRadius = nSigmaXY * std::sqrt(std::fabs(Gauss3D->X()[0]) + std::fabs(Gauss3D->X()[1])) / 2.;
317  maxLongLength = nSigmaZ * std::sqrt(std::fabs(Gauss3D->X()[2]));
318 
319  Gauss3D->Minimize();
320  goodData = Gauss3D->Status();
321  edm = Gauss3D->Edm();
322 
323  if (counterVx < minNentries) goodData = -2;
324  else if (edm::isNotFinite(edm) == true) goodData = -1;
325  else for (unsigned int j = 0; j < nParams; j++) if (edm::isNotFinite(Gauss3D->Errors()[j]) == true) { goodData = -1; break; }
326  if (goodData == 0)
327  {
328  covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
329  covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
330 
331  det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
332  Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
333  covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
334  if (det < 0.) { goodData = -1; if (internalDebug == true) cout << "Negative determinant !" << endl; }
335  }
336 
337  if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementX = i; }
338  }
339  if (internalDebug == true) cout << "Found bestMovementX --> " << bestMovementX << endl;
340 
341  // @@@ Fit at Y-deltaMean | Y | Y+deltaMean @@@
342  bestEdm = 1.;
343  for (int i = 0; i < 3; i++)
344  {
345  deltaMean = (double(i)-1.)*std::sqrt((*(it+1))*varFactor);
346  if (internalDebug == true)
347  {
348  cout << "deltaMean --> " << deltaMean << endl;
349  cout << "deltaMean X --> " << (double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor) << endl;
350  }
351 
352  Gauss3D->Clear();
353 
354  Gauss3D->SetVariable(0,"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY);
355  Gauss3D->SetVariable(1,"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY);
356  Gauss3D->SetVariable(2,"var z ", *(it+2), parDistanceZ*parDistanceZ);
357  Gauss3D->SetVariable(3,"cov xy", *(it+3), parDistanceCxy);
358  Gauss3D->SetVariable(4,"dydz ", *(it+4), parDistanceddZ);
359  Gauss3D->SetVariable(5,"dxdz ", *(it+5), parDistanceddZ);
360  Gauss3D->SetVariable(6,"mean x", *(it+6)+(double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor), parDistanceXY);
361  Gauss3D->SetVariable(7,"mean y", *(it+7)+deltaMean, parDistanceXY);
362  Gauss3D->SetVariable(8,"mean z", *(it+8), parDistanceZ);
363 
364  // Set the central positions of the centroid for vertex rejection
365  xPos = Gauss3D->X()[6];
366  yPos = Gauss3D->X()[7];
367  zPos = Gauss3D->X()[8];
368 
369  // Set dimensions of the centroid for vertex rejection
370  maxTransRadius = nSigmaXY * std::sqrt(std::fabs(Gauss3D->X()[0]) + std::fabs(Gauss3D->X()[1])) / 2.;
371  maxLongLength = nSigmaZ * std::sqrt(std::fabs(Gauss3D->X()[2]));
372 
373  Gauss3D->Minimize();
374  goodData = Gauss3D->Status();
375  edm = Gauss3D->Edm();
376 
377  if (counterVx < minNentries) goodData = -2;
378  else if (edm::isNotFinite(edm) == true) goodData = -1;
379  else for (unsigned int j = 0; j < nParams; j++) if (edm::isNotFinite(Gauss3D->Errors()[j]) == true) { goodData = -1; break; }
380  if (goodData == 0)
381  {
382  covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
383  covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
384 
385  det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
386  Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
387  covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
388  if (det < 0.) { goodData = -1; if (internalDebug == true) cout << "Negative determinant !" << endl; }
389  }
390 
391  if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementY = i; }
392  }
393  if (internalDebug == true) cout << "Found bestMovementY --> " << bestMovementY << endl;
394 
395  // @@@ Fit at Z-deltaMean | Z | Z+deltaMean @@@
396  bestEdm = 1.;
397  for (int i = 0; i < 3; i++)
398  {
399  deltaMean = (double(i)-1.)*std::sqrt(*(it+2));
400  if (internalDebug == true)
401  {
402  cout << "deltaMean --> " << deltaMean << endl;
403  cout << "deltaMean X --> " << (double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor) << endl;
404  cout << "deltaMean Y --> " << (double(bestMovementY)-1.)*std::sqrt((*(it+1))*varFactor) << endl;
405  }
406 
407  Gauss3D->Clear();
408 
409  Gauss3D->SetVariable(0,"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY);
410  Gauss3D->SetVariable(1,"var y ", *(it+1)*varFactor, 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))*varFactor), parDistanceXY);
416  Gauss3D->SetVariable(7,"mean y", *(it+7)+(double(bestMovementY)-1.)*std::sqrt((*(it+1))*varFactor), parDistanceXY);
417  Gauss3D->SetVariable(8,"mean z", *(it+8)+deltaMean, 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 (edm::isNotFinite(edm) == true) goodData = -1;
434  else for (unsigned int j = 0; j < nParams; j++) if (edm::isNotFinite(Gauss3D->Errors()[j]) == true) { goodData = -1; break; }
435  if (goodData == 0)
436  {
437  covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
438  covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
439 
440  det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
441  Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
442  covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
443  if (det < 0.) { goodData = -1; if (internalDebug == true) cout << "Negative determinant !" << endl; }
444  }
445 
446  if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementZ = i; }
447  }
448  if (internalDebug == true) cout << "Found bestMovementZ --> " << bestMovementZ << endl;
449 
450  Gauss3D->Clear();
451 
452  // @@@ FINAL FIT @@@
453  Gauss3D->SetVariable(0,"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY);
454  Gauss3D->SetVariable(1,"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY);
455  Gauss3D->SetVariable(2,"var z ", *(it+2), parDistanceZ*parDistanceZ);
456  Gauss3D->SetVariable(3,"cov xy", *(it+3), parDistanceCxy);
457  Gauss3D->SetVariable(4,"dydz ", *(it+4), parDistanceddZ);
458  Gauss3D->SetVariable(5,"dxdz ", *(it+5), parDistanceddZ);
459  Gauss3D->SetVariable(6,"mean x", *(it+6)+(double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor), parDistanceXY);
460  Gauss3D->SetVariable(7,"mean y", *(it+7)+(double(bestMovementY)-1.)*std::sqrt((*(it+1))*varFactor), parDistanceXY);
461  Gauss3D->SetVariable(8,"mean z", *(it+8)+(double(bestMovementZ)-1.)*std::sqrt(*(it+2)), parDistanceZ);
462 
463  // Set the central positions of the centroid for vertex rejection
464  xPos = Gauss3D->X()[6];
465  yPos = Gauss3D->X()[7];
466  zPos = Gauss3D->X()[8];
467 
468  // Set dimensions of the centroid for vertex rejection
469  maxTransRadius = nSigmaXY * std::sqrt(std::fabs(Gauss3D->X()[0]) + std::fabs(Gauss3D->X()[1])) / 2.;
470  maxLongLength = nSigmaZ * std::sqrt(std::fabs(Gauss3D->X()[2]));
471 
472  Gauss3D->Minimize();
473  goodData = Gauss3D->Status();
474  edm = Gauss3D->Edm();
475 
476  if (counterVx < minNentries) goodData = -2;
477  else if (edm::isNotFinite(edm) == true) goodData = -1;
478  else for (unsigned int j = 0; j < nParams; j++) if (edm::isNotFinite(Gauss3D->Errors()[j]) == true) { goodData = -1; break; }
479  if (goodData == 0)
480  {
481  covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
482  covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
483 
484  det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
485  Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
486  covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
487  if (det < 0.) { goodData = -1; if (internalDebug == true) cout << "Negative determinant !" << endl; }
488  }
489 
490  // @@@ FIT WITH DIFFERENT PARAMETER DISTANCES @@@
491  for (unsigned int i = 0; i < trials; i++)
492  {
493  if ((goodData != 0) && (goodData != -2))
494  {
495  Gauss3D->Clear();
496 
497  if (internalDebug == true) cout << "FIT WITH DIFFERENT PARAMETER DISTANCES - STEP " << i+1 << endl;
498 
499  Gauss3D->SetVariable(0,"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY * largerDist[i]);
500  Gauss3D->SetVariable(1,"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY * largerDist[i]);
501  Gauss3D->SetVariable(2,"var z ", *(it+2), parDistanceZ*parDistanceZ * largerDist[i]);
502  Gauss3D->SetVariable(3,"cov xy", *(it+3), parDistanceCxy * largerDist[i]);
503  Gauss3D->SetVariable(4,"dydz ", *(it+4), parDistanceddZ * largerDist[i]);
504  Gauss3D->SetVariable(5,"dxdz ", *(it+5), parDistanceddZ * largerDist[i]);
505  Gauss3D->SetVariable(6,"mean x", *(it+6)+(double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor), parDistanceXY * largerDist[i]);
506  Gauss3D->SetVariable(7,"mean y", *(it+7)+(double(bestMovementY)-1.)*std::sqrt((*(it+1))*varFactor), parDistanceXY * largerDist[i]);
507  Gauss3D->SetVariable(8,"mean z", *(it+8)+(double(bestMovementZ)-1.)*std::sqrt(*(it+2)), parDistanceZ * largerDist[i]);
508 
509  // Set the central positions of the centroid for vertex rejection
510  xPos = Gauss3D->X()[6];
511  yPos = Gauss3D->X()[7];
512  zPos = Gauss3D->X()[8];
513 
514  // Set dimensions of the centroid for vertex rejection
515  maxTransRadius = nSigmaXY * std::sqrt(std::fabs(Gauss3D->X()[0]) + std::fabs(Gauss3D->X()[1])) / 2.;
516  maxLongLength = nSigmaZ * std::sqrt(std::fabs(Gauss3D->X()[2]));
517 
518  Gauss3D->Minimize();
519  goodData = Gauss3D->Status();
520  edm = Gauss3D->Edm();
521 
522  if (counterVx < minNentries) goodData = -2;
523  else if (edm::isNotFinite(edm) == true) goodData = -1;
524  else for (unsigned int j = 0; j < nParams; j++) if (edm::isNotFinite(Gauss3D->Errors()[j]) == true) { goodData = -1; break; }
525  if (goodData == 0)
526  {
527  covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
528  covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
529 
530  det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
531  Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
532  covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
533  if (det < 0.) { goodData = -1; if (internalDebug == true) cout << "Negative determinant !" << endl; }
534  }
535  } else break;
536  }
537 
538  if (goodData == 0)
539  for (unsigned int i = 0; i < nParams; i++)
540  {
541  vals->operator[](i) = Gauss3D->X()[i];
542  vals->operator[](i+nParams) = Gauss3D->Errors()[i];
543  }
544 
545  delete Gauss3D;
546  return goodData;
547  }
548 
549  return -1;
550 }
551 
552 
553 void Vx3DHLTAnalyzer::reset(string ResetType)
554 {
555  if (ResetType.compare("scratch") == 0)
556  {
557  runNumber = 0;
558  numberGoodFits = 0;
559  numberFits = 0;
560  lastLumiOfFit = 0;
561 
562  Vx_X->Reset();
563  Vx_Y->Reset();
564  Vx_Z->Reset();
565 
566  Vx_ZX->Reset();
567  Vx_ZY->Reset();
568  Vx_XY->Reset();
569 
570  mXlumi->Reset();
571  mYlumi->Reset();
572  mZlumi->Reset();
573 
574  sXlumi->Reset();
575  sYlumi->Reset();
576  sZlumi->Reset();
577 
578  dxdzlumi->Reset();
579  dydzlumi->Reset();
580 
581  hitCounter->Reset();
582  hitCountHistory->Reset();
583  goodVxCounter->Reset();
584  goodVxCountHistory->Reset();
585  fitResults->Reset();
586 
587  reportSummary->Fill(0.);
588  reportSummaryMap->Fill(0.5, 0.5, 0.);
589 
590  Vertices.clear();
591 
592  lumiCounter = 0;
593  lumiCounterHisto = 0;
594  totalHits = 0;
595  beginTimeOfFit = 0;
596  endTimeOfFit = 0;
597  beginLumiOfFit = 0;
598  endLumiOfFit = 0;
599  }
600  else if (ResetType.compare("whole") == 0)
601  {
602  Vx_X->Reset();
603  Vx_Y->Reset();
604  Vx_Z->Reset();
605 
606  Vx_ZX->Reset();
607  Vx_ZY->Reset();
608  Vx_XY->Reset();
609 
610  Vertices.clear();
611 
612  lumiCounter = 0;
613  lumiCounterHisto = 0;
614  totalHits = 0;
615  beginTimeOfFit = 0;
616  endTimeOfFit = 0;
617  beginLumiOfFit = 0;
618  endLumiOfFit = 0;
619  }
620  else if (ResetType.compare("partial") == 0)
621  {
622  Vx_X->Reset();
623  Vx_Y->Reset();
624  Vx_Z->Reset();
625 
626  Vertices.clear();
627 
628  lumiCounter = 0;
629  totalHits = 0;
630  beginTimeOfFit = 0;
631  endTimeOfFit = 0;
632  beginLumiOfFit = 0;
633  endLumiOfFit = 0;
634  }
635  else if (ResetType.compare("nohisto") == 0)
636  {
637  Vertices.clear();
638 
639  lumiCounter = 0;
640  lumiCounterHisto = 0;
641  totalHits = 0;
642  beginTimeOfFit = 0;
643  endTimeOfFit = 0;
644  beginLumiOfFit = 0;
645  endLumiOfFit = 0;
646  }
647  else if (ResetType.compare("hitCounter") == 0)
648  totalHits = 0;
649 }
650 
651 
652 void Vx3DHLTAnalyzer::writeToFile(vector<double>* vals,
653  edm::TimeValue_t BeginTimeOfFit,
654  edm::TimeValue_t EndTimeOfFit,
655  unsigned int BeginLumiOfFit,
656  unsigned int EndLumiOfFit,
657  int dataType)
658 {
659  stringstream BufferString;
660  BufferString.precision(5);
661 
662  outputFile.open(fileName.c_str(), ios::out);
663 
664  if ((outputFile.is_open() == true) && (vals != NULL) && (vals->size() == 8*2))
665  {
666  vector<double>::const_iterator it = vals->begin();
667 
668  outputFile << "Runnumber " << runNumber << endl;
669  outputFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
670  outputFile << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " " << (endTimeOfFit >> 32) << endl;
671  outputFile << "LumiRange " << beginLumiOfFit << " - " << endLumiOfFit << endl;
672  outputFile << "Type " << dataType << endl;
673  // 3D Vertexing with Pixel Tracks:
674  // Good data = Type 3
675  // Bad data = Type -1
676 
677  BufferString << *(it+0);
678  outputFile << "X0 " << BufferString.str().c_str() << endl;
679  BufferString.str("");
680 
681  BufferString << *(it+1);
682  outputFile << "Y0 " << BufferString.str().c_str() << endl;
683  BufferString.str("");
684 
685  BufferString << *(it+2);
686  outputFile << "Z0 " << BufferString.str().c_str() << endl;
687  BufferString.str("");
688 
689  BufferString << *(it+3);
690  outputFile << "sigmaZ0 " << BufferString.str().c_str() << endl;
691  BufferString.str("");
692 
693  BufferString << *(it+4);
694  outputFile << "dxdz " << BufferString.str().c_str() << endl;
695  BufferString.str("");
696 
697  BufferString << *(it+5);
698  outputFile << "dydz " << BufferString.str().c_str() << endl;
699  BufferString.str("");
700 
701  BufferString << *(it+6);
702  outputFile << "BeamWidthX " << BufferString.str().c_str() << endl;
703  BufferString.str("");
704 
705  BufferString << *(it+7);
706  outputFile << "BeamWidthY " << BufferString.str().c_str() << endl;
707  BufferString.str("");
708 
709  outputFile << "Cov(0,j) " << *(it+8) << " 0.0 0.0 0.0 0.0 0.0 0.0" << endl;
710  outputFile << "Cov(1,j) 0.0 " << *(it+9) << " 0.0 0.0 0.0 0.0 0.0" << endl;
711  outputFile << "Cov(2,j) 0.0 0.0 " << *(it+10) << " 0.0 0.0 0.0 0.0" << endl;
712  outputFile << "Cov(3,j) 0.0 0.0 0.0 " << *(it+11) << " 0.0 0.0 0.0" << endl;
713  outputFile << "Cov(4,j) 0.0 0.0 0.0 0.0 " << *(it+12) << " 0.0 0.0" << endl;
714  outputFile << "Cov(5,j) 0.0 0.0 0.0 0.0 0.0 " << *(it+13) << " 0.0" << endl;
715  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;
716 
717  outputFile << "EmittanceX 0.0" << endl;
718  outputFile << "EmittanceY 0.0" << endl;
719  outputFile << "BetaStar 0.0" << endl;
720  }
721  outputFile.close();
722 
723  if ((debugMode == true) && (outputDebugFile.is_open() == true) && (vals != NULL) && (vals->size() == 8*2))
724  {
725  vector<double>::const_iterator it = vals->begin();
726 
727  outputDebugFile << "Runnumber " << runNumber << endl;
728  outputDebugFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
729  outputDebugFile << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " " << (endTimeOfFit >> 32) << endl;
730  outputDebugFile << "LumiRange " << beginLumiOfFit << " - " << endLumiOfFit << endl;
731  outputDebugFile << "Type " << dataType << endl;
732  // 3D Vertexing with Pixel Tracks:
733  // Good data = Type 3
734  // Bad data = Type -1
735 
736  BufferString << *(it+0);
737  outputDebugFile << "X0 " << BufferString.str().c_str() << endl;
738  BufferString.str("");
739 
740  BufferString << *(it+1);
741  outputDebugFile << "Y0 " << BufferString.str().c_str() << endl;
742  BufferString.str("");
743 
744  BufferString << *(it+2);
745  outputDebugFile << "Z0 " << BufferString.str().c_str() << endl;
746  BufferString.str("");
747 
748  BufferString << *(it+3);
749  outputDebugFile << "sigmaZ0 " << BufferString.str().c_str() << endl;
750  BufferString.str("");
751 
752  BufferString << *(it+4);
753  outputDebugFile << "dxdz " << BufferString.str().c_str() << endl;
754  BufferString.str("");
755 
756  BufferString << *(it+5);
757  outputDebugFile << "dydz " << BufferString.str().c_str() << endl;
758  BufferString.str("");
759 
760  BufferString << *(it+6);
761  outputDebugFile << "BeamWidthX " << BufferString.str().c_str() << endl;
762  BufferString.str("");
763 
764  BufferString << *(it+7);
765  outputDebugFile << "BeamWidthY " << BufferString.str().c_str() << endl;
766  BufferString.str("");
767 
768  outputDebugFile << "Cov(0,j) " << *(it+8) << " 0.0 0.0 0.0 0.0 0.0 0.0" << endl;
769  outputDebugFile << "Cov(1,j) 0.0 " << *(it+9) << " 0.0 0.0 0.0 0.0 0.0" << endl;
770  outputDebugFile << "Cov(2,j) 0.0 0.0 " << *(it+10) << " 0.0 0.0 0.0 0.0" << endl;
771  outputDebugFile << "Cov(3,j) 0.0 0.0 0.0 " << *(it+11) << " 0.0 0.0 0.0" << endl;
772  outputDebugFile << "Cov(4,j) 0.0 0.0 0.0 0.0 " << *(it+12) << " 0.0 0.0" << endl;
773  outputDebugFile << "Cov(5,j) 0.0 0.0 0.0 0.0 0.0 " << *(it+13) << " 0.0" << endl;
774  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;
775 
776  outputDebugFile << "EmittanceX 0.0" << endl;
777  outputDebugFile << "EmittanceY 0.0" << endl;
778  outputDebugFile << "BetaStar 0.0" << endl;
779  }
780 }
781 
782 
784  const EventSetup& iSetup)
785 {
786  if ((lumiCounter == 0) && (lumiBlock.luminosityBlock() > lastLumiOfFit))
787  {
788  beginTimeOfFit = lumiBlock.beginTime().value();
789  beginLumiOfFit = lumiBlock.luminosityBlock();
790  lumiCounter++;
791  lumiCounterHisto++;
792  }
793  else if ((lumiCounter != 0) && (lumiBlock.luminosityBlock() >= (beginLumiOfFit+lumiCounter))) { lumiCounter++; lumiCounterHisto++; }
794 }
795 
796 
798  const EventSetup& iSetup)
799 {
800  stringstream histTitle;
801  int goodData;
802  unsigned int nParams = 9;
803 
804  if ((lumiCounter%nLumiReset == 0) && (nLumiReset != 0) && (beginTimeOfFit != 0) && (runNumber != 0))
805  {
806  endTimeOfFit = lumiBlock.endTime().value();
807  endLumiOfFit = lumiBlock.luminosityBlock();
808  lastLumiOfFit = endLumiOfFit;
809  vector<double> vals;
810 
811  hitCounter->ShiftFillLast((double)totalHits, std::sqrt((double)totalHits), nLumiReset);
812 
813  if (lastLumiOfFit % prescaleHistory == 0)
814  {
815  hitCountHistory->getTH1()->SetBinContent(lastLumiOfFit, (double)totalHits);
816  hitCountHistory->getTH1()->SetBinError(lastLumiOfFit, std::sqrt((double)totalHits));
817  }
818 
819  if (dataFromFit == true)
820  {
821  vector<double> fitResults;
822 
823  fitResults.push_back(Vx_X->getTH1()->GetRMS()*Vx_X->getTH1()->GetRMS());
824  fitResults.push_back(Vx_Y->getTH1()->GetRMS()*Vx_Y->getTH1()->GetRMS());
825  fitResults.push_back(Vx_Z->getTH1()->GetRMS()*Vx_Z->getTH1()->GetRMS());
826  fitResults.push_back(0.0);
827  fitResults.push_back(0.0);
828  fitResults.push_back(0.0);
829  fitResults.push_back(Vx_X->getTH1()->GetMean());
830  fitResults.push_back(Vx_Y->getTH1()->GetMean());
831  fitResults.push_back(Vx_Z->getTH1()->GetMean());
832  for (unsigned int i = 0; i < nParams; i++) fitResults.push_back(0.0);
833 
834  goodData = MyFit(&fitResults);
835 
836  if (internalDebug == true)
837  {
838  cout << "goodData --> " << goodData << endl;
839  cout << "Used vertices --> " << counterVx << endl;
840  cout << "var x --> " << fitResults[0] << " +/- " << fitResults[0+nParams] << endl;
841  cout << "var y --> " << fitResults[1] << " +/- " << fitResults[1+nParams] << endl;
842  cout << "var z --> " << fitResults[2] << " +/- " << fitResults[2+nParams] << endl;
843  cout << "cov xy --> " << fitResults[3] << " +/- " << fitResults[3+nParams] << endl;
844  cout << "dydz --> " << fitResults[4] << " +/- " << fitResults[4+nParams] << endl;
845  cout << "dxdz --> " << fitResults[5] << " +/- " << fitResults[5+nParams] << endl;
846  cout << "mean x --> " << fitResults[6] << " +/- " << fitResults[6+nParams] << endl;
847  cout << "mean y --> " << fitResults[7] << " +/- " << fitResults[7+nParams] << endl;
848  cout << "mean z --> " << fitResults[8] << " +/- " << fitResults[8+nParams] << endl;
849  }
850 
851  if (goodData == 0)
852  {
853  vals.push_back(fitResults[6]);
854  vals.push_back(fitResults[7]);
855  vals.push_back(fitResults[8]);
856  vals.push_back(std::sqrt(std::fabs(fitResults[2])));
857  vals.push_back(fitResults[5]);
858  vals.push_back(fitResults[4]);
859  vals.push_back(std::sqrt(std::fabs(fitResults[0])));
860  vals.push_back(std::sqrt(std::fabs(fitResults[1])));
861 
862  vals.push_back(std::pow(fitResults[6+nParams],2.));
863  vals.push_back(std::pow(fitResults[7+nParams],2.));
864  vals.push_back(std::pow(fitResults[8+nParams],2.));
865  vals.push_back(std::pow(std::fabs(fitResults[2+nParams]) / (2.*std::sqrt(std::fabs(fitResults[2]))),2.));
866  vals.push_back(std::pow(fitResults[5+nParams],2.));
867  vals.push_back(std::pow(fitResults[4+nParams],2.));
868  vals.push_back(std::pow(std::fabs(fitResults[0+nParams]) / (2.*std::sqrt(std::fabs(fitResults[0]))),2.));
869  vals.push_back(std::pow(std::fabs(fitResults[1+nParams]) / (2.*std::sqrt(std::fabs(fitResults[1]))),2.));
870  }
871  else for (unsigned int i = 0; i < 8*2; i++) vals.push_back(0.0);
872 
873  fitResults.clear();
874  }
875  else
876  {
877  counterVx = Vx_X->getTH1F()->GetEntries();
878 
879  if (Vx_X->getTH1F()->GetEntries() >= minNentries)
880  {
881  goodData = 0;
882 
883  vals.push_back(Vx_X->getTH1F()->GetMean());
884  vals.push_back(Vx_Y->getTH1F()->GetMean());
885  vals.push_back(Vx_Z->getTH1F()->GetMean());
886  vals.push_back(Vx_Z->getTH1F()->GetRMS());
887  vals.push_back(0.0);
888  vals.push_back(0.0);
889  vals.push_back(Vx_X->getTH1F()->GetRMS());
890  vals.push_back(Vx_Y->getTH1F()->GetRMS());
891 
892  vals.push_back(std::pow(Vx_X->getTH1F()->GetMeanError(),2.));
893  vals.push_back(std::pow(Vx_Y->getTH1F()->GetMeanError(),2.));
894  vals.push_back(std::pow(Vx_Z->getTH1F()->GetMeanError(),2.));
895  vals.push_back(std::pow(Vx_Z->getTH1F()->GetRMSError(),2.));
896  vals.push_back(0.0);
897  vals.push_back(0.0);
898  vals.push_back(std::pow(Vx_X->getTH1F()->GetRMSError(),2.));
899  vals.push_back(std::pow(Vx_Y->getTH1F()->GetRMSError(),2.));
900  }
901  else
902  {
903  goodData = -2;
904  for (unsigned int i = 0; i < 8*2; i++) vals.push_back(0.0);
905  }
906  }
907 
908  // vals[0] = X0
909  // vals[1] = Y0
910  // vals[2] = Z0
911  // vals[3] = sigmaZ0
912  // vals[4] = dxdz
913  // vals[5] = dydz
914  // vals[6] = BeamWidthX
915  // vals[7] = BeamWidthY
916 
917  // vals[8] = err^2 X0
918  // vals[9] = err^2 Y0
919  // vals[10] = err^2 Z0
920  // vals[11] = err^2 sigmaZ0
921  // vals[12] = err^2 dxdz
922  // vals[13] = err^2 dydz
923  // vals[14] = err^2 BeamWidthX
924  // vals[15] = err^2 BeamWidthY
925 
926  // "goodData" CODE:
927  // 0 == OK --> Reset
928  // -2 == NO OK - not enough "minNentries" --> Wait for more lumisections
929  // Any other number == NO OK --> Reset
930 
931  numberFits++;
932  if (goodData == 0)
933  {
934  writeToFile(&vals, beginTimeOfFit, endTimeOfFit, beginLumiOfFit, endLumiOfFit, 3);
935  if ((internalDebug == true) && (outputDebugFile.is_open() == true)) outputDebugFile << "Used vertices: " << counterVx << endl;
936 
937  numberGoodFits++;
938 
939  histTitle << "Fitted Beam Spot [cm] (Lumi start: " << beginLumiOfFit << " - Lumi end: " << endLumiOfFit << ")";
940  if (lumiCounterHisto >= maxLumiIntegration) reset("whole");
941  else reset("partial");
942  }
943  else
944  {
945  writeToFile(&vals, beginTimeOfFit, endTimeOfFit, beginLumiOfFit, endLumiOfFit, -1);
946  if ((internalDebug == true) && (outputDebugFile.is_open() == true)) outputDebugFile << "Used vertices: " << counterVx << endl;
947 
948  if (goodData == -2)
949  {
950  histTitle << "Fitted Beam Spot [cm] (not enough statistics)";
951  if (lumiCounter >= maxLumiIntegration) reset("whole");
952  else reset("hitCounter");
953  }
954  else
955  {
956  histTitle << "Fitted Beam Spot [cm] (problems)";
957  if (lumiCounterHisto >= maxLumiIntegration) reset("whole");
958  else reset("partial");
959 
960  counterVx = 0;
961  }
962  }
963 
964  reportSummary->Fill(numberFits != 0 ? (double)numberGoodFits/(double)numberFits : 0.0);
965  reportSummaryMap->Fill(0.5, 0.5, numberFits != 0 ? (double)numberGoodFits/(double)numberFits : 0.0);
966 
967  fitResults->setAxisTitle(histTitle.str().c_str(), 1);
968 
969  fitResults->setBinContent(1, 9, vals[0]);
970  fitResults->setBinContent(1, 8, vals[1]);
971  fitResults->setBinContent(1, 7, vals[2]);
972  fitResults->setBinContent(1, 6, vals[3]);
973  fitResults->setBinContent(1, 5, vals[4]);
974  fitResults->setBinContent(1, 4, vals[5]);
975  fitResults->setBinContent(1, 3, vals[6]);
976  fitResults->setBinContent(1, 2, vals[7]);
977  fitResults->setBinContent(1, 1, counterVx);
978 
979  fitResults->setBinContent(2, 9, std::sqrt(vals[8]));
980  fitResults->setBinContent(2, 8, std::sqrt(vals[9]));
981  fitResults->setBinContent(2, 7, std::sqrt(vals[10]));
982  fitResults->setBinContent(2, 6, std::sqrt(vals[11]));
983  fitResults->setBinContent(2, 5, std::sqrt(vals[12]));
984  fitResults->setBinContent(2, 4, std::sqrt(vals[13]));
985  fitResults->setBinContent(2, 3, std::sqrt(vals[14]));
986  fitResults->setBinContent(2, 2, std::sqrt(vals[15]));
987  fitResults->setBinContent(2, 1, std::sqrt(counterVx));
988 
989  // Linear fit to the historical plots
990  TF1* myLinFit = new TF1("myLinFit", "[0] + [1]*x", mXlumi->getTH1()->GetXaxis()->GetXmin(), mXlumi->getTH1()->GetXaxis()->GetXmax());
991  myLinFit->SetLineColor(2);
992  myLinFit->SetLineWidth(2);
993  myLinFit->SetParName(0,"Intercept");
994  myLinFit->SetParName(1,"Slope");
995 
996  mXlumi->ShiftFillLast(vals[0], std::sqrt(vals[8]), nLumiReset);
997  myLinFit->SetParameter(0, mXlumi->getTH1()->GetMean(2));
998  myLinFit->SetParameter(1, 0.0);
999  mXlumi->getTH1()->Fit(myLinFit,"QR");
1000 
1001  mYlumi->ShiftFillLast(vals[1], std::sqrt(vals[9]), nLumiReset);
1002  myLinFit->SetParameter(0, mYlumi->getTH1()->GetMean(2));
1003  myLinFit->SetParameter(1, 0.0);
1004  mYlumi->getTH1()->Fit(myLinFit,"QR");
1005 
1006  mZlumi->ShiftFillLast(vals[2], std::sqrt(vals[10]), nLumiReset);
1007  myLinFit->SetParameter(0, mZlumi->getTH1()->GetMean(2));
1008  myLinFit->SetParameter(1, 0.0);
1009  mZlumi->getTH1()->Fit(myLinFit,"QR");
1010 
1011  sXlumi->ShiftFillLast(vals[6], std::sqrt(vals[14]), nLumiReset);
1012  myLinFit->SetParameter(0, sXlumi->getTH1()->GetMean(2));
1013  myLinFit->SetParameter(1, 0.0);
1014  sXlumi->getTH1()->Fit(myLinFit,"QR");
1015 
1016  sYlumi->ShiftFillLast(vals[7], std::sqrt(vals[15]), nLumiReset);
1017  myLinFit->SetParameter(0, sYlumi->getTH1()->GetMean(2));
1018  myLinFit->SetParameter(1, 0.0);
1019  sYlumi->getTH1()->Fit(myLinFit,"QR");
1020 
1021  sZlumi->ShiftFillLast(vals[3], std::sqrt(vals[11]), nLumiReset);
1022  myLinFit->SetParameter(0, sZlumi->getTH1()->GetMean(2));
1023  myLinFit->SetParameter(1, 0.0);
1024  sZlumi->getTH1()->Fit(myLinFit,"QR");
1025 
1026  dxdzlumi->ShiftFillLast(vals[4], std::sqrt(vals[12]), nLumiReset);
1027  myLinFit->SetParameter(0, dxdzlumi->getTH1()->GetMean(2));
1028  myLinFit->SetParameter(1, 0.0);
1029  dxdzlumi->getTH1()->Fit(myLinFit,"QR");
1030 
1031  dydzlumi->ShiftFillLast(vals[5], std::sqrt(vals[13]), nLumiReset);
1032  myLinFit->SetParameter(0, dydzlumi->getTH1()->GetMean(2));
1033  myLinFit->SetParameter(1, 0.0);
1034  dydzlumi->getTH1()->Fit(myLinFit,"QR");
1035 
1036  goodVxCounter->ShiftFillLast((double)counterVx, std::sqrt((double)counterVx), nLumiReset);
1037  myLinFit->SetParameter(0, goodVxCounter->getTH1()->GetMean(2));
1038  myLinFit->SetParameter(1, 0.0);
1039  goodVxCounter->getTH1()->Fit(myLinFit,"QR");
1040 
1041  if (lastLumiOfFit % prescaleHistory == 0)
1042  {
1043  goodVxCountHistory->getTH1()->SetBinContent(lastLumiOfFit, (double)counterVx);
1044  goodVxCountHistory->getTH1()->SetBinError(lastLumiOfFit, std::sqrt((double)counterVx));
1045  }
1046 
1047  delete myLinFit;
1048 
1049  vals.clear();
1050  }
1051  else if (nLumiReset == 0)
1052  {
1053  histTitle << "Fitted Beam Spot [cm] (no ongoing fits)";
1054  fitResults->setAxisTitle(histTitle.str().c_str(), 1);
1055  reportSummaryMap->Fill(0.5, 0.5, 1.0);
1056  hitCounter->ShiftFillLast(totalHits, std::sqrt(totalHits), 1);
1057  reset("nohisto");
1058  }
1059 }
1060 
1061 
1063 {
1064  // ### Set internal variables ###
1065  reset("scratch");
1066  prescaleHistory = 1; // Set the number of lumis to update historical plot
1067  maxLumiIntegration = 15; // If failing fits, this is the maximum number of integrated lumis after which a reset is issued
1068  minVxDoF = 10.; // Good vertex selection cut
1069  // For vertex fitter without track-weight: d.o.f. = 2*NTracks - 3
1070  // For vertex fitter with track-weight: d.o.f. = sum_NTracks(2*track_weight) - 3
1071  internalDebug = false;
1072  considerVxCovariance = true; // Deconvolute vertex covariance matrix
1073  pi = 3.141592653589793238;
1074  // ##############################
1075 }
1076 
1077 
1078 void Vx3DHLTAnalyzer::endJob() { reset("scratch"); }
1079 
1080 
1082 {
1083  DQMStore* dbe = 0;
1084  dbe = Service<DQMStore>().operator->();
1085 
1086  // ### Set internal variables ###
1087  nBinsHistoricalPlot = 80;
1088  nBinsWholeHistory = 3000; // Correspond to about 20h of data taking: 20h * 60min * 60s / 23s per lumi-block = 3130
1089  // ##############################
1090 
1091  if ( dbe )
1092  {
1093  dbe->setCurrentFolder("BeamPixel");
1094 
1095  Vx_X = dbe->book1D("vertex x", "Primary Vertex X Coordinate Distribution", int(rint(xRange/xStep)), -xRange/2., xRange/2.);
1096  Vx_Y = dbe->book1D("vertex y", "Primary Vertex Y Coordinate Distribution", int(rint(yRange/yStep)), -yRange/2., yRange/2.);
1097  Vx_Z = dbe->book1D("vertex z", "Primary Vertex Z Coordinate Distribution", int(rint(zRange/zStep)), -zRange/2., zRange/2.);
1098  Vx_X->setAxisTitle("Primary Vertices X [cm]",1);
1099  Vx_X->setAxisTitle("Entries [#]",2);
1100  Vx_Y->setAxisTitle("Primary Vertices Y [cm]",1);
1101  Vx_Y->setAxisTitle("Entries [#]",2);
1102  Vx_Z->setAxisTitle("Primary Vertices Z [cm]",1);
1103  Vx_Z->setAxisTitle("Entries [#]",2);
1104 
1105  mXlumi = dbe->book1D("muX vs lumi", "\\mu_{x} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
1106  mYlumi = dbe->book1D("muY vs lumi", "\\mu_{y} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
1107  mZlumi = dbe->book1D("muZ vs lumi", "\\mu_{z} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
1108  mXlumi->setAxisTitle("Lumisection [#]",1);
1109  mXlumi->setAxisTitle("\\mu_{x} [cm]",2);
1110  mXlumi->getTH1()->SetOption("E1");
1111  mYlumi->setAxisTitle("Lumisection [#]",1);
1112  mYlumi->setAxisTitle("\\mu_{y} [cm]",2);
1113  mYlumi->getTH1()->SetOption("E1");
1114  mZlumi->setAxisTitle("Lumisection [#]",1);
1115  mZlumi->setAxisTitle("\\mu_{z} [cm]",2);
1116  mZlumi->getTH1()->SetOption("E1");
1117 
1118  sXlumi = dbe->book1D("sigmaX vs lumi", "\\sigma_{x} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
1119  sYlumi = dbe->book1D("sigmaY vs lumi", "\\sigma_{y} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
1120  sZlumi = dbe->book1D("sigmaZ vs lumi", "\\sigma_{z} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
1121  sXlumi->setAxisTitle("Lumisection [#]",1);
1122  sXlumi->setAxisTitle("\\sigma_{x} [cm]",2);
1123  sXlumi->getTH1()->SetOption("E1");
1124  sYlumi->setAxisTitle("Lumisection [#]",1);
1125  sYlumi->setAxisTitle("\\sigma_{y} [cm]",2);
1126  sYlumi->getTH1()->SetOption("E1");
1127  sZlumi->setAxisTitle("Lumisection [#]",1);
1128  sZlumi->setAxisTitle("\\sigma_{z} [cm]",2);
1129  sZlumi->getTH1()->SetOption("E1");
1130 
1131  dxdzlumi = dbe->book1D("dxdz vs lumi", "dX/dZ vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
1132  dydzlumi = dbe->book1D("dydz vs lumi", "dY/dZ vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
1133  dxdzlumi->setAxisTitle("Lumisection [#]",1);
1134  dxdzlumi->setAxisTitle("dX/dZ [rad]",2);
1135  dxdzlumi->getTH1()->SetOption("E1");
1136  dydzlumi->setAxisTitle("Lumisection [#]",1);
1137  dydzlumi->setAxisTitle("dY/dZ [rad]",2);
1138  dydzlumi->getTH1()->SetOption("E1");
1139 
1140  Vx_ZX = dbe->book2D("vertex zx", "Primary Vertex ZX Coordinate Distribution", int(rint(zRange/zStep/5.)), -zRange/2., zRange/2., int(rint(xRange/xStep/5.)), -xRange/2., xRange/2.);
1141  Vx_ZY = dbe->book2D("vertex zy", "Primary Vertex ZY Coordinate Distribution", int(rint(zRange/zStep/5.)), -zRange/2., zRange/2., int(rint(yRange/yStep/5.)), -yRange/2., yRange/2.);
1142  Vx_XY = dbe->book2D("vertex xy", "Primary Vertex XY Coordinate Distribution", int(rint(xRange/xStep/5.)), -xRange/2., xRange/2., int(rint(yRange/yStep/5.)), -yRange/2., yRange/2.);
1143  Vx_ZX->setAxisTitle("Primary Vertices Z [cm]",1);
1144  Vx_ZX->setAxisTitle("Primary Vertices X [cm]",2);
1145  Vx_ZX->setAxisTitle("Entries [#]",3);
1146  Vx_ZY->setAxisTitle("Primary Vertices Z [cm]",1);
1147  Vx_ZY->setAxisTitle("Primary Vertices Y [cm]",2);
1148  Vx_ZY->setAxisTitle("Entries [#]",3);
1149  Vx_XY->setAxisTitle("Primary Vertices X [cm]",1);
1150  Vx_XY->setAxisTitle("Primary Vertices Y [cm]",2);
1151  Vx_XY->setAxisTitle("Entries [#]",3);
1152 
1153  hitCounter = dbe->book1D("pixelHits vs lumi", "# Pixel-Hits vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
1154  hitCounter->setAxisTitle("Lumisection [#]",1);
1155  hitCounter->setAxisTitle("Pixel-Hits [#]",2);
1156  hitCounter->getTH1()->SetOption("E1");
1157 
1158  hitCountHistory = dbe->book1D("hist pixelHits vs lumi", "History: # Pixel-Hits vs. Lumi", nBinsWholeHistory, 0.5, (double)nBinsWholeHistory+0.5);
1159  hitCountHistory->setAxisTitle("Lumisection [#]",1);
1160  hitCountHistory->setAxisTitle("Pixel-Hits [#]",2);
1161  hitCountHistory->getTH1()->SetOption("E1");
1162 
1163  goodVxCounter = dbe->book1D("good vertices vs lumi", "# Good vertices vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
1164  goodVxCounter->setAxisTitle("Lumisection [#]",1);
1165  goodVxCounter->setAxisTitle("Good vertices [#]",2);
1166  goodVxCounter->getTH1()->SetOption("E1");
1167 
1168  goodVxCountHistory = dbe->book1D("hist good vx vs lumi", "History: # Good vx vs. Lumi", nBinsWholeHistory, 0.5, (double)nBinsWholeHistory+0.5);
1169  goodVxCountHistory->setAxisTitle("Lumisection [#]",1);
1170  goodVxCountHistory->setAxisTitle("Good vertices [#]",2);
1171  goodVxCountHistory->getTH1()->SetOption("E1");
1172 
1173  fitResults = dbe->book2D("fit results","Results of Beam Spot Fit", 2, 0., 2., 9, 0., 9.);
1174  fitResults->setAxisTitle("Fitted Beam Spot [cm]", 1);
1175  fitResults->setBinLabel(9, "X", 2);
1176  fitResults->setBinLabel(8, "Y", 2);
1177  fitResults->setBinLabel(7, "Z", 2);
1178  fitResults->setBinLabel(6, "\\sigma_{Z}", 2);
1179  fitResults->setBinLabel(5, "#frac{dX}{dZ}[rad]", 2);
1180  fitResults->setBinLabel(4, "#frac{dY}{dZ}[rad]", 2);
1181  fitResults->setBinLabel(3, "\\sigma_{X}", 2);
1182  fitResults->setBinLabel(2, "\\sigma_{Y}", 2);
1183  fitResults->setBinLabel(1, "Vertices", 2);
1184  fitResults->setBinLabel(1, "Value", 1);
1185  fitResults->setBinLabel(2, "Stat. Error", 1);
1186  fitResults->getTH1()->SetOption("text");
1187 
1188  dbe->setCurrentFolder("BeamPixel/EventInfo");
1189  reportSummary = dbe->bookFloat("reportSummary");
1190  reportSummary->Fill(0.);
1191  reportSummaryMap = dbe->book2D("reportSummaryMap","Pixel-Vertices Beam Spot: % Good Fits", 1, 0., 1., 1, 0., 1.);
1192  reportSummaryMap->Fill(0.5, 0.5, 0.);
1193  dbe->setCurrentFolder("BeamPixel/EventInfo/reportSummaryContents");
1194 
1195  // Convention for reportSummary and reportSummaryMap:
1196  // - 0% at the moment of creation of the histogram
1197  // - n% numberGoodFits / numberFits
1198  }
1199 }
1200 
1201 
1202 // Define this as a plug-in
RunNumber_t run() const
Definition: EventID.h:39
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
boost::transform_iterator< IterHelp, const_IdIter > const_iterator
string reportSummary
double maxLongLength
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:446
double zPos
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
virtual int MyFit(std::vector< double > *vals)
#define NULL
Definition: scimark2.h:8
data_type const * const_iterator
Definition: DetSetNew.h:30
Timestamp const & beginTime() const
virtual unsigned int HitCounter(const edm::Event &iEvent)
tuple vertexCollection
float float float z
virtual void writeToFile(std::vector< double > *vals, edm::TimeValue_t BeginTimeOfFit, edm::TimeValue_t EndTimeOfFit, unsigned int BeginLumiOfFit, unsigned int EndLumiOfFit, int dataType)
const Double_t pi
double maxTransRadius
LuminosityBlockNumber_t luminosityBlock() const
double Gauss3DFunc(const double *par)
bool considerVxCovariance
double xPos
virtual void endLuminosityBlock(const edm::LuminosityBlock &lumiBlock, const edm::EventSetup &iSetup)
static ELstring formatTime(const time_t t)
Definition: ELoutput.cc:100
Vx3DHLTAnalyzer(const edm::ParameterSet &)
int iEvent
Definition: GenABIO.cc:230
bool isNotFinite(T x)
Definition: isFinite.h:10
virtual void beginJob()
virtual void beginRun()
T sqrt(T t)
Definition: SSEVec.h:48
Timestamp const & endTime() const
virtual void analyze(const edm::Event &, const edm::EventSetup &)
int j
Definition: DBlmapReader.cc:9
double yPos
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
LuminosityBlock const & getLuminosityBlock() const
Definition: Event.h:77
virtual void reset(std::string ResetType)
unsigned long long TimeValue_t
Definition: Timestamp.h:28
double Covariance[DIM][DIM]
virtual void beginLuminosityBlock(const edm::LuminosityBlock &lumiBlock, const edm::EventSetup &iSetup)
tuple out
Definition: dbtoconf.py:99
#define DIM
virtual std::string formatTime(const time_t &t)
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
edm::EventID id() const
Definition: EventBase.h:56
static std::atomic< unsigned int > counter
std::vector< VertexType > Vertices
virtual void endJob()
unsigned int counterVx
tuple cout
Definition: gather_cfg.py:121
double VxErrCorr
Definition: DDAxes.h:10
void reset(double vett[256])
Definition: TPedValues.cc:11
TimeValue_t value() const
Definition: Timestamp.h:56
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40