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