Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-08-17 00:01:51 

0001 //

0002 //  SiPixelTemplate2D.cc  Version 2.65

0003 //

0004 //  Full 2-D templates for cluster splitting, simulated cluster reweighting, and improved cluster probability

0005 //

0006 // Created by Morris Swartz on 12/01/09.

0007 // 2009 __TheJohnsHopkinsUniversity__.

0008 //

0009 // V1.01 - fix qavg_ filling

0010 // V1.02 - Add locBz to test if FPix use is out of range

0011 // V1.03 - Fix edge checking on final template to increase template size and to properly truncate cluster

0012 // v2.00 - Major changes to accommodate 2D reconstruction

0013 // v2.10 - Change chi2 and error scaling information to work with partially reconstructed clusters

0014 // v2.20 - Add cluster charge probability information, side loading for template generation

0015 // v2.21 - Double derivative interval [improves fit convergence]

0016 // v2.25 - Resize template store to accommodate FPix Templates

0017 // v2.30 - Fix bug found by P. Shuetze that compromises sqlite file loading

0018 // v2.35 - Add directory path selection to the ascii pushfile method

0019 // v2.50 - Change template storage to dynamically allocated 2D arrays of SiPixelTemplateEntry2D structs

0020 // v2.51 - Ensure that the derivative arrays are correctly zeroed between calls

0021 // v2.52 - Improve cosmetics for increased style points from judges

0022 // v2.60 - Fix FPix multiframe lookup problem [takes +-cotalpha and +-cotbeta]

0023 // v2.61a - Code 2.60 fix correctly

0024 // v2.65 - change double pixel flags to work with new shifted reco code definition

0025 //

0026 
0027 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0028 #else
0029 #include <math.h>
0030 #endif
0031 #include <algorithm>
0032 #include <vector>
0033 //#include "boost/multi_array.hpp"

0034 #include <iostream>
0035 #include <iomanip>
0036 #include <sstream>
0037 #include <fstream>
0038 
0039 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0040 #include "CondFormats/SiPixelTransient/interface/SiPixelTemplate2D.h"
0041 #include "FWCore/Utilities/interface/FileInPath.h"
0042 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0043 #include "FWCore/Utilities/interface/isFinite.h"
0044 #define LOGERROR(x) LogError(x)
0045 #define LOGINFO(x) LogInfo(x)
0046 #define ENDL " "
0047 #include "FWCore/Utilities/interface/Exception.h"
0048 using namespace edm;
0049 #else
0050 #include "SiPixelTemplate2D.h"
0051 #define LOGERROR(x) std::cout << x << ": "
0052 #define LOGINFO(x) std::cout << x << ": "
0053 #define ENDL std::endl
0054 #endif
0055 
0056 //****************************************************************

0057 //! This routine initializes the global template structures from

0058 //! an external file template_summary_zpNNNN where NNNN are four

0059 //! digits of filenum.

0060 //! \param filenum - an integer NNNN used in the filename template_summary_zpNNNN

0061 //****************************************************************

0062 bool SiPixelTemplate2D::pushfile(int filenum, std::vector<SiPixelTemplateStore2D>& pixelTemp, std::string dir) {
0063   // Add template stored in external file numbered filenum to theTemplateStore

0064 
0065   // Local variables

0066   const int code_version = {21};
0067 
0068   //  Create a filename for this run

0069   std::string tempfile = std::to_string(filenum);
0070 
0071   //  Create different path in CMSSW than standalone

0072 
0073 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0074   // If integer filenum has less than 4 digits, prepend 0's until we have four numerical characters, e.g. "0292"

0075   int nzeros = 4 - tempfile.length();
0076   if (nzeros > 0)
0077     tempfile = std::string(nzeros, '0') + tempfile;
0078   /// Alt implementation: for (unsigned cnt=4-tempfile.length(); cnt > 0; cnt-- ){ tempfile = "0" + tempfile; }

0079 
0080   tempfile = dir + "template_summary2D_zp" + tempfile + ".out";
0081   edm::FileInPath file(tempfile);  // Find the file in CMSSW

0082   tempfile = file.fullPath();      // Put it back with the whole path.

0083 
0084 #else
0085   // This is the same as above, but more elegant.

0086   std::ostringstream tout;
0087   tout << "template_summary2D_zp" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
0088   tempfile = tout.str();
0089 
0090 #endif
0091 
0092   //  Open the template file

0093   //

0094   std::ifstream in_file(tempfile);
0095   if (in_file.is_open() && in_file.good()) {
0096     // Create a local template storage entry

0097     SiPixelTemplateStore2D theCurrentTemp;
0098 
0099     // Read-in a header string first and print it

0100     char c;
0101     for (int i = 0; (c = in_file.get()) != '\n'; ++i) {
0102       if (i < 79) {
0103         theCurrentTemp.head.title[i] = c;
0104       } else {
0105         theCurrentTemp.head.title[79] = '\0';
0106       }
0107     }
0108     LOGINFO("SiPixelTemplate2D") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
0109 
0110     // next, the header information

0111     in_file >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >>
0112         theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx >> theCurrentTemp.head.Dtype >>
0113         theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >>
0114         theCurrentTemp.head.qscale >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >>
0115         theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >>
0116         theCurrentTemp.head.zsize;
0117 
0118     if (in_file.fail()) {
0119       LOGERROR("SiPixelTemplate2D") << "Error reading file 0A, no template load" << ENDL;
0120       return false;
0121     }
0122 
0123     if (theCurrentTemp.head.templ_version > 17) {
0124       in_file >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >> theCurrentTemp.head.lorxbias >>
0125           theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >> theCurrentTemp.head.fbin[2];
0126 
0127       if (in_file.fail()) {
0128         LOGERROR("SiPixelTemplate2D") << "Error reading file 0B, no template load" << ENDL;
0129         return false;
0130       }
0131     } else {
0132       // This is for older [legacy] payloads

0133       theCurrentTemp.head.ss50 = theCurrentTemp.head.s50;
0134       theCurrentTemp.head.lorybias = theCurrentTemp.head.lorywidth / 2.f;
0135       theCurrentTemp.head.lorxbias = theCurrentTemp.head.lorxwidth / 2.f;
0136       theCurrentTemp.head.fbin[0] = 1.5f;
0137       theCurrentTemp.head.fbin[1] = 1.00f;
0138       theCurrentTemp.head.fbin[2] = 0.85f;
0139     }
0140 
0141     LOGINFO("SiPixelTemplate2D") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version "
0142                                  << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
0143                                  << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx
0144                                  << ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
0145                                  << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
0146                                  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence
0147                                  << ", Q-scaling factor " << theCurrentTemp.head.qscale << ", 1/2 multi dcol threshold "
0148                                  << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold "
0149                                  << theCurrentTemp.head.ss50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth
0150                                  << ", y Lorentz Bias " << theCurrentTemp.head.lorybias << ", x Lorentz width "
0151                                  << theCurrentTemp.head.lorxwidth << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
0152                                  << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] << ", "
0153                                  << theCurrentTemp.head.fbin[1] << ", " << theCurrentTemp.head.fbin[2]
0154                                  << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size "
0155                                  << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
0156 
0157     if (theCurrentTemp.head.templ_version < code_version) {
0158       LOGERROR("SiPixelTemplate2D") << "code expects version " << code_version << ", no template load" << ENDL;
0159       return false;
0160     }
0161 
0162     if (theCurrentTemp.head.NTy != 0) {
0163       LOGERROR("SiPixelTemplate2D")
0164           << "Trying to load 1-d template info into the 2-d template object, check your DB/global tag!" << ENDL;
0165       return false;
0166     }
0167 
0168     // next, layout the 2-d structure needed to store template

0169 
0170     theCurrentTemp.entry.resize(theCurrentTemp.head.NTyx);
0171     for (auto& item : theCurrentTemp.entry)
0172       item.resize(theCurrentTemp.head.NTxx);
0173 
0174     // Read in the file info

0175 
0176     for (int iy = 0; iy < theCurrentTemp.head.NTyx; ++iy) {
0177       for (int jx = 0; jx < theCurrentTemp.head.NTxx; ++jx) {
0178         in_file >> theCurrentTemp.entry[iy][jx].runnum >> theCurrentTemp.entry[iy][jx].costrk[0] >>
0179             theCurrentTemp.entry[iy][jx].costrk[1] >> theCurrentTemp.entry[iy][jx].costrk[2];
0180 
0181         if (in_file.fail()) {
0182           LOGERROR("SiPixelTemplate2D") << "Error reading file 1, no template load, run # "
0183                                         << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0184           return false;
0185         }
0186 
0187         // Calculate cot(alpha) and cot(beta) for this entry

0188 
0189         theCurrentTemp.entry[iy][jx].cotalpha =
0190             theCurrentTemp.entry[iy][jx].costrk[0] / theCurrentTemp.entry[iy][jx].costrk[2];
0191 
0192         theCurrentTemp.entry[iy][jx].cotbeta =
0193             theCurrentTemp.entry[iy][jx].costrk[1] / theCurrentTemp.entry[iy][jx].costrk[2];
0194 
0195         in_file >> theCurrentTemp.entry[iy][jx].qavg >> theCurrentTemp.entry[iy][jx].pixmax >>
0196             theCurrentTemp.entry[iy][jx].sxymax >> theCurrentTemp.entry[iy][jx].iymin >>
0197             theCurrentTemp.entry[iy][jx].iymax >> theCurrentTemp.entry[iy][jx].jxmin >>
0198             theCurrentTemp.entry[iy][jx].jxmax;
0199 
0200         if (in_file.fail()) {
0201           LOGERROR("SiPixelTemplate2D") << "Error reading file 2, no template load, run # "
0202                                         << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0203           return false;
0204         }
0205 
0206         for (int k = 0; k < 2; ++k) {
0207           in_file >> theCurrentTemp.entry[iy][jx].xypar[k][0] >> theCurrentTemp.entry[iy][jx].xypar[k][1] >>
0208               theCurrentTemp.entry[iy][jx].xypar[k][2] >> theCurrentTemp.entry[iy][jx].xypar[k][3] >>
0209               theCurrentTemp.entry[iy][jx].xypar[k][4];
0210 
0211           if (in_file.fail()) {
0212             LOGERROR("SiPixelTemplate2D")
0213                 << "Error reading file 3, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0214             return false;
0215           }
0216         }
0217 
0218         for (int k = 0; k < 2; ++k) {
0219           in_file >> theCurrentTemp.entry[iy][jx].lanpar[k][0] >> theCurrentTemp.entry[iy][jx].lanpar[k][1] >>
0220               theCurrentTemp.entry[iy][jx].lanpar[k][2] >> theCurrentTemp.entry[iy][jx].lanpar[k][3] >>
0221               theCurrentTemp.entry[iy][jx].lanpar[k][4];
0222 
0223           if (in_file.fail()) {
0224             LOGERROR("SiPixelTemplate2D")
0225                 << "Error reading file 4, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0226             return false;
0227           }
0228         }
0229 
0230         //  Read the 2D template entries as floats [they are formatted that way] and cast to short ints

0231 
0232         float dummy[T2YSIZE];
0233         for (int l = 0; l < 7; ++l) {
0234           for (int k = 0; k < 7; ++k) {
0235             for (int j = 0; j < T2XSIZE; ++j) {
0236               for (int i = 0; i < T2YSIZE; ++i) {
0237                 in_file >> dummy[i];
0238               }
0239               if (in_file.fail()) {
0240                 LOGERROR("SiPixelTemplate2D")
0241                     << "Error reading file 5, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0242                 return false;
0243               }
0244               for (int i = 0; i < T2YSIZE; ++i) {
0245                 theCurrentTemp.entry[iy][jx].xytemp[k][l][i][j] = (short int)dummy[i];
0246               }
0247             }
0248           }
0249         }
0250 
0251         in_file >> theCurrentTemp.entry[iy][jx].chi2ppix >> theCurrentTemp.entry[iy][jx].chi2scale >>
0252             theCurrentTemp.entry[iy][jx].offsetx[0] >> theCurrentTemp.entry[iy][jx].offsetx[1] >>
0253             theCurrentTemp.entry[iy][jx].offsetx[2] >> theCurrentTemp.entry[iy][jx].offsetx[3] >>
0254             theCurrentTemp.entry[iy][jx].offsety[0] >> theCurrentTemp.entry[iy][jx].offsety[1] >>
0255             theCurrentTemp.entry[iy][jx].offsety[2] >> theCurrentTemp.entry[iy][jx].offsety[3];
0256 
0257         if (in_file.fail()) {
0258           LOGERROR("SiPixelTemplate2D") << "Error reading file 6, no template load, run # "
0259                                         << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0260           return false;
0261         }
0262 
0263         in_file >> theCurrentTemp.entry[iy][jx].clsleny >> theCurrentTemp.entry[iy][jx].clslenx >>
0264             theCurrentTemp.entry[iy][jx].mpvvav >> theCurrentTemp.entry[iy][jx].sigmavav >>
0265             theCurrentTemp.entry[iy][jx].kappavav >> theCurrentTemp.entry[iy][jx].scalexavg >>
0266             theCurrentTemp.entry[iy][jx].scaleyavg >> theCurrentTemp.entry[iy][jx].delyavg >>
0267             theCurrentTemp.entry[iy][jx].delysig >> theCurrentTemp.entry[iy][jx].spare[0];
0268 
0269         if (in_file.fail()) {
0270           LOGERROR("SiPixelTemplate2D") << "Error reading file 7, no template load, run # "
0271                                         << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0272           return false;
0273         }
0274 
0275         in_file >> theCurrentTemp.entry[iy][jx].scalex[0] >> theCurrentTemp.entry[iy][jx].scalex[1] >>
0276             theCurrentTemp.entry[iy][jx].scalex[2] >> theCurrentTemp.entry[iy][jx].scalex[3] >>
0277             theCurrentTemp.entry[iy][jx].scaley[0] >> theCurrentTemp.entry[iy][jx].scaley[1] >>
0278             theCurrentTemp.entry[iy][jx].scaley[2] >> theCurrentTemp.entry[iy][jx].scaley[3] >>
0279             theCurrentTemp.entry[iy][jx].spare[1] >> theCurrentTemp.entry[iy][jx].spare[2];
0280 
0281         if (in_file.fail()) {
0282           LOGERROR("SiPixelTemplate2D") << "Error reading file 8, no template load, run # "
0283                                         << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0284           return false;
0285         }
0286       }
0287     }
0288 
0289     in_file.close();
0290 
0291     // Add this template to the store

0292 
0293     pixelTemp.push_back(theCurrentTemp);
0294 
0295     return true;
0296 
0297   } else {
0298     // If file didn't open, report this

0299 
0300     LOGERROR("SiPixelTemplate2D") << "Error opening File" << tempfile << ENDL;
0301     return false;
0302   }
0303 
0304 }  // TempInit

