253 unsigned int nParams = 9;
257 double nSigmaXY = 100.;
258 double nSigmaZ = 100.;
259 double varFactor = 4./25.;
260 double parDistanceXY = 0.005;
261 double parDistanceZ = 0.5;
262 double parDistanceddZ = 1
e-3;
263 double parDistanceCxy = 1
e-5;
264 double bestEdm = 1
e-1;
266 const unsigned int trials = 4;
267 double largerDist[trials] = {0.1, 5., 10., 100.};
269 double covxz,covyz,det;
271 int bestMovementX = 1;
272 int bestMovementY = 1;
273 int bestMovementZ = 1;
277 double amin,errdef,edm;
280 vector<double>::const_iterator it =
vals->begin();
282 TFitterMinuit* Gauss3D =
new TFitterMinuit(nParams);
284 else Gauss3D->SetPrintLevel(0);
293 for (
int i = 0;
i < 3;
i++)
295 deltaMean = (double(
i)-1.)*
std::sqrt((*(it+0))*varFactor);
302 Gauss3D->SetParameter(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
303 Gauss3D->SetParameter(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
304 Gauss3D->SetParameter(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ, 0., 0.);
305 Gauss3D->SetParameter(3,
"cov xy", *(it+3), parDistanceCxy, 0., 0.);
306 Gauss3D->SetParameter(4,
"dydz ", *(it+4), parDistanceddZ, 0., 0.);
307 Gauss3D->SetParameter(5,
"dxdz ", *(it+5), parDistanceddZ, 0., 0.);
308 Gauss3D->SetParameter(6,
"mean x", *(it+6)+deltaMean, parDistanceXY, 0., 0.);
309 Gauss3D->SetParameter(7,
"mean y", *(it+7), parDistanceXY, 0., 0.);
310 Gauss3D->SetParameter(8,
"mean z", *(it+8), parDistanceZ, 0., 0.);
313 xPos = Gauss3D->GetParameter(6);
314 yPos = Gauss3D->GetParameter(7);
315 zPos = Gauss3D->GetParameter(8);
318 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(Gauss3D->GetParameter(0)) + std::fabs(Gauss3D->GetParameter(1))) / 2.;
321 goodData = Gauss3D->ExecuteCommand(
"MIGRAD",arglist,2);
322 Gauss3D->GetStats(amin, edm, errdef, nvpar, nparx);
326 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
edm::isNotFinite(Gauss3D->GetParError(
j)) ==
true) { goodData = -1;
break; }
329 covyz = Gauss3D->GetParameter(4)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(1))) - Gauss3D->GetParameter(5)*Gauss3D->GetParameter(3);
330 covxz = Gauss3D->GetParameter(5)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(0))) - Gauss3D->GetParameter(4)*Gauss3D->GetParameter(3);
332 det = std::fabs(Gauss3D->GetParameter(0)) * (std::fabs(Gauss3D->GetParameter(1))*std::fabs(Gauss3D->GetParameter(2)) - covyz*covyz) -
333 Gauss3D->GetParameter(3) * (Gauss3D->GetParameter(3)*std::fabs(Gauss3D->GetParameter(2)) - covxz*covyz) +
334 covxz * (Gauss3D->GetParameter(3)*covyz - covxz*std::fabs(Gauss3D->GetParameter(1)));
335 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
338 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementX =
i; }
340 if (
internalDebug ==
true)
cout <<
"Found bestMovementX --> " << bestMovementX << endl;
344 for (
int i = 0;
i < 3;
i++)
346 deltaMean = (double(
i)-1.)*
std::sqrt((*(it+1))*varFactor);
349 cout <<
"deltaMean --> " << deltaMean << endl;
350 cout <<
"deltaMean X --> " << (double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor) << endl;
357 Gauss3D->SetParameter(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
358 Gauss3D->SetParameter(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
359 Gauss3D->SetParameter(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ, 0., 0.);
360 Gauss3D->SetParameter(3,
"cov xy", *(it+3), parDistanceCxy, 0., 0.);
361 Gauss3D->SetParameter(4,
"dydz ", *(it+4), parDistanceddZ, 0., 0.);
362 Gauss3D->SetParameter(5,
"dxdz ", *(it+5), parDistanceddZ, 0., 0.);
363 Gauss3D->SetParameter(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY, 0., 0.);
364 Gauss3D->SetParameter(7,
"mean y", *(it+7)+deltaMean, parDistanceXY, 0., 0.);
365 Gauss3D->SetParameter(8,
"mean z", *(it+8), parDistanceZ, 0., 0.);
368 xPos = Gauss3D->GetParameter(6);
369 yPos = Gauss3D->GetParameter(7);
370 zPos = Gauss3D->GetParameter(8);
373 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(Gauss3D->GetParameter(0)) + std::fabs(Gauss3D->GetParameter(1))) / 2.;
376 goodData = Gauss3D->ExecuteCommand(
"MIGRAD",arglist,2);
377 Gauss3D->GetStats(amin, edm, errdef, nvpar, nparx);
381 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
edm::isNotFinite(Gauss3D->GetParError(
j)) ==
true) { goodData = -1;
break; }
384 covyz = Gauss3D->GetParameter(4)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(1))) - Gauss3D->GetParameter(5)*Gauss3D->GetParameter(3);
385 covxz = Gauss3D->GetParameter(5)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(0))) - Gauss3D->GetParameter(4)*Gauss3D->GetParameter(3);
387 det = std::fabs(Gauss3D->GetParameter(0)) * (std::fabs(Gauss3D->GetParameter(1))*std::fabs(Gauss3D->GetParameter(2)) - covyz*covyz) -
388 Gauss3D->GetParameter(3) * (Gauss3D->GetParameter(3)*std::fabs(Gauss3D->GetParameter(2)) - covxz*covyz) +
389 covxz * (Gauss3D->GetParameter(3)*covyz - covxz*std::fabs(Gauss3D->GetParameter(1)));
390 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
393 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementY =
i; }
395 if (
internalDebug ==
true)
cout <<
"Found bestMovementY --> " << bestMovementY << endl;
399 for (
int i = 0;
i < 3;
i++)
401 deltaMean = (double(
i)-1.)*
std::sqrt(*(it+2));
404 cout <<
"deltaMean --> " << deltaMean << endl;
405 cout <<
"deltaMean X --> " << (double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor) << endl;
406 cout <<
"deltaMean Y --> " << (double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor) << endl;
413 Gauss3D->SetParameter(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
414 Gauss3D->SetParameter(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
415 Gauss3D->SetParameter(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ, 0., 0.);
416 Gauss3D->SetParameter(3,
"cov xy", *(it+3), parDistanceCxy, 0., 0.);
417 Gauss3D->SetParameter(4,
"dydz ", *(it+4), parDistanceddZ, 0., 0.);
418 Gauss3D->SetParameter(5,
"dxdz ", *(it+5), parDistanceddZ, 0., 0.);
419 Gauss3D->SetParameter(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY, 0., 0.);
420 Gauss3D->SetParameter(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor), parDistanceXY, 0., 0.);
421 Gauss3D->SetParameter(8,
"mean z", *(it+8)+deltaMean, parDistanceZ, 0., 0.);
424 xPos = Gauss3D->GetParameter(6);
425 yPos = Gauss3D->GetParameter(7);
426 zPos = Gauss3D->GetParameter(8);
429 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(Gauss3D->GetParameter(0)) + std::fabs(Gauss3D->GetParameter(1))) / 2.;
432 goodData = Gauss3D->ExecuteCommand(
"MIGRAD",arglist,2);
433 Gauss3D->GetStats(amin, edm, errdef, nvpar, nparx);
437 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
edm::isNotFinite(Gauss3D->GetParError(
j)) ==
true) { goodData = -1;
break; }
440 covyz = Gauss3D->GetParameter(4)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(1))) - Gauss3D->GetParameter(5)*Gauss3D->GetParameter(3);
441 covxz = Gauss3D->GetParameter(5)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(0))) - Gauss3D->GetParameter(4)*Gauss3D->GetParameter(3);
443 det = std::fabs(Gauss3D->GetParameter(0)) * (std::fabs(Gauss3D->GetParameter(1))*std::fabs(Gauss3D->GetParameter(2)) - covyz*covyz) -
444 Gauss3D->GetParameter(3) * (Gauss3D->GetParameter(3)*std::fabs(Gauss3D->GetParameter(2)) - covxz*covyz) +
445 covxz * (Gauss3D->GetParameter(3)*covyz - covxz*std::fabs(Gauss3D->GetParameter(1)));
446 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
449 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementZ =
i; }
451 if (
internalDebug ==
true)
cout <<
"Found bestMovementZ --> " << bestMovementZ << endl;
458 Gauss3D->SetParameter(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
459 Gauss3D->SetParameter(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
460 Gauss3D->SetParameter(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ, 0., 0.);
461 Gauss3D->SetParameter(3,
"cov xy", *(it+3), parDistanceCxy, 0., 0.);
462 Gauss3D->SetParameter(4,
"dydz ", *(it+4), parDistanceddZ, 0., 0.);
463 Gauss3D->SetParameter(5,
"dxdz ", *(it+5), parDistanceddZ, 0., 0.);
464 Gauss3D->SetParameter(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY, 0., 0.);
465 Gauss3D->SetParameter(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor), parDistanceXY, 0., 0.);
466 Gauss3D->SetParameter(8,
"mean z", *(it+8)+(
double(bestMovementZ)-1.)*
std::sqrt(*(it+2)), parDistanceZ, 0., 0.);
469 xPos = Gauss3D->GetParameter(6);
470 yPos = Gauss3D->GetParameter(7);
471 zPos = Gauss3D->GetParameter(8);
474 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(Gauss3D->GetParameter(0)) + std::fabs(Gauss3D->GetParameter(1))) / 2.;
477 goodData = Gauss3D->ExecuteCommand(
"MIGRAD",arglist,2);
478 Gauss3D->GetStats(amin, edm, errdef, nvpar, nparx);
482 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
edm::isNotFinite(Gauss3D->GetParError(
j)) ==
true) { goodData = -1;
break; }
485 covyz = Gauss3D->GetParameter(4)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(1))) - Gauss3D->GetParameter(5)*Gauss3D->GetParameter(3);
486 covxz = Gauss3D->GetParameter(5)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(0))) - Gauss3D->GetParameter(4)*Gauss3D->GetParameter(3);
488 det = std::fabs(Gauss3D->GetParameter(0)) * (std::fabs(Gauss3D->GetParameter(1))*std::fabs(Gauss3D->GetParameter(2)) - covyz*covyz) -
489 Gauss3D->GetParameter(3) * (Gauss3D->GetParameter(3)*std::fabs(Gauss3D->GetParameter(2)) - covxz*covyz) +
490 covxz * (Gauss3D->GetParameter(3)*covyz - covxz*std::fabs(Gauss3D->GetParameter(1)));
491 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
497 for (
unsigned int i = 0;
i < trials;
i++)
499 if ((goodData != 0) && (goodData != -2))
503 if (
internalDebug ==
true)
cout <<
"FIT WITH DIFFERENT PARAMETER DISTANCES - STEP " <<
i+1 << endl;
505 Gauss3D->SetParameter(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY * largerDist[
i], 0, 0);
506 Gauss3D->SetParameter(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY * largerDist[i], 0, 0);
507 Gauss3D->SetParameter(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ * largerDist[i], 0, 0);
508 Gauss3D->SetParameter(3,
"cov xy", *(it+3), parDistanceCxy * largerDist[i], 0, 0);
509 Gauss3D->SetParameter(4,
"dydz ", *(it+4), parDistanceddZ * largerDist[i], 0, 0);
510 Gauss3D->SetParameter(5,
"dxdz ", *(it+5), parDistanceddZ * largerDist[i], 0, 0);
511 Gauss3D->SetParameter(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY * largerDist[i], 0, 0);
512 Gauss3D->SetParameter(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor), parDistanceXY * largerDist[i], 0, 0);
513 Gauss3D->SetParameter(8,
"mean z", *(it+8)+(
double(bestMovementZ)-1.)*
std::sqrt(*(it+2)), parDistanceZ * largerDist[i], 0, 0);
516 xPos = Gauss3D->GetParameter(6);
517 yPos = Gauss3D->GetParameter(7);
518 zPos = Gauss3D->GetParameter(8);
521 maxTransRadius = nSigmaXY *
std::sqrt(std::fabs(Gauss3D->GetParameter(0)) + std::fabs(Gauss3D->GetParameter(1))) / 2.;
524 goodData = Gauss3D->ExecuteCommand(
"MIGRAD",arglist,2);
525 Gauss3D->GetStats(amin, edm, errdef, nvpar, nparx);
529 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
edm::isNotFinite(Gauss3D->GetParError(
j)) ==
true) { goodData = -1;
break; }
532 covyz = Gauss3D->GetParameter(4)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(1))) - Gauss3D->GetParameter(5)*Gauss3D->GetParameter(3);
533 covxz = Gauss3D->GetParameter(5)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(0))) - Gauss3D->GetParameter(4)*Gauss3D->GetParameter(3);
535 det = std::fabs(Gauss3D->GetParameter(0)) * (std::fabs(Gauss3D->GetParameter(1))*std::fabs(Gauss3D->GetParameter(2)) - covyz*covyz) -
536 Gauss3D->GetParameter(3) * (Gauss3D->GetParameter(3)*std::fabs(Gauss3D->GetParameter(2)) - covxz*covyz) +
537 covxz * (Gauss3D->GetParameter(3)*covyz - covxz*std::fabs(Gauss3D->GetParameter(1)));
538 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
544 for (
unsigned int i = 0; i < nParams; i++)
546 vals->operator[](
i) = Gauss3D->GetParameter(i);
547 vals->operator[](i+nParams) = Gauss3D->GetParError(i);
void Gauss3DFunc(int &, double *, double &fval, double *par, int)