Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-06-04 06:14:14

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 }