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