CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Functions
SiStripTemplateSplit Namespace Reference

Functions

int StripTempSplit (int id, float cotalpha, float cotbeta, float locBy, int speed, std::vector< float > &cluster, SiStripTemplate &templ, float &xrec1, float &xrec2, float &sigmax, float &prob2x, int &q2bin, float &prob2Q)
 
int StripTempSplit (int id, float cotalpha, float cotbeta, float locBy, std::vector< float > &cluster, SiStripTemplate &templ, float &xrec1, float &xrec2, float &sigmax, float &prob2x, int &q2bin, float &prob2Q)
 

Function Documentation

int SiStripTemplateSplit::StripTempSplit ( int  id,
float  cotalpha,
float  cotbeta,
float  locBy,
int  speed,
std::vector< float > &  cluster,
SiStripTemplate &  templ,
float &  xrec1,
float &  xrec2,
float &  sigmax,
float &  prob2x,
int &  q2bin,
float &  prob2Q 
)

Reconstruct the best estimate of the hit positions for pixel clusters.

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
speed- (input) switch (-1->2) trading speed vs robustness -1 totally bombproof, searches the entire ranges of template bins, calculates Q probability w/ VVIObj 0 totally bombproof, searches the entire template bin range at full density (no Qprob) 1 faster, searches same range as 0 but at 1/2 density 2 fastest, searches same range as 1 but at 1/4 density (no big pix) and 1/2 density (big pix in cluster)
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].
templ- (input) the template used in the reconstruction
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 xrec1 and xrec2 in microns
prob2x- (output) probability describing goodness-of-fit to a merged cluster hypothesis for x-reco
q2bin- (output) index (0-4) describing the charge of the cluster assuming a merged 2-hit cluster hypothesis [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]
prob2Q- (output) probability that the cluster charge is compatible with a 2-hit merging

Definition at line 75 of file SiStripTemplateSplit.cc.

References BSHX, BSXM2, BSXSIZE, SiStripTemplate::chi2xavgc2m(), SiStripTemplate::chi2xavgone(), SiStripTemplate::chi2xminone(), SiStripTemplate::cxtemp(), delta, SiStripTemplate::dxone(), ENDL, edm::hlt::Exception, f, sistripvvi::VVIObj::fcn(), Gamma, i, SiStripTemplate::interpolate(), j, gen::k, LOGDEBUG, LOGERROR, max(), maxpix, bookConverter::min, SiStripTemplate::qavg(), SiStripTemplate::qmin(), SiStripTemplate::qscale(), SiStripTemplate::s50(), SiStripTemplate::sxmax(), SiStripTemplate::sxone(), theVerboseLevel, TSXSIZE, SiStripTemplate::vavilov2_pars(), SiStripTemplate::xavgc2m(), SiStripTemplate::xrmsc2m(), SiStripTemplate::xsigma2(), SiStripTemplate::xsize(), SiStripTemplate::xtemp3d(), and SiStripTemplate::xtemp3d_int().

Referenced by TrackClusterSplitter::splitCluster(), and StripTempSplit().

