SiPixelTemplate.h

Go to the documentation of this file.

//
//  SiPixelTemplate.h (v8.13)
//
//  Add goodness-of-fit info and spare entries to templates, version number in template header, more error checking
//  Add correction for (Q_F-Q_L)/(Q_F+Q_L) bias
//  Add cot(beta) reflection to reduce y-entries and more sophisticated x-interpolation
//  Fix small index searching bug in interpolate method
//  Change interpolation indexing to avoid complier complaining about possible un-initialized variables
//  Replace containers with static arrays in calls to ysigma2 and xsigma2
//  Add external threshold to calls to ysigma2 and xsigma2, fix parameter signal max for xsigma2
//  Return to 5 pixel spanning but adjust boundaries to use only when needed
//  Implement improved (faster) chi2min search that depends on pixel types
//  Fill template arrays in single calls to this object
//  Add qmin to the template
//  Add qscale to match charge scales
//  Small improvement to x-chisquare interpolation
//  Enlarge SiPixelTemplateStore to accommodate larger templates and increased alpha acceptance (reduce PT threshold to ~200 MeV)
//  Store x and y cluster sizes in fractional pixels to facilitate cluster splitting
//  Keep interpolated central 9 template bins in private storage and expand/shift in the getter functions (faster for speed=2/3) and easier to build 3d templates
//  Store error and bias information for the simple chi^2 min position analysis (no interpolation or Q_{FB} corrections) to use in cluster splitting
//  To save time, the gaussian centers and sigma are not interpolated right now (they aren't currently used).  They can be restored by un-commenting lines in the interpolate method.
//  Add a new method to calculate qbin for input cotbeta and cluster charge.  To be used for error estimation of merged clusters in PixelCPEGeneric.
//  Add bias info for Barrel and FPix separately in the header
//  Improve the charge estimation for larger cot(alpha) tracks
//  Change interpolate method to return false boolean if track angles are outside of range
//  Add template info and method for truncation information
//  Change to allow template sizes to be changed at compile time
//  Fix bug in track angle checking
//  Accommodate Dave's new DB pushfile which overloads the old method (file input)
//  Add CPEGeneric error information and expand qbin method to access useful info for PixelCPEGeneric
//  Fix large cot(alpha) bug in qmin interpolation
//  Add second qmin to allow a qbin=5 state
//  Use interpolated chi^2 info for one-pixel clusters
//  Separate BPix and FPix charge scales and thresholds
//  Fix DB pushfile version number checking bug.
//  Remove assert from qbin method
//  Replace asserts with exceptions in CMSSW
//  Change calling sequence to interpolate method to handle cot(beta)<0 for FPix cosmics
//  Add getter for pixelav Lorentz width estimates to qbin method
//  Add check on template size to interpolate and qbin methods
//  Add qbin population information, charge distribution information
//
//  V7.00 - Decouple BPix and FPix information into separate templates
//  Add methods to facilitate improved cluster splitting
//  Fix small charge scaling bug (affects FPix only)
//  Change y-slice used for the x-template to be closer to the actual cotalpha-cotbeta point 
//  (there is some weak breakdown of x-y factorization in the FPix after irradiation)
//
//  V8.00 - Add method to calculate a simple 2D template
//  Reorganize the interpolate method to extract header info only once per ID
//  V8.01 - Improve simple template normalization
//  V8.05 - Change qbin normalization to work better after irradiation
//  V8.10 - Add Vavilov distribution interpolation
//  V8.11 - Renormalize the x-templates for Guofan's cluster size calculation
//  V8.12 - Technical fix to qavg issue.
//  V8.13 - Fix qbin and fastsim interpolaters to avoid changing class variables
//
// Created by Morris Swartz on 10/27/06.
// Copyright 2006 __TheJohnsHopkinsUniversity__. All rights reserved.
//
//
 
// Build the template storage structure from several pieces 

#ifndef SiPixelTemplate_h
#define SiPixelTemplate_h 1

#include "SiPixelTemplateDefs.h"

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

#ifndef SI_PIXEL_TEMPLATE_STANDALONE
#include "CondFormats/SiPixelObjects/interface/SiPixelTemplateDBObject.h"
#include "FWCore/Utilities/interface/Exception.h"
#endif

typedef boost::multi_array<float, 3> array_3d;

