Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:02:36

0001 //
0002 //  SiPixelTemplate2D.h (v2.65)
0003 //
0004 //  Full 2-D templates for cluster splitting, simulated cluster reweighting, and improved cluster probability
0005 //
0006 // Created by Morris Swartz on 12/01/09.
0007 // V1.01 - fix qavg_ filling
0008 // V1.02 - Add locBz to test if FPix use is out of range
0009 // V1.03 - Fix edge checking on final template to increase template size and to properly truncate cluster
0010 // v2.00 - Major changes to accommodate 2D reconstruction
0011 // v2.10 - Change chi2 and error scaling information to work with partially reconstructed clusters
0012 // v2.20 - Add cluster charge probability information, side loading for template generation
0013 // v2.21 - Double derivative interval [improves fit convergence]
0014 // v2.25 - Resize template store to accommodate FPix Templates
0015 // v2.30 - Fix bug found by P. Shuetze that compromises sqlite file loading
0016 // v2.35 - Add directory path selection to the ascii pushfile method
0017 // v2.50 - Change template storage to dynamically allocated 2D arrays of SiPixelTemplateEntry2D structs
0018 // v2.51 - Ensure that the derivative arrays are correctly zeroed between calls
0019 // v2.52 - Improve cosmetics for increased style points from judges
0020 // v2.60 - Fix FPix multiframe lookup problem [takes +-cotalpha and +-cotbeta]
0021 // v2.61a - Code 2.60 fix correctly
0022 // v2.65 - change double pixel flags to work with new shifted reco code definition
0023 //
0024 
0025 // Build the template storage structure from several pieces
0026 
0027 #ifndef SiPixelTemplate2D_h
0028 #define SiPixelTemplate2D_h 1
0029 
0030 #include <vector>
0031 #include <cassert>
0032 #include "boost/multi_array.hpp"
0033 
0034 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0035 #include "CondFormats/SiPixelObjects/interface/SiPixel2DTemplateDBObject.h"
0036 #include "FWCore/Utilities/interface/Exception.h"
0037 #include "CondFormats/SiPixelTransient/interface/SiPixelTemplateDefs.h"
0038 #else
0039 #include "SiPixelTemplateDefs.h"
0040 #endif
0041 
0042 struct SiPixelTemplateEntry2D {  //!< Basic template entry corresponding to a single set of track angles
0043   int runnum;                    //!< number of pixelav run used to generate this entry
0044   float cotalpha;                //!< cot(alpha) is proportional to cluster length in x and is basis of interpolation
0045   float cotbeta;                 //!< cot(beta) is proportional to cluster length in y and is basis of interpolation
0046   float costrk[3];               //!< direction cosines of tracks used to generate this entry
0047   float qavg;                    //!< average cluster charge for this set of track angles
0048   float pixmax;                  //!< maximum charge for individual pixels in cluster
0049   float sxymax;                  //!< average pixel signal for use of the error parameterization
0050   int iymin;                     //!< the minimum nonzero pixel yindex in template (saves time during interpolation)
0051   int iymax;                     //!< the maximum nonzero pixel yindex in template (saves time during interpolation)
0052   int jxmin;                     //!< the minimum nonzero pixel xindex in template (saves time during interpolation)
0053   int jxmax;                     //!< the maximum nonzero pixel xindex in template (saves time during interpolation)
0054   float xypar[2][5];             //!< pixel uncertainty parameterization
0055   float lanpar[2][5];            //!< pixel landau distribution parameters
0056   short int xytemp[7][7][T2YSIZE][T2XSIZE];  //!< templates for y-reconstruction (binned over 1 central pixel)
0057   float chi2ppix;                            //!< average chi^2 per pixel
0058   float chi2scale;                           //!< scale factor for the chi2 distribution
0059   float chi2avgone;                          //!< average y chi^2 for 1 pixel clusters
0060   float chi2minone;                          //!< minimum of y chi^2 for 1 pixel clusters
0061   float clsleny;                             //!< cluster y-length in pixels at signal height symax/2
0062   float clslenx;                             //!< cluster x-length in pixels at signal height sxmax/2
0063   float mpvvav;      //!< most probable charge in Vavilov distribution (not actually for larger kappa)
0064   float sigmavav;    //!< "sigma" scale fctor for Vavilov distribution
0065   float kappavav;    //!< kappa parameter for Vavilov distribution
0066   float scalexavg;   //!< average x-error scale factor
0067   float scaleyavg;   //!< average y-error scale factor
0068   float delyavg;     //!< average length difference between template and cluster
0069   float delysig;     //!< rms of length difference between template and cluster
0070   float scalex[4];   //!< x-error scale factor in 4 charge bins
0071   float scaley[4];   //!< y-error scale factor in 4 charge bins
0072   float offsetx[4];  //!< x-offset in 4 charge bins
0073   float offsety[4];  //!< y-offset in 4 charge bins
0074   float spare[3];
0075 };
0076 
0077 struct SiPixelTemplateHeader2D {  //!< template header structure
0078   int ID;                         //!< template ID number
0079   int NTy;                        //!< number of Template y entries
0080   int NTyx;                       //!< number of Template y-slices of x entries
0081   int NTxx;                       //!< number of Template x-entries in each slice
0082   int Dtype;                      //!< detector type (0=BPix, 1=FPix)
0083   float qscale;                   //!< Charge scaling to match cmssw and pixelav
0084   float lorywidth;                //!< estimate of y-lorentz width for optimal resolution
0085   float lorxwidth;                //!< estimate of x-lorentz width for optimal resolution
0086   float lorybias;                 //!< estimate of y-lorentz bias
0087   float lorxbias;                 //!< estimate of x-lorentz bias
0088   float Vbias;                    //!< detector bias potential in Volts
0089   float temperature;              //!< detector temperature in deg K
0090   float fluence;                  //!< radiation fluence in n_eq/cm^2
0091   float s50;                      //!< 1/2 of the multihit dcol threshold in electrons
0092   float ss50;                     //!< 1/2 of the single hit dcol threshold in electrons
0093   char title[80];                 //!< template title
0094   int templ_version;              //!< Version number of the template to ensure code compatibility
0095   float Bfield;                   //!< Bfield in Tesla
0096   float fbin[3];                  //!< The QBin definitions in Q_clus/Q_avg
0097   float xsize;                    //!< pixel size (for future use in upgraded geometry)
0098   float ysize;                    //!< pixel size (for future use in upgraded geometry)
0099   float zsize;                    //!< pixel size (for future use in upgraded geometry)
0100 };
0101 
0102 struct SiPixelTemplateStore2D {  //!< template storage structure
0103 
0104   void resize(int ny, int nx) {
0105     entry.resize(ny);
0106     store.resize(nx * ny);
0107     int off = 0;
0108     for (int i = 0; i < ny; ++i) {
0109       entry[i] = store.data() + off;
0110       off += nx;
0111     }
0112     assert(nx * ny == off);
0113   }
0114 
0115   SiPixelTemplateHeader2D head;  //!< Header information
0116 
0117   //!< use 2d entry to store BPix and FPix entries [dynamically allocated
0118   std::vector<SiPixelTemplateEntry2D*> entry;
0119 
0120   std::vector<SiPixelTemplateEntry2D> store;
0121 };
0122 
0123 // ******************************************************************************************
0124 //! \class SiPixelTemplate2D
0125 //!
0126 //!  A template management class.  SiPixelTemplate contains thePixelTemp
0127 //!  (a std::vector  of SiPixelTemplateStore, each of which is a collection of many
0128 //!  SiPixelTemplateEntries).  Each SiPixelTemplateStore corresponds to a given detector
0129 //!  condition, and is valid for a range of runs.  We allow more than one Store since the
0130 //!  may change over time.
0131 //!
0132 //!  This class reads templates from files via pushfile() method.
0133 //!
0134 //!  The main functionality of SiPixelTemplate is xytemp(), which produces a template
0135 //!  on the fly, given a specific track's alpha and beta.  The results are kept in data
0136 //!  members and accessed via inline getters.
0137 //!
0138 //!  The resulting template is then used by PixelTempReco2D() (a global function) which
0139 //!  get the reference for SiPixelTemplate & templ and uses the current template to
0140 //!  reconstruct the SiPixelRecHit.
0141 // ******************************************************************************************
0142 class SiPixelTemplate2D {
0143 public:
0144   SiPixelTemplate2D(const std::vector<SiPixelTemplateStore2D>& thePixelTemp) : thePixelTemp_(thePixelTemp) {
0145     id_current_ = -1;
0146     index_id_ = -1;
0147     cota_current_ = 0.;
0148     cotb_current_ = 0.;
0149   }  //!< Default constructor
0150 
0151   // load the private store with info from the
0152   // file with the index (int) filenum ${dir}template_summary_zp${filenum}.out
0153 #ifdef SI_PIXEL_TEMPLATE_STANDALONE
0154   static bool pushfile(int filenum, std::vector<SiPixelTemplateStore2D>& pixelTemp, std::string dir = "");
0155 
0156   // For calibrations only: load precalculated values -- no interpolation.
0157   void sideload(SiPixelTemplateEntry2D* entry,
0158                 int iDtype,
0159                 float locBx,
0160                 float locBz,
0161                 float lorwdy,
0162                 float lorwdx,
0163                 float q50,
0164                 float fbin[3],
0165                 float xsize,
0166                 float ysize,
0167                 float zsize);
0168 
0169 #else
0170   static bool pushfile(int filenum,
0171                        std::vector<SiPixelTemplateStore2D>& pixelTemp,
0172                        std::string dir = "CalibTracker/SiPixelESProducers/data/");
0173 
0174   // Load from the DB (the default in CMSSW):
0175   static bool pushfile(const SiPixel2DTemplateDBObject& dbobject, std::vector<SiPixelTemplateStore2D>& pixelTemp);
0176 
0177 #endif
0178 
0179   //  Initialize things before interpolating
0180   bool getid(int id);
0181 
0182   bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx);
0183 
0184   // Interpolate input alpha and beta angles to produce a working template for each individual hit.
0185 
0186   // Works with Phase 0+1
0187   bool xytemp(float xhit,
0188               float yhit,
0189               bool ydouble[BYM2],
0190               bool xdouble[BXM2],
0191               float template2d[BXM2][BYM2],
0192               bool dervatives,
0193               float dpdx2d[2][BXM2][BYM2],
0194               float& QTemplate);
0195 
0196   // Overload for backward compatibility
0197 
0198   bool xytemp(float xhit, float yhit, bool ydouble[BYM2], bool xdouble[BXM2], float template2d[BXM2][BYM2]);
0199 
0200   // Overload for backward compatibility with re-weighting code
0201 
0202   bool xytemp(int id,
0203               float cotalpha,
0204               float cotbeta,
0205               float xhit,
0206               float yhit,
0207               std::vector<bool>& ydouble,
0208               std::vector<bool>& xdouble,
0209               float template2d[BXM2][BYM2]);
0210 
0211   void xysigma2(float qpixel, int index, float& xysig2);
0212 
0213   // Get the interpolated Landau distribution parameters
0214 
0215   void landau_par(float lanpar[2][5]);
0216 
0217   float qavg() { return qavg_; }      //!< average cluster charge for this set of track angles
0218   float pixmax() { return pixmax_; }  //!< maximum pixel charge
0219   float qscale() { return qscale_; }  //!< charge scaling factor
0220   float s50() { return s50_; }        //!< 1/2 of the pixel threshold signal in adc units
0221   float sxymax() { return sxymax_; }  //!< max pixel signal for pixel error calculation
0222   float scalex(int i) {
0223     if (checkIllegalIndex("scalex", 3, i)) {
0224       return scalex_[i];
0225     } else {
0226       return 0.f;
0227     }
0228   }  //!< x-error scale factor in 4 charge bins
0229   float scaley(int i) {
0230     if (checkIllegalIndex("scaley", 3, i)) {
0231       return scaley_[i];
0232     } else {
0233       return 0.f;
0234     }
0235   }  //!< y-error scale factor in 4 charge bins
0236   float offsetx(int i) {
0237     if (checkIllegalIndex("offsetx", 3, i)) {
0238       return offsetx_[i];
0239     } else {
0240       return 0.f;
0241     }
0242   }  //!< x-offset in 4 charge bins
0243   float offsety(int i) {
0244     if (checkIllegalIndex("offsety", 3, i)) {
0245       return offsety_[i];
0246     } else {
0247       return 0.f;
0248     }
0249   }  //!< y-offset in 4 charge bins
0250   float fbin(int i) {
0251     if (checkIllegalIndex("fbin", 2, i)) {
0252       return fbin_[i];
0253     } else {
0254       return 0.f;
0255     }
0256   }                                           //!< Return lower bound of Qbin definition
0257   float sizex() { return clslenx_; }          //!< return x size of template cluster
0258   float sizey() { return clsleny_; }          //!< return y size of template cluster
0259   float chi2ppix() { return chi2ppix_; }      //!< average chi^2 per struck pixel
0260   float chi2scale() { return chi2scale_; }    //!< scale factor for chi^2 distribution
0261   float chi2avgone() { return chi2avgone_; }  //!< average y chi^2 for 1 pixel clusters
0262   float chi2minone() { return chi2minone_; }  //!< minimum of y chi^2 for 1 pixel clusters
0263   float mpvvav() { return mpvvav_; }          //!< most probable Q in Vavilov distribution
0264   float sigmavav() { return sigmavav_; }      //!< scale factor in Vavilov distribution
0265   float kappavav() { return kappavav_; }      //!< kappa parameter in Vavilov distribution
0266   float lorydrift() { return lorydrift_; }    //!< signed lorentz y-width (microns)
0267   float lorxdrift() { return lorxdrift_; }    //!< signed lorentz x-width (microns)
0268   float clsleny() { return clsleny_; }        //!< cluster y-size
0269   float clslenx() { return clslenx_; }        //!< cluster x-size
0270   float scaleyavg() { return scaleyavg_; }    //!< y-reco error scaling factor
0271   float scalexavg() { return scalexavg_; }    //!< x-reco error scaling factor
0272   float delyavg() {
0273     return delyavg_;
0274   }  //!< average difference between clsleny_ and cluster length [with threshold effects]
0275   float delysig() { return delysig_; }  //!< rms difference between clsleny_ and cluster length [with threshold effects]
0276   float xsize() { return xsize_; }      //!< pixel x-size (microns)
0277   float ysize() { return ysize_; }      //!< pixel y-size (microns)
0278   float zsize() { return zsize_; }      //!< pixel z-size or thickness (microns)
0279   int storesize() {
0280     return (int)thePixelTemp_.size();
0281   }  //!< return the size of the template store (the number of stored IDs
0282 
0283 private:
0284   bool checkIllegalIndex(const std::string whichMethod, int indMax, int i) {
0285 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0286     if (i < 0 || i > indMax) {
0287       throw cms::Exception("DataCorrupt")
0288           << "SiPixelTemplate2D::" << whichMethod << " called with illegal index = " << i << std::endl;
0289     }
0290 #else
0291     assert(i >= 0 && i < indMax + 1);
0292 
0293 #endif
0294     return true;
0295   }
0296 
0297   // Keep current template interpolaion parameters
0298 
0299   int id_current_;      //!< current id
0300   int index_id_;        //!< current index
0301   float cota_current_;  //!< current cot alpha
0302   float cotb_current_;  //!< current cot beta
0303   int Nyx_;             //!< number of cot(beta)-entries (columns) in template
0304   int Nxx_;             //!< number of cot(alpha)-entries (rows) in template
0305   int Dtype_;           //!< flags BPix (=0) or FPix (=1)
0306   float cotbeta0_;      //!< minimum cot(beta) covered
0307   float cotbeta1_;      //!< maximum cot(beta) covered
0308   float deltacotb_;     //!< cot(beta) bin size
0309   float cotalpha0_;     //!< minimum cot(alpha) covered
0310   float cotalpha1_;     //!< maximum cot(alpha) covered
0311   float deltacota_;     //!< cot(alpha) bin size
0312   int iy0_;             //!< index of nearest cot(beta) bin
0313   int iy1_;             //!< index of next-nearest cot(beta) bin
0314   float adcotb_;        //!< fractional pixel distance of cot(beta) from iy0_
0315   int jx0_;             //!< index of nearest cot(alpha) bin
0316   int jx1_;             //!< index of next-nearest cot(alpha) bin
0317   float adcota_;        //!< fractional pixel distance of cot(alpha) from jx0_
0318   int imin_;            //!< min y index of templated cluster
0319   int imax_;            //!< max y index of templated cluster
0320   int jmin_;            //!< min x index of templated cluster
0321   int jmax_;            //!< max x index of templated cluster
0322   bool flip_y_;         //!< flip y sign-sensitive quantities
0323   bool flip_x_;         //!< flip x sign-sensitive quantities
0324   bool success_;        //!< true if cotalpha, cotbeta are inside of the acceptance (dynamically loaded)
0325 
0326   // Keep results of last interpolation to return through member functions
0327 
0328   float qavg_;                //!< average cluster charge for this set of track angles
0329   float pixmax_;              //!< maximum pixel charge
0330   float qscale_;              //!< charge scaling factor
0331   float s50_;                 //!< 1/2 of the pixel threshold signal in adc units
0332   float sxymax_;              //!< average pixel signal for y-projection of cluster
0333   float xytemp_[BXM2][BYM2];  //!< template for xy-reconstruction
0334   float xypary0x0_[2][5];     //!< Polynomial error parameterization at ix0,iy0
0335   float xypary1x0_[2][5];     //!< Polynomial error parameterization at ix0,iy1
0336   float xypary0x1_[2][5];     //!< Polynomial error parameterization at ix1,iy0
0337   float lanpar_[2][5];        //!< Interpolated Landau parameters
0338   float chi2ppix_;            //!< average chi^2 per struck pixel
0339   float chi2scale_;           //!< scale factor for chi2 distribution
0340   float chi2avgone_;          //!< average chi^2 for 1 pixel clusters
0341   float chi2minone_;          //!< minimum of chi^2 for 1 pixel clusters
0342   float clsleny_;             //!< projected y-length of cluster
0343   float clslenx_;             //!< projected x-length of cluster
0344   float scalexavg_;           //!< average x-error scale factor
0345   float scaleyavg_;           //!< average y-error scale factor
0346   float delyavg_;             //!< average difference between clsleny_ and cluster length [with threshold effects]
0347   float delysig_;             //!< rms of difference between clsleny_ and cluster length [with threshold effects]
0348   float scalex_[4];           //!< x-error scale factor in charge bins
0349   float scaley_[4];           //!< y-error scale factor in charge bins
0350   float offsetx_[4];          //!< x-offset in charge bins
0351   float offsety_[4];          //!< y-offset in charge bins
0352   float mpvvav_;              //!< most probable Q in Vavilov distribution
0353   float sigmavav_;            //!< scale factor in Vavilov distribution
0354   float kappavav_;            //!< kappa parameter in Vavilov distribution
0355   float lorywidth_;           //!< Lorentz y-width (sign corrected for fpix frame)
0356   float lorxwidth_;           //!< Lorentz x-width
0357   float lorydrift_;           //!< Lorentz y-drift
0358   float lorxdrift_;           //!< Lorentz x-drift
0359   float xsize_;               //!< Pixel x-size
0360   float ysize_;               //!< Pixel y-size
0361   float zsize_;               //!< Pixel z-size (thickness)
0362   float fbin_[3];             //!< The QBin definitions in Q_clus/Q_avg
0363   const SiPixelTemplateEntry2D* entry00_;  // Pointer to presently interpolated point [iy,ix]
0364   const SiPixelTemplateEntry2D* entry10_;  // Pointer to presently interpolated point [iy+1,ix]
0365   const SiPixelTemplateEntry2D* entry01_;  // Pointer to presently interpolated point [iy,ix+1]
0366 
0367   // The actual template store is a std::vector container
0368   const std::vector<SiPixelTemplateStore2D>& thePixelTemp_;
0369 };
0370 
0371 #endif