CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Functions
SiStripTemplateReco Namespace Reference

Functions

int StripTempReco1D (int id, float cotalpha, float cotbeta, float locBy, std::vector< float > &cluster, SiStripTemplate &templ, float &xrec, float &sigmax, float &probx, int &qbin, int speed, float &probQ)
 

Function Documentation

int SiStripTemplateReco::StripTempReco1D ( int  id,
float  cotalpha,
float  cotbeta,
float  locBy,
std::vector< float > &  cluster,
SiStripTemplate templ,
float &  xrec,
float &  sigmax,
float &  probx,
int &  qbin,
int  speed,
float &  probQ 
)

Reconstruct the best estimate of the hit position for strip clusters, includes autoswitching to barycenter when that technique is more accurate.

Parameters
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)
locBy- (input) the sign of the local B_y field to specify the Lorentz drift direction
cluster- (input) boost multi_array container array of 13 pixel signals, origin of local coords (0,0) at center of pixel cluster[0].
templ- (input) the template used in the reconstruction
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)
probQ- (output) the Vavilov-distribution-based cluster charge probability

Definition at line 84 of file SiStripTemplateReco.cc.

References assert(), BSHX, BSXM2, BSXSIZE, SiStripTemplate::chi2xavg(), SiStripTemplate::chi2xavgone(), SiStripTemplate::chi2xmin(), SiStripTemplate::chi2xminone(), SiStripTemplate::cxtemp(), delta, SiStripTemplate::dxone(), ENDL, Exception, f, sistripvvi::VVIObj::fcn(), Gamma, i, SiStripTemplate::interpolate(), j, kappa, LOGDEBUG, LOGERROR, SiStripTemplate::lorxwidth(), maxpix, SiStripTemplate::qavg(), SiStripTemplate::qmin(), SiStripTemplate::qscale(), SiStripTemplate::s50(), SiStripTemplate::sxmax(), SiStripTemplate::sxone(), theVerboseLevel, TSXSIZE, SiStripTemplate::vavilov_pars(), SiStripTemplate::xavg(), SiStripTemplate::xavgbcn(), SiStripTemplate::xflcorr(), SiStripTemplate::xrms(), SiStripTemplate::xrmsbcn(), SiStripTemplate::xsigma2(), xsize, SiStripTemplate::xsize(), and SiStripTemplate::xtemp().

Referenced by StripCPEfromTemplate::localParameters().

