CMS 3D CMS Logo
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, Exception, f, sistripvvi::VVIObj::fcn(), Gamma, mps_fire::i, createfilelist::int, SiStripTemplate::interpolate(), gen::k, kappa, LOGDEBUG, LOGERROR, maxpix, min(), SiStripTemplate::qavg(), SiStripTemplate::qmin(), SiStripTemplate::qscale(), SiStripTemplate::s50(), ctppsDiamondLocalTracks_cfi::sigma, SiStripTemplate::sxmax(), SiStripTemplate::sxone(), theVerboseLevel, TSXSIZE, SiStripTemplate::vavilov2_pars(), SiStripTemplate::xavgc2m(), SiStripTemplate::xrmsc2m(), SiStripTemplate::xsigma2(), xsize, 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  xrec1 = (0.125f*(minbink-xcbin)+BSHX-(float)shiftx+originx)*xsize - bias;
399  xrec2 = (0.125f*(minbinj-xcbin)+BSHX-(float)shiftx+originx)*xsize - bias;
400  sigmax = sqrt2x*templ.xrmsc2m(binqerr);
401 
402 
403 // Do goodness of fit test in y
404 
405  if(chi2xmin < 0.0) {chi2xmin = 0.0;}
406  meanx = templ.chi2xavgc2m(binq);
407  if(meanx < 0.01) {meanx = 0.01;}
408  hchi2 = chi2xmin/2.; hndof = meanx/2.;
409  prob2x = 1. - TMath::Gamma(hndof, hchi2);
410  }
411 
412  // Don't return exact zeros for the probability
413 
414  if(prob2x < probmin) {prob2x = probmin;}
415 
416 
417  return 0;
418 } // StripTempSplit
delta
dbl * delta
Definition: mlp_gen.cc:36
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
mps_fire.i
i
Definition: mps_fire.py:156
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
SiStripTemplate::xsigma2
void xsigma2(int fxstrp, int lxstrp, float sxthr, float xsum[17+4], float xsig2[17+4])
Definition: SiStripTemplate.cc:1017
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()
//!< minimum of x chi^2 for 1 strip clusters
Definition: SiStripTemplate.h:317
maxpix
static const int maxpix
Definition: SiPixelLorentzAngle.h:49
SiStripTemplate::qmin
float qmin()
minimum cluster charge for valid hit (keeps 99.9% of simulated hits)
Definition: SiStripTemplate.h:192
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
SiStripTemplate::dxone
float dxone()
mean offset/correction for one strip x-clusters
Definition: SiStripTemplate.h:190
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
SiStripTemplate::qscale
float qscale()
charge scaling factor
Definition: SiStripTemplate.h:187
min
T min(T a, T b)
Definition: MathUtil.h:58
PixelCPETemplateReco_cfi.speed
speed
Definition: PixelCPETemplateReco_cfi.py:5
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
SiStripTemplate::chi2xavgone
float chi2xavgone()
//!< average x chi^2 for 1 strip clusters
Definition: SiStripTemplate.h:316
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
ctppsDiamondLocalTracks_cfi.sigma
sigma
Definition: ctppsDiamondLocalTracks_cfi.py:10
Gamma
dbl * Gamma
Definition: mlp_gen.cc:38
LOGDEBUG
#define LOGDEBUG(x)
Definition: SiStripTemplateSplit.cc:36
kappa
static const G4double kappa
Definition: UrbanMscModel93.cc:34
theVerboseLevel
int theVerboseLevel
Definition: SiStripTemplateSplit.cc:37
xsize
const Int_t xsize
Definition: trackSplitPlot.h:42
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 438 of file SiStripTemplateSplit.cc.

References PixelCPETemplateReco_cfi::speed, and StripTempSplit().

442 {
443  // Local variables
444  const int speed = 1;
445 
446  return SiStripTemplateSplit::StripTempSplit(id, cotalpha, cotbeta, locBy, speed, cluster, templ, xrec1, xrec2, sigmax, prob2x, q2bin, prob2Q);
447 } // StripTempSplit
PixelCPETemplateReco_cfi.speed
speed
Definition: PixelCPETemplateReco_cfi.py:5
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.11