0305 
0306 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0307 
0308 //****************************************************************

0309 //! This routine initializes the global template structures from an

0310 //! external file template_summary_zpNNNN where NNNN are four digits

0311 //! \param dbobject - db storing multiple template calibrations

0312 //****************************************************************

0313 bool SiPixelTemplate2D::pushfile(const SiPixel2DTemplateDBObject& dbobject,
0314                                  std::vector<SiPixelTemplateStore2D>& pixelTemp) {
0315   // Add template stored in external dbobject to theTemplateStore

0316 
0317   const int code_version = {21};
0318 
0319   // We must create a new object because dbobject must be a const and our stream must not be

0320   SiPixel2DTemplateDBObject db = dbobject;
0321 
0322   // Create a local template storage entry

0323   SiPixelTemplateStore2D theCurrentTemp;
0324 
0325   // Fill the template storage for each template calibration stored in the db

0326   for (int m = 0; m < db.numOfTempl(); ++m) {
0327     // Read-in a header string first and print it

0328 
0329     SiPixel2DTemplateDBObject::char2float temp;
0330     for (int i = 0; i < 20; ++i) {
0331       temp.f = db.sVector()[db.index()];
0332       theCurrentTemp.head.title[4 * i] = temp.c[0];
0333       theCurrentTemp.head.title[4 * i + 1] = temp.c[1];
0334       theCurrentTemp.head.title[4 * i + 2] = temp.c[2];
0335       theCurrentTemp.head.title[4 * i + 3] = temp.c[3];
0336       db.incrementIndex(1);
0337     }
0338     theCurrentTemp.head.title[79] = '\0';
0339     LOGINFO("SiPixelTemplate2D") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
0340 
0341     // next, the header information

0342 
0343     db >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >>
0344         theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx >> theCurrentTemp.head.Dtype >>
0345         theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >>
0346         theCurrentTemp.head.qscale >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >>
0347         theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >>
0348         theCurrentTemp.head.zsize;
0349 
0350     if (db.fail()) {
0351       LOGERROR("SiPixelTemplate2D") << "Error reading file 0A, no template load" << ENDL;
0352       return false;
0353     }
0354 
0355     LOGINFO("SiPixelTemplate2D") << "Loading Pixel Template File - " << theCurrentTemp.head.title
0356                                  << " code version = " << code_version << " object version "
0357                                  << theCurrentTemp.head.templ_version << ENDL;
0358 
0359     if (theCurrentTemp.head.templ_version > 17) {
0360       db >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >> theCurrentTemp.head.lorxbias >>
0361           theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >> theCurrentTemp.head.fbin[2];
0362 
0363       if (db.fail()) {
0364         LOGERROR("SiPixelTemplate2D") << "Error reading file 0B, no template load" << ENDL;
0365         return false;
0366       }
0367     } else {
0368       // This is for older [legacy] payloads and the numbers are indeed magic [they are part of the payload for v>17]

0369       theCurrentTemp.head.ss50 = theCurrentTemp.head.s50;
0370       theCurrentTemp.head.lorybias = theCurrentTemp.head.lorywidth / 2.f;
0371       theCurrentTemp.head.lorxbias = theCurrentTemp.head.lorxwidth / 2.f;
0372       theCurrentTemp.head.fbin[0] = 1.50f;
0373       theCurrentTemp.head.fbin[1] = 1.00f;
0374       theCurrentTemp.head.fbin[2] = 0.85f;
0375     }
0376 
0377     LOGINFO("SiPixelTemplate2D") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version "
0378                                  << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
0379                                  << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx
0380                                  << ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
0381                                  << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
0382                                  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence
0383                                  << ", Q-scaling factor " << theCurrentTemp.head.qscale << ", 1/2 multi dcol threshold "
0384                                  << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold "
0385                                  << theCurrentTemp.head.ss50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth
0386                                  << ", y Lorentz Bias " << theCurrentTemp.head.lorybias << ", x Lorentz width "
0387                                  << theCurrentTemp.head.lorxwidth << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
0388                                  << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] << ", "
0389                                  << theCurrentTemp.head.fbin[1] << ", " << theCurrentTemp.head.fbin[2]
0390                                  << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size "
0391                                  << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
0392 
0393     if (theCurrentTemp.head.templ_version < code_version) {
0394       LOGINFO("SiPixelTemplate2D") << "code expects version " << code_version << " finds "
0395                                    << theCurrentTemp.head.templ_version << ", load anyway " << ENDL;
0396     }
0397 
0398     if (theCurrentTemp.head.NTy != 0) {
0399       LOGERROR("SiPixelTemplate2D")
0400           << "Trying to load 1-d template info into the 2-d template object, check your DB/global tag!" << ENDL;
0401       return false;
0402     }
0403 
0404     // next, layout the 2-d structure needed to store template

0405 
0406     theCurrentTemp.entry.resize(theCurrentTemp.head.NTyx);
0407     for (auto& item : theCurrentTemp.entry)
0408       item.resize(theCurrentTemp.head.NTxx);
0409 
0410     // Read in the file info

0411 
0412     for (int iy = 0; iy < theCurrentTemp.head.NTyx; ++iy) {
0413       for (int jx = 0; jx < theCurrentTemp.head.NTxx; ++jx) {
0414         db >> theCurrentTemp.entry[iy][jx].runnum >> theCurrentTemp.entry[iy][jx].costrk[0] >>
0415             theCurrentTemp.entry[iy][jx].costrk[1] >> theCurrentTemp.entry[iy][jx].costrk[2];
0416 
0417         if (db.fail()) {
0418           LOGERROR("SiPixelTemplate2D") << "Error reading file 1, no template load, run # "
0419                                         << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0420           return false;
0421         }
0422 
0423         // Calculate cot(alpha) and cot(beta) for this entry

0424 
0425         theCurrentTemp.entry[iy][jx].cotalpha =
0426             theCurrentTemp.entry[iy][jx].costrk[0] / theCurrentTemp.entry[iy][jx].costrk[2];
0427 
0428         theCurrentTemp.entry[iy][jx].cotbeta =
0429             theCurrentTemp.entry[iy][jx].costrk[1] / theCurrentTemp.entry[iy][jx].costrk[2];
0430 
0431         db >> theCurrentTemp.entry[iy][jx].qavg >> theCurrentTemp.entry[iy][jx].pixmax >>
0432             theCurrentTemp.entry[iy][jx].sxymax >> theCurrentTemp.entry[iy][jx].iymin >>
0433             theCurrentTemp.entry[iy][jx].iymax >> theCurrentTemp.entry[iy][jx].jxmin >>
0434             theCurrentTemp.entry[iy][jx].jxmax;
0435 
0436         if (db.fail()) {
0437           LOGERROR("SiPixelTemplate2D") << "Error reading file 2, no template load, run # "
0438                                         << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0439           return false;
0440         }
0441 
0442         for (int k = 0; k < 2; ++k) {
0443           db >> theCurrentTemp.entry[iy][jx].xypar[k][0] >> theCurrentTemp.entry[iy][jx].xypar[k][1] >>
0444               theCurrentTemp.entry[iy][jx].xypar[k][2] >> theCurrentTemp.entry[iy][jx].xypar[k][3] >>
0445               theCurrentTemp.entry[iy][jx].xypar[k][4];
0446 
0447           if (db.fail()) {
0448             LOGERROR("SiPixelTemplate2D")
0449                 << "Error reading file 3, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0450             return false;
0451           }
0452         }
0453 
0454         for (int k = 0; k < 2; ++k) {
0455           db >> theCurrentTemp.entry[iy][jx].lanpar[k][0] >> theCurrentTemp.entry[iy][jx].lanpar[k][1] >>
0456               theCurrentTemp.entry[iy][jx].lanpar[k][2] >> theCurrentTemp.entry[iy][jx].lanpar[k][3] >>
0457               theCurrentTemp.entry[iy][jx].lanpar[k][4];
0458 
0459           if (db.fail()) {
0460             LOGERROR("SiPixelTemplate2D")
0461                 << "Error reading file 4, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0462             return false;
0463           }
0464         }
0465 
0466         //  Read the 2D template entries as floats [they are formatted that way] and cast to short ints

0467 
0468         float dummy[T2YSIZE];
0469         for (int l = 0; l < 7; ++l) {
0470           for (int k = 0; k < 7; ++k) {
0471             for (int j = 0; j < T2XSIZE; ++j) {
0472               for (int i = 0; i < T2YSIZE; ++i) {
0473                 db >> dummy[i];
0474               }
0475               if (db.fail()) {
0476                 LOGERROR("SiPixelTemplate2D")
0477                     << "Error reading file 5, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0478                 return false;
0479               }
0480               for (int i = 0; i < T2YSIZE; ++i) {
0481                 theCurrentTemp.entry[iy][jx].xytemp[k][l][i][j] = (short int)dummy[i];
0482               }
0483             }
0484           }
0485         }
0486 
0487         db >> theCurrentTemp.entry[iy][jx].chi2ppix >> theCurrentTemp.entry[iy][jx].chi2scale >>
0488             theCurrentTemp.entry[iy][jx].offsetx[0] >> theCurrentTemp.entry[iy][jx].offsetx[1] >>
0489             theCurrentTemp.entry[iy][jx].offsetx[2] >> theCurrentTemp.entry[iy][jx].offsetx[3] >>
0490             theCurrentTemp.entry[iy][jx].offsety[0] >> theCurrentTemp.entry[iy][jx].offsety[1] >>
0491             theCurrentTemp.entry[iy][jx].offsety[2] >> theCurrentTemp.entry[iy][jx].offsety[3];
0492 
0493         if (db.fail()) {
0494           LOGERROR("SiPixelTemplate2D") << "Error reading file 6, no template load, run # "
0495                                         << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0496           return false;
0497         }
0498 
0499         db >> theCurrentTemp.entry[iy][jx].clsleny >> theCurrentTemp.entry[iy][jx].clslenx >>
0500             theCurrentTemp.entry[iy][jx].mpvvav >> theCurrentTemp.entry[iy][jx].sigmavav >>
0501             theCurrentTemp.entry[iy][jx].kappavav >> theCurrentTemp.entry[iy][jx].scalexavg >>
0502             theCurrentTemp.entry[iy][jx].scaleyavg >> theCurrentTemp.entry[iy][jx].delyavg >>
0503             theCurrentTemp.entry[iy][jx].delysig >> theCurrentTemp.entry[iy][jx].spare[0];
0504 
0505         if (db.fail()) {
0506           LOGERROR("SiPixelTemplate2D") << "Error reading file 7, no template load, run # "
0507                                         << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0508           return false;
0509         }
0510 
0511         db >> theCurrentTemp.entry[iy][jx].scalex[0] >> theCurrentTemp.entry[iy][jx].scalex[1] >>
0512             theCurrentTemp.entry[iy][jx].scalex[2] >> theCurrentTemp.entry[iy][jx].scalex[3] >>
0513             theCurrentTemp.entry[iy][jx].scaley[0] >> theCurrentTemp.entry[iy][jx].scaley[1] >>
0514             theCurrentTemp.entry[iy][jx].scaley[2] >> theCurrentTemp.entry[iy][jx].scaley[3] >>
0515             theCurrentTemp.entry[iy][jx].spare[1] >> theCurrentTemp.entry[iy][jx].spare[2];
0516 
0517         if (db.fail()) {
0518           LOGERROR("SiPixelTemplate2D") << "Error reading file 8, no template load, run # "
0519                                         << theCurrentTemp.entry[iy][jx].runnum << ENDL;
0520           return false;
0521         }
0522       }
0523     }
0524 
0525     // Add this template to the store

0526 
0527     pixelTemp.push_back(theCurrentTemp);
0528   }
0529 
0530   return true;
0531 
0532 }  // TempInit

