251 unsigned int nParams = 9;
255 double nSigmaXY = 100.;
256 double nSigmaZ = 100.;
257 double varFactor = 4./25.;
258 double parDistanceXY = 0.005;
259 double parDistanceZ = 0.5;
260 double parDistanceddZ = 1
e-3;
261 double parDistanceCxy = 1
e-5;
262 double bestEdm = 1
e-1;
264 const unsigned int trials = 4;
265 double largerDist[trials] = {0.1, 5., 10., 100.};
267 double covxz,covyz,det;
269 int bestMovementX = 1;
270 int bestMovementY = 1;
271 int bestMovementZ = 1;
276 vector<double>::const_iterator it =
vals->begin();
278 ROOT::Math::Minimizer* Gauss3D = ROOT::Math::Factory::CreateMinimizer(
"Minuit2",
"Migrad");
279 Gauss3D->SetMaxFunctionCalls(1e4);
280 Gauss3D->SetTolerance(1
e-9);
282 else Gauss3D->SetPrintLevel(0);
284 ROOT::Math::Functor _Gauss3DFunc(&
Gauss3DFunc,nParams);
285 Gauss3D->SetFunction(_Gauss3DFunc);
291 for (
int i = 0;
i < 3;
i++)
293 deltaMean = (double(
i)-1.)*
std::sqrt((*(it+0))*varFactor);
298 Gauss3D->SetVariable(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY);
299 Gauss3D->SetVariable(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY);
300 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ);
301 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
302 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
303 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
304 Gauss3D->SetVariable(6,
"mean x", *(it+6)+deltaMean, parDistanceXY);
305 Gauss3D->SetVariable(7,
"mean y", *(it+7), parDistanceXY);
306 Gauss3D->SetVariable(8,
"mean z", *(it+8), parDistanceZ);
309 xPos = Gauss3D->X()[6];
310 yPos = Gauss3D->X()[7];
311 zPos = Gauss3D->X()[8];
318 goodData = Gauss3D->Status();
319 edm = Gauss3D->Edm();
323 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
edm::isNotFinite(Gauss3D->Errors()[
j]) ==
true) { goodData = -1;
break; }
326 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
327 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
329 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
330 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
331 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
332 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
335 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementX =
i; }
337 if (
internalDebug ==
true)
cout <<
"Found bestMovementX --> " << bestMovementX << endl;
341 for (
int i = 0;
i < 3;
i++)
343 deltaMean = (double(
i)-1.)*
std::sqrt((*(it+1))*varFactor);
346 cout <<
"deltaMean --> " << deltaMean << endl;
347 cout <<
"deltaMean X --> " << (double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor) << endl;
352 Gauss3D->SetVariable(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY);
353 Gauss3D->SetVariable(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY);
354 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ);
355 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
356 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
357 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
358 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY);
359 Gauss3D->SetVariable(7,
"mean y", *(it+7)+deltaMean, parDistanceXY);
360 Gauss3D->SetVariable(8,
"mean z", *(it+8), parDistanceZ);
363 xPos = Gauss3D->X()[6];
364 yPos = Gauss3D->X()[7];
365 zPos = Gauss3D->X()[8];
372 goodData = Gauss3D->Status();
373 edm = Gauss3D->Edm();
377 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
edm::isNotFinite(Gauss3D->Errors()[
j]) ==
true) { goodData = -1;
break; }
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 = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
389 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementY =
i; }
391 if (
internalDebug ==
true)
cout <<
"Found bestMovementY --> " << bestMovementY << endl;
395 for (
int i = 0;
i < 3;
i++)
397 deltaMean = (double(
i)-1.)*
std::sqrt(*(it+2));
400 cout <<
"deltaMean --> " << deltaMean << endl;
401 cout <<
"deltaMean X --> " << (double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor) << endl;
402 cout <<
"deltaMean Y --> " << (double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor) << endl;
407 Gauss3D->SetVariable(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY);
408 Gauss3D->SetVariable(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY);
409 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ);
410 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
411 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
412 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
413 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY);
414 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor), parDistanceXY);
415 Gauss3D->SetVariable(8,
"mean z", *(it+8)+deltaMean, parDistanceZ);
418 xPos = Gauss3D->X()[6];
419 yPos = Gauss3D->X()[7];
420 zPos = Gauss3D->X()[8];
427 goodData = Gauss3D->Status();
428 edm = Gauss3D->Edm();
432 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
edm::isNotFinite(Gauss3D->Errors()[
j]) ==
true) { goodData = -1;
break; }
435 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
436 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
438 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
439 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
440 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
441 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
444 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementZ =
i; }
446 if (
internalDebug ==
true)
cout <<
"Found bestMovementZ --> " << bestMovementZ << endl;
451 Gauss3D->SetVariable(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY);
452 Gauss3D->SetVariable(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY);
453 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ);
454 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
455 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
456 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
457 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY);
458 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor), parDistanceXY);
459 Gauss3D->SetVariable(8,
"mean z", *(it+8)+(
double(bestMovementZ)-1.)*
std::sqrt(*(it+2)), parDistanceZ);
462 xPos = Gauss3D->X()[6];
463 yPos = Gauss3D->X()[7];
464 zPos = Gauss3D->X()[8];
471 goodData = Gauss3D->Status();
472 edm = Gauss3D->Edm();
476 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
edm::isNotFinite(Gauss3D->Errors()[
j]) ==
true) { goodData = -1;
break; }
479 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
480 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
482 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
483 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
484 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
485 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
489 for (
unsigned int i = 0;
i < trials;
i++)
491 if ((goodData != 0) && (goodData != -2))
495 if (
internalDebug ==
true)
cout <<
"FIT WITH DIFFERENT PARAMETER DISTANCES - STEP " <<
i+1 << endl;
497 Gauss3D->SetVariable(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY * largerDist[
i]);
498 Gauss3D->SetVariable(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY * largerDist[i]);
499 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ * largerDist[i]);
500 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy * largerDist[i]);
501 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ * largerDist[i]);
502 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ * largerDist[i]);
503 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY * largerDist[i]);
504 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor), parDistanceXY * largerDist[i]);
505 Gauss3D->SetVariable(8,
"mean z", *(it+8)+(
double(bestMovementZ)-1.)*
std::sqrt(*(it+2)), parDistanceZ * largerDist[i]);
508 xPos = Gauss3D->X()[6];
509 yPos = Gauss3D->X()[7];
510 zPos = Gauss3D->X()[8];
517 goodData = Gauss3D->Status();
518 edm = Gauss3D->Edm();
522 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
edm::isNotFinite(Gauss3D->Errors()[
j]) ==
true) { goodData = -1;
break; }
525 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
526 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
528 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
529 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
530 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
531 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
537 for (
unsigned int i = 0; i < nParams; i++)
539 vals->operator[](
i) = Gauss3D->X()[
i];
540 vals->operator[](i+nParams) = Gauss3D->Errors()[
i];
double Gauss3DFunc(const double *par)