15 #include <Math/Minimizer.h> 16 #include <Math/Factory.h> 17 #include <Math/Functor.h> 69 internalDebug =
false;
70 considerVxCovariance =
true;
71 pi = 3.141592653589793238;
90 stringstream debugFile;
93 if (outputDebugFile.is_open() ==
true)
94 outputDebugFile.close();
96 debugFile <<
tmp.c_str() <<
"_Run" <<
iEvent.id().run() <<
".txt";
97 outputDebugFile.open(debugFile.str().c_str(),
ios::out);
98 outputDebugFile.close();
99 outputDebugFile.open(debugFile.str().c_str(), ios::app);
102 dqmBeginLuminosityBlock(
iEvent.getLuminosityBlock(), iSetup);
103 }
else if (beginTimeOfFit != 0) {
104 totalHits += HitCounter(
iEvent);
106 if (internalDebug ==
true) {
107 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tI found " << totalHits <<
" pixel hits until now";
108 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tIn this event there are " << Vx3DCollection->size() <<
" vertex cadidates";
111 for (vector<Vertex>::const_iterator it3DVx = Vx3DCollection->begin(); it3DVx != Vx3DCollection->end(); it3DVx++) {
112 if (internalDebug ==
true) {
113 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tVertex selections:";
115 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tVertex number = " << it3DVx - Vx3DCollection->begin();
116 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tisValid = " << it3DVx->isValid();
117 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tisFake = " << it3DVx->isFake();
118 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tnodof = " << it3DVx->ndof();
119 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\ttracksSize = " << it3DVx->tracksSize();
122 if ((it3DVx->isValid() ==
true) && (it3DVx->isFake() ==
false) && (it3DVx->ndof() >=
minVxDoF) &&
123 (it3DVx->tracksSize() > 0) && ((it3DVx->ndof() + 3.) / ((double)it3DVx->tracksSize()) >= 2. *
minVxWgt)) {
124 for (
i = 0;
i <
DIM;
i++) {
125 for (
j = 0;
j <
DIM;
j++) {
144 if ((
i ==
DIM) && (det > 0.)) {
145 if (internalDebug ==
true)
146 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tVertex accepted !";
148 MyVertex.
x = it3DVx->x();
149 MyVertex.
y = it3DVx->y();
150 MyVertex.
z = it3DVx->z();
153 Vx_X->Fill(it3DVx->x());
154 Vx_Y->Fill(it3DVx->y());
155 Vx_Z->Fill(it3DVx->z());
157 Vx_ZX->Fill(it3DVx->z(), it3DVx->x());
158 Vx_ZY->Fill(it3DVx->z(), it3DVx->y());
159 Vx_XY->Fill(it3DVx->x(), it3DVx->y());
161 Vx_X_Cum->Fill(it3DVx->x());
162 Vx_Y_Cum->Fill(it3DVx->y());
163 Vx_Z_Cum->Fill(it3DVx->z());
165 Vx_ZX_Cum->Fill(it3DVx->z(), it3DVx->x());
166 Vx_ZY_Cum->Fill(it3DVx->z(), it3DVx->y());
167 Vx_XY_Cum->Fill(it3DVx->x(), it3DVx->y());
168 }
else if (internalDebug ==
true) {
169 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tVertex discarded !";
175 }
else if (internalDebug ==
true)
176 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tVertex discarded !";
196 strftime(ts,
sizeof(ts),
"%Y.%m.%d %H:%M:%S %Z", gmtime(&
t));
198 string ts_string(ts);
220 for (
unsigned int i = 0;
i <
Vertices.size();
i++) {
223 (std::fabs(
Vertices[
i].z - zPos) <= maxLongLength)) {
224 if (considerVxCovariance ==
true) {
229 K[1][2] = K[2][1] = par[4] * (std::fabs(par[2]) - std::fabs(par[1])) - par[5] * par[3] +
231 K[0][2] = K[2][0] = par[5] * (std::fabs(par[2]) - std::fabs(par[0])) - par[4] * par[3] +
234 K[0][0] = std::fabs(par[0]);
235 K[1][1] = std::fabs(par[1]);
236 K[2][2] = std::fabs(par[2]);
237 K[0][1] = K[1][0] = par[3];
238 K[1][2] = K[2][1] = par[4] * (std::fabs(par[2]) - std::fabs(par[1])) - par[5] * par[3];
239 K[0][2] = K[2][0] = par[5] * (std::fabs(par[2]) - std::fabs(par[0])) - par[4] * par[3];
242 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]) +
243 K[0][2] * (K[0][1] * K[1][2] - K[0][2] * K[1][1]);
245 M[0][0] = (K[1][1] * K[2][2] - K[1][2] * K[1][2]) / det;
246 M[1][1] = (K[0][0] * K[2][2] - K[0][2] * K[0][2]) / det;
247 M[2][2] = (K[0][0] * K[1][1] - K[0][1] * K[0][1]) / det;
248 M[0][1] = M[1][0] = (K[0][2] * K[1][2] - K[0][1] * K[2][2]) / det;
249 M[1][2] = M[2][1] = (K[0][2] * K[0][1] - K[1][2] * K[0][0]) / det;
250 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 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\t@@@ START FITTING @@@";
318 for (
int i = 0;
i < 3;
i++) {
320 if (internalDebug ==
true)
321 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tdeltaMean --> " << deltaMean;
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);
336 xPos = *(
it + 6) + deltaMean;
341 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(*(
it + 0)) + std::fabs(*(
it + 1)) / 2.);
346 goodData = Gauss3D->Status();
347 edm = Gauss3D->Edm();
349 edm::LogError(
"Vx3DHLTAnalyzer") <<
"\tCaught Minuit2 exception @ Fit at X: \n" << er.
what();
359 if (internalDebug ==
true)
360 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNot finite edm !";
361 }
else if (goodData != -6)
362 for (
unsigned int j = 0;
j < nParams;
j++)
365 if (internalDebug ==
true)
366 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNot finite errors !";
370 covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
371 Gauss3D->X()[5] * Gauss3D->X()[3];
372 covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
373 Gauss3D->X()[4] * Gauss3D->X()[3];
375 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
376 Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
377 covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
380 if (internalDebug ==
true)
381 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNegative determinant !";
385 if ((goodData == 0) && (std::fabs(
edm) < bestEdm)) {
390 if (internalDebug ==
true)
391 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tFound bestMovementX --> " << bestMovementX;
395 for (
int i = 0;
i < 3;
i++) {
397 if (internalDebug ==
true) {
398 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tdeltaMean --> " << deltaMean;
399 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tdeltaMean X --> " << (double(bestMovementX) - 1.) *
std::sqrt(*(
it + 0));
404 Gauss3D->SetVariable(0,
"var x ", *(
it + 0), parDistanceXY * parDistanceXY);
405 Gauss3D->SetVariable(1,
"var y ", *(
it + 1), parDistanceXY * parDistanceXY);
406 Gauss3D->SetVariable(2,
"var z ", *(
it + 2), parDistanceZ * parDistanceZ);
407 Gauss3D->SetVariable(3,
"cov xy", *(
it + 3), parDistanceCxy);
408 Gauss3D->SetVariable(4,
"dydz ", *(
it + 4), parDistanceddZ);
409 Gauss3D->SetVariable(5,
"dxdz ", *(
it + 5), parDistanceddZ);
410 Gauss3D->SetVariable(6,
"mean x", *(
it + 6) + (
double(bestMovementX) - 1.) *
std::sqrt(*(
it + 0)), parDistanceXY);
411 Gauss3D->SetVariable(7,
"mean y", *(
it + 7) + deltaMean, parDistanceXY);
412 Gauss3D->SetVariable(8,
"mean z", *(
it + 8), parDistanceZ);
415 xPos = *(
it + 6) + (
double(bestMovementX) - 1.) *
std::sqrt(*(
it + 0));
416 yPos = *(
it + 7) + deltaMean;
420 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(*(
it + 0)) + std::fabs(*(
it + 1)) / 2.);
425 goodData = Gauss3D->Status();
426 edm = Gauss3D->Edm();
428 edm::LogError(
"Vx3DHLTAnalyzer") <<
"\tCaught Minuit2 exception @ Fit at Y: \n" << er.
what();
438 if (internalDebug ==
true)
439 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNot finite edm !";
440 }
else if (goodData != -6)
441 for (
unsigned int j = 0;
j < nParams;
j++)
444 if (internalDebug ==
true)
445 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNot finite errors !";
449 covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
450 Gauss3D->X()[5] * Gauss3D->X()[3];
451 covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
452 Gauss3D->X()[4] * Gauss3D->X()[3];
454 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
455 Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
456 covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
459 if (internalDebug ==
true)
460 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNegative determinant !";
464 if ((goodData == 0) && (std::fabs(
edm) < bestEdm)) {
469 if (internalDebug ==
true)
470 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tFound bestMovementY --> " << bestMovementY;
474 for (
int i = 0;
i < 3;
i++) {
476 if (internalDebug ==
true) {
477 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tdeltaMean --> " << deltaMean;
478 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tdeltaMean X --> " << (double(bestMovementX) - 1.) *
std::sqrt(*(
it + 0));
479 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tdeltaMean Y --> " << (double(bestMovementY) - 1.) *
std::sqrt(*(
it + 1));
484 Gauss3D->SetVariable(0,
"var x ", *(
it + 0), parDistanceXY * parDistanceXY);
485 Gauss3D->SetVariable(1,
"var y ", *(
it + 1), parDistanceXY * parDistanceXY);
486 Gauss3D->SetVariable(2,
"var z ", *(
it + 2), parDistanceZ * parDistanceZ);
487 Gauss3D->SetVariable(3,
"cov xy", *(
it + 3), parDistanceCxy);
488 Gauss3D->SetVariable(4,
"dydz ", *(
it + 4), parDistanceddZ);
489 Gauss3D->SetVariable(5,
"dxdz ", *(
it + 5), parDistanceddZ);
490 Gauss3D->SetVariable(6,
"mean x", *(
it + 6) + (
double(bestMovementX) - 1.) *
std::sqrt(*(
it + 0)), parDistanceXY);
491 Gauss3D->SetVariable(7,
"mean y", *(
it + 7) + (
double(bestMovementY) - 1.) *
std::sqrt(*(
it + 1)), parDistanceXY);
492 Gauss3D->SetVariable(8,
"mean z", *(
it + 8) + deltaMean, parDistanceZ);
495 xPos = *(
it + 6) + (
double(bestMovementX) - 1.) *
std::sqrt(*(
it + 0));
496 yPos = *(
it + 7) + (
double(bestMovementY) - 1.) *
std::sqrt(*(
it + 1));
497 zPos = *(
it + 8) + deltaMean;
500 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(*(
it + 0)) + std::fabs(*(
it + 1)) / 2.);
505 goodData = Gauss3D->Status();
506 edm = Gauss3D->Edm();
508 edm::LogError(
"Vx3DHLTAnalyzer") <<
"\tCaught Minuit2 exception @ Fit at Z: \n" << er.
what();
518 if (internalDebug ==
true)
519 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNot finite edm !";
520 }
else if (goodData != -6)
521 for (
unsigned int j = 0;
j < nParams;
j++)
524 if (internalDebug ==
true)
525 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNot finite errors !";
529 covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
530 Gauss3D->X()[5] * Gauss3D->X()[3];
531 covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
532 Gauss3D->X()[4] * Gauss3D->X()[3];
534 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
535 Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
536 covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
539 if (internalDebug ==
true)
540 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNegative determinant !";
544 if ((goodData == 0) && (std::fabs(
edm) < bestEdm)) {
549 if (internalDebug ==
true)
550 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tFound bestMovementZ --> " << bestMovementZ;
555 Gauss3D->SetVariable(0,
"var x ", *(
it + 0), parDistanceXY * parDistanceXY);
556 Gauss3D->SetVariable(1,
"var y ", *(
it + 1), parDistanceXY * parDistanceXY);
557 Gauss3D->SetVariable(2,
"var z ", *(
it + 2), parDistanceZ * parDistanceZ);
558 Gauss3D->SetVariable(3,
"cov xy", *(
it + 3), parDistanceCxy);
559 Gauss3D->SetVariable(4,
"dydz ", *(
it + 4), parDistanceddZ);
560 Gauss3D->SetVariable(5,
"dxdz ", *(
it + 5), parDistanceddZ);
561 Gauss3D->SetVariable(6,
"mean x", *(
it + 6) + (
double(bestMovementX) - 1.) *
std::sqrt(*(
it + 0)), parDistanceXY);
562 Gauss3D->SetVariable(7,
"mean y", *(
it + 7) + (
double(bestMovementY) - 1.) *
std::sqrt(*(
it + 1)), parDistanceXY);
563 Gauss3D->SetVariable(8,
"mean z", *(
it + 8) + (
double(bestMovementZ) - 1.) *
std::sqrt(*(
it + 2)), parDistanceZ);
566 xPos = *(
it + 6) + (
double(bestMovementX) - 1.) *
std::sqrt(*(
it + 0));
567 yPos = *(
it + 7) + (
double(bestMovementY) - 1.) *
std::sqrt(*(
it + 1));
568 zPos = *(
it + 8) + (
double(bestMovementZ) - 1.) *
std::sqrt(*(
it + 2));
571 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(*(
it + 0)) + std::fabs(*(
it + 1)) / 2.);
576 goodData = Gauss3D->Status();
577 edm = Gauss3D->Edm();
579 edm::LogError(
"Vx3DHLTAnalyzer") <<
"\tCaught Minuit2 exception @ Final fit: \n" << er.
what();
588 if (internalDebug ==
true)
589 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNot finite edm !";
590 }
else if (goodData != -6)
591 for (
unsigned int j = 0;
j < nParams;
j++)
594 if (internalDebug ==
true)
595 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNot finite errors !";
599 covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
600 Gauss3D->X()[5] * Gauss3D->X()[3];
601 covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
602 Gauss3D->X()[4] * Gauss3D->X()[3];
604 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
605 Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
606 covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
609 if (internalDebug ==
true)
610 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNegative determinant !";
615 for (
unsigned int i = 0;
i < trials;
i++) {
616 if ((goodData != 0) && (goodData != -2)) {
619 if (internalDebug ==
true)
620 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tFIT WITH DIFFERENT PARAMETER DISTANCES - STEP " <<
i + 1;
622 Gauss3D->SetVariable(0,
"var x ", *(
it + 0), parDistanceXY * parDistanceXY * largerDist[
i]);
623 Gauss3D->SetVariable(1,
"var y ", *(
it + 1), parDistanceXY * parDistanceXY * largerDist[
i]);
624 Gauss3D->SetVariable(2,
"var z ", *(
it + 2), parDistanceZ * parDistanceZ * largerDist[
i]);
625 Gauss3D->SetVariable(3,
"cov xy", *(
it + 3), parDistanceCxy * largerDist[
i]);
626 Gauss3D->SetVariable(4,
"dydz ", *(
it + 4), parDistanceddZ * largerDist[
i]);
627 Gauss3D->SetVariable(5,
"dxdz ", *(
it + 5), parDistanceddZ * largerDist[
i]);
628 Gauss3D->SetVariable(6,
630 *(
it + 6) + (
double(bestMovementX) - 1.) *
std::sqrt(*(
it + 0)),
631 parDistanceXY * largerDist[
i]);
632 Gauss3D->SetVariable(7,
634 *(
it + 7) + (
double(bestMovementY) - 1.) *
std::sqrt(*(
it + 1)),
635 parDistanceXY * largerDist[
i]);
636 Gauss3D->SetVariable(
637 8,
"mean z", *(
it + 8) + (
double(bestMovementZ) - 1.) *
std::sqrt(*(
it + 2)), parDistanceZ * largerDist[
i]);
640 xPos = *(
it + 6) + (
double(bestMovementX) - 1.) *
std::sqrt(*(
it + 0));
641 yPos = *(
it + 7) + (
double(bestMovementY) - 1.) *
std::sqrt(*(
it + 1));
642 zPos = *(
it + 8) + (
double(bestMovementZ) - 1.) *
std::sqrt(*(
it + 2));
645 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(*(
it + 0)) + std::fabs(*(
it + 1)) / 2.);
650 goodData = Gauss3D->Status();
651 edm = Gauss3D->Edm();
653 edm::LogError(
"Vx3DHLTAnalyzer") <<
"\tCaught Minuit2 exception @ Fit with different distances: \n" 664 if (internalDebug ==
true)
665 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNot finite edm !";
666 }
else if (goodData != -6)
667 for (
unsigned int j = 0;
j < nParams;
j++)
670 if (internalDebug ==
true)
671 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNot finite errors !";
675 covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
676 Gauss3D->X()[5] * Gauss3D->X()[3];
677 covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
678 Gauss3D->X()[4] * Gauss3D->X()[3];
680 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
681 Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
682 covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
685 if (internalDebug ==
true)
686 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNegative determinant !";
694 for (
unsigned int i = 0;
i < nParams;
i++) {
695 vals->operator[](
i) = Gauss3D->X()[
i];
696 vals->operator[](
i + nParams) = Gauss3D->Errors()[
i];
707 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true)) {
708 outputDebugFile <<
"Runnumber " <<
runNumber << endl;
709 outputDebugFile <<
"BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) <<
" " << (beginTimeOfFit >> 32) << endl;
710 outputDebugFile <<
"BeginLumiRange " << beginLumiOfFit << endl;
711 outputDebugFile <<
"EndTimeOfFit " << formatTime(endTimeOfFit >> 32) <<
" " << (endTimeOfFit >> 32) << endl;
712 outputDebugFile <<
"EndLumiRange " << endLumiOfFit << endl;
713 outputDebugFile <<
"LumiCounter " << lumiCounter << endl;
714 outputDebugFile <<
"LastLumiOfFit " << lastLumiOfFit << endl;
717 if (ResetType ==
"scratch") {
755 goodVxCounter->Reset();
756 statusCounter->Reset();
759 reportSummary->Fill(-1);
760 reportSummaryMap->getTH1()->SetBinContent(1, 1, -1);
771 if (internalDebug ==
true)
772 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tReset issued: scratch";
773 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true))
774 outputDebugFile <<
"Reset -scratch- issued\n" << endl;
775 }
else if (ResetType ==
"whole") {
793 if (internalDebug ==
true)
794 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tReset issued: whole";
795 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true))
796 outputDebugFile <<
"Reset -whole- issued\n" << endl;
797 }
else if (ResetType ==
"fit") {
802 if (internalDebug ==
true)
803 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tReset issued: fit";
804 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true))
805 outputDebugFile <<
"Reset -fit- issued\n" << endl;
806 }
else if (ResetType ==
"hitCounter") {
809 if (internalDebug ==
true)
810 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tReset issued: hitCounter";
811 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true))
812 outputDebugFile <<
"Reset -hitCounter- issued\n" << endl;
819 unsigned int BeginLumiOfFit,
820 unsigned int EndLumiOfFit,
822 stringstream BufferString;
823 BufferString.precision(5);
827 if ((
outputFile.is_open() ==
true) && (vals !=
nullptr) && (vals->size() == (nParams - 1) * 2)) {
828 vector<double>::const_iterator
it = vals->begin();
831 outputFile <<
"BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) <<
" " << (beginTimeOfFit >> 32) << endl;
832 outputFile <<
"EndTimeOfFit " << formatTime(endTimeOfFit >> 32) <<
" " << (endTimeOfFit >> 32) << endl;
833 outputFile <<
"LumiRange " << beginLumiOfFit <<
" - " << endLumiOfFit << endl;
839 BufferString << *(
it + 0);
840 outputFile <<
"X0 " << BufferString.str().c_str() << endl;
841 BufferString.str(
"");
843 BufferString << *(
it + 1);
844 outputFile <<
"Y0 " << BufferString.str().c_str() << endl;
845 BufferString.str(
"");
847 BufferString << *(
it + 2);
848 outputFile <<
"Z0 " << BufferString.str().c_str() << endl;
849 BufferString.str(
"");
851 BufferString << *(
it + 5);
852 outputFile <<
"sigmaZ0 " << BufferString.str().c_str() << endl;
853 BufferString.str(
"");
855 BufferString << *(
it + 6);
856 outputFile <<
"dxdz " << BufferString.str().c_str() << endl;
857 BufferString.str(
"");
859 BufferString << *(
it + 7);
860 outputFile <<
"dydz " << BufferString.str().c_str() << endl;
861 BufferString.str(
"");
863 BufferString << *(
it + 3);
864 outputFile <<
"BeamWidthX " << BufferString.str().c_str() << endl;
865 BufferString.str(
"");
867 BufferString << *(
it + 4);
868 outputFile <<
"BeamWidthY " << BufferString.str().c_str() << endl;
869 BufferString.str(
"");
871 outputFile <<
"Cov(0,j) " << *(
it + 8) <<
" 0 0 0 0 0 0" << endl;
872 outputFile <<
"Cov(1,j) 0 " << *(
it + 9) <<
" 0 0 0 0 0" << endl;
873 outputFile <<
"Cov(2,j) 0 0 " << *(
it + 10) <<
" 0 0 0 0" << endl;
874 outputFile <<
"Cov(3,j) 0 0 0 " << *(
it + 13) <<
" 0 0 0" << endl;
875 outputFile <<
"Cov(4,j) 0 0 0 0 " << *(
it + 14) <<
" 0 0" << endl;
876 outputFile <<
"Cov(5,j) 0 0 0 0 0 " << *(
it + 15) <<
" 0" << endl;
878 << ((*(
it + 11)) + (*(
it + 12)) + 2. *
std::sqrt((*(
it + 11)) * (*(
it + 12)))) / 4. << endl;
886 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true) && (vals !=
nullptr) &&
887 (vals->size() == (nParams - 1) * 2)) {
888 vector<double>::const_iterator
it = vals->begin();
890 outputDebugFile <<
"Runnumber " <<
runNumber << endl;
891 outputDebugFile <<
"BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) <<
" " << (beginTimeOfFit >> 32) << endl;
892 outputDebugFile <<
"EndTimeOfFit " << formatTime(endTimeOfFit >> 32) <<
" " << (endTimeOfFit >> 32) << endl;
893 outputDebugFile <<
"LumiRange " << beginLumiOfFit <<
" - " << endLumiOfFit << endl;
894 outputDebugFile <<
"Type " <<
dataType << endl;
899 BufferString << *(
it + 0);
900 outputDebugFile <<
"X0 " << BufferString.str().c_str() << endl;
901 BufferString.str(
"");
903 BufferString << *(
it + 1);
904 outputDebugFile <<
"Y0 " << BufferString.str().c_str() << endl;
905 BufferString.str(
"");
907 BufferString << *(
it + 2);
908 outputDebugFile <<
"Z0 " << BufferString.str().c_str() << endl;
909 BufferString.str(
"");
911 BufferString << *(
it + 5);
912 outputDebugFile <<
"sigmaZ0 " << BufferString.str().c_str() << endl;
913 BufferString.str(
"");
915 BufferString << *(
it + 6);
916 outputDebugFile <<
"dxdz " << BufferString.str().c_str() << endl;
917 BufferString.str(
"");
919 BufferString << *(
it + 7);
920 outputDebugFile <<
"dydz " << BufferString.str().c_str() << endl;
921 BufferString.str(
"");
923 BufferString << *(
it + 3);
924 outputDebugFile <<
"BeamWidthX " << BufferString.str().c_str() << endl;
925 BufferString.str(
"");
927 BufferString << *(
it + 4);
928 outputDebugFile <<
"BeamWidthY " << BufferString.str().c_str() << endl;
929 BufferString.str(
"");
931 outputDebugFile <<
"Cov(0,j) " << *(
it + 8) <<
" 0 0 0 0 0 0" << endl;
932 outputDebugFile <<
"Cov(1,j) 0 " << *(
it + 9) <<
" 0 0 0 0 0" << endl;
933 outputDebugFile <<
"Cov(2,j) 0 0 " << *(
it + 10) <<
" 0 0 0 0" << endl;
934 outputDebugFile <<
"Cov(3,j) 0 0 0 " << *(
it + 13) <<
" 0 0 0" << endl;
935 outputDebugFile <<
"Cov(4,j) 0 0 0 0 " << *(
it + 14) <<
" 0 0" << endl;
936 outputDebugFile <<
"Cov(5,j) 0 0 0 0 0 " << *(
it + 15) <<
" 0" << endl;
937 outputDebugFile <<
"Cov(6,j) 0 0 0 0 0 0 " 938 << ((*(
it + 11)) + (*(
it + 12)) + 2. *
std::sqrt((*(
it + 11)) * (*(
it + 12)))) / 4. << endl;
940 outputDebugFile <<
"EmittanceX 0" << endl;
941 outputDebugFile <<
"EmittanceY 0" << endl;
942 outputDebugFile <<
"BetaStar 0" << endl;
944 outputDebugFile <<
"\nUsed vertices: " << counterVx <<
"\n" << endl;
949 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"var x --> " << fitResults[0] <<
" +/- " << fitResults[0 + nParams];
950 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"var y --> " << fitResults[1] <<
" +/- " << fitResults[1 + nParams];
951 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"var z --> " << fitResults[2] <<
" +/- " << fitResults[2 + nParams];
952 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"cov xy --> " << fitResults[3] <<
" +/- " << fitResults[3 + nParams];
953 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"dydz --> " << fitResults[4] <<
" +/- " << fitResults[4 + nParams];
954 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"dxdz --> " << fitResults[5] <<
" +/- " << fitResults[5 + nParams];
955 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"mean x --> " << fitResults[6] <<
" +/- " << fitResults[6 + nParams];
956 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"mean y --> " << fitResults[7] <<
" +/- " << fitResults[7 + nParams];
957 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"mean z --> " << fitResults[8] <<
" +/- " << fitResults[8 + nParams];
962 if ((lumiCounter == 0) && (lumiBlock.
luminosityBlock() > lastLumiOfFit)) {
966 }
else if ((lumiCounter != 0) && (lumiBlock.
luminosityBlock() >= (beginLumiOfFit + lumiCounter)))
974 double minXfit, maxXfit;
980 lastLumiOfFit = endLumiOfFit;
983 hitCounter->getTH1()->SetBinContent(lastLumiOfFit, (
double)totalHits);
984 hitCounter->getTH1()->SetBinError(
990 vector<double> fitResults;
992 fitResults.push_back(Vx_X->getTH1()->GetRMS() * Vx_X->getTH1()->GetRMS());
993 fitResults.push_back(Vx_Y->getTH1()->GetRMS() * Vx_Y->getTH1()->GetRMS());
994 fitResults.push_back(Vx_Z->getTH1()->GetRMS() * Vx_Z->getTH1()->GetRMS());
995 fitResults.push_back(0.0);
996 fitResults.push_back(0.0);
997 fitResults.push_back(0.0);
998 fitResults.push_back(Vx_X->getTH1()->GetMean());
999 fitResults.push_back(Vx_Y->getTH1()->GetMean());
1000 fitResults.push_back(Vx_Z->getTH1()->GetMean());
1001 for (
unsigned int i = 0;
i < nParams;
i++)
1002 fitResults.push_back(0.0);
1004 if (internalDebug ==
true) {
1005 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\t@@@ Beam Spot parameters - prefit @@@";
1007 printFitParams(fitResults);
1010 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) <<
" " 1011 << (beginTimeOfFit >> 32);
1012 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"EndTimeOfFit " << formatTime(endTimeOfFit >> 32) <<
" " 1013 << (endTimeOfFit >> 32);
1014 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"LumiRange " << beginLumiOfFit <<
" - " << endLumiOfFit;
1017 goodData = MyFit(&fitResults);
1019 if (internalDebug ==
true) {
1020 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\t@@@ Beam Spot parameters - postfit @@@";
1022 printFitParams(fitResults);
1024 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"goodData --> " << goodData;
1025 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"Used vertices --> " << counterVx;
1028 if (goodData == 0) {
1029 vals.push_back(fitResults[6]);
1030 vals.push_back(fitResults[7]);
1031 vals.push_back(fitResults[8]);
1032 vals.push_back(
std::sqrt(std::fabs(fitResults[0])));
1033 vals.push_back(
std::sqrt(std::fabs(fitResults[1])));
1034 vals.push_back(
std::sqrt(std::fabs(fitResults[2])));
1035 vals.push_back(fitResults[5]);
1036 vals.push_back(fitResults[4]);
1038 vals.push_back(
std::pow(fitResults[6 + nParams], 2.));
1039 vals.push_back(
std::pow(fitResults[7 + nParams], 2.));
1040 vals.push_back(
std::pow(fitResults[8 + nParams], 2.));
1041 vals.push_back(
std::pow(std::fabs(fitResults[0 + nParams]) / (2. *
std::sqrt(std::fabs(fitResults[0]))), 2.));
1042 vals.push_back(
std::pow(std::fabs(fitResults[1 + nParams]) / (2. *
std::sqrt(std::fabs(fitResults[1]))), 2.));
1043 vals.push_back(
std::pow(std::fabs(fitResults[2 + nParams]) / (2. *
std::sqrt(std::fabs(fitResults[2]))), 2.));
1044 vals.push_back(
std::pow(fitResults[5 + nParams], 2.));
1045 vals.push_back(
std::pow(fitResults[4 + nParams], 2.));
1047 for (
unsigned int i = 0;
i < (nParams - 1) * 2;
i++)
1048 vals.push_back(0.0);
1052 counterVx = Vx_X->getTH1F()->GetEntries();
1054 if (Vx_X->getTH1F()->GetEntries() >=
minNentries) {
1057 vals.push_back(Vx_X->getTH1F()->GetMean());
1058 vals.push_back(Vx_Y->getTH1F()->GetMean());
1059 vals.push_back(Vx_Z->getTH1F()->GetMean());
1060 vals.push_back(Vx_X->getTH1F()->GetRMS());
1061 vals.push_back(Vx_Y->getTH1F()->GetRMS());
1062 vals.push_back(Vx_Z->getTH1F()->GetRMS());
1063 vals.push_back(0.0);
1064 vals.push_back(0.0);
1066 vals.push_back(
std::pow(Vx_X->getTH1F()->GetMeanError(), 2.));
1067 vals.push_back(
std::pow(Vx_Y->getTH1F()->GetMeanError(), 2.));
1068 vals.push_back(
std::pow(Vx_Z->getTH1F()->GetMeanError(), 2.));
1069 vals.push_back(
std::pow(Vx_X->getTH1F()->GetRMSError(), 2.));
1070 vals.push_back(
std::pow(Vx_Y->getTH1F()->GetRMSError(), 2.));
1071 vals.push_back(
std::pow(Vx_Z->getTH1F()->GetRMSError(), 2.));
1072 vals.push_back(0.0);
1073 vals.push_back(0.0);
1076 for (
unsigned int i = 0;
i < (nParams - 1) * 2;
i++)
1077 vals.push_back(0.0);
1100 writeToFile(&vals, beginTimeOfFit, endTimeOfFit, beginLumiOfFit, endLumiOfFit, 3);
1101 if (internalDebug ==
true)
1102 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tUsed vertices: " << counterVx;
1104 statusCounter->getTH1()->SetBinContent(lastLumiOfFit, (
double)goodData);
1105 statusCounter->getTH1()->SetBinError(lastLumiOfFit, 1
e-3);
1115 for (
int i = 0;
i < Vx_X_Fit->getTH1()->GetNbinsX();
i++) {
1116 Vx_X_Fit->getTH1()->SetBinContent(
1117 i + 1, Vx_X_Fit->getTH1()->GetBinContent(
i + 1) + Vx_X->getTH1()->GetBinContent(
i + 1));
1118 Vx_X_Fit->getTH1()->SetBinError(
i + 1,
sqrt(Vx_X_Fit->getTH1()->GetBinContent(
i + 1)));
1121 for (
int i = 0;
i < Vx_Y_Fit->getTH1()->GetNbinsX();
i++) {
1122 Vx_Y_Fit->getTH1()->SetBinContent(
1123 i + 1, Vx_Y_Fit->getTH1()->GetBinContent(
i + 1) + Vx_Y->getTH1()->GetBinContent(
i + 1));
1124 Vx_Y_Fit->getTH1()->SetBinError(
i + 1,
sqrt(Vx_Y_Fit->getTH1()->GetBinContent(
i + 1)));
1127 for (
int i = 0;
i < Vx_Z_Fit->getTH1()->GetNbinsX();
i++) {
1128 Vx_Z_Fit->getTH1()->SetBinContent(
1129 i + 1, Vx_Z_Fit->getTH1()->GetBinContent(
i + 1) + Vx_Z->getTH1()->GetBinContent(
i + 1));
1130 Vx_Z_Fit->getTH1()->SetBinError(
i + 1,
sqrt(Vx_Z_Fit->getTH1()->GetBinContent(
i + 1)));
1134 if (goodData == 0) {
1137 histTitle <<
"Ongoing: fitted lumis " << beginLumiOfFit <<
" - " << endLumiOfFit;
1146 statusCounter->getTH1()->SetBinContent(lastLumiOfFit, -5);
1147 statusCounter->getTH1()->SetBinError(lastLumiOfFit, 1
e-3);
1149 reset(
"hitCounter");
1152 reportSummary->Fill((numberFits != 0 ? ((
double)numberGoodFits) / ((
double)numberFits) : -1));
1153 reportSummaryMap->getTH1()->SetBinContent(
1154 1, 1, (numberFits != 0 ? ((
double)numberGoodFits) / ((
double)numberFits) : -1));
1156 fitResults->setAxisTitle(
histTitle.str(), 1);
1158 fitResults->setBinContent(1, 9, vals[0]);
1159 fitResults->setBinContent(1, 8, vals[1]);
1160 fitResults->setBinContent(1, 7, vals[2]);
1161 fitResults->setBinContent(1, 6, vals[3]);
1162 fitResults->setBinContent(1, 5, vals[4]);
1163 fitResults->setBinContent(1, 4, vals[5]);
1164 fitResults->setBinContent(1, 3, vals[6]);
1165 fitResults->setBinContent(1, 2, vals[7]);
1166 fitResults->setBinContent(1, 1, counterVx);
1168 fitResults->setBinContent(2, 9,
std::sqrt(vals[8]));
1169 fitResults->setBinContent(2, 8,
std::sqrt(vals[9]));
1170 fitResults->setBinContent(2, 7,
std::sqrt(vals[10]));
1171 fitResults->setBinContent(2, 6,
std::sqrt(vals[11]));
1172 fitResults->setBinContent(2, 5,
std::sqrt(vals[12]));
1173 fitResults->setBinContent(2, 4,
std::sqrt(vals[13]));
1174 fitResults->setBinContent(2, 3,
std::sqrt(vals[14]));
1175 fitResults->setBinContent(2, 2,
std::sqrt(vals[15]));
1176 fitResults->setBinContent(2, 1,
std::sqrt(counterVx));
1179 TF1* myLinFit =
new TF1(
1180 "myLinFit",
"[0] + [1]*x", mXlumi->getTH1()->GetXaxis()->GetXmin(), mXlumi->getTH1()->GetXaxis()->GetXmax());
1181 myLinFit->SetLineColor(2);
1182 myLinFit->SetLineWidth(2);
1183 myLinFit->SetParName(0,
"Inter.");
1184 myLinFit->SetParName(1,
"Slope");
1186 mXlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[0]);
1187 mXlumi->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt(vals[8]));
1188 myLinFit->SetParameter(0, mXlumi->getTH1()->GetMean(2));
1189 myLinFit->SetParameter(1, 0.0);
1190 mXlumi->getTH1()->Fit(myLinFit,
"QR");
1192 mYlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[1]);
1193 mYlumi->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt(vals[9]));
1194 myLinFit->SetParameter(0, mYlumi->getTH1()->GetMean(2));
1195 myLinFit->SetParameter(1, 0.0);
1196 mYlumi->getTH1()->Fit(myLinFit,
"QR");
1198 mZlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[2]);
1199 mZlumi->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt(vals[10]));
1200 myLinFit->SetParameter(0, mZlumi->getTH1()->GetMean(2));
1201 myLinFit->SetParameter(1, 0.0);
1202 mZlumi->getTH1()->Fit(myLinFit,
"QR");
1204 sXlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[3]);
1205 sXlumi->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt(vals[11]));
1206 myLinFit->SetParameter(0, sXlumi->getTH1()->GetMean(2));
1207 myLinFit->SetParameter(1, 0.0);
1208 sXlumi->getTH1()->Fit(myLinFit,
"QR");
1210 sYlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[4]);
1211 sYlumi->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt(vals[12]));
1212 myLinFit->SetParameter(0, sYlumi->getTH1()->GetMean(2));
1213 myLinFit->SetParameter(1, 0.0);
1214 sYlumi->getTH1()->Fit(myLinFit,
"QR");
1216 sZlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[5]);
1217 sZlumi->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt(vals[13]));
1218 myLinFit->SetParameter(0, sZlumi->getTH1()->GetMean(2));
1219 myLinFit->SetParameter(1, 0.0);
1220 sZlumi->getTH1()->Fit(myLinFit,
"QR");
1222 dxdzlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[6]);
1223 dxdzlumi->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt(vals[14]));
1224 myLinFit->SetParameter(0, dxdzlumi->getTH1()->GetMean(2));
1225 myLinFit->SetParameter(1, 0.0);
1226 dxdzlumi->getTH1()->Fit(myLinFit,
"QR");
1228 dydzlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[7]);
1229 dydzlumi->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt(vals[15]));
1230 myLinFit->SetParameter(0, dydzlumi->getTH1()->GetMean(2));
1231 myLinFit->SetParameter(1, 0.0);
1232 dydzlumi->getTH1()->Fit(myLinFit,
"QR");
1234 myLinFit->SetParameter(0, hitCounter->getTH1()->GetMean(2));
1235 myLinFit->SetParameter(1, 0.0);
1236 hitCounter->getTH1()->Fit(myLinFit,
"QR");
1238 goodVxCounter->getTH1()->SetBinContent(lastLumiOfFit, (
double)counterVx);
1239 goodVxCounter->getTH1()->SetBinError(
1241 (counterVx != 0 ? 1.
1243 myLinFit->SetParameter(0, goodVxCounter->getTH1()->GetMean(2));
1244 myLinFit->SetParameter(1, 0.0);
1245 goodVxCounter->getTH1()->Fit(myLinFit,
"QR");
1251 TF1* myGaussFit =
new TF1(
"myGaussFit",
1252 "[0]*exp(-(x-[1])*(x-[1])/(2*[2]*[2]))",
1253 Vx_Z_Fit->getTH1()->GetXaxis()->GetXmin(),
1254 Vx_Z_Fit->getTH1()->GetXaxis()->GetXmax());
1255 myGaussFit->SetLineColor(2);
1256 myGaussFit->SetLineWidth(2);
1257 myGaussFit->SetParName(0,
"Ampl.");
1258 myGaussFit->SetParName(1,
"#mu");
1259 myGaussFit->SetParName(2,
"#sigma");
1261 myGaussFit->SetParameter(0, Vx_X_Fit->getTH1()->GetMaximum());
1262 myGaussFit->SetParameter(1, Vx_X_Fit->getTH1()->GetMean());
1263 myGaussFit->SetParameter(2, Vx_X_Fit->getTH1()->GetRMS());
1264 minXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(1);
1265 for (
int i = 0;
i < Vx_X_Fit->getTH1()->GetNbinsX();
i++) {
1266 if (Vx_X_Fit->getTH1()->GetBinContent(
i + 1) > 0) {
1267 minXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(
i + 1);
1271 maxXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(Vx_X_Fit->getTH1()->GetNbinsX());
1272 for (
int i = Vx_X_Fit->getTH1()->GetNbinsX();
i > 0;
i--) {
1273 if (Vx_X_Fit->getTH1()->GetBinContent(
i) > 0) {
1274 maxXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(
i);
1278 myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1279 if (Vx_X_Fit->getTH1()->GetEntries() > 0)
1280 Vx_X_Fit->getTH1()->Fit(myGaussFit,
"QRL");
1282 myGaussFit->SetParameter(0, Vx_Y_Fit->getTH1()->GetMaximum());
1283 myGaussFit->SetParameter(1, Vx_Y_Fit->getTH1()->GetMean());
1284 myGaussFit->SetParameter(2, Vx_Y_Fit->getTH1()->GetRMS());
1285 minXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(1);
1286 for (
int i = 0;
i < Vx_Y_Fit->getTH1()->GetNbinsX();
i++) {
1287 if (Vx_Y_Fit->getTH1()->GetBinContent(
i + 1) > 0) {
1288 minXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(
i + 1);
1292 maxXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(Vx_Y_Fit->getTH1()->GetNbinsX());
1293 for (
int i = Vx_Y_Fit->getTH1()->GetNbinsX();
i > 0;
i--) {
1294 if (Vx_Y_Fit->getTH1()->GetBinContent(
i) > 0) {
1295 maxXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(
i);
1299 myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1300 if (Vx_Y_Fit->getTH1()->GetEntries() > 0)
1301 Vx_Y_Fit->getTH1()->Fit(myGaussFit,
"QRL");
1303 myGaussFit->SetParameter(0, Vx_Z_Fit->getTH1()->GetMaximum());
1304 myGaussFit->SetParameter(1, Vx_Z_Fit->getTH1()->GetMean());
1305 myGaussFit->SetParameter(2, Vx_Z_Fit->getTH1()->GetRMS());
1306 minXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(1);
1307 for (
int i = 0;
i < Vx_Z_Fit->getTH1()->GetNbinsX();
i++) {
1308 if (Vx_Z_Fit->getTH1()->GetBinContent(
i + 1) > 0) {
1309 minXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(
i + 1);
1313 maxXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(Vx_Z_Fit->getTH1()->GetNbinsX());
1314 for (
int i = Vx_Z_Fit->getTH1()->GetNbinsX();
i > 0;
i--) {
1315 if (Vx_Z_Fit->getTH1()->GetBinContent(
i) > 0) {
1316 maxXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(
i);
1320 myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1321 if (Vx_Z_Fit->getTH1()->GetEntries() > 0)
1322 Vx_Z_Fit->getTH1()->Fit(myGaussFit,
"QRL");
1328 fitResults->setAxisTitle(
histTitle.str(), 1);
1329 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true)) {
1330 outputDebugFile <<
"\n" 1332 outputDebugFile <<
"BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) <<
" " << (beginTimeOfFit >> 32) << endl;
1333 outputDebugFile <<
"BeginLumiRange " << beginLumiOfFit << endl;
1334 outputDebugFile <<
histTitle.str().c_str() <<
"\n" << endl;
1337 histTitle <<
"Ongoing: no ongoing fits";
1338 fitResults->setAxisTitle(
histTitle.str(), 1);
1339 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true))
1340 outputDebugFile <<
histTitle.str().c_str() <<
"\n" << endl;
1344 hitCounter->getTH1()->SetBinContent(endLumiOfFit, (
double)totalHits);
1345 hitCounter->getTH1()->SetBinError(endLumiOfFit,
std::sqrt((
double)totalHits));
1350 if (internalDebug ==
true)
1364 Vx_X->setAxisTitle(
"Entries [#]", 2);
1365 Vx_Y->setAxisTitle(
"Primary Vertices Y [cm]", 1);
1366 Vx_Y->setAxisTitle(
"Entries [#]", 2);
1367 Vx_Z->setAxisTitle(
"Primary Vertices Z [cm]", 1);
1368 Vx_Z->setAxisTitle(
"Entries [#]", 2);
1370 Vx_X_Fit = ibooker.
book1D(
"G - vertex x fit",
1371 "Primary Vertex X Distribution (For Fit)",
1375 Vx_Y_Fit = ibooker.
book1D(
"G - vertex y fit",
1376 "Primary Vertex Y Distribution (For Fit)",
1380 Vx_Z_Fit = ibooker.
book1D(
"G - vertex z fit",
1381 "Primary Vertex Z Distribution (For Fit)",
1386 Vx_X_Fit->setAxisTitle(
"Entries [#]", 2);
1387 Vx_Y_Fit->setAxisTitle(
"Primary Vertices Y [cm]", 1);
1388 Vx_Y_Fit->setAxisTitle(
"Entries [#]", 2);
1389 Vx_Z_Fit->setAxisTitle(
"Primary Vertices Z [cm]", 1);
1390 Vx_Z_Fit->setAxisTitle(
"Entries [#]", 2);
1392 Vx_X_Cum = ibooker.
book1D(
"I - vertex x cum",
1393 "Primary Vertex X Distribution (Cumulative)",
1397 Vx_Y_Cum = ibooker.
book1D(
"I - vertex y cum",
1398 "Primary Vertex Y Distribution (Cumulative)",
1402 Vx_Z_Cum = ibooker.
book1D(
"I - vertex z cum",
1403 "Primary Vertex Z Distribution (Cumulative)",
1408 Vx_X_Cum->setAxisTitle(
"Entries [#]", 2);
1409 Vx_Y_Cum->setAxisTitle(
"Primary Vertices Y [cm]", 1);
1410 Vx_Y_Cum->setAxisTitle(
"Entries [#]", 2);
1411 Vx_Z_Cum->setAxisTitle(
"Primary Vertices Z [cm]", 1);
1412 Vx_Z_Cum->setAxisTitle(
"Entries [#]", 2);
1421 mXlumi->setAxisTitle(
"#mu_{x} [cm]", 2);
1422 mXlumi->getTH1()->SetOption(
"E1");
1423 mYlumi->setAxisTitle(
"Lumisection [#]", 1);
1424 mYlumi->setAxisTitle(
"#mu_{y} [cm]", 2);
1425 mYlumi->getTH1()->SetOption(
"E1");
1426 mZlumi->setAxisTitle(
"Lumisection [#]", 1);
1427 mZlumi->setAxisTitle(
"#mu_{z} [cm]", 2);
1428 mZlumi->getTH1()->SetOption(
"E1");
1437 sXlumi->setAxisTitle(
"#sigma_{x} [cm]", 2);
1438 sXlumi->getTH1()->SetOption(
"E1");
1439 sYlumi->setAxisTitle(
"Lumisection [#]", 1);
1440 sYlumi->setAxisTitle(
"#sigma_{y} [cm]", 2);
1441 sYlumi->getTH1()->SetOption(
"E1");
1442 sZlumi->setAxisTitle(
"Lumisection [#]", 1);
1443 sZlumi->setAxisTitle(
"#sigma_{z} [cm]", 2);
1444 sZlumi->getTH1()->SetOption(
"E1");
1446 dxdzlumi = ibooker.
book1D(
1448 dydzlumi = ibooker.
book1D(
1451 dxdzlumi->setAxisTitle(
"dX/dZ [rad]", 2);
1452 dxdzlumi->getTH1()->SetOption(
"E1");
1453 dydzlumi->setAxisTitle(
"Lumisection [#]", 1);
1454 dydzlumi->setAxisTitle(
"dY/dZ [rad]", 2);
1455 dydzlumi->getTH1()->SetOption(
"E1");
1457 Vx_ZX = ibooker.
book2D(
"E - vertex zx",
1458 "Primary Vertex ZX Distribution",
1465 Vx_ZY = ibooker.
book2D(
"E - vertex zy",
1466 "Primary Vertex ZY Distribution",
1473 Vx_XY = ibooker.
book2D(
"E - vertex xy",
1474 "Primary Vertex XY Distribution",
1482 Vx_ZX->setAxisTitle(
"Primary Vertices X [cm]", 2);
1483 Vx_ZX->setAxisTitle(
"Entries [#]", 3);
1484 Vx_ZY->setAxisTitle(
"Primary Vertices Z [cm]", 1);
1485 Vx_ZY->setAxisTitle(
"Primary Vertices Y [cm]", 2);
1486 Vx_ZY->setAxisTitle(
"Entries [#]", 3);
1487 Vx_XY->setAxisTitle(
"Primary Vertices X [cm]", 1);
1488 Vx_XY->setAxisTitle(
"Primary Vertices Y [cm]", 2);
1489 Vx_XY->setAxisTitle(
"Entries [#]", 3);
1491 Vx_ZX_Cum = ibooker.
book2D(
"H - vertex zx cum",
1492 "Primary Vertex ZX Distribution (Cumulative)",
1499 Vx_ZY_Cum = ibooker.
book2D(
"H - vertex zy cum",
1500 "Primary Vertex ZY Distribution (Cumulative)",
1507 Vx_XY_Cum = ibooker.
book2D(
"H - vertex xy cum",
1508 "Primary Vertex XY Distribution (Cumulative)",
1516 Vx_ZX_Cum->setAxisTitle(
"Primary Vertices X [cm]", 2);
1517 Vx_ZX_Cum->setAxisTitle(
"Entries [#]", 3);
1518 Vx_ZY_Cum->setAxisTitle(
"Primary Vertices Z [cm]", 1);
1519 Vx_ZY_Cum->setAxisTitle(
"Primary Vertices Y [cm]", 2);
1520 Vx_ZY_Cum->setAxisTitle(
"Entries [#]", 3);
1521 Vx_XY_Cum->setAxisTitle(
"Primary Vertices X [cm]", 1);
1522 Vx_XY_Cum->setAxisTitle(
"Primary Vertices Y [cm]", 2);
1523 Vx_XY_Cum->setAxisTitle(
"Entries [#]", 3);
1525 hitCounter = ibooker.
book1D(
1528 hitCounter->setAxisTitle(
"Pixel-Hits [#]", 2);
1529 hitCounter->getTH1()->SetOption(
"E1");
1531 goodVxCounter = ibooker.
book1D(
"K - good vertices vs lumi",
1532 "# Good vertices vs. Lumisection",
1537 goodVxCounter->setAxisTitle(
"Good vertices [#]", 2);
1538 goodVxCounter->getTH1()->SetOption(
"E1");
1540 statusCounter = ibooker.
book1D(
1543 statusCounter->getTH1()->SetOption(
"E1");
1544 statusCounter->getTH1()->GetYaxis()->Set(12, -6.5, 5.5);
1545 statusCounter->getTH1()->GetYaxis()->SetBinLabel(1,
"FatalRootError");
1546 statusCounter->getTH1()->GetYaxis()->SetBinLabel(2,
"Max Lumi.");
1547 statusCounter->getTH1()->GetYaxis()->SetBinLabel(3,
"Neg. det.");
1548 statusCounter->getTH1()->GetYaxis()->SetBinLabel(4,
"Infinite err.");
1549 statusCounter->getTH1()->GetYaxis()->SetBinLabel(5,
"No vtx.");
1550 statusCounter->getTH1()->GetYaxis()->SetBinLabel(6,
"Infinite EDM");
1551 statusCounter->getTH1()->GetYaxis()->SetBinLabel(7,
"OK");
1552 statusCounter->getTH1()->GetYaxis()->SetBinLabel(8,
"MINUIT stat.");
1553 statusCounter->getTH1()->GetYaxis()->SetBinLabel(9,
"MINUIT stat.");
1554 statusCounter->getTH1()->GetYaxis()->SetBinLabel(10,
"MINUIT stat.");
1555 statusCounter->getTH1()->GetYaxis()->SetBinLabel(11,
"MINUIT stat.");
1556 statusCounter->getTH1()->GetYaxis()->SetBinLabel(12,
"MINUIT stat.");
1558 fitResults = ibooker.
book2D(
"A - fit results",
"Results of Beam Spot Fit", 2, 0., 2., 9, 0., 9.);
1560 fitResults->setBinLabel(9,
"X[cm]", 2);
1561 fitResults->setBinLabel(8,
"Y[cm]", 2);
1562 fitResults->setBinLabel(7,
"Z[cm]", 2);
1563 fitResults->setBinLabel(6,
"#sigma_{X}[cm]", 2);
1564 fitResults->setBinLabel(5,
"#sigma_{Y}[cm]", 2);
1565 fitResults->setBinLabel(4,
"#sigma_{Z}[cm]", 2);
1566 fitResults->setBinLabel(3,
"#frac{dX}{dZ}[rad]", 2);
1567 fitResults->setBinLabel(2,
"#frac{dY}{dZ}[rad]", 2);
1568 fitResults->setBinLabel(1,
"Vtx[#]", 2);
1569 fitResults->setBinLabel(1,
"Value", 1);
1570 fitResults->setBinLabel(2,
"Error (stat)", 1);
1571 fitResults->getTH1()->SetOption(
"text");
1575 reportSummary = ibooker.
bookFloat(
"reportSummary");
1576 reportSummary->
Fill(-1);
1577 reportSummaryMap = ibooker.
book2D(
"reportSummaryMap",
"Pixel-Vertices Beam Spot: % Good Fits", 1, 0., 1., 1, 0., 1.);
1578 reportSummaryMap->
getTH1()->SetBinContent(1, 1, -1);
T getParameter(std::string const &) const
MonitorElement * bookFloat(TString const &name, FUNC onbooking=NOOP())
virtual void setCurrentFolder(std::string const &fullpath)
constexpr bool isNotFinite(T x)
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
Log< level::Error, false > LogError
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)
const_iterator end(bool update=false) const
T getUntrackedParameter(std::string const &, T const &) const
Vx3DHLTAnalyzer(const edm::ParameterSet &)
std::string formatTime(const time_t &t)
#define DEFINE_FWK_MODULE(type)
void reset(std::string ResetType)
unsigned long long TimeValue_t
Log< level::Info, false > LogInfo
boost::transform_iterator< IterHelp, const_IdIter > const_iterator
const_iterator begin(bool update=false) const
TimeValue_t value() const
Timestamp const & beginTime() const
Timestamp const & endTime() const
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, FUNC onbooking=NOOP())
~Vx3DHLTAnalyzer() override
virtual TH1 * getTH1() const
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())
char const * what() const noexcept override
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
LuminosityBlockNumber_t luminosityBlock() const
void reset(double vett[256])
Power< A, B >::type pow(const A &a, const B &b)
virtual void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)