0533 
0534 #endif
0535 
0536 bool SiPixelTemplate2D::getid(int id) {
0537   if (id != id_current_) {
0538     // Find the index corresponding to id

0539 
0540     index_id_ = -1;
0541     for (int i = 0; i < (int)thePixelTemp_.size(); ++i) {
0542       if (id == thePixelTemp_[i].head.ID) {
0543         index_id_ = i;
0544         id_current_ = id;
0545 
0546         // Copy the detector type to the private variable

0547 
0548         Dtype_ = thePixelTemp_[index_id_].head.Dtype;
0549 
0550         // Copy the charge scaling factor to the private variable

0551 
0552         qscale_ = thePixelTemp_[index_id_].head.qscale;
0553 
0554         // Copy the pseudopixel signal size to the private variable

0555 
0556         s50_ = thePixelTemp_[index_id_].head.s50;
0557 
0558         // Copy Qbinning info to private variables

0559 
0560         for (int j = 0; j < 3; ++j) {
0561           fbin_[j] = thePixelTemp_[index_id_].head.fbin[j];
0562         }
0563 
0564         // Copy the Lorentz widths to private variables

0565 
0566         lorywidth_ = thePixelTemp_[index_id_].head.lorywidth;
0567         lorxwidth_ = thePixelTemp_[index_id_].head.lorxwidth;
0568 
0569         // Copy the pixel sizes private variables

0570 
0571         xsize_ = thePixelTemp_[index_id_].head.xsize;
0572         ysize_ = thePixelTemp_[index_id_].head.ysize;
0573         zsize_ = thePixelTemp_[index_id_].head.zsize;
0574 
0575         // Determine the size of this template

0576 
0577         Nyx_ = thePixelTemp_[index_id_].head.NTyx;
0578         Nxx_ = thePixelTemp_[index_id_].head.NTxx;
0579 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0580         if (Nyx_ < 2 || Nxx_ < 2) {
0581           throw cms::Exception("DataCorrupt") << "template ID = " << id_current_
0582                                               << "has too few entries: Nyx/Nxx = " << Nyx_ << "/" << Nxx_ << std::endl;
0583         }
0584 #else
0585         assert(Nyx_ > 1 && Nxx_ > 1);
0586 #endif
0587         int imidx = Nxx_ / 2;
0588 
0589         cotalpha0_ = thePixelTemp_[index_id_].entry[0][0].cotalpha;
0590         cotalpha1_ = thePixelTemp_[index_id_].entry[0][Nxx_ - 1].cotalpha;
0591         deltacota_ = (cotalpha1_ - cotalpha0_) / (float)(Nxx_ - 1);
0592 
0593         cotbeta0_ = thePixelTemp_[index_id_].entry[0][imidx].cotbeta;
0594         cotbeta1_ = thePixelTemp_[index_id_].entry[Nyx_ - 1][imidx].cotbeta;
0595         deltacotb_ = (cotbeta1_ - cotbeta0_) / (float)(Nyx_ - 1);
0596 
0597         break;
0598       }
0599     }
0600   }
0601 
0602 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0603   if (index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
0604     throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::interpolate can't find needed template ID = " << id
0605                                         << ", Are you using the correct global tag?" << std::endl;
0606   }
0607 #else
0608   assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
0609 #endif
0610   return true;
0611 }
0612 
0613 // *************************************************************************************************************************************

