File indexing completed on 2024-04-06 12:02:37
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0091 #include <cmath>
0092 #else
0093 #include <math.h>
0094 #endif
0095 #include <algorithm>
0096 #include <vector>
0097 #include "boost/multi_array.hpp"
0098 #include <iostream>
0099 #include <iomanip>
0100 #include <sstream>
0101 #include <fstream>
0102 #include <list>
0103
0104 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0105 #include "CondFormats/SiPixelTransient/interface/SiPixelTemplate.h"
0106 #include "CondFormats/SiPixelTransient/interface/SimplePixel.h"
0107 #include "FWCore/Utilities/interface/FileInPath.h"
0108 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0109 #include "FWCore/Utilities/interface/isFinite.h"
0110 #define LOGERROR(x) LogError(x)
0111 #define LOGINFO(x) LogInfo(x)
0112 #define LOGWARNING(x) LogWarning(x)
0113 #define ENDL " "
0114 #include "FWCore/Utilities/interface/Exception.h"
0115 using namespace edm;
0116 #else
0117 #include "SiPixelTemplate.h"
0118 #include "SimplePixel.h"
0119 #define LOGERROR(x) std::cout << x << ": "
0120 #define LOGINFO(x) std::cout << x << ": "
0121 #define ENDL std::endl
0122 #endif
0123
0124
0125
0126
0127
0128
0129
0130 bool SiPixelTemplate::pushfile(int filenum, std::vector<SiPixelTemplateStore>& pixelTemp, std::string dir) {
0131
0132
0133
0134 int i, j, k, l;
0135 float qavg_avg;
0136 char c;
0137 const int code_version = {17};
0138
0139
0140 std::string tempfile = std::to_string(filenum);
0141
0142 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0143
0144 int nzeros = 4 - tempfile.length();
0145 if (nzeros > 0)
0146 tempfile = std::string(nzeros, '0') + tempfile;
0147
0148
0149 tempfile = dir + "template_summary_zp" + tempfile + ".out";
0150 edm::FileInPath file(tempfile);
0151 tempfile = file.fullPath();
0152
0153 #else
0154
0155 std::ostringstream tout;
0156 tout << "template_summary_zp" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
0157 tempfile = tout.str();
0158
0159 #endif
0160
0161
0162
0163 std::ifstream in_file(tempfile);
0164 if (in_file.is_open() && in_file.good()) {
0165
0166
0167 SiPixelTemplateStore theCurrentTemp;
0168
0169
0170
0171 for (i = 0; (c = in_file.get()) != '\n'; ++i) {
0172 if (i < 79) {
0173 theCurrentTemp.head.title[i] = c;
0174 }
0175 }
0176 if (i > 78) {
0177 i = 78;
0178 }
0179 theCurrentTemp.head.title[i + 1] = '\0';
0180 LOGINFO("SiPixelTemplate") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
0181
0182
0183
0184 in_file >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >>
0185 theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx >> theCurrentTemp.head.Dtype >>
0186 theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >>
0187 theCurrentTemp.head.qscale >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >>
0188 theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >>
0189 theCurrentTemp.head.zsize;
0190
0191 if (in_file.fail()) {
0192 LOGERROR("SiPixelTemplate") << "Error reading file 0A, no template load" << ENDL;
0193 return false;
0194 }
0195
0196 if (theCurrentTemp.head.templ_version > 17) {
0197 in_file >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >> theCurrentTemp.head.lorxbias >>
0198 theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >> theCurrentTemp.head.fbin[2];
0199
0200 if (in_file.fail()) {
0201 LOGERROR("SiPixelTemplate") << "Error reading file 0B, no template load" << ENDL;
0202 return false;
0203 }
0204 } else {
0205 theCurrentTemp.head.ss50 = theCurrentTemp.head.s50;
0206 theCurrentTemp.head.lorybias = theCurrentTemp.head.lorywidth / 2.f;
0207 theCurrentTemp.head.lorxbias = theCurrentTemp.head.lorxwidth / 2.f;
0208 theCurrentTemp.head.fbin[0] = 1.5f;
0209 theCurrentTemp.head.fbin[1] = 1.00f;
0210 theCurrentTemp.head.fbin[2] = 0.85f;
0211 }
0212
0213 LOGINFO("SiPixelTemplate") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version "
0214 << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
0215 << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx
0216 << ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
0217 << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
0218 << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence
0219 << ", Q-scaling factor " << theCurrentTemp.head.qscale << ", 1/2 multi dcol threshold "
0220 << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold " << theCurrentTemp.head.ss50
0221 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", y Lorentz Bias "
0222 << theCurrentTemp.head.lorybias << ", x Lorentz width " << theCurrentTemp.head.lorxwidth
0223 << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
0224 << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] << ", "
0225 << theCurrentTemp.head.fbin[1] << ", " << theCurrentTemp.head.fbin[2]
0226 << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size "
0227 << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
0228
0229 if (theCurrentTemp.head.templ_version < code_version) {
0230 LOGERROR("SiPixelTemplate") << "code expects version " << code_version << " finds "
0231 << theCurrentTemp.head.templ_version << ", no template load" << ENDL;
0232 return false;
0233 }
0234
0235 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST
0236
0237
0238
0239 theCurrentTemp.cotbetaY = std::vector<float>(theCurrentTemp.head.NTy);
0240 theCurrentTemp.cotbetaX = std::vector<float>(theCurrentTemp.head.NTyx);
0241 theCurrentTemp.cotalphaX = std::vector<float>(theCurrentTemp.head.NTxx);
0242
0243 theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
0244
0245 theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
0246
0247 #endif
0248
0249
0250
0251 for (i = 0; i < theCurrentTemp.head.NTy; ++i) {
0252 in_file >> theCurrentTemp.enty[i].runnum >> theCurrentTemp.enty[i].costrk[0] >>
0253 theCurrentTemp.enty[i].costrk[1] >> theCurrentTemp.enty[i].costrk[2];
0254
0255 if (in_file.fail()) {
0256 LOGERROR("SiPixelTemplate") << "Error reading file 1, no template load, run # " << theCurrentTemp.enty[i].runnum
0257 << ENDL;
0258 return false;
0259 }
0260
0261
0262
0263 theCurrentTemp.enty[i].alpha =
0264 static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[0]));
0265
0266 theCurrentTemp.enty[i].cotalpha = theCurrentTemp.enty[i].costrk[0] / theCurrentTemp.enty[i].costrk[2];
0267
0268 theCurrentTemp.enty[i].beta =
0269 static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[1]));
0270
0271 theCurrentTemp.enty[i].cotbeta = theCurrentTemp.enty[i].costrk[1] / theCurrentTemp.enty[i].costrk[2];
0272
0273 in_file >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].pixmax >> theCurrentTemp.enty[i].symax >>
0274 theCurrentTemp.enty[i].dyone >> theCurrentTemp.enty[i].syone >> theCurrentTemp.enty[i].sxmax >>
0275 theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone;
0276
0277 if (in_file.fail()) {
0278 LOGERROR("SiPixelTemplate") << "Error reading file 2, no template load, run # " << theCurrentTemp.enty[i].runnum
0279 << ENDL;
0280 return false;
0281 }
0282
0283 in_file >> theCurrentTemp.enty[i].dytwo >> theCurrentTemp.enty[i].sytwo >> theCurrentTemp.enty[i].dxtwo >>
0284 theCurrentTemp.enty[i].sxtwo >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].clsleny >>
0285 theCurrentTemp.enty[i].clslenx;
0286
0287 if (in_file.fail()) {
0288 LOGERROR("SiPixelTemplate") << "Error reading file 3, no template load, run # " << theCurrentTemp.enty[i].runnum
0289 << ENDL;
0290 return false;
0291 }
0292
0293 if (theCurrentTemp.enty[i].qmin <= 0.) {
0294 LOGERROR("SiPixelTemplate") << "Error in template ID " << theCurrentTemp.head.ID
0295 << " qmin = " << theCurrentTemp.enty[i].qmin << ", run # "
0296 << theCurrentTemp.enty[i].runnum << ENDL;
0297 return false;
0298 }
0299
0300 for (j = 0; j < 2; ++j) {
0301 in_file >> theCurrentTemp.enty[i].ypar[j][0] >> theCurrentTemp.enty[i].ypar[j][1] >>
0302 theCurrentTemp.enty[i].ypar[j][2] >> theCurrentTemp.enty[i].ypar[j][3] >> theCurrentTemp.enty[i].ypar[j][4];
0303
0304 if (in_file.fail()) {
0305 LOGERROR("SiPixelTemplate") << "Error reading file 4, no template load, run # "
0306 << theCurrentTemp.enty[i].runnum << ENDL;
0307 return false;
0308 }
0309 }
0310
0311 for (j = 0; j < 9; ++j) {
0312 for (k = 0; k < TYSIZE; ++k) {
0313 in_file >> theCurrentTemp.enty[i].ytemp[j][k];
0314 }
0315
0316 if (in_file.fail()) {
0317 LOGERROR("SiPixelTemplate") << "Error reading file 5, no template load, run # "
0318 << theCurrentTemp.enty[i].runnum << ENDL;
0319 return false;
0320 }
0321 }
0322
0323 for (j = 0; j < 2; ++j) {
0324 in_file >> theCurrentTemp.enty[i].xpar[j][0] >> theCurrentTemp.enty[i].xpar[j][1] >>
0325 theCurrentTemp.enty[i].xpar[j][2] >> theCurrentTemp.enty[i].xpar[j][3] >> theCurrentTemp.enty[i].xpar[j][4];
0326
0327 if (in_file.fail()) {
0328 LOGERROR("SiPixelTemplate") << "Error reading file 6, no template load, run # "
0329 << theCurrentTemp.enty[i].runnum << ENDL;
0330 return false;
0331 }
0332 }
0333
0334 qavg_avg = 0.f;
0335 for (j = 0; j < 9; ++j) {
0336 for (k = 0; k < TXSIZE; ++k) {
0337 in_file >> theCurrentTemp.enty[i].xtemp[j][k];
0338 qavg_avg += theCurrentTemp.enty[i].xtemp[j][k];
0339 }
0340
0341 if (in_file.fail()) {
0342 LOGERROR("SiPixelTemplate") << "Error reading file 7, no template load, run # "
0343 << theCurrentTemp.enty[i].runnum << ENDL;
0344 return false;
0345 }
0346 }
0347 theCurrentTemp.enty[i].qavg_avg = qavg_avg / 9.;
0348
0349 for (j = 0; j < 4; ++j) {
0350 in_file >> theCurrentTemp.enty[i].yavg[j] >> theCurrentTemp.enty[i].yrms[j] >> theCurrentTemp.enty[i].ygx0[j] >>
0351 theCurrentTemp.enty[i].ygsig[j];
0352
0353 if (in_file.fail()) {
0354 LOGERROR("SiPixelTemplate") << "Error reading file 8, no template load, run # "
0355 << theCurrentTemp.enty[i].runnum << ENDL;
0356 return false;
0357 }
0358 }
0359
0360 for (j = 0; j < 4; ++j) {
0361 in_file >> theCurrentTemp.enty[i].yflpar[j][0] >> theCurrentTemp.enty[i].yflpar[j][1] >>
0362 theCurrentTemp.enty[i].yflpar[j][2] >> theCurrentTemp.enty[i].yflpar[j][3] >>
0363 theCurrentTemp.enty[i].yflpar[j][4] >> theCurrentTemp.enty[i].yflpar[j][5];
0364
0365 if (in_file.fail()) {
0366 LOGERROR("SiPixelTemplate") << "Error reading file 9, no template load, run # "
0367 << theCurrentTemp.enty[i].runnum << ENDL;
0368 return false;
0369 }
0370 }
0371
0372 for (j = 0; j < 4; ++j) {
0373 in_file >> theCurrentTemp.enty[i].xavg[j] >> theCurrentTemp.enty[i].xrms[j] >> theCurrentTemp.enty[i].xgx0[j] >>
0374 theCurrentTemp.enty[i].xgsig[j];
0375
0376 if (in_file.fail()) {
0377 LOGERROR("SiPixelTemplate") << "Error reading file 10, no template load, run # "
0378 << theCurrentTemp.enty[i].runnum << ENDL;
0379 return false;
0380 }
0381 }
0382
0383 for (j = 0; j < 4; ++j) {
0384 in_file >> theCurrentTemp.enty[i].xflpar[j][0] >> theCurrentTemp.enty[i].xflpar[j][1] >>
0385 theCurrentTemp.enty[i].xflpar[j][2] >> theCurrentTemp.enty[i].xflpar[j][3] >>
0386 theCurrentTemp.enty[i].xflpar[j][4] >> theCurrentTemp.enty[i].xflpar[j][5];
0387
0388 if (in_file.fail()) {
0389 LOGERROR("SiPixelTemplate") << "Error reading file 11, no template load, run # "
0390 << theCurrentTemp.enty[i].runnum << ENDL;
0391 return false;
0392 }
0393 }
0394
0395 for (j = 0; j < 4; ++j) {
0396 in_file >> theCurrentTemp.enty[i].chi2yavg[j] >> theCurrentTemp.enty[i].chi2ymin[j] >>
0397 theCurrentTemp.enty[i].chi2xavg[j] >> theCurrentTemp.enty[i].chi2xmin[j];
0398
0399 if (in_file.fail()) {
0400 LOGERROR("SiPixelTemplate") << "Error reading file 12, no template load, run # "
0401 << theCurrentTemp.enty[i].runnum << ENDL;
0402 return false;
0403 }
0404 }
0405
0406 for (j = 0; j < 4; ++j) {
0407 in_file >> theCurrentTemp.enty[i].yavgc2m[j] >> theCurrentTemp.enty[i].yrmsc2m[j] >>
0408 theCurrentTemp.enty[i].chi2yavgc2m[j] >> theCurrentTemp.enty[i].chi2yminc2m[j];
0409
0410 if (in_file.fail()) {
0411 LOGERROR("SiPixelTemplate") << "Error reading file 13, no template load, run # "
0412 << theCurrentTemp.enty[i].runnum << ENDL;
0413 return false;
0414 }
0415 }
0416
0417 for (j = 0; j < 4; ++j) {
0418 in_file >> theCurrentTemp.enty[i].xavgc2m[j] >> theCurrentTemp.enty[i].xrmsc2m[j] >>
0419 theCurrentTemp.enty[i].chi2xavgc2m[j] >> theCurrentTemp.enty[i].chi2xminc2m[j];
0420
0421 if (in_file.fail()) {
0422 LOGERROR("SiPixelTemplate") << "Error reading file 14, no template load, run # "
0423 << theCurrentTemp.enty[i].runnum << ENDL;
0424 return false;
0425 }
0426 }
0427
0428 for (j = 0; j < 4; ++j) {
0429 in_file >> theCurrentTemp.enty[i].yavggen[j] >> theCurrentTemp.enty[i].yrmsgen[j] >>
0430 theCurrentTemp.enty[i].ygx0gen[j] >> theCurrentTemp.enty[i].ygsiggen[j];
0431
0432 if (in_file.fail()) {
0433 LOGERROR("SiPixelTemplate") << "Error reading file 14a, no template load, run # "
0434 << theCurrentTemp.enty[i].runnum << ENDL;
0435 return false;
0436 }
0437 }
0438
0439 for (j = 0; j < 4; ++j) {
0440 in_file >> theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j] >>
0441 theCurrentTemp.enty[i].xgx0gen[j] >> theCurrentTemp.enty[i].xgsiggen[j];
0442
0443 if (in_file.fail()) {
0444 LOGERROR("SiPixelTemplate") << "Error reading file 14b, no template load, run # "
0445 << theCurrentTemp.enty[i].runnum << ENDL;
0446 return false;
0447 }
0448 }
0449
0450 in_file >> theCurrentTemp.enty[i].chi2yavgone >> theCurrentTemp.enty[i].chi2yminone >>
0451 theCurrentTemp.enty[i].chi2xavgone >> theCurrentTemp.enty[i].chi2xminone >> theCurrentTemp.enty[i].qmin2 >>
0452 theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav >>
0453 theCurrentTemp.enty[i].r_qMeas_qTrue >> theCurrentTemp.enty[i].spare[0];
0454
0455 if (in_file.fail()) {
0456 LOGERROR("SiPixelTemplate") << "Error reading file 15, no template load, run # "
0457 << theCurrentTemp.enty[i].runnum << ENDL;
0458 return false;
0459 }
0460
0461 in_file >> theCurrentTemp.enty[i].mpvvav2 >> theCurrentTemp.enty[i].sigmavav2 >>
0462 theCurrentTemp.enty[i].kappavav2 >> theCurrentTemp.enty[i].qbfrac[0] >> theCurrentTemp.enty[i].qbfrac[1] >>
0463 theCurrentTemp.enty[i].qbfrac[2] >> theCurrentTemp.enty[i].fracyone >> theCurrentTemp.enty[i].fracxone >>
0464 theCurrentTemp.enty[i].fracytwo >> theCurrentTemp.enty[i].fracxtwo;
0465
0466
0467 if (in_file.fail()) {
0468 LOGERROR("SiPixelTemplate") << "Error reading file 16, no template load, run # "
0469 << theCurrentTemp.enty[i].runnum << ENDL;
0470 return false;
0471 }
0472 }
0473
0474
0475
0476 for (k = 0; k < theCurrentTemp.head.NTyx; ++k) {
0477 for (i = 0; i < theCurrentTemp.head.NTxx; ++i) {
0478 in_file >> theCurrentTemp.entx[k][i].runnum >> theCurrentTemp.entx[k][i].costrk[0] >>
0479 theCurrentTemp.entx[k][i].costrk[1] >> theCurrentTemp.entx[k][i].costrk[2];
0480
0481 if (in_file.fail()) {
0482 LOGERROR("SiPixelTemplate") << "Error reading file 17, no template load, run # "
0483 << theCurrentTemp.entx[k][i].runnum << ENDL;
0484 return false;
0485 }
0486
0487
0488
0489 theCurrentTemp.entx[k][i].alpha = static_cast<float>(
0490 atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[0]));
0491
0492 theCurrentTemp.entx[k][i].cotalpha = theCurrentTemp.entx[k][i].costrk[0] / theCurrentTemp.entx[k][i].costrk[2];
0493
0494 theCurrentTemp.entx[k][i].beta = static_cast<float>(
0495 atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[1]));
0496
0497 theCurrentTemp.entx[k][i].cotbeta = theCurrentTemp.entx[k][i].costrk[1] / theCurrentTemp.entx[k][i].costrk[2];
0498
0499 in_file >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].pixmax >>
0500 theCurrentTemp.entx[k][i].symax >> theCurrentTemp.entx[k][i].dyone >> theCurrentTemp.entx[k][i].syone >>
0501 theCurrentTemp.entx[k][i].sxmax >> theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone;
0502
0503 if (in_file.fail()) {
0504 LOGERROR("SiPixelTemplate") << "Error reading file 18, no template load, run # "
0505 << theCurrentTemp.entx[k][i].runnum << ENDL;
0506 return false;
0507 }
0508
0509 in_file >> theCurrentTemp.entx[k][i].dytwo >> theCurrentTemp.entx[k][i].sytwo >>
0510 theCurrentTemp.entx[k][i].dxtwo >> theCurrentTemp.entx[k][i].sxtwo >> theCurrentTemp.entx[k][i].qmin >>
0511 theCurrentTemp.entx[k][i].clsleny >> theCurrentTemp.entx[k][i].clslenx;
0512
0513
0514 if (in_file.fail()) {
0515 LOGERROR("SiPixelTemplate") << "Error reading file 19, no template load, run # "
0516 << theCurrentTemp.entx[k][i].runnum << ENDL;
0517 return false;
0518 }
0519
0520 for (j = 0; j < 2; ++j) {
0521 in_file >> theCurrentTemp.entx[k][i].ypar[j][0] >> theCurrentTemp.entx[k][i].ypar[j][1] >>
0522 theCurrentTemp.entx[k][i].ypar[j][2] >> theCurrentTemp.entx[k][i].ypar[j][3] >>
0523 theCurrentTemp.entx[k][i].ypar[j][4];
0524
0525 if (in_file.fail()) {
0526 LOGERROR("SiPixelTemplate") << "Error reading file 20, no template load, run # "
0527 << theCurrentTemp.entx[k][i].runnum << ENDL;
0528 return false;
0529 }
0530 }
0531
0532 for (j = 0; j < 9; ++j) {
0533 for (l = 0; l < TYSIZE; ++l) {
0534 in_file >> theCurrentTemp.entx[k][i].ytemp[j][l];
0535 }
0536
0537 if (in_file.fail()) {
0538 LOGERROR("SiPixelTemplate") << "Error reading file 21, no template load, run # "
0539 << theCurrentTemp.entx[k][i].runnum << ENDL;
0540 return false;
0541 }
0542 }
0543
0544 for (j = 0; j < 2; ++j) {
0545 in_file >> theCurrentTemp.entx[k][i].xpar[j][0] >> theCurrentTemp.entx[k][i].xpar[j][1] >>
0546 theCurrentTemp.entx[k][i].xpar[j][2] >> theCurrentTemp.entx[k][i].xpar[j][3] >>
0547 theCurrentTemp.entx[k][i].xpar[j][4];
0548
0549 if (in_file.fail()) {
0550 LOGERROR("SiPixelTemplate") << "Error reading file 22, no template load, run # "
0551 << theCurrentTemp.entx[k][i].runnum << ENDL;
0552 return false;
0553 }
0554 }
0555
0556 qavg_avg = 0.f;
0557 for (j = 0; j < 9; ++j) {
0558 for (l = 0; l < TXSIZE; ++l) {
0559 in_file >> theCurrentTemp.entx[k][i].xtemp[j][l];
0560 qavg_avg += theCurrentTemp.entx[k][i].xtemp[j][l];
0561 }
0562
0563 if (in_file.fail()) {
0564 LOGERROR("SiPixelTemplate") << "Error reading file 23, no template load, run # "
0565 << theCurrentTemp.entx[k][i].runnum << ENDL;
0566 return false;
0567 }
0568 }
0569 theCurrentTemp.entx[k][i].qavg_avg = qavg_avg / 9.;
0570
0571 for (j = 0; j < 4; ++j) {
0572 in_file >> theCurrentTemp.entx[k][i].yavg[j] >> theCurrentTemp.entx[k][i].yrms[j] >>
0573 theCurrentTemp.entx[k][i].ygx0[j] >> theCurrentTemp.entx[k][i].ygsig[j];
0574
0575 if (in_file.fail()) {
0576 LOGERROR("SiPixelTemplate") << "Error reading file 24, no template load, run # "
0577 << theCurrentTemp.entx[k][i].runnum << ENDL;
0578 return false;
0579 }
0580 }
0581
0582 for (j = 0; j < 4; ++j) {
0583 in_file >> theCurrentTemp.entx[k][i].yflpar[j][0] >> theCurrentTemp.entx[k][i].yflpar[j][1] >>
0584 theCurrentTemp.entx[k][i].yflpar[j][2] >> theCurrentTemp.entx[k][i].yflpar[j][3] >>
0585 theCurrentTemp.entx[k][i].yflpar[j][4] >> theCurrentTemp.entx[k][i].yflpar[j][5];
0586
0587 if (in_file.fail()) {
0588 LOGERROR("SiPixelTemplate") << "Error reading file 25, no template load, run # "
0589 << theCurrentTemp.entx[k][i].runnum << ENDL;
0590 return false;
0591 }
0592 }
0593
0594 for (j = 0; j < 4; ++j) {
0595 in_file >> theCurrentTemp.entx[k][i].xavg[j] >> theCurrentTemp.entx[k][i].xrms[j] >>
0596 theCurrentTemp.entx[k][i].xgx0[j] >> theCurrentTemp.entx[k][i].xgsig[j];
0597
0598 if (in_file.fail()) {
0599 LOGERROR("SiPixelTemplate") << "Error reading file 26, no template load, run # "
0600 << theCurrentTemp.entx[k][i].runnum << ENDL;
0601 return false;
0602 }
0603 }
0604
0605 for (j = 0; j < 4; ++j) {
0606 in_file >> theCurrentTemp.entx[k][i].xflpar[j][0] >> theCurrentTemp.entx[k][i].xflpar[j][1] >>
0607 theCurrentTemp.entx[k][i].xflpar[j][2] >> theCurrentTemp.entx[k][i].xflpar[j][3] >>
0608 theCurrentTemp.entx[k][i].xflpar[j][4] >> theCurrentTemp.entx[k][i].xflpar[j][5];
0609
0610 if (in_file.fail()) {
0611 LOGERROR("SiPixelTemplate") << "Error reading file 27, no template load, run # "
0612 << theCurrentTemp.entx[k][i].runnum << ENDL;
0613 return false;
0614 }
0615 }
0616
0617 for (j = 0; j < 4; ++j) {
0618 in_file >> theCurrentTemp.entx[k][i].chi2yavg[j] >> theCurrentTemp.entx[k][i].chi2ymin[j] >>
0619 theCurrentTemp.entx[k][i].chi2xavg[j] >> theCurrentTemp.entx[k][i].chi2xmin[j];
0620
0621 if (in_file.fail()) {
0622 LOGERROR("SiPixelTemplate") << "Error reading file 28, no template load, run # "
0623 << theCurrentTemp.entx[k][i].runnum << ENDL;
0624 return false;
0625 }
0626 }
0627
0628 for (j = 0; j < 4; ++j) {
0629 in_file >> theCurrentTemp.entx[k][i].yavgc2m[j] >> theCurrentTemp.entx[k][i].yrmsc2m[j] >>
0630 theCurrentTemp.entx[k][i].chi2yavgc2m[j] >> theCurrentTemp.entx[k][i].chi2yminc2m[j];
0631
0632 if (in_file.fail()) {
0633 LOGERROR("SiPixelTemplate") << "Error reading file 29, no template load, run # "
0634 << theCurrentTemp.entx[k][i].runnum << ENDL;
0635 return false;
0636 }
0637 }
0638
0639 for (j = 0; j < 4; ++j) {
0640 in_file >> theCurrentTemp.entx[k][i].xavgc2m[j] >> theCurrentTemp.entx[k][i].xrmsc2m[j] >>
0641 theCurrentTemp.entx[k][i].chi2xavgc2m[j] >> theCurrentTemp.entx[k][i].chi2xminc2m[j];
0642
0643 if (in_file.fail()) {
0644 LOGERROR("SiPixelTemplate") << "Error reading file 30, no template load, run # "
0645 << theCurrentTemp.entx[k][i].runnum << ENDL;
0646 return false;
0647 }
0648 }
0649
0650 for (j = 0; j < 4; ++j) {
0651 in_file >> theCurrentTemp.entx[k][i].yavggen[j] >> theCurrentTemp.entx[k][i].yrmsgen[j] >>
0652 theCurrentTemp.entx[k][i].ygx0gen[j] >> theCurrentTemp.entx[k][i].ygsiggen[j];
0653
0654 if (in_file.fail()) {
0655 LOGERROR("SiPixelTemplate") << "Error reading file 30a, no template load, run # "
0656 << theCurrentTemp.entx[k][i].runnum << ENDL;
0657 return false;
0658 }
0659 }
0660
0661 for (j = 0; j < 4; ++j) {
0662 in_file >> theCurrentTemp.entx[k][i].xavggen[j] >> theCurrentTemp.entx[k][i].xrmsgen[j] >>
0663 theCurrentTemp.entx[k][i].xgx0gen[j] >> theCurrentTemp.entx[k][i].xgsiggen[j];
0664
0665 if (in_file.fail()) {
0666 LOGERROR("SiPixelTemplate") << "Error reading file 30b, no template load, run # "
0667 << theCurrentTemp.entx[k][i].runnum << ENDL;
0668 return false;
0669 }
0670 }
0671
0672 in_file >> theCurrentTemp.entx[k][i].chi2yavgone >> theCurrentTemp.entx[k][i].chi2yminone >>
0673 theCurrentTemp.entx[k][i].chi2xavgone >> theCurrentTemp.entx[k][i].chi2xminone >>
0674 theCurrentTemp.entx[k][i].qmin2 >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >>
0675 theCurrentTemp.entx[k][i].kappavav >> theCurrentTemp.entx[k][i].r_qMeas_qTrue >>
0676 theCurrentTemp.entx[k][i].spare[0];
0677
0678 if (in_file.fail()) {
0679 LOGERROR("SiPixelTemplate") << "Error reading file 31, no template load, run # "
0680 << theCurrentTemp.entx[k][i].runnum << ENDL;
0681 return false;
0682 }
0683
0684 in_file >> theCurrentTemp.entx[k][i].mpvvav2 >> theCurrentTemp.entx[k][i].sigmavav2 >>
0685 theCurrentTemp.entx[k][i].kappavav2 >> theCurrentTemp.entx[k][i].qbfrac[0] >>
0686 theCurrentTemp.entx[k][i].qbfrac[1] >> theCurrentTemp.entx[k][i].qbfrac[2] >>
0687 theCurrentTemp.entx[k][i].fracyone >> theCurrentTemp.entx[k][i].fracxone >>
0688 theCurrentTemp.entx[k][i].fracytwo >> theCurrentTemp.entx[k][i].fracxtwo;
0689
0690
0691 if (in_file.fail()) {
0692 LOGERROR("SiPixelTemplate") << "Error reading file 32, no template load, run # "
0693 << theCurrentTemp.entx[k][i].runnum << ENDL;
0694 return false;
0695 }
0696 }
0697 }
0698
0699 in_file.close();
0700
0701
0702
0703 pixelTemp.push_back(theCurrentTemp);
0704
0705 postInit(pixelTemp);
0706 return true;
0707
0708 } else {
0709
0710
0711 LOGERROR("SiPixelTemplate") << "Error opening File" << tempfile << ENDL;
0712 return false;
0713 }
0714
0715 }
0716
0717 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0718
0719
0720
0721
0722
0723
0724 bool SiPixelTemplate::pushfile(const SiPixelTemplateDBObject& dbobject, std::vector<SiPixelTemplateStore>& pixelTemp) {
0725
0726
0727
0728 int i, j, k, l;
0729 float qavg_avg;
0730 const int code_version = {17};
0731
0732
0733 auto db(dbobject.reader());
0734
0735
0736
0737 auto tmpPtr = std::make_unique<SiPixelTemplateStore>();
0738 auto& theCurrentTemp = *tmpPtr;
0739
0740
0741 for (int m = 0; m < db.numOfTempl(); ++m) {
0742
0743
0744 SiPixelTemplateDBObject::char2float temp;
0745 for (i = 0; i < 20; ++i) {
0746 temp.f = db.sVector()[db.index()];
0747 theCurrentTemp.head.title[4 * i] = temp.c[0];
0748 theCurrentTemp.head.title[4 * i + 1] = temp.c[1];
0749 theCurrentTemp.head.title[4 * i + 2] = temp.c[2];
0750 theCurrentTemp.head.title[4 * i + 3] = temp.c[3];
0751 db.incrementIndex(1);
0752 }
0753 theCurrentTemp.head.title[79] = '\0';
0754 LOGINFO("SiPixelTemplate") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
0755
0756
0757
0758 db >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >>
0759 theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx >> theCurrentTemp.head.Dtype >>
0760 theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >>
0761 theCurrentTemp.head.qscale >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >>
0762 theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >>
0763 theCurrentTemp.head.zsize;
0764
0765 if (db.fail()) {
0766 LOGERROR("SiPixelTemplate") << "Error reading file 0A, no template load" << ENDL;
0767 return false;
0768 }
0769
0770 LOGINFO("SiPixelTemplate") << "Loading Pixel Template File - " << theCurrentTemp.head.title
0771 << " code version = " << code_version << " object version "
0772 << theCurrentTemp.head.templ_version << ENDL;
0773
0774 if (theCurrentTemp.head.templ_version > 17) {
0775 db >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >> theCurrentTemp.head.lorxbias >>
0776 theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >> theCurrentTemp.head.fbin[2];
0777
0778 if (db.fail()) {
0779 LOGERROR("SiPixelTemplate") << "Error reading file 0B, no template load" << ENDL;
0780 return false;
0781 }
0782 } else {
0783 theCurrentTemp.head.ss50 = theCurrentTemp.head.s50;
0784 theCurrentTemp.head.lorybias = theCurrentTemp.head.lorywidth / 2.f;
0785 theCurrentTemp.head.lorxbias = theCurrentTemp.head.lorxwidth / 2.f;
0786 theCurrentTemp.head.fbin[0] = 1.50f;
0787 theCurrentTemp.head.fbin[1] = 1.00f;
0788 theCurrentTemp.head.fbin[2] = 0.85f;
0789
0790
0791 }
0792
0793 LOGINFO("SiPixelTemplate") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version "
0794 << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
0795 << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx
0796 << ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
0797 << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
0798 << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence
0799 << ", Q-scaling factor " << theCurrentTemp.head.qscale << ", 1/2 multi dcol threshold "
0800 << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold " << theCurrentTemp.head.ss50
0801 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", y Lorentz Bias "
0802 << theCurrentTemp.head.lorybias << ", x Lorentz width " << theCurrentTemp.head.lorxwidth
0803 << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
0804 << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] << ", "
0805 << theCurrentTemp.head.fbin[1] << ", " << theCurrentTemp.head.fbin[2]
0806 << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size "
0807 << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
0808
0809 if (theCurrentTemp.head.templ_version < code_version) {
0810 LOGINFO("SiPixelTemplate") << "code expects version " << code_version << " finds "
0811 << theCurrentTemp.head.templ_version << ", load anyway " << ENDL;
0812
0813 }
0814
0815 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST
0816
0817
0818 theCurrentTemp.cotbetaY = std::vector<float>(theCurrentTemp.head.NTy);
0819 theCurrentTemp.cotbetaX = std::vector<float>(theCurrentTemp.head.NTyx);
0820 theCurrentTemp.cotalphaX = std::vector<float>(theCurrentTemp.head.NTxx);
0821 theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
0822 theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
0823
0824 #endif
0825
0826
0827
0828 for (i = 0; i < theCurrentTemp.head.NTy; ++i) {
0829 db >> theCurrentTemp.enty[i].runnum >> theCurrentTemp.enty[i].costrk[0] >> theCurrentTemp.enty[i].costrk[1] >>
0830 theCurrentTemp.enty[i].costrk[2];
0831
0832 if (db.fail()) {
0833 LOGERROR("SiPixelTemplate") << "Error reading file 1, no template load, run # " << theCurrentTemp.enty[i].runnum
0834 << ENDL;
0835 return false;
0836 }
0837
0838
0839
0840 theCurrentTemp.enty[i].alpha =
0841 static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[0]));
0842
0843 theCurrentTemp.enty[i].cotalpha = theCurrentTemp.enty[i].costrk[0] / theCurrentTemp.enty[i].costrk[2];
0844
0845 theCurrentTemp.enty[i].beta =
0846 static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[1]));
0847
0848 theCurrentTemp.enty[i].cotbeta = theCurrentTemp.enty[i].costrk[1] / theCurrentTemp.enty[i].costrk[2];
0849
0850 db >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].pixmax >> theCurrentTemp.enty[i].symax >>
0851 theCurrentTemp.enty[i].dyone >> theCurrentTemp.enty[i].syone >> theCurrentTemp.enty[i].sxmax >>
0852 theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone;
0853
0854 if (db.fail()) {
0855 LOGERROR("SiPixelTemplate") << "Error reading file 2, no template load, run # " << theCurrentTemp.enty[i].runnum
0856 << ENDL;
0857 return false;
0858 }
0859
0860 db >> theCurrentTemp.enty[i].dytwo >> theCurrentTemp.enty[i].sytwo >> theCurrentTemp.enty[i].dxtwo >>
0861 theCurrentTemp.enty[i].sxtwo >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].clsleny >>
0862 theCurrentTemp.enty[i].clslenx;
0863
0864
0865 if (db.fail()) {
0866 LOGERROR("SiPixelTemplate") << "Error reading file 3, no template load, run # " << theCurrentTemp.enty[i].runnum
0867 << ENDL;
0868 return false;
0869 }
0870
0871 if (theCurrentTemp.enty[i].qmin <= 0.) {
0872 LOGERROR("SiPixelTemplate") << "Error in template ID " << theCurrentTemp.head.ID
0873 << " qmin = " << theCurrentTemp.enty[i].qmin << ", run # "
0874 << theCurrentTemp.enty[i].runnum << ENDL;
0875 return false;
0876 }
0877
0878 for (j = 0; j < 2; ++j) {
0879 db >> theCurrentTemp.enty[i].ypar[j][0] >> theCurrentTemp.enty[i].ypar[j][1] >>
0880 theCurrentTemp.enty[i].ypar[j][2] >> theCurrentTemp.enty[i].ypar[j][3] >> theCurrentTemp.enty[i].ypar[j][4];
0881
0882 if (db.fail()) {
0883 LOGERROR("SiPixelTemplate") << "Error reading file 4, no template load, run # "
0884 << theCurrentTemp.enty[i].runnum << ENDL;
0885 return false;
0886 }
0887 }
0888
0889 for (j = 0; j < 9; ++j) {
0890 for (k = 0; k < TYSIZE; ++k) {
0891 db >> theCurrentTemp.enty[i].ytemp[j][k];
0892 }
0893
0894 if (db.fail()) {
0895 LOGERROR("SiPixelTemplate") << "Error reading file 5, no template load, run # "
0896 << theCurrentTemp.enty[i].runnum << ENDL;
0897 return false;
0898 }
0899 }
0900
0901 for (j = 0; j < 2; ++j) {
0902 db >> theCurrentTemp.enty[i].xpar[j][0] >> theCurrentTemp.enty[i].xpar[j][1] >>
0903 theCurrentTemp.enty[i].xpar[j][2] >> theCurrentTemp.enty[i].xpar[j][3] >> theCurrentTemp.enty[i].xpar[j][4];
0904
0905 if (db.fail()) {
0906 LOGERROR("SiPixelTemplate") << "Error reading file 6, no template load, run # "
0907 << theCurrentTemp.enty[i].runnum << ENDL;
0908 return false;
0909 }
0910 }
0911
0912 qavg_avg = 0.f;
0913 for (j = 0; j < 9; ++j) {
0914 for (k = 0; k < TXSIZE; ++k) {
0915 db >> theCurrentTemp.enty[i].xtemp[j][k];
0916 qavg_avg += theCurrentTemp.enty[i].xtemp[j][k];
0917 }
0918
0919 if (db.fail()) {
0920 LOGERROR("SiPixelTemplate") << "Error reading file 7, no template load, run # "
0921 << theCurrentTemp.enty[i].runnum << ENDL;
0922 return false;
0923 }
0924 }
0925 theCurrentTemp.enty[i].qavg_avg = qavg_avg / 9.;
0926
0927 for (j = 0; j < 4; ++j) {
0928 db >> theCurrentTemp.enty[i].yavg[j] >> theCurrentTemp.enty[i].yrms[j] >> theCurrentTemp.enty[i].ygx0[j] >>
0929 theCurrentTemp.enty[i].ygsig[j];
0930
0931 if (db.fail()) {
0932 LOGERROR("SiPixelTemplate") << "Error reading file 8, no template load, run # "
0933 << theCurrentTemp.enty[i].runnum << ENDL;
0934 return false;
0935 }
0936 }
0937
0938 for (j = 0; j < 4; ++j) {
0939 db >> theCurrentTemp.enty[i].yflpar[j][0] >> theCurrentTemp.enty[i].yflpar[j][1] >>
0940 theCurrentTemp.enty[i].yflpar[j][2] >> theCurrentTemp.enty[i].yflpar[j][3] >>
0941 theCurrentTemp.enty[i].yflpar[j][4] >> theCurrentTemp.enty[i].yflpar[j][5];
0942
0943 if (db.fail()) {
0944 LOGERROR("SiPixelTemplate") << "Error reading file 9, no template load, run # "
0945 << theCurrentTemp.enty[i].runnum << ENDL;
0946 return false;
0947 }
0948 }
0949
0950 for (j = 0; j < 4; ++j) {
0951 db >> theCurrentTemp.enty[i].xavg[j] >> theCurrentTemp.enty[i].xrms[j] >> theCurrentTemp.enty[i].xgx0[j] >>
0952 theCurrentTemp.enty[i].xgsig[j];
0953
0954 if (db.fail()) {
0955 LOGERROR("SiPixelTemplate") << "Error reading file 10, no template load, run # "
0956 << theCurrentTemp.enty[i].runnum << ENDL;
0957 return false;
0958 }
0959 }
0960
0961 for (j = 0; j < 4; ++j) {
0962 db >> theCurrentTemp.enty[i].xflpar[j][0] >> theCurrentTemp.enty[i].xflpar[j][1] >>
0963 theCurrentTemp.enty[i].xflpar[j][2] >> theCurrentTemp.enty[i].xflpar[j][3] >>
0964 theCurrentTemp.enty[i].xflpar[j][4] >> theCurrentTemp.enty[i].xflpar[j][5];
0965
0966 if (db.fail()) {
0967 LOGERROR("SiPixelTemplate") << "Error reading file 11, no template load, run # "
0968 << theCurrentTemp.enty[i].runnum << ENDL;
0969 return false;
0970 }
0971 }
0972
0973 for (j = 0; j < 4; ++j) {
0974 db >> theCurrentTemp.enty[i].chi2yavg[j] >> theCurrentTemp.enty[i].chi2ymin[j] >>
0975 theCurrentTemp.enty[i].chi2xavg[j] >> theCurrentTemp.enty[i].chi2xmin[j];
0976
0977 if (db.fail()) {
0978 LOGERROR("SiPixelTemplate") << "Error reading file 12, no template load, run # "
0979 << theCurrentTemp.enty[i].runnum << ENDL;
0980 return false;
0981 }
0982 }
0983
0984 for (j = 0; j < 4; ++j) {
0985 db >> theCurrentTemp.enty[i].yavgc2m[j] >> theCurrentTemp.enty[i].yrmsc2m[j] >>
0986 theCurrentTemp.enty[i].chi2yavgc2m[j] >> theCurrentTemp.enty[i].chi2yminc2m[j];
0987
0988 if (db.fail()) {
0989 LOGERROR("SiPixelTemplate") << "Error reading file 13, no template load, run # "
0990 << theCurrentTemp.enty[i].runnum << ENDL;
0991 return false;
0992 }
0993 }
0994
0995 for (j = 0; j < 4; ++j) {
0996 db >> theCurrentTemp.enty[i].xavgc2m[j] >> theCurrentTemp.enty[i].xrmsc2m[j] >>
0997 theCurrentTemp.enty[i].chi2xavgc2m[j] >> theCurrentTemp.enty[i].chi2xminc2m[j];
0998
0999 if (db.fail()) {
1000 LOGERROR("SiPixelTemplate") << "Error reading file 14, no template load, run # "
1001 << theCurrentTemp.enty[i].runnum << ENDL;
1002 return false;
1003 }
1004 }
1005
1006 for (j = 0; j < 4; ++j) {
1007 db >> theCurrentTemp.enty[i].yavggen[j] >> theCurrentTemp.enty[i].yrmsgen[j] >>
1008 theCurrentTemp.enty[i].ygx0gen[j] >> theCurrentTemp.enty[i].ygsiggen[j];
1009
1010 if (db.fail()) {
1011 LOGERROR("SiPixelTemplate") << "Error reading file 14a, no template load, run # "
1012 << theCurrentTemp.enty[i].runnum << ENDL;
1013 return false;
1014 }
1015 }
1016
1017 for (j = 0; j < 4; ++j) {
1018 db >> theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j] >>
1019 theCurrentTemp.enty[i].xgx0gen[j] >> theCurrentTemp.enty[i].xgsiggen[j];
1020
1021 if (db.fail()) {
1022 LOGERROR("SiPixelTemplate") << "Error reading file 14b, no template load, run # "
1023 << theCurrentTemp.enty[i].runnum << ENDL;
1024 return false;
1025 }
1026 }
1027
1028 db >> theCurrentTemp.enty[i].chi2yavgone >> theCurrentTemp.enty[i].chi2yminone >>
1029 theCurrentTemp.enty[i].chi2xavgone >> theCurrentTemp.enty[i].chi2xminone >> theCurrentTemp.enty[i].qmin2 >>
1030 theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav >>
1031 theCurrentTemp.enty[i].r_qMeas_qTrue >> theCurrentTemp.enty[i].spare[0];
1032
1033 if (db.fail()) {
1034 LOGERROR("SiPixelTemplate") << "Error reading file 15, no template load, run # "
1035 << theCurrentTemp.enty[i].runnum << ENDL;
1036 return false;
1037 }
1038
1039 db >> theCurrentTemp.enty[i].mpvvav2 >> theCurrentTemp.enty[i].sigmavav2 >> theCurrentTemp.enty[i].kappavav2 >>
1040 theCurrentTemp.enty[i].qbfrac[0] >> theCurrentTemp.enty[i].qbfrac[1] >> theCurrentTemp.enty[i].qbfrac[2] >>
1041 theCurrentTemp.enty[i].fracyone >> theCurrentTemp.enty[i].fracxone >> theCurrentTemp.enty[i].fracytwo >>
1042 theCurrentTemp.enty[i].fracxtwo;
1043
1044
1045 if (db.fail()) {
1046 LOGERROR("SiPixelTemplate") << "Error reading file 16, no template load, run # "
1047 << theCurrentTemp.enty[i].runnum << ENDL;
1048 return false;
1049 }
1050 }
1051
1052
1053
1054 for (k = 0; k < theCurrentTemp.head.NTyx; ++k) {
1055 for (i = 0; i < theCurrentTemp.head.NTxx; ++i) {
1056 db >> theCurrentTemp.entx[k][i].runnum >> theCurrentTemp.entx[k][i].costrk[0] >>
1057 theCurrentTemp.entx[k][i].costrk[1] >> theCurrentTemp.entx[k][i].costrk[2];
1058
1059 if (db.fail()) {
1060 LOGERROR("SiPixelTemplate") << "Error reading file 17, no template load, run # "
1061 << theCurrentTemp.entx[k][i].runnum << ENDL;
1062 return false;
1063 }
1064
1065
1066
1067 theCurrentTemp.entx[k][i].alpha = static_cast<float>(
1068 atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[0]));
1069
1070 theCurrentTemp.entx[k][i].cotalpha = theCurrentTemp.entx[k][i].costrk[0] / theCurrentTemp.entx[k][i].costrk[2];
1071
1072 theCurrentTemp.entx[k][i].beta = static_cast<float>(
1073 atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[1]));
1074
1075 theCurrentTemp.entx[k][i].cotbeta = theCurrentTemp.entx[k][i].costrk[1] / theCurrentTemp.entx[k][i].costrk[2];
1076
1077 db >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].pixmax >> theCurrentTemp.entx[k][i].symax >>
1078 theCurrentTemp.entx[k][i].dyone >> theCurrentTemp.entx[k][i].syone >> theCurrentTemp.entx[k][i].sxmax >>
1079 theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone;
1080
1081 if (db.fail()) {
1082 LOGERROR("SiPixelTemplate") << "Error reading file 18, no template load, run # "
1083 << theCurrentTemp.entx[k][i].runnum << ENDL;
1084 return false;
1085 }
1086
1087 db >> theCurrentTemp.entx[k][i].dytwo >> theCurrentTemp.entx[k][i].sytwo >> theCurrentTemp.entx[k][i].dxtwo >>
1088 theCurrentTemp.entx[k][i].sxtwo >> theCurrentTemp.entx[k][i].qmin >> theCurrentTemp.entx[k][i].clsleny >>
1089 theCurrentTemp.entx[k][i].clslenx;
1090
1091
1092 if (db.fail()) {
1093 LOGERROR("SiPixelTemplate") << "Error reading file 19, no template load, run # "
1094 << theCurrentTemp.entx[k][i].runnum << ENDL;
1095 return false;
1096 }
1097
1098 for (j = 0; j < 2; ++j) {
1099 db >> theCurrentTemp.entx[k][i].ypar[j][0] >> theCurrentTemp.entx[k][i].ypar[j][1] >>
1100 theCurrentTemp.entx[k][i].ypar[j][2] >> theCurrentTemp.entx[k][i].ypar[j][3] >>
1101 theCurrentTemp.entx[k][i].ypar[j][4];
1102
1103 if (db.fail()) {
1104 LOGERROR("SiPixelTemplate") << "Error reading file 20, no template load, run # "
1105 << theCurrentTemp.entx[k][i].runnum << ENDL;
1106 return false;
1107 }
1108 }
1109
1110 for (j = 0; j < 9; ++j) {
1111 for (l = 0; l < TYSIZE; ++l) {
1112 db >> theCurrentTemp.entx[k][i].ytemp[j][l];
1113 }
1114
1115 if (db.fail()) {
1116 LOGERROR("SiPixelTemplate") << "Error reading file 21, no template load, run # "
1117 << theCurrentTemp.entx[k][i].runnum << ENDL;
1118 return false;
1119 }
1120 }
1121
1122 for (j = 0; j < 2; ++j) {
1123 db >> theCurrentTemp.entx[k][i].xpar[j][0] >> theCurrentTemp.entx[k][i].xpar[j][1] >>
1124 theCurrentTemp.entx[k][i].xpar[j][2] >> theCurrentTemp.entx[k][i].xpar[j][3] >>
1125 theCurrentTemp.entx[k][i].xpar[j][4];
1126
1127 if (db.fail()) {
1128 LOGERROR("SiPixelTemplate") << "Error reading file 22, no template load, run # "
1129 << theCurrentTemp.entx[k][i].runnum << ENDL;
1130 return false;
1131 }
1132 }
1133
1134 qavg_avg = 0.f;
1135 for (j = 0; j < 9; ++j) {
1136 for (l = 0; l < TXSIZE; ++l) {
1137 db >> theCurrentTemp.entx[k][i].xtemp[j][l];
1138 qavg_avg += theCurrentTemp.entx[k][i].xtemp[j][l];
1139 }
1140
1141 if (db.fail()) {
1142 LOGERROR("SiPixelTemplate") << "Error reading file 23, no template load, run # "
1143 << theCurrentTemp.entx[k][i].runnum << ENDL;
1144 return false;
1145 }
1146 }
1147 theCurrentTemp.entx[k][i].qavg_avg = qavg_avg / 9.;
1148
1149 for (j = 0; j < 4; ++j) {
1150 db >> theCurrentTemp.entx[k][i].yavg[j] >> theCurrentTemp.entx[k][i].yrms[j] >>
1151 theCurrentTemp.entx[k][i].ygx0[j] >> theCurrentTemp.entx[k][i].ygsig[j];
1152
1153 if (db.fail()) {
1154 LOGERROR("SiPixelTemplate") << "Error reading file 24, no template load, run # "
1155 << theCurrentTemp.entx[k][i].runnum << ENDL;
1156 return false;
1157 }
1158 }
1159
1160 for (j = 0; j < 4; ++j) {
1161 db >> theCurrentTemp.entx[k][i].yflpar[j][0] >> theCurrentTemp.entx[k][i].yflpar[j][1] >>
1162 theCurrentTemp.entx[k][i].yflpar[j][2] >> theCurrentTemp.entx[k][i].yflpar[j][3] >>
1163 theCurrentTemp.entx[k][i].yflpar[j][4] >> theCurrentTemp.entx[k][i].yflpar[j][5];
1164
1165 if (db.fail()) {
1166 LOGERROR("SiPixelTemplate") << "Error reading file 25, no template load, run # "
1167 << theCurrentTemp.entx[k][i].runnum << ENDL;
1168 return false;
1169 }
1170 }
1171
1172 for (j = 0; j < 4; ++j) {
1173 db >> theCurrentTemp.entx[k][i].xavg[j] >> theCurrentTemp.entx[k][i].xrms[j] >>
1174 theCurrentTemp.entx[k][i].xgx0[j] >> theCurrentTemp.entx[k][i].xgsig[j];
1175
1176 if (db.fail()) {
1177 LOGERROR("SiPixelTemplate") << "Error reading file 26, no template load, run # "
1178 << theCurrentTemp.entx[k][i].runnum << ENDL;
1179 return false;
1180 }
1181 }
1182
1183 for (j = 0; j < 4; ++j) {
1184 db >> theCurrentTemp.entx[k][i].xflpar[j][0] >> theCurrentTemp.entx[k][i].xflpar[j][1] >>
1185 theCurrentTemp.entx[k][i].xflpar[j][2] >> theCurrentTemp.entx[k][i].xflpar[j][3] >>
1186 theCurrentTemp.entx[k][i].xflpar[j][4] >> theCurrentTemp.entx[k][i].xflpar[j][5];
1187
1188 if (db.fail()) {
1189 LOGERROR("SiPixelTemplate") << "Error reading file 27, no template load, run # "
1190 << theCurrentTemp.entx[k][i].runnum << ENDL;
1191 return false;
1192 }
1193 }
1194
1195 for (j = 0; j < 4; ++j) {
1196 db >> theCurrentTemp.entx[k][i].chi2yavg[j] >> theCurrentTemp.entx[k][i].chi2ymin[j] >>
1197 theCurrentTemp.entx[k][i].chi2xavg[j] >> theCurrentTemp.entx[k][i].chi2xmin[j];
1198
1199 if (db.fail()) {
1200 LOGERROR("SiPixelTemplate") << "Error reading file 28, no template load, run # "
1201 << theCurrentTemp.entx[k][i].runnum << ENDL;
1202 return false;
1203 }
1204 }
1205
1206 for (j = 0; j < 4; ++j) {
1207 db >> theCurrentTemp.entx[k][i].yavgc2m[j] >> theCurrentTemp.entx[k][i].yrmsc2m[j] >>
1208 theCurrentTemp.entx[k][i].chi2yavgc2m[j] >> theCurrentTemp.entx[k][i].chi2yminc2m[j];
1209
1210 if (db.fail()) {
1211 LOGERROR("SiPixelTemplate") << "Error reading file 29, no template load, run # "
1212 << theCurrentTemp.entx[k][i].runnum << ENDL;
1213 return false;
1214 }
1215 }
1216
1217 for (j = 0; j < 4; ++j) {
1218 db >> theCurrentTemp.entx[k][i].xavgc2m[j] >> theCurrentTemp.entx[k][i].xrmsc2m[j] >>
1219 theCurrentTemp.entx[k][i].chi2xavgc2m[j] >> theCurrentTemp.entx[k][i].chi2xminc2m[j];
1220
1221 if (db.fail()) {
1222 LOGERROR("SiPixelTemplate") << "Error reading file 30, no template load, run # "
1223 << theCurrentTemp.entx[k][i].runnum << ENDL;
1224 return false;
1225 }
1226 }
1227
1228 for (j = 0; j < 4; ++j) {
1229 db >> theCurrentTemp.entx[k][i].yavggen[j] >> theCurrentTemp.entx[k][i].yrmsgen[j] >>
1230 theCurrentTemp.entx[k][i].ygx0gen[j] >> theCurrentTemp.entx[k][i].ygsiggen[j];
1231
1232 if (db.fail()) {
1233 LOGERROR("SiPixelTemplate") << "Error reading file 30a, no template load, run # "
1234 << theCurrentTemp.entx[k][i].runnum << ENDL;
1235 return false;
1236 }
1237 }
1238
1239 for (j = 0; j < 4; ++j) {
1240 db >> theCurrentTemp.entx[k][i].xavggen[j] >> theCurrentTemp.entx[k][i].xrmsgen[j] >>
1241 theCurrentTemp.entx[k][i].xgx0gen[j] >> theCurrentTemp.entx[k][i].xgsiggen[j];
1242
1243 if (db.fail()) {
1244 LOGERROR("SiPixelTemplate") << "Error reading file 30b, no template load, run # "
1245 << theCurrentTemp.entx[k][i].runnum << ENDL;
1246 return false;
1247 }
1248 }
1249
1250 db >> theCurrentTemp.entx[k][i].chi2yavgone >> theCurrentTemp.entx[k][i].chi2yminone >>
1251 theCurrentTemp.entx[k][i].chi2xavgone >> theCurrentTemp.entx[k][i].chi2xminone >>
1252 theCurrentTemp.entx[k][i].qmin2 >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >>
1253 theCurrentTemp.entx[k][i].kappavav >> theCurrentTemp.entx[k][i].r_qMeas_qTrue >>
1254 theCurrentTemp.entx[k][i].spare[0];
1255
1256 if (db.fail()) {
1257 LOGERROR("SiPixelTemplate") << "Error reading file 31, no template load, run # "
1258 << theCurrentTemp.entx[k][i].runnum << ENDL;
1259 return false;
1260 }
1261
1262 db >> theCurrentTemp.entx[k][i].mpvvav2 >> theCurrentTemp.entx[k][i].sigmavav2 >>
1263 theCurrentTemp.entx[k][i].kappavav2 >> theCurrentTemp.entx[k][i].qbfrac[0] >>
1264 theCurrentTemp.entx[k][i].qbfrac[1] >> theCurrentTemp.entx[k][i].qbfrac[2] >>
1265 theCurrentTemp.entx[k][i].fracyone >> theCurrentTemp.entx[k][i].fracxone >>
1266 theCurrentTemp.entx[k][i].fracytwo >> theCurrentTemp.entx[k][i].fracxtwo;
1267
1268
1269 if (db.fail()) {
1270 LOGERROR("SiPixelTemplate") << "Error reading file 32, no template load, run # "
1271 << theCurrentTemp.entx[k][i].runnum << ENDL;
1272 return false;
1273 }
1274 }
1275 }
1276
1277
1278
1279 pixelTemp.push_back(theCurrentTemp);
1280 }
1281 postInit(pixelTemp);
1282 return true;
1283
1284 }
1285
1286 #endif
1287
1288 void SiPixelTemplate::postInit(std::vector<SiPixelTemplateStore>& thePixelTemp_) {
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306 for (auto& templ : thePixelTemp_) {
1307 for (auto iy = 0; iy < templ.head.NTy; ++iy)
1308 templ.cotbetaY[iy] = templ.enty[iy].cotbeta;
1309 for (auto iy = 0; iy < templ.head.NTyx; ++iy)
1310 templ.cotbetaX[iy] = templ.entx[iy][0].cotbeta;
1311 for (auto ix = 0; ix < templ.head.NTxx; ++ix)
1312 templ.cotalphaX[ix] = templ.entx[0][ix].cotalpha;
1313 }
1314 }
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328 bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx) {
1329
1330
1331
1332 if (!edm::isFinite(cotalpha) || !edm::isFinite(cotbeta)) {
1333 success_ = false;
1334 return success_;
1335 }
1336
1337
1338 int i, j;
1339 int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, imidy, imaxx;
1340 float yxratio, xxratio, sxmax, qcorrect, qxtempcor, symax, chi2xavgone, chi2xminone, cota, cotb, cotalpha0, cotbeta0;
1341 float chi2xavg[4], chi2xmin[4], chi2xavgc2m[4], chi2xminc2m[4];
1342
1343
1344
1345 if (entry_sideloaded_ != nullptr) {
1346 success_ = true;
1347 return success_;
1348 }
1349
1350
1351 if (id != id_current_ || cotalpha != cota_current_ || cotbeta != cotb_current_) {
1352 cota_current_ = cotalpha;
1353 cotb_current_ = cotbeta;
1354 success_ = true;
1355
1356 if (id != id_current_) {
1357
1358
1359 index_id_ = -1;
1360 for (i = 0; i < (int)thePixelTemp_.size(); ++i) {
1361
1362
1363 if (id == thePixelTemp_[i].head.ID) {
1364 index_id_ = i;
1365 id_current_ = id;
1366
1367
1368
1369 dtype_ = thePixelTemp_[index_id_].head.Dtype;
1370
1371
1372
1373 qscale_ = thePixelTemp_[index_id_].head.qscale;
1374
1375
1376
1377 s50_ = thePixelTemp_[index_id_].head.s50;
1378
1379
1380
1381 for (j = 0; j < 3; ++j) {
1382 fbin_[j] = thePixelTemp_[index_id_].head.fbin[j];
1383 }
1384
1385
1386
1387
1388 xsize_ = thePixelTemp_[index_id_].head.xsize;
1389 ysize_ = thePixelTemp_[index_id_].head.ysize;
1390 zsize_ = thePixelTemp_[index_id_].head.zsize;
1391
1392 break;
1393 }
1394 }
1395 }
1396
1397 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1398 if (index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
1399 throw cms::Exception("DataCorrupt")
1400 << "SiPixelTemplate::interpolate can't find needed template ID = " << id << std::endl;
1401 }
1402 #else
1403 assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
1404 #endif
1405
1406
1407
1408 cotalpha0 = thePixelTemp_[index_id_].enty[0].cotalpha;
1409 qcorrect =
1410 std::sqrt((1.f + cotbeta * cotbeta + cotalpha * cotalpha) / (1.f + cotbeta * cotbeta + cotalpha0 * cotalpha0));
1411
1412
1413 cota = cotalpha;
1414 cotb = abs_cotb_ = std::abs(cotbeta);
1415 flip_x_ = false;
1416 flip_y_ = false;
1417 switch (dtype_) {
1418 case 0:
1419 if (cotbeta < 0.f) {
1420 flip_y_ = true;
1421 }
1422 break;
1423 case 1:
1424 if (locBz < 0.f) {
1425 cotb = cotbeta;
1426 } else {
1427 cotb = -cotbeta;
1428 flip_y_ = true;
1429 }
1430 break;
1431 case 2:
1432 case 3:
1433 case 4:
1434 case 5:
1435 if (locBx * locBz < 0.f) {
1436 cota = -cotalpha;
1437 flip_x_ = true;
1438 }
1439 if (locBx > 0.f) {
1440 cotb = cotbeta;
1441 } else {
1442 cotb = -cotbeta;
1443 flip_y_ = true;
1444 }
1445 break;
1446 default:
1447 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1448 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::illegal subdetector ID = " << dtype_ << std::endl;
1449 #else
1450 std::cout << "SiPixelTemplate::illegal subdetector ID = " << dtype_ << std::endl;
1451 #endif
1452 }
1453
1454 Ny = thePixelTemp_[index_id_].head.NTy;
1455 Nyx = thePixelTemp_[index_id_].head.NTyx;
1456 Nxx = thePixelTemp_[index_id_].head.NTxx;
1457
1458 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1459 if (Ny < 2 || Nyx < 1 || Nxx < 2) {
1460 throw cms::Exception("DataCorrupt")
1461 << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx
1462 << std::endl;
1463 }
1464 #else
1465 assert(Ny > 1 && Nyx > 0 && Nxx > 1);
1466 #endif
1467 imaxx = Nyx - 1;
1468 imidy = Nxx / 2;
1469
1470
1471
1472 ilow = 0;
1473 yratio_ = 0.f;
1474
1475 if (cotb >= thePixelTemp_[index_id_].enty[Ny - 1].cotbeta) {
1476 ilow = Ny - 2;
1477 yratio_ = 1.;
1478 success_ = false;
1479
1480 } else {
1481 if (cotb >= thePixelTemp_[index_id_].enty[0].cotbeta) {
1482 for (i = 0; i < Ny - 1; ++i) {
1483 if (thePixelTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index_id_].enty[i + 1].cotbeta) {
1484 ilow = i;
1485 yratio_ = (cotb - thePixelTemp_[index_id_].enty[i].cotbeta) /
1486 (thePixelTemp_[index_id_].enty[i + 1].cotbeta - thePixelTemp_[index_id_].enty[i].cotbeta);
1487 break;
1488 }
1489 }
1490 } else {
1491 success_ = false;
1492 }
1493 }
1494
1495 ihigh = ilow + 1;
1496
1497
1498
1499 enty0_ = &thePixelTemp_[index_id_].enty[ilow];
1500 enty1_ = &thePixelTemp_[index_id_].enty[ihigh];
1501
1502
1503
1504 qavg_ = (1.f - yratio_) * enty0_->qavg + yratio_ * enty1_->qavg;
1505 qavg_ *= qcorrect;
1506 symax = (1.f - yratio_) * enty0_->symax + yratio_ * enty1_->symax;
1507 syparmax_ = symax;
1508 sxmax = (1.f - yratio_) * enty0_->sxmax + yratio_ * enty1_->sxmax;
1509 dyone_ = (1.f - yratio_) * enty0_->dyone + yratio_ * enty1_->dyone;
1510 if (flip_y_) {
1511 dyone_ = -dyone_;
1512 }
1513 syone_ = (1.f - yratio_) * enty0_->syone + yratio_ * enty1_->syone;
1514 dytwo_ = (1.f - yratio_) * enty0_->dytwo + yratio_ * enty1_->dytwo;
1515 if (flip_y_) {
1516 dytwo_ = -dytwo_;
1517 }
1518 sytwo_ = (1.f - yratio_) * enty0_->sytwo + yratio_ * enty1_->sytwo;
1519 qmin_ = (1.f - yratio_) * enty0_->qmin + yratio_ * enty1_->qmin;
1520 qmin_ *= qcorrect;
1521 qmin2_ = (1.f - yratio_) * enty0_->qmin2 + yratio_ * enty1_->qmin2;
1522 qmin2_ *= qcorrect;
1523 mpvvav_ = (1.f - yratio_) * enty0_->mpvvav + yratio_ * enty1_->mpvvav;
1524 mpvvav_ *= qcorrect;
1525 sigmavav_ = (1.f - yratio_) * enty0_->sigmavav + yratio_ * enty1_->sigmavav;
1526 kappavav_ = (1.f - yratio_) * enty0_->kappavav + yratio_ * enty1_->kappavav;
1527 mpvvav2_ = (1.f - yratio_) * enty0_->mpvvav2 + yratio_ * enty1_->mpvvav2;
1528 mpvvav2_ *= qcorrect;
1529 sigmavav2_ = (1.f - yratio_) * enty0_->sigmavav2 + yratio_ * enty1_->sigmavav2;
1530 kappavav2_ = (1.f - yratio_) * enty0_->kappavav2 + yratio_ * enty1_->kappavav2;
1531 clsleny_ = fminf(enty0_->clsleny, enty1_->clsleny);
1532 qavg_avg_ = (1.f - yratio_) * enty0_->qavg_avg + yratio_ * enty1_->qavg_avg;
1533 qavg_avg_ *= qcorrect;
1534 for (i = 0; i < 2; ++i) {
1535 for (j = 0; j < 5; ++j) {
1536
1537
1538 if (flip_y_) {
1539 yparl_[1 - i][j] = enty0_->ypar[i][j];
1540 yparh_[1 - i][j] = enty1_->ypar[i][j];
1541 } else {
1542 yparl_[i][j] = enty0_->ypar[i][j];
1543 yparh_[i][j] = enty1_->ypar[i][j];
1544 }
1545 if (flip_x_) {
1546 xparly0_[1 - i][j] = enty0_->xpar[i][j];
1547 xparhy0_[1 - i][j] = enty1_->xpar[i][j];
1548 } else {
1549 xparly0_[i][j] = enty0_->xpar[i][j];
1550 xparhy0_[i][j] = enty1_->xpar[i][j];
1551 }
1552 }
1553 }
1554
1555 for (i = 0; i < 4; ++i) {
1556 yavg_[i] = (1.f - yratio_) * thePixelTemp_[index_id_].enty[ilow].yavg[i] +
1557 yratio_ * thePixelTemp_[index_id_].enty[ihigh].yavg[i];
1558 if (flip_y_) {
1559 yavg_[i] = -yavg_[i];
1560 }
1561 yavg_[i] = (1.f - yratio_) * enty0_->yavg[i] + yratio_ * enty1_->yavg[i];
1562 if (flip_y_) {
1563 yavg_[i] = -yavg_[i];
1564 }
1565 yrms_[i] = (1.f - yratio_) * enty0_->yrms[i] + yratio_ * enty1_->yrms[i];
1566 chi2yavg_[i] = (1.f - yratio_) * enty0_->chi2yavg[i] + yratio_ * enty1_->chi2yavg[i];
1567 chi2ymin_[i] = (1.f - yratio_) * enty0_->chi2ymin[i] + yratio_ * enty1_->chi2ymin[i];
1568 chi2xavg[i] = (1.f - yratio_) * enty0_->chi2xavg[i] + yratio_ * enty1_->chi2xavg[i];
1569 chi2xmin[i] = (1.f - yratio_) * enty0_->chi2xmin[i] + yratio_ * enty1_->chi2xmin[i];
1570 yavgc2m_[i] = (1.f - yratio_) * enty0_->yavgc2m[i] + yratio_ * enty1_->yavgc2m[i];
1571 if (flip_y_) {
1572 yavgc2m_[i] = -yavgc2m_[i];
1573 }
1574 yrmsc2m_[i] = (1.f - yratio_) * enty0_->yrmsc2m[i] + yratio_ * enty1_->yrmsc2m[i];
1575 chi2yavgc2m_[i] = (1.f - yratio_) * enty0_->chi2yavgc2m[i] + yratio_ * enty1_->chi2yavgc2m[i];
1576
1577 chi2yminc2m_[i] = (1.f - yratio_) * enty0_->chi2yminc2m[i] + yratio_ * enty1_->chi2yminc2m[i];
1578
1579 chi2xavgc2m[i] = (1.f - yratio_) * enty0_->chi2xavgc2m[i] + yratio_ * enty1_->chi2xavgc2m[i];
1580 chi2xminc2m[i] = (1.f - yratio_) * enty0_->chi2xminc2m[i] + yratio_ * enty1_->chi2xminc2m[i];
1581 for (j = 0; j < 6; ++j) {
1582 yflparl_[i][j] = enty0_->yflpar[i][j];
1583 yflparh_[i][j] = enty1_->yflpar[i][j];
1584
1585
1586
1587 if (flip_y_ && (j == 0 || j == 2 || j == 4)) {
1588 yflparl_[i][j] = -yflparl_[i][j];
1589 yflparh_[i][j] = -yflparh_[i][j];
1590 }
1591 }
1592 }
1593
1594
1595
1596 chi2yavgone_ = (1.f - yratio_) * enty0_->chi2yavgone + yratio_ * enty1_->chi2yavgone;
1597 chi2yminone_ = (1.f - yratio_) * enty0_->chi2yminone + yratio_ * enty1_->chi2yminone;
1598 chi2xavgone = (1.f - yratio_) * enty0_->chi2xavgone + yratio_ * enty1_->chi2xavgone;
1599 chi2xminone = (1.f - yratio_) * enty0_->chi2xminone + yratio_ * enty1_->chi2xminone;
1600
1601 fracyone_ = (1.f - yratio_) * enty0_->fracyone + yratio_ * enty1_->fracyone;
1602 fracytwo_ = (1.f - yratio_) * enty0_->fracytwo + yratio_ * enty1_->fracytwo;
1603
1604
1605
1606
1607
1608
1609
1610 for (i = 0; i < 9; ++i) {
1611 ytemp_[i][0] = 0.f;
1612 ytemp_[i][1] = 0.f;
1613 ytemp_[i][BYM2] = 0.f;
1614 ytemp_[i][BYM1] = 0.f;
1615 for (j = 0; j < TYSIZE; ++j) {
1616
1617
1618 if (flip_y_) {
1619 ytemp_[8 - i][BYM3 - j] = (1.f - yratio_) * enty0_->ytemp[i][j] + yratio_ * enty1_->ytemp[i][j];
1620 } else {
1621 ytemp_[i][j + 2] = (1.f - yratio_) * enty0_->ytemp[i][j] + yratio_ * enty1_->ytemp[i][j];
1622 }
1623 }
1624 }
1625
1626
1627
1628 iylow = 0;
1629 yxratio = 0.f;
1630
1631 if (abs_cotb_ >= thePixelTemp_[index_id_].entx[Nyx - 1][0].cotbeta) {
1632 iylow = Nyx - 2;
1633 yxratio = 1.f;
1634
1635 } else if (abs_cotb_ >= thePixelTemp_[index_id_].entx[0][0].cotbeta) {
1636 for (i = 0; i < Nyx - 1; ++i) {
1637 if (thePixelTemp_[index_id_].entx[i][0].cotbeta <= abs_cotb_ &&
1638 abs_cotb_ < thePixelTemp_[index_id_].entx[i + 1][0].cotbeta) {
1639 iylow = i;
1640 yxratio = (abs_cotb_ - thePixelTemp_[index_id_].entx[i][0].cotbeta) /
1641 (thePixelTemp_[index_id_].entx[i + 1][0].cotbeta - thePixelTemp_[index_id_].entx[i][0].cotbeta);
1642 break;
1643 }
1644 }
1645 }
1646
1647 iyhigh = iylow + 1;
1648
1649 ilow = 0;
1650 xxratio = 0.f;
1651
1652 if (cota >= thePixelTemp_[index_id_].entx[0][Nxx - 1].cotalpha) {
1653 ilow = Nxx - 2;
1654 xxratio = 1.f;
1655 success_ = false;
1656
1657 } else {
1658 if (cota >= thePixelTemp_[index_id_].entx[0][0].cotalpha) {
1659 for (i = 0; i < Nxx - 1; ++i) {
1660 if (thePixelTemp_[index_id_].entx[0][i].cotalpha <= cota &&
1661 cota < thePixelTemp_[index_id_].entx[0][i + 1].cotalpha) {
1662 ilow = i;
1663 xxratio = (cota - thePixelTemp_[index_id_].entx[0][i].cotalpha) /
1664 (thePixelTemp_[index_id_].entx[0][i + 1].cotalpha - thePixelTemp_[index_id_].entx[0][i].cotalpha);
1665 break;
1666 }
1667 }
1668 } else {
1669 success_ = false;
1670 }
1671 }
1672
1673 ihigh = ilow + 1;
1674
1675
1676
1677 yxratio_ = yxratio;
1678 xxratio_ = xxratio;
1679
1680
1681
1682 sxparmax_ = (1.f - xxratio) * thePixelTemp_[index_id_].entx[imaxx][ilow].sxmax +
1683 xxratio * thePixelTemp_[index_id_].entx[imaxx][ihigh].sxmax;
1684 sxmax_ = sxparmax_;
1685 if (thePixelTemp_[index_id_].entx[imaxx][imidy].sxmax != 0.f) {
1686 sxmax_ = sxmax_ / thePixelTemp_[index_id_].entx[imaxx][imidy].sxmax * sxmax;
1687 }
1688 symax_ = (1.f - xxratio) * thePixelTemp_[index_id_].entx[imaxx][ilow].symax +
1689 xxratio * thePixelTemp_[index_id_].entx[imaxx][ihigh].symax;
1690 if (thePixelTemp_[index_id_].entx[imaxx][imidy].symax != 0.f) {
1691 symax_ = symax_ / thePixelTemp_[index_id_].entx[imaxx][imidy].symax * symax;
1692 }
1693 dxone_ = (1.f - xxratio) * thePixelTemp_[index_id_].entx[0][ilow].dxone +
1694 xxratio * thePixelTemp_[index_id_].entx[0][ihigh].dxone;
1695 if (flip_x_) {
1696 dxone_ = -dxone_;
1697 }
1698 sxone_ = (1.f - xxratio) * thePixelTemp_[index_id_].entx[0][ilow].sxone +
1699 xxratio * thePixelTemp_[index_id_].entx[0][ihigh].sxone;
1700 dxtwo_ = (1.f - xxratio) * thePixelTemp_[index_id_].entx[0][ilow].dxtwo +
1701 xxratio * thePixelTemp_[index_id_].entx[0][ihigh].dxtwo;
1702 if (flip_x_) {
1703 dxtwo_ = -dxtwo_;
1704 }
1705 sxtwo_ = (1.f - xxratio) * thePixelTemp_[index_id_].entx[0][ilow].sxtwo +
1706 xxratio * thePixelTemp_[index_id_].entx[0][ihigh].sxtwo;
1707 clslenx_ = fminf(thePixelTemp_[index_id_].entx[0][ilow].clslenx, thePixelTemp_[index_id_].entx[0][ihigh].clslenx);
1708
1709 for (i = 0; i < 2; ++i) {
1710 for (j = 0; j < 5; ++j) {
1711
1712 if (flip_x_) {
1713 xpar0_[1 - i][j] = thePixelTemp_[index_id_].entx[imaxx][imidy].xpar[i][j];
1714 xparl_[1 - i][j] = thePixelTemp_[index_id_].entx[imaxx][ilow].xpar[i][j];
1715 xparh_[1 - i][j] = thePixelTemp_[index_id_].entx[imaxx][ihigh].xpar[i][j];
1716 } else {
1717 xpar0_[i][j] = thePixelTemp_[index_id_].entx[imaxx][imidy].xpar[i][j];
1718 xparl_[i][j] = thePixelTemp_[index_id_].entx[imaxx][ilow].xpar[i][j];
1719 xparh_[i][j] = thePixelTemp_[index_id_].entx[imaxx][ihigh].xpar[i][j];
1720 }
1721 }
1722 }
1723
1724 entx00_ = &thePixelTemp_[index_id_].entx[iylow][ilow];
1725 entx02_ = &thePixelTemp_[index_id_].entx[iylow][ihigh];
1726 entx20_ = &thePixelTemp_[index_id_].entx[iyhigh][ilow];
1727 entx22_ = &thePixelTemp_[index_id_].entx[iyhigh][ihigh];
1728 entx21_ = &thePixelTemp_[index_id_].entx[iyhigh][imidy];
1729
1730
1731 pixmax_ = (1.f - yxratio) * ((1.f - xxratio) * entx00_->pixmax + xxratio * entx02_->pixmax) +
1732 yxratio * ((1.f - xxratio) * entx20_->pixmax + xxratio * entx22_->pixmax);
1733
1734 r_qMeas_qTrue_ = (1.f - yxratio) * ((1.f - xxratio) * entx00_->r_qMeas_qTrue + xxratio * entx02_->r_qMeas_qTrue) +
1735 yxratio * ((1.f - xxratio) * entx20_->r_qMeas_qTrue + xxratio * entx22_->r_qMeas_qTrue);
1736
1737 for (i = 0; i < 4; ++i) {
1738 xavg_[i] = (1.f - yxratio) * ((1.f - xxratio) * entx00_->xavg[i] + xxratio * entx02_->xavg[i]) +
1739 yxratio * ((1.f - xxratio) * entx20_->xavg[i] + xxratio * entx22_->xavg[i]);
1740 if (flip_x_) {
1741 xavg_[i] = -xavg_[i];
1742 }
1743
1744 xrms_[i] = (1.f - yxratio) * ((1.f - xxratio) * entx00_->xrms[i] + xxratio * entx02_->xrms[i]) +
1745 yxratio * ((1.f - xxratio) * entx20_->xrms[i] + xxratio * entx22_->xrms[i]);
1746
1747 xavgc2m_[i] = (1.f - yxratio) * ((1.f - xxratio) * entx00_->xavgc2m[i] + xxratio * entx02_->xavgc2m[i]) +
1748 yxratio * ((1.f - xxratio) * entx20_->xavgc2m[i] + xxratio * entx22_->xavgc2m[i]);
1749 if (flip_x_) {
1750 xavgc2m_[i] = -xavgc2m_[i];
1751 }
1752
1753 xrmsc2m_[i] = (1.f - yxratio) * ((1.f - xxratio) * entx00_->xrmsc2m[i] + xxratio * entx02_->xrmsc2m[i]) +
1754 yxratio * ((1.f - xxratio) * entx20_->xrmsc2m[i] + xxratio * entx22_->xrmsc2m[i]);
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771 chi2xavg_[i] = ((1.f - xxratio) * entx20_->chi2xavg[i] + xxratio * entx22_->chi2xavg[i]);
1772 if (entx21_->chi2xavg[i] != 0.f) {
1773 chi2xavg_[i] = chi2xavg_[i] / entx21_->chi2xavg[i] * chi2xavg[i];
1774 }
1775
1776 chi2xmin_[i] = ((1.f - xxratio) * entx20_->chi2xmin[i] + xxratio * entx22_->chi2xmin[i]);
1777 if (entx21_->chi2xmin[i] != 0.f) {
1778 chi2xmin_[i] = chi2xmin_[i] / entx21_->chi2xmin[i] * chi2xmin[i];
1779 }
1780
1781 chi2xavgc2m_[i] = ((1.f - xxratio) * entx20_->chi2xavgc2m[i] + xxratio * entx22_->chi2xavgc2m[i]);
1782 if (entx21_->chi2xavgc2m[i] != 0.f) {
1783 chi2xavgc2m_[i] = chi2xavgc2m_[i] / entx21_->chi2xavgc2m[i] * chi2xavgc2m[i];
1784 }
1785
1786 chi2xminc2m_[i] = ((1.f - xxratio) * entx20_->chi2xminc2m[i] + xxratio * entx22_->chi2xminc2m[i]);
1787 if (entx21_->chi2xminc2m[i] != 0.f) {
1788 chi2xminc2m_[i] = chi2xminc2m_[i] / entx21_->chi2xminc2m[i] * chi2xminc2m[i];
1789 }
1790
1791 for (j = 0; j < 6; ++j) {
1792 xflparll_[i][j] = entx00_->xflpar[i][j];
1793 xflparlh_[i][j] = entx02_->xflpar[i][j];
1794 xflparhl_[i][j] = entx20_->xflpar[i][j];
1795 xflparhh_[i][j] = entx22_->xflpar[i][j];
1796
1797 if (flip_x_ && (j == 0 || j == 2 || j == 4)) {
1798 xflparll_[i][j] = -xflparll_[i][j];
1799 xflparlh_[i][j] = -xflparlh_[i][j];
1800 xflparhl_[i][j] = -xflparhl_[i][j];
1801 xflparhh_[i][j] = -xflparhh_[i][j];
1802 }
1803 }
1804 }
1805
1806
1807
1808 chi2xavgone_ = ((1.f - xxratio) * entx20_->chi2xavgone + xxratio * entx22_->chi2xavgone);
1809 if (entx21_->chi2xavgone != 0.f) {
1810 chi2xavgone_ = chi2xavgone_ / entx21_->chi2xavgone * chi2xavgone;
1811 }
1812
1813 chi2xminone_ = ((1.f - xxratio) * entx20_->chi2xminone + xxratio * entx22_->chi2xminone);
1814 if (entx21_->chi2xminone != 0.f) {
1815 chi2xminone_ = chi2xminone_ / entx21_->chi2xminone * chi2xminone;
1816 }
1817
1818 fracxone_ = (1.f - yxratio) * ((1.f - xxratio) * entx00_->fracxone + xxratio * entx02_->fracxone) +
1819 yxratio * ((1.f - xxratio) * entx20_->fracxone + xxratio * entx22_->fracxone);
1820 fracxtwo_ = (1.f - yxratio) * ((1.f - xxratio) * entx00_->fracxtwo + xxratio * entx02_->fracxtwo) +
1821 yxratio * ((1.f - xxratio) * entx20_->fracxtwo + xxratio * entx22_->fracxtwo);
1822
1823
1824
1825
1826
1827
1828
1829
1830
1831
1832
1833 cotbeta0 = thePixelTemp_[index_id_].entx[iyhigh][0].cotbeta;
1834 qxtempcor =
1835 std::sqrt((1.f + cotbeta * cotbeta + cotalpha * cotalpha) / (1.f + cotbeta0 * cotbeta0 + cotalpha * cotalpha));
1836
1837 for (i = 0; i < 9; ++i) {
1838 xtemp_[i][0] = 0.f;
1839 xtemp_[i][1] = 0.f;
1840 xtemp_[i][BXM2] = 0.f;
1841 xtemp_[i][BXM1] = 0.f;
1842 for (j = 0; j < TXSIZE; ++j) {
1843
1844
1845
1846 if (flip_x_) {
1847 xtemp_[8 - i][BXM3 - j] =
1848 qxtempcor * ((1.f - xxratio) * entx20_->xtemp[i][j] + xxratio * entx22_->xtemp[i][j]);
1849 } else {
1850 xtemp_[i][j + 2] = qxtempcor * ((1.f - xxratio) * entx20_->xtemp[i][j] + xxratio * entx22_->xtemp[i][j]);
1851 }
1852 }
1853 }
1854
1855 lorywidth_ = thePixelTemp_[index_id_].head.lorywidth;
1856 lorxwidth_ = thePixelTemp_[index_id_].head.lorxwidth;
1857 lorybias_ = thePixelTemp_[index_id_].head.lorybias;
1858 lorxbias_ = thePixelTemp_[index_id_].head.lorxbias;
1859 if (flip_x_) {
1860 lorxwidth_ = -lorxwidth_;
1861 lorxbias_ = -lorxbias_;
1862 }
1863 if (flip_y_) {
1864 lorywidth_ = -lorywidth_;
1865 lorybias_ = -lorybias_;
1866 }
1867 }
1868
1869 return success_;
1870 }
1871
1872
1873
1874
1875
1876
1877
1878
1879 bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta) {
1880
1881
1882
1883 float locBx = 1.f;
1884 if (cotbeta < 0.f) {
1885 locBx = -1.f;
1886 }
1887 float locBz = locBx;
1888 if (cotalpha < 0.f) {
1889 locBz = -locBx;
1890 }
1891 return SiPixelTemplate::interpolate(id, cotalpha, cotbeta, locBz, locBx);
1892 }
1893
1894
1895
1896
1897
1898
1899
1900
1901 bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta, float locBz) {
1902
1903
1904
1905 float locBx = 1.f;
1906 return SiPixelTemplate::interpolate(id, cotalpha, cotbeta, locBz, locBx);
1907 }
1908
1909
1910
1911
1912
1913
1914
1915
1916 #ifdef SI_PIXEL_TEMPLATE_STANDALONE
1917 void SiPixelTemplate::sideload(SiPixelTemplateEntry* entry,
1918 int iDtype,
1919 float locBx,
1920 float locBz,
1921 float lorwdy,
1922 float lorwdx,
1923 float q50,
1924 float fbin[3],
1925 float xsize,
1926 float ysize,
1927 float zsize) {
1928
1929 entry_sideloaded_ = entry;
1930
1931 enty1_ = entry;
1932 enty0_ = entry;
1933
1934 entx00_ = entry;
1935 entx02_ = entry;
1936 entx20_ = entry;
1937 entx22_ = entry;
1938 entx21_ = entry;
1939
1940 dtype_ = iDtype;
1941 lorywidth_ = lorwdy;
1942 lorxwidth_ = lorwdx;
1943 xsize_ = xsize;
1944 ysize_ = ysize;
1945 zsize_ = zsize;
1946 s50_ = q50;
1947 qscale_ = 1.f;
1948 for (int i = 0; i < 3; ++i) {
1949 fbin_[i] = fbin[i];
1950 }
1951
1952
1953
1954 yratio_ = 0.f;
1955 yxratio_ = 0.f;
1956 xxratio_ = 0.f;
1957
1958 qavg_ = entry->qavg;
1959 qmin_ = 0.f;
1960 qmin2_ = 0.f;
1961
1962 pixmax_ = entry->pixmax;
1963 sxmax_ = entry->sxmax;
1964 symax_ = entry->symax;
1965 clsleny_ = entry->clsleny;
1966 clslenx_ = entry->clslenx;
1967
1968 scaleyavg_ = 1.f;
1969 scalexavg_ = 1.f;
1970 delyavg_ = 0.f;
1971 delysig_ = 0.f;
1972
1973 dyone_ = entry->dyone;
1974 syone_ = entry->syone;
1975 dytwo_ = entry->dytwo;
1976 sytwo_ = entry->sytwo;
1977
1978 dxone_ = entry->dxone;
1979 sxone_ = entry->sxone;
1980 dxtwo_ = entry->dxtwo;
1981 sxtwo_ = entry->sxtwo;
1982
1983 chi2yminone_ = 0.f;
1984 chi2yavgone_ = 0.1;
1985 chi2xminone_ = 0.f;
1986 chi2xavgone_ = 0.1;
1987
1988 for (int i = 0; i < 4; ++i) {
1989 scalex_[i] = 1.f;
1990 scaley_[i] = 1.f;
1991 offsetx_[i] = 0.f;
1992 offsety_[i] = 0.f;
1993 xrms_[i] = 0.f;
1994 yrms_[i] = 0.f;
1995 xavg_[i] = 0.f;
1996 yavg_[i] = 0.f;
1997
1998 chi2yavg_[i] = 0.1;
1999 chi2xavg_[i] = 0.1;
2000
2001 chi2xmin_[i] = 0.f;
2002 chi2ymin_[i] = 0.f;
2003
2004 for (int j = 0; j < 6; j++) {
2005 yflparl_[i][j] = yflparh_[i][j] = entry->yflpar[i][j];
2006 xflparhh_[i][j] = xflparhl_[i][j] = xflparll_[i][j] = xflparlh_[i][j] = entry->xflpar[i][j];
2007 }
2008 }
2009
2010 sxparmax_ = entry->sxmax;
2011 syparmax_ = entry->symax;
2012
2013
2014
2015 flip_x_ = false;
2016 flip_y_ = false;
2017 float cotbeta = entry->cotbeta;
2018 switch (dtype_) {
2019 case 0:
2020 if (cotbeta < 0.f) {
2021 flip_y_ = true;
2022 }
2023 break;
2024 case 1:
2025 if (locBz > 0.f) {
2026 flip_y_ = true;
2027 }
2028 break;
2029 case 2:
2030 case 3:
2031 case 4:
2032 case 5:
2033 if (locBx * locBz < 0.f) {
2034 flip_x_ = true;
2035 }
2036 if (locBx < 0.f) {
2037 flip_y_ = true;
2038 }
2039 break;
2040 default:
2041 std::cout << "SiPixelTemplate:illegal subdetector ID = " << dtype_ << std::endl;
2042 }
2043
2044
2045
2046
2047
2048
2049 float qxtempcor = 1.f;
2050
2051 for (int i = 0; i < 9; ++i) {
2052 xtemp_[i][0] = 0.f;
2053 xtemp_[i][1] = 0.f;
2054 xtemp_[i][BXM2] = 0.f;
2055 xtemp_[i][BXM1] = 0.f;
2056 for (int j = 0; j < TXSIZE; ++j) {
2057 if (flip_x_) {
2058 xtemp_[8 - i][BXM3 - j] = qxtempcor * entry_sideloaded_->xtemp[i][j];
2059 } else {
2060 xtemp_[i][j + 2] = qxtempcor * entry_sideloaded_->xtemp[i][j];
2061 }
2062 }
2063 }
2064
2065 for (int i = 0; i < 9; ++i) {
2066 ytemp_[i][0] = 0.f;
2067 ytemp_[i][1] = 0.f;
2068 ytemp_[i][BYM2] = 0.f;
2069 ytemp_[i][BYM1] = 0.f;
2070 for (int j = 0; j < TYSIZE; ++j) {
2071
2072
2073 if (flip_y_) {
2074 ytemp_[8 - i][BYM3 - j] = entry_sideloaded_->ytemp[i][j];
2075 } else {
2076 ytemp_[i][j + 2] = entry_sideloaded_->ytemp[i][j];
2077 }
2078 }
2079 }
2080
2081
2082 for (int i = 0; i < 2; i++) {
2083 for (int j = 0; j < 5; j++) {
2084 if (flip_y_) {
2085 yparl_[1 - i][j] = yparh_[1 - i][j] = entry->ypar[i][j];
2086 } else {
2087 yparl_[i][j] = yparh_[i][j] = entry->ypar[i][j];
2088 }
2089
2090 if (flip_x_) {
2091 xpar0_[1 - i][j] = xparly0_[1 - i][j] = xparhy0_[1 - i][j] = xparl_[1 - i][j] = xparh_[1 - i][j] =
2092 entry->xpar[i][j];
2093 } else {
2094 xpar0_[i][j] = xparly0_[i][j] = xparhy0_[i][j] = xparl_[i][j] = xparh_[i][j] = entry->xpar[i][j];
2095 }
2096 }
2097 }
2098
2099 return;
2100 }
2101 #endif
2102
2103
2104
2105
2106
2107
2108
2109
2110
2111
2112 void SiPixelTemplate::ysigma2(int fypix, int lypix, float sythr, float ysum[25], float ysig2[25])
2113
2114 {
2115
2116
2117
2118 int i;
2119 float sigi, sigi2, sigi3, sigi4, symax, qscale, s25;
2120
2121
2122
2123 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2124 if (fypix < 2 || fypix >= BYM2) {
2125 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ysigma2 called with fypix = " << fypix << std::endl;
2126 }
2127 #else
2128 assert(fypix > 1 && fypix < BYM2);
2129 #endif
2130 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2131 if (lypix < fypix || lypix >= BYM2) {
2132 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ysigma2 called with lypix/fypix = " << lypix << "/"
2133 << fypix << std::endl;
2134 }
2135 #else
2136 assert(lypix >= fypix && lypix < BYM2);
2137 #endif
2138
2139
2140
2141 symax = symax_;
2142 s25 = 0.5f * s50_;
2143 if (symax_ > syparmax_) {
2144 symax = syparmax_;
2145 }
2146
2147
2148
2149 for (i = fypix - 2; i <= lypix + 2; ++i) {
2150 if (i < fypix || i > lypix) {
2151
2152
2153 ysig2[i] = s50_ * s50_;
2154 } else {
2155 if (ysum[i] < symax) {
2156 sigi = ysum[i];
2157 qscale = 1.f;
2158 if (sigi < s25)
2159 sigi = s25;
2160 } else {
2161 sigi = symax;
2162 qscale = ysum[i] / symax;
2163 }
2164 sigi2 = sigi * sigi;
2165 sigi3 = sigi2 * sigi;
2166 sigi4 = sigi3 * sigi;
2167 if (i <= BHY) {
2168 ysig2[i] = (1.f - yratio_) * (yparl_[0][0] + yparl_[0][1] * sigi + yparl_[0][2] * sigi2 + yparl_[0][3] * sigi3 +
2169 yparl_[0][4] * sigi4) +
2170 yratio_ * (yparh_[0][0] + yparh_[0][1] * sigi + yparh_[0][2] * sigi2 + yparh_[0][3] * sigi3 +
2171 yparh_[0][4] * sigi4);
2172 } else {
2173 ysig2[i] = (1.f - yratio_) * (yparl_[1][0] + yparl_[1][1] * sigi + yparl_[1][2] * sigi2 + yparl_[1][3] * sigi3 +
2174 yparl_[1][4] * sigi4) +
2175 yratio_ * (yparh_[1][0] + yparh_[1][1] * sigi + yparh_[1][2] * sigi2 + yparh_[1][3] * sigi3 +
2176 yparh_[1][4] * sigi4);
2177 }
2178 ysig2[i] *= qscale;
2179 if (ysum[i] > sythr) {
2180 ysig2[i] = 1.e8;
2181 }
2182 if (ysig2[i] <= 0.f) {
2183 LOGERROR("SiPixelTemplate") << "neg y-error-squared, id = " << id_current_ << ", index = " << index_id_
2184 << ", cot(alpha) = " << cota_current_ << ", cot(beta) = " << cotb_current_
2185 << ", sigi = " << sigi << ENDL;
2186 }
2187 }
2188 }
2189
2190 return;
2191
2192 }
2193
2194
2195
2196
2197
2198
2199
2200
2201 void SiPixelTemplate::ysigma2(float qpixel, int index, float& ysig2)
2202
2203 {
2204
2205
2206
2207 float sigi, sigi2, sigi3, sigi4, symax, qscale, err2, s25;
2208
2209
2210
2211 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2212 if (index < 2 || index >= BYM2) {
2213 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ysigma2 called with index = " << index << std::endl;
2214 }
2215 #else
2216 assert(index > 1 && index < BYM2);
2217 #endif
2218
2219
2220
2221 symax = symax_;
2222 s25 = 0.5f * s50_;
2223 if (symax_ > syparmax_) {
2224 symax = syparmax_;
2225 }
2226
2227
2228
2229 if (qpixel < symax) {
2230 sigi = qpixel;
2231 qscale = 1.f;
2232 if (sigi < s25)
2233 sigi = s25;
2234 } else {
2235 sigi = symax;
2236 qscale = qpixel / symax;
2237 }
2238 sigi2 = sigi * sigi;
2239 sigi3 = sigi2 * sigi;
2240 sigi4 = sigi3 * sigi;
2241 if (index <= BHY) {
2242 err2 =
2243 (1.f - yratio_) *
2244 (yparl_[0][0] + yparl_[0][1] * sigi + yparl_[0][2] * sigi2 + yparl_[0][3] * sigi3 + yparl_[0][4] * sigi4) +
2245 yratio_ *
2246 (yparh_[0][0] + yparh_[0][1] * sigi + yparh_[0][2] * sigi2 + yparh_[0][3] * sigi3 + yparh_[0][4] * sigi4);
2247 } else {
2248 err2 =
2249 (1.f - yratio_) *
2250 (yparl_[1][0] + yparl_[1][1] * sigi + yparl_[1][2] * sigi2 + yparl_[1][3] * sigi3 + yparl_[1][4] * sigi4) +
2251 yratio_ *
2252 (yparh_[1][0] + yparh_[1][1] * sigi + yparh_[1][2] * sigi2 + yparh_[1][3] * sigi3 + yparh_[1][4] * sigi4);
2253 }
2254 ysig2 = qscale * err2;
2255 if (ysig2 <= 0.f) {
2256 LOGERROR("SiPixelTemplate") << "neg y-error-squared, id = " << id_current_ << ", index = " << index_id_
2257 << ", cot(alpha) = " << cota_current_ << ", cot(beta) = " << cotb_current_
2258 << ", sigi = " << sigi << ENDL;
2259 }
2260
2261 return;
2262
2263 }
2264
2265
2266
2267
2268
2269
2270
2271
2272
2273
2274 void SiPixelTemplate::xsigma2(int fxpix, int lxpix, float sxthr, float xsum[BXSIZE], float xsig2[BXSIZE])
2275
2276 {
2277
2278
2279
2280 int i;
2281 float sigi, sigi2, sigi3, sigi4, yint, sxmax, x0, qscale, s25;
2282
2283
2284
2285 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2286 if (fxpix < 2 || fxpix >= BXM2) {
2287 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xsigma2 called with fxpix = " << fxpix << std::endl;
2288 }
2289 #else
2290 assert(fxpix > 1 && fxpix < BXM2);
2291 #endif
2292 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2293 if (lxpix < fxpix || lxpix >= BXM2) {
2294 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xsigma2 called with lxpix/fxpix = " << lxpix << "/"
2295 << fxpix << std::endl;
2296 }
2297 #else
2298 assert(lxpix >= fxpix && lxpix < BXM2);
2299 #endif
2300
2301
2302
2303 sxmax = sxmax_;
2304 s25 = 0.5f * s50_;
2305 if (sxmax_ > sxparmax_) {
2306 sxmax = sxparmax_;
2307 }
2308
2309
2310
2311 for (i = fxpix - 2; i <= lxpix + 2; ++i) {
2312 if (i < fxpix || i > lxpix) {
2313
2314
2315 xsig2[i] = s50_ * s50_;
2316 } else {
2317 if (xsum[i] < sxmax) {
2318 sigi = xsum[i];
2319 qscale = 1.f;
2320 if (sigi < s25)
2321 sigi = s25;
2322 } else {
2323 sigi = sxmax;
2324 qscale = xsum[i] / sxmax;
2325 }
2326 sigi2 = sigi * sigi;
2327 sigi3 = sigi2 * sigi;
2328 sigi4 = sigi3 * sigi;
2329
2330
2331
2332 if (i <= BHX) {
2333 yint = (1.f - yratio_) * (xparly0_[0][0] + xparly0_[0][1] * sigi + xparly0_[0][2] * sigi2 +
2334 xparly0_[0][3] * sigi3 + xparly0_[0][4] * sigi4) +
2335 yratio_ * (xparhy0_[0][0] + xparhy0_[0][1] * sigi + xparhy0_[0][2] * sigi2 + xparhy0_[0][3] * sigi3 +
2336 xparhy0_[0][4] * sigi4);
2337 } else {
2338 yint = (1.f - yratio_) * (xparly0_[1][0] + xparly0_[1][1] * sigi + xparly0_[1][2] * sigi2 +
2339 xparly0_[1][3] * sigi3 + xparly0_[1][4] * sigi4) +
2340 yratio_ * (xparhy0_[1][0] + xparhy0_[1][1] * sigi + xparhy0_[1][2] * sigi2 + xparhy0_[1][3] * sigi3 +
2341 xparhy0_[1][4] * sigi4);
2342 }
2343
2344
2345
2346 if (i <= BHX) {
2347 xsig2[i] = (1.f - xxratio_) * (xparl_[0][0] + xparl_[0][1] * sigi + xparl_[0][2] * sigi2 +
2348 xparl_[0][3] * sigi3 + xparl_[0][4] * sigi4) +
2349 xxratio_ * (xparh_[0][0] + xparh_[0][1] * sigi + xparh_[0][2] * sigi2 + xparh_[0][3] * sigi3 +
2350 xparh_[0][4] * sigi4);
2351 } else {
2352 xsig2[i] = (1.f - xxratio_) * (xparl_[1][0] + xparl_[1][1] * sigi + xparl_[1][2] * sigi2 +
2353 xparl_[1][3] * sigi3 + xparl_[1][4] * sigi4) +
2354 xxratio_ * (xparh_[1][0] + xparh_[1][1] * sigi + xparh_[1][2] * sigi2 + xparh_[1][3] * sigi3 +
2355 xparh_[1][4] * sigi4);
2356 }
2357
2358
2359
2360 if (i <= BHX) {
2361 x0 = xpar0_[0][0] + xpar0_[0][1] * sigi + xpar0_[0][2] * sigi2 + xpar0_[0][3] * sigi3 + xpar0_[0][4] * sigi4;
2362 } else {
2363 x0 = xpar0_[1][0] + xpar0_[1][1] * sigi + xpar0_[1][2] * sigi2 + xpar0_[1][3] * sigi3 + xpar0_[1][4] * sigi4;
2364 }
2365
2366
2367
2368 if (x0 != 0.f) {
2369 xsig2[i] = xsig2[i] / x0 * yint;
2370 }
2371 xsig2[i] *= qscale;
2372 if (xsum[i] > sxthr) {
2373 xsig2[i] = 1.e8;
2374 }
2375 if (xsig2[i] <= 0.f) {
2376 LOGERROR("SiPixelTemplate") << "neg x-error-squared, id = " << id_current_ << ", index = " << index_id_
2377 << ", cot(alpha) = " << cota_current_ << ", cot(beta) = " << cotb_current_
2378 << ", sigi = " << sigi << ENDL;
2379 }
2380 }
2381 }
2382
2383 return;
2384
2385 }
2386
2387
2388
2389
2390
2391
2392 float SiPixelTemplate::yflcorr(int binq, float qfly)
2393
2394 {
2395
2396
2397
2398 float qfl, qfl2, qfl3, qfl4, qfl5, dy;
2399
2400
2401
2402 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2403 if (binq < 0 || binq > 3) {
2404 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::yflcorr called with binq = " << binq << std::endl;
2405 }
2406 #else
2407 assert(binq >= 0 && binq < 4);
2408 #endif
2409 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2410 if (fabs((double)qfly) > 1.) {
2411 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::yflcorr called with qfly = " << qfly << std::endl;
2412 }
2413 #else
2414 assert(fabs((double)qfly) <= 1.);
2415 #endif
2416
2417
2418
2419 qfl = qfly;
2420
2421 if (qfl < -0.9f) {
2422 qfl = -0.9f;
2423 }
2424 if (qfl > 0.9f) {
2425 qfl = 0.9f;
2426 }
2427
2428
2429
2430 qfl2 = qfl * qfl;
2431 qfl3 = qfl2 * qfl;
2432 qfl4 = qfl3 * qfl;
2433 qfl5 = qfl4 * qfl;
2434 dy = (1.f - yratio_) * (yflparl_[binq][0] + yflparl_[binq][1] * qfl + yflparl_[binq][2] * qfl2 +
2435 yflparl_[binq][3] * qfl3 + yflparl_[binq][4] * qfl4 + yflparl_[binq][5] * qfl5) +
2436 yratio_ * (yflparh_[binq][0] + yflparh_[binq][1] * qfl + yflparh_[binq][2] * qfl2 + yflparh_[binq][3] * qfl3 +
2437 yflparh_[binq][4] * qfl4 + yflparh_[binq][5] * qfl5);
2438
2439 return dy;
2440
2441 }
2442
2443
2444
2445
2446
2447
2448 float SiPixelTemplate::xflcorr(int binq, float qflx)
2449
2450 {
2451
2452
2453
2454 float qfl, qfl2, qfl3, qfl4, qfl5, dx;
2455
2456
2457
2458 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2459 if (binq < 0 || binq > 3) {
2460 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xflcorr called with binq = " << binq << std::endl;
2461 }
2462 #else
2463 assert(binq >= 0 && binq < 4);
2464 #endif
2465 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2466 if (fabs((double)qflx) > 1.) {
2467 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xflcorr called with qflx = " << qflx << std::endl;
2468 }
2469 #else
2470 assert(fabs((double)qflx) <= 1.);
2471 #endif
2472
2473
2474
2475 qfl = qflx;
2476
2477 if (qfl < -0.9f) {
2478 qfl = -0.9f;
2479 }
2480 if (qfl > 0.9f) {
2481 qfl = 0.9f;
2482 }
2483
2484
2485
2486 qfl2 = qfl * qfl;
2487 qfl3 = qfl2 * qfl;
2488 qfl4 = qfl3 * qfl;
2489 qfl5 = qfl4 * qfl;
2490 dx = (1.f - yxratio_) *
2491 ((1.f - xxratio_) * (xflparll_[binq][0] + xflparll_[binq][1] * qfl + xflparll_[binq][2] * qfl2 +
2492 xflparll_[binq][3] * qfl3 + xflparll_[binq][4] * qfl4 + xflparll_[binq][5] * qfl5) +
2493 xxratio_ * (xflparlh_[binq][0] + xflparlh_[binq][1] * qfl + xflparlh_[binq][2] * qfl2 +
2494 xflparlh_[binq][3] * qfl3 + xflparlh_[binq][4] * qfl4 + xflparlh_[binq][5] * qfl5)) +
2495 yxratio_ *
2496 ((1.f - xxratio_) * (xflparhl_[binq][0] + xflparhl_[binq][1] * qfl + xflparhl_[binq][2] * qfl2 +
2497 xflparhl_[binq][3] * qfl3 + xflparhl_[binq][4] * qfl4 + xflparhl_[binq][5] * qfl5) +
2498 xxratio_ * (xflparhh_[binq][0] + xflparhh_[binq][1] * qfl + xflparhh_[binq][2] * qfl2 +
2499 xflparhh_[binq][3] * qfl3 + xflparhh_[binq][4] * qfl4 + xflparhh_[binq][5] * qfl5));
2500
2501 return dx;
2502
2503 }
2504
2505
2506
2507
2508
2509
2510
2511 void SiPixelTemplate::ytemp(int fybin, int lybin, float ytemplate[41][BYSIZE])
2512
2513 {
2514
2515
2516
2517 int i, j;
2518
2519
2520
2521 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2522 if (fybin < 0 || fybin > 40) {
2523 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ytemp called with fybin = " << fybin << std::endl;
2524 }
2525 #else
2526 assert(fybin >= 0 && fybin < 41);
2527 #endif
2528 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2529 if (lybin < 0 || lybin > 40) {
2530 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ytemp called with lybin = " << lybin << std::endl;
2531 }
2532 #else
2533 assert(lybin >= 0 && lybin < 41);
2534 #endif
2535
2536
2537
2538 for (i = 0; i < 9; ++i) {
2539 for (j = 0; j < BYSIZE; ++j) {
2540 ytemplate[i + 16][j] = ytemp_[i][j];
2541 }
2542 }
2543 for (i = 0; i < 8; ++i) {
2544 ytemplate[i + 8][BYM1] = 0.f;
2545 for (j = 0; j < BYM1; ++j) {
2546 ytemplate[i + 8][j] = ytemp_[i][j + 1];
2547 }
2548 }
2549 for (i = 1; i < 9; ++i) {
2550 ytemplate[i + 24][0] = 0.f;
2551 for (j = 0; j < BYM1; ++j) {
2552 ytemplate[i + 24][j + 1] = ytemp_[i][j];
2553 }
2554 }
2555
2556
2557
2558 if (fybin < 8) {
2559 for (i = 0; i < 8; ++i) {
2560 ytemplate[i][BYM2] = 0.f;
2561 ytemplate[i][BYM1] = 0.f;
2562 for (j = 0; j < BYM2; ++j) {
2563 ytemplate[i][j] = ytemp_[i][j + 2];
2564 }
2565 }
2566 }
2567 if (lybin > 32) {
2568 for (i = 1; i < 9; ++i) {
2569 ytemplate[i + 32][0] = 0.f;
2570 ytemplate[i + 32][1] = 0.f;
2571 for (j = 0; j < BYM2; ++j) {
2572 ytemplate[i + 32][j + 2] = ytemp_[i][j];
2573 }
2574 }
2575 }
2576
2577 return;
2578
2579 }
2580
2581
2582
2583
2584
2585
2586
2587 void SiPixelTemplate::xtemp(int fxbin, int lxbin, float xtemplate[41][BXSIZE]) {
2588
2589
2590
2591 int i, j;
2592
2593
2594
2595 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2596 if (fxbin < 0 || fxbin > 40) {
2597 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp called with fxbin = " << fxbin << std::endl;
2598 }
2599 #else
2600 assert(fxbin >= 0 && fxbin < 41);
2601 #endif
2602 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2603 if (lxbin < 0 || lxbin > 40) {
2604 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp called with lxbin = " << lxbin << std::endl;
2605 }
2606 #else
2607 assert(lxbin >= 0 && lxbin < 41);
2608 #endif
2609
2610
2611
2612 for (i = 0; i < 9; ++i) {
2613 for (j = 0; j < BXSIZE; ++j) {
2614 xtemplate[i + 16][j] = xtemp_[i][j];
2615 }
2616 }
2617 for (i = 0; i < 8; ++i) {
2618 xtemplate[i + 8][BXM1] = 0.f;
2619 for (j = 0; j < BXM1; ++j) {
2620 xtemplate[i + 8][j] = xtemp_[i][j + 1];
2621 }
2622 }
2623 for (i = 1; i < 9; ++i) {
2624 xtemplate[i + 24][0] = 0.f;
2625 for (j = 0; j < BXM1; ++j) {
2626 xtemplate[i + 24][j + 1] = xtemp_[i][j];
2627 }
2628 }
2629
2630 if (fxbin < 8) {
2631 for (i = 0; i < 8; ++i) {
2632 xtemplate[i][BXM2] = 0.f;
2633 xtemplate[i][BXM1] = 0.f;
2634 for (j = 0; j < BXM2; ++j) {
2635 xtemplate[i][j] = xtemp_[i][j + 2];
2636 }
2637 }
2638 }
2639 if (lxbin > 32) {
2640 for (i = 1; i < 9; ++i) {
2641 xtemplate[i + 32][0] = 0.f;
2642 xtemplate[i + 32][1] = 0.f;
2643 for (j = 0; j < BXM2; ++j) {
2644 xtemplate[i + 32][j + 2] = xtemp_[i][j];
2645 }
2646 }
2647 }
2648
2649 return;
2650
2651 }
2652
2653
2654
2655
2656 int SiPixelTemplate::cytemp()
2657
2658 {
2659
2660
2661
2662 int j;
2663
2664
2665
2666
2667 float sigmax = 0.f;
2668 float qedge = 2. * s50_;
2669 int jmax = -1;
2670
2671 for (j = 0; j < BYSIZE; ++j) {
2672 if (ytemp_[4][j] > sigmax) {
2673 sigmax = ytemp_[4][j];
2674 jmax = j;
2675 }
2676 }
2677 if (sigmax < qedge) {
2678 qedge = s50_;
2679 }
2680 if (sigmax < qedge || jmax < 1 || jmax > BYM2) {
2681 return -1;
2682 }
2683
2684
2685
2686 int jend = jmax;
2687
2688 for (j = jmax + 1; j < BYM1; ++j) {
2689 if (ytemp_[4][j] < qedge)
2690 break;
2691 jend = j;
2692 }
2693
2694 int jbeg = jmax;
2695
2696 for (j = jmax - 1; j > 0; --j) {
2697 if (ytemp_[4][j] < qedge)
2698 break;
2699 jbeg = j;
2700 }
2701
2702 return (jbeg + jend) / 2;
2703
2704 }
2705
2706
2707
2708
2709 int SiPixelTemplate::cxtemp()
2710
2711 {
2712
2713
2714
2715 int j;
2716
2717
2718
2719
2720 float sigmax = 0.f;
2721 float qedge = 2. * s50_;
2722 int jmax = -1;
2723
2724 for (j = 0; j < BXSIZE; ++j) {
2725 if (xtemp_[4][j] > sigmax) {
2726 sigmax = xtemp_[4][j];
2727 jmax = j;
2728 }
2729 }
2730 if (sigmax < qedge) {
2731 qedge = s50_;
2732 }
2733 if (sigmax < qedge || jmax < 1 || jmax > BXM2) {
2734 return -1;
2735 }
2736
2737
2738
2739 int jend = jmax;
2740
2741 for (j = jmax + 1; j < BXM1; ++j) {
2742 if (xtemp_[4][j] < qedge)
2743 break;
2744 jend = j;
2745 }
2746
2747 int jbeg = jmax;
2748
2749 for (j = jmax - 1; j > 0; --j) {
2750 if (xtemp_[4][j] < qedge)
2751 break;
2752 jbeg = j;
2753 }
2754
2755 return (jbeg + jend) / 2;
2756
2757 }
2758
2759
2760
2761
2762
2763
2764 void SiPixelTemplate::ytemp3d_int(int nypix, int& nybins)
2765
2766 {
2767
2768
2769
2770 int i, j, k;
2771 int ioff0, ioffp, ioffm;
2772
2773
2774
2775 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2776 if (nypix < 1 || nypix >= BYM3) {
2777 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ytemp3d called with nypix = " << nypix << std::endl;
2778 }
2779 #else
2780 assert(nypix > 0 && nypix < BYM3);
2781 #endif
2782
2783
2784
2785 float diff = fabsf(nypix - clsleny_) / 2. + 1.f;
2786 int nshift = (int)diff;
2787 if ((diff - nshift) > 0.5f) {
2788 ++nshift;
2789 }
2790
2791
2792
2793 nybins_ = 9 + 16 * nshift;
2794
2795
2796
2797 temp2dy_.resize(boost::extents[nybins_][BYSIZE]);
2798
2799
2800
2801 ioff0 = 8 * nshift;
2802
2803 for (i = 0; i < 9; ++i) {
2804 for (j = 0; j < BYSIZE; ++j) {
2805 temp2dy_[i + ioff0][j] = ytemp_[i][j];
2806 }
2807 }
2808
2809
2810
2811 for (k = 1; k <= nshift; ++k) {
2812 ioffm = ioff0 - k * 8;
2813 for (i = 0; i < 8; ++i) {
2814 for (j = 0; j < k; ++j) {
2815 temp2dy_[i + ioffm][BYM1 - j] = 0.f;
2816 }
2817 for (j = 0; j < BYSIZE - k; ++j) {
2818 temp2dy_[i + ioffm][j] = ytemp_[i][j + k];
2819 }
2820 }
2821 ioffp = ioff0 + k * 8;
2822 for (i = 1; i < 9; ++i) {
2823 for (j = 0; j < k; ++j) {
2824 temp2dy_[i + ioffp][j] = 0.f;
2825 }
2826 for (j = 0; j < BYSIZE - k; ++j) {
2827 temp2dy_[i + ioffp][j + k] = ytemp_[i][j];
2828 }
2829 }
2830 }
2831
2832 nybins = nybins_;
2833 return;
2834
2835 }
2836
2837
2838
2839
2840
2841
2842 void SiPixelTemplate::ytemp3d(int i, int j, std::vector<float>& ytemplate)
2843
2844 {
2845
2846 if (i >= 0 && i < nybins_ && j <= i) {
2847 for (int k = 0; k < BYSIZE; ++k) {
2848 ytemplate[k] = temp2dy_[i][k] + temp2dy_[j][k];
2849 }
2850 } else {
2851 for (int k = 0; k < BYSIZE; ++k) {
2852 ytemplate[k] = 0.;
2853 }
2854 }
2855
2856 return;
2857
2858 }
2859
2860
2861
2862
2863
2864
2865 void SiPixelTemplate::xtemp3d_int(int nxpix, int& nxbins)
2866
2867 {
2868
2869
2870
2871 int i, j, k;
2872 int ioff0, ioffp, ioffm;
2873
2874
2875
2876 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2877 if (nxpix < 1 || nxpix >= BXM3) {
2878 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp3d called with nxpix = " << nxpix << std::endl;
2879 }
2880 #else
2881 assert(nxpix > 0 && nxpix < BXM3);
2882 #endif
2883
2884
2885
2886 float diff = fabsf(nxpix - clslenx_) / 2.f + 1.f;
2887 int nshift = (int)diff;
2888 if ((diff - nshift) > 0.5f) {
2889 ++nshift;
2890 }
2891
2892
2893
2894 nxbins_ = 9 + 16 * nshift;
2895
2896
2897
2898 temp2dx_.resize(boost::extents[nxbins_][BXSIZE]);
2899
2900
2901
2902 ioff0 = 8 * nshift;
2903
2904 for (i = 0; i < 9; ++i) {
2905 for (j = 0; j < BXSIZE; ++j) {
2906 temp2dx_[i + ioff0][j] = xtemp_[i][j];
2907 }
2908 }
2909
2910
2911
2912 for (k = 1; k <= nshift; ++k) {
2913 ioffm = ioff0 - k * 8;
2914 for (i = 0; i < 8; ++i) {
2915 for (j = 0; j < k; ++j) {
2916 temp2dx_[i + ioffm][BXM1 - j] = 0.f;
2917 }
2918 for (j = 0; j < BXSIZE - k; ++j) {
2919 temp2dx_[i + ioffm][j] = xtemp_[i][j + k];
2920 }
2921 }
2922 ioffp = ioff0 + k * 8;
2923 for (i = 1; i < 9; ++i) {
2924 for (j = 0; j < k; ++j) {
2925 temp2dx_[i + ioffp][j] = 0.f;
2926 }
2927 for (j = 0; j < BXSIZE - k; ++j) {
2928 temp2dx_[i + ioffp][j + k] = xtemp_[i][j];
2929 }
2930 }
2931 }
2932
2933 nxbins = nxbins_;
2934
2935 return;
2936
2937 }
2938
2939
2940
2941
2942
2943
2944 void SiPixelTemplate::xtemp3d(int i, int j, std::vector<float>& xtemplate)
2945
2946 {
2947
2948 if (i >= 0 && i < nxbins_ && j <= i) {
2949 for (int k = 0; k < BXSIZE; ++k) {
2950 xtemplate[k] = temp2dx_[i][k] + temp2dx_[j][k];
2951 }
2952 } else {
2953 for (int k = 0; k < BXSIZE; ++k) {
2954 xtemplate[k] = 0.;
2955 }
2956 }
2957
2958 return;
2959
2960 }
2961
2962
2963
2964
2965
2966
2967
2968
2969
2970
2971
2972
2973
2974
2975
2976
2977
2978
2979
2980
2981
2982
2983
2984
2985
2986
2987
2988
2989
2990
2991 int SiPixelTemplate::qbin(int id,
2992 float cotalpha,
2993 float cotbeta,
2994 float locBz,
2995 float locBx,
2996 float qclus,
2997 float& pixmx,
2998 float& sigmay,
2999 float& deltay,
3000 float& sigmax,
3001 float& deltax,
3002 float& sy1,
3003 float& dy1,
3004 float& sy2,
3005 float& dy2,
3006 float& sx1,
3007 float& dx1,
3008 float& sx2,
3009 float& dx2)
3010 {
3011
3012
3013
3014
3015
3016
3017 int index = -1;
3018 for (int i = 0; i < (int)thePixelTemp_.size(); ++i) {
3019 if (id == thePixelTemp_[i].head.ID) {
3020 index = i;
3021 break;
3022 }
3023 }
3024
3025 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3026 if (index < 0 || index >= (int)thePixelTemp_.size()) {
3027 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::qbin can't find needed template ID = " << id << std::endl;
3028 }
3029 #else
3030 assert(index >= 0 && index < (int)thePixelTemp_.size());
3031 #endif
3032
3033
3034
3035 auto const& templ = thePixelTemp_[index];
3036
3037
3038
3039 auto acotb = std::abs(cotbeta);
3040
3041
3042
3043 auto cotalpha0 = thePixelTemp_[index].enty[0].cotalpha;
3044 auto qcorrect =
3045 std::sqrt((1.f + cotbeta * cotbeta + cotalpha * cotalpha) / (1.f + cotbeta * cotbeta + cotalpha0 * cotalpha0));
3046
3047
3048
3049 float cota = cotalpha;
3050 bool flip_x = false;
3051
3052 switch (dtype_) {
3053 case 0:
3054 case 1:
3055 case 2:
3056 case 3:
3057 case 4:
3058 case 5:
3059 if (locBx * locBz < 0.f) {
3060 cota = -cotalpha;
3061 flip_x = true;
3062 }
3063 break;
3064 default:
3065 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3066 throw cms::Exception("DataCorrupt")
3067 << "SiPixelTemplate::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
3068 #else
3069 std::cout << "SiPixelTemplate::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
3070 #endif
3071 }
3072
3073
3074
3075 auto qscale = thePixelTemp_[index].head.qscale;
3076
3077
3078
3079
3080
3081
3082
3083 auto Ny = thePixelTemp_[index].head.NTy;
3084 auto Nyx = thePixelTemp_[index].head.NTyx;
3085 auto Nxx = thePixelTemp_[index].head.NTxx;
3086
3087 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3088 if (Ny < 2 || Nyx < 1 || Nxx < 2) {
3089 throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny
3090 << "/" << Nyx << "/" << Nxx << std::endl;
3091 }
3092 #else
3093 assert(Ny > 1 && Nyx > 0 && Nxx > 1);
3094 #endif
3095
3096
3097
3098 dy1 = (1.f - yratio_) * enty0_->dyone + yratio_ * enty1_->dyone;
3099 if (flip_y_) {
3100 dy1 = -dy1;
3101 }
3102 sy1 = (1.f - yratio_) * enty0_->syone + yratio_ * enty1_->syone;
3103 dy2 = (1.f - yratio_) * enty0_->dytwo + yratio_ * enty1_->dytwo;
3104 if (flip_y_) {
3105 dy2 = -dy2;
3106 }
3107 sy2 = (1.f - yratio_) * enty0_->sytwo + yratio_ * enty1_->sytwo;
3108
3109 auto qavg = (1.f - yratio_) * enty0_->qavg + yratio_ * enty1_->qavg;
3110 qavg *= qcorrect;
3111 auto qmin = (1.f - yratio_) * enty0_->qmin + yratio_ * enty1_->qmin;
3112 qmin *= qcorrect;
3113 auto qmin2 = (1.f - yratio_) * enty0_->qmin2 + yratio_ * enty1_->qmin2;
3114 qmin2 *= qcorrect;
3115
3116 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3117 if (qavg <= 0.f || qmin <= 0.f) {
3118 throw cms::Exception("DataCorrupt")
3119 << "SiPixelTemplate::qbin, qavg or qmin <= 0,"
3120 << " Probably someone called the generic pixel reconstruction with an illegal trajectory state" << std::endl;
3121 }
3122 #else
3123 assert(qavg > 0.f && qmin > 0.f);
3124 #endif
3125
3126
3127
3128 auto qtotal = qscale * qclus;
3129
3130
3131 auto fq = qtotal / qavg;
3132 int binq;
3133
3134 if (fq > fbin_[0]) {
3135 binq = 0;
3136
3137
3138 } else {
3139 if (fq > fbin_[1]) {
3140 binq = 1;
3141 } else {
3142 if (fq > fbin_[2]) {
3143 binq = 2;
3144 } else {
3145 binq = 3;
3146 }
3147 }
3148 }
3149
3150 auto yavggen = (1.f - yratio_) * enty0_->yavggen[binq] + yratio_ * enty1_->yavggen[binq];
3151 if (flip_y_) {
3152 yavggen = -yavggen;
3153 }
3154 auto yrmsgen = (1.f - yratio_) * enty0_->yrmsgen[binq] + yratio_ * enty1_->yrmsgen[binq];
3155
3156
3157
3158 auto iylow = 0;
3159 auto iyhigh = 0;
3160 auto yxratio = 0.f;
3161
3162 {
3163 auto j = std::lower_bound(templ.cotbetaX.begin(), templ.cotbetaX.begin() + Nyx, acotb);
3164 if (j == templ.cotbetaX.begin() + Nyx) {
3165 --j;
3166 yxratio = 1.f;
3167 } else if (j == templ.cotbetaX.begin()) {
3168 ++j;
3169 yxratio = 0.f;
3170 } else {
3171 yxratio = (acotb - (*(j - 1))) / ((*j) - (*(j - 1)));
3172 }
3173
3174 iyhigh = j - templ.cotbetaX.begin();
3175 iylow = iyhigh - 1;
3176 }
3177
3178 auto ilow = 0;
3179 auto ihigh = 0;
3180 auto xxratio = 0.f;
3181
3182 {
3183 auto j = std::lower_bound(templ.cotalphaX.begin(), templ.cotalphaX.begin() + Nxx, cota);
3184 if (j == templ.cotalphaX.begin() + Nxx) {
3185 --j;
3186 xxratio = 1.f;
3187 } else if (j == templ.cotalphaX.begin()) {
3188 ++j;
3189 xxratio = 0.f;
3190 } else {
3191 xxratio = (cota - (*(j - 1))) / ((*j) - (*(j - 1)));
3192 }
3193
3194 ihigh = j - templ.cotalphaX.begin();
3195 ilow = ihigh - 1;
3196 }
3197
3198 dx1 =
3199 (1.f - xxratio) * thePixelTemp_[index].entx[0][ilow].dxone + xxratio * thePixelTemp_[index].entx[0][ihigh].dxone;
3200 if (flip_x) {
3201 dx1 = -dx1;
3202 }
3203 sx1 =
3204 (1.f - xxratio) * thePixelTemp_[index].entx[0][ilow].sxone + xxratio * thePixelTemp_[index].entx[0][ihigh].sxone;
3205 dx2 =
3206 (1.f - xxratio) * thePixelTemp_[index].entx[0][ilow].dxtwo + xxratio * thePixelTemp_[index].entx[0][ihigh].dxtwo;
3207 if (flip_x) {
3208 dx2 = -dx2;
3209 }
3210 sx2 =
3211 (1.f - xxratio) * thePixelTemp_[index].entx[0][ilow].sxtwo + xxratio * thePixelTemp_[index].entx[0][ihigh].sxtwo;
3212
3213
3214
3215 pixmx = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].pixmax +
3216 xxratio * thePixelTemp_[index].entx[iylow][ihigh].pixmax) +
3217 yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].pixmax +
3218 xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].pixmax);
3219
3220 auto xavggen = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].xavggen[binq] +
3221 xxratio * thePixelTemp_[index].entx[iylow][ihigh].xavggen[binq]) +
3222 yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].xavggen[binq] +
3223 xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].xavggen[binq]);
3224 if (flip_x) {
3225 xavggen = -xavggen;
3226 }
3227
3228 auto xrmsgen = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].xrmsgen[binq] +
3229 xxratio * thePixelTemp_[index].entx[iylow][ihigh].xrmsgen[binq]) +
3230 yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].xrmsgen[binq] +
3231 xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].xrmsgen[binq]);
3232
3233
3234
3235 sigmay = yrmsgen;
3236 deltay = yavggen;
3237
3238 sigmax = xrmsgen;
3239 deltax = xavggen;
3240
3241
3242
3243 if (qtotal < 0.95f * qmin) {
3244 binq = 5;
3245 } else {
3246 if (qtotal < 0.95f * qmin2) {
3247 binq = 4;
3248 }
3249 }
3250
3251 return binq;
3252
3253 }
3254
3255
3256
3257
3258
3259
3260
3261
3262
3263
3264
3265
3266
3267
3268
3269
3270
3271
3272
3273
3274
3275
3276
3277
3278
3279
3280
3281
3282
3283
3284 int SiPixelTemplate::qbin(int id,
3285 float cotalpha,
3286 float cotbeta,
3287 float locBz,
3288 float qclus,
3289 float& pixmx,
3290 float& sigmay,
3291 float& deltay,
3292 float& sigmax,
3293 float& deltax,
3294 float& sy1,
3295 float& dy1,
3296 float& sy2,
3297 float& dy2,
3298 float& sx1,
3299 float& dx1,
3300 float& sx2,
3301 float& dx2)
3302 {
3303
3304
3305
3306 float locBx = 1.f;
3307
3308 return SiPixelTemplate::qbin(id,
3309 cotalpha,
3310 cotbeta,
3311 locBz,
3312 locBx,
3313 qclus,
3314 pixmx,
3315 sigmay,
3316 deltay,
3317 sigmax,
3318 deltax,
3319 sy1,
3320 dy1,
3321 sy2,
3322 dy2,
3323 sx1,
3324 dx1,
3325 sx2,
3326 dx2);
3327
3328 }
3329
3330
3331
3332
3333
3334
3335
3336
3337
3338
3339
3340
3341
3342
3343
3344
3345
3346
3347
3348
3349
3350
3351
3352
3353
3354
3355
3356
3357
3358
3359
3360
3361
3362
3363
3364
3365
3366
3367
3368
3369
3370
3371
3372
3373 int SiPixelTemplate::qbin(int id, float cotalpha, float cotbeta, float qclus) {
3374
3375
3376
3377 float pixmx, sigmay, deltay, sigmax, deltax, sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, locBz;
3378
3379 float locBx = 1.f;
3380 if (cotbeta < 0.f) {
3381 locBx = -1.f;
3382 }
3383 locBz = locBx;
3384 if (cotalpha < 0.f) {
3385 locBz = -locBx;
3386 }
3387
3388 return SiPixelTemplate::qbin(id,
3389 cotalpha,
3390 cotbeta,
3391 locBz,
3392 locBx,
3393 qclus,
3394 pixmx,
3395 sigmay,
3396 deltay,
3397 sigmax,
3398 deltax,
3399 sy1,
3400 dy1,
3401 sy2,
3402 dy2,
3403 sx1,
3404 dx1,
3405 sx2,
3406 dx2);
3407
3408 }
3409
3410
3411
3412
3413
3414
3415
3416
3417 int SiPixelTemplate::qbin(int id, float cotbeta, float qclus) {
3418
3419
3420
3421 float pixmx, sigmay, deltay, sigmax, deltax, sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, locBz,
3422 locBx;
3423 const float cotalpha = 0.f;
3424 locBx = 1.f;
3425 if (cotbeta < 0.f) {
3426 locBx = -1.f;
3427 }
3428 locBz = locBx;
3429 if (cotalpha < 0.f) {
3430 locBz = -locBx;
3431 }
3432 return SiPixelTemplate::qbin(id,
3433 cotalpha,
3434 cotbeta,
3435 locBz,
3436 locBx,
3437 qclus,
3438 pixmx,
3439 sigmay,
3440 deltay,
3441 sigmax,
3442 deltax,
3443 sy1,
3444 dy1,
3445 sy2,
3446 dy2,
3447 sx1,
3448 dx1,
3449 sx2,
3450 dx2);
3451
3452 }
3453
3454
3455
3456
3457
3458
3459
3460
3461
3462
3463
3464
3465
3466
3467 void SiPixelTemplate::temperrors(int id,
3468 float cotalpha,
3469 float cotbeta,
3470 int qBin,
3471 float& sigmay,
3472 float& sigmax,
3473 float& sy1,
3474 float& sy2,
3475 float& sx1,
3476 float& sx2)
3477
3478 {
3479
3480
3481
3482 int i;
3483 int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, index;
3484 float yxratio, xxratio;
3485 float acotb;
3486 float yrms, xrms;
3487
3488
3489
3490
3491 index = -1;
3492 for (i = 0; i < (int)thePixelTemp_.size(); ++i) {
3493 if (id == thePixelTemp_[i].head.ID) {
3494 index = i;
3495 break;
3496 }
3497 }
3498
3499 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3500 if (index < 0 || index >= (int)thePixelTemp_.size()) {
3501 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::temperrors can't find needed template ID = " << id
3502 << std::endl;
3503 }
3504 #else
3505 assert(index >= 0 && index < (int)thePixelTemp_.size());
3506 #endif
3507
3508 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3509 if (qBin < 0 || qBin > 5) {
3510 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::temperrors called with illegal qBin = " << qBin
3511 << std::endl;
3512 }
3513 #else
3514 assert(qBin >= 0 && qBin < 6);
3515 #endif
3516
3517
3518
3519 if (qBin > 3) {
3520 qBin = 3;
3521 }
3522
3523
3524
3525
3526 acotb = std::abs(cotbeta);
3527
3528
3529
3530 Ny = thePixelTemp_[index].head.NTy;
3531 Nyx = thePixelTemp_[index].head.NTyx;
3532 Nxx = thePixelTemp_[index].head.NTxx;
3533
3534 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3535 if (Ny < 2 || Nyx < 1 || Nxx < 2) {
3536 throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny
3537 << "/" << Nyx << "/" << Nxx << std::endl;
3538 }
3539 #else
3540 assert(Ny > 1 && Nyx > 0 && Nxx > 1);
3541 #endif
3542
3543
3544
3545
3546
3547 sy1 = (1.f - yratio_) * enty0_->syone + yratio_ * enty1_->syone;
3548 sy2 = (1.f - yratio_) * enty0_->sytwo + yratio_ * enty1_->sytwo;
3549 yrms = (1.f - yratio_) * enty0_->yrms[qBin] + yratio_ * enty1_->yrms[qBin];
3550
3551
3552
3553 iylow = 0;
3554 yxratio = 0.f;
3555
3556 if (acotb >= thePixelTemp_[index].entx[Nyx - 1][0].cotbeta) {
3557 iylow = Nyx - 2;
3558 yxratio = 1.f;
3559
3560 } else if (acotb >= thePixelTemp_[index].entx[0][0].cotbeta) {
3561 for (i = 0; i < Nyx - 1; ++i) {
3562 if (thePixelTemp_[index].entx[i][0].cotbeta <= acotb && acotb < thePixelTemp_[index].entx[i + 1][0].cotbeta) {
3563 iylow = i;
3564 yxratio = (acotb - thePixelTemp_[index].entx[i][0].cotbeta) /
3565 (thePixelTemp_[index].entx[i + 1][0].cotbeta - thePixelTemp_[index].entx[i][0].cotbeta);
3566 break;
3567 }
3568 }
3569 }
3570
3571 iyhigh = iylow + 1;
3572
3573 ilow = 0;
3574 xxratio = 0.f;
3575
3576 if (cotalpha >= thePixelTemp_[index].entx[0][Nxx - 1].cotalpha) {
3577 ilow = Nxx - 2;
3578 xxratio = 1.f;
3579
3580 } else {
3581 if (cotalpha >= thePixelTemp_[index].entx[0][0].cotalpha) {
3582 for (i = 0; i < Nxx - 1; ++i) {
3583 if (thePixelTemp_[index].entx[0][i].cotalpha <= cotalpha &&
3584 cotalpha < thePixelTemp_[index].entx[0][i + 1].cotalpha) {
3585 ilow = i;
3586 xxratio = (cotalpha - thePixelTemp_[index].entx[0][i].cotalpha) /
3587 (thePixelTemp_[index].entx[0][i + 1].cotalpha - thePixelTemp_[index].entx[0][i].cotalpha);
3588 break;
3589 }
3590 }
3591 }
3592 }
3593
3594 ihigh = ilow + 1;
3595
3596 sx1 =
3597 (1.f - xxratio) * thePixelTemp_[index].entx[0][ilow].sxone + xxratio * thePixelTemp_[index].entx[0][ihigh].sxone;
3598 sx2 =
3599 (1.f - xxratio) * thePixelTemp_[index].entx[0][ilow].sxtwo + xxratio * thePixelTemp_[index].entx[0][ihigh].sxtwo;
3600
3601 xrms = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].xrms[qBin] +
3602 xxratio * thePixelTemp_[index].entx[iylow][ihigh].xrms[qBin]) +
3603 yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].xrms[qBin] +
3604 xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].xrms[qBin]);
3605
3606
3607
3608 sigmay = yrms;
3609
3610 sigmax = xrms;
3611
3612 return;
3613
3614 }
3615
3616
3617
3618
3619
3620
3621
3622
3623
3624
3625
3626
3627 void SiPixelTemplate::qbin_dist(int id,
3628 float cotalpha,
3629 float cotbeta,
3630 float qbin_frac[4],
3631 float& ny1_frac,
3632 float& ny2_frac,
3633 float& nx1_frac,
3634 float& nx2_frac)
3635
3636 {
3637
3638
3639
3640 int i;
3641 int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, index;
3642 float yxratio, xxratio;
3643 float acotb;
3644 float qfrac[4];
3645
3646
3647
3648
3649 index = -1;
3650 for (i = 0; i < (int)thePixelTemp_.size(); ++i) {
3651 if (id == thePixelTemp_[i].head.ID) {
3652 index = i;
3653
3654 break;
3655 }
3656 }
3657
3658 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3659 if (index < 0 || index >= (int)thePixelTemp_.size()) {
3660 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::temperrors can't find needed template ID = " << id
3661 << std::endl;
3662 }
3663 #else
3664 assert(index >= 0 && index < (int)thePixelTemp_.size());
3665 #endif
3666
3667
3668
3669
3670
3671 acotb = fabs((double)cotbeta);
3672
3673
3674
3675 Ny = thePixelTemp_[index].head.NTy;
3676 Nyx = thePixelTemp_[index].head.NTyx;
3677 Nxx = thePixelTemp_[index].head.NTxx;
3678
3679 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3680 if (Ny < 2 || Nyx < 1 || Nxx < 2) {
3681 throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny
3682 << "/" << Nyx << "/" << Nxx << std::endl;
3683 }
3684 #else
3685 assert(Ny > 1 && Nyx > 0 && Nxx > 1);
3686 #endif
3687
3688
3689 ny1_frac = (1.f - yratio_) * enty0_->fracyone + yratio_ * enty1_->fracyone;
3690 ny2_frac = (1.f - yratio_) * enty0_->fracytwo + yratio_ * enty1_->fracytwo;
3691
3692
3693
3694 iylow = 0;
3695 yxratio = 0.f;
3696
3697 if (acotb >= thePixelTemp_[index].entx[Nyx - 1][0].cotbeta) {
3698 iylow = Nyx - 2;
3699 yxratio = 1.f;
3700
3701 } else if (acotb >= thePixelTemp_[index].entx[0][0].cotbeta) {
3702 for (i = 0; i < Nyx - 1; ++i) {
3703 if (thePixelTemp_[index].entx[i][0].cotbeta <= acotb && acotb < thePixelTemp_[index].entx[i + 1][0].cotbeta) {
3704 iylow = i;
3705 yxratio = (acotb - thePixelTemp_[index].entx[i][0].cotbeta) /
3706 (thePixelTemp_[index].entx[i + 1][0].cotbeta - thePixelTemp_[index].entx[i][0].cotbeta);
3707 break;
3708 }
3709 }
3710 }
3711
3712 iyhigh = iylow + 1;
3713
3714 ilow = 0;
3715 xxratio = 0.f;
3716
3717 if (cotalpha >= thePixelTemp_[index].entx[0][Nxx - 1].cotalpha) {
3718 ilow = Nxx - 2;
3719 xxratio = 1.f;
3720
3721 } else {
3722 if (cotalpha >= thePixelTemp_[index].entx[0][0].cotalpha) {
3723 for (i = 0; i < Nxx - 1; ++i) {
3724 if (thePixelTemp_[index].entx[0][i].cotalpha <= cotalpha &&
3725 cotalpha < thePixelTemp_[index].entx[0][i + 1].cotalpha) {
3726 ilow = i;
3727 xxratio = (cotalpha - thePixelTemp_[index].entx[0][i].cotalpha) /
3728 (thePixelTemp_[index].entx[0][i + 1].cotalpha - thePixelTemp_[index].entx[0][i].cotalpha);
3729 break;
3730 }
3731 }
3732 }
3733 }
3734
3735 ihigh = ilow + 1;
3736
3737 for (i = 0; i < 3; ++i) {
3738 qfrac[i] = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].qbfrac[i] +
3739 xxratio * thePixelTemp_[index].entx[iylow][ihigh].qbfrac[i]) +
3740 yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].qbfrac[i] +
3741 xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].qbfrac[i]);
3742 }
3743 nx1_frac = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].fracxone +
3744 xxratio * thePixelTemp_[index].entx[iylow][ihigh].fracxone) +
3745 yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].fracxone +
3746 xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].fracxone);
3747 nx2_frac = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].fracxtwo +
3748 xxratio * thePixelTemp_[index].entx[iylow][ihigh].fracxtwo) +
3749 yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].fracxtwo +
3750 xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].fracxtwo);
3751
3752 qbin_frac[0] = qfrac[0];
3753 qbin_frac[1] = qbin_frac[0] + qfrac[1];
3754 qbin_frac[2] = qbin_frac[1] + qfrac[2];
3755 qbin_frac[3] = 1.f;
3756 return;
3757
3758 }
3759
3760
3761
3762
3763
3764
3765
3766
3767
3768
3769
3770 bool SiPixelTemplate::simpletemplate2D(
3771 float xhit, float yhit, std::vector<bool>& ydouble, std::vector<bool>& xdouble, float template2d[BXM2][BYM2]) {
3772
3773
3774 float x0, y0, xf, yf, xi, yi, sf, si, s0, qpix, slopey, slopex, ds;
3775 int i, j, jpix0, ipix0, jpixf, ipixf, jpix, ipix, nx, ny, anx, any, jmax, imax;
3776 float qtotal;
3777
3778 std::list<SimplePixel> list;
3779 std::list<SimplePixel>::iterator listIter, listEnd;
3780
3781
3782
3783 x0 = xhit - 0.5 * zsize_ * cota_current_;
3784 y0 = yhit - 0.5 * zsize_ * cotb_current_;
3785
3786 jpix0 = floor(x0 / xsize_) + 1;
3787 ipix0 = floor(y0 / ysize_) + 1;
3788
3789 if (jpix0 < 0 || jpix0 > BXM3) {
3790 return false;
3791 }
3792 if (ipix0 < 0 || ipix0 > BYM3) {
3793 return false;
3794 }
3795
3796 xf = xhit + 0.5 * zsize_ * cota_current_ + lorxwidth_;
3797 yf = yhit + 0.5 * zsize_ * cotb_current_ + lorywidth_;
3798
3799 jpixf = floor(xf / xsize_) + 1;
3800 ipixf = floor(yf / ysize_) + 1;
3801
3802 if (jpixf < 0 || jpixf > BXM3) {
3803 return false;
3804 }
3805 if (ipixf < 0 || ipixf > BYM3) {
3806 return false;
3807 }
3808
3809
3810
3811 sf = std::sqrt((xf - x0) * (xf - x0) + (yf - y0) * (yf - y0));
3812 if ((xf - x0) != 0.f) {
3813 slopey = (yf - y0) / (xf - x0);
3814 } else {
3815 slopey = 1.e10;
3816 }
3817 if ((yf - y0) != 0.f) {
3818 slopex = (xf - x0) / (yf - y0);
3819 } else {
3820 slopex = 1.e10;
3821 }
3822
3823
3824
3825 qtotal = qavg_avg_;
3826
3827 SimplePixel element;
3828 element.s = sf;
3829 element.x = xf;
3830 element.y = yf;
3831 element.i = ipixf;
3832 element.j = jpixf;
3833 element.btype = 0;
3834 list.push_back(element);
3835
3836
3837
3838 nx = jpixf - jpix0;
3839 anx = abs(nx);
3840 if (anx > 0) {
3841 if (nx > 0) {
3842 for (j = jpix0; j < jpixf; ++j) {
3843 xi = xsize_ * j;
3844 yi = slopey * (xi - x0) + y0;
3845 ipix = (int)(yi / ysize_) + 1;
3846 si = std::sqrt((xi - x0) * (xi - x0) + (yi - y0) * (yi - y0));
3847 element.s = si;
3848 element.x = xi;
3849 element.y = yi;
3850 element.i = ipix;
3851 element.j = j;
3852 element.btype = 1;
3853 list.push_back(element);
3854 }
3855 } else {
3856 for (j = jpix0; j > jpixf; --j) {
3857 xi = xsize_ * (j - 1);
3858 yi = slopey * (xi - x0) + y0;
3859 ipix = (int)(yi / ysize_) + 1;
3860 si = std::sqrt((xi - x0) * (xi - x0) + (yi - y0) * (yi - y0));
3861 element.s = si;
3862 element.x = xi;
3863 element.y = yi;
3864 element.i = ipix;
3865 element.j = j;
3866 element.btype = 1;
3867 list.push_back(element);
3868 }
3869 }
3870 }
3871
3872 ny = ipixf - ipix0;
3873 any = abs(ny);
3874 if (any > 0) {
3875 if (ny > 0) {
3876 for (i = ipix0; i < ipixf; ++i) {
3877 yi = ysize_ * i;
3878 xi = slopex * (yi - y0) + x0;
3879 jpix = (int)(xi / xsize_) + 1;
3880 si = std::sqrt((xi - x0) * (xi - x0) + (yi - y0) * (yi - y0));
3881 element.s = si;
3882 element.x = xi;
3883 element.y = yi;
3884 element.i = i;
3885 element.j = jpix;
3886 element.btype = 2;
3887 list.push_back(element);
3888 }
3889 } else {
3890 for (i = ipix0; i > ipixf; --i) {
3891 yi = ysize_ * (i - 1);
3892 xi = slopex * (yi - y0) + x0;
3893 jpix = (int)(xi / xsize_) + 1;
3894 si = std::sqrt((xi - x0) * (xi - x0) + (yi - y0) * (yi - y0));
3895 element.s = si;
3896 element.x = xi;
3897 element.y = yi;
3898 element.i = i;
3899 element.j = jpix;
3900 element.btype = 2;
3901 list.push_back(element);
3902 }
3903 }
3904 }
3905
3906 imax = std::max(ipix0, ipixf);
3907 jmax = std::max(jpix0, jpixf);
3908
3909
3910
3911 list.sort();
3912
3913
3914
3915 for (i = 1; i < imax; ++i) {
3916 if (ydouble[i - 1]) {
3917 listIter = list.begin();
3918 if (ny > 0) {
3919 while (listIter != list.end()) {
3920 if (listIter->i == i && listIter->btype == 2) {
3921 listIter = list.erase(listIter);
3922 continue;
3923 }
3924 if (listIter->i > i) {
3925 --(listIter->i);
3926 }
3927 ++listIter;
3928 }
3929 } else {
3930 while (listIter != list.end()) {
3931 if (listIter->i == i + 1 && listIter->btype == 2) {
3932 listIter = list.erase(listIter);
3933 continue;
3934 }
3935 if (listIter->i > i + 1) {
3936 --(listIter->i);
3937 }
3938 ++listIter;
3939 }
3940 }
3941 }
3942 }
3943
3944 for (j = 1; j < jmax; ++j) {
3945 if (xdouble[j - 1]) {
3946 listIter = list.begin();
3947 if (nx > 0) {
3948 while (listIter != list.end()) {
3949 if (listIter->j == j && listIter->btype == 1) {
3950 listIter = list.erase(listIter);
3951 continue;
3952 }
3953 if (listIter->j > j) {
3954 --(listIter->j);
3955 }
3956 ++listIter;
3957 }
3958 } else {
3959 while (listIter != list.end()) {
3960 if (listIter->j == j + 1 && listIter->btype == 1) {
3961 listIter = list.erase(listIter);
3962 continue;
3963 }
3964 if (listIter->j > j + 1) {
3965 --(listIter->j);
3966 }
3967 ++listIter;
3968 }
3969 }
3970 }
3971 }
3972
3973
3974
3975 s0 = 0.f;
3976 listIter = list.begin();
3977 listEnd = list.end();
3978 for (; listIter != listEnd; ++listIter) {
3979 si = listIter->s;
3980 ds = si - s0;
3981 s0 = si;
3982 j = listIter->j;
3983 i = listIter->i;
3984 if (sf > 0.f) {
3985 qpix = qtotal * ds / sf;
3986 } else {
3987 qpix = qtotal;
3988 }
3989 template2d[j][i] += qpix;
3990 }
3991
3992 return true;
3993
3994 }
3995
3996
3997
3998
3999
4000
4001
4002 void SiPixelTemplate::vavilov_pars(double& mpv, double& sigma, double& kappa)
4003
4004 {
4005
4006 int i;
4007 int ilow, ihigh, Ny;
4008 float yratio, cotb, cotalpha0, arg;
4009
4010
4011
4012 cotalpha0 = thePixelTemp_[index_id_].enty[0].cotalpha;
4013 arg = cotb_current_ * cotb_current_ + cota_current_ * cota_current_ - cotalpha0 * cotalpha0;
4014 if (arg < 0.f)
4015 arg = 0.f;
4016 cotb = std::sqrt(arg);
4017
4018
4019
4020 Ny = thePixelTemp_[index_id_].head.NTy;
4021
4022 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
4023 if (Ny < 2) {
4024 throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny = " << Ny
4025 << std::endl;
4026 }
4027 #else
4028 assert(Ny > 1);
4029 #endif
4030
4031
4032
4033 ilow = 0;
4034 yratio = 0.f;
4035
4036 if (cotb >= thePixelTemp_[index_id_].enty[Ny - 1].cotbeta) {
4037 ilow = Ny - 2;
4038 yratio = 1.f;
4039
4040 } else {
4041 if (cotb >= thePixelTemp_[index_id_].enty[0].cotbeta) {
4042 for (i = 0; i < Ny - 1; ++i) {
4043 if (thePixelTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index_id_].enty[i + 1].cotbeta) {
4044 ilow = i;
4045 yratio = (cotb - thePixelTemp_[index_id_].enty[i].cotbeta) /
4046 (thePixelTemp_[index_id_].enty[i + 1].cotbeta - thePixelTemp_[index_id_].enty[i].cotbeta);
4047 break;
4048 }
4049 }
4050 }
4051 }
4052
4053 ihigh = ilow + 1;
4054
4055
4056
4057 mpvvav_ = (1.f - yratio) * thePixelTemp_[index_id_].enty[ilow].mpvvav +
4058 yratio * thePixelTemp_[index_id_].enty[ihigh].mpvvav;
4059 sigmavav_ = (1.f - yratio) * thePixelTemp_[index_id_].enty[ilow].sigmavav +
4060 yratio * thePixelTemp_[index_id_].enty[ihigh].sigmavav;
4061 kappavav_ = (1.f - yratio) * thePixelTemp_[index_id_].enty[ilow].kappavav +
4062 yratio * thePixelTemp_[index_id_].enty[ihigh].kappavav;
4063
4064
4065
4066
4067 if (kappavav_ <= 0.01f) {
4068 LOGERROR("SiPixelTemplate") << "Vavilov kappa value is " << kappavav_ << " changing it to be above 0.01" << ENDL;
4069 kappavav_ = 0.01f + 0.0000001f;
4070 }
4071
4072 mpv = (double)mpvvav_;
4073 sigma = (double)sigmavav_;
4074 kappa = (double)kappavav_;
4075
4076 return;
4077
4078 }
4079
4080
4081
4082
4083
4084
4085
4086 void SiPixelTemplate::vavilov2_pars(double& mpv, double& sigma, double& kappa)
4087
4088 {
4089
4090 int i;
4091 int ilow, ihigh, Ny;
4092 float yratio, cotb, cotalpha0, arg;
4093
4094
4095
4096 cotalpha0 = thePixelTemp_[index_id_].enty[0].cotalpha;
4097 arg = cotb_current_ * cotb_current_ + cota_current_ * cota_current_ - cotalpha0 * cotalpha0;
4098 if (arg < 0.f)
4099 arg = 0.f;
4100 cotb = std::sqrt(arg);
4101
4102
4103
4104 Ny = thePixelTemp_[index_id_].head.NTy;
4105
4106 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
4107 if (Ny < 2) {
4108 throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny = " << Ny
4109 << std::endl;
4110 }
4111 #else
4112 assert(Ny > 1);
4113 #endif
4114
4115
4116
4117 ilow = 0;
4118 yratio = 0.f;
4119
4120 if (cotb >= thePixelTemp_[index_id_].enty[Ny - 1].cotbeta) {
4121 ilow = Ny - 2;
4122 yratio = 1.f;
4123
4124 } else {
4125 if (cotb >= thePixelTemp_[index_id_].enty[0].cotbeta) {
4126 for (i = 0; i < Ny - 1; ++i) {
4127 if (thePixelTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index_id_].enty[i + 1].cotbeta) {
4128 ilow = i;
4129 yratio = (cotb - thePixelTemp_[index_id_].enty[i].cotbeta) /
4130 (thePixelTemp_[index_id_].enty[i + 1].cotbeta - thePixelTemp_[index_id_].enty[i].cotbeta);
4131 break;
4132 }
4133 }
4134 }
4135 }
4136
4137 ihigh = ilow + 1;
4138
4139
4140
4141 mpvvav2_ = (1.f - yratio) * thePixelTemp_[index_id_].enty[ilow].mpvvav2 +
4142 yratio * thePixelTemp_[index_id_].enty[ihigh].mpvvav2;
4143 sigmavav2_ = (1.f - yratio) * thePixelTemp_[index_id_].enty[ilow].sigmavav2 +
4144 yratio * thePixelTemp_[index_id_].enty[ihigh].sigmavav2;
4145 kappavav2_ = (1.f - yratio) * thePixelTemp_[index_id_].enty[ilow].kappavav2 +
4146 yratio * thePixelTemp_[index_id_].enty[ihigh].kappavav2;
4147
4148
4149
4150 mpv = (double)mpvvav2_;
4151 sigma = (double)sigmavav2_;
4152 kappa = (double)kappavav2_;
4153
4154 return;
4155
4156 }