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