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"
35 //_____________________________________________________________________________
36 AliTRDseed::AliTRDseed()
57 // Default constructor
60 for (Int_t i = 0; i < 25; i++) {
61 fX[i] = 0; // x position
62 fY[i] = 0; // y position
63 fZ[i] = 0; // z position
64 fIndexes[i] = 0; // Indexes
65 fClusters[i] = 0x0; // Clusters
66 fUsable[i] = 0; // Indication - usable cluster
69 for (Int_t i = 0; i < 2; i++) {
70 fYref[i] = 0; // Reference y
71 fZref[i] = 0; // Reference z
72 fYfit[i] = 0; // Y fit position +derivation
73 fYfitR[i] = 0; // Y fit position +derivation
74 fZfit[i] = 0; // Z fit position
75 fZfitR[i] = 0; // Z fit position
76 fLabels[i] = 0; // Labels
81 //_____________________________________________________________________________
82 AliTRDseed::AliTRDseed(const AliTRDseed &s)
85 ,fPadLength(s.fPadLength)
106 for (Int_t i = 0; i < 25; i++) {
107 fX[i] = s.fX[i]; // x position
108 fY[i] = s.fY[i]; // y position
109 fZ[i] = s.fZ[i]; // z position
110 fIndexes[i] = s.fIndexes[i]; // Indexes
111 fClusters[i] = s.fClusters[i]; // Clusters
112 fUsable[i] = s.fUsable[i]; // Indication - usable cluster
115 for (Int_t i = 0; i < 2; i++) {
116 fYref[i] = s.fYref[i]; // Reference y
117 fZref[i] = s.fZref[i]; // Reference z
118 fYfit[i] = s.fYfit[i]; // Y fit position +derivation
119 fYfitR[i] = s.fYfitR[i]; // Y fit position +derivation
120 fZfit[i] = s.fZfit[i]; // Z fit position
121 fZfitR[i] = s.fZfitR[i]; // Z fit position
122 fLabels[i] = s.fLabels[i]; // Labels
127 //_____________________________________________________________________________
128 void AliTRDseed::Reset()
134 for (Int_t i = 0; i < 25; i++) {
135 fX[i] = 0; // X position
136 fY[i] = 0; // Y position
137 fZ[i] = 0; // Z position
138 fIndexes[i] = 0; // Indexes
139 fClusters[i] = 0; // Clusters
143 for (Int_t i = 0; i < 2; i++) {
144 fYref[i] = 0; // Reference y
145 fZref[i] = 0; // Reference z
146 fYfit[i] = 0; // Y fit position +derivation
147 fYfitR[i] = 0; // Y fit position +derivation
148 fZfit[i] = 0; // Z fit position
149 fZfitR[i] = 0; // Z fit position
150 fLabels[i] = -1; // Labels
152 fSigmaY = 0; // "Robust" sigma in y
153 fSigmaY2 = 0; // "Robust" sigma in y
154 fMeanz = 0; // Mean vaue of z
155 fZProb = 0; // Max probbable z
157 fN = 0; // Number of associated clusters
158 fN2 = 0; // Number of not crossed
159 fNUsed = 0; // Number of used clusters
160 fNChange = 0; // Change z counter
164 //_____________________________________________________________________________
165 void AliTRDseed::CookLabels()
168 // Cook 2 labels for seed
175 for (Int_t i = 0; i < 25; i++) {
176 if (!fClusters[i]) continue;
177 for (Int_t ilab = 0; ilab < 3; ilab++) {
178 if (fClusters[i]->GetLabel(ilab) >= 0) {
179 labels[nlab] = fClusters[i]->GetLabel(ilab);
185 Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE);
194 //_____________________________________________________________________________
195 void AliTRDseed::UseClusters()
201 for (Int_t i = 0; i < 25; i++) {
202 if (!fClusters[i]) continue;
203 if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
208 //_____________________________________________________________________________
209 void AliTRDseed::Update()
215 const Float_t kRatio = 0.8;
216 const Int_t kClmin = 6;
217 const Float_t kmaxtan = 2;
219 if (TMath::Abs(fYref[1]) > kmaxtan) return; // Track inclined too much
221 Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
222 Float_t ycrosscor = fPadLength * fTilt * 0.5; // Y correction for crossing
233 Int_t zints[25]; // Histograming of the z coordinate
234 // Get 1 and second max probable coodinates in z
236 Float_t allowedz[25]; // Allowed z for given time bin
237 Float_t yres[25]; // Residuals from reference
238 Float_t anglecor = fTilt * fZref[1]; // Correction to the angle
243 for (Int_t i = 0; i < 25; i++) {
245 if (!fClusters[i]) continue;
246 yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i]; // Residual y
247 zints[fN] = Int_t(fZ[i]);
251 if (fN < kClmin) return;
252 Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
254 if (nz <= 1) zouts[3] = 0;
255 if (zouts[1] + zouts[3] < kClmin) return;
257 // Z distance bigger than pad - length
258 if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) {
262 Int_t breaktime = -1;
263 Bool_t mbefore = kFALSE;
265 Int_t counts[2] = { 0, 0 };
270 // Find the break time allowing one chage on pad-rows
271 // with maximal numebr of accepted clusters
274 for (Int_t i = 0; i < 25; i++) {
275 cumul[i][0] = counts[0];
276 cumul[i][1] = counts[1];
277 if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
278 if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
281 for (Int_t i = 0; i < 24; i++) {
282 Int_t after = cumul[24][0] - cumul[i][0];
283 Int_t before = cumul[i][1];
284 if (after + before > maxcount) {
285 maxcount = after + before;
289 after = cumul[24][1] - cumul[i][1];
290 before = cumul[i][0];
291 if (after + before > maxcount) {
292 maxcount = after + before;
302 for (Int_t i = 0; i < 25; i++) {
303 if (i > breaktime) allowedz[i] = mbefore ? zouts[2] : zouts[0];
304 if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
307 if (((allowedz[0] > allowedz[24]) && (fZref[1] < 0)) ||
308 ((allowedz[0] < allowedz[24]) && (fZref[1] > 0))) {
310 // Tracklet z-direction not in correspondance with track z direction
313 for (Int_t i = 0; i < 25; i++) {
314 allowedz[i] = zouts[0]; // Only longest taken
320 // Cross pad -row tracklet - take the step change into account
322 for (Int_t i = 0; i < 25; i++) {
323 if (!fClusters[i]) continue;
324 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
325 yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i]; // Residual y
326 if (TMath::Abs(fZ[i] - fZProb) > 2) {
327 if (fZ[i] > fZProb) yres[i] += fTilt * fPadLength;
328 if (fZ[i] < fZProb) yres[i] -= fTilt * fPadLength;
336 for (Int_t i = 0; i < 25; i++) {
337 if (!fClusters[i]) continue;
338 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
339 yres2[fN2] = yres[i];
346 AliMathBase::EvaluateUni(fN2,yres2,mean,sigma,Int_t(fN2*kRatio-2));
347 if (sigma < sigmaexp * 0.8) {
365 for (Int_t i = 0; i < 25; i++) {
368 if (!fClusters[i]) continue;
369 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
370 if (TMath::Abs(yres[i] - mean) > 4.0 * sigma) continue;
373 fMPads += fClusters[i]->GetNPads();
374 Float_t weight = 1.0;
375 if (fClusters[i]->GetNPads() > 4) weight = 0.5;
376 if (fClusters[i]->GetNPads() > 5) weight = 0.2;
381 sumwx2 += x*x * weight;
382 sumwy += weight * yres[i];
383 sumwxy += weight * (yres[i]) * x;
384 sumwz += weight * fZ[i];
385 sumwxz += weight * fZ[i] * x;
393 fMeanz = sumwz / sumw;
394 Float_t correction = 0;
396 // Tracklet on boundary
397 if (fMeanz < fZProb) correction = ycrosscor;
398 if (fMeanz > fZProb) correction = -ycrosscor;
401 Double_t det = sumw * sumwx2 - sumwx * sumwx;
402 fYfitR[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det;
403 fYfitR[1] = (sumw * sumwxy - sumwx * sumwy) / det;
406 for (Int_t i = 0; i < 25; i++) {
407 if (!fUsable[i]) continue;
408 Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
409 fSigmaY2 += delta*delta;
411 fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2));
413 fZfitR[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
414 fZfitR[1] = (sumw * sumwxz - sumwx * sumwz) / det;
415 fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
416 fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det;
417 fYfitR[0] += fYref[0] + correction;
418 fYfitR[1] += fYref[1];
419 fYfit[0] = fYfitR[0];
420 fYfit[1] = fYfitR[1];
426 //_____________________________________________________________________________
427 void AliTRDseed::UpdateUsed()
434 for (Int_t i = 0; i < 25; i++) {
438 if ((fClusters[i]->IsUsed())) {
445 //_____________________________________________________________________________
446 Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
449 // Fit the Rieman tilt
452 // Fitting with tilting pads - kz not fixed
453 TLinearFitter fitterT2(4,"hyp4");
454 fitterT2.StoreData(kTRUE);
455 Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
458 fitterT2.ClearPoints();
460 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
462 if (!cseed[iLayer].IsOK()) continue;
463 Double_t tilt = cseed[iLayer].fTilt;
465 for (Int_t itime = 0; itime < 25; itime++) {
467 if (!cseed[iLayer].fUsable[itime]) continue;
468 // x relative to the midle chamber
469 Double_t x = cseed[iLayer].fX[itime] + cseed[iLayer].fX0 - xref2;
470 Double_t y = cseed[iLayer].fY[itime];
471 Double_t z = cseed[iLayer].fZ[itime];
477 Double_t x2 = cseed[iLayer].fX[itime] + cseed[iLayer].fX0; // Global x
478 Double_t t = 1.0 / (x2*x2 + y*y);
480 uvt[0] = 2.0 * x2 * uvt[1];
481 uvt[2] = 2.0 * tilt * uvt[1];
482 uvt[3] = 2.0 * tilt *uvt[1] * x;
483 uvt[4] = 2.0 * (y + tilt * z) * uvt[1];
485 Double_t error = 2.0 * uvt[1];
487 error *= cseed[iLayer].fSigmaY;
490 error *= 0.2; //Default error
492 fitterT2.AddPoint(uvt,uvt[4],error);
500 Double_t rpolz0 = fitterT2.GetParameter(3);
501 Double_t rpolz1 = fitterT2.GetParameter(4);
504 // Linear fitter - not possible to make boundaries
505 // non accept non possible z and dzdx combination
507 Bool_t acceptablez = kTRUE;
508 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
509 if (cseed[iLayer].IsOK()) {
510 Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2);
511 if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) {
512 acceptablez = kFALSE;
517 Double_t zmf = cseed[2].fZref[0] + cseed[2].fZref[1] * (xref2 - cseed[2].fX0);
518 Double_t dzmf = (cseed[2].fZref[1] + cseed[3].fZref[1]) * 0.5;
519 fitterT2.FixParameter(3,zmf);
520 fitterT2.FixParameter(4,dzmf);
522 fitterT2.ReleaseParameter(3);
523 fitterT2.ReleaseParameter(4);
524 rpolz0 = fitterT2.GetParameter(3);
525 rpolz1 = fitterT2.GetParameter(4);
528 Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
530 params[0] = fitterT2.GetParameter(0);
531 params[1] = fitterT2.GetParameter(1);
532 params[2] = fitterT2.GetParameter(2);
533 Double_t curvature = 1.0 + params[1] * params[1] - params[2] * params[0];
535 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
537 Double_t x = cseed[iLayer].fX0;
544 Double_t res2 = (x * params[0] + params[1]);
546 res2 = 1.0 - params[2]*params[0] + params[1]*params[1] - res2;
548 res2 = TMath::Sqrt(res2);
549 y = (1.0 - res2) / params[0];
553 Double_t x0 = -params[1] / params[0];
554 if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) {
555 Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1);
556 if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) {
557 Double_t res = (x - x0) / TMath::Sqrt(1.0 / (rm1*rm1) - (x-x0)*(x-x0));
558 if (params[0] < 0) res *= -1.0;
562 z = rpolz0 + rpolz1 * (x - xref2);
564 cseed[iLayer].fYref[0] = y;
565 cseed[iLayer].fYref[1] = dy;
566 cseed[iLayer].fZref[0] = z;
567 cseed[iLayer].fZref[1] = dz;
568 cseed[iLayer].fC = curvature;