25 #include <Math/Minimizer.h>
26 #include <Math/Factory.h>
27 #include <Math/Functor.h>
90 stringstream debugFile;
93 if (outputDebugFile.is_open() ==
true) outputDebugFile.close();
94 tmp.erase(strlen(
fileName.c_str())-4,4);
95 debugFile << tmp.c_str() <<
"_Run" << iEvent.
id().
run() <<
".txt";
96 outputDebugFile.open(debugFile.str().c_str(),
ios::out);
97 outputDebugFile.close();
98 outputDebugFile.open(debugFile.str().c_str(), ios::app);
103 else if (beginTimeOfFit != 0)
105 totalHits += HitCounter(iEvent);
107 for (vector<Vertex>::const_iterator it3DVx = Vx3DCollection->begin(); it3DVx != Vx3DCollection->end(); it3DVx++)
109 if ((it3DVx->isValid() ==
true) &&
110 (it3DVx->isFake() ==
false) &&
111 (it3DVx->ndof() >= minVxDoF))
113 for (i = 0; i <
DIM; i++)
115 for (j = 0; j <
DIM; j++)
125 if ((i == DIM) && (det > 0.))
127 MyVertex.
x = it3DVx->x();
128 MyVertex.
y = it3DVx->y();
129 MyVertex.
z = it3DVx->z();
132 else if (internalDebug ==
true)
134 cout <<
"Vertex discarded !" << endl;
135 for (i = 0; i <
DIM; i++)
136 for (j = 0; j <
DIM; j++)
137 cout <<
"(i,j) --> " << i <<
"," << j <<
" --> " << MyVertex.
Covariance[i][j] << endl;
140 Vx_X->Fill(it3DVx->x());
141 Vx_Y->Fill(it3DVx->y());
142 Vx_Z->Fill(it3DVx->z());
144 Vx_ZX->Fill(it3DVx->z(), it3DVx->x());
145 Vx_ZY->Fill(it3DVx->z(), it3DVx->y());
146 Vx_XY->Fill(it3DVx->x(), it3DVx->y());
156 iEvent.
getByToken(pixelHitCollection, rechitspixel);
170 strftime(ts,
sizeof(ts),
"%Y.%m.%d %H:%M:%S %Z", gmtime(&t));
205 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];
206 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];
210 K[0][0] = std::fabs(par[0]);
211 K[1][1] = std::fabs(par[1]);
212 K[2][2] = std::fabs(par[2]);
213 K[0][1] = K[1][0] = par[3];
214 K[1][2] = K[2][1] = par[4]*(std::fabs(par[2])-std::fabs(par[1])) - par[5]*par[3];
215 K[0][2] = K[2][0] = par[5]*(std::fabs(par[2])-std::fabs(par[0])) - par[4]*par[3];
218 det = K[0][0]*(K[1][1]*K[2][2] - K[1][2]*K[1][2]) -
219 K[0][1]*(K[0][1]*K[2][2] - K[0][2]*K[1][2]) +
220 K[0][2]*(K[0][1]*K[1][2] - K[0][2]*K[1][1]);
222 M[0][0] = (K[1][1]*K[2][2] - K[1][2]*K[1][2]) / det;
223 M[1][1] = (K[0][0]*K[2][2] - K[0][2]*K[0][2]) / det;
224 M[2][2] = (K[0][0]*K[1][1] - K[0][1]*K[0][1]) / det;
225 M[0][1] = M[1][0] = (K[0][2]*K[1][2] - K[0][1]*K[2][2]) / det;
226 M[1][2] = M[2][1] = (K[0][2]*K[0][1] - K[1][2]*K[0][0]) / det;
227 M[0][2] = M[2][0] = (K[0][1]*K[1][2] - K[0][2]*K[1][1]) / det;
251 unsigned int nParams = 9;
253 if ((vals !=
NULL) && (vals->size() == nParams*2))
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);
281 if (internalDebug ==
true) Gauss3D->SetPrintLevel(3);
282 else Gauss3D->SetPrintLevel(0);
284 ROOT::Math::Functor _Gauss3DFunc(&
Gauss3DFunc,nParams);
285 Gauss3D->SetFunction(_Gauss3DFunc);
287 if (internalDebug ==
true)
cout <<
"\n@@@ START FITTING @@@" << endl;
291 for (
int i = 0;
i < 3;
i++)
293 deltaMean = (double(
i)-1.)*
std::sqrt((*(it+0))*varFactor);
294 if (internalDebug ==
true)
cout <<
"deltaMean --> " << deltaMean << endl;
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);
344 if (internalDebug ==
true)
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));
398 if (internalDebug ==
true)
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];
553 if (ResetType.compare(
"scratch") == 0)
580 hitCountHistory->Reset();
581 goodVxCounter->Reset();
582 goodVxCountHistory->Reset();
586 reportSummaryMap->Fill(0.5, 0.5, 0.);
591 lumiCounterHisto = 0;
598 else if (ResetType.compare(
"whole") == 0)
611 lumiCounterHisto = 0;
618 else if (ResetType.compare(
"partial") == 0)
633 else if (ResetType.compare(
"nohisto") == 0)
638 lumiCounterHisto = 0;
645 else if (ResetType.compare(
"hitCounter") == 0)
653 unsigned int BeginLumiOfFit,
654 unsigned int EndLumiOfFit,
657 stringstream BufferString;
658 BufferString.precision(5);
662 if ((
outputFile.is_open() ==
true) && (vals !=
NULL) && (vals->size() == 8*2))
664 vector<double>::const_iterator it = vals->begin();
667 outputFile <<
"BeginTimeOfFit " <<
formatTime(beginTimeOfFit >> 32) <<
" " << (beginTimeOfFit >> 32) << endl;
668 outputFile <<
"EndTimeOfFit " <<
formatTime(endTimeOfFit >> 32) <<
" " << (endTimeOfFit >> 32) << endl;
669 outputFile <<
"LumiRange " << beginLumiOfFit <<
" - " << endLumiOfFit << endl;
675 BufferString << *(it+0);
676 outputFile <<
"X0 " << BufferString.str().c_str() << endl;
677 BufferString.str(
"");
679 BufferString << *(it+1);
680 outputFile <<
"Y0 " << BufferString.str().c_str() << endl;
681 BufferString.str(
"");
683 BufferString << *(it+2);
684 outputFile <<
"Z0 " << BufferString.str().c_str() << endl;
685 BufferString.str(
"");
687 BufferString << *(it+3);
688 outputFile <<
"sigmaZ0 " << BufferString.str().c_str() << endl;
689 BufferString.str(
"");
691 BufferString << *(it+4);
692 outputFile <<
"dxdz " << BufferString.str().c_str() << endl;
693 BufferString.str(
"");
695 BufferString << *(it+5);
696 outputFile <<
"dydz " << BufferString.str().c_str() << endl;
697 BufferString.str(
"");
699 BufferString << *(it+6);
700 outputFile <<
"BeamWidthX " << BufferString.str().c_str() << endl;
701 BufferString.str(
"");
703 BufferString << *(it+7);
704 outputFile <<
"BeamWidthY " << BufferString.str().c_str() << endl;
705 BufferString.str(
"");
707 outputFile <<
"Cov(0,j) " << *(it+8) <<
" 0.0 0.0 0.0 0.0 0.0 0.0" << endl;
708 outputFile <<
"Cov(1,j) 0.0 " << *(it+9) <<
" 0.0 0.0 0.0 0.0 0.0" << endl;
709 outputFile <<
"Cov(2,j) 0.0 0.0 " << *(it+10) <<
" 0.0 0.0 0.0 0.0" << endl;
710 outputFile <<
"Cov(3,j) 0.0 0.0 0.0 " << *(it+11) <<
" 0.0 0.0 0.0" << endl;
711 outputFile <<
"Cov(4,j) 0.0 0.0 0.0 0.0 " << *(it+12) <<
" 0.0 0.0" << endl;
712 outputFile <<
"Cov(5,j) 0.0 0.0 0.0 0.0 0.0 " << *(it+13) <<
" 0.0" << endl;
713 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;
721 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true) && (vals !=
NULL) && (vals->size() == 8*2))
723 vector<double>::const_iterator it = vals->begin();
725 outputDebugFile <<
"Runnumber " <<
runNumber << endl;
726 outputDebugFile <<
"BeginTimeOfFit " <<
formatTime(beginTimeOfFit >> 32) <<
" " << (beginTimeOfFit >> 32) << endl;
727 outputDebugFile <<
"EndTimeOfFit " <<
formatTime(endTimeOfFit >> 32) <<
" " << (endTimeOfFit >> 32) << endl;
728 outputDebugFile <<
"LumiRange " << beginLumiOfFit <<
" - " << endLumiOfFit << endl;
729 outputDebugFile <<
"Type " << dataType << endl;
734 BufferString << *(it+0);
735 outputDebugFile <<
"X0 " << BufferString.str().c_str() << endl;
736 BufferString.str(
"");
738 BufferString << *(it+1);
739 outputDebugFile <<
"Y0 " << BufferString.str().c_str() << endl;
740 BufferString.str(
"");
742 BufferString << *(it+2);
743 outputDebugFile <<
"Z0 " << BufferString.str().c_str() << endl;
744 BufferString.str(
"");
746 BufferString << *(it+3);
747 outputDebugFile <<
"sigmaZ0 " << BufferString.str().c_str() << endl;
748 BufferString.str(
"");
750 BufferString << *(it+4);
751 outputDebugFile <<
"dxdz " << BufferString.str().c_str() << endl;
752 BufferString.str(
"");
754 BufferString << *(it+5);
755 outputDebugFile <<
"dydz " << BufferString.str().c_str() << endl;
756 BufferString.str(
"");
758 BufferString << *(it+6);
759 outputDebugFile <<
"BeamWidthX " << BufferString.str().c_str() << endl;
760 BufferString.str(
"");
762 BufferString << *(it+7);
763 outputDebugFile <<
"BeamWidthY " << BufferString.str().c_str() << endl;
764 BufferString.str(
"");
766 outputDebugFile <<
"Cov(0,j) " << *(it+8) <<
" 0.0 0.0 0.0 0.0 0.0 0.0" << endl;
767 outputDebugFile <<
"Cov(1,j) 0.0 " << *(it+9) <<
" 0.0 0.0 0.0 0.0 0.0" << endl;
768 outputDebugFile <<
"Cov(2,j) 0.0 0.0 " << *(it+10) <<
" 0.0 0.0 0.0 0.0" << endl;
769 outputDebugFile <<
"Cov(3,j) 0.0 0.0 0.0 " << *(it+11) <<
" 0.0 0.0 0.0" << endl;
770 outputDebugFile <<
"Cov(4,j) 0.0 0.0 0.0 0.0 " << *(it+12) <<
" 0.0 0.0" << endl;
771 outputDebugFile <<
"Cov(5,j) 0.0 0.0 0.0 0.0 0.0 " << *(it+13) <<
" 0.0" << endl;
772 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;
774 outputDebugFile <<
"EmittanceX 0.0" << endl;
775 outputDebugFile <<
"EmittanceY 0.0" << endl;
776 outputDebugFile <<
"BetaStar 0.0" << endl;
784 if ((lumiCounter == 0) && (lumiBlock.
luminosityBlock() > lastLumiOfFit))
791 else if ((lumiCounter != 0) && (lumiBlock.
luminosityBlock() >= (beginLumiOfFit+lumiCounter))) { lumiCounter++; lumiCounterHisto++; }
798 stringstream histTitle;
800 unsigned int nParams = 9;
806 lastLumiOfFit = endLumiOfFit;
811 if (lastLumiOfFit % prescaleHistory == 0)
813 hitCountHistory->getTH1()->SetBinContent(lastLumiOfFit, (
double)totalHits);
814 hitCountHistory->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt((
double)totalHits));
819 vector<double> fitResults;
821 fitResults.push_back(Vx_X->getTH1()->GetRMS()*Vx_X->getTH1()->GetRMS());
822 fitResults.push_back(Vx_Y->getTH1()->GetRMS()*Vx_Y->getTH1()->GetRMS());
823 fitResults.push_back(Vx_Z->getTH1()->GetRMS()*Vx_Z->getTH1()->GetRMS());
824 fitResults.push_back(0.0);
825 fitResults.push_back(0.0);
826 fitResults.push_back(0.0);
827 fitResults.push_back(Vx_X->getTH1()->GetMean());
828 fitResults.push_back(Vx_Y->getTH1()->GetMean());
829 fitResults.push_back(Vx_Z->getTH1()->GetMean());
830 for (
unsigned int i = 0;
i < nParams;
i++) fitResults.push_back(0.0);
832 goodData = MyFit(&fitResults);
834 if (internalDebug ==
true)
836 cout <<
"goodData --> " << goodData << endl;
838 cout <<
"var x --> " << fitResults[0] <<
" +/- " << fitResults[0+nParams] << endl;
839 cout <<
"var y --> " << fitResults[1] <<
" +/- " << fitResults[1+nParams] << endl;
840 cout <<
"var z --> " << fitResults[2] <<
" +/- " << fitResults[2+nParams] << endl;
841 cout <<
"cov xy --> " << fitResults[3] <<
" +/- " << fitResults[3+nParams] << endl;
842 cout <<
"dydz --> " << fitResults[4] <<
" +/- " << fitResults[4+nParams] << endl;
843 cout <<
"dxdz --> " << fitResults[5] <<
" +/- " << fitResults[5+nParams] << endl;
844 cout <<
"mean x --> " << fitResults[6] <<
" +/- " << fitResults[6+nParams] << endl;
845 cout <<
"mean y --> " << fitResults[7] <<
" +/- " << fitResults[7+nParams] << endl;
846 cout <<
"mean z --> " << fitResults[8] <<
" +/- " << fitResults[8+nParams] << endl;
851 vals.push_back(fitResults[6]);
852 vals.push_back(fitResults[7]);
853 vals.push_back(fitResults[8]);
854 vals.push_back(
std::sqrt(std::fabs(fitResults[2])));
855 vals.push_back(fitResults[5]);
856 vals.push_back(fitResults[4]);
857 vals.push_back(
std::sqrt(std::fabs(fitResults[0])));
858 vals.push_back(
std::sqrt(std::fabs(fitResults[1])));
860 vals.push_back(
std::pow(fitResults[6+nParams],2.));
861 vals.push_back(
std::pow(fitResults[7+nParams],2.));
862 vals.push_back(
std::pow(fitResults[8+nParams],2.));
863 vals.push_back(
std::pow(std::fabs(fitResults[2+nParams]) / (2.*
std::sqrt(std::fabs(fitResults[2]))),2.));
864 vals.push_back(
std::pow(fitResults[5+nParams],2.));
865 vals.push_back(
std::pow(fitResults[4+nParams],2.));
866 vals.push_back(
std::pow(std::fabs(fitResults[0+nParams]) / (2.*
std::sqrt(std::fabs(fitResults[0]))),2.));
867 vals.push_back(
std::pow(std::fabs(fitResults[1+nParams]) / (2.*
std::sqrt(std::fabs(fitResults[1]))),2.));
869 else for (
unsigned int i = 0;
i < 8*2;
i++) vals.push_back(0.0);
875 counterVx = Vx_X->getTH1F()->GetEntries();
881 vals.push_back(Vx_X->getTH1F()->GetMean());
882 vals.push_back(Vx_Y->getTH1F()->GetMean());
883 vals.push_back(Vx_Z->getTH1F()->GetMean());
884 vals.push_back(Vx_Z->getTH1F()->GetRMS());
887 vals.push_back(Vx_X->getTH1F()->GetRMS());
888 vals.push_back(Vx_Y->getTH1F()->GetRMS());
890 vals.push_back(
std::pow(Vx_X->getTH1F()->GetMeanError(),2.));
891 vals.push_back(
std::pow(Vx_Y->getTH1F()->GetMeanError(),2.));
892 vals.push_back(
std::pow(Vx_Z->getTH1F()->GetMeanError(),2.));
893 vals.push_back(
std::pow(Vx_Z->getTH1F()->GetRMSError(),2.));
896 vals.push_back(
std::pow(Vx_X->getTH1F()->GetRMSError(),2.));
897 vals.push_back(
std::pow(Vx_Y->getTH1F()->GetRMSError(),2.));
902 for (
unsigned int i = 0;
i < 8*2;
i++) vals.push_back(0.0);
932 writeToFile(&vals, beginTimeOfFit, endTimeOfFit, beginLumiOfFit, endLumiOfFit, 3);
933 if ((internalDebug ==
true) && (outputDebugFile.is_open() ==
true)) outputDebugFile <<
"Used vertices: " <<
counterVx << endl;
937 histTitle <<
"Fitted Beam Spot [cm] (Lumi start: " << beginLumiOfFit <<
" - Lumi end: " << endLumiOfFit <<
")";
938 if (lumiCounterHisto >= maxLumiIntegration)
reset(
"whole");
939 else reset(
"partial");
943 writeToFile(&vals, beginTimeOfFit, endTimeOfFit, beginLumiOfFit, endLumiOfFit, -1);
944 if ((internalDebug ==
true) && (outputDebugFile.is_open() ==
true)) outputDebugFile <<
"Used vertices: " <<
counterVx << endl;
948 histTitle <<
"Fitted Beam Spot [cm] (not enough statistics)";
949 if (lumiCounter >= maxLumiIntegration)
reset(
"whole");
950 else reset(
"hitCounter");
954 histTitle <<
"Fitted Beam Spot [cm] (problems)";
955 if (lumiCounterHisto >= maxLumiIntegration)
reset(
"whole");
956 else reset(
"partial");
962 reportSummary->Fill(numberFits != 0 ? (
double)numberGoodFits/(
double)numberFits : 0.0);
963 reportSummaryMap->Fill(0.5, 0.5, numberFits != 0 ? (
double)numberGoodFits/(
double)numberFits : 0.0);
965 fitResults->setAxisTitle(histTitle.str().c_str(), 1);
967 fitResults->setBinContent(1, 9, vals[0]);
968 fitResults->setBinContent(1, 8, vals[1]);
969 fitResults->setBinContent(1, 7, vals[2]);
970 fitResults->setBinContent(1, 6, vals[3]);
971 fitResults->setBinContent(1, 5, vals[4]);
972 fitResults->setBinContent(1, 4, vals[5]);
973 fitResults->setBinContent(1, 3, vals[6]);
974 fitResults->setBinContent(1, 2, vals[7]);
975 fitResults->setBinContent(1, 1,
counterVx);
977 fitResults->setBinContent(2, 9,
std::sqrt(vals[8]));
978 fitResults->setBinContent(2, 8,
std::sqrt(vals[9]));
979 fitResults->setBinContent(2, 7,
std::sqrt(vals[10]));
980 fitResults->setBinContent(2, 6,
std::sqrt(vals[11]));
981 fitResults->setBinContent(2, 5,
std::sqrt(vals[12]));
982 fitResults->setBinContent(2, 4,
std::sqrt(vals[13]));
983 fitResults->setBinContent(2, 3,
std::sqrt(vals[14]));
984 fitResults->setBinContent(2, 2,
std::sqrt(vals[15]));
988 TF1* myLinFit =
new TF1(
"myLinFit",
"[0] + [1]*x", mXlumi->getTH1()->GetXaxis()->GetXmin(), mXlumi->getTH1()->GetXaxis()->GetXmax());
989 myLinFit->SetLineColor(2);
990 myLinFit->SetLineWidth(2);
991 myLinFit->SetParName(0,
"Intercept");
992 myLinFit->SetParName(1,
"Slope");
995 myLinFit->SetParameter(0, mXlumi->getTH1()->GetMean(2));
996 myLinFit->SetParameter(1, 0.0);
997 mXlumi->getTH1()->Fit(
"myLinFit",
"QR");
1000 myLinFit->SetParameter(0, mYlumi->getTH1()->GetMean(2));
1001 myLinFit->SetParameter(1, 0.0);
1002 mYlumi->getTH1()->Fit(
"myLinFit",
"QR");
1005 myLinFit->SetParameter(0, mZlumi->getTH1()->GetMean(2));
1006 myLinFit->SetParameter(1, 0.0);
1007 mZlumi->getTH1()->Fit(
"myLinFit",
"QR");
1010 myLinFit->SetParameter(0, sXlumi->getTH1()->GetMean(2));
1011 myLinFit->SetParameter(1, 0.0);
1012 sXlumi->getTH1()->Fit(
"myLinFit",
"QR");
1015 myLinFit->SetParameter(0, sYlumi->getTH1()->GetMean(2));
1016 myLinFit->SetParameter(1, 0.0);
1017 sYlumi->getTH1()->Fit(
"myLinFit",
"QR");
1020 myLinFit->SetParameter(0, sZlumi->getTH1()->GetMean(2));
1021 myLinFit->SetParameter(1, 0.0);
1022 sZlumi->getTH1()->Fit(
"myLinFit",
"QR");
1025 myLinFit->SetParameter(0, dxdzlumi->getTH1()->GetMean(2));
1026 myLinFit->SetParameter(1, 0.0);
1027 dxdzlumi->getTH1()->Fit(
"myLinFit",
"QR");
1030 myLinFit->SetParameter(0, dydzlumi->getTH1()->GetMean(2));
1031 myLinFit->SetParameter(1, 0.0);
1032 dydzlumi->getTH1()->Fit(
"myLinFit",
"QR");
1035 myLinFit->SetParameter(0, goodVxCounter->getTH1()->GetMean(2));
1036 myLinFit->SetParameter(1, 0.0);
1037 goodVxCounter->getTH1()->Fit(
"myLinFit",
"QR");
1039 if (lastLumiOfFit % prescaleHistory == 0)
1041 goodVxCountHistory->getTH1()->SetBinContent(lastLumiOfFit, (
double)counterVx);
1042 goodVxCountHistory->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt((
double)counterVx));
1051 histTitle <<
"Fitted Beam Spot [cm] (no ongoing fits)";
1052 fitResults->setAxisTitle(histTitle.str().c_str(), 1);
1053 reportSummaryMap->Fill(0.5, 0.5, 1.0);
1054 hitCounter->ShiftFillLast(totalHits,
std::sqrt(totalHits), 1);
1064 prescaleHistory = 1;
1065 maxLumiIntegration = 15;
1069 internalDebug =
false;
1071 pi = 3.141592653589793238;
1085 nBinsHistoricalPlot = 80;
1086 nBinsWholeHistory = 3000;
1097 Vx_X->setAxisTitle(
"Entries [#]",2);
1098 Vx_Y->setAxisTitle(
"Primary Vertices Y [cm]",1);
1099 Vx_Y->setAxisTitle(
"Entries [#]",2);
1100 Vx_Z->setAxisTitle(
"Primary Vertices Z [cm]",1);
1101 Vx_Z->setAxisTitle(
"Entries [#]",2);
1103 mXlumi = dbe->
book1D(
"muX vs lumi",
"\\mu_{x} vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1104 mYlumi = dbe->
book1D(
"muY vs lumi",
"\\mu_{y} vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1105 mZlumi = dbe->
book1D(
"muZ vs lumi",
"\\mu_{z} vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1107 mXlumi->setAxisTitle(
"\\mu_{x} [cm]",2);
1108 mXlumi->getTH1()->SetOption(
"E1");
1109 mYlumi->setAxisTitle(
"Lumisection [#]",1);
1110 mYlumi->setAxisTitle(
"\\mu_{y} [cm]",2);
1111 mYlumi->getTH1()->SetOption(
"E1");
1112 mZlumi->setAxisTitle(
"Lumisection [#]",1);
1113 mZlumi->setAxisTitle(
"\\mu_{z} [cm]",2);
1114 mZlumi->getTH1()->SetOption(
"E1");
1116 sXlumi = dbe->
book1D(
"sigmaX vs lumi",
"\\sigma_{x} vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1117 sYlumi = dbe->
book1D(
"sigmaY vs lumi",
"\\sigma_{y} vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1118 sZlumi = dbe->
book1D(
"sigmaZ vs lumi",
"\\sigma_{z} vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1120 sXlumi->setAxisTitle(
"\\sigma_{x} [cm]",2);
1121 sXlumi->getTH1()->SetOption(
"E1");
1122 sYlumi->setAxisTitle(
"Lumisection [#]",1);
1123 sYlumi->setAxisTitle(
"\\sigma_{y} [cm]",2);
1124 sYlumi->getTH1()->SetOption(
"E1");
1125 sZlumi->setAxisTitle(
"Lumisection [#]",1);
1126 sZlumi->setAxisTitle(
"\\sigma_{z} [cm]",2);
1127 sZlumi->getTH1()->SetOption(
"E1");
1129 dxdzlumi = dbe->
book1D(
"dxdz vs lumi",
"dX/dZ vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1130 dydzlumi = dbe->
book1D(
"dydz vs lumi",
"dY/dZ vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1132 dxdzlumi->setAxisTitle(
"dX/dZ [rad]",2);
1133 dxdzlumi->getTH1()->SetOption(
"E1");
1134 dydzlumi->setAxisTitle(
"Lumisection [#]",1);
1135 dydzlumi->setAxisTitle(
"dY/dZ [rad]",2);
1136 dydzlumi->getTH1()->SetOption(
"E1");
1142 Vx_ZX->setAxisTitle(
"Primary Vertices X [cm]",2);
1143 Vx_ZX->setAxisTitle(
"Entries [#]",3);
1144 Vx_ZY->setAxisTitle(
"Primary Vertices Z [cm]",1);
1145 Vx_ZY->setAxisTitle(
"Primary Vertices Y [cm]",2);
1146 Vx_ZY->setAxisTitle(
"Entries [#]",3);
1147 Vx_XY->setAxisTitle(
"Primary Vertices X [cm]",1);
1148 Vx_XY->setAxisTitle(
"Primary Vertices Y [cm]",2);
1149 Vx_XY->setAxisTitle(
"Entries [#]",3);
1151 hitCounter = dbe->
book1D(
"pixelHits vs lumi",
"# Pixel-Hits vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1153 hitCounter->setAxisTitle(
"Pixel-Hits [#]",2);
1154 hitCounter->getTH1()->SetOption(
"E1");
1156 hitCountHistory = dbe->
book1D(
"hist pixelHits vs lumi",
"History: # Pixel-Hits vs. Lumi", nBinsWholeHistory, 0.5, (
double)nBinsWholeHistory+0.5);
1158 hitCountHistory->setAxisTitle(
"Pixel-Hits [#]",2);
1159 hitCountHistory->getTH1()->SetOption(
"E1");
1161 goodVxCounter = dbe->
book1D(
"good vertices vs lumi",
"# Good vertices vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1163 goodVxCounter->setAxisTitle(
"Good vertices [#]",2);
1164 goodVxCounter->getTH1()->SetOption(
"E1");
1166 goodVxCountHistory = dbe->
book1D(
"hist good vx vs lumi",
"History: # Good vx vs. Lumi", nBinsWholeHistory, 0.5, (
double)nBinsWholeHistory+0.5);
1168 goodVxCountHistory->setAxisTitle(
"Good vertices [#]",2);
1169 goodVxCountHistory->getTH1()->SetOption(
"E1");
1171 fitResults = dbe->
book2D(
"fit results",
"Results of Beam Spot Fit", 2, 0., 2., 9, 0., 9.);
1173 fitResults->setBinLabel(9,
"X", 2);
1174 fitResults->setBinLabel(8,
"Y", 2);
1175 fitResults->setBinLabel(7,
"Z", 2);
1176 fitResults->setBinLabel(6,
"\\sigma_{Z}", 2);
1177 fitResults->setBinLabel(5,
"#frac{dX}{dZ}[rad]", 2);
1178 fitResults->setBinLabel(4,
"#frac{dY}{dZ}[rad]", 2);
1179 fitResults->setBinLabel(3,
"\\sigma_{X}", 2);
1180 fitResults->setBinLabel(2,
"\\sigma_{Y}", 2);
1181 fitResults->setBinLabel(1,
"Vertices", 2);
1182 fitResults->setBinLabel(1,
"Value", 1);
1183 fitResults->setBinLabel(2,
"Stat. Error", 1);
1184 fitResults->getTH1()->SetOption(
"text");
1189 reportSummaryMap = dbe->
book2D(
"reportSummaryMap",
"Pixel-Vertices Beam Spot: % Good Fits", 1, 0., 1., 1, 0., 1.);
1190 reportSummaryMap->
Fill(0.5, 0.5, 0.);
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
boost::transform_iterator< IterHelp, const_IdIter > const_iterator
virtual char * formatTime(const time_t &t)
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
#define DEFINE_FWK_MODULE(type)
virtual int MyFit(std::vector< double > *vals)
data_type const * const_iterator
Timestamp const & beginTime() const
virtual unsigned int HitCounter(const edm::Event &iEvent)
virtual void writeToFile(std::vector< double > *vals, edm::TimeValue_t BeginTimeOfFit, edm::TimeValue_t EndTimeOfFit, unsigned int BeginLumiOfFit, unsigned int EndLumiOfFit, int dataType)
MonitorElement * bookFloat(const char *name)
Book float.
LuminosityBlockNumber_t luminosityBlock() const
double Gauss3DFunc(const double *par)
bool considerVxCovariance
virtual void endLuminosityBlock(const edm::LuminosityBlock &lumiBlock, const edm::EventSetup &iSetup)
static ELstring formatTime(const time_t t)
Vx3DHLTAnalyzer(const edm::ParameterSet &)
Timestamp const & endTime() const
virtual void analyze(const edm::Event &, const edm::EventSetup &)
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
LuminosityBlock const & getLuminosityBlock() const
virtual void reset(std::string ResetType)
unsigned long long TimeValue_t
double Covariance[DIM][DIM]
virtual void beginLuminosityBlock(const edm::LuminosityBlock &lumiBlock, const edm::EventSetup &iSetup)
std::vector< std::vector< double > > tmp
static std::atomic< unsigned int > counter
std::vector< VertexType > Vertices
MonitorElement * book2D(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Book 2D histogram.
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
void reset(double vett[256])
TimeValue_t value() const
Power< A, B >::type pow(const A &a, const B &b)
void setCurrentFolder(const std::string &fullpath)