0614 //! Interpolate stored 2-D information for input angles

0615 //! \param         id - (input) the id of the template

0616 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)

0617 //! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)

0618 //! \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)

0619 //!                    for Phase 0 FPix IP-related tracks, locBz < 0 for cot(beta) > 0 and locBz > 0 for cot(beta) < 0

0620 //!                    for Phase 1 FPix IP-related tracks, see next comment

0621 //! \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)

0622 //!                    for Phase 1 FPix IP-related tracks, locBx/locBz > 0 for cot(alpha) > 0 and locBx/locBz < 0 for cot(alpha) < 0

0623 //!                    for Phase 1 FPix IP-related tracks, locBx > 0 for cot(beta) > 0 and locBx < 0 for cot(beta) < 0

0624 // *************************************************************************************************************************************

0625 
0626 bool SiPixelTemplate2D::interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx) {
0627   // Interpolate for a new set of track angles

0628 
0629   // Local variables

0630 
0631   float acotb, dcota, dcotb;
0632 
0633   // Check to see if interpolation is valid

0634 
0635   if (id != id_current_ || cotalpha != cota_current_ || cotbeta != cotb_current_) {
0636     cota_current_ = cotalpha;
0637     cotb_current_ = cotbeta;
0638     // Try to find the correct template.  Fill the class variable index_id_ .

0639     success_ = getid(id);
0640   }
0641 
0642 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0643   if (index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
0644     throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::interpolate can't find needed template ID = " << id
0645                                         << ", Are you using the correct global tag?" << std::endl;
0646   }
0647 #else
0648   assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
0649 #endif
0650 
0651   // Check angle limits and et up interpolation parameters

0652 
0653   float cota = cotalpha;
0654   flip_x_ = false;
0655   flip_y_ = false;
0656   switch (Dtype_) {
0657     case 0:
0658       if (cotbeta < 0.f) {
0659         flip_y_ = true;
0660       }
0661       break;
0662     case 1:
0663       if (locBz > 0.f) {
0664         flip_y_ = true;
0665       }
0666       break;
0667     case 2:
0668     case 3:
0669     case 4:
0670     case 5:
0671       if (locBx * locBz < 0.f) {
0672         cota = std::abs(cotalpha);
0673         flip_x_ = true;
0674       }
0675       if (locBx < 0.f) {
0676         flip_y_ = true;
0677       }
0678       break;
0679     default:
0680 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0681       throw cms::Exception("DataCorrupt")
0682           << "SiPixelTemplate2D::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
0683 
0684       //check for nan's

0685       if (!edm::isFinite(cotalpha) || !edm::isFinite(cotbeta)) {
0686         success_ = false;
0687         return success_;
0688       }
0689 #else
0690       std::cout << "SiPixelTemplate:2D:illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
0691 #endif
0692   }
0693 
0694   if (cota < cotalpha0_) {
0695     success_ = false;
0696     jx0_ = 0;
0697     jx1_ = 1;
0698     adcota_ = 0.f;
0699   } else if (cota > cotalpha1_) {
0700     success_ = false;
0701     jx0_ = Nxx_ - 1;
0702     jx1_ = jx0_ - 1;
0703     adcota_ = 0.f;
0704   } else {
0705     jx0_ = (int)((cota - cotalpha0_) / deltacota_ + 0.5f);
0706     dcota = (cota - (cotalpha0_ + jx0_ * deltacota_)) / deltacota_;
0707     adcota_ = fabs(dcota);
0708     if (dcota > 0.f) {
0709       jx1_ = jx0_ + 1;
0710       if (jx1_ > Nxx_ - 1)
0711         jx1_ = jx0_ - 1;
0712     } else {
0713       jx1_ = jx0_ - 1;
0714       if (jx1_ < 0)
0715         jx1_ = jx0_ + 1;
0716     }
0717   }
0718 
0719   // Interpolate the absolute value of cot(beta)

0720 
0721   acotb = std::abs(cotbeta);
0722 
0723   if (acotb < cotbeta0_) {
0724     success_ = false;
0725     iy0_ = 0;
0726     iy1_ = 1;
0727     adcotb_ = 0.f;
0728   } else if (acotb > cotbeta1_) {
0729     success_ = false;
0730     iy0_ = Nyx_ - 1;
0731     iy1_ = iy0_ - 1;
0732     adcotb_ = 0.f;
0733   } else {
0734     iy0_ = (int)((acotb - cotbeta0_) / deltacotb_ + 0.5f);
0735     dcotb = (acotb - (cotbeta0_ + iy0_ * deltacotb_)) / deltacotb_;
0736     adcotb_ = fabs(dcotb);
0737     if (dcotb > 0.f) {
0738       iy1_ = iy0_ + 1;
0739       if (iy1_ > Nyx_ - 1)
0740         iy1_ = iy0_ - 1;
0741     } else {
0742       iy1_ = iy0_ - 1;
0743       if (iy1_ < 0)
0744         iy1_ = iy0_ + 1;
0745     }
0746   }
0747 
0748   //  Calculate signed quantities

0749 
0750   lorydrift_ = lorywidth_ / 2.;
0751   if (flip_y_)
0752     lorydrift_ = -lorydrift_;
0753   lorxdrift_ = lorxwidth_ / 2.;
0754   if (flip_x_)
0755     lorxdrift_ = -lorxdrift_;
0756 
0757   // Use pointers to the three angle pairs used in the interpolation

0758 
0759   entry00_ = &thePixelTemp_[index_id_].entry[iy0_][jx0_];
0760   entry10_ = &thePixelTemp_[index_id_].entry[iy1_][jx0_];
0761   entry01_ = &thePixelTemp_[index_id_].entry[iy0_][jx1_];
0762 
0763   // Interpolate things in cot(alpha)-cot(beta)

0764 
0765   qavg_ = entry00_->qavg + adcota_ * (entry01_->qavg - entry00_->qavg) + adcotb_ * (entry10_->qavg - entry00_->qavg);
0766 
0767   pixmax_ = entry00_->pixmax + adcota_ * (entry01_->pixmax - entry00_->pixmax) +
0768             adcotb_ * (entry10_->pixmax - entry00_->pixmax);
0769 
0770   sxymax_ = entry00_->sxymax + adcota_ * (entry01_->sxymax - entry00_->sxymax) +
0771             adcotb_ * (entry10_->sxymax - entry00_->sxymax);
0772 
0773   chi2avgone_ = entry00_->chi2avgone + adcota_ * (entry01_->chi2avgone - entry00_->chi2avgone) +
0774                 adcotb_ * (entry10_->chi2avgone - entry00_->chi2avgone);
0775 
0776   chi2minone_ = entry00_->chi2minone + adcota_ * (entry01_->chi2minone - entry00_->chi2minone) +
0777                 adcotb_ * (entry10_->chi2minone - entry00_->chi2minone);
0778 
0779   clsleny_ = entry00_->clsleny + adcota_ * (entry01_->clsleny - entry00_->clsleny) +
0780              adcotb_ * (entry10_->clsleny - entry00_->clsleny);
0781 
0782   clslenx_ = entry00_->clslenx + adcota_ * (entry01_->clslenx - entry00_->clslenx) +
0783              adcotb_ * (entry10_->clslenx - entry00_->clslenx);
0784 
0785   chi2ppix_ = entry00_->chi2ppix + adcota_ * (entry01_->chi2ppix - entry00_->chi2ppix) +
0786               adcotb_ * (entry10_->chi2ppix - entry00_->chi2ppix);
0787 
0788   chi2scale_ = entry00_->chi2scale + adcota_ * (entry01_->chi2scale - entry00_->chi2scale) +
0789                adcotb_ * (entry10_->chi2scale - entry00_->chi2scale);
0790 
0791   scaleyavg_ = entry00_->scaleyavg + adcota_ * (entry01_->scaleyavg - entry00_->scaleyavg) +
0792                adcotb_ * (entry10_->scaleyavg - entry00_->scaleyavg);
0793 
0794   scalexavg_ = entry00_->scalexavg + adcota_ * (entry01_->scalexavg - entry00_->scalexavg) +
0795                adcotb_ * (entry10_->scalexavg - entry00_->scalexavg);
0796 
0797   delyavg_ = entry00_->delyavg + adcota_ * (entry01_->delyavg - entry00_->delyavg) +
0798              adcotb_ * (entry10_->delyavg - entry00_->delyavg);
0799 
0800   delysig_ = entry00_->delysig + adcota_ * (entry01_->delysig - entry00_->delysig) +
0801              adcotb_ * (entry10_->delysig - entry00_->delysig);
0802 
0803   mpvvav_ = entry00_->mpvvav + adcota_ * (entry01_->mpvvav - entry00_->mpvvav) +
0804             adcotb_ * (entry10_->mpvvav - entry00_->mpvvav);
0805 
0806   sigmavav_ = entry00_->sigmavav + adcota_ * (entry01_->sigmavav - entry00_->sigmavav) +
0807               adcotb_ * (entry10_->sigmavav - entry00_->sigmavav);
0808 
0809   kappavav_ = entry00_->kappavav + adcota_ * (entry01_->kappavav - entry00_->kappavav) +
0810               adcotb_ * (entry10_->kappavav - entry00_->kappavav);
0811 
0812   for (int i = 0; i < 4; ++i) {
0813     scalex_[i] = entry00_->scalex[i] + adcota_ * (entry01_->scalex[i] - entry00_->scalex[i]) +
0814                  adcotb_ * (entry10_->scalex[i] - entry00_->scalex[i]);
0815 
0816     scaley_[i] = entry00_->scaley[i] + adcota_ * (entry01_->scaley[i] - entry00_->scaley[i]) +
0817                  adcotb_ * (entry10_->scaley[i] - entry00_->scaley[i]);
0818 
0819     offsetx_[i] = entry00_->offsetx[i] + adcota_ * (entry01_->offsetx[i] - entry00_->offsetx[i]) +
0820                   adcotb_ * (entry10_->offsetx[i] - entry00_->offsetx[i]);
0821     if (flip_x_)
0822       offsetx_[i] = -offsetx_[i];
0823 
0824     offsety_[i] = entry00_->offsety[i] + adcota_ * (entry01_->offsety[i] - entry00_->offsety[i]) +
0825                   adcotb_ * (entry10_->offsety[i] - entry00_->offsety[i]);
0826     if (flip_y_)
0827       offsety_[i] = -offsety_[i];
0828   }
0829 
0830   for (int i = 0; i < 2; ++i) {
0831     for (int j = 0; j < 5; ++j) {
0832       // Charge loss switches sides when cot(beta) changes sign

0833       if (flip_y_) {
0834         xypary0x0_[1 - i][j] = (float)entry00_->xypar[i][j];
0835         xypary1x0_[1 - i][j] = (float)entry10_->xypar[i][j];
0836         xypary0x1_[1 - i][j] = (float)entry01_->xypar[i][j];
0837         lanpar_[1 - i][j] = entry00_->lanpar[i][j] + adcota_ * (entry01_->lanpar[i][j] - entry00_->lanpar[i][j]) +
0838                             adcotb_ * (entry10_->lanpar[i][j] - entry00_->lanpar[i][j]);
0839       } else {
0840         xypary0x0_[i][j] = (float)entry00_->xypar[i][j];
0841         xypary1x0_[i][j] = (float)entry10_->xypar[i][j];
0842         xypary0x1_[i][j] = (float)entry01_->xypar[i][j];
0843         lanpar_[i][j] = entry00_->lanpar[i][j] + adcota_ * (entry01_->lanpar[i][j] - entry00_->lanpar[i][j]) +
0844                         adcotb_ * (entry10_->lanpar[i][j] - entry00_->lanpar[i][j]);
0845       }
0846     }
0847   }
0848 
0849   return success_;
0850 }  // interpolate

