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();
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);
103 dqmBeginLuminosityBlock(
iEvent.getLuminosityBlock(), iSetup);
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();
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;
174 cout <<
"(i,j) --> " <<
i <<
"," <<
j <<
" --> " << MyVertex.
Covariance[
i][
j] << endl;
176 }
else if (internalDebug ==
true)
177 cout <<
"[Vx3DHLTAnalyzer]::\tVertex discarded !" << endl;
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++) {
224 (std::fabs(
Vertices[
i].z - zPos) <= maxLongLength)) {
225 if (considerVxCovariance ==
true) {
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;
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.;
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.;
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.;
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.;
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.;
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)))
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;
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 Vx_X_Fit->getTH1()->Fit(myGaussFit,
"QRL");
1244 myGaussFit->SetParameter(0, Vx_Y_Fit->getTH1()->GetMaximum());
1245 myGaussFit->SetParameter(1, Vx_Y_Fit->getTH1()->GetMean());
1246 myGaussFit->SetParameter(2, Vx_Y_Fit->getTH1()->GetRMS());
1247 minXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(1);
1248 for (
int i = 0;
i < Vx_Y_Fit->getTH1()->GetNbinsX();
i++) {
1249 if (Vx_Y_Fit->getTH1()->GetBinContent(
i + 1) > 0) {
1250 minXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(
i + 1);
1254 maxXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(Vx_Y_Fit->getTH1()->GetNbinsX());
1255 for (
int i = Vx_Y_Fit->getTH1()->GetNbinsX();
i > 0;
i--) {
1256 if (Vx_Y_Fit->getTH1()->GetBinContent(
i) > 0) {
1257 maxXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(
i);
1261 myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1262 Vx_Y_Fit->getTH1()->Fit(myGaussFit,
"QRL");
1264 myGaussFit->SetParameter(0, Vx_Z_Fit->getTH1()->GetMaximum());
1265 myGaussFit->SetParameter(1, Vx_Z_Fit->getTH1()->GetMean());
1266 myGaussFit->SetParameter(2, Vx_Z_Fit->getTH1()->GetRMS());
1267 minXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(1);
1268 for (
int i = 0;
i < Vx_Z_Fit->getTH1()->GetNbinsX();
i++) {
1269 if (Vx_Z_Fit->getTH1()->GetBinContent(
i + 1) > 0) {
1270 minXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(
i + 1);
1274 maxXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(Vx_Z_Fit->getTH1()->GetNbinsX());
1275 for (
int i = Vx_Z_Fit->getTH1()->GetNbinsX();
i > 0;
i--) {
1276 if (Vx_Z_Fit->getTH1()->GetBinContent(
i) > 0) {
1277 maxXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(
i);
1281 myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1282 Vx_Z_Fit->getTH1()->Fit(myGaussFit,
"QRL");
1288 fitResults->setAxisTitle(
histTitle.str(), 1);
1289 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true)) {
1290 outputDebugFile <<
"\n"
1292 outputDebugFile <<
"BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) <<
" " << (beginTimeOfFit >> 32) << endl;
1293 outputDebugFile <<
"BeginLumiRange " << beginLumiOfFit << endl;
1294 outputDebugFile <<
histTitle.str().c_str() <<
"\n" << endl;
1297 histTitle <<
"Ongoing: no ongoing fits";
1298 fitResults->setAxisTitle(
histTitle.str(), 1);
1299 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true))
1300 outputDebugFile <<
histTitle.str().c_str() <<
"\n" << endl;
1304 hitCounter->getTH1()->SetBinContent(endLumiOfFit, (
double)totalHits);
1305 hitCounter->getTH1()->SetBinError(endLumiOfFit,
std::sqrt((
double)totalHits));
1310 if (internalDebug ==
true)
1311 cout <<
"[Vx3DHLTAnalyzer]::\tHistogram title: " <<
histTitle.str() << endl;
1324 Vx_X->setAxisTitle(
"Entries [#]", 2);
1325 Vx_Y->setAxisTitle(
"Primary Vertices Y [cm]", 1);
1326 Vx_Y->setAxisTitle(
"Entries [#]", 2);
1327 Vx_Z->setAxisTitle(
"Primary Vertices Z [cm]", 1);
1328 Vx_Z->setAxisTitle(
"Entries [#]", 2);
1330 Vx_X_Fit = ibooker.
book1D(
"G - vertex x fit",
1331 "Primary Vertex X Distribution (For Fit)",
1335 Vx_Y_Fit = ibooker.
book1D(
"G - vertex y fit",
1336 "Primary Vertex Y Distribution (For Fit)",
1340 Vx_Z_Fit = ibooker.
book1D(
"G - vertex z fit",
1341 "Primary Vertex Z Distribution (For Fit)",
1346 Vx_X_Fit->setAxisTitle(
"Entries [#]", 2);
1347 Vx_Y_Fit->setAxisTitle(
"Primary Vertices Y [cm]", 1);
1348 Vx_Y_Fit->setAxisTitle(
"Entries [#]", 2);
1349 Vx_Z_Fit->setAxisTitle(
"Primary Vertices Z [cm]", 1);
1350 Vx_Z_Fit->setAxisTitle(
"Entries [#]", 2);
1352 Vx_X_Cum = ibooker.
book1D(
"I - vertex x cum",
1353 "Primary Vertex X Distribution (Cumulative)",
1357 Vx_Y_Cum = ibooker.
book1D(
"I - vertex y cum",
1358 "Primary Vertex Y Distribution (Cumulative)",
1362 Vx_Z_Cum = ibooker.
book1D(
"I - vertex z cum",
1363 "Primary Vertex Z Distribution (Cumulative)",
1368 Vx_X_Cum->setAxisTitle(
"Entries [#]", 2);
1369 Vx_Y_Cum->setAxisTitle(
"Primary Vertices Y [cm]", 1);
1370 Vx_Y_Cum->setAxisTitle(
"Entries [#]", 2);
1371 Vx_Z_Cum->setAxisTitle(
"Primary Vertices Z [cm]", 1);
1372 Vx_Z_Cum->setAxisTitle(
"Entries [#]", 2);
1381 mXlumi->setAxisTitle(
"#mu_{x} [cm]", 2);
1382 mXlumi->getTH1()->SetOption(
"E1");
1383 mYlumi->setAxisTitle(
"Lumisection [#]", 1);
1384 mYlumi->setAxisTitle(
"#mu_{y} [cm]", 2);
1385 mYlumi->getTH1()->SetOption(
"E1");
1386 mZlumi->setAxisTitle(
"Lumisection [#]", 1);
1387 mZlumi->setAxisTitle(
"#mu_{z} [cm]", 2);
1388 mZlumi->getTH1()->SetOption(
"E1");
1397 sXlumi->setAxisTitle(
"#sigma_{x} [cm]", 2);
1398 sXlumi->getTH1()->SetOption(
"E1");
1399 sYlumi->setAxisTitle(
"Lumisection [#]", 1);
1400 sYlumi->setAxisTitle(
"#sigma_{y} [cm]", 2);
1401 sYlumi->getTH1()->SetOption(
"E1");
1402 sZlumi->setAxisTitle(
"Lumisection [#]", 1);
1403 sZlumi->setAxisTitle(
"#sigma_{z} [cm]", 2);
1404 sZlumi->getTH1()->SetOption(
"E1");
1406 dxdzlumi = ibooker.
book1D(
1408 dydzlumi = ibooker.
book1D(
1411 dxdzlumi->setAxisTitle(
"dX/dZ [rad]", 2);
1412 dxdzlumi->getTH1()->SetOption(
"E1");
1413 dydzlumi->setAxisTitle(
"Lumisection [#]", 1);
1414 dydzlumi->setAxisTitle(
"dY/dZ [rad]", 2);
1415 dydzlumi->getTH1()->SetOption(
"E1");
1417 Vx_ZX = ibooker.
book2D(
"E - vertex zx",
1418 "Primary Vertex ZX Distribution",
1425 Vx_ZY = ibooker.
book2D(
"E - vertex zy",
1426 "Primary Vertex ZY Distribution",
1433 Vx_XY = ibooker.
book2D(
"E - vertex xy",
1434 "Primary Vertex XY Distribution",
1442 Vx_ZX->setAxisTitle(
"Primary Vertices X [cm]", 2);
1443 Vx_ZX->setAxisTitle(
"Entries [#]", 3);
1444 Vx_ZY->setAxisTitle(
"Primary Vertices Z [cm]", 1);
1445 Vx_ZY->setAxisTitle(
"Primary Vertices Y [cm]", 2);
1446 Vx_ZY->setAxisTitle(
"Entries [#]", 3);
1447 Vx_XY->setAxisTitle(
"Primary Vertices X [cm]", 1);
1448 Vx_XY->setAxisTitle(
"Primary Vertices Y [cm]", 2);
1449 Vx_XY->setAxisTitle(
"Entries [#]", 3);
1451 Vx_ZX_Cum = ibooker.
book2D(
"H - vertex zx cum",
1452 "Primary Vertex ZX Distribution (Cumulative)",
1459 Vx_ZY_Cum = ibooker.
book2D(
"H - vertex zy cum",
1460 "Primary Vertex ZY Distribution (Cumulative)",
1467 Vx_XY_Cum = ibooker.
book2D(
"H - vertex xy cum",
1468 "Primary Vertex XY Distribution (Cumulative)",
1476 Vx_ZX_Cum->setAxisTitle(
"Primary Vertices X [cm]", 2);
1477 Vx_ZX_Cum->setAxisTitle(
"Entries [#]", 3);
1478 Vx_ZY_Cum->setAxisTitle(
"Primary Vertices Z [cm]", 1);
1479 Vx_ZY_Cum->setAxisTitle(
"Primary Vertices Y [cm]", 2);
1480 Vx_ZY_Cum->setAxisTitle(
"Entries [#]", 3);
1481 Vx_XY_Cum->setAxisTitle(
"Primary Vertices X [cm]", 1);
1482 Vx_XY_Cum->setAxisTitle(
"Primary Vertices Y [cm]", 2);
1483 Vx_XY_Cum->setAxisTitle(
"Entries [#]", 3);
1485 hitCounter = ibooker.
book1D(
1488 hitCounter->setAxisTitle(
"Pixel-Hits [#]", 2);
1489 hitCounter->getTH1()->SetOption(
"E1");
1491 goodVxCounter = ibooker.
book1D(
"K - good vertices vs lumi",
1492 "# Good vertices vs. Lumisection",
1497 goodVxCounter->setAxisTitle(
"Good vertices [#]", 2);
1498 goodVxCounter->getTH1()->SetOption(
"E1");
1500 statusCounter = ibooker.
book1D(
1503 statusCounter->getTH1()->SetOption(
"E1");
1504 statusCounter->getTH1()->GetYaxis()->Set(11, -5.5, 5.5);
1505 statusCounter->getTH1()->GetYaxis()->SetBinLabel(1,
"Max Lumi.");
1506 statusCounter->getTH1()->GetYaxis()->SetBinLabel(2,
"Neg. det.");
1507 statusCounter->getTH1()->GetYaxis()->SetBinLabel(3,
"Infinite err.");
1508 statusCounter->getTH1()->GetYaxis()->SetBinLabel(4,
"No vtx.");
1509 statusCounter->getTH1()->GetYaxis()->SetBinLabel(5,
"Infinite EDM");
1510 statusCounter->getTH1()->GetYaxis()->SetBinLabel(6,
"OK");
1511 statusCounter->getTH1()->GetYaxis()->SetBinLabel(7,
"MINUIT stat.");
1512 statusCounter->getTH1()->GetYaxis()->SetBinLabel(8,
"MINUIT stat.");
1513 statusCounter->getTH1()->GetYaxis()->SetBinLabel(9,
"MINUIT stat.");
1514 statusCounter->getTH1()->GetYaxis()->SetBinLabel(10,
"MINUIT stat.");
1515 statusCounter->getTH1()->GetYaxis()->SetBinLabel(11,
"MINUIT stat.");
1517 fitResults = ibooker.
book2D(
"A - fit results",
"Results of Beam Spot Fit", 2, 0., 2., 9, 0., 9.);
1519 fitResults->setBinLabel(9,
"X[cm]", 2);
1520 fitResults->setBinLabel(8,
"Y[cm]", 2);
1521 fitResults->setBinLabel(7,
"Z[cm]", 2);
1522 fitResults->setBinLabel(6,
"#sigma_{X}[cm]", 2);
1523 fitResults->setBinLabel(5,
"#sigma_{Y}[cm]", 2);
1524 fitResults->setBinLabel(4,
"#sigma_{Z}[cm]", 2);
1525 fitResults->setBinLabel(3,
"#frac{dX}{dZ}[rad]", 2);
1526 fitResults->setBinLabel(2,
"#frac{dY}{dZ}[rad]", 2);
1527 fitResults->setBinLabel(1,
"Vtx[#]", 2);
1528 fitResults->setBinLabel(1,
"Value", 1);
1529 fitResults->setBinLabel(2,
"Error (stat)", 1);
1530 fitResults->getTH1()->SetOption(
"text");
1534 reportSummary = ibooker.
bookFloat(
"reportSummary");
1535 reportSummary->
Fill(-1);
1536 reportSummaryMap = ibooker.
book2D(
"reportSummaryMap",
"Pixel-Vertices Beam Spot: % Good Fits", 1, 0., 1., 1, 0., 1.);
1537 reportSummaryMap->
getTH1()->SetBinContent(1, 1, -1);