252 unsigned int nParams = 9;
256 double nSigmaXY = 100.;
257 double nSigmaZ = 100.;
258 double varFactor = 4./25.;
259 double parDistanceXY = 0.005;
260 double parDistanceZ = 0.5;
261 double parDistanceddZ = 1
e-3;
262 double parDistanceCxy = 1
e-5;
263 double bestEdm = 1
e-1;
265 const unsigned int trials = 4;
266 double largerDist[trials] = {0.1, 5., 10., 100.};
268 double covxz,covyz,det;
270 int bestMovementX = 1;
271 int bestMovementY = 1;
272 int bestMovementZ = 1;
276 double amin,errdef,edm;
279 vector<double>::const_iterator it =
vals->begin();
281 TFitterMinuit* Gauss3D =
new TFitterMinuit(nParams);
283 else Gauss3D->SetPrintLevel(0);
292 for (
int i = 0;
i < 3;
i++)
294 deltaMean = (double(
i)-1.)*
std::sqrt((*(it+0))*varFactor);
301 Gauss3D->SetParameter(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
302 Gauss3D->SetParameter(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
303 Gauss3D->SetParameter(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ, 0., 0.);
304 Gauss3D->SetParameter(3,
"cov xy", *(it+3), parDistanceCxy, 0., 0.);
305 Gauss3D->SetParameter(4,
"dydz ", *(it+4), parDistanceddZ, 0., 0.);
306 Gauss3D->SetParameter(5,
"dxdz ", *(it+5), parDistanceddZ, 0., 0.);
307 Gauss3D->SetParameter(6,
"mean x", *(it+6)+deltaMean, parDistanceXY, 0., 0.);
308 Gauss3D->SetParameter(7,
"mean y", *(it+7), parDistanceXY, 0., 0.);
309 Gauss3D->SetParameter(8,
"mean z", *(it+8), parDistanceZ, 0., 0.);
312 xPos = Gauss3D->GetParameter(6);
313 yPos = Gauss3D->GetParameter(7);
314 zPos = Gauss3D->GetParameter(8);
317 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(Gauss3D->GetParameter(0)) + std::fabs(Gauss3D->GetParameter(1))) / 2.;
320 goodData = Gauss3D->ExecuteCommand(
"MIGRAD",arglist,2);
321 Gauss3D->GetStats(amin, edm, errdef, nvpar, nparx);
324 else if (
std::isnan(edm) ==
true) goodData = -1;
325 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
std::isnan(Gauss3D->GetParError(
j)) ==
true) { goodData = -1;
break; }
328 covyz = Gauss3D->GetParameter(4)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(1))) - Gauss3D->GetParameter(5)*Gauss3D->GetParameter(3);
329 covxz = Gauss3D->GetParameter(5)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(0))) - Gauss3D->GetParameter(4)*Gauss3D->GetParameter(3);
331 det = std::fabs(Gauss3D->GetParameter(0)) * (std::fabs(Gauss3D->GetParameter(1))*std::fabs(Gauss3D->GetParameter(2)) - covyz*covyz) -
332 Gauss3D->GetParameter(3) * (Gauss3D->GetParameter(3)*std::fabs(Gauss3D->GetParameter(2)) - covxz*covyz) +
333 covxz * (Gauss3D->GetParameter(3)*covyz - covxz*std::fabs(Gauss3D->GetParameter(1)));
334 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
337 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementX =
i; }
339 if (
internalDebug ==
true)
cout <<
"Found bestMovementX --> " << bestMovementX << endl;
343 for (
int i = 0;
i < 3;
i++)
345 deltaMean = (double(
i)-1.)*
std::sqrt((*(it+1))*varFactor);
348 cout <<
"deltaMean --> " << deltaMean << endl;
349 cout <<
"deltaMean X --> " << (double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor) << endl;
356 Gauss3D->SetParameter(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
357 Gauss3D->SetParameter(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
358 Gauss3D->SetParameter(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ, 0., 0.);
359 Gauss3D->SetParameter(3,
"cov xy", *(it+3), parDistanceCxy, 0., 0.);
360 Gauss3D->SetParameter(4,
"dydz ", *(it+4), parDistanceddZ, 0., 0.);
361 Gauss3D->SetParameter(5,
"dxdz ", *(it+5), parDistanceddZ, 0., 0.);
362 Gauss3D->SetParameter(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY, 0., 0.);
363 Gauss3D->SetParameter(7,
"mean y", *(it+7)+deltaMean, parDistanceXY, 0., 0.);
364 Gauss3D->SetParameter(8,
"mean z", *(it+8), parDistanceZ, 0., 0.);
367 xPos = Gauss3D->GetParameter(6);
368 yPos = Gauss3D->GetParameter(7);
369 zPos = Gauss3D->GetParameter(8);
372 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(Gauss3D->GetParameter(0)) + std::fabs(Gauss3D->GetParameter(1))) / 2.;
375 goodData = Gauss3D->ExecuteCommand(
"MIGRAD",arglist,2);
376 Gauss3D->GetStats(amin, edm, errdef, nvpar, nparx);
379 else if (
std::isnan(edm) ==
true) goodData = -1;
380 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
std::isnan(Gauss3D->GetParError(
j)) ==
true) { goodData = -1;
break; }
383 covyz = Gauss3D->GetParameter(4)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(1))) - Gauss3D->GetParameter(5)*Gauss3D->GetParameter(3);
384 covxz = Gauss3D->GetParameter(5)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(0))) - Gauss3D->GetParameter(4)*Gauss3D->GetParameter(3);
386 det = std::fabs(Gauss3D->GetParameter(0)) * (std::fabs(Gauss3D->GetParameter(1))*std::fabs(Gauss3D->GetParameter(2)) - covyz*covyz) -
387 Gauss3D->GetParameter(3) * (Gauss3D->GetParameter(3)*std::fabs(Gauss3D->GetParameter(2)) - covxz*covyz) +
388 covxz * (Gauss3D->GetParameter(3)*covyz - covxz*std::fabs(Gauss3D->GetParameter(1)));
389 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
392 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementY =
i; }
394 if (
internalDebug ==
true)
cout <<
"Found bestMovementY --> " << bestMovementY << endl;
398 for (
int i = 0;
i < 3;
i++)
400 deltaMean = (double(
i)-1.)*
std::sqrt(*(it+2));
403 cout <<
"deltaMean --> " << deltaMean << endl;
404 cout <<
"deltaMean X --> " << (double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor) << endl;
405 cout <<
"deltaMean Y --> " << (double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor) << endl;
412 Gauss3D->SetParameter(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
413 Gauss3D->SetParameter(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
414 Gauss3D->SetParameter(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ, 0., 0.);
415 Gauss3D->SetParameter(3,
"cov xy", *(it+3), parDistanceCxy, 0., 0.);
416 Gauss3D->SetParameter(4,
"dydz ", *(it+4), parDistanceddZ, 0., 0.);
417 Gauss3D->SetParameter(5,
"dxdz ", *(it+5), parDistanceddZ, 0., 0.);
418 Gauss3D->SetParameter(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY, 0., 0.);
419 Gauss3D->SetParameter(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor), parDistanceXY, 0., 0.);
420 Gauss3D->SetParameter(8,
"mean z", *(it+8)+deltaMean, parDistanceZ, 0., 0.);
423 xPos = Gauss3D->GetParameter(6);
424 yPos = Gauss3D->GetParameter(7);
425 zPos = Gauss3D->GetParameter(8);
428 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(Gauss3D->GetParameter(0)) + std::fabs(Gauss3D->GetParameter(1))) / 2.;
431 goodData = Gauss3D->ExecuteCommand(
"MIGRAD",arglist,2);
432 Gauss3D->GetStats(amin, edm, errdef, nvpar, nparx);
435 else if (
std::isnan(edm) ==
true) goodData = -1;
436 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
std::isnan(Gauss3D->GetParError(
j)) ==
true) { goodData = -1;
break; }
439 covyz = Gauss3D->GetParameter(4)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(1))) - Gauss3D->GetParameter(5)*Gauss3D->GetParameter(3);
440 covxz = Gauss3D->GetParameter(5)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(0))) - Gauss3D->GetParameter(4)*Gauss3D->GetParameter(3);
442 det = std::fabs(Gauss3D->GetParameter(0)) * (std::fabs(Gauss3D->GetParameter(1))*std::fabs(Gauss3D->GetParameter(2)) - covyz*covyz) -
443 Gauss3D->GetParameter(3) * (Gauss3D->GetParameter(3)*std::fabs(Gauss3D->GetParameter(2)) - covxz*covyz) +
444 covxz * (Gauss3D->GetParameter(3)*covyz - covxz*std::fabs(Gauss3D->GetParameter(1)));
445 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
448 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementZ =
i; }
450 if (
internalDebug ==
true)
cout <<
"Found bestMovementZ --> " << bestMovementZ << endl;
457 Gauss3D->SetParameter(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
458 Gauss3D->SetParameter(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
459 Gauss3D->SetParameter(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ, 0., 0.);
460 Gauss3D->SetParameter(3,
"cov xy", *(it+3), parDistanceCxy, 0., 0.);
461 Gauss3D->SetParameter(4,
"dydz ", *(it+4), parDistanceddZ, 0., 0.);
462 Gauss3D->SetParameter(5,
"dxdz ", *(it+5), parDistanceddZ, 0., 0.);
463 Gauss3D->SetParameter(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY, 0., 0.);
464 Gauss3D->SetParameter(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor), parDistanceXY, 0., 0.);
465 Gauss3D->SetParameter(8,
"mean z", *(it+8)+(
double(bestMovementZ)-1.)*
std::sqrt(*(it+2)), parDistanceZ, 0., 0.);
468 xPos = Gauss3D->GetParameter(6);
469 yPos = Gauss3D->GetParameter(7);
470 zPos = Gauss3D->GetParameter(8);
473 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(Gauss3D->GetParameter(0)) + std::fabs(Gauss3D->GetParameter(1))) / 2.;
476 goodData = Gauss3D->ExecuteCommand(
"MIGRAD",arglist,2);
477 Gauss3D->GetStats(amin, edm, errdef, nvpar, nparx);
480 else if (
std::isnan(edm) ==
true) goodData = -1;
481 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
std::isnan(Gauss3D->GetParError(
j)) ==
true) { goodData = -1;
break; }
484 covyz = Gauss3D->GetParameter(4)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(1))) - Gauss3D->GetParameter(5)*Gauss3D->GetParameter(3);
485 covxz = Gauss3D->GetParameter(5)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(0))) - Gauss3D->GetParameter(4)*Gauss3D->GetParameter(3);
487 det = std::fabs(Gauss3D->GetParameter(0)) * (std::fabs(Gauss3D->GetParameter(1))*std::fabs(Gauss3D->GetParameter(2)) - covyz*covyz) -
488 Gauss3D->GetParameter(3) * (Gauss3D->GetParameter(3)*std::fabs(Gauss3D->GetParameter(2)) - covxz*covyz) +
489 covxz * (Gauss3D->GetParameter(3)*covyz - covxz*std::fabs(Gauss3D->GetParameter(1)));
490 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
496 for (
unsigned int i = 0;
i < trials;
i++)
498 if ((goodData != 0) && (goodData != -2))
502 if (
internalDebug ==
true)
cout <<
"FIT WITH DIFFERENT PARAMETER DISTANCES - STEP " <<
i+1 << endl;
504 Gauss3D->SetParameter(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY * largerDist[
i], 0, 0);
505 Gauss3D->SetParameter(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY * largerDist[i], 0, 0);
506 Gauss3D->SetParameter(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ * largerDist[i], 0, 0);
507 Gauss3D->SetParameter(3,
"cov xy", *(it+3), parDistanceCxy * largerDist[i], 0, 0);
508 Gauss3D->SetParameter(4,
"dydz ", *(it+4), parDistanceddZ * largerDist[i], 0, 0);
509 Gauss3D->SetParameter(5,
"dxdz ", *(it+5), parDistanceddZ * largerDist[i], 0, 0);
510 Gauss3D->SetParameter(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY * largerDist[i], 0, 0);
511 Gauss3D->SetParameter(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor), parDistanceXY * largerDist[i], 0, 0);
512 Gauss3D->SetParameter(8,
"mean z", *(it+8)+(
double(bestMovementZ)-1.)*
std::sqrt(*(it+2)), parDistanceZ * largerDist[i], 0, 0);
515 xPos = Gauss3D->GetParameter(6);
516 yPos = Gauss3D->GetParameter(7);
517 zPos = Gauss3D->GetParameter(8);
520 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(Gauss3D->GetParameter(0)) + std::fabs(Gauss3D->GetParameter(1))) / 2.;
523 goodData = Gauss3D->ExecuteCommand(
"MIGRAD",arglist,2);
524 Gauss3D->GetStats(amin, edm, errdef, nvpar, nparx);
527 else if (
std::isnan(edm) ==
true) goodData = -1;
528 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
std::isnan(Gauss3D->GetParError(
j)) ==
true) { goodData = -1;
break; }
531 covyz = Gauss3D->GetParameter(4)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(1))) - Gauss3D->GetParameter(5)*Gauss3D->GetParameter(3);
532 covxz = Gauss3D->GetParameter(5)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(0))) - Gauss3D->GetParameter(4)*Gauss3D->GetParameter(3);
534 det = std::fabs(Gauss3D->GetParameter(0)) * (std::fabs(Gauss3D->GetParameter(1))*std::fabs(Gauss3D->GetParameter(2)) - covyz*covyz) -
535 Gauss3D->GetParameter(3) * (Gauss3D->GetParameter(3)*std::fabs(Gauss3D->GetParameter(2)) - covxz*covyz) +
536 covxz * (Gauss3D->GetParameter(3)*covyz - covxz*std::fabs(Gauss3D->GetParameter(1)));
537 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
543 for (
unsigned int i = 0; i < nParams; i++)
545 vals->operator[](
i) = Gauss3D->GetParameter(i);
546 vals->operator[](i+nParams) = Gauss3D->GetParError(i);
void Gauss3DFunc(int &, double *, double &fval, double *par, int)