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 //
22 ///////////////////////////////////////////////////////////////////////////////
25 #include "TLinearFitter.h"
27 #include "AliMathBase.h"
29 #include "AliTRDseed.h"
30 #include "AliTRDcluster.h"
31 #include "AliTRDtracker.h"
32 #include "AliTRDtrackerV1.h"
36 //_____________________________________________________________________________
37 AliTRDseed::AliTRDseed()
58 // Default constructor
61 for (Int_t i = 0; i < knTimebins; i++) {
62 fX[i] = 0; // x position
63 fY[i] = 0; // y position
64 fZ[i] = 0; // z position
65 fIndexes[i] = 0; // Indexes
66 fClusters[i] = NULL; // Clusters
67 fUsable[i] = 0; // Indication - usable cluster
70 for (Int_t i = 0; i < 2; i++) {
71 fYref[i] = 0; // Reference y
72 fZref[i] = 0; // Reference z
73 fYfit[i] = 0; // Y fit position +derivation
74 fYfitR[i] = 0; // Y fit position +derivation
75 fZfit[i] = 0; // Z fit position
76 fZfitR[i] = 0; // Z fit position
77 fLabels[i] = 0; // Labels
82 //_____________________________________________________________________________
83 AliTRDseed::AliTRDseed(const AliTRDseed &s)
86 ,fPadLength(s.fPadLength)
107 for (Int_t i = 0; i < knTimebins; i++) {
108 fX[i] = s.fX[i]; // x position
109 fY[i] = s.fY[i]; // y position
110 fZ[i] = s.fZ[i]; // z position
111 fIndexes[i] = s.fIndexes[i]; // Indexes
112 fClusters[i] = s.fClusters[i]; // Clusters
113 fUsable[i] = s.fUsable[i]; // Indication - usable cluster
116 for (Int_t i = 0; i < 2; i++) {
117 fYref[i] = s.fYref[i]; // Reference y
118 fZref[i] = s.fZref[i]; // Reference z
119 fYfit[i] = s.fYfit[i]; // Y fit position +derivation
120 fYfitR[i] = s.fYfitR[i]; // Y fit position +derivation
121 fZfit[i] = s.fZfit[i]; // Z fit position
122 fZfitR[i] = s.fZfitR[i]; // Z fit position
123 fLabels[i] = s.fLabels[i]; // Labels
128 //_____________________________________________________________________________
129 void AliTRDseed::Copy(TObject &o) const
131 //printf("AliTRDseed::Copy()\n");
133 AliTRDseed &seed = (AliTRDseed &)o;
136 seed.fPadLength = fPadLength;
138 seed.fSigmaY = fSigmaY;
139 seed.fSigmaY2 = fSigmaY2;
140 seed.fMeanz = fMeanz;
141 seed.fZProb = fZProb;
144 seed.fNUsed = fNUsed;
146 seed.fNChange = fNChange;
147 seed.fMPads = fMPads;
151 seed.fChi2Z = fChi2Z;
152 for (Int_t i = 0; i < knTimebins; i++) {
156 seed.fIndexes[i] = fIndexes[i];
157 seed.fClusters[i] = fClusters[i];
158 seed.fUsable[i] = fUsable[i];
161 for (Int_t i = 0; i < 2; i++) {
162 seed.fYref[i] = fYref[i];
163 seed.fZref[i] = fZref[i];
164 seed.fYfit[i] = fYfit[i];
165 seed.fYfitR[i] = fYfitR[i];
166 seed.fZfit[i] = fZfit[i];
167 seed.fZfitR[i] = fZfitR[i];
168 seed.fLabels[i] = fLabels[i];
175 //_____________________________________________________________________________
176 void AliTRDseed::Reset()
182 for (Int_t i = 0; i < knTimebins; i++) {
183 fX[i] = 0; // X position
184 fY[i] = 0; // Y position
185 fZ[i] = 0; // Z position
186 fIndexes[i] = 0; // Indexes
187 fClusters[i] = NULL; // Clusters
191 for (Int_t i = 0; i < 2; i++) {
192 fYref[i] = 0; // Reference y
193 fZref[i] = 0; // Reference z
194 fYfit[i] = 0; // Y fit position +derivation
195 fYfitR[i] = 0; // Y fit position +derivation
196 fZfit[i] = 0; // Z fit position
197 fZfitR[i] = 0; // Z fit position
198 fLabels[i] = -1; // Labels
200 fSigmaY = 0; // "Robust" sigma in y
201 fSigmaY2 = 0; // "Robust" sigma in y
202 fMeanz = 0; // Mean vaue of z
203 fZProb = 0; // Max probbable z
205 fN = 0; // Number of associated clusters
206 fN2 = 0; // Number of not crossed
207 fNUsed = 0; // Number of used clusters
208 fNChange = 0; // Change z counter
212 //_____________________________________________________________________________
213 void AliTRDseed::CookLabels()
216 // Cook 2 labels for seed
223 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
224 if (!fClusters[i]) continue;
225 for (Int_t ilab = 0; ilab < 3; ilab++) {
226 if (fClusters[i]->GetLabel(ilab) >= 0) {
227 labels[nlab] = fClusters[i]->GetLabel(ilab);
233 Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE);
242 //_____________________________________________________________________________
243 void AliTRDseed::UseClusters()
249 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
250 if (!fClusters[i]) continue;
251 if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
256 //_____________________________________________________________________________
257 void AliTRDseed::Update()
265 // linear fit on the y direction with respect to the reference direction.
266 // The residuals for each x (x = xc - x0) are deduced from:
268 // the tilting correction is written :
269 // y = yc + h*(zc-zt) (2)
270 // yt = y0+dy/dx*x (3)
271 // zt = z0+dz/dx*x (4)
272 // from (1),(2),(3) and (4)
273 // dy = yc - y0 - (dy/dx + h*dz/dx)*x + h*(zc-z0)
274 // the last term introduces the correction on y direction due to tilting pads. There are 2 ways to account for this:
275 // 1. use tilting correction for calculating the y
276 // 2. neglect tilting correction here and account for it in the error parametrization of the tracklet.
278 const Float_t kRatio = 0.8;
279 const Int_t kClmin = 5;
280 const Float_t kmaxtan = 2;
282 if (TMath::Abs(fYref[1]) > kmaxtan){
283 //printf("Exit: Abs(fYref[1]) = %3.3f, kmaxtan = %3.3f\n", TMath::Abs(fYref[1]), kmaxtan);
284 return; // Track inclined too much
287 Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
288 Float_t ycrosscor = fPadLength * fTilt * 0.5; // Y correction for crossing
299 // Buffering: Leave it constant fot Performance issues
300 Int_t zints[knTimebins]; // Histograming of the z coordinate
301 // Get 1 and second max probable coodinates in z
302 Int_t zouts[2*knTimebins];
303 Float_t allowedz[knTimebins]; // Allowed z for given time bin
304 Float_t yres[knTimebins]; // Residuals from reference
305 //Float_t anglecor = fTilt * fZref[1]; // Correction to the angle
310 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
312 if (!fClusters[i]) continue;
313 if(!fClusters[i]->IsInChamber()) continue;
315 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + fTilt*(fZ[i] - fZref[0]);
316 yres[i] = fY[i] - fTilt*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
317 zints[fN] = Int_t(fZ[i]);
322 //printf("Exit fN < kClmin: fN = %d\n", fN);
325 Int_t nz = AliTRDtracker::Freq(fN, zints, zouts, kFALSE);
327 if (nz <= 1) zouts[3] = 0;
328 if (zouts[1] + zouts[3] < kClmin) {
329 //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]);
333 // Z distance bigger than pad - length
334 if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) zouts[3] = 0;
336 Int_t breaktime = -1;
337 Bool_t mbefore = kFALSE;
338 Int_t cumul[knTimebins][2];
339 Int_t counts[2] = { 0, 0 };
344 // Find the break time allowing one chage on pad-rows
345 // with maximal number of accepted clusters
348 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
349 cumul[i][0] = counts[0];
350 cumul[i][1] = counts[1];
351 if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
352 if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
355 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
356 Int_t after = cumul[AliTRDtrackerV1::GetNTimeBins()][0] - cumul[i][0];
357 Int_t before = cumul[i][1];
358 if (after + before > maxcount) {
359 maxcount = after + before;
363 after = cumul[AliTRDtrackerV1::GetNTimeBins()-1][1] - cumul[i][1];
364 before = cumul[i][0];
365 if (after + before > maxcount) {
366 maxcount = after + before;
376 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
377 if (i > breaktime) allowedz[i] = mbefore ? zouts[2] : zouts[0];
378 if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
381 if (((allowedz[0] > allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] < 0)) ||
382 ((allowedz[0] < allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] > 0))) {
384 // Tracklet z-direction not in correspondance with track z direction
387 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
388 allowedz[i] = zouts[0]; // Only longest taken
394 // Cross pad -row tracklet - take the step change into account
396 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
397 if (!fClusters[i]) continue;
398 if(!fClusters[i]->IsInChamber()) continue;
399 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
401 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] /*+ fTilt*(fZ[i] - fZref[0])*/;
402 yres[i] = fY[i] - fTilt*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
403 /* if (TMath::Abs(fZ[i] - fZProb) > 2) {
404 if (fZ[i] > fZProb) yres[i] += fTilt * fPadLength;
405 if (fZ[i] < fZProb) yres[i] -= fTilt * fPadLength;
410 Double_t yres2[knTimebins];
413 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
414 if (!fClusters[i]) continue;
415 if(!fClusters[i]->IsInChamber()) continue;
416 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
417 yres2[fN2] = yres[i];
421 //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
425 AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
426 if (sigma < sigmaexp * 0.8) {
444 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
447 if (!fClusters[i]) continue;
448 if (!fClusters[i]->IsInChamber()) continue;
449 if (TMath::Abs(fZ[i] - allowedz[i]) > 2){fClusters[i] = NULL; continue;}
450 if (TMath::Abs(yres[i] - mean) > 4.0 * sigma){fClusters[i] = NULL; continue;}
453 fMPads += fClusters[i]->GetNPads();
454 Float_t weight = 1.0;
455 if (fClusters[i]->GetNPads() > 4) weight = 0.5;
456 if (fClusters[i]->GetNPads() > 5) weight = 0.2;
460 //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]);
464 sumwx2 += x*x * weight;
465 sumwy += weight * yres[i];
466 sumwxy += weight * (yres[i]) * x;
467 sumwz += weight * fZ[i];
468 sumwxz += weight * fZ[i] * x;
473 //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
477 fMeanz = sumwz / sumw;
478 Float_t correction = 0;
480 // Tracklet on boundary
481 if (fMeanz < fZProb) correction = ycrosscor;
482 if (fMeanz > fZProb) correction = -ycrosscor;
485 Double_t det = sumw * sumwx2 - sumwx * sumwx;
486 fYfitR[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det;
487 fYfitR[1] = (sumw * sumwxy - sumwx * sumwy) / det;
490 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
491 if (!fUsable[i]) continue;
492 Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
493 fSigmaY2 += delta*delta;
495 fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2));
496 // TEMPORARY UNTIL covariance properly calculated
497 fSigmaY2 = TMath::Max(fSigmaY2, Float_t(.1));
499 fZfitR[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
500 fZfitR[1] = (sumw * sumwxz - sumwx * sumwz) / det;
501 fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
502 fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det;
503 // fYfitR[0] += fYref[0] + correction;
504 // fYfitR[1] += fYref[1];
505 fYfit[0] = fYfitR[0];
506 fYfit[1] = -fYfitR[1];
508 //printf("y0 = %7.3f tgy = %7.3f z0 = %7.3f tgz = %7.3f \n", fYfitR[0], fYfitR[1], fZfitR[0], fZfitR[1]);
514 //_____________________________________________________________________________
515 void AliTRDseed::UpdateUsed()
522 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
523 if (!fClusters[i]) continue;
524 if(!fUsable[i]) continue;
525 if ((fClusters[i]->IsUsed())) fNUsed++;
530 //_____________________________________________________________________________
531 Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
534 // Fit the Rieman tilt
537 // Fitting with tilting pads - kz not fixed
538 TLinearFitter fitterT2(4,"hyp4");
539 fitterT2.StoreData(kTRUE);
541 Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
544 fitterT2.ClearPoints();
546 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
548 if (!cseed[iLayer].IsOK()) continue;
549 Double_t tilt = cseed[iLayer].fTilt;
551 for (Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins()+1; itime++) {
553 if (!cseed[iLayer].fUsable[itime]) continue;
554 // x relative to the midle chamber
555 Double_t x = cseed[iLayer].fX[itime] + cseed[iLayer].fX0 - xref2;
556 Double_t y = cseed[iLayer].fY[itime];
557 Double_t z = cseed[iLayer].fZ[itime];
563 Double_t x2 = cseed[iLayer].fX[itime] + cseed[iLayer].fX0; // Global x
564 Double_t t = 1.0 / (x2*x2 + y*y);
566 uvt[0] = 2.0 * x2 * uvt[1];
567 uvt[2] = 2.0 * tilt * uvt[1];
568 uvt[3] = 2.0 * tilt *uvt[1] * x;
569 uvt[4] = 2.0 * (y + tilt * z) * uvt[1];
571 Double_t error = 2.0 * uvt[1];
573 error *= cseed[iLayer].fSigmaY;
576 error *= 0.2; //Default error
578 fitterT2.AddPoint(uvt,uvt[4],error);
586 Double_t rpolz0 = fitterT2.GetParameter(3);
587 Double_t rpolz1 = fitterT2.GetParameter(4);
590 // Linear fitter - not possible to make boundaries
591 // non accept non possible z and dzdx combination
593 Bool_t acceptablez = kTRUE;
594 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
595 if (!cseed[iLayer].IsOK()) continue;
596 Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2);
597 if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) acceptablez = kFALSE;
600 Double_t zmf = cseed[2].fZref[0] + cseed[2].fZref[1] * (xref2 - cseed[2].fX0);
601 Double_t dzmf = (cseed[2].fZref[1] + cseed[3].fZref[1]) * 0.5;
602 fitterT2.FixParameter(3,zmf);
603 fitterT2.FixParameter(4,dzmf);
605 fitterT2.ReleaseParameter(3);
606 fitterT2.ReleaseParameter(4);
607 rpolz0 = fitterT2.GetParameter(3);
608 rpolz1 = fitterT2.GetParameter(4);
611 Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
613 params[0] = fitterT2.GetParameter(0);
614 params[1] = fitterT2.GetParameter(1);
615 params[2] = fitterT2.GetParameter(2);
616 Double_t curvature = 1.0 + params[1] * params[1] - params[2] * params[0];
619 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
620 Double_t x = cseed[iLayer].fX0;
627 Double_t res2 = (x * params[0] + params[1]);
629 res2 = 1.0 - params[2]*params[0] + params[1]*params[1] - res2;
631 res2 = TMath::Sqrt(res2);
632 y = (1.0 - res2) / params[0];
636 Double_t x0 = -params[1] / params[0];
637 if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) {
638 Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1);
639 if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) {
640 Double_t res = (x - x0) / TMath::Sqrt((1./rm1-(x-x0))*(1./rm1+(x-x0)));
641 if (params[0] < 0) res *= -1.0;
645 z = rpolz0 + rpolz1 * (x - xref2);
647 cseed[iLayer].fYref[0] = y;
648 cseed[iLayer].fYref[1] = dy;
649 cseed[iLayer].fZref[0] = z;
650 cseed[iLayer].fZref[1] = dz;
651 cseed[iLayer].fC = curvature;