Reference Manual

Go to the documentation of this file.
00001 //
00002 //  SiPixelTemplate.cc  Version 8.30 
00003 //
00004 //  Add goodness-of-fit info and spare entries to templates, version number in template header, more error checking
00005 //  Add correction for (Q_F-Q_L)/(Q_F+Q_L) bias
00006 //  Add cot(beta) reflection to reduce y-entries and more sophisticated x-interpolation
00007 //  Fix small index searching bug in interpolate method
00008 //  Change interpolation indexing to avoid complier complaining about possible un-initialized variables
00009 //  Replace containers with static arrays in calls to ysigma2 and xsigma2
00010 //  Add external threshold to calls to ysigma2 and xsigma2, fix parameter signal max for xsigma2
00011 //  Return to 5 pixel spanning but adjust boundaries to use only when needed
00012 //  Implement improved (faster) chi2min search that depends on pixel types
00013 //  Fill template arrays in single calls to this object
00014 //  Add qmin to the template
00015 //  Add qscale to match charge scales
00016 //  Small improvement to x-chisquare interpolation
00017 //  Enlarge SiPixelTemplateStore to accommodate larger templates and increased alpha acceptance (reduce PT threshold to ~200 MeV)
00018 //  Change output streams to conform to CMSSW info and error logging.
00019 //  Store x and y cluster sizes in fractional pixels to facilitate cluster splitting
00020 //  Add methods to return 3-d templates needed for cluster splitting
00021 //  Keep interpolated central 9 template bins in private storage and expand/shift in the getter functions (faster for speed=2/3) and easier to build 3d templates
00022 //  Store error and bias information for the simple chi^2 min position analysis (no interpolation or Q_{FB} corrections) to use in cluster splitting
00023 //  To save time, the gaussian centers and sigma are not interpolated right now (they aren't currently used).  They can be restored by un-commenting lines in the interpolate method.
00024 //  Add a new method to calculate qbin for input cotbeta and cluster charge.  To be used for error estimation of merged clusters in PixelCPEGeneric.
00025 //  Improve the charge estimation for larger cot(alpha) tracks
00026 //  Change interpolate method to return false boolean if track angles are outside of range
00027 //  Add template info and method for truncation information
00028 //  Change to allow template sizes to be changed at compile time
00029 //  Fix bug in track angle checking
00030 //  Accommodate Dave's new DB pushfile which overloads the old method (file input)
00031 //  Add CPEGeneric error information and expand qbin method to access useful info for PixelCPEGeneric
00032 //  Fix large cot(alpha) bug in qmin interpolation
00033 //  Add second qmin to allow a qbin=5 state
00034 //  Use interpolated chi^2 info for one-pixel clusters
00035 //  Fix DB pushfile version number checking bug.
00036 //  Remove assert from qbin method
00037 //  Replace asserts with exceptions in CMSSW
00038 //  Change calling sequence to interpolate method to handle cot(beta)<0 for FPix cosmics
00039 //  Add getter for pixelav Lorentz width estimates to qbin method
00040 //  Add check on template size to interpolate and qbin methods
00041 //  Add qbin population information, charge distribution information
00042 //
00043 //
00044 //  V7.00 - Decouple BPix and FPix information into separate templates
00045 //  Add methods to facilitate improved cluster splitting
00046 //  Fix small charge scaling bug (affects FPix only)
00047 //  Change y-slice used for the x-template to be closer to the actual cotalpha-cotbeta point 
00048 //  (there is some weak breakdown of x-y factorization in the FPix after irradiation)
00049 //
00050 //
00051 //  V8.00 - Add method to calculate a simple 2D template
00052 //  Reorganize the interpolate method to extract header info only once per ID
00053 //  V8.01 - Improve simple template normalization
00054 //  V8.05 - Change qbin normalization to work better after irradiation
00055 //  V8.10 - Add Vavilov distribution interpolation
00056 //  V8.11 - Renormalize the x-templates for Guofan's cluster size calculation
00057 //  V8.12 - Technical fix to qavg issue.
00058 //  V8.13 - Fix qbin and fastsim interpolaters to avoid changing class variables
00059 //  V8.20 - Add methods to identify the central pixels in the x- and y-templates (to help align templates with clusters in radiation damaged detectors)
00060 //          Rename class variables from pxxxx (private xxxx) to xxxx_ to follow standard convention.
00061 //          Add compiler option to store the template entries in BOOST multiarrays of structs instead of simple c arrays
00062 //          (allows dynamic resizing template storage and has bounds checking but costs ~10% in execution time).
00063 //  V8.21 - Add new qbin method to use in cluster splitting
00064 //  V8.23 - Replace chi-min position errors with merged cluster chi2 probability info
00065 //  V8.25 - Incorporate VI's speed changes into the current version
00066 //  V8.26 - Modify the Vavilov lookups to work better with the FPix (offset y-templates)
00067 //  V8.30 - Change the splitting template generation and access to improve speed and eliminate triple index boost::multiarray
00068 
00069 //  Created by Morris Swartz on 10/27/06.
00070 //
00071 //
00072 
00073 //#include <stdlib.h> 
00074 //#include <stdio.h>
00075 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00076 //#include <cmath.h>
00077 #else
00078 #include <math.h>
00079 #endif
00080 #include <algorithm>
00081 #include <vector>
00082 //#include "boost/multi_array.hpp"
00083 #include <iostream>
00084 #include <iomanip>
00085 #include <sstream>
00086 #include <fstream>
00087 #include <list>
00088 
00089 
00090 
00091 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00092 #include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplate.h"
00093 #include "RecoLocalTracker/SiPixelRecHits/interface/SimplePixel.h"
00094 #include "FWCore/ParameterSet/interface/FileInPath.h"
00095 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00096 #define LOGERROR(x) LogError(x)
00097 #define LOGINFO(x) LogInfo(x)
00098 #define ENDL " "
00099 #include "FWCore/Utilities/interface/Exception.h"
00100 using namespace edm;
00101 #else
00102 #include "SiPixelTemplate.h"
00103 #include "SimplePixel.h"
00104 #define LOGERROR(x) std::cout << x << ": "
00105 #define LOGINFO(x) std::cout << x << ": "
00106 #define ENDL std::endl
00107 #endif
00108 
00109 //**************************************************************** 
00114 //**************************************************************** 
00115 bool SiPixelTemplate::pushfile(int filenum)
00116 {
00117     // Add template stored in external file numbered filenum to theTemplateStore
00118     
00119     // Local variables 
00120     int i, j, k, l;
00121         float qavg_avg;
00122         const char *tempfile;
00123         //      char title[80]; remove this
00124     char c;
00125         const int code_version={17};
00126         
00127         
00128         
00129         //  Create a filename for this run 
00130         
00131         std::ostringstream tout;
00132         
00133         //  Create different path in CMSSW than standalone
00134         
00135 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00136         tout << "CalibTracker/SiPixelESProducers/data/template_summary_zp" 
00137         << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
00138         std::string tempf = tout.str();
00139         edm::FileInPath file( tempf.c_str() );
00140         tempfile = (file.fullPath()).c_str();
00141 #else
00142         tout << "template_summary_zp" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
00143         std::string tempf = tout.str();
00144         tempfile = tempf.c_str();
00145 #endif
00146         
00147         //  open the template file 
00148         
00149         std::ifstream in_file(tempfile, std::ios::in);
00150         
00151         if(in_file.is_open()) {
00152                 
00153                 // Create a local template storage entry
00154                 
00155                 SiPixelTemplateStore theCurrentTemp;
00156                 
00157                 // Read-in a header string first and print it    
00158                 
00159                 for (i=0; (c=in_file.get()) != '\n'; ++i) {
00160                         if(i < 79) {theCurrentTemp.head.title[i] = c;}
00161                 }
00162                 if(i > 78) {i=78;}
00163                 theCurrentTemp.head.title[i+1] ='\0';
00164                 LOGINFO("SiPixelTemplate") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
00165                 
00166                 // next, the header information     
00167                 
00168                 in_file >> theCurrentTemp.head.ID  >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx
00169                 >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
00170                 >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >> theCurrentTemp.head.zsize;
00171                 
00172                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file, no template load" << ENDL; return false;}
00173                 
00174                 LOGINFO("SiPixelTemplate") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version " << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield 
00175                 << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
00176                 << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
00177                 << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
00178                 << ", 1/2 threshold " << theCurrentTemp.head.s50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", x Lorentz width " << theCurrentTemp.head.lorxwidth    
00179                 << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
00180                 
00181                 if(theCurrentTemp.head.templ_version < code_version) {LOGERROR("SiPixelTemplate") << "code expects version " << code_version << ", no template load" << ENDL; return false;}
00182                 
00183 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST 
00184                 
00185 // next, layout the 1-d/2-d structures needed to store template
00186                                 
00187                 theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
00188 
00189                 theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
00190                 
00191 #endif
00192                 
00193 // next, loop over all y-angle entries   
00194                 
00195                 for (i=0; i < theCurrentTemp.head.NTy; ++i) {     
00196                         
00197                         in_file >> theCurrentTemp.enty[i].runnum >> theCurrentTemp.enty[i].costrk[0] 
00198                         >> theCurrentTemp.enty[i].costrk[1] >> theCurrentTemp.enty[i].costrk[2]; 
00199                         
00200                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 1, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00201                         
00202                         // Calculate the alpha, beta, and cot(beta) for this entry 
00203                         
00204                         theCurrentTemp.enty[i].alpha = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[0]));
00205                         
00206                         theCurrentTemp.enty[i].cotalpha = theCurrentTemp.enty[i].costrk[0]/theCurrentTemp.enty[i].costrk[2];
00207                         
00208                         theCurrentTemp.enty[i].beta = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[1]));
00209                         
00210                         theCurrentTemp.enty[i].cotbeta = theCurrentTemp.enty[i].costrk[1]/theCurrentTemp.enty[i].costrk[2];
00211                         
00212                         in_file >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].pixmax >> theCurrentTemp.enty[i].symax >> theCurrentTemp.enty[i].dyone
00213                         >> theCurrentTemp.enty[i].syone >> theCurrentTemp.enty[i].sxmax >> theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone;
00214                         
00215                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 2, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00216                         
00217                         in_file >> theCurrentTemp.enty[i].dytwo >> theCurrentTemp.enty[i].sytwo >> theCurrentTemp.enty[i].dxtwo 
00218                         >> theCurrentTemp.enty[i].sxtwo >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].clsleny >> theCurrentTemp.enty[i].clslenx;
00219                         
00220                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 3, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00221                         
00222                         for (j=0; j<2; ++j) {
00223                                 
00224                                 in_file >> theCurrentTemp.enty[i].ypar[j][0] >> theCurrentTemp.enty[i].ypar[j][1] 
00225                                 >> theCurrentTemp.enty[i].ypar[j][2] >> theCurrentTemp.enty[i].ypar[j][3] >> theCurrentTemp.enty[i].ypar[j][4];
00226                                 
00227                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 4, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00228                                 
00229                         }
00230                         
00231                         for (j=0; j<9; ++j) {
00232                                 
00233                                 for (k=0; k<TYSIZE; ++k) {in_file >> theCurrentTemp.enty[i].ytemp[j][k];}
00234                                 
00235                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 5, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00236                         }
00237                         
00238                         for (j=0; j<2; ++j) {
00239                                 
00240                                 in_file >> theCurrentTemp.enty[i].xpar[j][0] >> theCurrentTemp.enty[i].xpar[j][1] 
00241                                 >> theCurrentTemp.enty[i].xpar[j][2] >> theCurrentTemp.enty[i].xpar[j][3] >> theCurrentTemp.enty[i].xpar[j][4];
00242                                 
00243                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 6, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00244                                 
00245                         }
00246                         
00247                         qavg_avg = 0.f;
00248                         for (j=0; j<9; ++j) {
00249                                 
00250                                 for (k=0; k<TXSIZE; ++k) {in_file >> theCurrentTemp.enty[i].xtemp[j][k]; qavg_avg += theCurrentTemp.enty[i].xtemp[j][k];} 
00251                                 
00252                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 7, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00253                         }
00254                         theCurrentTemp.enty[i].qavg_avg = qavg_avg/9.;
00255                         
00256                         for (j=0; j<4; ++j) {
00257                                 
00258                                 in_file >> theCurrentTemp.enty[i].yavg[j] >> theCurrentTemp.enty[i].yrms[j] >> theCurrentTemp.enty[i].ygx0[j] >> theCurrentTemp.enty[i].ygsig[j];
00259                                 
00260                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 8, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00261                         }
00262                         
00263                         for (j=0; j<4; ++j) {
00264                                 
00265                                 in_file >> theCurrentTemp.enty[i].yflpar[j][0] >> theCurrentTemp.enty[i].yflpar[j][1] >> theCurrentTemp.enty[i].yflpar[j][2] 
00266                                 >> theCurrentTemp.enty[i].yflpar[j][3] >> theCurrentTemp.enty[i].yflpar[j][4] >> theCurrentTemp.enty[i].yflpar[j][5];
00267                                 
00268                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 9, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00269                         }
00270                         
00271                         for (j=0; j<4; ++j) {
00272                                 
00273                                 in_file >> theCurrentTemp.enty[i].xavg[j] >> theCurrentTemp.enty[i].xrms[j] >> theCurrentTemp.enty[i].xgx0[j] >> theCurrentTemp.enty[i].xgsig[j];
00274                                 
00275                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 10, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00276                         }
00277                         
00278                         for (j=0; j<4; ++j) {
00279                                 
00280                                 in_file >> theCurrentTemp.enty[i].xflpar[j][0] >> theCurrentTemp.enty[i].xflpar[j][1] >> theCurrentTemp.enty[i].xflpar[j][2] 
00281                                 >> theCurrentTemp.enty[i].xflpar[j][3] >> theCurrentTemp.enty[i].xflpar[j][4] >> theCurrentTemp.enty[i].xflpar[j][5];
00282                                 
00283                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 11, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00284                         }
00285                         
00286                         for (j=0; j<4; ++j) {
00287                                 
00288                                 in_file >> theCurrentTemp.enty[i].chi2yavg[j] >> theCurrentTemp.enty[i].chi2ymin[j] >> theCurrentTemp.enty[i].chi2xavg[j] >> theCurrentTemp.enty[i].chi2xmin[j];
00289                                 
00290                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 12, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00291                         }
00292                         
00293                         for (j=0; j<4; ++j) {
00294                                 
00295                                 in_file >> theCurrentTemp.enty[i].yavgc2m[j] >> theCurrentTemp.enty[i].yrmsc2m[j] >> theCurrentTemp.enty[i].chi2yavgc2m[j] >> theCurrentTemp.enty[i].chi2yminc2m[j];
00296                                 
00297                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 13, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00298                         }
00299                         
00300                         for (j=0; j<4; ++j) {
00301                                 
00302                                 in_file >> theCurrentTemp.enty[i].xavgc2m[j] >> theCurrentTemp.enty[i].xrmsc2m[j] >> theCurrentTemp.enty[i].chi2xavgc2m[j] >> theCurrentTemp.enty[i].chi2xminc2m[j];
00303                                 
00304                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00305                         } 
00306                         
00307                         for (j=0; j<4; ++j) {
00308                                 
00309                                 in_file >> theCurrentTemp.enty[i].yavggen[j] >> theCurrentTemp.enty[i].yrmsgen[j] >> theCurrentTemp.enty[i].ygx0gen[j] >> theCurrentTemp.enty[i].ygsiggen[j];
00310                                 
00311                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14a, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00312                         }
00313                         
00314                         for (j=0; j<4; ++j) {
00315                                 
00316                                 in_file >> theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j] >> theCurrentTemp.enty[i].xgx0gen[j] >> theCurrentTemp.enty[i].xgsiggen[j];
00317                                 
00318                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14b, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00319                         } 
00320                         
00321                         in_file >> theCurrentTemp.enty[i].chi2yavgone >> theCurrentTemp.enty[i].chi2yminone >> theCurrentTemp.enty[i].chi2xavgone >> theCurrentTemp.enty[i].chi2xminone >> theCurrentTemp.enty[i].qmin2
00322                         >> theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav >> theCurrentTemp.enty[i].qavg_spare >> theCurrentTemp.enty[i].spare[0];
00323                         
00324                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 15, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00325                         
00326                         in_file >> theCurrentTemp.enty[i].mpvvav2 >> theCurrentTemp.enty[i].sigmavav2 >> theCurrentTemp.enty[i].kappavav2 >> theCurrentTemp.enty[i].qbfrac[0] >> theCurrentTemp.enty[i].qbfrac[1]
00327                         >> theCurrentTemp.enty[i].qbfrac[2] >> theCurrentTemp.enty[i].fracyone >> theCurrentTemp.enty[i].fracxone >> theCurrentTemp.enty[i].fracytwo >> theCurrentTemp.enty[i].fracxtwo;
00328                         //              theCurrentTemp.enty[i].qbfrac[3] = 1. - theCurrentTemp.enty[i].qbfrac[0] - theCurrentTemp.enty[i].qbfrac[1] - theCurrentTemp.enty[i].qbfrac[2];
00329                         
00330                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 16, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00331                         
00332                 }
00333                 
00334                 // next, loop over all barrel x-angle entries   
00335                 
00336                 for (k=0; k < theCurrentTemp.head.NTyx; ++k) { 
00337                         
00338                         for (i=0; i < theCurrentTemp.head.NTxx; ++i) { 
00339                                 
00340                                 in_file >> theCurrentTemp.entx[k][i].runnum >> theCurrentTemp.entx[k][i].costrk[0] 
00341                                 >> theCurrentTemp.entx[k][i].costrk[1] >> theCurrentTemp.entx[k][i].costrk[2]; 
00342                                 
00343                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 17, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00344                                 
00345                                 // Calculate the alpha, beta, and cot(beta) for this entry 
00346                                 
00347                                 theCurrentTemp.entx[k][i].alpha = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[0]));
00348                                 
00349                                 theCurrentTemp.entx[k][i].cotalpha = theCurrentTemp.entx[k][i].costrk[0]/theCurrentTemp.entx[k][i].costrk[2];
00350                                 
00351                                 theCurrentTemp.entx[k][i].beta = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[1]));
00352                                 
00353                                 theCurrentTemp.entx[k][i].cotbeta = theCurrentTemp.entx[k][i].costrk[1]/theCurrentTemp.entx[k][i].costrk[2];
00354                                 
00355                                 in_file >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].pixmax >> theCurrentTemp.entx[k][i].symax >> theCurrentTemp.entx[k][i].dyone
00356                                 >> theCurrentTemp.entx[k][i].syone >> theCurrentTemp.entx[k][i].sxmax >> theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone;
00357                                 
00358                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 18, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00359                                 
00360                                 in_file >> theCurrentTemp.entx[k][i].dytwo >> theCurrentTemp.entx[k][i].sytwo >> theCurrentTemp.entx[k][i].dxtwo 
00361                                 >> theCurrentTemp.entx[k][i].sxtwo >> theCurrentTemp.entx[k][i].qmin >> theCurrentTemp.entx[k][i].clsleny >> theCurrentTemp.entx[k][i].clslenx;
00362                                 //                         >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav;
00363                                 
00364                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 19, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00365                                 
00366                                 for (j=0; j<2; ++j) {
00367                                         
00368                                         in_file >> theCurrentTemp.entx[k][i].ypar[j][0] >> theCurrentTemp.entx[k][i].ypar[j][1] 
00369                                         >> theCurrentTemp.entx[k][i].ypar[j][2] >> theCurrentTemp.entx[k][i].ypar[j][3] >> theCurrentTemp.entx[k][i].ypar[j][4];
00370                                         
00371                                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 20, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00372                                 }
00373                                 
00374                                 for (j=0; j<9; ++j) {
00375                                         
00376                                         for (l=0; l<TYSIZE; ++l) {in_file >> theCurrentTemp.entx[k][i].ytemp[j][l];} 
00377                                         
00378                                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 21, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00379                                 }
00380                                 
00381                                 for (j=0; j<2; ++j) {
00382                                         
00383                                         in_file >> theCurrentTemp.entx[k][i].xpar[j][0] >> theCurrentTemp.entx[k][i].xpar[j][1] 
00384                                         >> theCurrentTemp.entx[k][i].xpar[j][2] >> theCurrentTemp.entx[k][i].xpar[j][3] >> theCurrentTemp.entx[k][i].xpar[j][4];
00385                                         
00386                                         
00387                                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 22, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00388                                 }
00389                                 
00390                                 qavg_avg = 0.f;
00391                                 for (j=0; j<9; ++j) {
00392                                         
00393                                         for (l=0; l<TXSIZE; ++l) {in_file >> theCurrentTemp.entx[k][i].xtemp[j][l]; qavg_avg += theCurrentTemp.entx[k][i].xtemp[j][l];} 
00394                                         
00395                                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 23, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00396                                 }
00397                                 theCurrentTemp.entx[k][i].qavg_avg = qavg_avg/9.;
00398                                 
00399                                 for (j=0; j<4; ++j) {
00400                                         
00401                                         in_file >> theCurrentTemp.entx[k][i].yavg[j] >> theCurrentTemp.entx[k][i].yrms[j] >> theCurrentTemp.entx[k][i].ygx0[j] >> theCurrentTemp.entx[k][i].ygsig[j];
00402                                         
00403                                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 24, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00404                                 }
00405                                 
00406                                 for (j=0; j<4; ++j) {
00407                                         
00408                                         in_file >> theCurrentTemp.entx[k][i].yflpar[j][0] >> theCurrentTemp.entx[k][i].yflpar[j][1] >> theCurrentTemp.entx[k][i].yflpar[j][2] 
00409                                         >> theCurrentTemp.entx[k][i].yflpar[j][3] >> theCurrentTemp.entx[k][i].yflpar[j][4] >> theCurrentTemp.entx[k][i].yflpar[j][5];
00410                                         
00411                                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 25, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00412                                 }
00413                                 
00414                                 for (j=0; j<4; ++j) {
00415                                         
00416                                         in_file >> theCurrentTemp.entx[k][i].xavg[j] >> theCurrentTemp.entx[k][i].xrms[j] >> theCurrentTemp.entx[k][i].xgx0[j] >> theCurrentTemp.entx[k][i].xgsig[j];
00417                                         
00418                                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 26, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00419                                 }
00420                                 
00421                                 for (j=0; j<4; ++j) {
00422                                         
00423                                         in_file >> theCurrentTemp.entx[k][i].xflpar[j][0] >> theCurrentTemp.entx[k][i].xflpar[j][1] >> theCurrentTemp.entx[k][i].xflpar[j][2] 
00424                                         >> theCurrentTemp.entx[k][i].xflpar[j][3] >> theCurrentTemp.entx[k][i].xflpar[j][4] >> theCurrentTemp.entx[k][i].xflpar[j][5];
00425                                         
00426                                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 27, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00427                                 }
00428                                 
00429                                 for (j=0; j<4; ++j) {
00430                                         
00431                                         in_file >> theCurrentTemp.entx[k][i].chi2yavg[j] >> theCurrentTemp.entx[k][i].chi2ymin[j] >> theCurrentTemp.entx[k][i].chi2xavg[j] >> theCurrentTemp.entx[k][i].chi2xmin[j];
00432                                         
00433                                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 28, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00434                                 }
00435                                 
00436                                 for (j=0; j<4; ++j) {
00437                                         
00438                                         in_file >> theCurrentTemp.entx[k][i].yavgc2m[j] >> theCurrentTemp.entx[k][i].yrmsc2m[j] >> theCurrentTemp.entx[k][i].chi2yavgc2m[j] >> theCurrentTemp.entx[k][i].chi2yminc2m[j];
00439                                         
00440                                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 29, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00441                                 }
00442                                 
00443                                 for (j=0; j<4; ++j) {
00444                                         
00445                                         in_file >> theCurrentTemp.entx[k][i].xavgc2m[j] >> theCurrentTemp.entx[k][i].xrmsc2m[j] >> theCurrentTemp.entx[k][i].chi2xavgc2m[j] >> theCurrentTemp.entx[k][i].chi2xminc2m[j];
00446                                         
00447                                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00448                                 }
00449                                 
00450                                 for (j=0; j<4; ++j) {
00451                                         
00452                                         in_file >> theCurrentTemp.entx[k][i].yavggen[j] >> theCurrentTemp.entx[k][i].yrmsgen[j] >> theCurrentTemp.entx[k][i].ygx0gen[j] >> theCurrentTemp.entx[k][i].ygsiggen[j];
00453                                         
00454                                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30a, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00455                                 }
00456                                 
00457                                 for (j=0; j<4; ++j) {
00458                                         
00459                                         in_file >> theCurrentTemp.entx[k][i].xavggen[j] >> theCurrentTemp.entx[k][i].xrmsgen[j] >> theCurrentTemp.entx[k][i].xgx0gen[j] >> theCurrentTemp.entx[k][i].xgsiggen[j];
00460                                         
00461                                         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30b, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00462                                 }
00463                                 
00464                                 in_file >> theCurrentTemp.entx[k][i].chi2yavgone >> theCurrentTemp.entx[k][i].chi2yminone >> theCurrentTemp.entx[k][i].chi2xavgone >> theCurrentTemp.entx[k][i].chi2xminone >> theCurrentTemp.entx[k][i].qmin2
00465                                 >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav >> theCurrentTemp.entx[k][i].qavg_spare >> theCurrentTemp.entx[k][i].spare[0];
00466                                 
00467                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 31, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00468                                 
00469                                 in_file >> theCurrentTemp.entx[k][i].mpvvav2 >> theCurrentTemp.entx[k][i].sigmavav2 >> theCurrentTemp.entx[k][i].kappavav2 >> theCurrentTemp.entx[k][i].qbfrac[0] >> theCurrentTemp.entx[k][i].qbfrac[1]
00470                                 >> theCurrentTemp.entx[k][i].qbfrac[2] >> theCurrentTemp.entx[k][i].fracyone >> theCurrentTemp.entx[k][i].fracxone >> theCurrentTemp.entx[k][i].fracytwo >> theCurrentTemp.entx[k][i].fracxtwo;
00471                                 //              theCurrentTemp.entx[k][i].qbfrac[3] = 1. - theCurrentTemp.entx[k][i].qbfrac[0] - theCurrentTemp.entx[k][i].qbfrac[1] - theCurrentTemp.entx[k][i].qbfrac[2];
00472                                 
00473                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 32, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00474                                 
00475                         }
00476                 }       
00477                 
00478                 
00479                 in_file.close();
00480                 
00481                 // Add this template to the store
00482                 
00483                 thePixelTemp_.push_back(theCurrentTemp);
00484                 
00485                 return true;
00486                 
00487         } else {
00488                 
00489                 // If file didn't open, report this
00490                 
00491                 LOGERROR("SiPixelTemplate") << "Error opening File" << tempfile << ENDL;
00492                 return false;
00493                 
00494         }
00495         
00496 } // TempInit 
00497 
00498 
00499 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00500 
00501 //**************************************************************** 
00505 //**************************************************************** 
00506 bool SiPixelTemplate::pushfile(const SiPixelTemplateDBObject& dbobject)
00507 {
00508         // Add template stored in external dbobject to theTemplateStore
00509     
00510         // Local variables 
00511         int i, j, k, l;
00512         float qavg_avg;
00513         //      const char *tempfile;
00514         const int code_version={16};
00515         
00516         // We must create a new object because dbobject must be a const and our stream must not be
00517         SiPixelTemplateDBObject db = dbobject;
00518         
00519         // Create a local template storage entry
00520         SiPixelTemplateStore theCurrentTemp;
00521         
00522         // Fill the template storage for each template calibration stored in the db
00523         for(int m=0; m<db.numOfTempl(); ++m)
00524         {
00525                 
00526                 // Read-in a header string first and print it    
00527                 
00528                 SiPixelTemplateDBObject::char2float temp;
00529                 for (i=0; i<20; ++i) {
00530                         temp.f = db.sVector()[db.index()];
00531                         theCurrentTemp.head.title[4*i] = temp.c[0];
00532                         theCurrentTemp.head.title[4*i+1] = temp.c[1];
00533                         theCurrentTemp.head.title[4*i+2] = temp.c[2];
00534                         theCurrentTemp.head.title[4*i+3] = temp.c[3];
00535                         db.incrementIndex(1);
00536                 }
00537                 theCurrentTemp.head.title[79] = '\0';
00538                 LOGINFO("SiPixelTemplate") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
00539                 
00540                 // next, the header information     
00541                 
00542                 db >> theCurrentTemp.head.ID  >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx
00543                 >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
00544                 >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >> theCurrentTemp.head.zsize;
00545                 
00546                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file, no template load" << ENDL; return false;}
00547                 
00548                 LOGINFO("SiPixelTemplate") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version " << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield 
00549                 << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
00550                 << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
00551                 << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
00552                 << ", 1/2 threshold " << theCurrentTemp.head.s50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", x Lorentz width " << theCurrentTemp.head.lorxwidth    
00553                 << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
00554                 
00555                 if(theCurrentTemp.head.templ_version < code_version) {LOGERROR("SiPixelTemplate") << "code expects version " << code_version << ", no template load" << ENDL; return false;}
00556                 
00557                 
00558 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST 
00559                 
00560 // next, layout the 1-d/2-d structures needed to store template
00561                 
00562                 theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
00563                 
00564                 theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
00565                 
00566 #endif
00567                                 
00568 // next, loop over all barrel y-angle entries   
00569                 
00570                 for (i=0; i < theCurrentTemp.head.NTy; ++i) {     
00571                         
00572                         db >> theCurrentTemp.enty[i].runnum >> theCurrentTemp.enty[i].costrk[0] 
00573                         >> theCurrentTemp.enty[i].costrk[1] >> theCurrentTemp.enty[i].costrk[2]; 
00574                         
00575                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 1, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00576                         
00577 // Calculate the alpha, beta, and cot(beta) for this entry 
00578                         
00579                         theCurrentTemp.enty[i].alpha = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[0]));
00580                         
00581                         theCurrentTemp.enty[i].cotalpha = theCurrentTemp.enty[i].costrk[0]/theCurrentTemp.enty[i].costrk[2];
00582                         
00583                         theCurrentTemp.enty[i].beta = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[1]));
00584                         
00585                         theCurrentTemp.enty[i].cotbeta = theCurrentTemp.enty[i].costrk[1]/theCurrentTemp.enty[i].costrk[2];
00586                         
00587                         db >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].pixmax >> theCurrentTemp.enty[i].symax >> theCurrentTemp.enty[i].dyone
00588                         >> theCurrentTemp.enty[i].syone >> theCurrentTemp.enty[i].sxmax >> theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone;
00589                         
00590                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 2, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00591                         
00592                         db >> theCurrentTemp.enty[i].dytwo >> theCurrentTemp.enty[i].sytwo >> theCurrentTemp.enty[i].dxtwo 
00593                         >> theCurrentTemp.enty[i].sxtwo >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].clsleny >> theCurrentTemp.enty[i].clslenx;
00594                         //                           >> theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav;
00595                         
00596                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 3, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00597                         
00598                         for (j=0; j<2; ++j) {
00599                                 
00600                                 db >> theCurrentTemp.enty[i].ypar[j][0] >> theCurrentTemp.enty[i].ypar[j][1] 
00601                                 >> theCurrentTemp.enty[i].ypar[j][2] >> theCurrentTemp.enty[i].ypar[j][3] >> theCurrentTemp.enty[i].ypar[j][4];
00602                                 
00603                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 4, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00604                                 
00605                         }
00606                         
00607                         for (j=0; j<9; ++j) {
00608                                 
00609                                 for (k=0; k<TYSIZE; ++k) {db >> theCurrentTemp.enty[i].ytemp[j][k];}
00610                                 
00611                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 5, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00612                         }
00613                         
00614                         for (j=0; j<2; ++j) {
00615                                 
00616                                 db >> theCurrentTemp.enty[i].xpar[j][0] >> theCurrentTemp.enty[i].xpar[j][1] 
00617                                 >> theCurrentTemp.enty[i].xpar[j][2] >> theCurrentTemp.enty[i].xpar[j][3] >> theCurrentTemp.enty[i].xpar[j][4];
00618                                 
00619                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 6, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00620                                 
00621                         }
00622                         
00623                         qavg_avg = 0.f;
00624                         for (j=0; j<9; ++j) {
00625                                 
00626                                 for (k=0; k<TXSIZE; ++k) {db >> theCurrentTemp.enty[i].xtemp[j][k]; qavg_avg += theCurrentTemp.enty[i].xtemp[j][k];} 
00627                                 
00628                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 7, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00629                         }
00630                         theCurrentTemp.enty[i].qavg_avg = qavg_avg/9.;
00631                         
00632                         for (j=0; j<4; ++j) {
00633                                 
00634                                 db >> theCurrentTemp.enty[i].yavg[j] >> theCurrentTemp.enty[i].yrms[j] >> theCurrentTemp.enty[i].ygx0[j] >> theCurrentTemp.enty[i].ygsig[j];
00635                                 
00636                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 8, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00637                         }
00638                         
00639                         for (j=0; j<4; ++j) {
00640                                 
00641                                 db >> theCurrentTemp.enty[i].yflpar[j][0] >> theCurrentTemp.enty[i].yflpar[j][1] >> theCurrentTemp.enty[i].yflpar[j][2] 
00642                                 >> theCurrentTemp.enty[i].yflpar[j][3] >> theCurrentTemp.enty[i].yflpar[j][4] >> theCurrentTemp.enty[i].yflpar[j][5];
00643                                 
00644                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 9, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00645                         }
00646                         
00647                         for (j=0; j<4; ++j) {
00648                                 
00649                                 db >> theCurrentTemp.enty[i].xavg[j] >> theCurrentTemp.enty[i].xrms[j] >> theCurrentTemp.enty[i].xgx0[j] >> theCurrentTemp.enty[i].xgsig[j];
00650                                 
00651                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 10, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00652                         }
00653                         
00654                         for (j=0; j<4; ++j) {
00655                                 
00656                                 db >> theCurrentTemp.enty[i].xflpar[j][0] >> theCurrentTemp.enty[i].xflpar[j][1] >> theCurrentTemp.enty[i].xflpar[j][2] 
00657                                 >> theCurrentTemp.enty[i].xflpar[j][3] >> theCurrentTemp.enty[i].xflpar[j][4] >> theCurrentTemp.enty[i].xflpar[j][5];
00658                                 
00659                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 11, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00660                         }
00661                         
00662                         for (j=0; j<4; ++j) {
00663                                 
00664                                 db >> theCurrentTemp.enty[i].chi2yavg[j] >> theCurrentTemp.enty[i].chi2ymin[j] >> theCurrentTemp.enty[i].chi2xavg[j] >> theCurrentTemp.enty[i].chi2xmin[j];
00665                                 
00666                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 12, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00667                         }
00668                         
00669                         for (j=0; j<4; ++j) {
00670                                 
00671                                 db >> theCurrentTemp.enty[i].yavgc2m[j] >> theCurrentTemp.enty[i].yrmsc2m[j] >> theCurrentTemp.enty[i].chi2yavgc2m[j] >> theCurrentTemp.enty[i].chi2yminc2m[j];
00672                                 
00673                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 13, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00674                         }
00675                         
00676                         for (j=0; j<4; ++j) {
00677                                 
00678                                 db >> theCurrentTemp.enty[i].xavgc2m[j] >> theCurrentTemp.enty[i].xrmsc2m[j] >> theCurrentTemp.enty[i].chi2xavgc2m[j] >> theCurrentTemp.enty[i].chi2xminc2m[j];
00679                                 
00680                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00681                         } 
00682                         
00683                         for (j=0; j<4; ++j) {
00684                                 
00685                                 db >> theCurrentTemp.enty[i].yavggen[j] >> theCurrentTemp.enty[i].yrmsgen[j] >> theCurrentTemp.enty[i].ygx0gen[j] >> theCurrentTemp.enty[i].ygsiggen[j];
00686                                 
00687                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14a, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00688                         }
00689                         
00690                         for (j=0; j<4; ++j) {
00691                                 
00692                                 db >> theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j] >> theCurrentTemp.enty[i].xgx0gen[j] >> theCurrentTemp.enty[i].xgsiggen[j];
00693                                 
00694                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14b, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00695                         } 
00696                         
00697                         
00698                         db >> theCurrentTemp.enty[i].chi2yavgone >> theCurrentTemp.enty[i].chi2yminone >> theCurrentTemp.enty[i].chi2xavgone >> theCurrentTemp.enty[i].chi2xminone >> theCurrentTemp.enty[i].qmin2
00699                         >> theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav >> theCurrentTemp.enty[i].qavg_spare >> theCurrentTemp.enty[i].spare[0];
00700                         
00701                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 15, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00702                         
00703                         db >> theCurrentTemp.enty[i].mpvvav2 >> theCurrentTemp.enty[i].sigmavav2 >> theCurrentTemp.enty[i].kappavav2 >> theCurrentTemp.enty[i].qbfrac[0] >> theCurrentTemp.enty[i].qbfrac[1]
00704                         >> theCurrentTemp.enty[i].qbfrac[2] >> theCurrentTemp.enty[i].fracyone >> theCurrentTemp.enty[i].fracxone >> theCurrentTemp.enty[i].fracytwo >> theCurrentTemp.enty[i].fracxtwo;
00705                         //                      theCurrentTemp.enty[i].qbfrac[3] = 1. - theCurrentTemp.enty[i].qbfrac[0] - theCurrentTemp.enty[i].qbfrac[1] - theCurrentTemp.enty[i].qbfrac[2];
00706                         
00707                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 16, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00708                         
00709                 }
00710                 
00711                 // next, loop over all barrel x-angle entries   
00712                 
00713                 for (k=0; k < theCurrentTemp.head.NTyx; ++k) { 
00714                         
00715                         for (i=0; i < theCurrentTemp.head.NTxx; ++i) { 
00716                                 
00717                                 db >> theCurrentTemp.entx[k][i].runnum >> theCurrentTemp.entx[k][i].costrk[0] 
00718                                 >> theCurrentTemp.entx[k][i].costrk[1] >> theCurrentTemp.entx[k][i].costrk[2]; 
00719                                 
00720                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 17, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00721                                 
00722                                 // Calculate the alpha, beta, and cot(beta) for this entry 
00723                                 
00724                                 theCurrentTemp.entx[k][i].alpha = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[0]));
00725                                 
00726                                 theCurrentTemp.entx[k][i].cotalpha = theCurrentTemp.entx[k][i].costrk[0]/theCurrentTemp.entx[k][i].costrk[2];
00727                                 
00728                                 theCurrentTemp.entx[k][i].beta = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[1]));
00729                                 
00730                                 theCurrentTemp.entx[k][i].cotbeta = theCurrentTemp.entx[k][i].costrk[1]/theCurrentTemp.entx[k][i].costrk[2];
00731                                 
00732                                 db >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].pixmax >> theCurrentTemp.entx[k][i].symax >> theCurrentTemp.entx[k][i].dyone
00733                                 >> theCurrentTemp.entx[k][i].syone >> theCurrentTemp.entx[k][i].sxmax >> theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone;
00734                                 
00735                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 18, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00736                                 
00737                                 db >> theCurrentTemp.entx[k][i].dytwo >> theCurrentTemp.entx[k][i].sytwo >> theCurrentTemp.entx[k][i].dxtwo 
00738                                 >> theCurrentTemp.entx[k][i].sxtwo >> theCurrentTemp.entx[k][i].qmin >> theCurrentTemp.entx[k][i].clsleny >> theCurrentTemp.entx[k][i].clslenx;
00739                                 //                     >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav;
00740                                 
00741                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 19, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00742                                 
00743                                 for (j=0; j<2; ++j) {
00744                                         
00745                                         db >> theCurrentTemp.entx[k][i].ypar[j][0] >> theCurrentTemp.entx[k][i].ypar[j][1] 
00746                                         >> theCurrentTemp.entx[k][i].ypar[j][2] >> theCurrentTemp.entx[k][i].ypar[j][3] >> theCurrentTemp.entx[k][i].ypar[j][4];
00747                                         
00748                                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 20, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00749                                 }
00750                                 
00751                                 for (j=0; j<9; ++j) {
00752                                         
00753                                         for (l=0; l<TYSIZE; ++l) {db >> theCurrentTemp.entx[k][i].ytemp[j][l];} 
00754                                         
00755                                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 21, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00756                                 }
00757                                 
00758                                 for (j=0; j<2; ++j) {
00759                                         
00760                                         db >> theCurrentTemp.entx[k][i].xpar[j][0] >> theCurrentTemp.entx[k][i].xpar[j][1] 
00761                                         >> theCurrentTemp.entx[k][i].xpar[j][2] >> theCurrentTemp.entx[k][i].xpar[j][3] >> theCurrentTemp.entx[k][i].xpar[j][4];
00762                                         
00763                                         
00764                                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 22, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00765                                 }
00766                                 
00767                                 qavg_avg = 0.f;
00768                                 for (j=0; j<9; ++j) {
00769                                         
00770                                         for (l=0; l<TXSIZE; ++l) {db >> theCurrentTemp.entx[k][i].xtemp[j][l]; qavg_avg += theCurrentTemp.entx[k][i].xtemp[j][l];} 
00771                                         
00772                                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 23, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00773                                 }
00774                                 theCurrentTemp.entx[k][i].qavg_avg = qavg_avg/9.;
00775                                 
00776                                 for (j=0; j<4; ++j) {
00777                                         
00778                                         db >> theCurrentTemp.entx[k][i].yavg[j] >> theCurrentTemp.entx[k][i].yrms[j] >> theCurrentTemp.entx[k][i].ygx0[j] >> theCurrentTemp.entx[k][i].ygsig[j];
00779                                         
00780                                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 24, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00781                                 }
00782                                 
00783                                 for (j=0; j<4; ++j) {
00784                                         
00785                                         db >> theCurrentTemp.entx[k][i].yflpar[j][0] >> theCurrentTemp.entx[k][i].yflpar[j][1] >> theCurrentTemp.entx[k][i].yflpar[j][2] 
00786                                         >> theCurrentTemp.entx[k][i].yflpar[j][3] >> theCurrentTemp.entx[k][i].yflpar[j][4] >> theCurrentTemp.entx[k][i].yflpar[j][5];
00787                                         
00788                                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 25, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00789                                 }
00790                                 
00791                                 for (j=0; j<4; ++j) {
00792                                         
00793                                         db >> theCurrentTemp.entx[k][i].xavg[j] >> theCurrentTemp.entx[k][i].xrms[j] >> theCurrentTemp.entx[k][i].xgx0[j] >> theCurrentTemp.entx[k][i].xgsig[j];
00794                                         
00795                                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 26, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00796                                 }
00797                                 
00798                                 for (j=0; j<4; ++j) {
00799                                         
00800                                         db >> theCurrentTemp.entx[k][i].xflpar[j][0] >> theCurrentTemp.entx[k][i].xflpar[j][1] >> theCurrentTemp.entx[k][i].xflpar[j][2] 
00801                                         >> theCurrentTemp.entx[k][i].xflpar[j][3] >> theCurrentTemp.entx[k][i].xflpar[j][4] >> theCurrentTemp.entx[k][i].xflpar[j][5];
00802                                         
00803                                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 27, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00804                                 }
00805                                 
00806                                 for (j=0; j<4; ++j) {
00807                                         
00808                                         db >> theCurrentTemp.entx[k][i].chi2yavg[j] >> theCurrentTemp.entx[k][i].chi2ymin[j] >> theCurrentTemp.entx[k][i].chi2xavg[j] >> theCurrentTemp.entx[k][i].chi2xmin[j];
00809                                         
00810                                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 28, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00811                                 }
00812                                 
00813                                 for (j=0; j<4; ++j) {
00814                                         
00815                                         db >> theCurrentTemp.entx[k][i].yavgc2m[j] >> theCurrentTemp.entx[k][i].yrmsc2m[j] >> theCurrentTemp.entx[k][i].chi2yavgc2m[j] >> theCurrentTemp.entx[k][i].chi2yminc2m[j];
00816                                         
00817                                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 29, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00818                                 }
00819                                 
00820                                 for (j=0; j<4; ++j) {
00821                                         
00822                                         db >> theCurrentTemp.entx[k][i].xavgc2m[j] >> theCurrentTemp.entx[k][i].xrmsc2m[j] >> theCurrentTemp.entx[k][i].chi2xavgc2m[j] >> theCurrentTemp.entx[k][i].chi2xminc2m[j];
00823                                         
00824                                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00825                                 }
00826                                 
00827                                 for (j=0; j<4; ++j) {
00828                                         
00829                                         db >> theCurrentTemp.entx[k][i].yavggen[j] >> theCurrentTemp.entx[k][i].yrmsgen[j] >> theCurrentTemp.entx[k][i].ygx0gen[j] >> theCurrentTemp.entx[k][i].ygsiggen[j];
00830                                         
00831                                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30a, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00832                                 }
00833                                 
00834                                 for (j=0; j<4; ++j) {
00835                                         
00836                                         db >> theCurrentTemp.entx[k][i].xavggen[j] >> theCurrentTemp.entx[k][i].xrmsgen[j] >> theCurrentTemp.entx[k][i].xgx0gen[j] >> theCurrentTemp.entx[k][i].xgsiggen[j];
00837                                         
00838                                         if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30b, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00839                                 }
00840                                 
00841                                 
00842                                 db >> theCurrentTemp.entx[k][i].chi2yavgone >> theCurrentTemp.entx[k][i].chi2yminone >> theCurrentTemp.entx[k][i].chi2xavgone >> theCurrentTemp.entx[k][i].chi2xminone >> theCurrentTemp.entx[k][i].qmin2
00843                                 >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav >> theCurrentTemp.entx[k][i].qavg_spare >> theCurrentTemp.entx[k][i].spare[0];
00844                                 
00845                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 31, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00846                                 
00847                                 db >> theCurrentTemp.entx[k][i].mpvvav2 >> theCurrentTemp.entx[k][i].sigmavav2 >> theCurrentTemp.entx[k][i].kappavav2 >> theCurrentTemp.entx[k][i].qbfrac[0] >> theCurrentTemp.entx[k][i].qbfrac[1]
00848                                 >> theCurrentTemp.entx[k][i].qbfrac[2] >> theCurrentTemp.entx[k][i].fracyone >> theCurrentTemp.entx[k][i].fracxone >> theCurrentTemp.entx[k][i].fracytwo >> theCurrentTemp.entx[k][i].fracxtwo;
00849                                 //                              theCurrentTemp.entx[k][i].qbfrac[3] = 1. - theCurrentTemp.entx[k][i].qbfrac[0] - theCurrentTemp.entx[k][i].qbfrac[1] - theCurrentTemp.entx[k][i].qbfrac[2];
00850                                 
00851                                 if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 32, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00852                                 
00853                         }
00854                 }       
00855                 
00856                 
00857                 // Add this template to the store
00858                 
00859                 thePixelTemp_.push_back(theCurrentTemp);
00860                 
00861         }
00862         return true;
00863         
00864 } // TempInit 
00865 
00866 #endif
00867 
00868 
00869 // ************************************************************************************************************ 
00876 // ************************************************************************************************************ 
00877 bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta, float locBz){
00878     // Interpolate for a new set of track angles 
00879     
00880     // Local variables 
00881     int i, j;
00882         int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, imidy, imaxx;
00883         float yratio, yxratio, xxratio, sxmax, qcorrect, qxtempcor, symax, chi2xavgone, chi2xminone, cotb, cotalpha0, cotbeta0;
00884         bool flip_y;
00885 //      std::vector <float> xrms(4), xgsig(4), xrmsc2m(4);
00886         float chi2xavg[4], chi2xmin[4], chi2xavgc2m[4], chi2xminc2m[4];
00887 
00888 
00889 // Check to see if interpolation is valid     
00890 
00891 if(id != id_current_ || cotalpha != cota_current_ || cotbeta != cotb_current_) {
00892 
00893         cota_current_ = cotalpha; cotb_current_ = cotbeta; success_ = true;
00894         
00895         if(id != id_current_) {
00896 
00897 // Find the index corresponding to id
00898 
00899        index_id_ = -1;
00900        for(i=0; i<(int)thePixelTemp_.size(); ++i) {
00901         
00902               if(id == thePixelTemp_[i].head.ID) {
00903            
00904                  index_id_ = i;
00905                       id_current_ = id;
00906                                 
00907 // Copy the charge scaling factor to the private variable     
00908                                 
00909                                 qscale_ = thePixelTemp_[index_id_].head.qscale;
00910                                 
00911 // Copy the pseudopixel signal size to the private variable     
00912                                 
00913                                 s50_ = thePixelTemp_[index_id_].head.s50;
00914                         
00915 // Pixel sizes to the private variables     
00916                                 
00917                                 xsize_ = thePixelTemp_[index_id_].head.xsize;
00918                                 ysize_ = thePixelTemp_[index_id_].head.ysize;
00919                                 zsize_ = thePixelTemp_[index_id_].head.zsize;
00920                                 
00921                                 break;
00922           }
00923             }
00924      }
00925          
00926 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00927     if(index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
00928                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::interpolate can't find needed template ID = " << id << std::endl;
00929         }
00930 #else
00931         assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
00932 #endif
00933          
00934 // Interpolate the absolute value of cot(beta)     
00935     
00936         abs_cotb_ = std::abs(cotbeta);
00937         
00938 //      qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)       
00939 
00940         cotalpha0 =  thePixelTemp_[index_id_].enty[0].cotalpha;
00941         qcorrect=std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta*cotbeta+cotalpha0*cotalpha0));
00942         
00943 // for some cosmics, the ususal gymnastics are incorrect   
00944         if(thePixelTemp_[index_id_].head.Dtype == 0) {
00945                 cotb = abs_cotb_;
00946                 flip_y = false;
00947                 if(cotbeta < 0.f) {flip_y = true;}
00948         } else {
00949             if(locBz < 0.f) {
00950                         cotb = cotbeta;
00951                         flip_y = false;
00952                 } else {
00953                         cotb = -cotbeta;
00954                         flip_y = true;
00955                 }       
00956         }
00957                 
00958         Ny = thePixelTemp_[index_id_].head.NTy;
00959         Nyx = thePixelTemp_[index_id_].head.NTyx;
00960         Nxx = thePixelTemp_[index_id_].head.NTxx;
00961                 
00962 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00963         if(Ny < 2 || Nyx < 1 || Nxx < 2) {
00964                 throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
00965         }
00966 #else
00967         assert(Ny > 1 && Nyx > 0 && Nxx > 1);
00968 #endif
00969         imaxx = Nyx - 1;
00970         imidy = Nxx/2;
00971         
00972 // next, loop over all y-angle entries   
00973 
00974         ilow = 0;
00975         yratio = 0.f;
00976 
00977         if(cotb >= thePixelTemp_[index_id_].enty[Ny-1].cotbeta) {
00978         
00979                 ilow = Ny-2;
00980                 yratio = 1.;
00981                 success_ = false;
00982                 
00983         } else {
00984            
00985                 if(cotb >= thePixelTemp_[index_id_].enty[0].cotbeta) {
00986 
00987                         for (i=0; i<Ny-1; ++i) { 
00988     
00989                         if( thePixelTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index_id_].enty[i+1].cotbeta) {
00990                   
00991                                 ilow = i;
00992                                 yratio = (cotb - thePixelTemp_[index_id_].enty[i].cotbeta)/(thePixelTemp_[index_id_].enty[i+1].cotbeta - thePixelTemp_[index_id_].enty[i].cotbeta);
00993                                 break;                   
00994                         }
00995                 }
00996         } else { success_ = false; }
00997   }
00998         
00999         ihigh=ilow + 1;
01000                           
01001 // Interpolate/store all y-related quantities (flip displacements when flip_y)
01002 
01003         yratio_ = yratio;
01004         qavg_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].qavg + yratio*thePixelTemp_[index_id_].enty[ihigh].qavg;
01005         qavg_ *= qcorrect;
01006         symax = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].symax + yratio*thePixelTemp_[index_id_].enty[ihigh].symax;
01007         syparmax_ = symax;
01008         sxmax = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].sxmax + yratio*thePixelTemp_[index_id_].enty[ihigh].sxmax;
01009         dyone_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].dyone + yratio*thePixelTemp_[index_id_].enty[ihigh].dyone;
01010         if(flip_y) {dyone_ = -dyone_;}
01011         syone_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].syone + yratio*thePixelTemp_[index_id_].enty[ihigh].syone;
01012         dytwo_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].dytwo + yratio*thePixelTemp_[index_id_].enty[ihigh].dytwo;
01013         if(flip_y) {dytwo_ = -dytwo_;}
01014         sytwo_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].sytwo + yratio*thePixelTemp_[index_id_].enty[ihigh].sytwo;
01015         qmin_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].qmin + yratio*thePixelTemp_[index_id_].enty[ihigh].qmin;
01016         qmin_ *= qcorrect;
01017         qmin2_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].qmin2 + yratio*thePixelTemp_[index_id_].enty[ihigh].qmin2;
01018         qmin2_ *= qcorrect;
01019         mpvvav_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].mpvvav + yratio*thePixelTemp_[index_id_].enty[ihigh].mpvvav;
01020         mpvvav_ *= qcorrect;
01021         sigmavav_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].sigmavav + yratio*thePixelTemp_[index_id_].enty[ihigh].sigmavav;
01022         kappavav_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].kappavav + yratio*thePixelTemp_[index_id_].enty[ihigh].kappavav;
01023         mpvvav2_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].mpvvav2 + yratio*thePixelTemp_[index_id_].enty[ihigh].mpvvav2;
01024         mpvvav2_ *= qcorrect;
01025         sigmavav2_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].sigmavav2 + yratio*thePixelTemp_[index_id_].enty[ihigh].sigmavav2;
01026         kappavav2_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].kappavav2 + yratio*thePixelTemp_[index_id_].enty[ihigh].kappavav2;
01027         clsleny_ = fminf(thePixelTemp_[index_id_].enty[ilow].clsleny, thePixelTemp_[index_id_].enty[ihigh].clsleny);
01028         qavg_avg_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].qavg_avg + yratio*thePixelTemp_[index_id_].enty[ihigh].qavg_avg;
01029         qavg_avg_ *= qcorrect;
01030         for(i=0; i<2 ; ++i) {
01031                 for(j=0; j<5 ; ++j) {
01032 // Charge loss switches sides when cot(beta) changes sign
01033                         if(flip_y) {
01034                     yparl_[1-i][j] = thePixelTemp_[index_id_].enty[ilow].ypar[i][j];
01035                     yparh_[1-i][j] = thePixelTemp_[index_id_].enty[ihigh].ypar[i][j];
01036                         } else {
01037                     yparl_[i][j] = thePixelTemp_[index_id_].enty[ilow].ypar[i][j];
01038                     yparh_[i][j] = thePixelTemp_[index_id_].enty[ihigh].ypar[i][j];
01039                         }
01040                         xparly0_[i][j] = thePixelTemp_[index_id_].enty[ilow].xpar[i][j];
01041                         xparhy0_[i][j] = thePixelTemp_[index_id_].enty[ihigh].xpar[i][j];
01042                 }
01043         }
01044         for(i=0; i<4; ++i) {
01045                 yavg_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].yavg[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].yavg[i];
01046                 if(flip_y) {yavg_[i] = -yavg_[i];}
01047                 yrms_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].yrms[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].yrms[i];
01048 //            ygx0_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].ygx0[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].ygx0[i];
01049 //            if(flip_y) {ygx0_[i] = -ygx0_[i];}
01050 //            ygsig_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].ygsig[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].ygsig[i];
01051 //            xrms[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].xrms[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].xrms[i];
01052 //            xgsig[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].xgsig[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].xgsig[i];
01053                 chi2yavg_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2yavg[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2yavg[i];
01054                 chi2ymin_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2ymin[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2ymin[i];
01055                 chi2xavg[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2xavg[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2xavg[i];
01056                 chi2xmin[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2xmin[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2xmin[i];
01057                 yavgc2m_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].yavgc2m[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].yavgc2m[i];
01058                 if(flip_y) {yavgc2m_[i] = -yavgc2m_[i];}
01059                 yrmsc2m_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].yrmsc2m[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].yrmsc2m[i];
01060                 chi2yavgc2m_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2yavgc2m[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2yavgc2m[i];
01061 //            if(flip_y) {chi2yavgc2m_[i] = -chi2yavgc2m_[i];}
01062                 chi2yminc2m_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2yminc2m[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2yminc2m[i];
01063 //            xrmsc2m[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].xrmsc2m[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].xrmsc2m[i];
01064                 chi2xavgc2m[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2xavgc2m[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2xavgc2m[i];
01065                 chi2xminc2m[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2xminc2m[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2xminc2m[i];
01066                 for(j=0; j<6 ; ++j) {
01067                         yflparl_[i][j] = thePixelTemp_[index_id_].enty[ilow].yflpar[i][j];
01068                         yflparh_[i][j] = thePixelTemp_[index_id_].enty[ihigh].yflpar[i][j];
01069                          
01070 // Since Q_fl is odd under cotbeta, it flips qutomatically, change only even terms
01071 
01072                         if(flip_y && (j == 0 || j == 2 || j == 4)) {
01073                            yflparl_[i][j] = - yflparl_[i][j];
01074                            yflparh_[i][j] = - yflparh_[i][j];
01075                         }
01076                 }
01077         }
01078            
01080 
01081         chi2yavgone_=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2yavgone + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2yavgone;
01082         chi2yminone_=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2yminone + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2yminone;
01083         chi2xavgone=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2xavgone + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2xavgone;
01084         chi2xminone=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2xminone + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2xminone;
01085                 //       for(i=0; i<10; ++i) {
01086 //                  pyspare[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].yspare[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].yspare[i];
01087 //       }
01088                           
01089 // Interpolate and build the y-template 
01090         
01091         for(i=0; i<9; ++i) {
01092                 ytemp_[i][0] = 0.f;
01093                 ytemp_[i][1] = 0.f;
01094                 ytemp_[i][BYM2] = 0.f;
01095                 ytemp_[i][BYM1] = 0.f;
01096                 for(j=0; j<TYSIZE; ++j) {
01097                   
01098 // Flip the basic y-template when the cotbeta is negative
01099 
01100                         if(flip_y) {
01101                            ytemp_[8-i][BYM3-j]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].ytemp[i][j] + yratio*thePixelTemp_[index_id_].enty[ihigh].ytemp[i][j];
01102                         } else {
01103                            ytemp_[i][j+2]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].ytemp[i][j] + yratio*thePixelTemp_[index_id_].enty[ihigh].ytemp[i][j];
01104                         }
01105                 }
01106         }
01107         
01108 // next, loop over all x-angle entries, first, find relevant y-slices   
01109         
01110         iylow = 0;
01111         yxratio = 0.f;
01112 
01113         if(abs_cotb_ >= thePixelTemp_[index_id_].entx[Nyx-1][0].cotbeta) {
01114         
01115                 iylow = Nyx-2;
01116                 yxratio = 1.f;
01117                 
01118         } else if(abs_cotb_ >= thePixelTemp_[index_id_].entx[0][0].cotbeta) {
01119 
01120                 for (i=0; i<Nyx-1; ++i) { 
01121     
01122                         if( thePixelTemp_[index_id_].entx[i][0].cotbeta <= abs_cotb_ && abs_cotb_ < thePixelTemp_[index_id_].entx[i+1][0].cotbeta) {
01123                   
01124                            iylow = i;
01125                            yxratio = (abs_cotb_ - thePixelTemp_[index_id_].entx[i][0].cotbeta)/(thePixelTemp_[index_id_].entx[i+1][0].cotbeta - thePixelTemp_[index_id_].entx[i][0].cotbeta);
01126                            break;                        
01127                         }
01128                 }
01129         }
01130         
01131         iyhigh=iylow + 1;
01132 
01133         ilow = 0;
01134         xxratio = 0.f;
01135 
01136         if(cotalpha >= thePixelTemp_[index_id_].entx[0][Nxx-1].cotalpha) {
01137         
01138                 ilow = Nxx-2;
01139                 xxratio = 1.f;
01140                 success_ = false;
01141                 
01142         } else {
01143            
01144                 if(cotalpha >= thePixelTemp_[index_id_].entx[0][0].cotalpha) {
01145 
01146                         for (i=0; i<Nxx-1; ++i) { 
01147     
01148                            if( thePixelTemp_[index_id_].entx[0][i].cotalpha <= cotalpha && cotalpha < thePixelTemp_[index_id_].entx[0][i+1].cotalpha) {
01149                   
01150                                   ilow = i;
01151                                   xxratio = (cotalpha - thePixelTemp_[index_id_].entx[0][i].cotalpha)/(thePixelTemp_[index_id_].entx[0][i+1].cotalpha - thePixelTemp_[index_id_].entx[0][i].cotalpha);
01152                                   break;
01153                         }
01154                   }
01155                 } else { success_ = false; }
01156         }
01157         
01158         ihigh=ilow + 1;
01159                           
01160 // Interpolate/store all x-related quantities 
01161 
01162         yxratio_ = yxratio;
01163         xxratio_ = xxratio;             
01164                           
01165 // sxparmax defines the maximum charge for which the parameters xpar are defined (not rescaled by cotbeta) 
01166 
01167         sxparmax_ = (1.f - xxratio)*thePixelTemp_[index_id_].entx[imaxx][ilow].sxmax + xxratio*thePixelTemp_[index_id_].entx[imaxx][ihigh].sxmax;
01168         sxmax_ = sxparmax_;
01169         if(thePixelTemp_[index_id_].entx[imaxx][imidy].sxmax != 0.f) {sxmax_=sxmax_/thePixelTemp_[index_id_].entx[imaxx][imidy].sxmax*sxmax;}
01170         symax_ = (1.f - xxratio)*thePixelTemp_[index_id_].entx[imaxx][ilow].symax + xxratio*thePixelTemp_[index_id_].entx[imaxx][ihigh].symax;
01171         if(thePixelTemp_[index_id_].entx[imaxx][imidy].symax != 0.f) {symax_=symax_/thePixelTemp_[index_id_].entx[imaxx][imidy].symax*symax;}
01172         dxone_ = (1.f - xxratio)*thePixelTemp_[index_id_].entx[0][ilow].dxone + xxratio*thePixelTemp_[index_id_].entx[0][ihigh].dxone;
01173         sxone_ = (1.f - xxratio)*thePixelTemp_[index_id_].entx[0][ilow].sxone + xxratio*thePixelTemp_[index_id_].entx[0][ihigh].sxone;
01174         dxtwo_ = (1.f - xxratio)*thePixelTemp_[index_id_].entx[0][ilow].dxtwo + xxratio*thePixelTemp_[index_id_].entx[0][ihigh].dxtwo;
01175         sxtwo_ = (1.f - xxratio)*thePixelTemp_[index_id_].entx[0][ilow].sxtwo + xxratio*thePixelTemp_[index_id_].entx[0][ihigh].sxtwo;
01176         clslenx_ = fminf(thePixelTemp_[index_id_].entx[0][ilow].clslenx, thePixelTemp_[index_id_].entx[0][ihigh].clslenx);
01177         for(i=0; i<2 ; ++i) {
01178                 for(j=0; j<5 ; ++j) {
01179                  xpar0_[i][j] = thePixelTemp_[index_id_].entx[imaxx][imidy].xpar[i][j];
01180                  xparl_[i][j] = thePixelTemp_[index_id_].entx[imaxx][ilow].xpar[i][j];
01181                  xparh_[i][j] = thePixelTemp_[index_id_].entx[imaxx][ihigh].xpar[i][j];
01182                 }
01183         }
01184                           
01185 // pixmax is the maximum allowed pixel charge (used for truncation)
01186 
01187         pixmax_=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].pixmax + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].pixmax)
01188                         +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].pixmax + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].pixmax);
01189                           
01190         for(i=0; i<4; ++i) {
01191                 xavg_[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].xavg[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].xavg[i])
01192                                 +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xavg[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xavg[i]);
01193                   
01194                 xrms_[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].xrms[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].xrms[i])
01195                                 +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xrms[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xrms[i]);
01196                   
01197 //            xgx0_[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].xgx0[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].xgx0[i])
01198 //                        +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xgx0[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xgx0[i]);
01199                                                         
01200 //            xgsig_[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].xgsig[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].xgsig[i])
01201 //                        +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xgsig[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xgsig[i]);
01202                                   
01203                 xavgc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].xavgc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].xavgc2m[i])
01204                                 +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xavgc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xavgc2m[i]);
01205                   
01206                 xrmsc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].xrmsc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].xrmsc2m[i])
01207                                 +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xrmsc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xrmsc2m[i]);
01208                   
01209 //            chi2xavgc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].chi2xavgc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].chi2xavgc2m[i])
01210 //                        +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].chi2xavgc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].chi2xavgc2m[i]);
01211                                                         
01212 //            chi2xminc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].chi2xminc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].chi2xminc2m[i])
01213 //                        +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].chi2xminc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].chi2xminc2m[i]);
01214 //
01215 //  Try new interpolation scheme
01216 //                                                                                                                      
01217 //            chi2xavg_[i]=((1.f - xxratio)*thePixelTemp_[index_id_].entx[imaxx][ilow].chi2xavg[i] + xxratio*thePixelTemp_[index_id_].entx[imaxx][ihigh].chi2xavg[i]);
01218 //                if(thePixelTemp_[index_id_].entx[imaxx][imidy].chi2xavg[i] != 0.f) {chi2xavg_[i]=chi2xavg_[i]/thePixelTemp_[index_id_].entx[imaxx][imidy].chi2xavg[i]*chi2xavg[i];}
01219 //                                                      
01220 //            chi2xmin_[i]=((1.f - xxratio)*thePixelTemp_[index_id_].entx[imaxx][ilow].chi2xmin[i] + xxratio*thePixelTemp_[index_id_].entx[imaxx][ihigh].chi2xmin[i]);
01221 //                if(thePixelTemp_[index_id_].entx[imaxx][imidy].chi2xmin[i] != 0.f) {chi2xmin_[i]=chi2xmin_[i]/thePixelTemp_[index_id_].entx[imaxx][imidy].chi2xmin[i]*chi2xmin[i];}
01222 //                
01223                 chi2xavg_[i]=((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].chi2xavg[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].chi2xavg[i]);
01224                 if(thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xavg[i] != 0.f) {chi2xavg_[i]=chi2xavg_[i]/thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xavg[i]*chi2xavg[i];}
01225                                                         
01226                 chi2xmin_[i]=((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].chi2xmin[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].chi2xmin[i]);
01227                 if(thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xmin[i] != 0.f) {chi2xmin_[i]=chi2xmin_[i]/thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xmin[i]*chi2xmin[i];}
01228                 
01229                 chi2xavgc2m_[i]=((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].chi2xavgc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].chi2xavgc2m[i]);
01230                 if(thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xavgc2m[i] != 0.f) {chi2xavgc2m_[i]=chi2xavgc2m_[i]/thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xavgc2m[i]*chi2xavgc2m[i];}
01231                 
01232                 chi2xminc2m_[i]=((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].chi2xminc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].chi2xminc2m[i]);
01233                 if(thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xminc2m[i] != 0.f) {chi2xminc2m_[i]=chi2xminc2m_[i]/thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xminc2m[i]*chi2xminc2m[i];}
01234                 
01235                 for(j=0; j<6 ; ++j) {
01236                  xflparll_[i][j] = thePixelTemp_[index_id_].entx[iylow][ilow].xflpar[i][j];
01237                  xflparlh_[i][j] = thePixelTemp_[index_id_].entx[iylow][ihigh].xflpar[i][j];
01238                  xflparhl_[i][j] = thePixelTemp_[index_id_].entx[iyhigh][ilow].xflpar[i][j];
01239                  xflparhh_[i][j] = thePixelTemp_[index_id_].entx[iyhigh][ihigh].xflpar[i][j];
01240                 }
01241         }
01242            
01243 // Do the spares next
01244 
01245         chi2xavgone_=((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].chi2xavgone + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].chi2xavgone);
01246         if(thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xavgone != 0.f) {chi2xavgone_=chi2xavgone_/thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xavgone*chi2xavgone;}
01247                 
01248         chi2xminone_=((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].chi2xminone + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].chi2xminone);
01249         if(thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xminone != 0.f) {chi2xminone_=chi2xminone_/thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xminone*chi2xminone;}
01250                 //       for(i=0; i<10; ++i) {
01251 //            pxspare[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].xspare[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].xspare[i])
01252 //                        +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xspare[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xspare[i]);
01253 //       }
01254                           
01255 // Interpolate and build the x-template 
01256         
01257 //      qxtempcor corrects the total charge to the actual track angles (not actually needed for the template fits, but useful for Guofan)
01258         
01259         cotbeta0 =  thePixelTemp_[index_id_].entx[iyhigh][0].cotbeta;
01260         qxtempcor=std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta0*cotbeta0+cotalpha*cotalpha));
01261         
01262         for(i=0; i<9; ++i) {
01263                 xtemp_[i][0] = 0.f;
01264                 xtemp_[i][1] = 0.f;
01265                 xtemp_[i][BXM2] = 0.f;
01266                 xtemp_[i][BXM1] = 0.f;
01267                 for(j=0; j<TXSIZE; ++j) {
01268 //  Take next largest x-slice for the x-template (it reduces bias in the forward direction after irradiation)
01269 //                 xtemp_[i][j+2]=(1.f - xxratio)*thePixelTemp_[index_id_].entx[imaxx][ilow].xtemp[i][j] + xxratio*thePixelTemp_[index_id_].entx[imaxx][ihigh].xtemp[i][j];
01270 //                 xtemp_[i][j+2]=(1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xtemp[i][j] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xtemp[i][j];
01271          xtemp_[i][j+2]=qxtempcor*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xtemp[i][j] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xtemp[i][j]);
01272                 }
01273         }
01274         
01275         lorywidth_ = thePixelTemp_[index_id_].head.lorywidth;
01276         if(locBz > 0.f) {lorywidth_ = -lorywidth_;}
01277         lorxwidth_ = thePixelTemp_[index_id_].head.lorxwidth;
01278         
01279   }
01280         
01281   return success_;
01282 } // interpolate
01283 
01284 
01285 
01286 
01287 
01288 // ************************************************************************************************************ 
01293 // ************************************************************************************************************ 
01294 bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta)
01295 {
01296     // Interpolate for a new set of track angles 
01297     
01298     // Local variables 
01299     float locBz;
01300         locBz = -1.f;
01301         if(cotbeta < 0.f) {locBz = -locBz;}
01302     return SiPixelTemplate::interpolate(id, cotalpha, cotbeta, locBz);
01303 }
01304 
01305 
01306 
01307 
01308 // ************************************************************************************************************ 
01316 // ************************************************************************************************************ 
01317   void SiPixelTemplate::ysigma2(int fypix, int lypix, float sythr, float ysum[25], float ysig2[25])
01318   
01319 {
01320     // Interpolate using quantities already stored in the private variables
01321     
01322     // Local variables 
01323     int i;
01324         float sigi, sigi2, sigi3, sigi4, symax, qscale;
01325         
01326     // Make sure that input is OK
01327     
01328 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01329     if(fypix < 2 || fypix >= BYM2) {
01330                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ysigma2 called with fypix = " << fypix << std::endl;
01331         }
01332 #else
01333         assert(fypix > 1 && fypix < BYM2);
01334 #endif
01335 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01336            if(lypix < fypix || lypix >= BYM2) {
01337                   throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ysigma2 called with lypix/fypix = " << lypix << "/" << fypix << std::endl;
01338                 }
01339 #else
01340                 assert(lypix >= fypix && lypix < BYM2);
01341 #endif
01342                      
01343 // Define the maximum signal to use in the parameterization 
01344 
01345        symax = symax_;
01346            if(symax_ > syparmax_) {symax = syparmax_;}
01347            
01348 // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis 
01349 
01350            for(i=fypix-2; i<=lypix+2; ++i) {
01351                   if(i < fypix || i > lypix) {
01352            
01353 // Nearest pseudopixels have uncertainties of 50% of threshold, next-nearest have 10% of threshold
01354 
01355                          ysig2[i] = s50_*s50_;
01356                   } else {
01357                          if(ysum[i] < symax) {
01358                                 sigi = ysum[i];
01359                                 qscale = 1.f;
01360                          } else {
01361                                 sigi = symax;
01362                                 qscale = ysum[i]/symax;
01363                          }
01364                          sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
01365                          if(i <= BHY) {
01366                                 ysig2[i] = (1.f-yratio_)*
01367                                 (yparl_[0][0]+yparl_[0][1]*sigi+yparl_[0][2]*sigi2+yparl_[0][3]*sigi3+yparl_[0][4]*sigi4)
01368                                 + yratio_*
01369                                 (yparh_[0][0]+yparh_[0][1]*sigi+yparh_[0][2]*sigi2+yparh_[0][3]*sigi3+yparh_[0][4]*sigi4);
01370                          } else {
01371                                 ysig2[i] = (1.f-yratio_)*
01372                                 (yparl_[1][0]+yparl_[1][1]*sigi+yparl_[1][2]*sigi2+yparl_[1][3]*sigi3+yparl_[1][4]*sigi4)
01373                                 + yratio_*
01374                             (yparh_[1][0]+yparh_[1][1]*sigi+yparh_[1][2]*sigi2+yparh_[1][3]*sigi3+yparh_[1][4]*sigi4);
01375                          }
01376                          ysig2[i] *=qscale;
01377                      if(ysum[i] > sythr) {ysig2[i] = 1.e8;}
01378                          if(ysig2[i] <= 0.f) {LOGERROR("SiPixelTemplate") << "neg y-error-squared, id = " << id_current_ << ", index = " << index_id_ << 
01379                          ", cot(alpha) = " << cota_current_ << ", cot(beta) = " << cotb_current_ <<  ", sigi = " << sigi << ENDL;}
01380               }
01381            }
01382         
01383         return;
01384         
01385 } // End ysigma2
01386 
01387 
01388 // ************************************************************************************************************ 
01394 // ************************************************************************************************************ 
01395 void SiPixelTemplate::ysigma2(float qpixel, int index, float& ysig2)
01396 
01397 {
01398     // Interpolate using quantities already stored in the private variables
01399     
01400     // Local variables 
01401         float sigi, sigi2, sigi3, sigi4, symax, qscale, err2;
01402         
01403     // Make sure that input is OK
01404     
01405 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01406     if(index < 2 || index >= BYM2) {
01407                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ysigma2 called with index = " << index << std::endl;
01408         }
01409 #else
01410         assert(index > 1 && index < BYM2);
01411 #endif
01412         
01413         // Define the maximum signal to use in the parameterization 
01414         
01415         symax = symax_;
01416         if(symax_ > syparmax_) {symax = syparmax_;}
01417         
01418         // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis 
01419         
01420                         if(qpixel < symax) {
01421                                 sigi = qpixel;
01422                                 qscale = 1.f;
01423                         } else {
01424                                 sigi = symax;
01425                                 qscale = qpixel/symax;
01426                         }
01427                         sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
01428                         if(index <= BHY) {
01429                                 err2 = (1.f-yratio_)*
01430                                 (yparl_[0][0]+yparl_[0][1]*sigi+yparl_[0][2]*sigi2+yparl_[0][3]*sigi3+yparl_[0][4]*sigi4)
01431                                 + yratio_*
01432                                 (yparh_[0][0]+yparh_[0][1]*sigi+yparh_[0][2]*sigi2+yparh_[0][3]*sigi3+yparh_[0][4]*sigi4);
01433                         } else {
01434                                 err2 = (1.f-yratio_)*
01435                                 (yparl_[1][0]+yparl_[1][1]*sigi+yparl_[1][2]*sigi2+yparl_[1][3]*sigi3+yparl_[1][4]*sigi4)
01436                                 + yratio_*
01437                             (yparh_[1][0]+yparh_[1][1]*sigi+yparh_[1][2]*sigi2+yparh_[1][3]*sigi3+yparh_[1][4]*sigi4);
01438                         }
01439                         ysig2 =qscale*err2;
01440                         if(ysig2 <= 0.f) {LOGERROR("SiPixelTemplate") << "neg y-error-squared, id = " << id_current_ << ", index = " << index_id_ << 
01441                         ", cot(alpha) = " << cota_current_ << ", cot(beta) = " << cotb_current_ <<  ", sigi = " << sigi << ENDL;}
01442         
01443         return;
01444         
01445 } // End ysigma2
01446 
01447 
01448 
01449 
01450 
01451 // ************************************************************************************************************ 
01459 // ************************************************************************************************************ 
01460   void SiPixelTemplate::xsigma2(int fxpix, int lxpix, float sxthr, float xsum[11], float xsig2[11])
01461   
01462 {
01463     // Interpolate using quantities already stored in the private variables
01464     
01465     // Local variables 
01466     int i;
01467         float sigi, sigi2, sigi3, sigi4, yint, sxmax, x0, qscale;
01468         
01469     // Make sure that input is OK
01470     
01471 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01472                   if(fxpix < 2 || fxpix >= BXM2) {
01473                          throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xsigma2 called with fxpix = " << fxpix << std::endl;
01474                    }
01475 #else
01476                    assert(fxpix > 1 && fxpix < BXM2);
01477 #endif
01478 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01479                          if(lxpix < fxpix || lxpix >= BXM2) {
01480                                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xsigma2 called with lxpix/fxpix = " << lxpix << "/" << fxpix << std::endl;
01481                          }
01482 #else
01483                          assert(lxpix >= fxpix && lxpix < BXM2);
01484 #endif
01485                      
01486 // Define the maximum signal to use in the parameterization 
01487 
01488        sxmax = sxmax_;
01489            if(sxmax_ > sxparmax_) {sxmax = sxparmax_;}
01490            
01491 // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis 
01492 
01493            for(i=fxpix-2; i<=lxpix+2; ++i) {
01494                   if(i < fxpix || i > lxpix) {
01495            
01496 // Nearest pseudopixels have uncertainties of 50% of threshold, next-nearest have 10% of threshold
01497 
01498                          xsig2[i] = s50_*s50_;
01499                   } else {
01500                          if(xsum[i] < sxmax) {
01501                                 sigi = xsum[i];
01502                                 qscale = 1.f;
01503                          } else {
01504                                 sigi = sxmax;
01505                                 qscale = xsum[i]/sxmax;
01506                          }
01507                          sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
01508                          
01509 // First, do the cotbeta interpolation                   
01510                          
01511                          if(i <= BHX) {
01512                                 yint = (1.f-yratio_)*
01513                                 (xparly0_[0][0]+xparly0_[0][1]*sigi+xparly0_[0][2]*sigi2+xparly0_[0][3]*sigi3+xparly0_[0][4]*sigi4)
01514                                 + yratio_*
01515                                 (xparhy0_[0][0]+xparhy0_[0][1]*sigi+xparhy0_[0][2]*sigi2+xparhy0_[0][3]*sigi3+xparhy0_[0][4]*sigi4);
01516                          } else {
01517                                 yint = (1.f-yratio_)*
01518                                 (xparly0_[1][0]+xparly0_[1][1]*sigi+xparly0_[1][2]*sigi2+xparly0_[1][3]*sigi3+xparly0_[1][4]*sigi4)
01519                                 + yratio_*
01520                             (xparhy0_[1][0]+xparhy0_[1][1]*sigi+xparhy0_[1][2]*sigi2+xparhy0_[1][3]*sigi3+xparhy0_[1][4]*sigi4);
01521                          }
01522                          
01523 // Next, do the cotalpha interpolation                   
01524                          
01525                          if(i <= BHX) {
01526                                 xsig2[i] = (1.f-xxratio_)*
01527                                 (xparl_[0][0]+xparl_[0][1]*sigi+xparl_[0][2]*sigi2+xparl_[0][3]*sigi3+xparl_[0][4]*sigi4)
01528                                 + xxratio_*
01529                                 (xparh_[0][0]+xparh_[0][1]*sigi+xparh_[0][2]*sigi2+xparh_[0][3]*sigi3+xparh_[0][4]*sigi4);
01530                          } else {
01531                                 xsig2[i] = (1.f-xxratio_)*
01532                                 (xparl_[1][0]+xparl_[1][1]*sigi+xparl_[1][2]*sigi2+xparl_[1][3]*sigi3+xparl_[1][4]*sigi4)
01533                                 + xxratio_*
01534                             (xparh_[1][0]+xparh_[1][1]*sigi+xparh_[1][2]*sigi2+xparh_[1][3]*sigi3+xparh_[1][4]*sigi4);
01535                          }
01536                          
01537 // Finally, get the mid-point value of the cotalpha function                     
01538                          
01539                          if(i <= BHX) {
01540                                 x0 = xpar0_[0][0]+xpar0_[0][1]*sigi+xpar0_[0][2]*sigi2+xpar0_[0][3]*sigi3+xpar0_[0][4]*sigi4;
01541                          } else {
01542                                 x0 = xpar0_[1][0]+xpar0_[1][1]*sigi+xpar0_[1][2]*sigi2+xpar0_[1][3]*sigi3+xpar0_[1][4]*sigi4;
01543                          }
01544                          
01545 // Finally, rescale the yint value for cotalpha variation                        
01546                          
01547                          if(x0 != 0.f) {xsig2[i] = xsig2[i]/x0 * yint;}
01548                          xsig2[i] *=qscale;
01549                      if(xsum[i] > sxthr) {xsig2[i] = 1.e8;}
01550                          if(xsig2[i] <= 0.f) {LOGERROR("SiPixelTemplate") << "neg x-error-squared, id = " << id_current_ << ", index = " << index_id_ << 
01551                          ", cot(alpha) = " << cota_current_ << ", cot(beta) = " << cotb_current_  << ", sigi = " << sigi << ENDL;}
01552               }
01553            }
01554         
01555         return;
01556         
01557 } // End xsigma2
01558 
01559 
01560 
01561 
01562 
01563 
01564 // ************************************************************************************************************ 
01568 // ************************************************************************************************************ 
01569   float SiPixelTemplate::yflcorr(int binq, float qfly)
01570   
01571 {
01572     // Interpolate using quantities already stored in the private variables
01573     
01574     // Local variables 
01575         float qfl, qfl2, qfl3, qfl4, qfl5, dy;
01576         
01577     // Make sure that input is OK
01578     
01579 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01580         if(binq < 0 || binq > 3) {
01581            throw cms::Exception("DataCorrupt") << "SiPixelTemplate::yflcorr called with binq = " << binq << std::endl;
01582         }
01583 #else
01584          assert(binq >= 0 && binq < 4);
01585 #endif
01586 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01587          if(fabs((double)qfly) > 1.) {
01588                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::yflcorr called with qfly = " << qfly << std::endl;
01589          }
01590 #else
01591          assert(fabs((double)qfly) <= 1.);
01592 #endif
01593                      
01594 // Define the maximum signal to allow before de-weighting a pixel 
01595 
01596        qfl = qfly;
01597 
01598        if(qfl < -0.9f) {qfl = -0.9f;}
01599            if(qfl > 0.9f) {qfl = 0.9f;}
01600            
01601 // Interpolate between the two polynomials
01602 
01603            qfl2 = qfl*qfl; qfl3 = qfl2*qfl; qfl4 = qfl3*qfl; qfl5 = qfl4*qfl;
01604            dy = (1.f-yratio_)*(yflparl_[binq][0]+yflparl_[binq][1]*qfl+yflparl_[binq][2]*qfl2+yflparl_[binq][3]*qfl3+yflparl_[binq][4]*qfl4+yflparl_[binq][5]*qfl5)
01605                   + yratio_*(yflparh_[binq][0]+yflparh_[binq][1]*qfl+yflparh_[binq][2]*qfl2+yflparh_[binq][3]*qfl3+yflparh_[binq][4]*qfl4+yflparh_[binq][5]*qfl5);
01606         
01607         return dy;
01608         
01609 } // End yflcorr
01610 
01611 
01612 
01613 
01614 
01615 
01616 // ************************************************************************************************************ 
01620 // ************************************************************************************************************ 
01621   float SiPixelTemplate::xflcorr(int binq, float qflx)
01622   
01623 {
01624     // Interpolate using quantities already stored in the private variables
01625     
01626     // Local variables 
01627         float qfl, qfl2, qfl3, qfl4, qfl5, dx;
01628         
01629     // Make sure that input is OK
01630     
01631 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01632         if(binq < 0 || binq > 3) {
01633                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xflcorr called with binq = " << binq << std::endl;
01634         }
01635 #else
01636         assert(binq >= 0 && binq < 4);
01637 #endif
01638 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01639         if(fabs((double)qflx) > 1.) {
01640                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xflcorr called with qflx = " << qflx << std::endl;
01641         }
01642 #else
01643         assert(fabs((double)qflx) <= 1.);
01644 #endif
01645                      
01646 // Define the maximum signal to allow before de-weighting a pixel 
01647 
01648        qfl = qflx;
01649 
01650        if(qfl < -0.9f) {qfl = -0.9f;}
01651            if(qfl > 0.9f) {qfl = 0.9f;}
01652            
01653 // Interpolate between the two polynomials
01654 
01655            qfl2 = qfl*qfl; qfl3 = qfl2*qfl; qfl4 = qfl3*qfl; qfl5 = qfl4*qfl;
01656            dx = (1.f - yxratio_)*((1.f-xxratio_)*(xflparll_[binq][0]+xflparll_[binq][1]*qfl+xflparll_[binq][2]*qfl2+xflparll_[binq][3]*qfl3+xflparll_[binq][4]*qfl4+xflparll_[binq][5]*qfl5)
01657                   + xxratio_*(xflparlh_[binq][0]+xflparlh_[binq][1]*qfl+xflparlh_[binq][2]*qfl2+xflparlh_[binq][3]*qfl3+xflparlh_[binq][4]*qfl4+xflparlh_[binq][5]*qfl5))
01658               + yxratio_*((1.f-xxratio_)*(xflparhl_[binq][0]+xflparhl_[binq][1]*qfl+xflparhl_[binq][2]*qfl2+xflparhl_[binq][3]*qfl3+xflparhl_[binq][4]*qfl4+xflparhl_[binq][5]*qfl5)
01659                   + xxratio_*(xflparhh_[binq][0]+xflparhh_[binq][1]*qfl+xflparhh_[binq][2]*qfl2+xflparhh_[binq][3]*qfl3+xflparhh_[binq][4]*qfl4+xflparhh_[binq][5]*qfl5));
01660         
01661         return dx;
01662         
01663 } // End xflcorr
01664 
01665 
01666 
01667 // ************************************************************************************************************ 
01672 // ************************************************************************************************************ 
01673   void SiPixelTemplate::ytemp(int fybin, int lybin, float ytemplate[41][BYSIZE])
01674   
01675 {
01676     // Retrieve already interpolated quantities
01677     
01678     // Local variables 
01679     int i, j;
01680 
01681    // Verify that input parameters are in valid range
01682 
01683 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01684         if(fybin < 0 || fybin > 40) {
01685                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ytemp called with fybin = " << fybin << std::endl;
01686         }
01687 #else
01688         assert(fybin >= 0 && fybin < 41);
01689 #endif
01690 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01691         if(lybin < 0 || lybin > 40) {
01692                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ytemp called with lybin = " << lybin << std::endl;
01693         }
01694 #else
01695         assert(lybin >= 0 && lybin < 41);
01696 #endif
01697 
01698 // Build the y-template, the central 25 bins are here in all cases
01699         
01700         for(i=0; i<9; ++i) {
01701            for(j=0; j<BYSIZE; ++j) {
01702                         ytemplate[i+16][j]=ytemp_[i][j];
01703            }
01704         }
01705         for(i=0; i<8; ++i) {
01706            ytemplate[i+8][BYM1] = 0.f;
01707            for(j=0; j<BYM1; ++j) {
01708               ytemplate[i+8][j]=ytemp_[i][j+1];
01709            }
01710         }
01711         for(i=1; i<9; ++i) {
01712            ytemplate[i+24][0] = 0.f;
01713            for(j=0; j<BYM1; ++j) {
01714               ytemplate[i+24][j+1]=ytemp_[i][j];
01715            }
01716         }
01717         
01718 //  Add more bins if needed
01719 
01720         if(fybin < 8) {
01721            for(i=0; i<8; ++i) {
01722               ytemplate[i][BYM2] = 0.f;
01723               ytemplate[i][BYM1] = 0.f;
01724               for(j=0; j<BYM2; ++j) {
01725                 ytemplate[i][j]=ytemp_[i][j+2];
01726               }
01727            }
01728         }
01729         if(lybin > 32) {
01730            for(i=1; i<9; ++i) {
01731           ytemplate[i+32][0] = 0.f;
01732               ytemplate[i+32][1] = 0.f;
01733               for(j=0; j<BYM2; ++j) {
01734                  ytemplate[i+32][j+2]=ytemp_[i][j];
01735               }
01736            }
01737         }
01738         
01739         return;
01740         
01741 } // End ytemp
01742 
01743 
01744 
01745 // ************************************************************************************************************ 
01750 // ************************************************************************************************************ 
01751   void SiPixelTemplate::xtemp(int fxbin, int lxbin, float xtemplate[41][BXSIZE])
01752   
01753 {
01754     // Retrieve already interpolated quantities
01755     
01756     // Local variables 
01757     int i, j;
01758 
01759    // Verify that input parameters are in valid range
01760 
01761 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01762         if(fxbin < 0 || fxbin > 40) {
01763                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp called with fxbin = " << fxbin << std::endl;
01764         }
01765 #else
01766         assert(fxbin >= 0 && fxbin < 41);
01767 #endif
01768 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01769         if(lxbin < 0 || lxbin > 40) {
01770                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp called with lxbin = " << lxbin << std::endl;
01771         }
01772 #else
01773         assert(lxbin >= 0 && lxbin < 41);
01774 #endif
01775 
01776 // Build the x-template, the central 25 bins are here in all cases
01777         
01778         for(i=0; i<9; ++i) {
01779            for(j=0; j<BXSIZE; ++j) {
01780               xtemplate[i+16][j]=xtemp_[i][j];
01781            }
01782         }
01783         for(i=0; i<8; ++i) {
01784            xtemplate[i+8][BXM1] = 0.f;
01785            for(j=0; j<BXM1; ++j) {
01786               xtemplate[i+8][j]=xtemp_[i][j+1];
01787            }
01788         }
01789         for(i=1; i<9; ++i) {
01790            xtemplate[i+24][0] = 0.f;
01791            for(j=0; j<BXM1; ++j) {
01792               xtemplate[i+24][j+1]=xtemp_[i][j];
01793            }
01794         }
01795         
01796 //  Add more bins if needed     
01797         
01798         if(fxbin < 8) {
01799            for(i=0; i<8; ++i) {
01800           xtemplate[i][BXM2] = 0.f;
01801               xtemplate[i][BXM1] = 0.f;
01802               for(j=0; j<BXM2; ++j) {
01803                 xtemplate[i][j]=xtemp_[i][j+2];
01804               }
01805            }
01806         }
01807         if(lxbin > 32) {
01808            for(i=1; i<9; ++i) {
01809           xtemplate[i+32][0] = 0.f;
01810               xtemplate[i+32][1] = 0.f;
01811               for(j=0; j<BXM2; ++j) {
01812                 xtemplate[i+32][j+2]=xtemp_[i][j];
01813               }
01814            }
01815         }
01816         
01817         return;
01818         
01819 } // End xtemp
01820 
01821 
01822 // ************************************************************************************************************ 
01824 // ************************************************************************************************************ 
01825 int SiPixelTemplate::cytemp()
01826 
01827 {
01828         // Retrieve already interpolated quantities
01829         
01830         // Local variables 
01831         int j;
01832         
01833 // Analyze only pixels along the central entry
01834 // First, find the maximum signal and then work out to the edges
01835         
01836         float sigmax = 0.f;
01837         float qedge = 2.*s50_;
01838         int jmax = -1;
01839         
01840         for(j=0; j<BYSIZE; ++j) {
01841                 if(ytemp_[4][j] > sigmax) {
01842                         sigmax = ytemp_[4][j];
01843                         jmax = j;
01844                 }       
01845         }
01846         if(sigmax < qedge) {qedge = s50_;}
01847         if(sigmax < qedge || jmax<1 || jmax>BYM2) {return -1;}
01848         
01849 //  Now search forward and backward
01850         
01851         int jend = jmax;
01852         
01853         for(j=jmax+1; j<BYM1; ++j) {
01854                 if(ytemp_[4][j] < qedge) break;
01855            jend = j;
01856         }
01857 
01858    int jbeg = jmax;
01859 
01860 for(j=jmax-1; j>0; --j) {
01861                 if(ytemp_[4][j] < qedge) break;
01862                 jbeg = j;
01863         }
01864                 
01865         return (jbeg+jend)/2;
01866         
01867 } // End cytemp
01868 
01869 
01870 
01871 // ************************************************************************************************************ 
01873 // ************************************************************************************************************ 
01874 int SiPixelTemplate::cxtemp()
01875 
01876 {
01877         // Retrieve already interpolated quantities
01878         
01879         // Local variables 
01880         int j;
01881         
01882         // Analyze only pixels along the central entry
01883         // First, find the maximum signal and then work out to the edges
01884         
01885         float sigmax = 0.f;
01886         float qedge = 2.*s50_;
01887         int jmax = -1;
01888         
01889         for(j=0; j<BXSIZE; ++j) {
01890                 if(xtemp_[4][j] > sigmax) {
01891                         sigmax = xtemp_[4][j];
01892                         jmax = j;
01893                 }       
01894         }
01895         if(sigmax < qedge) {qedge = s50_;}
01896         if(sigmax < qedge || jmax<1 || jmax>BXM2) {return -1;}
01897         
01898         //  Now search forward and backward
01899         
01900         int jend = jmax;
01901         
01902         for(j=jmax+1; j<BXM1; ++j) {
01903                 if(xtemp_[4][j] < qedge) break;
01904            jend = j;
01905    }
01906 
01907    int jbeg = jmax;
01908 
01909    for(j=jmax-1; j>0; --j) {
01910            if(xtemp_[4][j] < qedge) break;
01911            jbeg = j;
01912    }
01913 
01914    return (jbeg+jend)/2;
01915 
01916 } // End cxtemp
01917 
01918 
01919 // ************************************************************************************************************ 
01923 // ************************************************************************************************************ 
01924 void SiPixelTemplate::ytemp3d_int(int nypix, int& nybins)
01925 
01926 {
01927         
01928     // Retrieve already interpolated quantities
01929     
01930     // Local variables 
01931     int i, j, k;
01932         int ioff0, ioffp, ioffm;
01933     
01934     // Verify that input parameters are in valid range
01935     
01936 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01937         if(nypix < 1 || nypix >= BYM3) {
01938                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ytemp3d called with nypix = " << nypix << std::endl;
01939         }
01940 #else
01941         assert(nypix > 0 && nypix < BYM3);
01942 #endif
01943         
01944     // Calculate the size of the shift in pixels needed to span the entire cluster
01945     
01946     float diff = fabsf(nypix - clsleny_)/2. + 1.f;
01947         int nshift = (int)diff;
01948         if((diff - nshift) > 0.5f) {++nshift;}
01949     
01950     // Calculate the number of bins needed to specify each hit range
01951     
01952     nybins_ = 9 + 16*nshift;
01953         
01954     // Create a 2-d working template with the correct size
01955     
01956         temp2dy_.resize(boost::extents[nybins_][BYSIZE]);
01957         
01958     //  The 9 central bins are copied from the interpolated private store
01959     
01960         ioff0 = 8*nshift;
01961         
01962         for(i=0; i<9; ++i) {
01963         for(j=0; j<BYSIZE; ++j) {
01964             temp2dy_[i+ioff0][j]=ytemp_[i][j];
01965         }
01966         }
01967         
01968     // Add the +- shifted templates     
01969     
01970         for(k=1; k<=nshift; ++k) {
01971         ioffm=ioff0-k*8;
01972         for(i=0; i<8; ++i) {
01973             for(j=0; j<k; ++j) {
01974                 temp2dy_[i+ioffm][BYM1-j] = 0.f;
01975             }
01976             for(j=0; j<BYSIZE-k; ++j) {
01977                 temp2dy_[i+ioffm][j]=ytemp_[i][j+k];
01978             }
01979         }
01980         ioffp=ioff0+k*8;
01981         for(i=1; i<9; ++i) {
01982             for(j=0; j<k; ++j) {
01983                 temp2dy_[i+ioffp][j] = 0.f;
01984             }
01985             for(j=0; j<BYSIZE-k; ++j) {
01986                 temp2dy_[i+ioffp][j+k]=ytemp_[i][j];
01987             }
01988         }
01989         }
01990     
01991     nybins = nybins_;   
01992         return;
01993         
01994 } // End ytemp3d_int
01995 
01996 
01997 // ************************************************************************************************************ 
02001 // ************************************************************************************************************ 
02002 void SiPixelTemplate::ytemp3d(int i, int j, std::vector<float>& ytemplate)
02003 
02004 {       
02005     // Sum two 2-d templates to make the 3-d template
02006         if(i >= 0 && i < nybins_ && j <= i) {
02007        for(int k=0; k<BYSIZE; ++k) {
02008           ytemplate[k]=temp2dy_[i][k]+temp2dy_[j][k];   
02009        }
02010     } else {
02011         for(int k=0; k<BYSIZE; ++k) {
02012             ytemplate[k]=0.;    
02013         } 
02014     }
02015         
02016         return;
02017         
02018 } // End ytemp3d
02019 
02020 
02021 // ************************************************************************************************************ 
02025 // ************************************************************************************************************ 
02026   void SiPixelTemplate::xtemp3d_int(int nxpix, int& nxbins)
02027   
02028 {       
02029     // Retrieve already interpolated quantities
02030     
02031     // Local variables 
02032     int i, j, k;
02033         int ioff0, ioffp, ioffm;
02034 
02035    // Verify that input parameters are in valid range
02036 
02037 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
02038         if(nxpix < 1 || nxpix >= BXM3) {
02039            throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp3d called with nxpix = " << nxpix << std::endl;
02040         }
02041 #else
02042         assert(nxpix > 0 && nxpix < BXM3);
02043 #endif
02044         
02045 // Calculate the size of the shift in pixels needed to span the entire cluster
02046 
02047     float diff = fabsf(nxpix - clslenx_)/2. + 1.f;
02048         int nshift = (int)diff;
02049         if((diff - nshift) > 0.5f) {++nshift;}
02050 
02051 // Calculate the number of bins needed to specify each hit range
02052 
02053     nxbins_ = 9 + 16*nshift;
02054         
02055 // Create a 2-d working template with the correct size
02056 
02057         temp2dx_.resize(boost::extents[nxbins_][BXSIZE]);
02058         
02059 //  The 9 central bins are copied from the interpolated private store
02060 
02061         ioff0 = 8*nshift;
02062         
02063         for(i=0; i<9; ++i) {
02064            for(j=0; j<BXSIZE; ++j) {
02065           temp2dx_[i+ioff0][j]=xtemp_[i][j];
02066            }
02067         }
02068         
02069 // Add the +- shifted templates 
02070 
02071         for(k=1; k<=nshift; ++k) {
02072           ioffm=ioff0-k*8;
02073           for(i=0; i<8; ++i) {
02074              for(j=0; j<k; ++j) {
02075                 temp2dx_[i+ioffm][BXM1-j] = 0.f;
02076                  }
02077              for(j=0; j<BXSIZE-k; ++j) {
02078                 temp2dx_[i+ioffm][j]=xtemp_[i][j+k];
02079              }
02080            }
02081            ioffp=ioff0+k*8;
02082            for(i=1; i<9; ++i) {
02083               for(j=0; j<k; ++j) {
02084                  temp2dx_[i+ioffp][j] = 0.f;
02085               }
02086               for(j=0; j<BXSIZE-k; ++j) {
02087                  temp2dx_[i+ioffp][j+k]=xtemp_[i][j];
02088               }
02089            }
02090         }
02091     
02092     nxbins = nxbins_;                                   
02093         
02094         return;
02095         
02096 } // End xtemp3d_int
02097 
02098 
02099 
02100 // ************************************************************************************************************ 
02104 // ************************************************************************************************************ 
02105 void SiPixelTemplate::xtemp3d(int i, int j, std::vector<float>& xtemplate)
02106 
02107 {       
02108     // Sum two 2-d templates to make the 3-d template
02109         if(i >= 0 && i < nxbins_ && j <= i) {
02110         for(int k=0; k<BXSIZE; ++k) {
02111             xtemplate[k]=temp2dx_[i][k]+temp2dx_[j][k]; 
02112         }
02113     } else {
02114         for(int k=0; k<BXSIZE; ++k) {
02115             xtemplate[k]=0.;    
02116         } 
02117     }
02118         
02119         return;
02120         
02121 } // End xtemp3d
02122 
02123 
02124 // ************************************************************************************************************ 
02148 // ************************************************************************************************************ 
02149 int SiPixelTemplate::qbin(int id, float cotalpha, float cotbeta, float locBz, float qclus, float& pixmx, float& sigmay, float& deltay, float& sigmax, float& deltax, 
02150                           float& sy1, float& dy1, float& sy2, float& dy2, float& sx1, float& dx1, float& sx2, float& dx2, float& lorywidth, float& lorxwidth)
02151                  
02152 {
02153     // Interpolate for a new set of track angles 
02154     
02155     // Local variables 
02156     int i, binq;
02157         int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, index;
02158         float yratio, yxratio, xxratio;
02159         float acotb, qscale, qavg, qmin, qmin2, fq, qtotal, qcorrect, cotb, cotalpha0;
02160         float yavggen[4], yrmsgen[4], xavggen[4], xrmsgen[4];
02161         bool flip_y;
02162         
02163 
02164 // Find the index corresponding to id
02165 
02166        index = -1;
02167        for(i=0; i<(int)thePixelTemp_.size(); ++i) {
02168         
02169               if(id == thePixelTemp_[i].head.ID) {
02170            
02171                  index = i;
02172                      break;
02173           }
02174             }
02175          
02176 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
02177            if(index < 0 || index >= (int)thePixelTemp_.size()) {
02178               throw cms::Exception("DataCorrupt") << "SiPixelTemplate::qbin can't find needed template ID = " << id << std::endl;
02179            }
02180 #else
02181            assert(index >= 0 && index < (int)thePixelTemp_.size());
02182 #endif
02183          
02184 //              
02185 
02186 // Interpolate the absolute value of cot(beta)     
02187     
02188     acotb = fabs((double)cotbeta);
02189         
02190 //      qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)       
02191 
02192         //      qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)       
02193         
02194         cotalpha0 =  thePixelTemp_[index].enty[0].cotalpha;
02195         qcorrect=std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta*cotbeta+cotalpha0*cotalpha0));
02196                                 
02197         // for some cosmics, the ususal gymnastics are incorrect   
02198         
02199         if(thePixelTemp_[index].head.Dtype == 0) {
02200                 cotb = acotb;
02201                 flip_y = false;
02202                 if(cotbeta < 0.f) {flip_y = true;}
02203         } else {
02204             if(locBz < 0.f) {
02205                         cotb = cotbeta;
02206                         flip_y = false;
02207                 } else {
02208                         cotb = -cotbeta;
02209                         flip_y = true;
02210                 }       
02211         }
02212         
02213         // Copy the charge scaling factor to the private variable     
02214                 
02215            qscale = thePixelTemp_[index].head.qscale;
02216                 
02217        Ny = thePixelTemp_[index].head.NTy;
02218        Nyx = thePixelTemp_[index].head.NTyx;
02219        Nxx = thePixelTemp_[index].head.NTxx;
02220                 
02221 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
02222                 if(Ny < 2 || Nyx < 1 || Nxx < 2) {
02223                         throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
02224                 }
02225 #else
02226                 assert(Ny > 1 && Nyx > 0 && Nxx > 1);
02227 #endif
02228         
02229 // next, loop over all y-angle entries   
02230 
02231            ilow = 0;
02232            yratio = 0.f;
02233 
02234            if(cotb >= thePixelTemp_[index].enty[Ny-1].cotbeta) {
02235         
02236                ilow = Ny-2;
02237                    yratio = 1.f;
02238                 
02239            } else {
02240            
02241               if(cotb >= thePixelTemp_[index].enty[0].cotbeta) {
02242 
02243              for (i=0; i<Ny-1; ++i) { 
02244     
02245                 if( thePixelTemp_[index].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index].enty[i+1].cotbeta) {
02246                   
02247                        ilow = i;
02248                            yratio = (cotb - thePixelTemp_[index].enty[i].cotbeta)/(thePixelTemp_[index].enty[i+1].cotbeta - thePixelTemp_[index].enty[i].cotbeta);
02249                            break;                        
02250                         }
02251                  }
02252                   } 
02253            }
02254         
02255            ihigh=ilow + 1;
02256                           
02257 // Interpolate/store all y-related quantities (flip displacements when flip_y)
02258 
02259            qavg = (1.f - yratio)*thePixelTemp_[index].enty[ilow].qavg + yratio*thePixelTemp_[index].enty[ihigh].qavg;
02260            qavg *= qcorrect;
02261            dy1 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].dyone + yratio*thePixelTemp_[index].enty[ihigh].dyone;
02262            if(flip_y) {dy1 = -dy1;}
02263            sy1 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].syone + yratio*thePixelTemp_[index].enty[ihigh].syone;
02264            dy2 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].dytwo + yratio*thePixelTemp_[index].enty[ihigh].dytwo;
02265            if(flip_y) {dy2 = -dy2;}
02266            sy2 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].sytwo + yratio*thePixelTemp_[index].enty[ihigh].sytwo;
02267            qmin = (1.f - yratio)*thePixelTemp_[index].enty[ilow].qmin + yratio*thePixelTemp_[index].enty[ihigh].qmin;
02268            qmin *= qcorrect;
02269            qmin2 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].qmin2 + yratio*thePixelTemp_[index].enty[ihigh].qmin2;
02270            qmin2 *= qcorrect;
02271            for(i=0; i<4; ++i) {
02272               yavggen[i]=(1.f - yratio)*thePixelTemp_[index].enty[ilow].yavggen[i] + yratio*thePixelTemp_[index].enty[ihigh].yavggen[i];
02273               if(flip_y) {yavggen[i] = -yavggen[i];}
02274               yrmsgen[i]=(1.f - yratio)*thePixelTemp_[index].enty[ilow].yrmsgen[i] + yratio*thePixelTemp_[index].enty[ihigh].yrmsgen[i];
02275            }
02276            
02277         
02278 // next, loop over all x-angle entries, first, find relevant y-slices   
02279         
02280            iylow = 0;
02281            yxratio = 0.f;
02282 
02283            if(acotb >= thePixelTemp_[index].entx[Nyx-1][0].cotbeta) {
02284         
02285                iylow = Nyx-2;
02286                    yxratio = 1.f;
02287                 
02288            } else if(acotb >= thePixelTemp_[index].entx[0][0].cotbeta) {
02289 
02290           for (i=0; i<Nyx-1; ++i) { 
02291     
02292              if( thePixelTemp_[index].entx[i][0].cotbeta <= acotb && acotb < thePixelTemp_[index].entx[i+1][0].cotbeta) {
02293                   
02294                     iylow = i;
02295                         yxratio = (acotb - thePixelTemp_[index].entx[i][0].cotbeta)/(thePixelTemp_[index].entx[i+1][0].cotbeta - thePixelTemp_[index].entx[i][0].cotbeta);
02296                         break;                   
02297                      }
02298               }
02299            }
02300         
02301            iyhigh=iylow + 1;
02302 
02303            ilow = 0;
02304            xxratio = 0.f;
02305 
02306            if(cotalpha >= thePixelTemp_[index].entx[0][Nxx-1].cotalpha) {
02307         
02308                ilow = Nxx-2;
02309                    xxratio = 1.f;
02310                 
02311            } else {
02312            
02313               if(cotalpha >= thePixelTemp_[index].entx[0][0].cotalpha) {
02314 
02315              for (i=0; i<Nxx-1; ++i) { 
02316     
02317                 if( thePixelTemp_[index].entx[0][i].cotalpha <= cotalpha && cotalpha < thePixelTemp_[index].entx[0][i+1].cotalpha) {
02318                   
02319                        ilow = i;
02320                            xxratio = (cotalpha - thePixelTemp_[index].entx[0][i].cotalpha)/(thePixelTemp_[index].entx[0][i+1].cotalpha - thePixelTemp_[index].entx[0][i].cotalpha);
02321                            break;
02322                             }
02323                      }
02324                   } 
02325            }
02326         
02327            ihigh=ilow + 1;
02328                           
02329            dx1 = (1.f - xxratio)*thePixelTemp_[index].entx[0][ilow].dxone + xxratio*thePixelTemp_[index].entx[0][ihigh].dxone;
02330            sx1 = (1.f - xxratio)*thePixelTemp_[index].entx[0][ilow].sxone + xxratio*thePixelTemp_[index].entx[0][ihigh].sxone;
02331            dx2 = (1.f - xxratio)*thePixelTemp_[index].entx[0][ilow].dxtwo + xxratio*thePixelTemp_[index].entx[0][ihigh].dxtwo;
02332            sx2 = (1.f - xxratio)*thePixelTemp_[index].entx[0][ilow].sxtwo + xxratio*thePixelTemp_[index].entx[0][ihigh].sxtwo;
02333                           
02334 // pixmax is the maximum allowed pixel charge (used for truncation)
02335 
02336            pixmx=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].pixmax + xxratio*thePixelTemp_[index].entx[iylow][ihigh].pixmax)
02337                           +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].pixmax + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].pixmax);
02338                           
02339            for(i=0; i<4; ++i) {
02340                                   
02341               xavggen[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].xavggen[i] + xxratio*thePixelTemp_[index].entx[iylow][ihigh].xavggen[i])
02342                           +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].xavggen[i] + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].xavggen[i]);
02343                   
02344               xrmsgen[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].xrmsgen[i] + xxratio*thePixelTemp_[index].entx[iylow][ihigh].xrmsgen[i])
02345                           +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].xrmsgen[i] + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].xrmsgen[i]);            
02346            }
02347            
02348                 lorywidth = thePixelTemp_[index].head.lorywidth;
02349             if(locBz > 0.f) {lorywidth = -lorywidth;}
02350                 lorxwidth = thePixelTemp_[index].head.lorxwidth;
02351                 
02352         
02353         
02354 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
02355         if(qavg <= 0.f || qmin <= 0.f) {
02356                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::qbin, qavg or qmin <= 0," 
02357                 << " Probably someone called the generic pixel reconstruction with an illegal trajectory state" << std::endl;
02358         }
02359 #else
02360         assert(qavg > 0.f && qmin > 0.f);
02361 #endif
02362         
02363 //  Scale the input charge to account for differences between pixelav and CMSSW simulation or data      
02364         
02365         qtotal = qscale*qclus;
02366         
02367 // uncertainty and final corrections depend upon total charge bin          
02368            
02369         fq = qtotal/qavg;
02370         if(fq > 1.5f) {
02371            binq=0;
02372         } else {
02373            if(fq > 1.0f) {
02374               binq=1;
02375            } else {
02376                   if(fq > 0.85f) {
02377                          binq=2;
02378                   } else {
02379                          binq=3;
02380                   }
02381            }
02382         }
02383         
02384 //  Take the errors and bias from the correct charge bin
02385         
02386         sigmay = yrmsgen[binq]; deltay = yavggen[binq];
02387         
02388         sigmax = xrmsgen[binq]; deltax = xavggen[binq];
02389         
02390 // If the charge is too small (then flag it)
02391         
02392         if(qtotal < 0.95f*qmin) {binq = 5;} else {if(qtotal < 0.95f*qmin2) {binq = 4;}}
02393                 
02394     return binq;
02395   
02396 } // qbin
02397 
02398 // ************************************************************************************************************ 
02420 // ************************************************************************************************************ 
02421 int SiPixelTemplate::qbin(int id, float cotalpha, float cotbeta, float locBz, float qclus, float& pixmx, float& sigmay, float& deltay, float& sigmax, float& deltax, 
02422                           float& sy1, float& dy1, float& sy2, float& dy2, float& sx1, float& dx1, float& sx2, float& dx2)
02423 
02424 {
02425         float lorywidth, lorxwidth;
02426         return SiPixelTemplate::qbin(id, cotalpha, cotbeta, locBz, qclus, pixmx, sigmay, deltay, sigmax, deltax, 
02427                                                                  sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, lorywidth, lorxwidth);
02428         
02429 } // qbin
02430 
02431 // ************************************************************************************************************ 
02438 // ************************************************************************************************************ 
02439 int SiPixelTemplate::qbin(int id, float cotalpha, float cotbeta, float qclus)
02440 {
02441         // Interpolate for a new set of track angles 
02442         
02443         // Local variables 
02444         float pixmx, sigmay, deltay, sigmax, deltax, sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, locBz, lorywidth, lorxwidth;
02445         locBz = -1.f;
02446         if(cotbeta < 0.f) {locBz = -locBz;}
02447         return SiPixelTemplate::qbin(id, cotalpha, cotbeta, locBz, qclus, pixmx, sigmay, deltay, sigmax, deltax, 
02448                                                                                   sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, lorywidth, lorxwidth);
02449         
02450 } // qbin
02451 
02452 // ************************************************************************************************************ 
02458 // ************************************************************************************************************ 
02459 int SiPixelTemplate::qbin(int id, float cotbeta, float qclus)
02460 {
02461 // Interpolate for a new set of track angles 
02462                                 
02463 // Local variables 
02464         float pixmx, sigmay, deltay, sigmax, deltax, sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, locBz, lorywidth, lorxwidth;
02465         const float cotalpha = 0.f;
02466         locBz = -1.f;
02467         if(cotbeta < 0.f) {locBz = -locBz;}
02468         return SiPixelTemplate::qbin(id, cotalpha, cotbeta, locBz, qclus, pixmx, sigmay, deltay, sigmax, deltax, 
02469                                                                 sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, lorywidth, lorxwidth);
02470                                 
02471 } // qbin
02472                                 
02473 
02474 
02475 // ************************************************************************************************************ 
02487 // ************************************************************************************************************ 
02488 void SiPixelTemplate::temperrors(int id, float cotalpha, float cotbeta, int qBin, float& sigmay, float& sigmax, float& sy1, float& sy2, float& sx1, float& sx2)
02489 
02490 {
02491     // Interpolate for a new set of track angles 
02492     
02493     // Local variables 
02494     int i;
02495         int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, index;
02496         float yratio, yxratio, xxratio;
02497         float acotb, cotb;
02498         float yrms, xrms;
02499         //bool flip_y;
02500         
02501                 
02502                 // Find the index corresponding to id
02503                 
02504                 index = -1;
02505                 for(i=0; i<(int)thePixelTemp_.size(); ++i) {
02506                         
02507                         if(id == thePixelTemp_[i].head.ID) {
02508                                 
02509                                 index = i;
02510                                 break;
02511                         }
02512             }
02513         
02514 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
02515         if(index < 0 || index >= (int)thePixelTemp_.size()) {
02516                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::temperrors can't find needed template ID = " << id << std::endl;
02517         }
02518 #else
02519         assert(index >= 0 && index < (int)thePixelTemp_.size());
02520 #endif
02521         
02522         
02523 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
02524         if(qBin < 0 || qBin > 5) {
02525                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::temperrors called with illegal qBin = " << qBin << std::endl;
02526         }
02527 #else
02528         assert(qBin >= 0 && qBin < 6);
02529 #endif
02530         
02531 // The error information for qBin > 3 is taken to be the same as qBin=3 
02532         
02533         if(qBin > 3) {qBin = 3;}
02534         //              
02535         
02536         // Interpolate the absolute value of cot(beta)     
02537     
02538     acotb = fabs((double)cotbeta);
02539         cotb = cotbeta;
02540         
02541         // for some cosmics, the ususal gymnastics are incorrect   
02542         
02543 //      if(thePixelTemp_[index].head.Dtype == 0) {
02544                 cotb = acotb;
02545                 //flip_y = false;
02546                 //if(cotbeta < 0.f) {flip_y = true;}
02547 //      } else {
02548 //          if(locBz < 0.f) {
02549 //                      cotb = cotbeta;
02550 //                      flip_y = false;
02551 //              } else {
02552 //                      cotb = -cotbeta;
02553 //                      flip_y = true;
02554 //              }       
02555 //      }
02556                                 
02557                 // Copy the charge scaling factor to the private variable     
02558                 
02559                 Ny = thePixelTemp_[index].head.NTy;
02560                 Nyx = thePixelTemp_[index].head.NTyx;
02561                 Nxx = thePixelTemp_[index].head.NTxx;
02562                 
02563 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
02564                 if(Ny < 2 || Nyx < 1 || Nxx < 2) {
02565                         throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
02566                 }
02567 #else
02568                 assert(Ny > 1 && Nyx > 0 && Nxx > 1);
02569 #endif
02570         
02571                 // next, loop over all y-angle entries   
02572                 
02573                 ilow = 0;
02574                 yratio = 0.f;
02575                 
02576                 if(cotb >= thePixelTemp_[index].enty[Ny-1].cotbeta) {
02577                         
02578                         ilow = Ny-2;
02579                         yratio = 1.f;
02580                         
02581                 } else {
02582                         
02583                         if(cotb >= thePixelTemp_[index].enty[0].cotbeta) {
02584                                 
02585                                 for (i=0; i<Ny-1; ++i) { 
02586                                         
02587                                         if( thePixelTemp_[index].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index].enty[i+1].cotbeta) {
02588                                                 
02589                                                 ilow = i;
02590                                                 yratio = (cotb - thePixelTemp_[index].enty[i].cotbeta)/(thePixelTemp_[index].enty[i+1].cotbeta - thePixelTemp_[index].enty[i].cotbeta);
02591                                                 break;                   
02592                                         }
02593                                 }
02594                         } 
02595                 }
02596                 
02597                 ihigh=ilow + 1;
02598                 
02599                 // Interpolate/store all y-related quantities (flip displacements when flip_y)
02600                 
02601                 sy1 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].syone + yratio*thePixelTemp_[index].enty[ihigh].syone;
02602                 sy2 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].sytwo + yratio*thePixelTemp_[index].enty[ihigh].sytwo;
02603                 yrms=(1.f - yratio)*thePixelTemp_[index].enty[ilow].yrms[qBin] + yratio*thePixelTemp_[index].enty[ihigh].yrms[qBin];
02604                 
02605                 
02606                 // next, loop over all x-angle entries, first, find relevant y-slices   
02607                 
02608                 iylow = 0;
02609                 yxratio = 0.f;
02610                 
02611                 if(acotb >= thePixelTemp_[index].entx[Nyx-1][0].cotbeta) {
02612                         
02613                         iylow = Nyx-2;
02614                         yxratio = 1.f;
02615                         
02616                 } else if(acotb >= thePixelTemp_[index].entx[0][0].cotbeta) {
02617                         
02618                         for (i=0; i<Nyx-1; ++i) { 
02619                                 
02620                                 if( thePixelTemp_[index].entx[i][0].cotbeta <= acotb && acotb < thePixelTemp_[index].entx[i+1][0].cotbeta) {
02621                                         
02622                                         iylow = i;
02623                                         yxratio = (acotb - thePixelTemp_[index].entx[i][0].cotbeta)/(thePixelTemp_[index].entx[i+1][0].cotbeta - thePixelTemp_[index].entx[i][0].cotbeta);
02624                                         break;                   
02625                                 }
02626                         }
02627                 }
02628                 
02629                 iyhigh=iylow + 1;
02630                 
02631                 ilow = 0;
02632                 xxratio = 0.f;
02633                 
02634                 if(cotalpha >= thePixelTemp_[index].entx[0][Nxx-1].cotalpha) {
02635                         
02636                         ilow = Nxx-2;
02637                         xxratio = 1.f;
02638                         
02639                 } else {
02640                         
02641                         if(cotalpha >= thePixelTemp_[index].entx[0][0].cotalpha) {
02642                                 
02643                                 for (i=0; i<Nxx-1; ++i) { 
02644                                         
02645                                         if( thePixelTemp_[index].entx[0][i].cotalpha <= cotalpha && cotalpha < thePixelTemp_[index].entx[0][i+1].cotalpha) {
02646                                                 
02647                                                 ilow = i;
02648                                                 xxratio = (cotalpha - thePixelTemp_[index].entx[0][i].cotalpha)/(thePixelTemp_[index].entx[0][i+1].cotalpha - thePixelTemp_[index].entx[0][i].cotalpha);
02649                                                 break;
02650                                         }
02651                                 }
02652                         } 
02653                 }
02654                 
02655                 ihigh=ilow + 1;
02656                 
02657                 sx1 = (1.f - xxratio)*thePixelTemp_[index].entx[0][ilow].sxone + xxratio*thePixelTemp_[index].entx[0][ihigh].sxone;
02658                 sx2 = (1.f - xxratio)*thePixelTemp_[index].entx[0][ilow].sxtwo + xxratio*thePixelTemp_[index].entx[0][ihigh].sxtwo;
02659                 
02660                 xrms=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].xrms[qBin] + xxratio*thePixelTemp_[index].entx[iylow][ihigh].xrms[qBin])
02661                         +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].xrms[qBin] + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].xrms[qBin]);              
02662                 
02663         
02664         
02665         
02666         //  Take the errors and bias from the correct charge bin
02667         
02668         sigmay = yrms;
02669         
02670         sigmax = xrms;
02671                 
02672     return;
02673         
02674 } // temperrors
02675 
02676 // ************************************************************************************************************ 
02686 // ************************************************************************************************************ 
02687 void SiPixelTemplate::qbin_dist(int id, float cotalpha, float cotbeta, float qbin_frac[4], float& ny1_frac, float& ny2_frac, float& nx1_frac, float& nx2_frac)
02688 
02689 {
02690     // Interpolate for a new set of track angles 
02691     
02692     // Local variables 
02693     int i;
02694         int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, index;
02695         float yratio, yxratio, xxratio;
02696         float acotb, cotb;
02697         float qfrac[4];
02698         //bool flip_y;
02699                         
02700                 // Find the index corresponding to id
02701                 
02702                 index = -1;
02703                 for(i=0; i<(int)thePixelTemp_.size(); ++i) {
02704                         
02705                         if(id == thePixelTemp_[i].head.ID) {
02706                                 
02707                                 index = i;
02708 //                              id_current_ = id;
02709                                 break;
02710                         }
02711             }
02712         
02713 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
02714         if(index < 0 || index >= (int)thePixelTemp_.size()) {
02715                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate::temperrors can't find needed template ID = " << id << std::endl;
02716         }
02717 #else
02718         assert(index >= 0 && index < (int)thePixelTemp_.size());
02719 #endif
02720         
02721         //              
02722         
02723         // Interpolate the absolute value of cot(beta)     
02724     
02725     acotb = fabs((double)cotbeta);
02726         cotb = cotbeta;
02727         
02728         
02729         // for some cosmics, the ususal gymnastics are incorrect   
02730         
02731 //      if(thePixelTemp_[index].head.Dtype == 0) {
02732             cotb = acotb;
02733                 //flip_y = false;
02734                 //if(cotbeta < 0.f) {flip_y = true;}
02735 //      } else {
02736 //          if(locBz < 0.f) {
02737 //                      cotb = cotbeta;
02738 //                      flip_y = false;
02739 //              } else {
02740 //                      cotb = -cotbeta;
02741 //                      flip_y = true;
02742 //              }       
02743 //      }
02744                 
02745                 // Copy the charge scaling factor to the private variable     
02746                 
02747                 Ny = thePixelTemp_[index].head.NTy;
02748                 Nyx = thePixelTemp_[index].head.NTyx;
02749                 Nxx = thePixelTemp_[index].head.NTxx;
02750                 
02751 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
02752                 if(Ny < 2 || Nyx < 1 || Nxx < 2) {
02753                         throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
02754                 }
02755 #else
02756                 assert(Ny > 1 && Nyx > 0 && Nxx > 1);
02757 #endif
02758         
02759                 // next, loop over all y-angle entries   
02760                 
02761                 ilow = 0;
02762                 yratio = 0.f;
02763                 
02764                 if(cotb >= thePixelTemp_[index].enty[Ny-1].cotbeta) {
02765                         
02766                         ilow = Ny-2;
02767                         yratio = 1.f;
02768                         
02769                 } else {
02770                         
02771                         if(cotb >= thePixelTemp_[index].enty[0].cotbeta) {
02772                                 
02773                                 for (i=0; i<Ny-1; ++i) { 
02774                                         
02775                                         if( thePixelTemp_[index].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index].enty[i+1].cotbeta) {
02776                                                 
02777                                                 ilow = i;
02778                                                 yratio = (cotb - thePixelTemp_[index].enty[i].cotbeta)/(thePixelTemp_[index].enty[i+1].cotbeta - thePixelTemp_[index].enty[i].cotbeta);
02779                                                 break;                   
02780                                         }
02781                                 }
02782                         } 
02783                 }
02784                 
02785                 ihigh=ilow + 1;
02786                 
02787                 // Interpolate/store all y-related quantities (flip displacements when flip_y)
02788                 ny1_frac = (1.f - yratio)*thePixelTemp_[index].enty[ilow].fracyone + yratio*thePixelTemp_[index].enty[ihigh].fracyone;
02789             ny2_frac = (1.f - yratio)*thePixelTemp_[index].enty[ilow].fracytwo + yratio*thePixelTemp_[index].enty[ihigh].fracytwo;
02790                 
02791                 // next, loop over all x-angle entries, first, find relevant y-slices   
02792                 
02793                 iylow = 0;
02794                 yxratio = 0.f;
02795                 
02796                 if(acotb >= thePixelTemp_[index].entx[Nyx-1][0].cotbeta) {
02797                         
02798                         iylow = Nyx-2;
02799                         yxratio = 1.f;
02800                         
02801                 } else if(acotb >= thePixelTemp_[index].entx[0][0].cotbeta) {
02802                         
02803                         for (i=0; i<Nyx-1; ++i) { 
02804                                 
02805                                 if( thePixelTemp_[index].entx[i][0].cotbeta <= acotb && acotb < thePixelTemp_[index].entx[i+1][0].cotbeta) {
02806                                         
02807                                         iylow = i;
02808                                         yxratio = (acotb - thePixelTemp_[index].entx[i][0].cotbeta)/(thePixelTemp_[index].entx[i+1][0].cotbeta - thePixelTemp_[index].entx[i][0].cotbeta);
02809                                         break;                   
02810                                 }
02811                         }
02812                 }
02813                 
02814                 iyhigh=iylow + 1;
02815                 
02816                 ilow = 0;
02817                 xxratio = 0.f;
02818                 
02819                 if(cotalpha >= thePixelTemp_[index].entx[0][Nxx-1].cotalpha) {
02820                         
02821                         ilow = Nxx-2;
02822                         xxratio = 1.f;
02823                         
02824                 } else {
02825                         
02826                         if(cotalpha >= thePixelTemp_[index].entx[0][0].cotalpha) {
02827                                 
02828                                 for (i=0; i<Nxx-1; ++i) { 
02829                                         
02830                                         if( thePixelTemp_[index].entx[0][i].cotalpha <= cotalpha && cotalpha < thePixelTemp_[index].entx[0][i+1].cotalpha) {
02831                                                 
02832                                                 ilow = i;
02833                                                 xxratio = (cotalpha - thePixelTemp_[index].entx[0][i].cotalpha)/(thePixelTemp_[index].entx[0][i+1].cotalpha - thePixelTemp_[index].entx[0][i].cotalpha);
02834                                                 break;
02835                                         }
02836                                 }
02837                         } 
02838                 }
02839                 
02840                 ihigh=ilow + 1;
02841                 
02842                 for(i=0; i<3; ++i) {
02843                    qfrac[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].qbfrac[i] + xxratio*thePixelTemp_[index].entx[iylow][ihigh].qbfrac[i])
02844                    +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].qbfrac[i] + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].qbfrac[i]);             
02845                 }
02846                 nx1_frac = (1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].fracxone + xxratio*thePixelTemp_[index].entx[iylow][ihigh].fracxone)
02847                    +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].fracxone + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].fracxone);               
02848                 nx2_frac = (1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].fracxtwo + xxratio*thePixelTemp_[index].entx[iylow][ihigh].fracxtwo)
02849                    +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].fracxtwo + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].fracxtwo);               
02850                 
02851         
02852 
02853         qbin_frac[0] = qfrac[0];
02854         qbin_frac[1] = qbin_frac[0] + qfrac[1];
02855         qbin_frac[2] = qbin_frac[1] + qfrac[2];
02856         qbin_frac[3] = 1.f;
02857     return;
02858         
02859 } // qbin
02860 
02861 // *************************************************************************************************************************************
02863 
02869 // *************************************************************************************************************************************
02870 
02871 bool SiPixelTemplate::simpletemplate2D(float xhit, float yhit, std::vector<bool>& ydouble, std::vector<bool>& xdouble, float template2d[BXM2][BYM2])
02872 {
02873         
02874         // Local variables
02875         
02876         float x0, y0, xf, yf, xi, yi, sf, si, s0, qpix, slopey, slopex, ds;
02877         int i, j, jpix0, ipix0, jpixf, ipixf, jpix, ipix, nx, ny, anx, any, jmax, imax;
02878         float qtotal;
02879         //      double path;
02880         std::list<SimplePixel> list;
02881         std::list<SimplePixel>::iterator listIter, listEnd;     
02882         
02883         // Calculate the entry and exit points for the line charge from the track
02884         
02885         x0 = xhit - 0.5*zsize_*cota_current_;
02886         y0 = yhit - 0.5*zsize_*cotb_current_;
02887         
02888         jpix0 = floor(x0/xsize_)+1;
02889         ipix0 = floor(y0/ysize_)+1;
02890         
02891         if(jpix0 < 0 || jpix0 > BXM3) {return false;}
02892         if(ipix0 < 0 || ipix0 > BYM3) {return false;}
02893         
02894         xf = xhit + 0.5*zsize_*cota_current_ + lorxwidth_;
02895         yf = yhit + 0.5*zsize_*cotb_current_ + lorywidth_;
02896         
02897         jpixf = floor(xf/xsize_)+1;
02898         ipixf = floor(yf/ysize_)+1;
02899         
02900         if(jpixf < 0 || jpixf > BXM3) {return false;}
02901         if(ipixf < 0 || ipixf > BYM3) {return false;}
02902         
02903 // total charge length 
02904         
02905         sf = std::sqrt((xf-x0)*(xf-x0) + (yf-y0)*(yf-y0));
02906         if((xf-x0) != 0.f) {slopey = (yf-y0)/(xf-x0);} else { slopey = 1.e10;}
02907         if((yf-y0) != 0.f) {slopex = (xf-x0)/(yf-y0);} else { slopex = 1.e10;}
02908         
02909 // use average charge in this direction
02910         
02911         qtotal = qavg_avg_;
02912         
02913         SimplePixel element;
02914         element.s = sf;
02915         element.x = xf;
02916         element.y = yf;
02917         element.i = ipixf;
02918         element.j = jpixf;
02919         element.btype = 0;
02920         list.push_back(element);
02921         
02922         //  nx is the number of x interfaces crossed by the line charge 
02923         
02924         nx = jpixf - jpix0;
02925         anx = abs(nx);
02926         if(anx > 0) {
02927                 if(nx > 0) {
02928                         for(j=jpix0; j<jpixf; ++j) {
02929                                 xi = xsize_*j;
02930                                 yi = slopey*(xi-x0) + y0;
02931                                 ipix = (int)(yi/ysize_)+1;
02932                                 si = std::sqrt((xi-x0)*(xi-x0) + (yi-y0)*(yi-y0));
02933                                 element.s = si;
02934                                 element.x = xi;
02935                                 element.y = yi;
02936                                 element.i = ipix;
02937                                 element.j = j;
02938                                 element.btype = 1;
02939                                 list.push_back(element);
02940                         }
02941                 } else {
02942                         for(j=jpix0; j>jpixf; --j) {
02943                                 xi = xsize_*(j-1);
02944                                 yi = slopey*(xi-x0) + y0;
02945                                 ipix = (int)(yi/ysize_)+1;
02946                                 si = std::sqrt((xi-x0)*(xi-x0) + (yi-y0)*(yi-y0));
02947                                 element.s = si;
02948                                 element.x = xi;
02949                                 element.y = yi;
02950                                 element.i = ipix;
02951                                 element.j = j;
02952                                 element.btype = 1;
02953                                 list.push_back(element);
02954                         }
02955                 }
02956         }
02957         
02958         ny = ipixf - ipix0;
02959         any = abs(ny);
02960         if(any > 0) {
02961                 if(ny > 0) {
02962                         for(i=ipix0; i<ipixf; ++i) {
02963                                 yi = ysize_*i;
02964                                 xi = slopex*(yi-y0) + x0;
02965                                 jpix = (int)(xi/xsize_)+1;
02966                                 si = std::sqrt((xi-x0)*(xi-x0) + (yi-y0)*(yi-y0));
02967                                 element.s = si;
02968                                 element.x = xi;
02969                                 element.y = yi;
02970                                 element.i = i;
02971                                 element.j = jpix;
02972                                 element.btype = 2;
02973                                 list.push_back(element);
02974                         }
02975                 } else {
02976                         for(i=ipix0; i>ipixf; --i) {
02977                                 yi = ysize_*(i-1);
02978                                 xi = slopex*(yi-y0) + x0;
02979                                 jpix = (int)(xi/xsize_)+1;
02980                                 si = std::sqrt((xi-x0)*(xi-x0) + (yi-y0)*(yi-y0));
02981                                 element.s = si;
02982                                 element.x = xi;
02983                                 element.y = yi;
02984                                 element.i = i;
02985                                 element.j = jpix;
02986                                 element.btype = 2;
02987                                 list.push_back(element);
02988                         }
02989                 }
02990         }
02991         
02992         imax = std::max(ipix0, ipixf);
02993         jmax = std::max(jpix0, jpixf);
02994         
02995         // Sort the list according to the distance from the initial point
02996         
02997         list.sort();
02998         
02999         // Look for double pixels and adjust the list appropriately
03000         
03001         for(i=1; i<imax; ++i) {
03002                 if(ydouble[i-1]) {
03003                         listIter = list.begin();
03004                         if(ny > 0) {
03005                                 while(listIter != list.end()) {
03006                                         if(listIter->i == i && listIter->btype == 2) {
03007                                                 listIter = list.erase(listIter);
03008                                                 continue;
03009                                         }
03010                                         if(listIter->i > i) {
03011                                                 --(listIter->i);
03012                                         }
03013                                         ++listIter;
03014                                 }
03015                         } else {
03016                                 while(listIter != list.end()) {
03017                                         if(listIter->i == i+1 && listIter->btype == 2) {
03018                                                 listIter = list.erase(listIter);
03019                                                 continue;
03020                                         }
03021                                         if(listIter->i > i+1) {
03022                                                 --(listIter->i);
03023                                         }
03024                                         ++listIter;
03025                                 }                               
03026                         }
03027                 }
03028         }
03029         
03030         for(j=1; j<jmax; ++j) {
03031                 if(xdouble[j-1]) {
03032                         listIter = list.begin();
03033                         if(nx > 0) {
03034                                 while(listIter != list.end()) {
03035                                         if(listIter->j == j && listIter->btype == 1) {
03036                                                 listIter = list.erase(listIter);
03037                                                 continue;
03038                                         }
03039                                         if(listIter->j > j) {
03040                                                 --(listIter->j);
03041                                         }
03042                                         ++listIter;
03043                                 }
03044                         } else {
03045                                 while(listIter != list.end()) {
03046                                         if(listIter->j == j+1 && listIter->btype == 1) {
03047                                                 listIter = list.erase(listIter);
03048                                                 continue;
03049                                         }
03050                                         if(listIter->j > j+1) {
03051                                                 --(listIter->j);
03052                                         }
03053                                         ++listIter;
03054                                 }                               
03055                         }
03056                 }
03057         }
03058         
03059         // The list now contains the path lengths of the line charge in each pixel from (x0,y0).  Cacluate the lengths of the segments and the charge. 
03060         
03061         s0 = 0.f;
03062         listIter = list.begin();
03063         listEnd = list.end();
03064         for( ;listIter != listEnd; ++listIter) {
03065                 si = listIter->s;
03066                 ds = si - s0;
03067                 s0 = si;
03068                 j = listIter->j;
03069                 i = listIter->i;
03070                 if(sf > 0.f) { qpix = qtotal*ds/sf;} else {qpix = qtotal;}
03071                 template2d[j][i] += qpix;
03072         }
03073         
03074         return true;
03075         
03076 }  // simpletemplate2D
03077 
03078 
03079 // ************************************************************************************************************ 
03084 // ************************************************************************************************************ 
03085 void SiPixelTemplate::vavilov_pars(double& mpv, double& sigma, double& kappa)
03086 
03087 {
03088         // Local variables 
03089         int i;
03090         int ilow, ihigh, Ny;
03091         float yratio, cotb, cotalpha0, arg;
03092         
03093 // Interpolate in cotbeta only for the correct total path length (converts cotalpha, cotbeta into an effective cotbeta) 
03094         
03095         cotalpha0 =  thePixelTemp_[index_id_].enty[0].cotalpha;
03096     arg = cotb_current_*cotb_current_ + cota_current_*cota_current_ - cotalpha0*cotalpha0;
03097     if(arg < 0.f) arg = 0.f;
03098         cotb = std::sqrt(arg);
03099                 
03100 // Copy the charge scaling factor to the private variable     
03101         
03102         Ny = thePixelTemp_[index_id_].head.NTy;
03103         
03104 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
03105         if(Ny < 2) {
03106                 throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny = " << Ny << std::endl;
03107         }
03108 #else
03109         assert(Ny > 1);
03110 #endif
03111         
03112 // next, loop over all y-angle entries   
03113         
03114         ilow = 0;
03115         yratio = 0.f;
03116         
03117         if(cotb >= thePixelTemp_[index_id_].enty[Ny-1].cotbeta) {
03118                 
03119                 ilow = Ny-2;
03120                 yratio = 1.f;
03121                 
03122         } else {
03123                 
03124                 if(cotb >= thePixelTemp_[index_id_].enty[0].cotbeta) {
03125                         
03126                         for (i=0; i<Ny-1; ++i) { 
03127                                 
03128                                 if( thePixelTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index_id_].enty[i+1].cotbeta) {
03129                                         
03130                                         ilow = i;
03131                                         yratio = (cotb - thePixelTemp_[index_id_].enty[i].cotbeta)/(thePixelTemp_[index_id_].enty[i+1].cotbeta - thePixelTemp_[index_id_].enty[i].cotbeta);
03132                                         break;                   
03133                                 }
03134                         }
03135                 } 
03136         }
03137         
03138         ihigh=ilow + 1;
03139         
03140 // Interpolate Vavilov parameters
03141         
03142         mpvvav_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].mpvvav + yratio*thePixelTemp_[index_id_].enty[ihigh].mpvvav;
03143         sigmavav_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].sigmavav + yratio*thePixelTemp_[index_id_].enty[ihigh].sigmavav;
03144         kappavav_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].kappavav + yratio*thePixelTemp_[index_id_].enty[ihigh].kappavav;
03145         
03146 // Copy to parameter list
03147         
03148         
03149         mpv = (double)mpvvav_;
03150         sigma = (double)sigmavav_;
03151         kappa = (double)kappavav_;
03152         
03153         return;
03154         
03155 } // vavilov_pars
03156 
03157 // ************************************************************************************************************ 
03162 // ************************************************************************************************************ 
03163 void SiPixelTemplate::vavilov2_pars(double& mpv, double& sigma, double& kappa)
03164 
03165 {
03166         // Local variables 
03167         int i;
03168         int ilow, ihigh, Ny;
03169         float yratio, cotb, cotalpha0, arg;
03170         
03171     // Interpolate in cotbeta only for the correct total path length (converts cotalpha, cotbeta into an effective cotbeta) 
03172         
03173         cotalpha0 =  thePixelTemp_[index_id_].enty[0].cotalpha;
03174     arg = cotb_current_*cotb_current_ + cota_current_*cota_current_ - cotalpha0*cotalpha0;
03175     if(arg < 0.f) arg = 0.f;
03176         cotb = std::sqrt(arg);
03177         
03178         // Copy the charge scaling factor to the private variable     
03179         
03180         Ny = thePixelTemp_[index_id_].head.NTy;
03181         
03182 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
03183         if(Ny < 2) {
03184                 throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny = " << Ny << std::endl;
03185         }
03186 #else
03187         assert(Ny > 1);
03188 #endif
03189         
03190         // next, loop over all y-angle entries   
03191         
03192         ilow = 0;
03193         yratio = 0.f;
03194         
03195         if(cotb >= thePixelTemp_[index_id_].enty[Ny-1].cotbeta) {
03196                 
03197                 ilow = Ny-2;
03198                 yratio = 1.f;
03199                 
03200         } else {
03201                 
03202                 if(cotb >= thePixelTemp_[index_id_].enty[0].cotbeta) {
03203                         
03204                         for (i=0; i<Ny-1; ++i) { 
03205                                 
03206                                 if( thePixelTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index_id_].enty[i+1].cotbeta) {
03207                                         
03208                                         ilow = i;
03209                                         yratio = (cotb - thePixelTemp_[index_id_].enty[i].cotbeta)/(thePixelTemp_[index_id_].enty[i+1].cotbeta - thePixelTemp_[index_id_].enty[i].cotbeta);
03210                                         break;                   
03211                                 }
03212                         }
03213                 } 
03214         }
03215         
03216         ihigh=ilow + 1;
03217         
03218         // Interpolate Vavilov parameters
03219         
03220         mpvvav2_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].mpvvav2 + yratio*thePixelTemp_[index_id_].enty[ihigh].mpvvav2;
03221         sigmavav2_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].sigmavav2 + yratio*thePixelTemp_[index_id_].enty[ihigh].sigmavav2;
03222         kappavav2_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].kappavav2 + yratio*thePixelTemp_[index_id_].enty[ihigh].kappavav2;
03223         
03224         // Copy to parameter list
03225                 
03226         mpv = (double)mpvvav2_;
03227         sigma = (double)sigmavav2_;
03228         kappa = (double)kappavav2_;
03229         
03230         return;
03231         
03232 } // vavilov2_pars
03233 
03234
/data/refman/pasoursint/CMSSW_6_1_2_SLHC2_patch1/src/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplate.cc
