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