79 {
80  // Local variables
81  int i, j, k, binq, binqerr, midpix;
82  int fxpix, nxpix, lxpix, logxpx, shiftx;
83  int nclusx;
84  int nxbin, xcbin, minbinj, minbink;
85  int deltaj, jmin, jmax, kmin, kmax, km, fxbin, lxbin, djx;
86  float sxthr, delta, sigma, pseudopix, xsize, qscale;
87  float ss2, ssa, sa2, rat, fq, qtotal, qavg;
88  float originx, bias, maxpix;
89  double chi2x, meanx, chi2xmin, chi21max;
90  double hchi2, hndof;
91  double prvav, mpv, sigmaQ, kappa, xvav, beta2;
92  float xsum[BSXSIZE];
93  float xsig2[BSXSIZE];
94  float xw2[BSXSIZE], xsw[BSXSIZE];
95  const float sqrt2x={2.00000};
96  const float sqrt12={3.4641};
97  const float probmin={1.110223e-16};
98  const float prob2Qmin={1.e-5};
99 
100 // bool SiStripTemplateSplit::SimpleTemplate2D(float cotalpha, float cotbeta, float xhit, float yhit, float thick, float lorxwidth, float lorywidth,
101 // float qavg, std::vector<bool> ydouble, std::vector<bool> xdouble, float template2d[BSXM2][BYM2]);
102 
103 // The minimum chi2 for a valid one pixel cluster = pseudopixel contribution only
104 
105  const double mean1pix={0.100}, chi21min={0.160};
106 
107 // First, interpolate the template needed to analyze this cluster
108 // check to see of the track direction is in the physical range of the loaded template
109 
110  if(!templ.interpolate(id, cotalpha, cotbeta, locBy)) {
111  LOGDEBUG("SiStripTemplateReco") << "input cluster direction cot(alpha) = " << cotalpha << ", cot(beta) = " << cotbeta
112  << " is not within the acceptance of template ID = " << id << ", no reconstruction performed" << ENDL;
113  return 20;
114  }
115 
116  // Get pixel dimensions from the template (to allow multiple detectors in the future)
117 
118  xsize = templ.xsize();
119 
120  // Define size of pseudopixel
121 
122  pseudopix = templ.s50();
123 
124  // Get charge scaling factor
125 
126  qscale = templ.qscale();
127 
128  // enforce maximum size
129 
130  nclusx = (int)cluster.size();
131 
132  if(nclusx > TSXSIZE) {nclusx = TSXSIZE;}
133 
134  // First, rescale all strip charges, sum them and trunate the strip charges
135 
136  qtotal = 0.f;
137  for(i=0; i<BSXSIZE; ++i) {xsum[i] = 0.f;}
138  maxpix = 2.f*templ.sxmax();
139  for(j=0; j<nclusx; ++j) {
140  xsum[j] = qscale*cluster[j];
141  qtotal += xsum[j];
142  if(xsum[j] > maxpix) {xsum[j] = maxpix;}
143  }
144 
145  // next, identify the x-cluster ends, count total pixels, nxpix, and logical pixels, logxpx
146 
147  fxpix = -1;
148  nxpix=0;
149  lxpix=0;
150  logxpx=0;
151  for(i=0; i<BSXSIZE; ++i) {
152  if(xsum[i] > 0.f) {
153  if(fxpix == -1) {fxpix = i;}
154  ++logxpx;
155  ++nxpix;
156  lxpix = i;
157  }
158  }
159 
160 
161  // dlengthx = (float)nxpix - templ.clslenx();
162 
163  // Make sure cluster is continuous
164 
165  if((lxpix-fxpix+1) != nxpix) {
166 
167  LOGDEBUG("SiStripTemplateReco") << "x-length of pixel cluster doesn't agree with number of pixels above threshold" << ENDL;
168  if (theVerboseLevel > 2) {
169  LOGDEBUG("SiStripTemplateReco") << "xsum[] = ";
170  for(i=0; i<BSXSIZE-1; ++i) {LOGDEBUG("SiStripTemplateReco") << xsum[i] << ", ";}
171  LOGDEBUG("SiStripTemplateReco") << ENDL;
172  }
173 
174  return 2;
175  }
176 
177  // If cluster is longer than max template size, technique fails
178 
179  if(nxpix > TSXSIZE) {
180 
181  LOGDEBUG("SiStripTemplateReco") << "x-length of pixel cluster is larger than maximum template size" << ENDL;
182  if (theVerboseLevel > 2) {
183  LOGDEBUG("SiStripTemplateReco") << "xsum[] = ";
184  for(i=0; i<BSXSIZE-1; ++i) {LOGDEBUG("SiStripTemplateReco") << xsum[i] << ", ";}
185  LOGDEBUG("SiStripTemplateReco") << ENDL;
186  }
187 
188  return 7;
189  }
190 
191  // next, center the cluster on template center if necessary
192 
193  midpix = (fxpix+lxpix)/2;
194  shiftx = templ.cxtemp() - midpix;
195  if(shiftx > 0) {
196  for(i=lxpix; i>=fxpix; --i) {
197  xsum[i+shiftx] = xsum[i];
198  xsum[i] = 0.f;
199  }
200  } else if (shiftx < 0) {
201  for(i=fxpix; i<=lxpix; ++i) {
202  xsum[i+shiftx] = xsum[i];
203  xsum[i] = 0.f;
204  }
205  }
206  lxpix +=shiftx;
207  fxpix +=shiftx;
208 
209  // If the cluster boundaries are OK, add pesudopixels, otherwise quit
210 
211  if(fxpix > 1 && fxpix <BSXM2) {
212  xsum[fxpix-1] = pseudopix;
213  xsum[fxpix-2] = 0.2f*pseudopix;
214  } else {return 9;}
215  if(lxpix > 1 && lxpix < BSXM2) {
216  xsum[lxpix+1] = pseudopix;
217  xsum[lxpix+2] = 0.2f*pseudopix;
218  } else {return 9;}
219 
220  // finally, determine if pixel[0] is a double pixel and make an origin correction if it is
221 
222  originx = 0.f;
223 
224 // uncertainty and final corrections depend upon total charge bin
225 
226  qavg = templ.qavg();
227  fq = qtotal/qavg;
228  if(fq > 3.0f) {
229  binq=0;
230  } else {
231  if(fq > 2.0f) {
232  binq=1;
233  } else {
234  if(fq > 1.70f) {
235  binq=2;
236  } else {
237  binq=3;
238  }
239  }
240  }
241 
242  // Return the charge bin via the parameter list unless the charge is too small (then flag it)
243 
244  q2bin = binq;
245  if(qtotal < 1.9f*templ.qmin()) {q2bin = 5;} else {
246  if(qtotal < 1.9f*templ.qmin(1)) {q2bin = 4;}
247  }
248  if (theVerboseLevel > 9) {
249  LOGDEBUG("SiStripTemplateReco") <<
250  "ID = " << id <<
251  " cot(alpha) = " << cotalpha << " cot(beta) = " << cotbeta <<
252  " nclusx = " << nclusx << ENDL;
253  }
254 
255 // binqerr is the charge bin for error estimation
256 
257  binqerr = binq;
258  if(binqerr > 2) binqerr = 2;
259 
260 // Calculate the Vavilov probability that the cluster charge is consistent with a merged cluster
261 
262  if(speed < 0) {
263  templ.vavilov2_pars(mpv, sigmaQ, kappa);
264 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
265  if((sigmaQ <=0.) || (mpv <= 0.) || (kappa < 0.01) || (kappa > 9.9)) {
266  throw cms::Exception("DataCorrupt") << "SiStripTemplateSplit::Vavilov parameters mpv/sigmaQ/kappa = " << mpv << "/" << sigmaQ << "/" << kappa << std::endl;
267  }
268 #else
269  assert((sigmaQ > 0.) && (mpv > 0.) && (kappa > 0.01) && (kappa < 10.));
270 #endif
271  xvav = ((double)qtotal-mpv)/sigmaQ;
272  beta2 = 1.;
273 // VVIObj is a private port of CERNLIB VVIDIS
274  sistripvvi::VVIObj vvidist(kappa, beta2, 1);
275  prvav = vvidist.fcn(xvav);
276  prob2Q = 1. - prvav;
277  if(prob2Q < prob2Qmin) {prob2Q = prob2Qmin;}
278  } else {
279  prob2Q = -1.f;
280  }
281 
282 
283 // Next, generate the 3d x-template
284 
285  templ.xtemp3d_int(nxpix, nxbin);
286 
287 // retrieve the number of x-bins
288 
289  xcbin = nxbin/2;
290 
291 // Next, decide on chi^2 min search parameters
292 
293 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
294  if(speed < -1 || speed > 2) {
295  throw cms::Exception("DataCorrupt") << "SiStripTemplateReco::PixelTempReco2D called with illegal speed = " << speed << std::endl;
296  }
297 #else
298  assert(speed >= -1 && speed < 3);
299 #endif
300  fxbin = 0; lxbin = nxbin-1; djx = 1;
301  if(speed > 0) {
302  djx = 2;
303  if(speed > 1) {djx = 4;}
304  }
305 
306 // Define the maximum signal to allow before de-weighting a pixel
307 
308  sxthr = 1.1f*maxpix;
309 
310 // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
311 
312 // for(i=0; i<BYSIZE; ++i) { xsig2[i] = 0.; }
313  templ.xsigma2(fxpix, lxpix, sxthr, xsum, xsig2);
314 
315 // Find the template bin that minimizes the Chi^2
316 
317  chi2xmin = 1.e15;
318  ss2 = 0.f;
319  for(i=fxpix-2; i<=lxpix+2; ++i) {
320  xw2[i] = 1.f/xsig2[i];
321  xsw[i] = xsum[i]*xw2[i];
322  ss2 += xsum[i]*xsw[i];
323  }
324  minbinj = -1;
325  minbink = -1;
326  deltaj = djx;
327  jmin = fxbin;
328  jmax = lxbin;
329  kmin = fxbin;
330  kmax = lxbin;
331  std::vector<float> xtemp(BSXSIZE);
332  while(deltaj > 0) {
333  for(j=jmin; j<jmax; j+=deltaj) {
334  km = std::min(kmax, j);
335  for(k=kmin; k<=km; k+=deltaj) {
336 
337  // Get the template for this set of indices
338 
339  templ.xtemp3d(j, k, xtemp);
340 
341  ssa = 0.f;
342  sa2 = 0.f;
343  for(i=fxpix-2; i<=lxpix+2; ++i) {
344  ssa += xsw[i]*xtemp[i];
345  sa2 += xtemp[i]*xtemp[i]*xw2[i];
346  }
347  rat=ssa/ss2;
348  if(rat <= 0.f) {LOGERROR("SiStripTemplateSplit") << "illegal chi2xmin normalization = " << rat << ENDL; rat = 1.;}
349  chi2x=ss2-2.f*ssa/rat+sa2/(rat*rat);
350  if(chi2x < chi2xmin) {
351  chi2xmin = chi2x;
352  minbinj = j;
353  minbink = k;
354  }
355  }
356  }
357  deltaj /= 2;
358  if(minbinj > fxbin) {jmin = minbinj - deltaj;} else {jmin = fxbin;}
359  if(minbinj < lxbin) {jmax = minbinj + deltaj;} else {jmax = lxbin;}
360  if(minbink > fxbin) {kmin = minbink - deltaj;} else {kmin = fxbin;}
361  if(minbink < lxbin) {kmax = minbink + deltaj;} else {kmax = lxbin;}
362  }
363 
364  if (theVerboseLevel > 9) {
365  LOGDEBUG("SiStripTemplateReco") <<
366  "minbinj/minbink " << minbinj<< "/" << minbink << " chi2xmin = " << chi2xmin << ENDL;
367  }
368 
369 // Do not apply final template pass to 1-pixel clusters (use calibrated offset)
370 
371  if(logxpx == 1) {
372 
373  delta = templ.dxone();
374  sigma = templ.sxone();
375  xrec1 = 0.5f*(fxpix+lxpix-2*shiftx+2.f*originx)*xsize-delta;
376  xrec2 = xrec1;
377  if(sigma <= 0.f) {
378  sigmax = xsize/sqrt12;
379  } else {
380  sigmax = sigma;
381  }
382 
383 // Do probability calculation for one-pixel clusters
384 
385  chi21max = fmax(chi21min, (double)templ.chi2xminone());
386  chi2xmin -=chi21max;
387  if(chi2xmin < 0.) {chi2xmin = 0.;}
388  meanx = fmax(mean1pix, (double)templ.chi2xavgone());
389  hchi2 = chi2xmin/2.; hndof = meanx/2.;
390  prob2x = 1. - TMath::Gamma(hndof, hchi2);
391 
392  } else {
393 
394 
395 // uncertainty and final correction depend upon charge bin
396 
397  bias = templ.xavgc2m(binq);
398  k = std::min(minbink, minbinj);
399  j = std::max(minbink, minbinj);
400  xrec1 = (0.125f*(minbink-xcbin)+BSHX-(float)shiftx+originx)*xsize - bias;
401  xrec2 = (0.125f*(minbinj-xcbin)+BSHX-(float)shiftx+originx)*xsize - bias;
402  sigmax = sqrt2x*templ.xrmsc2m(binqerr);
403 
404 
405 // Do goodness of fit test in y
406 
407  if(chi2xmin < 0.0) {chi2xmin = 0.0;}
408  meanx = templ.chi2xavgc2m(binq);
409  if(meanx < 0.01) {meanx = 0.01;}
410  hchi2 = chi2xmin/2.; hndof = meanx/2.;
411  prob2x = 1. - TMath::Gamma(hndof, hchi2);
412  }
413 
414  // Don't return exact zeros for the probability
415 
416  if(prob2x < probmin) {prob2x = probmin;}
417 
418 
419  return 0;
420 } // StripTempSplit
delta
dbl * delta
Definition: mlp_gen.cc:36
i
int i
Definition: DBlmapReader.cc:9
SiStripTemplate::xavgc2m
float xavgc2m(int i)
1st pass chi2 min search: average x-bias of reconstruction binned in 4 charge bins ...
Definition: SiStripTemplate.h:246
SiStripTemplate::chi2xavgc2m
float chi2xavgc2m(int i)
1st pass chi2 min search: average x-chisq for merged clusters
Definition: SiStripTemplate.h:274
SiStripTemplate::xsigma2
void xsigma2(int fxstrp, int lxstrp, float sxthr, float xsum[BSXSIZE], float xsig2[BSXSIZE])
Definition: SiStripTemplate.cc:1017
SiStripTemplate::vavilov2_pars
void vavilov2_pars(double &mpv, double &sigma, double &kappa)
Definition: SiStripTemplate.cc:1691
BSXSIZE
#define BSXSIZE
Definition: SiStripTemplateDefs.h:25
SiStripTemplate::xsize
float xsize()
strip x-size (microns)
Definition: SiStripTemplate.h:322
SiStripTemplate::xtemp3d
void xtemp3d(int j, int k, std::vector< float > &xtemplate)
Definition: SiStripTemplate.cc:1441
BSHX
#define BSHX
Definition: SiStripTemplateDefs.h:26
SiStripTemplate::qavg
float qavg()
average cluster charge for this set of track angles
Definition: SiStripTemplate.h:186
SiStripTemplate::chi2xminone
float chi2xminone()
//!&lt; minimum of x chi^2 for 1 strip clusters
Definition: SiStripTemplate.h:317
SiStripTemplate::qmin
float qmin()
minimum cluster charge for valid hit (keeps 99.9% of simulated hits)
Definition: SiStripTemplate.h:192
theVerboseLevel
static int theVerboseLevel
Definition: SiPixelTemplateReco.cc:70
SiStripTemplate::xtemp3d_int
void xtemp3d_int(int nxpix, int &nxbins)
Definition: SiStripTemplate.cc:1362
SiStripTemplate::cxtemp
int cxtemp()
Return central pixel of x-template pixels above readout threshold.
Definition: SiStripTemplate.cc:1314
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
SiStripTemplate::dxone
float dxone()
mean offset/correction for one strip x-clusters
Definition: SiStripTemplate.h:190
j
int j
Definition: DBlmapReader.cc:9
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
SiStripTemplate::qscale
float qscale()
charge scaling factor
Definition: SiStripTemplate.h:187
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:81
SiStripTemplate::chi2xavgone
float chi2xavgone()
//!&lt; average x chi^2 for 1 strip clusters
Definition: SiStripTemplate.h:316
maxpix
static const int maxpix
Definition: SiPixelLorentzAngle.h:49
SiStripTemplate::s50
float s50()
1/2 of the strip threshold signal in electrons
Definition: SiStripTemplate.h:188
SiStripTemplate::xrmsc2m
float xrmsc2m(int i)
1st pass chi2 min search: average x-rms of reconstruction binned in 4 charge bins ...
Definition: SiStripTemplate.h:253
SiStripTemplate::sxone
float sxone()
rms for one strip x-clusters
Definition: SiStripTemplate.h:191
LOGERROR
#define LOGERROR(x)
Definition: SiStripTemplateSplit.cc:35
SiStripTemplate::sxmax
float sxmax()
average strip signal for x-projection of cluster
Definition: SiStripTemplate.h:189
Gamma
dbl * Gamma
Definition: mlp_gen.cc:38
LOGDEBUG
#define LOGDEBUG(x)
Definition: SiStripTemplateSplit.cc:36
SiStripTemplate::interpolate
bool interpolate(int id, float cotalpha, float cotbeta, float locBy)
Definition: SiStripTemplate.cc:641
TSXSIZE
#define TSXSIZE
Definition: SiStripTemplateDefs.h:22
BSXM2
#define BSXM2
Definition: SiStripTemplateDefs.h:28
ENDL
#define ENDL
Definition: SiStripTemplateSplit.cc:38
int SiStripTemplateSplit::StripTempSplit ( int  id,
float  cotalpha,
float  cotbeta,
float  locBy,
std::vector< float > &  cluster,
SiStripTemplate &  templ,
float &  xrec1,
float &  xrec2,
float &  sigmax,
float &  prob2x,
int &  q2bin,
float &  prob2Q 
)

