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