CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_2_9_HLT1_bphpatch4/src/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplate.cc

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