SiPixelRecHits/src/SiPixelTemplateReco2D.cc

0001 //

0002 //  SiPixelTemplateReco2D.cc (Version 2.90)

0003 //  Updated to work with the 2D template generation code

0004 //  Include all bells and whistles for edge clusters

0005 //  2.10 - Add y-lorentz drift to estimate starting point [for FPix]

0006 //  2.10 - Remove >1 pixel requirement

0007 //  2.20 - Fix major bug, change chi2 scan to 9x5 [YxX]

0008 //  2.30 - Allow one pixel clusters, improve cosmetics for increased style points from judges

0009 //  2.50 - Add variable cluster shifting to make the position parameter space more symmetric,

0010 //         also fix potential problems with variable size input clusters and double pixel flags

0011 //  2.55 - Fix another double pixel flag problem and a small pseudopixel problem in the edgegflagy = 3 case.

0012 //  2.60 - Modify the algorithm to return the point with the best chi2 from the starting point scan when

0013 //         the iterative procedure does not converge [eg 1 pixel clusters]

0014 //  2.70 - Change convergence criterion to require it in both planes [it was either]

0015 //  2.80 - Change 3D to 2D

0016 //  2.90 - Fix divide by zero for separate 1D convergence branch

0017 //  3.00 - Change VVIObjF so it only reads kappa

0018 //

0019 //

0020 //

0021 //  Created by Morris Swartz on 7/13/17.

0022 //  Simplification of VVIObjF by Tamas Vami

0023 //

0024
0025 #ifdef SI_PIXEL_TEMPLATE_STANDALONE
0026 #include <math.h>
0027 #endif
0028 #include <algorithm>
0029 #include <vector>
0030 #include <utility>
0031 #include <iostream>
0032 // ROOT::Math has a c++ function that does the probability calc, but only in v5.12 and later

0033 #include "TMath.h"
0034 #include "Math/DistFunc.h"
0035
0036 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0037 #include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplateReco2D.h"
0038 #include "RecoLocalTracker/SiPixelRecHits/interface/VVIObjF.h"
0039 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0040 #define LOGERROR(x) edm::LogError(x)
0041 #define LOGDEBUG(x) LogDebug(x)
0042 #define ENDL " "
0043 #include "FWCore/Utilities/interface/Exception.h"
0044 #else
0045 #include "SiPixelTemplateReco2D.h"
0046 #include "VVIObjF.h"
0047 #define LOGERROR(x) std::cout << x << ": "
0048 #define LOGDEBUG(x) std::cout << x << ": "
0049 #define ENDL std::endl
0050 #endif
0051
0052 using namespace SiPixelTemplateReco2D;
0053
0054 // *************************************************************************************************************************************

0055 //! Reconstruct the best estimate of the hit position for pixel clusters.

0056 //! \param         id - (input) identifier of the template to use

0057 //! \param   cotalpha - (input) the cotangent of the alpha track angle (see CMS IN 2004/014)

0058 //! \param    cotbeta - (input) the cotangent of the beta track angle (see CMS IN 2004/014)

0059 //! \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)

0060 //!                    for Phase 0 FPix IP-related tracks, locBz < 0 for cot(beta) > 0 and locBz > 0 for cot(beta) < 0

0061 //!                    for Phase 1 FPix IP-related tracks, see next comment

0062 //! \param locBx - (input) the sign of this quantity is used to determine whether to flip cot(alpha/beta)<0 quantities from cot(alpha/beta)>0 (FPix only)

0063 //!                    for Phase 1 FPix IP-related tracks, locBx/locBz > 0 for cot(alpha) > 0 and locBx/locBz < 0 for cot(alpha) < 0

0064 //!                    for Phase 1 FPix IP-related tracks, locBx > 0 for cot(beta) > 0 and locBx < 0 for cot(beta) < 0

0065 //! \param   edgeflagy - (input) flag to indicate the present of edges in y: 0-none (or interior gap), 1-edge at small y, 2-edge at large y,

0066 //!                                                                          3-edge at either end

0067 //! \param   edgeflagx - (input) flag to indicate the present of edges in x: 0-none, 1-edge at small x, 2-edge at large x

