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