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