1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // The seed of a local TRD track //
23 ///////////////////////////////////////////////////////////////////////////////
25 #include "AliTRDseed.h"
26 #include "AliTRDcluster.h"
27 #include "AliTRDtracker.h"
30 #include "TLinearFitter.h"
34 //_____________________________________________________________________________
35 AliTRDseed::AliTRDseed() :
56 // Default constructor
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
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
79 //_____________________________________________________________________________
80 void AliTRDseed::Reset()
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
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
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
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
115 //_____________________________________________________________________________
116 void AliTRDseed::CookLabels()
119 // Cook 2 labels for seed
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);
134 Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE);
136 if (nlab2>1 && out[3]>1) fLabels[1] =out[2];
140 //_____________________________________________________________________________
141 void AliTRDseed::UseClusters()
147 for (Int_t i=0;i<25;i++){
148 if (!fClusters[i]) continue;
149 if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
154 //_____________________________________________________________________________
155 void AliTRDseed::Update()
161 const Float_t kRatio = 0.8;
162 const Int_t kClmin = 6;
163 const Float_t kmaxtan = 2;
165 if (TMath::Abs(fYref[1])>kmaxtan) return; // too much inclined track
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
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
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
182 for (Int_t i=0;i<25;i++){
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]);
190 if (fN<kClmin) return;
191 Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
193 if (nz<=1) zouts[3]=0;
194 if (zouts[1]+zouts[3]<kClmin) return;
196 if (TMath::Abs(zouts[0]-zouts[2])>12.) zouts[3]=0; // z distance bigger than pad - length
198 Int_t breaktime = -1;
199 Bool_t mbefore = kFALSE;
201 Int_t counts[2]={0,0};
205 // find the break time allowing one chage on pad-rows with maximal numebr of accepted clusters
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]++;
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;
223 after = cumul[24][1]-cumul[i][1];
224 before = cumul[i][0];
225 if (after+before>maxcount) {
226 maxcount=after+before;
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];
237 if ( (allowedz[0]>allowedz[24] && fZref[1]<0) || (allowedz[0]<allowedz[24] && fZref[1]>0)){
239 // tracklet z-direction not in correspondance with track z direction
242 for (Int_t i=0;i<25;i++){
243 allowedz[i] = zouts[0]; //only longest taken
249 // cross pad -row tracklet - take the step change into account
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;
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];
274 EvaluateUni(fN2,yres2,mean,sigma,Int_t(fN2*kRatio-2));
275 if (sigma<sigmaexp*0.8) sigma=sigmaexp;
280 sumw=0; sumwx=0; sumwx2=0;
281 sumwy=0; sumwxy=0; sumwz=0;sumwxz=0;
286 for (Int_t i=0;i<25;i++){
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;
293 fMPads+=fClusters[i]->GetNPads();
295 if (fClusters[i]->GetNPads()>4) weight=0.5;
296 if (fClusters[i]->GetNPads()>5) weight=0.2;
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;
308 Float_t correction =0;
310 // tracklet on boundary
311 if (fMeanz<fZProb) correction = ycrosscor;
312 if (fMeanz>fZProb) correction = -ycrosscor;
314 Double_t det = sumw*sumwx2-sumwx*sumwx;
315 fYfitR[0] = (sumwx2*sumwy-sumwx*sumwxy)/det;
316 fYfitR[1] = (sumw*sumwxy-sumwx*sumwy)/det;
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;
324 fSigmaY2 = TMath::Sqrt(fSigmaY2/Float_t(fN2-2));
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];
339 //_____________________________________________________________________________
340 void AliTRDseed::UpdateUsed()
347 for (Int_t i=0;i<25;i++){
348 if (!fClusters[i]) continue;
349 if ((fClusters[i]->IsUsed())) fNUsed++;
354 //_____________________________________________________________________________
355 void AliTRDseed::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean
356 , Double_t &sigma, Int_t hh)
359 // Robust estimator in 1D case MI version
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
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);
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];
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]];
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){
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]];
402 Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
409 //_____________________________________________________________________________
410 Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
413 // Fit the Rieman tilt
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
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;
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];
437 Double_t x2 = cseed[iLayer].fX[itime]+cseed[iLayer].fX0; // global x
438 Double_t t = 1./(x2*x2+y*y);
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];
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);
453 Double_t rpolz0 = fitterT2.GetParameter(3);
454 Double_t rpolz1 = fitterT2.GetParameter(4);
456 // linear fitter - not possible to make boundaries
457 // non accept non possible z and dzdx combination
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;
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);
473 fitterT2.ReleaseParameter(3);
474 fitterT2.ReleaseParameter(4);
475 rpolz0 = fitterT2.GetParameter(3);
476 rpolz1 = fitterT2.GetParameter(4);
479 Double_t chi2TR = fitterT2.GetChisquare()/Float_t(npointsT);
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;
489 Double_t res2 = (x*params[0]+params[1]);
491 res2 = 1.-params[2]*params[0]+params[1]*params[1]-res2;
493 res2 = TMath::Sqrt(res2);
494 y = (1-res2)/params[0];
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.;
506 z = rpolz0+rpolz1*(x-xref2);
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;