Reference Manual

Go to the documentation of this file.
00001 //
00002 //  SiStripTemplate.cc  Version 1.00 (based on SiPixelTemplate v8.20)
00003 //
00004 //  V1.05 - add VI optimizations from pixel template object
00005 //  V1.06 - increase angular acceptance (and structure size)
00006 //  V2.00 - add barycenter interpolation and getters, fix calculation for charge deposition to accommodate cota-offsets in the central cotb entries.
00007 //  V2.01 - fix problem with number of spare entries
00008 //
00009 
00010 //  Created by Morris Swartz on 2/2/11.
00011 //
00012 //
00013 
00014 //#include <stdlib.h> 
00015 //#include <stdio.h>
00016 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00017 //#include <cmath.h>
00018 #else
00019 #include <math.h>
00020 #endif
00021 #include <algorithm>
00022 #include <vector>
00023 //#include "boost/multi_array.hpp"
00024 #include <iostream>
00025 #include <iomanip>
00026 #include <sstream>
00027 #include <fstream>
00028 #include <list>
00029 
00030 
00031 
00032 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00033 #include "RecoLocalTracker/SiStripRecHitConverter/interface/SiStripTemplate.h"
00034 #include "FWCore/ParameterSet/interface/FileInPath.h"
00035 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00036 #define LOGERROR(x) LogError(x)
00037 #define LOGINFO(x) LogInfo(x)
00038 #define ENDL " "
00039 #include "FWCore/Utilities/interface/Exception.h"
00040 using namespace edm;
00041 #else
00042 #include "SiStripTemplate.h"
00043 #define LOGERROR(x) std::cout << x << ": "
00044 #define LOGINFO(x) std::cout << x << ": "
00045 #define ENDL std::endl
00046 #endif
00047 
00048 //**************************************************************** 
00053 //**************************************************************** 
00054 bool SiStripTemplate::pushfile(int filenum)
00055 {
00056     // Add template stored in external file numbered filenum to theTemplateStore
00057     
00058     // Local variables 
00059     int i, j, k, l;
00060         float qavg_avg;
00061         const char *tempfile;
00062         //      char title[80]; remove this
00063     char c;
00064         const int code_version={18};
00065         
00066         
00067         
00068         //  Create a filename for this run 
00069         
00070         std::ostringstream tout;
00071         
00072         //  Create different path in CMSSW than standalone
00073         
00074 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00075         tout << "CalibTracker/SiPixelESProducers/data/stemplate_summary_p" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
00076         std::string tempf = tout.str();
00077         edm::FileInPath file( tempf.c_str() );
00078         tempfile = (file.fullPath()).c_str();
00079 #else
00080         tout << "stemplate_summary_p" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
00081         std::string tempf = tout.str();
00082         tempfile = tempf.c_str();
00083 #endif
00084         
00085         //  open the template file 
00086         
00087         std::ifstream in_file(tempfile, std::ios::in);
00088         
00089         if(in_file.is_open()) {
00090                 
00091                 // Create a local template storage entry
00092                 
00093                 SiStripTemplateStore theCurrentTemp;
00094                 
00095                 // Read-in a header string first and print it    
00096                 
00097                 for (i=0; (c=in_file.get()) != '\n'; ++i) {
00098                         if(i < 79) {theCurrentTemp.head.title[i] = c;}
00099                 }
00100                 if(i > 78) {i=78;}
00101                 theCurrentTemp.head.title[i+1] ='\0';
00102                 LOGINFO("SiStripTemplate") << "Loading Strip Template File - " << theCurrentTemp.head.title << ENDL;
00103                 
00104                 // next, the header information     
00105                 
00106                 in_file >> theCurrentTemp.head.ID  >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx
00107                 >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
00108                 >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >> theCurrentTemp.head.zsize;
00109                 
00110                 if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file, no template load" << ENDL; return false;}
00111                 
00112                 LOGINFO("SiStripTemplate") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version " << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield 
00113                 << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
00114                 << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
00115                 << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
00116                 << ", 1/2 threshold " << theCurrentTemp.head.s50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", x Lorentz width " << theCurrentTemp.head.lorxwidth    
00117                 << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
00118                 
00119                 if(theCurrentTemp.head.templ_version < code_version) {LOGERROR("SiStripTemplate") << "code expects version " << code_version << ", no template load" << ENDL; return false;}
00120                 
00121 #ifdef SI_STRIP_TEMPLATE_USE_BOOST 
00122                 
00123 // next, layout the 1-d/2-d structures needed to store template
00124                                 
00125                 theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
00126 
00127                 theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
00128                 
00129 #endif
00130                 
00131 // next, loop over all y-angle entries   
00132                 
00133                 for (i=0; i < theCurrentTemp.head.NTy; ++i) {     
00134                         
00135                         in_file >> theCurrentTemp.enty[i].runnum >> theCurrentTemp.enty[i].costrk[0] 
00136                         >> theCurrentTemp.enty[i].costrk[1] >> theCurrentTemp.enty[i].costrk[2]; 
00137                         
00138                         if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 1, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00139                         
00140                         // Calculate the alpha, beta, and cot(beta) for this entry 
00141                         
00142                         theCurrentTemp.enty[i].alpha = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[0]));
00143                         
00144                         theCurrentTemp.enty[i].cotalpha = theCurrentTemp.enty[i].costrk[0]/theCurrentTemp.enty[i].costrk[2];
00145                         
00146                         theCurrentTemp.enty[i].beta = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[1]));
00147                         
00148                         theCurrentTemp.enty[i].cotbeta = theCurrentTemp.enty[i].costrk[1]/theCurrentTemp.enty[i].costrk[2];
00149                         
00150                         in_file >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].sxmax >> theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].clslenx;
00151                         
00152                         if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 2, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00153                         for (j=0; j<2; ++j) {
00154                                 
00155                                 in_file >> theCurrentTemp.enty[i].xpar[j][0] >> theCurrentTemp.enty[i].xpar[j][1] 
00156                                 >> theCurrentTemp.enty[i].xpar[j][2] >> theCurrentTemp.enty[i].xpar[j][3] >> theCurrentTemp.enty[i].xpar[j][4];
00157                                 
00158                                 if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 6, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00159                                 
00160                         }
00161                         
00162                         qavg_avg = 0.f;
00163                         for (j=0; j<9; ++j) {
00164                                 
00165                                 for (k=0; k<TSXSIZE; ++k) {in_file >> theCurrentTemp.enty[i].xtemp[j][k]; qavg_avg += theCurrentTemp.enty[i].xtemp[j][k];} 
00166                                 
00167                                 if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 7, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00168                         }
00169                         theCurrentTemp.enty[i].qavg_avg = qavg_avg/9.;
00170                         
00171                         for (j=0; j<4; ++j) {
00172                                 
00173                                 in_file >> theCurrentTemp.enty[i].xavg[j] >> theCurrentTemp.enty[i].xrms[j] >> theCurrentTemp.enty[i].xgx0[j] >> theCurrentTemp.enty[i].xgsig[j];
00174                                 
00175                                 if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 10, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00176                         }
00177                         
00178                         for (j=0; j<4; ++j) {
00179                                 
00180                                 in_file >> theCurrentTemp.enty[i].xflpar[j][0] >> theCurrentTemp.enty[i].xflpar[j][1] >> theCurrentTemp.enty[i].xflpar[j][2] 
00181                                 >> theCurrentTemp.enty[i].xflpar[j][3] >> theCurrentTemp.enty[i].xflpar[j][4] >> theCurrentTemp.enty[i].xflpar[j][5];
00182                                 
00183                                 if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 11, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00184                         }
00185                         
00186                         for (j=0; j<4; ++j) {
00187                                 
00188                                 in_file >> theCurrentTemp.enty[i].chi2xavg[j] >> theCurrentTemp.enty[i].chi2xmin[j] >> theCurrentTemp.enty[i].chi2xavgc2m[j] >> theCurrentTemp.enty[i].chi2xminc2m[j];
00189                                 
00190                                 if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 12, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00191                         }
00192                         
00193                                 for (j=0; j<4; ++j) {
00194                                 
00195                                 in_file >> theCurrentTemp.enty[i].xavgc2m[j] >> theCurrentTemp.enty[i].xrmsc2m[j] >> theCurrentTemp.enty[i].xgx0c2m[j] >> theCurrentTemp.enty[i].xgsigc2m[j];
00196                                 
00197                                 if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 14, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00198                         } 
00199                         
00200                                 for (j=0; j<4; ++j) {
00201                                 
00202                                 in_file >> theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j] >> theCurrentTemp.enty[i].xgx0gen[j] >> theCurrentTemp.enty[i].xgsiggen[j];
00203                                 
00204                                 if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 14b, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00205                         } 
00206                         
00207                         for (j=0; j<4; ++j) {
00208                                 
00209                                 in_file >> theCurrentTemp.enty[i].xavgbcn[j] >> theCurrentTemp.enty[i].xrmsbcn[j] >> theCurrentTemp.enty[i].xgx0bcn[j] >> theCurrentTemp.enty[i].xgsigbcn[j];
00210                                 
00211                                 if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 14c, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00212                         } 
00213                         
00214                         in_file >> theCurrentTemp.enty[i].chi2xavgone >> theCurrentTemp.enty[i].chi2xminone >> theCurrentTemp.enty[i].qmin2
00215                         >> theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav 
00216                         >> theCurrentTemp.enty[i].mpvvav2 >> theCurrentTemp.enty[i].sigmavav2 >> theCurrentTemp.enty[i].kappavav2 >> theCurrentTemp.enty[i].spare[0];
00217                         
00218                         if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 15, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00219                         
00220                         in_file >> theCurrentTemp.enty[i].qbfrac[0] >> theCurrentTemp.enty[i].qbfrac[1] >> theCurrentTemp.enty[i].qbfrac[2] >> theCurrentTemp.enty[i].fracxone
00221                         >> theCurrentTemp.enty[i].spare[1] >> theCurrentTemp.enty[i].spare[2] >> theCurrentTemp.enty[i].spare[3] >> theCurrentTemp.enty[i].spare[4]
00222                         >> theCurrentTemp.enty[i].spare[5] >> theCurrentTemp.enty[i].spare[6];
00223                         
00224                         if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 16, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00225                         
00226                 }
00227                 
00228                 // next, loop over all barrel x-angle entries   
00229                 
00230                 for (k=0; k < theCurrentTemp.head.NTyx; ++k) { 
00231                         
00232                         for (i=0; i < theCurrentTemp.head.NTxx; ++i) { 
00233                                 
00234                                 in_file >> theCurrentTemp.entx[k][i].runnum >> theCurrentTemp.entx[k][i].costrk[0] 
00235                                 >> theCurrentTemp.entx[k][i].costrk[1] >> theCurrentTemp.entx[k][i].costrk[2]; 
00236                                 
00237                                 if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 17, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00238                                 
00239                                 // Calculate the alpha, beta, and cot(beta) for this entry 
00240                                 
00241                                 theCurrentTemp.entx[k][i].alpha = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[0]));
00242                                 
00243                                 theCurrentTemp.entx[k][i].cotalpha = theCurrentTemp.entx[k][i].costrk[0]/theCurrentTemp.entx[k][i].costrk[2];
00244                                 
00245                                 theCurrentTemp.entx[k][i].beta = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[1]));
00246                                 
00247                                 theCurrentTemp.entx[k][i].cotbeta = theCurrentTemp.entx[k][i].costrk[1]/theCurrentTemp.entx[k][i].costrk[2];
00248                                 
00249                                 in_file >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].sxmax >> theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone >> theCurrentTemp.entx[k][i].qmin >> theCurrentTemp.entx[k][i].clslenx;
00250                                 
00251                                 if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 18, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00252                                 
00253                                 for (j=0; j<2; ++j) {
00254                                         
00255                                         in_file >> theCurrentTemp.entx[k][i].xpar[j][0] >> theCurrentTemp.entx[k][i].xpar[j][1] 
00256                                         >> theCurrentTemp.entx[k][i].xpar[j][2] >> theCurrentTemp.entx[k][i].xpar[j][3] >> theCurrentTemp.entx[k][i].xpar[j][4];
00257                                         
00258                                         if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 19, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00259                                         
00260                                 }
00261                                 
00262                                 qavg_avg = 0.f;
00263                                 for (j=0; j<9; ++j) {
00264                                         
00265                                         for (l=0; l<TSXSIZE; ++l) {in_file >> theCurrentTemp.entx[k][i].xtemp[j][l]; qavg_avg += theCurrentTemp.entx[k][i].xtemp[j][l];} 
00266                                         
00267                                         if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 20, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00268                                 }
00269                                 theCurrentTemp.entx[k][i].qavg_avg = qavg_avg/9.;
00270                                 
00271                                 for (j=0; j<4; ++j) {
00272                                         
00273                                         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];
00274                                         
00275                                         if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 21, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00276                                 }
00277                                 
00278                                 for (j=0; j<4; ++j) {
00279                                         
00280                                         in_file >> theCurrentTemp.entx[k][i].xflpar[j][0] >> theCurrentTemp.entx[k][i].xflpar[j][1] >> theCurrentTemp.entx[k][i].xflpar[j][2] 
00281                                         >> theCurrentTemp.entx[k][i].xflpar[j][3] >> theCurrentTemp.entx[k][i].xflpar[j][4] >> theCurrentTemp.entx[k][i].xflpar[j][5];
00282                                         
00283                                         if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 22, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00284                                 }
00285                                 
00286                                 for (j=0; j<4; ++j) {
00287                                         
00288                                         in_file >> theCurrentTemp.entx[k][i].chi2xavg[j] >> theCurrentTemp.entx[k][i].chi2xmin[j] >> theCurrentTemp.entx[k][i].chi2xavgc2m[j] >> theCurrentTemp.entx[k][i].chi2xminc2m[j];
00289                                         
00290                                         if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 23, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00291                                 }
00292                                 
00293                                 for (j=0; j<4; ++j) {
00294                                         
00295                                         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];
00296                                         
00297                                         if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 24, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00298                                 } 
00299                                 
00300                                 for (j=0; j<4; ++j) {
00301                                         
00302                                         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];
00303                                         
00304                                         if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 25, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00305                                 } 
00306                                 
00307                                 for (j=0; j<4; ++j) {
00308                                         
00309                                         in_file >> theCurrentTemp.entx[k][i].xavgbcn[j] >> theCurrentTemp.entx[k][i].xrmsbcn[j] >> theCurrentTemp.entx[k][i].xgx0bcn[j] >> theCurrentTemp.entx[k][i].xgsigbcn[j];
00310                                         
00311                                         if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 26, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00312                                 } 
00313                                 
00314                                 in_file >> theCurrentTemp.entx[k][i].chi2xavgone >> theCurrentTemp.entx[k][i].chi2xminone >> theCurrentTemp.entx[k][i].qmin2
00315                                 >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav 
00316                                 >> theCurrentTemp.entx[k][i].mpvvav2 >> theCurrentTemp.entx[k][i].sigmavav2 >> theCurrentTemp.entx[k][i].kappavav2 >> theCurrentTemp.entx[k][i].spare[0];
00317                                 
00318                                 if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 27, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00319                                 
00320                                 in_file >> theCurrentTemp.entx[k][i].qbfrac[0] >> theCurrentTemp.entx[k][i].qbfrac[1] >> theCurrentTemp.entx[k][i].qbfrac[2] >> theCurrentTemp.entx[k][i].fracxone
00321                                 >> theCurrentTemp.entx[k][i].spare[1] >> theCurrentTemp.entx[k][i].spare[2] >> theCurrentTemp.entx[k][i].spare[3] >> theCurrentTemp.entx[k][i].spare[4]
00322                                 >> theCurrentTemp.entx[k][i].spare[5] >> theCurrentTemp.entx[k][i].spare[6];
00323                                 
00324                                 if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 28, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00325                                 
00326 
00327                         }
00328                 }       
00329                 
00330                 
00331                 in_file.close();
00332                 
00333                 // Add this template to the store
00334                 
00335                 theStripTemp_.push_back(theCurrentTemp);
00336                 
00337                 return true;
00338                 
00339         } else {
00340                 
00341                 // If file didn't open, report this
00342                 
00343                 LOGERROR("SiStripTemplate") << "Error opening File " << tempfile << ENDL;
00344                 return false;
00345                 
00346         }
00347         
00348 } // TempInit 
00349 
00350 
00351 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00352 
00353 //**************************************************************** 
00357 //**************************************************************** 
00358 bool SiStripTemplate::pushfile(const SiPixelTemplateDBObject& dbobject)
00359 {
00360         // Add template stored in external dbobject to theTemplateStore
00361     
00362         // Local variables 
00363         int i, j, k, l;
00364         float qavg_avg;
00365         //      const char *tempfile;
00366         const int code_version={17};
00367         
00368         // We must create a new object because dbobject must be a const and our stream must not be
00369         SiPixelTemplateDBObject db = dbobject;
00370         
00371         // Create a local template storage entry
00372         SiStripTemplateStore theCurrentTemp;
00373         
00374         // Fill the template storage for each template calibration stored in the db
00375         for(int m=0; m<db.numOfTempl(); ++m)
00376         {
00377                 
00378                 // Read-in a header string first and print it    
00379                 
00380                 SiPixelTemplateDBObject::char2float temp;
00381                 for (i=0; i<20; ++i) {
00382                         temp.f = db.sVector()[db.index()];
00383                         theCurrentTemp.head.title[4*i] = temp.c[0];
00384                         theCurrentTemp.head.title[4*i+1] = temp.c[1];
00385                         theCurrentTemp.head.title[4*i+2] = temp.c[2];
00386                         theCurrentTemp.head.title[4*i+3] = temp.c[3];
00387                         db.incrementIndex(1);
00388                 }
00389                 theCurrentTemp.head.title[79] = '\0';
00390                 LOGINFO("SiStripTemplate") << "Loading Strip Template File - " << theCurrentTemp.head.title << ENDL;
00391                 
00392                 // next, the header information     
00393                 
00394                 db >> theCurrentTemp.head.ID  >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx
00395                 >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
00396                 >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >> theCurrentTemp.head.zsize;
00397                 
00398                 if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file, no template load" << ENDL; return false;}
00399                 
00400                 LOGINFO("SiStripTemplate") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version " << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield 
00401                 << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
00402                 << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
00403                 << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
00404                 << ", 1/2 threshold " << theCurrentTemp.head.s50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", x Lorentz width " << theCurrentTemp.head.lorxwidth    
00405                 << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
00406                 
00407                 if(theCurrentTemp.head.templ_version < code_version) {LOGERROR("SiStripTemplate") << "code expects version " << code_version << ", no template load" << ENDL; return false;}
00408                 
00409                 
00410 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST 
00411                 
00412 // next, layout the 1-d/2-d structures needed to store template
00413                 
00414                 theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
00415                 
00416                 theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
00417                 
00418 #endif
00419                                 
00420                 // next, loop over all y-angle entries   
00421                 
00422                 for (i=0; i < theCurrentTemp.head.NTy; ++i) {     
00423                         
00424                         db >> theCurrentTemp.enty[i].runnum >> theCurrentTemp.enty[i].costrk[0] 
00425                         >> theCurrentTemp.enty[i].costrk[1] >> theCurrentTemp.enty[i].costrk[2]; 
00426                         
00427                         if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 1, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00428                         
00429                         // Calculate the alpha, beta, and cot(beta) for this entry 
00430                         
00431                         theCurrentTemp.enty[i].alpha = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[0]));
00432                         
00433                         theCurrentTemp.enty[i].cotalpha = theCurrentTemp.enty[i].costrk[0]/theCurrentTemp.enty[i].costrk[2];
00434                         
00435                         theCurrentTemp.enty[i].beta = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[1]));
00436                         
00437                         theCurrentTemp.enty[i].cotbeta = theCurrentTemp.enty[i].costrk[1]/theCurrentTemp.enty[i].costrk[2];
00438                         
00439                         db >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].sxmax >> theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].clslenx;
00440                         
00441                         if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 2, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00442                         
00443                         for (j=0; j<2; ++j) {
00444                                 
00445                                 db >> theCurrentTemp.enty[i].xpar[j][0] >> theCurrentTemp.enty[i].xpar[j][1] 
00446                                 >> theCurrentTemp.enty[i].xpar[j][2] >> theCurrentTemp.enty[i].xpar[j][3] >> theCurrentTemp.enty[i].xpar[j][4];
00447                                 
00448                                 if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 6, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00449                                 
00450                         }
00451                         
00452                         qavg_avg = 0.f;
00453                         for (j=0; j<9; ++j) {
00454                                 
00455                                 for (k=0; k<TSXSIZE; ++k) {db >> theCurrentTemp.enty[i].xtemp[j][k]; qavg_avg += theCurrentTemp.enty[i].xtemp[j][k];} 
00456                                 
00457                                 if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 7, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00458                         }
00459                         theCurrentTemp.enty[i].qavg_avg = qavg_avg/9.;
00460                         
00461                         for (j=0; j<4; ++j) {
00462                                 
00463                                 db >> theCurrentTemp.enty[i].xavg[j] >> theCurrentTemp.enty[i].xrms[j] >> theCurrentTemp.enty[i].xgx0[j] >> theCurrentTemp.enty[i].xgsig[j];
00464                                 
00465                                 if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 10, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00466                         }
00467                         
00468                         for (j=0; j<4; ++j) {
00469                                 
00470                                 db >> theCurrentTemp.enty[i].xflpar[j][0] >> theCurrentTemp.enty[i].xflpar[j][1] >> theCurrentTemp.enty[i].xflpar[j][2] 
00471                                 >> theCurrentTemp.enty[i].xflpar[j][3] >> theCurrentTemp.enty[i].xflpar[j][4] >> theCurrentTemp.enty[i].xflpar[j][5];
00472                                 
00473                                 if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 11, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00474                         }
00475                         
00476                         for (j=0; j<4; ++j) {
00477                                 
00478                                 db >> theCurrentTemp.enty[i].chi2xavg[j] >> theCurrentTemp.enty[i].chi2xmin[j] >> theCurrentTemp.enty[i].chi2xavgc2m[j] >> theCurrentTemp.enty[i].chi2xminc2m[j];
00479                                 
00480                                 if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 12, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00481                         }
00482                         
00483                         for (j=0; j<4; ++j) {
00484                                 
00485                                 db >> theCurrentTemp.enty[i].xavgc2m[j] >> theCurrentTemp.enty[i].xrmsc2m[j] >> theCurrentTemp.enty[i].xgx0c2m[j] >> theCurrentTemp.enty[i].xgsigc2m[j];
00486                                 
00487                                 if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 14, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00488                         } 
00489                         
00490                         for (j=0; j<4; ++j) {
00491                                 
00492                                 db >> theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j] >> theCurrentTemp.enty[i].xgx0gen[j] >> theCurrentTemp.enty[i].xgsiggen[j];
00493                                 
00494                                 if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 14b, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00495                         } 
00496                         
00497                         for (j=0; j<4; ++j) {
00498                                 
00499                                 db >> theCurrentTemp.enty[i].xavgbcn[j] >> theCurrentTemp.enty[i].xrmsbcn[j] >> theCurrentTemp.enty[i].xgx0bcn[j] >> theCurrentTemp.enty[i].xgsigbcn[j];
00500                                 
00501                                 if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 14c, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00502                         } 
00503                         
00504                         db >> theCurrentTemp.enty[i].chi2xavgone >> theCurrentTemp.enty[i].chi2xminone >> theCurrentTemp.enty[i].qmin2
00505                         >> theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav 
00506                         >> theCurrentTemp.enty[i].mpvvav2 >> theCurrentTemp.enty[i].sigmavav2 >> theCurrentTemp.enty[i].kappavav2 >> theCurrentTemp.enty[i].spare[0];
00507                         
00508                         if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 15, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00509                         
00510                         db >> theCurrentTemp.enty[i].qbfrac[0] >> theCurrentTemp.enty[i].qbfrac[1] >> theCurrentTemp.enty[i].qbfrac[2] >> theCurrentTemp.enty[i].fracxone
00511                         >> theCurrentTemp.enty[i].spare[1] >> theCurrentTemp.enty[i].spare[2] >> theCurrentTemp.enty[i].spare[3] >> theCurrentTemp.enty[i].spare[4]
00512                         >> theCurrentTemp.enty[i].spare[5] >> theCurrentTemp.enty[i].spare[6];
00513                         
00514                         if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 16, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
00515                         
00516                 }
00517                 
00518                 // next, loop over all barrel x-angle entries   
00519                 
00520                 for (k=0; k < theCurrentTemp.head.NTyx; ++k) { 
00521                         
00522                         for (i=0; i < theCurrentTemp.head.NTxx; ++i) { 
00523                                 
00524                                 db >> theCurrentTemp.entx[k][i].runnum >> theCurrentTemp.entx[k][i].costrk[0] 
00525                                 >> theCurrentTemp.entx[k][i].costrk[1] >> theCurrentTemp.entx[k][i].costrk[2]; 
00526                                 
00527                                 if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 17, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00528                                 
00529                                 // Calculate the alpha, beta, and cot(beta) for this entry 
00530                                 
00531                                 theCurrentTemp.entx[k][i].alpha = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[0]));
00532                                 
00533                                 theCurrentTemp.entx[k][i].cotalpha = theCurrentTemp.entx[k][i].costrk[0]/theCurrentTemp.entx[k][i].costrk[2];
00534                                 
00535                                 theCurrentTemp.entx[k][i].beta = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[1]));
00536                                 
00537                                 theCurrentTemp.entx[k][i].cotbeta = theCurrentTemp.entx[k][i].costrk[1]/theCurrentTemp.entx[k][i].costrk[2];
00538                                 
00539                                 db >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].sxmax >> theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone >> theCurrentTemp.entx[k][i].qmin >> theCurrentTemp.entx[k][i].clslenx;
00540                                 
00541                                 if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 18, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00542                                 
00543                                 for (j=0; j<2; ++j) {
00544                                         
00545                                         db >> theCurrentTemp.entx[k][i].xpar[j][0] >> theCurrentTemp.entx[k][i].xpar[j][1] 
00546                                         >> theCurrentTemp.entx[k][i].xpar[j][2] >> theCurrentTemp.entx[k][i].xpar[j][3] >> theCurrentTemp.entx[k][i].xpar[j][4];
00547                                         
00548                                         if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 19, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00549                                         
00550                                 }
00551                                 
00552                                 qavg_avg = 0.f;
00553                                 for (j=0; j<9; ++j) {
00554                                         
00555                                         for (l=0; l<TSXSIZE; ++k) {db >> theCurrentTemp.entx[k][i].xtemp[j][l]; qavg_avg += theCurrentTemp.entx[k][i].xtemp[j][l];} 
00556                                         
00557                                         if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 20, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00558                                 }
00559                                 theCurrentTemp.entx[k][i].qavg_avg = qavg_avg/9.;
00560                                 
00561                                 for (j=0; j<4; ++j) {
00562                                         
00563                                         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];
00564                                         
00565                                         if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 21, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00566                                 }
00567                                 
00568                                 for (j=0; j<4; ++j) {
00569                                         
00570                                         db >> theCurrentTemp.entx[k][i].xflpar[j][0] >> theCurrentTemp.entx[k][i].xflpar[j][1] >> theCurrentTemp.entx[k][i].xflpar[j][2] 
00571                                         >> theCurrentTemp.entx[k][i].xflpar[j][3] >> theCurrentTemp.entx[k][i].xflpar[j][4] >> theCurrentTemp.entx[k][i].xflpar[j][5];
00572                                         
00573                                         if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 22, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00574                                 }
00575                                 
00576                                 for (j=0; j<4; ++j) {
00577                                         
00578                                         db >> theCurrentTemp.entx[k][i].chi2xavg[j] >> theCurrentTemp.entx[k][i].chi2xmin[j] >> theCurrentTemp.entx[k][i].chi2xavgc2m[j] >> theCurrentTemp.entx[k][i].chi2xminc2m[j];
00579                                         
00580                                         if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 23, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00581                                 }
00582                                 
00583                                 for (j=0; j<4; ++j) {
00584                                         
00585                                         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];
00586                                         
00587                                         if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 24, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00588                                 } 
00589                                 
00590                                 for (j=0; j<4; ++j) {
00591                                         
00592                                         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];
00593                                         
00594                                         if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 25, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00595                                 } 
00596                                 
00597                                 for (j=0; j<4; ++j) {
00598                                         
00599                                         db >> theCurrentTemp.entx[k][i].xavgbcn[j] >> theCurrentTemp.entx[k][i].xrmsbcn[j] >> theCurrentTemp.entx[k][i].xgx0bcn[j] >> theCurrentTemp.entx[k][i].xgsigbcn[j];
00600                                         
00601                                         if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 26, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00602                                 } 
00603                                 
00604                                 db >> theCurrentTemp.entx[k][i].chi2xavgone >> theCurrentTemp.entx[k][i].chi2xminone >> theCurrentTemp.entx[k][i].qmin2
00605                                 >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav 
00606                                 >> theCurrentTemp.entx[k][i].mpvvav2 >> theCurrentTemp.entx[k][i].sigmavav2 >> theCurrentTemp.entx[k][i].kappavav2 >> theCurrentTemp.entx[k][i].spare[0];
00607                                 
00608                                 if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 27, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00609                                 
00610                                 db >> theCurrentTemp.entx[k][i].qbfrac[0] >> theCurrentTemp.entx[k][i].qbfrac[1] >> theCurrentTemp.entx[k][i].qbfrac[2] >> theCurrentTemp.entx[k][i].fracxone
00611                                 >> theCurrentTemp.entx[k][i].spare[1] >> theCurrentTemp.entx[k][i].spare[2] >> theCurrentTemp.entx[k][i].spare[3] >> theCurrentTemp.entx[k][i].spare[4]
00612                                 >> theCurrentTemp.entx[k][i].spare[5] >> theCurrentTemp.entx[k][i].spare[6];
00613                                 
00614                                 if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 28, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
00615                                 
00616                                 
00617                                 
00618                         }
00619                 }       
00620                 
00621                 
00622                 // Add this template to the store
00623                 
00624                 theStripTemp_.push_back(theCurrentTemp);
00625                 
00626         }
00627         return true;
00628         
00629 } // TempInit 
00630 
00631 #endif
00632 
00633 
00634 // ************************************************************************************************************ 
00640 // ************************************************************************************************************ 
00641 bool SiStripTemplate::interpolate(int id, float cotalpha, float cotbeta, float locBy)
00642 {
00643     // Interpolate for a new set of track angles 
00644     
00645     // Local variables 
00646     int i, j;
00647         int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, imidy, imaxx;
00648         float yratio, yxratio, xxratio, sxmax, qcorrect, qxtempcor, chi2xavgone, chi2xminone, cota, cotb, cotalpha0, cotbeta0;
00649         bool flip_x;
00650 //      std::vector <float> xrms(4), xgsig(4), xrmsc2m(4), xgsigc2m(4);
00651         std::vector <float> chi2xavg(4), chi2xmin(4), chi2xavgc2m(4), chi2xminc2m(4);
00652 
00653 
00654 // Check to see if interpolation is valid     
00655 
00656 if(id != id_current_ || cotalpha != cota_current_ || cotbeta != cotb_current_) {
00657 
00658         cota_current_ = cotalpha; cotb_current_ = cotbeta; success_ = true;
00659         
00660         if(id != id_current_) {
00661 
00662 // Find the index corresponding to id
00663 
00664        index_id_ = -1;
00665        for(i=0; i<(int)theStripTemp_.size(); ++i) {
00666         
00667               if(id == theStripTemp_[i].head.ID) {
00668            
00669                  index_id_ = i;
00670                       id_current_ = id;
00671                                 
00672 // Copy the charge scaling factor to the private variable     
00673                                 
00674                                 qscale_ = theStripTemp_[index_id_].head.qscale;
00675                                 
00676 // Copy the pseudopixel signal size to the private variable     
00677                                 
00678                                 s50_ = theStripTemp_[index_id_].head.s50;
00679                         
00680 // Pixel sizes to the private variables     
00681                                 
00682                                 xsize_ = theStripTemp_[index_id_].head.xsize;
00683                                 ysize_ = theStripTemp_[index_id_].head.ysize;
00684                                 zsize_ = theStripTemp_[index_id_].head.zsize;
00685                                 
00686                                 break;
00687           }
00688             }
00689      }
00690          
00691 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00692     if(index_id_ < 0 || index_id_ >= (int)theStripTemp_.size()) {
00693                 throw cms::Exception("DataCorrupt") << "SiStripTemplate::interpolate can't find needed template ID = " << id << std::endl;
00694         }
00695 #else
00696         assert(index_id_ >= 0 && index_id_ < (int)theStripTemp_.size());
00697 #endif
00698          
00699 // Interpolate the absolute value of cot(beta)     
00700     
00701         abs_cotb_ = std::abs(cotbeta);
00702         cotb = abs_cotb_; 
00703         
00704 //      qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)       
00705 
00706         cotalpha0 =  theStripTemp_[index_id_].enty[0].cotalpha;
00707         qcorrect=std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta*cotbeta+cotalpha0*cotalpha0));  
00708 // flip quantities when the magnetic field in in the positive y local direction  
00709         if(locBy > 0.f) {
00710                 flip_x = true;
00711         } else {
00712                 flip_x = false;
00713         }
00714                 
00715         Ny = theStripTemp_[index_id_].head.NTy;
00716         Nyx = theStripTemp_[index_id_].head.NTyx;
00717         Nxx = theStripTemp_[index_id_].head.NTxx;
00718                 
00719 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00720         if(Ny < 2 || Nyx < 1 || Nxx < 2) {
00721                 throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
00722         }
00723 #else
00724         assert(Ny > 1 && Nyx > 0 && Nxx > 1);
00725 #endif
00726         imaxx = Nyx - 1;
00727         imidy = Nxx/2;
00728         
00729 // next, loop over all y-angle entries   
00730 
00731         ilow = 0;
00732         yratio = 0.f;
00733 
00734         if(cotb >= theStripTemp_[index_id_].enty[Ny-1].cotbeta) {
00735         
00736                 ilow = Ny-2;
00737                 yratio = 1.f;
00738                 success_ = false;
00739                 
00740         } else {
00741            
00742                 if(cotb >= theStripTemp_[index_id_].enty[0].cotbeta) {
00743 
00744                         for (i=0; i<Ny-1; ++i) { 
00745     
00746                         if( theStripTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < theStripTemp_[index_id_].enty[i+1].cotbeta) {
00747                   
00748                                 ilow = i;
00749                                 yratio = (cotb - theStripTemp_[index_id_].enty[i].cotbeta)/(theStripTemp_[index_id_].enty[i+1].cotbeta - theStripTemp_[index_id_].enty[i].cotbeta);
00750                                 break;                   
00751                         }
00752                 }
00753         } else { success_ = false; }
00754   }
00755         
00756         ihigh=ilow + 1;
00757                           
00758 // Interpolate/store all y-related quantities (flip displacements when flip_y)
00759 
00760         yratio_ = yratio;
00761         qavg_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].qavg + yratio*theStripTemp_[index_id_].enty[ihigh].qavg;
00762         qavg_ *= qcorrect;
00763         sxmax = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].sxmax + yratio*theStripTemp_[index_id_].enty[ihigh].sxmax;
00764         syparmax_ = sxmax;
00765         qmin_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].qmin + yratio*theStripTemp_[index_id_].enty[ihigh].qmin;
00766         qmin_ *= qcorrect;
00767         qmin2_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].qmin2 + yratio*theStripTemp_[index_id_].enty[ihigh].qmin2;
00768         qmin2_ *= qcorrect;
00769         mpvvav_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].mpvvav + yratio*theStripTemp_[index_id_].enty[ihigh].mpvvav;
00770         mpvvav_ *= qcorrect;
00771         sigmavav_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].sigmavav + yratio*theStripTemp_[index_id_].enty[ihigh].sigmavav;
00772         kappavav_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].kappavav + yratio*theStripTemp_[index_id_].enty[ihigh].kappavav;
00773         mpvvav2_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].mpvvav2 + yratio*theStripTemp_[index_id_].enty[ihigh].mpvvav2;
00774         mpvvav2_ *= qcorrect;
00775         sigmavav2_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].sigmavav2 + yratio*theStripTemp_[index_id_].enty[ihigh].sigmavav2;
00776         kappavav2_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].kappavav2 + yratio*theStripTemp_[index_id_].enty[ihigh].kappavav2;
00777         qavg_avg_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].qavg_avg + yratio*theStripTemp_[index_id_].enty[ihigh].qavg_avg;
00778         qavg_avg_ *= qcorrect;
00779         for(i=0; i<2 ; ++i) {
00780                 for(j=0; j<5 ; ++j) {
00781 // Charge loss switches sides when cot(alpha) changes sign
00782                         if(flip_x) {
00783                                 xparly0_[1-i][j] = theStripTemp_[index_id_].enty[ilow].xpar[i][j];
00784                                 xparhy0_[1-i][j] = theStripTemp_[index_id_].enty[ihigh].xpar[i][j];
00785                         } else {
00786                                 xparly0_[i][j] = theStripTemp_[index_id_].enty[ilow].xpar[i][j];
00787                                 xparhy0_[i][j] = theStripTemp_[index_id_].enty[ihigh].xpar[i][j];
00788                         }
00789                 }
00790         }
00791         for(i=0; i<4; ++i) {
00792                 chi2xavg[i]=(1.f - yratio)*theStripTemp_[index_id_].enty[ilow].chi2xavg[i] + yratio*theStripTemp_[index_id_].enty[ihigh].chi2xavg[i];
00793                 chi2xmin[i]=(1.f - yratio)*theStripTemp_[index_id_].enty[ilow].chi2xmin[i] + yratio*theStripTemp_[index_id_].enty[ihigh].chi2xmin[i];
00794                 chi2xavgc2m[i]=(1.f - yratio)*theStripTemp_[index_id_].enty[ilow].chi2xavgc2m[i] + yratio*theStripTemp_[index_id_].enty[ihigh].chi2xavgc2m[i];
00795                 chi2xminc2m[i]=(1.f - yratio)*theStripTemp_[index_id_].enty[ilow].chi2xminc2m[i] + yratio*theStripTemp_[index_id_].enty[ihigh].chi2xminc2m[i];
00796         }
00797            
00799 
00800         chi2xavgone=(1.f - yratio)*theStripTemp_[index_id_].enty[ilow].chi2xavgone + yratio*theStripTemp_[index_id_].enty[ihigh].chi2xavgone;
00801         chi2xminone=(1.f - yratio)*theStripTemp_[index_id_].enty[ilow].chi2xminone + yratio*theStripTemp_[index_id_].enty[ihigh].chi2xminone;
00802                 //       for(i=0; i<10; ++i) {
00803 //                  pyspare[i]=(1.f - yratio)*theStripTemp_[index_id_].enty[ilow].yspare[i] + yratio*theStripTemp_[index_id_].enty[ihigh].yspare[i];
00804 //       }
00805                           
00806         
00807 // next, loop over all x-angle entries, first, find relevant y-slices   
00808         
00809         iylow = 0;
00810         yxratio = 0.f;
00811 
00812         if(abs_cotb_ >= theStripTemp_[index_id_].entx[Nyx-1][0].cotbeta) {
00813         
00814                 iylow = Nyx-2;
00815                 yxratio = 1.f;
00816                 
00817         } else if(abs_cotb_ >= theStripTemp_[index_id_].entx[0][0].cotbeta) {
00818 
00819                 for (i=0; i<Nyx-1; ++i) { 
00820     
00821                         if( theStripTemp_[index_id_].entx[i][0].cotbeta <= abs_cotb_ && abs_cotb_ < theStripTemp_[index_id_].entx[i+1][0].cotbeta) {
00822                   
00823                            iylow = i;
00824                            yxratio = (abs_cotb_ - theStripTemp_[index_id_].entx[i][0].cotbeta)/(theStripTemp_[index_id_].entx[i+1][0].cotbeta - theStripTemp_[index_id_].entx[i][0].cotbeta);
00825                            break;                        
00826                         }
00827                 }
00828         }
00829         
00830         iyhigh=iylow + 1;
00831 
00832         ilow = 0;
00833         xxratio = 0.f;
00834         if(flip_x) {cota = -cotalpha;} else {cota = cotalpha;}
00835 
00836         if(cota >= theStripTemp_[index_id_].entx[0][Nxx-1].cotalpha) {
00837         
00838                 ilow = Nxx-2;
00839                 xxratio = 1.f;
00840                 success_ = false;
00841                 
00842         } else {
00843            
00844                 if(cota >= theStripTemp_[index_id_].entx[0][0].cotalpha) {
00845 
00846                         for (i=0; i<Nxx-1; ++i) { 
00847     
00848                            if( theStripTemp_[index_id_].entx[0][i].cotalpha <= cota && cota < theStripTemp_[index_id_].entx[0][i+1].cotalpha) {
00849                   
00850                                   ilow = i;
00851                                   xxratio = (cota - theStripTemp_[index_id_].entx[0][i].cotalpha)/(theStripTemp_[index_id_].entx[0][i+1].cotalpha - theStripTemp_[index_id_].entx[0][i].cotalpha);
00852                                   break;
00853                         }
00854                   }
00855                 } else { success_ = false; }
00856         }
00857         
00858         ihigh=ilow + 1;
00859                           
00860 // Interpolate/store all x-related quantities 
00861 
00862         yxratio_ = yxratio;
00863         xxratio_ = xxratio;             
00864                           
00865 // sxparmax defines the maximum charge for which the parameters xpar are defined (not rescaled by cotbeta) 
00866 
00867         sxparmax_ = (1.f - xxratio)*theStripTemp_[index_id_].entx[imaxx][ilow].sxmax + xxratio*theStripTemp_[index_id_].entx[imaxx][ihigh].sxmax;
00868         sxmax_ = sxparmax_;
00869         if(theStripTemp_[index_id_].entx[imaxx][imidy].sxmax != 0.f) {sxmax_=sxmax_/theStripTemp_[index_id_].entx[imaxx][imidy].sxmax*sxmax;}
00870         dxone_ = (1.f - xxratio)*theStripTemp_[index_id_].entx[0][ilow].dxone + xxratio*theStripTemp_[index_id_].entx[0][ihigh].dxone;
00871         if(flip_x) {dxone_ = -dxone_;}
00872         sxone_ = (1.f - xxratio)*theStripTemp_[index_id_].entx[0][ilow].sxone + xxratio*theStripTemp_[index_id_].entx[0][ihigh].sxone;
00873         clslenx_ = fminf(theStripTemp_[index_id_].entx[0][ilow].clslenx, theStripTemp_[index_id_].entx[0][ihigh].clslenx);
00874         for(i=0; i<2 ; ++i) {
00875                 for(j=0; j<5 ; ++j) {
00876                         if(flip_x) {
00877                  xpar0_[1-i][j] = theStripTemp_[index_id_].entx[imaxx][imidy].xpar[i][j];
00878                  xparl_[1-i][j] = theStripTemp_[index_id_].entx[imaxx][ilow].xpar[i][j];
00879                  xparh_[1-i][j] = theStripTemp_[index_id_].entx[imaxx][ihigh].xpar[i][j];
00880                         } else {
00881                  xpar0_[i][j] = theStripTemp_[index_id_].entx[imaxx][imidy].xpar[i][j];
00882                  xparl_[i][j] = theStripTemp_[index_id_].entx[imaxx][ilow].xpar[i][j];
00883                  xparh_[i][j] = theStripTemp_[index_id_].entx[imaxx][ihigh].xpar[i][j];
00884                         }
00885                 }
00886         }
00887                           
00888 // sxmax is the maximum allowed strip charge (used for truncation)
00889 
00890         sxmax_=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].sxmax + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].sxmax)
00891                         +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].sxmax + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].sxmax);
00892                           
00893         for(i=0; i<4; ++i) {
00894                 xavg_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xavg[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xavg[i])
00895                                 +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xavg[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xavg[i]);
00896                 if(flip_x) {xavg_[i] = -xavg_[i];}
00897                   
00898                 xrms_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xrms[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xrms[i])
00899                                 +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xrms[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xrms[i]);
00900                   
00901 //            xgx0_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xgx0[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xgx0[i])
00902 //                        +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xgx0[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xgx0[i]);
00903                                                         
00904 //            xgsig_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xgsig[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xgsig[i])
00905 //                        +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xgsig[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xgsig[i]);
00906                                   
00907                 xavgc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xavgc2m[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xavgc2m[i])
00908                                 +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xavgc2m[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xavgc2m[i]);
00909                 if(flip_x) {xavgc2m_[i] = -xavgc2m_[i];}
00910                   
00911                 xrmsc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xrmsc2m[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xrmsc2m[i])
00912                                 +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xrmsc2m[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xrmsc2m[i]);
00913                   
00914                 xavgbcn_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xavgbcn[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xavgbcn[i])
00915         +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xavgbcn[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xavgbcn[i]);
00916                 if(flip_x) {xavgbcn_[i] = -xavgbcn_[i];}
00917         
00918                 xrmsbcn_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xrmsbcn[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xrmsbcn[i])
00919         +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xrmsbcn[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xrmsbcn[i]);
00920         
00921 //            xgx0c2m_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xgx0c2m[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xgx0c2m[i])
00922 //                        +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xgx0c2m[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xgx0c2m[i]);
00923                                                         
00924 //            xgsigc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xgsigc2m[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xgsigc2m[i])
00925 //                        +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xgsigc2m[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xgsigc2m[i]);
00926 //
00927 //  Try new interpolation scheme
00928 //                                                                                                                      
00929 //            chi2xavg_[i]=((1.f - xxratio)*theStripTemp_[index_id_].entx[imaxx][ilow].chi2xavg[i] + xxratio*theStripTemp_[index_id_].entx[imaxx][ihigh].chi2xavg[i]);
00930 //                if(theStripTemp_[index_id_].entx[imaxx][imidy].chi2xavg[i] != 0.f) {chi2xavg_[i]=chi2xavg_[i]/theStripTemp_[index_id_].entx[imaxx][imidy].chi2xavg[i]*chi2xavg[i];}
00931 //                                                      
00932 //            chi2xmin_[i]=((1.f - xxratio)*theStripTemp_[index_id_].entx[imaxx][ilow].chi2xmin[i] + xxratio*theStripTemp_[index_id_].entx[imaxx][ihigh].chi2xmin[i]);
00933 //                if(theStripTemp_[index_id_].entx[imaxx][imidy].chi2xmin[i] != 0.f) {chi2xmin_[i]=chi2xmin_[i]/theStripTemp_[index_id_].entx[imaxx][imidy].chi2xmin[i]*chi2xmin[i];}
00934 //                
00935                 chi2xavg_[i]=((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].chi2xavg[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].chi2xavg[i]);
00936                 if(theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xavg[i] != 0.f) {chi2xavg_[i]=chi2xavg_[i]/theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xavg[i]*chi2xavg[i];}
00937                                                         
00938                 chi2xmin_[i]=((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].chi2xmin[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].chi2xmin[i]);
00939                 if(theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xmin[i] != 0.f) {chi2xmin_[i]=chi2xmin_[i]/theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xmin[i]*chi2xmin[i];}
00940 
00941                 chi2xavgc2m_[i]=((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].chi2xavgc2m[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].chi2xavgc2m[i]);
00942                 if(theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xavgc2m[i] != 0.f) {chi2xavgc2m_[i]=chi2xavgc2m_[i]/theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xavgc2m[i]*chi2xavgc2m[i];}
00943                 
00944                 chi2xminc2m_[i]=((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].chi2xminc2m[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].chi2xminc2m[i]);
00945                 if(theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xminc2m[i] != 0.f) {chi2xminc2m_[i]=chi2xminc2m_[i]/theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xminc2m[i]*chi2xminc2m[i];}    
00946                 
00947                 for(j=0; j<6 ; ++j) {
00948                  xflparll_[i][j] = theStripTemp_[index_id_].entx[iylow][ilow].xflpar[i][j];
00949                  xflparlh_[i][j] = theStripTemp_[index_id_].entx[iylow][ihigh].xflpar[i][j];
00950                  xflparhl_[i][j] = theStripTemp_[index_id_].entx[iyhigh][ilow].xflpar[i][j];
00951                  xflparhh_[i][j] = theStripTemp_[index_id_].entx[iyhigh][ihigh].xflpar[i][j];
00952                         // Since Q_fl is odd under cotbeta, it flips qutomatically, change only even terms
00953                         
00954                         if(flip_x && (j == 0 || j == 2 || j == 4)) {
00955                  xflparll_[i][j] = -xflparll_[i][j];
00956                  xflparlh_[i][j] = -xflparlh_[i][j];
00957                  xflparhl_[i][j] = -xflparhl_[i][j];
00958                  xflparhh_[i][j] = -xflparhh_[i][j];
00959                         }
00960                 }
00961         }
00962            
00963 // Do the spares next
00964 
00965         chi2xavgone_=((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].chi2xavgone + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].chi2xavgone);
00966         if(theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xavgone != 0.f) {chi2xavgone_=chi2xavgone_/theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xavgone*chi2xavgone;}
00967                 
00968         chi2xminone_=((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].chi2xminone + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].chi2xminone);
00969         if(theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xminone != 0.f) {chi2xminone_=chi2xminone_/theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xminone*chi2xminone;}
00970                 //       for(i=0; i<10; ++i) {
00971 //            pxspare[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xspare[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xspare[i])
00972 //                        +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xspare[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xspare[i]);
00973 //       }
00974                           
00975 // Interpolate and build the x-template 
00976         
00977 //      qxtempcor corrects the total charge to the actual track angles (not actually needed for the template fits, but useful for Guofan)
00978         
00979         cotbeta0 =  theStripTemp_[index_id_].entx[iyhigh][0].cotbeta;
00980         qxtempcor=std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta0*cotbeta0+cotalpha*cotalpha));
00981         
00982         for(i=0; i<9; ++i) {
00983                 xtemp_[i][0] = 0.f;
00984                 xtemp_[i][1] = 0.f;
00985                 xtemp_[i][BSXM2] = 0.f;
00986                 xtemp_[i][BSXM1] = 0.f;
00987                 for(j=0; j<TSXSIZE; ++j) {
00988 //  Take next largest x-slice for the x-template (it reduces bias in the forward direction after irradiation)
00989 //                 xtemp_[i][j+2]=(1.f - xxratio)*theStripTemp_[index_id_].entx[imaxx][ilow].xtemp[i][j] + xxratio*theStripTemp_[index_id_].entx[imaxx][ihigh].xtemp[i][j];
00990          if(flip_x) {
00991                                 xtemp_[8-i][BSXM3-j]=qxtempcor*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xtemp[i][j] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xtemp[i][j]);
00992                         } else {
00993             xtemp_[i][j+2]=qxtempcor*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xtemp[i][j] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xtemp[i][j]);
00994                         }
00995                 }
00996         }
00997         
00998         lorxwidth_ = theStripTemp_[index_id_].head.lorxwidth;
00999         if(locBy > 0.f) {lorxwidth_ = -lorxwidth_;}
01000         
01001   }
01002         
01003   return success_;
01004 } // interpolate
01005 
01006 
01007 
01008 // ************************************************************************************************************ 
01016 // ************************************************************************************************************ 
01017   void SiStripTemplate::xsigma2(int fxpix, int lxpix, float sxthr, float xsum[11], float xsig2[11])
01018   
01019 {
01020     // Interpolate using quantities already stored in the private variables
01021     
01022     // Local variables 
01023     int i;
01024         float sigi, sigi2, sigi3, sigi4, yint, sxmax, x0, qscale;
01025         float sigiy, sigiy2, sigiy3, sigiy4;
01026         
01027     // Make sure that input is OK
01028     
01029 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01030                   if(fxpix < 2 || fxpix >= BSXM2) {
01031                          throw cms::Exception("DataCorrupt") << "SiStripTemplate::xsigma2 called with fxpix = " << fxpix << std::endl;
01032                    }
01033 #else
01034                    assert(fxpix > 1 && fxpix < BSXM2);
01035 #endif
01036 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01037                          if(lxpix < fxpix || lxpix >= BSXM2) {
01038                                 throw cms::Exception("DataCorrupt") << "SiStripTemplate::xsigma2 called with lxpix/fxpix = " << lxpix << "/" << fxpix << std::endl;
01039                          }
01040 #else
01041                          assert(lxpix >= fxpix && lxpix < BSXM2);
01042 #endif
01043                      
01044 // Define the maximum signal to use in the parameterization 
01045 
01046        sxmax = sxmax_;
01047            if(sxmax_ > sxparmax_) {sxmax = sxparmax_;}
01048            
01049 // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis 
01050 
01051            for(i=fxpix-2; i<=lxpix+2; ++i) {
01052                   if(i < fxpix || i > lxpix) {
01053            
01054 // Nearest pseudopixels have uncertainties of 50% of threshold, next-nearest have 10% of threshold
01055 
01056                          xsig2[i] = s50_*s50_;
01057                   } else {
01058                          if(xsum[i] < sxmax) {
01059                                 sigi = xsum[i];
01060                                 qscale = 1.f;
01061                          } else {
01062                                 sigi = sxmax;
01063                                 qscale = xsum[i]/sxmax;
01064                          }
01065                          sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
01066                          
01067                           if(xsum[i] < syparmax_) {
01068                                   sigiy = xsum[i];
01069                           } else {
01070                                   sigiy = syparmax_;
01071                           }
01072                           sigiy2 = sigiy*sigiy; sigiy3 = sigiy2*sigiy; sigiy4 = sigiy3*sigiy;
01073                           
01074 // First, do the cotbeta interpolation                   
01075                          
01076                          if(i <= BSHX) {
01077                                 yint = (1.f-yratio_)*
01078                                 (xparly0_[0][0]+xparly0_[0][1]*sigiy+xparly0_[0][2]*sigiy2+xparly0_[0][3]*sigiy3+xparly0_[0][4]*sigiy4)
01079                                 + yratio_*
01080                                 (xparhy0_[0][0]+xparhy0_[0][1]*sigiy+xparhy0_[0][2]*sigiy2+xparhy0_[0][3]*sigiy3+xparhy0_[0][4]*sigiy4);
01081                          } else {
01082                                 yint = (1.f-yratio_)*
01083                                 (xparly0_[1][0]+xparly0_[1][1]*sigiy+xparly0_[1][2]*sigiy2+xparly0_[1][3]*sigiy3+xparly0_[1][4]*sigiy4)
01084                                 + yratio_*
01085                                 (xparhy0_[1][0]+xparhy0_[1][1]*sigiy+xparhy0_[1][2]*sigiy2+xparhy0_[1][3]*sigiy3+xparhy0_[1][4]*sigiy4);
01086                          }
01087                          
01088 // Next, do the cotalpha interpolation                   
01089                          
01090                          if(i <= BSHX) {
01091                                 xsig2[i] = (1.f-xxratio_)*
01092                                 (xparl_[0][0]+xparl_[0][1]*sigi+xparl_[0][2]*sigi2+xparl_[0][3]*sigi3+xparl_[0][4]*sigi4)
01093                                 + xxratio_*
01094                                 (xparh_[0][0]+xparh_[0][1]*sigi+xparh_[0][2]*sigi2+xparh_[0][3]*sigi3+xparh_[0][4]*sigi4);
01095                          } else {
01096                                 xsig2[i] = (1.f-xxratio_)*
01097                                 (xparl_[1][0]+xparl_[1][1]*sigi+xparl_[1][2]*sigi2+xparl_[1][3]*sigi3+xparl_[1][4]*sigi4)
01098                                 + xxratio_*
01099                             (xparh_[1][0]+xparh_[1][1]*sigi+xparh_[1][2]*sigi2+xparh_[1][3]*sigi3+xparh_[1][4]*sigi4);
01100                          }
01101                          
01102 // Finally, get the mid-point value of the cotalpha function                     
01103                          
01104                          if(i <= BSHX) {
01105                                 x0 = xpar0_[0][0]+xpar0_[0][1]*sigi+xpar0_[0][2]*sigi2+xpar0_[0][3]*sigi3+xpar0_[0][4]*sigi4;
01106                          } else {
01107                                 x0 = xpar0_[1][0]+xpar0_[1][1]*sigi+xpar0_[1][2]*sigi2+xpar0_[1][3]*sigi3+xpar0_[1][4]*sigi4;
01108                          }
01109                          
01110 // Finally, rescale the yint value for cotalpha variation                        
01111                          
01112                          if(x0 != 0.f) {xsig2[i] = xsig2[i]/x0 * yint;}
01113                          xsig2[i] *=qscale;
01114                      if(xsum[i] > sxthr) {xsig2[i] = 1.e8f;}
01115                          if(xsig2[i] <= 0.f) {LOGERROR("SiStripTemplate") << "neg x-error-squared = " << xsig2[i] << ", id = " << id_current_ << ", index = " << index_id_ << 
01116                          ", cot(alpha) = " << cota_current_ << ", cot(beta) = " << cotb_current_  << ", sigi = " << sigi << ", sxparmax = " << sxparmax_ << ", sxmax = " << sxmax_ << ENDL;}
01117               }
01118            }
01119         
01120         return;
01121         
01122 } // End xsigma2
01123 
01124 
01125 
01126 
01127 // ************************************************************************************************************ 
01131 // ************************************************************************************************************ 
01132   float SiStripTemplate::xflcorr(int binq, float qflx)
01133   
01134 {
01135     // Interpolate using quantities already stored in the private variables
01136     
01137     // Local variables 
01138         float qfl, qfl2, qfl3, qfl4, qfl5, dx;
01139         
01140     // Make sure that input is OK
01141     
01142 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01143         if(binq < 0 || binq > 3) {
01144                 throw cms::Exception("DataCorrupt") << "SiStripTemplate::xflcorr called with binq = " << binq << std::endl;
01145         }
01146 #else
01147         assert(binq >= 0 && binq < 4);
01148 #endif
01149 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01150         if(fabs((double)qflx) > 1.) {
01151                 throw cms::Exception("DataCorrupt") << "SiStripTemplate::xflcorr called with qflx = " << qflx << std::endl;
01152         }
01153 #else
01154         assert(fabs((double)qflx) <= 1.);
01155 #endif
01156                      
01157 // Define the maximum signal to allow before de-weighting a pixel 
01158 
01159        qfl = qflx;
01160 
01161        if(qfl < -0.9f) {qfl = -0.9f;}
01162            if(qfl > 0.9f) {qfl = 0.9f;}
01163            
01164 // Interpolate between the two polynomials
01165 
01166            qfl2 = qfl*qfl; qfl3 = qfl2*qfl; qfl4 = qfl3*qfl; qfl5 = qfl4*qfl;
01167            dx = (1.f - yxratio_)*((1.f-xxratio_)*(xflparll_[binq][0]+xflparll_[binq][1]*qfl+xflparll_[binq][2]*qfl2+xflparll_[binq][3]*qfl3+xflparll_[binq][4]*qfl4+xflparll_[binq][5]*qfl5)
01168                   + xxratio_*(xflparlh_[binq][0]+xflparlh_[binq][1]*qfl+xflparlh_[binq][2]*qfl2+xflparlh_[binq][3]*qfl3+xflparlh_[binq][4]*qfl4+xflparlh_[binq][5]*qfl5))
01169               + yxratio_*((1.f-xxratio_)*(xflparhl_[binq][0]+xflparhl_[binq][1]*qfl+xflparhl_[binq][2]*qfl2+xflparhl_[binq][3]*qfl3+xflparhl_[binq][4]*qfl4+xflparhl_[binq][5]*qfl5)
01170                   + xxratio_*(xflparhh_[binq][0]+xflparhh_[binq][1]*qfl+xflparhh_[binq][2]*qfl2+xflparhh_[binq][3]*qfl3+xflparhh_[binq][4]*qfl4+xflparhh_[binq][5]*qfl5));
01171         
01172         return dx;
01173         
01174 } // End xflcorr
01175 
01176 
01177 // ************************************************************************************************************ 
01182 // ************************************************************************************************************ 
01183   void SiStripTemplate::xtemp(int fxbin, int lxbin, float xtemplate[41][BSXSIZE])
01184   
01185 {
01186     // Retrieve already interpolated quantities
01187     
01188     // Local variables 
01189     int i, j;
01190 
01191    // Verify that input parameters are in valid range
01192 
01193 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01194         if(fxbin < 0 || fxbin > 40) {
01195                 throw cms::Exception("DataCorrupt") << "SiStripTemplate::xtemp called with fxbin = " << fxbin << std::endl;
01196         }
01197 #else
01198         assert(fxbin >= 0 && fxbin < 41);
01199 #endif
01200 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01201         if(lxbin < 0 || lxbin > 40) {
01202                 throw cms::Exception("DataCorrupt") << "SiStripTemplate::xtemp called with lxbin = " << lxbin << std::endl;
01203         }
01204 #else
01205         assert(lxbin >= 0 && lxbin < 41);
01206 #endif
01207 
01208 // Build the x-template, the central 25 bins are here in all cases
01209         
01210         for(i=0; i<9; ++i) {
01211            for(j=0; j<BSXSIZE; ++j) {
01212               xtemplate[i+16][j]=xtemp_[i][j];
01213            }
01214         }
01215         for(i=0; i<8; ++i) {
01216            xtemplate[i+8][BSXM1] = 0.f;
01217            for(j=0; j<BSXM1; ++j) {
01218               xtemplate[i+8][j]=xtemp_[i][j+1];
01219            }
01220         }
01221         for(i=1; i<9; ++i) {
01222            xtemplate[i+24][0] = 0.f;
01223            for(j=0; j<BSXM1; ++j) {
01224               xtemplate[i+24][j+1]=xtemp_[i][j];
01225            }
01226         }
01227         
01228 //  Add more bins if needed     
01229         
01230         if(fxbin < 8) {
01231            for(i=0; i<8; ++i) {
01232           xtemplate[i][BSXM2] = 0.f;
01233               xtemplate[i][BSXM1] = 0.f;
01234               for(j=0; j<BSXM2; ++j) {
01235                 xtemplate[i][j]=xtemp_[i][j+2];
01236               }
01237            }
01238         }
01239         if(lxbin > 32) {
01240            for(i=1; i<9; ++i) {
01241           xtemplate[i+32][0] = 0.f;
01242               xtemplate[i+32][1] = 0.f;
01243               for(j=0; j<BSXM2; ++j) {
01244                 xtemplate[i+32][j+2]=xtemp_[i][j];
01245               }
01246            }
01247         }
01248         
01249         return;
01250         
01251 } // End xtemp
01252 
01253 
01254 // ************************************************************************************************************ 
01258 // ************************************************************************************************************ 
01259 void SiStripTemplate::sxtemp(float xhit, std::vector<float>& cluster)
01260 
01261 {
01262         // Retrieve already interpolated quantities
01263         
01264         // Local variables 
01265         int i, j;
01266         
01267         // Extract x template based upon the hit position 
01268    
01269         float xpix = xhit/xsize_ + 0.5f;
01270 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01271         if(xpix < 0.f) {
01272                 throw cms::Exception("DataCorrupt") << "SiStripTemplate::2xtemp called with xhit = " << xhit << std::endl;
01273         }
01274 #else
01275         assert(xpix >= 0.f);
01276 #endif
01277         
01278 // cpix is struck pixel(strip) of the cluster   
01279         
01280         int cpix = (int)xpix;
01281         int shift = BSHX - cpix;
01282         
01283 // xbin the floating bin number and cbin is the bin number (between 0 and 7) of the interpolated template
01284         
01285         float xbin = 8.*(xpix-(float)cpix);
01286         int cbin = (int)xbin;
01287         
01288         float xfrac = xbin-(float)cbin;
01289         
01290         int sizex = std::min((int)cluster.size(), BSXSIZE);
01291         
01292 // shift and interpolate the correct cluster shape      
01293         
01294         for(i=0; i<sizex; ++i) {
01295                 j = i+shift;
01296                 if(j < 0 || j > sizex-1) {cluster[i] = 0.f;} else {
01297               cluster[i]=(1.f-xfrac)*xtemp_[cbin][j]+xfrac*xtemp_[cbin+1][j];
01298                         if(cluster[i] < s50_) cluster[i] = 0.f;
01299                 }
01300                 
01301 // Return cluster in same charge units          
01302                 
01303                 cluster[i] /= qscale_;
01304         }
01305 
01306         
01307         return;
01308         
01309 } // End sxtemp
01310 
01311 // ************************************************************************************************************ 
01313 // ************************************************************************************************************ 
01314 int SiStripTemplate::cxtemp()
01315 
01316 {
01317         // Retrieve already interpolated quantities
01318         
01319         // Local variables 
01320         int j;
01321         
01322         // Analyze only pixels along the central entry
01323         // First, find the maximum signal and then work out to the edges
01324         
01325         float sigmax = 0.f;
01326         int jmax = -1;
01327         
01328         for(j=0; j<BSXSIZE; ++j) {
01329                 if(xtemp_[4][j] > sigmax) {
01330                         sigmax = xtemp_[4][j];
01331                         jmax = j;
01332                 }       
01333         }
01334         if(sigmax < 2.*s50_ || jmax<1 || jmax>BSXM2) {return -1;}
01335         
01336         //  Now search forward and backward
01337         
01338         int jend = jmax;
01339         
01340         for(j=jmax+1; j<BSXM1; ++j) {
01341                 if(xtemp_[4][j] < 2.*s50_) break;
01342            jend = j;
01343    }
01344 
01345    int jbeg = jmax;
01346 
01347    for(j=jmax-1; j>0; --j) {
01348            if(xtemp_[4][j] < 2.*s50_) break;
01349            jbeg = j;
01350    }
01351 
01352    return (jbeg+jend)/2;
01353 
01354 } // End cxtemp
01355 
01356 
01357 // ************************************************************************************************************ 
01361 // ************************************************************************************************************ 
01362 void SiStripTemplate::xtemp3d_int(int nxpix, int& nxbins)
01363 
01364 {       
01365     // Retrieve already interpolated quantities
01366     
01367     // Local variables 
01368     int i, j, k;
01369         int ioff0, ioffp, ioffm;
01370     
01371     // Verify that input parameters are in valid range
01372     
01373 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01374         if(nxpix < 1 || nxpix >= BSXM3) {
01375         throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp3d called with nxpix = " << nxpix << std::endl;
01376         }
01377 #else
01378         assert(nxpix > 0 && nxpix < BSXM3);
01379 #endif
01380         
01381     // Calculate the size of the shift in pixels needed to span the entire cluster
01382     
01383     float diff = fabsf(nxpix - clslenx_)/2. + 1.f;
01384         int nshift = (int)diff;
01385         if((diff - nshift) > 0.5f) {++nshift;}
01386     
01387     // Calculate the number of bins needed to specify each hit range
01388     
01389     nxbins_ = 9 + 16*nshift;
01390         
01391     // Create a 2-d working template with the correct size
01392     
01393         temp2dx_.resize(boost::extents[nxbins_][BSXSIZE]);
01394         
01395     //  The 9 central bins are copied from the interpolated private store
01396     
01397         ioff0 = 8*nshift;
01398         
01399         for(i=0; i<9; ++i) {
01400         for(j=0; j<BSXSIZE; ++j) {
01401             temp2dx_[i+ioff0][j]=xtemp_[i][j];
01402         }
01403         }
01404         
01405     // Add the +- shifted templates     
01406     
01407         for(k=1; k<=nshift; ++k) {
01408         ioffm=ioff0-k*8;
01409         for(i=0; i<8; ++i) {
01410             for(j=0; j<k; ++j) {
01411                 temp2dx_[i+ioffm][BSXM1-j] = 0.f;
01412             }
01413             for(j=0; j<BSXSIZE-k; ++j) {
01414                 temp2dx_[i+ioffm][j]=xtemp_[i][j+k];
01415             }
01416         }
01417         ioffp=ioff0+k*8;
01418         for(i=1; i<9; ++i) {
01419             for(j=0; j<k; ++j) {
01420                 temp2dx_[i+ioffp][j] = 0.f;
01421             }
01422             for(j=0; j<BSXSIZE-k; ++j) {
01423                 temp2dx_[i+ioffp][j+k]=xtemp_[i][j];
01424             }
01425         }
01426         }
01427     
01428     nxbins = nxbins_;                                   
01429         
01430         return;
01431         
01432 } // End xtemp3d_int
01433 
01434 
01435 
01436 // ************************************************************************************************************ 
01440 // ************************************************************************************************************ 
01441 void SiStripTemplate::xtemp3d(int i, int j, std::vector<float>& xtemplate)
01442 
01443 {       
01444     // Sum two 2-d templates to make the 3-d template
01445         if(i >= 0 && i < nxbins_ && j <= i) {
01446         for(int k=0; k<BSXSIZE; ++k) {
01447             xtemplate[k]=temp2dx_[i][k]+temp2dx_[j][k]; 
01448         }
01449     } else {
01450         for(int k=0; k<BSXSIZE; ++k) {
01451             xtemplate[k]=0.;    
01452         } 
01453     }
01454         
01455         return;
01456         
01457 } // End xtemp3d
01458 
01459 
01460 
01461 // ************************************************************************************************************ 
01468 // ************************************************************************************************************ 
01469 int SiStripTemplate::qbin(int id, float cotalpha, float cotbeta, float qclus)
01470                  
01471 {
01472     // Interpolate for a new set of track angles 
01473     
01474     // Local variables 
01475     int i, binq;
01476         int ilow, ihigh, Ny, Nxx, Nyx, index;
01477         float yratio;
01478         float acotb, qscale, qavg, qmin, qmin2, fq, qtotal, qcorrect, cotalpha0;
01479         
01480 
01481 // Find the index corresponding to id
01482 
01483        index = -1;
01484        for(i=0; i<(int)theStripTemp_.size(); ++i) {
01485         
01486               if(id == theStripTemp_[i].head.ID) {
01487            
01488                  index = i;
01489                      break;
01490           }
01491             }
01492          
01493 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01494            if(index < 0 || index >= (int)theStripTemp_.size()) {
01495               throw cms::Exception("DataCorrupt") << "SiStripTemplate::qbin can't find needed template ID = " << id << std::endl;
01496            }
01497 #else
01498            assert(index >= 0 && index < (int)theStripTemp_.size());
01499 #endif
01500          
01501 //              
01502 
01503 // Interpolate the absolute value of cot(beta)     
01504     
01505     acotb = fabs((double)cotbeta);
01506         
01507 //      qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)       
01508 
01509         //      qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)       
01510         
01511         cotalpha0 =  theStripTemp_[index].enty[0].cotalpha;
01512         qcorrect=std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta*cotbeta+cotalpha0*cotalpha0));
01513                                         
01514         // Copy the charge scaling factor to the private variable     
01515                 
01516            qscale = theStripTemp_[index].head.qscale;
01517                 
01518        Ny = theStripTemp_[index].head.NTy;
01519        Nyx = theStripTemp_[index].head.NTyx;
01520        Nxx = theStripTemp_[index].head.NTxx;
01521                 
01522 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01523                 if(Ny < 2 || Nyx < 1 || Nxx < 2) {
01524                         throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
01525                 }
01526 #else
01527                 assert(Ny > 1 && Nyx > 0 && Nxx > 1);
01528 #endif
01529         
01530 // next, loop over all y-angle entries   
01531 
01532            ilow = 0;
01533            yratio = 0.f;
01534 
01535            if(acotb >= theStripTemp_[index].enty[Ny-1].cotbeta) {
01536         
01537                ilow = Ny-2;
01538                    yratio = 1.f;
01539                 
01540            } else {
01541            
01542               if(acotb >= theStripTemp_[index].enty[0].cotbeta) {
01543 
01544              for (i=0; i<Ny-1; ++i) { 
01545     
01546                 if( theStripTemp_[index].enty[i].cotbeta <= acotb && acotb < theStripTemp_[index].enty[i+1].cotbeta) {
01547                   
01548                        ilow = i;
01549                            yratio = (acotb - theStripTemp_[index].enty[i].cotbeta)/(theStripTemp_[index].enty[i+1].cotbeta - theStripTemp_[index].enty[i].cotbeta);
01550                            break;                        
01551                         }
01552                  }
01553                   } 
01554            }
01555         
01556            ihigh=ilow + 1;
01557                           
01558 // Interpolate/store all y-related quantities (flip displacements when flip_y)
01559 
01560            qavg = (1.f - yratio)*theStripTemp_[index].enty[ilow].qavg + yratio*theStripTemp_[index].enty[ihigh].qavg;
01561            qavg *= qcorrect;
01562            qmin = (1.f - yratio)*theStripTemp_[index].enty[ilow].qmin + yratio*theStripTemp_[index].enty[ihigh].qmin;
01563            qmin *= qcorrect;
01564            qmin2 = (1.f - yratio)*theStripTemp_[index].enty[ilow].qmin2 + yratio*theStripTemp_[index].enty[ihigh].qmin2;
01565            qmin2 *= qcorrect;
01566            
01567         
01568 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01569         if(qavg <= 0.f || qmin <= 0.f) {
01570                 throw cms::Exception("DataCorrupt") << "SiStripTemplate::qbin, qavg or qmin <= 0," 
01571                 << " Probably someone called the generic pixel reconstruction with an illegal trajectory state" << std::endl;
01572         }
01573 #else
01574         assert(qavg > 0.f && qmin > 0.f);
01575 #endif
01576         
01577 //  Scale the input charge to account for differences between pixelav and CMSSW simulation or data      
01578         
01579         qtotal = qscale*qclus;
01580         
01581 // uncertainty and final corrections depend upon total charge bin          
01582            
01583         fq = qtotal/qavg;
01584         if(fq > 1.5f) {
01585            binq=0;
01586         } else {
01587            if(fq > 1.0f) {
01588               binq=1;
01589            } else {
01590                   if(fq > 0.85f) {
01591                          binq=2;
01592                   } else {
01593                          binq=3;
01594                   }
01595            }
01596         }
01597         
01598 // If the charge is too small (then flag it)
01599         
01600         if(qtotal < 0.95f*qmin) {binq = 5;} else {if(qtotal < 0.95f*qmin2) {binq = 4;}}
01601                 
01602     return binq;
01603   
01604 } // qbin
01605 
01606 
01607 // ************************************************************************************************************ 
01612 // ************************************************************************************************************ 
01613 void SiStripTemplate::vavilov_pars(double& mpv, double& sigma, double& kappa)
01614 
01615 {
01616         // Local variables 
01617         int i;
01618         int ilow, ihigh, Ny;
01619         float yratio, cotb, cotalpha0, arg;
01620         
01621 // Interpolate in cotbeta only for the correct total path length (converts cotalpha, cotbeta into an effective cotbeta) 
01622         
01623         cotalpha0 =  theStripTemp_[index_id_].enty[0].cotalpha;
01624     arg = cotb_current_*cotb_current_ + cota_current_*cota_current_ - cotalpha0*cotalpha0;
01625     if(arg < 0.f) arg = 0.f;
01626         cotb = std::sqrt(arg);
01627                 
01628 // Copy the charge scaling factor to the private variable     
01629         
01630         Ny = theStripTemp_[index_id_].head.NTy;
01631         
01632 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01633         if(Ny < 2) {
01634                 throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny = " << Ny << std::endl;
01635         }
01636 #else
01637         assert(Ny > 1);
01638 #endif
01639         
01640 // next, loop over all y-angle entries   
01641         
01642         ilow = 0;
01643         yratio = 0.f;
01644         
01645         if(cotb >= theStripTemp_[index_id_].enty[Ny-1].cotbeta) {
01646                 
01647                 ilow = Ny-2;
01648                 yratio = 1.f;
01649                 
01650         } else {
01651                 
01652                 if(cotb >= theStripTemp_[index_id_].enty[0].cotbeta) {
01653                         
01654                         for (i=0; i<Ny-1; ++i) { 
01655                                 
01656                                 if( theStripTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < theStripTemp_[index_id_].enty[i+1].cotbeta) {
01657                                         
01658                                         ilow = i;
01659                                         yratio = (cotb - theStripTemp_[index_id_].enty[i].cotbeta)/(theStripTemp_[index_id_].enty[i+1].cotbeta - theStripTemp_[index_id_].enty[i].cotbeta);
01660                                         break;                   
01661                                 }
01662                         }
01663                 } 
01664         }
01665         
01666         ihigh=ilow + 1;
01667         
01668 // Interpolate Vavilov parameters
01669         
01670         mpvvav_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].mpvvav + yratio*theStripTemp_[index_id_].enty[ihigh].mpvvav;
01671         sigmavav_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].sigmavav + yratio*theStripTemp_[index_id_].enty[ihigh].sigmavav;
01672         kappavav_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].kappavav + yratio*theStripTemp_[index_id_].enty[ihigh].kappavav;
01673         
01674 // Copy to parameter list
01675         
01676         
01677         mpv = (double)mpvvav_;
01678         sigma = (double)sigmavav_;
01679         kappa = (double)kappavav_;
01680         
01681         return;
01682         
01683 } // vavilov_pars
01684 
01685 // ************************************************************************************************************ 
01690 // ************************************************************************************************************ 
01691 void SiStripTemplate::vavilov2_pars(double& mpv, double& sigma, double& kappa)
01692 
01693 {
01694         // Local variables 
01695         int i;
01696         int ilow, ihigh, Ny;
01697         float yratio, cotb, cotalpha0, arg;
01698         
01699 // Interpolate in cotbeta only for the correct total path length (converts cotalpha, cotbeta into an effective cotbeta) 
01700         
01701         cotalpha0 =  theStripTemp_[index_id_].enty[0].cotalpha;
01702     arg = cotb_current_*cotb_current_ + cota_current_*cota_current_ - cotalpha0*cotalpha0;
01703     if(arg < 0.f) arg = 0.f;
01704         cotb = std::sqrt(arg);
01705         
01706         // Copy the charge scaling factor to the private variable     
01707         
01708         Ny = theStripTemp_[index_id_].head.NTy;
01709         
01710 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
01711         if(Ny < 2) {
01712                 throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny = " << Ny << std::endl;
01713         }
01714 #else
01715         assert(Ny > 1);
01716 #endif
01717         
01718         // next, loop over all y-angle entries   
01719         
01720         ilow = 0;
01721         yratio = 0.f;
01722         
01723         if(cotb >= theStripTemp_[index_id_].enty[Ny-1].cotbeta) {
01724                 
01725                 ilow = Ny-2;
01726                 yratio = 1.f;
01727                 
01728         } else {
01729                 
01730                 if(cotb >= theStripTemp_[index_id_].enty[0].cotbeta) {
01731                         
01732                         for (i=0; i<Ny-1; ++i) { 
01733                                 
01734                                 if( theStripTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < theStripTemp_[index_id_].enty[i+1].cotbeta) {
01735                                         
01736                                         ilow = i;
01737                                         yratio = (cotb - theStripTemp_[index_id_].enty[i].cotbeta)/(theStripTemp_[index_id_].enty[i+1].cotbeta - theStripTemp_[index_id_].enty[i].cotbeta);
01738                                         break;                   
01739                                 }
01740                         }
01741                 } 
01742         }
01743         
01744         ihigh=ilow + 1;
01745         
01746         // Interpolate Vavilov parameters
01747         
01748         mpvvav2_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].mpvvav2 + yratio*theStripTemp_[index_id_].enty[ihigh].mpvvav2;
01749         sigmavav2_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].sigmavav2 + yratio*theStripTemp_[index_id_].enty[ihigh].sigmavav2;
01750         kappavav2_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].kappavav2 + yratio*theStripTemp_[index_id_].enty[ihigh].kappavav2;
01751         
01752         // Copy to parameter list
01753         
01754         mpv = (double)mpvvav2_;
01755         sigma = (double)sigmavav2_;
01756         kappa = (double)kappavav2_;
01757         
01758         return;
01759         
01760 } // vavilov2_pars
01761 
01762 
01763
/data/refman/pasoursint/CMSSW_5_3_10_patch1/src/RecoLocalTracker/SiStripRecHitConverter/src/SiStripTemplate.cc
