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);
306 Gauss3D->SetPrintLevel(3);
308 Gauss3D->SetPrintLevel(0);
311 Gauss3D->SetFunction(_Gauss3DFunc);
314 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\t@@@ START FITTING @@@";
318 for (
int i = 0;
i < 3;
i++) {
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;
346 goodData = Gauss3D->Status();
347 edm = Gauss3D->Edm();
349 edm::LogError(
"Vx3DHLTAnalyzer") <<
"\tCaught Minuit2 exception @ Fit at X: \n" << er.
what();
360 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNot finite edm !";
361 }
else if (goodData != -6)
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]));
381 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNegative determinant !";
385 if ((goodData == 0) && (std::fabs(
edm) < bestEdm)) {
391 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tFound bestMovementX --> " << bestMovementX;
395 for (
int i = 0;
i < 3;
i++) {
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);
416 yPos = *(
it + 7) + deltaMean;
425 goodData = Gauss3D->Status();
426 edm = Gauss3D->Edm();
428 edm::LogError(
"Vx3DHLTAnalyzer") <<
"\tCaught Minuit2 exception @ Fit at Y: \n" << er.
what();
439 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNot finite edm !";
440 }
else if (goodData != -6)
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]));
460 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNegative determinant !";
464 if ((goodData == 0) && (std::fabs(
edm) < bestEdm)) {
470 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tFound bestMovementY --> " << bestMovementY;
474 for (
int i = 0;
i < 3;
i++) {
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);
497 zPos = *(
it + 8) + deltaMean;
505 goodData = Gauss3D->Status();
506 edm = Gauss3D->Edm();
508 edm::LogError(
"Vx3DHLTAnalyzer") <<
"\tCaught Minuit2 exception @ Fit at Z: \n" << er.
what();
519 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNot finite edm !";
520 }
else if (goodData != -6)
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]));
540 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNegative determinant !";
544 if ((goodData == 0) && (std::fabs(
edm) < bestEdm)) {
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);
576 goodData = Gauss3D->Status();
577 edm = Gauss3D->Edm();
579 edm::LogError(
"Vx3DHLTAnalyzer") <<
"\tCaught Minuit2 exception @ Final fit: \n" << er.
what();
589 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNot finite edm !";
590 }
else if (goodData != -6)
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]));
610 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNegative determinant !";
615 for (
unsigned int i = 0;
i < trials;
i++) {
616 if ((goodData != 0) && (goodData != -2)) {
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]);
650 goodData = Gauss3D->Status();
651 edm = Gauss3D->Edm();
653 edm::LogError(
"Vx3DHLTAnalyzer") <<
"\tCaught Minuit2 exception @ Fit with different distances: \n" 665 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNot finite edm !";
666 }
else if (goodData != -6)
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]));
686 edm::LogInfo(
"Vx3DHLTAnalyzer") <<
"\tNegative determinant !";
695 vals->operator[](
i) = Gauss3D->X()[
i];
696 vals->operator[](
i +
nParams) = Gauss3D->Errors()[
i];
constexpr bool isNotFinite(T x)
Log< level::Error, false > LogError
double Gauss3DFunc(const double *par)
Log< level::Info, false > LogInfo
char const * what() const noexcept override