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