|
|
|
Typedefs | |
| typedef boost::multi_array < float, 2 > | array_2d |
| typedef boost::multi_array < float, 3 > | array_3d |
Functions | |
| int | PixelTempReco2D (int id, float cotalpha, float cotbeta, float locBz, array_2d &cluster, std::vector< bool > &ydouble, std::vector< bool > &xdouble, SiPixelTemplate &templ, float &yrec, float &sigmay, float &proby, float &xrec, float &sigmax, float &probx, int &qbin, int speed, bool deadpix, std::vector< std::pair< int, int > > &zeropix, float &probQ) |
| int | PixelTempReco2D (int id, float cotalpha, float cotbeta, float locBz, array_2d &cluster, std::vector< bool > &ydouble, std::vector< bool > &xdouble, SiPixelTemplate &templ, float &yrec, float &sigmay, float &proby, float &xrec, float &sigmax, float &probx, int &qbin, int speed, float &probQ) |
| int | PixelTempReco2D (int id, float cotalpha, float cotbeta, array_2d &cluster, std::vector< bool > &ydouble, std::vector< bool > &xdouble, SiPixelTemplate &templ, float &yrec, float &sigmay, float &proby, float &xrec, float &sigmax, float &probx, int &qbin, int speed) |
| int | PixelTempReco2D (int id, float cotalpha, float cotbeta, array_2d &cluster, std::vector< bool > &ydouble, std::vector< bool > &xdouble, SiPixelTemplate &templ, float &yrec, float &sigmay, float &proby, float &xrec, float &sigmax, float &probx, int &qbin, int speed, float &probQ) |
| int | PixelTempSplit (int id, bool fpix, float cotalpha, float cotbeta, array_2d cluster, std::vector< bool > ydouble, std::vector< bool > xdouble, SiPixelTemplate &templ, float &yrec1, float &yrec2, float &sigmay, float &proby, float &xrec1, float &xrec2, float &sigmax, float &probx, int &qbin) |
| int | PixelTempSplit (int id, bool fpix, float cotalpha, float cotbeta, array_2d cluster, std::vector< bool > ydouble, std::vector< bool > xdouble, SiPixelTemplate &templ, float &yrec1, float &yrec2, float &sigmay, float &proby, float &xrec1, float &xrec2, float &sigmax, float &probx, int &qbin, bool deadpix, std::vector< std::pair< int, int > > zeropix) |
| typedef boost::multi_array< float, 2 > SiPixelTemplateReco::array_2d |
Definition at line 68 of file SiPixelTemplateReco.h.
| typedef boost::multi_array<float, 3> SiPixelTemplateReco::array_3d |
Definition at line 45 of file SiPixelTemplateSplit.h.
| int SiPixelTemplateReco::PixelTempReco2D | ( | int | id, |
| float | cotalpha, | ||
| float | cotbeta, | ||
| float | locBz, | ||
| array_2d & | clust, | ||
| std::vector< bool > & | ydouble, | ||
| std::vector< bool > & | xdouble, | ||
| SiPixelTemplate & | templ, | ||
| float & | yrec, | ||
| float & | sigmay, | ||
| float & | proby, | ||
| float & | xrec, | ||
| float & | sigmax, | ||
| float & | probx, | ||
| int & | qbin, | ||
| int | speed, | ||
| bool | deadpix, | ||
| std::vector< std::pair< int, int > > & | zeropix, | ||
| float & | probQ | ||
| ) |
Reconstruct the best estimate of the hit position for pixel clusters.
| id | - (input) identifier of the template to use |
| cotalpha | - (input) the cotangent of the alpha track angle (see CMS IN 2004/014) |
| cotbeta | - (input) the cotangent of the beta track angle (see CMS IN 2004/014) |
| locBz | - (input) the sign of the local B_z field for FPix (usually B_z<0 when cot(beta)>0 and B_z>0 when cot(beta)<0 |
| cluster | - (input) boost multi_array container of 7x21 array of pixel signals, origin of local coords (0,0) at center of pixel cluster[0][0]. Set dead pixels to small non-zero values (0.1 e). |
| ydouble | - (input) STL vector of 21 element array to flag a double-pixel |
| xdouble | - (input) STL vector of 7 element array to flag a double-pixel |
| templ | - (input) the template used in the reconstruction |
| yrec | - (output) best estimate of y-coordinate of hit in microns |
| sigmay | - (output) best estimate of uncertainty on yrec in microns |
| proby | - (output) probability describing goodness-of-fit for y-reco |
| xrec | - (output) best estimate of x-coordinate of hit in microns |
| sigmax | - (output) best estimate of uncertainty on xrec in microns |
| probx | - (output) probability describing goodness-of-fit for x-reco |
| qbin | - (output) index (0-4) describing the charge of the cluster qbin = 0 Q/Q_avg > 1.5 [few % of all hits] 1 1.5 > Q/Q_avg > 1.0 [~30% of all hits] 2 1.0 > Q/Q_avg > 0.85 [~30% of all hits] 3 0.85 > Q/Q_avg > min1 [~30% of all hits] 4 min1 > Q/Q_avg > min2 [~0.1% of all hits] 5 min2 > Q/Q_avg [~0.1% of all hits] |
| speed | - (input) switch (-2->5) trading speed vs robustness -2 totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4), calculates Q probability w/ VVIObj (better but slower) -1 totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4), calculates Q probability w/ TMath::VavilovI (poorer but faster) 0 totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4) 1 faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density 2 faster yet, searches same range as 1 but at 1/2 density 3 fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster) 4 fastest w/ Q prob, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster), calculates Q probability w/ VVIObj (better but slower) 5 fastest w/ Q prob, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster), calculates Q probability w/ TMath::VavilovI (poorer but faster) |
| deadpix | - (input) bool to indicate that there are dead pixels to be included in the analysis |
| zeropix | - (input) vector of index pairs pointing to the dead pixels |
| probQ | - (output) the Vavilov-distribution-based cluster charge probability |
Definition at line 120 of file SiPixelTemplateReco.cc.
References BHX, BHY, BXM1, BXM2, BXSIZE, BYM1, BYM2, BYM3, BYSIZE, SiPixelTemplate::chi2xavg(), SiPixelTemplate::chi2xavgone(), SiPixelTemplate::chi2xmin(), SiPixelTemplate::chi2xminone(), SiPixelTemplate::chi2yavg(), SiPixelTemplate::chi2yavgone(), SiPixelTemplate::chi2ymin(), SiPixelTemplate::chi2yminone(), delta, SiPixelTemplate::dxone(), SiPixelTemplate::dxtwo(), SiPixelTemplate::dyone(), SiPixelTemplate::dytwo(), ENDL, Exception, f, VVIObj::fcn(), Gamma, i, SiPixelTemplate::interpolate(), j, gen::k, LOGDEBUG, LOGERROR, maxpix, SiPixelTemplate::pixmax(), SiPixelTemplate::qavg(), SiPixelTemplate::qmin(), SiPixelTemplate::qscale(), SiPixelTemplate::s50(), python::multivaluedict::sort(), SiPixelTemplate::sxmax(), SiPixelTemplate::sxone(), SiPixelTemplate::sxtwo(), SiPixelTemplate::symax(), SiPixelTemplate::syone(), SiPixelTemplate::sytwo(), theVerboseLevel, TXSIZE, TYSIZE, SiPixelTemplate::vavilov_pars(), SiPixelTemplate::xavg(), SiPixelTemplate::xflcorr(), SiPixelTemplate::xrms(), SiPixelTemplate::xsigma2(), SiPixelTemplate::xsize(), SiPixelTemplate::xtemp(), SiPixelTemplate::yavg(), SiPixelTemplate::yflcorr(), SiPixelTemplate::yrms(), SiPixelTemplate::ysigma2(), SiPixelTemplate::ysize(), and SiPixelTemplate::ytemp().
Referenced by PixelCPETemplateReco::localPosition(), and PixelTempReco2D().
{
// Local variables
int i, j, k, minbin, binl, binh, binq, midpix, fypix, nypix, lypix, logypx;
int fxpix, nxpix, lxpix, logxpx, shifty, shiftx, nyzero[TYSIZE];
int nclusx, nclusy;
int deltaj, jmin, jmax, fxbin, lxbin, fybin, lybin, djy, djx;
int fypix2D, lypix2D, fxpix2D, lxpix2D;
float sythr, sxthr, rnorm, delta, sigma, sigavg, pseudopix, qscale, q50;
float ss2, ssa, sa2, ssba, saba, sba2, rat, fq, qtotal, qpixel;
float originx, originy, qfy, qly, qfx, qlx, bias, maxpix, minmax;
double chi2x, meanx, chi2y, meany, chi2ymin, chi2xmin, chi21max;
double hchi2, hndof, prvav, mpv, sigmaQ, kappa, xvav, beta2;
float ytemp[41][BYSIZE], xtemp[41][BXSIZE], ysum[BYSIZE], xsum[BXSIZE], ysort[BYSIZE], xsort[BXSIZE];
float chi2ybin[41], chi2xbin[41], ysig2[BYSIZE], xsig2[BXSIZE];
float yw2[BYSIZE], xw2[BXSIZE], ysw[BYSIZE], xsw[BXSIZE];
bool yd[BYSIZE], xd[BXSIZE], anyyd, anyxd, calc_probQ, use_VVIObj;
float ysize, xsize;
const float probmin={1.110223e-16};
const float probQmin={1.e-5};
// The minimum chi2 for a valid one pixel cluster = pseudopixel contribution only
const double mean1pix={0.100}, chi21min={0.160};
// First, interpolate the template needed to analyze this cluster
// check to see of the track direction is in the physical range of the loaded template
if(!templ.interpolate(id, cotalpha, cotbeta, locBz)) {
if (theVerboseLevel > 2) {LOGDEBUG("SiPixelTemplateReco") << "input cluster direction cot(alpha) = " << cotalpha << ", cot(beta) = " << cotbeta<< ", local B_z = " << locBz << ", template ID = " << id << ", no reconstruction performed" << ENDL;}
return 20;
}
// Make a local copy of the cluster container so that we can muck with it
array_2d cluster = clust;
// Check to see if Q probability is selected
calc_probQ = false;
use_VVIObj = false;
if(speed < 0) {
calc_probQ = true;
if(speed < -1) use_VVIObj = true;
speed = 0;
}
if(speed > 3) {
calc_probQ = true;
if(speed < 5) use_VVIObj = true;
speed = 3;
}
// Get pixel dimensions from the template (to allow multiple detectors in the future)
xsize = templ.xsize();
ysize = templ.ysize();
// Define size of pseudopixel
q50 = templ.s50();
pseudopix = 0.2*q50;
// Get charge scaling factor
qscale = templ.qscale();
// Check that the cluster container is (up to) a 7x21 matrix and matches the dimensions of the double pixel flags
if(cluster.num_dimensions() != 2) {
LOGERROR("SiPixelTemplateReco") << "input cluster container (BOOST Multiarray) has wrong number of dimensions" << ENDL;
return 3;
}
nclusx = (int)cluster.shape()[0];
nclusy = (int)cluster.shape()[1];
if(nclusx != (int)xdouble.size()) {
LOGERROR("SiPixelTemplateReco") << "input cluster container x-size is not equal to double pixel flag container size" << ENDL;
return 4;
}
if(nclusy != (int)ydouble.size()) {
LOGERROR("SiPixelTemplateReco") << "input cluster container y-size is not equal to double pixel flag container size" << ENDL;
return 5;
}
// enforce maximum size
if(nclusx > TXSIZE) {nclusx = TXSIZE;}
if(nclusy > TYSIZE) {nclusy = TYSIZE;}
// First, rescale all pixel charges
for(j=0; j<nclusx; ++j)
for(i=0; i<nclusy; ++i)
if(cluster[j][i] > 0) {cluster[j][i] *= qscale;}
// Next, sum the total charge and "decapitate" big pixels
qtotal = 0.;
minmax = templ.pixmax();
for(i=0; i<nclusy; ++i) {
maxpix = minmax;
if(ydouble[i]) {maxpix *=2.;}
for(j=0; j<nclusx; ++j) {
qtotal += cluster[j][i];
if(cluster[j][i] > maxpix) {cluster[j][i] = maxpix;}
}
}
// Do the cluster repair here
if(deadpix) {
fypix = BYM3; lypix = -1;
for(i=0; i<nclusy; ++i) {
ysum[i] = 0.; nyzero[i] = 0;
// Do preliminary cluster projection in y
for(j=0; j<nclusx; ++j) {
ysum[i] += cluster[j][i];
}
if(ysum[i] > 0.) {
// identify ends of cluster to determine what the missing charge should be
if(i < fypix) {fypix = i;}
if(i > lypix) {lypix = i;}
}
}
// Now loop over dead pixel list and "fix" everything
//First see if the cluster ends are redefined and that we have only one dead pixel per column
std::vector<std::pair<int, int> >::const_iterator zeroIter = zeropix.begin(), zeroEnd = zeropix.end();
for ( ; zeroIter != zeroEnd; ++zeroIter ) {
i = zeroIter->second;
if(i<0 || i>TYSIZE-1) {LOGERROR("SiPixelTemplateReco") << "dead pixel column y-index " << i << ", no reconstruction performed" << ENDL;
return 11;}
// count the number of dead pixels in each column
++nyzero[i];
// allow them to redefine the cluster ends
if(i < fypix) {fypix = i;}
if(i > lypix) {lypix = i;}
}
nypix = lypix-fypix+1;
// Now adjust the charge in the dead pixels to sum to 0.5*truncation value in the end columns and the truncation value in the interior columns
for (zeroIter = zeropix.begin(); zeroIter != zeroEnd; ++zeroIter ) {
i = zeroIter->second; j = zeroIter->first;
if(j<0 || j>TXSIZE-1) {LOGERROR("SiPixelTemplateReco") << "dead pixel column x-index " << j << ", no reconstruction performed" << ENDL;
return 12;}
if((i == fypix || i == lypix) && nypix > 1) {maxpix = templ.symax()/2.;} else {maxpix = templ.symax();}
if(ydouble[i]) {maxpix *=2.;}
if(nyzero[i] > 0 && nyzero[i] < 3) {qpixel = (maxpix - ysum[i])/(float)nyzero[i];} else {qpixel = 1.;}
if(qpixel < 1.) {qpixel = 1.;}
cluster[j][i] = qpixel;
// Adjust the total cluster charge to reflect the charge of the "repaired" cluster
qtotal += qpixel;
}
// End of cluster repair section
}
// Next, make y-projection of the cluster and copy the double pixel flags into a 25 element container
for(i=0; i<BYSIZE; ++i) { ysum[i] = 0.; yd[i] = false;}
k=0;
anyyd = false;
for(i=0; i<nclusy; ++i) {
for(j=0; j<nclusx; ++j) {
ysum[k] += cluster[j][i];
}
// If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
if(ydouble[i]) {
ysum[k] /= 2.;
ysum[k+1] = ysum[k];
yd[k] = true;
yd[k+1] = false;
k=k+2;
anyyd = true;
} else {
yd[k] = false;
++k;
}
if(k > BYM1) {break;}
}
// Next, make x-projection of the cluster and copy the double pixel flags into an 11 element container
for(i=0; i<BXSIZE; ++i) { xsum[i] = 0.; xd[i] = false;}
k=0;
anyxd = false;
for(j=0; j<nclusx; ++j) {
for(i=0; i<nclusy; ++i) {
xsum[k] += cluster[j][i];
}
// If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
if(xdouble[j]) {
xsum[k] /= 2.;
xsum[k+1] = xsum[k];
xd[k]=true;
xd[k+1]=false;
k=k+2;
anyxd = true;
} else {
xd[k]=false;
++k;
}
if(k > BXM1) {break;}
}
// next, identify the y-cluster ends, count total pixels, nypix, and logical pixels, logypx
fypix=-1;
nypix=0;
lypix=0;
logypx=0;
for(i=0; i<BYSIZE; ++i) {
if(ysum[i] > 0.) {
if(fypix == -1) {fypix = i;}
if(!yd[i]) {
ysort[logypx] = ysum[i];
++logypx;
}
++nypix;
lypix = i;
}
}
// dlengthy = (float)nypix - templ.clsleny();
// Make sure cluster is continuous
if((lypix-fypix+1) != nypix || nypix == 0) {
LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster doesn't agree with number of pixels above threshold" << ENDL;
if (theVerboseLevel > 2) {
LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
for(i=0; i<BYSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";}
LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE-1] << ENDL;
}
return 1;
}
// If cluster is longer than max template size, technique fails
if(nypix > TYSIZE) {
LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster is larger than maximum template size" << ENDL;
if (theVerboseLevel > 2) {
LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
for(i=0; i<BYSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";}
LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE-1] << ENDL;
}
return 6;
}
// Remember these numbers for later
fypix2D = fypix;
lypix2D = lypix;
// next, center the cluster on pixel 12 if necessary
midpix = (fypix+lypix)/2;
shifty = BHY - midpix;
if(shifty > 0) {
for(i=lypix; i>=fypix; --i) {
ysum[i+shifty] = ysum[i];
ysum[i] = 0.;
yd[i+shifty] = yd[i];
yd[i] = false;
}
} else if (shifty < 0) {
for(i=fypix; i<=lypix; ++i) {
ysum[i+shifty] = ysum[i];
ysum[i] = 0.;
yd[i+shifty] = yd[i];
yd[i] = false;
}
}
lypix +=shifty;
fypix +=shifty;
// If the cluster boundaries are OK, add pesudopixels, otherwise quit
if(fypix > 1 && fypix < BYM2) {
ysum[fypix-1] = pseudopix;
ysum[fypix-2] = pseudopix;
} else {return 8;}
if(lypix > 1 && lypix < BYM2) {
ysum[lypix+1] = pseudopix;
ysum[lypix+2] = pseudopix;
} else {return 8;}
// finally, determine if pixel[0] is a double pixel and make an origin correction if it is
if(ydouble[0]) {
originy = -0.5;
} else {
originy = 0.;
}
// next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx
fxpix=-1;
nxpix=0;
lxpix=0;
logxpx=0;
for(i=0; i<BXSIZE; ++i) {
if(xsum[i] > 0.) {
if(fxpix == -1) {fxpix = i;}
if(!xd[i]) {
xsort[logxpx] = xsum[i];
++logxpx;
}
++nxpix;
lxpix = i;
}
}
// dlengthx = (float)nxpix - templ.clslenx();
// Make sure cluster is continuous
if((lxpix-fxpix+1) != nxpix) {
LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster doesn't agree with number of pixels above threshold" << ENDL;
if (theVerboseLevel > 2) {
LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
for(i=0; i<BXSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";}
LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE-1] << ENDL;
}
return 2;
}
// If cluster is longer than max template size, technique fails
if(nxpix > TXSIZE) {
LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster is larger than maximum template size" << ENDL;
if (theVerboseLevel > 2) {
LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
for(i=0; i<BXSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";}
LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE-1] << ENDL;
}
return 7;
}
// Remember these numbers for later
fxpix2D = fxpix;
lxpix2D = lxpix;
// next, center the cluster on pixel 5 if necessary
midpix = (fxpix+lxpix)/2;
shiftx = BHX - midpix;
if(shiftx > 0) {
for(i=lxpix; i>=fxpix; --i) {
xsum[i+shiftx] = xsum[i];
xsum[i] = 0.;
xd[i+shiftx] = xd[i];
xd[i] = false;
}
} else if (shiftx < 0) {
for(i=fxpix; i<=lxpix; ++i) {
xsum[i+shiftx] = xsum[i];
xsum[i] = 0.;
xd[i+shiftx] = xd[i];
xd[i] = false;
}
}
lxpix +=shiftx;
fxpix +=shiftx;
// If the cluster boundaries are OK, add pesudopixels, otherwise quit
if(fxpix > 1 && fxpix < BXM2) {
xsum[fxpix-1] = pseudopix;
xsum[fxpix-2] = pseudopix;
} else {return 9;}
if(lxpix > 1 && lxpix < BXM2) {
xsum[lxpix+1] = pseudopix;
xsum[lxpix+2] = pseudopix;
} else {return 9;}
// finally, determine if pixel[0] is a double pixel and make an origin correction if it is
if(xdouble[0]) {
originx = -0.5;
} else {
originx = 0.;
}
// uncertainty and final corrections depend upon total charge bin
fq = qtotal/templ.qavg();
if(fq > 1.5) {
binq=0;
} else {
if(fq > 1.0) {
binq=1;
} else {
if(fq > 0.85) {
binq=2;
} else {
binq=3;
}
}
}
// Return the charge bin via the parameter list unless the charge is too small (then flag it)
qbin = binq;
if(!deadpix && qtotal < 0.95*templ.qmin()) {qbin = 5;} else {
if(!deadpix && qtotal < 0.95*templ.qmin(1)) {qbin = 4;}
}
if (theVerboseLevel > 9) {
LOGDEBUG("SiPixelTemplateReco") <<
"ID = " << id <<
" cot(alpha) = " << cotalpha << " cot(beta) = " << cotbeta <<
" nclusx = " << nclusx << " nclusy = " << nclusy << ENDL;
}
// Next, copy the y- and x-templates to local arrays
// First, decide on chi^2 min search parameters
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
if(speed < 0 || speed > 3) {
throw cms::Exception("DataCorrupt") << "SiPixelTemplateReco::PixelTempReco2D called with illegal speed = " << speed << std::endl;
}
#else
assert(speed >= 0 && speed < 4);
#endif
fybin = 2; lybin = 38; fxbin = 2; lxbin = 38; djy = 1; djx = 1;
if(speed > 0) {
fybin = 8; lybin = 32;
if(yd[fypix]) {fybin = 4; lybin = 36;}
if(lypix > fypix) {
if(yd[lypix-1]) {fybin = 4; lybin = 36;}
}
fxbin = 8; lxbin = 32;
if(xd[fxpix]) {fxbin = 4; lxbin = 36;}
if(lxpix > fxpix) {
if(xd[lxpix-1]) {fxbin = 4; lxbin = 36;}
}
}
if(speed > 1) {
djy = 2; djx = 2;
if(speed > 2) {
if(!anyyd) {djy = 4;}
if(!anyxd) {djx = 4;}
}
}
if (theVerboseLevel > 9) {
LOGDEBUG("SiPixelTemplateReco") <<
"fypix " << fypix << " lypix = " << lypix <<
" fybin = " << fybin << " lybin = " << lybin <<
" djy = " << djy << " logypx = " << logypx << ENDL;
LOGDEBUG("SiPixelTemplateReco") <<
"fxpix " << fxpix << " lxpix = " << lxpix <<
" fxbin = " << fxbin << " lxbin = " << lxbin <<
" djx = " << djx << " logxpx = " << logxpx << ENDL;
}
// Now do the copies
templ.ytemp(fybin, lybin, ytemp);
templ.xtemp(fxbin, lxbin, xtemp);
// Do the y-reconstruction first
// Apply the first-pass template algorithm to all clusters
// Modify the template if double pixels are present
if(nypix > logypx) {
i=fypix;
while(i < lypix) {
if(yd[i] && !yd[i+1]) {
for(j=fybin; j<=lybin; ++j) {
// Sum the adjacent cells and put the average signal in both
sigavg = (ytemp[j][i] + ytemp[j][i+1])/2.;
ytemp[j][i] = sigavg;
ytemp[j][i+1] = sigavg;
}
i += 2;
} else {
++i;
}
}
}
// Define the maximum signal to allow before de-weighting a pixel
sythr = 1.1*(templ.symax());
// Make sure that there will be at least two pixels that are not de-weighted
std::sort(&ysort[0], &ysort[logypx]);
if(logypx == 1) {sythr = 1.01f*ysort[0];} else {
if (ysort[1] > sythr) { sythr = 1.01f*ysort[1]; }
}
// Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
// for(i=0; i<BYSIZE; ++i) { ysig2[i] = 0.;}
templ.ysigma2(fypix, lypix, sythr, ysum, ysig2);
// Find the template bin that minimizes the Chi^2
chi2ymin = 1.e15;
for(i=fybin; i<=lybin; ++i) {
chi2ybin[i] = -1.e15f;
}
ss2 = 0.f;
for(i=fypix-2; i<=lypix+2; ++i) {
yw2[i] = 1.f/ysig2[i];
ysw[i] = ysum[i]*yw2[i];
ss2 += ysum[i]*ysw[i];
}
minbin = -1;
deltaj = djy;
jmin = fybin;
jmax = lybin;
while(deltaj > 0) {
for(j=jmin; j<=jmax; j+=deltaj) {
if(chi2ybin[j] < -100.f) {
ssa = 0.f;
sa2 = 0.f;
for(i=fypix-2; i<=lypix+2; ++i) {
ssa += ysw[i]*ytemp[j][i];
sa2 += ytemp[j][i]*ytemp[j][i]*yw2[i];
}
rat=ssa/ss2;
if(rat <= 0.f) {LOGERROR("SiPixelTemplateReco") << "illegal chi2ymin normalization (1) = " << rat << ENDL; rat = 1.;}
chi2ybin[j]=ss2-2.f*ssa/rat+sa2/(rat*rat);
}
if(chi2ybin[j] < chi2ymin) {
chi2ymin = chi2ybin[j];
minbin = j;
}
}
deltaj /= 2;
if(minbin > fybin) {jmin = minbin - deltaj;} else {jmin = fybin;}
if(minbin < lybin) {jmax = minbin + deltaj;} else {jmax = lybin;}
}
if (theVerboseLevel > 9) {
LOGDEBUG("SiPixelTemplateReco") <<
"minbin " << minbin << " chi2ymin = " << chi2ymin << ENDL;
}
// Do not apply final template pass to 1-pixel clusters (use calibrated offset)
if(logypx == 1) {
if(nypix ==1) {
delta = templ.dyone();
sigma = templ.syone();
} else {
delta = templ.dytwo();
sigma = templ.sytwo();
}
yrec = 0.5f*(fypix+lypix-2*shifty+2.f*originy)*ysize-delta;
if(sigma <= 0.) {
sigmay = 43.3f;
} else {
sigmay = sigma;
}
// Do probability calculation for one-pixel clusters
chi21max = fmax(chi21min, (double)templ.chi2yminone());
chi2ymin -=chi21max;
if(chi2ymin < 0.) {chi2ymin = 0.;}
// proby = gsl_cdf_chisq_Q(chi2ymin, mean1pix);
meany = fmax(mean1pix, (double)templ.chi2yavgone());
hchi2 = chi2ymin/2.; hndof = meany/2.;
proby = 1. - TMath::Gamma(hndof, hchi2);
} else {
// For cluster > 1 pix, make the second, interpolating pass with the templates
binl = minbin - 1;
binh = binl + 2;
if(binl < fybin) { binl = fybin;}
if(binh > lybin) { binh = lybin;}
ssa = 0.;
sa2 = 0.;
ssba = 0.;
saba = 0.;
sba2 = 0.;
for(i=fypix-2; i<=lypix+2; ++i) {
ssa += ysw[i]*ytemp[binl][i];
sa2 += ytemp[binl][i]*ytemp[binl][i]*yw2[i];
ssba += ysw[i]*(ytemp[binh][i] - ytemp[binl][i]);
saba += ytemp[binl][i]*(ytemp[binh][i] - ytemp[binl][i])*yw2[i];
sba2 += (ytemp[binh][i] - ytemp[binl][i])*(ytemp[binh][i] - ytemp[binl][i])*yw2[i];
}
// rat is the fraction of the "distance" from template a to template b
rat=(ssba*ssa-ss2*saba)/(ss2*sba2-ssba*ssba);
if(rat < 0.) {rat=0.;}
if(rat > 1.) {rat=1.0;}
rnorm = (ssa+rat*ssba)/ss2;
// Calculate the charges in the first and last pixels
qfy = ysum[fypix];
if(yd[fypix]) {qfy+=ysum[fypix+1];}
if(logypx > 1) {
qly=ysum[lypix];
if(yd[lypix-1]) {qly+=ysum[lypix-1];}
} else {
qly = qfy;
}
// Now calculate the mean bias correction and uncertainties
float qyfrac = (qfy-qly)/(qfy+qly);
bias = templ.yflcorr(binq,qyfrac)+templ.yavg(binq);
// uncertainty and final correction depend upon charge bin
yrec = (0.125f*binl+BHY-2.5f+rat*(binh-binl)*0.125f-(float)shifty+originy)*ysize - bias;
sigmay = templ.yrms(binq);
// Do goodness of fit test in y
if(rnorm <= 0.) {LOGERROR("SiPixelTemplateReco") << "illegal chi2y normalization (2) = " << rnorm << ENDL; rnorm = 1.;}
chi2y=ss2-2./rnorm*ssa-2./rnorm*rat*ssba+(sa2+2.*rat*saba+rat*rat*sba2)/(rnorm*rnorm)-templ.chi2ymin(binq);
if(chi2y < 0.0) {chi2y = 0.0;}
meany = templ.chi2yavg(binq);
if(meany < 0.01) {meany = 0.01;}
// gsl function that calculates the chi^2 tail prob for non-integral dof
// proby = gsl_cdf_chisq_Q(chi2y, meany);
// proby = ROOT::Math::chisquared_cdf_c(chi2y, meany);
hchi2 = chi2y/2.; hndof = meany/2.;
proby = 1. - TMath::Gamma(hndof, hchi2);
}
// Do the x-reconstruction next
// Apply the first-pass template algorithm to all clusters
// Modify the template if double pixels are present
if(nxpix > logxpx) {
i=fxpix;
while(i < lxpix) {
if(xd[i] && !xd[i+1]) {
for(j=fxbin; j<=lxbin; ++j) {
// Sum the adjacent cells and put the average signal in both
sigavg = (xtemp[j][i] + xtemp[j][i+1])/2.;
xtemp[j][i] = sigavg;
xtemp[j][i+1] = sigavg;
}
i += 2;
} else {
++i;
}
}
}
// Define the maximum signal to allow before de-weighting a pixel
sxthr = 1.1f*templ.sxmax();
// Make sure that there will be at least two pixels that are not de-weighted
std::sort(&xsort[0], &xsort[logxpx]);
if(logxpx == 1) {sxthr = 1.01f*xsort[0];} else {
if (xsort[1] > sxthr) { sxthr = 1.01f*xsort[1]; }
}
// Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
// for(i=0; i<BXSIZE; ++i) { xsig2[i] = 0.; }
templ.xsigma2(fxpix, lxpix, sxthr, xsum, xsig2);
// Find the template bin that minimizes the Chi^2
chi2xmin = 1.e15;
for(i=fxbin; i<=lxbin; ++i) { chi2xbin[i] = -1.e15f;}
ss2 = 0.f;
for(i=fxpix-2; i<=lxpix+2; ++i) {
xw2[i] = 1.f/xsig2[i];
xsw[i] = xsum[i]*xw2[i];
ss2 += xsum[i]*xsw[i];
}
minbin = -1;
deltaj = djx;
jmin = fxbin;
jmax = lxbin;
while(deltaj > 0) {
for(j=jmin; j<=jmax; j+=deltaj) {
if(chi2xbin[j] < -100.f) {
ssa = 0.f;
sa2 = 0.f;
for(i=fxpix-2; i<=lxpix+2; ++i) {
ssa += xsw[i]*xtemp[j][i];
sa2 += xtemp[j][i]*xtemp[j][i]*xw2[i];
}
rat=ssa/ss2;
if(rat <= 0.f) {LOGERROR("SiPixelTemplateReco") << "illegal chi2xmin normalization (1) = " << rat << ENDL; rat = 1.;}
chi2xbin[j]=ss2-2.f*ssa/rat+sa2/(rat*rat);
}
if(chi2xbin[j] < chi2xmin) {
chi2xmin = chi2xbin[j];
minbin = j;
}
}
deltaj /= 2;
if(minbin > fxbin) {jmin = minbin - deltaj;} else {jmin = fxbin;}
if(minbin < lxbin) {jmax = minbin + deltaj;} else {jmax = lxbin;}
}
if (theVerboseLevel > 9) {
LOGDEBUG("SiPixelTemplateReco") <<
"minbin " << minbin << " chi2xmin = " << chi2xmin << ENDL;
}
// Do not apply final template pass to 1-pixel clusters (use calibrated offset)
if(logxpx == 1) {
if(nxpix ==1) {
delta = templ.dxone();
sigma = templ.sxone();
} else {
delta = templ.dxtwo();
sigma = templ.sxtwo();
}
xrec = 0.5*(fxpix+lxpix-2*shiftx+2.*originx)*xsize-delta;
if(sigma <= 0.) {
sigmax = 28.9;
} else {
sigmax = sigma;
}
// Do probability calculation for one-pixel clusters
chi21max = fmax(chi21min, (double)templ.chi2xminone());
chi2xmin -=chi21max;
if(chi2xmin < 0.) {chi2xmin = 0.;}
meanx = fmax(mean1pix, (double)templ.chi2xavgone());
hchi2 = chi2xmin/2.; hndof = meanx/2.;
probx = 1. - TMath::Gamma(hndof, hchi2);
} else {
// Now make the second, interpolating pass with the templates
binl = minbin - 1;
binh = binl + 2;
if(binl < fxbin) { binl = fxbin;}
if(binh > lxbin) { binh = lxbin;}
ssa = 0.;
sa2 = 0.;
ssba = 0.;
saba = 0.;
sba2 = 0.;
for(i=fxpix-2; i<=lxpix+2; ++i) {
ssa += xsw[i]*xtemp[binl][i];
sa2 += xtemp[binl][i]*xtemp[binl][i]*xw2[i];
ssba += xsw[i]*(xtemp[binh][i] - xtemp[binl][i]);
saba += xtemp[binl][i]*(xtemp[binh][i] - xtemp[binl][i])*xw2[i];
sba2 += (xtemp[binh][i] - xtemp[binl][i])*(xtemp[binh][i] - xtemp[binl][i])*xw2[i];
}
// rat is the fraction of the "distance" from template a to template b
rat=(ssba*ssa-ss2*saba)/(ss2*sba2-ssba*ssba);
if(rat < 0.f) {rat=0.f;}
if(rat > 1.f) {rat=1.0f;}
rnorm = (ssa+rat*ssba)/ss2;
// Calculate the charges in the first and last pixels
qfx = xsum[fxpix];
if(xd[fxpix]) {qfx+=xsum[fxpix+1];}
if(logxpx > 1) {
qlx=xsum[lxpix];
if(xd[lxpix-1]) {qlx+=xsum[lxpix-1];}
} else {
qlx = qfx;
}
// Now calculate the mean bias correction and uncertainties
float qxfrac = (qfx-qlx)/(qfx+qlx);
bias = templ.xflcorr(binq,qxfrac)+templ.xavg(binq);
// uncertainty and final correction depend upon charge bin
xrec = (0.125f*binl+BHX-2.5f+rat*(binh-binl)*0.125f-(float)shiftx+originx)*xsize - bias;
sigmax = templ.xrms(binq);
// Do goodness of fit test in x
if(rnorm <= 0.) {LOGERROR("SiPixelTemplateReco") << "illegal chi2x normalization (2) = " << rnorm << ENDL; rnorm = 1.;}
chi2x=ss2-2./rnorm*ssa-2./rnorm*rat*ssba+(sa2+2.*rat*saba+rat*rat*sba2)/(rnorm*rnorm)-templ.chi2xmin(binq);
if(chi2x < 0.0) {chi2x = 0.0;}
meanx = templ.chi2xavg(binq);
if(meanx < 0.01) {meanx = 0.01;}
// gsl function that calculates the chi^2 tail prob for non-integral dof
// probx = gsl_cdf_chisq_Q(chi2x, meanx);
// probx = ROOT::Math::chisquared_cdf_c(chi2x, meanx, trx0);
hchi2 = chi2x/2.; hndof = meanx/2.;
probx = 1. - TMath::Gamma(hndof, hchi2);
}
// Don't return exact zeros for the probability
if(proby < probmin) {proby = probmin;}
if(probx < probmin) {probx = probmin;}
// Decide whether to generate a cluster charge probability
if(calc_probQ) {
// Calculate the Vavilov probability that the cluster charge is OK
templ.vavilov_pars(mpv, sigmaQ, kappa);
#ifndef SI_PIXEL_TEMPLATE_STANDALONE
if((sigmaQ <=0.) || (mpv <= 0.) || (kappa < 0.01) || (kappa > 9.9)) {
throw cms::Exception("DataCorrupt") << "SiPixelTemplateReco::Vavilov parameters mpv/sigmaQ/kappa = " << mpv << "/" << sigmaQ << "/" << kappa << std::endl;
}
#else
assert((sigmaQ > 0.) && (mpv > 0.) && (kappa > 0.01) && (kappa < 10.));
#endif
xvav = ((double)qtotal-mpv)/sigmaQ;
beta2 = 1.;
if(use_VVIObj) {
// VVIObj is a private port of CERNLIB VVIDIS
VVIObj vvidist(kappa, beta2, 1);
prvav = vvidist.fcn(xvav);
// std::cout << "vav: " << kappa << " " << xvav << " " << prvav
// << " " << TMath::VavilovI(xvav, kappa, beta2) << std::endl;
} else {
// Use faster but less accurate TMath Vavilov distribution function
prvav = TMath::VavilovI(xvav, kappa, beta2);
}
// Change to upper tail probability
// if(prvav > 0.5) prvav = 1. - prvav;
// probQ = (float)(2.*prvav);
probQ = 1. - prvav;
if(probQ < probQmin) {probQ = probQmin;}
} else {
probQ = -1;
}
return 0;
} // PixelTempReco2D
| int SiPixelTemplateReco::PixelTempReco2D | ( | int | id, |
| float | cotalpha, | ||
| float | cotbeta, | ||
| float | locBz, | ||
| array_2d & | cluster, | ||
| std::vector< bool > & | ydouble, | ||
| std::vector< bool > & | xdouble, | ||
| SiPixelTemplate & | templ, | ||
| float & | yrec, | ||
| float & | sigmay, | ||
| float & | proby, | ||
| float & | xrec, | ||
| float & | sigmax, | ||
| float & | probx, | ||
| int & | qbin, | ||
| int | speed, | ||
| float & | probQ | ||
| ) |
Reconstruct the best estimate of the hit position for pixel clusters.
| id | - (input) identifier of the template to use |
| cotalpha | - (input) the cotangent of the alpha track angle (see CMS IN 2004/014) |
| cotbeta | - (input) the cotangent of the beta track angle (see CMS IN 2004/014) |
| locBz | - (input) the sign of the local B_z field for FPix (usually B_z<0 when cot(beta)>0 and B_z>0 when cot(beta)<0 |
| cluster | - (input) boost multi_array container of 7x21 array of pixel signals, origin of local coords (0,0) at center of pixel cluster[0][0]. |
| ydouble | - (input) STL vector of 21 element array to flag a double-pixel |
| xdouble | - (input) STL vector of 7 element array to flag a double-pixel |
| templ | - (input) the template used in the reconstruction |
| yrec | - (output) best estimate of y-coordinate of hit in microns |
| sigmay | - (output) best estimate of uncertainty on yrec in microns |
| proby | - (output) probability describing goodness-of-fit for y-reco |
| xrec | - (output) best estimate of x-coordinate of hit in microns |
| sigmax | - (output) best estimate of uncertainty on xrec in microns |
| probx | - (output) probability describing goodness-of-fit for x-reco |
| qbin | - (output) index (0-4) describing the charge 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] |
| speed | - (input) switch (-1-4) trading speed vs robustness -1 totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4), calculates Q probability 0 totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4) 1 faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density 2 faster yet, searches same range as 1 but at 1/2 density 3 fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster) 4 fastest w/ Q prob, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster), calculates Q probability |
| probQ | - (output) the Vavilov-distribution-based cluster charge probability |
Definition at line 1030 of file SiPixelTemplateReco.cc.
References PixelTempReco2D().
{
// Local variables
const bool deadpix = false;
std::vector<std::pair<int, int> > zeropix;
return SiPixelTemplateReco::PixelTempReco2D(id, cotalpha, cotbeta, locBz, cluster, ydouble, xdouble, templ,
yrec, sigmay, proby, xrec, sigmax, probx, qbin, speed, deadpix, zeropix, probQ);
} // PixelTempReco2D
| int SiPixelTemplateReco::PixelTempReco2D | ( | int | id, |
| float | cotalpha, | ||
| float | cotbeta, | ||
| array_2d & | cluster, | ||
| std::vector< bool > & | ydouble, | ||
| std::vector< bool > & | xdouble, | ||
| SiPixelTemplate & | templ, | ||
| float & | yrec, | ||
| float & | sigmay, | ||
| float & | proby, | ||
| float & | xrec, | ||
| float & | sigmax, | ||
| float & | probx, | ||
| int & | qbin, | ||
| int | speed | ||
| ) |
Reconstruct the best estimate of the hit position for pixel clusters.
| id | - (input) identifier of the template to use |
| cotalpha | - (input) the cotangent of the alpha track angle (see CMS IN 2004/014) |
| cotbeta | - (input) the cotangent of the beta track angle (see CMS IN 2004/014) |
| cluster | - (input) boost multi_array container of 7x21 array of pixel signals, origin of local coords (0,0) at center of pixel cluster[0][0]. |
| ydouble | - (input) STL vector of 21 element array to flag a double-pixel |
| xdouble | - (input) STL vector of 7 element array to flag a double-pixel |
| templ | - (input) the template used in the reconstruction |
| yrec | - (output) best estimate of y-coordinate of hit in microns |
| sigmay | - (output) best estimate of uncertainty on yrec in microns |
| proby | - (output) probability describing goodness-of-fit for y-reco |
| xrec | - (output) best estimate of x-coordinate of hit in microns |
| sigmax | - (output) best estimate of uncertainty on xrec in microns |
| probx | - (output) probability describing goodness-of-fit for x-reco |
| qbin | - (output) index (0-4) describing the charge 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] |
| speed | - (input) switch (0-3) trading speed vs robustness 0 totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4) 1 faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density 2 faster yet, searches same range as 1 but at 1/2 density 3 fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster) |
Definition at line 1118 of file SiPixelTemplateReco.cc.
References PixelTempReco2D().
{
// Local variables
const bool deadpix = false;
std::vector<std::pair<int, int> > zeropix;
float locBz = -1.;
if(cotbeta < 0.) {locBz = -locBz;}
float probQ;
if(speed < 0) speed = 0;
if(speed > 3) speed = 3;
return SiPixelTemplateReco::PixelTempReco2D(id, cotalpha, cotbeta, locBz, cluster, ydouble, xdouble, templ,
yrec, sigmay, proby, xrec, sigmax, probx, qbin, speed, deadpix, zeropix, probQ);
} // PixelTempReco2D
| int SiPixelTemplateReco::PixelTempReco2D | ( | int | id, |
| float | cotalpha, | ||
| float | cotbeta, | ||
| array_2d & | cluster, | ||
| std::vector< bool > & | ydouble, | ||
| std::vector< bool > & | xdouble, | ||
| SiPixelTemplate & | templ, | ||
| float & | yrec, | ||
| float & | sigmay, | ||
| float & | proby, | ||
| float & | xrec, | ||
| float & | sigmax, | ||
| float & | probx, | ||
| int & | qbin, | ||
| int | speed, | ||
| float & | probQ | ||
| ) |
Reconstruct the best estimate of the hit position for pixel clusters.
| id | - (input) identifier of the template to use |
| cotalpha | - (input) the cotangent of the alpha track angle (see CMS IN 2004/014) |
| cotbeta | - (input) the cotangent of the beta track angle (see CMS IN 2004/014) |
| cluster | - (input) boost multi_array container of 7x21 array of pixel signals, origin of local coords (0,0) at center of pixel cluster[0][0]. |
| ydouble | - (input) STL vector of 21 element array to flag a double-pixel |
| xdouble | - (input) STL vector of 7 element array to flag a double-pixel |
| templ | - (input) the template used in the reconstruction |
| yrec | - (output) best estimate of y-coordinate of hit in microns |
| sigmay | - (output) best estimate of uncertainty on yrec in microns |
| proby | - (output) probability describing goodness-of-fit for y-reco |
| xrec | - (output) best estimate of x-coordinate of hit in microns |
| sigmax | - (output) best estimate of uncertainty on xrec in microns |
| probx | - (output) probability describing goodness-of-fit for x-reco |
| qbin | - (output) index (0-4) describing the charge 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] |
| speed | - (input) switch (-1-4) trading speed vs robustness -1 totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4), calculates Q probability 0 totally bombproof, searches the entire 41 bin range at full density (equiv to V2_4) 1 faster, searches reduced 25 bin range (no big pix) + 33 bins (big pix at ends) at full density 2 faster yet, searches same range as 1 but at 1/2 density 3 fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster) 4 fastest w/ Q prob, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster), calculates Q probability |
| probQ | - (output) the Vavilov-distribution-based cluster charge probability |
Definition at line 1074 of file SiPixelTemplateReco.cc.
References PixelTempReco2D().
{
// Local variables
const bool deadpix = false;
std::vector<std::pair<int, int> > zeropix;
float locBz = -1.;
if(cotbeta < 0.) {locBz = -locBz;}
return SiPixelTemplateReco::PixelTempReco2D(id, cotalpha, cotbeta, locBz, cluster, ydouble, xdouble, templ,
yrec, sigmay, proby, xrec, sigmax, probx, qbin, speed, deadpix, zeropix, probQ);
} // PixelTempReco2D
| int SiPixelTemplateReco::PixelTempSplit | ( | int | id, |
| bool | fpix, | ||
| float | cotalpha, | ||
| float | cotbeta, | ||
| array_2d | cluster, | ||
| std::vector< bool > | ydouble, | ||
| std::vector< bool > | xdouble, | ||
| SiPixelTemplate & | templ, | ||
| float & | yrec1, | ||
| float & | yrec2, | ||
| float & | sigmay, | ||
| float & | proby, | ||
| float & | xrec1, | ||
| float & | xrec2, | ||
| float & | sigmax, | ||
| float & | probx, | ||
| int & | qbin | ||
| ) |
Reconstruct the best estimate of the hit positions for pixel clusters.
| id | - (input) identifier of the template to use |
| fpix | - (input) logical input indicating whether to use FPix templates (true) or Barrel templates (false) |
| cotalpha | - (input) the cotangent of the alpha track angle (see CMS IN 2004/014) |
| cotbeta | - (input) the cotangent of the beta track angle (see CMS IN 2004/014) |
| cluster | - (input) boost multi_array container of 7x21 array of pixel signals, origin of local coords (0,0) at center of pixel cluster[0][0]. |
| ydouble | - (input) STL vector of 21 element array to flag a double-pixel |
| xdouble | - (input) STL vector of 7 element array to flag a double-pixel |
| templ | - (input) the template used in the reconstruction |
| yrec1 | - (output) best estimate of first y-coordinate of hit in microns |
| yrec2 | - (output) best estimate of second y-coordinate of hit in microns |
| sigmay | - (output) best estimate of uncertainty on yrec in microns |
| proby | - (output) probability describing goodness-of-fit for y-reco |
| xrec1 | - (output) best estimate of first x-coordinate of hit in microns |
| xrec2 | - (output) best estimate of second x-coordinate of hit in microns |
| sigmax | - (output) best estimate of uncertainty on xrec in microns |
| probx | - (output) probability describing goodness-of-fit for x-reco |
| qbin | - (output) index (0-4) describing the charge 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] |
Definition at line 800 of file SiPixelTemplateSplit.cc.
References PixelTempSplit().
{
// Local variables
const bool deadpix = false;
std::vector<std::pair<int, int> > zeropix;
return SiPixelTemplateReco::PixelTempSplit(id, fpix, cotalpha, cotbeta, cluster, ydouble, xdouble, templ,
yrec1, yrec2, sigmay, proby, xrec1, xrec2, sigmax, probx, qbin, deadpix, zeropix);
} // PixelTempSplit
| int SiPixelTemplateReco::PixelTempSplit | ( | int | id, |
| bool | fpix, | ||
| float | cotalpha, | ||
| float | cotbeta, | ||
| array_2d | cluster, | ||
| std::vector< bool > | ydouble, | ||
| std::vector< bool > | xdouble, | ||
| SiPixelTemplate & | templ, | ||
| float & | yrec1, | ||
| float & | yrec2, | ||
| float & | sigmay, | ||
| float & | proby, | ||
| float & | xrec1, | ||
| float & | xrec2, | ||
| float & | sigmax, | ||
| float & | probx, | ||
| int & | qbin, | ||
| bool | deadpix, | ||
| std::vector< std::pair< int, int > > | zeropix | ||
| ) |
Reconstruct the best estimate of the hit positions for pixel clusters.
| id | - (input) identifier of the template to use |
| fpix | - (input) logical input indicating whether to use FPix templates (true) or Barrel templates (false) |
| cotalpha | - (input) the cotangent of the alpha track angle (see CMS IN 2004/014) |
| cotbeta | - (input) the cotangent of the beta track angle (see CMS IN 2004/014) |
| cluster | - (input) boost multi_array container of 7x21 array of pixel signals, origin of local coords (0,0) at center of pixel cluster[0][0]. |
| ydouble | - (input) STL vector of 21 element array to flag a double-pixel |
| xdouble | - (input) STL vector of 7 element array to flag a double-pixel |
| templ | - (input) the template used in the reconstruction |
| yrec1 | - (output) best estimate of first y-coordinate of hit in microns |
| yrec2 | - (output) best estimate of second y-coordinate of hit in microns |
| sigmay | - (output) best estimate of uncertainty on yrec in microns |
| proby | - (output) probability describing goodness-of-fit for y-reco |
| xrec1 | - (output) best estimate of first x-coordinate of hit in microns |
| xrec2 | - (output) best estimate of second x-coordinate of hit in microns |
| sigmax | - (output) best estimate of uncertainty on xrec in microns |
| probx | - (output) probability describing goodness-of-fit for x-reco |
| qbin | - (output) index (0-4) describing the charge 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] |
| deadpix | - (input) bool to indicate that there are dead pixels to be included in the analysis |
| zeropix | - (input) vector of index pairs pointing to the dead pixels |
Definition at line 70 of file SiPixelTemplateSplit.cc.
References BHX, BHY, BXM1, BXM2, BXSIZE, BYM1, BYM2, BYM3, BYSIZE, SiPixelTemplate::chi2xavg(), SiPixelTemplate::chi2xavgone(), SiPixelTemplate::chi2xminone(), SiPixelTemplate::chi2yavg(), SiPixelTemplate::chi2yavgone(), SiPixelTemplate::chi2yminone(), delta, SiPixelTemplate::dxone(), SiPixelTemplate::dxtwo(), SiPixelTemplate::dyone(), SiPixelTemplate::dytwo(), ENDL, Gamma, i, SiPixelTemplate::interpolate(), j, gen::k, LOGDEBUG, LOGERROR, max(), maxpix, min, SiPixelTemplate::pixmax(), SiPixelTemplate::qavg(), SiPixelTemplate::qmin(), SiPixelTemplate::qscale(), SiPixelTemplate::s50(), python::multivaluedict::sort(), SiPixelTemplate::sxmax(), SiPixelTemplate::sxone(), SiPixelTemplate::sxtwo(), SiPixelTemplate::symax(), SiPixelTemplate::syone(), SiPixelTemplate::sytwo(), theVerboseLevel, TXSIZE, TYSIZE, SiPixelTemplate::xavgc2m(), SiPixelTemplate::xrmsc2m(), SiPixelTemplate::xsigma2(), SiPixelTemplate::xtemp3d(), SiPixelTemplate::yavgc2m(), SiPixelTemplate::yrmsc2m(), SiPixelTemplate::ysigma2(), and SiPixelTemplate::ytemp3d().
Referenced by PixelCPETemplateReco::localPosition(), and PixelTempSplit().
{
// Local variables
static int i, j, k, minbin, binq, midpix, fypix, nypix, lypix, logypx;
static int fxpix, nxpix, lxpix, logxpx, shifty, shiftx, nyzero[TYSIZE];
static int nclusx, nclusy;
static int nybin, ycbin, nxbin, xcbin, minbinj, minbink;
static float sythr, sxthr, delta, sigma, sigavg, pseudopix, qscale;
static float ss2, ssa, sa2, rat, fq, qtotal, qpixel;
static float originx, originy, bias, maxpix, minmax;
static double chi2x, meanx, chi2y, meany, chi2ymin, chi2xmin, chi21max;
static double hchi2, hndof;
static float ysum[BYSIZE], xsum[BXSIZE], ysort[BYSIZE], xsort[BXSIZE];
static float ysig2[BYSIZE], xsig2[BXSIZE];
static bool yd[BYSIZE], xd[BXSIZE], anyyd, anyxd;
const float ysize={150.}, xsize={100.}, sqrt2={2.};
const float probmin={1.110223e-16};
// const float sqrt2={1.41421356};
// The minimum chi2 for a valid one pixel cluster = pseudopixel contribution only
const double mean1pix={0.100}, chi21min={0.160};
// First, interpolate the template needed to analyze this cluster
// check to see of the track direction is in the physical range of the loaded template
if(!templ.interpolate(id, fpix, cotalpha, cotbeta)) {
LOGDEBUG("SiPixelTemplateReco") << "input cluster direction cot(alpha) = " << cotalpha << ", cot(beta) = " << cotbeta << " is not within the acceptance of fpix = "
<< fpix << ", template ID = " << id << ", no reconstruction performed" << ENDL;
return 20;
}
// Define size of pseudopixel
pseudopix = templ.s50();
// Get charge scaling factor
qscale = templ.qscale();
// Check that the cluster container is (up to) a 7x21 matrix and matches the dimensions of the double pixel flags
if(cluster.num_dimensions() != 2) {
LOGERROR("SiPixelTemplateReco") << "input cluster container (BOOST Multiarray) has wrong number of dimensions" << ENDL;
return 3;
}
nclusx = (int)cluster.shape()[0];
nclusy = (int)cluster.shape()[1];
if(nclusx != (int)xdouble.size()) {
LOGERROR("SiPixelTemplateReco") << "input cluster container x-size is not equal to double pixel flag container size" << ENDL;
return 4;
}
if(nclusy != (int)ydouble.size()) {
LOGERROR("SiPixelTemplateReco") << "input cluster container y-size is not equal to double pixel flag container size" << ENDL;
return 5;
}
// enforce maximum size
if(nclusx > TXSIZE) {nclusx = TXSIZE;}
if(nclusy > TYSIZE) {nclusy = TYSIZE;}
// First, rescale all pixel charges
for(i=0; i<nclusy; ++i) {
for(j=0; j<nclusx; ++j) {
if(cluster[j][i] > 0) {cluster[j][i] *= qscale;}
}
}
// Next, sum the total charge and "decapitate" big pixels
qtotal = 0.;
minmax = 2.*templ.pixmax();
for(i=0; i<nclusy; ++i) {
maxpix = minmax;
if(ydouble[i]) {maxpix *=2.;}
for(j=0; j<nclusx; ++j) {
qtotal += cluster[j][i];
if(cluster[j][i] > maxpix) {cluster[j][i] = maxpix;}
}
}
// Do the cluster repair here
if(deadpix) {
fypix = BYM3; lypix = -1;
for(i=0; i<nclusy; ++i) {
ysum[i] = 0.; nyzero[i] = 0;
// Do preliminary cluster projection in y
for(j=0; j<nclusx; ++j) {
ysum[i] += cluster[j][i];
}
if(ysum[i] > 0.) {
// identify ends of cluster to determine what the missing charge should be
if(i < fypix) {fypix = i;}
if(i > lypix) {lypix = i;}
}
}
// Now loop over dead pixel list and "fix" everything
//First see if the cluster ends are redefined and that we have only one dead pixel per column
std::vector<std::pair<int, int> >::const_iterator zeroIter = zeropix.begin(), zeroEnd = zeropix.end();
for ( ; zeroIter != zeroEnd; ++zeroIter ) {
i = zeroIter->second;
if(i<0 || i>TYSIZE-1) {LOGERROR("SiPixelTemplateReco") << "dead pixel column y-index " << i << ", no reconstruction performed" << ENDL;
return 11;}
// count the number of dead pixels in each column
++nyzero[i];
// allow them to redefine the cluster ends
if(i < fypix) {fypix = i;}
if(i > lypix) {lypix = i;}
}
nypix = lypix-fypix+1;
// Now adjust the charge in the dead pixels to sum to 0.5*truncation value in the end columns and the truncation value in the interior columns
for (zeroIter = zeropix.begin(); zeroIter != zeroEnd; ++zeroIter ) {
i = zeroIter->second; j = zeroIter->first;
if(j<0 || j>TXSIZE-1) {LOGERROR("SiPixelTemplateReco") << "dead pixel column x-index " << j << ", no reconstruction performed" << ENDL;
return 12;}
if((i == fypix || i == lypix) && nypix > 1) {maxpix = templ.symax()/2.;} else {maxpix = templ.symax();}
if(ydouble[i]) {maxpix *=2.;}
if(nyzero[i] > 0 && nyzero[i] < 3) {qpixel = (maxpix - ysum[i])/(float)nyzero[i];} else {qpixel = 1.;}
if(qpixel < 1.) {qpixel = 1.;}
cluster[j][i] = qpixel;
// Adjust the total cluster charge to reflect the charge of the "repaired" cluster
qtotal += qpixel;
}
// End of cluster repair section
}
// Next, make y-projection of the cluster and copy the double pixel flags into a 25 element container
for(i=0; i<BYSIZE; ++i) { ysum[i] = 0.; yd[i] = false;}
k=0;
anyyd = false;
for(i=0; i<nclusy; ++i) {
for(j=0; j<nclusx; ++j) {
ysum[k] += cluster[j][i];
}
// If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
if(ydouble[i]) {
ysum[k] /= 2.;
ysum[k+1] = ysum[k];
yd[k] = true;
yd[k+1] = false;
k=k+2;
anyyd = true;
} else {
yd[k] = false;
++k;
}
if(k > BYM1) {break;}
}
// Next, make x-projection of the cluster and copy the double pixel flags into an 11 element container
for(i=0; i<BXSIZE; ++i) { xsum[i] = 0.; xd[i] = false;}
k=0;
anyxd = false;
for(j=0; j<nclusx; ++j) {
for(i=0; i<nclusy; ++i) {
xsum[k] += cluster[j][i];
}
// If this is a double pixel, put 1/2 of the charge in 2 consective single pixels
if(xdouble[j]) {
xsum[k] /= 2.;
xsum[k+1] = xsum[k];
xd[k]=true;
xd[k+1]=false;
k=k+2;
anyxd = true;
} else {
xd[k]=false;
++k;
}
if(k > BXM1) {break;}
}
// next, identify the y-cluster ends, count total pixels, nypix, and logical pixels, logypx
fypix=-1;
nypix=0;
lypix=0;
logypx=0;
for(i=0; i<BYSIZE; ++i) {
if(ysum[i] > 0.) {
if(fypix == -1) {fypix = i;}
if(!yd[i]) {
ysort[logypx] = ysum[i];
++logypx;
}
++nypix;
lypix = i;
}
}
// Make sure cluster is continuous
if((lypix-fypix+1) != nypix || nypix == 0) {
LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster doesn't agree with number of pixels above threshold" << ENDL;
if (theVerboseLevel > 2) {
LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
for(i=0; i<BYSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";}
LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE-1] << ENDL;
}
return 1;
}
// If cluster is longer than max template size, technique fails
if(nypix > TYSIZE) {
LOGDEBUG("SiPixelTemplateReco") << "y-length of pixel cluster is larger than maximum template size" << ENDL;
if (theVerboseLevel > 2) {
LOGDEBUG("SiPixelTemplateReco") << "ysum[] = ";
for(i=0; i<BYSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << ysum[i] << ", ";}
LOGDEBUG("SiPixelTemplateReco") << ysum[BYSIZE-1] << ENDL;
}
return 6;
}
// next, center the cluster on pixel 12 if necessary
midpix = (fypix+lypix)/2;
shifty = BHY - midpix;
if(shifty > 0) {
for(i=lypix; i>=fypix; --i) {
ysum[i+shifty] = ysum[i];
ysum[i] = 0.;
yd[i+shifty] = yd[i];
yd[i] = false;
}
} else if (shifty < 0) {
for(i=fypix; i<=lypix; ++i) {
ysum[i+shifty] = ysum[i];
ysum[i] = 0.;
yd[i+shifty] = yd[i];
yd[i] = false;
}
}
lypix +=shifty;
fypix +=shifty;
// If the cluster boundaries are OK, add pesudopixels, otherwise quit
if(fypix > 1 && fypix < BYM2) {
ysum[fypix-1] = pseudopix;
ysum[fypix-2] = 0.2*pseudopix;
} else {return 8;}
if(lypix > 1 && lypix < BYM2) {
ysum[lypix+1] = pseudopix;
ysum[lypix+2] = 0.2*pseudopix;
} else {return 8;}
// finally, determine if pixel[0] is a double pixel and make an origin correction if it is
if(ydouble[0]) {
originy = -0.5;
} else {
originy = 0.;
}
// next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx
fxpix=-1;
nxpix=0;
lxpix=0;
logxpx=0;
for(i=0; i<BXSIZE; ++i) {
if(xsum[i] > 0.) {
if(fxpix == -1) {fxpix = i;}
if(!xd[i]) {
xsort[logxpx] = xsum[i];
++logxpx;
}
++nxpix;
lxpix = i;
}
}
// Make sure cluster is continuous
if((lxpix-fxpix+1) != nxpix) {
LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster doesn't agree with number of pixels above threshold" << ENDL;
if (theVerboseLevel > 2) {
LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
for(i=0; i<BXSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";}
LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE-1] << ENDL;
}
return 2;
}
// If cluster is longer than max template size, technique fails
if(nxpix > TXSIZE) {
LOGDEBUG("SiPixelTemplateReco") << "x-length of pixel cluster is larger than maximum template size" << ENDL;
if (theVerboseLevel > 2) {
LOGDEBUG("SiPixelTemplateReco") << "xsum[] = ";
for(i=0; i<BXSIZE-1; ++i) {LOGDEBUG("SiPixelTemplateReco") << xsum[i] << ", ";}
LOGDEBUG("SiPixelTemplateReco") << ysum[BXSIZE-1] << ENDL;
}
return 7;
}
// next, center the cluster on pixel 5 if necessary
midpix = (fxpix+lxpix)/2;
shiftx = BHX - midpix;
if(shiftx > 0) {
for(i=lxpix; i>=fxpix; --i) {
xsum[i+shiftx] = xsum[i];
xsum[i] = 0.;
xd[i+shiftx] = xd[i];
xd[i] = false;
}
} else if (shiftx < 0) {
for(i=fxpix; i<=lxpix; ++i) {
xsum[i+shiftx] = xsum[i];
xsum[i] = 0.;
xd[i+shiftx] = xd[i];
xd[i] = false;
}
}
lxpix +=shiftx;
fxpix +=shiftx;
// If the cluster boundaries are OK, add pesudopixels, otherwise quit
if(fxpix > 1 && fxpix <BXM2) {
xsum[fxpix-1] = pseudopix;
xsum[fxpix-2] = 0.2*pseudopix;
} else {return 9;}
if(lxpix > 1 && lxpix < BXM2) {
xsum[lxpix+1] = pseudopix;
xsum[lxpix+2] = 0.2*pseudopix;
} else {return 9;}
// finally, determine if pixel[0] is a double pixel and make an origin correction if it is
if(xdouble[0]) {
originx = -0.5;
} else {
originx = 0.;
}
// uncertainty and final corrections depend upon total charge bin
fq = qtotal/templ.qavg();
if(fq > 3.0) {
binq=0;
} else {
if(fq > 2.0) {
binq=1;
} else {
if(fq > 1.70) {
binq=2;
} else {
binq=3;
}
}
}
// Return the charge bin via the parameter list unless the charge is too small (then flag it)
qbin = binq;
if(!deadpix && qtotal < 1.9*templ.qmin()) {qbin = 5;} else {
if(!deadpix && qtotal < 1.9*templ.qmin(1)) {qbin = 4;}
}
if (theVerboseLevel > 9) {
LOGDEBUG("SiPixelTemplateReco") <<
"ID = " << id << " FPix = " << fpix <<
" cot(alpha) = " << cotalpha << " cot(beta) = " << cotbeta <<
" nclusx = " << nclusx << " nclusy = " << nclusy << ENDL;
}
// Next, generate the 3d y- and x-templates
array_3d ytemp;
templ.ytemp3d(nypix, ytemp);
array_3d xtemp;
templ.xtemp3d(nxpix, xtemp);
// Do the y-reconstruction first
// retrieve the number of y-bins
nybin = ytemp.shape()[0]; ycbin = nybin/2;
// Modify the template if double pixels are present
if(nypix > logypx) {
i=fypix;
while(i < lypix) {
if(yd[i] && !yd[i+1]) {
for(j=0; j<nybin; ++j) {
for(k=0; k<=j; ++k) {
// Sum the adjacent cells and put the average signal in both
sigavg = (ytemp[j][k][i] + ytemp[j][k][i+1])/2.;
ytemp[j][k][i] = sigavg;
ytemp[j][k][i+1] = sigavg;
}
}
i += 2;
} else {
++i;
}
}
}
// Define the maximum signal to allow before de-weighting a pixel
sythr = 2.1*(templ.symax());
// Make sure that there will be at least two pixels that are not de-weighted
std::sort(&ysort[0], &ysort[logypx]);
if(logypx == 1) {sythr = 1.01*ysort[0];} else {
if (ysort[1] > sythr) { sythr = 1.01*ysort[1]; }
}
// Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
for(i=0; i<BYSIZE; ++i) { ysig2[i] = 0.;}
templ.ysigma2(fypix, lypix, sythr, ysum, ysig2);
// Find the template bin that minimizes the Chi^2
chi2ymin = 1.e15;
minbinj = -1;
minbink = -1;
for(j=0; j<nybin; ++j) {
for(k=0; k<=j; ++k) {
ss2 = 0.;
ssa = 0.;
sa2 = 0.;
for(i=fypix-2; i<=lypix+2; ++i) {
ss2 += ysum[i]*ysum[i]/ysig2[i];
ssa += ysum[i]*ytemp[j][k][i]/ysig2[i];
sa2 += ytemp[j][k][i]*ytemp[j][k][i]/ysig2[i];
}
rat=ssa/ss2;
if(rat <= 0.) {LOGERROR("SiPixelTemplateReco") << "illegal chi2ymin normalization = " << rat << ENDL; rat = 1.;}
chi2y=ss2-2.*ssa/rat+sa2/(rat*rat);
if(chi2y < chi2ymin) {
chi2ymin = chi2y;
minbinj = j;
minbink = k;
}
}
}
if (theVerboseLevel > 9) {
LOGDEBUG("SiPixelTemplateReco") <<
"minbins " << minbinj << "," << minbink << " chi2ymin = " << chi2ymin << ENDL;
}
// Do not apply final template pass to 1-pixel clusters (use calibrated offset)
if(logypx == 1) {
if(nypix ==1) {
delta = templ.dyone();
sigma = templ.syone();
} else {
delta = templ.dytwo();
sigma = templ.sytwo();
}
yrec1 = 0.5*(fypix+lypix-2*shifty+2.*originy)*ysize-delta;
yrec2 = yrec1;
if(sigma <= 0.) {
sigmay = 43.3;
} else {
sigmay = sigma;
}
// Do probability calculation for one-pixel clusters
chi21max = fmax(chi21min, (double)templ.chi2yminone());
chi2ymin -=chi21max;
if(chi2ymin < 0.) {chi2ymin = 0.;}
// proby = gsl_cdf_chisq_Q(chi2ymin, mean1pix);
meany = fmax(mean1pix, (double)templ.chi2yavgone());
hchi2 = chi2ymin/2.; hndof = meany/2.;
proby = 1. - TMath::Gamma(hndof, hchi2);
} else {
// For cluster > 1 pix, use chi^2 minimm to recontruct the two y-positions
// at small eta, the templates won't actually work on two pixel y-clusters so just return the pixel centers
if(logypx == 2 && fabsf(cotbeta) < 0.25) {
switch(nypix) {
case 2:
// Both pixels are small
yrec1 = (fypix-shifty+originy)*ysize;
yrec2 = (lypix-shifty+originy)*ysize;
sigmay = 43.3;
break;
case 3:
// One big pixel and one small pixel
if(yd[fypix]) {
yrec1 = (fypix+0.5-shifty+originy)*ysize;
yrec2 = (lypix-shifty+originy)*ysize;
sigmay = 43.3;
} else {
yrec1 = (fypix-shifty+originy)*ysize;
yrec2 = (lypix-0.5-shifty+originy)*ysize;
sigmay = 65.;
}
break;
case 4:
// Two big pixels
yrec1 = (fypix+0.5-shifty+originy)*ysize;
yrec2 = (lypix-0.5-shifty+originy)*ysize;
sigmay = 86.6;
break;
default:
// Something is screwy ...
LOGERROR("SiPixelTemplateReco") << "weird problem: logical y-pixels = " << logypx << ", total ysize in normal pixels = " << nypix << ENDL;
return 10;
}
} else {
// uncertainty and final correction depend upon charge bin
bias = templ.yavgc2m(binq);
yrec1 = (0.125*(minbink-ycbin)+BHY-(float)shifty+originy)*ysize - bias;
yrec2 = (0.125*(minbinj-ycbin)+BHY-(float)shifty+originy)*ysize - bias;
sigmay = sqrt2*templ.yrmsc2m(binq);
}
// Do goodness of fit test in y
if(chi2ymin < 0.0) {chi2ymin = 0.0;}
meany = 2.*templ.chi2yavg(binq);
if(meany < 0.01) {meany = 0.01;}
// gsl function that calculates the chi^2 tail prob for non-integral dof
// proby = gsl_cdf_chisq_Q(chi2y, meany);
// proby = ROOT::Math::chisquared_cdf_c(chi2y, meany);
hchi2 = chi2ymin/2.; hndof = meany/2.;
proby = 1. - TMath::Gamma(hndof, hchi2);
}
// Do the x-reconstruction next
// retrieve the number of x-bins
nxbin = xtemp.shape()[0]; xcbin = nxbin/2;
// Modify the template if double pixels are present
if(nxpix > logxpx) {
i=fxpix;
while(i < lxpix) {
if(xd[i] && !xd[i+1]) {
for(j=0; j<nxbin; ++j) {
for(k=0; k<=j; ++k) {
// Sum the adjacent cells and put the average signal in both
sigavg = (xtemp[j][k][i] + xtemp[j][k][i+1])/2.;
xtemp[j][k][i] = sigavg;
xtemp[j][k][i+1] = sigavg;
}
}
i += 2;
} else {
++i;
}
}
}
// Define the maximum signal to allow before de-weighting a pixel
sxthr = 2.1*templ.sxmax();
// Make sure that there will be at least two pixels that are not de-weighted
std::sort(&xsort[0], &xsort[logxpx]);
if(logxpx == 1) {sxthr = 1.01*xsort[0];} else {
if (xsort[1] > sxthr) { sxthr = 1.01*xsort[1]; }
}
// Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
for(i=0; i<BYSIZE; ++i) { xsig2[i] = 0.; }
templ.xsigma2(fxpix, lxpix, sxthr, xsum, xsig2);
// Find the template bin that minimizes the Chi^2
chi2xmin = 1.e15;
minbinj = -1;
minbink = -1;
for(j=0; j<nxbin; ++j) {
for(k=0; k<=j; ++k) {
ss2 = 0.;
ssa = 0.;
sa2 = 0.;
for(i=fxpix-2; i<=lxpix+2; ++i) {
ss2 += xsum[i]*xsum[i]/xsig2[i];
ssa += xsum[i]*xtemp[j][k][i]/xsig2[i];
sa2 += xtemp[j][k][i]*xtemp[j][k][i]/xsig2[i];
}
rat=ssa/ss2;
if(rat <= 0.) {LOGERROR("SiPixelTemplateReco") << "illegal chi2ymin normalization = " << rat << ENDL; rat = 1.;}
chi2x=ss2-2.*ssa/rat+sa2/(rat*rat);
if(chi2x < chi2xmin) {
chi2xmin = chi2x;
minbinj = j;
minbink = k;
}
}
}
if (theVerboseLevel > 9) {
LOGDEBUG("SiPixelTemplateReco") <<
"minbin " << minbin << " chi2xmin = " << chi2xmin << ENDL;
}
// Do not apply final template pass to 1-pixel clusters (use calibrated offset)
if(logxpx == 1) {
if(nxpix ==1) {
delta = templ.dxone();
sigma = templ.sxone();
} else {
delta = templ.dxtwo();
sigma = templ.sxtwo();
}
xrec1 = 0.5*(fxpix+lxpix-2*shiftx+2.*originx)*xsize-delta;
xrec2 = xrec1;
if(sigma <= 0.) {
sigmax = 28.9;
} else {
sigmax = sigma;
}
// Do probability calculation for one-pixel clusters
chi21max = fmax(chi21min, (double)templ.chi2xminone());
chi2xmin -=chi21max;
if(chi2xmin < 0.) {chi2xmin = 0.;}
meanx = fmax(mean1pix, (double)templ.chi2xavgone());
hchi2 = chi2xmin/2.; hndof = meanx/2.;
probx = 1. - TMath::Gamma(hndof, hchi2);
} else {
// For cluster > 1 pix, use chi^2 minimm to recontruct the two x-positions
// uncertainty and final correction depend upon charge bin
bias = templ.xavgc2m(binq);
k = std::min(minbink, minbinj);
j = std::max(minbink, minbinj);
xrec1 = (0.125*(minbink-xcbin)+BHX-(float)shiftx+originx)*xsize - bias;
xrec2 = (0.125*(minbinj-xcbin)+BHX-(float)shiftx+originx)*xsize - bias;
sigmax = sqrt2*templ.xrmsc2m(binq);
// Do goodness of fit test in y
if(chi2xmin < 0.0) {chi2xmin = 0.0;}
meanx = 2.*templ.chi2xavg(binq);
if(meanx < 0.01) {meanx = 0.01;}
hchi2 = chi2xmin/2.; hndof = meanx/2.;
probx = 1. - TMath::Gamma(hndof, hchi2);
}
// Don't return exact zeros for the probability
if(proby < probmin) {proby = probmin;}
if(probx < probmin) {probx = probmin;}
return 0;
} // PixelTempSplit
1.7.3