305 double nSigmaXY = 10.;
307 double parDistanceXY = 1
e-3;
308 double parDistanceZ = 1
e-2;
309 double parDistanceddZ = 1
e-3;
310 double parDistanceCxy = 1
e-5;
313 const unsigned int trials = 4;
314 double largerDist[trials] = {0.1, 5., 10., 100.};
316 double covxz,covyz,det;
318 int bestMovementX = 1;
319 int bestMovementY = 1;
320 int bestMovementZ = 1;
325 vector<double>::const_iterator it =
vals->begin();
327 ROOT::Math::Minimizer* Gauss3D = ROOT::Math::Factory::CreateMinimizer(
"Minuit2",
"Migrad");
328 Gauss3D->SetErrorDef(1.0);
330 else Gauss3D->SetPrintLevel(0);
333 Gauss3D->SetFunction(_Gauss3DFunc);
335 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\t@@@ START FITTING @@@" << endl;
339 for (
int i = 0;
i < 3;
i++)
341 deltaMean = (double(
i)-1.)*
std::sqrt(*(it+0));
342 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl;
346 Gauss3D->SetVariable(0,
"var x ", *(it+0), parDistanceXY * parDistanceXY);
347 Gauss3D->SetVariable(1,
"var y ", *(it+1), parDistanceXY * parDistanceXY);
348 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ * parDistanceZ);
349 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
350 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
351 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
352 Gauss3D->SetVariable(6,
"mean x", *(it+6)+deltaMean, parDistanceXY);
353 Gauss3D->SetVariable(7,
"mean y", *(it+7), parDistanceXY);
354 Gauss3D->SetVariable(8,
"mean z", *(it+8), parDistanceZ);
357 xPos = Gauss3D->X()[6];
358 yPos = Gauss3D->X()[7];
359 zPos = Gauss3D->X()[8];
366 goodData = Gauss3D->Status();
367 edm = Gauss3D->Edm();
370 else if (
isNotFinite(edm) ==
true) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; }
371 else for (
unsigned int j = 0; j <
nParams; j++)
375 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
380 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
381 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
383 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
384 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
385 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
386 if (det < 0.) { goodData = -4;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; }
389 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementX =
i; }
391 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tFound bestMovementX --> " << bestMovementX << endl;
395 for (
int i = 0;
i < 3;
i++)
397 deltaMean = (double(
i)-1.)*
std::sqrt(*(it+1));
400 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl;
401 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean X --> " << (double(bestMovementX)-1.)*
std::sqrt(*(it+0)) << endl;
406 Gauss3D->SetVariable(0,
"var x ", *(it+0), parDistanceXY * parDistanceXY);
407 Gauss3D->SetVariable(1,
"var y ", *(it+1), parDistanceXY * parDistanceXY);
408 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ * parDistanceZ);
409 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
410 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
411 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
412 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt(*(it+0)), parDistanceXY);
413 Gauss3D->SetVariable(7,
"mean y", *(it+7)+deltaMean, parDistanceXY);
414 Gauss3D->SetVariable(8,
"mean z", *(it+8), parDistanceZ);
417 xPos = Gauss3D->X()[6];
418 yPos = Gauss3D->X()[7];
419 zPos = Gauss3D->X()[8];
426 goodData = Gauss3D->Status();
427 edm = Gauss3D->Edm();
430 else if (
isNotFinite(edm) ==
true) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; }
431 else for (
unsigned int j = 0; j <
nParams; j++)
435 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
440 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
441 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
443 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
444 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
445 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
446 if (det < 0.) { goodData = -4;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; }
449 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementY =
i; }
451 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tFound bestMovementY --> " << bestMovementY << endl;
455 for (
int i = 0;
i < 3;
i++)
457 deltaMean = (double(
i)-1.)*
std::sqrt(*(it+2));
460 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl;
461 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean X --> " << (double(bestMovementX)-1.)*
std::sqrt(*(it+0)) << endl;
462 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean Y --> " << (double(bestMovementY)-1.)*
std::sqrt(*(it+1)) << endl;
467 Gauss3D->SetVariable(0,
"var x ", *(it+0), parDistanceXY * parDistanceXY);
468 Gauss3D->SetVariable(1,
"var y ", *(it+1), parDistanceXY * parDistanceXY);
469 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ * parDistanceZ);
470 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
471 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
472 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
473 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt(*(it+0)), parDistanceXY);
474 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt(*(it+1)), parDistanceXY);
475 Gauss3D->SetVariable(8,
"mean z", *(it+8)+deltaMean, parDistanceZ);
478 xPos = Gauss3D->X()[6];
479 yPos = Gauss3D->X()[7];
480 zPos = Gauss3D->X()[8];
487 goodData = Gauss3D->Status();
488 edm = Gauss3D->Edm();
491 else if (
isNotFinite(edm) ==
true) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; }
492 else for (
unsigned int j = 0; j <
nParams; j++)
496 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
501 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
502 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
504 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
505 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
506 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
507 if (det < 0.) { goodData = -4;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; }
510 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementZ =
i; }
512 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tFound bestMovementZ --> " << bestMovementZ << endl;
517 Gauss3D->SetVariable(0,
"var x ", *(it+0), parDistanceXY * parDistanceXY);
518 Gauss3D->SetVariable(1,
"var y ", *(it+1), parDistanceXY * parDistanceXY);
519 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ * parDistanceZ);
520 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
521 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
522 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
523 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt(*(it+0)), parDistanceXY);
524 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt(*(it+1)), parDistanceXY);
525 Gauss3D->SetVariable(8,
"mean z", *(it+8)+(
double(bestMovementZ)-1.)*
std::sqrt(*(it+2)), parDistanceZ);
528 xPos = Gauss3D->X()[6];
529 yPos = Gauss3D->X()[7];
530 zPos = Gauss3D->X()[8];
537 goodData = Gauss3D->Status();
538 edm = Gauss3D->Edm();
541 else if (
isNotFinite(edm) ==
true) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; }
542 else for (
unsigned int j = 0; j <
nParams; j++)
546 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
551 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
552 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
554 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
555 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
556 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
557 if (det < 0.) { goodData = -4;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; }
561 for (
unsigned int i = 0;
i < trials;
i++)
563 if ((goodData != 0) && (goodData != -2))
567 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tFIT WITH DIFFERENT PARAMETER DISTANCES - STEP " <<
i+1 << endl;
569 Gauss3D->SetVariable(0,
"var x ", *(it+0), parDistanceXY * parDistanceXY * largerDist[
i]);
570 Gauss3D->SetVariable(1,
"var y ", *(it+1), parDistanceXY * parDistanceXY * largerDist[i]);
571 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ * parDistanceZ * largerDist[i]);
572 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy * largerDist[i]);
573 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ * largerDist[i]);
574 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ * largerDist[i]);
575 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt(*(it+0)), parDistanceXY * largerDist[i]);
576 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt(*(it+1)), parDistanceXY * largerDist[i]);
577 Gauss3D->SetVariable(8,
"mean z", *(it+8)+(
double(bestMovementZ)-1.)*
std::sqrt(*(it+2)), parDistanceZ * largerDist[i]);
580 xPos = Gauss3D->X()[6];
581 yPos = Gauss3D->X()[7];
582 zPos = Gauss3D->X()[8];
589 goodData = Gauss3D->Status();
590 edm = Gauss3D->Edm();
593 else if (
isNotFinite(edm) ==
true) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; }
594 else for (
unsigned int j = 0; j <
nParams; j++)
598 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
603 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
604 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
606 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
607 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
608 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
609 if (det < 0.) { goodData = -4;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; }
615 for (
unsigned int i = 0; i <
nParams; i++)
617 vals->operator[](
i) = Gauss3D->X()[
i];
double Gauss3DFunc(const double *par)