struct SiPixelTemplateEntry { 
  int runnum;              
  float alpha;             
  float cotalpha;          
  float beta;              
  float cotbeta;           
  float costrk[3];            
  float qavg;              
  float pixmax;            
  float symax;             
  float dyone;             
  float syone;             
  float sxmax;             
  float dxone;             
  float sxone;             
  float dytwo;             
  float sytwo;             
  float dxtwo;             
  float sxtwo;             
  float qmin;              
  float qmin2;             
  float clsleny;           
  float clslenx;           
  float mpvvav;            
  float sigmavav;          
  float kappavav;          
  float ypar[2][5];        
  float ytemp[9][TYSIZE];  
  float xpar[2][5];        
  float xtemp[9][TXSIZE];  
  float yavg[4];           
  float yrms[4];           
  float ygx0[4];           
  float ygsig[4];          
  float yflpar[4][6];      
  float xavg[4];           
  float xrms[4];           
  float xgx0[4];           
  float xgsig[4];          
  float xflpar[4][6];      
  float chi2yavg[4];       
  float chi2ymin[4];       
  float chi2xavg[4];       
  float chi2xmin[4];       
  float chi2yavgone;       
  float chi2yminone;       
  float chi2xavgone;       
  float chi2xminone;       
  float yavgc2m[4];        
  float yrmsc2m[4];        
  float ygx0c2m[4];        
  float ygsigc2m[4];       
  float xavgc2m[4];        
  float xrmsc2m[4];        
  float xgx0c2m[4];        
  float xgsigc2m[4];       
  float yavggen[4];        
  float yrmsgen[4];        
  float ygx0gen[4];        
  float ygsiggen[4];       
  float xavggen[4];        
  float xrmsgen[4];        
  float xgx0gen[4];        
  float xgsiggen[4];       
  float qbfrac[3];         
  float fracyone;          
  float fracxone;          
  float fracytwo;          
  float fracxtwo;          
  float qavg_avg;          
  float qavg_spare;        
  float spare[4];
} ;




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



struct SiPixelTemplateStore { 
  SiPixelTemplateHeader head;
  SiPixelTemplateEntry enty[60];     
  SiPixelTemplateEntry entx[5][29];  
} ;


// ******************************************************************************************
// ******************************************************************************************
class SiPixelTemplate {
 public:
  SiPixelTemplate() {id_current = -1; index_id = -1; cota_current = 0.; cotb_current = 0.;} 
  bool pushfile(int filenum);     // load the private store with info from the 
                                  // file with the index (int) filenum
                                  
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
  bool pushfile(const SiPixelTemplateDBObject& dbobject);     // load the private store with info from db
#endif
  
    
// Interpolate input alpha and beta angles to produce a working template for each individual hit. 
  bool interpolate(int id, float cotalpha, float cotbeta, float locBz);
    
// overload for compatibility. 
  bool interpolate(int id, float cotalpha, float cotbeta);
  
// retreive interpolated templates. 
  void ytemp(int fybin, int lybin, float ytemplate[41][BYSIZE]);
  
  void xtemp(int fxbin, int lxbin, float xtemplate[41][BXSIZE]);
  
// new methods to build templates from two interpolated clusters (for splitting) 
  void ytemp3d(int nypix, array_3d& ytemplate);
  
  void xtemp3d(int nxpix, array_3d& xtemplate);
  
// Convert vector of projected signals into uncertainties for fitting. 
  void ysigma2(int fypix, int lypix, float sythr, float ysum[BYSIZE], float ysig2[BYSIZE]);
    
  void ysigma2(float qpixel, int index, float& ysig2);