0851 
0852 // *************************************************************************************************************************************

0853 //! Load template info for single angle point to invoke template reco for template generation

0854 //! \param      entry - (input) pointer to template entry

0855 //! \param      sizex - (input) pixel x-size

0856 //! \param      sizey - (input) pixel y-size

0857 //! \param      sizez - (input) pixel z-size

0858 // *************************************************************************************************************************************

0859 
0860 #ifdef SI_PIXEL_TEMPLATE_STANDALONE
0861 void SiPixelTemplate2D::sideload(SiPixelTemplateEntry2D* entry,
0862                                  int iDtype,
0863                                  float locBx,
0864                                  float locBz,
0865                                  float lorwdy,
0866                                  float lorwdx,
0867                                  float q50,
0868                                  float fbin[3],
0869                                  float xsize,
0870                                  float ysize,
0871                                  float zsize) {
0872   // Set class variables to the input parameters

0873 
0874   entry00_ = entry;
0875   entry01_ = entry;
0876   entry10_ = entry;
0877   Dtype_ = iDtype;
0878   lorywidth_ = lorwdy;
0879   lorxwidth_ = lorwdx;
0880   xsize_ = xsize;
0881   ysize_ = ysize;
0882   zsize_ = zsize;
0883   s50_ = q50;
0884   qscale_ = 1.f;
0885   for (int i = 0; i < 3; ++i) {
0886     fbin_[i] = fbin[i];
0887   }
0888 
0889   // Set other class variables

0890 
0891   adcota_ = 0.f;
0892   adcotb_ = 0.f;
0893 
0894   // Interpolate things in cot(alpha)-cot(beta)

0895 
0896   qavg_ = entry00_->qavg;
0897 
0898   pixmax_ = entry00_->pixmax;
0899 
0900   sxymax_ = entry00_->sxymax;
0901 
0902   clsleny_ = entry00_->clsleny;
0903 
0904   clslenx_ = entry00_->clslenx;
0905 
0906   scaleyavg_ = 1.f;
0907 
0908   scalexavg_ = 1.f;
0909 
0910   delyavg_ = 0.f;
0911 
0912   delysig_ = 0.f;
0913 
0914   for (int i = 0; i < 4; ++i) {
0915     scalex_[i] = 1.f;
0916     scaley_[i] = 1.f;
0917     offsetx_[i] = 0.f;
0918     offsety_[i] = 0.f;
0919   }
0920 
0921   // This works only for IP-related tracks

0922 
0923   flip_x_ = false;
0924   flip_y_ = false;
0925   float cotbeta = entry00_->cotbeta;
0926   switch (Dtype_) {
0927     case 0:
0928       if (cotbeta < 0.f) {
0929         flip_y_ = true;
0930       }
0931       break;
0932     case 1:
0933       if (locBz > 0.f) {
0934         flip_y_ = true;
0935       }
0936       break;
0937     case 2:
0938     case 3:
0939     case 4:
0940     case 5:
0941       if (locBx * locBz < 0.f) {
0942         flip_x_ = true;
0943       }
0944       if (locBx < 0.f) {
0945         flip_y_ = true;
0946       }
0947       break;
0948     default:
0949       std::cout << "SiPixelTemplate:2D:illegal subdetector ID = " << iDtype << std::endl;
0950   }
0951 
0952   //  Calculate signed quantities

0953 
0954   lorydrift_ = lorywidth_ / 2.;
0955   if (flip_y_)
0956     lorydrift_ = -lorydrift_;
0957   lorxdrift_ = lorxwidth_ / 2.;
0958   if (flip_x_)
0959     lorxdrift_ = -lorxdrift_;
0960 
0961   for (int i = 0; i < 2; ++i) {
0962     for (int j = 0; j < 5; ++j) {
0963       // Charge loss switches sides when cot(beta) changes sign

0964       if (flip_y_) {
0965         xypary0x0_[1 - i][j] = (float)entry00_->xypar[i][j];
0966         xypary1x0_[1 - i][j] = (float)entry00_->xypar[i][j];
0967         xypary0x1_[1 - i][j] = (float)entry00_->xypar[i][j];
0968         lanpar_[1 - i][j] = entry00_->lanpar[i][j];
0969       } else {
0970         xypary0x0_[i][j] = (float)entry00_->xypar[i][j];
0971         xypary1x0_[i][j] = (float)entry00_->xypar[i][j];
0972         xypary0x1_[i][j] = (float)entry00_->xypar[i][j];
0973         lanpar_[i][j] = entry00_->lanpar[i][j];
0974       }
0975     }
0976   }
0977   return;
0978 }
0979 #endif
0980 
0981 // *************************************************************************************************************************************

0982 //! \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)

0983 //! \param       yhit - (input) y-position of hit relative to the lower left corner of pixel[1][1]

0984 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel starting at cluster[1][1]

0985 //! \param    xdouble - (input) STL vector of 11 element array to flag a double-pixel starting at cluster[1][1]

0986 //! \param template2d - (output) 2d template of size matched to the cluster.  Input must be zeroed since charge is added only.

0987 // *************************************************************************************************************************************

