310 double nSigmaXY = 10.;
312 double parDistanceXY = 1
e-3;
313 double parDistanceZ = 1
e-2;
314 double parDistanceddZ = 1
e-3;
315 double parDistanceCxy = 1
e-5;
318 const unsigned int trials = 4;
319 double largerDist[trials] = {0.1, 5., 10., 100.};
321 double covxz,covyz,det;
323 int bestMovementX = 1;
324 int bestMovementY = 1;
325 int bestMovementZ = 1;
330 vector<double>::const_iterator it =
vals->begin();
332 ROOT::Math::Minimizer* Gauss3D = ROOT::Math::Factory::CreateMinimizer(
"Minuit2",
"Migrad");
333 Gauss3D->SetErrorDef(1.0);
335 else Gauss3D->SetPrintLevel(0);
338 Gauss3D->SetFunction(_Gauss3DFunc);
340 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\t@@@ START FITTING @@@" << endl;
344 for (
int i = 0;
i < 3;
i++)
346 deltaMean = (double(
i)-1.)*
std::sqrt(*(it+0));
347 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl;
351 Gauss3D->SetVariable(0,
"var x ", *(it+0), parDistanceXY * parDistanceXY);
352 Gauss3D->SetVariable(1,
"var y ", *(it+1), parDistanceXY * parDistanceXY);
353 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ * parDistanceZ);
354 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
355 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
356 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
357 Gauss3D->SetVariable(6,
"mean x", *(it+6)+deltaMean, parDistanceXY);
358 Gauss3D->SetVariable(7,
"mean y", *(it+7), parDistanceXY);
359 Gauss3D->SetVariable(8,
"mean z", *(it+8), parDistanceZ);
362 xPos = Gauss3D->X()[6];
363 yPos = Gauss3D->X()[7];
364 zPos = Gauss3D->X()[8];
371 goodData = Gauss3D->Status();
372 edm = Gauss3D->Edm();
375 else if (
isNotFinite(edm) ==
true) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; }
376 else for (
unsigned int j = 0;
j <
nParams;
j++)
380 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
385 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
386 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
388 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
389 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
390 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
391 if (det < 0.) { goodData = -4;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; }
394 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementX =
i; }
396 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tFound bestMovementX --> " << bestMovementX << endl;
400 for (
int i = 0;
i < 3;
i++)
402 deltaMean = (double(
i)-1.)*
std::sqrt(*(it+1));
405 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl;
406 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean X --> " << (double(bestMovementX)-1.)*
std::sqrt(*(it+0)) << endl;
411 Gauss3D->SetVariable(0,
"var x ", *(it+0), parDistanceXY * parDistanceXY);
412 Gauss3D->SetVariable(1,
"var y ", *(it+1), parDistanceXY * parDistanceXY);
413 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ * parDistanceZ);
414 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
415 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
416 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
417 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt(*(it+0)), parDistanceXY);
418 Gauss3D->SetVariable(7,
"mean y", *(it+7)+deltaMean, parDistanceXY);
419 Gauss3D->SetVariable(8,
"mean z", *(it+8), parDistanceZ);
422 xPos = Gauss3D->X()[6];
423 yPos = Gauss3D->X()[7];
424 zPos = Gauss3D->X()[8];
431 goodData = Gauss3D->Status();
432 edm = Gauss3D->Edm();
435 else if (
isNotFinite(edm) ==
true) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; }
436 else for (
unsigned int j = 0;
j <
nParams;
j++)
440 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
445 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
446 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
448 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
449 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
450 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
451 if (det < 0.) { goodData = -4;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; }
454 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementY =
i; }
456 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tFound bestMovementY --> " << bestMovementY << endl;
460 for (
int i = 0;
i < 3;
i++)
462 deltaMean = (double(
i)-1.)*
std::sqrt(*(it+2));
465 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl;
466 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean X --> " << (double(bestMovementX)-1.)*
std::sqrt(*(it+0)) << endl;
467 cout <<
"[Vx3DHLTAnalyzer]::\tdeltaMean Y --> " << (double(bestMovementY)-1.)*
std::sqrt(*(it+1)) << endl;
472 Gauss3D->SetVariable(0,
"var x ", *(it+0), parDistanceXY * parDistanceXY);
473 Gauss3D->SetVariable(1,
"var y ", *(it+1), parDistanceXY * parDistanceXY);
474 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ * parDistanceZ);
475 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
476 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
477 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
478 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt(*(it+0)), parDistanceXY);
479 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt(*(it+1)), parDistanceXY);
480 Gauss3D->SetVariable(8,
"mean z", *(it+8)+deltaMean, parDistanceZ);
483 xPos = Gauss3D->X()[6];
484 yPos = Gauss3D->X()[7];
485 zPos = Gauss3D->X()[8];
492 goodData = Gauss3D->Status();
493 edm = Gauss3D->Edm();
496 else if (
isNotFinite(edm) ==
true) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; }
497 else for (
unsigned int j = 0;
j <
nParams;
j++)
501 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
506 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
507 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
509 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
510 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
511 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
512 if (det < 0.) { goodData = -4;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; }
515 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementZ =
i; }
517 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tFound bestMovementZ --> " << bestMovementZ << endl;
522 Gauss3D->SetVariable(0,
"var x ", *(it+0), parDistanceXY * parDistanceXY);
523 Gauss3D->SetVariable(1,
"var y ", *(it+1), parDistanceXY * parDistanceXY);
524 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ * parDistanceZ);
525 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
526 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
527 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
528 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt(*(it+0)), parDistanceXY);
529 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt(*(it+1)), parDistanceXY);
530 Gauss3D->SetVariable(8,
"mean z", *(it+8)+(
double(bestMovementZ)-1.)*
std::sqrt(*(it+2)), parDistanceZ);
533 xPos = Gauss3D->X()[6];
534 yPos = Gauss3D->X()[7];
535 zPos = Gauss3D->X()[8];
542 goodData = Gauss3D->Status();
543 edm = Gauss3D->Edm();
546 else if (
isNotFinite(edm) ==
true) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; }
547 else for (
unsigned int j = 0;
j <
nParams;
j++)
551 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
556 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
557 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
559 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
560 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
561 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
562 if (det < 0.) { goodData = -4;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; }
566 for (
unsigned int i = 0;
i < trials;
i++)
568 if ((goodData != 0) && (goodData != -2))
572 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tFIT WITH DIFFERENT PARAMETER DISTANCES - STEP " <<
i+1 << endl;
574 Gauss3D->SetVariable(0,
"var x ", *(it+0), parDistanceXY * parDistanceXY * largerDist[
i]);
575 Gauss3D->SetVariable(1,
"var y ", *(it+1), parDistanceXY * parDistanceXY * largerDist[i]);
576 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ * parDistanceZ * largerDist[i]);
577 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy * largerDist[i]);
578 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ * largerDist[i]);
579 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ * largerDist[i]);
580 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt(*(it+0)), parDistanceXY * largerDist[i]);
581 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt(*(it+1)), parDistanceXY * largerDist[i]);
582 Gauss3D->SetVariable(8,
"mean z", *(it+8)+(
double(bestMovementZ)-1.)*
std::sqrt(*(it+2)), parDistanceZ * largerDist[i]);
585 xPos = Gauss3D->X()[6];
586 yPos = Gauss3D->X()[7];
587 zPos = Gauss3D->X()[8];
594 goodData = Gauss3D->Status();
595 edm = Gauss3D->Edm();
598 else if (
isNotFinite(edm) ==
true) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; }
599 else for (
unsigned int j = 0;
j <
nParams;
j++)
603 if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
608 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
609 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
611 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
612 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
613 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
614 if (det < 0.) { goodData = -4;
if (
internalDebug ==
true)
cout <<
"[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; }
620 for (
unsigned int i = 0; i <
nParams; i++)
622 vals->operator[](
i) = Gauss3D->X()[
i];
double Gauss3DFunc(const double *par)