]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDseed.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Removing DISPLAY
[u/mrichter/AliRoot.git] / TRD / AliTRDseed.cxx
CommitLineData
75bd7f81 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18///////////////////////////////////////////////////////////////////////////////
19// //
20// The seed of a local TRD track //
21// //
22// //
23///////////////////////////////////////////////////////////////////////////////
24
25#include "AliTRDseed.h"
26#include "AliTRDcluster.h"
27#include "AliTRDtracker.h"
28
29#include "TMath.h"
30#include "TLinearFitter.h"
31
32ClassImp(AliTRDseed)
33
34//_____________________________________________________________________________
35AliTRDseed::AliTRDseed() :
36 TObject(),
37 fTilt(0),
38 fPadLength(0),
39 fX0(0),
40 fSigmaY(0),
41 fSigmaY2(0),
42 fMeanz(0),
43 fZProb(0),
44 fN(0),
45 fN2(0),
46 fNUsed(0),
47 fFreq(0),
48 fNChange(0),
49 fMPads(0),
50 fC(0),
51 fCC(0),
52 fChi2(1e10),
53 fChi2Z(1e10)
54{
55 //
56 // Default constructor
57 //
58
59 for (Int_t i = 0; i < 25; i++) {
60 fX[i] = 0; // x position
61 fY[i] = 0; // y position
62 fZ[i] = 0; // z position
63 fIndexes[i] = 0; // indexes
64 fClusters[i] = 0x0; // clusters
65 fUsable[i] = 0; // indication - usable cluster
66 }
67 for (Int_t i = 0; i < 2; i++) {
68 fYref[i] = 0; // reference y
69 fZref[i] = 0; // reference z
70 fYfit[i] = 0; // y fit position +derivation
71 fYfitR[i] = 0; // y fit position +derivation
72 fZfit[i] = 0; // z fit position
73 fZfitR[i] = 0; // z fit position
74 fLabels[i] = 0; // labels
75 }
76
77}
78
79//_____________________________________________________________________________
80void AliTRDseed::Reset()
81{
82 //
83 // Reset seed
84 //
85
86 for (Int_t i = 0; i < 25; i++) {
87 fX[i] = 0; // !x position
88 fY[i] = 0; // !y position
89 fZ[i] = 0; // !z position
90 fIndexes[i] = 0; // !indexes
91 fClusters[i] = 0; // !clusters
92 fUsable[i] = kFALSE;
93 }
94 for (Int_t i = 0; i < 2; i++) {
95 fYref[i] = 0; // reference y
96 fZref[i] = 0; // reference z
97 fYfit[i] = 0; // y fit position +derivation
98 fYfitR[i] = 0; // y fit position +derivation
99 fZfit[i] = 0; // z fit position
100 fZfitR[i] = 0; // z fit position
101 fLabels[i] =-1; // labels
102 }
103 fSigmaY = 0; //"robust" sigma in y
104 fSigmaY2 = 0; //"robust" sigma in y
105 fMeanz = 0; // mean vaue of z
106 fZProb = 0; // max probbable z
107 fMPads = 0;
108 fN = 0; // number of associated clusters
109 fN2 = 0; // number of not crossed
110 fNUsed = 0; // number of used clusters
111 fNChange = 0; // change z counter
112
113}
114
115//_____________________________________________________________________________
116void AliTRDseed::CookLabels()
117{
118 //
119 // Cook 2 labels for seed
120 //
121
122 Int_t labels[200];
123 Int_t out[200];
124 Int_t nlab =0;
125 for (Int_t i=0;i<25;i++){
126 if (!fClusters[i]) continue;
127 for (Int_t ilab=0;ilab<3;ilab++){
128 if (fClusters[i]->GetLabel(ilab)>=0){
129 labels[nlab] = fClusters[i]->GetLabel(ilab);
130 nlab++;
131 }
132 }
133 }
134 Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE);
135 fLabels[0] = out[0];
136 if (nlab2>1 && out[3]>1) fLabels[1] =out[2];
137
138}
139
140//_____________________________________________________________________________
141void AliTRDseed::UseClusters()
142{
143 //
144 // Use clusters
145 //
146
147 for (Int_t i=0;i<25;i++){
148 if (!fClusters[i]) continue;
149 if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
150 }
151
152}
153
154//_____________________________________________________________________________
155void AliTRDseed::Update()
156{
157 //
158 // Update the seed.
159 //
160
161 const Float_t kRatio = 0.8;
162 const Int_t kClmin = 6;
163 const Float_t kmaxtan = 2;
164
165 if (TMath::Abs(fYref[1])>kmaxtan) return; // too much inclined track
166
167 Float_t sigmaexp = 0.05+TMath::Abs(fYref[1]*0.25); // expected r.m.s in y direction
168 Float_t ycrosscor = fPadLength*fTilt*0.5; // y correction for crossing
169 fNChange =0;
170
171 Double_t sumw, sumwx,sumwx2;
172 Double_t sumwy, sumwxy, sumwz,sumwxz;
173 Int_t zints[25]; // histograming of the z coordinate
174 // get 1 and second max probable coodinates in z
175 Int_t zouts[50]; //
176 Float_t allowedz[25]; // allowed z for given time bin
177 Float_t yres[25]; // residuals from reference
178 Float_t anglecor = fTilt*fZref[1]; //correction to the angle
179
180
181 fN=0; fN2 =0;
182 for (Int_t i=0;i<25;i++){
183 yres[i] =10000;
184 if (!fClusters[i]) continue;
185 yres[i] = fY[i]-fYref[0]-(fYref[1]+anglecor)*fX[i]; // residual y
186 zints[fN] = Int_t(fZ[i]);
187 fN++;
188 }
189
190 if (fN<kClmin) return;
191 Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
192 fZProb = zouts[0];
193 if (nz<=1) zouts[3]=0;
194 if (zouts[1]+zouts[3]<kClmin) return;
195
196 if (TMath::Abs(zouts[0]-zouts[2])>12.) zouts[3]=0; // z distance bigger than pad - length
197
198 Int_t breaktime = -1;
199 Bool_t mbefore = kFALSE;
200 Int_t cumul[25][2];
201 Int_t counts[2]={0,0};
202
203 if (zouts[3]>=3){
204 //
205 // find the break time allowing one chage on pad-rows with maximal numebr of accepted clusters
206 //
207 fNChange=1;
208 for (Int_t i=0;i<25;i++){
209 cumul[i][0] = counts[0];
210 cumul[i][1] = counts[1];
211 if (TMath::Abs(fZ[i]-zouts[0])<2) counts[0]++;
212 if (TMath::Abs(fZ[i]-zouts[2])<2) counts[1]++;
213 }
214 Int_t maxcount = 0;
215 for (Int_t i=0;i<24;i++) {
216 Int_t after = cumul[24][0]-cumul[i][0];
217 Int_t before = cumul[i][1];
218 if (after+before>maxcount) {
219 maxcount=after+before;
220 breaktime=i;
221 mbefore=kFALSE;
222 }
223 after = cumul[24][1]-cumul[i][1];
224 before = cumul[i][0];
225 if (after+before>maxcount) {
226 maxcount=after+before;
227 breaktime=i;
228 mbefore=kTRUE;
229 }
230 }
231 breaktime-=1;
232 }
233 for (Int_t i=0;i<25;i++){
234 if (i>breaktime) allowedz[i] = mbefore ? zouts[2]:zouts[0];
235 if (i<=breaktime) allowedz[i] = (!mbefore) ? zouts[2]:zouts[0];
236 }
237 if ( (allowedz[0]>allowedz[24] && fZref[1]<0) || (allowedz[0]<allowedz[24] && fZref[1]>0)){
238 //
239 // tracklet z-direction not in correspondance with track z direction
240 //
241 fNChange =0;
242 for (Int_t i=0;i<25;i++){
243 allowedz[i] = zouts[0]; //only longest taken
244 }
245 }
246
247 if (fNChange>0){
248 //
249 // cross pad -row tracklet - take the step change into account
250 //
251 for (Int_t i=0;i<25;i++){
252 if (!fClusters[i]) continue;
253 if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
254 yres[i] = fY[i]-fYref[0]-(fYref[1]+anglecor)*fX[i]; // residual y
255 if (TMath::Abs(fZ[i]-fZProb)>2){
256 if (fZ[i]>fZProb) yres[i]+=fTilt*fPadLength;
257 if (fZ[i]<fZProb) yres[i]-=fTilt*fPadLength;
258 }
259 }
260 }
261
262 Double_t yres2[25];
263 Double_t mean,sigma;
264 for (Int_t i=0;i<25;i++){
265 if (!fClusters[i]) continue;
266 if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
267 yres2[fN2] = yres[i];
268 fN2++;
269 }
270 if (fN2<kClmin){
271 fN2 = 0;
272 return;
273 }
274 EvaluateUni(fN2,yres2,mean,sigma,Int_t(fN2*kRatio-2));
275 if (sigma<sigmaexp*0.8) sigma=sigmaexp;
276 fSigmaY = sigma;
277 //
278 //
279 // reset sums
280 sumw=0; sumwx=0; sumwx2=0;
281 sumwy=0; sumwxy=0; sumwz=0;sumwxz=0;
282 fN2 =0;
283 fMeanz =0;
284 fMPads =0;
285 //
286 for (Int_t i=0;i<25;i++){
287 fUsable[i]=kFALSE;
288 if (!fClusters[i]) continue;
289 if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
290 if (TMath::Abs(yres[i]-mean)>4.*sigma) continue;
291 fUsable[i] = kTRUE;
292 fN2++;
293 fMPads+=fClusters[i]->GetNPads();
294 Float_t weight =1;
295 if (fClusters[i]->GetNPads()>4) weight=0.5;
296 if (fClusters[i]->GetNPads()>5) weight=0.2;
297 //
298 Double_t x = fX[i];
299 sumw+=weight; sumwx+=x*weight; sumwx2+=x*x*weight;
300 sumwy+=weight*yres[i]; sumwxy+=weight*(yres[i])*x;
301 sumwz+=weight*fZ[i]; sumwxz+=weight*fZ[i]*x;
302 }
303 if (fN2<kClmin){
304 fN2 = 0;
305 return;
306 }
307 fMeanz = sumwz/sumw;
308 Float_t correction =0;
309 if (fNChange>0){
310 // tracklet on boundary
311 if (fMeanz<fZProb) correction = ycrosscor;
312 if (fMeanz>fZProb) correction = -ycrosscor;
313 }
314 Double_t det = sumw*sumwx2-sumwx*sumwx;
315 fYfitR[0] = (sumwx2*sumwy-sumwx*sumwxy)/det;
316 fYfitR[1] = (sumw*sumwxy-sumwx*sumwy)/det;
317
318 fSigmaY2 =0;
319 for (Int_t i=0;i<25;i++){
320 if (!fUsable[i]) continue;
321 Float_t delta = yres[i]-fYfitR[0]-fYfitR[1]*fX[i];
322 fSigmaY2+=delta*delta;
323 }
324 fSigmaY2 = TMath::Sqrt(fSigmaY2/Float_t(fN2-2));
325 //
326 fZfitR[0] = (sumwx2*sumwz-sumwx*sumwxz)/det;
327 fZfitR[1] = (sumw*sumwxz-sumwx*sumwz)/det;
328 fZfit[0] = (sumwx2*sumwz-sumwx*sumwxz)/det;
329 fZfit[1] = (sumw*sumwxz-sumwx*sumwz)/det;
330 fYfitR[0] += fYref[0]+correction;
331 fYfitR[1] += fYref[1];
332 fYfit[0] = fYfitR[0];
333 fYfit[1] = fYfitR[1];
334
335 UpdateUsed();
336
337}
338
339//_____________________________________________________________________________
340void AliTRDseed::UpdateUsed()
341{
342 //
343 // Update used seed
344 //
345
346 fNUsed = 0;
347 for (Int_t i=0;i<25;i++){
348 if (!fClusters[i]) continue;
349 if ((fClusters[i]->IsUsed())) fNUsed++;
350 }
351
352}
353
354//_____________________________________________________________________________
355void AliTRDseed::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean
356 , Double_t &sigma, Int_t hh)
357{
358 //
359 // Robust estimator in 1D case MI version
360 //
361 // For the univariate case
362 // estimates of location and scatter are returned in mean and sigma parameters
363 // the algorithm works on the same principle as in multivariate case -
364 // it finds a subset of size hh with smallest sigma, and then returns mean and
365 // sigma of this subset
366 //
367
368 if (hh==0)
369 hh=(nvectors+2)/2;
370 Double_t faclts[]={2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660,1.879,1.7973,1.7203,1.6473};
371 Int_t *index=new Int_t[nvectors];
372 TMath::Sort(nvectors, data, index, kFALSE);
373
374 Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
375 Double_t factor = faclts[nquant-1];
376
377 Double_t sumx =0;
378 Double_t sumx2 =0;
379 Int_t bestindex = -1;
380 Double_t bestmean = 0;
381 Double_t bestsigma = data[index[nvectors-1]]-data[index[0]]; // maximal possible sigma
382 for (Int_t i=0; i<hh; i++){
383 sumx += data[index[i]];
384 sumx2 += data[index[i]]*data[index[i]];
385 }
386
387 Double_t norm = 1./Double_t(hh);
388 Double_t norm2 = 1./Double_t(hh-1);
389 for (Int_t i=hh; i<nvectors; i++){
390 Double_t cmean = sumx*norm;
391 Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2;
392 if (csigma<bestsigma){
393 bestmean = cmean;
394 bestsigma = csigma;
395 bestindex = i-hh;
396 }
397
398 sumx += data[index[i]]-data[index[i-hh]];
399 sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
400 }
401
402 Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
403 mean = bestmean;
404 sigma = bstd;
405 delete [] index;
406
407}
408
409//_____________________________________________________________________________
410Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
411{
412 //
413 // Fit the Rieman tilt
414 //
415
416 TLinearFitter fitterT2(4,"hyp4"); // fitting with tilting pads - kz not fixed
417 fitterT2.StoreData(kTRUE);
418 Float_t xref2 = (cseed[2].fX0+cseed[3].fX0)*0.5; // reference x0 for z
419 //
420 Int_t npointsT =0;
421 fitterT2.ClearPoints();
422 for (Int_t iLayer=0; iLayer<6;iLayer++){
423 if (!cseed[iLayer].IsOK()) continue;
424 Double_t tilt = cseed[iLayer].fTilt;
425
426 for (Int_t itime=0;itime<25;itime++){
427 if (!cseed[iLayer].fUsable[itime]) continue;
428 // x relative to the midle chamber
429 Double_t x = cseed[iLayer].fX[itime]+cseed[iLayer].fX0-xref2;
430 Double_t y = cseed[iLayer].fY[itime];
431 Double_t z = cseed[iLayer].fZ[itime];
432
433 //
434 // tilted rieman
435 //
436 Double_t uvt[6];
437 Double_t x2 = cseed[iLayer].fX[itime]+cseed[iLayer].fX0; // global x
438 Double_t t = 1./(x2*x2+y*y);
439 uvt[1] = t; // t
440 uvt[0] = 2.*x2*uvt[1]; // u
441 uvt[2] = 2.0*tilt*uvt[1];
442 uvt[3] = 2.0*tilt*x*uvt[1];
443 uvt[4] = 2.0*(y+tilt*z)*uvt[1];
444 //
445 Double_t error = 2*uvt[1];
446 if (terror) error*=cseed[iLayer].fSigmaY;
447 else {error *=0.2;} //default error
448 fitterT2.AddPoint(uvt,uvt[4],error);
449 npointsT++;
450 }
451 }
452 fitterT2.Eval();
453 Double_t rpolz0 = fitterT2.GetParameter(3);
454 Double_t rpolz1 = fitterT2.GetParameter(4);
455 //
456 // linear fitter - not possible to make boundaries
457 // non accept non possible z and dzdx combination
458 //
459 Bool_t acceptablez =kTRUE;
460 for (Int_t iLayer=0; iLayer<6;iLayer++){
461 if (cseed[iLayer].IsOK()){
462 Double_t zT2 = rpolz0+rpolz1*(cseed[iLayer].fX0 - xref2);
463 if (TMath::Abs(cseed[iLayer].fZProb-zT2)>cseed[iLayer].fPadLength*0.5+1)
464 acceptablez = kFALSE;
465 }
466 }
467 if (!acceptablez){
468 Double_t zmf = cseed[2].fZref[0]+cseed[2].fZref[1]*(xref2-cseed[2].fX0);
469 Double_t dzmf = (cseed[2].fZref[1]+ cseed[3].fZref[1])*0.5;
470 fitterT2.FixParameter(3,zmf);
471 fitterT2.FixParameter(4,dzmf);
472 fitterT2.Eval();
473 fitterT2.ReleaseParameter(3);
474 fitterT2.ReleaseParameter(4);
475 rpolz0 = fitterT2.GetParameter(3);
476 rpolz1 = fitterT2.GetParameter(4);
477 }
478
479 Double_t chi2TR = fitterT2.GetChisquare()/Float_t(npointsT);
480 Double_t params[3];
481 params[0] = fitterT2.GetParameter(0);
482 params[1] = fitterT2.GetParameter(1);
483 params[2] = fitterT2.GetParameter(2);
484 Double_t curvature = 1+params[1]*params[1]-params[2]*params[0];
485 for (Int_t iLayer = 0; iLayer<6;iLayer++){
486 Double_t x = cseed[iLayer].fX0;
487 Double_t y=0,dy=0, z=0, dz=0;
488 // y
489 Double_t res2 = (x*params[0]+params[1]);
490 res2*=res2;
491 res2 = 1.-params[2]*params[0]+params[1]*params[1]-res2;
492 if (res2>=0){
493 res2 = TMath::Sqrt(res2);
494 y = (1-res2)/params[0];
495 }
496 //dy
497 Double_t x0 = -params[1]/params[0];
498 if (-params[2]*params[0]+params[1]*params[1]+1>0){
499 Double_t rm1 = params[0]/TMath::Sqrt(-params[2]*params[0]+params[1]*params[1]+1);
500 if ( 1./(rm1*rm1)-(x-x0)*(x-x0)>0){
501 Double_t res = (x-x0)/TMath::Sqrt(1./(rm1*rm1)-(x-x0)*(x-x0));
502 if (params[0]<0) res*=-1.;
503 dy = res;
504 }
505 }
506 z = rpolz0+rpolz1*(x-xref2);
507 dz = rpolz1;
508 cseed[iLayer].fYref[0] = y;
509 cseed[iLayer].fYref[1] = dy;
510 cseed[iLayer].fZref[0] = z;
511 cseed[iLayer].fZref[1] = dz;
512 cseed[iLayer].fC = curvature;
513
514 }
515
516 return chi2TR;
517
518}