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