88 {
89  // Local variables
90  int i, j, minbin, binl, binh, binq, midpix;
91  int fxpix, nxpix, lxpix, logxpx, shiftx, ftpix, ltpix;
92  int nclusx;
93  int deltaj, jmin, jmax, fxbin, lxbin, djx;
94  float sxthr, rnorm, delta, sigma, pseudopix, qscale, q50, q100;
95  float ss2, ssa, sa2, ssba, saba, sba2, rat, fq, qtotal, barycenter, sigmaxbcn;
96  float originx, qfx, qlx, bias, biasbcn, maxpix;
97  double chi2x, meanx, chi2xmin, chi21max;
98  double hchi2, hndof, prvav, mpv, sigmaQ, kappa, xvav, beta2;
99  float xtemp[41][BSXSIZE], xsum[BSXSIZE];
100  float chi2xbin[41], xsig2[BSXSIZE];
101  float xw2[BSXSIZE], xsw[BSXSIZE];
102  bool calc_probQ, use_VVIObj;
103  float xsize;
104  const float probmin={1.110223e-16};
105  const float probQmin={1.e-5};
106 
107 // The minimum chi2 for a valid one pixel cluster = pseudopixel contribution only
108 
109  const double mean1pix={0.100}, chi21min={0.160};
110 
111 // First, interpolate the template needed to analyze this cluster
112 // check to see of the track direction is in the physical range of the loaded template
113 
114  if(!templ.interpolate(id, cotalpha, cotbeta, locBy)) {
115  if (theVerboseLevel > 2) {LOGDEBUG("SiStripTemplateReco") << "input cluster direction cot(alpha) = " << cotalpha << ", cot(beta) = " << cotbeta << ", local B_y = " << locBy << ", template ID = " << id << ", no reconstruction performed" << ENDL;}
116  return 20;
117  }
118 
119 // Check to see if Q probability is selected
120 
121  calc_probQ = false;
122  use_VVIObj = false;
123  if(speed < 0) {
124  calc_probQ = true;
125  if(speed < -1) use_VVIObj = true;
126  speed = 0;
127  }
128 
129  if(speed > 3) {
130  calc_probQ = true;
131  if(speed < 5) use_VVIObj = true;
132  speed = 3;
133  }
134 
135 // Get pixel dimensions from the template (to allow multiple detectors in the future)
136 
137  xsize = templ.xsize();
138 
139 // Define size of pseudopixel
140 
141  q50 = templ.s50();
142  q100 = 2.f * q50;
143  pseudopix = q50;
144 
145 // Get charge scaling factor
146 
147  qscale = templ.qscale();
148 
149 // enforce maximum size
150 
151  nclusx = (int)cluster.size();
152 
153  if(nclusx > TSXSIZE) {nclusx = TSXSIZE;}
154 
155 // First, rescale all strip charges, sum them and trunate the strip charges
156 
157  qtotal = 0.;
158  for(i=0; i<BSXSIZE; ++i) {xsum[i] = 0.f;}
159  maxpix = templ.sxmax();
160  barycenter = 0.f;
161  for(j=0; j<nclusx; ++j) {
162  xsum[j] = qscale*cluster[j];
163  qtotal += xsum[j];
164  barycenter += j*xsize*xsum[j];
165  if(xsum[j] > maxpix) {xsum[j] = maxpix;}
166  }
167 
168  barycenter = barycenter/qtotal - 0.5f*templ.lorxwidth();
169 
170 // next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx
171 
172  fxpix = -1;
173  ftpix = -1;
174  nxpix=0;
175  lxpix=0;
176  ltpix=0;
177  logxpx=0;
178  for(i=0; i<BSXSIZE; ++i) {
179  if(xsum[i] > 0.f) {
180  if(fxpix == -1) {fxpix = i;}
181  ++logxpx;
182  ++nxpix;
183  lxpix = i;
184  if(xsum[i] > q100) {
185  if(ftpix == -1) {ftpix = i;}
186  ltpix = i;
187  }
188  }
189  }
190 
191 
192 // dlengthx = (float)nxpix - templ.clslenx();
193 
194 // Make sure cluster is continuous
195 
196  if((lxpix-fxpix+1) != nxpix) {
197 
198  LOGDEBUG("SiStripTemplateReco") << "x-length of pixel cluster doesn't agree with number of pixels above threshold" << ENDL;
199  if (theVerboseLevel > 2) {
200  LOGDEBUG("SiStripTemplateReco") << "xsum[] = ";
201  for(i=0; i<BSXSIZE-1; ++i) {LOGDEBUG("SiStripTemplateReco") << xsum[i] << ", ";}
202  LOGDEBUG("SiStripTemplateReco") << ENDL;
203  }
204 
205  return 2;
206  }
207 
208 // If cluster is longer than max template size, technique fails
209 
210  if(nxpix > TSXSIZE) {
211 
212  LOGDEBUG("SiStripTemplateReco") << "x-length of pixel cluster is larger than maximum template size" << ENDL;
213  if (theVerboseLevel > 2) {
214  LOGDEBUG("SiStripTemplateReco") << "xsum[] = ";
215  for(i=0; i<BSXSIZE-1; ++i) {LOGDEBUG("SiStripTemplateReco") << xsum[i] << ", ";}
216  LOGDEBUG("SiStripTemplateReco") << ENDL;
217  }
218 
219  return 7;
220  }
221 
222 // next, center the cluster on template center if necessary
223 
224  midpix = (ftpix+ltpix)/2;
225  shiftx = templ.cxtemp() - midpix;
226 
227  if(shiftx > 0) {
228  for(i=lxpix; i>=fxpix; --i) {
229  xsum[i+shiftx] = xsum[i];
230  xsum[i] = 0.f;
231  }
232  } else if (shiftx < 0) {
233  for(i=fxpix; i<=lxpix; ++i) {
234  xsum[i+shiftx] = xsum[i];
235  xsum[i] = 0.f;
236  }
237  }
238  lxpix +=shiftx;
239  fxpix +=shiftx;
240 
241 // If the cluster boundaries are OK, add pesudopixels, otherwise quit
242 
243  if(fxpix > 1 && fxpix < BSXM2) {
244  xsum[fxpix-1] = pseudopix;
245  xsum[fxpix-2] = 0.2f*pseudopix;
246  } else {return 9;}
247  if(lxpix > 1 && lxpix < BSXM2) {
248  xsum[lxpix+1] = pseudopix;
249  xsum[lxpix+2] = 0.2f*pseudopix;
250  } else {return 9;}
251 
252 // finally, determine if pixel[0] is a double pixel and make an origin correction if it is
253 
254  originx = 0.f;
255 
256 // uncertainty and final corrections depend upon total charge bin
257 
258  fq = qtotal/templ.qavg();
259  if(fq > 1.5f) {
260  binq=0;
261  } else {
262  if(fq > 1.0f) {
263  binq=1;
264  } else {
265  if(fq > 0.85f) {
266  binq=2;
267  } else {
268  binq=3;
269  }
270  }
271  }
272 
273 // Return the charge bin via the parameter list unless the charge is too small (then flag it)
274 
275  qbin = binq;
276  if(qtotal < 0.95f*templ.qmin()) {qbin = 5;} else {
277  if(qtotal < 0.95f*templ.qmin(1)) {qbin = 4;}
278  }
279  if (theVerboseLevel > 9) {
280  LOGDEBUG("SiStripTemplateReco") <<
281  "ID = " << id <<
282  " cot(alpha) = " << cotalpha << " cot(beta) = " << cotbeta <<
283  " nclusx = " << nclusx << ENDL;
284  }
285 
286 
287 // Next, copy the y- and x-templates to local arrays
288 
289 // First, decide on chi^2 min search parameters
290 
291 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
292  if(speed < 0 || speed > 3) {
293  throw cms::Exception("DataCorrupt") << "SiStripTemplateReco::StripTempReco2D called with illegal speed = " << speed << std::endl;
294  }
295 #else
296  assert(speed >= 0 && speed < 4);
297 #endif
298  fxbin = 2; lxbin = 38; djx = 1;
299  if(speed > 0) {
300  fxbin = 8; lxbin = 32;
301  }
302 
303  if(speed > 1) {
304  djx = 2;
305  if(speed > 2) {
306  djx = 4;
307  }
308  }
309 
310  if (theVerboseLevel > 9) {
311  LOGDEBUG("SiStripTemplateReco") <<
312  "fxpix " << fxpix << " lxpix = " << lxpix <<
313  " fxbin = " << fxbin << " lxbin = " << lxbin <<
314  " djx = " << djx << " logxpx = " << logxpx << ENDL;
315  }
316 
317 // Now do the copies
318 
319  templ.xtemp(fxbin, lxbin, xtemp);
320 
321 // Do the x-reconstruction next
322 
323 // Apply the first-pass template algorithm to all clusters
324 
325 // Modify the template if double pixels are present
326 
327 // Define the maximum signal to allow before de-weighting a pixel
328 
329  sxthr = 1.1f*maxpix;
330 
331 // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
332 
333 // for(i=0; i<BSXSIZE; ++i) { xsig2[i] = 0.; }
334  templ.xsigma2(fxpix, lxpix, sxthr, xsum, xsig2);
335 
336 // Find the template bin that minimizes the Chi^2
337 
338  chi2xmin = 1.e15;
339  for(i=fxbin; i<=lxbin; ++i) { chi2xbin[i] = -1.e15f;}
340  ss2 = 0.f;
341  for(i=fxpix-2; i<=lxpix+2; ++i) {
342  xw2[i] = 1.f/xsig2[i];
343  xsw[i] = xsum[i]*xw2[i];
344  ss2 += xsum[i]*xsw[i];
345  }
346  minbin = -1;
347  deltaj = djx;
348  jmin = fxbin;
349  jmax = lxbin;
350  while(deltaj > 0) {
351  for(j=jmin; j<=jmax; j+=deltaj) {
352  if(chi2xbin[j] < -100.f) {
353  ssa = 0.f;
354  sa2 = 0.f;
355  for(i=fxpix-2; i<=lxpix+2; ++i) {
356  ssa += xsw[i]*xtemp[j][i];
357  sa2 += xtemp[j][i]*xtemp[j][i]*xw2[i];
358  }
359  rat=ssa/ss2;
360  if(rat <= 0.f) {LOGERROR("SiStripTemplateReco") << "illegal chi2xmin normalization (1) = " << rat << ENDL; rat = 1.;}
361  chi2xbin[j]=ss2-2.*ssa/rat+sa2/(rat*rat);
362  }
363  if(chi2xbin[j] < chi2xmin) {
364  chi2xmin = chi2xbin[j];
365  minbin = j;
366  }
367  }
368  deltaj /= 2;
369  if(minbin > fxbin) {jmin = minbin - deltaj;} else {jmin = fxbin;}
370  if(minbin < lxbin) {jmax = minbin + deltaj;} else {jmax = lxbin;}
371  }
372 
373  if (theVerboseLevel > 9) {
374  LOGDEBUG("SiStripTemplateReco") <<
375  "minbin " << minbin << " chi2xmin = " << chi2xmin << ENDL;
376  }
377 
378 // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
379 
380  if(nxpix == 1) {
381 
382  delta = templ.dxone();
383  sigma = templ.sxone();
384  xrec = 0.5f*(fxpix+lxpix-2*shiftx+2.f*originx)*xsize-delta;
385  if(sigma <= 0.f) {
386  sigmax = 28.9f;
387  } else {
388  sigmax = sigma;
389  }
390 
391 // Do probability calculation for one-pixel clusters
392 
393  chi21max = fmax(chi21min, (double)templ.chi2xminone());
394  chi2xmin -=chi21max;
395  if(chi2xmin < 0.) {chi2xmin = 0.;}
396  meanx = fmax(mean1pix, (double)templ.chi2xavgone());
397  hchi2 = chi2xmin/2.; hndof = meanx/2.;
398  probx = 1. - TMath::Gamma(hndof, hchi2);
399 
400  } else {
401 
402 // Now make the second, interpolating pass with the templates
403 
404  binl = minbin - 1;
405  binh = binl + 2;
406  if(binl < fxbin) { binl = fxbin;}
407  if(binh > lxbin) { binh = lxbin;}
408  ssa = 0.;
409  sa2 = 0.;
410  ssba = 0.;
411  saba = 0.;
412  sba2 = 0.;
413  for(i=fxpix-2; i<=lxpix+2; ++i) {
414  ssa += xsw[i]*xtemp[binl][i];
415  sa2 += xtemp[binl][i]*xtemp[binl][i]*xw2[i];
416  ssba += xsw[i]*(xtemp[binh][i] - xtemp[binl][i]);
417  saba += xtemp[binl][i]*(xtemp[binh][i] - xtemp[binl][i])*xw2[i];
418  sba2 += (xtemp[binh][i] - xtemp[binl][i])*(xtemp[binh][i] - xtemp[binl][i])*xw2[i];
419  }
420 
421 // rat is the fraction of the "distance" from template a to template b
422 
423  rat=(ssba*ssa-ss2*saba)/(ss2*sba2-ssba*ssba);
424  if(rat < 0.f) {rat=0.f;}
425  if(rat > 1.f) {rat=1.0f;}
426  rnorm = (ssa+rat*ssba)/ss2;
427 
428 // Calculate the charges in the first and last pixels
429 
430  qfx = xsum[fxpix];
431  if(logxpx > 1) {
432  qlx=xsum[lxpix];
433  } else {
434  qlx = qfx;
435  }
436 
437 // Now calculate the mean bias correction and uncertainties
438 
439  float qxfrac = (qfx-qlx)/(qfx+qlx);
440  bias = templ.xflcorr(binq,qxfrac)+templ.xavg(binq);
441 
442 // uncertainty and final correction depend upon charge bin
443 
444  xrec = (0.125f*binl+BSHX-2.5f+rat*(binh-binl)*0.125f-(float)shiftx+originx)*xsize - bias;
445  sigmax = templ.xrms(binq);
446 
447 // Do goodness of fit test in x
448 
449  if(rnorm <= 0.f) {LOGERROR("SiStripTemplateReco") << "illegal chi2x normalization (2) = " << rnorm << ENDL; rnorm = 1.;}
450  chi2x=ss2-2.f/rnorm*ssa-2.f/rnorm*rat*ssba+(sa2+2.f*rat*saba+rat*rat*sba2)/(rnorm*rnorm)-templ.chi2xmin(binq);
451  if(chi2x < 0.0) {chi2x = 0.0;}
452  meanx = templ.chi2xavg(binq);
453  if(meanx < 0.01) {meanx = 0.01;}
454  // gsl function that calculates the chi^2 tail prob for non-integral dof
455  // probx = gsl_cdf_chisq_Q(chi2x, meanx);
456  // probx = ROOT::Math::chisquared_cdf_c(chi2x, meanx, trx0);
457  hchi2 = chi2x/2.; hndof = meanx/2.;
458  probx = 1. - TMath::Gamma(hndof, hchi2);
459 
460 // Now choose the better result
461 
462  bias = templ.xavg(binq);
463  biasbcn = templ.xavgbcn(binq);
464  sigmaxbcn = templ.xrmsbcn(binq);
465 
466  if((bias*bias+sigmax*sigmax) > (biasbcn*biasbcn+sigmaxbcn*sigmaxbcn)) {
467 
468  xrec = barycenter - biasbcn;
469  sigmax = sigmaxbcn;
470 
471  }
472 
473  }
474 
475 // Don't return exact zeros for the probability
476 
477  if(probx < probmin) {probx = probmin;}
478 
479 // Decide whether to generate a cluster charge probability
480 
481  if(calc_probQ) {
482 
483 // Calculate the Vavilov probability that the cluster charge is OK
484 
485  templ.vavilov_pars(mpv, sigmaQ, kappa);
486 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
487  if((sigmaQ <=0.) || (mpv <= 0.) || (kappa < 0.01) || (kappa > 9.9)) {
488  throw cms::Exception("DataCorrupt") << "SiStripTemplateReco::Vavilov parameters mpv/sigmaQ/kappa = " << mpv << "/" << sigmaQ << "/" << kappa << std::endl;
489  }
490 #else
491  assert((sigmaQ > 0.) && (mpv > 0.) && (kappa > 0.01) && (kappa < 10.));
492 #endif
493  xvav = ((double)qtotal-mpv)/sigmaQ;
494  beta2 = 1.;
495  if(use_VVIObj) {
496 // VVIObj is a private port of CERNLIB VVIDIS
497  sistripvvi::VVIObj vvidist(kappa, beta2, 1);
498  prvav = vvidist.fcn(xvav);
499  } else {
500 // Use faster but less accurate TMath Vavilov distribution function
501  prvav = TMath::VavilovI(xvav, kappa, beta2);
502  }
503 // Change to upper tail probability
504 // if(prvav > 0.5) prvav = 1. - prvav;
505 // probQ = (float)(2.*prvav);
506  probQ = 1. - prvav;
507  if(probQ < probQmin) {probQ = probQmin;}
508  } else {
509  probQ = -1;
510  }
511 
512  return 0;
513 } // StripTempReco2D
dbl * delta
Definition: mlp_gen.cc:36
int i
Definition: DBlmapReader.cc:9
float chi2xavg(int i)
averaage x chi^2 in 4 charge bins
#define BSXSIZE
float xsize()
strip x-size (microns)
void xtemp(int fxbin, int lxbin, float xtemplate[41][17+4])
void vavilov_pars(double &mpv, double &sigma, double &kappa)
assert(m_qm.get())
void xsigma2(int fxstrp, int lxstrp, float sxthr, float xsum[17+4], float xsig2[17+4])
#define BSHX
float qavg()
average cluster charge for this set of track angles
float chi2xminone()
//!&lt; minimum of x chi^2 for 1 strip clusters
static const int maxpix
#define ENDL
float xrmsbcn(int i)
1st pass chi2 min search: average x-rms of reconstruction binned in 4 charge bins ...
float lorxwidth()
signed lorentz x-width (microns)
float xflcorr(int binq, float qflx)
float qmin()
minimum cluster charge for valid hit (keeps 99.9% of simulated hits)
int cxtemp()
Return central pixel of x-template pixels above readout threshold.
float xavg(int i)
average x-bias of reconstruction binned in 4 charge bins
float dxone()
mean offset/correction for one strip x-clusters
int j
Definition: DBlmapReader.cc:9
if(c.getParameter< edm::InputTag >("puppiValueMap").label().size()!=0)
double f[11][100]
float qscale()
charge scaling factor
float xavgbcn(int i)
1st pass chi2 min search: average x-bias of reconstruction binned in 4 charge bins ...
float chi2xavgone()
//!&lt; average x chi^2 for 1 strip clusters
float chi2xmin(int i)
minimum y chi^2 in 4 charge bins
float s50()
1/2 of the strip threshold signal in electrons
static const int theVerboseLevel
float sxone()
rms for one strip x-clusters
float sxmax()
average strip signal for x-projection of cluster
dbl * Gamma
Definition: mlp_gen.cc:38
static const G4double kappa
#define LOGERROR(x)
const Int_t xsize
bool interpolate(int id, float cotalpha, float cotbeta, float locBy)
#define LOGDEBUG(x)
#define TSXSIZE
#define BSXM2