0068 //! \param    cluster - (input) pixel array struct with double pixel flags and bounds

0069 //!           origin of local coords (0,0) at center of pixel cluster[0][0].

0070 //! \param    templ2D - (input) the 2D template used in the reconstruction

0071 //! \param       yrec - (output) best estimate of y-coordinate of hit in microns

0072 //! \param     sigmay - (output) best estimate of uncertainty on yrec in microns

0073 //! \param       xrec - (output) best estimate of x-coordinate of hit in microns

0074 //! \param     sigmax - (output) best estimate of uncertainty on xrec in microns

0075 //! \param     probxy - (output) probability describing goodness-of-fit

0076 //! \param     probQ  - (output) probability describing upper cluster charge tail

0077 //! \param       qbin - (output) index (0-4) describing the charge of the cluster

0078 //!                       qbin = 0        Q/Q_avg > 1.5   [few % of all hits]

0079 //!                              1  1.5 > Q/Q_avg > 1.0   [~30% of all hits]

0080 //!                              2  1.0 > Q/Q_avg > 0.85  [~30% of all hits]

0081 //!                              3 0.85 > Q/Q_avg > min1  [~30% of all hits]

0082 //! \param     deltay - (output) template y-length - cluster length [when > 0, possibly missing end]

0083 // *************************************************************************************************************************************

