302 double nSigmaXY = 10.;
304 double parDistanceXY = 1
e-3;
305 double parDistanceZ = 1
e-2;
306 double parDistanceddZ = 1
e-3;
307 double parDistanceCxy = 1
e-5;
310 const unsigned int trials = 4;
311 double largerDist[trials] = {0.1, 5., 10., 100.};
313 double covxz,covyz,det;
315 int bestMovementX = 1;
316 int bestMovementY = 1;
317 int bestMovementZ = 1;
322 vector<double>::const_iterator it =
vals->begin();
324 ROOT::Math::Minimizer* Gauss3D = ROOT::Math::Factory::CreateMinimizer(
"Minuit2",
"Migrad");
325 Gauss3D->SetErrorDef(1.0);
327 else Gauss3D->SetPrintLevel(0);
330 Gauss3D->SetFunction(_Gauss3DFunc);
332 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\t@@@ START FITTING @@@" << endl;
336 for (
int i = 0;
i < 3;
i++)
338 deltaMean = (double(
i)-1.)*
std::sqrt(*(it+0));
339 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl;
343 Gauss3D->SetVariable(0,
"var x ", *(it+0), parDistanceXY * parDistanceXY);
344 Gauss3D->SetVariable(1,
"var y ", *(it+1), parDistanceXY * parDistanceXY);
345 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ * parDistanceZ);
346 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
347 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
348 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
349 Gauss3D->SetVariable(6,
"mean x", *(it+6)+deltaMean, parDistanceXY);
350 Gauss3D->SetVariable(7,
"mean y", *(it+7), parDistanceXY);
351 Gauss3D->SetVariable(8,
"mean z", *(it+8), parDistanceZ);
354 xPos = Gauss3D->X()[6];
355 yPos = Gauss3D->X()[7];
356 zPos = Gauss3D->X()[8];
363 goodData = Gauss3D->Status();
364 edm = Gauss3D->Edm();
367 else if (
isNotFinite(edm) ==
true) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; }
368 else for (
unsigned int j = 0;
j <
nParams;
j++)
372 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
377 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
378 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
380 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
381 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
382 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
383 if (det < 0.) { goodData = -4;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; }
386 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementX =
i; }
388 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tFound bestMovementX --> " << bestMovementX << endl;
392 for (
int i = 0;
i < 3;
i++)
394 deltaMean = (double(
i)-1.)*
std::sqrt(*(it+1));
397 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl;
398 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean X --> " << (double(bestMovementX)-1.)*
std::sqrt(*(it+0)) << endl;
403 Gauss3D->SetVariable(0,
"var x ", *(it+0), parDistanceXY * parDistanceXY);
404 Gauss3D->SetVariable(1,
"var y ", *(it+1), parDistanceXY * parDistanceXY);
405 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ * parDistanceZ);
406 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
407 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
408 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
409 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt(*(it+0)), parDistanceXY);
410 Gauss3D->SetVariable(7,
"mean y", *(it+7)+deltaMean, parDistanceXY);
411 Gauss3D->SetVariable(8,
"mean z", *(it+8), parDistanceZ);
414 xPos = Gauss3D->X()[6];
415 yPos = Gauss3D->X()[7];
416 zPos = Gauss3D->X()[8];
423 goodData = Gauss3D->Status();
424 edm = Gauss3D->Edm();
427 else if (
isNotFinite(edm) ==
true) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; }
428 else for (
unsigned int j = 0;
j <
nParams;
j++)
432 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
437 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
438 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - 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]));
443 if (det < 0.) { goodData = -4;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; }
446 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementY =
i; }
448 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tFound bestMovementY --> " << bestMovementY << endl;
452 for (
int i = 0;
i < 3;
i++)
454 deltaMean = (double(
i)-1.)*
std::sqrt(*(it+2));
457 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl;
458 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean X --> " << (double(bestMovementX)-1.)*
std::sqrt(*(it+0)) << endl;
459 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean Y --> " << (double(bestMovementY)-1.)*
std::sqrt(*(it+1)) << endl;
464 Gauss3D->SetVariable(0,
"var x ", *(it+0), parDistanceXY * parDistanceXY);
465 Gauss3D->SetVariable(1,
"var y ", *(it+1), parDistanceXY * parDistanceXY);
466 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ * parDistanceZ);
467 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
468 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
469 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
470 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt(*(it+0)), parDistanceXY);
471 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt(*(it+1)), parDistanceXY);
472 Gauss3D->SetVariable(8,
"mean z", *(it+8)+deltaMean, parDistanceZ);
475 xPos = Gauss3D->X()[6];
476 yPos = Gauss3D->X()[7];
477 zPos = Gauss3D->X()[8];
484 goodData = Gauss3D->Status();
485 edm = Gauss3D->Edm();
488 else if (
isNotFinite(edm) ==
true) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; }
489 else for (
unsigned int j = 0;
j <
nParams;
j++)
493 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
498 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
499 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
501 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
502 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
503 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
504 if (det < 0.) { goodData = -4;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; }
507 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementZ =
i; }
509 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tFound bestMovementZ --> " << bestMovementZ << endl;
514 Gauss3D->SetVariable(0,
"var x ", *(it+0), parDistanceXY * parDistanceXY);
515 Gauss3D->SetVariable(1,
"var y ", *(it+1), parDistanceXY * parDistanceXY);
516 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ * parDistanceZ);
517 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
518 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
519 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
520 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt(*(it+0)), parDistanceXY);
521 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt(*(it+1)), parDistanceXY);
522 Gauss3D->SetVariable(8,
"mean z", *(it+8)+(
double(bestMovementZ)-1.)*
std::sqrt(*(it+2)), parDistanceZ);
525 xPos = Gauss3D->X()[6];
526 yPos = Gauss3D->X()[7];
527 zPos = Gauss3D->X()[8];
534 goodData = Gauss3D->Status();
535 edm = Gauss3D->Edm();
538 else if (
isNotFinite(edm) ==
true) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; }
539 else for (
unsigned int j = 0;
j <
nParams;
j++)
543 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
548 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
549 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
551 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
552 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
553 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
554 if (det < 0.) { goodData = -4;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; }
558 for (
unsigned int i = 0;
i < trials;
i++)
560 if ((goodData != 0) && (goodData != -2))
564 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tFIT WITH DIFFERENT PARAMETER DISTANCES - STEP " <<
i+1 << endl;
566 Gauss3D->SetVariable(0,
"var x ", *(it+0), parDistanceXY * parDistanceXY * largerDist[
i]);
567 Gauss3D->SetVariable(1,
"var y ", *(it+1), parDistanceXY * parDistanceXY * largerDist[i]);
568 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ * parDistanceZ * largerDist[i]);
569 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy * largerDist[i]);
570 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ * largerDist[i]);
571 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ * largerDist[i]);
572 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt(*(it+0)), parDistanceXY * largerDist[i]);
573 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt(*(it+1)), parDistanceXY * largerDist[i]);
574 Gauss3D->SetVariable(8,
"mean z", *(it+8)+(
double(bestMovementZ)-1.)*
std::sqrt(*(it+2)), parDistanceZ * largerDist[i]);
577 xPos = Gauss3D->X()[6];
578 yPos = Gauss3D->X()[7];
579 zPos = Gauss3D->X()[8];
586 goodData = Gauss3D->Status();
587 edm = Gauss3D->Edm();
590 else if (
isNotFinite(edm) ==
true) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; }
591 else for (
unsigned int j = 0;
j <
nParams;
j++)
595 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
600 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
601 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
603 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
604 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
605 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
606 if (det < 0.) { goodData = -4;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; }
612 for (
unsigned int i = 0; i <
nParams; i++)
614 vals->operator[](
i) = Gauss3D->X()[
i];
double Gauss3DFunc(const double *par)