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));
207 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];
208 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];
212 K[0][0] = std::fabs(par[0]);
213 K[1][1] = std::fabs(par[1]);
214 K[2][2] = std::fabs(par[2]);
215 K[0][1] = K[1][0] = par[3];
216 K[1][2] = K[2][1] = par[4]*(std::fabs(par[2])-std::fabs(par[1])) - par[5]*par[3];
217 K[0][2] = K[2][0] = par[5]*(std::fabs(par[2])-std::fabs(par[0])) - par[4]*par[3];
220 det = K[0][0]*(K[1][1]*K[2][2] - K[1][2]*K[1][2]) -
221 K[0][1]*(K[0][1]*K[2][2] - K[0][2]*K[1][2]) +
222 K[0][2]*(K[0][1]*K[1][2] - K[0][2]*K[1][1]);
224 M[0][0] = (K[1][1]*K[2][2] - K[1][2]*K[1][2]) / det;
225 M[1][1] = (K[0][0]*K[2][2] - K[0][2]*K[0][2]) / det;
226 M[2][2] = (K[0][0]*K[1][1] - K[0][1]*K[0][1]) / det;
227 M[0][1] = M[1][0] = (K[0][2]*K[1][2] - K[0][1]*K[2][2]) / det;
228 M[1][2] = M[2][1] = (K[0][2]*K[0][1] - K[1][2]*K[0][0]) / det;
229 M[0][2] = M[2][0] = (K[0][1]*K[1][2] - K[0][2]*K[1][1]) / det;
253 unsigned int nParams = 9;
255 if ((vals !=
NULL) && (vals->size() == nParams*2))
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;
278 vector<double>::const_iterator it = vals->begin();
280 ROOT::Math::Minimizer* Gauss3D = ROOT::Math::Factory::CreateMinimizer(
"Minuit2",
"Migrad");
281 Gauss3D->SetMaxFunctionCalls(1e4);
282 Gauss3D->SetTolerance(1
e-9);
283 if (internalDebug ==
true) Gauss3D->SetPrintLevel(3);
284 else Gauss3D->SetPrintLevel(0);
286 ROOT::Math::Functor _Gauss3DFunc(&
Gauss3DFunc,nParams);
287 Gauss3D->SetFunction(_Gauss3DFunc);
289 if (internalDebug ==
true)
cout <<
"\n@@@ START FITTING @@@" << endl;
293 for (
int i = 0;
i < 3;
i++)
295 deltaMean = (double(
i)-1.)*
std::sqrt((*(it+0))*varFactor);
296 if (internalDebug ==
true)
cout <<
"deltaMean --> " << deltaMean << endl;
300 Gauss3D->SetVariable(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY);
301 Gauss3D->SetVariable(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY);
302 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ);
303 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
304 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
305 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
306 Gauss3D->SetVariable(6,
"mean x", *(it+6)+deltaMean, parDistanceXY);
307 Gauss3D->SetVariable(7,
"mean y", *(it+7), parDistanceXY);
308 Gauss3D->SetVariable(8,
"mean z", *(it+8), parDistanceZ);
311 xPos = Gauss3D->X()[6];
312 yPos = Gauss3D->X()[7];
313 zPos = Gauss3D->X()[8];
320 goodData = Gauss3D->Status();
321 edm = Gauss3D->Edm();
325 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
edm::isNotFinite(Gauss3D->Errors()[
j]) ==
true) { goodData = -1;
break; }
328 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
329 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
331 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
332 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
333 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[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);
346 if (internalDebug ==
true)
348 cout <<
"deltaMean --> " << deltaMean << endl;
349 cout <<
"deltaMean X --> " << (double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor) << endl;
354 Gauss3D->SetVariable(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY);
355 Gauss3D->SetVariable(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY);
356 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ);
357 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
358 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
359 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
360 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY);
361 Gauss3D->SetVariable(7,
"mean y", *(it+7)+deltaMean, parDistanceXY);
362 Gauss3D->SetVariable(8,
"mean z", *(it+8), parDistanceZ);
365 xPos = Gauss3D->X()[6];
366 yPos = Gauss3D->X()[7];
367 zPos = Gauss3D->X()[8];
374 goodData = Gauss3D->Status();
375 edm = Gauss3D->Edm();
379 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
edm::isNotFinite(Gauss3D->Errors()[
j]) ==
true) { goodData = -1;
break; }
382 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
383 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
385 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
386 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
387 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
388 if (det < 0.) { goodData = -1;
if (internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
391 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementY =
i; }
393 if (internalDebug ==
true)
cout <<
"Found bestMovementY --> " << bestMovementY << endl;
397 for (
int i = 0;
i < 3;
i++)
399 deltaMean = (double(
i)-1.)*
std::sqrt(*(it+2));
400 if (internalDebug ==
true)
402 cout <<
"deltaMean --> " << deltaMean << endl;
403 cout <<
"deltaMean X --> " << (double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor) << endl;
404 cout <<
"deltaMean Y --> " << (double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor) << endl;
409 Gauss3D->SetVariable(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY);
410 Gauss3D->SetVariable(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY);
411 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ);
412 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
413 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
414 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
415 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY);
416 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor), parDistanceXY);
417 Gauss3D->SetVariable(8,
"mean z", *(it+8)+deltaMean, parDistanceZ);
420 xPos = Gauss3D->X()[6];
421 yPos = Gauss3D->X()[7];
422 zPos = Gauss3D->X()[8];
429 goodData = Gauss3D->Status();
430 edm = Gauss3D->Edm();
434 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
edm::isNotFinite(Gauss3D->Errors()[
j]) ==
true) { goodData = -1;
break; }
437 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
438 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
440 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
441 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
442 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
443 if (det < 0.) { goodData = -1;
if (internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
446 if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementZ =
i; }
448 if (internalDebug ==
true)
cout <<
"Found bestMovementZ --> " << bestMovementZ << endl;
453 Gauss3D->SetVariable(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY);
454 Gauss3D->SetVariable(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY);
455 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ);
456 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy);
457 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ);
458 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ);
459 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*
std::sqrt((*(it+0))*varFactor), parDistanceXY);
460 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*
std::sqrt((*(it+1))*varFactor), parDistanceXY);
461 Gauss3D->SetVariable(8,
"mean z", *(it+8)+(
double(bestMovementZ)-1.)*
std::sqrt(*(it+2)), parDistanceZ);
464 xPos = Gauss3D->X()[6];
465 yPos = Gauss3D->X()[7];
466 zPos = Gauss3D->X()[8];
473 goodData = Gauss3D->Status();
474 edm = Gauss3D->Edm();
478 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
edm::isNotFinite(Gauss3D->Errors()[
j]) ==
true) { goodData = -1;
break; }
481 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
482 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
484 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
485 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
486 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
487 if (det < 0.) { goodData = -1;
if (internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
491 for (
unsigned int i = 0;
i < trials;
i++)
493 if ((goodData != 0) && (goodData != -2))
497 if (internalDebug ==
true)
cout <<
"FIT WITH DIFFERENT PARAMETER DISTANCES - STEP " <<
i+1 << endl;
499 Gauss3D->SetVariable(0,
"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY * largerDist[
i]);
500 Gauss3D->SetVariable(1,
"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY * largerDist[i]);
501 Gauss3D->SetVariable(2,
"var z ", *(it+2), parDistanceZ*parDistanceZ * largerDist[i]);
502 Gauss3D->SetVariable(3,
"cov xy", *(it+3), parDistanceCxy * largerDist[i]);
503 Gauss3D->SetVariable(4,
"dydz ", *(it+4), parDistanceddZ * largerDist[i]);
504 Gauss3D->SetVariable(5,
"dxdz ", *(it+5), parDistanceddZ * largerDist[i]);
505 Gauss3D->SetVariable(6,
"mean x", *(it+6)+(
double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor), parDistanceXY * largerDist[i]);
506 Gauss3D->SetVariable(7,
"mean y", *(it+7)+(
double(bestMovementY)-1.)*std::sqrt((*(it+1))*varFactor), parDistanceXY * largerDist[i]);
507 Gauss3D->SetVariable(8,
"mean z", *(it+8)+(
double(bestMovementZ)-1.)*std::sqrt(*(it+2)), parDistanceZ * largerDist[i]);
510 xPos = Gauss3D->X()[6];
511 yPos = Gauss3D->X()[7];
512 zPos = Gauss3D->X()[8];
519 goodData = Gauss3D->Status();
520 edm = Gauss3D->Edm();
524 else for (
unsigned int j = 0;
j < nParams;
j++)
if (
edm::isNotFinite(Gauss3D->Errors()[
j]) ==
true) { goodData = -1;
break; }
527 covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3];
528 covxz = Gauss3D->X()[5]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[0])) - Gauss3D->X()[4]*Gauss3D->X()[3];
530 det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) -
531 Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) +
532 covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1]));
533 if (det < 0.) { goodData = -1;
if (internalDebug ==
true)
cout <<
"Negative determinant !" << endl; }
539 for (
unsigned int i = 0;
i < nParams;
i++)
541 vals->operator[](
i) = Gauss3D->X()[
i];
542 vals->operator[](
i+nParams) = Gauss3D->Errors()[
i];
555 if (ResetType.compare(
"scratch") == 0)
582 hitCountHistory->Reset();
583 goodVxCounter->Reset();
584 goodVxCountHistory->Reset();
588 reportSummaryMap->Fill(0.5, 0.5, 0.);
593 lumiCounterHisto = 0;
600 else if (ResetType.compare(
"whole") == 0)
613 lumiCounterHisto = 0;
620 else if (ResetType.compare(
"partial") == 0)
635 else if (ResetType.compare(
"nohisto") == 0)
640 lumiCounterHisto = 0;
647 else if (ResetType.compare(
"hitCounter") == 0)
655 unsigned int BeginLumiOfFit,
656 unsigned int EndLumiOfFit,
659 stringstream BufferString;
660 BufferString.precision(5);
664 if ((
outputFile.is_open() ==
true) && (vals !=
NULL) && (vals->size() == 8*2))
666 vector<double>::const_iterator it = vals->begin();
669 outputFile <<
"BeginTimeOfFit " <<
formatTime(beginTimeOfFit >> 32) <<
" " << (beginTimeOfFit >> 32) << endl;
670 outputFile <<
"EndTimeOfFit " <<
formatTime(endTimeOfFit >> 32) <<
" " << (endTimeOfFit >> 32) << endl;
671 outputFile <<
"LumiRange " << beginLumiOfFit <<
" - " << endLumiOfFit << endl;
677 BufferString << *(it+0);
678 outputFile <<
"X0 " << BufferString.str().c_str() << endl;
679 BufferString.str(
"");
681 BufferString << *(it+1);
682 outputFile <<
"Y0 " << BufferString.str().c_str() << endl;
683 BufferString.str(
"");
685 BufferString << *(it+2);
686 outputFile <<
"Z0 " << BufferString.str().c_str() << endl;
687 BufferString.str(
"");
689 BufferString << *(it+3);
690 outputFile <<
"sigmaZ0 " << BufferString.str().c_str() << endl;
691 BufferString.str(
"");
693 BufferString << *(it+4);
694 outputFile <<
"dxdz " << BufferString.str().c_str() << endl;
695 BufferString.str(
"");
697 BufferString << *(it+5);
698 outputFile <<
"dydz " << BufferString.str().c_str() << endl;
699 BufferString.str(
"");
701 BufferString << *(it+6);
702 outputFile <<
"BeamWidthX " << BufferString.str().c_str() << endl;
703 BufferString.str(
"");
705 BufferString << *(it+7);
706 outputFile <<
"BeamWidthY " << BufferString.str().c_str() << endl;
707 BufferString.str(
"");
709 outputFile <<
"Cov(0,j) " << *(it+8) <<
" 0.0 0.0 0.0 0.0 0.0 0.0" << endl;
710 outputFile <<
"Cov(1,j) 0.0 " << *(it+9) <<
" 0.0 0.0 0.0 0.0 0.0" << endl;
711 outputFile <<
"Cov(2,j) 0.0 0.0 " << *(it+10) <<
" 0.0 0.0 0.0 0.0" << endl;
712 outputFile <<
"Cov(3,j) 0.0 0.0 0.0 " << *(it+11) <<
" 0.0 0.0 0.0" << endl;
713 outputFile <<
"Cov(4,j) 0.0 0.0 0.0 0.0 " << *(it+12) <<
" 0.0 0.0" << endl;
714 outputFile <<
"Cov(5,j) 0.0 0.0 0.0 0.0 0.0 " << *(it+13) <<
" 0.0" << endl;
715 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;
723 if ((
debugMode ==
true) && (outputDebugFile.is_open() ==
true) && (vals !=
NULL) && (vals->size() == 8*2))
725 vector<double>::const_iterator it = vals->begin();
727 outputDebugFile <<
"Runnumber " <<
runNumber << endl;
728 outputDebugFile <<
"BeginTimeOfFit " <<
formatTime(beginTimeOfFit >> 32) <<
" " << (beginTimeOfFit >> 32) << endl;
729 outputDebugFile <<
"EndTimeOfFit " <<
formatTime(endTimeOfFit >> 32) <<
" " << (endTimeOfFit >> 32) << endl;
730 outputDebugFile <<
"LumiRange " << beginLumiOfFit <<
" - " << endLumiOfFit << endl;
731 outputDebugFile <<
"Type " << dataType << endl;
736 BufferString << *(it+0);
737 outputDebugFile <<
"X0 " << BufferString.str().c_str() << endl;
738 BufferString.str(
"");
740 BufferString << *(it+1);
741 outputDebugFile <<
"Y0 " << BufferString.str().c_str() << endl;
742 BufferString.str(
"");
744 BufferString << *(it+2);
745 outputDebugFile <<
"Z0 " << BufferString.str().c_str() << endl;
746 BufferString.str(
"");
748 BufferString << *(it+3);
749 outputDebugFile <<
"sigmaZ0 " << BufferString.str().c_str() << endl;
750 BufferString.str(
"");
752 BufferString << *(it+4);
753 outputDebugFile <<
"dxdz " << BufferString.str().c_str() << endl;
754 BufferString.str(
"");
756 BufferString << *(it+5);
757 outputDebugFile <<
"dydz " << BufferString.str().c_str() << endl;
758 BufferString.str(
"");
760 BufferString << *(it+6);
761 outputDebugFile <<
"BeamWidthX " << BufferString.str().c_str() << endl;
762 BufferString.str(
"");
764 BufferString << *(it+7);
765 outputDebugFile <<
"BeamWidthY " << BufferString.str().c_str() << endl;
766 BufferString.str(
"");
768 outputDebugFile <<
"Cov(0,j) " << *(it+8) <<
" 0.0 0.0 0.0 0.0 0.0 0.0" << endl;
769 outputDebugFile <<
"Cov(1,j) 0.0 " << *(it+9) <<
" 0.0 0.0 0.0 0.0 0.0" << endl;
770 outputDebugFile <<
"Cov(2,j) 0.0 0.0 " << *(it+10) <<
" 0.0 0.0 0.0 0.0" << endl;
771 outputDebugFile <<
"Cov(3,j) 0.0 0.0 0.0 " << *(it+11) <<
" 0.0 0.0 0.0" << endl;
772 outputDebugFile <<
"Cov(4,j) 0.0 0.0 0.0 0.0 " << *(it+12) <<
" 0.0 0.0" << endl;
773 outputDebugFile <<
"Cov(5,j) 0.0 0.0 0.0 0.0 0.0 " << *(it+13) <<
" 0.0" << endl;
774 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;
776 outputDebugFile <<
"EmittanceX 0.0" << endl;
777 outputDebugFile <<
"EmittanceY 0.0" << endl;
778 outputDebugFile <<
"BetaStar 0.0" << endl;
786 if ((lumiCounter == 0) && (lumiBlock.
luminosityBlock() > lastLumiOfFit))
793 else if ((lumiCounter != 0) && (lumiBlock.
luminosityBlock() >= (beginLumiOfFit+lumiCounter))) { lumiCounter++; lumiCounterHisto++; }
800 stringstream histTitle;
802 unsigned int nParams = 9;
808 lastLumiOfFit = endLumiOfFit;
813 if (lastLumiOfFit % prescaleHistory == 0)
815 hitCountHistory->getTH1()->SetBinContent(lastLumiOfFit, (
double)totalHits);
816 hitCountHistory->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt((
double)totalHits));
821 vector<double> fitResults;
823 fitResults.push_back(Vx_X->getTH1()->GetRMS()*Vx_X->getTH1()->GetRMS());
824 fitResults.push_back(Vx_Y->getTH1()->GetRMS()*Vx_Y->getTH1()->GetRMS());
825 fitResults.push_back(Vx_Z->getTH1()->GetRMS()*Vx_Z->getTH1()->GetRMS());
826 fitResults.push_back(0.0);
827 fitResults.push_back(0.0);
828 fitResults.push_back(0.0);
829 fitResults.push_back(Vx_X->getTH1()->GetMean());
830 fitResults.push_back(Vx_Y->getTH1()->GetMean());
831 fitResults.push_back(Vx_Z->getTH1()->GetMean());
832 for (
unsigned int i = 0;
i < nParams;
i++) fitResults.push_back(0.0);
834 goodData = MyFit(&fitResults);
836 if (internalDebug ==
true)
838 cout <<
"goodData --> " << goodData << endl;
840 cout <<
"var x --> " << fitResults[0] <<
" +/- " << fitResults[0+nParams] << endl;
841 cout <<
"var y --> " << fitResults[1] <<
" +/- " << fitResults[1+nParams] << endl;
842 cout <<
"var z --> " << fitResults[2] <<
" +/- " << fitResults[2+nParams] << endl;
843 cout <<
"cov xy --> " << fitResults[3] <<
" +/- " << fitResults[3+nParams] << endl;
844 cout <<
"dydz --> " << fitResults[4] <<
" +/- " << fitResults[4+nParams] << endl;
845 cout <<
"dxdz --> " << fitResults[5] <<
" +/- " << fitResults[5+nParams] << endl;
846 cout <<
"mean x --> " << fitResults[6] <<
" +/- " << fitResults[6+nParams] << endl;
847 cout <<
"mean y --> " << fitResults[7] <<
" +/- " << fitResults[7+nParams] << endl;
848 cout <<
"mean z --> " << fitResults[8] <<
" +/- " << fitResults[8+nParams] << endl;
853 vals.push_back(fitResults[6]);
854 vals.push_back(fitResults[7]);
855 vals.push_back(fitResults[8]);
856 vals.push_back(
std::sqrt(std::fabs(fitResults[2])));
857 vals.push_back(fitResults[5]);
858 vals.push_back(fitResults[4]);
859 vals.push_back(
std::sqrt(std::fabs(fitResults[0])));
860 vals.push_back(
std::sqrt(std::fabs(fitResults[1])));
862 vals.push_back(
std::pow(fitResults[6+nParams],2.));
863 vals.push_back(
std::pow(fitResults[7+nParams],2.));
864 vals.push_back(
std::pow(fitResults[8+nParams],2.));
865 vals.push_back(
std::pow(std::fabs(fitResults[2+nParams]) / (2.*
std::sqrt(std::fabs(fitResults[2]))),2.));
866 vals.push_back(
std::pow(fitResults[5+nParams],2.));
867 vals.push_back(
std::pow(fitResults[4+nParams],2.));
868 vals.push_back(
std::pow(std::fabs(fitResults[0+nParams]) / (2.*
std::sqrt(std::fabs(fitResults[0]))),2.));
869 vals.push_back(
std::pow(std::fabs(fitResults[1+nParams]) / (2.*
std::sqrt(std::fabs(fitResults[1]))),2.));
871 else for (
unsigned int i = 0;
i < 8*2;
i++) vals.push_back(0.0);
877 counterVx = Vx_X->getTH1F()->GetEntries();
883 vals.push_back(Vx_X->getTH1F()->GetMean());
884 vals.push_back(Vx_Y->getTH1F()->GetMean());
885 vals.push_back(Vx_Z->getTH1F()->GetMean());
886 vals.push_back(Vx_Z->getTH1F()->GetRMS());
889 vals.push_back(Vx_X->getTH1F()->GetRMS());
890 vals.push_back(Vx_Y->getTH1F()->GetRMS());
892 vals.push_back(
std::pow(Vx_X->getTH1F()->GetMeanError(),2.));
893 vals.push_back(
std::pow(Vx_Y->getTH1F()->GetMeanError(),2.));
894 vals.push_back(
std::pow(Vx_Z->getTH1F()->GetMeanError(),2.));
895 vals.push_back(
std::pow(Vx_Z->getTH1F()->GetRMSError(),2.));
898 vals.push_back(
std::pow(Vx_X->getTH1F()->GetRMSError(),2.));
899 vals.push_back(
std::pow(Vx_Y->getTH1F()->GetRMSError(),2.));
904 for (
unsigned int i = 0;
i < 8*2;
i++) vals.push_back(0.0);
934 writeToFile(&vals, beginTimeOfFit, endTimeOfFit, beginLumiOfFit, endLumiOfFit, 3);
935 if ((internalDebug ==
true) && (outputDebugFile.is_open() ==
true)) outputDebugFile <<
"Used vertices: " <<
counterVx << endl;
939 histTitle <<
"Fitted Beam Spot [cm] (Lumi start: " << beginLumiOfFit <<
" - Lumi end: " << endLumiOfFit <<
")";
940 if (lumiCounterHisto >= maxLumiIntegration)
reset(
"whole");
941 else reset(
"partial");
945 writeToFile(&vals, beginTimeOfFit, endTimeOfFit, beginLumiOfFit, endLumiOfFit, -1);
946 if ((internalDebug ==
true) && (outputDebugFile.is_open() ==
true)) outputDebugFile <<
"Used vertices: " <<
counterVx << endl;
950 histTitle <<
"Fitted Beam Spot [cm] (not enough statistics)";
951 if (lumiCounter >= maxLumiIntegration)
reset(
"whole");
952 else reset(
"hitCounter");
956 histTitle <<
"Fitted Beam Spot [cm] (problems)";
957 if (lumiCounterHisto >= maxLumiIntegration)
reset(
"whole");
958 else reset(
"partial");
964 reportSummary->Fill(numberFits != 0 ? (
double)numberGoodFits/(
double)numberFits : 0.0);
965 reportSummaryMap->Fill(0.5, 0.5, numberFits != 0 ? (
double)numberGoodFits/(
double)numberFits : 0.0);
967 fitResults->setAxisTitle(histTitle.str().c_str(), 1);
969 fitResults->setBinContent(1, 9, vals[0]);
970 fitResults->setBinContent(1, 8, vals[1]);
971 fitResults->setBinContent(1, 7, vals[2]);
972 fitResults->setBinContent(1, 6, vals[3]);
973 fitResults->setBinContent(1, 5, vals[4]);
974 fitResults->setBinContent(1, 4, vals[5]);
975 fitResults->setBinContent(1, 3, vals[6]);
976 fitResults->setBinContent(1, 2, vals[7]);
977 fitResults->setBinContent(1, 1,
counterVx);
979 fitResults->setBinContent(2, 9,
std::sqrt(vals[8]));
980 fitResults->setBinContent(2, 8,
std::sqrt(vals[9]));
981 fitResults->setBinContent(2, 7,
std::sqrt(vals[10]));
982 fitResults->setBinContent(2, 6,
std::sqrt(vals[11]));
983 fitResults->setBinContent(2, 5,
std::sqrt(vals[12]));
984 fitResults->setBinContent(2, 4,
std::sqrt(vals[13]));
985 fitResults->setBinContent(2, 3,
std::sqrt(vals[14]));
986 fitResults->setBinContent(2, 2,
std::sqrt(vals[15]));
990 TF1* myLinFit =
new TF1(
"myLinFit",
"[0] + [1]*x", mXlumi->getTH1()->GetXaxis()->GetXmin(), mXlumi->getTH1()->GetXaxis()->GetXmax());
991 myLinFit->SetLineColor(2);
992 myLinFit->SetLineWidth(2);
993 myLinFit->SetParName(0,
"Intercept");
994 myLinFit->SetParName(1,
"Slope");
997 myLinFit->SetParameter(0, mXlumi->getTH1()->GetMean(2));
998 myLinFit->SetParameter(1, 0.0);
999 mXlumi->getTH1()->Fit(myLinFit,
"QR");
1002 myLinFit->SetParameter(0, mYlumi->getTH1()->GetMean(2));
1003 myLinFit->SetParameter(1, 0.0);
1004 mYlumi->getTH1()->Fit(myLinFit,
"QR");
1007 myLinFit->SetParameter(0, mZlumi->getTH1()->GetMean(2));
1008 myLinFit->SetParameter(1, 0.0);
1009 mZlumi->getTH1()->Fit(myLinFit,
"QR");
1012 myLinFit->SetParameter(0, sXlumi->getTH1()->GetMean(2));
1013 myLinFit->SetParameter(1, 0.0);
1014 sXlumi->getTH1()->Fit(myLinFit,
"QR");
1017 myLinFit->SetParameter(0, sYlumi->getTH1()->GetMean(2));
1018 myLinFit->SetParameter(1, 0.0);
1019 sYlumi->getTH1()->Fit(myLinFit,
"QR");
1022 myLinFit->SetParameter(0, sZlumi->getTH1()->GetMean(2));
1023 myLinFit->SetParameter(1, 0.0);
1024 sZlumi->getTH1()->Fit(myLinFit,
"QR");
1027 myLinFit->SetParameter(0, dxdzlumi->getTH1()->GetMean(2));
1028 myLinFit->SetParameter(1, 0.0);
1029 dxdzlumi->getTH1()->Fit(myLinFit,
"QR");
1032 myLinFit->SetParameter(0, dydzlumi->getTH1()->GetMean(2));
1033 myLinFit->SetParameter(1, 0.0);
1034 dydzlumi->getTH1()->Fit(myLinFit,
"QR");
1037 myLinFit->SetParameter(0, goodVxCounter->getTH1()->GetMean(2));
1038 myLinFit->SetParameter(1, 0.0);
1039 goodVxCounter->getTH1()->Fit(myLinFit,
"QR");
1041 if (lastLumiOfFit % prescaleHistory == 0)
1043 goodVxCountHistory->getTH1()->SetBinContent(lastLumiOfFit, (
double)counterVx);
1044 goodVxCountHistory->getTH1()->SetBinError(lastLumiOfFit,
std::sqrt((
double)counterVx));
1053 histTitle <<
"Fitted Beam Spot [cm] (no ongoing fits)";
1054 fitResults->setAxisTitle(histTitle.str().c_str(), 1);
1055 reportSummaryMap->Fill(0.5, 0.5, 1.0);
1056 hitCounter->ShiftFillLast(totalHits,
std::sqrt(totalHits), 1);
1066 prescaleHistory = 1;
1067 maxLumiIntegration = 15;
1071 internalDebug =
false;
1073 pi = 3.141592653589793238;
1087 nBinsHistoricalPlot = 80;
1088 nBinsWholeHistory = 3000;
1093 dbe->setCurrentFolder(
"BeamPixel");
1095 Vx_X = dbe->book1D(
"vertex x",
"Primary Vertex X Coordinate Distribution",
int(rint(
xRange/
xStep)), -
xRange/2.,
xRange/2.);
1096 Vx_Y = dbe->book1D(
"vertex y",
"Primary Vertex Y Coordinate Distribution",
int(rint(
yRange/
yStep)), -
yRange/2.,
yRange/2.);
1097 Vx_Z = dbe->book1D(
"vertex z",
"Primary Vertex Z Coordinate Distribution",
int(rint(
zRange/
zStep)), -
zRange/2.,
zRange/2.);
1098 Vx_X->setAxisTitle(
"Primary Vertices X [cm]",1);
1099 Vx_X->setAxisTitle(
"Entries [#]",2);
1100 Vx_Y->setAxisTitle(
"Primary Vertices Y [cm]",1);
1101 Vx_Y->setAxisTitle(
"Entries [#]",2);
1102 Vx_Z->setAxisTitle(
"Primary Vertices Z [cm]",1);
1103 Vx_Z->setAxisTitle(
"Entries [#]",2);
1105 mXlumi = dbe->book1D(
"muX vs lumi",
"\\mu_{x} vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1106 mYlumi = dbe->book1D(
"muY vs lumi",
"\\mu_{y} vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1107 mZlumi = dbe->book1D(
"muZ vs lumi",
"\\mu_{z} vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1108 mXlumi->setAxisTitle(
"Lumisection [#]",1);
1109 mXlumi->setAxisTitle(
"\\mu_{x} [cm]",2);
1110 mXlumi->getTH1()->SetOption(
"E1");
1111 mYlumi->setAxisTitle(
"Lumisection [#]",1);
1112 mYlumi->setAxisTitle(
"\\mu_{y} [cm]",2);
1113 mYlumi->getTH1()->SetOption(
"E1");
1114 mZlumi->setAxisTitle(
"Lumisection [#]",1);
1115 mZlumi->setAxisTitle(
"\\mu_{z} [cm]",2);
1116 mZlumi->getTH1()->SetOption(
"E1");
1118 sXlumi = dbe->book1D(
"sigmaX vs lumi",
"\\sigma_{x} vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1119 sYlumi = dbe->book1D(
"sigmaY vs lumi",
"\\sigma_{y} vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1120 sZlumi = dbe->book1D(
"sigmaZ vs lumi",
"\\sigma_{z} vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1121 sXlumi->setAxisTitle(
"Lumisection [#]",1);
1122 sXlumi->setAxisTitle(
"\\sigma_{x} [cm]",2);
1123 sXlumi->getTH1()->SetOption(
"E1");
1124 sYlumi->setAxisTitle(
"Lumisection [#]",1);
1125 sYlumi->setAxisTitle(
"\\sigma_{y} [cm]",2);
1126 sYlumi->getTH1()->SetOption(
"E1");
1127 sZlumi->setAxisTitle(
"Lumisection [#]",1);
1128 sZlumi->setAxisTitle(
"\\sigma_{z} [cm]",2);
1129 sZlumi->getTH1()->SetOption(
"E1");
1131 dxdzlumi = dbe->book1D(
"dxdz vs lumi",
"dX/dZ vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1132 dydzlumi = dbe->book1D(
"dydz vs lumi",
"dY/dZ vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1133 dxdzlumi->setAxisTitle(
"Lumisection [#]",1);
1134 dxdzlumi->setAxisTitle(
"dX/dZ [rad]",2);
1135 dxdzlumi->getTH1()->SetOption(
"E1");
1136 dydzlumi->setAxisTitle(
"Lumisection [#]",1);
1137 dydzlumi->setAxisTitle(
"dY/dZ [rad]",2);
1138 dydzlumi->getTH1()->SetOption(
"E1");
1143 Vx_ZX->setAxisTitle(
"Primary Vertices Z [cm]",1);
1144 Vx_ZX->setAxisTitle(
"Primary Vertices X [cm]",2);
1145 Vx_ZX->setAxisTitle(
"Entries [#]",3);
1146 Vx_ZY->setAxisTitle(
"Primary Vertices Z [cm]",1);
1147 Vx_ZY->setAxisTitle(
"Primary Vertices Y [cm]",2);
1148 Vx_ZY->setAxisTitle(
"Entries [#]",3);
1149 Vx_XY->setAxisTitle(
"Primary Vertices X [cm]",1);
1150 Vx_XY->setAxisTitle(
"Primary Vertices Y [cm]",2);
1151 Vx_XY->setAxisTitle(
"Entries [#]",3);
1153 hitCounter = dbe->book1D(
"pixelHits vs lumi",
"# Pixel-Hits vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1154 hitCounter->setAxisTitle(
"Lumisection [#]",1);
1155 hitCounter->setAxisTitle(
"Pixel-Hits [#]",2);
1156 hitCounter->getTH1()->SetOption(
"E1");
1158 hitCountHistory = dbe->book1D(
"hist pixelHits vs lumi",
"History: # Pixel-Hits vs. Lumi", nBinsWholeHistory, 0.5, (
double)nBinsWholeHistory+0.5);
1159 hitCountHistory->setAxisTitle(
"Lumisection [#]",1);
1160 hitCountHistory->setAxisTitle(
"Pixel-Hits [#]",2);
1161 hitCountHistory->getTH1()->SetOption(
"E1");
1163 goodVxCounter = dbe->book1D(
"good vertices vs lumi",
"# Good vertices vs. Lumisection", nBinsHistoricalPlot, 0.5, (
double)nBinsHistoricalPlot+0.5);
1164 goodVxCounter->setAxisTitle(
"Lumisection [#]",1);
1165 goodVxCounter->setAxisTitle(
"Good vertices [#]",2);
1166 goodVxCounter->getTH1()->SetOption(
"E1");
1168 goodVxCountHistory = dbe->book1D(
"hist good vx vs lumi",
"History: # Good vx vs. Lumi", nBinsWholeHistory, 0.5, (
double)nBinsWholeHistory+0.5);
1169 goodVxCountHistory->setAxisTitle(
"Lumisection [#]",1);
1170 goodVxCountHistory->setAxisTitle(
"Good vertices [#]",2);
1171 goodVxCountHistory->getTH1()->SetOption(
"E1");
1173 fitResults = dbe->book2D(
"fit results",
"Results of Beam Spot Fit", 2, 0., 2., 9, 0., 9.);
1174 fitResults->setAxisTitle(
"Fitted Beam Spot [cm]", 1);
1175 fitResults->setBinLabel(9,
"X", 2);
1176 fitResults->setBinLabel(8,
"Y", 2);
1177 fitResults->setBinLabel(7,
"Z", 2);
1178 fitResults->setBinLabel(6,
"\\sigma_{Z}", 2);
1179 fitResults->setBinLabel(5,
"#frac{dX}{dZ}[rad]", 2);
1180 fitResults->setBinLabel(4,
"#frac{dY}{dZ}[rad]", 2);
1181 fitResults->setBinLabel(3,
"\\sigma_{X}", 2);
1182 fitResults->setBinLabel(2,
"\\sigma_{Y}", 2);
1183 fitResults->setBinLabel(1,
"Vertices", 2);
1184 fitResults->setBinLabel(1,
"Value", 1);
1185 fitResults->setBinLabel(2,
"Stat. Error", 1);
1186 fitResults->getTH1()->SetOption(
"text");
1188 dbe->setCurrentFolder(
"BeamPixel/EventInfo");
1191 reportSummaryMap = dbe->book2D(
"reportSummaryMap",
"Pixel-Vertices Beam Spot: % Good Fits", 1, 0., 1., 1, 0., 1.);
1192 reportSummaryMap->Fill(0.5, 0.5, 0.);
1193 dbe->setCurrentFolder(
"BeamPixel/EventInfo/reportSummaryContents");
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
boost::transform_iterator< IterHelp, const_IdIter > const_iterator
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)
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)
virtual std::string formatTime(const time_t &t)
std::vector< std::vector< double > > tmp
static std::atomic< unsigned int > counter
std::vector< VertexType > Vertices
void reset(double vett[256])
TimeValue_t value() const
Power< A, B >::type pow(const A &a, const B &b)