0084 int SiPixelTemplateReco2D::PixelTempReco2D(int id,
0085                                            float cotalpha,
0086                                            float cotbeta,
0087                                            float locBz,
0088                                            float locBx,
0089                                            int edgeflagy,
0090                                            int edgeflagx,
0091                                            ClusMatrix& cluster,
0092                                            SiPixelTemplate2D& templ2D,
0093                                            float& yrec,
0094                                            float& sigmay,
0095                                            float& xrec,
0096                                            float& sigmax,
0097                                            float& probxy,
0098                                            float& probQ,
0099                                            int& qbin,
0100                                            float& deltay,
0101                                            int& npixels)
0102
0103 {
0104   // Local variables

0105
0106   float template2d[BXM2][BYM2], dpdx2d[2][BXM2][BYM2], fbin[3];
0107
0108   // fraction of truncation signal to measure the cluster ends

0109   const float fracpix = 0.45f;
0110
0111   const int nilist = 9, njlist = 5;
0112   const float ilist[nilist] = {0.f, -1.f, -0.75f, -0.5f, -0.25f, 0.25f, 0.5f, 0.75f, 1.f};
0113   const float jlist[njlist] = {0.f, -0.5f, -0.25f, 0.25f, 0.50f};
0114
0115   // Extract some relevant info from the 2D template

0116
0117   if (id >= 0) {  // if id < 0, bypass interpolation (used in calibration)

0118     if (!templ2D.interpolate(id, cotalpha, cotbeta, locBz, locBx))
0119       return 4;
0120   }
0121   float xsize = templ2D.xsize();
0122   float ysize = templ2D.ysize();
0123
0124   // Allow Qbin Q/Q_avg fractions to vary to optimize error estimation

0125
0126   for (int i = 0; i < 3; ++i) {
0127     fbin[i] = templ2D.fbin(i);
0128   }
0129
0130   float q50 = templ2D.s50();
0131   float pseudopix = 0.2f * q50;
0132   float pseudopix2 = q50 * q50;
0133
0134   // Get charge scaling factor

0135
0136   float qscale = templ2D.qscale();
0137
0138   // Check that the cluster container is (up to) a 7x21 matrix and matches the dimensions of the double pixel flags

0139
0140   int nclusx = cluster.mrow;
0141   int nclusy = (int)cluster.mcol;
0142   bool* xdouble = cluster.xdouble;
0143   bool* ydouble = cluster.ydouble;
0144
0145   // First, rescale all pixel charges and compute total charge

0146   float qtotal = 0.f;
0147   int imin = BYM2, imax = 0, jmin = BXM2, jmax = 0;
0148   for (int k = 0; k < nclusx * nclusy; ++k) {
0149     cluster.matrix[k] *= qscale;
0150     qtotal += cluster.matrix[k];
0151     int j = k / nclusy;
0152     int i = k - j * nclusy;
0153     if (cluster.matrix[k] > 0.f) {
0154       if (j < jmin) {
0155         jmin = j;
0156       }
0157       if (j > jmax) {
0158         jmax = j;
0159       }
0160       if (i < imin) {
0161         imin = i;
0162       }
0163       if (i > imax) {
0164         imax = i;
0165       }
0166     }
0167   }
0168
0169   //  Calculate the shifts needed to center the cluster in the buffer

0170
0171   int shiftx = THXP1 - (jmin + jmax) / 2;
0172   int shifty = THYP1 - (imin + imax) / 2;
0173   //  Always shift by at least one pixel

0174   if (shiftx < 1)
0175     shiftx = 1;
0176   if (shifty < 1)
0177     shifty = 1;
0178   //  Shift the min maxes too

0179   jmin += shiftx;
0180   jmax += shiftx;
0181   imin += shifty;
0182   imax += shifty;
0183
0184   // uncertainty and final corrections depend upon total charge bin

0185
0186   float fq0 = qtotal / templ2D.qavg();
0187
0188   // Next, copy the cluster and double pixel flags into a shifted workspace to

0189   // allow for significant zeros [pseudopixels] to be added

0190
0191   float clusxy[BXM2][BYM2];
0192   for (int j = 0; j < BXM2; ++j) {
0193     for (int i = 0; i < BYM2; ++i) {
0194       clusxy[j][i] = 0.f;
0195     }
0196   }
0197
0198   const unsigned int NPIXMAX = 200;
0199
0200   int indexxy[2][NPIXMAX];
0201   float pixel[NPIXMAX];
0202   float sigma2[NPIXMAX];
0203   float minmax = templ2D.pixmax();
0204   float ylow0 = 0.f, yhigh0 = 0.f, xlow0 = 0.f, xhigh0 = 0.f;
0205   int npixel = 0;
0206   float ysum[BYM2], ye[BYM2 + 1], ypos[BYM2], xpos[BXM2], xe[BXM2 + 1];
0207   bool yd[BYM2], xd[BXM2];
0208
0209   //  Copy double pixel flags to shifted arrays

0210
0211   for (int i = 0; i < BYM2; ++i) {
0212     ysum[i] = 0.f;
0213     yd[i] = false;
0214   }
0215   for (int j = 0; j < BXM2; ++j) {
0216     xd[j] = false;
0217   }
0218   for (int i = 0; i < nclusy; ++i) {
0219     if (ydouble[i]) {
0220       int iy = i + shifty;
0221       if (iy > -1 && iy < BYM2)
0222         yd[iy] = true;
0223     }
0224   }
0225   for (int j = 0; j < nclusx; ++j) {
0226     if (xdouble[j]) {
0227       int jx = j + shiftx;
0228       if (jx > -1 && jx < BXM2)
0229         xd[jx] = true;
0230     }
0231   }
0232   // Work out the positions of the rows+columns relative to the lower left edge of pixel[1,1]

0233   xe[0] = -xsize;
0234   ye[0] = -ysize;
0235   for (int i = 0; i < BYM2; ++i) {
0236     float ypitch = ysize;
0237     if (yd[i]) {
0238       ypitch += ysize;
0239     }
0240     ye[i + 1] = ye[i] + ypitch;
0241     ypos[i] = ye[i] + ypitch / 2.;
0242   }
0243   for (int j = 0; j < BXM2; ++j) {
0244     float xpitch = xsize;
0245     if (xd[j]) {
0246       xpitch += xsize;
0247     }
0248     xe[j + 1] = xe[j] + xpitch;
0249     xpos[j] = xe[j] + xpitch / 2.;
0250   }
0251   // Shift the cluster center to the central pixel of the array, truncate big signals

0252   for (int i = 0; i < nclusy; ++i) {
0253     int iy = i + shifty;
0254     float maxpix = minmax;
0255     if (ydouble[i]) {
0256       maxpix *= 2.f;
0257     }
0258     if (iy > -1 && iy < BYM2) {
0259       for (int j = 0; j < nclusx; ++j) {
0260         int jx = j + shiftx;
0261         if (jx > -1 && jx < BXM2) {
0262           if (cluster(j, i) > maxpix) {
0263             clusxy[jx][iy] = maxpix;
0264           } else {
0265             clusxy[jx][iy] = cluster(j, i);
0266           }
0267           if (clusxy[jx][iy] > 0.f) {
0268             ysum[iy] += clusxy[jx][iy];
0269             indexxy[0][npixel] = jx;
0270             indexxy[1][npixel] = iy;
0271             pixel[npixel] = clusxy[jx][iy];
0272             ++npixel;
0273           }
0274         }
0275       }
0276     }
0277   }
0278
0279   // Make sure that we find at least one pixel

0280   if (npixel < 1) {
0281 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0282     throw cms::Exception("DataCorrupt") << "PixelTemplateReco2D::number of pixels above threshold = " << npixel
0283                                         << std::endl;
0284 #else
0285     std::cout << "PixelTemplateReco2D::number of pixels above threshold = " << npixel << std::endl;
0286 #endif
0287     return 1;
0288   }
0289
0290   // Get the shifted coordinates of the cluster ends

0291   xlow0 = xe[jmin];
0292   xhigh0 = xe[jmax + 1];
0293   ylow0 = ye[imin];
0294   yhigh0 = ye[imax + 1];
0295
0296   // Next, calculate the error^2 [need to know which is the middle y pixel of the cluster]

0297
0298   int ypixoff = T2HYP1 - (imin + imax) / 2;
0299   for (int k = 0; k < npixel; ++k) {
0300     int ypixeff = ypixoff + indexxy[1][k];
0301     templ2D.xysigma2(pixel[k], ypixeff, sigma2[k]);
0302   }
0303
0304   // Next, find the half-height edges of the y-projection and identify any

0305   // missing columns to remove from the fit

0306
0307   int imisslow = -1, imisshigh = -1, jmisslow = -1, jmisshigh = -1;
0308   float ylow = -1.f, yhigh = -1.f;
0309   float hmaxpix = fracpix * templ2D.sxymax();
0310   for (int i = imin; i <= imax; ++i) {
0311     if (ysum[i] > hmaxpix && ysum[i - 1] < hmaxpix && ylow < 0.f) {
0312       ylow = ypos[i - 1] + (ypos[i] - ypos[i - 1]) * (hmaxpix - ysum[i - 1]) / (ysum[i] - ysum[i - 1]);
0313     }
0314     if (ysum[i] < q50) {
0315       if (imisslow < 0)
0316         imisslow = i;
0317       imisshigh = i;
0318     }
0319     if (ysum[i] > hmaxpix && ysum[i + 1] < hmaxpix) {
0320       yhigh = ypos[i] + (ypos[i + 1] - ypos[i]) * (ysum[i] - hmaxpix) / (ysum[i] - ysum[i + 1]);
0321     }
0322   }
0323   if (ylow < 0.f || yhigh < 0.f) {
0324     ylow = ylow0;
0325     yhigh = yhigh0;
0326   }
0327
0328   float templeny = templ2D.clsleny();
0329   deltay = templeny - (yhigh - ylow) / ysize;
0330
0331   //  x0 and y0 are best guess seeds for the fit

0332
0333   float x0 = 0.5f * (xlow0 + xhigh0) - templ2D.lorxdrift();
0334   float y0 = 0.5f * (ylow + yhigh) - templ2D.lorydrift();
0335   //   float y1 = yhigh - halfy - templ2D.lorydrift();

0336   //   printf("y0 = %f, y1 = %f \n", y0, y1);

0337   //   float y0 = 0.5f*(ylow + yhigh);

0338
0339   // If there are missing edge columns, set up missing column flags and number

0340   // of minimization passes

0341
0342   int npass = 1;
0343
0344   switch (edgeflagy) {
0345     case 0:
0346       break;
0347     case 1:
0348       imisshigh = imin - 1;
0349       imisslow = -1;
0350       break;
0351     case 2:
0352       imisshigh = -1;
0353       imisslow = imax + 1;
0354       break;
0355     case 3:
0356       imisshigh = imin - 1;
0357       imisslow = -1;
0358       npass = 2;
0359       break;
0360     default:
0361 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0362       throw cms::Exception("DataCorrupt") << "PixelTemplateReco2D::illegal edgeflagy = " << edgeflagy << std::endl;
0363 #else
0364       std::cout << "PixelTemplate:2D:illegal edgeflagy = " << edgeflagy << std::endl;
0365 #endif
0366   }
0367
0368   switch (edgeflagx) {
0369     case 0:
0370       break;
0371     case 1:
0372       jmisshigh = jmin - 1;
0373       jmisslow = -1;
0374       break;
0375     case 2:
0376       jmisshigh = -1;
0377       jmisslow = jmax + 1;
0378       break;
0379     default:
0380 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0381       throw cms::Exception("DataCorrupt") << "PixelTemplateReco2D::illegal edgeflagx = " << edgeflagx << std::endl;
0382 #else
0383       std::cout << "PixelTemplate:2D:illegal edgeflagx = " << edgeflagx << std::endl;
0384 #endif
0385   }
0386
0387   // Define quantities to be saved for each pass

0388
0389   float chi2min[2], xerr2[2], yerr2[2];
0390   float x2D0[2], y2D0[2], qtfrac0[2];
0391   int ipass, tpixel;
0392   //int niter0[2];

0393
0394   for (ipass = 0; ipass < npass; ++ipass) {
0395     if (ipass == 1) {
0396       // Now try again if both edges are possible

0397
0398       imisshigh = -1;
0399       imisslow = imax + 1;
0400       // erase pseudo pixels from previous pass

0401       for (int k = npixel; k < tpixel; ++k) {
0402         int j = indexxy[0][k];
0403         int i = indexxy[1][k];
0404         clusxy[j][i] = 0.f;
0405       }
0406     }
0407
0408     // Next, add pseudo pixels around the periphery of the cluster

0409
0410     tpixel = npixel;
0411     for (int k = 0; k < npixel; ++k) {
0412       int j = indexxy[0][k];
0413       int i = indexxy[1][k];
0414       if ((j - 1) != jmisshigh) {
0415         if (clusxy[j - 1][i] < pseudopix) {
0416           indexxy[0][tpixel] = j - 1;
0417           indexxy[1][tpixel] = i;
0418           clusxy[j - 1][i] = pseudopix;
0419           pixel[tpixel] = pseudopix;
0420           sigma2[tpixel] = pseudopix2;
0421           ++tpixel;
0422         }
0423       }
0424       if ((j + 1) != jmisslow) {
0425         if (clusxy[j + 1][i] < pseudopix) {
0426           indexxy[0][tpixel] = j + 1;
0427           indexxy[1][tpixel] = i;
0428           clusxy[j + 1][i] = pseudopix;
0429           pixel[tpixel] = pseudopix;
0430           sigma2[tpixel] = pseudopix2;
0431           ++tpixel;
0432         }
0433       }
0434       // Don't add them if this is a dead column

0435       if ((i + 1) != imisslow) {
0436         if ((j - 1) != jmisshigh) {
0437           if (clusxy[j - 1][i + 1] < pseudopix) {
0438             indexxy[0][tpixel] = j - 1;
0439             indexxy[1][tpixel] = i + 1;
0440             clusxy[j - 1][i + 1] = pseudopix;
0441             pixel[tpixel] = pseudopix;
0442             sigma2[tpixel] = pseudopix2;
0443             ++tpixel;
0444           }
0445         }
0446         if (clusxy[j][i + 1] < pseudopix) {
0447           indexxy[0][tpixel] = j;
0448           indexxy[1][tpixel] = i + 1;
0449           clusxy[j][i + 1] = pseudopix;
0450           pixel[tpixel] = pseudopix;
0451           sigma2[tpixel] = pseudopix2;
0452           ++tpixel;
0453         }
0454         if ((j + 1) != jmisslow) {
0455           if (clusxy[j + 1][i + 1] < pseudopix) {
0456             indexxy[0][tpixel] = j + 1;
0457             indexxy[1][tpixel] = i + 1;
0458             clusxy[j + 1][i + 1] = pseudopix;
0459             pixel[tpixel] = pseudopix;
0460             sigma2[tpixel] = pseudopix2;
0461             ++tpixel;
0462           }
0463         }
0464       }
0465       // Don't add them if this is a dead column

0466       if ((i - 1) != imisshigh) {
0467         if ((j - 1) != jmisshigh) {
0468           if (clusxy[j - 1][i - 1] < pseudopix) {
0469             indexxy[0][tpixel] = j - 1;
0470             indexxy[1][tpixel] = i - 1;
0471             clusxy[j - 1][i - 1] = pseudopix;
0472             pixel[tpixel] = pseudopix;
0473             sigma2[tpixel] = pseudopix2;
0474             ++tpixel;
0475           }
0476         }
0477         if (clusxy[j][i - 1] < pseudopix) {
0478           indexxy[0][tpixel] = j;
0479           indexxy[1][tpixel] = i - 1;
0480           clusxy[j][i - 1] = pseudopix;
0481           pixel[tpixel] = pseudopix;
0482           sigma2[tpixel] = pseudopix2;
0483           ++tpixel;
0484         }
0485         if ((j + 1) != jmisslow) {
0486           if (clusxy[j + 1][i - 1] < pseudopix) {
0487             indexxy[0][tpixel] = j + 1;
0488             indexxy[1][tpixel] = i - 1;
0489             clusxy[j + 1][i - 1] = pseudopix;
0490             pixel[tpixel] = pseudopix;
0491             sigma2[tpixel] = pseudopix2;
0492             ++tpixel;
0493           }
0494         }
0495       }
0496     }
0497
0498     // Calculate chi2 over a grid of several seeds and choose the smallest

0499
0500     chi2min[ipass] = 1000000.f;
0501     float chi2, qtemplate, qactive, qtfrac = 0.f, x2D = 0.f, y2D = 0.f;
0502     //  Scale the y search grid for long clusters [longer than 7 pixels]

0503     float ygridscale = 0.271 * cotbeta;
0504     if (ygridscale < 1.f)
0505       ygridscale = 1.f;
0506     for (int is = 0; is < nilist; ++is) {
0507       for (int js = 0; js < njlist; ++js) {
0508         float xtry = x0 + jlist[js] * xsize;
0509         float ytry = y0 + ilist[is] * ygridscale * ysize;
0510         chi2 = 0.f;
0511         qactive = 0.f;
0512         for (int j = 0; j < BXM2; ++j) {
0513           for (int i = 0; i < BYM2; ++i) {
0514             template2d[j][i] = 0.f;
0515           }
0516         }
0517         templ2D.xytemp(xtry, ytry, yd, xd, template2d, false, dpdx2d, qtemplate);
0518         for (int k = 0; k < tpixel; ++k) {
0519           int jpix = indexxy[0][k];
0520           int ipix = indexxy[1][k];
0521           float qtpixel = template2d[jpix][ipix];
0522           float pt = pixel[k] - qtpixel;
0523           chi2 += pt * pt / sigma2[k];
0524           if (k < npixel) {
0525             qactive += qtpixel;
0526           }
0527         }
0528         if (chi2 < chi2min[ipass]) {
0529           chi2min[ipass] = chi2;
0530           x2D = xtry;
0531           y2D = ytry;
0532           qtfrac = qactive / qtemplate;
0533         }
0534       }
0535     }
0536
0537     int niter = 0;
0538     float xstep = 1.0f, ystep = 1.0f;
0539     float minv11 = 1000.f, minv12 = 1000.f, minv22 = 1000.f;
0540     chi2 = chi2min[ipass];
0541     while (chi2 <= chi2min[ipass] && niter < 15 && (niter < 2 || (std::abs(xstep) > 0.2 || std::abs(ystep) > 0.2))) {
0542       // Remember the present parameters

0543       x2D0[ipass] = x2D;
0544       y2D0[ipass] = y2D;
0545       qtfrac0[ipass] = qtfrac;
0546       xerr2[ipass] = minv11;
0547       yerr2[ipass] = minv22;
0548       chi2min[ipass] = chi2;
0549       //niter0[ipass] = niter;

0550
0551       // Calculate the initial template which also allows the error calculation for the struck pixels

0552
0553       for (int j = 0; j < BXM2; ++j) {
0554         for (int i = 0; i < BYM2; ++i) {
0555           template2d[j][i] = 0.f;
0556         }
0557       }
0558       templ2D.xytemp(x2D, y2D, yd, xd, template2d, true, dpdx2d, qtemplate);
0559
0560       float sumptdt1 = 0., sumptdt2 = 0.;
0561       float sumdtdt11 = 0., sumdtdt12 = 0., sumdtdt22 = 0.;
0562       chi2 = 0.f;
0563       qactive = 0.f;
0564       // Loop over all pixels and calculate matrices

0565
0566       for (int k = 0; k < tpixel; ++k) {
0567         int jpix = indexxy[0][k];
0568         int ipix = indexxy[1][k];
0569         float qtpixel = template2d[jpix][ipix];
0570         float pt = pixel[k] - qtpixel;
0571         float dtdx = dpdx2d[0][jpix][ipix];
0572         float dtdy = dpdx2d[1][jpix][ipix];
0573         chi2 += pt * pt / sigma2[k];
0574         sumptdt1 += pt * dtdx / sigma2[k];
0575         sumptdt2 += pt * dtdy / sigma2[k];
0576         sumdtdt11 += dtdx * dtdx / sigma2[k];
0577         sumdtdt12 += dtdx * dtdy / sigma2[k];
0578         sumdtdt22 += dtdy * dtdy / sigma2[k];
0579         if (k < npixel) {
0580           qactive += qtpixel;
0581         }
0582       }
0583
0584       // invert the parameter covariance matrix

0585
0586       float D = sumdtdt11 * sumdtdt22 - sumdtdt12 * sumdtdt12;
0587
0588       // If the matrix is non-singular invert and solve

0589
0590       if (std::abs(D) > 1.e-3) {
0591         minv11 = sumdtdt22 / D;
0592         minv12 = -sumdtdt12 / D;
0593         minv22 = sumdtdt11 / D;
0594
0595         qtfrac = qactive / qtemplate;
0596
0597         //Calculate the step size

0598
0599         xstep = minv11 * sumptdt1 + minv12 * sumptdt2;
0600         ystep = minv12 * sumptdt1 + minv22 * sumptdt2;
0601
0602       } else {
0603         //  Assume alternately that ystep = 0 and then xstep = 0

0604         if (sumdtdt11 > 0.0001f) {
0605           xstep = sumptdt1 / sumdtdt11;
0606         } else {
0607           xstep = 0.f;
0608         }
0609         if (sumdtdt22 > 0.0001f) {
0610           ystep = sumptdt2 / sumdtdt22;
0611         } else {
0612           ystep = 0.f;
0613         }
0614       }
0615       xstep *= 0.9f;
0616       ystep *= 0.9f;
0617       if (std::abs(xstep) > 2. * xsize || std::abs(ystep) > 2. * ysize)
0618         break;
0619       x2D += xstep;
0620       y2D += ystep;
0621       ++niter;
0622     }
0623   }
0624
0625   ipass = 0;
0626
0627   if (npass > 1) {
0628     // two passes, take smaller chisqared

0629     if (chi2min[1] < chi2min[0]) {
0630       ipass = 1;
0631     }
0632   }
0633
0634   // Correct the charge ratio for missing pixels

0635
0636   float fq;
0637   if (qtfrac0[ipass] < 0.10f || qtfrac0[ipass] > 1.f) {
0638     qtfrac0[ipass] = 1.f;
0639   }
0640   fq = fq0 / qtfrac0[ipass];
0641
0642   //   printf("qtfrac0 = %f \n", qtfrac0);

0643
0644   if (fq > fbin[0]) {
0645     qbin = 0;
0646   } else {
0647     if (fq > fbin[1]) {
0648       qbin = 1;
0649     } else {
0650       if (fq > fbin[2]) {
0651         qbin = 2;
0652       } else {
0653         qbin = 3;
0654       }
0655     }
0656   }
0657
0658   // Get charge related quantities

0659
0660   float scalex = templ2D.scalex(qbin);
0661   float scaley = templ2D.scaley(qbin);
0662   float offsetx = templ2D.offsetx(qbin);
0663   float offsety = templ2D.offsety(qbin);
0664
0665   // This 2D code has the origin (0,0) at the lower left edge of the input cluster

0666   // That is now pixel [shiftx,shifty] and the template reco convention is the middle

0667   // of that pixel, so we need to correct

0668
0669   xrec = x2D0[ipass] - xpos[shiftx] - offsetx;
0670   yrec = y2D0[ipass] - ypos[shifty] - offsety;
0671   if (xerr2[ipass] > 0.f) {
0672     sigmax = scalex * sqrt(xerr2[ipass]);
0673     if (sigmax < 3.f)
0674       sigmax = 3.f;
0675   } else {
0676     sigmax = 10000.f;
0677   }
0678   if (yerr2[ipass] > 0.f) {
0679     sigmay = scaley * sqrt(yerr2[ipass]);
0680     if (sigmay < 3.f)
0681       sigmay = 3.f;
0682   } else {
0683     sigmay = 10000.f;
0684   }
0685   if (id >= 0) {  //if id < 0 bypass interpolation (used in calibration)

0686     // Do chi2 probability calculation

0687     double meanxy = (double)(npixel * templ2D.chi2ppix());
0688     double chi2scale = (double)templ2D.chi2scale();
0689     if (meanxy < 0.01) {
0690       meanxy = 0.01;
0691     }
0692     double hndof = meanxy / 2.f;
0693     double hchi2 = chi2scale * chi2min[ipass] / 2.f;
0694     probxy = (float)(1. - TMath::Gamma(hndof, hchi2));
0695     // Do the charge probability

0696     float mpv = templ2D.mpvvav();
0697     float sigmaQ = templ2D.sigmavav();
0698     float kappa = templ2D.kappavav();
0699     float xvav = (qtotal / qtfrac0[ipass] - mpv) / sigmaQ;
0700     //  VVIObj is a private port of CERNLIB VVIDIS

0701     VVIObjF vvidist(kappa);
0702     float prvav = vvidist.fcn(xvav);
0703     probQ = 1.f - prvav;
0704   } else {
0705     probxy = chi2min[ipass];
0706     npixels = npixel;
0707     probQ = 0.f;
0708   }
0709
0710   return 0;
0711 }  // PixelTempReco2D

