16 #include <Math/Minimizer.h>
17 #include <Math/Factory.h>
18 #include <Math/Functor.h>
70 internalDebug =
false;
71 considerVxCovariance =
true;
72 pi = 3.141592653589793238;
91 stringstream debugFile;
94 if (outputDebugFile.is_open() ==
true)
95 outputDebugFile.close();
96 tmp.erase(strlen(
fileName.c_str()) - 4, 4);
97 debugFile << tmp.c_str() <<
"_Run" << iEvent.
id().
run() <<
".txt";
98 outputDebugFile.open(debugFile.str().c_str(),
ios::out);
99 outputDebugFile.close();
100 outputDebugFile.open(debugFile.str().c_str(), ios::app);
104 }
else if (beginTimeOfFit != 0) {
105 totalHits += HitCounter(iEvent);
107 if (internalDebug ==
true) {
108 cout <<
"[Vx3DHLTAnalyzer]::\tI found " << totalHits <<
" pixel hits until now" << endl;
109 cout <<
"[Vx3DHLTAnalyzer]::\tIn this event there are " << Vx3DCollection->size() <<
" vertex cadidates" << endl;
112 for (vector<Vertex>::const_iterator it3DVx = Vx3DCollection->begin(); it3DVx != Vx3DCollection->end(); it3DVx++) {
113 if (internalDebug ==
true) {
114 cout <<
"[Vx3DHLTAnalyzer]::\tVertex selections:" << endl;
115 cout <<
"[Vx3DHLTAnalyzer]::\tEvent ID = " << iEvent.
id() << endl;
116 cout <<
"[Vx3DHLTAnalyzer]::\tVertex number = " << it3DVx - Vx3DCollection->begin() << endl;
117 cout <<
"[Vx3DHLTAnalyzer]::\tisValid = " << it3DVx->isValid() << endl;
118 cout <<
"[Vx3DHLTAnalyzer]::\tisFake = " << it3DVx->isFake() << endl;
119 cout <<
"[Vx3DHLTAnalyzer]::\tnodof = " << it3DVx->ndof() << endl;
120 cout <<
"[Vx3DHLTAnalyzer]::\ttracksSize = " << it3DVx->tracksSize() << endl;
123 if ((it3DVx->isValid() ==
true) && (it3DVx->isFake() ==
false) && (it3DVx->ndof() >=
minVxDoF) &&
124 (it3DVx->tracksSize() > 0) && ((it3DVx->ndof() + 3.) / ((double)it3DVx->tracksSize()) >= 2. *
minVxWgt)) {
125 for (i = 0; i <
DIM; i++) {
126 for (j = 0; j <
DIM; j++) {
145 if ((i == DIM) && (det > 0.)) {
146 if (internalDebug ==
true)
147 cout <<
"[Vx3DHLTAnalyzer]::\tVertex accepted !" << endl;
149 MyVertex.
x = it3DVx->x();
150 MyVertex.
y = it3DVx->y();
151 MyVertex.
z = it3DVx->z();
152 Vertices.push_back(MyVertex);
154 Vx_X->Fill(it3DVx->x());
155 Vx_Y->Fill(it3DVx->y());
156 Vx_Z->Fill(it3DVx->z());
158 Vx_ZX->Fill(it3DVx->z(), it3DVx->x());
159 Vx_ZY->Fill(it3DVx->z(), it3DVx->y());
160 Vx_XY->Fill(it3DVx->x(), it3DVx->y());
162 Vx_X_Cum->Fill(it3DVx->x());
163 Vx_Y_Cum->Fill(it3DVx->y());
164 Vx_Z_Cum->Fill(it3DVx->z());
166 Vx_ZX_Cum->Fill(it3DVx->z(), it3DVx->x());
167 Vx_ZY_Cum->Fill(it3DVx->z(), it3DVx->y());
168 Vx_XY_Cum->Fill(it3DVx->x(), it3DVx->y());
169 }
else if (internalDebug ==
true) {
170 cout <<
"[Vx3DHLTAnalyzer]::\tVertex discarded !" << endl;
172 for (i = 0; i <
DIM; i++)
173 for (j = 0; j <
DIM; j++)
174 cout <<
"(i,j) --> " << i <<
"," << j <<
" --> " << MyVertex.
Covariance[i][j] << endl;
176 }
else if (internalDebug ==
true)
177 cout <<
"[Vx3DHLTAnalyzer]::\tVertex discarded !" << endl;
190 counter +=
h->cluster()->size();
197 strftime(ts,
sizeof(ts),
"%Y.%m.%d %H:%M:%S %Z", gmtime(&t));
199 string ts_string(ts);
221 for (
unsigned int i = 0;
i < Vertices.size();
i++) {
222 if ((
std::sqrt((Vertices[
i].x - xPos) * (Vertices[
i].x - xPos) + (Vertices[
i].y - yPos) * (Vertices[
i].y - yPos)) <=
224 (std::fabs(Vertices[
i].z - zPos) <= maxLongLength)) {
225 if (considerVxCovariance ==
true) {
227 K[1][1] = std::fabs(par[1]) +
VxErrCorr *
VxErrCorr * std::fabs(Vertices[
i].Covariance[1][1]);
228 K[2][2] = std::fabs(par[2]) +
VxErrCorr *
VxErrCorr * std::fabs(Vertices[
i].Covariance[2][2]);
230 K[1][2] = K[2][1] = par[4] * (std::fabs(par[2]) - std::fabs(par[1])) - par[5] * par[3] +
232 K[0][2] = K[2][0] = par[5] * (std::fabs(par[2]) - std::fabs(par[0])) - par[4] * par[3] +
235 K[0][0] = std::fabs(par[0]);
236 K[1][1] = std::fabs(par[1]);
237 K[2][2] = std::fabs(par[2]);
238 K[0][1] = K[1][0] = par[3];
239 K[1][2] = K[2][1] = par[4] * (std::fabs(par[2]) - std::fabs(par[1])) - par[5] * par[3];
240 K[0][2] = K[2][0] = par[5] * (std::fabs(par[2]) - std::fabs(par[0])) - par[4] * par[3];
243 det = K[0][0] * (K[1][1] * K[2][2] - K[1][2] * K[1][2]) - K[0][1] * (K[0][1] * K[2][2] - K[0][2] * K[1][2]) +
244 K[0][2] * (K[0][1] * K[1][2] - K[0][2] * K[1][1]);
246 M[0][0] = (K[1][1] * K[2][2] - K[1][2] * K[1][2]) / det;
247 M[1][1] = (K[0][0] * K[2][2] - K[0][2] * K[0][2]) / det;
248 M[2][2] = (K[0][0] * K[1][1] - K[0][1] * K[0][1]) / det;
249 M[0][1] = M[1][0] = (K[0][2] * K[1][2] - K[0][1] * K[2][2]) / det;
250 M[1][2] = M[2][1] = (K[0][2] * K[0][1] - K[1][2] * K[0][0]) / det;
251 M[0][2] = M[2][0] = (K[0][1] * K[1][2] - K[0][2] * K[1][1]) / det;
254 (M[0][0] * (Vertices[
i].x - par[6]) * (Vertices[
i].x - par[6]) +
255 M[1][1] * (Vertices[
i].y - par[7]) * (Vertices[
i].y - par[7]) +
256 M[2][2] * (Vertices[
i].z - par[8]) * (Vertices[
i].z - par[8]) +
257 2. * M[0][1] * (Vertices[
i].x - par[6]) * (Vertices[
i].y - par[7]) +
258 2. * M[1][2] * (Vertices[
i].y - par[7]) * (Vertices[
i].z - par[8]) +
259 2. * M[0][2] * (Vertices[
i].x - par[6]) * (Vertices[
i].z - par[8]));
280 if ((vals !=
nullptr) && (vals->size() == nParams * 2)) {
281 double nSigmaXY = 10.;
283 double parDistanceXY = 1
e-3;
284 double parDistanceZ = 1
e-2;
285 double parDistanceddZ = 1
e-3;
286 double parDistanceCxy = 1
e-5;
289 const unsigned int trials = 4;
290 double largerDist[trials] = {0.1, 5., 10., 100.};
292 double covxz, covyz, det;
294 int bestMovementX = 1;
295 int bestMovementY = 1;
296 int bestMovementZ = 1;
301 vector<double>::const_iterator it = vals->begin();
303 ROOT::Math::Minimizer* Gauss3D = ROOT::Math::Factory::CreateMinimizer(
"Minuit2",
"Migrad");
304 Gauss3D->SetErrorDef(1.0);
305 if (internalDebug ==
true)
306 Gauss3D->SetPrintLevel(3);
308 Gauss3D->SetPrintLevel(0);
311 Gauss3D->SetFunction(_Gauss3DFunc);
313 if (internalDebug ==
true)
314 cout <<
"[Vx3DHLTAnalyzer]::\t@@@ START FITTING @@@" << endl;
318 for (
int i = 0;
i < 3;
i++) {
319 deltaMean = (double(
i) - 1.) *
std::sqrt(*(it + 0));
320 if (internalDebug ==
true)
321 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl;
325 Gauss3D->SetVariable(0,
"var x ", *(it + 0), parDistanceXY * parDistanceXY);
326 Gauss3D->SetVariable(1,
"var y ", *(it + 1), parDistanceXY * parDistanceXY);
327 Gauss3D->SetVariable(2,
"var z ", *(it + 2), parDistanceZ * parDistanceZ);
328 Gauss3D->SetVariable(3,
"cov xy", *(it + 3), parDistanceCxy);
329 Gauss3D->SetVariable(4,
"dydz ", *(it + 4), parDistanceddZ);
330 Gauss3D->SetVariable(5,
"dxdz ", *(it + 5), parDistanceddZ);
331 Gauss3D->SetVariable(6,
"mean x", *(it + 6) + deltaMean, parDistanceXY);
332 Gauss3D->SetVariable(7,
"mean y", *(it + 7), parDistanceXY);
333 Gauss3D->SetVariable(8,
"mean z", *(it + 8), parDistanceZ);
341 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs((*vals)[0]) + std::fabs((*vals)[1])) / 2.;
342 maxLongLength = nSigmaZ *
std::sqrt(std::fabs((*vals)[2]));
345 goodData = Gauss3D->Status();
346 edm = Gauss3D->Edm();
352 if (internalDebug ==
true)
353 cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl;
355 for (
unsigned int j = 0;
j < nParams;
j++)
358 if (internalDebug ==
true)
359 cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
363 covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
364 Gauss3D->X()[5] * Gauss3D->X()[3];
365 covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
366 Gauss3D->X()[4] * Gauss3D->X()[3];
368 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
369 Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
370 covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
373 if (internalDebug ==
true)
374 cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl;
378 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) {
383 if (internalDebug ==
true)
384 cout <<
"[Vx3DHLTAnalyzer]::\tFound bestMovementX --> " << bestMovementX << endl;
388 for (
int i = 0;
i < 3;
i++) {
389 deltaMean = (double(
i) - 1.) *
std::sqrt(*(it + 1));
390 if (internalDebug ==
true) {
391 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl;
392 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean X --> " << (double(bestMovementX) - 1.) *
std::sqrt(*(it + 0)) << endl;
397 Gauss3D->SetVariable(0,
"var x ", *(it + 0), parDistanceXY * parDistanceXY);
398 Gauss3D->SetVariable(1,
"var y ", *(it + 1), parDistanceXY * parDistanceXY);
399 Gauss3D->SetVariable(2,
"var z ", *(it + 2), parDistanceZ * parDistanceZ);
400 Gauss3D->SetVariable(3,
"cov xy", *(it + 3), parDistanceCxy);
401 Gauss3D->SetVariable(4,
"dydz ", *(it + 4), parDistanceddZ);
402 Gauss3D->SetVariable(5,
"dxdz ", *(it + 5), parDistanceddZ);
403 Gauss3D->SetVariable(6,
"mean x", *(it + 6) + (
double(bestMovementX) - 1.) *
std::sqrt(*(it + 0)), parDistanceXY);
404 Gauss3D->SetVariable(7,
"mean y", *(it + 7) + deltaMean, parDistanceXY);
405 Gauss3D->SetVariable(8,
"mean z", *(it + 8), parDistanceZ);
408 xPos = Gauss3D->X()[6];
409 yPos = Gauss3D->X()[7];
410 zPos = Gauss3D->X()[8];
413 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(Gauss3D->X()[0]) + std::fabs(Gauss3D->X()[1])) / 2.;
414 maxLongLength = nSigmaZ *
std::sqrt(std::fabs(Gauss3D->X()[2]));
417 goodData = Gauss3D->Status();
418 edm = Gauss3D->Edm();
424 if (internalDebug ==
true)
425 cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl;
427 for (
unsigned int j = 0;
j < nParams;
j++)
430 if (internalDebug ==
true)
431 cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
435 covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
436 Gauss3D->X()[5] * Gauss3D->X()[3];
437 covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
438 Gauss3D->X()[4] * Gauss3D->X()[3];
440 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
441 Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
442 covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
445 if (internalDebug ==
true)
446 cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl;
450 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) {
455 if (internalDebug ==
true)
456 cout <<
"[Vx3DHLTAnalyzer]::\tFound bestMovementY --> " << bestMovementY << endl;
460 for (
int i = 0;
i < 3;
i++) {
461 deltaMean = (double(
i) - 1.) *
std::sqrt(*(it + 2));
462 if (internalDebug ==
true) {
463 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl;
464 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean X --> " << (double(bestMovementX) - 1.) *
std::sqrt(*(it + 0)) << endl;
465 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean Y --> " << (double(bestMovementY) - 1.) *
std::sqrt(*(it + 1)) << endl;
470 Gauss3D->SetVariable(0,
"var x ", *(it + 0), parDistanceXY * parDistanceXY);
471 Gauss3D->SetVariable(1,
"var y ", *(it + 1), parDistanceXY * parDistanceXY);
472 Gauss3D->SetVariable(2,
"var z ", *(it + 2), parDistanceZ * parDistanceZ);
473 Gauss3D->SetVariable(3,
"cov xy", *(it + 3), parDistanceCxy);
474 Gauss3D->SetVariable(4,
"dydz ", *(it + 4), parDistanceddZ);
475 Gauss3D->SetVariable(5,
"dxdz ", *(it + 5), parDistanceddZ);
476 Gauss3D->SetVariable(6,
"mean x", *(it + 6) + (
double(bestMovementX) - 1.) *
std::sqrt(*(it + 0)), parDistanceXY);
477 Gauss3D->SetVariable(7,
"mean y", *(it + 7) + (
double(bestMovementY) - 1.) *
std::sqrt(*(it + 1)), parDistanceXY);
478 Gauss3D->SetVariable(8,
"mean z", *(it + 8) + deltaMean, parDistanceZ);
481 xPos = Gauss3D->X()[6];
482 yPos = Gauss3D->X()[7];
483 zPos = Gauss3D->X()[8];
486 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(Gauss3D->X()[0]) + std::fabs(Gauss3D->X()[1])) / 2.;
487 maxLongLength = nSigmaZ *
std::sqrt(std::fabs(Gauss3D->X()[2]));
490 goodData = Gauss3D->Status();
491 edm = Gauss3D->Edm();
497 if (internalDebug ==
true)
498 cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl;
500 for (
unsigned int j = 0;
j < nParams;
j++)
503 if (internalDebug ==
true)
504 cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
508 covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
509 Gauss3D->X()[5] * Gauss3D->X()[3];
510 covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
511 Gauss3D->X()[4] * Gauss3D->X()[3];
513 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
514 Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
515 covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
518 if (internalDebug ==
true)
519 cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl;
523 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) {
528 if (internalDebug ==
true)
529 cout <<
"[Vx3DHLTAnalyzer]::\tFound bestMovementZ --> " << bestMovementZ << endl;
534 Gauss3D->SetVariable(0,
"var x ", *(it + 0), parDistanceXY * parDistanceXY);
535 Gauss3D->SetVariable(1,
"var y ", *(it + 1), parDistanceXY * parDistanceXY);
536 Gauss3D->SetVariable(2,
"var z ", *(it + 2), parDistanceZ * parDistanceZ);
537 Gauss3D->SetVariable(3,
"cov xy", *(it + 3), parDistanceCxy);
538 Gauss3D->SetVariable(4,
"dydz ", *(it + 4), parDistanceddZ);
539 Gauss3D->SetVariable(5,
"dxdz ", *(it + 5), parDistanceddZ);
540 Gauss3D->SetVariable(6,
"mean x", *(it + 6) + (
double(bestMovementX) - 1.) *
std::sqrt(*(it + 0)), parDistanceXY);
541 Gauss3D->SetVariable(7,
"mean y", *(it + 7) + (
double(bestMovementY) - 1.) *
std::sqrt(*(it + 1)), parDistanceXY);
542 Gauss3D->SetVariable(8,
"mean z", *(it + 8) + (
double(bestMovementZ) - 1.) *
std::sqrt(*(it + 2)), parDistanceZ);
550 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs((*vals)[0]) + std::fabs((*vals)[1])) / 2.;
551 maxLongLength = nSigmaZ *
std::sqrt(std::fabs((*vals)[2]));
554 goodData = Gauss3D->Status();
555 edm = Gauss3D->Edm();
561 if (internalDebug ==
true)
562 cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl;
564 for (
unsigned int j = 0;
j < nParams;
j++)
567 if (internalDebug ==
true)
568 cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
572 covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
573 Gauss3D->X()[5] * Gauss3D->X()[3];
574 covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
575 Gauss3D->X()[4] * Gauss3D->X()[3];
577 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
578 Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
579 covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
582 if (internalDebug ==
true)
583 cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl;
588 for (
unsigned int i = 0;
i < trials;
i++) {
589 if ((goodData != 0) && (goodData != -2)) {
592 if (internalDebug ==
true)
593 cout <<
"[Vx3DHLTAnalyzer]::\tFIT WITH DIFFERENT PARAMETER DISTANCES - STEP " <<
i + 1 << endl;
595 Gauss3D->SetVariable(0,
"var x ", *(it + 0), parDistanceXY * parDistanceXY * largerDist[
i]);
596 Gauss3D->SetVariable(1,
"var y ", *(it + 1), parDistanceXY * parDistanceXY * largerDist[i]);
597 Gauss3D->SetVariable(2,
"var z ", *(it + 2), parDistanceZ * parDistanceZ * largerDist[i]);
598 Gauss3D->SetVariable(3,
"cov xy", *(it + 3), parDistanceCxy * largerDist[i]);
599 Gauss3D->SetVariable(4,
"dydz ", *(it + 4), parDistanceddZ * largerDist[i]);
600 Gauss3D->SetVariable(5,
"dxdz ", *(it + 5), parDistanceddZ * largerDist[i]);
601 Gauss3D->SetVariable(6,
603 *(it + 6) + (
double(bestMovementX) - 1.) * std::sqrt(*(it + 0)),
604 parDistanceXY * largerDist[i]);
605 Gauss3D->SetVariable(7,
607 *(it + 7) + (
double(bestMovementY) - 1.) * std::sqrt(*(it + 1)),
608 parDistanceXY * largerDist[i]);
609 Gauss3D->SetVariable(
610 8,
"mean z", *(it + 8) + (
double(bestMovementZ) - 1.) * std::sqrt(*(it + 2)), parDistanceZ * largerDist[i]);
613 xPos = Gauss3D->X()[6];
614 yPos = Gauss3D->X()[7];
615 zPos = Gauss3D->X()[8];
618 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(Gauss3D->X()[0]) + std::fabs(Gauss3D->X()[1])) / 2.;
619 maxLongLength = nSigmaZ *
std::sqrt(std::fabs(Gauss3D->X()[2]));
622 goodData = Gauss3D->Status();
623 edm = Gauss3D->Edm();
629 if (internalDebug ==
true)
630 cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl;
632 for (
unsigned int j = 0;
j < nParams;
j++)
635 if (internalDebug ==
true)
636 cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
640 covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
641 Gauss3D->X()[5] * Gauss3D->X()[3];
642 covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
643 Gauss3D->X()[4] * Gauss3D->X()[3];
645 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
646 Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
647 covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
650 if (internalDebug ==
true)
651 cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl;
659 for (
unsigned int i = 0;
i < nParams;
i++) {
660 vals->operator[](
i) = Gauss3D->X()[
i];
661 vals->operator[](
i + nParams) = Gauss3D->Errors()[
i];
672 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true)) {
673 outputDebugFile <<
"Runnumber " <<
runNumber << endl;
674 outputDebugFile <<
"BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) <<
" " << (beginTimeOfFit >> 32) << endl;
675 outputDebugFile <<
"BeginLumiRange " << beginLumiOfFit << endl;
676 outputDebugFile <<
"EndTimeOfFit " << formatTime(endTimeOfFit >> 32) <<
" " << (endTimeOfFit >> 32) << endl;
677 outputDebugFile <<
"EndLumiRange " << endLumiOfFit << endl;
678 outputDebugFile <<
"LumiCounter " << lumiCounter << endl;
679 outputDebugFile <<
"LastLumiOfFit " << lastLumiOfFit << endl;
682 if (ResetType ==
"scratch") {
720 goodVxCounter->Reset();
721 statusCounter->Reset();
724 reportSummary->Fill(-1);
725 reportSummaryMap->getTH1()->SetBinContent(1, 1, -1);
736 if (internalDebug ==
true)
737 cout <<
"[Vx3DHLTAnalyzer]::\tReset issued: scratch" << endl;
738 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true))
739 outputDebugFile <<
"Reset -scratch- issued\n" << endl;
740 }
else if (ResetType ==
"whole") {
758 if (internalDebug ==
true)
759 cout <<
"[Vx3DHLTAnalyzer]::\tReset issued: whole" << endl;
760 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true))
761 outputDebugFile <<
"Reset -whole- issued\n" << endl;
762 }
else if (ResetType ==
"fit") {
767 if (internalDebug ==
true)
768 cout <<
"[Vx3DHLTAnalyzer]::\tReset issued: fit" << endl;
769 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true))
770 outputDebugFile <<
"Reset -fit- issued\n" << endl;
771 }
else if (ResetType ==
"hitCounter") {
774 if (internalDebug ==
true)
775 cout <<
"[Vx3DHLTAnalyzer]::\tReset issued: hitCounter" << endl;
776 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true))
777 outputDebugFile <<
"Reset -hitCounter- issued\n" << endl;
784 unsigned int BeginLumiOfFit,
785 unsigned int EndLumiOfFit,
787 stringstream BufferString;
788 BufferString.precision(5);
792 if ((
outputFile.is_open() ==
true) && (vals !=
nullptr) && (vals->size() == (nParams - 1) * 2)) {
793 vector<double>::const_iterator it = vals->begin();
796 outputFile <<
"BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) <<
" " << (beginTimeOfFit >> 32) << endl;
797 outputFile <<
"EndTimeOfFit " << formatTime(endTimeOfFit >> 32) <<
" " << (endTimeOfFit >> 32) << endl;
798 outputFile <<
"LumiRange " << beginLumiOfFit <<
" - " << endLumiOfFit << endl;
804 BufferString << *(it + 0);
805 outputFile <<
"X0 " << BufferString.str().c_str() << endl;
806 BufferString.str(
"");
808 BufferString << *(it + 1);
809 outputFile <<
"Y0 " << BufferString.str().c_str() << endl;
810 BufferString.str(
"");
812 BufferString << *(it + 2);
813 outputFile <<
"Z0 " << BufferString.str().c_str() << endl;
814 BufferString.str(
"");
816 BufferString << *(it + 5);
817 outputFile <<
"sigmaZ0 " << BufferString.str().c_str() << endl;
818 BufferString.str(
"");
820 BufferString << *(it + 6);
821 outputFile <<
"dxdz " << BufferString.str().c_str() << endl;
822 BufferString.str(
"");
824 BufferString << *(it + 7);
825 outputFile <<
"dydz " << BufferString.str().c_str() << endl;
826 BufferString.str(
"");
828 BufferString << *(it + 3);
829 outputFile <<
"BeamWidthX " << BufferString.str().c_str() << endl;
830 BufferString.str(
"");
832 BufferString << *(it + 4);
833 outputFile <<
"BeamWidthY " << BufferString.str().c_str() << endl;
834 BufferString.str(
"");
836 outputFile <<
"Cov(0,j) " << *(it + 8) <<
" 0 0 0 0 0 0" << endl;
837 outputFile <<
"Cov(1,j) 0 " << *(it + 9) <<
" 0 0 0 0 0" << endl;
838 outputFile <<
"Cov(2,j) 0 0 " << *(it + 10) <<
" 0 0 0 0" << endl;
839 outputFile <<
"Cov(3,j) 0 0 0 " << *(it + 13) <<
" 0 0 0" << endl;
840 outputFile <<
"Cov(4,j) 0 0 0 0 " << *(it + 14) <<
" 0 0" << endl;
841 outputFile <<
"Cov(5,j) 0 0 0 0 0 " << *(it + 15) <<
" 0" << endl;
843 << ((*(it + 11)) + (*(it + 12)) + 2. *
std::sqrt((*(it + 11)) * (*(it + 12)))) / 4. << endl;
851 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true) && (vals !=
nullptr) &&
852 (vals->size() == (nParams - 1) * 2)) {
853 vector<double>::const_iterator it = vals->begin();
855 outputDebugFile <<
"Runnumber " <<
runNumber << endl;
856 outputDebugFile <<
"BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) <<
" " << (beginTimeOfFit >> 32) << endl;
857 outputDebugFile <<
"EndTimeOfFit " << formatTime(endTimeOfFit >> 32) <<
" " << (endTimeOfFit >> 32) << endl;
858 outputDebugFile <<
"LumiRange " << beginLumiOfFit <<
" - " << endLumiOfFit << endl;
859 outputDebugFile <<
"Type " << dataType << endl;
864 BufferString << *(it + 0);
865 outputDebugFile <<
"X0 " << BufferString.str().c_str() << endl;
866 BufferString.str(
"");
868 BufferString << *(it + 1);
869 outputDebugFile <<
"Y0 " << BufferString.str().c_str() << endl;
870 BufferString.str(
"");
872 BufferString << *(it + 2);
873 outputDebugFile <<
"Z0 " << BufferString.str().c_str() << endl;
874 BufferString.str(
"");
876 BufferString << *(it + 5);
877 outputDebugFile <<
"sigmaZ0 " << BufferString.str().c_str() << endl;
878 BufferString.str(
"");
880 BufferString << *(it + 6);
881 outputDebugFile <<
"dxdz " << BufferString.str().c_str() << endl;
882 BufferString.str(
"");
884 BufferString << *(it + 7);
885 outputDebugFile <<
"dydz " << BufferString.str().c_str() << endl;
886 BufferString.str(
"");
888 BufferString << *(it + 3);
889 outputDebugFile <<
"BeamWidthX " << BufferString.str().c_str() << endl;
890 BufferString.str(
"");
892 BufferString << *(it + 4);
893 outputDebugFile <<
"BeamWidthY " << BufferString.str().c_str() << endl;
894 BufferString.str(
"");
896 outputDebugFile <<
"Cov(0,j) " << *(it + 8) <<
" 0 0 0 0 0 0" << endl;
897 outputDebugFile <<
"Cov(1,j) 0 " << *(it + 9) <<
" 0 0 0 0 0" << endl;
898 outputDebugFile <<
"Cov(2,j) 0 0 " << *(it + 10) <<
" 0 0 0 0" << endl;
899 outputDebugFile <<
"Cov(3,j) 0 0 0 " << *(it + 13) <<
" 0 0 0" << endl;
900 outputDebugFile <<
"Cov(4,j) 0 0 0 0 " << *(it + 14) <<
" 0 0" << endl;
901 outputDebugFile <<
"Cov(5,j) 0 0 0 0 0 " << *(it + 15) <<
" 0" << endl;
902 outputDebugFile <<
"Cov(6,j) 0 0 0 0 0 0 "
903 << ((*(it + 11)) + (*(it + 12)) + 2. *
std::sqrt((*(it + 11)) * (*(it + 12)))) / 4. << endl;
905 outputDebugFile <<
"EmittanceX 0" << endl;
906 outputDebugFile <<
"EmittanceY 0" << endl;
907 outputDebugFile <<
"BetaStar 0" << endl;
909 outputDebugFile <<
"\nUsed vertices: " << counterVx <<
"\n" << endl;
914 cout <<
"var x --> " << fitResults[0] <<
" +/- " << fitResults[0 + nParams] << endl;
915 cout <<
"var y --> " << fitResults[1] <<
" +/- " << fitResults[1 + nParams] << endl;
916 cout <<
"var z --> " << fitResults[2] <<
" +/- " << fitResults[2 + nParams] << endl;
917 cout <<
"cov xy --> " << fitResults[3] <<
" +/- " << fitResults[3 + nParams] << endl;
918 cout <<
"dydz --> " << fitResults[4] <<
" +/- " << fitResults[4 + nParams] << endl;
919 cout <<
"dxdz --> " << fitResults[5] <<
" +/- " << fitResults[5 + nParams] << endl;
920 cout <<
"mean x --> " << fitResults[6] <<
" +/- " << fitResults[6 + nParams] << endl;
921 cout <<
"mean y --> " << fitResults[7] <<
" +/- " << fitResults[7 + nParams] << endl;
922 cout <<
"mean z --> " << fitResults[8] <<
" +/- " << fitResults[8 + nParams] << endl;
927 if ((lumiCounter == 0) && (lumiBlock.
luminosityBlock() > lastLumiOfFit)) {
931 }
else if ((lumiCounter != 0) && (lumiBlock.
luminosityBlock() >= (beginLumiOfFit + lumiCounter)))
938 stringstream histTitle;
939 double minXfit, maxXfit;
945 lastLumiOfFit = endLumiOfFit;
948 hitCounter->getTH1()->SetBinContent(lastLumiOfFit, (
double)totalHits);
949 hitCounter->getTH1()->SetBinError(
955 vector<double> fitResults;
957 fitResults.push_back(Vx_X->getTH1()->GetRMS() * Vx_X->getTH1()->GetRMS());
958 fitResults.push_back(Vx_Y->getTH1()->GetRMS() * Vx_Y->getTH1()->GetRMS());
959 fitResults.push_back(Vx_Z->getTH1()->GetRMS() * Vx_Z->getTH1()->GetRMS());
960 fitResults.push_back(0.0);
961 fitResults.push_back(0.0);
962 fitResults.push_back(0.0);
963 fitResults.push_back(Vx_X->getTH1()->GetMean());
964 fitResults.push_back(Vx_Y->getTH1()->GetMean());
965 fitResults.push_back(Vx_Z->getTH1()->GetMean());
966 for (
unsigned int i = 0;
i < nParams;
i++)
967 fitResults.push_back(0.0);
969 if (internalDebug ==
true) {
970 cout <<
"[Vx3DHLTAnalyzer]::\t@@@ Beam Spot parameters - prefit @@@" << endl;
972 printFitParams(fitResults);
975 cout <<
"BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) <<
" " << (beginTimeOfFit >> 32) << endl;
976 cout <<
"EndTimeOfFit " << formatTime(endTimeOfFit >> 32) <<
" " << (endTimeOfFit >> 32) << endl;
977 cout <<
"LumiRange " << beginLumiOfFit <<
" - " << endLumiOfFit << endl;
980 goodData = MyFit(&fitResults);
982 if (internalDebug ==
true) {
983 cout <<
"[Vx3DHLTAnalyzer]::\t@@@ Beam Spot parameters - postfit @@@" << endl;
985 printFitParams(fitResults);
987 cout <<
"goodData --> " << goodData << endl;
988 cout <<
"Used vertices --> " << counterVx << endl;
992 vals.push_back(fitResults[6]);
993 vals.push_back(fitResults[7]);
994 vals.push_back(fitResults[8]);
995 vals.push_back(
std::sqrt(std::fabs(fitResults[0])));
996 vals.push_back(
std::sqrt(std::fabs(fitResults[1])));
997 vals.push_back(
std::sqrt(std::fabs(fitResults[2])));
998 vals.push_back(fitResults[5]);
999 vals.push_back(fitResults[4]);
1001 vals.push_back(
std::pow(fitResults[6 + nParams], 2.));
1002 vals.push_back(
std::pow(fitResults[7 + nParams], 2.));
1003 vals.push_back(
std::pow(fitResults[8 + nParams], 2.));
1004 vals.push_back(
std::pow(std::fabs(fitResults[0 + nParams]) / (2. *
std::sqrt(std::fabs(fitResults[0]))), 2.));
1005 vals.push_back(
std::pow(std::fabs(fitResults[1 + nParams]) / (2. *
std::sqrt(std::fabs(fitResults[1]))), 2.));
1006 vals.push_back(
std::pow(std::fabs(fitResults[2 + nParams]) / (2. *
std::sqrt(std::fabs(fitResults[2]))), 2.));
1007 vals.push_back(
std::pow(fitResults[5 + nParams], 2.));
1008 vals.push_back(
std::pow(fitResults[4 + nParams], 2.));
1010 for (
unsigned int i = 0;
i < (nParams - 1) * 2;
i++)
1011 vals.push_back(0.0);
1015 counterVx = Vx_X->getTH1F()->GetEntries();
1017 if (Vx_X->getTH1F()->GetEntries() >=
minNentries) {
1020 vals.push_back(Vx_X->getTH1F()->GetMean());
1021 vals.push_back(Vx_Y->getTH1F()->GetMean());
1022 vals.push_back(Vx_Z->getTH1F()->GetMean());
1023 vals.push_back(Vx_X->getTH1F()->GetRMS());
1024 vals.push_back(Vx_Y->getTH1F()->GetRMS());
1025 vals.push_back(Vx_Z->getTH1F()->GetRMS());
1026 vals.push_back(0.0);
1027 vals.push_back(0.0);
1029 vals.push_back(
std::pow(Vx_X->getTH1F()->GetMeanError(), 2.));
1030 vals.push_back(
std::pow(Vx_Y->getTH1F()->GetMeanError(), 2.));
1031 vals.push_back(
std::pow(Vx_Z->getTH1F()->GetMeanError(), 2.));
1032 vals.push_back(
std::pow(Vx_X->getTH1F()->GetRMSError(), 2.));
1033 vals.push_back(
std::pow(Vx_Y->getTH1F()->GetRMSError(), 2.));
1034 vals.push_back(
std::pow(Vx_Z->getTH1F()->GetRMSError(), 2.));
1035 vals.push_back(0.0);
1036 vals.push_back(0.0);
1039 for (
unsigned int i = 0;
i < (nParams - 1) * 2;
i++)
1040 vals.push_back(0.0);
1063 writeToFile(&vals, beginTimeOfFit, endTimeOfFit, beginLumiOfFit, endLumiOfFit, 3);
1064 if (internalDebug ==
true)
1065 cout <<
"[Vx3DHLTAnalyzer]::\tUsed vertices: " << counterVx << endl;
1067 statusCounter->getTH1()->SetBinContent(lastLumiOfFit, (
double)goodData);
1068 statusCounter->getTH1()->SetBinError(lastLumiOfFit, 1
e-3);
1078 for (
int i = 0;
i < Vx_X_Fit->getTH1()->GetNbinsX();
i++) {
1079 Vx_X_Fit->getTH1()->SetBinContent(
1080 i + 1, Vx_X_Fit->getTH1()->GetBinContent(
i + 1) + Vx_X->getTH1()->GetBinContent(
i + 1));
1081 Vx_X_Fit->getTH1()->SetBinError(
i + 1,
sqrt(Vx_X_Fit->getTH1()->GetBinContent(
i + 1)));
1084 for (
int i = 0;
i < Vx_Y_Fit->getTH1()->GetNbinsX();
i++) {
1085 Vx_Y_Fit->getTH1()->SetBinContent(
1086 i + 1, Vx_Y_Fit->getTH1()->GetBinContent(
i + 1) + Vx_Y->getTH1()->GetBinContent(
i + 1));
1087 Vx_Y_Fit->getTH1()->SetBinError(
i + 1,
sqrt(Vx_Y_Fit->getTH1()->GetBinContent(
i + 1)));
1090 for (
int i = 0;
i < Vx_Z_Fit->getTH1()->GetNbinsX();
i++) {
1091 Vx_Z_Fit->getTH1()->SetBinContent(
1092 i + 1, Vx_Z_Fit->getTH1()->GetBinContent(
i + 1) + Vx_Z->getTH1()->GetBinContent(
i + 1));
1093 Vx_Z_Fit->getTH1()->SetBinError(
i + 1,
sqrt(Vx_Z_Fit->getTH1()->GetBinContent(
i + 1)));
1097 if (goodData == 0) {
1100 histTitle <<
"Ongoing: fitted lumis " << beginLumiOfFit <<
" - " << endLumiOfFit;
1104 histTitle <<
"Ongoing: not enough evts (" << lumiCounter <<
" - " <<
maxLumiIntegration <<
" lumis)";
1106 histTitle <<
"Ongoing: temporary problems (" << lumiCounter <<
" - " << maxLumiIntegration <<
" lumis)";
1108 if (lumiCounter >= maxLumiIntegration) {
1109 statusCounter->getTH1()->SetBinContent(lastLumiOfFit, -5);
1110 statusCounter->getTH1()->SetBinError(lastLumiOfFit, 1
e-3);
1112 reset(
"hitCounter");
1115 reportSummary->Fill((numberFits != 0 ? ((
double)numberGoodFits) / ((
double)numberFits) : -1));
1116 reportSummaryMap->getTH1()->SetBinContent(
1117 1, 1, (numberFits != 0 ? ((
double)numberGoodFits) / ((
double)numberFits) : -1));
1119 fitResults->setAxisTitle(histTitle.str(), 1);
1121 fitResults->setBinContent(1, 9, vals[0]);
1122 fitResults->setBinContent(1, 8, vals[1]);
1123 fitResults->setBinContent(1, 7, vals[2]);
1124 fitResults->setBinContent(1, 6, vals[3]);
1125 fitResults->setBinContent(1, 5, vals[4]);
1126 fitResults->setBinContent(1, 4, vals[5]);
1127 fitResults->setBinContent(1, 3, vals[6]);
1128 fitResults->setBinContent(1, 2, vals[7]);
1129 fitResults->setBinContent(1, 1, counterVx);
1131 fitResults->setBinContent(2, 9,
std::sqrt(vals[8]));
1132 fitResults->setBinContent(2, 8,
std::sqrt(vals[9]));
1133 fitResults->setBinContent(2, 7,
std::sqrt(vals[10]));
1134 fitResults->setBinContent(2, 6,
std::sqrt(vals[11]));
1135 fitResults->setBinContent(2, 5,
std::sqrt(vals[12]));
1136 fitResults->setBinContent(2, 4,
std::sqrt(vals[13]));
1137 fitResults->setBinContent(2, 3,
std::sqrt(vals[14]));
1138 fitResults->setBinContent(2, 2,
std::sqrt(vals[15]));
1139 fitResults->setBinContent(2, 1,
std::sqrt(counterVx));
1142 TF1* myLinFit =
new TF1(
1143 "myLinFit",
"[0] + [1]*x", mXlumi->getTH1()->GetXaxis()->GetXmin(), mXlumi->getTH1()->GetXaxis()->GetXmax());
1144 myLinFit->SetLineColor(2);
1145 myLinFit->SetLineWidth(2);
1146 myLinFit->SetParName(0,
"Inter.");
1147 myLinFit->SetParName(1,
"Slope");
1149 mXlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[0]);
1150 mXlumi->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt(vals[8]));
1151 myLinFit->SetParameter(0, mXlumi->getTH1()->GetMean(2));
1152 myLinFit->SetParameter(1, 0.0);
1153 mXlumi->getTH1()->Fit(myLinFit,
"QR");
1155 mYlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[1]);
1156 mYlumi->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt(vals[9]));
1157 myLinFit->SetParameter(0, mYlumi->getTH1()->GetMean(2));
1158 myLinFit->SetParameter(1, 0.0);
1159 mYlumi->getTH1()->Fit(myLinFit,
"QR");
1161 mZlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[2]);
1162 mZlumi->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt(vals[10]));
1163 myLinFit->SetParameter(0, mZlumi->getTH1()->GetMean(2));
1164 myLinFit->SetParameter(1, 0.0);
1165 mZlumi->getTH1()->Fit(myLinFit,
"QR");
1167 sXlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[3]);
1168 sXlumi->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt(vals[11]));
1169 myLinFit->SetParameter(0, sXlumi->getTH1()->GetMean(2));
1170 myLinFit->SetParameter(1, 0.0);
1171 sXlumi->getTH1()->Fit(myLinFit,
"QR");
1173 sYlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[4]);
1174 sYlumi->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt(vals[12]));
1175 myLinFit->SetParameter(0, sYlumi->getTH1()->GetMean(2));
1176 myLinFit->SetParameter(1, 0.0);
1177 sYlumi->getTH1()->Fit(myLinFit,
"QR");
1179 sZlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[5]);
1180 sZlumi->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt(vals[13]));
1181 myLinFit->SetParameter(0, sZlumi->getTH1()->GetMean(2));
1182 myLinFit->SetParameter(1, 0.0);
1183 sZlumi->getTH1()->Fit(myLinFit,
"QR");
1185 dxdzlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[6]);
1186 dxdzlumi->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt(vals[14]));
1187 myLinFit->SetParameter(0, dxdzlumi->getTH1()->GetMean(2));
1188 myLinFit->SetParameter(1, 0.0);
1189 dxdzlumi->getTH1()->Fit(myLinFit,
"QR");
1191 dydzlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[7]);
1192 dydzlumi->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt(vals[15]));
1193 myLinFit->SetParameter(0, dydzlumi->getTH1()->GetMean(2));
1194 myLinFit->SetParameter(1, 0.0);
1195 dydzlumi->getTH1()->Fit(myLinFit,
"QR");
1197 myLinFit->SetParameter(0, hitCounter->getTH1()->GetMean(2));
1198 myLinFit->SetParameter(1, 0.0);
1199 hitCounter->getTH1()->Fit(myLinFit,
"QR");
1201 goodVxCounter->getTH1()->SetBinContent(lastLumiOfFit, (
double)counterVx);
1202 goodVxCounter->getTH1()->SetBinError(
1204 (counterVx != 0 ? 1.
1206 myLinFit->SetParameter(0, goodVxCounter->getTH1()->GetMean(2));
1207 myLinFit->SetParameter(1, 0.0);
1208 goodVxCounter->getTH1()->Fit(myLinFit,
"QR");
1214 TF1* myGaussFit =
new TF1(
"myGaussFit",
1215 "[0]*exp(-(x-[1])*(x-[1])/(2*[2]*[2]))",
1216 Vx_Z_Fit->getTH1()->GetXaxis()->GetXmin(),
1217 Vx_Z_Fit->getTH1()->GetXaxis()->GetXmax());
1218 myGaussFit->SetLineColor(2);
1219 myGaussFit->SetLineWidth(2);
1220 myGaussFit->SetParName(0,
"Ampl.");
1221 myGaussFit->SetParName(1,
"#mu");
1222 myGaussFit->SetParName(2,
"#sigma");
1224 myGaussFit->SetParameter(0, Vx_X_Fit->getTH1()->GetMaximum());
1225 myGaussFit->SetParameter(1, Vx_X_Fit->getTH1()->GetMean());
1226 myGaussFit->SetParameter(2, Vx_X_Fit->getTH1()->GetRMS());
1227 minXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(1);
1228 for (
int i = 0;
i < Vx_X_Fit->getTH1()->GetNbinsX();
i++) {
1229 if (Vx_X_Fit->getTH1()->GetBinContent(
i + 1) > 0) {
1230 minXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(
i + 1);
1234 maxXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(Vx_X_Fit->getTH1()->GetNbinsX());
1235 for (
int i = Vx_X_Fit->getTH1()->GetNbinsX();
i > 0;
i--) {
1236 if (Vx_X_Fit->getTH1()->GetBinContent(
i) > 0) {
1237 maxXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(
i);
1241 myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1242 if (Vx_X_Fit->getTH1()->GetEntries() > 0)
1243 Vx_X_Fit->getTH1()->Fit(myGaussFit,
"QRL");
1245 myGaussFit->SetParameter(0, Vx_Y_Fit->getTH1()->GetMaximum());
1246 myGaussFit->SetParameter(1, Vx_Y_Fit->getTH1()->GetMean());
1247 myGaussFit->SetParameter(2, Vx_Y_Fit->getTH1()->GetRMS());
1248 minXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(1);
1249 for (
int i = 0;
i < Vx_Y_Fit->getTH1()->GetNbinsX();
i++) {
1250 if (Vx_Y_Fit->getTH1()->GetBinContent(
i + 1) > 0) {
1251 minXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(
i + 1);
1255 maxXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(Vx_Y_Fit->getTH1()->GetNbinsX());
1256 for (
int i = Vx_Y_Fit->getTH1()->GetNbinsX();
i > 0;
i--) {
1257 if (Vx_Y_Fit->getTH1()->GetBinContent(
i) > 0) {
1258 maxXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(
i);
1262 myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1263 if (Vx_Y_Fit->getTH1()->GetEntries() > 0)
1264 Vx_Y_Fit->getTH1()->Fit(myGaussFit,
"QRL");
1266 myGaussFit->SetParameter(0, Vx_Z_Fit->getTH1()->GetMaximum());
1267 myGaussFit->SetParameter(1, Vx_Z_Fit->getTH1()->GetMean());
1268 myGaussFit->SetParameter(2, Vx_Z_Fit->getTH1()->GetRMS());
1269 minXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(1);
1270 for (
int i = 0;
i < Vx_Z_Fit->getTH1()->GetNbinsX();
i++) {
1271 if (Vx_Z_Fit->getTH1()->GetBinContent(
i + 1) > 0) {
1272 minXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(
i + 1);
1276 maxXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(Vx_Z_Fit->getTH1()->GetNbinsX());
1277 for (
int i = Vx_Z_Fit->getTH1()->GetNbinsX();
i > 0;
i--) {
1278 if (Vx_Z_Fit->getTH1()->GetBinContent(
i) > 0) {
1279 maxXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(
i);
1283 myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1284 if (Vx_Z_Fit->getTH1()->GetEntries() > 0)
1285 Vx_Z_Fit->getTH1()->Fit(myGaussFit,
"QRL");
1289 histTitle <<
"Ongoing: accumulating evts (" << lumiCounter %
nLumiFit <<
" - " <<
nLumiFit <<
" in " << lumiCounter
1291 fitResults->setAxisTitle(histTitle.str(), 1);
1292 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true)) {
1293 outputDebugFile <<
"\n"
1295 outputDebugFile <<
"BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) <<
" " << (beginTimeOfFit >> 32) << endl;
1296 outputDebugFile <<
"BeginLumiRange " << beginLumiOfFit << endl;
1297 outputDebugFile << histTitle.str().c_str() <<
"\n" << endl;
1300 histTitle <<
"Ongoing: no ongoing fits";
1301 fitResults->setAxisTitle(histTitle.str(), 1);
1302 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true))
1303 outputDebugFile << histTitle.str().c_str() <<
"\n" << endl;
1307 hitCounter->getTH1()->SetBinContent(endLumiOfFit, (
double)totalHits);
1308 hitCounter->getTH1()->SetBinError(endLumiOfFit,
std::sqrt((
double)totalHits));
1313 if (internalDebug ==
true)
1314 cout <<
"[Vx3DHLTAnalyzer]::\tHistogram title: " << histTitle.str() << endl;
1327 Vx_X->setAxisTitle(
"Entries [#]", 2);
1328 Vx_Y->setAxisTitle(
"Primary Vertices Y [cm]", 1);
1329 Vx_Y->setAxisTitle(
"Entries [#]", 2);
1330 Vx_Z->setAxisTitle(
"Primary Vertices Z [cm]", 1);
1331 Vx_Z->setAxisTitle(
"Entries [#]", 2);
1333 Vx_X_Fit = ibooker.
book1D(
"G - vertex x fit",
1334 "Primary Vertex X Distribution (For Fit)",
1338 Vx_Y_Fit = ibooker.
book1D(
"G - vertex y fit",
1339 "Primary Vertex Y Distribution (For Fit)",
1343 Vx_Z_Fit = ibooker.
book1D(
"G - vertex z fit",
1344 "Primary Vertex Z Distribution (For Fit)",
1349 Vx_X_Fit->setAxisTitle(
"Entries [#]", 2);
1350 Vx_Y_Fit->setAxisTitle(
"Primary Vertices Y [cm]", 1);
1351 Vx_Y_Fit->setAxisTitle(
"Entries [#]", 2);
1352 Vx_Z_Fit->setAxisTitle(
"Primary Vertices Z [cm]", 1);
1353 Vx_Z_Fit->setAxisTitle(
"Entries [#]", 2);
1355 Vx_X_Cum = ibooker.
book1D(
"I - vertex x cum",
1356 "Primary Vertex X Distribution (Cumulative)",
1360 Vx_Y_Cum = ibooker.
book1D(
"I - vertex y cum",
1361 "Primary Vertex Y Distribution (Cumulative)",
1365 Vx_Z_Cum = ibooker.
book1D(
"I - vertex z cum",
1366 "Primary Vertex Z Distribution (Cumulative)",
1371 Vx_X_Cum->setAxisTitle(
"Entries [#]", 2);
1372 Vx_Y_Cum->setAxisTitle(
"Primary Vertices Y [cm]", 1);
1373 Vx_Y_Cum->setAxisTitle(
"Entries [#]", 2);
1374 Vx_Z_Cum->setAxisTitle(
"Primary Vertices Z [cm]", 1);
1375 Vx_Z_Cum->setAxisTitle(
"Entries [#]", 2);
1380 "B - muY vs lumi",
"#mu_{y} vs. Lumisection", nLumiXaxisRange, 0.5, ((
double)nLumiXaxisRange) + 0.5);
1382 "B - muZ vs lumi",
"#mu_{z} vs. Lumisection", nLumiXaxisRange, 0.5, ((
double)nLumiXaxisRange) + 0.5);
1384 mXlumi->setAxisTitle(
"#mu_{x} [cm]", 2);
1385 mXlumi->getTH1()->SetOption(
"E1");
1386 mYlumi->setAxisTitle(
"Lumisection [#]", 1);
1387 mYlumi->setAxisTitle(
"#mu_{y} [cm]", 2);
1388 mYlumi->getTH1()->SetOption(
"E1");
1389 mZlumi->setAxisTitle(
"Lumisection [#]", 1);
1390 mZlumi->setAxisTitle(
"#mu_{z} [cm]", 2);
1391 mZlumi->getTH1()->SetOption(
"E1");
1394 "C - sigmaX vs lumi",
"#sigma_{x} vs. Lumisection", nLumiXaxisRange, 0.5, ((
double)nLumiXaxisRange) + 0.5);
1396 "C - sigmaY vs lumi",
"#sigma_{y} vs. Lumisection", nLumiXaxisRange, 0.5, ((
double)nLumiXaxisRange) + 0.5);
1398 "C - sigmaZ vs lumi",
"#sigma_{z} vs. Lumisection", nLumiXaxisRange, 0.5, ((
double)nLumiXaxisRange) + 0.5);
1400 sXlumi->setAxisTitle(
"#sigma_{x} [cm]", 2);
1401 sXlumi->getTH1()->SetOption(
"E1");
1402 sYlumi->setAxisTitle(
"Lumisection [#]", 1);
1403 sYlumi->setAxisTitle(
"#sigma_{y} [cm]", 2);
1404 sYlumi->getTH1()->SetOption(
"E1");
1405 sZlumi->setAxisTitle(
"Lumisection [#]", 1);
1406 sZlumi->setAxisTitle(
"#sigma_{z} [cm]", 2);
1407 sZlumi->getTH1()->SetOption(
"E1");
1409 dxdzlumi = ibooker.
book1D(
1410 "D - dxdz vs lumi",
"dX/dZ vs. Lumisection", nLumiXaxisRange, 0.5, ((
double)nLumiXaxisRange) + 0.5);
1411 dydzlumi = ibooker.
book1D(
1412 "D - dydz vs lumi",
"dY/dZ vs. Lumisection", nLumiXaxisRange, 0.5, ((
double)nLumiXaxisRange) + 0.5);
1414 dxdzlumi->setAxisTitle(
"dX/dZ [rad]", 2);
1415 dxdzlumi->getTH1()->SetOption(
"E1");
1416 dydzlumi->setAxisTitle(
"Lumisection [#]", 1);
1417 dydzlumi->setAxisTitle(
"dY/dZ [rad]", 2);
1418 dydzlumi->getTH1()->SetOption(
"E1");
1420 Vx_ZX = ibooker.
book2D(
"E - vertex zx",
1421 "Primary Vertex ZX Distribution",
1428 Vx_ZY = ibooker.
book2D(
"E - vertex zy",
1429 "Primary Vertex ZY Distribution",
1436 Vx_XY = ibooker.
book2D(
"E - vertex xy",
1437 "Primary Vertex XY Distribution",
1445 Vx_ZX->setAxisTitle(
"Primary Vertices X [cm]", 2);
1446 Vx_ZX->setAxisTitle(
"Entries [#]", 3);
1447 Vx_ZY->setAxisTitle(
"Primary Vertices Z [cm]", 1);
1448 Vx_ZY->setAxisTitle(
"Primary Vertices Y [cm]", 2);
1449 Vx_ZY->setAxisTitle(
"Entries [#]", 3);
1450 Vx_XY->setAxisTitle(
"Primary Vertices X [cm]", 1);
1451 Vx_XY->setAxisTitle(
"Primary Vertices Y [cm]", 2);
1452 Vx_XY->setAxisTitle(
"Entries [#]", 3);
1454 Vx_ZX_Cum = ibooker.
book2D(
"H - vertex zx cum",
1455 "Primary Vertex ZX Distribution (Cumulative)",
1462 Vx_ZY_Cum = ibooker.
book2D(
"H - vertex zy cum",
1463 "Primary Vertex ZY Distribution (Cumulative)",
1470 Vx_XY_Cum = ibooker.
book2D(
"H - vertex xy cum",
1471 "Primary Vertex XY Distribution (Cumulative)",
1479 Vx_ZX_Cum->setAxisTitle(
"Primary Vertices X [cm]", 2);
1480 Vx_ZX_Cum->setAxisTitle(
"Entries [#]", 3);
1481 Vx_ZY_Cum->setAxisTitle(
"Primary Vertices Z [cm]", 1);
1482 Vx_ZY_Cum->setAxisTitle(
"Primary Vertices Y [cm]", 2);
1483 Vx_ZY_Cum->setAxisTitle(
"Entries [#]", 3);
1484 Vx_XY_Cum->setAxisTitle(
"Primary Vertices X [cm]", 1);
1485 Vx_XY_Cum->setAxisTitle(
"Primary Vertices Y [cm]", 2);
1486 Vx_XY_Cum->setAxisTitle(
"Entries [#]", 3);
1488 hitCounter = ibooker.
book1D(
1489 "J - pixelHits vs lumi",
"# Pixel-Hits vs. Lumisection", nLumiXaxisRange, 0.5, ((
double)nLumiXaxisRange) + 0.5);
1491 hitCounter->setAxisTitle(
"Pixel-Hits [#]", 2);
1492 hitCounter->getTH1()->SetOption(
"E1");
1494 goodVxCounter = ibooker.
book1D(
"K - good vertices vs lumi",
1495 "# Good vertices vs. Lumisection",
1498 ((
double)nLumiXaxisRange) + 0.5);
1500 goodVxCounter->setAxisTitle(
"Good vertices [#]", 2);
1501 goodVxCounter->getTH1()->SetOption(
"E1");
1503 statusCounter = ibooker.
book1D(
1504 "L - status vs lumi",
"App. Status vs. Lumisection", nLumiXaxisRange, 0.5, ((
double)nLumiXaxisRange) + 0.5);
1506 statusCounter->getTH1()->SetOption(
"E1");
1507 statusCounter->getTH1()->GetYaxis()->Set(11, -5.5, 5.5);
1508 statusCounter->getTH1()->GetYaxis()->SetBinLabel(1,
"Max Lumi.");
1509 statusCounter->getTH1()->GetYaxis()->SetBinLabel(2,
"Neg. det.");
1510 statusCounter->getTH1()->GetYaxis()->SetBinLabel(3,
"Infinite err.");
1511 statusCounter->getTH1()->GetYaxis()->SetBinLabel(4,
"No vtx.");
1512 statusCounter->getTH1()->GetYaxis()->SetBinLabel(5,
"Infinite EDM");
1513 statusCounter->getTH1()->GetYaxis()->SetBinLabel(6,
"OK");
1514 statusCounter->getTH1()->GetYaxis()->SetBinLabel(7,
"MINUIT stat.");
1515 statusCounter->getTH1()->GetYaxis()->SetBinLabel(8,
"MINUIT stat.");
1516 statusCounter->getTH1()->GetYaxis()->SetBinLabel(9,
"MINUIT stat.");
1517 statusCounter->getTH1()->GetYaxis()->SetBinLabel(10,
"MINUIT stat.");
1518 statusCounter->getTH1()->GetYaxis()->SetBinLabel(11,
"MINUIT stat.");
1520 fitResults = ibooker.
book2D(
"A - fit results",
"Results of Beam Spot Fit", 2, 0., 2., 9, 0., 9.);
1522 fitResults->setBinLabel(9,
"X[cm]", 2);
1523 fitResults->setBinLabel(8,
"Y[cm]", 2);
1524 fitResults->setBinLabel(7,
"Z[cm]", 2);
1525 fitResults->setBinLabel(6,
"#sigma_{X}[cm]", 2);
1526 fitResults->setBinLabel(5,
"#sigma_{Y}[cm]", 2);
1527 fitResults->setBinLabel(4,
"#sigma_{Z}[cm]", 2);
1528 fitResults->setBinLabel(3,
"#frac{dX}{dZ}[rad]", 2);
1529 fitResults->setBinLabel(2,
"#frac{dY}{dZ}[rad]", 2);
1530 fitResults->setBinLabel(1,
"Vtx[#]", 2);
1531 fitResults->setBinLabel(1,
"Value", 1);
1532 fitResults->setBinLabel(2,
"Error (stat)", 1);
1533 fitResults->getTH1()->SetOption(
"text");
1537 reportSummary = ibooker.
bookFloat(
"reportSummary");
1538 reportSummary->
Fill(-1);
1539 reportSummaryMap = ibooker.
book2D(
"reportSummaryMap",
"Pixel-Vertices Beam Spot: % Good Fits", 1, 0., 1., 1, 0., 1.);
1540 reportSummaryMap->
getTH1()->SetBinContent(1, 1, -1);
T getUntrackedParameter(std::string const &, T const &) const
MonitorElement * bookFloat(TString const &name, FUNC onbooking=NOOP())
static std::vector< std::string > checklist log
virtual void setCurrentFolder(std::string const &fullpath)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
constexpr bool isNotFinite(T x)
#define DEFINE_FWK_MODULE(type)
int MyFit(std::vector< double > *vals)
void analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) override
void dqmBeginLuminosityBlock(const edm::LuminosityBlock &lumiBlock, const edm::EventSetup &iSetup) override
void dqmEndLuminosityBlock(const edm::LuminosityBlock &lumiBlock, const edm::EventSetup &iSetup) override
data_type const * const_iterator
Timestamp const & beginTime() const
unsigned int HitCounter(const edm::Event &iEvent)
double Gauss3DFunc(const double *par)
void writeToFile(std::vector< double > *vals, edm::TimeValue_t BeginTimeOfFit, edm::TimeValue_t EndTimeOfFit, unsigned int BeginLumiOfFit, unsigned int EndLumiOfFit, int dataType)
LuminosityBlockNumber_t luminosityBlock() const
Vx3DHLTAnalyzer(const edm::ParameterSet &)
std::string formatTime(const time_t &t)
Timestamp const & endTime() const
LuminosityBlock const & getLuminosityBlock() const
void reset(std::string ResetType)
unsigned long long TimeValue_t
boost::transform_iterator< IterHelp, const_IdIter > const_iterator
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, FUNC onbooking=NOOP())
T getParameter(std::string const &) const
~Vx3DHLTAnalyzer() override
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
void printFitParams(const std::vector< double > &fitResults)
static std::atomic< unsigned int > counter
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
void reset(double vett[256])
TimeValue_t value() const
virtual TH1 * getTH1() const
Power< A, B >::type pow(const A &a, const B &b)
void writeToFile(const CaloCluster &c, FILE *file)
virtual void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)