0988 
0989 bool SiPixelTemplate2D::xytemp(float xhit,
0990                                float yhit,
0991                                bool ydouble[BYM2],
0992                                bool xdouble[BXM2],
0993                                float template2d[BXM2][BYM2],
0994                                bool derivatives,
0995                                float dpdx2d[2][BXM2][BYM2],
0996                                float& QTemplate) {
0997   // Interpolate for a new set of track angles

0998 
0999   // Local variables

1000   int pixx, pixy, k0, k1, l0, l1, deltax, deltay, iflipy, jflipx, imin, imax, jmin, jmax;
1001   int m, n;
1002   float dx, dy, ddx, ddy, adx, ady;
1003   //   const float deltaxy[2] = {8.33f, 12.5f};

1004   const float deltaxy[2] = {16.67f, 25.0f};
1005 
1006   // Check to see if interpolation is valid

1007 
1008   // next, determine the indices of the closest point in k (y-displacement), l (x-displacement)

1009   // pixy and pixx are the indices of the struck pixel in the (Ty,Tx) system

1010   // k0,k1 are the k-indices of the closest and next closest point

1011   // l0,l1 are the l-indices of the closest and next closest point

1012 
1013   pixy = (int)floorf(yhit / ysize_);
1014   dy = yhit - (pixy + 0.5f) * ysize_;
1015   if (flip_y_) {
1016     dy = -dy;
1017   }
1018   k0 = (int)(dy / ysize_ * 6.f + 3.5f);
1019   if (k0 < 0)
1020     k0 = 0;
1021   if (k0 > 6)
1022     k0 = 6;
1023   ddy = 6.f * dy / ysize_ - (k0 - 3);
1024   ady = fabs(ddy);
1025   if (ddy > 0.f) {
1026     k1 = k0 + 1;
1027     if (k1 > 6)
1028       k1 = k0 - 1;
1029   } else {
1030     k1 = k0 - 1;
1031     if (k1 < 0)
1032       k1 = k0 + 1;
1033   }
1034   pixx = (int)floorf(xhit / xsize_);
1035   dx = xhit - (pixx + 0.5f) * xsize_;
1036   if (flip_x_) {
1037     dx = -dx;
1038   }
1039   l0 = (int)(dx / xsize_ * 6.f + 3.5f);
1040   if (l0 < 0)
1041     l0 = 0;
1042   if (l0 > 6)
1043     l0 = 6;
1044   ddx = 6.f * dx / xsize_ - (l0 - 3);
1045   adx = fabs(ddx);
1046   if (ddx > 0.f) {
1047     l1 = l0 + 1;
1048     if (l1 > 6)
1049       l1 = l0 - 1;
1050   } else {
1051     l1 = l0 - 1;
1052     if (l1 < 0)
1053       l1 = l0 + 1;
1054   }
1055 
1056   // OK, lets do the template interpolation.

1057 
1058   // First find the limits of the indices for non-zero pixels

1059 
1060   imin = std::min(entry00_->iymin, entry10_->iymin);
1061   imin_ = std::min(imin, entry01_->iymin);
1062 
1063   jmin = std::min(entry00_->jxmin, entry10_->jxmin);
1064   jmin_ = std::min(jmin, entry01_->jxmin);
1065 
1066   imax = std::max(entry00_->iymax, entry10_->iymax);
1067   imax_ = std::max(imax, entry01_->iymax);
1068 
1069   jmax = std::max(entry00_->jxmax, entry10_->jxmax);
1070   jmax_ = std::max(jmax, entry01_->jxmax);
1071 
1072   // Calculate the x and y offsets to make the new template

1073 
1074   // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system

1075 
1076   ++pixy;
1077   ++pixx;
1078 
1079   // In the template store, the struck pixel is always (THy,THx)

1080 
1081   deltax = pixx - T2HX;
1082   deltay = pixy - T2HY;
1083 
1084   //  First zero the local 2-d template

1085 
1086   for (int j = 0; j < BXM2; ++j) {
1087     for (int i = 0; i < BYM2; ++i) {
1088       xytemp_[j][i] = 0.f;
1089     }
1090   }
1091 
1092   // Loop over the non-zero part of the template index space and interpolate

1093 
1094   for (int j = jmin_; j <= jmax_; ++j) {
1095     // Flip indices as needed

1096     if (flip_x_) {
1097       jflipx = T2XSIZE - 1 - j;
1098       m = deltax + jflipx;
1099     } else {
1100       m = deltax + j;
1101     }
1102     for (int i = imin_; i <= imax_; ++i) {
1103       if (flip_y_) {
1104         iflipy = T2YSIZE - 1 - i;
1105         n = deltay + iflipy;
1106       } else {
1107         n = deltay + i;
1108       }
1109       if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1110         xytemp_[m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1111                         adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1112                         ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1113                         adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1114                         adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1115       }
1116     }
1117   }
1118 
1119   //combine rows and columns to simulate double pixels

1120 
1121   for (int n = 1; n < BYM3; ++n) {
1122     if (ydouble[n]) {
1123       //  Combine the y-columns

1124       for (int m = 1; m < BXM3; ++m) {
1125         xytemp_[m][n] += xytemp_[m][n + 1];
1126       }
1127       //  Now shift the remaining pixels over by one column

1128       for (int i = n + 1; i < BYM3; ++i) {
1129         for (int m = 1; m < BXM3; ++m) {
1130           xytemp_[m][i] = xytemp_[m][i + 1];
1131         }
1132       }
1133     }
1134   }
1135 
1136   //combine rows and columns to simulate double pixels

1137 
1138   for (int m = 1; m < BXM3; ++m) {
1139     if (xdouble[m]) {
1140       //  Combine the x-rows

1141       for (int n = 1; n < BYM3; ++n) {
1142         xytemp_[m][n] += xytemp_[m + 1][n];
1143       }
1144       //  Now shift the remaining pixels over by one row

1145       for (int j = m + 1; j < BXM3; ++j) {
1146         for (n = 1; n < BYM3; ++n) {
1147           xytemp_[j][n] = xytemp_[j + 1][n];
1148         }
1149       }
1150     }
1151   }
1152 
1153   //  Finally, loop through and increment the external template

1154 
1155   float qtemptot = 0.f;
1156 
1157   for (int n = 1; n < BYM3; ++n) {
1158     for (int m = 1; m < BXM3; ++m) {
1159       if (xytemp_[m][n] != 0.f) {
1160         template2d[m][n] += xytemp_[m][n];
1161         qtemptot += xytemp_[m][n];
1162       }
1163     }
1164   }
1165 
1166   QTemplate = qtemptot;
1167 
1168   if (derivatives) {
1169     float dxytempdx[2][BXM2][BYM2], dxytempdy[2][BXM2][BYM2];
1170 
1171     for (int k = 0; k < 2; ++k) {
1172       for (int i = 0; i < BXM2; ++i) {
1173         for (int j = 0; j < BYM2; ++j) {
1174           dxytempdx[k][i][j] = 0.f;
1175           dxytempdy[k][i][j] = 0.f;
1176           dpdx2d[k][i][j] = 0.f;
1177         }
1178       }
1179     }
1180 
1181     // First do shifted +x template

1182 
1183     pixx = (int)floorf((xhit + deltaxy[0]) / xsize_);
1184     dx = (xhit + deltaxy[0]) - (pixx + 0.5f) * xsize_;
1185     if (flip_x_) {
1186       dx = -dx;
1187     }
1188     l0 = (int)(dx / xsize_ * 6.f + 3.5f);
1189     if (l0 < 0)
1190       l0 = 0;
1191     if (l0 > 6)
1192       l0 = 6;
1193     ddx = 6.f * dx / xsize_ - (l0 - 3);
1194     adx = fabs(ddx);
1195     if (ddx > 0.f) {
1196       l1 = l0 + 1;
1197       if (l1 > 6)
1198         l1 = l0 - 1;
1199     } else {
1200       l1 = l0 - 1;
1201       if (l1 < 0)
1202         l1 = l0 + 1;
1203     }
1204 
1205     // OK, lets do the template interpolation.

1206 
1207     // Calculate the x and y offsets to make the new template

1208 
1209     // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system

1210 
1211     ++pixx;
1212 
1213     // In the template store, the struck pixel is always (THy,THx)

1214 
1215     deltax = pixx - T2HX;
1216 
1217     // Loop over the non-zero part of the template index space and interpolate

1218 
1219     for (int j = jmin_; j <= jmax_; ++j) {
1220       // Flip indices as needed

1221       if (flip_x_) {
1222         jflipx = T2XSIZE - 1 - j;
1223         m = deltax + jflipx;
1224       } else {
1225         m = deltax + j;
1226       }
1227       for (int i = imin_; i <= imax_; ++i) {
1228         if (flip_y_) {
1229           iflipy = T2YSIZE - 1 - i;
1230           n = deltay + iflipy;
1231         } else {
1232           n = deltay + i;
1233         }
1234         if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1235           dxytempdx[1][m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1236                                adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1237                                ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1238                                adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1239                                adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1240         }
1241       }
1242     }
1243 
1244     //combine rows and columns to simulate double pixels

1245 
1246     for (int n = 1; n < BYM3; ++n) {
1247       if (ydouble[n]) {
1248         //  Combine the y-columns

1249         for (int m = 1; m < BXM3; ++m) {
1250           dxytempdx[1][m][n] += dxytempdx[1][m][n + 1];
1251         }
1252         //  Now shift the remaining pixels over by one column

1253         for (int i = n + 1; i < BYM3; ++i) {
1254           for (int m = 1; m < BXM3; ++m) {
1255             dxytempdx[1][m][i] = dxytempdx[1][m][i + 1];
1256           }
1257         }
1258       }
1259     }
1260 
1261     //combine rows and columns to simulate double pixels

1262 
1263     for (int m = 1; m < BXM3; ++m) {
1264       if (xdouble[m]) {
1265         //  Combine the x-rows

1266         for (int n = 1; n < BYM3; ++n) {
1267           dxytempdx[1][m][n] += dxytempdx[1][m + 1][n];
1268         }
1269         //  Now shift the remaining pixels over by one row

1270         for (int j = m + 1; j < BXM3; ++j) {
1271           for (int n = 1; n < BYM3; ++n) {
1272             dxytempdx[1][j][n] = dxytempdx[1][j + 1][n];
1273           }
1274         }
1275       }
1276     }
1277 
1278     // Next do shifted -x template

1279 
1280     pixx = (int)floorf((xhit - deltaxy[0]) / xsize_);
1281     dx = (xhit - deltaxy[0]) - (pixx + 0.5f) * xsize_;
1282     if (flip_x_) {
1283       dx = -dx;
1284     }
1285     l0 = (int)(dx / xsize_ * 6.f + 3.5f);
1286     if (l0 < 0)
1287       l0 = 0;
1288     if (l0 > 6)
1289       l0 = 6;
1290     ddx = 6.f * dx / xsize_ - (l0 - 3);
1291     adx = fabs(ddx);
1292     if (ddx > 0.f) {
1293       l1 = l0 + 1;
1294       if (l1 > 6)
1295         l1 = l0 - 1;
1296     } else {
1297       l1 = l0 - 1;
1298       if (l1 < 0)
1299         l1 = l0 + 1;
1300     }
1301 
1302     // OK, lets do the template interpolation.

1303 
1304     // Calculate the x and y offsets to make the new template

1305 
1306     // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system

1307 
1308     ++pixx;
1309 
1310     // In the template store, the struck pixel is always (THy,THx)

1311 
1312     deltax = pixx - T2HX;
1313 
1314     // Loop over the non-zero part of the template index space and interpolate

1315 
1316     for (int j = jmin_; j <= jmax_; ++j) {
1317       // Flip indices as needed

1318       if (flip_x_) {
1319         jflipx = T2XSIZE - 1 - j;
1320         m = deltax + jflipx;
1321       } else {
1322         m = deltax + j;
1323       }
1324       for (int i = imin_; i <= imax_; ++i) {
1325         if (flip_y_) {
1326           iflipy = T2YSIZE - 1 - i;
1327           n = deltay + iflipy;
1328         } else {
1329           n = deltay + i;
1330         }
1331         if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1332           dxytempdx[0][m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1333                                adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1334                                ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1335                                adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1336                                adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1337         }
1338       }
1339     }
1340 
1341     //combine rows and columns to simulate double pixels

1342 
1343     for (int n = 1; n < BYM3; ++n) {
1344       if (ydouble[n]) {
1345         //  Combine the y-columns

1346         for (int m = 1; m < BXM3; ++m) {
1347           dxytempdx[0][m][n] += dxytempdx[0][m][n + 1];
1348         }
1349         //  Now shift the remaining pixels over by one column

1350         for (int i = n + 1; i < BYM3; ++i) {
1351           for (int m = 1; m < BXM3; ++m) {
1352             dxytempdx[0][m][i] = dxytempdx[0][m][i + 1];
1353           }
1354         }
1355       }
1356     }
1357 
1358     //combine rows and columns to simulate double pixels

1359 
1360     for (int m = 1; m < BXM3; ++m) {
1361       if (xdouble[m]) {
1362         //  Combine the x-rows

1363         for (int n = 1; n < BYM3; ++n) {
1364           dxytempdx[0][m][n] += dxytempdx[0][m + 1][n];
1365         }
1366         //  Now shift the remaining pixels over by one row

1367         for (int j = m + 1; j < BXM3; ++j) {
1368           for (int n = 1; n < BYM3; ++n) {
1369             dxytempdx[0][j][n] = dxytempdx[0][j + 1][n];
1370           }
1371         }
1372       }
1373     }
1374 
1375     //  Finally, normalize the derivatives and copy the results to the output array

1376 
1377     for (int n = 1; n < BYM3; ++n) {
1378       for (int m = 1; m < BXM3; ++m) {
1379         dpdx2d[0][m][n] = (dxytempdx[1][m][n] - dxytempdx[0][m][n]) / (2. * deltaxy[0]);
1380       }
1381     }
1382 
1383     // Next, do shifted y template

1384 
1385     pixy = (int)floorf((yhit + deltaxy[1]) / ysize_);
1386     dy = (yhit + deltaxy[1]) - (pixy + 0.5f) * ysize_;
1387     if (flip_y_) {
1388       dy = -dy;
1389     }
1390     k0 = (int)(dy / ysize_ * 6.f + 3.5f);
1391     if (k0 < 0)
1392       k0 = 0;
1393     if (k0 > 6)
1394       k0 = 6;
1395     ddy = 6.f * dy / ysize_ - (k0 - 3);
1396     ady = fabs(ddy);
1397     if (ddy > 0.f) {
1398       k1 = k0 + 1;
1399       if (k1 > 6)
1400         k1 = k0 - 1;
1401     } else {
1402       k1 = k0 - 1;
1403       if (k1 < 0)
1404         k1 = k0 + 1;
1405     }
1406     pixx = (int)floorf(xhit / xsize_);
1407     dx = xhit - (pixx + 0.5f) * xsize_;
1408     if (flip_x_) {
1409       dx = -dx;
1410     }
1411     l0 = (int)(dx / xsize_ * 6.f + 3.5f);
1412     if (l0 < 0)
1413       l0 = 0;
1414     if (l0 > 6)
1415       l0 = 6;
1416     ddx = 6.f * dx / xsize_ - (l0 - 3);
1417     adx = fabs(ddx);
1418     if (ddx > 0.f) {
1419       l1 = l0 + 1;
1420       if (l1 > 6)
1421         l1 = l0 - 1;
1422     } else {
1423       l1 = l0 - 1;
1424       if (l1 < 0)
1425         l1 = l0 + 1;
1426     }
1427 
1428     // OK, lets do the template interpolation.

1429 
1430     // Calculate the x and y offsets to make the new template

1431 
1432     // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system

1433 
1434     ++pixy;
1435     ++pixx;
1436 
1437     // In the template store, the struck pixel is always (THy,THx)

1438 
1439     deltax = pixx - T2HX;
1440     deltay = pixy - T2HY;
1441 
1442     // Loop over the non-zero part of the template index space and interpolate

1443 
1444     for (int j = jmin_; j <= jmax_; ++j) {
1445       // Flip indices as needed

1446       if (flip_x_) {
1447         jflipx = T2XSIZE - 1 - j;
1448         m = deltax + jflipx;
1449       } else {
1450         m = deltax + j;
1451       }
1452       for (int i = imin_; i <= imax_; ++i) {
1453         if (flip_y_) {
1454           iflipy = T2YSIZE - 1 - i;
1455           n = deltay + iflipy;
1456         } else {
1457           n = deltay + i;
1458         }
1459         if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1460           dxytempdy[1][m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1461                                adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1462                                ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1463                                adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1464                                adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1465         }
1466       }
1467     }
1468 
1469     //combine rows and columns to simulate double pixels

1470 
1471     for (int n = 1; n < BYM3; ++n) {
1472       if (ydouble[n]) {
1473         //  Combine the y-columns

1474         for (int m = 1; m < BXM3; ++m) {
1475           dxytempdy[1][m][n] += dxytempdy[1][m][n + 1];
1476         }
1477         //  Now shift the remaining pixels over by one column

1478         for (int i = n + 1; i < BYM3; ++i) {
1479           for (int m = 1; m < BXM3; ++m) {
1480             dxytempdy[1][m][i] = dxytempdy[1][m][i + 1];
1481           }
1482         }
1483       }
1484     }
1485 
1486     //combine rows and columns to simulate double pixels

1487 
1488     for (int m = 1; m < BXM3; ++m) {
1489       if (xdouble[m]) {
1490         //  Combine the x-rows

1491         for (int n = 1; n < BYM3; ++n) {
1492           dxytempdy[1][m][n] += dxytempdy[1][m + 1][n];
1493         }
1494         //  Now shift the remaining pixels over by one row

1495         for (int j = m + 1; j < BXM3; ++j) {
1496           for (int n = 1; n < BYM3; ++n) {
1497             dxytempdy[1][j][n] = dxytempdy[1][j + 1][n];
1498           }
1499         }
1500       }
1501     }
1502 
1503     // Next, do shifted -y template

1504 
1505     pixy = (int)floorf((yhit - deltaxy[1]) / ysize_);
1506     dy = (yhit - deltaxy[1]) - (pixy + 0.5f) * ysize_;
1507     if (flip_y_) {
1508       dy = -dy;
1509     }
1510     k0 = (int)(dy / ysize_ * 6.f + 3.5f);
1511     if (k0 < 0)
1512       k0 = 0;
1513     if (k0 > 6)
1514       k0 = 6;
1515     ddy = 6.f * dy / ysize_ - (k0 - 3);
1516     ady = fabs(ddy);
1517     if (ddy > 0.f) {
1518       k1 = k0 + 1;
1519       if (k1 > 6)
1520         k1 = k0 - 1;
1521     } else {
1522       k1 = k0 - 1;
1523       if (k1 < 0)
1524         k1 = k0 + 1;
1525     }
1526 
1527     // OK, lets do the template interpolation.

1528 
1529     // Calculate the x and y offsets to make the new template

1530 
1531     // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system

1532 
1533     ++pixy;
1534 
1535     // In the template store, the struck pixel is always (THy,THx)

1536 
1537     deltay = pixy - T2HY;
1538 
1539     // Loop over the non-zero part of the template index space and interpolate

1540 
1541     for (int j = jmin_; j <= jmax_; ++j) {
1542       // Flip indices as needed

1543       if (flip_x_) {
1544         jflipx = T2XSIZE - 1 - j;
1545         m = deltax + jflipx;
1546       } else {
1547         m = deltax + j;
1548       }
1549       for (int i = imin_; i <= imax_; ++i) {
1550         if (flip_y_) {
1551           iflipy = T2YSIZE - 1 - i;
1552           n = deltay + iflipy;
1553         } else {
1554           n = deltay + i;
1555         }
1556         if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1557           dxytempdy[0][m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1558                                adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1559                                ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1560                                adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1561                                adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1562         }
1563       }
1564     }
1565 
1566     //combine rows and columns to simulate double pixels

1567 
1568     for (int n = 1; n < BYM3; ++n) {
1569       if (ydouble[n]) {
1570         //  Combine the y-columns

1571         for (int m = 1; m < BXM3; ++m) {
1572           dxytempdy[0][m][n] += dxytempdy[0][m][n + 1];
1573         }
1574         //  Now shift the remaining pixels over by one column

1575         for (int i = n + 1; i < BYM3; ++i) {
1576           for (int m = 1; m < BXM3; ++m) {
1577             dxytempdy[0][m][i] = dxytempdy[0][m][i + 1];
1578           }
1579         }
1580       }
1581     }
1582 
1583     //combine rows and columns to simulate double pixels

1584 
1585     for (int m = 1; m < BXM3; ++m) {
1586       if (xdouble[m]) {
1587         //  Combine the x-rows

1588         for (int n = 1; n < BYM3; ++n) {
1589           dxytempdy[0][m][n] += dxytempdy[0][m + 1][n];
1590         }
1591         //  Now shift the remaining pixels over by one row

1592         for (int j = m + 1; j < BXM3; ++j) {
1593           for (int n = 1; n < BYM3; ++n) {
1594             dxytempdy[0][j][n] = dxytempdy[0][j + 1][n];
1595           }
1596         }
1597       }
1598     }
1599 
1600     //  Finally, normalize the derivatives and copy the results to the output array

1601 
1602     for (int n = 1; n < BYM3; ++n) {
1603       for (int m = 1; m < BXM3; ++m) {
1604         dpdx2d[1][m][n] = (dxytempdy[1][m][n] - dxytempdy[0][m][n]) / (2. * deltaxy[1]);
1605       }
1606     }
1607   }
1608 
1609   return success_;
1610 }  // xytemp