  void xsigma2(int fxpix, int lxpix, float sxthr, float xsum[BXSIZE], float xsig2[BXSIZE]);
  
// Interpolate qfl correction in y. 
  float yflcorr(int binq, float qfly);
  
// Interpolate qfl correction in x. 
  float xflcorr(int binq, float qflx);
  
// Interpolate input beta angle to estimate the average charge. return qbin flag for input cluster charge, and estimate y/x errors and biases for the Generic Algorithm. 
  int qbin(int id, float cotalpha, float cotbeta, float locBz, float qclus, float& pixmx, float& sigmay, float& deltay, float& sigmax, float& deltax, 
           float& sy1, float& dy1, float& sy2, float& dy2, float& sx1, float& dx1, float& sx2, float& dx2, float& lorywidth, float& lorxwidth);
    
// Overload for backward compatibility. 
    int qbin(int id, float cotalpha, float cotbeta, float locBz, float qclus, float& pixmx, float& sigmay, float& deltay, float& sigmax, float& deltax, 
             float& sy1, float& dy1, float& sy2, float& dy2, float& sx1, float& dx1, float& sx2, float& dx2);
    
// Overload to keep legacy interface 
  int qbin(int id, float cotbeta, float qclus);
    
// Method to return template errors for fastsim
  void temperrors(int id, float cotalpha, float cotbeta, int qBin, float& sigmay, float& sigmax, float& sy1, float& sy2, float& sx1, float& sx2);
    
//Method to return qbin and size probabilities for fastsim
  void qbin_dist(int id, float cotalpha, float cotbeta, float qbin_frac[4], float& ny1_frac, float& ny2_frac, float& nx1_frac, float& nx2_frac);
    
//Method to calculate simple 2D templates 
  bool simpletemplate2D(float xhitp, float yhitp, std::vector<bool>& ydouble, std::vector<bool>& xdouble, float template2d[BXM2][BYM2]);
    
//Method to interpolate Vavilov distribution parameters
  void vavilov_pars(double& mpv, double& sigma, double& kappa);
    
   
  float qavg() {return pqavg;}        
  float pixmax() {return ppixmax;}         
  float qscale() {return pqscale;}         
  float s50() {return ps50;}               
  float symax() {return psymax;}             
  float dyone() {return pdyone;}             
  float syone() {return psyone;}             
  float dytwo() {return pdytwo;}             
  float sytwo() {return psytwo;}            
  float sxmax() {return psxmax;}            
  float dxone() {return pdxone;}             
  float sxone() {return psxone;}             
  float dxtwo() {return pdxtwo;}             
  float sxtwo() {return psxtwo;}             
  float qmin() {return pqmin;}               
  float qmin(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 1) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::qmin called with illegal index = " << i << std::endl;}
#else
      assert(i>=0 && i<2); 
#endif
     if(i==0){return pqmin;}else{return pqmin2;}} 
  float clsleny() {return pclsleny;}         
  float clslenx() {return pclslenx;}         
  float yratio() {return pyratio;}            
  float yxratio() {return pyxratio;}           
  float xxratio() {return pxxratio;}           
  float yavg(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::yavg called with illegal index = " << i << std::endl;}
#else
      assert(i>=0 && i<4); 
#endif
     return pyavg[i];}         
  float yrms(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::yrms called with illegal index = " << i << std::endl;}
#else
      assert(i>=0 && i<4); 
#endif
     return pyrms[i];}         
  float ygx0(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ygx0 called with illegal index = " << i << std::endl;}
#else
      assert(i>=0 && i<4); 
#endif
     return pygx0[i];}         
  float ygsig(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ygsig called with illegal index = " << i << std::endl;}
#else     
     assert(i>=0 && i<4); 
#endif
     return pygsig[i];}       
  float xavg(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xavg called with illegal index = " << i << std::endl;}
#else     
      assert(i>=0 && i<4); 
#endif
     return pxavg[i];}         
  float xrms(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xrms called with illegal index = " << i << std::endl;}
#else     
      assert(i>=0 && i<4); 
#endif
     return pxrms[i];}         
  float xgx0(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xgx0 called with illegal index = " << i << std::endl;}
#else     
      assert(i>=0 && i<4); 
#endif
     return pxgx0[i];}         
  float xgsig(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xgsig called with illegal index = " << i << std::endl;}
#else     
      assert(i>=0 && i<4); 
#endif
     return pxgsig[i];}       
  float chi2yavg(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::chi2yavg called with illegal index = " << i << std::endl;}
#else     
      assert(i>=0 && i<4); 
#endif
     return pchi2yavg[i];} 
  float chi2ymin(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::chi2ymin called with illegal index = " << i << std::endl;}
#else         
     assert(i>=0 && i<4); 
#endif
     return pchi2ymin[i];} 
  float chi2xavg(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::chi2xavg called with illegal index = " << i << std::endl;}
#else     
      assert(i>=0 && i<4); 
#endif
     return pchi2xavg[i];} 
  float chi2xmin(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::chi2xmin called with illegal index = " << i << std::endl;}
#else     
      assert(i>=0 && i<4);
#endif
     return pchi2xmin[i];} 
  float yavgc2m(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::yavgc2m called with illegal index = " << i << std::endl;}
#else     
      assert(i>=0 && i<4); 
#endif
     return pyavgc2m[i];}   
  float yrmsc2m(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::yrmsc2m called with illegal index = " << i << std::endl;}
#else         
     assert(i>=0 && i<4); 
#endif
     return pyrmsc2m[i];}   
  float ygx0c2m(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ygx0c2m called with illegal index = " << i << std::endl;}
#else     
      assert(i>=0 && i<4); 
#endif
     return pygx0c2m[i];}   
  float ygsigc2m(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ygsigc2m called with illegal index = " << i << std::endl;}
#else     
      assert(i>=0 && i<4); 
#endif
     return pygsigc2m[i];} 
  float xavgc2m(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xavgc2m called with illegal index = " << i << std::endl;}
#else     
      assert(i>=0 && i<4); 
#endif
     return pxavgc2m[i];}   
  float xrmsc2m(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xrmsc2m called with illegal index = " << i << std::endl;}
#else     
      assert(i>=0 && i<4); 
#endif
     return pxrmsc2m[i];}   
  float xgx0c2m(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xgx0cm2 called with illegal index = " << i << std::endl;}
#else     
      assert(i>=0 && i<4); 
#endif
     return pxgx0c2m[i];}   
  float xgsigc2m(int i) {
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
      if(i < 0 || i > 3) {throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xgsigc2m called with illegal index = " << i << std::endl;}
#else     
      assert(i>=0 && i<4); 
#endif
     return pxgsigc2m[i];} 
  float chi2yavgone() {return pchi2yavgone;}                        
  float chi2yminone() {return pchi2yminone;}                        
  float chi2xavgone() {return pchi2xavgone;}                        
  float chi2xminone() {return pchi2xminone;}                        
  float lorywidth() {return plorywidth;}                            
  float lorxwidth() {return plorxwidth;}                            
  float mpvvav() {return pmpvvav;}                                  
  float sigmavav() {return psigmavav;}                              
  float kappavav() {return pkappavav;}                              
  float xsize() {return pxsize;}                                    
  float ysize() {return pysize;}                                    
  float zsize() {return pzsize;}                                    
//  float yspare(int i) {assert(i>=0 && i<5); return pyspare[i];}    //!< vector of 5 spares interpolated in beta only
//  float xspare(int i) {assert(i>=0 && i<10); return pxspare[i];}    //!< vector of 10 spares interpolated in alpha and beta
  
  
 private:
  
  // Keep current template interpolaion parameters  
  
  int id_current;           
  int index_id;             
  float cota_current;       
  float cotb_current;       
  float abs_cotb;           
  bool success;             
  
  
  // Keep results of last interpolation to return through member functions
  
  float pqavg;              
  float ppixmax;            
  float pqscale;            
  float ps50;               
  float psymax;             
  float psyparmax;          
  float pdyone;             
  float psyone;             
  float pdytwo;             
  float psytwo;             
  float psxmax;             
  float psxparmax;          
  float pdxone;             
  float psxone;             
  float pdxtwo;             
  float psxtwo;             
  float pqmin;              
  float pclsleny;           
  float pclslenx;           
  float pyratio;            
  float pyparl[2][5];       
  float pyparh[2][5];       
  float pxparly0[2][5];     
  float pxparhy0[2][5];     
  float pytemp[9][BYSIZE];  
  float pyxratio;           
  float pxxratio;           
  float pxpar0[2][5];       
  float pxparl[2][5];       
  float pxparh[2][5];       
  float pxtemp[9][BXSIZE];  
  float pyavg[4];           
  float pyrms[4];           
  float pygx0[4];           
  float pygsig[4];          
  float pyflparl[4][6];     
  float pyflparh[4][6];     
  float pxavg[4];           
  float pxrms[4];           
  float pxgx0[4];           
  float pxgsig[4];           
  float pxflparll[4][6];    
  float pxflparlh[4][6];    
  float pxflparhl[4][6];    
  float pxflparhh[4][6];    
  float pchi2yavg[4];       
  float pchi2ymin[4];       
  float pchi2xavg[4];       
  float pchi2xmin[4];       
  float pyavgc2m[4];        
  float pyrmsc2m[4];        
  float pygx0c2m[4];        
  float pygsigc2m[4];       
  float pxavgc2m[4];        
  float pxrmsc2m[4];        
  float pxgx0c2m[4];        
  float pxgsigc2m[4];       
  float pchi2yavgone;       
  float pchi2yminone;       
  float pchi2xavgone;       
  float pchi2xminone;       
  float pqmin2;             
  float pmpvvav;            
  float psigmavav;          
  float pkappavav;          
  float plorywidth;         
  float plorxwidth;         
  float pxsize;             
  float pysize;             
  float pzsize;             
  float pqavg_avg;          

  
  // The actual template store is a std::vector container

  std::vector< SiPixelTemplateStore > thePixelTemp;
} ;


#endif