Reconstruct the best estimate of the hit positions for pixel clusters.

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 of 7x21 array of pixel signals, origin of local coords (0,0) at center of pixel cluster[0][0].
templ- (input) the template used in the reconstruction
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 xrec1 and xrec2 in microns
prob2x- (output) probability describing goodness-of-fit to a merged cluster hypothesis for x-reco
q2bin- (output) index (0-4) describing the charge of the cluster assuming a merged 2-hit cluster hypothesis [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]
prob2Q- (output) probability that the cluster charge is compatible with a 2-hit merging

Definition at line 440 of file SiStripTemplateSplit.cc.

References StripTempSplit().

444 {
445  // Local variables
446  const int speed = 1;
447 
448  return SiStripTemplateSplit::StripTempSplit(id, cotalpha, cotbeta, locBy, speed, cluster, templ, xrec1, xrec2, sigmax, prob2x, q2bin, prob2Q);
449 } // StripTempSplit
SiStripTemplateSplit::StripTempSplit
int StripTempSplit(int id, float cotalpha, float cotbeta, float locBy, int speed, std::vector< float > &cluster, SiStripTemplate &templ, float &xrec1, float &xrec2, float &sigmax, float &prob2x, int &q2bin, float &prob2Q)
Definition: SiStripTemplateSplit.cc:75
Generated for CMSSW Reference Manual by doxygen 1.8.5