1611 
1612 // *************************************************************************************************************************************

1613 //! Interpolate stored 2-D information for input angles and hit position to make a 2-D template

1614 //! \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)

1615 //! \param       yhit - (input) y-position of hit relative to the lower left corner of pixel[1][1]

1616 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel starting at cluster[1][1]

1617 //! \param    xdouble - (input) STL vector of 11 element array to flag a double-pixel starting at cluster[1][1]

1618 //! \param template2d - (output) 2d template of size matched to the cluster.  Input must be zeroed since charge is added only.

1619 // *************************************************************************************************************************************

1620 
1621 bool SiPixelTemplate2D::xytemp(
1622     float xhit, float yhit, bool ydouble[BYM2], bool xdouble[BXM2], float template2d[BXM2][BYM2]) {
1623   // Interpolate for a new set of track angles

1624 
1625   bool derivatives = false;
1626   float dpdx2d[2][BXM2][BYM2];
1627   float QTemplate;
1628 
1629   return SiPixelTemplate2D::xytemp(xhit, yhit, ydouble, xdouble, template2d, derivatives, dpdx2d, QTemplate);
1630 
1631 }  // xytemp

1632 
1633 // *************************************************************************************************************************************

1634 //! Interpolate stored 2-D information for input angles and hit position to make a 2-D template

