SiPixelTransient/src/SiPixelGenError.cc

0001 //
0002 //  SiPixelGenError.cc  Version 2.30
0003 //
0004 //  Object stores Lorentz widths, bias corrections, and errors for the Generic Algorithm
0005 //
0006 //  Created by Morris Swartz on 10/27/06.
0007 //  Add some debugging messages. d.k. 5/14
0008 //  Update for Phase 1 FPix, M.S. 1/15/17
0009 //  V2.01 - Allow subdetector ID=5 for FPix R2P2, Fix error message
0010 //  V2.10 - Update the variable size [SI_PIXEL_TEMPLATE_USE_BOOST] option so that it works with VI's enhancements
0011 //  V2.20 - Add directory path selection to the ascii pushfile method
0012 //  V2.21 - Move templateStore to the heap, fix variable name in pushfile()
0013 //  V2.30 - Fix interpolation of IrradiationBias corrections
0014
0015 //#include <stdlib.h>
0016 //#include <stdio.h>
0017 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0018 #include <cmath>
0019 #else
0020 #include <math.h>
0021 #endif
0022 #include <algorithm>
0023 #include <vector>
0024 //#include "boost/multi_array.hpp"
0025 #include <iostream>
0026 #include <iomanip>
0027 #include <sstream>
0028 #include <fstream>
0029 #include <list>
0030
0031 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0032 #include "CondFormats/SiPixelTransient/interface/SiPixelGenError.h"
0033 #include "FWCore/Utilities/interface/FileInPath.h"
0034 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0035 #define LOGERROR(x) LogError(x)
0036 #define LOGWARNING(x) LogWarning(x)
0037 #define LOGINFO(x) LogInfo(x)
0038 #define ENDL " "
0039 #include "FWCore/Utilities/interface/Exception.h"
0040 using namespace edm;
0041 #else
0042 #include "SiPixelGenError.h"
0043 #define LOGERROR(x) std::cout << x << ": "
0044 #define LOGINFO(x) std::cout << x << ": "
0045 #define LOGWARNING(x) std::cout << x << ": "
0046 #define ENDL std::endl
0047 #endif
0048
0049 //****************************************************************
0050 //! This routine initializes the global GenError structures from
0051 //! an external file generror_summary_zpNNNN where NNNN are four
0052 //! digits of filenum.
0053 //! \param filenum - an integer NNNN used in the filename generror_summary_zpNNNN
0054 //****************************************************************
0055 bool SiPixelGenError::pushfile(int filenum, std::vector<SiPixelGenErrorStore>& pixelTemp, std::string dir) {
0056   // Add info stored in external file numbered filenum to theGenErrorStore
0057
0058   // Local variables
0059   int i, j, k;
0060   float costrk[3] = {0, 0, 0};
0061   const char* tempfile;
0062   //    char title[80]; remove this
0063   char c;
0064   const int code_version = {1};
0065
0066   //  Create a filename for this run
0067
0068   std::ostringstream tout;
0069
0070   //  Create different path in CMSSW than standalone
0071
0072 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0073   tout << dir << "generror_summary_zp" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out"
0074        << std::ends;
0075   std::string tempf = tout.str();
0076   edm::FileInPath file(tempf.c_str());
0077   tempfile = (file.fullPath()).c_str();
0078 #else
0079   tout << "generror_summary_zp" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
0080   std::string tempf = tout.str();
0081   tempfile = tempf.c_str();
0082 #endif
0083
0084   //  open the Generic file
0085
0086   std::ifstream in_file(tempfile, std::ios::in);
0087
0088   if (in_file.is_open()) {
0089     // Create a local GenError storage entry
0090
0091     SiPixelGenErrorStore theCurrentTemp;
0092
0093     // Read-in a header string first and print it
0094
0095     for (i = 0; (c = in_file.get()) != '\n'; ++i) {
0096       if (i < 79) {
0097         theCurrentTemp.head.title[i] = c;
0098       }
0099     }
0100     if (i > 78) {
0101       i = 78;
0102     }
0103     theCurrentTemp.head.title[i + 1] = '\0';
0104     LOGINFO("SiPixelGenError") << "Loading Pixel GenError File - " << theCurrentTemp.head.title << ENDL;
0105
0106     // next, the header information
0107
0108     in_file >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >>
0109         theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx >> theCurrentTemp.head.Dtype >>
0110         theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >>
0111         theCurrentTemp.head.qscale >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >>
0112         theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >>
0113         theCurrentTemp.head.zsize >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >>
0114         theCurrentTemp.head.lorxbias >> theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >>
0115         theCurrentTemp.head.fbin[2];
0116
0117     if (in_file.fail()) {
0118       LOGERROR("SiPixelGenError") << "Error reading file, no GenError load" << ENDL;
0119       return false;
0120     }
0121
0122     LOGINFO("SiPixelGenError") << "GenError ID = " << theCurrentTemp.head.ID << ", GenError Version "
0123                                << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
0124                                << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx
0125                                << ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
0126                                << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
0127                                << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence
0128                                << ", Q-scaling factor " << theCurrentTemp.head.qscale << ", 1/2 multi dcol threshold "
0129                                << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold " << theCurrentTemp.head.ss50
0130                                << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", y Lorentz Bias "
0131                                << theCurrentTemp.head.lorybias << ", x Lorentz width " << theCurrentTemp.head.lorxwidth
0132                                << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
0133                                << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] << ", "
0134                                << theCurrentTemp.head.fbin[1] << ", " << theCurrentTemp.head.fbin[2]
0135                                << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size "
0136                                << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
0137
0138     if (theCurrentTemp.head.templ_version < code_version) {
0139       LOGERROR("SiPixelGenError") << "code expects version " << code_version << ", no GenError load" << ENDL;
0140       return false;
0141     }
0142
0143 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST
0144
0145     // next, layout the 1-d/2-d structures needed to store GenError info
0146
0147     theCurrentTemp.cotbetaY = new float[theCurrentTemp.head.NTy];
0148     theCurrentTemp.cotbetaX = new float[theCurrentTemp.head.NTyx];
0149     theCurrentTemp.cotalphaX = new float[theCurrentTemp.head.NTxx];
0150
0151     theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
0152
0153     theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
0154
0155 #endif
0156
0157     // next, loop over all y-angle entries
0158
0159     for (i = 0; i < theCurrentTemp.head.NTy; ++i) {
0160       in_file >> theCurrentTemp.enty[i].runnum >> costrk[0] >> costrk[1] >> costrk[2];
0161
0162       if (in_file.fail()) {
0163         LOGERROR("SiPixelGenError") << "Error reading file 1, no GenError load, run # " << theCurrentTemp.enty[i].runnum
0164                                     << ENDL;
0165         return false;
0166       }
0167
0168       // Calculate the alpha, beta, and cot(beta) for this entry
0169
0170       theCurrentTemp.enty[i].cotalpha = costrk[0] / costrk[2];
0171
0172       theCurrentTemp.enty[i].cotbeta = costrk[1] / costrk[2];
0173
0174       in_file >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].pixmax >> theCurrentTemp.enty[i].dyone >>
0175           theCurrentTemp.enty[i].syone >> theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone;
0176
0177       if (in_file.fail()) {
0178         LOGERROR("SiPixelGenError") << "Error reading file 2, no GenError load, run # " << theCurrentTemp.enty[i].runnum
0179                                     << ENDL;
0180         return false;
0181       }
0182
0183       in_file >> theCurrentTemp.enty[i].dytwo >> theCurrentTemp.enty[i].sytwo >> theCurrentTemp.enty[i].dxtwo >>
0184           theCurrentTemp.enty[i].sxtwo >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].qmin2;
0185
0186       if (in_file.fail()) {
0187         LOGERROR("SiPixelGenError") << "Error reading file 3, no GenError load, run # " << theCurrentTemp.enty[i].runnum
0188                                     << ENDL;
0189         return false;
0190       }
0191
0192       if (theCurrentTemp.enty[i].qmin <= 0.) {
0193         LOGERROR("SiPixelGenError") << "Error in GenError ID " << theCurrentTemp.head.ID
0194                                     << " qmin = " << theCurrentTemp.enty[i].qmin << ", run # "
0195                                     << theCurrentTemp.enty[i].runnum << ENDL;
0196         return false;
0197       }
0198
0199       for (j = 0; j < 4; ++j) {
0200         in_file >> theCurrentTemp.enty[i].yavggen[j] >> theCurrentTemp.enty[i].yrmsgen[j] >>
0201             theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j];
0202
0203         if (in_file.fail()) {
0204           LOGERROR("SiPixelGenError") << "Error reading file 14a, no GenError load, run # "
0205                                       << theCurrentTemp.enty[i].runnum << ENDL;
0206           return false;
0207         }
0208       }
0209     }
0210
0211     // next, loop over all barrel x-angle entries
0212
0213     for (k = 0; k < theCurrentTemp.head.NTyx; ++k) {
0214       for (i = 0; i < theCurrentTemp.head.NTxx; ++i) {
0215         in_file >> theCurrentTemp.entx[k][i].runnum >> costrk[0] >> costrk[1] >> costrk[2];
0216
0217         if (in_file.fail()) {
0218           LOGERROR("SiPixelGenError") << "Error reading file 17, no GenError load, run # "
0219                                       << theCurrentTemp.entx[k][i].runnum << ENDL;
0220           return false;
0221         }
0222
0223         // Calculate the alpha, beta, and cot(beta) for this entry
0224
0225         theCurrentTemp.entx[k][i].cotalpha = costrk[0] / costrk[2];
0226
0227         theCurrentTemp.entx[k][i].cotbeta = costrk[1] / costrk[2];
0228
0229         in_file >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].pixmax >>
0230             theCurrentTemp.entx[k][i].dyone >> theCurrentTemp.entx[k][i].syone >> theCurrentTemp.entx[k][i].dxone >>
0231             theCurrentTemp.entx[k][i].sxone;
0232
0233         if (in_file.fail()) {
0234           LOGERROR("SiPixelGenError") << "Error reading file 18, no GenError load, run # "
0235                                       << theCurrentTemp.entx[k][i].runnum << ENDL;
0236           return false;
0237         }
0238
0239         in_file >> theCurrentTemp.entx[k][i].dytwo >> theCurrentTemp.entx[k][i].sytwo >>
0240             theCurrentTemp.entx[k][i].dxtwo >> theCurrentTemp.entx[k][i].sxtwo >> theCurrentTemp.entx[k][i].qmin >>
0241             theCurrentTemp.entx[k][i].qmin2;
0242         //             >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav;
0243
0244         if (in_file.fail()) {
0245           LOGERROR("SiPixelGenError") << "Error reading file 19, no GenError load, run # "
0246                                       << theCurrentTemp.entx[k][i].runnum << ENDL;
0247           return false;
0248         }
0249
0250         for (j = 0; j < 4; ++j) {
0251           in_file >> theCurrentTemp.entx[k][i].yavggen[j] >> theCurrentTemp.entx[k][i].yrmsgen[j] >>
0252               theCurrentTemp.entx[k][i].xavggen[j] >> theCurrentTemp.entx[k][i].xrmsgen[j];
0253
0254           if (in_file.fail()) {
0255             LOGERROR("SiPixelGenError") << "Error reading file 30a, no GenError load, run # "
0256                                         << theCurrentTemp.entx[k][i].runnum << ENDL;
0257             return false;
0258           }
0259         }
0260       }
0261     }
0262
0263     in_file.close();
0264
0265     // Add this info to the store
0266
0267     pixelTemp.push_back(theCurrentTemp);
0268
0269     postInit(pixelTemp);
0270
0271     return true;
0272
0273   } else {
0274     // If file didn't open, report this
0275
0276     LOGERROR("SiPixelGenError") << "Error opening File" << tempfile << ENDL;
0277     return false;
0278   }
0279
0280 }  // TempInit
0281
0282 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0283
0284 //****************************************************************
0285 //! This routine initializes the global GenError structures from an
0286 //! SiPixelGenErrorDBObject
0287 //! \param dbobject - db storing multiple generic calibrations
0288 //****************************************************************
0289 bool SiPixelGenError::pushfile(const SiPixelGenErrorDBObject& dbobject, std::vector<SiPixelGenErrorStore>& pixelTemp) {
0290   // Add GenError stored in external dbobject to theGenErrorStore
0291
0292   // Local variables
0293   int i, j, k;
0294   float costrk[3] = {0, 0, 0};
0295   //    const char *tempfile;
0296   const int code_version = {1};
0297
0298   // We must create a new object because dbobject must be a const and our stream must not be
0299   SiPixelGenErrorDBObject db = dbobject;
0300
0301   // Create a local GenError storage entry
0302   /// SiPixelGenErrorStore theCurrentTemp;    // not on the stack...
0303   auto tmpPtr = std::make_unique<SiPixelGenErrorStore>();  // must be allocated on the heap instead
0304   auto& theCurrentTemp = *tmpPtr;
0305
0306   // Fill the GenError storage for each GenError calibration stored in the db
0307   for (int m = 0; m < db.numOfTempl(); ++m) {
0308     // Read-in a header string first and print it
0309
0310     SiPixelGenErrorDBObject::char2float temp;
0311     for (i = 0; i < 20; ++i) {
0312       temp.f = db.sVector()[db.index()];
0313       theCurrentTemp.head.title[4 * i] = temp.c[0];
0314       theCurrentTemp.head.title[4 * i + 1] = temp.c[1];
0315       theCurrentTemp.head.title[4 * i + 2] = temp.c[2];
0316       theCurrentTemp.head.title[4 * i + 3] = temp.c[3];
0317       db.incrementIndex(1);
0318     }
0319
0320     theCurrentTemp.head.title[79] = '\0';
0321     LOGINFO("SiPixelGenError") << "Loading Pixel GenError File - " << theCurrentTemp.head.title << ENDL;
0322
0323     // next, the header information
0324
0325     db >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >>
0326         theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx >> theCurrentTemp.head.Dtype >>
0327         theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >>
0328         theCurrentTemp.head.qscale >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >>
0329         theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >>
0330         theCurrentTemp.head.zsize >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >>
0331         theCurrentTemp.head.lorxbias >> theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >>
0332         theCurrentTemp.head.fbin[2];
0333
0334     if (db.fail()) {
0335       LOGERROR("SiPixelGenError") << "Error reading file, no GenError load" << ENDL;
0336       return false;
0337     }
0338
0339     LOGINFO("SiPixelGenError") << "GenError ID = " << theCurrentTemp.head.ID << ", GenError Version "
0340                                << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
0341                                << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx
0342                                << ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
0343                                << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
0344                                << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence
0345                                << ", Q-scaling factor " << theCurrentTemp.head.qscale << ", 1/2 multi dcol threshold "
0346                                << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold " << theCurrentTemp.head.ss50
0347                                << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", y Lorentz Bias "
0348                                << theCurrentTemp.head.lorybias << ", x Lorentz width " << theCurrentTemp.head.lorxwidth
0349                                << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
0350                                << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] << ", "
0351                                << theCurrentTemp.head.fbin[1] << ", " << theCurrentTemp.head.fbin[2]
0352                                << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size "
0353                                << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
0354
0355     LOGINFO("SiPixelGenError") << "Loading Pixel GenError - " << theCurrentTemp.head.title << " version "
0356                                << theCurrentTemp.head.templ_version << " code v." << code_version << ENDL;
0357     if (theCurrentTemp.head.templ_version < code_version) {
0358       LOGERROR("SiPixelGenError") << "code expects version " << code_version << ", no GenError load" << ENDL;
0359       return false;
0360     }
0361
0362 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST
0363
0364     // next, layout the 1-d/2-d structures needed to store GenError
0365
0366     // &&&  Who is going to delete these?  Are we leaking memory?
0367     theCurrentTemp.cotbetaY = new float[theCurrentTemp.head.NTy];
0368     theCurrentTemp.cotbetaX = new float[theCurrentTemp.head.NTyx];
0369     theCurrentTemp.cotalphaX = new float[theCurrentTemp.head.NTxx];
0370
0371     theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
0372
0373     theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
0374
0375 #endif
0376
0377     // next, loop over all barrel y-angle entries
0378
0379     for (i = 0; i < theCurrentTemp.head.NTy; ++i) {
0380       db >> theCurrentTemp.enty[i].runnum >> costrk[0] >> costrk[1] >> costrk[2];
0381
0382       if (db.fail()) {
0383         LOGERROR("SiPixelGenError") << "Error reading file 1, no GenError load, run # " << theCurrentTemp.enty[i].runnum
0384                                     << ENDL;
0385         return false;
0386       }
0387
0388       // Calculate the alpha, beta, and cot(beta) for this entry
0389
0390       theCurrentTemp.enty[i].cotalpha = costrk[0] / costrk[2];
0391
0392       theCurrentTemp.enty[i].cotbeta = costrk[1] / costrk[2];
0393
0394       db >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].pixmax >> theCurrentTemp.enty[i].dyone >>
0395           theCurrentTemp.enty[i].syone >> theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone;
0396
0397       if (db.fail()) {
0398         LOGERROR("SiPixelGenError") << "Error reading file 2, no GenError load, run # " << theCurrentTemp.enty[i].runnum
0399                                     << ENDL;
0400         return false;
0401       }
0402
0403       db >> theCurrentTemp.enty[i].dytwo >> theCurrentTemp.enty[i].sytwo >> theCurrentTemp.enty[i].dxtwo >>
0404           theCurrentTemp.enty[i].sxtwo >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].qmin2;
0405
0406       for (j = 0; j < 4; ++j) {
0407         db >> theCurrentTemp.enty[i].yavggen[j] >> theCurrentTemp.enty[i].yrmsgen[j] >>
0408             theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j];
0409
0410         if (db.fail()) {
0411           LOGERROR("SiPixelGenError") << "Error reading file 14a, no GenError load, run # "
0412                                       << theCurrentTemp.enty[i].runnum << ENDL;
0413           return false;
0414         }
0415       }
0416     }
0417
0418     // next, loop over all barrel x-angle entries
0419
0420     for (k = 0; k < theCurrentTemp.head.NTyx; ++k) {
0421       for (i = 0; i < theCurrentTemp.head.NTxx; ++i) {
0422         db >> theCurrentTemp.entx[k][i].runnum >> costrk[0] >> costrk[1] >> costrk[2];
0423
0424         if (db.fail()) {
0425           LOGERROR("SiPixelGenError") << "Error reading file 17, no GenError load, run # "
0426                                       << theCurrentTemp.entx[k][i].runnum << ENDL;
0427           return false;
0428         }
0429
0430         // Calculate the alpha, beta, and cot(beta) for this entry
0431
0432         theCurrentTemp.entx[k][i].cotalpha = costrk[0] / costrk[2];
0433
0434         theCurrentTemp.entx[k][i].cotbeta = costrk[1] / costrk[2];
0435
0436         db >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].pixmax >> theCurrentTemp.entx[k][i].dyone >>
0437             theCurrentTemp.entx[k][i].syone >> theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone;
0438
0439         if (db.fail()) {
0440           LOGERROR("SiPixelGenError") << "Error reading file 18, no GenError load, run # "
0441                                       << theCurrentTemp.entx[k][i].runnum << ENDL;
0442           return false;
0443         }
0444
0445         db >> theCurrentTemp.entx[k][i].dytwo >> theCurrentTemp.entx[k][i].sytwo >> theCurrentTemp.entx[k][i].dxtwo >>
0446             theCurrentTemp.entx[k][i].sxtwo >> theCurrentTemp.entx[k][i].qmin >> theCurrentTemp.entx[k][i].qmin2;
0447
0448         for (j = 0; j < 4; ++j) {
0449           db >> theCurrentTemp.entx[k][i].yavggen[j] >> theCurrentTemp.entx[k][i].yrmsgen[j] >>
0450               theCurrentTemp.entx[k][i].xavggen[j] >> theCurrentTemp.entx[k][i].xrmsgen[j];
0451
0452           if (db.fail()) {
0453             LOGERROR("SiPixelGenError") << "Error reading file 30a, no GenError load, run # "
0454                                         << theCurrentTemp.entx[k][i].runnum << ENDL;
0455             return false;
0456           }
0457         }
0458       }
0459     }
0460
0461     // Add this GenError to the store
0462
0463     pixelTemp.push_back(theCurrentTemp);
0464
0465     postInit(pixelTemp);
0466   }
0467   return true;
0468
0469 }  // TempInit
0470
0471 #endif
0472
0473 void SiPixelGenError::postInit(std::vector<SiPixelGenErrorStore>& thePixelTemp_) {
0474   for (auto& templ : thePixelTemp_) {
0475     for (auto iy = 0; iy < templ.head.NTy; ++iy)
0476       templ.cotbetaY[iy] = templ.enty[iy].cotbeta;
0477     for (auto iy = 0; iy < templ.head.NTyx; ++iy)
0478       templ.cotbetaX[iy] = templ.entx[iy][0].cotbeta;
0479     for (auto ix = 0; ix < templ.head.NTxx; ++ix)
0480       templ.cotalphaX[ix] = templ.entx[0][ix].cotalpha;
0481   }
0482 }
0483
0484 // ************************************************************************************************************
0485 //! Interpolate beta/alpha angles to produce an expected average charge. Return int (0-4) describing the charge
0486 //! of the cluster [0: 1.5<Q/Qavg, 1: 1<Q/Qavg<1.5, 2: 0.85<Q/Qavg<1, 3: 0.95Qmin<Q<0.85Qavg, 4: Q<0.95Qmin].
0487 //! \param id - (input) index of the GenError to use
0488 //! \param cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)
0489 //! \param cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)
0490 //! \param locBz - (input) the sign of this quantity is used to determine whether to flip cot(beta)<0 quantities from cot(beta)>0 (FPix only)
0491 //!                    for FPix IP-related tracks, locBz < 0 for cot(beta) > 0 and locBz > 0 for cot(beta) < 0
0492 //! \param qclus - (input) the cluster charge in electrons
0493 //! \param pixmax - (output) the maximum pixel charge in electrons (truncation value)
0494 //! \param sigmay - (output) the estimated y-error for CPEGeneric in microns
0495 //! \param deltay - (output) the estimated y-bias for CPEGeneric in microns
0496 //! \param sigmax - (output) the estimated x-error for CPEGeneric in microns
0497 //! \param deltax - (output) the estimated x-bias for CPEGeneric in microns
0498 //! \param sy1 - (output) the estimated y-error for 1 single-pixel clusters in microns
0499 //! \param dy1 - (output) the estimated y-bias for 1 single-pixel clusters in microns
0500 //! \param sy2 - (output) the estimated y-error for 1 double-pixel clusters in microns
0501 //! \param dy2 - (output) the estimated y-bias for 1 double-pixel clusters in microns
0502 //! \param sx1 - (output) the estimated x-error for 1 single-pixel clusters in microns
0503 //! \param dx1 - (output) the estimated x-bias for 1 single-pixel clusters in microns
0504 //! \param sx2 - (output) the estimated x-error for 1 double-pixel clusters in microns
0505 //! \param dx2 - (output) the estimated x-bias for 1 double-pixel clusters in microns
0506 // ************************************************************************************************************
0507 // a simpler method just to return the LA
0508 int SiPixelGenError::qbin(int id) {
0509   // Find the index corresponding to id
0510
0511   if (id != id_current_) {
0512     index_id_ = -1;
0513     for (int i = 0; i < (int)thePixelTemp_.size(); ++i) {
0514       if (id == thePixelTemp_[i].head.ID) {
0515         index_id_ = i;
0516         id_current_ = id;
0517         //
0518         lorywidth_ = thePixelTemp_[i].head.lorywidth;
0519         lorxwidth_ = thePixelTemp_[i].head.lorxwidth;
0520         lorybias_ = thePixelTemp_[i].head.lorybias;
0521         lorxbias_ = thePixelTemp_[i].head.lorxbias;
0522
0523         //for(int j=0; j<3; ++j) {fbin_[j] = thePixelTemp_[i].head.fbin[j];}
0524
0525         // Pixel sizes to the private variables
0526         xsize_ = thePixelTemp_[i].head.xsize;
0527         ysize_ = thePixelTemp_[i].head.ysize;
0528         zsize_ = thePixelTemp_[i].head.zsize;
0529
0530         break;
0531       }
0532     }
0533   }
0534   return index_id_;
0535 }
0536 //-----------------------------------------------------------------------
0537 // Full method
0538 int SiPixelGenError::qbin(int id,
0539                           float cotalpha,
0540                           float cotbeta,
0541                           float locBz,
0542                           float locBx,
0543                           float qclus,
0544                           bool irradiationCorrections,
0545                           int& pixmx,
0546                           float& sigmay,
0547                           float& deltay,
0548                           float& sigmax,
0549                           float& deltax,
0550                           float& sy1,
0551                           float& dy1,
0552                           float& sy2,
0553                           float& dy2,
0554                           float& sx1,
0555                           float& dx1,
0556                           float& sx2,
0557                           float& dx2) {
0558   // Interpolate for a new set of track angles
0559
0560   // Find the index corresponding to id
0561
0562   if (id != id_current_) {
0563     index_id_ = -1;
0564     for (int i = 0; i < (int)thePixelTemp_.size(); ++i) {
0565       if (id == thePixelTemp_[i].head.ID) {
0566         index_id_ = i;
0567         id_current_ = id;
0568         lorywidth_ = thePixelTemp_[i].head.lorywidth;
0569         lorxwidth_ = thePixelTemp_[i].head.lorxwidth;
0570         lorybias_ = thePixelTemp_[i].head.lorybias;
0571         lorxbias_ = thePixelTemp_[i].head.lorxbias;
0572         for (int j = 0; j < 3; ++j) {
0573           fbin_[j] = thePixelTemp_[i].head.fbin[j];
0574         }
0575
0576         // Pixel sizes to the private variables
0577
0578         xsize_ = thePixelTemp_[i].head.xsize;
0579         ysize_ = thePixelTemp_[i].head.ysize;
0580         zsize_ = thePixelTemp_[i].head.zsize;
0581
0582         break;
0583       }
0584     }
0585   }
0586
0587   int index = index_id_;
0588
0589 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0590   if (index < 0 || index >= (int)thePixelTemp_.size()) {
0591     throw cms::Exception("DataCorrupt") << "SiPixelGenError::qbin can't find needed GenError ID = " << id << std::endl;
0592   }
0593 #else
0594   assert(index >= 0 && index < (int)thePixelTemp_.size());
0595 #endif
0596
0597   //
0598
0599   auto const& templ = thePixelTemp_[index];
0600
0601   // Interpolate the absolute value of cot(beta)
0602
0603   auto acotb = std::abs(cotbeta);
0604
0605   //    qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)
0606
0607   auto cotalpha0 = thePixelTemp_[index].enty[0].cotalpha;
0608   auto qcorrect =
0609       std::sqrt((1.f + cotbeta * cotbeta + cotalpha * cotalpha) / (1.f + cotbeta * cotbeta + cotalpha0 * cotalpha0));
0610
0611   // for some cosmics, the ususal gymnastics are incorrect
0612
0613   float cota = cotalpha;
0614   float cotb = acotb;
0615   bool flip_y;
0616   bool flip_x;
0617   // for some cosmics, the ususal gymnastics are incorrect
0618   flip_x = false;
0619   flip_y = false;
0620   switch (thePixelTemp_[index_id_].head.Dtype) {
0621     case 0:
0622       if (cotbeta < 0.f) {
0623         flip_y = true;
0624       }
0625       break;
0626     case 1:
0627       if (locBz < 0.f) {
0628         cotb = cotbeta;
0629       } else {
0630         cotb = -cotbeta;
0631         flip_y = true;
0632       }
0633       break;
0634     case 2:
0635     case 3:
0636     case 4:
0637     case 5:
0638       if (locBx * locBz < 0.f) {
0639         cota = -cotalpha;
0640         flip_x = true;
0641       }
0642       if (locBx > 0.f) {
0643         cotb = cotbeta;
0644       } else {
0645         cotb = -cotbeta;
0646         flip_y = true;
0647       }
0648       break;
0649     default:
0650 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0651       throw cms::Exception("DataCorrupt")
0652           << "SiPixelGenError::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
0653 #else
0654       std::cout << "SiPixelGenError::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
0655 #endif
0656   }
0657
0658   // Copy the charge scaling factor to the private variable
0659
0660   if (flip_y) {
0661     lorybias_ = -lorybias_;
0662     lorywidth_ = -lorywidth_;
0663   }
0664   if (flip_x) {
0665     lorxbias_ = -lorxbias_;
0666     lorxwidth_ = -lorxwidth_;
0667   }
0668
0669   auto qscale = thePixelTemp_[index].head.qscale;
0670
0671   /*
0672     lorywidth = thePixelTemp_[index].head.lorywidth;
0673     if(locBz > 0.f) {lorywidth = -lorywidth;}
0674     lorxwidth = thePixelTemp_[index].head.lorxwidth;
0675     */
0676
0677   auto Ny = thePixelTemp_[index].head.NTy;
0678   auto Nyx = thePixelTemp_[index].head.NTyx;
0679   auto Nxx = thePixelTemp_[index].head.NTxx;
0680
0681 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0682   if (Ny < 2 || Nyx < 1 || Nxx < 2) {
0683     throw cms::Exception("DataCorrupt") << "GenError ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny
0684                                         << "/" << Nyx << "/" << Nxx << std::endl;
0685   }
0686 #else
0687   assert(Ny > 1 && Nyx > 0 && Nxx > 1);
0688 #endif
0689
0690   // next, loop over all y-angle entries
0691
0692   auto ilow = 0;
0693   auto ihigh = 0;
0694   auto yratio = 0.f;
0695
0696   {
0697     auto j = std::lower_bound(templ.cotbetaY, templ.cotbetaY + Ny, cotb);
0698     if (j == templ.cotbetaY + Ny) {
0699       --j;
0700       yratio = 1.f;
0701     } else if (j == templ.cotbetaY) {
0702       ++j;
0703       yratio = 0.f;
0704     } else {
0705       yratio = (cotb - (*(j - 1))) / ((*j) - (*(j - 1)));
0706     }
0707
0708     ihigh = j - templ.cotbetaY;
0709     ilow = ihigh - 1;
0710   }
0711
0712   // Interpolate/store all y-related quantities (flip displacements when flip_y)
0713
0714   auto qavg = (1.f - yratio) * thePixelTemp_[index].enty[ilow].qavg + yratio * thePixelTemp_[index].enty[ihigh].qavg;
0715   qavg *= qcorrect;
0716   auto qmin = (1.f - yratio) * thePixelTemp_[index].enty[ilow].qmin + yratio * thePixelTemp_[index].enty[ihigh].qmin;
0717   qmin *= qcorrect;
0718   auto qmin2 = (1.f - yratio) * thePixelTemp_[index].enty[ilow].qmin2 + yratio * thePixelTemp_[index].enty[ihigh].qmin2;
0719   qmin2 *= qcorrect;
0720
0721 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0722   if (qavg <= 0.f || qmin <= 0.f) {
0723     throw cms::Exception("DataCorrupt")
0724         << "SiPixelGenError::qbin, qavg or qmin <= 0,"
0725         << " Probably someone called the generic pixel reconstruction with an illegal trajectory state" << std::endl;
0726   }
0727 #else
0728   assert(qavg > 0.f && qmin > 0.f);
0729 #endif
0730
0731   //  Scale the input charge to account for differences between pixelav and CMSSW simulation or data
0732
0733   auto qtotal = qscale * qclus;
0734
0735   // uncertainty and final corrections depend upon total charge bin
0736
0737   auto fq = qtotal / qavg;
0738   int binq;
0739   if (fq > fbin_[0]) {
0740     binq = 0;
0741   } else {
0742     if (fq > fbin_[1]) {
0743       binq = 1;
0744     } else {
0745       if (fq > fbin_[2]) {
0746         binq = 2;
0747       } else {
0748         binq = 3;
0749       }
0750     }
0751   }
0752
0753   auto yrmsgen = (1.f - yratio) * thePixelTemp_[index].enty[ilow].yrmsgen[binq] +
0754                  yratio * thePixelTemp_[index].enty[ihigh].yrmsgen[binq];
0755   sy1 = (1.f - yratio) * thePixelTemp_[index].enty[ilow].syone + yratio * thePixelTemp_[index].enty[ihigh].syone;
0756   sy2 = (1.f - yratio) * thePixelTemp_[index].enty[ilow].sytwo + yratio * thePixelTemp_[index].enty[ihigh].sytwo;
0757
0758   if (irradiationCorrections) {
0759     auto yavggen = (1.f - yratio) * thePixelTemp_[index].enty[ilow].yavggen[binq] +
0760                    yratio * thePixelTemp_[index].enty[ihigh].yavggen[binq];
0761     if (flip_y) {
0762       yavggen = -yavggen;
0763     }
0764     deltay = yavggen;
0765     dy1 = (1.f - yratio) * thePixelTemp_[index].enty[ilow].dyone + yratio * thePixelTemp_[index].enty[ihigh].dyone;
0766     if (flip_y) {
0767       dy1 = -dy1;
0768     }
0769     dy2 = (1.f - yratio) * thePixelTemp_[index].enty[ilow].dytwo + yratio * thePixelTemp_[index].enty[ihigh].dytwo;
0770     if (flip_y) {
0771       dy2 = -dy2;
0772     }
0773   }
0774
0775   // next, loop over all x-angle entries, first, find relevant y-slices
0776
0777   auto iylow = 0;
0778   auto iyhigh = 0;
0779   auto yxratio = 0.f;
0780
0781   {
0782     auto j = std::lower_bound(templ.cotbetaX, templ.cotbetaX + Nyx, acotb);
0783     if (j == templ.cotbetaX + Nyx) {
0784       --j;
0785       yxratio = 1.f;
0786     } else if (j == templ.cotbetaX) {
0787       ++j;
0788       yxratio = 0.f;
0789     } else {
0790       yxratio = (acotb - (*(j - 1))) / ((*j) - (*(j - 1)));
0791     }
0792
0793     iyhigh = j - templ.cotbetaX;
0794     iylow = iyhigh - 1;
0795   }
0796
0797   auto xxratio = 0.f;
0798
0799   {
0800     auto j = std::lower_bound(templ.cotalphaX, templ.cotalphaX + Nxx, cota);
0801     if (j == templ.cotalphaX + Nxx) {
0802       --j;
0803       xxratio = 1.f;
0804     } else if (j == templ.cotalphaX) {
0805       ++j;
0806       xxratio = 0.f;
0807     } else {
0808       xxratio = (cota - (*(j - 1))) / ((*j) - (*(j - 1)));
0809     }
0810
0811     ihigh = j - templ.cotalphaX;
0812     ilow = ihigh - 1;
0813   }
0814
0815   sx1 =
0816       (1.f - xxratio) * thePixelTemp_[index].entx[0][ilow].sxone + xxratio * thePixelTemp_[index].entx[0][ihigh].sxone;
0817   sx2 =
0818       (1.f - xxratio) * thePixelTemp_[index].entx[0][ilow].sxtwo + xxratio * thePixelTemp_[index].entx[0][ihigh].sxtwo;
0819
0820   // pixmax is the maximum allowed pixel charge (used for truncation)
0821
0822   pixmx = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].pixmax +
0823                              xxratio * thePixelTemp_[index].entx[iylow][ihigh].pixmax) +
0824           yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].pixmax +
0825                      xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].pixmax);
0826
0827   auto xrmsgen = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].xrmsgen[binq] +
0828                                     xxratio * thePixelTemp_[index].entx[iylow][ihigh].xrmsgen[binq]) +
0829                  yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].xrmsgen[binq] +
0830                             xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].xrmsgen[binq]);
0831
0832   if (irradiationCorrections) {
0833     auto xavggen = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].xavggen[binq] +
0834                                       xxratio * thePixelTemp_[index].entx[iylow][ihigh].xavggen[binq]) +
0835                    yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].xavggen[binq] +
0836                               xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].xavggen[binq]);
0837     if (flip_x) {
0838       xavggen = -xavggen;
0839     }
0840     deltax = xavggen;
0841     dx1 = (1.f - xxratio) * thePixelTemp_[index].entx[0][ilow].dxone +
0842           xxratio * thePixelTemp_[index].entx[0][ihigh].dxone;
0843     if (flip_x) {
0844       dx1 = -dx1;
0845     }
0846     dx2 = (1.f - xxratio) * thePixelTemp_[index].entx[0][ilow].dxtwo +
0847           xxratio * thePixelTemp_[index].entx[0][ihigh].dxtwo;
0848     if (flip_x) {
0849       dx2 = -dx2;
0850     }
0851   }
0852
0853   //  Take the errors and bias from the correct charge bin
0854
0855   sigmay = yrmsgen;
0856
0857   sigmax = xrmsgen;
0858
0859   // If the charge is too small (then flag it)
0860
0861   if (qtotal < 0.95f * qmin) {
0862     binq = 5;
0863   } else {
0864     if (qtotal < 0.95f * qmin2) {
0865       binq = 4;
0866     }
0867   }
0868
0869   return binq;
0870
0871 }  // qbin
0872
0873 int SiPixelGenError::qbin(int id,
0874                           float cotalpha,
0875                           float cotbeta,
0876                           float locBz,
0877                           float locBx,
0878                           float qclus,
0879                           float& pixmx,
0880                           float& sigmay,
0881                           float& deltay,
0882                           float& sigmax,
0883                           float& deltax,
0884                           float& sy1,
0885                           float& dy1,
0886                           float& sy2,
0887                           float& dy2,
0888                           float& sx1,
0889                           float& dx1,
0890                           float& sx2,
0891                           float& dx2) {
0892   // Interpolate for a new set of track angles
0893
0894   bool irradiationCorrections = true;
0895   int ipixmx, ibin;
0896
0897   ibin = SiPixelGenError::qbin(id,
0898                                cotalpha,
0899                                cotbeta,
0900                                locBz,
0901                                locBx,
0902                                qclus,
0903                                irradiationCorrections,
0904                                ipixmx,
0905                                sigmay,
0906                                deltay,
0907                                sigmax,
0908                                deltax,
0909                                sy1,
0910                                dy1,
0911                                sy2,
0912                                dy2,
0913                                sx1,
0914                                dx1,
0915                                sx2,
0916                                dx2);
0917
0918   pixmx = (float)ipixmx;
0919
0920   return ibin;
0921 }