21 #include <Math/Minimizer.h>
22 #include <Math/Factory.h>
23 #include <Math/Functor.h>
80 stringstream debugFile;
84 tmp.erase(strlen(
fileName.c_str())-4,4);
85 debugFile << tmp.c_str() <<
"_Run" << iEvent.
id().
run() <<
".txt";
97 for (vector<Vertex>::const_iterator it3DVx = Vx3DCollection->begin(); it3DVx != Vx3DCollection->end(); it3DVx++)
99 if ((it3DVx->isValid() ==
true) &&
100 (it3DVx->isFake() ==
false) &&
103 for (i = 0; i <
DIM; i++)
105 for (j = 0; j <
DIM; j++)
115 if ((i == DIM) && (det > 0.))
117 MyVertex.
x = it3DVx->x();
118 MyVertex.
y = it3DVx->y();
119 MyVertex.
z = it3DVx->z();
124 cout <<
"Vertex discarded !" << endl;
125 for (i = 0; i <
DIM; i++)
126 for (j = 0; j <
DIM; j++)
127 cout <<
"(i,j) --> " << i <<
"," << j <<
" --> " << MyVertex.
Covariance[i][j] << endl;
160 strftime(ts,
sizeof(ts),
"%Y.%m.%d %H:%M:%S %Z", gmtime(&t));
162 string ts_string(ts);
197 K[1][2] = K[2][1] = par[4]*(std::fabs(par[2])-std::fabs(par[1])) - par[5]*par[3] +
VxErrCorr*
VxErrCorr *
Vertices[
i].Covariance[1][2];
198 K[0][2] = K[2][0] = par[5]*(std::fabs(par[2])-std::fabs(par[0])) - par[4]*par[3] +
VxErrCorr*
VxErrCorr *
Vertices[
i].Covariance[0][2];
202 K[0][0] = std::fabs(par[0]);
203 K[1][1] = std::fabs(par[1]);
204 K[2][2] = std::fabs(par[2]);
205 K[0][1] = K[1][0] = par[3];
206 K[1][2] = K[2][1] = par[4]*(std::fabs(par[2])-std::fabs(par[1])) - par[5]*par[3];
207 K[0][2] = K[2][0] = par[5]*(std::fabs(par[2])-std::fabs(par[0])) - par[4]*par[3];
210 det = K[0][0]*(K[1][1]*K[2][2] - K[1][2]*K[1][2]) -
211 K[0][1]*(K[0][1]*K[2][2] - K[0][2]*K[1][2]) +
212 K[0][2]*(K[0][1]*K[1][2] - K[0][2]*K[1][1]);
214 M[0][0] = (K[1][1]*K[2][2] - K[1][2]*K[1][2]) / det;
215 M[1][1] = (K[0][0]*K[2][2] - K[0][2]*K[0][2]) / det;
216 M[2][2] = (K[0][0]*K[1][1] - K[0][1]*K[0][1]) / det;
217 M[0][1] = M[1][0] = (K[0][2]*K[1][2] - K[0][1]*K[2][2]) / det;
218 M[1][2] = M[2][1] = (K[0][2]*K[0][1] - K[1][2]*K[0][0]) / det;
219 M[0][2] = M[2][0] = (K[0][1]*K[1][2] - K[0][2]*K[1][1]) / det;
243 unsigned int nParams = 9;
245 if ((vals !=
NULL) && (vals->size() == nParams*2))
247 double nSigmaXY = 100.;
248 double nSigmaZ = 100.;
249 double varFactor = 4./25.;
250 double parDistanceXY = 0.005;
251 double parDistanceZ = 0.5;
252 double parDistanceddZ = 1
e-3;
253 double parDistanceCxy = 1
e-5;
254 double bestEdm = 1
e-1;
256 const unsigned int trials = 4;
257 double largerDist[trials] = {0.1, 5., 10., 100.};
259 double covxz,covyz,det;
261 int bestMovementX = 1;
262 int bestMovementY = 1;
263 int bestMovementZ = 1;
268 vector<double>::const_iterator it = vals->begin();
270 ROOT::Math::Minimizer* Gauss3D = ROOT::Math::Factory::CreateMinimizer(
"Minuit2",
"Migrad");
271 Gauss3D->SetMaxFunctionCalls(1e4);
272 Gauss3D->SetTolerance(1
e-9);
274 else Gauss3D->SetPrintLevel(0);
276 ROOT::Math::Functor _Gauss3DFunc(&
Gauss3DFunc,nParams);
277 Gauss3D->SetFunction(_Gauss3DFunc);
283 for (
int i = 0;
i < 3;
i++)
285 deltaMean = (double(
i)-1.)*
std::sqrt((*(it+0))*varFactor);
290 Gauss3D->SetVariable(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY);
291 Gauss3D->SetVariable(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY);
292 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ);
293 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
294 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
295 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
296 Gauss3D->SetVariable(6,
"mean x", *(it+6)+deltaMean, parDistanceXY);
297 Gauss3D->SetVariable(7,
"mean y", *(it+7), parDistanceXY);
298 Gauss3D->SetVariable(8,
"mean z", *(it+8), parDistanceZ);
301 xPos = Gauss3D->X()[6];
302 yPos = Gauss3D->X()[7];
303 zPos = Gauss3D->X()[8];
310 goodData = Gauss3D->Status();
311 edm = Gauss3D->Edm();
315 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
isNotFinite(Gauss3D->Errors()[
j]) ==
true) { goodData = -1;
break; }
318 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
319 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
321 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
322 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
323 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
324 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
327 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementX =
i; }
329 if (
internalDebug ==
true)
cout <<
"Found bestMovementX --> " << bestMovementX << endl;
333 for (
int i = 0;
i < 3;
i++)
335 deltaMean = (double(
i)-1.)*
std::sqrt((*(it+1))*varFactor);
338 cout <<
"deltaMean --> " << deltaMean << endl;
339 cout <<
"deltaMean X --> " << (double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor) << endl;
344 Gauss3D->SetVariable(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY);
345 Gauss3D->SetVariable(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY);
346 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ);
347 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
348 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
349 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
350 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY);
351 Gauss3D->SetVariable(7,
"mean y", *(it+7)+deltaMean, parDistanceXY);
352 Gauss3D->SetVariable(8,
"mean z", *(it+8), parDistanceZ);
355 xPos = Gauss3D->X()[6];
356 yPos = Gauss3D->X()[7];
357 zPos = Gauss3D->X()[8];
364 goodData = Gauss3D->Status();
365 edm = Gauss3D->Edm();
369 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
isNotFinite(Gauss3D->Errors()[
j]) ==
true) { goodData = -1;
break; }
372 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
373 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
375 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
376 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
377 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
378 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
381 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementY =
i; }
383 if (
internalDebug ==
true)
cout <<
"Found bestMovementY --> " << bestMovementY << endl;
387 for (
int i = 0;
i < 3;
i++)
389 deltaMean = (double(
i)-1.)*
std::sqrt(*(it+2));
392 cout <<
"deltaMean --> " << deltaMean << endl;
393 cout <<
"deltaMean X --> " << (double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor) << endl;
394 cout <<
"deltaMean Y --> " << (double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor) << endl;
399 Gauss3D->SetVariable(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY);
400 Gauss3D->SetVariable(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY);
401 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ);
402 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
403 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
404 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
405 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY);
406 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor), parDistanceXY);
407 Gauss3D->SetVariable(8,
"mean z", *(it+8)+deltaMean, parDistanceZ);
410 xPos = Gauss3D->X()[6];
411 yPos = Gauss3D->X()[7];
412 zPos = Gauss3D->X()[8];
419 goodData = Gauss3D->Status();
420 edm = Gauss3D->Edm();
424 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
isNotFinite(Gauss3D->Errors()[
j]) ==
true) { goodData = -1;
break; }
427 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
428 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
430 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
431 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
432 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
433 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
436 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementZ =
i; }
438 if (
internalDebug ==
true)
cout <<
"Found bestMovementZ --> " << bestMovementZ << endl;
443 Gauss3D->SetVariable(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY);
444 Gauss3D->SetVariable(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY);
445 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ);
446 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
447 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
448 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
449 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY);
450 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor), parDistanceXY);
451 Gauss3D->SetVariable(8,
"mean z", *(it+8)+(
double(bestMovementZ)-1.)*
std::sqrt(*(it+2)), parDistanceZ);
454 xPos = Gauss3D->X()[6];
455 yPos = Gauss3D->X()[7];
456 zPos = Gauss3D->X()[8];
463 goodData = Gauss3D->Status();
464 edm = Gauss3D->Edm();
468 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
isNotFinite(Gauss3D->Errors()[
j]) ==
true) { goodData = -1;
break; }
471 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
472 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
474 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
475 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
476 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
477 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
481 for (
unsigned int i = 0;
i < trials;
i++)
483 if ((goodData != 0) && (goodData != -2))
487 if (
internalDebug ==
true)
cout <<
"FIT WITH DIFFERENT PARAMETER DISTANCES - STEP " <<
i+1 << endl;
489 Gauss3D->SetVariable(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY * largerDist[
i]);
490 Gauss3D->SetVariable(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY * largerDist[i]);
491 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ * largerDist[i]);
492 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy * largerDist[i]);
493 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ * largerDist[i]);
494 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ * largerDist[i]);
495 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor), parDistanceXY * largerDist[i]);
496 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*std::sqrt((*(it+1))*varFactor), parDistanceXY * largerDist[i]);
497 Gauss3D->SetVariable(8,
"mean z", *(it+8)+(
double(bestMovementZ)-1.)*std::sqrt(*(it+2)), parDistanceZ * largerDist[i]);
500 xPos = Gauss3D->X()[6];
501 yPos = Gauss3D->X()[7];
502 zPos = Gauss3D->X()[8];
509 goodData = Gauss3D->Status();
510 edm = Gauss3D->Edm();
514 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
isNotFinite(Gauss3D->Errors()[
j]) ==
true) { goodData = -1;
break; }
517 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
518 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
520 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
521 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
522 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
523 if (det < 0.) { goodData = -1;
if (
internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
529 for (
unsigned int i = 0;
i < nParams;
i++)
531 vals->operator[](
i) = Gauss3D->X()[
i];
532 vals->operator[](
i+nParams) = Gauss3D->Errors()[
i];
545 if (ResetType.compare(
"scratch") == 0)
590 else if (ResetType.compare(
"whole") == 0)
610 else if (ResetType.compare(
"partial") == 0)
625 else if (ResetType.compare(
"nohisto") == 0)
637 else if (ResetType.compare(
"hitCounter") == 0)
645 unsigned int BeginLumiOfFit,
646 unsigned int EndLumiOfFit,
649 stringstream BufferString;
650 BufferString.precision(5);
654 if ((
outputFile.is_open() ==
true) && (vals !=
NULL) && (vals->size() == 8*2))
656 vector<double>::const_iterator it = vals->begin();
667 BufferString << *(it+0);
668 outputFile <<
"X0 " << BufferString.str().c_str() << endl;
669 BufferString.str(
"");
671 BufferString << *(it+1);
672 outputFile <<
"Y0 " << BufferString.str().c_str() << endl;
673 BufferString.str(
"");
675 BufferString << *(it+2);
676 outputFile <<
"Z0 " << BufferString.str().c_str() << endl;
677 BufferString.str(
"");
679 BufferString << *(it+3);
680 outputFile <<
"sigmaZ0 " << BufferString.str().c_str() << endl;
681 BufferString.str(
"");
683 BufferString << *(it+4);
684 outputFile <<
"dxdz " << BufferString.str().c_str() << endl;
685 BufferString.str(
"");
687 BufferString << *(it+5);
688 outputFile <<
"dydz " << BufferString.str().c_str() << endl;
689 BufferString.str(
"");
691 BufferString << *(it+6);
692 outputFile <<
"BeamWidthX " << BufferString.str().c_str() << endl;
693 BufferString.str(
"");
695 BufferString << *(it+7);
696 outputFile <<
"BeamWidthY " << BufferString.str().c_str() << endl;
697 BufferString.str(
"");
699 outputFile <<
"Cov(0,j) " << *(it+8) <<
" 0.0 0.0 0.0 0.0 0.0 0.0" << endl;
700 outputFile <<
"Cov(1,j) 0.0 " << *(it+9) <<
" 0.0 0.0 0.0 0.0 0.0" << endl;
701 outputFile <<
"Cov(2,j) 0.0 0.0 " << *(it+10) <<
" 0.0 0.0 0.0 0.0" << endl;
702 outputFile <<
"Cov(3,j) 0.0 0.0 0.0 " << *(it+11) <<
" 0.0 0.0 0.0" << endl;
703 outputFile <<
"Cov(4,j) 0.0 0.0 0.0 0.0 " << *(it+12) <<
" 0.0 0.0" << endl;
704 outputFile <<
"Cov(5,j) 0.0 0.0 0.0 0.0 0.0 " << *(it+13) <<
" 0.0" << endl;
705 outputFile <<
"Cov(6,j) 0.0 0.0 0.0 0.0 0.0 0.0 " << ((*(it+14)) + (*(it+15)) + 2.*
std::sqrt((*(it+14))*(*(it+15)))) / 4. << endl;
715 vector<double>::const_iterator it = vals->begin();
726 BufferString << *(it+0);
728 BufferString.str(
"");
730 BufferString << *(it+1);
732 BufferString.str(
"");
734 BufferString << *(it+2);
736 BufferString.str(
"");
738 BufferString << *(it+3);
740 BufferString.str(
"");
742 BufferString << *(it+4);
744 BufferString.str(
"");
746 BufferString << *(it+5);
748 BufferString.str(
"");
750 BufferString << *(it+6);
751 outputDebugFile <<
"BeamWidthX " << BufferString.str().c_str() << endl;
752 BufferString.str(
"");
754 BufferString << *(it+7);
755 outputDebugFile <<
"BeamWidthY " << BufferString.str().c_str() << endl;
756 BufferString.str(
"");
758 outputDebugFile <<
"Cov(0,j) " << *(it+8) <<
" 0.0 0.0 0.0 0.0 0.0 0.0" << endl;
759 outputDebugFile <<
"Cov(1,j) 0.0 " << *(it+9) <<
" 0.0 0.0 0.0 0.0 0.0" << endl;
760 outputDebugFile <<
"Cov(2,j) 0.0 0.0 " << *(it+10) <<
" 0.0 0.0 0.0 0.0" << endl;
761 outputDebugFile <<
"Cov(3,j) 0.0 0.0 0.0 " << *(it+11) <<
" 0.0 0.0 0.0" << endl;
762 outputDebugFile <<
"Cov(4,j) 0.0 0.0 0.0 0.0 " << *(it+12) <<
" 0.0 0.0" << endl;
763 outputDebugFile <<
"Cov(5,j) 0.0 0.0 0.0 0.0 0.0 " << *(it+13) <<
" 0.0" << endl;
764 outputDebugFile <<
"Cov(6,j) 0.0 0.0 0.0 0.0 0.0 0.0 " << ((*(it+14)) + (*(it+15)) + 2.*
std::sqrt((*(it+14))*(*(it+15)))) / 4. << endl;
788 stringstream histTitle;
790 unsigned int nParams = 9;
814 fitResults.push_back(0.0);
815 fitResults.push_back(0.0);
816 fitResults.push_back(0.0);
817 fitResults.push_back(
Vx_X->
getTH1()->GetMean());
818 fitResults.push_back(
Vx_Y->
getTH1()->GetMean());
819 fitResults.push_back(
Vx_Z->
getTH1()->GetMean());
820 for (
unsigned int i = 0;
i < nParams;
i++) fitResults.push_back(0.0);
822 goodData =
MyFit(&fitResults);
826 cout <<
"goodData --> " << goodData << endl;
828 cout <<
"var x --> " << fitResults[0] <<
" +/- " << fitResults[0+nParams] << endl;
829 cout <<
"var y --> " << fitResults[1] <<
" +/- " << fitResults[1+nParams] << endl;
830 cout <<
"var z --> " << fitResults[2] <<
" +/- " << fitResults[2+nParams] << endl;
831 cout <<
"cov xy --> " << fitResults[3] <<
" +/- " << fitResults[3+nParams] << endl;
832 cout <<
"dydz --> " << fitResults[4] <<
" +/- " << fitResults[4+nParams] << endl;
833 cout <<
"dxdz --> " << fitResults[5] <<
" +/- " << fitResults[5+nParams] << endl;
834 cout <<
"mean x --> " << fitResults[6] <<
" +/- " << fitResults[6+nParams] << endl;
835 cout <<
"mean y --> " << fitResults[7] <<
" +/- " << fitResults[7+nParams] << endl;
836 cout <<
"mean z --> " << fitResults[8] <<
" +/- " << fitResults[8+nParams] << endl;
841 vals.push_back(fitResults[6]);
842 vals.push_back(fitResults[7]);
843 vals.push_back(fitResults[8]);
844 vals.push_back(
std::sqrt(std::fabs(fitResults[2])));
845 vals.push_back(fitResults[5]);
846 vals.push_back(fitResults[4]);
847 vals.push_back(
std::sqrt(std::fabs(fitResults[0])));
848 vals.push_back(
std::sqrt(std::fabs(fitResults[1])));
850 vals.push_back(
std::pow(fitResults[6+nParams],2.));
851 vals.push_back(
std::pow(fitResults[7+nParams],2.));
852 vals.push_back(
std::pow(fitResults[8+nParams],2.));
853 vals.push_back(
std::pow(std::fabs(fitResults[2+nParams]) / (2.*
std::sqrt(std::fabs(fitResults[2]))),2.));
854 vals.push_back(
std::pow(fitResults[5+nParams],2.));
855 vals.push_back(
std::pow(fitResults[4+nParams],2.));
856 vals.push_back(
std::pow(std::fabs(fitResults[0+nParams]) / (2.*
std::sqrt(std::fabs(fitResults[0]))),2.));
857 vals.push_back(
std::pow(std::fabs(fitResults[1+nParams]) / (2.*
std::sqrt(std::fabs(fitResults[1]))),2.));
859 else for (
unsigned int i = 0;
i < 8*2;
i++) vals.push_back(0.0);
892 for (
unsigned int i = 0;
i < 8*2;
i++) vals.push_back(0.0);
929 else reset(
"partial");
938 histTitle <<
"Fitted Beam Spot [cm] (not enough statistics)";
940 else reset(
"hitCounter");
944 histTitle <<
"Fitted Beam Spot [cm] (problems)";
946 else reset(
"partial");
978 TF1* myLinFit =
new TF1(
"myLinFit",
"[0] + [1]*x",
mXlumi->
getTH1()->GetXaxis()->GetXmin(),
mXlumi->
getTH1()->GetXaxis()->GetXmax());
979 myLinFit->SetLineColor(2);
980 myLinFit->SetLineWidth(2);
981 myLinFit->SetParName(0,
"Intercept");
982 myLinFit->SetParName(1,
"Slope");
985 myLinFit->SetParameter(0,
mXlumi->
getTH1()->GetMean(2));
986 myLinFit->SetParameter(1, 0.0);
990 myLinFit->SetParameter(0,
mYlumi->
getTH1()->GetMean(2));
991 myLinFit->SetParameter(1, 0.0);
995 myLinFit->SetParameter(0,
mZlumi->
getTH1()->GetMean(2));
996 myLinFit->SetParameter(1, 0.0);
1000 myLinFit->SetParameter(0,
sXlumi->
getTH1()->GetMean(2));
1001 myLinFit->SetParameter(1, 0.0);
1005 myLinFit->SetParameter(0,
sYlumi->
getTH1()->GetMean(2));
1006 myLinFit->SetParameter(1, 0.0);
1010 myLinFit->SetParameter(0,
sZlumi->
getTH1()->GetMean(2));
1011 myLinFit->SetParameter(1, 0.0);
1016 myLinFit->SetParameter(1, 0.0);
1021 myLinFit->SetParameter(1, 0.0);
1026 myLinFit->SetParameter(1, 0.0);
1033 myExpFit->SetLineColor(2);
1034 myExpFit->SetLineWidth(2);
1035 myExpFit->SetParName(0,
"Amplitude");
1036 myExpFit->SetParName(1,
"#tau");
1058 histTitle <<
"Fitted Beam Spot [cm] (no ongoing fits)";
1077 pi = 3.141592653589793238;
1097 dbe->setCurrentFolder(
"BeamPixel");
1177 fitResults = dbe->book2D(
"fit results",
"Results of Beam Spot Fit", 2, 0., 2., 9, 0., 9.);
1192 dbe->setCurrentFolder(
"BeamPixel/EventInfo");
1195 reportSummaryMap = dbe->book2D(
"reportSummaryMap",
"Pixel-Vertices Beam Spot: % Good Fits", 1, 0., 1., 1, 0., 1.);
1197 dbe->setCurrentFolder(
"BeamPixel/EventInfo/reportSummaryContents");
EDGetTokenT< SiPixelRecHitCollection > pixelHitCollection
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
MonitorElement * dxdzlumi
boost::transform_iterator< IterHelp, const_IdIter > const_iterator
void setBinContent(int binx, double content)
set content of bin (1-D)
virtual unsigned int HitCounter(const Event &iEvent)
unsigned int maxLumiIntegration
bool getByToken(EDGetToken token, Handle< PROD > &result) const
vector< VertexType > Vertices
#define DEFINE_FWK_MODULE(type)
unsigned int beginLumiOfFit
MonitorElement * hitCountHistory
void setBinLabel(int bin, const std::string &label, int axis=1)
set bin label for x, y or z axis (axis=1, 2, 3 respectively)
virtual void beginRun(const Run &iRun, const EventSetup &iSetup)
Vx3DHLTAnalyzer(const ParameterSet &)
data_type const * const_iterator
Timestamp const & beginTime() const
virtual void beginLuminosityBlock(const LuminosityBlock &lumiBlock, const EventSetup &iSetup)
MonitorElement * reportSummaryMap
virtual void endLuminosityBlock(const LuminosityBlock &lumiBlock, const EventSetup &iSetup)
unsigned int nBinsHistoricalPlot
LuminosityBlockNumber_t luminosityBlock() const
double Gauss3DFunc(const double *par)
bool considerVxCovariance
void ShiftFillLast(double y, double ye=0., int32_t xscale=1)
MonitorElement * fitResults
virtual string formatTime(const time_t &t)
unsigned int endLumiOfFit
unsigned int numberGoodFits
MonitorElement * reportSummary
Timestamp const & endTime() const
MonitorElement * goodVxCounter
virtual int MyFit(vector< double > *vals)
MonitorElement * dydzlumi
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
LuminosityBlock const & getLuminosityBlock() const
unsigned long long TimeValue_t
unsigned int lastLumiOfFit
unsigned int nBinsWholeHistory
unsigned int lumiCounterHisto
TimeValue_t beginTimeOfFit
virtual void analyze(const Event &, const EventSetup &)
TH1F * getTH1F(void) const
EDGetTokenT< VertexCollection > vertexCollection
std::vector< std::vector< double > > tmp
static std::atomic< unsigned int > counter
virtual void writeToFile(vector< double > *vals, TimeValue_t BeginTimeOfFit, TimeValue_t EndTimeOfFit, unsigned int BeginLumiOfFit, unsigned int EndLumiOfFit, int dataType)
unsigned int prescaleHistory
MonitorElement * goodVxCountHistory
MonitorElement * hitCounter
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
void Reset(void)
reset ME (ie. contents, errors, etc)
TimeValue_t value() const
virtual void reset(string ResetType)
Power< A, B >::type pow(const A &a, const B &b)