CMSSW Reference Manual

00001 //
00002 //  SiPixelTemplate.cc  Version 5.00 
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 //  Add bias info for Barrel and FPix separately in the header
00026 //  Improve the charge estimation for larger cot(alpha) tracks
00027 //  Change interpolate method to return false boolean if track angles are outside of range
00028 //  Add template info and method for truncation information
00029 //  Change to allow template sizes to be changed at compile time
00030 //
00031 //  Created by Morris Swartz on 10/27/06.
00032 //  Copyright 2006 __TheJohnsHopkinsUniversity__. All rights reserved.
00033 //
00034 //
00035 
00036 //#include <stdlib.h> 
00037 //#include <stdio.h>
00038 #include <math.h>
00039 #include <algorithm>
00040 #include <vector>
00041 //#include "boost/multi_array.hpp"
00042 #include <iostream>
00043 #include <iomanip>
00044 #include <sstream>
00045 #include <fstream>
00046 
00047 
00048 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00049 #include "CondFormats/SiPixelObjects/interface/SiPixelTemplate.h"
00050 #include "FWCore/ParameterSet/interface/FileInPath.h"
00051 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00052 #define LOGERROR(x) LogError(x)
00053 #define LOGINFO(x) LogInfo(x)
00054 #define ENDL " "
00055 using namespace edm;
00056 #else
00057 #include "SiPixelTemplate.h"
00058 #define LOGERROR(x) std::cout << x << ": "
00059 #define LOGINFO(x) std::cout << x << ": "
00060 #define ENDL std::endl
00061 #endif
00062 
00063 //**************************************************************** 
00068 //**************************************************************** 
00069 bool SiPixelTemplate::pushfile(int filenum)
00070 {
00071     // Add template stored in external file numbered filenum to theTemplateStore
00072     
00073     // Local variables 
00074     int i, j, k, l;
00075         const char *tempfile;
00076         char title[80];
00077     char c;
00078         const int code_version={11};
00079         
00080 
00081 
00082 //  Create a filename for this run 
00083 
00084         //std::ostringstream tout;
00085         //tout << "template_summary_zp" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
00086         //std::string tempf = tout.str();
00087         //tempfile = tempf.c_str();
00088         
00089  std::ostringstream tout;
00090  tout << "RecoLocalTracker/SiPixelRecHits/data/template_summary_zp" 
00091       << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
00092  std::string tempf = tout.str();
00093  
00094  // std::cout << "tempf = " << tempf << std::endl;
00095  edm::FileInPath file( tempf.c_str() );
00096  tempfile = (file.fullPath()).c_str();
00097  // std::cout << "tempfile = " << tempfile << std::endl;
00098 
00099 
00100 //  open the template file 
00101 
00102  std::ifstream in_file(tempfile, std::ios::in);
00103  
00104  if(in_file.is_open()) {
00105  
00106  // Create a local template storage entry
00107         
00108         SiPixelTemplateStore theCurrentTemp;
00109         
00110 // Read-in a header string first and print it    
00111     
00112     for (i=0; (c=in_file.get()) != '\n'; ++i) {
00113        if(i < 79) {theCurrentTemp.head.title[i] = c;}
00114     }
00115         if(i > 78) {i=78;}
00116         theCurrentTemp.head.title[i+1] ='\0';
00117     LOGINFO("SiPixelTemplate") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
00118     
00119 // next, the header information     
00120     
00121     in_file >> theCurrentTemp.head.ID >> theCurrentTemp.head.NBy >> theCurrentTemp.head.NByx >> theCurrentTemp.head.NBxx
00122                 >> theCurrentTemp.head.NFy >> theCurrentTemp.head.NFyx >> theCurrentTemp.head.NFxx >> theCurrentTemp.head.Bbias 
00123                          >> theCurrentTemp.head.Fbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
00124                         >> theCurrentTemp.head.s50 >> theCurrentTemp.head.templ_version;
00125                         
00126         if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file, no template load" << ENDL; return false;}
00127         
00128     LOGINFO("SiPixelTemplate") << "Template ID = " << theCurrentTemp.head.ID << ", NBy = " << theCurrentTemp.head.NBy << ", NByx = " << theCurrentTemp.head.NByx 
00129                  << ", NBxx = " << theCurrentTemp.head.NBxx << ", NFy = " << theCurrentTemp.head.NFy << ", NFyx = " << theCurrentTemp.head.NFyx
00130                  << ", NFxx = " << theCurrentTemp.head.NFxx << ", Barrel bias voltage " << theCurrentTemp.head.Bbias << ", FPix bias voltage " << theCurrentTemp.head.Fbias << ", temperature "
00131                  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
00132                  << ", 1/2 threshold " << theCurrentTemp.head.s50 << ", Template Version " << theCurrentTemp.head.templ_version << ENDL;    
00133                         
00134         if(theCurrentTemp.head.templ_version != code_version) {LOGERROR("SiPixelTemplate") << "code expects version " << code_version << ", no template load" << ENDL; return false;}
00135                  
00136 // next, loop over all barrel y-angle entries   
00137 
00138     for (i=0; i < theCurrentTemp.head.NBy; ++i) {     
00139     
00140        in_file >> theCurrentTemp.entby[i].runnum >> theCurrentTemp.entby[i].costrk[0] 
00141                    >> theCurrentTemp.entby[i].costrk[1] >> theCurrentTemp.entby[i].costrk[2]; 
00142                         
00143            if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 1, no template load, run # " << theCurrentTemp.entby[i].runnum << ENDL; return false;}
00144                           
00145 // Calculate the alpha, beta, and cot(beta) for this entry 
00146 
00147        theCurrentTemp.entby[i].alpha = static_cast<float>(atan2((double)theCurrentTemp.entby[i].costrk[2], (double)theCurrentTemp.entby[i].costrk[0]));
00148            
00149            theCurrentTemp.entby[i].cotalpha = theCurrentTemp.entby[i].costrk[0]/theCurrentTemp.entby[i].costrk[2];
00150 
00151        theCurrentTemp.entby[i].beta = static_cast<float>(atan2((double)theCurrentTemp.entby[i].costrk[2], (double)theCurrentTemp.entby[i].costrk[1]));
00152            
00153            theCurrentTemp.entby[i].cotbeta = theCurrentTemp.entby[i].costrk[1]/theCurrentTemp.entby[i].costrk[2];
00154     
00155        in_file >> theCurrentTemp.entby[i].qavg >> theCurrentTemp.entby[i].pixmax >> theCurrentTemp.entby[i].symax >> theCurrentTemp.entby[i].dyone
00156                    >> theCurrentTemp.entby[i].syone >> theCurrentTemp.entby[i].sxmax >> theCurrentTemp.entby[i].dxone >> theCurrentTemp.entby[i].sxone;
00157                         
00158        if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 2, no template load, run # " << theCurrentTemp.entby[i].runnum << ENDL; return false;}
00159     
00160        in_file >> theCurrentTemp.entby[i].dytwo >> theCurrentTemp.entby[i].sytwo >> theCurrentTemp.entby[i].dxtwo 
00161                    >> theCurrentTemp.entby[i].sxtwo >> theCurrentTemp.entby[i].qmin >> theCurrentTemp.entby[i].clsleny >> theCurrentTemp.entby[i].clslenx;
00162                         
00163        if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 3, no template load, run # " << theCurrentTemp.entby[i].runnum << ENDL; return false;}
00164                           
00165            for (j=0; j<2; ++j) {
00166     
00167           in_file >> theCurrentTemp.entby[i].ypar[j][0] >> theCurrentTemp.entby[i].ypar[j][1] 
00168                       >> theCurrentTemp.entby[i].ypar[j][2] >> theCurrentTemp.entby[i].ypar[j][3] >> theCurrentTemp.entby[i].ypar[j][4];
00169                         
00170           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 4, no template load, run # " << theCurrentTemp.entby[i].runnum << ENDL; return false;}
00171                           
00172            }
00173                           
00174            for (j=0; j<9; ++j) {
00175     
00176           for (k=0; k<TYSIZE; ++k) {in_file >> theCurrentTemp.entby[i].ytemp[j][k];}
00177                         
00178           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 5, no template load, run # " << theCurrentTemp.entby[i].runnum << ENDL; return false;}
00179            }
00180                           
00181            for (j=0; j<2; ++j) {
00182     
00183                   in_file >> theCurrentTemp.entby[i].xpar[j][0] >> theCurrentTemp.entby[i].xpar[j][1] 
00184                       >> theCurrentTemp.entby[i].xpar[j][2] >> theCurrentTemp.entby[i].xpar[j][3] >> theCurrentTemp.entby[i].xpar[j][4];
00185                         
00186           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 6, no template load, run # " << theCurrentTemp.entby[i].runnum << ENDL; return false;}
00187                           
00188            }
00189                           
00190            for (j=0; j<9; ++j) {
00191     
00192               for (k=0; k<TXSIZE; ++k) {in_file >> theCurrentTemp.entby[i].xtemp[j][k];} 
00193                                         
00194           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 7, no template load, run # " << theCurrentTemp.entby[i].runnum << ENDL; return false;}
00195            }
00196            
00197            for (j=0; j<4; ++j) {
00198     
00199           in_file >> theCurrentTemp.entby[i].yavg[j] >> theCurrentTemp.entby[i].yrms[j] >> theCurrentTemp.entby[i].ygx0[j] >> theCurrentTemp.entby[i].ygsig[j];
00200                         
00201           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 8, no template load, run # " << theCurrentTemp.entby[i].runnum << ENDL; return false;}
00202            }
00203                                   
00204            for (j=0; j<4; ++j) {
00205     
00206           in_file >> theCurrentTemp.entby[i].yflpar[j][0] >> theCurrentTemp.entby[i].yflpar[j][1] >> theCurrentTemp.entby[i].yflpar[j][2] 
00207                                   >> theCurrentTemp.entby[i].yflpar[j][3] >> theCurrentTemp.entby[i].yflpar[j][4] >> theCurrentTemp.entby[i].yflpar[j][5];
00208                         
00209           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 9, no template load, run # " << theCurrentTemp.entby[i].runnum << ENDL; return false;}
00210            }
00211            
00212                    for (j=0; j<4; ++j) {
00213     
00214           in_file >> theCurrentTemp.entby[i].xavg[j] >> theCurrentTemp.entby[i].xrms[j] >> theCurrentTemp.entby[i].xgx0[j] >> theCurrentTemp.entby[i].xgsig[j];
00215                         
00216           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 10, no template load, run # " << theCurrentTemp.entby[i].runnum << ENDL; return false;}
00217            }
00218                           
00219            for (j=0; j<4; ++j) {
00220     
00221           in_file >> theCurrentTemp.entby[i].xflpar[j][0] >> theCurrentTemp.entby[i].xflpar[j][1] >> theCurrentTemp.entby[i].xflpar[j][2] 
00222                           >> theCurrentTemp.entby[i].xflpar[j][3] >> theCurrentTemp.entby[i].xflpar[j][4] >> theCurrentTemp.entby[i].xflpar[j][5];
00223                         
00224           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 11, no template load, run # " << theCurrentTemp.entby[i].runnum << ENDL; return false;}
00225            }
00226                           
00227            for (j=0; j<4; ++j) {
00228     
00229           in_file >> theCurrentTemp.entby[i].chi2yavg[j] >> theCurrentTemp.entby[i].chi2ymin[j] >> theCurrentTemp.entby[i].chi2xavg[j] >> theCurrentTemp.entby[i].chi2xmin[j];
00230 
00231           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 12, no template load, run # " << theCurrentTemp.entby[i].runnum << ENDL; return false;}
00232            }
00233            
00234            for (j=0; j<4; ++j) {
00235     
00236           in_file >> theCurrentTemp.entby[i].yavgc2m[j] >> theCurrentTemp.entby[i].yrmsc2m[j] >> theCurrentTemp.entby[i].ygx0c2m[j] >> theCurrentTemp.entby[i].ygsigc2m[j];
00237                         
00238           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 13, no template load, run # " << theCurrentTemp.entby[i].runnum << ENDL; return false;}
00239            }
00240            
00241                    for (j=0; j<4; ++j) {
00242     
00243           in_file >> theCurrentTemp.entby[i].xavgc2m[j] >> theCurrentTemp.entby[i].xrmsc2m[j] >> theCurrentTemp.entby[i].xgx0c2m[j] >> theCurrentTemp.entby[i].xgsigc2m[j];
00244                         
00245           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14, no template load, run # " << theCurrentTemp.entby[i].runnum << ENDL; return false;}
00246            } 
00247            
00248            in_file >> theCurrentTemp.entby[i].yspare[0] >> theCurrentTemp.entby[i].yspare[1] >> theCurrentTemp.entby[i].yspare[2] >> theCurrentTemp.entby[i].yspare[3] >> theCurrentTemp.entby[i].yspare[4]
00249             >> theCurrentTemp.entby[i].yspare[5] >> theCurrentTemp.entby[i].yspare[6] >> theCurrentTemp.entby[i].yspare[7] >> theCurrentTemp.entby[i].yspare[8] >> theCurrentTemp.entby[i].yspare[9];
00250 
00251            if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 15, no template load, run # " << theCurrentTemp.entby[i].runnum << ENDL; return false;}
00252 
00253            in_file >> theCurrentTemp.entby[i].xspare[0] >> theCurrentTemp.entby[i].xspare[1] >> theCurrentTemp.entby[i].xspare[2] >> theCurrentTemp.entby[i].xspare[3] >> theCurrentTemp.entby[i].xspare[4]
00254             >> theCurrentTemp.entby[i].xspare[5] >> theCurrentTemp.entby[i].xspare[6] >> theCurrentTemp.entby[i].xspare[7] >> theCurrentTemp.entby[i].xspare[8] >> theCurrentTemp.entby[i].xspare[9];
00255 
00256            if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 16, no template load, run # " << theCurrentTemp.entby[i].runnum << ENDL; return false;}
00257            
00258         }
00259         
00260 // next, loop over all barrel x-angle entries   
00261 
00262   for (k=0; k < theCurrentTemp.head.NByx; ++k) { 
00263 
00264     for (i=0; i < theCurrentTemp.head.NBxx; ++i) { 
00265         
00266        in_file >> theCurrentTemp.entbx[k][i].runnum >> theCurrentTemp.entbx[k][i].costrk[0] 
00267                    >> theCurrentTemp.entbx[k][i].costrk[1] >> theCurrentTemp.entbx[k][i].costrk[2]; 
00268                         
00269        if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 17, no template load, run # " << theCurrentTemp.entbx[k][i].runnum << ENDL; return false;}
00270                           
00271 // Calculate the alpha, beta, and cot(beta) for this entry 
00272 
00273        theCurrentTemp.entbx[k][i].alpha = static_cast<float>(atan2((double)theCurrentTemp.entbx[k][i].costrk[2], (double)theCurrentTemp.entbx[k][i].costrk[0]));
00274            
00275            theCurrentTemp.entbx[k][i].cotalpha = theCurrentTemp.entbx[k][i].costrk[0]/theCurrentTemp.entbx[k][i].costrk[2];
00276 
00277        theCurrentTemp.entbx[k][i].beta = static_cast<float>(atan2((double)theCurrentTemp.entbx[k][i].costrk[2], (double)theCurrentTemp.entbx[k][i].costrk[1]));
00278            
00279            theCurrentTemp.entbx[k][i].cotbeta = theCurrentTemp.entbx[k][i].costrk[1]/theCurrentTemp.entbx[k][i].costrk[2];
00280     
00281        in_file >> theCurrentTemp.entbx[k][i].qavg >> theCurrentTemp.entbx[k][i].pixmax >> theCurrentTemp.entbx[k][i].symax >> theCurrentTemp.entbx[k][i].dyone
00282                    >> theCurrentTemp.entbx[k][i].syone >> theCurrentTemp.entbx[k][i].sxmax >> theCurrentTemp.entbx[k][i].dxone >> theCurrentTemp.entbx[k][i].sxone;
00283                         
00284        if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 18, no template load, run # " << theCurrentTemp.entbx[k][i].runnum << ENDL; return false;}
00285     
00286        in_file >> theCurrentTemp.entbx[k][i].dytwo >> theCurrentTemp.entbx[k][i].sytwo >> theCurrentTemp.entbx[k][i].dxtwo 
00287                    >> theCurrentTemp.entbx[k][i].sxtwo >> theCurrentTemp.entbx[k][i].qmin >> theCurrentTemp.entbx[k][i].clsleny >> theCurrentTemp.entbx[k][i].clslenx;
00288                         
00289        if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 19, no template load, run # " << theCurrentTemp.entbx[k][i].runnum << ENDL; return false;}
00290                           
00291            for (j=0; j<2; ++j) {
00292     
00293           in_file >> theCurrentTemp.entbx[k][i].ypar[j][0] >> theCurrentTemp.entbx[k][i].ypar[j][1] 
00294                       >> theCurrentTemp.entbx[k][i].ypar[j][2] >> theCurrentTemp.entbx[k][i].ypar[j][3] >> theCurrentTemp.entbx[k][i].ypar[j][4];
00295                                                 
00296           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 20, no template load, run # " << theCurrentTemp.entbx[k][i].runnum << ENDL; return false;}
00297            }
00298                           
00299            for (j=0; j<9; ++j) {
00300     
00301               for (l=0; l<TYSIZE; ++l) {in_file >> theCurrentTemp.entbx[k][i].ytemp[j][l];} 
00302                                         
00303                   if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 21, no template load, run # " << theCurrentTemp.entbx[k][i].runnum << ENDL; return false;}
00304            }
00305                           
00306            for (j=0; j<2; ++j) {
00307     
00308                   in_file >> theCurrentTemp.entbx[k][i].xpar[j][0] >> theCurrentTemp.entbx[k][i].xpar[j][1] 
00309                       >> theCurrentTemp.entbx[k][i].xpar[j][2] >> theCurrentTemp.entbx[k][i].xpar[j][3] >> theCurrentTemp.entbx[k][i].xpar[j][4];
00310                           
00311                         
00312           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 22, no template load, run # " << theCurrentTemp.entbx[k][i].runnum << ENDL; return false;}
00313            }
00314                           
00315            for (j=0; j<9; ++j) {
00316     
00317           for (l=0; l<TXSIZE; ++l) {in_file >> theCurrentTemp.entbx[k][i].xtemp[j][l];} 
00318                                         
00319           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 23, no template load, run # " << theCurrentTemp.entbx[k][i].runnum << ENDL; return false;}
00320            }
00321            
00322            for (j=0; j<4; ++j) {
00323     
00324           in_file >> theCurrentTemp.entbx[k][i].yavg[j] >> theCurrentTemp.entbx[k][i].yrms[j] >> theCurrentTemp.entbx[k][i].ygx0[j] >> theCurrentTemp.entbx[k][i].ygsig[j];
00325                         
00326           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 24, no template load, run # " << theCurrentTemp.entbx[k][i].runnum << ENDL; return false;}
00327            }
00328                                   
00329            for (j=0; j<4; ++j) {
00330     
00331           in_file >> theCurrentTemp.entbx[k][i].yflpar[j][0] >> theCurrentTemp.entbx[k][i].yflpar[j][1] >> theCurrentTemp.entbx[k][i].yflpar[j][2] 
00332                                   >> theCurrentTemp.entbx[k][i].yflpar[j][3] >> theCurrentTemp.entbx[k][i].yflpar[j][4] >> theCurrentTemp.entbx[k][i].yflpar[j][5];
00333                         
00334           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 25, no template load, run # " << theCurrentTemp.entbx[k][i].runnum << ENDL; return false;}
00335            }
00336                                   
00337            for (j=0; j<4; ++j) {
00338     
00339           in_file >> theCurrentTemp.entbx[k][i].xavg[j] >> theCurrentTemp.entbx[k][i].xrms[j] >> theCurrentTemp.entbx[k][i].xgx0[j] >> theCurrentTemp.entbx[k][i].xgsig[j];
00340                         
00341           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 26, no template load, run # " << theCurrentTemp.entbx[k][i].runnum << ENDL; return false;}
00342            }
00343                           
00344            for (j=0; j<4; ++j) {
00345     
00346           in_file >> theCurrentTemp.entbx[k][i].xflpar[j][0] >> theCurrentTemp.entbx[k][i].xflpar[j][1] >> theCurrentTemp.entbx[k][i].xflpar[j][2] 
00347                           >> theCurrentTemp.entbx[k][i].xflpar[j][3] >> theCurrentTemp.entbx[k][i].xflpar[j][4] >> theCurrentTemp.entbx[k][i].xflpar[j][5];
00348                         
00349           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 27, no template load, run # " << theCurrentTemp.entbx[k][i].runnum << ENDL; return false;}
00350            }
00351                           
00352            for (j=0; j<4; ++j) {
00353     
00354           in_file >> theCurrentTemp.entbx[k][i].chi2yavg[j] >> theCurrentTemp.entbx[k][i].chi2ymin[j] >> theCurrentTemp.entbx[k][i].chi2xavg[j] >> theCurrentTemp.entbx[k][i].chi2xmin[j];
00355                         
00356           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 28, no template load, run # " << theCurrentTemp.entbx[k][i].runnum << ENDL; return false;}
00357            }
00358            
00359            for (j=0; j<4; ++j) {
00360     
00361           in_file >> theCurrentTemp.entbx[k][i].yavgc2m[j] >> theCurrentTemp.entbx[k][i].yrmsc2m[j] >> theCurrentTemp.entbx[k][i].ygx0c2m[j] >> theCurrentTemp.entbx[k][i].ygsigc2m[j];
00362                         
00363           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 29, no template load, run # " << theCurrentTemp.entbx[k][i].runnum << ENDL; return false;}
00364            }
00365            
00366                    for (j=0; j<4; ++j) {
00367     
00368           in_file >> theCurrentTemp.entbx[k][i].xavgc2m[j] >> theCurrentTemp.entbx[k][i].xrmsc2m[j] >> theCurrentTemp.entbx[k][i].xgx0c2m[j] >> theCurrentTemp.entbx[k][i].xgsigc2m[j];
00369                         
00370           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30, no template load, run # " << theCurrentTemp.entbx[k][i].runnum << ENDL; return false;}
00371            }
00372            
00373            in_file >> theCurrentTemp.entbx[k][i].yspare[0] >> theCurrentTemp.entbx[k][i].yspare[1] >> theCurrentTemp.entbx[k][i].yspare[2] >> theCurrentTemp.entbx[k][i].yspare[3] >> theCurrentTemp.entbx[k][i].yspare[4]
00374             >> theCurrentTemp.entbx[k][i].yspare[5] >> theCurrentTemp.entbx[k][i].yspare[6] >> theCurrentTemp.entbx[k][i].yspare[7] >> theCurrentTemp.entbx[k][i].yspare[8] >> theCurrentTemp.entbx[k][i].yspare[9];
00375                         
00376            if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 31, no template load, run # " << theCurrentTemp.entbx[k][i].runnum << ENDL; return false;}
00377 
00378            in_file >> theCurrentTemp.entbx[k][i].xspare[0] >> theCurrentTemp.entbx[k][i].xspare[1] >> theCurrentTemp.entbx[k][i].xspare[2] >> theCurrentTemp.entbx[k][i].xspare[3] >> theCurrentTemp.entbx[k][i].xspare[4]
00379             >> theCurrentTemp.entbx[k][i].xspare[5] >> theCurrentTemp.entbx[k][i].xspare[6] >> theCurrentTemp.entbx[k][i].xspare[7] >> theCurrentTemp.entbx[k][i].xspare[8] >> theCurrentTemp.entbx[k][i].xspare[9];
00380                         
00381            if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 32, no template load, run # " << theCurrentTemp.entbx[k][i].runnum << ENDL; return false;}
00382            
00383         }
00384   }     
00385     
00386 // next, loop over all forward y-angle entries   
00387 
00388     for (i=0; i < theCurrentTemp.head.NFy; ++i) {     
00389     
00390        in_file >> theCurrentTemp.entfy[i].runnum >> theCurrentTemp.entfy[i].costrk[0] 
00391                    >> theCurrentTemp.entfy[i].costrk[1] >> theCurrentTemp.entfy[i].costrk[2]; 
00392                         
00393        if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 33, no template load, run # " << theCurrentTemp.entfy[i].runnum << ENDL; return false;}
00394                           
00395 // Calculate the alpha, beta, and cot(beta) for this entry 
00396 
00397        theCurrentTemp.entfy[i].alpha = static_cast<float>(atan2((double)theCurrentTemp.entfy[i].costrk[2], (double)theCurrentTemp.entfy[i].costrk[0]));
00398            
00399            theCurrentTemp.entfy[i].cotalpha = theCurrentTemp.entfy[i].costrk[0]/theCurrentTemp.entfy[i].costrk[2];
00400 
00401        theCurrentTemp.entfy[i].beta = static_cast<float>(atan2((double)theCurrentTemp.entfy[i].costrk[2], (double)theCurrentTemp.entfy[i].costrk[1]));
00402            
00403            theCurrentTemp.entfy[i].cotbeta = theCurrentTemp.entfy[i].costrk[1]/theCurrentTemp.entfy[i].costrk[2];
00404     
00405        in_file >> theCurrentTemp.entfy[i].qavg >> theCurrentTemp.entfy[i].pixmax >> theCurrentTemp.entfy[i].symax >> theCurrentTemp.entfy[i].dyone
00406                    >> theCurrentTemp.entfy[i].syone >> theCurrentTemp.entfy[i].sxmax >> theCurrentTemp.entfy[i].dxone >> theCurrentTemp.entfy[i].sxone;
00407                         
00408        if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 34, no template load, run # " << theCurrentTemp.entfy[i].runnum << ENDL; return false;}
00409            
00410        in_file >> theCurrentTemp.entfy[i].dytwo >> theCurrentTemp.entfy[i].sytwo >> theCurrentTemp.entfy[i].dxtwo 
00411                    >> theCurrentTemp.entfy[i].sxtwo >> theCurrentTemp.entfy[i].qmin >> theCurrentTemp.entfy[i].clsleny >> theCurrentTemp.entfy[i].clslenx;
00412                         
00413        if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 35, no template load, run # " << theCurrentTemp.entfy[i].runnum << ENDL; return false;}
00414                           
00415            for (j=0; j<2; ++j) {
00416     
00417           in_file >> theCurrentTemp.entfy[i].ypar[j][0] >> theCurrentTemp.entfy[i].ypar[j][1] 
00418                       >> theCurrentTemp.entfy[i].ypar[j][2] >> theCurrentTemp.entfy[i].ypar[j][3] >> theCurrentTemp.entfy[i].ypar[j][4];
00419                         
00420           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 36, no template load, run # " << theCurrentTemp.entfy[i].runnum << ENDL; return false;}
00421                           
00422            }
00423                           
00424            for (j=0; j<9; ++j) {
00425     
00426           for (l=0; l<TYSIZE; ++l) {in_file >> theCurrentTemp.entfy[i].ytemp[j][l];} 
00427                                         
00428           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 37, no template load, run # " << theCurrentTemp.entfy[i].runnum << ENDL; return false;}
00429            }
00430                           
00431            for (j=0; j<2; ++j) {
00432     
00433                   in_file >> theCurrentTemp.entfy[i].xpar[j][0] >> theCurrentTemp.entfy[i].xpar[j][1] 
00434                       >> theCurrentTemp.entfy[i].xpar[j][2] >> theCurrentTemp.entfy[i].xpar[j][3] >> theCurrentTemp.entfy[i].xpar[j][4];
00435                           
00436                         
00437           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 38, no template load, run # " << theCurrentTemp.entfy[i].runnum << ENDL; return false;}
00438            }
00439                           
00440            for (j=0; j<9; ++j) {
00441     
00442           for (l=0; l<TXSIZE; ++l) {in_file >> theCurrentTemp.entfy[i].xtemp[j][l];} 
00443                                         
00444           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 39, no template load, run # " << theCurrentTemp.entfy[i].runnum << ENDL; return false;}
00445            }
00446            
00447            for (j=0; j<4; ++j) {
00448     
00449           in_file >> theCurrentTemp.entfy[i].yavg[j] >> theCurrentTemp.entfy[i].yrms[j] >> theCurrentTemp.entfy[i].ygx0[j] >> theCurrentTemp.entfy[i].ygsig[j];
00450                         
00451           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 40, no template load, run # " << theCurrentTemp.entfy[i].runnum << ENDL; return false;}
00452            }
00453                                   
00454            for (j=0; j<4; ++j) {
00455     
00456           in_file >> theCurrentTemp.entfy[i].yflpar[j][0] >> theCurrentTemp.entfy[i].yflpar[j][1] >> theCurrentTemp.entfy[i].yflpar[j][2]
00457                                   >> theCurrentTemp.entfy[i].yflpar[j][3] >> theCurrentTemp.entfy[i].yflpar[j][4] >> theCurrentTemp.entfy[i].yflpar[j][5];
00458                         
00459           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 41, no template load, run # " << theCurrentTemp.entfy[i].runnum << ENDL; return false;}
00460            }
00461            
00462            for (j=0; j<4; ++j) {
00463     
00464           in_file >> theCurrentTemp.entfy[i].xavg[j] >> theCurrentTemp.entfy[i].xrms[j] >> theCurrentTemp.entfy[i].xgx0[j] >> theCurrentTemp.entfy[i].xgsig[j];
00465                         
00466           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 42, no template load, run # " << theCurrentTemp.entfy[i].runnum << ENDL; return false;}
00467            }
00468                           
00469            for (j=0; j<4; ++j) {
00470     
00471           in_file >> theCurrentTemp.entfy[i].xflpar[j][0] >> theCurrentTemp.entfy[i].xflpar[j][1] >> theCurrentTemp.entfy[i].xflpar[j][2] 
00472                           >> theCurrentTemp.entfy[i].xflpar[j][3] >> theCurrentTemp.entfy[i].xflpar[j][4] >> theCurrentTemp.entfy[i].xflpar[j][5];
00473                         
00474           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 43, no template load, run # " << theCurrentTemp.entfy[i].runnum << ENDL; return false;}
00475            }
00476                           
00477            for (j=0; j<4; ++j) {
00478     
00479           in_file >> theCurrentTemp.entfy[i].chi2yavg[j] >> theCurrentTemp.entfy[i].chi2ymin[j] >> theCurrentTemp.entfy[i].chi2xavg[j] >> theCurrentTemp.entfy[i].chi2xmin[j];
00480                         
00481           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 44, no template load, run # " << theCurrentTemp.entfy[i].runnum << ENDL; return false;}
00482            }
00483            
00484            for (j=0; j<4; ++j) {
00485     
00486           in_file >> theCurrentTemp.entfy[i].yavgc2m[j] >> theCurrentTemp.entfy[i].yrmsc2m[j] >> theCurrentTemp.entfy[i].ygx0c2m[j] >> theCurrentTemp.entfy[i].ygsigc2m[j];
00487                         
00488           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 45, no template load, run # " << theCurrentTemp.entfy[i].runnum << ENDL; return false;}
00489            }
00490            
00491                    for (j=0; j<4; ++j) {
00492     
00493           in_file >> theCurrentTemp.entfy[i].xavgc2m[j] >> theCurrentTemp.entfy[i].xrmsc2m[j] >> theCurrentTemp.entfy[i].xgx0c2m[j] >> theCurrentTemp.entfy[i].xgsigc2m[j];
00494                         
00495           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 46, no template load, run # " << theCurrentTemp.entfy[i].runnum << ENDL; return false;}
00496            }
00497            
00498            in_file >> theCurrentTemp.entfy[i].yspare[0] >> theCurrentTemp.entfy[i].yspare[1] >> theCurrentTemp.entfy[i].yspare[2] >> theCurrentTemp.entfy[i].yspare[3] >> theCurrentTemp.entfy[i].yspare[4]
00499             >> theCurrentTemp.entfy[i].yspare[5] >> theCurrentTemp.entfy[i].yspare[6] >> theCurrentTemp.entfy[i].yspare[7] >> theCurrentTemp.entfy[i].yspare[8] >> theCurrentTemp.entfy[i].yspare[9];
00500                         
00501            if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 47, no template load, run # " << theCurrentTemp.entfy[i].runnum << ENDL; return false;}
00502 
00503            in_file >> theCurrentTemp.entfy[i].xspare[0] >> theCurrentTemp.entfy[i].xspare[1] >> theCurrentTemp.entfy[i].xspare[2] >> theCurrentTemp.entfy[i].xspare[3] >> theCurrentTemp.entfy[i].xspare[4]
00504             >> theCurrentTemp.entfy[i].xspare[5] >> theCurrentTemp.entfy[i].xspare[6] >> theCurrentTemp.entfy[i].xspare[7] >> theCurrentTemp.entfy[i].xspare[8] >> theCurrentTemp.entfy[i].xspare[9];
00505                         
00506            if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 48, no template load, run # " << theCurrentTemp.entfy[i].runnum << ENDL; return false;}
00507            
00508         }
00509         
00510 // next, loop over all forward x-angle entries   
00511 
00512   for (k=0; k < theCurrentTemp.head.NFyx; ++k) { 
00513   
00514     for (i=0; i < theCurrentTemp.head.NFxx; ++i) {     
00515     
00516        in_file >> theCurrentTemp.entfx[k][i].runnum >> theCurrentTemp.entfx[k][i].costrk[0] 
00517                    >> theCurrentTemp.entfx[k][i].costrk[1] >> theCurrentTemp.entfx[k][i].costrk[2]; 
00518                         
00519        if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 49, no template load, run # " << theCurrentTemp.entfx[k][i].runnum << ENDL; return false;}
00520                            
00521 // Calculate the alpha, beta, and cot(beta) for this entry 
00522 
00523        theCurrentTemp.entfx[k][i].alpha = static_cast<float>(atan2((double)theCurrentTemp.entfx[k][i].costrk[2], (double)theCurrentTemp.entfx[k][i].costrk[0]));
00524            
00525            theCurrentTemp.entfx[k][i].cotalpha = theCurrentTemp.entfx[k][i].costrk[0]/theCurrentTemp.entfx[k][i].costrk[2];
00526 
00527        theCurrentTemp.entfx[k][i].beta = static_cast<float>(atan2((double)theCurrentTemp.entfx[k][i].costrk[2], (double)theCurrentTemp.entfx[k][i].costrk[1]));
00528            
00529            theCurrentTemp.entfx[k][i].cotbeta = theCurrentTemp.entfx[k][i].costrk[1]/theCurrentTemp.entfx[k][i].costrk[2];
00530     
00531        in_file >> theCurrentTemp.entfx[k][i].qavg >> theCurrentTemp.entfx[k][i].pixmax >> theCurrentTemp.entfx[k][i].symax >> theCurrentTemp.entfx[k][i].dyone
00532                    >> theCurrentTemp.entfx[k][i].syone >> theCurrentTemp.entfx[k][i].sxmax >> theCurrentTemp.entfx[k][i].dxone >> theCurrentTemp.entfx[k][i].sxone;
00533                         
00534        if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 50, no template load, run # " << theCurrentTemp.entfx[k][i].runnum << ENDL; return false;}
00535     
00536        in_file >> theCurrentTemp.entfx[k][i].dytwo >> theCurrentTemp.entfx[k][i].sytwo >> theCurrentTemp.entfx[k][i].dxtwo 
00537                    >> theCurrentTemp.entfx[k][i].sxtwo >> theCurrentTemp.entfx[k][i].qmin >> theCurrentTemp.entfx[k][i].clsleny >> theCurrentTemp.entfx[k][i].clslenx;
00538                         
00539        if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 51, no template load, run # " << theCurrentTemp.entfx[k][i].runnum << ENDL; return false;}
00540                           
00541            for (j=0; j<2; ++j) {
00542     
00543           in_file >> theCurrentTemp.entfx[k][i].ypar[j][0] >> theCurrentTemp.entfx[k][i].ypar[j][1] 
00544                       >> theCurrentTemp.entfx[k][i].ypar[j][2] >> theCurrentTemp.entfx[k][i].ypar[j][3] >> theCurrentTemp.entfx[k][i].ypar[j][4];
00545                         
00546           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 52, no template load, run # " << theCurrentTemp.entfx[k][i].runnum << ENDL; return false;}
00547                           
00548            }
00549                           
00550            for (j=0; j<9; ++j) {
00551     
00552           for (l=0; l<TYSIZE; ++l) {in_file >> theCurrentTemp.entfx[k][i].ytemp[j][l];} 
00553                                         
00554           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 53, no template load, run # " << theCurrentTemp.entfx[k][i].runnum << ENDL; return false;}
00555            }
00556                           
00557            for (j=0; j<2; ++j) {
00558     
00559                   in_file >> theCurrentTemp.entfx[k][i].xpar[j][0] >> theCurrentTemp.entfx[k][i].xpar[j][1] 
00560                       >> theCurrentTemp.entfx[k][i].xpar[j][2] >> theCurrentTemp.entfx[k][i].xpar[j][3] >> theCurrentTemp.entfx[k][i].xpar[j][4];
00561                           
00562                         
00563           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 54, no template load, run # " << theCurrentTemp.entfx[k][i].runnum << ENDL; return false;}
00564            }
00565                           
00566            for (j=0; j<9; ++j) {
00567     
00568           for (l=0; l<TXSIZE; ++l) {in_file >> theCurrentTemp.entfx[k][i].xtemp[j][l];} 
00569                                         
00570           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 55, no template load, run # " << theCurrentTemp.entfx[k][i].runnum << ENDL; return false;}
00571            }
00572            
00573            for (j=0; j<4; ++j) {
00574     
00575           in_file >> theCurrentTemp.entfx[k][i].yavg[j] >> theCurrentTemp.entfx[k][i].yrms[j] >> theCurrentTemp.entfx[k][i].ygx0[j] >> theCurrentTemp.entfx[k][i].ygsig[j];
00576                         
00577           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 56, no template load, run # " << theCurrentTemp.entfx[k][i].runnum << ENDL; return false;}
00578            }
00579                                   
00580            for (j=0; j<4; ++j) {
00581     
00582           in_file >> theCurrentTemp.entfx[k][i].yflpar[j][0] >> theCurrentTemp.entfx[k][i].yflpar[j][1] >> theCurrentTemp.entfx[k][i].yflpar[j][2] 
00583                           >> theCurrentTemp.entfx[k][i].yflpar[j][3] >> theCurrentTemp.entfx[k][i].yflpar[j][4] >> theCurrentTemp.entfx[k][i].yflpar[j][5];
00584                         
00585           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 57, no template load, run # " << theCurrentTemp.entfx[k][i].runnum << ENDL; return false;}
00586            }
00587            
00588            for (j=0; j<4; ++j) {
00589     
00590           in_file >> theCurrentTemp.entfx[k][i].xavg[j] >> theCurrentTemp.entfx[k][i].xrms[j] >> theCurrentTemp.entfx[k][i].xgx0[j] >> theCurrentTemp.entfx[k][i].xgsig[j];
00591                         
00592           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 58, no template load, run # " << theCurrentTemp.entfx[k][i].runnum << ENDL; return false;}
00593            }
00594                           
00595            for (j=0; j<4; ++j) {
00596     
00597           in_file >> theCurrentTemp.entfx[k][i].xflpar[j][0] >> theCurrentTemp.entfx[k][i].xflpar[j][1] >> theCurrentTemp.entfx[k][i].xflpar[j][2] 
00598                           >> theCurrentTemp.entfx[k][i].xflpar[j][3] >> theCurrentTemp.entfx[k][i].xflpar[j][4] >> theCurrentTemp.entfx[k][i].xflpar[j][5];
00599                         
00600           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 59, no template load, run # " << theCurrentTemp.entfx[k][i].runnum << ENDL; return false;}
00601            }
00602                           
00603            for (j=0; j<4; ++j) {
00604     
00605           in_file >> theCurrentTemp.entfx[k][i].chi2yavg[j] >> theCurrentTemp.entfx[k][i].chi2ymin[j] >> theCurrentTemp.entfx[k][i].chi2xavg[j] >> theCurrentTemp.entfx[k][i].chi2xmin[j];
00606 
00607           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 60, no template load, run # " << theCurrentTemp.entfx[k][i].runnum << ENDL; return false;}
00608            }
00609            for (j=0; j<4; ++j) {
00610     
00611           in_file >> theCurrentTemp.entfx[k][i].yavgc2m[j] >> theCurrentTemp.entfx[k][i].yrmsc2m[j] >> theCurrentTemp.entfx[k][i].ygx0c2m[j] >> theCurrentTemp.entfx[k][i].ygsigc2m[j];
00612                         
00613           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 61, no template load, run # " << theCurrentTemp.entfx[k][i].runnum << ENDL; return false;}
00614            }
00615            
00616                    for (j=0; j<4; ++j) {
00617     
00618           in_file >> theCurrentTemp.entfx[k][i].xavgc2m[j] >> theCurrentTemp.entfx[k][i].xrmsc2m[j] >> theCurrentTemp.entfx[k][i].xgx0c2m[j] >> theCurrentTemp.entfx[k][i].xgsigc2m[j];
00619                         
00620           if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 62, no template load, run # " << theCurrentTemp.entfx[k][i].runnum << ENDL; return false;}
00621            }
00622            
00623            in_file >> theCurrentTemp.entfx[k][i].yspare[0] >> theCurrentTemp.entfx[k][i].yspare[1] >> theCurrentTemp.entfx[k][i].yspare[2] >> theCurrentTemp.entfx[k][i].yspare[3] >> theCurrentTemp.entfx[k][i].yspare[4]
00624             >> theCurrentTemp.entfx[k][i].yspare[5] >> theCurrentTemp.entfx[k][i].yspare[6] >> theCurrentTemp.entfx[k][i].yspare[7] >> theCurrentTemp.entfx[k][i].yspare[8] >> theCurrentTemp.entfx[k][i].yspare[9];
00625 
00626            if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 63, no template load, run # " << theCurrentTemp.entfx[k][i].runnum << ENDL; return false;}
00627 
00628            in_file >> theCurrentTemp.entfx[k][i].xspare[0] >> theCurrentTemp.entfx[k][i].xspare[1] >> theCurrentTemp.entfx[k][i].xspare[2] >> theCurrentTemp.entfx[k][i].xspare[3] >> theCurrentTemp.entfx[k][i].xspare[4]
00629             >> theCurrentTemp.entfx[k][i].xspare[5] >> theCurrentTemp.entfx[k][i].xspare[6] >> theCurrentTemp.entfx[k][i].xspare[7] >> theCurrentTemp.entfx[k][i].xspare[8] >> theCurrentTemp.entfx[k][i].xspare[9];
00630 
00631            if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 64, no template load, run # " << theCurrentTemp.entfx[k][i].runnum << ENDL; return false;}
00632            
00633         }       
00634   }
00635     
00636     in_file.close();
00637         
00638 // Add this template to the store
00639         
00640         thePixelTemp.push_back(theCurrentTemp);
00641         
00642         return true;
00643         
00644  } else {
00645  
00646  // If file didn't open, report this
00647  
00648     LOGERROR("SiPixelTemplate") << "Error opening File" << tempfile << ENDL;
00649         return false;
00650         
00651  }
00652         
00653 } // TempInit 
00654 
00655 
00656 
00657 
00658 
00659 
00660 
00661 
00662 
00663 
00664 
00665 // ************************************************************************************************************ 
00672 // ************************************************************************************************************ 
00673 bool SiPixelTemplate::interpolate(int id, bool fpix, float cotalpha, float cotbeta)
00674 {
00675     // Interpolate for a new set of track angles 
00676     
00677     // Local variables 
00678     int i, j, ind;
00679         int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, imidy, imaxx;
00680         float yratio, yxratio, xxratio, sxmax, qcorrect, symax;
00681         bool success;
00682 //      std::vector <float> xrms(4), xgsig(4), xrmsc2m(4), xgsigc2m(4);
00683         std::vector <float> chi2xavg(4), chi2xmin(4);
00684 
00685 
00686 // Check to see if interpolation is valid     
00687 
00688 if(id != id_current || fpix != fpix_current || cotalpha != cota_current || cotbeta != cotb_current) {
00689 
00690         fpix_current = fpix; cota_current = cotalpha; cotb_current = cotbeta;
00691         
00692         if(id != id_current) {
00693 
00694 // Find the index corresponding to id
00695 
00696        index_id = -1;
00697        for(i=0; i<thePixelTemp.size(); ++i) {
00698         
00699               if(id == thePixelTemp[i].head.ID) {
00700            
00701                  index_id = i;
00702                      id_current = id;
00703                      break;
00704           }
00705             }
00706      }
00707          
00708          assert(index_id >= 0 && index_id < thePixelTemp.size());
00709          
00710 // Interpolate the absolute value of cot(beta)     
00711     
00712     abs_cotb = fabs((double)cotbeta);
00713         
00714 //      qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)       
00715 
00716     qcorrect=(float)sqrt((double)((1.+cotbeta*cotbeta+cotalpha*cotalpha)/(1.+cotbeta*cotbeta)));
00717 
00718 // Copy the charge scaling factor to the private variable     
00719     
00720     pqscale = thePixelTemp[index_id].head.qscale;
00721 
00722 // Copy the pseudopixel signal size to the private variable     
00723     
00724     ps50 = thePixelTemp[index_id].head.s50;
00725 
00726 // success flags whether or not the track angles are inside the interpolation range
00727     
00728     success = true;
00729         
00730 // Decide which template (FPix or BPix) to use 
00731 
00732     if(fpix) {
00733     
00734 // Begin FPix section, make the index counters easier to use     
00735     
00736        Ny = thePixelTemp[index_id].head.NFy;
00737        Nyx = thePixelTemp[index_id].head.NFyx;
00738        Nxx = thePixelTemp[index_id].head.NFxx;
00739            imaxx = Nyx - 1;
00740            imidy = Nxx/2;
00741         
00742 // next, loop over all y-angle entries   
00743 
00744            ilow = 0;
00745            yratio = 0.;
00746 
00747            if(abs_cotb >= thePixelTemp[index_id].entfy[Ny-1].cotbeta) {
00748         
00749                ilow = Ny-2;
00750                    yratio = 1.;
00751                    success = false;
00752                 
00753            } else {
00754            
00755               if(abs_cotb >= thePixelTemp[index_id].entfy[0].cotbeta) {
00756 
00757              for (i=0; i<Ny-1; ++i) { 
00758     
00759                 if( thePixelTemp[index_id].entfy[i].cotbeta <= abs_cotb && abs_cotb < thePixelTemp[index_id].entfy[i+1].cotbeta) {
00760                   
00761                        ilow = i;
00762                            yratio = (abs_cotb - thePixelTemp[index_id].entfy[i].cotbeta)/(thePixelTemp[index_id].entfy[i+1].cotbeta - thePixelTemp[index_id].entfy[i].cotbeta);
00763                            break;                        
00764                         }
00765                  }
00766                   } else { success = false; }
00767            }
00768         
00769            ihigh=ilow + 1;
00770                           
00771 // Interpolate/store all y-related quantities (flip displacements when cotbeta < 0)
00772 
00773        pyratio = yratio;
00774            pqavg = (1. - yratio)*thePixelTemp[index_id].entfy[ilow].qavg + yratio*thePixelTemp[index_id].entfy[ihigh].qavg;
00775            pqavg *= qcorrect;
00776            symax = (1. - yratio)*thePixelTemp[index_id].entfy[ilow].symax + yratio*thePixelTemp[index_id].entfy[ihigh].symax;
00777            psyparmax = symax;
00778            sxmax = (1. - yratio)*thePixelTemp[index_id].entfy[ilow].sxmax + yratio*thePixelTemp[index_id].entfy[ihigh].sxmax;
00779            pdyone = (1. - yratio)*thePixelTemp[index_id].entfy[ilow].dyone + yratio*thePixelTemp[index_id].entfy[ihigh].dyone;
00780            if(cotbeta < 0.) {pdyone = -pdyone;}
00781            psyone = (1. - yratio)*thePixelTemp[index_id].entfy[ilow].syone + yratio*thePixelTemp[index_id].entfy[ihigh].syone;
00782            pdytwo = (1. - yratio)*thePixelTemp[index_id].entfy[ilow].dytwo + yratio*thePixelTemp[index_id].entfy[ihigh].dytwo;
00783            if(cotbeta < 0.) {pdytwo = -pdytwo;}
00784            psytwo = (1. - yratio)*thePixelTemp[index_id].entfy[ilow].sytwo + yratio*thePixelTemp[index_id].entfy[ihigh].sytwo;
00785            pqmin = (1. - yratio)*thePixelTemp[index_id].entfy[ilow].qmin + yratio*thePixelTemp[index_id].entfy[ihigh].qmin;
00786            pclsleny = fminf(thePixelTemp[index_id].entfy[ilow].clsleny, thePixelTemp[index_id].entfy[ihigh].clsleny);
00787            for(i=0; i<2 ; ++i) {
00788               for(j=0; j<5 ; ++j) {
00789 // Charge loss switches sides when cot(beta) changes sign
00790                      if(cotbeta < 0) {
00791                     pyparl[1-i][j] = thePixelTemp[index_id].entfy[ilow].ypar[i][j];
00792                     pyparh[1-i][j] = thePixelTemp[index_id].entfy[ihigh].ypar[i][j];
00793                          } else {
00794                     pyparl[i][j] = thePixelTemp[index_id].entfy[ilow].ypar[i][j];
00795                     pyparh[i][j] = thePixelTemp[index_id].entfy[ihigh].ypar[i][j];
00796                          }
00797                  pxparly0[i][j] = thePixelTemp[index_id].entfy[ilow].xpar[i][j];
00798                  pxparhy0[i][j] = thePixelTemp[index_id].entfy[ihigh].xpar[i][j];
00799               }
00800            }
00801            for(i=0; i<4; ++i) {
00802               pyavg[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].yavg[i] + yratio*thePixelTemp[index_id].entfy[ihigh].yavg[i];
00803               if(cotbeta < 0.) {pyavg[i] = -pyavg[i];}
00804               pyrms[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].yrms[i] + yratio*thePixelTemp[index_id].entfy[ihigh].yrms[i];
00805 //            pygx0[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].ygx0[i] + yratio*thePixelTemp[index_id].entfy[ihigh].ygx0[i];
00806 //            if(cotbeta < 0.) {pygx0[i] = -pygx0[i];}
00807 //            pygsig[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].ygsig[i] + yratio*thePixelTemp[index_id].entfy[ihigh].ygsig[i];
00808 //            xrms[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].xrms[i] + yratio*thePixelTemp[index_id].entfy[ihigh].xrms[i];
00809 //            xgsig[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].xgsig[i] + yratio*thePixelTemp[index_id].entfy[ihigh].xgsig[i];
00810               pchi2yavg[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].chi2yavg[i] + yratio*thePixelTemp[index_id].entfy[ihigh].chi2yavg[i];
00811               pchi2ymin[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].chi2ymin[i] + yratio*thePixelTemp[index_id].entfy[ihigh].chi2ymin[i];
00812               chi2xavg[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].chi2xavg[i] + yratio*thePixelTemp[index_id].entfy[ihigh].chi2xavg[i];
00813               chi2xmin[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].chi2xmin[i] + yratio*thePixelTemp[index_id].entfy[ihigh].chi2xmin[i];
00814               pyavgc2m[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].yavgc2m[i] + yratio*thePixelTemp[index_id].entfy[ihigh].yavgc2m[i];
00815               if(cotbeta < 0.) {pyavgc2m[i] = -pyavgc2m[i];}
00816               pyrmsc2m[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].yrmsc2m[i] + yratio*thePixelTemp[index_id].entfy[ihigh].yrmsc2m[i];
00817 //            pygx0c2m[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].ygx0c2m[i] + yratio*thePixelTemp[index_id].entfy[ihigh].ygx0c2m[i];
00818 //            if(cotbeta < 0.) {pygx0c2m[i] = -pygx0c2m[i];}
00819 //            pygsigc2m[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].ygsigc2m[i] + yratio*thePixelTemp[index_id].entfy[ihigh].ygsigc2m[i];
00820 //            xrmsc2m[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].xrmsc2m[i] + yratio*thePixelTemp[index_id].entfy[ihigh].xrmsc2m[i];
00821 //            xgsigc2m[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].xgsigc2m[i] + yratio*thePixelTemp[index_id].entfy[ihigh].xgsigc2m[i];
00822                   for(j=0; j<6 ; ++j) {
00823                          pyflparl[i][j] = thePixelTemp[index_id].entfy[ilow].yflpar[i][j];
00824                          pyflparh[i][j] = thePixelTemp[index_id].entfy[ihigh].yflpar[i][j];
00825                          
00826 // Since Q_fl is odd under cotbeta, it flips qutomatically, change only even terms
00827 
00828                          if(cotbeta < 0. && (j == 0 || j == 2 || j == 4)) {
00829                             pyflparl[i][j] = - pyflparl[i][j];
00830                             pyflparh[i][j] = - pyflparh[i][j];
00831                          }
00832                   }
00833            }
00834            
00836 
00837 //       for(i=0; i<10; ++i) {
00838 //                  pyspare[i]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].yspare[i] + yratio*thePixelTemp[index_id].entfy[ihigh].yspare[i];
00839 //       }
00840                           
00841 // Interpolate and build the y-template 
00842         
00843            for(i=0; i<9; ++i) {
00844           pytemp[i][0] = 0.;
00845           pytemp[i][1] = 0.;
00846               pytemp[i][BYM2] = 0.;
00847               pytemp[i][BYM1] = 0.;
00848               for(j=0; j<TYSIZE; ++j) {
00849                   
00850 // Flip the basic y-template when the cotbeta is negative
00851 
00852                      if(cotbeta < 0.) {
00853                     pytemp[8-i][BYM3-j]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].ytemp[i][j] + yratio*thePixelTemp[index_id].entfy[ihigh].ytemp[i][j];
00854                          } else {
00855                     pytemp[i][j+2]=(1. - yratio)*thePixelTemp[index_id].entfy[ilow].ytemp[i][j] + yratio*thePixelTemp[index_id].entfy[ihigh].ytemp[i][j];
00856                          }
00857               }
00858            }
00859         
00860 // next, loop over all x-angle entries, first, find relevant y-slices   
00861         
00862            iylow = 0;
00863            yxratio = 0.;
00864 
00865            if(abs_cotb >= thePixelTemp[index_id].entfx[Nyx-1][0].cotbeta) {
00866         
00867                iylow = Nyx-2;
00868                    yxratio = 1.;
00869                 
00870            } else if(abs_cotb >= thePixelTemp[index_id].entfx[0][0].cotbeta) {
00871 
00872           for (i=0; i<Nyx-1; ++i) { 
00873     
00874              if( thePixelTemp[index_id].entfx[i][0].cotbeta <= abs_cotb && abs_cotb < thePixelTemp[index_id].entfx[i+1][0].cotbeta) {
00875                   
00876                     iylow = i;
00877                         yxratio = (abs_cotb - thePixelTemp[index_id].entfx[i][0].cotbeta)/(thePixelTemp[index_id].entfx[i+1][0].cotbeta - thePixelTemp[index_id].entfx[i][0].cotbeta);
00878                         break;                   
00879                      }
00880               }
00881            }
00882         
00883            iyhigh=iylow + 1;
00884 
00885            ilow = 0;
00886            xxratio = 0.;
00887 
00888            if(cotalpha >= thePixelTemp[index_id].entfx[0][Nxx-1].cotalpha) {
00889         
00890                ilow = Nxx-2;
00891                    xxratio = 1.;
00892                    success = false;
00893                 
00894            } else {
00895            
00896               if(cotalpha >= thePixelTemp[index_id].entfx[0][0].cotalpha) {
00897 
00898              for (i=0; i<Nxx-1; ++i) { 
00899     
00900                 if( thePixelTemp[index_id].entfx[0][i].cotalpha <= cotalpha && cotalpha < thePixelTemp[index_id].entfx[0][i+1].cotalpha) {
00901                   
00902                        ilow = i;
00903                            xxratio = (cotalpha - thePixelTemp[index_id].entfx[0][i].cotalpha)/(thePixelTemp[index_id].entfx[0][i+1].cotalpha - thePixelTemp[index_id].entfx[0][i].cotalpha);
00904                            break;
00905                             }
00906                      }
00907                   } else { success = false; }
00908            }
00909         
00910            ihigh=ilow + 1;
00911                           
00912 // Interpolate/store all x-related quantities 
00913 
00914        pyxratio = yxratio;
00915        pxxratio = xxratio;              
00916                           
00917 // sxparmax defines the maximum charge for which the parameters xpar are defined (not rescaled by cotbeta) 
00918 
00919            psxparmax = (1. - xxratio)*thePixelTemp[index_id].entfx[imaxx][ilow].sxmax + xxratio*thePixelTemp[index_id].entfx[imaxx][ihigh].sxmax;
00920            psxmax = psxparmax;
00921        if(thePixelTemp[index_id].entfx[imaxx][imidy].sxmax != 0.) {psxmax=psxmax/thePixelTemp[index_id].entfx[imaxx][imidy].sxmax*sxmax;}
00922            psymax = (1. - xxratio)*thePixelTemp[index_id].entfx[imaxx][ilow].symax + xxratio*thePixelTemp[index_id].entfx[imaxx][ihigh].symax;
00923        if(thePixelTemp[index_id].entfx[imaxx][imidy].symax != 0.) {psymax=psymax/thePixelTemp[index_id].entfx[imaxx][imidy].symax*symax;}
00924            pdxone = (1. - xxratio)*thePixelTemp[index_id].entfx[0][ilow].dxone + xxratio*thePixelTemp[index_id].entfx[0][ihigh].dxone;
00925            psxone = (1. - xxratio)*thePixelTemp[index_id].entfx[0][ilow].sxone + xxratio*thePixelTemp[index_id].entfx[0][ihigh].sxone;
00926            pdxtwo = (1. - xxratio)*thePixelTemp[index_id].entfx[0][ilow].dxtwo + xxratio*thePixelTemp[index_id].entfx[0][ihigh].dxtwo;
00927            psxtwo = (1. - xxratio)*thePixelTemp[index_id].entfx[0][ilow].sxtwo + xxratio*thePixelTemp[index_id].entfx[0][ihigh].sxtwo;
00928            pclslenx = fminf(thePixelTemp[index_id].entfx[0][ilow].clslenx, thePixelTemp[index_id].entfx[0][ihigh].clslenx);
00929            for(i=0; i<2 ; ++i) {
00930               for(j=0; j<5 ; ++j) {
00931                  pxpar0[i][j] = thePixelTemp[index_id].entfx[imaxx][imidy].xpar[i][j];
00932                  pxparl[i][j] = thePixelTemp[index_id].entfx[imaxx][ilow].xpar[i][j];
00933                  pxparh[i][j] = thePixelTemp[index_id].entfx[imaxx][ihigh].xpar[i][j];
00934               }
00935            }
00936                           
00937 // pixmax is the maximum allowed pixel charge (used for truncation)
00938 
00939            ppixmax=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entfx[iylow][ilow].pixmax + xxratio*thePixelTemp[index_id].entfx[iylow][ihigh].pixmax)
00940                           +yxratio*((1. - xxratio)*thePixelTemp[index_id].entfx[iyhigh][ilow].pixmax + xxratio*thePixelTemp[index_id].entfx[iyhigh][ihigh].pixmax);
00941                           
00942            for(i=0; i<4; ++i) {
00943               pxavg[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entfx[iylow][ilow].xavg[i] + xxratio*thePixelTemp[index_id].entfx[iylow][ihigh].xavg[i])
00944                           +yxratio*((1. - xxratio)*thePixelTemp[index_id].entfx[iyhigh][ilow].xavg[i] + xxratio*thePixelTemp[index_id].entfx[iyhigh][ihigh].xavg[i]);
00945                   
00946               pxrms[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entfx[iylow][ilow].xrms[i] + xxratio*thePixelTemp[index_id].entfx[iylow][ihigh].xrms[i])
00947                           +yxratio*((1. - xxratio)*thePixelTemp[index_id].entfx[iyhigh][ilow].xrms[i] + xxratio*thePixelTemp[index_id].entfx[iyhigh][ihigh].xrms[i]);
00948                   
00949 //            pxgx0[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entfx[iylow][ilow].xgx0[i] + xxratio*thePixelTemp[index_id].entfx[iylow][ihigh].xgx0[i])
00950 //                        +yxratio*((1. - xxratio)*thePixelTemp[index_id].entfx[iyhigh][ilow].xgx0[i] + xxratio*thePixelTemp[index_id].entfx[iyhigh][ihigh].xgx0[i]);
00951                                                         
00952 //            pxgsig[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entfx[iylow][ilow].xgsig[i] + xxratio*thePixelTemp[index_id].entfx[iylow][ihigh].xgsig[i])
00953 //                        +yxratio*((1. - xxratio)*thePixelTemp[index_id].entfx[iyhigh][ilow].xgsig[i] + xxratio*thePixelTemp[index_id].entfx[iyhigh][ihigh].xgsig[i]);
00954                                   
00955               pxavgc2m[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entfx[iylow][ilow].xavgc2m[i] + xxratio*thePixelTemp[index_id].entfx[iylow][ihigh].xavgc2m[i])
00956                           +yxratio*((1. - xxratio)*thePixelTemp[index_id].entfx[iyhigh][ilow].xavgc2m[i] + xxratio*thePixelTemp[index_id].entfx[iyhigh][ihigh].xavgc2m[i]);
00957                   
00958               pxrmsc2m[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entfx[iylow][ilow].xrmsc2m[i] + xxratio*thePixelTemp[index_id].entfx[iylow][ihigh].xrmsc2m[i])
00959                           +yxratio*((1. - xxratio)*thePixelTemp[index_id].entfx[iyhigh][ilow].xrmsc2m[i] + xxratio*thePixelTemp[index_id].entfx[iyhigh][ihigh].xrmsc2m[i]);
00960                   
00961 //            pxgx0c2m[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entfx[iylow][ilow].xgx0c2m[i] + xxratio*thePixelTemp[index_id].entfx[iylow][ihigh].xgx0c2m[i])
00962 //                        +yxratio*((1. - xxratio)*thePixelTemp[index_id].entfx[iyhigh][ilow].xgx0c2m[i] + xxratio*thePixelTemp[index_id].entfx[iyhigh][ihigh].xgx0c2m[i]);
00963                                                         
00964 //            pxgsigc2m[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entfx[iylow][ilow].xgsigc2m[i] + xxratio*thePixelTemp[index_id].entfx[iylow][ihigh].xgsigc2m[i])
00965 //                        +yxratio*((1. - xxratio)*thePixelTemp[index_id].entfx[iyhigh][ilow].xgsigc2m[i] + xxratio*thePixelTemp[index_id].entfx[iyhigh][ihigh].xgsigc2m[i]);
00966 //
00967 //  Try new interpolation scheme
00968 //                                                                                                                      
00969 //            pchi2xavg[i]=((1. - xxratio)*thePixelTemp[index_id].entfx[imaxx][ilow].chi2xavg[i] + xxratio*thePixelTemp[index_id].entfx[imaxx][ihigh].chi2xavg[i]);
00970 //                if(thePixelTemp[index_id].entfx[imaxx][imidy].chi2xavg[i] != 0.) {pchi2xavg[i]=pchi2xavg[i]/thePixelTemp[index_id].entfx[imaxx][imidy].chi2xavg[i]*chi2xavg[i];}
00971 //                                                      
00972 //            pchi2xmin[i]=((1. - xxratio)*thePixelTemp[index_id].entfx[imaxx][ilow].chi2xmin[i] + xxratio*thePixelTemp[index_id].entfx[imaxx][ihigh].chi2xmin[i]);
00973 //                if(thePixelTemp[index_id].entfx[imaxx][imidy].chi2xmin[i] != 0.) {pchi2xmin[i]=pchi2xmin[i]/thePixelTemp[index_id].entfx[imaxx][imidy].chi2xmin[i]*chi2xmin[i];}
00974 //                
00975               pchi2xavg[i]=((1. - xxratio)*thePixelTemp[index_id].entfx[iyhigh][ilow].chi2xavg[i] + xxratio*thePixelTemp[index_id].entfx[iyhigh][ihigh].chi2xavg[i]);
00976                   if(thePixelTemp[index_id].entfx[iyhigh][imidy].chi2xavg[i] != 0.) {pchi2xavg[i]=pchi2xavg[i]/thePixelTemp[index_id].entfx[iyhigh][imidy].chi2xavg[i]*chi2xavg[i];}
00977                                                         
00978               pchi2xmin[i]=((1. - xxratio)*thePixelTemp[index_id].entfx[iyhigh][ilow].chi2xmin[i] + xxratio*thePixelTemp[index_id].entfx[iyhigh][ihigh].chi2xmin[i]);
00979                   if(thePixelTemp[index_id].entfx[iyhigh][imidy].chi2xmin[i] != 0.) {pchi2xmin[i]=pchi2xmin[i]/thePixelTemp[index_id].entfx[iyhigh][imidy].chi2xmin[i]*chi2xmin[i];}
00980                   
00981               for(j=0; j<6 ; ++j) {
00982                  pxflparll[i][j] = thePixelTemp[index_id].entfx[iylow][ilow].xflpar[i][j];
00983                  pxflparlh[i][j] = thePixelTemp[index_id].entfx[iylow][ihigh].xflpar[i][j];
00984                  pxflparhl[i][j] = thePixelTemp[index_id].entfx[iyhigh][ilow].xflpar[i][j];
00985                  pxflparhh[i][j] = thePixelTemp[index_id].entfx[iyhigh][ihigh].xflpar[i][j];
00986                   }
00987            }
00988            
00989 // Do the spares next
00990 
00991 //       for(i=0; i<10; ++i) {
00992 //            pxspare[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entfx[iylow][ilow].xspare[i] + xxratio*thePixelTemp[index_id].entfx[iylow][ihigh].xspare[i])
00993 //                        +yxratio*((1. - xxratio)*thePixelTemp[index_id].entfx[iyhigh][ilow].xspare[i] + xxratio*thePixelTemp[index_id].entfx[iyhigh][ihigh].xspare[i]);
00994 //       }
00995                           
00996 // Interpolate and build the x-template 
00997         
00998            for(i=0; i<9; ++i) {
00999           pxtemp[i][0] = 0.;
01000           pxtemp[i][1] = 0.;
01001               pxtemp[i][BXM2] = 0.;
01002               pxtemp[i][BXM1] = 0.;
01003               for(j=0; j<TXSIZE; ++j) {
01004                 pxtemp[i][j+2]=(1. - xxratio)*thePixelTemp[index_id].entfx[imaxx][ilow].xtemp[i][j] + xxratio*thePixelTemp[index_id].entfx[imaxx][ihigh].xtemp[i][j];
01005               }
01006            }
01007 
01008         } else {
01009         
01010     
01011 // Begin BPix section, make the index counters easier to use     
01012     
01013        Ny = thePixelTemp[index_id].head.NBy;
01014        Nyx = thePixelTemp[index_id].head.NByx;
01015        Nxx = thePixelTemp[index_id].head.NBxx;
01016            imaxx = Nyx - 1;
01017            imidy = Nxx/2;
01018         
01019 // next, loop over all y-angle entries   
01020 
01021            ilow = 0;
01022            yratio = 0.;
01023 
01024            if(abs_cotb >= thePixelTemp[index_id].entby[Ny-1].cotbeta) {
01025         
01026                ilow = Ny-2;
01027                    yratio = 1.;
01028                    success = false;
01029                 
01030            } else {
01031            
01032               if(abs_cotb >= thePixelTemp[index_id].entby[0].cotbeta) {
01033 
01034              for (i=0; i<Ny-1; ++i) { 
01035     
01036                 if( thePixelTemp[index_id].entby[i].cotbeta <= abs_cotb && abs_cotb < thePixelTemp[index_id].entby[i+1].cotbeta) {
01037                   
01038                        ilow = i;
01039                            yratio = (abs_cotb - thePixelTemp[index_id].entby[i].cotbeta)/(thePixelTemp[index_id].entby[i+1].cotbeta - thePixelTemp[index_id].entby[i].cotbeta);
01040                            break;                        
01041                         }
01042                  }
01043                   } else { success = false;}
01044            }
01045         
01046            ihigh=ilow + 1;
01047                           
01048 // Interpolate/store all y-related quantities (flip displacements when cotbeta < 0)
01049 
01050        pyratio = yratio;
01051            pqavg = (1. - yratio)*thePixelTemp[index_id].entby[ilow].qavg + yratio*thePixelTemp[index_id].entby[ihigh].qavg;
01052            pqavg *= qcorrect;
01053            symax = (1. - yratio)*thePixelTemp[index_id].entby[ilow].symax + yratio*thePixelTemp[index_id].entby[ihigh].symax;
01054            psyparmax = symax;
01055            sxmax = (1. - yratio)*thePixelTemp[index_id].entby[ilow].sxmax + yratio*thePixelTemp[index_id].entby[ihigh].sxmax;
01056            pdyone = (1. - yratio)*thePixelTemp[index_id].entby[ilow].dyone + yratio*thePixelTemp[index_id].entby[ihigh].dyone;
01057            if(cotbeta < 0.) {pdyone = -pdyone;}
01058            psyone = (1. - yratio)*thePixelTemp[index_id].entby[ilow].syone + yratio*thePixelTemp[index_id].entby[ihigh].syone;
01059            pdytwo = (1. - yratio)*thePixelTemp[index_id].entby[ilow].dytwo + yratio*thePixelTemp[index_id].entby[ihigh].dytwo;
01060            if(cotbeta < 0.) {pdytwo = -pdytwo;}
01061            psytwo = (1. - yratio)*thePixelTemp[index_id].entby[ilow].sytwo + yratio*thePixelTemp[index_id].entby[ihigh].sytwo;
01062            pqmin = (1. - yratio)*thePixelTemp[index_id].entby[ilow].qmin + yratio*thePixelTemp[index_id].entby[ihigh].qmin;
01063            pclsleny = fminf(thePixelTemp[index_id].entby[ilow].clsleny, thePixelTemp[index_id].entby[ihigh].clsleny);
01064            for(i=0; i<2 ; ++i) {
01065               for(j=0; j<5 ; ++j) {
01066 // Charge loss switches sides when cot(beta) changes sign
01067                      if(cotbeta < 0) {
01068                     pyparl[1-i][j] = thePixelTemp[index_id].entby[ilow].ypar[i][j];
01069                     pyparh[1-i][j] = thePixelTemp[index_id].entby[ihigh].ypar[i][j];
01070                          } else {
01071                     pyparl[i][j] = thePixelTemp[index_id].entby[ilow].ypar[i][j];
01072                     pyparh[i][j] = thePixelTemp[index_id].entby[ihigh].ypar[i][j];
01073                          }
01074                  pxparly0[i][j] = thePixelTemp[index_id].entby[ilow].xpar[i][j];
01075                  pxparhy0[i][j] = thePixelTemp[index_id].entby[ihigh].xpar[i][j];
01076               }
01077            }
01078            for(i=0; i<4; ++i) {
01079               pyavg[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].yavg[i] + yratio*thePixelTemp[index_id].entby[ihigh].yavg[i];
01080               if(cotbeta < 0.) {pyavg[i] = -pyavg[i];}
01081               pyrms[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].yrms[i] + yratio*thePixelTemp[index_id].entby[ihigh].yrms[i];
01082 //            pygx0[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].ygx0[i] + yratio*thePixelTemp[index_id].entby[ihigh].ygx0[i];
01083 //            if(cotbeta < 0.) {pygx0[i] = -pygx0[i];}
01084 //            pygsig[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].ygsig[i] + yratio*thePixelTemp[index_id].entby[ihigh].ygsig[i];
01085 //            xrms[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].xrms[i] + yratio*thePixelTemp[index_id].entby[ihigh].xrms[i];
01086 //            xgsig[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].xgsig[i] + yratio*thePixelTemp[index_id].entby[ihigh].xgsig[i];
01087               pchi2yavg[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].chi2yavg[i] + yratio*thePixelTemp[index_id].entby[ihigh].chi2yavg[i];
01088               pchi2ymin[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].chi2ymin[i] + yratio*thePixelTemp[index_id].entby[ihigh].chi2ymin[i];
01089               chi2xavg[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].chi2xavg[i] + yratio*thePixelTemp[index_id].entby[ihigh].chi2xavg[i];
01090               chi2xmin[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].chi2xmin[i] + yratio*thePixelTemp[index_id].entby[ihigh].chi2xmin[i];
01091               pyavgc2m[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].yavgc2m[i] + yratio*thePixelTemp[index_id].entby[ihigh].yavgc2m[i];
01092               if(cotbeta < 0.) {pyavgc2m[i] = -pyavgc2m[i];}
01093               pyrmsc2m[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].yrmsc2m[i] + yratio*thePixelTemp[index_id].entby[ihigh].yrmsc2m[i];
01094 //            pygx0c2m[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].ygx0c2m[i] + yratio*thePixelTemp[index_id].entby[ihigh].ygx0c2m[i];
01095 //            if(cotbeta < 0.) {pygx0c2m[i] = -pygx0c2m[i];}
01096 //            pygsigc2m[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].ygsigc2m[i] + yratio*thePixelTemp[index_id].entby[ihigh].ygsigc2m[i];
01097 //            xrmsc2m[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].xrmsc2m[i] + yratio*thePixelTemp[index_id].entby[ihigh].xrmsc2m[i];
01098 //            xgsigc2m[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].xgsigc2m[i] + yratio*thePixelTemp[index_id].entby[ihigh].xgsigc2m[i];
01099                   for(j=0; j<6 ; ++j) {
01100                          pyflparl[i][j] = thePixelTemp[index_id].entby[ilow].yflpar[i][j];
01101                          pyflparh[i][j] = thePixelTemp[index_id].entby[ihigh].yflpar[i][j];
01102                          
01103 // Since Q_fl is odd under cotbeta, it flips qutomatically, change only even terms
01104 
01105                          if(cotbeta < 0. && (j == 0 || j == 2 || j == 4)) {
01106                             pyflparl[i][j] = - pyflparl[i][j];
01107                             pyflparh[i][j] = - pyflparh[i][j];
01108                          }
01109                   }
01110            }
01111            
01112 // Do the spares next
01113 
01114 //       for(i=0; i<10; ++i) {
01115 //                pyspare[i]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].yspare[i] + yratio*thePixelTemp[index_id].entby[ihigh].yspare[i];
01116 //       }
01117                           
01118 // Interpolate and build the y-template 
01119         
01120            for(i=0; i<9; ++i) {
01121           pytemp[i][0] = 0.;
01122           pytemp[i][1] = 0.;
01123               pytemp[i][BYM2] = 0.;
01124               pytemp[i][BYM1] = 0.;
01125               for(j=0; j<TYSIZE; ++j) {
01126                   
01127 // Flip the basic y-template when the cotbeta is negative
01128 
01129                      if(cotbeta < 0.) {
01130                     pytemp[8-i][BYM3-j]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].ytemp[i][j] + yratio*thePixelTemp[index_id].entby[ihigh].ytemp[i][j];
01131                          } else {
01132                     pytemp[i][j+2]=(1. - yratio)*thePixelTemp[index_id].entby[ilow].ytemp[i][j] + yratio*thePixelTemp[index_id].entby[ihigh].ytemp[i][j];
01133                          }
01134               }
01135            }
01136         
01137 // next, loop over all x-angle entries, first, find relevant y-slices   
01138 
01139            iylow = 0;
01140            yxratio = 0.;
01141 
01142            if(abs_cotb >= thePixelTemp[index_id].entbx[Nyx-1][0].cotbeta) {
01143         
01144                iylow = Nyx-2;
01145                    yxratio = 1.;
01146                 
01147            } else if(abs_cotb >= thePixelTemp[index_id].entbx[0][0].cotbeta) {
01148 
01149           for (i=0; i<Nyx-1; ++i) { 
01150     
01151              if( thePixelTemp[index_id].entbx[i][0].cotbeta <= abs_cotb && abs_cotb < thePixelTemp[index_id].entbx[i+1][0].cotbeta) {
01152                   
01153                     iylow = i;
01154                         yxratio = (abs_cotb - thePixelTemp[index_id].entbx[i][0].cotbeta)/(thePixelTemp[index_id].entbx[i+1][0].cotbeta - thePixelTemp[index_id].entbx[i][0].cotbeta);
01155                         break;                   
01156                      }
01157               }
01158            }
01159         
01160            iyhigh=iylow + 1;
01161 
01162            ilow = 0;
01163            xxratio = 0.;
01164 
01165            if(cotalpha >= thePixelTemp[index_id].entbx[0][Nxx-1].cotalpha) {
01166         
01167                ilow = Nxx-2;
01168                    xxratio = 1.;
01169                    success = false;
01170                 
01171            } else {
01172               if(cotalpha >= thePixelTemp[index_id].entbx[0][0].cotalpha) {
01173 
01174              for (i=0; i<Nxx-1; ++i) { 
01175     
01176                 if( thePixelTemp[index_id].entbx[0][i].cotalpha <= cotalpha && cotalpha < thePixelTemp[index_id].entbx[0][i+1].cotalpha) {
01177                   
01178                        ilow = i;
01179                            xxratio = (cotalpha - thePixelTemp[index_id].entbx[0][i].cotalpha)/(thePixelTemp[index_id].entbx[0][i+1].cotalpha - thePixelTemp[index_id].entbx[0][i].cotalpha);
01180                            break;
01181                             }
01182                      }
01183                   } else { success = false; }
01184            }
01185         
01186            ihigh=ilow + 1;
01187                           
01188 // Interpolate/store all x-related quantities 
01189 
01190        pyxratio = yxratio;
01191        pxxratio = xxratio;              
01192                           
01193 // sxparmax defines the maximum charge for which the parameters xpar are defined (not rescaled by cotbeta) 
01194 
01195            psxparmax = (1. - xxratio)*thePixelTemp[index_id].entbx[imaxx][ilow].sxmax + xxratio*thePixelTemp[index_id].entbx[imaxx][ihigh].sxmax;
01196            psxmax = psxparmax;
01197        if(thePixelTemp[index_id].entbx[imaxx][imidy].sxmax != 0.) {psxmax=psxmax/thePixelTemp[index_id].entbx[imaxx][imidy].sxmax*sxmax;}
01198            psymax = (1. - xxratio)*thePixelTemp[index_id].entbx[imaxx][ilow].symax + xxratio*thePixelTemp[index_id].entbx[imaxx][ihigh].symax;
01199        if(thePixelTemp[index_id].entbx[imaxx][imidy].symax != 0.) {psymax=psymax/thePixelTemp[index_id].entbx[imaxx][imidy].symax*symax;}
01200            pdxone = (1. - xxratio)*thePixelTemp[index_id].entbx[0][ilow].dxone + xxratio*thePixelTemp[index_id].entbx[0][ihigh].dxone;
01201            psxone = (1. - xxratio)*thePixelTemp[index_id].entbx[0][ilow].sxone + xxratio*thePixelTemp[index_id].entbx[0][ihigh].sxone;
01202            pdxtwo = (1. - xxratio)*thePixelTemp[index_id].entbx[0][ilow].dxtwo + xxratio*thePixelTemp[index_id].entbx[0][ihigh].dxtwo;
01203            psxtwo = (1. - xxratio)*thePixelTemp[index_id].entbx[0][ilow].sxtwo + xxratio*thePixelTemp[index_id].entbx[0][ihigh].sxtwo;
01204            pclslenx = fminf(thePixelTemp[index_id].entbx[0][ilow].clslenx, thePixelTemp[index_id].entbx[0][ihigh].clslenx);
01205            for(i=0; i<2 ; ++i) {
01206               for(j=0; j<5 ; ++j) {
01207                  pxpar0[i][j] = thePixelTemp[index_id].entbx[imaxx][imidy].xpar[i][j];
01208                  pxparl[i][j] = thePixelTemp[index_id].entbx[imaxx][ilow].xpar[i][j];
01209                  pxparh[i][j] = thePixelTemp[index_id].entbx[imaxx][ihigh].xpar[i][j];
01210               }
01211            }
01212                           
01213 // pixmax is the maximum allowed pixel charge (used for truncation)
01214 
01215            ppixmax=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entbx[iylow][ilow].pixmax + xxratio*thePixelTemp[index_id].entbx[iylow][ihigh].pixmax)
01216                           +yxratio*((1. - xxratio)*thePixelTemp[index_id].entbx[iyhigh][ilow].pixmax + xxratio*thePixelTemp[index_id].entbx[iyhigh][ihigh].pixmax);
01217                           
01218            for(i=0; i<4; ++i) {
01219               pxavg[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entbx[iylow][ilow].xavg[i] + xxratio*thePixelTemp[index_id].entbx[iylow][ihigh].xavg[i])
01220                           +yxratio*((1. - xxratio)*thePixelTemp[index_id].entbx[iyhigh][ilow].xavg[i] + xxratio*thePixelTemp[index_id].entbx[iyhigh][ihigh].xavg[i]);
01221                   
01222               pxrms[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entbx[iylow][ilow].xrms[i] + xxratio*thePixelTemp[index_id].entbx[iylow][ihigh].xrms[i])
01223                           +yxratio*((1. - xxratio)*thePixelTemp[index_id].entbx[iyhigh][ilow].xrms[i] + xxratio*thePixelTemp[index_id].entbx[iyhigh][ihigh].xrms[i]);
01224                   
01225 //            pxgx0[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entbx[iylow][ilow].xgx0[i] + xxratio*thePixelTemp[index_id].entbx[iylow][ihigh].xgx0[i])
01226 //                        +yxratio*((1. - xxratio)*thePixelTemp[index_id].entbx[iyhigh][ilow].xgx0[i] + xxratio*thePixelTemp[index_id].entbx[iyhigh][ihigh].xgx0[i]);
01227                                                         
01228 //            pxgsig[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entbx[iylow][ilow].xgsig[i] + xxratio*thePixelTemp[index_id].entbx[iylow][ihigh].xgsig[i])
01229 //                        +yxratio*((1. - xxratio)*thePixelTemp[index_id].entbx[iyhigh][ilow].xgsig[i] + xxratio*thePixelTemp[index_id].entbx[iyhigh][ihigh].xgsig[i]);
01230                                   
01231               pxavgc2m[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entbx[iylow][ilow].xavgc2m[i] + xxratio*thePixelTemp[index_id].entbx[iylow][ihigh].xavgc2m[i])
01232                           +yxratio*((1. - xxratio)*thePixelTemp[index_id].entbx[iyhigh][ilow].xavgc2m[i] + xxratio*thePixelTemp[index_id].entbx[iyhigh][ihigh].xavgc2m[i]);
01233                   
01234               pxrmsc2m[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entbx[iylow][ilow].xrmsc2m[i] + xxratio*thePixelTemp[index_id].entbx[iylow][ihigh].xrmsc2m[i])
01235                           +yxratio*((1. - xxratio)*thePixelTemp[index_id].entbx[iyhigh][ilow].xrmsc2m[i] + xxratio*thePixelTemp[index_id].entbx[iyhigh][ihigh].xrmsc2m[i]);
01236                   
01237 //            pxgx0c2m[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entbx[iylow][ilow].xgx0c2m[i] + xxratio*thePixelTemp[index_id].entbx[iylow][ihigh].xgx0c2m[i])
01238 //                        +yxratio*((1. - xxratio)*thePixelTemp[index_id].entbx[iyhigh][ilow].xgx0c2m[i] + xxratio*thePixelTemp[index_id].entbx[iyhigh][ihigh].xgx0c2m[i]);
01239                                                         
01240 //            pxgsigc2m[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entbx[iylow][ilow].xgsigc2m[i] + xxratio*thePixelTemp[index_id].entbx[iylow][ihigh].xgsigc2m[i])
01241 //                        +yxratio*((1. - xxratio)*thePixelTemp[index_id].entbx[iyhigh][ilow].xgsigc2m[i] + xxratio*thePixelTemp[index_id].entbx[iyhigh][ihigh].xgsigc2m[i]);
01242 //
01243 //  Try new interpolation scheme
01244 //                                                                                                                                                                                                                                              
01245 //            pchi2xavg[i]=((1. - xxratio)*thePixelTemp[index_id].entbx[imaxx][ilow].chi2xavg[i] + xxratio*thePixelTemp[index_id].entbx[imaxx][ihigh].chi2xavg[i]);
01246 //                if(thePixelTemp[index_id].entbx[imaxx][imidy].chi2xavg[i] != 0.) {pchi2xavg[i]=pchi2xavg[i]/thePixelTemp[index_id].entbx[imaxx][imidy].chi2xavg[i]*chi2xavg[i];}
01247 //                                                      
01248 //            pchi2xmin[i]=((1. - xxratio)*thePixelTemp[index_id].entbx[imaxx][ilow].chi2xmin[i] + xxratio*thePixelTemp[index_id].entbx[imaxx][ihigh].chi2xmin[i]);
01249 //                if(thePixelTemp[index_id].entbx[imaxx][imidy].chi2xmin[i] != 0.) {pchi2xmin[i]=pchi2xmin[i]/thePixelTemp[index_id].entbx[imaxx][imidy].chi2xmin[i]*chi2xmin[i];}
01250                   
01251               pchi2xavg[i]=((1. - xxratio)*thePixelTemp[index_id].entbx[iyhigh][ilow].chi2xavg[i] + xxratio*thePixelTemp[index_id].entbx[iyhigh][ihigh].chi2xavg[i]);
01252                   if(thePixelTemp[index_id].entbx[iyhigh][imidy].chi2xavg[i] != 0.) {pchi2xavg[i]=pchi2xavg[i]/thePixelTemp[index_id].entbx[iyhigh][imidy].chi2xavg[i]*chi2xavg[i];}
01253                                                         
01254               pchi2xmin[i]=((1. - xxratio)*thePixelTemp[index_id].entbx[iyhigh][ilow].chi2xmin[i] + xxratio*thePixelTemp[index_id].entbx[iyhigh][ihigh].chi2xmin[i]);
01255                   if(thePixelTemp[index_id].entbx[iyhigh][imidy].chi2xmin[i] != 0.) {pchi2xmin[i]=pchi2xmin[i]/thePixelTemp[index_id].entbx[iyhigh][imidy].chi2xmin[i]*chi2xmin[i];}
01256                   
01257               for(j=0; j<6 ; ++j) {
01258                  pxflparll[i][j] = thePixelTemp[index_id].entbx[iylow][ilow].xflpar[i][j];
01259                  pxflparlh[i][j] = thePixelTemp[index_id].entbx[iylow][ihigh].xflpar[i][j];
01260                  pxflparhl[i][j] = thePixelTemp[index_id].entbx[iyhigh][ilow].xflpar[i][j];
01261                  pxflparhh[i][j] = thePixelTemp[index_id].entbx[iyhigh][ihigh].xflpar[i][j];
01262                   }
01263            }
01264            
01265 // Do the spares next
01266 
01267 //       for(i=0; i<10; ++i) {
01268 //            pxspare[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entbx[iylow][ilow].xspare[i] + xxratio*thePixelTemp[index_id].entbx[iylow][ihigh].xspare[i])
01269 //                        +yxratio*((1. - xxratio)*thePixelTemp[index_id].entbx[iyhigh][ilow].xspare[i] + xxratio*thePixelTemp[index_id].entbx[iyhigh][ihigh].xspare[i]);
01270 //       }
01271                           
01272 // Interpolate and build the x-template 
01273         
01274            for(i=0; i<9; ++i) {
01275           pxtemp[i][0] = 0.;
01276           pxtemp[i][1] = 0.;
01277               pxtemp[i][BXM2] = 0.;
01278               pxtemp[i][BXM1] = 0.;
01279               for(j=0; j<TXSIZE; ++j) {
01280                 pxtemp[i][j+2]=(1. - xxratio)*thePixelTemp[index_id].entbx[imaxx][ilow].xtemp[i][j] + xxratio*thePixelTemp[index_id].entbx[imaxx][ihigh].xtemp[i][j];
01281               }
01282            }
01283         }
01284   }
01285   return success;
01286 } // interpolate
01287 
01288 
01289 
01290 
01291 
01292 // ************************************************************************************************************ 
01300 // ************************************************************************************************************ 
01301   void SiPixelTemplate::ysigma2(int fypix, int lypix, float sythr, float ysum[25], float ysig2[25])
01302   
01303 {
01304     // Interpolate using quantities already stored in the private variables
01305     
01306     // Local variables 
01307     int i;
01308         float sigi, sigi2, sigi3, sigi4, symax, qscale;
01309         
01310     // Make sure that input is OK
01311     
01312         assert(fypix > 1 && fypix < BYM2);
01313         assert(lypix >= fypix && lypix < BYM2);
01314                      
01315 // Define the maximum signal to use in the parameterization 
01316 
01317        symax = psymax;
01318            if(psymax > psyparmax) {symax = psyparmax;}
01319            
01320 // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis 
01321 
01322            for(i=fypix-2; i<=lypix+2; ++i) {
01323                   if(i < fypix || i > lypix) {
01324            
01325 // Nearest pseudopixels have uncertainties of 50% of threshold, next-nearest have 10% of threshold
01326 
01327                          ysig2[i] = ps50*ps50;
01328                   } else {
01329                          if(ysum[i] < symax) {
01330                                 sigi = ysum[i];
01331                                 qscale = 1.;
01332                          } else {
01333                                 sigi = symax;
01334                                 qscale = ysum[i]/symax;
01335                          }
01336                          sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
01337                          if(i <= BHY) {
01338                                 ysig2[i] = (1.-pyratio)*
01339                                 (pyparl[0][0]+pyparl[0][1]*sigi+pyparl[0][2]*sigi2+pyparl[0][3]*sigi3+pyparl[0][4]*sigi4)
01340                                 + pyratio*
01341                                 (pyparh[0][0]+pyparh[0][1]*sigi+pyparh[0][2]*sigi2+pyparh[0][3]*sigi3+pyparh[0][4]*sigi4);
01342                          } else {
01343                                 ysig2[i] = (1.-pyratio)*
01344                                 (pyparl[1][0]+pyparl[1][1]*sigi+pyparl[1][2]*sigi2+pyparl[1][3]*sigi3+pyparl[1][4]*sigi4)
01345                                 + pyratio*
01346                             (pyparh[1][0]+pyparh[1][1]*sigi+pyparh[1][2]*sigi2+pyparh[1][3]*sigi3+pyparh[1][4]*sigi4);
01347                          }
01348                          ysig2[i] *=qscale;
01349                      if(ysum[i] > sythr) {ysig2[i] = 1.e8;}
01350                          if(ysig2[i] <= 0.) {LOGERROR("SiPixelTemplate") << "neg y-error-squared, id = " << id_current << ", index = " << index_id << 
01351                          ", cot(alpha) = " << cota_current << ", cot(beta) = " << cotb_current << ", fpix = " << fpix_current << " sigi = " << sigi << ENDL;}
01352               }
01353            }
01354         
01355         return;
01356         
01357 } // End ysigma2
01358 
01359 
01360 
01361 
01362 
01363 
01364 
01365 // ************************************************************************************************************ 
01373 // ************************************************************************************************************ 
01374   void SiPixelTemplate::xsigma2(int fxpix, int lxpix, float sxthr, float xsum[11], float xsig2[11])
01375   
01376 {
01377     // Interpolate using quantities already stored in the private variables
01378     
01379     // Local variables 
01380     int i;
01381         float sigi, sigi2, sigi3, sigi4, yint, sxmax, x0, qscale;
01382         
01383     // Make sure that input is OK
01384     
01385         assert(fxpix > 1 && fxpix < BXM2);
01386         assert(lxpix >= fxpix && lxpix < BXM2);
01387                      
01388 // Define the maximum signal to use in the parameterization 
01389 
01390        sxmax = psxmax;
01391            if(psxmax > psxparmax) {sxmax = psxparmax;}
01392            
01393 // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis 
01394 
01395            for(i=fxpix-2; i<=lxpix+2; ++i) {
01396                   if(i < fxpix || i > lxpix) {
01397            
01398 // Nearest pseudopixels have uncertainties of 50% of threshold, next-nearest have 10% of threshold
01399 
01400                          xsig2[i] = ps50*ps50;
01401                   } else {
01402                          if(xsum[i] < sxmax) {
01403                                 sigi = xsum[i];
01404                                 qscale = 1.;
01405                          } else {
01406                                 sigi = sxmax;
01407                                 qscale = xsum[i]/sxmax;
01408                          }
01409                          sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
01410                          
01411 // First, do the cotbeta interpolation                   
01412                          
01413                          if(i <= BHX) {
01414                                 yint = (1.-pyratio)*
01415                                 (pxparly0[0][0]+pxparly0[0][1]*sigi+pxparly0[0][2]*sigi2+pxparly0[0][3]*sigi3+pxparly0[0][4]*sigi4)
01416                                 + pyratio*
01417                                 (pxparhy0[0][0]+pxparhy0[0][1]*sigi+pxparhy0[0][2]*sigi2+pxparhy0[0][3]*sigi3+pxparhy0[0][4]*sigi4);
01418                          } else {
01419                                 yint = (1.-pyratio)*
01420                                 (pxparly0[1][0]+pxparly0[1][1]*sigi+pxparly0[1][2]*sigi2+pxparly0[1][3]*sigi3+pxparly0[1][4]*sigi4)
01421                                 + pyratio*
01422                             (pxparhy0[1][0]+pxparhy0[1][1]*sigi+pxparhy0[1][2]*sigi2+pxparhy0[1][3]*sigi3+pxparhy0[1][4]*sigi4);
01423                          }
01424                          
01425 // Next, do the cotalpha interpolation                   
01426                          
01427                          if(i <= BHX) {
01428                                 xsig2[i] = (1.-pxxratio)*
01429                                 (pxparl[0][0]+pxparl[0][1]*sigi+pxparl[0][2]*sigi2+pxparl[0][3]*sigi3+pxparl[0][4]*sigi4)
01430                                 + pxxratio*
01431                                 (pxparh[0][0]+pxparh[0][1]*sigi+pxparh[0][2]*sigi2+pxparh[0][3]*sigi3+pxparh[0][4]*sigi4);
01432                          } else {
01433                                 xsig2[i] = (1.-pxxratio)*
01434                                 (pxparl[1][0]+pxparl[1][1]*sigi+pxparl[1][2]*sigi2+pxparl[1][3]*sigi3+pxparl[1][4]*sigi4)
01435                                 + pxxratio*
01436                             (pxparh[1][0]+pxparh[1][1]*sigi+pxparh[1][2]*sigi2+pxparh[1][3]*sigi3+pxparh[1][4]*sigi4);
01437                          }
01438                          
01439 // Finally, get the mid-point value of the cotalpha function                     
01440                          
01441                          if(i <= BHX) {
01442                                 x0 = pxpar0[0][0]+pxpar0[0][1]*sigi+pxpar0[0][2]*sigi2+pxpar0[0][3]*sigi3+pxpar0[0][4]*sigi4;
01443                          } else {
01444                                 x0 = pxpar0[1][0]+pxpar0[1][1]*sigi+pxpar0[1][2]*sigi2+pxpar0[1][3]*sigi3+pxpar0[1][4]*sigi4;
01445                          }
01446                          
01447 // Finally, rescale the yint value for cotalpha variation                        
01448                          
01449                          if(x0 != 0.) {xsig2[i] = xsig2[i]/x0 * yint;}
01450                          xsig2[i] *=qscale;
01451                      if(xsum[i] > sxthr) {xsig2[i] = 1.e8;}
01452                          if(xsig2[i] <= 0.) {LOGERROR("SiPixelTemplate") << "neg x-error-squared, id = " << id_current << ", index = " << index_id << 
01453                          ", cot(alpha) = " << cota_current << ", cot(beta) = " << cotb_current << ", fpix = " << fpix_current << " sigi = " << sigi << ENDL;}
01454               }
01455            }
01456         
01457         return;
01458         
01459 } // End xsigma2
01460 
01461 
01462 
01463 
01464 
01465 
01466 // ************************************************************************************************************ 
01470 // ************************************************************************************************************ 
01471   float SiPixelTemplate::yflcorr(int binq, float qfly)
01472   
01473 {
01474     // Interpolate using quantities already stored in the private variables
01475     
01476     // Local variables 
01477     int i;
01478         float qfl, qfl2, qfl3, qfl4, qfl5, dy;
01479         
01480     // Make sure that input is OK
01481     
01482         assert(binq >= 0 && binq < 4);
01483         assert(fabs((double)qfly) <= 1.);
01484                      
01485 // Define the maximum signal to allow before de-weighting a pixel 
01486 
01487        qfl = qfly;
01488 
01489        if(qfl < -0.9) {qfl = -0.9;}
01490            if(qfl > 0.9) {qfl = 0.9;}
01491            
01492 // Interpolate between the two polynomials
01493 
01494            qfl2 = qfl*qfl; qfl3 = qfl2*qfl; qfl4 = qfl3*qfl; qfl5 = qfl4*qfl;
01495            dy = (1.-pyratio)*(pyflparl[binq][0]+pyflparl[binq][1]*qfl+pyflparl[binq][2]*qfl2+pyflparl[binq][3]*qfl3+pyflparl[binq][4]*qfl4+pyflparl[binq][5]*qfl5)
01496                   + pyratio*(pyflparh[binq][0]+pyflparh[binq][1]*qfl+pyflparh[binq][2]*qfl2+pyflparh[binq][3]*qfl3+pyflparh[binq][4]*qfl4+pyflparh[binq][5]*qfl5);
01497         
01498         return dy;
01499         
01500 } // End yflcorr
01501 
01502 
01503 
01504 
01505 
01506 
01507 // ************************************************************************************************************ 
01511 // ************************************************************************************************************ 
01512   float SiPixelTemplate::xflcorr(int binq, float qflx)
01513   
01514 {
01515     // Interpolate using quantities already stored in the private variables
01516     
01517     // Local variables 
01518     int i;
01519         float qfl, qfl2, qfl3, qfl4, qfl5, dx;
01520         
01521     // Make sure that input is OK
01522     
01523         assert(binq >= 0 && binq < 4);
01524         assert(fabs((double)qflx) <= 1.);
01525                      
01526 // Define the maximum signal to allow before de-weighting a pixel 
01527 
01528        qfl = qflx;
01529 
01530        if(qfl < -0.9) {qfl = -0.9;}
01531            if(qfl > 0.9) {qfl = 0.9;}
01532            
01533 // Interpolate between the two polynomials
01534 
01535            qfl2 = qfl*qfl; qfl3 = qfl2*qfl; qfl4 = qfl3*qfl; qfl5 = qfl4*qfl;
01536            dx = (1. - pyxratio)*((1.-pxxratio)*(pxflparll[binq][0]+pxflparll[binq][1]*qfl+pxflparll[binq][2]*qfl2+pxflparll[binq][3]*qfl3+pxflparll[binq][4]*qfl4+pxflparll[binq][5]*qfl5)
01537                   + pxxratio*(pxflparlh[binq][0]+pxflparlh[binq][1]*qfl+pxflparlh[binq][2]*qfl2+pxflparlh[binq][3]*qfl3+pxflparlh[binq][4]*qfl4+pxflparlh[binq][5]*qfl5))
01538               + pyxratio*((1.-pxxratio)*(pxflparhl[binq][0]+pxflparhl[binq][1]*qfl+pxflparhl[binq][2]*qfl2+pxflparhl[binq][3]*qfl3+pxflparhl[binq][4]*qfl4+pxflparhl[binq][5]*qfl5)
01539                   + pxxratio*(pxflparhh[binq][0]+pxflparhh[binq][1]*qfl+pxflparhh[binq][2]*qfl2+pxflparhh[binq][3]*qfl3+pxflparhh[binq][4]*qfl4+pxflparhh[binq][5]*qfl5));
01540         
01541         return dx;
01542         
01543 } // End xflcorr
01544 
01545 
01546 
01547 // ************************************************************************************************************ 
01552 // ************************************************************************************************************ 
01553   void SiPixelTemplate::ytemp(int fybin, int lybin, float ytemplate[41][BYSIZE])
01554   
01555 {
01556     // Retrieve already interpolated quantities
01557     
01558     // Local variables 
01559     int i, j;
01560 
01561    // Verify that input parameters are in valid range
01562 
01563         assert(fybin >= 0 && fybin < 41);
01564         assert(lybin >= 0 && lybin < 41);
01565 
01566 // Build the y-template, the central 25 bins are here in all cases
01567         
01568         for(i=0; i<9; ++i) {
01569            for(j=0; j<BYSIZE; ++j) {
01570                         ytemplate[i+16][j]=pytemp[i][j];
01571            }
01572         }
01573         for(i=0; i<8; ++i) {
01574            ytemplate[i+8][BYM1] = 0.;
01575            for(j=0; j<BYM1; ++j) {
01576               ytemplate[i+8][j]=pytemp[i][j+1];
01577            }
01578         }
01579         for(i=1; i<9; ++i) {
01580            ytemplate[i+24][0] = 0.;
01581            for(j=0; j<BYM1; ++j) {
01582               ytemplate[i+24][j+1]=pytemp[i][j];
01583            }
01584         }
01585         
01586 //  Add more bins if needed
01587 
01588         if(fybin < 8) {
01589            for(i=0; i<8; ++i) {
01590               ytemplate[i][BYM2] = 0.;
01591               ytemplate[i][BYM1] = 0.;
01592               for(j=0; j<BYM2; ++j) {
01593                 ytemplate[i][j]=pytemp[i][j+2];
01594               }
01595            }
01596         }
01597         if(lybin > 32) {
01598            for(i=1; i<9; ++i) {
01599           ytemplate[i+32][0] = 0.;
01600               ytemplate[i+32][1] = 0.;
01601               for(j=0; j<BYM2; ++j) {
01602                  ytemplate[i+32][j+2]=pytemp[i][j];
01603               }
01604            }
01605         }
01606         
01607         return;
01608         
01609 } // End ytemp
01610 
01611 
01612 
01613 // ************************************************************************************************************ 
01618 // ************************************************************************************************************ 
01619   void SiPixelTemplate::xtemp(int fxbin, int lxbin, float xtemplate[41][BXSIZE])
01620   
01621 {
01622     // Retrieve already interpolated quantities
01623     
01624     // Local variables 
01625     int i, j;
01626 
01627    // Verify that input parameters are in valid range
01628 
01629         assert(fxbin >= 0 && fxbin < 41);
01630         assert(lxbin >= 0 && lxbin < 41);
01631 
01632 // Build the x-template, the central 25 bins are here in all cases
01633         
01634         for(i=0; i<9; ++i) {
01635            for(j=0; j<BXSIZE; ++j) {
01636               xtemplate[i+16][j]=pxtemp[i][j];
01637            }
01638         }
01639         for(i=0; i<8; ++i) {
01640            xtemplate[i+8][BXM1] = 0.;
01641            for(j=0; j<BXM1; ++j) {
01642               xtemplate[i+8][j]=pxtemp[i][j+1];
01643            }
01644         }
01645         for(i=1; i<9; ++i) {
01646            xtemplate[i+24][0] = 0.;
01647            for(j=0; j<BXM1; ++j) {
01648               xtemplate[i+24][j+1]=pxtemp[i][j];
01649            }
01650         }
01651         
01652 //  Add more bins if needed     
01653         
01654         if(fxbin < 8) {
01655            for(i=0; i<8; ++i) {
01656           xtemplate[i][BXM2] = 0.;
01657               xtemplate[i][BXM1] = 0.;
01658               for(j=0; j<BXM2; ++j) {
01659                 xtemplate[i][j]=pxtemp[i][j+2];
01660               }
01661            }
01662         }
01663         if(lxbin > 32) {
01664            for(i=1; i<9; ++i) {
01665           xtemplate[i+32][0] = 0.;
01666               xtemplate[i+32][1] = 0.;
01667               for(j=0; j<BXM2; ++j) {
01668                 xtemplate[i+32][j+2]=pxtemp[i][j];
01669               }
01670            }
01671         }
01672         
01673         return;
01674         
01675 } // End xtemp
01676 
01677 
01678 // ************************************************************************************************************ 
01682 // ************************************************************************************************************ 
01683   void SiPixelTemplate::ytemp3d(int nypix, array_3d& ytemplate)
01684   
01685 {
01686     typedef boost::multi_array<float, 2> array_2d;
01687         
01688     // Retrieve already interpolated quantities
01689     
01690     // Local variables 
01691     int i, j, k;
01692         int ioff0, ioffp, ioffm;
01693 
01694    // Verify that input parameters are in valid range
01695 
01696         assert(nypix > 0 && nypix < BYM3);
01697         
01698 // Calculate the size of the shift in pixels needed to span the entire cluster
01699 
01700     float diff = fabsf(nypix - pclsleny)/2. + 1.;
01701         int nshift = (int) diff;
01702         if((diff - nshift) > 0.5) {++nshift;}
01703 
01704 // Calculate the number of bins needed to specify each hit range
01705 
01706     int nbins = 9 + 16*nshift;
01707         
01708 // Create a 2-d working template with the correct size
01709 
01710         array_2d temp2d(boost::extents[nbins][BYSIZE]);
01711         
01712 //  The 9 central bins are copied from the interpolated private store
01713 
01714         ioff0 = 8*nshift;
01715         
01716         for(i=0; i<9; ++i) {
01717            for(j=0; j<BYSIZE; ++j) {
01718           temp2d[i+ioff0][j]=pytemp[i][j];
01719            }
01720         }
01721         
01722 // Add the +- shifted templates 
01723 
01724         for(k=1; k<=nshift; ++k) {
01725           ioffm=ioff0-k*8;
01726           for(i=0; i<8; ++i) {
01727              for(j=0; j<k; ++j) {
01728                 temp2d[i+ioffm][BYM1-j] = 0.;
01729                  }
01730              for(j=0; j<BYSIZE-k; ++j) {
01731                 temp2d[i+ioffm][j]=pytemp[i][j+k];
01732              }
01733            }
01734            ioffp=ioff0+k*8;
01735            for(i=1; i<9; ++i) {
01736               for(j=0; j<k; ++j) {
01737                  temp2d[i+ioffp][j] = 0.;
01738               }
01739               for(j=0; j<BYSIZE-k; ++j) {
01740                  temp2d[i+ioffp][j+k]=pytemp[i][j];
01741               }
01742            }
01743         }
01744                 
01745 // Resize the 3d template container
01746 
01747     ytemplate.resize(boost::extents[nbins][nbins][BYSIZE]);
01748         
01749 // Sum two 2-d templates to make the 3-d template
01750         
01751    for(i=0; i<nbins; ++i) {
01752       for(j=0; j<=i; ++j) {
01753          for(k=0; k<BYSIZE; ++k) {
01754             ytemplate[i][j][k]=temp2d[i][k]+temp2d[j][k];       
01755                  }
01756           }
01757    }
01758         
01759         return;
01760         
01761 } // End ytemp3d
01762 
01763 
01764 // ************************************************************************************************************ 
01768 // ************************************************************************************************************ 
01769   void SiPixelTemplate::xtemp3d(int nxpix, array_3d& xtemplate)
01770   
01771 {
01772     typedef boost::multi_array<float, 2> array_2d;
01773         
01774     // Retrieve already interpolated quantities
01775     
01776     // Local variables 
01777     int i, j, k;
01778         int ioff0, ioffp, ioffm;
01779 
01780    // Verify that input parameters are in valid range
01781 
01782         assert(nxpix > 0 && nxpix < BXM3);
01783         
01784 // Calculate the size of the shift in pixels needed to span the entire cluster
01785 
01786     float diff = fabsf(nxpix - pclslenx)/2. + 1.;
01787         int nshift = (int) diff;
01788         if((diff - nshift) > 0.5) {++nshift;}
01789 
01790 // Calculate the number of bins needed to specify each hit range
01791 
01792     int nbins = 9 + 16*nshift;
01793         
01794 // Create a 2-d working template with the correct size
01795 
01796         array_2d temp2d(boost::extents[nbins][BXSIZE]);
01797         
01798 //  The 9 central bins are copied from the interpolated private store
01799 
01800         ioff0 = 8*nshift;
01801         
01802         for(i=0; i<9; ++i) {
01803            for(j=0; j<BXSIZE; ++j) {
01804           temp2d[i+ioff0][j]=pxtemp[i][j];
01805            }
01806         }
01807         
01808 // Add the +- shifted templates 
01809 
01810         for(k=1; k<=nshift; ++k) {
01811           ioffm=ioff0-k*8;
01812           for(i=0; i<8; ++i) {
01813              for(j=0; j<k; ++j) {
01814                 temp2d[i+ioffm][BXM1-j] = 0.;
01815                  }
01816              for(j=0; j<BXSIZE-k; ++j) {
01817                 temp2d[i+ioffm][j]=pxtemp[i][j+k];
01818              }
01819            }
01820            ioffp=ioff0+k*8;
01821            for(i=1; i<9; ++i) {
01822               for(j=0; j<k; ++j) {
01823                  temp2d[i+ioffp][j] = 0.;
01824               }
01825               for(j=0; j<BXSIZE-k; ++j) {
01826                  temp2d[i+ioffp][j+k]=pxtemp[i][j];
01827               }
01828            }
01829         }
01830                                 
01831 // Resize the 3d template container
01832 
01833     xtemplate.resize(boost::extents[nbins][nbins][BXSIZE]);
01834         
01835 // Sum two 2-d templates to make the 3-d template
01836         
01837    for(i=0; i<nbins; ++i) {
01838       for(j=0; j<=i; ++j) {
01839          for(k=0; k<BXSIZE; ++k) {
01840             xtemplate[i][j][k]=temp2d[i][k]+temp2d[j][k];       
01841                  }
01842           }
01843    }
01844         
01845         return;
01846         
01847 } // End xtemp3d
01848 
01849 
01850 
01851 // ************************************************************************************************************ 
01859 // ************************************************************************************************************ 
01860 int SiPixelTemplate::qbin(int id, bool fpix, float cotbeta, float qclus)
01861 {
01862     // Interpolate for a new set of track angles 
01863     
01864     // Local variables 
01865     int i, j, binq;
01866         int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, imidy, imaxx;
01867         float yratio, yxratio, xxratio, sxmax;
01868         float acotb, qscale, qavg, qmin, fq, qtotal;
01869         
01870         if(id != id_current) {
01871 
01872 // Find the index corresponding to id
01873 
01874        index_id = -1;
01875        for(i=0; i<thePixelTemp.size(); ++i) {
01876         
01877               if(id == thePixelTemp[i].head.ID) {
01878            
01879                  index_id = i;
01880                      id_current = id;
01881                      break;
01882           }
01883             }
01884      }
01885          
01886          assert(index_id >= 0 && index_id < thePixelTemp.size());
01887          
01888 //              
01889 
01890 // Interpolate the absolute value of cot(beta)     
01891     
01892     acotb = fabs((double)cotbeta);
01893 
01894 // Copy the charge scaling factor to the private variable     
01895     
01896     qscale = thePixelTemp[index_id].head.qscale;
01897         
01898 // Decide which template (FPix or BPix) to use 
01899 
01900     if(fpix) {
01901     
01902 // Begin FPix section, make the index counters easier to use     
01903     
01904        Ny = thePixelTemp[index_id].head.NFy;
01905        Nyx = thePixelTemp[index_id].head.NFyx;
01906        Nxx = thePixelTemp[index_id].head.NFxx;
01907            imaxx = Nyx - 1;
01908            imidy = Nxx/2;
01909         
01910 // next, loop over all y-angle entries   
01911 
01912            ilow = 0;
01913            yratio = 0.;
01914 
01915            if(acotb >= thePixelTemp[index_id].entfy[Ny-1].cotbeta) {
01916         
01917                ilow = Ny-2;
01918                    yratio = 1.;
01919                 
01920            } else if(acotb >= thePixelTemp[index_id].entfy[0].cotbeta) {
01921 
01922           for (i=0; i<Ny-1; ++i) { 
01923     
01924              if( thePixelTemp[index_id].entfy[i].cotbeta <= acotb && acotb < thePixelTemp[index_id].entfy[i+1].cotbeta) {
01925                   
01926                     ilow = i;
01927                         yratio = (acotb - thePixelTemp[index_id].entfy[i].cotbeta)/(thePixelTemp[index_id].entfy[i+1].cotbeta - thePixelTemp[index_id].entfy[i].cotbeta);
01928                         break;                   
01929                      }
01930               }
01931            }
01932         
01933            ihigh=ilow + 1;
01934                           
01935 // Interpolate/store all y-related quantities (flip displacements when cotbeta < 0)
01936 
01937            qavg = (1. - yratio)*thePixelTemp[index_id].entfy[ilow].qavg + yratio*thePixelTemp[index_id].entfy[ihigh].qavg;
01938            qmin = (1. - yratio)*thePixelTemp[index_id].entfy[ilow].qmin + yratio*thePixelTemp[index_id].entfy[ihigh].qmin;
01939 
01940         } else {
01941         
01942     
01943 // Begin BPix section, make the index counters easier to use     
01944     
01945        Ny = thePixelTemp[index_id].head.NBy;
01946        Nyx = thePixelTemp[index_id].head.NByx;
01947        Nxx = thePixelTemp[index_id].head.NBxx;
01948            imaxx = Nyx - 1;
01949            imidy = Nxx/2;
01950         
01951 // next, loop over all y-angle entries   
01952 
01953            ilow = 0;
01954            yratio = 0.;
01955 
01956            if(acotb >= thePixelTemp[index_id].entby[Ny-1].cotbeta) {
01957         
01958                ilow = Ny-2;
01959                    yratio = 1.;
01960                 
01961            } else if(acotb >= thePixelTemp[index_id].entby[0].cotbeta) {
01962 
01963           for (i=0; i<Ny-1; ++i) { 
01964     
01965              if( thePixelTemp[index_id].entby[i].cotbeta <= acotb && abs_cotb < thePixelTemp[index_id].entby[i+1].cotbeta) {
01966                   
01967                     ilow = i;
01968                         yratio = (acotb - thePixelTemp[index_id].entby[i].cotbeta)/(thePixelTemp[index_id].entby[i+1].cotbeta - thePixelTemp[index_id].entby[i].cotbeta);
01969                         break;                   
01970                      }
01971               }
01972            }
01973         
01974            ihigh=ilow + 1;
01975                           
01976 // Interpolate/store all y-related quantities (flip displacements when cotbeta < 0)
01977 
01978            qavg = (1. - yratio)*thePixelTemp[index_id].entby[ilow].qavg + yratio*thePixelTemp[index_id].entby[ihigh].qavg;
01979            qmin = (1. - yratio)*thePixelTemp[index_id].entby[ilow].qmin + yratio*thePixelTemp[index_id].entby[ihigh].qmin;
01980         }
01981         
01982         assert(qavg > 0. && qmin > 0.);
01983         
01984 //  Scale the input charge to account for differences between pixelav and CMSSW simulation or data      
01985         
01986         qtotal = qscale*qclus;
01987         
01988 // uncertainty and final corrections depend upon total charge bin          
01989            
01990         fq = qtotal/qavg;
01991         if(fq > 1.5) {
01992            binq=0;
01993         } else {
01994            if(fq > 1.0) {
01995               binq=1;
01996            } else {
01997                   if(fq > 0.85) {
01998                          binq=2;
01999                   } else {
02000                          binq=3;
02001                   }
02002            }
02003         }
02004         
02005 // If the charge is too small (then flag it)
02006         
02007         if(qtotal < 0.95*qmin) {binq = 4;}
02008                 
02009     return binq;
02010   
02011 } // qbin
02012 
02013 
02014 
02015 
02016