SiPixelTemplate2D.h

Go to the documentation of this file.

//
//  SiPixelTemplate2D.h (v2.65)
//
//  Full 2-D templates for cluster splitting, simulated cluster reweighting, and improved cluster probability
//
// Created by Morris Swartz on 12/01/09.
// V1.01 - fix qavg_ filling
// V1.02 - Add locBz to test if FPix use is out of range
// V1.03 - Fix edge checking on final template to increase template size and to properly truncate cluster
// v2.00 - Major changes to accommodate 2D reconstruction
// v2.10 - Change chi2 and error scaling information to work with partially reconstructed clusters
// v2.20 - Add cluster charge probability information, side loading for template generation
// v2.21 - Double derivative interval [improves fit convergence]
// v2.25 - Resize template store to accommodate FPix Templates
// v2.30 - Fix bug found by P. Shuetze that compromises sqlite file loading
// v2.35 - Add directory path selection to the ascii pushfile method
// v2.50 - Change template storage to dynamically allocated 2D arrays of SiPixelTemplateEntry2D structs
// v2.51 - Ensure that the derivative arrays are correctly zeroed between calls
// v2.52 - Improve cosmetics for increased style points from judges
// v2.60 - Fix FPix multiframe lookup problem [takes +-cotalpha and +-cotbeta]
// v2.61a - Code 2.60 fix correctly
// v2.65 - change double pixel flags to work with new shifted reco code definition
//

// Build the template storage structure from several pieces

#ifndef SiPixelTemplate2D_h
#define SiPixelTemplate2D_h 1

#include <vector>
#include <cassert>
#include "boost/multi_array.hpp"

#ifndef SI_PIXEL_TEMPLATE_STANDALONE
#include "CondFormats/SiPixelObjects/interface/SiPixel2DTemplateDBObject.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "CondFormats/SiPixelTransient/interface/SiPixelTemplateDefs.h"
#else
#include "SiPixelTemplateDefs.h"
#endif

struct SiPixelTemplateEntry2D {  
  int runnum;                    
  float cotalpha;                
  float cotbeta;                 
  float costrk[3];               
  float qavg;                    
  float pixmax;                  
  float sxymax;                  
  int iymin;                     
  int iymax;                     
  int jxmin;                     
  int jxmax;                     
  float xypar[2][5];             
  float lanpar[2][5];            
  short int xytemp[7][7][T2YSIZE][T2XSIZE];  
  float chi2ppix;                            
  float chi2scale;                           
  float chi2avgone;                          
  float chi2minone;                          
  float clsleny;                             
  float clslenx;                             
  float mpvvav;      
  float sigmavav;    
  float kappavav;    
  float scalexavg;   
  float scaleyavg;   
  float delyavg;     
  float delysig;     
  float scalex[4];   
  float scaley[4];   
  float offsetx[4];  
  float offsety[4];  
  float spare[3];
};

struct SiPixelTemplateHeader2D {  
  int ID;                         
  int NTy;                        
  int NTyx;                       
  int NTxx;                       
  int Dtype;                      
  float qscale;                   
  float lorywidth;                
  float lorxwidth;                
  float lorybias;                 
  float lorxbias;                 
  float Vbias;                    
  float temperature;              
  float fluence;                  
  float s50;                      
  float ss50;                     
  char title[80];                 
  int templ_version;              
  float Bfield;                   
  float fbin[3];                  
  float xsize;                    
  float ysize;                    
  float zsize;                    
};

struct SiPixelTemplateStore2D {  

  void resize(int ny, int nx) {
    entry.resize(ny);
    store.resize(nx * ny);
    int off = 0;
    for (int i = 0; i < ny; ++i) {
      entry[i] = store.data() + off;
      off += nx;
    }
    assert(nx * ny == off);
  }

  SiPixelTemplateHeader2D head;  

  std::vector<SiPixelTemplateEntry2D*> entry;

  std::vector<SiPixelTemplateEntry2D> store;
};

// ******************************************************************************************
// ******************************************************************************************
class SiPixelTemplate2D {
public:
  SiPixelTemplate2D(const std::vector<SiPixelTemplateStore2D>& thePixelTemp) : thePixelTemp_(thePixelTemp) {
    id_current_ = -1;
    index_id_ = -1;
    cota_current_ = 0.;
    cotb_current_ = 0.;
  }  

  // load the private store with info from the
  // file with the index (int) filenum ${dir}template_summary_zp${filenum}.out
#ifdef SI_PIXEL_TEMPLATE_STANDALONE
  static bool pushfile(int filenum, std::vector<SiPixelTemplateStore2D>& pixelTemp, std::string dir = "");

  // For calibrations only: load precalculated values -- no interpolation.
  void sideload(SiPixelTemplateEntry2D* entry,
                int iDtype,
                float locBx,
                float locBz,
                float lorwdy,
                float lorwdx,
                float q50,
                float fbin[3],
                float xsize,
                float ysize,
                float zsize);

#else
  static bool pushfile(int filenum,
                       std::vector<SiPixelTemplateStore2D>& pixelTemp,
                       std::string dir = "CalibTracker/SiPixelESProducers/data/");

  // Load from the DB (the default in CMSSW):
  static bool pushfile(const SiPixel2DTemplateDBObject& dbobject, std::vector<SiPixelTemplateStore2D>& pixelTemp);

#endif

  //  Initialize things before interpolating
  bool getid(int id);

  bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx);

  // Interpolate input alpha and beta angles to produce a working template for each individual hit.

  // Works with Phase 0+1
  bool xytemp(float xhit,
              float yhit,
              bool ydouble[BYM2],
              bool xdouble[BXM2],
              float template2d[BXM2][BYM2],
              bool dervatives,
              float dpdx2d[2][BXM2][BYM2],
              float& QTemplate);

  // Overload for backward compatibility

  bool xytemp(float xhit, float yhit, bool ydouble[BYM2], bool xdouble[BXM2], float template2d[BXM2][BYM2]);

  // Overload for backward compatibility with re-weighting code

  bool xytemp(int id,
              float cotalpha,
              float cotbeta,
              float xhit,
              float yhit,
              std::vector<bool>& ydouble,
              std::vector<bool>& xdouble,
              float template2d[BXM2][BYM2]);

  void xysigma2(float qpixel, int index, float& xysig2);

  // Get the interpolated Landau distribution parameters

  void landau_par(float lanpar[2][5]);

  float qavg() { return qavg_; }      
  float pixmax() { return pixmax_; }  
  float qscale() { return qscale_; }  
  float s50() { return s50_; }        
  float sxymax() { return sxymax_; }  
  float scalex(int i) {
    if (checkIllegalIndex("scalex", 3, i)) {
      return scalex_[i];
    } else {
      return 0.f;
    }
  }  
  float scaley(int i) {
    if (checkIllegalIndex("scaley", 3, i)) {
      return scaley_[i];
    } else {
      return 0.f;
    }
  }  
  float offsetx(int i) {
    if (checkIllegalIndex("offsetx", 3, i)) {
      return offsetx_[i];
    } else {
      return 0.f;
    }
  }  
  float offsety(int i) {
    if (checkIllegalIndex("offsety", 3, i)) {
      return offsety_[i];
    } else {
      return 0.f;
    }
  }  
  float fbin(int i) {
    if (checkIllegalIndex("fbin", 2, i)) {
      return fbin_[i];
    } else {
      return 0.f;
    }
  }                                           
  float sizex() { return clslenx_; }          
  float sizey() { return clsleny_; }          
  float chi2ppix() { return chi2ppix_; }      
  float chi2scale() { return chi2scale_; }    
  float chi2avgone() { return chi2avgone_; }  
  float chi2minone() { return chi2minone_; }  
  float mpvvav() { return mpvvav_; }          
  float sigmavav() { return sigmavav_; }      
  float kappavav() { return kappavav_; }      
  float lorydrift() { return lorydrift_; }    
  float lorxdrift() { return lorxdrift_; }    
  float clsleny() { return clsleny_; }        
  float clslenx() { return clslenx_; }        
  float scaleyavg() { return scaleyavg_; }    
  float scalexavg() { return scalexavg_; }    
  float delyavg() {
    return delyavg_;
  }  
  float delysig() { return delysig_; }  
  float xsize() { return xsize_; }      
  float ysize() { return ysize_; }      
  float zsize() { return zsize_; }      
  int storesize() {
    return (int)thePixelTemp_.size();
  }  

private:
  bool checkIllegalIndex(const std::string whichMethod, int indMax, int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
    if (i < 0 || i > indMax) {
      throw cms::Exception("DataCorrupt")
          << "SiPixelTemplate2D::" << whichMethod << " called with illegal index = " << i << std::endl;
    }
#else
    assert(i >= 0 && i < indMax + 1);

#endif
    return true;
  }

  // Keep current template interpolaion parameters

  int id_current_;      
  int index_id_;        
  float cota_current_;  
  float cotb_current_;  
  int Nyx_;             
  int Nxx_;             
  int Dtype_;           
  float cotbeta0_;      
  float cotbeta1_;      
  float deltacotb_;     
  float cotalpha0_;     
  float cotalpha1_;     
  float deltacota_;     
  int iy0_;             
  int iy1_;             
  float adcotb_;        
  int jx0_;             
  int jx1_;             
  float adcota_;        
  int imin_;            
  int imax_;            
  int jmin_;            
  int jmax_;            
  bool flip_y_;         
  bool flip_x_;         
  bool success_;        

  // Keep results of last interpolation to return through member functions

  float qavg_;                
  float pixmax_;              
  float qscale_;              
  float s50_;                 
  float sxymax_;              
  float xytemp_[BXM2][BYM2];  
  float xypary0x0_[2][5];     
  float xypary1x0_[2][5];     
  float xypary0x1_[2][5];     
  float lanpar_[2][5];        
  float chi2ppix_;            
  float chi2scale_;           
  float chi2avgone_;          
  float chi2minone_;          
  float clsleny_;             
  float clslenx_;             
  float scalexavg_;           
  float scaleyavg_;           
  float delyavg_;             
  float delysig_;             
  float scalex_[4];           
  float scaley_[4];           
  float offsetx_[4];          
  float offsety_[4];          
  float mpvvav_;              
  float sigmavav_;            
  float kappavav_;            
  float lorywidth_;           
  float lorxwidth_;           
  float lorydrift_;           
  float lorxdrift_;           
  float xsize_;               
  float ysize_;               
  float zsize_;               
  float fbin_[3];             
  const SiPixelTemplateEntry2D* entry00_;  // Pointer to presently interpolated point [iy,ix]
  const SiPixelTemplateEntry2D* entry10_;  // Pointer to presently interpolated point [iy+1,ix]
  const SiPixelTemplateEntry2D* entry01_;  // Pointer to presently interpolated point [iy,ix+1]

  // The actual template store is a std::vector container
  const std::vector<SiPixelTemplateStore2D>& thePixelTemp_;
};

#endif