CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplate2D.cc

Go to the documentation of this file.
00001 //
00002 //  SiPixelTemplate2D.cc  Version 1.03 
00003 //
00004 //  Full 2-D templates for cluster splitting, simulated cluster reweighting, and improved cluster probability
00005 //
00006 // Created by Morris Swartz on 12/01/09.
00007 // Copyright 2009 __TheJohnsHopkinsUniversity__. All rights reserved.
00008 //
00009 // V1.01 - fix qavg_ filling
00010 // V1.02 - Add locBz to test if FPix use is out of range
00011 // V1.03 - Fix edge checking on final template to increase template size and to properly truncate cluster
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 
00029 
00030 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00031 #include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplate2D.h"
00032 #include "FWCore/ParameterSet/interface/FileInPath.h"
00033 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00034 #define LOGERROR(x) LogError(x)
00035 #define LOGINFO(x) LogInfo(x)
00036 #define ENDL " "
00037 #include "FWCore/Utilities/interface/Exception.h"
00038 using namespace edm;
00039 #else
00040 #include "SiPixelTemplate2D.h"
00041 #define LOGERROR(x) std::cout << x << ": "
00042 #define LOGINFO(x) std::cout << x << ": "
00043 #define ENDL std::endl
00044 #endif
00045 
00046 //**************************************************************** 
00051 //**************************************************************** 
00052 bool SiPixelTemplate2D::pushfile(int filenum)
00053 {
00054     // Add template stored in external file numbered filenum to theTemplateStore
00055     
00056     // Local variables 
00057     int i, j, k, l, iy, jx;
00058         const char *tempfile;
00059         //      char title[80]; remove this
00060     char c;
00061         const int code_version={16};
00062         
00063         
00064         
00065         //  Create a filename for this run 
00066         
00067         std::ostringstream tout;
00068         
00069         //  Create different path in CMSSW than standalone
00070         
00071 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00072         tout << "CalibTracker/SiPixelESProducers/data/template_summary2D_zp" 
00073         << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
00074         std::string tempf = tout.str();
00075         edm::FileInPath file( tempf.c_str() );
00076         tempfile = (file.fullPath()).c_str();
00077 #else
00078         tout << "template_summary2D_zp" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
00079         std::string tempf = tout.str();
00080         tempfile = tempf.c_str();
00081 #endif
00082         
00083         //  open the template file 
00084         
00085         std::ifstream in_file(tempfile, std::ios::in);
00086         
00087         if(in_file.is_open()) {
00088                 
00089                 // Create a local template storage entry
00090                 
00091                 SiPixelTemplateStore2D theCurrentTemp;
00092                 
00093                 // Read-in a header string first and print it    
00094                 
00095                 for (i=0; (c=in_file.get()) != '\n'; ++i) {
00096                         if(i < 79) {theCurrentTemp.head.title[i] = c;}
00097                 }
00098                 if(i > 78) {i=78;}
00099                 theCurrentTemp.head.title[i+1] ='\0';
00100                 LOGINFO("SiPixelTemplate2D") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
00101                 
00102                 // next, the header information     
00103                 
00104                 in_file >> theCurrentTemp.head.ID  >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx
00105                 >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
00106                 >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >> theCurrentTemp.head.zsize;
00107                 
00108                 if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file, no template load" << ENDL; return false;}
00109                 
00110                 LOGINFO("SiPixelTemplate2D") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version " << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield 
00111                 << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
00112                 << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
00113                 << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
00114                 << ", 1/2 threshold " << theCurrentTemp.head.s50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", x Lorentz width " << theCurrentTemp.head.lorxwidth    
00115                 << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
00116                 
00117                 if(theCurrentTemp.head.templ_version != code_version) {LOGERROR("SiPixelTemplate2D") << "code expects version " << code_version << ", no template load" << ENDL; return false;}
00118                 
00119                 if(theCurrentTemp.head.NTy != 0) {LOGERROR("SiPixelTemplate2D") << "Trying to load 1-d template info into the 2-d template object, check your DB/global tag!" << ENDL; return false;}
00120                 
00121                 // next, layout the 2-d structure needed to store template
00122                 
00123                 theCurrentTemp.entry.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
00124                 
00125                 // Read in the file info
00126                 
00127                 for (iy=0; iy < theCurrentTemp.head.NTyx; ++iy) {    
00128                         for(jx=0; jx < theCurrentTemp.head.NTxx; ++jx) {
00129                         
00130                            in_file >> theCurrentTemp.entry[iy][jx].runnum >> theCurrentTemp.entry[iy][jx].costrk[0] 
00131                            >> theCurrentTemp.entry[iy][jx].costrk[1] >> theCurrentTemp.entry[iy][jx].costrk[2]; 
00132                         
00133                            if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 1, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
00134                         
00135                         // Calculate cot(alpha) and cot(beta) for this entry 
00136                         
00137                            theCurrentTemp.entry[iy][jx].cotalpha = theCurrentTemp.entry[iy][jx].costrk[0]/theCurrentTemp.entry[iy][jx].costrk[2];
00138                         
00139                            theCurrentTemp.entry[iy][jx].cotbeta = theCurrentTemp.entry[iy][jx].costrk[1]/theCurrentTemp.entry[iy][jx].costrk[2];
00140                         
00141                            in_file >> theCurrentTemp.entry[iy][jx].qavg >> theCurrentTemp.entry[iy][jx].pixmax >> theCurrentTemp.entry[iy][jx].sxymax >> theCurrentTemp.entry[iy][jx].iymin
00142                            >> theCurrentTemp.entry[iy][jx].iymax >> theCurrentTemp.entry[iy][jx].jxmin >> theCurrentTemp.entry[iy][jx].jxmax;
00143                         
00144                            if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 2, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
00145                         
00146                            for (k=0; k<2; ++k) {
00147                                 
00148                                   in_file >> theCurrentTemp.entry[iy][jx].xypar[k][0] >> theCurrentTemp.entry[iy][jx].xypar[k][1] 
00149                                   >> theCurrentTemp.entry[iy][jx].xypar[k][2] >> theCurrentTemp.entry[iy][jx].xypar[k][3] >> theCurrentTemp.entry[iy][jx].xypar[k][4];
00150                                 
00151                                   if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 3, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
00152                                 
00153                            }
00154                                 
00155                            for (k=0; k<2; ++k) {
00156                                         
00157                                   in_file >> theCurrentTemp.entry[iy][jx].lanpar[k][0] >> theCurrentTemp.entry[iy][jx].lanpar[k][1] 
00158                                   >> theCurrentTemp.entry[iy][jx].lanpar[k][2] >> theCurrentTemp.entry[iy][jx].lanpar[k][3] >> theCurrentTemp.entry[iy][jx].lanpar[k][4];
00159                                         
00160                                   if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 4, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
00161                                         
00162                            }
00163                                 
00164                            for (l=0; l<7; ++l) {        
00165                               for (k=0; k<7; ++k) {
00166                                                 for (j=0; j<T2XSIZE; ++j) {
00167                                 
00168                                         for (i=0; i<T2YSIZE; ++i) {in_file >> theCurrentTemp.entry[iy][jx].xytemp[k][l][i][j];}
00169                                 
00170                                         if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 5, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
00171                                       }
00172                               }
00173                            }
00174                                                 
00175                                 
00176                            in_file >> theCurrentTemp.entry[iy][jx].chi2minone >> theCurrentTemp.entry[iy][jx].chi2avgone >> theCurrentTemp.entry[iy][jx].chi2min[0] >> theCurrentTemp.entry[iy][jx].chi2avg[0]
00177                                    >> theCurrentTemp.entry[iy][jx].chi2min[1] >> theCurrentTemp.entry[iy][jx].chi2avg[1]>> theCurrentTemp.entry[iy][jx].chi2min[2] >> theCurrentTemp.entry[iy][jx].chi2avg[2]
00178                                    >> theCurrentTemp.entry[iy][jx].chi2min[3] >> theCurrentTemp.entry[iy][jx].chi2avg[3];
00179                                 
00180                            if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 6, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
00181                         
00182                            in_file >> theCurrentTemp.entry[iy][jx].spare[0] >> theCurrentTemp.entry[iy][jx].spare[1] >> theCurrentTemp.entry[iy][jx].spare[2] >> theCurrentTemp.entry[iy][jx].spare[3]
00183                                     >> theCurrentTemp.entry[iy][jx].spare[4] >> theCurrentTemp.entry[iy][jx].spare[5] >> theCurrentTemp.entry[iy][jx].spare[6] >> theCurrentTemp.entry[iy][jx].spare[7]
00184                                     >> theCurrentTemp.entry[iy][jx].spare[8]  >> theCurrentTemp.entry[iy][jx].spare[9];
00185                                 
00186                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 7, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
00187                                         
00188                            in_file >> theCurrentTemp.entry[iy][jx].spare[10] >> theCurrentTemp.entry[iy][jx].spare[11] >> theCurrentTemp.entry[iy][jx].spare[12] >> theCurrentTemp.entry[iy][jx].spare[13]
00189                                 >> theCurrentTemp.entry[iy][jx].spare[14] >> theCurrentTemp.entry[iy][jx].spare[15] >> theCurrentTemp.entry[iy][jx].spare[16] >> theCurrentTemp.entry[iy][jx].spare[17]
00190                                 >> theCurrentTemp.entry[iy][jx].spare[18]  >> theCurrentTemp.entry[iy][jx].spare[19];
00191                                 
00192                                 if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 8, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
00193                    }
00194                         
00195                 }
00196                                 
00197                 
00198                 in_file.close();
00199                 
00200                 // Add this template to the store
00201                 
00202                 thePixelTemp_.push_back(theCurrentTemp);
00203                 
00204                 return true;
00205                 
00206         } else {
00207                 
00208                 // If file didn't open, report this
00209                 
00210                 LOGERROR("SiPixelTemplate2D") << "Error opening File" << tempfile << ENDL;
00211                 return false;
00212                 
00213    }
00214         
00215 } // TempInit 
00216 
00217 
00218 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00219 
00220 //**************************************************************** 
00224 //**************************************************************** 
00225 bool SiPixelTemplate2D::pushfile(const SiPixelTemplateDBObject& dbobject)
00226 {
00227         // Add template stored in external dbobject to theTemplateStore
00228     
00229         // Local variables 
00230         int i, j, k, l, iy, jx;
00231         //      const char *tempfile;
00232         const int code_version={16};
00233         
00234         // We must create a new object because dbobject must be a const and our stream must not be
00235         SiPixelTemplateDBObject db = dbobject;
00236         
00237         // Create a local template storage entry
00238         SiPixelTemplateStore2D theCurrentTemp;
00239         
00240         // Fill the template storage for each template calibration stored in the db
00241         for(int m=0; m<db.numOfTempl(); ++m)
00242         {
00243                 
00244                 // Read-in a header string first and print it    
00245                 
00246                 SiPixelTemplateDBObject::char2float temp;
00247                 for (i=0; i<20; ++i) {
00248                         temp.f = db.sVector()[db.index()];
00249                         theCurrentTemp.head.title[4*i] = temp.c[0];
00250                         theCurrentTemp.head.title[4*i+1] = temp.c[1];
00251                         theCurrentTemp.head.title[4*i+2] = temp.c[2];
00252                         theCurrentTemp.head.title[4*i+3] = temp.c[3];
00253                         db.incrementIndex(1);
00254                 }
00255                 theCurrentTemp.head.title[79] = '\0';
00256                 LOGINFO("SiPixelTemplate2D") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
00257                 
00258                 // next, the header information     
00259                 
00260                 db >> theCurrentTemp.head.ID  >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx
00261                 >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
00262                 >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >> theCurrentTemp.head.zsize;
00263                 
00264                 if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file, no template load" << ENDL; return false;}
00265                 
00266                 LOGINFO("SiPixelTemplate2D") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version " << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield 
00267                 << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
00268                 << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
00269                 << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
00270                 << ", 1/2 threshold " << theCurrentTemp.head.s50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", x Lorentz width " << theCurrentTemp.head.lorxwidth    
00271                 << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
00272                 
00273                 if(theCurrentTemp.head.templ_version != code_version) {LOGERROR("SiPixelTemplate2D") << "code expects version " << code_version << ", no template load" << ENDL; return false;}
00274                 
00275                 if(theCurrentTemp.head.NTy != 0) {LOGERROR("SiPixelTemplate2D") << "Trying to load 1-d template info into the 2-d template object, check your DB/global tag!" << ENDL; return false;}
00276                 
00277                 // next, layout the 2-d structure needed to store template
00278                 
00279                 theCurrentTemp.entry.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
00280                 
00281                 // Read in the file info
00282                 
00283                 for (iy=0; iy < theCurrentTemp.head.NTyx; ++iy) {    
00284                         for(jx=0; jx < theCurrentTemp.head.NTxx; ++jx) {
00285                                 
00286                                 db >> theCurrentTemp.entry[iy][jx].runnum >> theCurrentTemp.entry[iy][jx].costrk[0] 
00287                                 >> theCurrentTemp.entry[iy][jx].costrk[1] >> theCurrentTemp.entry[iy][jx].costrk[2]; 
00288                                 
00289                                 if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 1, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
00290                                 
00291                                 // Calculate cot(alpha) and cot(beta) for this entry 
00292                                 
00293                                 theCurrentTemp.entry[iy][jx].cotalpha = theCurrentTemp.entry[iy][jx].costrk[0]/theCurrentTemp.entry[iy][jx].costrk[2];
00294                                 
00295                                 theCurrentTemp.entry[iy][jx].cotbeta = theCurrentTemp.entry[iy][jx].costrk[1]/theCurrentTemp.entry[iy][jx].costrk[2];
00296                                 
00297                                 db >> theCurrentTemp.entry[iy][jx].qavg >> theCurrentTemp.entry[iy][jx].pixmax >> theCurrentTemp.entry[iy][jx].sxymax >> theCurrentTemp.entry[iy][jx].iymin
00298                                 >> theCurrentTemp.entry[iy][jx].iymax >> theCurrentTemp.entry[iy][jx].jxmin >> theCurrentTemp.entry[iy][jx].jxmax;
00299                                 
00300                                 if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 2, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
00301                                 
00302                                 for (k=0; k<2; ++k) {
00303                                         
00304                                         db >> theCurrentTemp.entry[iy][jx].xypar[k][0] >> theCurrentTemp.entry[iy][jx].xypar[k][1] 
00305                                         >> theCurrentTemp.entry[iy][jx].xypar[k][2] >> theCurrentTemp.entry[iy][jx].xypar[k][3] >> theCurrentTemp.entry[iy][jx].xypar[k][4];
00306                                         
00307                                         if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 3, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
00308                                         
00309                                 }
00310                                 
00311                                 for (k=0; k<2; ++k) {
00312                                         
00313                                         db >> theCurrentTemp.entry[iy][jx].lanpar[k][0] >> theCurrentTemp.entry[iy][jx].lanpar[k][1] 
00314                                         >> theCurrentTemp.entry[iy][jx].lanpar[k][2] >> theCurrentTemp.entry[iy][jx].lanpar[k][3] >> theCurrentTemp.entry[iy][jx].lanpar[k][4];
00315                                         
00316                                         if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 4, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
00317                                         
00318                                 }
00319                                 
00320                                 for (l=0; l<7; ++l) {   
00321                                         for (k=0; k<7; ++k) {
00322                                                 for (j=0; j<T2XSIZE; ++j) {
00323                                                         
00324                                                         for (i=0; i<T2YSIZE; ++i) {db >> theCurrentTemp.entry[iy][jx].xytemp[k][l][i][j];}
00325                                                         
00326                                                         if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 5, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
00327                                                 }
00328                                         }
00329                                 }
00330                                 
00331                                 
00332                                 db >> theCurrentTemp.entry[iy][jx].chi2minone >> theCurrentTemp.entry[iy][jx].chi2avgone >> theCurrentTemp.entry[iy][jx].chi2min[0] >> theCurrentTemp.entry[iy][jx].chi2avg[0]
00333                                 >> theCurrentTemp.entry[iy][jx].chi2min[1] >> theCurrentTemp.entry[iy][jx].chi2avg[1]>> theCurrentTemp.entry[iy][jx].chi2min[2] >> theCurrentTemp.entry[iy][jx].chi2avg[2]
00334                                 >> theCurrentTemp.entry[iy][jx].chi2min[3] >> theCurrentTemp.entry[iy][jx].chi2avg[3];
00335                                 
00336                                 if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 6, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
00337                                 
00338                            db >> theCurrentTemp.entry[iy][jx].spare[0] >> theCurrentTemp.entry[iy][jx].spare[1] >> theCurrentTemp.entry[iy][jx].spare[2] >> theCurrentTemp.entry[iy][jx].spare[3]
00339                                 >> theCurrentTemp.entry[iy][jx].spare[4] >> theCurrentTemp.entry[iy][jx].spare[5] >> theCurrentTemp.entry[iy][jx].spare[6] >> theCurrentTemp.entry[iy][jx].spare[7]
00340                                 >> theCurrentTemp.entry[iy][jx].spare[8]  >> theCurrentTemp.entry[iy][jx].spare[9];
00341                                 
00342                                 if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 7, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
00343                                 
00344                            db >> theCurrentTemp.entry[iy][jx].spare[10] >> theCurrentTemp.entry[iy][jx].spare[11] >> theCurrentTemp.entry[iy][jx].spare[12] >> theCurrentTemp.entry[iy][jx].spare[13]
00345                                 >> theCurrentTemp.entry[iy][jx].spare[14] >> theCurrentTemp.entry[iy][jx].spare[15] >> theCurrentTemp.entry[iy][jx].spare[16] >> theCurrentTemp.entry[iy][jx].spare[17]
00346                                 >> theCurrentTemp.entry[iy][jx].spare[18]  >> theCurrentTemp.entry[iy][jx].spare[19];
00347                                 
00348                                 if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 8, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
00349                                         
00350                                 }
00351                         }
00352                         
00353                 }
00354                 
00355                                 
00356 // Add this template to the store
00357                 
00358                 thePixelTemp_.push_back(theCurrentTemp);
00359                 
00360         return true;
00361         
00362 } // TempInit 
00363 
00364 #endif
00365 
00366 
00367 
00368 // *************************************************************************************************************************************
00380 // *************************************************************************************************************************************
00381 
00382 bool SiPixelTemplate2D::xytemp(int id, float cotalpha, float cotbeta, float locBz, float xhit, float yhit, std::vector<bool>& ydouble, std::vector<bool>& xdouble, float template2d[BXM2][BYM2])
00383 {
00384     // Interpolate for a new set of track angles 
00385     
00386     // Local variables 
00387         int i, j;
00388         int pixx, pixy, k0, k1, l0, l1, imidx, deltax, deltay, iflipy, imin, imax, jmin, jmax;
00389         int m, n;
00390         float acotb, dcota, dcotb, dx, dy, ddx, ddy, adx, ady, tmpxy;
00391         bool flip_y;
00392 //      std::vector <float> xrms(4), xgsig(4), xrmsc2m(4), xgsigc2m(4);
00393         std::vector <float> chi2xavg(4), chi2xmin(4);
00394 
00395 
00396 // Check to see if interpolation is valid     
00397 
00398         if(id != id_current_ || cotalpha != cota_current_ || cotbeta != cotb_current_) {
00399                 
00400                 cota_current_ = cotalpha; cotb_current_ = cotbeta; success_ = true;
00401                 
00402                 if(id != id_current_) {
00403                         
00404                         // Find the index corresponding to id
00405                         
00406                         index_id_ = -1;
00407                         for(i=0; i<(int)thePixelTemp_.size(); ++i) {
00408                                 
00409                                 if(id == thePixelTemp_[i].head.ID) {
00410                                         
00411                                         index_id_ = i;
00412                                         id_current_ = id;
00413                                         
00414 // Copy the charge scaling factor to the private variable     
00415                                         
00416                                         Dtype_ = thePixelTemp_[index_id_].head.Dtype;
00417                                         
00418 // Copy the charge scaling factor to the private variable     
00419                                         
00420                                         qscale_ = thePixelTemp_[index_id_].head.qscale;
00421                                         
00422 // Copy the pseudopixel signal size to the private variable     
00423                                         
00424                                         s50_ = thePixelTemp_[index_id_].head.s50;
00425                                         
00426 // Copy the Lorentz widths to private variables     
00427                                         
00428                                         lorywidth_ = thePixelTemp_[index_id_].head.lorywidth;
00429                                         lorxwidth_ = thePixelTemp_[index_id_].head.lorxwidth;
00430                                         
00431 // Copy the pixel sizes private variables    
00432                                         
00433                                         xsize_ = thePixelTemp_[index_id_].head.xsize;
00434                                         ysize_ = thePixelTemp_[index_id_].head.ysize;
00435                                         zsize_ = thePixelTemp_[index_id_].head.zsize;
00436                                         
00437 // Determine the size of this template    
00438                                         
00439                                         Nyx_ = thePixelTemp_[index_id_].head.NTyx;
00440                                         Nxx_ = thePixelTemp_[index_id_].head.NTxx;
00441 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00442                                         if(Nyx_ < 2 || Nxx_ < 2) {
00443                                                 throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Nyx/Nxx = " << Nyx_ << "/" << Nxx_ << std::endl;
00444                                         }
00445 #else
00446                                         assert(Nyx_ > 1 && Nxx_ > 1);
00447 #endif
00448                                         imidx = Nxx_/2;
00449                                         
00450                                         cotalpha0_ =  thePixelTemp_[index_id_].entry[0][0].cotalpha;
00451                                         cotalpha1_ =  thePixelTemp_[index_id_].entry[0][Nxx_-1].cotalpha;
00452                                         deltacota_ = (cotalpha1_-cotalpha0_)/(float)(Nxx_-1);
00453                                         
00454                                         cotbeta0_ =  thePixelTemp_[index_id_].entry[0][imidx].cotbeta;
00455                                         cotbeta1_ =  thePixelTemp_[index_id_].entry[Nyx_-1][imidx].cotbeta;
00456                                         deltacotb_ = (cotbeta1_-cotbeta0_)/(float)(Nyx_-1);
00457                                         
00458                                         break;
00459                                 }
00460                         }
00461                 }
00462         }
00463         
00464 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00465         if(index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
00466                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::interpolate can't find needed template ID = " << id 
00467                 << ", Are you using the correct global tag?" << std::endl;
00468         }
00469 #else
00470         assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
00471 #endif
00472         
00473 // Check angle limits and et up interpolation parameters  
00474         
00475         if(cotalpha < cotalpha0_) {
00476                 success_ = false;
00477                 jx0_ = 0;
00478                 jx1_ = 1;
00479                 adcota_ = 0.f;
00480         } else if(cotalpha > cotalpha1_) {
00481             success_ = false;
00482                 jx0_ = Nxx_ - 1;
00483                 jx1_ = jx0_ - 1;
00484                 adcota_ = 0.f;
00485         } else {
00486            jx0_ = (int)((cotalpha-cotalpha0_)/deltacota_+0.5f);
00487            dcota = (cotalpha - (cotalpha0_ + jx0_*deltacota_))/deltacota_;
00488            adcota_ = fabs(dcota);
00489            if(dcota > 0.f) {jx1_ = jx0_ + 1;if(jx1_ > Nxx_-1) jx1_ = jx0_-1;} else {jx1_ = jx0_ - 1; if(jx1_ < 0) jx1_ = jx0_+1;}               
00490         }
00491                 
00492 // Interpolate the absolute value of cot(beta)     
00493     
00494         acotb = std::abs(cotbeta);
00495         
00496         if(acotb < cotbeta0_) {
00497             success_ = false;
00498                 iy0_ = 0;
00499                 iy1_ = 1;
00500                 adcotb_ = 0.f;
00501         } else if(acotb > cotbeta1_) {
00502             success_ = false;
00503                 iy0_ = Nyx_ - 1;
00504                 iy1_ = iy0_ - 1;
00505                 adcotb_ = 0.f;
00506         } else {
00507                 iy0_ = (int)((acotb-cotbeta0_)/deltacotb_+0.5f);
00508                 dcotb = (acotb - (cotbeta0_ + iy0_*deltacotb_))/deltacotb_;
00509                 adcotb_ = fabs(dcotb);
00510                 if(dcotb > 0.f) {iy1_ = iy0_ + 1; if(iy1_ > Nyx_-1) iy1_ = iy0_-1;} else {iy1_ = iy0_ - 1; if(iy1_ < 0) iy1_ = iy0_+1;}         
00511     }   
00512         
00513 // This works only for IP-related tracks 
00514         
00515         flip_y = false;
00516         if(cotbeta < 0.f) {flip_y = true;}
00517         
00518 // If Fpix-related track has wrong cotbeta-field correlation, return false to trigger simple template instead   
00519         
00520         if(Dtype_ == 1) {
00521                 if(cotbeta*locBz > 0.f) success_ = false;
00522         }
00523         
00524 // Interpolate things in cot(alpha)-cot(beta)   
00525         
00526         qavg_ = thePixelTemp_[index_id_].entry[iy0_][jx0_].qavg
00527         +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].qavg - thePixelTemp_[index_id_].entry[iy0_][jx0_].qavg)
00528         +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].qavg - thePixelTemp_[index_id_].entry[iy0_][jx0_].qavg);
00529         
00530         pixmax_ = thePixelTemp_[index_id_].entry[iy0_][jx0_].pixmax
00531         +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].pixmax - thePixelTemp_[index_id_].entry[iy0_][jx0_].pixmax)
00532         +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].pixmax - thePixelTemp_[index_id_].entry[iy0_][jx0_].pixmax);
00533         
00534         sxymax_ = thePixelTemp_[index_id_].entry[iy0_][jx0_].sxymax
00535         +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].sxymax - thePixelTemp_[index_id_].entry[iy0_][jx0_].sxymax)
00536         +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].sxymax - thePixelTemp_[index_id_].entry[iy0_][jx0_].sxymax);
00537         
00538         chi2avgone_ = thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2avgone
00539         +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].chi2avgone - thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2avgone)
00540         +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].chi2avgone - thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2avgone);
00541         
00542         chi2minone_ = thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2minone
00543         +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].chi2minone - thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2minone)
00544         +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].chi2minone - thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2minone);
00545         
00546         for(i=0; i<4 ; ++i) {
00547                 chi2avg_[i] = thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2avg[i]
00548                 +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].chi2avg[i] - thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2avg[i])
00549                 +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].chi2avg[i] - thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2avg[i]);
00550                 
00551                 chi2min_[i] = thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2min[i]
00552                 +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].chi2min[i] - thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2min[i])
00553                 +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].chi2min[i] - thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2min[i]);
00554         }
00555         
00556         for(i=0; i<2 ; ++i) {
00557                 for(j=0; j<5 ; ++j) {
00558                         // Charge loss switches sides when cot(beta) changes sign
00559                         if(flip_y) {
00560                                 xypary0x0_[1-i][j] = thePixelTemp_[index_id_].entry[iy0_][jx0_].xypar[i][j];
00561                                 xypary1x0_[1-i][j] = thePixelTemp_[index_id_].entry[iy1_][jx0_].xypar[i][j];
00562                                 xypary0x1_[1-i][j] = thePixelTemp_[index_id_].entry[iy0_][jx1_].xypar[i][j];
00563                                 lanpar_[1-i][j] = thePixelTemp_[index_id_].entry[iy0_][jx0_].lanpar[i][j]
00564                                 +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].lanpar[i][j] - thePixelTemp_[index_id_].entry[iy0_][jx0_].lanpar[i][j])
00565                                 +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].lanpar[i][j] - thePixelTemp_[index_id_].entry[iy0_][jx0_].lanpar[i][j]);
00566                         } else {
00567                                 xypary0x0_[i][j] = thePixelTemp_[index_id_].entry[iy0_][jx0_].xypar[i][j];
00568                                 xypary1x0_[i][j] = thePixelTemp_[index_id_].entry[iy1_][jx0_].xypar[i][j];
00569                                 xypary0x1_[i][j] = thePixelTemp_[index_id_].entry[iy0_][jx1_].xypar[i][j];
00570                                 lanpar_[i][j] = thePixelTemp_[index_id_].entry[iy0_][jx0_].lanpar[i][j]
00571                                 +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].lanpar[i][j] - thePixelTemp_[index_id_].entry[iy0_][jx0_].lanpar[i][j])
00572                                 +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].lanpar[i][j] - thePixelTemp_[index_id_].entry[iy0_][jx0_].lanpar[i][j]);
00573                         }
00574                 }
00575         }
00576         
00577 // next, determine the indices of the closest point in k (y-displacement), l (x-displacement)
00578 // pixy and pixx are the indices of the struck pixel in the (Ty,Tx) system
00579 // k0,k1 are the k-indices of the closest and next closest point
00580 // l0,l1 are the l-indices of the closest and next closest point
00581 
00582         pixy = (int)floorf(yhit/ysize_);
00583         dy = yhit-(pixy+0.5f)*ysize_;
00584         if(flip_y) {dy = -dy;}
00585         k0 = (int)(dy/ysize_*6.f+3.5f);
00586         if(k0 < 0) k0 = 0;
00587         if(k0 > 6) k0 = 6;
00588         ddy = 6.f*dy/ysize_ - (k0-3);
00589         ady = fabs(ddy);
00590         if(ddy > 0.f) {k1 = k0 + 1; if(k1 > 6) k1 = k0-1;} else {k1 = k0 - 1; if(k1 < 0) k1 = k0+1;}
00591         pixx = (int)floorf(xhit/xsize_);
00592         dx = xhit-(pixx+0.5f)*xsize_;
00593         l0 = (int)(dx/xsize_*6.f+3.5f);
00594         if(l0 < 0) l0 = 0;
00595         if(l0 > 6) l0 = 6;
00596         ddx = 6.f*dx/xsize_ - (l0-3);
00597         adx = fabs(ddx);
00598         if(ddx > 0.f) {l1 = l0 + 1; if(l1 > 6) l1 = l0-1;} else {l1 = l0 - 1; if(l1 < 0) l1 = l0+1;}
00599         
00600 // OK, lets do the template interpolation.  
00601         
00602 // First find the limits of the indices for non-zero pixels
00603         
00604         imin = std::min(thePixelTemp_[index_id_].entry[iy0_][jx0_].iymin,thePixelTemp_[index_id_].entry[iy1_][jx0_].iymin);
00605         imin = std::min(imin,thePixelTemp_[index_id_].entry[iy0_][jx1_].iymin);
00606         
00607         jmin = std::min(thePixelTemp_[index_id_].entry[iy0_][jx0_].jxmin,thePixelTemp_[index_id_].entry[iy1_][jx0_].jxmin);
00608         jmin = std::min(jmin,thePixelTemp_[index_id_].entry[iy0_][jx1_].jxmin);
00609         
00610         imax = std::max(thePixelTemp_[index_id_].entry[iy0_][jx0_].iymax,thePixelTemp_[index_id_].entry[iy1_][jx0_].iymax);
00611         imax = std::max(imax,thePixelTemp_[index_id_].entry[iy0_][jx1_].iymax);
00612         
00613         jmax = std::max(thePixelTemp_[index_id_].entry[iy0_][jx0_].jxmax,thePixelTemp_[index_id_].entry[iy1_][jx0_].jxmax);
00614         jmax = std::max(jmax,thePixelTemp_[index_id_].entry[iy0_][jx1_].jxmax);
00615                 
00616 // Calculate the x and y offsets to make the new template
00617         
00618 // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
00619         
00620         ++pixy; ++pixx;
00621         
00622 // In the template store, the struck pixel is always (THy,THx)
00623         
00624         deltax = pixx - T2HX;
00625         deltay = pixy - T2HY;
00626         
00627 //  First zero the local 2-d template
00628         
00629         for(j=0; j<BXM2; ++j) {for(i=0; i<BYM2; ++i) {xytemp_[j][i] = 0.f;}}
00630         
00631 // Loop over the non-zero part of the template index space and interpolate
00632         
00633         for(j=jmin; j<=jmax; ++j) {
00634            for(i=imin; i<=imax; ++i) {
00635                         m = deltax+j;
00636                         
00637 // If cot(beta) < 0, we must flip the cluster, iflipy is the flipped y-index
00638                         
00639                         if(flip_y) {iflipy=T2YSIZE-1-i; n = deltay+iflipy;} else {n = deltay+i;}
00640                         if(m>=0 && m<=BXM3 && n>=0 && n<=BYM3) {
00641                            tmpxy = thePixelTemp_[index_id_].entry[iy0_][jx0_].xytemp[k0][l0][i][j] 
00642                                    + adx*(thePixelTemp_[index_id_].entry[iy0_][jx0_].xytemp[k0][l1][i][j] - thePixelTemp_[index_id_].entry[iy0_][jx0_].xytemp[k0][l0][i][j])
00643                                    + ady*(thePixelTemp_[index_id_].entry[iy0_][jx0_].xytemp[k1][l0][i][j] - thePixelTemp_[index_id_].entry[iy0_][jx0_].xytemp[k0][l0][i][j])
00644                                    + adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].xytemp[k0][l0][i][j] - thePixelTemp_[index_id_].entry[iy0_][jx0_].xytemp[k0][l0][i][j])
00645                                    + adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].xytemp[k0][l0][i][j] - thePixelTemp_[index_id_].entry[iy0_][jx0_].xytemp[k0][l0][i][j]);
00646                            if(tmpxy > 0.f) {xytemp_[m][n] = tmpxy;} else {xytemp_[m][n] = 0.f;}
00647                         }
00648                 }
00649         }
00650         
00651 //combine rows and columns to simulate double pixels
00652         
00653         for(n=1; n<BYM3; ++n) {
00654                 if(ydouble[n-1]) {
00655 //  Combine the y-columns
00656                         for(m=1; m<BXM3; ++m) {
00657                                 xytemp_[m][n] += xytemp_[m][n+1];
00658                         }
00659 //  Now shift the remaining pixels over by one column
00660                         for(i=n+1; i<BYM3; ++i) {
00661                            for(m=1; m<BXM3; ++m) {
00662                                    xytemp_[m][i] = xytemp_[m][i+1];
00663                            }
00664                         }
00665                 }
00666         }
00667 
00668 //combine rows and columns to simulate double pixels
00669         
00670         for(m=1; m<BXM3; ++m) {
00671                 if(xdouble[m-1]) {
00672 //  Combine the x-rows
00673                         for(n=1; n<BYM3; ++n) {
00674                                 xytemp_[m][n] += xytemp_[m+1][n];
00675                         }
00676 //  Now shift the remaining pixels over by one row
00677                         for(j=m+1; j<BXM3; ++j) {
00678                            for(n=1; n<BYM3; ++n) {
00679                                    xytemp_[j][n] = xytemp_[j+1][n];
00680                            }
00681                         }
00682                 }
00683         }
00684         
00685 //  Finally, loop through and increment the external template
00686 
00687         for(n=1; n<BYM3; ++n) {
00688                 for(m=1; m<BXM3; ++m) {
00689                         if(xytemp_[m][n] > 0.f) {template2d[m][n] += xytemp_[m][n];}
00690                 }
00691         }
00692         
00693   return success_;
00694 } // xytemp
00695 
00696 
00697 
00698 // *************************************************************************************************************************************
00708 // *************************************************************************************************************************************
00709 
00710 bool SiPixelTemplate2D::xytemp(int id, float cotalpha, float cotbeta, float xhit, float yhit, std::vector<bool>& ydouble, std::vector<bool>& xdouble, float template2d[BXM2][BYM2])
00711 {
00712         // Interpolate for a new set of track angles 
00713         
00714         // Local variables 
00715         float locBz = -1;
00716         if(cotbeta < 0.f) {locBz = -locBz;}
00717 
00718         return SiPixelTemplate2D::xytemp(id, cotalpha, cotbeta, locBz, xhit, yhit, ydouble, xdouble, template2d);
00719         
00720 } // xytemp
00721 
00722 
00723 
00724 // ************************************************************************************************************ 
00730 // ************************************************************************************************************ 
00731 void SiPixelTemplate2D::xysigma2(float qpixel, int index, float& xysig2)
00732 
00733 {
00734     // Interpolate using quantities already stored in the private variables
00735     
00736     // Local variables 
00737         float sigi, sigi2, sigi3, sigi4, qscale, err2, err00;
00738         
00739     // Make sure that input is OK
00740     
00741 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00742     if(index < 2 || index >= BYM2) {
00743                 throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::ysigma2 called with index = " << index << std::endl;
00744         }
00745 #else
00746         assert(index > 1 && index < BYM2);
00747 #endif
00748         
00749         // Define the maximum signal to use in the parameterization 
00750         
00751         // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis 
00752         
00753                         if(qpixel < sxymax_) {
00754                                 sigi = qpixel;
00755                                 qscale = 1.f;
00756                         } else {
00757                                 sigi = sxymax_;
00758                                 qscale = qpixel/sxymax_;
00759                         }
00760                         sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
00761                         if(index <= THXP1) {
00762                                 err00 = xypary0x0_[0][0]+xypary0x0_[0][1]*sigi+xypary0x0_[0][2]*sigi2+xypary0x0_[0][3]*sigi3+xypary0x0_[0][4]*sigi4;
00763                                 err2 = err00
00764                                      +adcota_*(xypary0x1_[0][0]+xypary0x1_[0][1]*sigi+xypary0x1_[0][2]*sigi2+xypary0x1_[0][3]*sigi3+xypary0x1_[0][4]*sigi4 - err00)
00765                                      +adcotb_*(xypary1x0_[0][0]+xypary1x0_[0][1]*sigi+xypary1x0_[0][2]*sigi2+xypary1x0_[0][3]*sigi3+xypary1x0_[0][4]*sigi4 - err00);
00766                         } else {
00767                                 err00 = xypary0x0_[1][0]+xypary0x0_[1][1]*sigi+xypary0x0_[1][2]*sigi2+xypary0x0_[1][3]*sigi3+xypary0x0_[1][4]*sigi4;
00768                                 err2 = err00
00769                                      +adcota_*(xypary0x1_[1][0]+xypary0x1_[1][1]*sigi+xypary0x1_[1][2]*sigi2+xypary0x1_[1][3]*sigi3+xypary0x1_[1][4]*sigi4 - err00)
00770                                      +adcotb_*(xypary1x0_[1][0]+xypary1x0_[1][1]*sigi+xypary1x0_[1][2]*sigi2+xypary1x0_[1][3]*sigi3+xypary1x0_[1][4]*sigi4 - err00);
00771                         }
00772                         xysig2 =qscale*err2;
00773                         if(xysig2 <= 0.f) {LOGERROR("SiPixelTemplate2D") << "neg y-error-squared, id = " << id_current_ << ", index = " << index_id_ << 
00774                         ", cot(alpha) = " << cota_current_ << ", cot(beta) = " << cotb_current_ <<  ", sigi = " << sigi << ENDL;}
00775         
00776         return;
00777         
00778 } // End xysigma2
00779 
00780 
00781 // ************************************************************************************************************ 
00783 // ************************************************************************************************************ 
00784 void SiPixelTemplate2D::landau_par(float lanpar[2][5])
00785 
00786 {
00787         // Interpolate using quantities already stored in the private variables
00788         
00789         // Local variables 
00790         int i,j;
00791         for(i=0; i<2; ++i) {
00792                 for(j=0; j<5; ++j) {
00793                         lanpar[i][j] = lanpar_[i][j];
00794                 }
00795         }
00796         return;
00797         
00798 } // End lan_par
00799 
00800 
00801