17 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
31 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
35 #define LOGERROR(x) LogError(x)
36 #define LOGWARNING(x) LogWarning(x)
37 #define LOGINFO(x) LogInfo(x)
43 #define LOGERROR(x) std::cout << x << ": "
44 #define LOGINFO(x) std::cout << x << ": "
45 #define LOGWARNING(x) std::cout << x << ": "
46 #define ENDL std::endl
60 float costrk[3] = {0, 0, 0};
64 const int code_version = {1};
68 std::ostringstream tout;
72 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
73 tout <<
dir <<
"generror_summary_zp" << std::setw(4) << std::setfill(
'0') << std::right << filenum <<
".out"
77 tempfile = (
file.fullPath()).c_str();
79 tout <<
"generror_summary_zp" << std::setw(4) << std::setfill(
'0') << std::right << filenum <<
".out" << std::ends;
81 tempfile = tempf.c_str();
88 if (in_file.is_open()) {
95 for (
i = 0; (
c = in_file.get()) !=
'\n'; ++
i) {
104 LOGINFO(
"SiPixelGenError") <<
"Loading Pixel GenError File - " << theCurrentTemp.
head.
title <<
ENDL;
117 if (in_file.fail()) {
118 LOGERROR(
"SiPixelGenError") <<
"Error reading file, no GenError load" <<
ENDL;
122 LOGINFO(
"SiPixelGenError") <<
"GenError ID = " << theCurrentTemp.
head.
ID <<
", GenError Version "
124 <<
", NTy = " << theCurrentTemp.
head.
NTy <<
", NTyx = " << theCurrentTemp.
head.
NTyx
125 <<
", NTxx = " << theCurrentTemp.
head.
NTxx <<
", Dtype = " << theCurrentTemp.
head.
Dtype
126 <<
", Bias voltage " << theCurrentTemp.
head.
Vbias <<
", temperature "
128 <<
", Q-scaling factor " << theCurrentTemp.
head.
qscale <<
", 1/2 multi dcol threshold "
129 << theCurrentTemp.
head.
s50 <<
", 1/2 single dcol threshold " << theCurrentTemp.
head.
ss50
130 <<
", y Lorentz Width " << theCurrentTemp.
head.
lorywidth <<
", y Lorentz Bias "
133 <<
", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.
head.
fbin[0] <<
", "
135 <<
", pixel x-size " << theCurrentTemp.
head.
xsize <<
", y-size "
139 LOGERROR(
"SiPixelGenError") <<
"code expects version " << code_version <<
", no GenError load" <<
ENDL;
143 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST
151 theCurrentTemp.
enty.resize(boost::extents[theCurrentTemp.
head.
NTy]);
159 for (
i = 0;
i < theCurrentTemp.
head.
NTy; ++
i) {
160 in_file >> theCurrentTemp.
enty[
i].
runnum >> costrk[0] >> costrk[1] >> costrk[2];
162 if (in_file.fail()) {
163 LOGERROR(
"SiPixelGenError") <<
"Error reading file 1, no GenError load, run # " << theCurrentTemp.
enty[
i].
runnum
177 if (in_file.fail()) {
178 LOGERROR(
"SiPixelGenError") <<
"Error reading file 2, no GenError load, run # " << theCurrentTemp.
enty[
i].
runnum
186 if (in_file.fail()) {
187 LOGERROR(
"SiPixelGenError") <<
"Error reading file 3, no GenError load, run # " << theCurrentTemp.
enty[
i].
runnum
193 LOGERROR(
"SiPixelGenError") <<
"Error in GenError ID " << theCurrentTemp.
head.
ID
194 <<
" qmin = " << theCurrentTemp.
enty[
i].
qmin <<
", run # "
199 for (
j = 0;
j < 4; ++
j) {
203 if (in_file.fail()) {
204 LOGERROR(
"SiPixelGenError") <<
"Error reading file 14a, no GenError load, run # "
215 in_file >> theCurrentTemp.
entx[
k][
i].
runnum >> costrk[0] >> costrk[1] >> costrk[2];
217 if (in_file.fail()) {
218 LOGERROR(
"SiPixelGenError") <<
"Error reading file 17, no GenError load, run # "
233 if (in_file.fail()) {
234 LOGERROR(
"SiPixelGenError") <<
"Error reading file 18, no GenError load, run # "
244 if (in_file.fail()) {
245 LOGERROR(
"SiPixelGenError") <<
"Error reading file 19, no GenError load, run # "
250 for (
j = 0;
j < 4; ++
j) {
254 if (in_file.fail()) {
255 LOGERROR(
"SiPixelGenError") <<
"Error reading file 30a, no GenError load, run # "
267 pixelTemp.push_back(theCurrentTemp);
276 LOGERROR(
"SiPixelGenError") <<
"Error opening File" << tempfile <<
ENDL;
282 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
294 float costrk[3] = {0, 0, 0};
296 const int code_version = {1};
303 auto tmpPtr = std::make_unique<SiPixelGenErrorStore>();
304 auto& theCurrentTemp = *tmpPtr;
307 for (
int m = 0;
m <
db.numOfTempl(); ++
m) {
311 for (
i = 0;
i < 20; ++
i) {
313 theCurrentTemp.head.title[4 *
i] =
temp.c[0];
314 theCurrentTemp.head.title[4 *
i + 1] =
temp.c[1];
315 theCurrentTemp.head.title[4 *
i + 2] =
temp.c[2];
316 theCurrentTemp.head.title[4 *
i + 3] =
temp.c[3];
317 db.incrementIndex(1);
320 theCurrentTemp.head.title[79] =
'\0';
321 LOGINFO(
"SiPixelGenError") <<
"Loading Pixel GenError File - " << theCurrentTemp.head.title <<
ENDL;
325 db >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >>
326 theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx >> theCurrentTemp.head.Dtype >>
327 theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >>
328 theCurrentTemp.head.qscale >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >>
329 theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >>
330 theCurrentTemp.head.zsize >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >>
331 theCurrentTemp.head.lorxbias >> theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >>
332 theCurrentTemp.head.fbin[2];
335 LOGERROR(
"SiPixelGenError") <<
"Error reading file, no GenError load" <<
ENDL;
339 LOGINFO(
"SiPixelGenError") <<
"GenError ID = " << theCurrentTemp.head.ID <<
", GenError Version "
340 << theCurrentTemp.head.templ_version <<
", Bfield = " << theCurrentTemp.head.Bfield
341 <<
", NTy = " << theCurrentTemp.head.NTy <<
", NTyx = " << theCurrentTemp.head.NTyx
342 <<
", NTxx = " << theCurrentTemp.head.NTxx <<
", Dtype = " << theCurrentTemp.head.Dtype
343 <<
", Bias voltage " << theCurrentTemp.head.Vbias <<
", temperature "
344 << theCurrentTemp.head.temperature <<
", fluence " << theCurrentTemp.head.fluence
345 <<
", Q-scaling factor " << theCurrentTemp.head.qscale <<
", 1/2 multi dcol threshold "
346 << theCurrentTemp.head.s50 <<
", 1/2 single dcol threshold " << theCurrentTemp.head.ss50
347 <<
", y Lorentz Width " << theCurrentTemp.head.lorywidth <<
", y Lorentz Bias "
348 << theCurrentTemp.head.lorybias <<
", x Lorentz width " << theCurrentTemp.head.lorxwidth
349 <<
", x Lorentz Bias " << theCurrentTemp.head.lorxbias
350 <<
", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] <<
", "
351 << theCurrentTemp.head.fbin[1] <<
", " << theCurrentTemp.head.fbin[2]
352 <<
", pixel x-size " << theCurrentTemp.head.xsize <<
", y-size "
353 << theCurrentTemp.head.ysize <<
", zsize " << theCurrentTemp.head.zsize <<
ENDL;
355 LOGINFO(
"SiPixelGenError") <<
"Loading Pixel GenError - " << theCurrentTemp.head.title <<
" version "
356 << theCurrentTemp.head.templ_version <<
" code v." << code_version <<
ENDL;
357 if (theCurrentTemp.head.templ_version < code_version) {
358 LOGERROR(
"SiPixelGenError") <<
"code expects version " << code_version <<
", no GenError load" <<
ENDL;
362 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST
367 theCurrentTemp.cotbetaY =
new float[theCurrentTemp.head.NTy];
368 theCurrentTemp.cotbetaX =
new float[theCurrentTemp.head.NTyx];
369 theCurrentTemp.cotalphaX =
new float[theCurrentTemp.head.NTxx];
371 theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
373 theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
379 for (
i = 0;
i < theCurrentTemp.head.NTy; ++
i) {
380 db >> theCurrentTemp.enty[
i].runnum >> costrk[0] >> costrk[1] >> costrk[2];
383 LOGERROR(
"SiPixelGenError") <<
"Error reading file 1, no GenError load, run # " << theCurrentTemp.enty[
i].runnum
390 theCurrentTemp.enty[
i].cotalpha = costrk[0] / costrk[2];
392 theCurrentTemp.enty[
i].cotbeta = costrk[1] / costrk[2];
394 db >> theCurrentTemp.enty[
i].qavg >> theCurrentTemp.enty[
i].pixmax >> theCurrentTemp.enty[
i].dyone >>
395 theCurrentTemp.enty[
i].syone >> theCurrentTemp.enty[
i].dxone >> theCurrentTemp.enty[
i].sxone;
398 LOGERROR(
"SiPixelGenError") <<
"Error reading file 2, no GenError load, run # " << theCurrentTemp.enty[
i].runnum
403 db >> theCurrentTemp.enty[
i].dytwo >> theCurrentTemp.enty[
i].sytwo >> theCurrentTemp.enty[
i].dxtwo >>
404 theCurrentTemp.enty[
i].sxtwo >> theCurrentTemp.enty[
i].qmin >> theCurrentTemp.enty[
i].qmin2;
406 for (
j = 0;
j < 4; ++
j) {
407 db >> theCurrentTemp.enty[
i].yavggen[
j] >> theCurrentTemp.enty[
i].yrmsgen[
j] >>
408 theCurrentTemp.enty[
i].xavggen[
j] >> theCurrentTemp.enty[
i].xrmsgen[
j];
411 LOGERROR(
"SiPixelGenError") <<
"Error reading file 14a, no GenError load, run # "
412 << theCurrentTemp.enty[
i].runnum <<
ENDL;
420 for (
k = 0;
k < theCurrentTemp.head.NTyx; ++
k) {
421 for (
i = 0;
i < theCurrentTemp.head.NTxx; ++
i) {
422 db >> theCurrentTemp.entx[
k][
i].runnum >> costrk[0] >> costrk[1] >> costrk[2];
425 LOGERROR(
"SiPixelGenError") <<
"Error reading file 17, no GenError load, run # "
426 << theCurrentTemp.entx[
k][
i].runnum <<
ENDL;
432 theCurrentTemp.entx[
k][
i].cotalpha = costrk[0] / costrk[2];
434 theCurrentTemp.entx[
k][
i].cotbeta = costrk[1] / costrk[2];
436 db >> theCurrentTemp.entx[
k][
i].qavg >> theCurrentTemp.entx[
k][
i].pixmax >> theCurrentTemp.entx[
k][
i].dyone >>
437 theCurrentTemp.entx[
k][
i].syone >> theCurrentTemp.entx[
k][
i].dxone >> theCurrentTemp.entx[
k][
i].sxone;
440 LOGERROR(
"SiPixelGenError") <<
"Error reading file 18, no GenError load, run # "
441 << theCurrentTemp.entx[
k][
i].runnum <<
ENDL;
445 db >> theCurrentTemp.entx[
k][
i].dytwo >> theCurrentTemp.entx[
k][
i].sytwo >> theCurrentTemp.entx[
k][
i].dxtwo >>
446 theCurrentTemp.entx[
k][
i].sxtwo >> theCurrentTemp.entx[
k][
i].qmin >> theCurrentTemp.entx[
k][
i].qmin2;
448 for (
j = 0;
j < 4; ++
j) {
449 db >> theCurrentTemp.entx[
k][
i].yavggen[
j] >> theCurrentTemp.entx[
k][
i].yrmsgen[
j] >>
450 theCurrentTemp.entx[
k][
i].xavggen[
j] >> theCurrentTemp.entx[
k][
i].xrmsgen[
j];
453 LOGERROR(
"SiPixelGenError") <<
"Error reading file 30a, no GenError load, run # "
454 << theCurrentTemp.entx[
k][
i].runnum <<
ENDL;
463 pixelTemp.push_back(theCurrentTemp);
474 for (
auto& templ : thePixelTemp_) {
475 for (
auto iy = 0; iy < templ.head.NTy; ++iy)
476 templ.cotbetaY[iy] = templ.enty[iy].cotbeta;
477 for (
auto iy = 0; iy < templ.head.NTyx; ++iy)
478 templ.cotbetaX[iy] = templ.entx[iy][0].cotbeta;
479 for (
auto ix = 0; ix < templ.head.NTxx; ++ix)
480 templ.cotalphaX[ix] = templ.entx[0][ix].cotalpha;
511 if (
id != id_current_) {
513 for (
int i = 0;
i < (
int)thePixelTemp_.size(); ++
i) {
514 if (
id == thePixelTemp_[
i].head.ID) {
518 lorywidth_ = thePixelTemp_[
i].head.lorywidth;
519 lorxwidth_ = thePixelTemp_[
i].head.lorxwidth;
520 lorybias_ = thePixelTemp_[
i].head.lorybias;
521 lorxbias_ = thePixelTemp_[
i].head.lorxbias;
526 xsize_ = thePixelTemp_[
i].head.xsize;
527 ysize_ = thePixelTemp_[
i].head.ysize;
528 zsize_ = thePixelTemp_[
i].head.zsize;
544 bool irradiationCorrections,
562 if (
id != id_current_) {
564 for (
int i = 0;
i < (
int)thePixelTemp_.size(); ++
i) {
565 if (
id == thePixelTemp_[
i].head.ID) {
568 lorywidth_ = thePixelTemp_[
i].head.lorywidth;
569 lorxwidth_ = thePixelTemp_[
i].head.lorxwidth;
570 lorybias_ = thePixelTemp_[
i].head.lorybias;
571 lorxbias_ = thePixelTemp_[
i].head.lorxbias;
572 for (
int j = 0;
j < 3; ++
j) {
573 fbin_[
j] = thePixelTemp_[
i].head.fbin[
j];
578 xsize_ = thePixelTemp_[
i].head.xsize;
579 ysize_ = thePixelTemp_[
i].head.ysize;
580 zsize_ = thePixelTemp_[
i].head.zsize;
587 int index = index_id_;
589 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
590 if (index < 0 || index >= (
int)thePixelTemp_.size()) {
591 throw cms::Exception(
"DataCorrupt") <<
"SiPixelGenError::qbin can't find needed GenError ID = " <<
id << std::endl;
599 auto const& templ = thePixelTemp_[
index];
607 auto cotalpha0 = thePixelTemp_[
index].enty[0].cotalpha;
609 std::sqrt((1.
f + cotbeta * cotbeta + cotalpha * cotalpha) / (1.
f + cotbeta * cotbeta + cotalpha0 * cotalpha0));
613 float cota = cotalpha;
620 switch (thePixelTemp_[index_id_].head.Dtype) {
638 if (locBx * locBz < 0.
f) {
650 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
652 <<
"SiPixelGenError::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
654 std::cout <<
"SiPixelGenError::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
661 lorybias_ = -lorybias_;
662 lorywidth_ = -lorywidth_;
665 lorxbias_ = -lorxbias_;
666 lorxwidth_ = -lorxwidth_;
669 auto qscale = thePixelTemp_[
index].head.qscale;
677 auto Ny = thePixelTemp_[
index].head.NTy;
678 auto Nyx = thePixelTemp_[
index].head.NTyx;
679 auto Nxx = thePixelTemp_[
index].head.NTxx;
681 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
682 if (
Ny < 2 || Nyx < 1 || Nxx < 2) {
683 throw cms::Exception(
"DataCorrupt") <<
"GenError ID = " << id_current_ <<
"has too few entries: Ny/Nyx/Nxx = " <<
Ny
684 <<
"/" << Nyx <<
"/" << Nxx << std::endl;
687 assert(
Ny > 1 && Nyx > 0 && Nxx > 1);
698 if (
j == templ.cotbetaY +
Ny) {
701 }
else if (
j == templ.cotbetaY) {
705 yratio = (cotb - (*(
j - 1))) / ((*
j) - (*(
j - 1)));
708 ihigh =
j - templ.cotbetaY;
714 auto qavg = (1.f - yratio) * thePixelTemp_[
index].enty[ilow].qavg + yratio * thePixelTemp_[
index].enty[ihigh].qavg;
716 auto qmin = (1.f - yratio) * thePixelTemp_[
index].enty[ilow].qmin + yratio * thePixelTemp_[
index].enty[ihigh].qmin;
718 auto qmin2 = (1.f - yratio) * thePixelTemp_[
index].enty[ilow].qmin2 + yratio * thePixelTemp_[
index].enty[ihigh].qmin2;
721 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
722 if (qavg <= 0.
f || qmin <= 0.
f) {
724 <<
"SiPixelGenError::qbin, qavg or qmin <= 0,"
725 <<
" Probably someone called the generic pixel reconstruction with an illegal trajectory state" << std::endl;
733 auto qtotal = qscale * qclus;
737 auto fq = qtotal / qavg;
753 auto yrmsgen = (1.f - yratio) * thePixelTemp_[
index].enty[ilow].yrmsgen[binq] +
754 yratio * thePixelTemp_[
index].enty[ihigh].yrmsgen[binq];
755 sy1 = (1.f - yratio) * thePixelTemp_[
index].enty[ilow].syone + yratio * thePixelTemp_[
index].enty[ihigh].syone;
756 sy2 = (1.f - yratio) * thePixelTemp_[
index].enty[ilow].sytwo + yratio * thePixelTemp_[
index].enty[ihigh].sytwo;
758 if (irradiationCorrections) {
759 auto yavggen = (1.f - yratio) * thePixelTemp_[
index].enty[ilow].yavggen[binq] +
760 yratio * thePixelTemp_[
index].enty[ihigh].yavggen[binq];
765 dy1 = (1.f - yratio) * thePixelTemp_[
index].enty[ilow].dyone + yratio * thePixelTemp_[
index].enty[ihigh].dyone;
769 dy2 = (1.f - yratio) * thePixelTemp_[
index].enty[ilow].dytwo + yratio * thePixelTemp_[
index].enty[ihigh].dytwo;
783 if (
j == templ.cotbetaX + Nyx) {
786 }
else if (
j == templ.cotbetaX) {
790 yxratio = (acotb - (*(
j - 1))) / ((*
j) - (*(
j - 1)));
793 iyhigh =
j - templ.cotbetaX;
801 if (
j == templ.cotalphaX + Nxx) {
804 }
else if (
j == templ.cotalphaX) {
808 xxratio = (cota - (*(
j - 1))) / ((*
j) - (*(
j - 1)));
811 ihigh =
j - templ.cotalphaX;
816 (1.f - xxratio) * thePixelTemp_[
index].entx[0][ilow].sxone + xxratio * thePixelTemp_[
index].entx[0][ihigh].sxone;
818 (1.f - xxratio) * thePixelTemp_[
index].entx[0][ilow].sxtwo + xxratio * thePixelTemp_[
index].entx[0][ihigh].sxtwo;
822 pixmx = (1.f - yxratio) * ((1.
f - xxratio) * thePixelTemp_[
index].entx[iylow][ilow].pixmax +
823 xxratio * thePixelTemp_[
index].entx[iylow][ihigh].pixmax) +
824 yxratio * ((1.
f - xxratio) * thePixelTemp_[
index].entx[iyhigh][ilow].pixmax +
825 xxratio * thePixelTemp_[
index].entx[iyhigh][ihigh].pixmax);
827 auto xrmsgen = (1.f - yxratio) * ((1.
f - xxratio) * thePixelTemp_[
index].entx[iylow][ilow].xrmsgen[binq] +
828 xxratio * thePixelTemp_[
index].entx[iylow][ihigh].xrmsgen[binq]) +
829 yxratio * ((1.
f - xxratio) * thePixelTemp_[
index].entx[iyhigh][ilow].xrmsgen[binq] +
830 xxratio * thePixelTemp_[
index].entx[iyhigh][ihigh].xrmsgen[binq]);
832 if (irradiationCorrections) {
833 auto xavggen = (1.f - yxratio) * ((1.
f - xxratio) * thePixelTemp_[
index].entx[iylow][ilow].xavggen[binq] +
834 xxratio * thePixelTemp_[
index].entx[iylow][ihigh].xavggen[binq]) +
835 yxratio * ((1.
f - xxratio) * thePixelTemp_[
index].entx[iyhigh][ilow].xavggen[binq] +
836 xxratio * thePixelTemp_[
index].entx[iyhigh][ihigh].xavggen[binq]);
841 dx1 = (1.f - xxratio) * thePixelTemp_[
index].entx[0][ilow].dxone +
842 xxratio * thePixelTemp_[
index].entx[0][ihigh].dxone;
846 dx2 = (1.f - xxratio) * thePixelTemp_[
index].entx[0][ilow].dxtwo +
847 xxratio * thePixelTemp_[
index].entx[0][ihigh].dxtwo;
861 if (qtotal < 0.95
f * qmin) {
864 if (qtotal < 0.95
f * qmin2) {
894 bool irradiationCorrections =
true;
903 irradiationCorrections,
918 pixmx = (
float)ipixmx;