File indexing completed on 2021-08-17 00:01:51
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 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0028 #else
0029 #include <math.h>
0030 #endif
0031 #include <algorithm>
0032 #include <vector>
0033
0034 #include <iostream>
0035 #include <iomanip>
0036 #include <sstream>
0037 #include <fstream>
0038
0039 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0040 #include "CondFormats/SiPixelTransient/interface/SiPixelTemplate2D.h"
0041 #include "FWCore/Utilities/interface/FileInPath.h"
0042 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0043 #include "FWCore/Utilities/interface/isFinite.h"
0044 #define LOGERROR(x) LogError(x)
0045 #define LOGINFO(x) LogInfo(x)
0046 #define ENDL " "
0047 #include "FWCore/Utilities/interface/Exception.h"
0048 using namespace edm;
0049 #else
0050 #include "SiPixelTemplate2D.h"
0051 #define LOGERROR(x) std::cout << x << ": "
0052 #define LOGINFO(x) std::cout << x << ": "
0053 #define ENDL std::endl
0054 #endif
0055
0056
0057
0058
0059
0060
0061
0062 bool SiPixelTemplate2D::pushfile(int filenum, std::vector<SiPixelTemplateStore2D>& pixelTemp, std::string dir) {
0063
0064
0065
0066 const int code_version = {21};
0067
0068
0069 std::string tempfile = std::to_string(filenum);
0070
0071
0072
0073 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0074
0075 int nzeros = 4 - tempfile.length();
0076 if (nzeros > 0)
0077 tempfile = std::string(nzeros, '0') + tempfile;
0078
0079
0080 tempfile = dir + "template_summary2D_zp" + tempfile + ".out";
0081 edm::FileInPath file(tempfile);
0082 tempfile = file.fullPath();
0083
0084 #else
0085
0086 std::ostringstream tout;
0087 tout << "template_summary2D_zp" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
0088 tempfile = tout.str();
0089
0090 #endif
0091
0092
0093
0094 std::ifstream in_file(tempfile);
0095 if (in_file.is_open() && in_file.good()) {
0096
0097 SiPixelTemplateStore2D theCurrentTemp;
0098
0099
0100 char c;
0101 for (int i = 0; (c = in_file.get()) != '\n'; ++i) {
0102 if (i < 79) {
0103 theCurrentTemp.head.title[i] = c;
0104 } else {
0105 theCurrentTemp.head.title[79] = '\0';
0106 }
0107 }
0108 LOGINFO("SiPixelTemplate2D") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
0109
0110
0111 in_file >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >>
0112 theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx >> theCurrentTemp.head.Dtype >>
0113 theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >>
0114 theCurrentTemp.head.qscale >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >>
0115 theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >>
0116 theCurrentTemp.head.zsize;
0117
0118 if (in_file.fail()) {
0119 LOGERROR("SiPixelTemplate2D") << "Error reading file 0A, no template load" << ENDL;
0120 return false;
0121 }
0122
0123 if (theCurrentTemp.head.templ_version > 17) {
0124 in_file >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >> theCurrentTemp.head.lorxbias >>
0125 theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >> theCurrentTemp.head.fbin[2];
0126
0127 if (in_file.fail()) {
0128 LOGERROR("SiPixelTemplate2D") << "Error reading file 0B, no template load" << ENDL;
0129 return false;
0130 }
0131 } else {
0132
0133 theCurrentTemp.head.ss50 = theCurrentTemp.head.s50;
0134 theCurrentTemp.head.lorybias = theCurrentTemp.head.lorywidth / 2.f;
0135 theCurrentTemp.head.lorxbias = theCurrentTemp.head.lorxwidth / 2.f;
0136 theCurrentTemp.head.fbin[0] = 1.5f;
0137 theCurrentTemp.head.fbin[1] = 1.00f;
0138 theCurrentTemp.head.fbin[2] = 0.85f;
0139 }
0140
0141 LOGINFO("SiPixelTemplate2D") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version "
0142 << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
0143 << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx
0144 << ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
0145 << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
0146 << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence
0147 << ", Q-scaling factor " << theCurrentTemp.head.qscale << ", 1/2 multi dcol threshold "
0148 << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold "
0149 << theCurrentTemp.head.ss50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth
0150 << ", y Lorentz Bias " << theCurrentTemp.head.lorybias << ", x Lorentz width "
0151 << theCurrentTemp.head.lorxwidth << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
0152 << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] << ", "
0153 << theCurrentTemp.head.fbin[1] << ", " << theCurrentTemp.head.fbin[2]
0154 << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size "
0155 << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
0156
0157 if (theCurrentTemp.head.templ_version < code_version) {
0158 LOGERROR("SiPixelTemplate2D") << "code expects version " << code_version << ", no template load" << ENDL;
0159 return false;
0160 }
0161
0162 if (theCurrentTemp.head.NTy != 0) {
0163 LOGERROR("SiPixelTemplate2D")
0164 << "Trying to load 1-d template info into the 2-d template object, check your DB/global tag!" << ENDL;
0165 return false;
0166 }
0167
0168
0169
0170 theCurrentTemp.entry.resize(theCurrentTemp.head.NTyx);
0171 for (auto& item : theCurrentTemp.entry)
0172 item.resize(theCurrentTemp.head.NTxx);
0173
0174
0175
0176 for (int iy = 0; iy < theCurrentTemp.head.NTyx; ++iy) {
0177 for (int jx = 0; jx < theCurrentTemp.head.NTxx; ++jx) {
0178 in_file >> theCurrentTemp.entry[iy][jx].runnum >> theCurrentTemp.entry[iy][jx].costrk[0] >>
0179 theCurrentTemp.entry[iy][jx].costrk[1] >> theCurrentTemp.entry[iy][jx].costrk[2];
0180
0181 if (in_file.fail()) {
0182 LOGERROR("SiPixelTemplate2D") << "Error reading file 1, no template load, run # "
0183 << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0184 return false;
0185 }
0186
0187
0188
0189 theCurrentTemp.entry[iy][jx].cotalpha =
0190 theCurrentTemp.entry[iy][jx].costrk[0] / theCurrentTemp.entry[iy][jx].costrk[2];
0191
0192 theCurrentTemp.entry[iy][jx].cotbeta =
0193 theCurrentTemp.entry[iy][jx].costrk[1] / theCurrentTemp.entry[iy][jx].costrk[2];
0194
0195 in_file >> theCurrentTemp.entry[iy][jx].qavg >> theCurrentTemp.entry[iy][jx].pixmax >>
0196 theCurrentTemp.entry[iy][jx].sxymax >> theCurrentTemp.entry[iy][jx].iymin >>
0197 theCurrentTemp.entry[iy][jx].iymax >> theCurrentTemp.entry[iy][jx].jxmin >>
0198 theCurrentTemp.entry[iy][jx].jxmax;
0199
0200 if (in_file.fail()) {
0201 LOGERROR("SiPixelTemplate2D") << "Error reading file 2, no template load, run # "
0202 << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0203 return false;
0204 }
0205
0206 for (int k = 0; k < 2; ++k) {
0207 in_file >> theCurrentTemp.entry[iy][jx].xypar[k][0] >> theCurrentTemp.entry[iy][jx].xypar[k][1] >>
0208 theCurrentTemp.entry[iy][jx].xypar[k][2] >> theCurrentTemp.entry[iy][jx].xypar[k][3] >>
0209 theCurrentTemp.entry[iy][jx].xypar[k][4];
0210
0211 if (in_file.fail()) {
0212 LOGERROR("SiPixelTemplate2D")
0213 << "Error reading file 3, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0214 return false;
0215 }
0216 }
0217
0218 for (int k = 0; k < 2; ++k) {
0219 in_file >> theCurrentTemp.entry[iy][jx].lanpar[k][0] >> theCurrentTemp.entry[iy][jx].lanpar[k][1] >>
0220 theCurrentTemp.entry[iy][jx].lanpar[k][2] >> theCurrentTemp.entry[iy][jx].lanpar[k][3] >>
0221 theCurrentTemp.entry[iy][jx].lanpar[k][4];
0222
0223 if (in_file.fail()) {
0224 LOGERROR("SiPixelTemplate2D")
0225 << "Error reading file 4, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0226 return false;
0227 }
0228 }
0229
0230
0231
0232 float dummy[T2YSIZE];
0233 for (int l = 0; l < 7; ++l) {
0234 for (int k = 0; k < 7; ++k) {
0235 for (int j = 0; j < T2XSIZE; ++j) {
0236 for (int i = 0; i < T2YSIZE; ++i) {
0237 in_file >> dummy[i];
0238 }
0239 if (in_file.fail()) {
0240 LOGERROR("SiPixelTemplate2D")
0241 << "Error reading file 5, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0242 return false;
0243 }
0244 for (int i = 0; i < T2YSIZE; ++i) {
0245 theCurrentTemp.entry[iy][jx].xytemp[k][l][i][j] = (short int)dummy[i];
0246 }
0247 }
0248 }
0249 }
0250
0251 in_file >> theCurrentTemp.entry[iy][jx].chi2ppix >> theCurrentTemp.entry[iy][jx].chi2scale >>
0252 theCurrentTemp.entry[iy][jx].offsetx[0] >> theCurrentTemp.entry[iy][jx].offsetx[1] >>
0253 theCurrentTemp.entry[iy][jx].offsetx[2] >> theCurrentTemp.entry[iy][jx].offsetx[3] >>
0254 theCurrentTemp.entry[iy][jx].offsety[0] >> theCurrentTemp.entry[iy][jx].offsety[1] >>
0255 theCurrentTemp.entry[iy][jx].offsety[2] >> theCurrentTemp.entry[iy][jx].offsety[3];
0256
0257 if (in_file.fail()) {
0258 LOGERROR("SiPixelTemplate2D") << "Error reading file 6, no template load, run # "
0259 << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0260 return false;
0261 }
0262
0263 in_file >> theCurrentTemp.entry[iy][jx].clsleny >> theCurrentTemp.entry[iy][jx].clslenx >>
0264 theCurrentTemp.entry[iy][jx].mpvvav >> theCurrentTemp.entry[iy][jx].sigmavav >>
0265 theCurrentTemp.entry[iy][jx].kappavav >> theCurrentTemp.entry[iy][jx].scalexavg >>
0266 theCurrentTemp.entry[iy][jx].scaleyavg >> theCurrentTemp.entry[iy][jx].delyavg >>
0267 theCurrentTemp.entry[iy][jx].delysig >> theCurrentTemp.entry[iy][jx].spare[0];
0268
0269 if (in_file.fail()) {
0270 LOGERROR("SiPixelTemplate2D") << "Error reading file 7, no template load, run # "
0271 << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0272 return false;
0273 }
0274
0275 in_file >> theCurrentTemp.entry[iy][jx].scalex[0] >> theCurrentTemp.entry[iy][jx].scalex[1] >>
0276 theCurrentTemp.entry[iy][jx].scalex[2] >> theCurrentTemp.entry[iy][jx].scalex[3] >>
0277 theCurrentTemp.entry[iy][jx].scaley[0] >> theCurrentTemp.entry[iy][jx].scaley[1] >>
0278 theCurrentTemp.entry[iy][jx].scaley[2] >> theCurrentTemp.entry[iy][jx].scaley[3] >>
0279 theCurrentTemp.entry[iy][jx].spare[1] >> theCurrentTemp.entry[iy][jx].spare[2];
0280
0281 if (in_file.fail()) {
0282 LOGERROR("SiPixelTemplate2D") << "Error reading file 8, no template load, run # "
0283 << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0284 return false;
0285 }
0286 }
0287 }
0288
0289 in_file.close();
0290
0291
0292
0293 pixelTemp.push_back(theCurrentTemp);
0294
0295 return true;
0296
0297 } else {
0298
0299
0300 LOGERROR("SiPixelTemplate2D") << "Error opening File" << tempfile << ENDL;
0301 return false;
0302 }
0303
0304 }
0305
0306 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0307
0308
0309
0310
0311
0312
0313 bool SiPixelTemplate2D::pushfile(const SiPixel2DTemplateDBObject& dbobject,
0314 std::vector<SiPixelTemplateStore2D>& pixelTemp) {
0315
0316
0317 const int code_version = {21};
0318
0319
0320 SiPixel2DTemplateDBObject db = dbobject;
0321
0322
0323 SiPixelTemplateStore2D theCurrentTemp;
0324
0325
0326 for (int m = 0; m < db.numOfTempl(); ++m) {
0327
0328
0329 SiPixel2DTemplateDBObject::char2float temp;
0330 for (int i = 0; i < 20; ++i) {
0331 temp.f = db.sVector()[db.index()];
0332 theCurrentTemp.head.title[4 * i] = temp.c[0];
0333 theCurrentTemp.head.title[4 * i + 1] = temp.c[1];
0334 theCurrentTemp.head.title[4 * i + 2] = temp.c[2];
0335 theCurrentTemp.head.title[4 * i + 3] = temp.c[3];
0336 db.incrementIndex(1);
0337 }
0338 theCurrentTemp.head.title[79] = '\0';
0339 LOGINFO("SiPixelTemplate2D") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
0340
0341
0342
0343 db >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >>
0344 theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx >> theCurrentTemp.head.Dtype >>
0345 theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >>
0346 theCurrentTemp.head.qscale >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >>
0347 theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >>
0348 theCurrentTemp.head.zsize;
0349
0350 if (db.fail()) {
0351 LOGERROR("SiPixelTemplate2D") << "Error reading file 0A, no template load" << ENDL;
0352 return false;
0353 }
0354
0355 LOGINFO("SiPixelTemplate2D") << "Loading Pixel Template File - " << theCurrentTemp.head.title
0356 << " code version = " << code_version << " object version "
0357 << theCurrentTemp.head.templ_version << ENDL;
0358
0359 if (theCurrentTemp.head.templ_version > 17) {
0360 db >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >> theCurrentTemp.head.lorxbias >>
0361 theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >> theCurrentTemp.head.fbin[2];
0362
0363 if (db.fail()) {
0364 LOGERROR("SiPixelTemplate2D") << "Error reading file 0B, no template load" << ENDL;
0365 return false;
0366 }
0367 } else {
0368
0369 theCurrentTemp.head.ss50 = theCurrentTemp.head.s50;
0370 theCurrentTemp.head.lorybias = theCurrentTemp.head.lorywidth / 2.f;
0371 theCurrentTemp.head.lorxbias = theCurrentTemp.head.lorxwidth / 2.f;
0372 theCurrentTemp.head.fbin[0] = 1.50f;
0373 theCurrentTemp.head.fbin[1] = 1.00f;
0374 theCurrentTemp.head.fbin[2] = 0.85f;
0375 }
0376
0377 LOGINFO("SiPixelTemplate2D") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version "
0378 << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
0379 << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx
0380 << ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
0381 << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
0382 << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence
0383 << ", Q-scaling factor " << theCurrentTemp.head.qscale << ", 1/2 multi dcol threshold "
0384 << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold "
0385 << theCurrentTemp.head.ss50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth
0386 << ", y Lorentz Bias " << theCurrentTemp.head.lorybias << ", x Lorentz width "
0387 << theCurrentTemp.head.lorxwidth << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
0388 << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] << ", "
0389 << theCurrentTemp.head.fbin[1] << ", " << theCurrentTemp.head.fbin[2]
0390 << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size "
0391 << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
0392
0393 if (theCurrentTemp.head.templ_version < code_version) {
0394 LOGINFO("SiPixelTemplate2D") << "code expects version " << code_version << " finds "
0395 << theCurrentTemp.head.templ_version << ", load anyway " << ENDL;
0396 }
0397
0398 if (theCurrentTemp.head.NTy != 0) {
0399 LOGERROR("SiPixelTemplate2D")
0400 << "Trying to load 1-d template info into the 2-d template object, check your DB/global tag!" << ENDL;
0401 return false;
0402 }
0403
0404
0405
0406 theCurrentTemp.entry.resize(theCurrentTemp.head.NTyx);
0407 for (auto& item : theCurrentTemp.entry)
0408 item.resize(theCurrentTemp.head.NTxx);
0409
0410
0411
0412 for (int iy = 0; iy < theCurrentTemp.head.NTyx; ++iy) {
0413 for (int jx = 0; jx < theCurrentTemp.head.NTxx; ++jx) {
0414 db >> theCurrentTemp.entry[iy][jx].runnum >> theCurrentTemp.entry[iy][jx].costrk[0] >>
0415 theCurrentTemp.entry[iy][jx].costrk[1] >> theCurrentTemp.entry[iy][jx].costrk[2];
0416
0417 if (db.fail()) {
0418 LOGERROR("SiPixelTemplate2D") << "Error reading file 1, no template load, run # "
0419 << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0420 return false;
0421 }
0422
0423
0424
0425 theCurrentTemp.entry[iy][jx].cotalpha =
0426 theCurrentTemp.entry[iy][jx].costrk[0] / theCurrentTemp.entry[iy][jx].costrk[2];
0427
0428 theCurrentTemp.entry[iy][jx].cotbeta =
0429 theCurrentTemp.entry[iy][jx].costrk[1] / theCurrentTemp.entry[iy][jx].costrk[2];
0430
0431 db >> theCurrentTemp.entry[iy][jx].qavg >> theCurrentTemp.entry[iy][jx].pixmax >>
0432 theCurrentTemp.entry[iy][jx].sxymax >> theCurrentTemp.entry[iy][jx].iymin >>
0433 theCurrentTemp.entry[iy][jx].iymax >> theCurrentTemp.entry[iy][jx].jxmin >>
0434 theCurrentTemp.entry[iy][jx].jxmax;
0435
0436 if (db.fail()) {
0437 LOGERROR("SiPixelTemplate2D") << "Error reading file 2, no template load, run # "
0438 << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0439 return false;
0440 }
0441
0442 for (int k = 0; k < 2; ++k) {
0443 db >> theCurrentTemp.entry[iy][jx].xypar[k][0] >> theCurrentTemp.entry[iy][jx].xypar[k][1] >>
0444 theCurrentTemp.entry[iy][jx].xypar[k][2] >> theCurrentTemp.entry[iy][jx].xypar[k][3] >>
0445 theCurrentTemp.entry[iy][jx].xypar[k][4];
0446
0447 if (db.fail()) {
0448 LOGERROR("SiPixelTemplate2D")
0449 << "Error reading file 3, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0450 return false;
0451 }
0452 }
0453
0454 for (int k = 0; k < 2; ++k) {
0455 db >> theCurrentTemp.entry[iy][jx].lanpar[k][0] >> theCurrentTemp.entry[iy][jx].lanpar[k][1] >>
0456 theCurrentTemp.entry[iy][jx].lanpar[k][2] >> theCurrentTemp.entry[iy][jx].lanpar[k][3] >>
0457 theCurrentTemp.entry[iy][jx].lanpar[k][4];
0458
0459 if (db.fail()) {
0460 LOGERROR("SiPixelTemplate2D")
0461 << "Error reading file 4, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0462 return false;
0463 }
0464 }
0465
0466
0467
0468 float dummy[T2YSIZE];
0469 for (int l = 0; l < 7; ++l) {
0470 for (int k = 0; k < 7; ++k) {
0471 for (int j = 0; j < T2XSIZE; ++j) {
0472 for (int i = 0; i < T2YSIZE; ++i) {
0473 db >> dummy[i];
0474 }
0475 if (db.fail()) {
0476 LOGERROR("SiPixelTemplate2D")
0477 << "Error reading file 5, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0478 return false;
0479 }
0480 for (int i = 0; i < T2YSIZE; ++i) {
0481 theCurrentTemp.entry[iy][jx].xytemp[k][l][i][j] = (short int)dummy[i];
0482 }
0483 }
0484 }
0485 }
0486
0487 db >> theCurrentTemp.entry[iy][jx].chi2ppix >> theCurrentTemp.entry[iy][jx].chi2scale >>
0488 theCurrentTemp.entry[iy][jx].offsetx[0] >> theCurrentTemp.entry[iy][jx].offsetx[1] >>
0489 theCurrentTemp.entry[iy][jx].offsetx[2] >> theCurrentTemp.entry[iy][jx].offsetx[3] >>
0490 theCurrentTemp.entry[iy][jx].offsety[0] >> theCurrentTemp.entry[iy][jx].offsety[1] >>
0491 theCurrentTemp.entry[iy][jx].offsety[2] >> theCurrentTemp.entry[iy][jx].offsety[3];
0492
0493 if (db.fail()) {
0494 LOGERROR("SiPixelTemplate2D") << "Error reading file 6, no template load, run # "
0495 << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0496 return false;
0497 }
0498
0499 db >> theCurrentTemp.entry[iy][jx].clsleny >> theCurrentTemp.entry[iy][jx].clslenx >>
0500 theCurrentTemp.entry[iy][jx].mpvvav >> theCurrentTemp.entry[iy][jx].sigmavav >>
0501 theCurrentTemp.entry[iy][jx].kappavav >> theCurrentTemp.entry[iy][jx].scalexavg >>
0502 theCurrentTemp.entry[iy][jx].scaleyavg >> theCurrentTemp.entry[iy][jx].delyavg >>
0503 theCurrentTemp.entry[iy][jx].delysig >> theCurrentTemp.entry[iy][jx].spare[0];
0504
0505 if (db.fail()) {
0506 LOGERROR("SiPixelTemplate2D") << "Error reading file 7, no template load, run # "
0507 << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0508 return false;
0509 }
0510
0511 db >> theCurrentTemp.entry[iy][jx].scalex[0] >> theCurrentTemp.entry[iy][jx].scalex[1] >>
0512 theCurrentTemp.entry[iy][jx].scalex[2] >> theCurrentTemp.entry[iy][jx].scalex[3] >>
0513 theCurrentTemp.entry[iy][jx].scaley[0] >> theCurrentTemp.entry[iy][jx].scaley[1] >>
0514 theCurrentTemp.entry[iy][jx].scaley[2] >> theCurrentTemp.entry[iy][jx].scaley[3] >>
0515 theCurrentTemp.entry[iy][jx].spare[1] >> theCurrentTemp.entry[iy][jx].spare[2];
0516
0517 if (db.fail()) {
0518 LOGERROR("SiPixelTemplate2D") << "Error reading file 8, no template load, run # "
0519 << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0520 return false;
0521 }
0522 }
0523 }
0524
0525
0526
0527 pixelTemp.push_back(theCurrentTemp);
0528 }
0529
0530 return true;
0531
0532 }
0533
0534 #endif
0535
0536 bool SiPixelTemplate2D::getid(int id) {
0537 if (id != id_current_) {
0538
0539
0540 index_id_ = -1;
0541 for (int i = 0; i < (int)thePixelTemp_.size(); ++i) {
0542 if (id == thePixelTemp_[i].head.ID) {
0543 index_id_ = i;
0544 id_current_ = id;
0545
0546
0547
0548 Dtype_ = thePixelTemp_[index_id_].head.Dtype;
0549
0550
0551
0552 qscale_ = thePixelTemp_[index_id_].head.qscale;
0553
0554
0555
0556 s50_ = thePixelTemp_[index_id_].head.s50;
0557
0558
0559
0560 for (int j = 0; j < 3; ++j) {
0561 fbin_[j] = thePixelTemp_[index_id_].head.fbin[j];
0562 }
0563
0564
0565
0566 lorywidth_ = thePixelTemp_[index_id_].head.lorywidth;
0567 lorxwidth_ = thePixelTemp_[index_id_].head.lorxwidth;
0568
0569
0570
0571 xsize_ = thePixelTemp_[index_id_].head.xsize;
0572 ysize_ = thePixelTemp_[index_id_].head.ysize;
0573 zsize_ = thePixelTemp_[index_id_].head.zsize;
0574
0575
0576
0577 Nyx_ = thePixelTemp_[index_id_].head.NTyx;
0578 Nxx_ = thePixelTemp_[index_id_].head.NTxx;
0579 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0580 if (Nyx_ < 2 || Nxx_ < 2) {
0581 throw cms::Exception("DataCorrupt") << "template ID = " << id_current_
0582 << "has too few entries: Nyx/Nxx = " << Nyx_ << "/" << Nxx_ << std::endl;
0583 }
0584 #else
0585 assert(Nyx_ > 1 && Nxx_ > 1);
0586 #endif
0587 int imidx = Nxx_ / 2;
0588
0589 cotalpha0_ = thePixelTemp_[index_id_].entry[0][0].cotalpha;
0590 cotalpha1_ = thePixelTemp_[index_id_].entry[0][Nxx_ - 1].cotalpha;
0591 deltacota_ = (cotalpha1_ - cotalpha0_) / (float)(Nxx_ - 1);
0592
0593 cotbeta0_ = thePixelTemp_[index_id_].entry[0][imidx].cotbeta;
0594 cotbeta1_ = thePixelTemp_[index_id_].entry[Nyx_ - 1][imidx].cotbeta;
0595 deltacotb_ = (cotbeta1_ - cotbeta0_) / (float)(Nyx_ - 1);
0596
0597 break;
0598 }
0599 }
0600 }
0601
0602 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0603 if (index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
0604 throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::interpolate can't find needed template ID = " << id
0605 << ", Are you using the correct global tag?" << std::endl;
0606 }
0607 #else
0608 assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
0609 #endif
0610 return true;
0611 }
0612
0613
0614
0615
0616
0617
0618
0619
0620
0621
0622
0623
0624
0625
0626 bool SiPixelTemplate2D::interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx) {
0627
0628
0629
0630
0631 float acotb, dcota, dcotb;
0632
0633
0634
0635 if (id != id_current_ || cotalpha != cota_current_ || cotbeta != cotb_current_) {
0636 cota_current_ = cotalpha;
0637 cotb_current_ = cotbeta;
0638
0639 success_ = getid(id);
0640 }
0641
0642 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0643 if (index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
0644 throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::interpolate can't find needed template ID = " << id
0645 << ", Are you using the correct global tag?" << std::endl;
0646 }
0647 #else
0648 assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
0649 #endif
0650
0651
0652
0653 float cota = cotalpha;
0654 flip_x_ = false;
0655 flip_y_ = false;
0656 switch (Dtype_) {
0657 case 0:
0658 if (cotbeta < 0.f) {
0659 flip_y_ = true;
0660 }
0661 break;
0662 case 1:
0663 if (locBz > 0.f) {
0664 flip_y_ = true;
0665 }
0666 break;
0667 case 2:
0668 case 3:
0669 case 4:
0670 case 5:
0671 if (locBx * locBz < 0.f) {
0672 cota = std::abs(cotalpha);
0673 flip_x_ = true;
0674 }
0675 if (locBx < 0.f) {
0676 flip_y_ = true;
0677 }
0678 break;
0679 default:
0680 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0681 throw cms::Exception("DataCorrupt")
0682 << "SiPixelTemplate2D::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
0683
0684
0685 if (!edm::isFinite(cotalpha) || !edm::isFinite(cotbeta)) {
0686 success_ = false;
0687 return success_;
0688 }
0689 #else
0690 std::cout << "SiPixelTemplate:2D:illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
0691 #endif
0692 }
0693
0694 if (cota < cotalpha0_) {
0695 success_ = false;
0696 jx0_ = 0;
0697 jx1_ = 1;
0698 adcota_ = 0.f;
0699 } else if (cota > cotalpha1_) {
0700 success_ = false;
0701 jx0_ = Nxx_ - 1;
0702 jx1_ = jx0_ - 1;
0703 adcota_ = 0.f;
0704 } else {
0705 jx0_ = (int)((cota - cotalpha0_) / deltacota_ + 0.5f);
0706 dcota = (cota - (cotalpha0_ + jx0_ * deltacota_)) / deltacota_;
0707 adcota_ = fabs(dcota);
0708 if (dcota > 0.f) {
0709 jx1_ = jx0_ + 1;
0710 if (jx1_ > Nxx_ - 1)
0711 jx1_ = jx0_ - 1;
0712 } else {
0713 jx1_ = jx0_ - 1;
0714 if (jx1_ < 0)
0715 jx1_ = jx0_ + 1;
0716 }
0717 }
0718
0719
0720
0721 acotb = std::abs(cotbeta);
0722
0723 if (acotb < cotbeta0_) {
0724 success_ = false;
0725 iy0_ = 0;
0726 iy1_ = 1;
0727 adcotb_ = 0.f;
0728 } else if (acotb > cotbeta1_) {
0729 success_ = false;
0730 iy0_ = Nyx_ - 1;
0731 iy1_ = iy0_ - 1;
0732 adcotb_ = 0.f;
0733 } else {
0734 iy0_ = (int)((acotb - cotbeta0_) / deltacotb_ + 0.5f);
0735 dcotb = (acotb - (cotbeta0_ + iy0_ * deltacotb_)) / deltacotb_;
0736 adcotb_ = fabs(dcotb);
0737 if (dcotb > 0.f) {
0738 iy1_ = iy0_ + 1;
0739 if (iy1_ > Nyx_ - 1)
0740 iy1_ = iy0_ - 1;
0741 } else {
0742 iy1_ = iy0_ - 1;
0743 if (iy1_ < 0)
0744 iy1_ = iy0_ + 1;
0745 }
0746 }
0747
0748
0749
0750 lorydrift_ = lorywidth_ / 2.;
0751 if (flip_y_)
0752 lorydrift_ = -lorydrift_;
0753 lorxdrift_ = lorxwidth_ / 2.;
0754 if (flip_x_)
0755 lorxdrift_ = -lorxdrift_;
0756
0757
0758
0759 entry00_ = &thePixelTemp_[index_id_].entry[iy0_][jx0_];
0760 entry10_ = &thePixelTemp_[index_id_].entry[iy1_][jx0_];
0761 entry01_ = &thePixelTemp_[index_id_].entry[iy0_][jx1_];
0762
0763
0764
0765 qavg_ = entry00_->qavg + adcota_ * (entry01_->qavg - entry00_->qavg) + adcotb_ * (entry10_->qavg - entry00_->qavg);
0766
0767 pixmax_ = entry00_->pixmax + adcota_ * (entry01_->pixmax - entry00_->pixmax) +
0768 adcotb_ * (entry10_->pixmax - entry00_->pixmax);
0769
0770 sxymax_ = entry00_->sxymax + adcota_ * (entry01_->sxymax - entry00_->sxymax) +
0771 adcotb_ * (entry10_->sxymax - entry00_->sxymax);
0772
0773 chi2avgone_ = entry00_->chi2avgone + adcota_ * (entry01_->chi2avgone - entry00_->chi2avgone) +
0774 adcotb_ * (entry10_->chi2avgone - entry00_->chi2avgone);
0775
0776 chi2minone_ = entry00_->chi2minone + adcota_ * (entry01_->chi2minone - entry00_->chi2minone) +
0777 adcotb_ * (entry10_->chi2minone - entry00_->chi2minone);
0778
0779 clsleny_ = entry00_->clsleny + adcota_ * (entry01_->clsleny - entry00_->clsleny) +
0780 adcotb_ * (entry10_->clsleny - entry00_->clsleny);
0781
0782 clslenx_ = entry00_->clslenx + adcota_ * (entry01_->clslenx - entry00_->clslenx) +
0783 adcotb_ * (entry10_->clslenx - entry00_->clslenx);
0784
0785 chi2ppix_ = entry00_->chi2ppix + adcota_ * (entry01_->chi2ppix - entry00_->chi2ppix) +
0786 adcotb_ * (entry10_->chi2ppix - entry00_->chi2ppix);
0787
0788 chi2scale_ = entry00_->chi2scale + adcota_ * (entry01_->chi2scale - entry00_->chi2scale) +
0789 adcotb_ * (entry10_->chi2scale - entry00_->chi2scale);
0790
0791 scaleyavg_ = entry00_->scaleyavg + adcota_ * (entry01_->scaleyavg - entry00_->scaleyavg) +
0792 adcotb_ * (entry10_->scaleyavg - entry00_->scaleyavg);
0793
0794 scalexavg_ = entry00_->scalexavg + adcota_ * (entry01_->scalexavg - entry00_->scalexavg) +
0795 adcotb_ * (entry10_->scalexavg - entry00_->scalexavg);
0796
0797 delyavg_ = entry00_->delyavg + adcota_ * (entry01_->delyavg - entry00_->delyavg) +
0798 adcotb_ * (entry10_->delyavg - entry00_->delyavg);
0799
0800 delysig_ = entry00_->delysig + adcota_ * (entry01_->delysig - entry00_->delysig) +
0801 adcotb_ * (entry10_->delysig - entry00_->delysig);
0802
0803 mpvvav_ = entry00_->mpvvav + adcota_ * (entry01_->mpvvav - entry00_->mpvvav) +
0804 adcotb_ * (entry10_->mpvvav - entry00_->mpvvav);
0805
0806 sigmavav_ = entry00_->sigmavav + adcota_ * (entry01_->sigmavav - entry00_->sigmavav) +
0807 adcotb_ * (entry10_->sigmavav - entry00_->sigmavav);
0808
0809 kappavav_ = entry00_->kappavav + adcota_ * (entry01_->kappavav - entry00_->kappavav) +
0810 adcotb_ * (entry10_->kappavav - entry00_->kappavav);
0811
0812 for (int i = 0; i < 4; ++i) {
0813 scalex_[i] = entry00_->scalex[i] + adcota_ * (entry01_->scalex[i] - entry00_->scalex[i]) +
0814 adcotb_ * (entry10_->scalex[i] - entry00_->scalex[i]);
0815
0816 scaley_[i] = entry00_->scaley[i] + adcota_ * (entry01_->scaley[i] - entry00_->scaley[i]) +
0817 adcotb_ * (entry10_->scaley[i] - entry00_->scaley[i]);
0818
0819 offsetx_[i] = entry00_->offsetx[i] + adcota_ * (entry01_->offsetx[i] - entry00_->offsetx[i]) +
0820 adcotb_ * (entry10_->offsetx[i] - entry00_->offsetx[i]);
0821 if (flip_x_)
0822 offsetx_[i] = -offsetx_[i];
0823
0824 offsety_[i] = entry00_->offsety[i] + adcota_ * (entry01_->offsety[i] - entry00_->offsety[i]) +
0825 adcotb_ * (entry10_->offsety[i] - entry00_->offsety[i]);
0826 if (flip_y_)
0827 offsety_[i] = -offsety_[i];
0828 }
0829
0830 for (int i = 0; i < 2; ++i) {
0831 for (int j = 0; j < 5; ++j) {
0832
0833 if (flip_y_) {
0834 xypary0x0_[1 - i][j] = (float)entry00_->xypar[i][j];
0835 xypary1x0_[1 - i][j] = (float)entry10_->xypar[i][j];
0836 xypary0x1_[1 - i][j] = (float)entry01_->xypar[i][j];
0837 lanpar_[1 - i][j] = entry00_->lanpar[i][j] + adcota_ * (entry01_->lanpar[i][j] - entry00_->lanpar[i][j]) +
0838 adcotb_ * (entry10_->lanpar[i][j] - entry00_->lanpar[i][j]);
0839 } else {
0840 xypary0x0_[i][j] = (float)entry00_->xypar[i][j];
0841 xypary1x0_[i][j] = (float)entry10_->xypar[i][j];
0842 xypary0x1_[i][j] = (float)entry01_->xypar[i][j];
0843 lanpar_[i][j] = entry00_->lanpar[i][j] + adcota_ * (entry01_->lanpar[i][j] - entry00_->lanpar[i][j]) +
0844 adcotb_ * (entry10_->lanpar[i][j] - entry00_->lanpar[i][j]);
0845 }
0846 }
0847 }
0848
0849 return success_;
0850 }
0851
0852
0853
0854
0855
0856
0857
0858
0859
0860 #ifdef SI_PIXEL_TEMPLATE_STANDALONE
0861 void SiPixelTemplate2D::sideload(SiPixelTemplateEntry2D* entry,
0862 int iDtype,
0863 float locBx,
0864 float locBz,
0865 float lorwdy,
0866 float lorwdx,
0867 float q50,
0868 float fbin[3],
0869 float xsize,
0870 float ysize,
0871 float zsize) {
0872
0873
0874 entry00_ = entry;
0875 entry01_ = entry;
0876 entry10_ = entry;
0877 Dtype_ = iDtype;
0878 lorywidth_ = lorwdy;
0879 lorxwidth_ = lorwdx;
0880 xsize_ = xsize;
0881 ysize_ = ysize;
0882 zsize_ = zsize;
0883 s50_ = q50;
0884 qscale_ = 1.f;
0885 for (int i = 0; i < 3; ++i) {
0886 fbin_[i] = fbin[i];
0887 }
0888
0889
0890
0891 adcota_ = 0.f;
0892 adcotb_ = 0.f;
0893
0894
0895
0896 qavg_ = entry00_->qavg;
0897
0898 pixmax_ = entry00_->pixmax;
0899
0900 sxymax_ = entry00_->sxymax;
0901
0902 clsleny_ = entry00_->clsleny;
0903
0904 clslenx_ = entry00_->clslenx;
0905
0906 scaleyavg_ = 1.f;
0907
0908 scalexavg_ = 1.f;
0909
0910 delyavg_ = 0.f;
0911
0912 delysig_ = 0.f;
0913
0914 for (int i = 0; i < 4; ++i) {
0915 scalex_[i] = 1.f;
0916 scaley_[i] = 1.f;
0917 offsetx_[i] = 0.f;
0918 offsety_[i] = 0.f;
0919 }
0920
0921
0922
0923 flip_x_ = false;
0924 flip_y_ = false;
0925 float cotbeta = entry00_->cotbeta;
0926 switch (Dtype_) {
0927 case 0:
0928 if (cotbeta < 0.f) {
0929 flip_y_ = true;
0930 }
0931 break;
0932 case 1:
0933 if (locBz > 0.f) {
0934 flip_y_ = true;
0935 }
0936 break;
0937 case 2:
0938 case 3:
0939 case 4:
0940 case 5:
0941 if (locBx * locBz < 0.f) {
0942 flip_x_ = true;
0943 }
0944 if (locBx < 0.f) {
0945 flip_y_ = true;
0946 }
0947 break;
0948 default:
0949 std::cout << "SiPixelTemplate:2D:illegal subdetector ID = " << iDtype << std::endl;
0950 }
0951
0952
0953
0954 lorydrift_ = lorywidth_ / 2.;
0955 if (flip_y_)
0956 lorydrift_ = -lorydrift_;
0957 lorxdrift_ = lorxwidth_ / 2.;
0958 if (flip_x_)
0959 lorxdrift_ = -lorxdrift_;
0960
0961 for (int i = 0; i < 2; ++i) {
0962 for (int j = 0; j < 5; ++j) {
0963
0964 if (flip_y_) {
0965 xypary0x0_[1 - i][j] = (float)entry00_->xypar[i][j];
0966 xypary1x0_[1 - i][j] = (float)entry00_->xypar[i][j];
0967 xypary0x1_[1 - i][j] = (float)entry00_->xypar[i][j];
0968 lanpar_[1 - i][j] = entry00_->lanpar[i][j];
0969 } else {
0970 xypary0x0_[i][j] = (float)entry00_->xypar[i][j];
0971 xypary1x0_[i][j] = (float)entry00_->xypar[i][j];
0972 xypary0x1_[i][j] = (float)entry00_->xypar[i][j];
0973 lanpar_[i][j] = entry00_->lanpar[i][j];
0974 }
0975 }
0976 }
0977 return;
0978 }
0979 #endif
0980
0981
0982
0983
0984
0985
0986
0987
0988
0989 bool SiPixelTemplate2D::xytemp(float xhit,
0990 float yhit,
0991 bool ydouble[BYM2],
0992 bool xdouble[BXM2],
0993 float template2d[BXM2][BYM2],
0994 bool derivatives,
0995 float dpdx2d[2][BXM2][BYM2],
0996 float& QTemplate) {
0997
0998
0999
1000 int pixx, pixy, k0, k1, l0, l1, deltax, deltay, iflipy, jflipx, imin, imax, jmin, jmax;
1001 int m, n;
1002 float dx, dy, ddx, ddy, adx, ady;
1003
1004 const float deltaxy[2] = {16.67f, 25.0f};
1005
1006
1007
1008
1009
1010
1011
1012
1013 pixy = (int)floorf(yhit / ysize_);
1014 dy = yhit - (pixy + 0.5f) * ysize_;
1015 if (flip_y_) {
1016 dy = -dy;
1017 }
1018 k0 = (int)(dy / ysize_ * 6.f + 3.5f);
1019 if (k0 < 0)
1020 k0 = 0;
1021 if (k0 > 6)
1022 k0 = 6;
1023 ddy = 6.f * dy / ysize_ - (k0 - 3);
1024 ady = fabs(ddy);
1025 if (ddy > 0.f) {
1026 k1 = k0 + 1;
1027 if (k1 > 6)
1028 k1 = k0 - 1;
1029 } else {
1030 k1 = k0 - 1;
1031 if (k1 < 0)
1032 k1 = k0 + 1;
1033 }
1034 pixx = (int)floorf(xhit / xsize_);
1035 dx = xhit - (pixx + 0.5f) * xsize_;
1036 if (flip_x_) {
1037 dx = -dx;
1038 }
1039 l0 = (int)(dx / xsize_ * 6.f + 3.5f);
1040 if (l0 < 0)
1041 l0 = 0;
1042 if (l0 > 6)
1043 l0 = 6;
1044 ddx = 6.f * dx / xsize_ - (l0 - 3);
1045 adx = fabs(ddx);
1046 if (ddx > 0.f) {
1047 l1 = l0 + 1;
1048 if (l1 > 6)
1049 l1 = l0 - 1;
1050 } else {
1051 l1 = l0 - 1;
1052 if (l1 < 0)
1053 l1 = l0 + 1;
1054 }
1055
1056
1057
1058
1059
1060 imin = std::min(entry00_->iymin, entry10_->iymin);
1061 imin_ = std::min(imin, entry01_->iymin);
1062
1063 jmin = std::min(entry00_->jxmin, entry10_->jxmin);
1064 jmin_ = std::min(jmin, entry01_->jxmin);
1065
1066 imax = std::max(entry00_->iymax, entry10_->iymax);
1067 imax_ = std::max(imax, entry01_->iymax);
1068
1069 jmax = std::max(entry00_->jxmax, entry10_->jxmax);
1070 jmax_ = std::max(jmax, entry01_->jxmax);
1071
1072
1073
1074
1075
1076 ++pixy;
1077 ++pixx;
1078
1079
1080
1081 deltax = pixx - T2HX;
1082 deltay = pixy - T2HY;
1083
1084
1085
1086 for (int j = 0; j < BXM2; ++j) {
1087 for (int i = 0; i < BYM2; ++i) {
1088 xytemp_[j][i] = 0.f;
1089 }
1090 }
1091
1092
1093
1094 for (int j = jmin_; j <= jmax_; ++j) {
1095
1096 if (flip_x_) {
1097 jflipx = T2XSIZE - 1 - j;
1098 m = deltax + jflipx;
1099 } else {
1100 m = deltax + j;
1101 }
1102 for (int i = imin_; i <= imax_; ++i) {
1103 if (flip_y_) {
1104 iflipy = T2YSIZE - 1 - i;
1105 n = deltay + iflipy;
1106 } else {
1107 n = deltay + i;
1108 }
1109 if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1110 xytemp_[m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1111 adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1112 ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1113 adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1114 adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1115 }
1116 }
1117 }
1118
1119
1120
1121 for (int n = 1; n < BYM3; ++n) {
1122 if (ydouble[n]) {
1123
1124 for (int m = 1; m < BXM3; ++m) {
1125 xytemp_[m][n] += xytemp_[m][n + 1];
1126 }
1127
1128 for (int i = n + 1; i < BYM3; ++i) {
1129 for (int m = 1; m < BXM3; ++m) {
1130 xytemp_[m][i] = xytemp_[m][i + 1];
1131 }
1132 }
1133 }
1134 }
1135
1136
1137
1138 for (int m = 1; m < BXM3; ++m) {
1139 if (xdouble[m]) {
1140
1141 for (int n = 1; n < BYM3; ++n) {
1142 xytemp_[m][n] += xytemp_[m + 1][n];
1143 }
1144
1145 for (int j = m + 1; j < BXM3; ++j) {
1146 for (n = 1; n < BYM3; ++n) {
1147 xytemp_[j][n] = xytemp_[j + 1][n];
1148 }
1149 }
1150 }
1151 }
1152
1153
1154
1155 float qtemptot = 0.f;
1156
1157 for (int n = 1; n < BYM3; ++n) {
1158 for (int m = 1; m < BXM3; ++m) {
1159 if (xytemp_[m][n] != 0.f) {
1160 template2d[m][n] += xytemp_[m][n];
1161 qtemptot += xytemp_[m][n];
1162 }
1163 }
1164 }
1165
1166 QTemplate = qtemptot;
1167
1168 if (derivatives) {
1169 float dxytempdx[2][BXM2][BYM2], dxytempdy[2][BXM2][BYM2];
1170
1171 for (int k = 0; k < 2; ++k) {
1172 for (int i = 0; i < BXM2; ++i) {
1173 for (int j = 0; j < BYM2; ++j) {
1174 dxytempdx[k][i][j] = 0.f;
1175 dxytempdy[k][i][j] = 0.f;
1176 dpdx2d[k][i][j] = 0.f;
1177 }
1178 }
1179 }
1180
1181
1182
1183 pixx = (int)floorf((xhit + deltaxy[0]) / xsize_);
1184 dx = (xhit + deltaxy[0]) - (pixx + 0.5f) * xsize_;
1185 if (flip_x_) {
1186 dx = -dx;
1187 }
1188 l0 = (int)(dx / xsize_ * 6.f + 3.5f);
1189 if (l0 < 0)
1190 l0 = 0;
1191 if (l0 > 6)
1192 l0 = 6;
1193 ddx = 6.f * dx / xsize_ - (l0 - 3);
1194 adx = fabs(ddx);
1195 if (ddx > 0.f) {
1196 l1 = l0 + 1;
1197 if (l1 > 6)
1198 l1 = l0 - 1;
1199 } else {
1200 l1 = l0 - 1;
1201 if (l1 < 0)
1202 l1 = l0 + 1;
1203 }
1204
1205
1206
1207
1208
1209
1210
1211 ++pixx;
1212
1213
1214
1215 deltax = pixx - T2HX;
1216
1217
1218
1219 for (int j = jmin_; j <= jmax_; ++j) {
1220
1221 if (flip_x_) {
1222 jflipx = T2XSIZE - 1 - j;
1223 m = deltax + jflipx;
1224 } else {
1225 m = deltax + j;
1226 }
1227 for (int i = imin_; i <= imax_; ++i) {
1228 if (flip_y_) {
1229 iflipy = T2YSIZE - 1 - i;
1230 n = deltay + iflipy;
1231 } else {
1232 n = deltay + i;
1233 }
1234 if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1235 dxytempdx[1][m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1236 adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1237 ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1238 adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1239 adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1240 }
1241 }
1242 }
1243
1244
1245
1246 for (int n = 1; n < BYM3; ++n) {
1247 if (ydouble[n]) {
1248
1249 for (int m = 1; m < BXM3; ++m) {
1250 dxytempdx[1][m][n] += dxytempdx[1][m][n + 1];
1251 }
1252
1253 for (int i = n + 1; i < BYM3; ++i) {
1254 for (int m = 1; m < BXM3; ++m) {
1255 dxytempdx[1][m][i] = dxytempdx[1][m][i + 1];
1256 }
1257 }
1258 }
1259 }
1260
1261
1262
1263 for (int m = 1; m < BXM3; ++m) {
1264 if (xdouble[m]) {
1265
1266 for (int n = 1; n < BYM3; ++n) {
1267 dxytempdx[1][m][n] += dxytempdx[1][m + 1][n];
1268 }
1269
1270 for (int j = m + 1; j < BXM3; ++j) {
1271 for (int n = 1; n < BYM3; ++n) {
1272 dxytempdx[1][j][n] = dxytempdx[1][j + 1][n];
1273 }
1274 }
1275 }
1276 }
1277
1278
1279
1280 pixx = (int)floorf((xhit - deltaxy[0]) / xsize_);
1281 dx = (xhit - deltaxy[0]) - (pixx + 0.5f) * xsize_;
1282 if (flip_x_) {
1283 dx = -dx;
1284 }
1285 l0 = (int)(dx / xsize_ * 6.f + 3.5f);
1286 if (l0 < 0)
1287 l0 = 0;
1288 if (l0 > 6)
1289 l0 = 6;
1290 ddx = 6.f * dx / xsize_ - (l0 - 3);
1291 adx = fabs(ddx);
1292 if (ddx > 0.f) {
1293 l1 = l0 + 1;
1294 if (l1 > 6)
1295 l1 = l0 - 1;
1296 } else {
1297 l1 = l0 - 1;
1298 if (l1 < 0)
1299 l1 = l0 + 1;
1300 }
1301
1302
1303
1304
1305
1306
1307
1308 ++pixx;
1309
1310
1311
1312 deltax = pixx - T2HX;
1313
1314
1315
1316 for (int j = jmin_; j <= jmax_; ++j) {
1317
1318 if (flip_x_) {
1319 jflipx = T2XSIZE - 1 - j;
1320 m = deltax + jflipx;
1321 } else {
1322 m = deltax + j;
1323 }
1324 for (int i = imin_; i <= imax_; ++i) {
1325 if (flip_y_) {
1326 iflipy = T2YSIZE - 1 - i;
1327 n = deltay + iflipy;
1328 } else {
1329 n = deltay + i;
1330 }
1331 if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1332 dxytempdx[0][m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1333 adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1334 ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1335 adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1336 adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1337 }
1338 }
1339 }
1340
1341
1342
1343 for (int n = 1; n < BYM3; ++n) {
1344 if (ydouble[n]) {
1345
1346 for (int m = 1; m < BXM3; ++m) {
1347 dxytempdx[0][m][n] += dxytempdx[0][m][n + 1];
1348 }
1349
1350 for (int i = n + 1; i < BYM3; ++i) {
1351 for (int m = 1; m < BXM3; ++m) {
1352 dxytempdx[0][m][i] = dxytempdx[0][m][i + 1];
1353 }
1354 }
1355 }
1356 }
1357
1358
1359
1360 for (int m = 1; m < BXM3; ++m) {
1361 if (xdouble[m]) {
1362
1363 for (int n = 1; n < BYM3; ++n) {
1364 dxytempdx[0][m][n] += dxytempdx[0][m + 1][n];
1365 }
1366
1367 for (int j = m + 1; j < BXM3; ++j) {
1368 for (int n = 1; n < BYM3; ++n) {
1369 dxytempdx[0][j][n] = dxytempdx[0][j + 1][n];
1370 }
1371 }
1372 }
1373 }
1374
1375
1376
1377 for (int n = 1; n < BYM3; ++n) {
1378 for (int m = 1; m < BXM3; ++m) {
1379 dpdx2d[0][m][n] = (dxytempdx[1][m][n] - dxytempdx[0][m][n]) / (2. * deltaxy[0]);
1380 }
1381 }
1382
1383
1384
1385 pixy = (int)floorf((yhit + deltaxy[1]) / ysize_);
1386 dy = (yhit + deltaxy[1]) - (pixy + 0.5f) * ysize_;
1387 if (flip_y_) {
1388 dy = -dy;
1389 }
1390 k0 = (int)(dy / ysize_ * 6.f + 3.5f);
1391 if (k0 < 0)
1392 k0 = 0;
1393 if (k0 > 6)
1394 k0 = 6;
1395 ddy = 6.f * dy / ysize_ - (k0 - 3);
1396 ady = fabs(ddy);
1397 if (ddy > 0.f) {
1398 k1 = k0 + 1;
1399 if (k1 > 6)
1400 k1 = k0 - 1;
1401 } else {
1402 k1 = k0 - 1;
1403 if (k1 < 0)
1404 k1 = k0 + 1;
1405 }
1406 pixx = (int)floorf(xhit / xsize_);
1407 dx = xhit - (pixx + 0.5f) * xsize_;
1408 if (flip_x_) {
1409 dx = -dx;
1410 }
1411 l0 = (int)(dx / xsize_ * 6.f + 3.5f);
1412 if (l0 < 0)
1413 l0 = 0;
1414 if (l0 > 6)
1415 l0 = 6;
1416 ddx = 6.f * dx / xsize_ - (l0 - 3);
1417 adx = fabs(ddx);
1418 if (ddx > 0.f) {
1419 l1 = l0 + 1;
1420 if (l1 > 6)
1421 l1 = l0 - 1;
1422 } else {
1423 l1 = l0 - 1;
1424 if (l1 < 0)
1425 l1 = l0 + 1;
1426 }
1427
1428
1429
1430
1431
1432
1433
1434 ++pixy;
1435 ++pixx;
1436
1437
1438
1439 deltax = pixx - T2HX;
1440 deltay = pixy - T2HY;
1441
1442
1443
1444 for (int j = jmin_; j <= jmax_; ++j) {
1445
1446 if (flip_x_) {
1447 jflipx = T2XSIZE - 1 - j;
1448 m = deltax + jflipx;
1449 } else {
1450 m = deltax + j;
1451 }
1452 for (int i = imin_; i <= imax_; ++i) {
1453 if (flip_y_) {
1454 iflipy = T2YSIZE - 1 - i;
1455 n = deltay + iflipy;
1456 } else {
1457 n = deltay + i;
1458 }
1459 if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1460 dxytempdy[1][m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1461 adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1462 ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1463 adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1464 adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1465 }
1466 }
1467 }
1468
1469
1470
1471 for (int n = 1; n < BYM3; ++n) {
1472 if (ydouble[n]) {
1473
1474 for (int m = 1; m < BXM3; ++m) {
1475 dxytempdy[1][m][n] += dxytempdy[1][m][n + 1];
1476 }
1477
1478 for (int i = n + 1; i < BYM3; ++i) {
1479 for (int m = 1; m < BXM3; ++m) {
1480 dxytempdy[1][m][i] = dxytempdy[1][m][i + 1];
1481 }
1482 }
1483 }
1484 }
1485
1486
1487
1488 for (int m = 1; m < BXM3; ++m) {
1489 if (xdouble[m]) {
1490
1491 for (int n = 1; n < BYM3; ++n) {
1492 dxytempdy[1][m][n] += dxytempdy[1][m + 1][n];
1493 }
1494
1495 for (int j = m + 1; j < BXM3; ++j) {
1496 for (int n = 1; n < BYM3; ++n) {
1497 dxytempdy[1][j][n] = dxytempdy[1][j + 1][n];
1498 }
1499 }
1500 }
1501 }
1502
1503
1504
1505 pixy = (int)floorf((yhit - deltaxy[1]) / ysize_);
1506 dy = (yhit - deltaxy[1]) - (pixy + 0.5f) * ysize_;
1507 if (flip_y_) {
1508 dy = -dy;
1509 }
1510 k0 = (int)(dy / ysize_ * 6.f + 3.5f);
1511 if (k0 < 0)
1512 k0 = 0;
1513 if (k0 > 6)
1514 k0 = 6;
1515 ddy = 6.f * dy / ysize_ - (k0 - 3);
1516 ady = fabs(ddy);
1517 if (ddy > 0.f) {
1518 k1 = k0 + 1;
1519 if (k1 > 6)
1520 k1 = k0 - 1;
1521 } else {
1522 k1 = k0 - 1;
1523 if (k1 < 0)
1524 k1 = k0 + 1;
1525 }
1526
1527
1528
1529
1530
1531
1532
1533 ++pixy;
1534
1535
1536
1537 deltay = pixy - T2HY;
1538
1539
1540
1541 for (int j = jmin_; j <= jmax_; ++j) {
1542
1543 if (flip_x_) {
1544 jflipx = T2XSIZE - 1 - j;
1545 m = deltax + jflipx;
1546 } else {
1547 m = deltax + j;
1548 }
1549 for (int i = imin_; i <= imax_; ++i) {
1550 if (flip_y_) {
1551 iflipy = T2YSIZE - 1 - i;
1552 n = deltay + iflipy;
1553 } else {
1554 n = deltay + i;
1555 }
1556 if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1557 dxytempdy[0][m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1558 adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1559 ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1560 adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1561 adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1562 }
1563 }
1564 }
1565
1566
1567
1568 for (int n = 1; n < BYM3; ++n) {
1569 if (ydouble[n]) {
1570
1571 for (int m = 1; m < BXM3; ++m) {
1572 dxytempdy[0][m][n] += dxytempdy[0][m][n + 1];
1573 }
1574
1575 for (int i = n + 1; i < BYM3; ++i) {
1576 for (int m = 1; m < BXM3; ++m) {
1577 dxytempdy[0][m][i] = dxytempdy[0][m][i + 1];
1578 }
1579 }
1580 }
1581 }
1582
1583
1584
1585 for (int m = 1; m < BXM3; ++m) {
1586 if (xdouble[m]) {
1587
1588 for (int n = 1; n < BYM3; ++n) {
1589 dxytempdy[0][m][n] += dxytempdy[0][m + 1][n];
1590 }
1591
1592 for (int j = m + 1; j < BXM3; ++j) {
1593 for (int n = 1; n < BYM3; ++n) {
1594 dxytempdy[0][j][n] = dxytempdy[0][j + 1][n];
1595 }
1596 }
1597 }
1598 }
1599
1600
1601
1602 for (int n = 1; n < BYM3; ++n) {
1603 for (int m = 1; m < BXM3; ++m) {
1604 dpdx2d[1][m][n] = (dxytempdy[1][m][n] - dxytempdy[0][m][n]) / (2. * deltaxy[1]);
1605 }
1606 }
1607 }
1608
1609 return success_;
1610 }
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621 bool SiPixelTemplate2D::xytemp(
1622 float xhit, float yhit, bool ydouble[BYM2], bool xdouble[BXM2], float template2d[BXM2][BYM2]) {
1623
1624
1625 bool derivatives = false;
1626 float dpdx2d[2][BXM2][BYM2];
1627 float QTemplate;
1628
1629 return SiPixelTemplate2D::xytemp(xhit, yhit, ydouble, xdouble, template2d, derivatives, dpdx2d, QTemplate);
1630
1631 }
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645 bool SiPixelTemplate2D::xytemp(int id,
1646 float cotalpha,
1647 float cotbeta,
1648 float xhit,
1649 float yhit,
1650 std::vector<bool>& ydouble,
1651 std::vector<bool>& xdouble,
1652 float template2d[BXM2][BYM2]) {
1653
1654
1655 bool derivatives = false;
1656 float dpdx2d[2][BXM2][BYM2];
1657 float QTemplate;
1658 float locBx = 1.f;
1659 if (cotbeta < 0.f) {
1660 locBx = -1.f;
1661 }
1662 float locBz = locBx;
1663 if (cotalpha < 0.f) {
1664 locBz = -locBx;
1665 }
1666
1667 bool yd[BYM2], xd[BXM2];
1668
1669 yd[0] = false;
1670 yd[BYM2 - 1] = false;
1671 for (int i = 0; i < TYSIZE; ++i) {
1672 yd[i + 1] = ydouble[i];
1673 }
1674 xd[0] = false;
1675 xd[BXM2 - 1] = false;
1676 for (int j = 0; j < TXSIZE; ++j) {
1677 xd[j + 1] = xdouble[j];
1678 }
1679
1680
1681
1682 if (SiPixelTemplate2D::interpolate(id, cotalpha, cotbeta, locBz, locBx)) {
1683 return SiPixelTemplate2D::xytemp(xhit, yhit, yd, xd, template2d, derivatives, dpdx2d, QTemplate);
1684 } else {
1685 return false;
1686 }
1687
1688 }
1689
1690
1691
1692
1693
1694
1695
1696
1697 void SiPixelTemplate2D::xysigma2(float qpixel, int index, float& xysig2)
1698
1699 {
1700
1701
1702
1703 float sigi, sigi2, sigi3, sigi4, qscale, err2, err00;
1704
1705
1706
1707 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1708 if (index < 1 || index >= BYM2) {
1709 throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::ysigma2 called with index = " << index << std::endl;
1710 }
1711 #else
1712 assert(index > 0 && index < BYM2);
1713 #endif
1714
1715
1716
1717
1718
1719 if (qpixel < sxymax_) {
1720 sigi = qpixel;
1721 qscale = 1.f;
1722 } else {
1723 sigi = sxymax_;
1724 qscale = qpixel / sxymax_;
1725 }
1726 sigi2 = sigi * sigi;
1727 sigi3 = sigi2 * sigi;
1728 sigi4 = sigi3 * sigi;
1729 if (index <= T2HYP1) {
1730 err00 = xypary0x0_[0][0] + xypary0x0_[0][1] * sigi + xypary0x0_[0][2] * sigi2 + xypary0x0_[0][3] * sigi3 +
1731 xypary0x0_[0][4] * sigi4;
1732 err2 = err00 +
1733 adcota_ * (xypary0x1_[0][0] + xypary0x1_[0][1] * sigi + xypary0x1_[0][2] * sigi2 + xypary0x1_[0][3] * sigi3 +
1734 xypary0x1_[0][4] * sigi4 - err00) +
1735 adcotb_ * (xypary1x0_[0][0] + xypary1x0_[0][1] * sigi + xypary1x0_[0][2] * sigi2 + xypary1x0_[0][3] * sigi3 +
1736 xypary1x0_[0][4] * sigi4 - err00);
1737 } else {
1738 err00 = xypary0x0_[1][0] + xypary0x0_[1][1] * sigi + xypary0x0_[1][2] * sigi2 + xypary0x0_[1][3] * sigi3 +
1739 xypary0x0_[1][4] * sigi4;
1740 err2 = err00 +
1741 adcota_ * (xypary0x1_[1][0] + xypary0x1_[1][1] * sigi + xypary0x1_[1][2] * sigi2 + xypary0x1_[1][3] * sigi3 +
1742 xypary0x1_[1][4] * sigi4 - err00) +
1743 adcotb_ * (xypary1x0_[1][0] + xypary1x0_[1][1] * sigi + xypary1x0_[1][2] * sigi2 + xypary1x0_[1][3] * sigi3 +
1744 xypary1x0_[1][4] * sigi4 - err00);
1745 }
1746 xysig2 = qscale * err2;
1747 if (xysig2 <= 0.f) {
1748 xysig2 = s50_ * s50_;
1749 }
1750
1751 return;
1752
1753 }
1754
1755
1756
1757
1758 void SiPixelTemplate2D::landau_par(float lanpar[2][5])
1759
1760 {
1761
1762
1763 for (int i = 0; i < 2; ++i) {
1764 for (int j = 0; j < 5; ++j) {
1765 lanpar[i][j] = lanpar_[i][j];
1766 }
1767 }
1768 return;
1769
1770 }