]> 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
Updated RecParam calibration object (correct subversion)
[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// //
75bd7f81 22///////////////////////////////////////////////////////////////////////////////
23
ad670fba 24#include "TMath.h"
25#include "TLinearFitter.h"
26
4fad09c9 27#include "AliMathBase.h"
28
75bd7f81 29#include "AliTRDseed.h"
30#include "AliTRDcluster.h"
31#include "AliTRDtracker.h"
32
75bd7f81 33ClassImp(AliTRDseed)
34
35//_____________________________________________________________________________
ad670fba 36AliTRDseed::AliTRDseed()
37 :TObject()
38 ,fTilt(0)
39 ,fPadLength(0)
40 ,fX0(0)
41 ,fSigmaY(0)
42 ,fSigmaY2(0)
43 ,fMeanz(0)
44 ,fZProb(0)
45 ,fN(0)
46 ,fN2(0)
47 ,fNUsed(0)
48 ,fFreq(0)
49 ,fNChange(0)
50 ,fMPads(0)
51 ,fC(0)
52 ,fCC(0)
53 ,fChi2(1.0e10)
54 ,fChi2Z(1.0e10)
75bd7f81 55{
56 //
57 // Default constructor
58 //
59
60 for (Int_t i = 0; i < 25; i++) {
ad670fba 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
67 }
68
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
77 }
78
79}
80
81//_____________________________________________________________________________
82AliTRDseed::AliTRDseed(const AliTRDseed &s)
83 :TObject(s)
4fad09c9 84 ,fTilt(s.fTilt)
85 ,fPadLength(s.fPadLength)
86 ,fX0(s.fX0)
87 ,fSigmaY(s.fSigmaY)
88 ,fSigmaY2(s.fSigmaY2)
89 ,fMeanz(s.fMeanz)
90 ,fZProb(s.fZProb)
91 ,fN(s.fN)
92 ,fN2(s.fN2)
93 ,fNUsed(s.fNUsed)
94 ,fFreq(s.fFreq)
95 ,fNChange(s.fNChange)
96 ,fMPads(s.fMPads)
97 ,fC(s.fC)
98 ,fCC(s.fCC)
99 ,fChi2(s.fChi2)
100 ,fChi2Z(s.fChi2Z)
ad670fba 101{
102 //
103 // Copy constructor
104 //
105
106 for (Int_t i = 0; i < 25; i++) {
4fad09c9 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
75bd7f81 113 }
ad670fba 114
75bd7f81 115 for (Int_t i = 0; i < 2; i++) {
4fad09c9 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
75bd7f81 123 }
124
125}
126
127//_____________________________________________________________________________
128void AliTRDseed::Reset()
129{
130 //
131 // Reset seed
132 //
133
134 for (Int_t i = 0; i < 25; i++) {
ad670fba 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
75bd7f81 140 fUsable[i] = kFALSE;
141 }
ad670fba 142
75bd7f81 143 for (Int_t i = 0; i < 2; i++) {
ad670fba 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
75bd7f81 151 }
ad670fba 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
75bd7f81 156 fMPads = 0;
ad670fba 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
75bd7f81 161
162}
163
164//_____________________________________________________________________________
165void AliTRDseed::CookLabels()
166{
167 //
168 // Cook 2 labels for seed
169 //
170
171 Int_t labels[200];
172 Int_t out[200];
ad670fba 173 Int_t nlab = 0;
174
175 for (Int_t i = 0; i < 25; i++) {
75bd7f81 176 if (!fClusters[i]) continue;
ad670fba 177 for (Int_t ilab = 0; ilab < 3; ilab++) {
178 if (fClusters[i]->GetLabel(ilab) >= 0) {
75bd7f81 179 labels[nlab] = fClusters[i]->GetLabel(ilab);
180 nlab++;
181 }
182 }
183 }
ad670fba 184
75bd7f81 185 Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE);
186 fLabels[0] = out[0];
ad670fba 187 if ((nlab2 > 1) &&
188 (out[3] > 1)) {
189 fLabels[1] = out[2];
190 }
75bd7f81 191
192}
193
194//_____________________________________________________________________________
195void AliTRDseed::UseClusters()
196{
197 //
198 // Use clusters
199 //
200
ad670fba 201 for (Int_t i = 0; i < 25; i++) {
75bd7f81 202 if (!fClusters[i]) continue;
203 if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
204 }
205
206}
207
208//_____________________________________________________________________________
209void AliTRDseed::Update()
210{
211 //
212 // Update the seed.
213 //
214
215 const Float_t kRatio = 0.8;
216 const Int_t kClmin = 6;
217 const Float_t kmaxtan = 2;
218
ad670fba 219 if (TMath::Abs(fYref[1]) > kmaxtan) return; // Track inclined too much
220
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
223 fNChange = 0;
224
225 Double_t sumw;
226 Double_t sumwx;
227 Double_t sumwx2;
228 Double_t sumwy;
229 Double_t sumwxy;
230 Double_t sumwz;
231 Double_t sumwxz;
232
233 Int_t zints[25]; // Histograming of the z coordinate
234 // Get 1 and second max probable coodinates in z
235 Int_t zouts[50];
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
75bd7f81 239
240
ad670fba 241 fN = 0;
242 fN2 = 0;
243 for (Int_t i = 0; i < 25; i++) {
244 yres[i] = 10000.0;
75bd7f81 245 if (!fClusters[i]) continue;
ad670fba 246 yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i]; // Residual y
75bd7f81 247 zints[fN] = Int_t(fZ[i]);
248 fN++;
249 }
250
ad670fba 251 if (fN < kClmin) return;
75bd7f81 252 Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
253 fZProb = zouts[0];
ad670fba 254 if (nz <= 1) zouts[3] = 0;
255 if (zouts[1] + zouts[3] < kClmin) return;
75bd7f81 256
ad670fba 257 // Z distance bigger than pad - length
258 if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) {
259 zouts[3]=0;
260 }
75bd7f81 261
262 Int_t breaktime = -1;
263 Bool_t mbefore = kFALSE;
264 Int_t cumul[25][2];
ad670fba 265 Int_t counts[2] = { 0, 0 };
75bd7f81 266
ad670fba 267 if (zouts[3] >= 3) {
268
75bd7f81 269 //
ad670fba 270 // Find the break time allowing one chage on pad-rows
271 // with maximal numebr of accepted clusters
75bd7f81 272 //
ad670fba 273 fNChange = 1;
274 for (Int_t i = 0; i < 25; i++) {
75bd7f81 275 cumul[i][0] = counts[0];
276 cumul[i][1] = counts[1];
ad670fba 277 if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
278 if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
75bd7f81 279 }
ad670fba 280 Int_t maxcount = 0;
281 for (Int_t i = 0; i < 24; i++) {
282 Int_t after = cumul[24][0] - cumul[i][0];
75bd7f81 283 Int_t before = cumul[i][1];
ad670fba 284 if (after + before > maxcount) {
285 maxcount = after + before;
286 breaktime = i;
287 mbefore = kFALSE;
75bd7f81 288 }
ad670fba 289 after = cumul[24][1] - cumul[i][1];
75bd7f81 290 before = cumul[i][0];
ad670fba 291 if (after + before > maxcount) {
292 maxcount = after + before;
293 breaktime = i;
294 mbefore = kTRUE;
75bd7f81 295 }
296 }
ad670fba 297
298 breaktime -= 1;
299
75bd7f81 300 }
ad670fba 301
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];
75bd7f81 305 }
ad670fba 306
307 if (((allowedz[0] > allowedz[24]) && (fZref[1] < 0)) ||
308 ((allowedz[0] < allowedz[24]) && (fZref[1] > 0))) {
75bd7f81 309 //
ad670fba 310 // Tracklet z-direction not in correspondance with track z direction
75bd7f81 311 //
ad670fba 312 fNChange = 0;
313 for (Int_t i = 0; i < 25; i++) {
314 allowedz[i] = zouts[0]; // Only longest taken
75bd7f81 315 }
316 }
317
ad670fba 318 if (fNChange > 0) {
75bd7f81 319 //
ad670fba 320 // Cross pad -row tracklet - take the step change into account
75bd7f81 321 //
ad670fba 322 for (Int_t i = 0; i < 25; i++) {
75bd7f81 323 if (!fClusters[i]) continue;
ad670fba 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;
75bd7f81 329 }
330 }
331 }
332
333 Double_t yres2[25];
ad670fba 334 Double_t mean;
335 Double_t sigma;
336 for (Int_t i = 0; i < 25; i++) {
75bd7f81 337 if (!fClusters[i]) continue;
ad670fba 338 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
339 yres2[fN2] = yres[i];
75bd7f81 340 fN2++;
341 }
ad670fba 342 if (fN2 < kClmin) {
75bd7f81 343 fN2 = 0;
344 return;
345 }
4fad09c9 346 AliMathBase::EvaluateUni(fN2,yres2,mean,sigma,Int_t(fN2*kRatio-2));
ad670fba 347 if (sigma < sigmaexp * 0.8) {
348 sigma = sigmaexp;
349 }
75bd7f81 350 fSigmaY = sigma;
ad670fba 351
352 // Reset sums
353 sumw = 0;
354 sumwx = 0;
355 sumwx2 = 0;
356 sumwy = 0;
357 sumwxy = 0;
358 sumwz = 0;
359 sumwxz = 0;
360
361 fN2 = 0;
362 fMeanz = 0;
363 fMPads = 0;
364
365 for (Int_t i = 0; i < 25; i++) {
366
367 fUsable[i] = kFALSE;
75bd7f81 368 if (!fClusters[i]) continue;
ad670fba 369 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
370 if (TMath::Abs(yres[i] - mean) > 4.0 * sigma) continue;
75bd7f81 371 fUsable[i] = kTRUE;
372 fN2++;
ad670fba 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;
377
75bd7f81 378 Double_t x = fX[i];
ad670fba 379 sumw += weight;
380 sumwx += x * weight;
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;
386
75bd7f81 387 }
ad670fba 388
389 if (fN2 < kClmin){
75bd7f81 390 fN2 = 0;
391 return;
392 }
ad670fba 393 fMeanz = sumwz / sumw;
394 Float_t correction = 0;
395 if (fNChange > 0) {
396 // Tracklet on boundary
397 if (fMeanz < fZProb) correction = ycrosscor;
398 if (fMeanz > fZProb) correction = -ycrosscor;
75bd7f81 399 }
ad670fba 400
401 Double_t det = sumw * sumwx2 - sumwx * sumwx;
402 fYfitR[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det;
403 fYfitR[1] = (sumw * sumwxy - sumwx * sumwy) / det;
75bd7f81 404
ad670fba 405 fSigmaY2 = 0;
406 for (Int_t i = 0; i < 25; i++) {
75bd7f81 407 if (!fUsable[i]) continue;
ad670fba 408 Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
409 fSigmaY2 += delta*delta;
75bd7f81 410 }
ad670fba 411 fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2));
412
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];
75bd7f81 421
422 UpdateUsed();
423
424}
425
426//_____________________________________________________________________________
427void AliTRDseed::UpdateUsed()
428{
429 //
430 // Update used seed
431 //
432
433 fNUsed = 0;
ad670fba 434 for (Int_t i = 0; i < 25; i++) {
95a5eb17 435 if (!fClusters[i]) {
436 continue;
437 }
438 if ((fClusters[i]->IsUsed())) {
439 fNUsed++;
440 }
75bd7f81 441 }
442
443}
444
75bd7f81 445//_____________________________________________________________________________
446Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
447{
448 //
449 // Fit the Rieman tilt
450 //
451
ad670fba 452 // Fitting with tilting pads - kz not fixed
453 TLinearFitter fitterT2(4,"hyp4");
75bd7f81 454 fitterT2.StoreData(kTRUE);
ad670fba 455 Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
456
457 Int_t npointsT = 0;
75bd7f81 458 fitterT2.ClearPoints();
ad670fba 459
460 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
461
75bd7f81 462 if (!cseed[iLayer].IsOK()) continue;
463 Double_t tilt = cseed[iLayer].fTilt;
464
ad670fba 465 for (Int_t itime = 0; itime < 25; itime++) {
466
75bd7f81 467 if (!cseed[iLayer].fUsable[itime]) continue;
468 // x relative to the midle chamber
ad670fba 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];
75bd7f81 472
473 //
ad670fba 474 // Tilted rieman
75bd7f81 475 //
476 Double_t uvt[6];
ad670fba 477 Double_t x2 = cseed[iLayer].fX[itime] + cseed[iLayer].fX0; // Global x
478 Double_t t = 1.0 / (x2*x2 + y*y);
479 uvt[1] = t;
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];
484
485 Double_t error = 2.0 * uvt[1];
486 if (terror) {
487 error *= cseed[iLayer].fSigmaY;
488 }
489 else {
490 error *= 0.2; //Default error
491 }
75bd7f81 492 fitterT2.AddPoint(uvt,uvt[4],error);
493 npointsT++;
ad670fba 494
75bd7f81 495 }
ad670fba 496
75bd7f81 497 }
ad670fba 498
75bd7f81 499 fitterT2.Eval();
500 Double_t rpolz0 = fitterT2.GetParameter(3);
501 Double_t rpolz1 = fitterT2.GetParameter(4);
ad670fba 502
75bd7f81 503 //
ad670fba 504 // Linear fitter - not possible to make boundaries
75bd7f81 505 // non accept non possible z and dzdx combination
506 //
ad670fba 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) {
75bd7f81 512 acceptablez = kFALSE;
ad670fba 513 }
75bd7f81 514 }
515 }
ad670fba 516 if (!acceptablez) {
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;
75bd7f81 519 fitterT2.FixParameter(3,zmf);
520 fitterT2.FixParameter(4,dzmf);
521 fitterT2.Eval();
522 fitterT2.ReleaseParameter(3);
523 fitterT2.ReleaseParameter(4);
524 rpolz0 = fitterT2.GetParameter(3);
525 rpolz1 = fitterT2.GetParameter(4);
526 }
527
ad670fba 528 Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
75bd7f81 529 Double_t params[3];
ad670fba 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];
534
535 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
536
537 Double_t x = cseed[iLayer].fX0;
538 Double_t y = 0;
539 Double_t dy = 0;
540 Double_t z = 0;
541 Double_t dz = 0;
542
75bd7f81 543 // y
ad670fba 544 Double_t res2 = (x * params[0] + params[1]);
545 res2 *= res2;
546 res2 = 1.0 - params[2]*params[0] + params[1]*params[1] - res2;
547 if (res2 >= 0) {
75bd7f81 548 res2 = TMath::Sqrt(res2);
ad670fba 549 y = (1.0 - res2) / params[0];
75bd7f81 550 }
ad670fba 551
75bd7f81 552 //dy
ad670fba 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;
75bd7f81 559 dy = res;
560 }
561 }
ad670fba 562 z = rpolz0 + rpolz1 * (x - xref2);
75bd7f81 563 dz = rpolz1;
564 cseed[iLayer].fYref[0] = y;
565 cseed[iLayer].fYref[1] = dy;
566 cseed[iLayer].fZref[0] = z;
567 cseed[iLayer].fZref[1] = dz;
ad670fba 568 cseed[iLayer].fC = curvature;
75bd7f81 569
570 }
571
572 return chi2TR;
573
574}