First version of the parallel TPC tracking (M.Ivanov)
[u/mrichter/AliRoot.git] / TPC / AliTPCtrackerMI.cxx
CommitLineData
1c53abe2 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
17
18
19
20
21/*
22 AliTPC parallel tracker -
23 How to use? -
24 run AliTPCFindClusters.C macro - clusters neccessary for tracker are founded
25 run AliTPCFindTracksMI.C macro - to find tracks
26 tracks are written to AliTPCtracks.root file
27 for comparison also seeds are written to the same file - to special branch
28*/
29
30//-------------------------------------------------------
31// Implementation of the TPC tracker
32//
33// Origin: Marian Ivanov Marian.Ivanov@cern.ch
34//
35//-------------------------------------------------------
36
37#include <TObjArray.h>
38#include <TFile.h>
39#include <TTree.h>
40#include <iostream.h>
41
42#include "AliTPCtrackerMI.h"
43#include "AliTPCclusterMI.h"
44#include "AliTPCParam.h"
45#include "AliTPCClustersRow.h"
46#include "AliComplexCluster.h"
47#include "TStopwatch.h"
48
49
50ClassImp(AliTPCseed)
51
52AliTPCclusterTracks::AliTPCclusterTracks(){
53 // class for storing overlaping info
54 fTrackIndex[0]=-1;
55 fTrackIndex[1]=-1;
56 fTrackIndex[2]=-1;
57 fDistance[0]=1000;
58 fDistance[1]=1000;
59 fDistance[2]=1000;
60}
61
62
63
64
65
66Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, AliTPCclusterMI* c, Double_t chi2, UInt_t i){
67
68 Int_t sec=(i&0xff000000)>>24;
69 Int_t row = (i&0x00ff0000)>>16;
70 track->fRow=(i&0x00ff0000)>>16;
71 track->fSector = sec;
72 // Int_t index = i&0xFFFF;
73 if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow();
74 track->fClusterIndex[track->fRow] = i;
75 track->fFirstPoint = row;
76 //
77
78 AliTPCTrackPoint *trpoint =track->GetTrackPoint(track->fRow);
79 Float_t angle2 = track->GetSnp()*track->GetSnp();
80 angle2 = TMath::Sqrt(angle2/(1-angle2));
81 //
82 //SET NEW Track Point
83 //
84 if (c!=0){
85 //if we have a cluster
86 trpoint->GetCPoint().SetY(c->GetY());
87 trpoint->GetCPoint().SetZ(c->GetZ());
88 //
89 trpoint->GetCPoint().SetSigmaY(c->GetSigmaY2()/(track->fCurrentSigmaY*track->fCurrentSigmaY));
90 trpoint->GetCPoint().SetSigmaZ(c->GetSigmaZ2()/(track->fCurrentSigmaZ*track->fCurrentSigmaZ));
91 //
92 trpoint->GetCPoint().SetType(c->GetType());
93 trpoint->GetCPoint().SetQ(c->GetQ());
94 trpoint->GetCPoint().SetMax(c->GetMax());
95 //
96 trpoint->GetCPoint().SetErrY(TMath::Sqrt(track->fErrorY2));
97 trpoint->GetCPoint().SetErrZ(TMath::Sqrt(track->fErrorZ2));
98 //
99 }
100 trpoint->GetTPoint().SetX(track->GetX());
101 trpoint->GetTPoint().SetY(track->GetY());
102 trpoint->GetTPoint().SetZ(track->GetZ());
103 //
104 trpoint->GetTPoint().SetAngleY(angle2);
105 trpoint->GetTPoint().SetAngleZ(track->GetTgl());
106
107
108
109 if (chi2>10){
110 // printf("suspicious chi2 %f\n",chi2);
111 }
112 // if (track->fIsSeeding){
113 track->fErrorY2 *= 1.2;
114 track->fErrorY2 += 0.0064;
115 track->fErrorZ2 *= 1.2;
116 track->fErrorY2 += 0.005;
117
118 //}
119
120 return track->Update(c,chi2,i);
121
122}
123//_____________________________________________________________________________
124AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par, Int_t eventn):
125AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
126{
127 //---------------------------------------------------------------------
128 // The main TPC tracker constructor
129 //---------------------------------------------------------------------
130 fInnerSec=new AliTPCSector[fkNIS];
131 fOuterSec=new AliTPCSector[fkNOS];
132
133 Int_t i;
134 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
135 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
136
137 fN=0; fSectors=0;
138
139 fClustersArray.Setup(par);
140 fClustersArray.SetClusterType("AliTPCclusterMI");
141
142 char cname[100];
143 if (eventn==-1) {
144 sprintf(cname,"TreeC_TPC");
145 }
146 else {
147 sprintf(cname,"TreeC_TPC_%d",eventn);
148 }
149
150 fClustersArray.ConnectTree(cname);
151
152 fEventN = eventn;
153 fSeeds=0;
154 fNtracks = 0;
155 fParam = par;
156}
157
158//_____________________________________________________________________________
159AliTPCtrackerMI::~AliTPCtrackerMI() {
160 //------------------------------------------------------------------
161 // TPC tracker destructor
162 //------------------------------------------------------------------
163 delete[] fInnerSec;
164 delete[] fOuterSec;
165 if (fSeeds) {
166 fSeeds->Delete();
167 delete fSeeds;
168 }
169}
170
171
172Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
173 //
174 //
175 Float_t snoise2;
176 Float_t z = fParam->GetZLength()-TMath::Abs(seed->GetZ());
177
178 //cluster "quality"
179 Float_t rsigmay = 1;
180 Int_t ctype = 0;
181
182 //standard if we don't have cluster - take MIP
183 const Float_t chmip = 50.;
184 Float_t amp = chmip/0.3;
185 Float_t nel;
186 Float_t nprim;
187 if (cl){
188 amp = cl->GetQ();
189 rsigmay = cl->GetSigmaY2()/(seed->fCurrentSigmaY*seed->fCurrentSigmaY);
190 ctype = cl->GetType();
191 }
192
193
194 Float_t landau=2 ; //landau fluctuation part
195 Float_t gg=2; // gg fluctuation part
196 Float_t padlength= fSectors->GetPadPitchLength(seed->GetX());
197
198
199 if (fSectors==fInnerSec){
200 snoise2 = 0.0004;
201 nel = 0.268*amp;
202 nprim = 0.155*amp;
203 gg = (2+0.0002*amp)/nel;
204 landau = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
205 if (landau>1) landau=1;
206 }
207 else {
208 snoise2 = 0.0004;
209 nel = 0.3*amp;
210 nprim = 0.133*amp;
211 gg = (2+0.0002*amp)/nel;
212 landau = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
213 if (landau>1) landau=1;
214 }
215
216
217 Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
218 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
219 angle2 = angle2/(1-angle2);
220 Float_t angular = landau*angle2*padlength*padlength/12.;
221 Float_t res = sdiff + angular;
222
223
224 if ((ctype==0) && (fSectors ==fOuterSec))
225 res *= 0.78 +TMath::Exp(7.4*(rsigmay-1.2));
226
227 if ((ctype==0) && (fSectors ==fInnerSec))
228 res *= 0.72 +TMath::Exp(3.36*(rsigmay-1.2));
229
230
231 if ((ctype>0))
232 res*= TMath::Power((rsigmay+0.5),1.5)+0.0064;
233
234 if (ctype<0)
235 res*=2.4; // overestimate error 2 times
236
237 res+= snoise2;
238
239 if (res<2*snoise2)
240 res = 2*snoise2;
241
242 seed->SetErrorY2(res);
243 return res;
244}
245
246
247
248
249Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
250 //
251 //
252 Float_t snoise2;
253 Float_t z = fParam->GetZLength()-TMath::Abs(seed->GetZ());
254 //signal quality
255 Float_t rsigmaz=1;
256 Int_t ctype =0;
257
258 const Float_t chmip = 50.;
259 Float_t amp = chmip/0.3;
260 Float_t nel;
261 Float_t nprim;
262 if (cl){
263 amp = cl->GetQ();
264 rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ*seed->fCurrentSigmaZ);
265 ctype = cl->GetType();
266 }
267
268 //
269 Float_t landau=2 ; //landau fluctuation part
270 Float_t gg=2; // gg fluctuation part
271 Float_t padlength= fSectors->GetPadPitchLength(seed->GetX());
272
273 if (fSectors==fInnerSec){
274 snoise2 = 0.0004;
275 nel = 0.268*amp;
276 nprim = 0.155*amp;
277 gg = (2+0.0002*amp)/nel;
278 landau = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
279 if (landau>1) landau=1;
280 }
281 else {
282 snoise2 = 0.0004;
283 nel = 0.3*amp;
284 nprim = 0.133*amp;
285 gg = (2+0.0002*amp)/nel;
286 landau = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
287 if (landau>1) landau=1;
288 }
289 Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
290
291 Float_t angle = seed->GetTgl();
292 Float_t angular = landau*angle*angle*padlength*padlength/12.;
293 Float_t res = sdiff + angular;
294
295 if ((ctype==0) && (fSectors ==fOuterSec))
296 res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
297
298 if ((ctype==0) && (fSectors ==fInnerSec))
299 res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
300 if ((ctype>0))
301 res*= TMath::Power(rsigmaz+0.5,1.5)+0.0064; //0.31+0.147*ctype;
302 if (ctype<0)
303 res*=1.3;
304 if ((ctype<0) &&amp<70)
305 res*=1.3;
306
307 res += snoise2;
308 if (res<2*snoise2)
309 res = 2*snoise2;
310
311 seed->SetErrorZ2(res);
312 return res;
313}
314
315
316
317Int_t AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
318{
319 //-----------------------------------------------------------------
320 // This function find proloncation of a track to a reference plane x=xk.
321 // doesn't change internal state of the track
322 //-----------------------------------------------------------------
323
324 Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
325 // Double_t y1=fP0, z1=fP1;
326 Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
327 Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
328
329 y = fP0;
330 z = fP1;
331 y += dx*(c1+c2)/(r1+r2);
332 z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
333 return 0;
334}
335
336
337//_____________________________________________________________________________
338Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const
339{
340 //-----------------------------------------------------------------
341 // This function calculates a predicted chi2 increment.
342 //-----------------------------------------------------------------
343 //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
344 Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
345 r00+=fC00; r01+=fC10; r11+=fC11;
346
347 Double_t det=r00*r11 - r01*r01;
348 if (TMath::Abs(det) < 1.e-10) {
349 Int_t n=GetNumberOfClusters();
350 if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
351 return 1e10;
352 }
353 Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
354
355 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
356
357 return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
358}
359
360
361//_________________________________________________________________________________________
362
363
364Int_t AliTPCseed::Compare(const TObject *o) const {
365 //-----------------------------------------------------------------
366 // This function compares tracks according to the sector - for given sector according z
367 //-----------------------------------------------------------------
368 AliTPCseed *t=(AliTPCseed*)o;
369 if (t->fSector>fSector) return -1;
370 if (t->fSector<fSector) return 1;
371
372 Double_t z2 = t->GetZ();
373 Double_t z1 = GetZ();
374 if (z2>z1) return 1;
375 if (z2<z1) return -1;
376 return 0;
377}
378
379
380
381
382
383
384//_____________________________________________________________________________
385Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t index) {
386 //-----------------------------------------------------------------
387 // This function associates a cluster with this track.
388 //-----------------------------------------------------------------
389 // Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
390 //Double_t r00=sigmay2, r01=0., r11=sigmaz2;
391 Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
392
393 r00+=fC00; r01+=fC10; r11+=fC11;
394 Double_t det=r00*r11 - r01*r01;
395 Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
396
397 Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
398 Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
399 Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
400 Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
401 Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
402
403 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
404 Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
405 if (TMath::Abs(cur*fX-eta) >= 0.99999) {
406 Int_t n=GetNumberOfClusters();
407 if (n>4) cerr<<n<<" AliTPCtrack warning: Filtering failed !\n";
408 return 0;
409 }
410
411 fP0 += k00*dy + k01*dz;
412 fP1 += k10*dy + k11*dz;
413 fP2 = eta;
414 fP3 += k30*dy + k31*dz;
415 fP4 = cur;
416
417 Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
418 Double_t c12=fC21, c13=fC31, c14=fC41;
419
420 fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
421 fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
422 fC40-=k00*c04+k01*c14;
423
424 fC11-=k10*c01+k11*fC11;
425 fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
426 fC41-=k10*c04+k11*c14;
427
428 fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
429 fC42-=k20*c04+k21*c14;
430
431 fC33-=k30*c03+k31*c13;
432 fC43-=k40*c03+k41*c13;
433
434 fC44-=k40*c04+k41*c14;
435
436 Int_t n=GetNumberOfClusters();
437 fIndex[n]=index;
438 SetNumberOfClusters(n+1);
439 SetChi2(GetChi2()+chisq);
440
441 return 1;
442}
443
444
445
446//_____________________________________________________________________________
447Double_t AliTPCtrackerMI::f1(Double_t x1,Double_t y1,
448 Double_t x2,Double_t y2,
449 Double_t x3,Double_t y3)
450{
451 //-----------------------------------------------------------------
452 // Initial approximation of the track curvature
453 //-----------------------------------------------------------------
454 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
455 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
456 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
457 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
458 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
459
460 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
461
462 return -xr*yr/sqrt(xr*xr+yr*yr);
463}
464
465
466//_____________________________________________________________________________
467Double_t AliTPCtrackerMI::f2(Double_t x1,Double_t y1,
468 Double_t x2,Double_t y2,
469 Double_t x3,Double_t y3)
470{
471 //-----------------------------------------------------------------
472 // Initial approximation of the track curvature times center of curvature
473 //-----------------------------------------------------------------
474 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
475 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
476 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
477 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
478 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
479
480 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
481
482 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
483}
484
485//_____________________________________________________________________________
486Double_t AliTPCtrackerMI::f3(Double_t x1,Double_t y1,
487 Double_t x2,Double_t y2,
488 Double_t z1,Double_t z2)
489{
490 //-----------------------------------------------------------------
491 // Initial approximation of the tangent of the track dip angle
492 //-----------------------------------------------------------------
493 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
494}
495
496
497void AliTPCtrackerMI::LoadClusters()
498{
499 //
500 // load clusters to the memory
501 Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
502 for (Int_t i=0; i<j; i++) {
503 fClustersArray.LoadEntry(i);
504 }
505}
506
507void AliTPCtrackerMI::UnloadClusters()
508{
509 //
510 // load clusters to the memory
511 Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
512 for (Int_t i=0; i<j; i++) {
513 fClustersArray.ClearSegment(i);
514 }
515}
516
517
518
519//_____________________________________________________________________________
520void AliTPCtrackerMI::LoadOuterSectors() {
521 //-----------------------------------------------------------------
522 // This function fills outer TPC sectors with clusters.
523 //-----------------------------------------------------------------
524 UInt_t index;
525 Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
526 for (Int_t i=0; i<j; i++) {
527 // AliSegmentID *s=fClustersArray.LoadEntry(i);
528 AliSegmentID *s= fClustersArray.At(i);
529 Int_t sec,row;
530 AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
531 par->AdjustSectorRow(s->GetID(),sec,row);
532 if (sec<fkNIS*2) continue;
533 AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
534 Int_t ncl=clrow->GetArray()->GetEntriesFast();
535 while (ncl--) {
536 AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl];
537 index=(((sec<<8)+row)<<16)+ncl;
538 fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index);
539 }
540 }
541 fN=fkNOS;
542 fSectors=fOuterSec;
543}
544
545
546//_____________________________________________________________________________
547void AliTPCtrackerMI::LoadInnerSectors() {
548 //-----------------------------------------------------------------
549 // This function fills inner TPC sectors with clusters.
550 //-----------------------------------------------------------------
551 UInt_t index;
552 Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
553 for (Int_t i=0; i<j; i++) {
554 // AliSegmentID *s=fClustersArray.LoadEntry(i);
555 AliSegmentID *s=fClustersArray.At(i);
556 Int_t sec,row;
557 AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
558 par->AdjustSectorRow(s->GetID(),sec,row);
559 if (sec>=fkNIS*2) continue;
560 AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
561 Int_t ncl=clrow->GetArray()->GetEntriesFast();
562 while (ncl--) {
563 AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl];
564 index=(((sec<<8)+row)<<16)+ncl;
565 fInnerSec[sec%fkNIS][row].InsertCluster(c,index);
566 }
567 }
568
569 fN=fkNIS;
570 fSectors=fInnerSec;
571}
572
573Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
574 //-----------------------------------------------------------------
575 // This function tries to find a track prolongation to next pad row
576 //-----------------------------------------------------------------
577 // Double_t xt=t.GetX();
578 // Int_t row = fSectors->GetRowNumber(xt)-1;
579 // if (row < nr) return 1; // don't prolongate if not information until now -
580 //
581 Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
582 if (!t.PropagateTo(x)) return 0;
583 // update current
584 t.fCurrentSigmaY = GetSigmaY(&t);
585 t.fCurrentSigmaZ = GetSigmaZ(&t);
586 //
587 AliTPCclusterMI *cl=0;
588 UInt_t index=0;
589 const AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
590 Double_t sy2=ErrY2(&t)*2;
591 Double_t sz2=ErrZ2(&t)*2;
592
593
594 Double_t roady =3.*sqrt(t.GetSigmaY2() + sy2);
595 Double_t roadz = 3 *sqrt(t.GetSigmaZ2() + sz2);
596 Double_t y=t.GetY(), z=t.GetZ();
597
598 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
599 t.fInDead = kTRUE;
600 return 0;
601 }
602 else
603 {
604 t.fNFoundable++;
605 }
606 //calculate
607 Float_t maxdistance = roady*roady + roadz*roadz;
608 if (krow) {
609 for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
610 AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
611 if (c->GetZ() > z+roadz) break;
612 if ( (c->GetY()-y) > roady ) continue;
613 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
614 if (maxdistance>distance) {
615 maxdistance = distance;
616 cl=c;
617 index=krow.GetIndex(i);
618 }
619 }
620 }
621 if (cl) {
622 //Double_t sy2=
623 ErrY2(&t,cl);
624 //Double_t sz2=
625 ErrZ2(&t,cl);
626 Double_t chi2= t.GetPredictedChi2(cl);
627 UpdateTrack(&t,cl,chi2,index);
628 } else {
629 if (y > ymax) {
630 t.fRelativeSector= (t.fRelativeSector+1) % fN;
631 if (!t.Rotate(fSectors->GetAlpha()))
632 return 0;
633 } else if (y <-ymax) {
634 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
635 if (!t.Rotate(-fSectors->GetAlpha()))
636 return 0;
637 }
638 }
639 return 1;
640}
641
642
643Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,Int_t trindex, Int_t nr) {
644 //-----------------------------------------------------------------
645 // This function tries to find a track prolongation to next pad row
646 //-----------------------------------------------------------------
647 t.fCurrentCluster = 0;
648 t.fCurrentClusterIndex1 = 0;
649 t.fCurrentClusterIndex2 = 0;
650
651 Double_t xt=t.GetX();
652 Int_t row = fSectors->GetRowNumber(xt)-1;
653 if (row < nr) return 1; // don't prolongate if not information until now -
654 Double_t x=fSectors->GetX(nr);
655 if (!t.PropagateTo(x)) return 0;
656 // update current
657 t.fCurrentSigmaY = GetSigmaY(&t);
658 t.fCurrentSigmaZ = GetSigmaZ(&t);
659
660 AliTPCclusterMI *cl=0;
661 UInt_t index=0;
662 AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
663 //
664 Double_t y=t.GetY(), z=t.GetZ();
665 Double_t roady = 3.* TMath::Sqrt(t.GetSigmaY2() + t.fCurrentSigmaY*t.fCurrentSigmaY);
666 Double_t roadz = 3.* TMath::Sqrt(t.GetSigmaZ2() + t.fCurrentSigmaZ*t.fCurrentSigmaZ);
667 //
668
669 Float_t maxdistance = 1000000;
670 if (krow) {
671 for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
672 AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
673 if (c->GetZ() > z+roadz) break;
674 if (TMath::Abs(c->GetY()-y)>roady) continue;
675
676 //krow.UpdateClusterTrack(i,trindex,&t);
677
678 Float_t dy2 = (c->GetY()- t.GetY());
679 dy2*=dy2;
680 Float_t dz2 = (c->GetZ()- t.GetZ());
681 dz2*=dz2;
682 //
683 Float_t distance = dy2+dz2;
684 //
685 if (distance > maxdistance) continue;
686 maxdistance = distance;
687 cl=c;
688 index=i;
689 }
690 }
691 t.fCurrentCluster = cl;
692 t.fCurrentClusterIndex1 = krow.GetIndex(index);
693 t.fCurrentClusterIndex2 = index;
694 return 1;
695}
696
697
698Int_t AliTPCtrackerMI::FollowToNextCluster(Int_t trindex, Int_t nr) {
699 //-----------------------------------------------------------------
700 // This function tries to find a track prolongation to next pad row
701 //-----------------------------------------------------------------
702 AliTPCseed & t = *((AliTPCseed*)(fSeeds->At(trindex)));
703 AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
704 // Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
705 Double_t y=t.GetY();
706 Double_t ymax=fSectors->GetMaxY(nr);
707
708 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
709 t.fInDead = kTRUE;
710 return 0;
711 }
712 else
713 {
714 t.fNFoundable++;
715 }
716
717 if (t.fCurrentCluster) {
718 // Float_t l=fSectors->GetPadPitchLength();
719 // AliTPCclusterTracks * cltrack = krow.GetClusterTracks(t.fCurrentClusterIndex1);
720
721 Double_t sy2=ErrY2(&t,t.fCurrentCluster);
722 Double_t sz2=ErrZ2(&t,t.fCurrentCluster);
723
724
725 Double_t sdistancey = TMath::Sqrt(sy2+t.GetSigmaY2());
726 Double_t sdistancez = TMath::Sqrt(sz2+t.GetSigmaZ2());
727
728 Double_t rdistancey = TMath::Abs(t.fCurrentCluster->GetY()-t.GetY());
729 Double_t rdistancez = TMath::Abs(t.fCurrentCluster->GetZ()-t.GetZ());
730
731 Double_t rdistance = TMath::Sqrt(TMath::Power(rdistancey/sdistancey,2)+TMath::Power(rdistancez/sdistancez,2));
732
733
734 // printf("\t%f\t%f\t%f\n",rdistancey/sdistancey,rdistancez/sdistancez,rdistance);
735 if (rdistance>4) return 0;
736
737 if ((rdistancey/sdistancey>2.5 || rdistancez/sdistancez>2.5) && t.fCurrentCluster->GetType()==0)
738 return 0; //suspisiouce - will be changed
739
740 if ((rdistancey/sdistancey>2. || rdistancez/sdistancez>2.0) && t.fCurrentCluster->GetType()>0)
741 // strict cut on overlaped cluster
742 return 0; //suspisiouce - will be changed
743
744 if ( (rdistancey/sdistancey>1. || rdistancez/sdistancez>2.5 ||t.fCurrentCluster->GetQ()<70 )
745 && t.fCurrentCluster->GetType()<0)
746 return 0;
747
748 // t.SetSampledEdx(0.3*t.fCurrentCluster->GetQ()/l,t.GetNumberOfClusters(), GetSigmaY(&t), GetSigmaZ(&t));
749 UpdateTrack(&t,t.fCurrentCluster,t.GetPredictedChi2(t.fCurrentCluster),t.fCurrentClusterIndex1);
750
751 } else {
752 if (y > ymax) {
753 t.fRelativeSector= (t.fRelativeSector+1) % fN;
754 if (!t.Rotate(fSectors->GetAlpha()))
755 return 0;
756 } else if (y <-ymax) {
757 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
758 if (!t.Rotate(-fSectors->GetAlpha()))
759 return 0;
760 }
761 }
762 return 1;
763}
764
765
766
767
768/*
769Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t step)
770{
771 //-----------------------------------------------------------------
772 // fast prolongation mathod -
773 // don't update track only after step clusters
774 //-----------------------------------------------------------------
775 Double_t xt=t.GetX();
776 //
777 Double_t alpha=t.GetAlpha();
778 alpha =- fSectors->GetAlphaShift();
779 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
780 if (alpha < 0. ) alpha += 2.*TMath::Pi();
781 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
782 Int_t row0 = fSectors->GetRowNumber(xt);
783 Double_t x = fSectors->GetX(row0);
784 Double_t ymax = fSectors->GetMaxY(row0);
785 //
786 Double_t sy2=ErrY2(&t)*2;
787 Double_t sz2=ErrZ2(&t)*2;
788 Double_t roady =3.*sqrt(t.GetSigmaY2() + sy2);
789 Double_t roadz = 3 *sqrt(t.GetSigmaZ2() + sz2);
790 Float_t maxdistance = roady*roady + roadz*roadz;
791 t.fCurrentSigmaY = GetSigmaY(&t);
792 t.fCurrentSigmaZ = GetSigmaZ(&t);
793 //
794 Int_t nclusters = 0;
795 Double_t y;
796 Double_t z;
797 Double_t yy[200]; //track prolongation
798 Double_t zz[200];
799 Double_t cy[200]; // founded cluster position
800 Double_t cz[200];
801 Double_t sy[200]; // founded cluster error
802 Double_t sz[200];
803 Bool_t hitted[200]; // indication of cluster presence
804 //
805
806 //
807 for (Int_t drow = step; drow>=0; drow--) {
808 Int_t row = row0-drow;
809 if (row<0) break;
810 Double_t x = fSectors->GetX(row);
811 Double_t ymax = fSectors->GetMaxY(row);
812 t.GetProlongation(x,y,z);
813 yy[drow] =y;
814 zz[drow] =z;
815 const AliTPCRow &krow=fSectors[t.fRelativeSector][row];
816 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
817 t.fInDead = kTRUE;
818 break;
819 }
820 else
821 {
822 t.fNFoundable++;
823 }
824
825 //find nearest cluster
826 AliTPCclusterMI *cl= 0;
827 if (krow) {
828 for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
829 AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
830 if (c->GetZ() > z+roadz) break;
831 if ( (c->GetY()-y) > roady ) continue;
832 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
833 if (maxdistance>distance) {
834 maxdistance = distance;
835 cl=c;
836 // index=krow.GetIndex(i);
837 }
838 }
839 } // end of seearch
840 //update cluster information
841 if (cl){
842 cy[drow] = cl->GetY();
843 cz[drow] = cl->GetZ();
844 sy[drow] = ErrY2(&t,cl);
845 sz[drow] = ErrZ2(&t,cl);
846 hitted[drow] = kTRUE;
847 nclusters++;
848 }
849 else
850 hitted[drow] = kFALSE;
851 }
852 //if we have information - update track
853 if (nclusters>0){
854 Float_t sumyw0 = 0;
855 Float_t sumzw0 = 0;
856 Float_t sumyw = 0;
857 Float_t sumzw = 0;
858 Float_t sumrow = 0;
859 for (Int_t i=0;i<step;i++)
860 if (hitted[i]){
861 sumyw0+= 1/sy[i];
862 sumzw0+= 1/sz[i];
863 //
864 sumyw+= (cy[i]-yy[i])/sy[i];
865 sumzw+= (cz[i]-zz[i])/sz[i];
866 sumrow+= i;
867 }
868 Float_t dy = sumyw/sumyw0;
869 Float_t dz = sumzw/sumzw0;
870 Float_t mrow = sumrow/nclusters+row0;
871 Float_t x = fSectors->GetX(mrow);
872 t.PropagateTo(x);
873 AliTPCclusterMI cvirtual;
874 cvirtual.SetZ(dz+t.GetZ());
875 cvirtual.SetY(dy+t.GetY());
876 t.SetErrorY2(1.2*t.fErrorY2/TMath::Sqrt(Float_t(nclusters)));
877 t.SetErrorZ2(1.2*t.fErrorZ2/TMath::Sqrt(Float_t(nclusters)));
878 Float_t chi2 = t.GetPredictedChi2(&cvirtual);
879 t.Update(&cvirtual,chi2,0);
880 Int_t ncl = t.GetNumberOfClusters();
881 ncl = ncl-1+nclusters;
882 t.SetN(ncl);
883 }
884 return 1;
885}
886*/
887
888//_____________________________________________________________________________
889Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf) {
890 //-----------------------------------------------------------------
891 // This function tries to find a track prolongation.
892 //-----------------------------------------------------------------
893 Double_t xt=t.GetX();
894 //
895 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
896 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
897 if (alpha < 0. ) alpha += 2.*TMath::Pi();
898 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
899
900 for (Int_t nr=fSectors->GetRowNumber(xt)-1; nr>=rf; nr--) {
901
902 if (FollowToNext(t,nr)==0) {
903 }
904 }
905 return 1;
906}
907
908
909//_____________________________________________________________________________
910Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
911 //-----------------------------------------------------------------
912 // This function tries to find a track prolongation.
913 //-----------------------------------------------------------------
914 Double_t xt=t.GetX();
915 //
916 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
917 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
918 if (alpha < 0. ) alpha += 2.*TMath::Pi();
919 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
920
921 for (Int_t nr=fSectors->GetRowNumber(xt)+1; nr<=rf; nr++) {
922 if (t.GetSnp()>0.2)
923 FollowToNext(t,nr);
924 }
925 return 1;
926}
927
928
929
930
931
932Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
933{
934 //
935 //
936 sum1=0;
937 sum2=0;
938 Int_t sum=0;
939 if (s1->fSector!=s2->fSector) return 0;
940 //
941 Float_t dz2 =(s1->GetZ() - s2->GetZ());
942 dz2*=dz2;
943 Float_t dy2 =(s1->GetY() - s2->GetY());
944 dy2*=dy2;
945 Float_t distance = TMath::Sqrt(dz2+dy2);
946 if (distance>5.) return 0; // if there are far away - not overlap - to reduce combinatorics
947
948 for (Int_t i=0;i<160;i++){
949 if (s1->fClusterIndex[i]>0) sum1++;
950 if (s2->fClusterIndex[i]>0) sum2++;
951 if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) {
952 sum++;
953 }
954 }
955
956 Float_t summin = TMath::Min(sum1,sum2);
957 Float_t ratio = sum/Float_t(summin);
958 return ratio;
959}
960
961void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
962{
963 //
964 //
965 if (s1->fSector!=s2->fSector) return;
966 //
967 Float_t dz2 =(s1->GetZ() - s2->GetZ());
968 dz2*=dz2;
969 Float_t dy2 =(s1->GetY() - s2->GetY());
970 dy2*=dy2;
971 Float_t distance = TMath::Sqrt(dz2+dy2);
972 if (distance>15.) return ; // if there are far away - not overlap - to reduce combinatorics
973 //trpoint = new (pointarray[track->fRow]) AliTPCTrackPoint;
974 // TClonesArray &pointarray1 = *(s1->fPoints);
975 //TClonesArray &pointarray2 = *(s2->fPoints);
976 //
977 for (Int_t i=0;i<160;i++){
978 if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) {
979 // AliTPCTrackPoint *p1 = (AliTPCTrackPoint *)(pointarray1.UncheckedAt(i));
980 //AliTPCTrackPoint *p2 = (AliTPCTrackPoint *)(pointarray2.UncheckedAt(i));
981 AliTPCTrackPoint *p1 = s1->GetTrackPoint(i);
982 AliTPCTrackPoint *p2 = s2->GetTrackPoint(i);;
983 p1->fIsShared = kTRUE;
984 p2->fIsShared = kTRUE;
985 }
986 }
987}
988
989
990
991
992void AliTPCtrackerMI::RemoveOverlap(TObjArray * arr, Float_t factor, Int_t removalindex , Bool_t shared=kFALSE){
993
994
995
996 //
997 // remove overlap - used removal factor - removal index stored in the track
998 arr->Sort(); // sorting according z
999 arr->Expand(arr->GetEntries());
1000 Int_t nseed=arr->GetEntriesFast();
1001 // printf("seeds \t%p \t%d\n",arr, nseed);
1002 // arr->Expand(arr->GetEntries()); //remove 0 pointers
1003 nseed = arr->GetEntriesFast();
1004 Int_t removed = 0;
1005 for (Int_t i=0; i<nseed; i++) {
1006 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1007 if (!pt) {
1008 continue;
1009 }
1010 if (!(pt->IsActive())) continue;
1011 for (Int_t j=i+1; j<nseed; j++){
1012 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1013 if ((pt2) && pt2->IsActive())
1014 if (pt->fSector == pt2->fSector)
1015 if (TMath::Abs(pt2->GetZ()-pt->GetZ())<2){
1016 Int_t sum1,sum2;
1017 Float_t ratio = OverlapFactor(pt,pt2,sum1,sum2);
1018 if (ratio>factor){
1019 // if (pt->GetChi2()<pt2->GetChi2()) pt2->Desactivate(removalindex); // arr->RemoveAt(j);
1020 Float_t ratio2 = (pt->GetChi2()*sum2)/(pt2->GetChi2()*sum1);
1021 Float_t ratio3 = Float_t(sum1-sum2)/Float_t(sum1+sum2);
1022 removed++;
1023 if (TMath::Abs(ratio3)>0.025){ // if much more points
1024 if (sum1>sum2) pt2->Desactivate(removalindex);
1025 else {
1026 pt->Desactivate(removalindex); // arr->RemoveAt(i);
1027 break;
1028 }
1029 }
1030 else{ //decide on mean chi2
1031 if (ratio2<1)
1032 pt2->Desactivate(removalindex);
1033 else {
1034 pt->Desactivate(removalindex); // arr->RemoveAt(i);
1035 break;
1036 }
1037 }
1038
1039 } // if suspicious ratio
1040 }
1041 else
1042 break;
1043 }
1044 }
1045 // printf("removed\t%d\n",removed);
1046 Int_t good =0;
1047 for (Int_t i=0; i<nseed; i++) {
1048 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1049 if (!pt) break;
1050 if (pt->GetNumberOfClusters() < pt->fNFoundable*0.5) {
1051 //desactivate tracks with small number of points
1052 // printf("%d\t%d\t%f\n", pt->GetNumberOfClusters(), pt->fNFoundable,pt->GetNumberOfClusters()/Float_t(pt->fNFoundable));
1053 pt->Desactivate(10); //desactivate - small muber of points
1054 }
1055 if (!(pt->IsActive())) continue;
1056 good++;
1057 }
1058
1059
1060 if (shared)
1061 for (Int_t i=0; i<nseed; i++) {
1062 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1063 if (!pt) continue;
1064 if (!(pt->IsActive())) continue;
1065 for (Int_t j=i+1; j<nseed; j++){
1066 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1067 if ((pt2) && pt2->IsActive()) {
1068 if ( TMath::Abs(pt->fSector-pt2->fSector)>1) break;
1069 SignShared(pt,pt2);
1070 }
1071 }
1072 }
1073 fNtracks = good;
1074 printf("\n*****\nNumber of good tracks after overlap removal\t%d\n",fNtracks);
1075
1076
1077}
1078void AliTPCtrackerMI::MakeSeedsAll()
1079{
1080 if (fSeeds == 0) fSeeds = new TObjArray;
1081 TObjArray * arr;
1082 for (Int_t sec=0;sec<fkNOS;sec+=3){
1083 arr = MakeSeedsSectors(sec,sec+3);
1084 Int_t nseed = arr->GetEntriesFast();
1085 for (Int_t i=0;i<nseed;i++)
1086 fSeeds->AddLast(arr->RemoveAt(i));
1087 }
1088 // fSeeds = MakeSeedsSectors(0,fkNOS);
1089}
1090
1091TObjArray * AliTPCtrackerMI::MakeSeedsSectors(Int_t sec1, Int_t sec2)
1092{
1093 //
1094 // loop over all sectors and make seed
1095 //find track seeds
1096 TStopwatch timer;
1097 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
1098 Int_t nrows=nlow+nup;
1099 Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
1100 // if (fSeeds==0) fSeeds = new TObjArray;
1101 TObjArray * arr = new TObjArray;
1102
1103 for (Int_t sec=sec1; sec<sec2;sec++){
1104 MakeSeeds(arr, sec, nup-1, nup-1-gap);
1105 MakeSeeds(arr, sec, nup-1-shift, nup-1-shift-gap);
1106 }
1107 Int_t nseed=arr->GetEntriesFast();
1108 gap=Int_t(0.3*nrows);
1109 // continue seeds
1110 Int_t i;
1111 for (i=0; i<nseed; i++) {
1112 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
1113 if (!pt) continue;
1114 if (FollowProlongation(t,nup-gap)) {
1115 pt->fIsSeeding =kFALSE;
1116 continue;
1117 }
1118 delete arr->RemoveAt(i);
1119 }
1120 //
1121 //remove seeds which overlaps
1122 RemoveOverlap(arr,0.6,1);
1123 //delete seeds - which were sign
1124 nseed=arr->GetEntriesFast();
1125 for (i=0; i<nseed; i++) {
1126 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1127 //, &t=*pt;
1128 if (!pt) continue;
1129 if ((pt->IsActive()) && pt->GetNumberOfClusters() > pt->fNFoundable*0.5 ) {
1130 //FollowBackProlongation(t,nup);
1131 continue;
1132 }
1133 delete arr->RemoveAt(i);
1134 }
1135 return arr;
1136}
1137
1138
1139
1140//_____________________________________________________________________________
1141void AliTPCtrackerMI::MakeSeeds(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2) {
1142 //-----------------------------------------------------------------
1143 // This function creates track seeds.
1144 //-----------------------------------------------------------------
1145 // if (fSeeds==0) fSeeds=new TObjArray(15000);
1146
1147 Double_t x[5], c[15];
1148
1149 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
1150 Double_t cs=cos(alpha), sn=sin(alpha);
1151
1152 Double_t x1 =fOuterSec->GetX(i1);
1153 Double_t xx2=fOuterSec->GetX(i2);
1154 //
1155 // for (Int_t ns=0; ns<fkNOS; ns++)
1156 Int_t ns =sec;
1157 {
1158 Int_t nl=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
1159 Int_t nm=fOuterSec[ns][i2];
1160 Int_t nu=fOuterSec[(ns+1)%fkNOS][i2];
1161 const AliTPCRow& kr1=fOuterSec[ns][i1];
1162 AliTPCRow& kr21 = fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
1163 AliTPCRow& kr22 = fOuterSec[(ns)%fkNOS][i2];
1164 AliTPCRow& kr23 = fOuterSec[(ns+1)%fkNOS][i2];
1165
1166 for (Int_t is=0; is < kr1; is++) {
1167 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
1168 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
1169
1170 Float_t anglez = (z1-z3)/(x1-x3);
1171 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
1172
1173 for (Int_t js=0; js < nl+nm+nu; js++) {
1174 const AliTPCclusterMI *kcl;
1175 Double_t x2, y2, z2;
1176 if (js<nl) {
1177 if (js==0) {
1178 js = kr21.Find(extraz-15.);
1179 if (js>=nl) continue;
1180 }
1181 kcl=kr21[js];
1182 z2=kcl->GetZ();
1183 if ((extraz-z2)>10) continue;
1184 if ((extraz-z2)<-10) {
1185 js = nl;
1186 continue;
1187 }
1188 y2=kcl->GetY();
1189 x2= xx2*cs+y2*sn;
1190 y2=-xx2*sn+y2*cs;
1191 } else
1192 if (js<nl+nm) {
1193 if (js==nl) {
1194 js = nl+kr22.Find(extraz-15.);
1195 if (js>=nl+nm) continue;
1196 }
1197 kcl=kr22[js-nl];
1198 z2=kcl->GetZ();
1199 if ((extraz-z2)>10) continue;
1200 if ((extraz-z2)<-10) {
1201 js = nl+nm;
1202 continue;
1203 }
1204 x2=xx2; y2=kcl->GetY();
1205 } else {
1206 //const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][i2];
1207 if (js==nl+nm) {
1208 js = nl+nm+kr23.Find(extraz-15.);
1209 if (js>=nl+nm+nu) break;
1210 }
1211 kcl=kr23[js-nl-nm];
1212 z2=kcl->GetZ();
1213 if ((extraz-z2)>10) continue;
1214 if ((extraz-z2)<-10) {
1215 break;
1216 }
1217 y2=kcl->GetY();
1218 x2=xx2*cs-y2*sn;
1219 y2=xx2*sn+y2*cs;
1220 }
1221
1222 Double_t zz=z1 - anglez*(x1-x2);
1223 if (TMath::Abs(zz-z2)>10.) continue;
1224
1225 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
1226 if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
1227
1228 x[0]=y1;
1229 x[1]=z1;
1230 x[4]=f1(x1,y1,x2,y2,x3,y3);
1231 if (TMath::Abs(x[4]) >= 0.0066) continue;
1232 x[2]=f2(x1,y1,x2,y2,x3,y3);
1233 //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
1234 x[3]=f3(x1,y1,x2,y2,z1,z2);
1235 if (TMath::Abs(x[3]) > 1.2) continue;
1236 Double_t a=asin(x[2]);
1237 Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
1238 if (TMath::Abs(zv-z3)>10.) continue;
1239
1240 Double_t sy1=kr1[is]->GetSigmaY2()*2, sz1=kr1[is]->GetSigmaZ2()*4;
1241 Double_t sy2=kcl->GetSigmaY2()*2, sz2=kcl->GetSigmaZ2()*4;
1242 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
1243 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
1244 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
1245
1246 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
1247 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
1248 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
1249 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
1250 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1251 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
1252 Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
1253 Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
1254 Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
1255 Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
1256
1257 c[0]=sy1;
1258 c[1]=0.; c[2]=sz1;
1259 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
1260 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
1261 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
1262 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
1263 c[13]=f30*sy1*f40+f32*sy2*f42;
1264 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
1265
1266 UInt_t index=kr1.GetIndex(is);
1267 AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
1268 track->fIsSeeding = kTRUE;
1269 Int_t rc=FollowProlongation(*track, i2);
1270 //FollowProlongationFast(*track, 5);
1271 //FollowProlongationFast(*track, 5);
1272 //FollowProlongationFast(*track, 5);
1273 //FollowProlongationFast(*track, 5);
1274 //Int_t rc = 1;
1275
1276 track->fLastPoint = i1; // first cluster in track position
1277 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/4 || track->GetNumberOfClusters() < track->fNFoundable/2. ) delete track;
1278 else arr->AddLast(track);
1279 }
1280 }
1281 }
1282}
1283
1284//_____________________________________________________________________________
1285Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
1286 //-----------------------------------------------------------------
1287 // This function reades track seeds.
1288 //-----------------------------------------------------------------
1289 TDirectory *savedir=gDirectory;
1290
1291 TFile *in=(TFile*)inp;
1292 if (!in->IsOpen()) {
1293 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
1294 return 1;
1295 }
1296
1297 in->cd();
1298 TTree *seedTree=(TTree*)in->Get("Seeds");
1299 if (!seedTree) {
1300 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
1301 cerr<<"can't get a tree with track seeds !\n";
1302 return 2;
1303 }
1304 AliTPCtrack *seed=new AliTPCtrack;
1305 seedTree->SetBranchAddress("tracks",&seed);
1306
1307 if (fSeeds==0) fSeeds=new TObjArray(15000);
1308
1309 Int_t n=(Int_t)seedTree->GetEntries();
1310 for (Int_t i=0; i<n; i++) {
1311 seedTree->GetEvent(i);
1312 fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
1313 }
1314
1315 delete seed;
1316 delete seedTree;
1317 savedir->cd();
1318 return 0;
1319}
1320
1321//_____________________________________________________________________________
1322Int_t AliTPCtrackerMI::Clusters2Tracks(const TFile *inp, TFile *out) {
1323 //-----------------------------------------------------------------
1324 // This is a track finder.
1325 //-----------------------------------------------------------------
1326 TDirectory *savedir=gDirectory;
1327
1328 if (inp) {
1329 TFile *in=(TFile*)inp;
1330 if (!in->IsOpen()) {
1331 cerr<<"AliTPCtrackerMI::Clusters2Tracks(): input file is not open !\n";
1332 return 1;
1333 }
1334 }
1335
1336 if (!out->IsOpen()) {
1337 cerr<<"AliTPCtrackerMI::Clusters2Tracks(): output file is not open !\n";
1338 return 2;
1339 }
1340
1341 out->cd();
1342
1343 char tname[100];
1344 sprintf(tname,"TreeT_TPC_%d",fEventN);
1345 TTree tracktree(tname,"Tree with TPC tracks");
1346 TTree seedtree("Seeds","Seeds");
1347 AliTPCtrack *iotrack=0;
1348 AliTPCseed *ioseed=0;
1349 tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
1350 TStopwatch timer;
1351
1352 printf("Loading clusters \n");
1353 LoadClusters();
1354 printf("Time for loading clusters: \t");timer.Print();timer.Start();
1355
1356 printf("Loading outer sectors\n");
1357 LoadOuterSectors();
1358 printf("Time for loading outer sectors: \t");timer.Print();timer.Start();
1359
1360 printf("Loading inner sectors\n");
1361 LoadInnerSectors();
1362 printf("Time for loading inner sectors: \t");timer.Print();timer.Start();
1363 fSectors = fOuterSec;
1364 fN=fkNOS;
1365
1366
1367 //find track seeds
1368 MakeSeedsAll();
1369 printf("Time for seeding: \t"); timer.Print();timer.Start();
1370 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
1371 Int_t nrows=nlow+nup;
1372
1373 Int_t gap=Int_t(0.3*nrows);
1374 Int_t i;
1375 Int_t nseed=fSeeds->GetEntriesFast();
1376 // outer sectors parallel tracking
1377 ParallelTracking(fSectors->GetNRows()-gap-1,0);
1378 printf("Time for parralel tracking outer sectors: \t"); timer.Print();timer.Start();
1379
1380 RemoveOverlap(fSeeds, 0.6,3);
1381 printf("Time for removal overlap- outer sectors: \t");timer.Print();timer.Start();
1382 //parallel tracking
1383 fSectors = fInnerSec;
1384 fN=fkNIS;
1385 ParallelTracking(fSectors->GetNRows()-1,0);
1386 printf("Number of tracks after inner tracking %d\n",fNtracks);
1387 printf("Time for parralel tracking inner sectors: \t"); timer.Print();timer.Start();
1388 //
1389 RemoveOverlap(fSeeds,0.6,5,kTRUE); // remove overlap - shared points signed
1390 printf("Time for removal overlap- inner sectors: \t"); timer.Print();timer.Start();
1391 //
1392 //
1393
1394 ioseed = (AliTPCseed*)(fSeeds->UncheckedAt(0));
1395 AliTPCseed * vseed = new AliTPCseed;
1396 vseed->fPoints = new TClonesArray("AliTPCTrackPoint",1);
1397 vseed->fEPoints = new TClonesArray("AliTPCExactPoint",1);
1398 vseed->fPoints->ExpandCreateFast(2);
1399
1400 TBranch * seedbranch = seedtree.Branch("seeds","AliTPCseed",&vseed,32000,99);
1401 //delete vseed;
1402 nseed=fSeeds->GetEntriesFast();
1403
1404 Int_t found = 0;
1405 for (i=0; i<nseed; i++) {
1406 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
1407 if (!pt) continue;
1408 Int_t nc=t.GetNumberOfClusters();
1409 t.CookdEdx(0.02,0.6);
1410 CookLabel(pt,0.1); //For comparison only
1411 // if ((pt->IsActive()) && (nc>Int_t(0.4*nrows))){
1412 if ((pt->IsActive()) && (nc>Int_t(0.5*t.fNFoundable) && (nc>Int_t(0.3*nrows)))){
1413 iotrack=pt;
1414 tracktree.Fill();
1415 cerr<<found++<<'\r';
1416 }
1417 else
1418 if ( (pt->IsActive())) fNtracks--;
1419 pt->RebuildSeed();
1420 seedbranch->SetAddress(&pt);
1421
1422 seedtree.Fill();
1423 for (Int_t j=0;j<160;j++){
1424 delete pt->fPoints->RemoveAt(j);
1425 }
1426 delete pt->fPoints;
1427 pt->fPoints =0;
1428 delete fSeeds->RemoveAt(i);
1429 }
1430 printf("Time for track writing and dedx cooking: \t"); timer.Print();timer.Start();
1431
1432 UnloadClusters();
1433 printf("Time for unloading cluster: \t"); timer.Print();timer.Start();
1434
1435 tracktree.Write();
1436 seedtree.Write();
1437 cerr<<"Number of found tracks : "<<fNtracks<<"\t"<<found<<endl;
1438
1439 savedir->cd();
1440
1441 return 0;
1442}
1443
1444
1445void AliTPCtrackerMI::ParallelTracking(Int_t rfirst, Int_t rlast)
1446{
1447 //
1448 // try to track in parralel
1449
1450 Int_t nseed=fSeeds->GetEntriesFast();
1451 //prepare seeds for tracking
1452 for (Int_t i=0; i<nseed; i++) {
1453 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
1454 if (!pt) continue;
1455 if (!t.IsActive()) continue;
1456 // follow prolongation to the first layer
1457 FollowProlongation(t, rfirst+1);
1458 }
1459
1460
1461 //
1462 for (Int_t nr=rfirst; nr>=rlast; nr--){
1463 // make indexes with the cluster tracks for given
1464 // for (Int_t i = 0;i<fN;i++)
1465 // fSectors[i][nr].MakeClusterTracks();
1466
1467 // find nearest cluster
1468 for (Int_t i=0; i<nseed; i++) {
1469 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
1470 if (!pt) continue;
1471 if (!pt->IsActive()) continue;
1472 UpdateClusters(t,i,nr);
1473 }
1474 // prolonagate to the nearest cluster - if founded
1475 for (Int_t i=0; i<nseed; i++) {
1476 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
1477 if (!pt) continue;
1478 if (!pt->IsActive()) continue;
1479 FollowToNextCluster(i,nr);
1480 }
1481 // for (Int_t i= 0;i<fN;i++)
1482 // fSectors[i][nr].ClearClusterTracks();
1483 }
1484
1485}
1486
1487Float_t AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed)
1488{
1489 //
1490 //
1491 Float_t sd2 = (fParam->GetZLength()-TMath::Abs(seed->GetZ()))*fParam->GetDiffL()*fParam->GetDiffL();
1492 Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
1493 Float_t sres = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
1494 Float_t angular = seed->GetSnp();
1495 angular = angular*angular/(1-angular*angular);
1496 // angular*=angular;
1497 //angular = TMath::Sqrt(angular/(1-angular));
1498 Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres);
1499 return res;
1500}
1501Float_t AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
1502{
1503 //
1504 //
1505 Float_t sd2 = (fParam->GetZLength()-TMath::Abs(seed->GetZ()))*fParam->GetDiffL()*fParam->GetDiffL();
1506 Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
1507 Float_t sres = fParam->GetZSigma();
1508 Float_t angular = seed->GetTgl();
1509 Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres);
1510 return res;
1511}
1512
1513
1514
1515//_________________________________________________________________________
1516AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1517 //--------------------------------------------------------------------
1518 // Return pointer to a given cluster
1519 //--------------------------------------------------------------------
1520 Int_t sec=(index&0xff000000)>>24;
1521 Int_t row=(index&0x00ff0000)>>16;
1522 Int_t ncl=(index&0x0000ffff)>>00;
1523
1524 AliTPCClustersRow *clrow=((AliTPCtrackerMI *) this)->fClustersArray.GetRow(sec,row);
1525 return (AliTPCclusterMI*)(*clrow)[ncl];
1526}
1527
1528//__________________________________________________________________________
1529void AliTPCtrackerMI::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
1530 //--------------------------------------------------------------------
1531 //This function "cooks" a track label. If label<0, this track is fake.
1532 //--------------------------------------------------------------------
1533 Int_t noc=t->GetNumberOfClusters();
1534 Int_t *lb=new Int_t[noc];
1535 Int_t *mx=new Int_t[noc];
1536 AliTPCclusterMI **clusters=new AliTPCclusterMI*[noc];
1537
1538 Int_t i;
1539 for (i=0; i<noc; i++) {
1540 lb[i]=mx[i]=0;
1541 Int_t index=t->GetClusterIndex(i);
1542 clusters[i]=GetClusterMI(index);
1543 }
1544
1545 Int_t lab=123456789;
1546 for (i=0; i<noc; i++) {
1547 AliTPCclusterMI *c=clusters[i];
1548 lab=TMath::Abs(c->GetLabel(0));
1549 Int_t j;
1550 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
1551 lb[j]=lab;
1552 (mx[j])++;
1553 }
1554
1555 Int_t max=0;
1556 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
1557
1558 for (i=0; i<noc; i++) {
1559 AliTPCclusterMI *c=clusters[i];
1560 if (TMath::Abs(c->GetLabel(1)) == lab ||
1561 TMath::Abs(c->GetLabel(2)) == lab ) max++;
1562 }
1563
1564 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
1565
1566 else {
1567 Int_t tail=Int_t(0.10*noc);
1568 max=0;
1569 for (i=1; i<=tail; i++) {
1570 AliTPCclusterMI *c=clusters[noc-i];
1571 if (lab == TMath::Abs(c->GetLabel(0)) ||
1572 lab == TMath::Abs(c->GetLabel(1)) ||
1573 lab == TMath::Abs(c->GetLabel(2))) max++;
1574 }
1575 if (max < Int_t(0.5*tail)) lab=-lab;
1576 }
1577
1578 t->SetLabel(lab);
1579
1580 delete[] lb;
1581 delete[] mx;
1582 delete[] clusters;
1583}
1584
1585//_________________________________________________________________________
1586void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
1587 //-----------------------------------------------------------------------
1588 // Setup inner sector
1589 //-----------------------------------------------------------------------
1590 if (f==0) {
1591 fAlpha=par->GetInnerAngle();
1592 fAlphaShift=par->GetInnerAngleShift();
1593 fPadPitchWidth=par->GetInnerPadPitchWidth();
1594 fPadPitchLength=par->GetInnerPadPitchLength();
1595 fN=par->GetNRowLow();
1596 fRow=new AliTPCRow[fN];
1597 for (Int_t i=0; i<fN; i++) {
1598 fRow[i].SetX(par->GetPadRowRadiiLow(i));
1599 fRow[i].fDeadZone =1.5; //1.5 cm of dead zone
1600 }
1601 } else {
1602 fAlpha=par->GetOuterAngle();
1603 fAlphaShift=par->GetOuterAngleShift();
1604 fPadPitchWidth = par->GetOuterPadPitchWidth();
1605 fPadPitchLength = par->GetOuter1PadPitchLength();
1606 f1PadPitchLength = par->GetOuter1PadPitchLength();
1607 f2PadPitchLength = par->GetOuter2PadPitchLength();
1608
1609 fN=par->GetNRowUp();
1610 fRow=new AliTPCRow[fN];
1611 for (Int_t i=0; i<fN; i++) {
1612 fRow[i].SetX(par->GetPadRowRadiiUp(i));
1613 fRow[i].fDeadZone =1.5; // 1.5 cm of dead zone
1614 }
1615 }
1616}
1617
1618
1619AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
1620 //
1621 if (fClusterTracks) delete [] fClusterTracks;
1622 fClusterTracks = 0;
1623}
1624
1625void AliTPCtrackerMI::AliTPCRow::MakeClusterTracks(){
1626 //create cluster tracks
1627 if (fN>0)
1628 fClusterTracks = new AliTPCclusterTracks[fN];
1629}
1630
1631void AliTPCtrackerMI::AliTPCRow::ClearClusterTracks(){
1632 if (fClusterTracks) delete[] fClusterTracks;
1633 fClusterTracks =0;
1634}
1635
1636
1637
1638void AliTPCtrackerMI::AliTPCRow::UpdateClusterTrack(Int_t clindex, Int_t trindex, AliTPCseed * seed){
1639 //
1640 //
1641 // update information of the cluster tracks - if track is nearer then other tracks to the
1642 // given track
1643 const AliTPCclusterMI * cl = (*this)[clindex];
1644 AliTPCclusterTracks * cltracks = GetClusterTracks(clindex);
1645 // find the distance of the cluster to the track
1646 Float_t dy2 = (cl->GetY()- seed->GetY());
1647 dy2*=dy2;
1648 Float_t dz2 = (cl->GetZ()- seed->GetZ());
1649 dz2*=dz2;
1650 //
1651 Float_t distance = TMath::Sqrt(dy2+dz2);
1652 if (distance> 3)
1653 return; // MI - to be changed - AliTPCtrackerParam
1654
1655 if ( distance < cltracks->fDistance[0]){
1656 cltracks->fDistance[2] =cltracks->fDistance[1];
1657 cltracks->fDistance[1] =cltracks->fDistance[0];
1658 cltracks->fDistance[0] =distance;
1659 cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1];
1660 cltracks->fTrackIndex[1] =cltracks->fTrackIndex[0];
1661 cltracks->fTrackIndex[0] =trindex;
1662 }
1663 else
1664 if ( distance < cltracks->fDistance[1]){
1665 cltracks->fDistance[2] =cltracks->fDistance[1];
1666 cltracks->fDistance[1] =distance;
1667 cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1];
1668 cltracks->fTrackIndex[1] =trindex;
1669 } else
1670 if (distance < cltracks->fDistance[2]){
1671 cltracks->fDistance[2] =distance;
1672 cltracks->fTrackIndex[2] =trindex;
1673 }
1674}
1675
1676
1677//_________________________________________________________________________
1678void
1679AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
1680 //-----------------------------------------------------------------------
1681 // Insert a cluster into this pad row in accordence with its y-coordinate
1682 //-----------------------------------------------------------------------
1683 if (fN==kMaxClusterPerRow) {
1684 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
1685 }
1686 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
1687 Int_t i=Find(c->GetZ());
1688 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
1689 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
1690 fIndex[i]=index; fClusters[i]=c; fN++;
1691}
1692
1693//___________________________________________________________________
1694Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
1695 //-----------------------------------------------------------------------
1696 // Return the index of the nearest cluster
1697 //-----------------------------------------------------------------------
1698 if (fN==0) return 0;
1699 if (z <= fClusters[0]->GetZ()) return 0;
1700 if (z > fClusters[fN-1]->GetZ()) return fN;
1701 Int_t b=0, e=fN-1, m=(b+e)/2;
1702 for (; b<e; m=(b+e)/2) {
1703 if (z > fClusters[m]->GetZ()) b=m+1;
1704 else e=m;
1705 }
1706 return m;
1707}
1708
1709
1710
1711AliTPCseed::AliTPCseed():AliTPCtrack(){
1712 //
1713 fRow=0;
1714 fRemoval =0;
1715 memset(fClusterIndex,0,sizeof(Int_t)*200);
1716 fPoints = 0;
1717 fEPoints = 0;
1718 fNFoundable =0;
1719 fNShared =0;
1720 fTrackPoints =0;
1721 fRemoval = 0;
1722}
1723
1724AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
1725 fPoints = 0;
1726 fEPoints = 0;
1727 fNShared =0;
1728 fTrackPoints =0;
1729 fRemoval =0;
1730}
1731
1732AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){
1733 fRow=0;
1734 memset(fClusterIndex,0,sizeof(Int_t)*200);
1735 fPoints = 0;
1736 fEPoints = 0;
1737 fNFoundable =0;
1738 fNShared =0;
1739 fTrackPoints =0;
1740 fRemoval =0;
1741}
1742
1743AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15],
1744 Double_t xr, Double_t alpha):
1745 AliTPCtrack(index, xx, cc, xr, alpha) {
1746 //
1747 //
1748 fRow =0;
1749 memset(fClusterIndex,0,sizeof(Int_t)*200);
1750 fPoints = 0;
1751 fEPoints = 0;
1752 fNFoundable =0;
1753 fNShared = 0;
1754 fTrackPoints =0;
1755 fRemoval =0;
1756}
1757
1758AliTPCseed::~AliTPCseed(){
1759 if (fPoints) delete fPoints;
1760 fPoints =0;
1761 fEPoints = 0;
1762 if (fTrackPoints){
1763 for (Int_t i=0;i<8;i++){
1764 delete [] fTrackPoints[i];
1765 }
1766 delete fTrackPoints;
1767 fTrackPoints =0;
1768 }
1769
1770}
1771
1772AliTPCTrackPoint * AliTPCseed::GetTrackPoint(Int_t i)
1773{
1774 //
1775 //
1776 if (!fTrackPoints) {
1777 fTrackPoints = new AliTPCTrackPoint*[8];
1778 for ( Int_t i=0;i<8;i++)
1779 fTrackPoints[i]=0;
1780 }
1781 Int_t index1 = i/20;
1782 if (!fTrackPoints[index1]) fTrackPoints[index1] = new AliTPCTrackPoint[20];
1783 return &(fTrackPoints[index1][i%20]);
1784}
1785
1786void AliTPCseed::RebuildSeed()
1787{
1788 //
1789 // rebuild seed to be ready for storing
1790 fPoints = new TClonesArray("AliTPCTrackPoint",160);
1791 fPoints->ExpandCreateFast(160);
1792 fEPoints = new TClonesArray("AliTPCExactPoint",1);
1793 for (Int_t i=0;i<160;i++){
1794 AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
1795 *trpoint = *(GetTrackPoint(i));
1796 }
1797
1798}
1799
1800//_____________________________________________________________________________
1801void AliTPCseed::CookdEdx(Double_t low, Double_t up) {
1802 //-----------------------------------------------------------------
1803 // This funtion calculates dE/dX within the "low" and "up" cuts.
1804 //-----------------------------------------------------------------
1805
1806 Float_t amp[200];
1807 Float_t angular[200];
1808 Float_t weight[200];
1809 Int_t index[200];
1810 //Int_t nc = 0;
1811 // TClonesArray & arr = *fPoints;
1812 Float_t meanlog = 100.;
1813
1814 Float_t mean[4] = {0,0,0,0};
1815 Float_t sigma[4] = {1000,1000,1000,1000};
1816 Int_t nc[4] = {0,0,0,0};
1817 Float_t norm[4] = {1000,1000,1000,1000};
1818 //
1819 //
1820 fNShared =0;
1821
1822 for (Int_t of =0; of<4; of++){
1823 for (Int_t i=of;i<160;i+=4)
1824 {
1825 //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
1826 AliTPCTrackPoint * point = GetTrackPoint(i);
1827 if (point==0) continue;
1828 if (point->fIsShared){
1829 fNShared++;
1830 continue;
1831 }
1832 if (point->GetCPoint().GetMax()<5) continue;
1833 Float_t angley = point->GetTPoint().GetAngleY();
1834 Float_t anglez = point->GetTPoint().GetAngleZ();
1835 Int_t type = point->GetCPoint().GetType();
1836 Float_t rsigmay = point->GetCPoint().GetSigmaY();
1837 Float_t rsigmaz = point->GetCPoint().GetSigmaZ();
1838 Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
1839
1840 Float_t ampc = 0; // normalization to the number of electrons
1841 if (i>64){
1842 ampc = 1.*point->GetCPoint().GetMax();
1843 //ampc = 1.*point->GetCPoint().GetQ();
1844 // AliTPCClusterPoint & p = point->GetCPoint();
1845 // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
1846 // Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
1847 //Float_t dz =
1848 // TMath::Abs( Int_t(iz) - iz + 0.5);
1849 //ampc *= 1.15*(1-0.3*dy);
1850 //ampc *= 1.15*(1-0.3*dz);
1851
1852 }
1853 else{
1854 ampc = 1.0*point->GetCPoint().GetMax();
1855 //ampc = 1.0*point->GetCPoint().GetQ();
1856 //AliTPCClusterPoint & p = point->GetCPoint();
1857 // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
1858 //Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
1859 //Float_t dz =
1860 // TMath::Abs( Int_t(iz) - iz + 0.5);
1861
1862 //ampc *= 1.15*(1-0.3*dy);
1863 //ampc *= 1.15*(1-0.3*dz);
1864
1865 }
1866 ampc *= 2.0; // put mean value to channel 50
1867 //ampc *= 0.58; // put mean value to channel 50
1868 Float_t w = 1.;
1869 // if (type>0) w = 1./(type/2.-0.5);
1870 if (i<64) {
1871 ampc /= 0.6;
1872 }
1873 if (i>128){
1874 ampc /=1.5;
1875 }
1876 if (type<0) { //amp at the border - lower weight
1877 // w*= 2.;
1878 continue;
1879 }
1880 if (rsigma>1.5) ampc/=1.3; // if big backround
1881 amp[nc[of]] = ampc;
1882 angular[nc[of]] = TMath::Sqrt(1.+angley*angley+anglez*anglez);
1883 weight[nc[of]] = w;
1884 nc[of]++;
1885 }
1886
1887 TMath::Sort(nc[of],amp,index,kFALSE);
1888 Float_t sumamp=0;
1889 Float_t sumamp2=0;
1890 Float_t sumw=0;
1891 //meanlog = amp[index[Int_t(nc[of]*0.33)]];
1892 meanlog = 200;
1893 for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
1894 Float_t ampl = amp[index[i]]/angular[index[i]];
1895 ampl = meanlog*TMath::Log(1.+ampl/meanlog);
1896 //
1897 sumw += weight[index[i]];
1898 sumamp += weight[index[i]]*ampl;
1899 sumamp2 += weight[index[i]]*ampl*ampl;
1900 norm[of] += angular[index[i]]*weight[index[i]];
1901 }
1902 if (sumw<1){
1903 SetdEdx(0);
1904 }
1905 else {
1906 norm[of] /= sumw;
1907 mean[of] = sumamp/sumw;
1908 sigma[of] = sumamp2/sumw-mean[of]*mean[of];
1909 if (sigma[of]>0.1)
1910 sigma[of] = TMath::Sqrt(sigma[of]);
1911 else
1912 sigma[of] = 1000;
1913
1914 mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
1915 //mean *=(1-0.02*(sigma/(mean*0.17)-1.));
1916 //mean *=(1-0.1*(norm-1.));
1917 }
1918 }
1919
1920 Float_t dedx =0;
1921 fSdEdx =0;
1922 fMAngular =0;
1923 // mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
1924 // mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
1925
1926
1927 // dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
1928 // ( TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
1929
1930 Int_t norm2 = 0;
1931 Int_t norm3 = 0;
1932 for (Int_t i =0;i<4;i++){
1933 if (nc[i]>2&&nc[i]<1000){
1934 dedx += mean[i] *nc[i];
1935 fSdEdx += sigma[i]*(nc[i]-2);
1936 fMAngular += norm[i] *nc[i];
1937 norm2 += nc[i];
1938 norm3 += nc[i]-2;
1939 }
1940 fDEDX[i] = mean[i];
1941 fSDEDX[i] = sigma[i];
1942 fNCDEDX[i]= nc[i];
1943 }
1944
1945 if (norm3>0){
1946 dedx /=norm2;
1947 fSdEdx /=norm3;
1948 fMAngular/=norm2;
1949 }
1950 else{
1951 SetdEdx(0);
1952 return;
1953 }
1954 // Float_t dedx1 =dedx;
1955
1956 dedx =0;
1957 for (Int_t i =0;i<4;i++){
1958 if (nc[i]>2&&nc[i]<1000){
1959 mean[i] = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
1960 dedx += mean[i] *nc[i];
1961 }
1962 fDEDX[i] = mean[i];
1963 }
1964 dedx /= norm2;
1965
1966
1967
1968 SetdEdx(dedx);
1969
1970 //mi deDX
1971
1972
1973
1974 //Very rough PID
1975 Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
1976
1977 if (p<0.6) {
1978 if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
1979 if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
1980 SetMass(0.93827); return;
1981 }
1982
1983 if (p<1.2) {
1984 if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
1985 SetMass(0.93827); return;
1986 }
1987
1988 SetMass(0.13957); return;
1989
1990}
1991
1992