Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-01-25 23:30:12 

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.resize(theCurrentTemp.head.NTyx, theCurrentTemp.head.NTxx);
0171 
0172     // Read in the file info

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

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

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

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

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

0303 
0304 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0305 
0306 //****************************************************************

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

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

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

0310 //****************************************************************

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

0314 
0315   const int code_version = {21};
0316 
0317   // We must use a reader because dbobject must be a const and its stream must not be

0318   SiPixel2DTemplateDBObject::Reader reader(dbobject);
0319 
0320   struct Failed {};
0321 
0322   // Fill the template storage for each template calibration stored in the db

0323   for (int m = 0; m < dbobject.numOfTempl(); ++m) {
0324     // Create a template storage entry

0325     // SiPixelTemplateStore2D

0326     pixelTemp.emplace_back();
0327     auto& theCurrentTemp = pixelTemp.back();
0328 
0329     try {
0330       // Read-in a header string first and print it

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

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

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

0410       theCurrentTemp.resize(theCurrentTemp.head.NTyx, theCurrentTemp.head.NTxx);
0411 
0412       // Read in the file info

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

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

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

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

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

0550 
0551         Dtype_ = thePixelTemp_[index_id_].head.Dtype;
0552 
0553         // Copy the charge scaling factor to the private variable

0554 
0555         qscale_ = thePixelTemp_[index_id_].head.qscale;
0556 
0557         // Copy the pseudopixel signal size to the private variable

0558 
0559         s50_ = thePixelTemp_[index_id_].head.s50;
0560 
0561         // Copy Qbinning info to private variables

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

0568 
0569         lorywidth_ = thePixelTemp_[index_id_].head.lorywidth;
0570         lorxwidth_ = thePixelTemp_[index_id_].head.lorxwidth;
0571 
0572         // Copy the pixel sizes private variables

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

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

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

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

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

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

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

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

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

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

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

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

0627 // *************************************************************************************************************************************

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

0631 
0632   //check for nan's

0633   if (!edm::isFinite(cotalpha) || !edm::isFinite(cotbeta)) {
0634     success_ = false;
0635     return success_;
0636   }
0637 
0638   // Local variables

0639 
0640   float acotb, dcota, dcotb;
0641 
0642   // Check to see if interpolation is valid

0643 
0644   if (id != id_current_ || cotalpha != cota_current_ || cotbeta != cotb_current_) {
0645     cota_current_ = cotalpha;
0646     cotb_current_ = cotbeta;
0647     // Try to find the correct template.  Fill the class variable index_id_ .

0648     success_ = getid(id);
0649   }
0650 
0651 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0652   if (index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
0653     throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::interpolate can't find needed template ID = " << id
0654                                         << ", Are you using the correct global tag?" << std::endl;
0655   }
0656 #else
0657   assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
0658 #endif
0659 
0660   // Check angle limits and et up interpolation parameters

0661 
0662   float cota = cotalpha;
0663   flip_x_ = false;
0664   flip_y_ = false;
0665   switch (Dtype_) {
0666     case 0:
0667       if (cotbeta < 0.f) {
0668         flip_y_ = true;
0669       }
0670       break;
0671     case 1:
0672       if (locBz > 0.f) {
0673         flip_y_ = true;
0674       }
0675       break;
0676     case 2:
0677     case 3:
0678     case 4:
0679     case 5:
0680       if (locBx * locBz < 0.f) {
0681         cota = std::abs(cotalpha);
0682         flip_x_ = true;
0683       }
0684       if (locBx < 0.f) {
0685         flip_y_ = true;
0686       }
0687       break;
0688     default:
0689 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0690       throw cms::Exception("DataCorrupt")
0691           << "SiPixelTemplate2D::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
0692 #else
0693       std::cout << "SiPixelTemplate:2D:illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
0694 #endif
0695   }
0696 
0697   if (cota < cotalpha0_) {
0698     success_ = false;
0699     jx0_ = 0;
0700     jx1_ = 1;
0701     adcota_ = 0.f;
0702   } else if (cota > cotalpha1_) {
0703     success_ = false;
0704     jx0_ = Nxx_ - 1;
0705     jx1_ = jx0_ - 1;
0706     adcota_ = 0.f;
0707   } else {
0708     jx0_ = (int)((cota - cotalpha0_) / deltacota_ + 0.5f);
0709     dcota = (cota - (cotalpha0_ + jx0_ * deltacota_)) / deltacota_;
0710     adcota_ = fabs(dcota);
0711     if (dcota > 0.f) {
0712       jx1_ = jx0_ + 1;
0713       if (jx1_ > Nxx_ - 1)
0714         jx1_ = jx0_ - 1;
0715     } else {
0716       jx1_ = jx0_ - 1;
0717       if (jx1_ < 0)
0718         jx1_ = jx0_ + 1;
0719     }
0720   }
0721 
0722   // Interpolate the absolute value of cot(beta)

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

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

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

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

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

0854 
0855 // *************************************************************************************************************************************

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

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

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

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

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

0861 // *************************************************************************************************************************************

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

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

0893 
0894   adcota_ = 0.f;
0895   adcotb_ = 0.f;
0896 
0897   // Interpolate things in cot(alpha)-cot(beta)

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

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

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

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

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

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

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

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

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

0990 // *************************************************************************************************************************************

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

1001 
1002   // Local variables

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

1007   const float deltaxy[2] = {16.67f, 25.0f};
1008 
1009   // Check to see if interpolation is valid

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

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

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

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

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

1060 
1061   // First find the limits of the indices for non-zero pixels

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

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

1078 
1079   ++pixy;
1080   ++pixx;
1081 
1082   // In the template store, the struck pixel is always (THy,THx)

1083 
1084   deltax = pixx - T2HX;
1085   deltay = pixy - T2HY;
1086 
1087   //  First zero the local 2-d template

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

1096 
1097   for (int j = jmin_; j <= jmax_; ++j) {
1098     // Flip indices as needed

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

1123 
1124   for (int n = 1; n < BYM3; ++n) {
1125     if (ydouble[n]) {
1126       //  Combine the y-columns

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

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

1140 
1141   for (int m = 1; m < BXM3; ++m) {
1142     if (xdouble[m]) {
1143       //  Combine the x-rows

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

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

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

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

1209 
1210     // Calculate the x and y offsets to make the new template

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

1213 
1214     ++pixx;
1215 
1216     // In the template store, the struck pixel is always (THy,THx)

1217 
1218     deltax = pixx - T2HX;
1219 
1220     // Loop over the non-zero part of the template index space and interpolate

1221 
1222     for (int j = jmin_; j <= jmax_; ++j) {
1223       // Flip indices as needed

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

1248 
1249     for (int n = 1; n < BYM3; ++n) {
1250       if (ydouble[n]) {
1251         //  Combine the y-columns

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

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

1265 
1266     for (int m = 1; m < BXM3; ++m) {
1267       if (xdouble[m]) {
1268         //  Combine the x-rows

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

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

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

1306 
1307     // Calculate the x and y offsets to make the new template

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

1310 
1311     ++pixx;
1312 
1313     // In the template store, the struck pixel is always (THy,THx)

1314 
1315     deltax = pixx - T2HX;
1316 
1317     // Loop over the non-zero part of the template index space and interpolate

1318 
1319     for (int j = jmin_; j <= jmax_; ++j) {
1320       // Flip indices as needed

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

1345 
1346     for (int n = 1; n < BYM3; ++n) {
1347       if (ydouble[n]) {
1348         //  Combine the y-columns

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

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

1362 
1363     for (int m = 1; m < BXM3; ++m) {
1364       if (xdouble[m]) {
1365         //  Combine the x-rows

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

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

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

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

1432 
1433     // Calculate the x and y offsets to make the new template

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

1436 
1437     ++pixy;
1438     ++pixx;
1439 
1440     // In the template store, the struck pixel is always (THy,THx)

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

1446 
1447     for (int j = jmin_; j <= jmax_; ++j) {
1448       // Flip indices as needed

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

1473 
1474     for (int n = 1; n < BYM3; ++n) {
1475       if (ydouble[n]) {
1476         //  Combine the y-columns

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

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

1490 
1491     for (int m = 1; m < BXM3; ++m) {
1492       if (xdouble[m]) {
1493         //  Combine the x-rows

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

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

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

1531 
1532     // Calculate the x and y offsets to make the new template

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

1535 
1536     ++pixy;
1537 
1538     // In the template store, the struck pixel is always (THy,THx)

1539 
1540     deltay = pixy - T2HY;
1541 
1542     // Loop over the non-zero part of the template index space and interpolate

1543 
1544     for (int j = jmin_; j <= jmax_; ++j) {
1545       // Flip indices as needed

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

1570 
1571     for (int n = 1; n < BYM3; ++n) {
1572       if (ydouble[n]) {
1573         //  Combine the y-columns

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

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

1587 
1588     for (int m = 1; m < BXM3; ++m) {
1589       if (xdouble[m]) {
1590         //  Combine the x-rows

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

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

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

1614 
1615 // *************************************************************************************************************************************

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

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

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

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

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

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

1622 // *************************************************************************************************************************************

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

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

1635 
1636 // *************************************************************************************************************************************

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

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

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

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

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

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

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

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

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

1646 // *************************************************************************************************************************************

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

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

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

1692 
1693 // ************************************************************************************************************

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

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

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

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

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

1699 // ************************************************************************************************************

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

1704 
1705   // Local variables

1706   float sigi, sigi2, sigi3, sigi4, qscale, err2, err00;
1707 
1708   // Make sure that input is OK

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

1719 
1720   // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis

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

1757 
1758 // ************************************************************************************************************

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

1760 // ************************************************************************************************************

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

1765 
1766   for (int i = 0; i < 2; ++i) {
1767     for (int j = 0; j < 5; ++j) {
1768       lanpar[i][j] = lanpar_[i][j];
1769     }
1770   }
1771   return;
1772 
1773 }  // 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!