00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00017
00018 #else
00019 #include <math.h>
00020 #endif
00021 #include <algorithm>
00022 #include <vector>
00023
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
00055
00056
00057 int i, j, k, l, iy, jx;
00058 const char *tempfile;
00059
00060 char c;
00061 const int code_version={16};
00062
00063
00064
00065
00066
00067 std::ostringstream tout;
00068
00069
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
00084
00085 std::ifstream in_file(tempfile, std::ios::in);
00086
00087 if(in_file.is_open()) {
00088
00089
00090
00091 SiPixelTemplateStore2D theCurrentTemp;
00092
00093
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
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
00122
00123 theCurrentTemp.entry.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
00124
00125
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
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
00201
00202 thePixelTemp_.push_back(theCurrentTemp);
00203
00204 return true;
00205
00206 } else {
00207
00208
00209
00210 LOGERROR("SiPixelTemplate2D") << "Error opening File" << tempfile << ENDL;
00211 return false;
00212
00213 }
00214
00215 }
00216
00217
00218 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
00219
00220
00224
00225 bool SiPixelTemplate2D::pushfile(const SiPixelTemplateDBObject& dbobject)
00226 {
00227
00228
00229
00230 int i, j, k, l, iy, jx;
00231
00232 const int code_version={16};
00233
00234
00235 SiPixelTemplateDBObject db = dbobject;
00236
00237
00238 SiPixelTemplateStore2D theCurrentTemp;
00239
00240
00241 for(int m=0; m<db.numOfTempl(); ++m)
00242 {
00243
00244
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
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
00278
00279 theCurrentTemp.entry.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
00280
00281
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
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
00357
00358 thePixelTemp_.push_back(theCurrentTemp);
00359
00360 return true;
00361
00362 }
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
00385
00386
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
00393 std::vector <float> chi2xavg(4), chi2xmin(4);
00394
00395
00396
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
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
00415
00416 Dtype_ = thePixelTemp_[index_id_].head.Dtype;
00417
00418
00419
00420 qscale_ = thePixelTemp_[index_id_].head.qscale;
00421
00422
00423
00424 s50_ = thePixelTemp_[index_id_].head.s50;
00425
00426
00427
00428 lorywidth_ = thePixelTemp_[index_id_].head.lorywidth;
00429 lorxwidth_ = thePixelTemp_[index_id_].head.lorxwidth;
00430
00431
00432
00433 xsize_ = thePixelTemp_[index_id_].head.xsize;
00434 ysize_ = thePixelTemp_[index_id_].head.ysize;
00435 zsize_ = thePixelTemp_[index_id_].head.zsize;
00436
00437
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
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
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
00514
00515 flip_y = false;
00516 if(cotbeta < 0.f) {flip_y = true;}
00517
00518
00519
00520 if(Dtype_ == 1) {
00521 if(cotbeta*locBz > 0.f) success_ = false;
00522 }
00523
00524
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
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
00578
00579
00580
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
00601
00602
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
00617
00618
00619
00620 ++pixy; ++pixx;
00621
00622
00623
00624 deltax = pixx - T2HX;
00625 deltay = pixy - T2HY;
00626
00627
00628
00629 for(j=0; j<BXM2; ++j) {for(i=0; i<BYM2; ++i) {xytemp_[j][i] = 0.f;}}
00630
00631
00632
00633 for(j=jmin; j<=jmax; ++j) {
00634 for(i=imin; i<=imax; ++i) {
00635 m = deltax+j;
00636
00637
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
00652
00653 for(n=1; n<BYM3; ++n) {
00654 if(ydouble[n-1]) {
00655
00656 for(m=1; m<BXM3; ++m) {
00657 xytemp_[m][n] += xytemp_[m][n+1];
00658 }
00659
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
00669
00670 for(m=1; m<BXM3; ++m) {
00671 if(xdouble[m-1]) {
00672
00673 for(n=1; n<BYM3; ++n) {
00674 xytemp_[m][n] += xytemp_[m+1][n];
00675 }
00676
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
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 }
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
00713
00714
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 }
00721
00722
00723
00724
00730
00731 void SiPixelTemplate2D::xysigma2(float qpixel, int index, float& xysig2)
00732
00733 {
00734
00735
00736
00737 float sigi, sigi2, sigi3, sigi4, qscale, err2, err00;
00738
00739
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
00750
00751
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 }
00779
00780
00781
00783
00784 void SiPixelTemplate2D::landau_par(float lanpar[2][5])
00785
00786 {
00787
00788
00789
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 }
00799
00800
00801