1635 //! \param         id - (input) the id of the template

1636 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)

1637 //! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)

1638 //! \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)

1639 //! \param       yhit - (input) y-position of hit relative to the lower left corner of pixel[1][1]

1640 //! \param    ydouble - (input) STL vector of 21 element array to flag a double-pixel starting at cluster[1][1]

1641 //! \param    xdouble - (input) STL vector of 11 element array to flag a double-pixel starting at cluster[1][1]

1642 //! \param template2d - (output) 2d template of size matched to the cluster.  Input must be zeroed since charge is added only.

1643 // *************************************************************************************************************************************

1644 
1645 bool SiPixelTemplate2D::xytemp(int id,
1646                                float cotalpha,
1647                                float cotbeta,
1648                                float xhit,
1649                                float yhit,
1650                                std::vector<bool>& ydouble,
1651                                std::vector<bool>& xdouble,
1652                                float template2d[BXM2][BYM2]) {
1653   // Local variables

1654 
1655   bool derivatives = false;
1656   float dpdx2d[2][BXM2][BYM2];
1657   float QTemplate;
1658   float locBx = 1.f;
1659   if (cotbeta < 0.f) {
1660     locBx = -1.f;
1661   }
1662   float locBz = locBx;
1663   if (cotalpha < 0.f) {
1664     locBz = -locBx;
1665   }
1666 
1667   bool yd[BYM2], xd[BXM2];
1668 
1669   yd[0] = false;
1670   yd[BYM2 - 1] = false;
1671   for (int i = 0; i < TYSIZE; ++i) {
1672     yd[i + 1] = ydouble[i];
1673   }
1674   xd[0] = false;
1675   xd[BXM2 - 1] = false;
1676   for (int j = 0; j < TXSIZE; ++j) {
1677     xd[j + 1] = xdouble[j];
1678   }
1679 
1680   // Interpolate for a new set of track angles

1681 
1682   if (SiPixelTemplate2D::interpolate(id, cotalpha, cotbeta, locBz, locBx)) {
1683     return SiPixelTemplate2D::xytemp(xhit, yhit, yd, xd, template2d, derivatives, dpdx2d, QTemplate);
1684   } else {
1685     return false;
1686   }
1687 
1688 }  // xytemp

1689 
1690 // ************************************************************************************************************

1691 //! Return y error (squared) for an input signal and yindex

1692 //! Add large Q scaling for use in cluster splitting.

1693 //! \param qpixel - (input) pixel charge

1694 //! \param index - (input) y-index index of pixel

1695 //! \param xysig2 - (output) square error

1696 // ************************************************************************************************************

1697 void SiPixelTemplate2D::xysigma2(float qpixel, int index, float& xysig2)
1698 
1699 {
1700   // Interpolate using quantities already stored in the private variables

1701 
1702   // Local variables

1703   float sigi, sigi2, sigi3, sigi4, qscale, err2, err00;
1704 
1705   // Make sure that input is OK

1706 
1707 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1708   if (index < 1 || index >= BYM2) {
1709     throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::ysigma2 called with index = " << index << std::endl;
1710   }
1711 #else
1712   assert(index > 0 && index < BYM2);
1713 #endif
1714 
1715   // Define the maximum signal to use in the parameterization

1716 
1717   // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis

1718 
1719   if (qpixel < sxymax_) {
1720     sigi = qpixel;
1721     qscale = 1.f;
1722   } else {
1723     sigi = sxymax_;
1724     qscale = qpixel / sxymax_;
1725   }
1726   sigi2 = sigi * sigi;
1727   sigi3 = sigi2 * sigi;
1728   sigi4 = sigi3 * sigi;
1729   if (index <= T2HYP1) {
1730     err00 = xypary0x0_[0][0] + xypary0x0_[0][1] * sigi + xypary0x0_[0][2] * sigi2 + xypary0x0_[0][3] * sigi3 +
1731             xypary0x0_[0][4] * sigi4;
1732     err2 = err00 +
1733            adcota_ * (xypary0x1_[0][0] + xypary0x1_[0][1] * sigi + xypary0x1_[0][2] * sigi2 + xypary0x1_[0][3] * sigi3 +
1734                       xypary0x1_[0][4] * sigi4 - err00) +
1735            adcotb_ * (xypary1x0_[0][0] + xypary1x0_[0][1] * sigi + xypary1x0_[0][2] * sigi2 + xypary1x0_[0][3] * sigi3 +
1736                       xypary1x0_[0][4] * sigi4 - err00);
1737   } else {
1738     err00 = xypary0x0_[1][0] + xypary0x0_[1][1] * sigi + xypary0x0_[1][2] * sigi2 + xypary0x0_[1][3] * sigi3 +
1739             xypary0x0_[1][4] * sigi4;
1740     err2 = err00 +
1741            adcota_ * (xypary0x1_[1][0] + xypary0x1_[1][1] * sigi + xypary0x1_[1][2] * sigi2 + xypary0x1_[1][3] * sigi3 +
1742                       xypary0x1_[1][4] * sigi4 - err00) +
1743            adcotb_ * (xypary1x0_[1][0] + xypary1x0_[1][1] * sigi + xypary1x0_[1][2] * sigi2 + xypary1x0_[1][3] * sigi3 +
1744                       xypary1x0_[1][4] * sigi4 - err00);
1745   }
1746   xysig2 = qscale * err2;
1747   if (xysig2 <= 0.f) {
1748     xysig2 = s50_ * s50_;
1749   }
1750 
1751   return;
1752 
1753 }  // End xysigma2

1754 
1755 // ************************************************************************************************************

1756 //! Return the Landau probability parameters for this set of cot(alpha, cot(beta)

1757 // ************************************************************************************************************

1758 void SiPixelTemplate2D::landau_par(float lanpar[2][5])
1759 
1760 {
1761   // Interpolate using quantities already stored in the private variables

1762 
1763   for (int i = 0; i < 2; ++i) {
1764     for (int j = 0; j < 5; ++j) {
1765       lanpar[i][j] = lanpar_[i][j];
1766     }
1767   }
1768   return;
1769 
1770 }  // End lan_par
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!