0712
0713 int SiPixelTemplateReco2D::PixelTempReco2D(int id,
0714                                            float cotalpha,
0715                                            float cotbeta,
0716                                            float locBz,
0717                                            float locBx,
0718                                            int edgeflagy,
0719                                            int edgeflagx,
0720                                            ClusMatrix& cluster,
0721                                            SiPixelTemplate2D& templ2D,
0722                                            float& yrec,
0723                                            float& sigmay,
0724                                            float& xrec,
0725                                            float& sigmax,
0726                                            float& probxy,
0727                                            float& probQ,
0728                                            int& qbin,
0729                                            float& deltay)
0730
0731 {
0732   // Local variables

0733   int npixels;
0734   return SiPixelTemplateReco2D::PixelTempReco2D(id,
0735                                                 cotalpha,
0736                                                 cotbeta,
0737                                                 locBz,
0738                                                 locBx,
0739                                                 edgeflagy,
0740                                                 edgeflagx,
0741                                                 cluster,
0742                                                 templ2D,
0743                                                 yrec,
0744                                                 sigmay,
0745                                                 xrec,
0746                                                 sigmax,
0747                                                 probxy,
0748                                                 probQ,
0749                                                 qbin,
0750                                                 deltay,
0751                                                 npixels);
0752 }  // PixelTempReco2D