SiPixelTransient/src/SiPixelTemplate.cc

0001 //
0002 //  SiPixelTemplate.cc  Version 10.24
0003 //
0004 //  Add goodness-of-fit info and spare entries to templates, version number in template header, more error checking
0005 //  Add correction for (Q_F-Q_L)/(Q_F+Q_L) bias
0006 //  Add cot(beta) reflection to reduce y-entries and more sophisticated x-interpolation
0007 //  Fix small index searching bug in interpolate method
0008 //  Change interpolation indexing to avoid complier complaining about possible un-initialized variables
0009 //  Replace containers with static arrays in calls to ysigma2 and xsigma2
0010 //  Add external threshold to calls to ysigma2 and xsigma2, fix parameter signal max for xsigma2
0011 //  Return to 5 pixel spanning but adjust boundaries to use only when needed
0012 //  Implement improved (faster) chi2min search that depends on pixel types
0013 //  Fill template arrays in single calls to this object
0014 //  Add qmin to the template
0015 //  Add qscale to match charge scales
0016 //  Small improvement to x-chisquare interpolation
0017 //  Enlarge SiPixelTemplateStore to accommodate larger templates and increased alpha acceptance (reduce PT threshold to ~200 MeV)
0018 //  Change output streams to conform to CMSSW info and error logging.
0019 //  Store x and y cluster sizes in fractional pixels to facilitate cluster splitting
0020 //  Add methods to return 3-d templates needed for cluster splitting
0021 //  Keep interpolated central 9 template bins in private storage and expand/shift in the getter functions (faster for speed=2/3) and easier to build 3d templates
0022 //  Store error and bias information for the simple chi^2 min position analysis (no interpolation or Q_{FB} corrections) to use in cluster splitting
0023 //  To save time, the gaussian centers and sigma are not interpolated right now (they aren't currently used).  They can be restored by un-commenting lines in the interpolate method.
0024 //  Add a new method to calculate qbin for input cotbeta and cluster charge.  To be used for error estimation of merged clusters in PixelCPEGeneric.
0025 //  Improve the charge estimation for larger cot(alpha) tracks
0026 //  Change interpolate method to return false boolean if track angles are outside of range
0027 //  Add template info and method for truncation information
0028 //  Change to allow template sizes to be changed at compile time
0029 //  Fix bug in track angle checking
0030 //  Accommodate Dave's new DB pushfile which overloads the old method (file input)
0031 //  Add CPEGeneric error information and expand qbin method to access useful info for PixelCPEGeneric
0032 //  Fix large cot(alpha) bug in qmin interpolation
0033 //  Add second qmin to allow a qbin=5 state
0034 //  Use interpolated chi^2 info for one-pixel clusters
0035 //  Fix DB pushfile version number checking bug.
0036 //  Remove assert from qbin method
0037 //  Replace asserts with exceptions in CMSSW
0038 //  Change calling sequence to interpolate method to handle cot(beta)<0 for FPix cosmics
0039 //  Add getter for pixelav Lorentz width estimates to qbin method
0040 //  Add check on template size to interpolate and qbin methods
0041 //  Add qbin population information, charge distribution information
0042 //
0043 //
0044 //  V7.00 - Decouple BPix and FPix information into separate templates
0045 //  Add methods to facilitate improved cluster splitting
0046 //  Fix small charge scaling bug (affects FPix only)
0047 //  Change y-slice used for the x-template to be closer to the actual cotalpha-cotbeta point
0048 //  (there is some weak breakdown of x-y factorization in the FPix after irradiation)
0049 //
0050 //
0051 //  V8.00 - Add method to calculate a simple 2D template
0052 //  Reorganize the interpolate method to extract header info only once per ID
0053 //  V8.01 - Improve simple template normalization
0054 //  V8.05 - Change qbin normalization to work better after irradiation
0055 //  V8.10 - Add Vavilov distribution interpolation
0056 //  V8.11 - Renormalize the x-templates for Guofan's cluster size calculation
0057 //  V8.12 - Technical fix to qavg issue.
0058 //  V8.13 - Fix qbin and fastsim interpolaters to avoid changing class variables
0059 //  V8.20 - Add methods to identify the central pixels in the x- and y-templates (to help align templates with clusters in radiation damaged detectors)
0060 //          Rename class variables from pxxxx (private xxxx) to xxxx_ to follow standard convention.
0061 //          Add compiler option to store the template entries in BOOST multiarrays of structs instead of simple c arrays
0062 //          (allows dynamic resizing template storage and has bounds checking but costs ~10% in execution time).
0063 //  V8.21 - Add new qbin method to use in cluster splitting
0064 //  V8.23 - Replace chi-min position errors with merged cluster chi2 probability info
0065 //  V8.25 - Incorporate VI's speed changes into the current version
0066 //  V8.26 - Modify the Vavilov lookups to work better with the FPix (offset y-templates)
0067 //  V8.30 - Change the splitting template generation and access to improve speed and eliminate triple index boost::multiarray
0068 //  V8.31 - Add correction factor: measured/true charge
0069 //  V8.31 - Fix version number bug in db object I/O (pushfile)
0070 //  V8.32 - Check for illegal qmin during loading
0071 //  V8.33 - Fix small type conversion warnings
0072 //  V8.40 - Incorporate V.I. optimizations
0073 //  V9.00 - Expand header to include multi and single dcol thresholds, LA biases, and (variable) Qbin definitions
0074 //  V9.01 - Protect against negative error squared
0075 //  V10.00 - Update to work with Phase 1 FPix.  Fix some code problems introduced by other maintainers.
0076 //  V10.01 - Fix initialization style as suggested by S. Krutelyov
0077 //  V10.10 - Add class variables and methods to be used to correctly calculate the probabilities of single pixel clusters
0078 //  V10.11 - Allow subdetector ID=5 for FPix R2P2 [allows better internal labeling of templates]
0079 //  V10.12 - Enforce minimum signal size in pixel charge uncertainty calculation
0080 //  V10.13 - Update the variable size [SI_PIXEL_TEMPLATE_USE_BOOST] option so that it works with VI's enhancements
0081 //  V10.20 - Add directory path selection to the ascii pushfile method
0082 //  V10.21 - Address runtime issues in pushfile() for gcc 7.X due to using tempfile as char string + misc. cleanup [Petar]
0083 //  V10.22 - Move templateStore to the heap, fix variable name in pushfile() [Petar]
0084 //  V10.24 - Add sideload() + associated gymnastics [Petar and Oz]
0085
0086 //  Created by Morris Swartz on 10/27/06.
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 //! This routine initializes the global template structures from
0126 //! an external file template_summary_zpNNNN where NNNN are four
0127 //! digits of filenum.
0128 //! \param filenum - an integer NNNN used in the filename template_summary_zpNNNN
0129 //****************************************************************
0130 bool SiPixelTemplate::pushfile(int filenum, std::vector<SiPixelTemplateStore>& pixelTemp, std::string dir) {
0131   // Add template stored in external file numbered filenum to theTemplateStore
0132
0133   // Local variables
0134   int i, j, k, l;
0135   float qavg_avg;
0136   char c;
0137   const int code_version = {17};
0138
0139   //  Create a filename for this run
0140   std::string tempfile = std::to_string(filenum);
0141
0142 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0143   // If integer filenum has less than 4 digits, prepend 0's until we have four numerical characters, e.g. "0292"
0144   int nzeros = 4 - tempfile.length();
0145   if (nzeros > 0)
0146     tempfile = std::string(nzeros, '0') + tempfile;
0147   /// Alt implementation: for (unsigned cnt=4-tempfile.length(); cnt > 0; cnt-- ){ tempfile = "0" + tempfile; }
0148
0149   tempfile = dir + "template_summary_zp" + tempfile + ".out";
0150   edm::FileInPath file(tempfile);  // Find the file in CMSSW
0151   tempfile = file.fullPath();      // Put it back with the whole path.
0152
0153 #else
0154   // This is the same as above, but more elegant.  (Elegance not allowed in CMSSW...)
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   //  Open the template file
0162   //
0163   std::ifstream in_file(tempfile);
0164   if (in_file.is_open() && in_file.good()) {
0165     // Create a local template storage entry
0166
0167     SiPixelTemplateStore theCurrentTemp;
0168
0169     // Read-in a header string first and print it
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     // next, the header information
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     // next, layout the 1-d/2-d structures needed to store template
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     // next, loop over all y-angle entries
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       // Calculate the alpha, beta, and cot(beta) for this entry
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       //        theCurrentTemp.enty[i].qbfrac[3] = 1. - theCurrentTemp.enty[i].qbfrac[0] - theCurrentTemp.enty[i].qbfrac[1] - theCurrentTemp.enty[i].qbfrac[2];
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     // next, loop over all barrel x-angle entries
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         // Calculate the alpha, beta, and cot(beta) for this entry
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         //             >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav;
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         //      theCurrentTemp.entx[k][i].qbfrac[3] = 1. - theCurrentTemp.entx[k][i].qbfrac[0] - theCurrentTemp.entx[k][i].qbfrac[1] - theCurrentTemp.entx[k][i].qbfrac[2];
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     // Add this template to the store
0702
0703     pixelTemp.push_back(theCurrentTemp);
0704
0705     postInit(pixelTemp);
0706     return true;
0707
0708   } else {
0709     // If file didn't open, report this
0710
0711     LOGERROR("SiPixelTemplate") << "Error opening File" << tempfile << ENDL;
0712     return false;
0713   }
0714
0715 }  // TempInit
0716
0717 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0718
0719 //****************************************************************
0720 //! This routine initializes the global template structures from an
0721 //! external file template_summary_zpNNNN where NNNN are four digits
0722 //! \param dbobject - db storing multiple template calibrations
0723 //****************************************************************
0724 bool SiPixelTemplate::pushfile(const SiPixelTemplateDBObject& dbobject, std::vector<SiPixelTemplateStore>& pixelTemp) {
0725   // Add template stored in external dbobject to theTemplateStore
0726
0727   // Local variables
0728   int i, j, k, l;
0729   float qavg_avg;
0730   const int code_version = {17};
0731
0732   // We must create a new object because dbobject must be a const and our stream must not be
0733   auto db(dbobject.reader());
0734
0735   // Create a local template storage entry
0736   /// SiPixelTemplateStore theCurrentTemp;    // large, don't allocate it on the stack
0737   auto tmpPtr = std::make_unique<SiPixelTemplateStore>();  // must be allocated on the heap instead
0738   auto& theCurrentTemp = *tmpPtr;
0739
0740   // Fill the template storage for each template calibration stored in the db
0741   for (int m = 0; m < db.numOfTempl(); ++m) {
0742     // Read-in a header string first and print it
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     // next, the header information
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       //std::cout<<" set fbin  "<< theCurrentTemp.head.fbin[0]<<" "<<theCurrentTemp.head.fbin[1]<<" "
0790       //         <<theCurrentTemp.head.fbin[2]<<std::endl;
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       //return false; // dk
0813     }
0814
0815 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST
0816
0817     // next, layout the 1-d/2-d structures needed to store template
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     // next, loop over all barrel y-angle entries
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       // Calculate the alpha, beta, and cot(beta) for this entry
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       //                 >> theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav;
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       //            theCurrentTemp.enty[i].qbfrac[3] = 1. - theCurrentTemp.enty[i].qbfrac[0] - theCurrentTemp.enty[i].qbfrac[1] - theCurrentTemp.enty[i].qbfrac[2];
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     // next, loop over all barrel x-angle entries
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         // Calculate the alpha, beta, and cot(beta) for this entry
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         //                     >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav;
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         //              theCurrentTemp.entx[k][i].qbfrac[3] = 1. - theCurrentTemp.entx[k][i].qbfrac[0] - theCurrentTemp.entx[k][i].qbfrac[1] - theCurrentTemp.entx[k][i].qbfrac[2];
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     // Add this template to the store
1278
1279     pixelTemp.push_back(theCurrentTemp);
1280   }
1281   postInit(pixelTemp);
1282   return true;
1283
1284 }  // TempInit
1285
1286 #endif
1287
1288 void SiPixelTemplate::postInit(std::vector<SiPixelTemplateStore>& thePixelTemp_) {
1289   /*
1290     std::cout << "SiPixelTemplate size " << thePixelTemp_.size() << std::endl;
1291     #ifndef SI_PIXEL_TEMPLATE_USE_BOOST
1292     std::cout <<"uses C arrays" << std::endl;
1293     #endif
1294
1295     int i=0;
1296     for (auto & templ : thePixelTemp_) {
1297     std::cout << i <<':' << templ.head.ID << ' ' << templ.head.NTy <<','<< templ.head.NTyx <<','<< templ.head.NTxx << std::endl;
1298     for ( auto iy=1; iy<templ.head.NTy; ++iy  ) { auto & ent = templ.enty[iy]; std::cout << ent.cotbeta <<',' << ent.cotbeta-templ.enty[iy-1].cotbeta << ' '; }
1299     std::cout << std::endl;
1300     for ( auto ix=1; ix<templ.head.NTxx; ++ix  ){ auto & ent = templ.entx[0][ix]; std::cout << ent.cotalpha <<','<< ent.cotalpha-templ.entx[0][ix-1].cotalpha << ' ';}
1301     std::cout << std::endl;
1302     ++i;
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 //! Interpolate input alpha and beta angles to produce a working template for each individual hit.
1318 //! \param id - (input) index of the template to use
1319 //! \param cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
1320 //! \param cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
1321 //! \param locBz - (input) the sign of this quantity is used to determine whether to flip cot(beta)<0 quantities from cot(beta)>0 (FPix only)
1322 //!                    for Phase 0 FPix IP-related tracks, locBz < 0 for cot(beta) > 0 and locBz > 0 for cot(beta) < 0
1323 //!                    for Phase 1 FPix IP-related tracks, see next comment
1324 //! \param locBx - (input) the sign of this quantity is used to determine whether to flip cot(alpha/beta)<0 quantities from cot(alpha/beta)>0 (FPix only)
1325 //!                    for Phase 1 FPix IP-related tracks, locBx/locBz > 0 for cot(alpha) > 0 and locBx/locBz < 0 for cot(alpha) < 0
1326 //!                    for Phase 1 FPix IP-related tracks, locBx > 0 for cot(beta) > 0 and locBx < 0 for cot(beta) < 0
1327 // ************************************************************************************************************
1328 bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx) {
1329   // Interpolate for a new set of track angles
1330
1331   //check for nan's
1332   if (!edm::isFinite(cotalpha) || !edm::isFinite(cotbeta)) {
1333     success_ = false;
1334     return success_;
1335   }
1336
1337   // Local variables
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   // If the sideloading is turned on, xtemp_ and ytemp_ are already set to what they need to be.
1344   // So we'll just exit.
1345   if (entry_sideloaded_ != nullptr) {
1346     success_ = true;
1347     return success_;
1348   }
1349
1350   // Check to see if interpolation is valid
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       // Find the index corresponding to id
1358
1359       index_id_ = -1;
1360       for (i = 0; i < (int)thePixelTemp_.size(); ++i) {
1361         //std::cout<<i<<" "<<id<<" "<<thePixelTemp_[i].head.ID<<std::endl;
1362
1363         if (id == thePixelTemp_[i].head.ID) {
1364           index_id_ = i;
1365           id_current_ = id;
1366
1367           // Copy the detector type to the private variable
1368
1369           dtype_ = thePixelTemp_[index_id_].head.Dtype;
1370
1371           // Copy the charge scaling factor to the private variable
1372
1373           qscale_ = thePixelTemp_[index_id_].head.qscale;
1374
1375           // Copy the pseudopixel signal size to the private variable
1376
1377           s50_ = thePixelTemp_[index_id_].head.s50;
1378
1379           // Copy Qbinning info to private variables
1380
1381           for (j = 0; j < 3; ++j) {
1382             fbin_[j] = thePixelTemp_[index_id_].head.fbin[j];
1383           }
1384           //std::cout<<" set fbin  "<< fbin_[0]<<" "<<fbin_[1]<<" "<<fbin_[2]<<std::endl;
1385
1386           // Pixel sizes to the private variables
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     //  qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)
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     // for some cosmics, the ususal gymnastics are incorrect
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     // next, loop over all y-angle entries
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     // Use pointers to the three angle pairs used in the interpolation
1498     //
1499     enty0_ = &thePixelTemp_[index_id_].enty[ilow];
1500     enty1_ = &thePixelTemp_[index_id_].enty[ihigh];
1501
1502     // Interpolate/store all y-related quantities (flip displacements when flip_y_)
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         // Charge loss switches sides when cot(beta) changes sign
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       //          if(flip_y_) {chi2yavgc2m_[i] = -chi2yavgc2m_[i];}
1577       chi2yminc2m_[i] = (1.f - yratio_) * enty0_->chi2yminc2m[i] + yratio_ * enty1_->chi2yminc2m[i];
1578       //          xrmsc2m[i]=(1.f - yratio_)*enty0_->xrmsc2m[i] + yratio_*enty1_->xrmsc2m[i];
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         // Since Q_fl is odd under cotbeta, it flips qutomatically, change only even terms
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     //// Single pixel cluster probabilities
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     //       If using y-spares
1604     //       for(i=0; i<10; ++i) {
1605     //          pyspare[i]=(1.f - yratio_)*enty0_->yspare[i] + yratio_*enty1_->yspare[i];
1606     //       }
1607
1608     // Interpolate and build the y-template
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         // Flip the basic y-template when the cotbeta is negative
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     // next, loop over all x-angle entries, first, find relevant y-slices
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     // Interpolate/store all x-related quantities
1676
1677     yxratio_ = yxratio;
1678     xxratio_ = xxratio;
1679
1680     // sxparmax defines the maximum charge for which the parameters xpar are defined (not rescaled by cotbeta)
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         // Charge loss switches sides when cot(alpha) changes sign
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     // Pointers to the currently interpolated point.
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     // pixmax is the maximum allowed pixel charge (used for truncation)
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       //  Try new interpolation scheme instead
1757       //
1758       //
1759       //          chi2xavgc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*entx00_->chi2xavgc2m[i] + xxratio*entx02_->chi2xavgc2m[i])
1760       //                  +yxratio*((1.f - xxratio)*entx20_->chi2xavgc2m[i] + xxratio*entx22_->chi2xavgc2m[i]);
1761
1762       //          chi2xminc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*entx00_->chi2xminc2m[i] + xxratio*entx02_->chi2xminc2m[i])
1763       //                  +yxratio*((1.f - xxratio)*entx20_->chi2xminc2m[i] + xxratio*entx22_->chi2xminc2m[i]);
1764       //
1765       //          chi2xavg_[i]=((1.f - xxratio)*thePixelTemp_[index_id_].entx[imaxx][ilow].chi2xavg[i] + xxratio*thePixelTemp_[index_id_].entx[imaxx][ihigh].chi2xavg[i]);
1766       //          if(thePixelTemp_[index_id_].entx[imaxx][imidy].chi2xavg[i] != 0.f) {chi2xavg_[i]=chi2xavg_[i]/thePixelTemp_[index_id_].entx[imaxx][imidy].chi2xavg[i]*chi2xavg[i];}
1767       //
1768       //          chi2xmin_[i]=((1.f - xxratio)*thePixelTemp_[index_id_].entx[imaxx][ilow].chi2xmin[i] + xxratio*thePixelTemp_[index_id_].entx[imaxx][ihigh].chi2xmin[i]);
1769       //          if(thePixelTemp_[index_id_].entx[imaxx][imidy].chi2xmin[i] != 0.f) {chi2xmin_[i]=chi2xmin_[i]/thePixelTemp_[index_id_].entx[imaxx][imidy].chi2xmin[i]*chi2xmin[i];}
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         // Since Q_fl is odd under cotalpha, it flips qutomatically, change only even terms
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     // Do the spares next
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     //       If using x-spares
1824     //       for(i=0; i<10; ++i) {
1825     //        pxspare[i]=(1.f - yxratio)*((1.f - xxratio)*entx00_->xspare[i] + xxratio*entx02_->xspare[i])
1826     //                +yxratio*((1.f - xxratio)*entx20_->xspare[i] + xxratio*entx22_->xspare[i]);
1827     //       }
1828
1829     // Interpolate and build the x-template
1830
1831     //  qxtempcor corrects the total charge to the actual track angles (not actually needed for the template fits, but useful for Guofan)
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         //  Take next largest x-slice for the x-template (it reduces bias in the forward direction after irradiation)
1844         //         xtemp_[i][j+2]=(1.f - xxratio)*thePixelTemp_[index_id_].entx[imaxx][ilow].xtemp[i][j] + xxratio*thePixelTemp_[index_id_].entx[imaxx][ihigh].xtemp[i][j];
1845         //         xtemp_[i][j+2]=(1.f - xxratio)*entx20_->xtemp[i][j] + xxratio*entx22_->xtemp[i][j];
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 }  // interpolate
1871
1872 // ************************************************************************************************************
1873 //! Interpolate input alpha and beta angles to produce a working template for each individual hit.
1874 //! \param id - (input) index of the template to use
1875 //! \param cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
1876 //! \param cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
1877 //! Use this for Phase 1, IP related hits
1878 // ************************************************************************************************************
1879 bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta) {
1880   // Interpolate for a new set of track angles
1881
1882   // Local variables
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 //! Interpolate input alpha and beta angles to produce a working template for each individual hit.
1896 //! \param id - (input) index of the template to use
1897 //! \param cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
1898 //! \param cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
1899 //! Use this for Phase 0, IP related hits
1900 // ************************************************************************************************************
1901 bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta, float locBz) {
1902   // Interpolate for a new set of track angles
1903
1904   // Local variables
1905   float locBx = 1.f;
1906   return SiPixelTemplate::interpolate(id, cotalpha, cotbeta, locBz, locBx);
1907 }
1908
1909 // *************************************************************************************************************************************
1910 //! Load template info for single angle point to invoke template reco for template generation
1911 //! \param      entry - (input) pointer to template entry
1912 //! \param      sizex - (input) pixel x-size
1913 //! \param      sizey - (input) pixel y-size
1914 //! \param      sizez - (input) pixel z-size
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   // Set class variables to the input parameters
1929   entry_sideloaded_ = entry;  // will bypass xtemp_[] and ytemp_[] and use those from this 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   // Set other class variables
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   // This works only for IP-related tracks
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   //  Calculate signed quantities
2045
2046   //    qxtempcor corrects the total charge to the actual track angles (not actually needed for the template fits, but useful for Guofan)
2047   // &&& What to do here?
2048   // qxtempcor=std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta0*cotbeta0+cotalpha*cotalpha));
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       // Flip the basic y-template when the cotbeta is negative
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   // Fitted errors params
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 }  // sideload
2101 #endif
2102
2103 // ************************************************************************************************************
2104 //! Return vector of y errors (squared) for an input vector of projected signals
2105 //! Add large Q scaling for use in cluster splitting.
2106 //! \param fypix - (input) index of the first real pixel in the projected cluster (doesn't include pseudopixels)
2107 //! \param lypix - (input) index of the last real pixel in the projected cluster (doesn't include pseudopixels)
2108 //! \param sythr - (input) maximum signal before de-weighting
2109 //! \param ysum - (input) 25-element vector of pixel signals
2110 //! \param ysig2 - (output) 25-element vector of y errors (squared)
2111 // ************************************************************************************************************
2112 void SiPixelTemplate::ysigma2(int fypix, int lypix, float sythr, float ysum[25], float ysig2[25])
2113
2114 {
2115   // Interpolate using quantities already stored in the private variables
2116
2117   // Local variables
2118   int i;
2119   float sigi, sigi2, sigi3, sigi4, symax, qscale, s25;
2120
2121   // Make sure that input is OK
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   // Define the maximum signal to use in the parameterization
2140
2141   symax = symax_;
2142   s25 = 0.5f * s50_;
2143   if (symax_ > syparmax_) {
2144     symax = syparmax_;
2145   }
2146
2147   // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
2148
2149   for (i = fypix - 2; i <= lypix + 2; ++i) {
2150     if (i < fypix || i > lypix) {
2151       // Nearest pseudopixels have uncertainties of 50% of threshold, next-nearest have 10% of threshold
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 }  // End ysigma2
2193
2194 // ************************************************************************************************************
2195 //! Return y error (squared) for an input signal and yindex
2196 //! Add large Q scaling for use in cluster splitting.
2197 //! \param qpixel - (input) pixel charge
2198 //! \param index - (input) y-index index of pixel
2199 //! \param ysig2 - (output) square error
2200 // ************************************************************************************************************
2201 void SiPixelTemplate::ysigma2(float qpixel, int index, float& ysig2)
2202
2203 {
2204   // Interpolate using quantities already stored in the private variables
2205
2206   // Local variables
2207   float sigi, sigi2, sigi3, sigi4, symax, qscale, err2, s25;
2208
2209   // Make sure that input is OK
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   // Define the maximum signal to use in the parameterization
2220
2221   symax = symax_;
2222   s25 = 0.5f * s50_;
2223   if (symax_ > syparmax_) {
2224     symax = syparmax_;
2225   }
2226
2227   // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
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 }  // End ysigma2
2264
2265 // ************************************************************************************************************
2266 //! Return vector of x errors (squared) for an input vector of projected signals
2267 //! Add large Q scaling for use in cluster splitting.
2268 //! \param fxpix - (input) index of the first real pixel in the projected cluster (doesn't include pseudopixels)
2269 //! \param lxpix - (input) index of the last real pixel in the projected cluster (doesn't include pseudopixels)
2270 //! \param sxthr - (input) maximum signal before de-weighting
2271 //! \param xsum - (input) 11-element vector of pixel signals
2272 //! \param xsig2 - (output) 11-element vector of x errors (squared)
2273 // ************************************************************************************************************
2274 void SiPixelTemplate::xsigma2(int fxpix, int lxpix, float sxthr, float xsum[11], float xsig2[11])
2275
2276 {
2277   // Interpolate using quantities already stored in the private variables
2278
2279   // Local variables
2280   int i;
2281   float sigi, sigi2, sigi3, sigi4, yint, sxmax, x0, qscale, s25;
2282
2283   // Make sure that input is OK
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   // Define the maximum signal to use in the parameterization
2302
2303   sxmax = sxmax_;
2304   s25 = 0.5f * s50_;
2305   if (sxmax_ > sxparmax_) {
2306     sxmax = sxparmax_;
2307   }
2308
2309   // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
2310
2311   for (i = fxpix - 2; i <= lxpix + 2; ++i) {
2312     if (i < fxpix || i > lxpix) {
2313       // Nearest pseudopixels have uncertainties of 50% of threshold, next-nearest have 10% of threshold
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       // First, do the cotbeta interpolation
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       // Next, do the cotalpha interpolation
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       // Finally, get the mid-point value of the cotalpha function
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       // Finally, rescale the yint value for cotalpha variation
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 }  // End xsigma2
2386
2387 // ************************************************************************************************************
2388 //! Return interpolated y-correction for input charge bin and qfly
2389 //! \param binq - (input) charge bin [0-3]
2390 //! \param qfly - (input) (Q_f-Q_l)/(Q_f+Q_l) for this cluster
2391 // ************************************************************************************************************
2392 float SiPixelTemplate::yflcorr(int binq, float qfly)
2393
2394 {
2395   // Interpolate using quantities already stored in the private variables
2396
2397   // Local variables
2398   float qfl, qfl2, qfl3, qfl4, qfl5, dy;
2399
2400   // Make sure that input is OK
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   // Define the maximum signal to allow before de-weighting a pixel
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   // Interpolate between the two polynomials
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 }  // End yflcorr
2442
2443 // ************************************************************************************************************
2444 //! Return interpolated x-correction for input charge bin and qflx
2445 //! \param binq - (input) charge bin [0-3]
2446 //! \param qflx - (input) (Q_f-Q_l)/(Q_f+Q_l) for this cluster
2447 // ************************************************************************************************************
2448 float SiPixelTemplate::xflcorr(int binq, float qflx)
2449
2450 {
2451   // Interpolate using quantities already stored in the private variables
2452
2453   // Local variables
2454   float qfl, qfl2, qfl3, qfl4, qfl5, dx;
2455
2456   // Make sure that input is OK
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   // Define the maximum signal to allow before de-weighting a pixel
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   // Interpolate between the two polynomials
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 }  // End xflcorr
2504
2505 // ************************************************************************************************************
2506 //! Return interpolated y-template in single call
2507 //! \param fybin - (input) index of first bin (0-40) to fill
2508 //! \param fybin - (input) index of last bin (0-40) to fill
2509 //! \param ytemplate - (output) a 41x25 output buffer
2510 // ************************************************************************************************************
2511 void SiPixelTemplate::ytemp(int fybin, int lybin, float ytemplate[41][BYSIZE])
2512
2513 {
2514   // Retrieve already interpolated quantities
2515
2516   // Local variables
2517   int i, j;
2518
2519   // Verify that input parameters are in valid range
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   // Build the y-template, the central 25 bins are here in all cases
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   //  Add   more bins if needed
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 }  // End ytemp
2580
2581 // ************************************************************************************************************
2582 //! Return interpolated y-template in single call
2583 //! \param fxbin - (input) index of first bin (0-40) to fill
2584 //! \param fxbin - (input) index of last bin (0-40) to fill
2585 //! \param xtemplate - (output) a 41x11 output buffer
2586 // ************************************************************************************************************
2587 void SiPixelTemplate::xtemp(int fxbin, int lxbin, float xtemplate[41][BXSIZE]) {
2588   // Retrieve already interpolated quantities
2589
2590   // Local variables
2591   int i, j;
2592
2593   // Verify that input parameters are in valid range
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   // Build the x-template, the central 25 bins are here in all cases
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   //  Add more bins if needed
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 }  // End xtemp
2652
2653 // ************************************************************************************************************
2654 //! Return central pixel of y template pixels above readout threshold
2655 // ************************************************************************************************************
2656 int SiPixelTemplate::cytemp()
2657
2658 {
2659   // Retrieve already interpolated quantities
2660
2661   // Local variables
2662   int j;
2663
2664   // Analyze only pixels along the central entry
2665   // First, find the maximum signal and then work out to the edges
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   //  Now search forward and backward
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 }  // End cytemp
2705
2706 // ************************************************************************************************************
2707 //! Return central pixel of x-template pixels above readout threshold
2708 // ************************************************************************************************************
2709 int SiPixelTemplate::cxtemp()
2710
2711 {
2712   // Retrieve already interpolated quantities
2713
2714   // Local variables
2715   int j;
2716
2717   // Analyze only pixels along the central entry
2718   // First, find the maximum signal and then work out to the edges
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   //  Now search forward and backward
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 }  // End cxtemp
2758
2759 // ************************************************************************************************************
2760 //! Make interpolated 3d y-template (stored as class variables)
2761 //! \param nypix - (input) number of pixels in cluster (needed to size template)
2762 //! \param nybins - (output) number of bins needed for each template projection
2763 // ************************************************************************************************************
2764 void SiPixelTemplate::ytemp3d_int(int nypix, int& nybins)
2765
2766 {
2767   // Retrieve already interpolated quantities
2768
2769   // Local variables
2770   int i, j, k;
2771   int ioff0, ioffp, ioffm;
2772
2773   // Verify that input parameters are in valid range
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   // Calculate the size of the shift in pixels needed to span the entire cluster
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   // Calculate the number of bins needed to specify each hit range
2792
2793   nybins_ = 9 + 16 * nshift;
2794
2795   // Create a 2-d working template with the correct size
2796
2797   temp2dy_.resize(boost::extents[nybins_][BYSIZE]);
2798
2799   //  The 9 central bins are copied from the interpolated private store
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   // Add the +- shifted templates
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 }  // End ytemp3d_int
2836
2837 // ************************************************************************************************************
2838 //! Return interpolated 3d y-template in single call
2839 //! \param i,j - (input) template indices
2840 //! \param ytemplate - (output) a boost 3d array containing two sets of temlate indices and the combined pixel signals
2841 // ************************************************************************************************************
2842 void SiPixelTemplate::ytemp3d(int i, int j, std::vector<float>& ytemplate)
2843
2844 {
2845   // Sum two 2-d templates to make the 3-d template
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 }  // End ytemp3d
2859
2860 // ************************************************************************************************************
2861 //! Make interpolated 3d x-template (stored as class variables)
2862 //! \param nxpix - (input) number of pixels in cluster (needed to size template)
2863 //! \param nxbins - (output) number of bins needed for each template projection
2864 // ************************************************************************************************************
2865 void SiPixelTemplate::xtemp3d_int(int nxpix, int& nxbins)
2866
2867 {
2868   // Retrieve already interpolated quantities
2869
2870   // Local variables
2871   int i, j, k;
2872   int ioff0, ioffp, ioffm;
2873
2874   // Verify that input parameters are in valid range
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   // Calculate the size of the shift in pixels needed to span the entire cluster
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   // Calculate the number of bins needed to specify each hit range
2893
2894   nxbins_ = 9 + 16 * nshift;
2895
2896   // Create a 2-d working template with the correct size
2897
2898   temp2dx_.resize(boost::extents[nxbins_][BXSIZE]);
2899
2900   //  The 9 central bins are copied from the interpolated private store
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   // Add the +- shifted templates
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 }  // End xtemp3d_int
2938
2939 // ************************************************************************************************************
2940 //! Return interpolated 3d x-template in single call
2941 //! \param i,j - (input) template indices
2942 //! \param xtemplate - (output) a boost 3d array containing two sets of temlate indices and the combined pixel signals
2943 // ************************************************************************************************************
2944 void SiPixelTemplate::xtemp3d(int i, int j, std::vector<float>& xtemplate)
2945
2946 {
2947   // Sum two 2-d templates to make the 3-d template
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 }  // End xtemp3d
2961
2962 // ************************************************************************************************************
2963 //! Interpolate beta/alpha angles to produce an expected average charge. Return int (0-4) describing the charge
2964 //! of the cluster [0: 1.5<Q/Qavg, 1: 1<Q/Qavg<1.5, 2: 0.85<Q/Qavg<1, 3: 0.95Qmin<Q<0.85Qavg, 4: Q<0.95Qmin].
2965 //! \param id - (input) index of the template to use
2966 //! \param cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
2967 //! \param cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
2968 //! \param locBz - (input) the sign of this quantity is used to determine whether to flip cot(beta)<0 quantities from cot(beta)>0 (FPix only)
2969 //!                    for Phase 0 FPix IP-related tracks, locBz < 0 for cot(beta) > 0 and locBz > 0 for cot(beta) < 0
2970 //!                    for Phase 1 FPix IP-related tracks, see next comment
2971 //! \param locBx - (input) the sign of this quantity is used to determine whether to flip cot(alpha/beta)<0 quantities from cot(alpha/beta)>0 (FPix only)
2972 //!                    for Phase 1 FPix IP-related tracks, locBx/locBz > 0 for cot(alpha) > 0 and locBx/locBz < 0 for cot(alpha) < 0
2973 //!                    for Phase 1 FPix IP-related tracks, locBx > 0 for cot(beta) > 0 and locBx < 0 for cot(beta) < 0//!
2974 //! \param qclus - (input) the cluster charge in electrons
2975 //! \param pixmax - (output) the maximum pixel charge in electrons (truncation value)
2976 //! \param sigmay - (output) the estimated y-error for CPEGeneric in microns
2977 //! \param deltay - (output) the estimated y-bias for CPEGeneric in microns
2978 //! \param sigmax - (output) the estimated x-error for CPEGeneric in microns
2979 //! \param deltax - (output) the estimated x-bias for CPEGeneric in microns
2980 //! \param sy1 - (output) the estimated y-error for 1 single-pixel clusters in microns
2981 //! \param dy1 - (output) the estimated y-bias for 1 single-pixel clusters in microns
2982 //! \param sy2 - (output) the estimated y-error for 1 double-pixel clusters in microns
2983 //! \param dy2 - (output) the estimated y-bias for 1 double-pixel clusters in microns
2984 //! \param sx1 - (output) the estimated x-error for 1 single-pixel clusters in microns
2985 //! \param dx1 - (output) the estimated x-bias for 1 single-pixel clusters in microns
2986 //! \param sx2 - (output) the estimated x-error for 1 double-pixel clusters in microns
2987 //! \param dx2 - (output) the estimated x-bias for 1 double-pixel clusters in microns
2988 //! \param lorywidth - (output) the estimated y Lorentz width
2989 //! \param lorxwidth - (output) the estimated x Lorentz width
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)  // float& lorywidth, float& lorxwidth)
3010 {
3011   // Interpolate for a new set of track angles
3012
3013   // &&& New approach: use cached pointers.
3014
3015   // Find the index corresponding to id
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   // Interpolate the absolute value of cot(beta)
3038
3039   auto acotb = std::abs(cotbeta);
3040
3041   //    qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)
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   // for some cosmics, the ususal gymnastics are incorrect
3048
3049   float cota = cotalpha;
3050   bool flip_x = false;
3051   //y flipping already taken care of by interpolate
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   // Copy the charge scaling factor to the private variable
3074
3075   auto qscale = thePixelTemp_[index].head.qscale;
3076
3077   /*
3078     lorywidth = thePixelTemp_[index].head.lorywidth;
3079     if(locBz > 0.f) {lorywidth = -lorywidth;}
3080     lorxwidth = thePixelTemp_[index].head.lorxwidth;
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   // Interpolate/store all y-related quantities (flip displacements when flip_y)
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   //  Scale the input charge to account for differences between pixelav and CMSSW simulation or data
3127
3128   auto qtotal = qscale * qclus;
3129
3130   // uncertainty and final corrections depend upon total charge bin
3131   auto fq = qtotal / qavg;
3132   int binq;
3133
3134   if (fq > fbin_[0]) {
3135     binq = 0;
3136     //std::cout<<" fq "<<fq<<" "<<qtotal<<" "<<qavg<<" "<<qclus<<" "<<qscale<<" "
3137     //       <<fbin_[0]<<" "<<fbin_[1]<<" "<<fbin_[2]<<std::endl;
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   // next, loop over all x-angle entries, first, find relevant y-slices
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   // pixmax is the maximum allowed pixel charge (used for truncation)
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   //  Take the errors and bias from the correct charge bin
3234
3235   sigmay = yrmsgen;
3236   deltay = yavggen;
3237
3238   sigmax = xrmsgen;
3239   deltax = xavggen;
3240
3241   // If the charge is too small (then flag it)
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 }  // qbin
3254
3255 // ************************************************************************************************************
3256 //! Interpolate beta/alpha angles to produce an expected average charge. Return int (0-4) describing the charge
3257 //! of the cluster [0: 1.5<Q/Qavg, 1: 1<Q/Qavg<1.5, 2: 0.85<Q/Qavg<1, 3: 0.95Qmin<Q<0.85Qavg, 4: Q<0.95Qmin].
3258 //! \param id - (input) index of the template to use
3259 //! \param cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
3260 //! \param cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
3261 //! \param locBz - (input) the sign of this quantity is used to determine whether to flip cot(beta)<0 quantities from cot(beta)>0 (FPix only)
3262 //!                    for Phase 0 FPix IP-related tracks, locBz < 0 for cot(beta) > 0 and locBz > 0 for cot(beta) < 0
3263 //!                    for Phase 1 FPix IP-related tracks, see next comment
3264 //! \param locBx - (input) the sign of this quantity is used to determine whether to flip cot(alpha/beta)<0 quantities from cot(alpha/beta)>0 (FPix only)
3265 //!                    for Phase 1 FPix IP-related tracks, locBx/locBz > 0 for cot(alpha) > 0 and locBx/locBz < 0 for cot(alpha) < 0
3266 //!                    for Phase 1 FPix IP-related tracks, locBx > 0 for cot(beta) > 0 and locBx < 0 for cot(beta) < 0//!
3267 //! \param qclus - (input) the cluster charge in electrons
3268 //! \param pixmax - (output) the maximum pixel charge in electrons (truncation value)
3269 //! \param sigmay - (output) the estimated y-error for CPEGeneric in microns
3270 //! \param deltay - (output) the estimated y-bias for CPEGeneric in microns
3271 //! \param sigmax - (output) the estimated x-error for CPEGeneric in microns
3272 //! \param deltax - (output) the estimated x-bias for CPEGeneric in microns
3273 //! \param sy1 - (output) the estimated y-error for 1 single-pixel clusters in microns
3274 //! \param dy1 - (output) the estimated y-bias for 1 single-pixel clusters in microns
3275 //! \param sy2 - (output) the estimated y-error for 1 double-pixel clusters in microns
3276 //! \param dy2 - (output) the estimated y-bias for 1 double-pixel clusters in microns
3277 //! \param sx1 - (output) the estimated x-error for 1 single-pixel clusters in microns
3278 //! \param dx1 - (output) the estimated x-bias for 1 single-pixel clusters in microns
3279 //! \param sx2 - (output) the estimated x-error for 1 double-pixel clusters in microns
3280 //! \param dx2 - (output) the estimated x-bias for 1 double-pixel clusters in microns
3281 //! \param lorywidth - (output) the estimated y Lorentz width
3282 //! \param lorxwidth - (output) the estimated x Lorentz width
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)  // float& lorywidth, float& lorxwidth)
3302 {
3303   // Interpolate for a new set of track angles
3304
3305   // Local variables
3306   float locBx = 1.f;  //  lorywidth, lorxwidth;
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);  // , lorywidth, lorxwidth);
3327
3328 }  // qbin
3329
3330 // ************************************************************************************************************
3331 //! Interpolate beta/alpha angles to produce an expected average charge. Return int (0-4) describing the charge
3332 //! of the cluster [0: 1.5<Q/Qavg, 1: 1<Q/Qavg<1.5, 2: 0.85<Q/Qavg<1, 3: 0.95Qmin<Q<0.85Qavg, 4: Q<0.95Qmin].
3333 //! \param id - (input) index of the template to use
3334 //! \param cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
3335 //! \param cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
3336 //! \param locBz - (input) the sign of this quantity is used to determine whether to flip cot(beta)<0 quantities from cot(beta)>0 (FPix only)
3337 //!                    for FPix IP-related tracks, locBz < 0 for cot(beta) > 0 and locBz > 0 for cot(beta) < 0
3338 //! \param qclus - (input) the cluster charge in electrons
3339 //! \param pixmax - (output) the maximum pixel charge in electrons (truncation value)
3340 //! \param sigmay - (output) the estimated y-error for CPEGeneric in microns
3341 //! \param deltay - (output) the estimated y-bias for CPEGeneric in microns
3342 //! \param sigmax - (output) the estimated x-error for CPEGeneric in microns
3343 //! \param deltax - (output) the estimated x-bias for CPEGeneric in microns
3344 //! \param sy1 - (output) the estimated y-error for 1 single-pixel clusters in microns
3345 //! \param dy1 - (output) the estimated y-bias for 1 single-pixel clusters in microns
3346 //! \param sy2 - (output) the estimated y-error for 1 double-pixel clusters in microns
3347 //! \param dy2 - (output) the estimated y-bias for 1 double-pixel clusters in microns
3348 //! \param sx1 - (output) the estimated x-error for 1 single-pixel clusters in microns
3349 //! \param dx1 - (output) the estimated x-bias for 1 single-pixel clusters in microns
3350 //! \param sx2 - (output) the estimated x-error for 1 double-pixel clusters in microns
3351 //! \param dx2 - (output) the estimated x-bias for 1 double-pixel clusters in microns
3352 // ************************************************************************************************************
3353 /*
3354  int SiPixelTemplate::qbin(int id, float cotalpha, float cotbeta, float locBz, float qclus, float& pixmx, float& sigmay, float& deltay, float& sigmax, float& deltax,
3355  float& sy1, float& dy1, float& sy2, float& dy2, float& sx1, float& dx1, float& sx2, float& dx2)
3356
3357  {
3358     float lorywidth, lorxwidth;
3359     return SiPixelTemplate::qbin(id, cotalpha, cotbeta, locBz, qclus, pixmx, sigmay, deltay, sigmax, deltax,
3360  sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, lorywidth, lorxwidth);
3361
3362  } // qbin
3363  */
3364
3365 // ************************************************************************************************************
3366 //! Interpolate beta/alpha angles to produce an expected average charge. Return int (0-4) describing the charge
3367 //! of the cluster [0: 1.5<Q/Qavg, 1: 1<Q/Qavg<1.5, 2: 0.85<Q/Qavg<1, 3: 0.95Qmin<Q<0.85Qavg, 4: Q<0.95Qmin].
3368 //! \param id - (input) index of the template to use
3369 //! \param cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
3370 //! \param cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
3371 //! \param qclus - (input) the cluster charge in electrons
3372 // ************************************************************************************************************
3373 int SiPixelTemplate::qbin(int id, float cotalpha, float cotbeta, float qclus) {
3374   // Interpolate for a new set of track angles
3375
3376   // Local variables
3377   float pixmx, sigmay, deltay, sigmax, deltax, sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, locBz;  //  lorywidth, lorxwidth;
3378   // Local variables
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);  // , lorywidth, lorxwidth);
3407
3408 }  // qbin
3409
3410 // ************************************************************************************************************
3411 //! Interpolate beta/alpha angles to produce an expected average charge. Return int (0-4) describing the charge
3412 //! of the cluster [0: 1.5<Q/Qavg, 1: 1<Q/Qavg<1.5, 2: 0.85<Q/Qavg<1, 3: 0.95Qmin<Q<0.85Qavg, 4: Q<0.95Qmin].
3413 //! \param id - (input) index of the template to use
3414 //! \param cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
3415 //! \param qclus - (input) the cluster charge in electrons
3416 // ************************************************************************************************************
3417 int SiPixelTemplate::qbin(int id, float cotbeta, float qclus) {
3418   // Interpolate for a new set of track angles
3419
3420   // Local variables
3421   float pixmx, sigmay, deltay, sigmax, deltax, sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, locBz,
3422       locBx;  //, lorywidth, lorxwidth;
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);  // , lorywidth, lorxwidth);
3451
3452 }  // qbin
3453
3454 // ************************************************************************************************************
3455 //! Interpolate beta/alpha angles to produce estimated errors for fastsim
3456 //! \param id - (input) index of the template to use
3457 //! \param cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
3458 //! \param cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
3459 //! \param qBin - (input) charge bin from 0-3
3460 //! \param sigmay - (output) the estimated y-error for CPETemplate in microns
3461 //! \param sigmax - (output) the estimated x-error for CPETemplate in microns
3462 //! \param sy1 - (output) the estimated y-error for 1 single-pixel clusters in microns
3463 //! \param sy2 - (output) the estimated y-error for 1 double-pixel clusters in microns
3464 //! \param sx1 - (output) the estimated x-error for 1 single-pixel clusters in microns
3465 //! \param sx2 - (output) the estimated x-error for 1 double-pixel clusters in microns
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   // Interpolate for a new set of track angles
3480
3481   // Local variables
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   //bool flip_y;
3488
3489   // Find the index corresponding to id
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   // The error information for qBin > 3 is taken to be the same as qBin=3
3518
3519   if (qBin > 3) {
3520     qBin = 3;
3521   }
3522   //
3523
3524   // Interpolate the absolute value of cot(beta)
3525
3526   acotb = std::abs(cotbeta);
3527
3528   // Copy the charge scaling factor to the private variable
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   // next, loop over all y-angle entries
3544
3545   // Interpolate/store all y-related quantities (flip displacements when flip_y)
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   // next, loop over all x-angle entries, first, find relevant y-slices
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   //  Take the errors and bias from the correct charge bin
3607
3608   sigmay = yrms;
3609
3610   sigmax = xrms;
3611
3612   return;
3613
3614 }  // temperrors
3615
3616 // ************************************************************************************************************
3617 //! Interpolate beta/alpha angles to produce estimated errors for fastsim
3618 //! \param id - (input) index of the template to use
3619 //! \param cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
3620 //! \param cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
3621 //! \param qbin_frac[4] - (output) the integrated probability for qbin=0, 0+1, 0+1+2, 0+1+2+3 (1.)
3622 //! \param ny1_frac - (output) the probability for ysize = 1 for a single-size pixel
3623 //! \param ny2_frac - (output) the probability for ysize = 1 for a double-size pixel
3624 //! \param nx1_frac - (output) the probability for xsize = 1 for a single-size pixel
3625 //! \param nx2_frac - (output) the probability for xsize = 1 for a double-size pixel
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   // Interpolate for a new set of track angles
3638
3639   // Local variables
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   //bool flip_y;
3646
3647   // Find the index corresponding to id
3648
3649   index = -1;
3650   for (i = 0; i < (int)thePixelTemp_.size(); ++i) {
3651     if (id == thePixelTemp_[i].head.ID) {
3652       index = i;
3653       //                id_current_ = id;
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   // Interpolate the absolute value of cot(beta)
3670
3671   acotb = fabs((double)cotbeta);
3672
3673   // Copy the charge scaling factor to the private variable
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   // Interpolate/store all y-related quantities (flip displacements when flip_y)
3689   ny1_frac = (1.f - yratio_) * enty0_->fracyone + yratio_ * enty1_->fracyone;
3690   ny2_frac = (1.f - yratio_) * enty0_->fracytwo + yratio_ * enty1_->fracytwo;
3691
3692   // next, loop over all x-angle entries, first, find relevant y-slices
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 }  // qbin_dist
3759
3760 // *************************************************************************************************************************************
3761 //! Make simple 2-D templates from track angles set in interpolate and hit position.
3762
3763 //! \param       xhit - (input) x-position of hit relative to the lower left corner of pixel[1][1] (to allow for the "padding" of the two-d clusters in the splitter)
3764 //! \param       yhit - (input) y-position of hit relative to the lower left corner of pixel[1][1]
3765 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel starting at cluster[1][1]
3766 //! \param    xdouble - (input) STL vector of 11 element array to flag a double-pixel starting at cluster[1][1]
3767 //! \param template2d - (output) 2d template of size matched to the cluster.  Input must be zeroed since charge is added only.
3768 // *************************************************************************************************************************************
3769
3770 bool SiPixelTemplate::simpletemplate2D(
3771     float xhit, float yhit, std::vector<bool>& ydouble, std::vector<bool>& xdouble, float template2d[BXM2][BYM2]) {
3772   // Local variables
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   //    double path;
3778   std::list<SimplePixel> list;
3779   std::list<SimplePixel>::iterator listIter, listEnd;
3780
3781   // Calculate the entry and exit points for the line charge from the track
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   // total charge length
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   // use average charge in this direction
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   //  nx is the number of x interfaces crossed by the line charge
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   // Sort the list according to the distance from the initial point
3910
3911   list.sort();
3912
3913   // Look for double pixels and adjust the list appropriately
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   // The list now contains the path lengths of the line charge in each pixel from (x0,y0).  Cacluate the lengths of the segments and the charge.
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 }  // simpletemplate2D
3995
3996 // ************************************************************************************************************
3997 //! Interpolate beta/alpha angles to produce Vavilov parameters for the charge distribution
3998 //! \param mpv   - (output) the Vavilov most probable charge (well, not really the most probable esp at large kappa)
3999 //! \param sigma - (output) the Vavilov sigma parameter
4000 //! \param kappa - (output) the Vavilov kappa parameter [0.01 (Landau-like) < kappa < 10 (Gaussian-like)
4001 // ************************************************************************************************************
4002 void SiPixelTemplate::vavilov_pars(double& mpv, double& sigma, double& kappa)
4003
4004 {
4005   // Local variables
4006   int i;
4007   int ilow, ihigh, Ny;
4008   float yratio, cotb, cotalpha0, arg;
4009
4010   // Interpolate in cotbeta only for the correct total path length (converts cotalpha, cotbeta into an effective cotbeta)
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   // Copy the charge scaling factor to the private variable
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   // next, loop over all y-angle entries
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   // Interpolate Vavilov parameters
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   // Copy to parameter list
4065
4066   //Avoid rounding difference between floats and doubles causing issues later
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 }  // vavilov_pars
4079
4080 // ************************************************************************************************************
4081 //! Interpolate beta/alpha angles to produce Vavilov parameters for the 2-cluster charge distribution
4082 //! \param mpv   - (output) the Vavilov most probable charge (well, not really the most probable esp at large kappa)
4083 //! \param sigma - (output) the Vavilov sigma parameter
4084 //! \param kappa - (output) the Vavilov kappa parameter [0.01 (Landau-like) < kappa < 10 (Gaussian-like)
4085 // ************************************************************************************************************
4086 void SiPixelTemplate::vavilov2_pars(double& mpv, double& sigma, double& kappa)
4087
4088 {
4089   // Local variables
4090   int i;
4091   int ilow, ihigh, Ny;
4092   float yratio, cotb, cotalpha0, arg;
4093
4094   // Interpolate in cotbeta only for the correct total path length (converts cotalpha, cotbeta into an effective cotbeta)
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   // Copy the charge scaling factor to the private variable
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   // next, loop over all y-angle entries
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   // Interpolate Vavilov parameters
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   // Copy to parameter list
4149
4150   mpv = (double)mpvvav2_;
4151   sigma = (double)sigmavav2_;
4152   kappa = (double)kappavav2_;
4153
4154   return;
4155
4156 }  // vavilov2_pars