]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCtracker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
18-feb-2003 NvE Memberfunction ShowCalorimeters() introduced in AliEvent to display...
[u/mrichter/AliRoot.git] / TPC / AliTPCtracker.cxx
CommitLineData
73042f01 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/*
17$Log$
e3f7105e 18Revision 1.25 2003/01/28 16:43:35 hristov
19Additional protection: to be discussed with the Root team (M.Ivanov)
20
646a38dc 21Revision 1.24 2002/11/19 16:13:24 hristov
22stdlib.h included to declare exit() on HP
23
1f9a65c4 24Revision 1.23 2002/11/19 11:50:08 hristov
25Removing CONTAINERS (Yu.Belikov)
26
b9d0a01d 27Revision 1.19 2002/07/19 07:31:40 kowal2
28Improvement in tracking by J. Belikov
29
024a7fe9 30Revision 1.18 2002/05/13 07:33:52 kowal2
31Added protection in Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const
32in the case of defined region of interests.
33
11c7b548 34Revision 1.17 2002/03/18 17:59:13 kowal2
35Chnges in the pad geometry - 3 pad lengths introduced.
36
f03e3423 37Revision 1.16 2001/11/08 16:39:03 hristov
38Additional protection (M.Masera)
39
1c8d9aad 40Revision 1.15 2001/11/08 16:36:33 hristov
41Updated V2 stream of tracking (Yu.Belikov). The new long waited features are: 1) Possibility to pass the primary vertex position to the trackers (both for the TPC and the ITS) 2) Possibility to specify the number of tracking passes together with applying (or not applying) the vertex constraint (ITS only) 3) Possibility to make some use of partial PID provided by the TPC when doing tracking in the ITS (ITS only) 4) V0 reconstruction with a helix minimisation of the DCA. (new macros: AliV0FindVertices.C and AliV0Comparison.C) 4a) ( Consequence of the 4) ) All the efficiencies and resolutions are from now on calculated including *secondary*particles* too. (Don't be surprised by the drop in efficiency etc)
42
7f6ddf58 43Revision 1.14 2001/10/21 19:04:55 hristov
44Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
45
f38c8ae5 46Revision 1.13 2001/07/23 12:01:30 hristov
47Initialisation added
48
7afb86d1 49Revision 1.12 2001/07/20 14:32:44 kowal2
50Processing of many events possible now
51
afc42102 52Revision 1.11 2001/05/23 08:50:10 hristov
53Weird inline removed
54
bdbd0f7a 55Revision 1.10 2001/05/16 14:57:25 alibrary
56New files for folders and Stack
57
9e1a0ddb 58Revision 1.9 2001/05/11 07:16:56 hristov
59Fix needed on Sun and Alpha
60
3d39dd8b 61Revision 1.8 2001/05/08 15:00:15 hristov
62Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
63
be9c5115 64Revision 1.5 2000/12/20 07:51:59 kowal2
65Changes suggested by Alessandra and Paolo to avoid overlapped
66data fields in encapsulated classes.
67
d4cf1daa 68Revision 1.4 2000/11/02 07:27:16 kowal2
69code corrections
70
98ab53f4 71Revision 1.2 2000/06/30 12:07:50 kowal2
72Updated from the TPC-PreRelease branch
73
73042f01 74Revision 1.1.2.1 2000/06/25 08:53:55 kowal2
75Splitted from AliTPCtracking
76
77*/
78
79//-------------------------------------------------------
80// Implementation of the TPC tracker
81//
82// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
83//-------------------------------------------------------
73042f01 84#include <TObjArray.h>
85#include <TFile.h>
be5b9287 86#include <TTree.h>
19364939 87#include <Riostream.h>
be5b9287 88
89#include "AliTPCtracker.h"
90#include "AliTPCcluster.h"
b9de75e1 91#include "AliTPCParam.h"
bcc04e2a 92#include "AliClusters.h"
73042f01 93
1f9a65c4 94#include <stdlib.h>
b9de75e1 95//_____________________________________________________________________________
bcc04e2a 96AliTPCtracker::AliTPCtracker(const AliTPCParam *par):
7f6ddf58 97AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
b9de75e1 98{
99 //---------------------------------------------------------------------
100 // The main TPC tracker constructor
101 //---------------------------------------------------------------------
102 fInnerSec=new AliTPCSector[fkNIS];
103 fOuterSec=new AliTPCSector[fkNOS];
104
105 Int_t i;
106 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
107 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
108
e3f7105e 109 fParam = (AliTPCParam*) par;
b9de75e1 110 fSeeds=0;
111}
112
113//_____________________________________________________________________________
114AliTPCtracker::~AliTPCtracker() {
115 //------------------------------------------------------------------
116 // TPC tracker destructor
117 //------------------------------------------------------------------
118 delete[] fInnerSec;
119 delete[] fOuterSec;
1c8d9aad 120 if (fSeeds) {
121 fSeeds->Delete();
122 delete fSeeds;
123 }
b9de75e1 124}
73042f01 125
126//_____________________________________________________________________________
127Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
128{
129 //
130 // Parametrised error of the cluster reconstruction (pad direction)
131 //
132 // Sigma rphi
133 const Float_t kArphi=0.41818e-2;
134 const Float_t kBrphi=0.17460e-4;
135 const Float_t kCrphi=0.30993e-2;
136 const Float_t kDrphi=0.41061e-3;
137
138 pt=TMath::Abs(pt)*1000.;
139 Double_t x=r/pt;
140 tgl=TMath::Abs(tgl);
141 Double_t s=kArphi - kBrphi*r*tgl + kCrphi*x*x + kDrphi*x;
142 if (s<0.4e-3) s=0.4e-3;
143 s*=1.3; //Iouri Belikov
144
145 return s;
146}
147
148//_____________________________________________________________________________
149Double_t SigmaZ2(Double_t r, Double_t tgl)
150{
151 //
152 // Parametrised error of the cluster reconstruction (drift direction)
153 //
154 // Sigma z
155 const Float_t kAz=0.39614e-2;
156 const Float_t kBz=0.22443e-4;
157 const Float_t kCz=0.51504e-1;
158
159
160 tgl=TMath::Abs(tgl);
161 Double_t s=kAz - kBz*r*tgl + kCz*tgl*tgl;
162 if (s<0.4e-3) s=0.4e-3;
163 s*=1.3; //Iouri Belikov
164
165 return s;
166}
167
168//_____________________________________________________________________________
bdbd0f7a 169Double_t f1(Double_t x1,Double_t y1,
73042f01 170 Double_t x2,Double_t y2,
171 Double_t x3,Double_t y3)
172{
173 //-----------------------------------------------------------------
174 // Initial approximation of the track curvature
175 //-----------------------------------------------------------------
176 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
177 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
178 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
179 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
180 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
181
182 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
183
184 return -xr*yr/sqrt(xr*xr+yr*yr);
185}
186
187
188//_____________________________________________________________________________
bdbd0f7a 189Double_t f2(Double_t x1,Double_t y1,
73042f01 190 Double_t x2,Double_t y2,
191 Double_t x3,Double_t y3)
192{
193 //-----------------------------------------------------------------
194 // Initial approximation of the track curvature times center of curvature
195 //-----------------------------------------------------------------
196 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
197 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
198 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
199 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
200 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
201
202 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
203
204 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
205}
206
207//_____________________________________________________________________________
bdbd0f7a 208Double_t f3(Double_t x1,Double_t y1,
73042f01 209 Double_t x2,Double_t y2,
210 Double_t z1,Double_t z2)
211{
212 //-----------------------------------------------------------------
213 // Initial approximation of the tangent of the track dip angle
214 //-----------------------------------------------------------------
215 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
216}
217
218//_____________________________________________________________________________
bcc04e2a 219void AliTPCtracker::LoadClusters() {
b9de75e1 220 //-----------------------------------------------------------------
bcc04e2a 221 // This function loads TPC clusters.
b9de75e1 222 //-----------------------------------------------------------------
bcc04e2a 223 if (!gFile->IsOpen()) {
224 cerr<<"AliTPCtracker::LoadClusters : "<<
225 "file with clusters has not been open !\n";
226 return;
b9de75e1 227 }
228
bcc04e2a 229 Char_t name[100];
230 sprintf(name,"TreeC_TPC_%d",GetEventNumber());
231 TTree *cTree=(TTree*)gFile->Get(name);
232 if (!cTree) {
233 cerr<<"AliTPCtracker::LoadClusters : "<<
234 "can't get the tree with TPC clusters !\n";
235 return;
b9de75e1 236 }
237
bcc04e2a 238 TBranch *branch=cTree->GetBranch("Segment");
239 if (!branch) {
240 cerr<<"AliTPCtracker::LoadClusters : "<<
241 "can't get the segment branch !\n";
242 return;
b9de75e1 243 }
646a38dc 244 AliClusters carray, *addr=&carray;
245 carray.SetClass("AliTPCcluster");
246 carray.SetArray(0);
bcc04e2a 247 branch->SetAddress(&addr);
248
249 Int_t nentr=(Int_t)cTree->GetEntries();
250
251 for (Int_t i=0; i<nentr; i++) {
252 cTree->GetEvent(i);
253
254 Int_t ncl=carray.GetArray()->GetEntriesFast();
255
256 Int_t nir=fInnerSec->GetNRows(), nor=fOuterSec->GetNRows();
257 Int_t id=carray.GetID();
258 if ((id<0) || (id>2*(fkNIS*nir + fkNOS*nor))) {
259 cerr<<"AliTPCtracker::LoadClusters : "<<
260 "wrong index !\n";
261 exit(1);
262 }
263 Int_t outindex = 2*fkNIS*nir;
264 if (id<outindex) {
265 Int_t sec = id/nir;
266 Int_t row = id - sec*nir;
267 sec %= fkNIS;
268 AliTPCRow &padrow=fInnerSec[sec][row];
269 while (ncl--) {
270 AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
271 padrow.InsertCluster(c,sec,row);
272 }
273 } else {
274 id -= outindex;
275 Int_t sec = id/nor;
276 Int_t row = id - sec*nor;
277 sec %= fkNOS;
278 AliTPCRow &padrow=fOuterSec[sec][row];
279 while (ncl--) {
280 AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
281 padrow.InsertCluster(c,sec+fkNIS,row);
282 }
283 }
b9de75e1 284
bcc04e2a 285 carray.GetArray()->Clear();
286 }
287 delete cTree;
b9de75e1 288}
289
290//_____________________________________________________________________________
bcc04e2a 291void AliTPCtracker::UnloadClusters() {
b9de75e1 292 //-----------------------------------------------------------------
bcc04e2a 293 // This function unloads TPC clusters.
b9de75e1 294 //-----------------------------------------------------------------
bcc04e2a 295 Int_t i;
296 for (i=0; i<fkNIS; i++) {
297 Int_t nr=fInnerSec->GetNRows();
298 for (Int_t n=0; n<nr; n++) fInnerSec[i][n].ResetClusters();
299 }
300 for (i=0; i<fkNOS; i++) {
301 Int_t nr=fOuterSec->GetNRows();
302 for (Int_t n=0; n<nr; n++) fOuterSec[i][n].ResetClusters();
b9de75e1 303 }
b9de75e1 304}
305
306//_____________________________________________________________________________
307Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
73042f01 308 //-----------------------------------------------------------------
309 // This function tries to find a track prolongation.
310 //-----------------------------------------------------------------
b9de75e1 311 Double_t xt=t.GetX();
312 const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip :
313 Int_t(0.5*fSectors->GetNRows());
73042f01 314 Int_t tryAgain=kSKIP;
73042f01 315
b9de75e1 316 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
317 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
318 if (alpha < 0. ) alpha += 2.*TMath::Pi();
319 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
320
f03e3423 321 Int_t nrows=fSectors->GetRowNumber(xt)-1;
322 for (Int_t nr=nrows; nr>=rf; nr--) {
b9de75e1 323 Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
73042f01 324 if (!t.PropagateTo(x)) return 0;
325
326 AliTPCcluster *cl=0;
327 UInt_t index=0;
b9de75e1 328 Double_t maxchi2=kMaxCHI2;
329 const AliTPCRow &krow=fSectors[s][nr];
330 Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
331 Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),pt);
73042f01 332 Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
333 Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
334
b9de75e1 335 if (road>kMaxROAD) {
73042f01 336 if (t.GetNumberOfClusters()>4)
337 cerr<<t.GetNumberOfClusters()
338 <<"FindProlongation warning: Too broad road !\n";
339 return 0;
340 }
341
342 if (krow) {
343 for (Int_t i=krow.Find(y-road); i<krow; i++) {
be5b9287 344 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
73042f01 345 if (c->GetY() > y+road) break;
346 if (c->IsUsed()) continue;
347 if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue;
348 Double_t chi2=t.GetPredictedChi2(c);
349 if (chi2 > maxchi2) continue;
350 maxchi2=chi2;
351 cl=c;
352 index=krow.GetIndex(i);
353 }
354 }
355 if (cl) {
b9de75e1 356 Float_t l=fSectors->GetPadPitchWidth();
024a7fe9 357 Float_t corr=1.; if (nr>63) corr=0.67; // new (third) pad response !
358 t.SetSampledEdx(cl->GetQ()/l*corr,t.GetNumberOfClusters());
be9c5115 359 if (!t.Update(cl,maxchi2,index)) {
360 if (!tryAgain--) return 0;
361 } else tryAgain=kSKIP;
73042f01 362 } else {
363 if (tryAgain==0) break;
364 if (y > ymax) {
b9de75e1 365 s = (s+1) % fN;
366 if (!t.Rotate(fSectors->GetAlpha())) return 0;
73042f01 367 } else if (y <-ymax) {
b9de75e1 368 s = (s-1+fN) % fN;
369 if (!t.Rotate(-fSectors->GetAlpha())) return 0;
73042f01 370 }
371 tryAgain--;
372 }
373 }
374
375 return 1;
b9de75e1 376}
377
378//_____________________________________________________________________________
379Int_t AliTPCtracker::FollowBackProlongation
380(AliTPCseed& seed, const AliTPCtrack &track) {
381 //-----------------------------------------------------------------
382 // This function propagates tracks back through the TPC
383 //-----------------------------------------------------------------
384 Double_t alpha=seed.GetAlpha() - fSectors->GetAlphaShift();
385 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
386 if (alpha < 0. ) alpha += 2.*TMath::Pi();
387 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
388
389 Int_t idx=-1, sec=-1, row=-1;
390 Int_t nc=seed.GetLabel(); //index of the cluster to start with
391 if (nc--) {
392 idx=track.GetClusterIndex(nc);
393 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
394 }
bcc04e2a 395 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
396 else { if (sec < fkNIS) row=-1; }
b9de75e1 397
398 Int_t nr=fSectors->GetNRows();
399 for (Int_t i=0; i<nr; i++) {
400 Double_t x=fSectors->GetX(i), ymax=fSectors->GetMaxY(i);
401
402 if (!seed.PropagateTo(x)) return 0;
403
404 Double_t y=seed.GetY();
405 if (y > ymax) {
406 s = (s+1) % fN;
407 if (!seed.Rotate(fSectors->GetAlpha())) return 0;
408 } else if (y <-ymax) {
409 s = (s-1+fN) % fN;
410 if (!seed.Rotate(-fSectors->GetAlpha())) return 0;
411 }
412
413 AliTPCcluster *cl=0;
414 Int_t index=0;
415 Double_t maxchi2=kMaxCHI2;
416 Double_t pt=seed.GetConvConst()/(100/0.299792458/0.2)/seed.Get1Pt();
417 Double_t sy2=SigmaY2(seed.GetX(),seed.GetTgl(),pt);
418 Double_t sz2=SigmaZ2(seed.GetX(),seed.GetTgl());
419 Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ();
420 if (road>kMaxROAD) {
421 cerr<<seed.GetNumberOfClusters()
422 <<"AliTPCtracker::FollowBackProlongation: Too broad road !\n";
423 return 0;
424 }
425
426
427 Int_t accepted=seed.GetNumberOfClusters();
428 if (row==i) {
429 //try to accept already found cluster
430 AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
431 Double_t chi2;
432 if ((chi2=seed.GetPredictedChi2(c))<maxchi2 || accepted<27) {
433 index=idx; cl=c; maxchi2=chi2;
434 } //else cerr<<"AliTPCtracker::FollowBackProlongation: oulier !\n";
435
436 if (nc--) {
437 idx=track.GetClusterIndex(nc);
438 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
439 }
bcc04e2a 440 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
441 else { if (sec < fkNIS) row=-1; }
73042f01 442
b9de75e1 443 }
444 if (!cl) {
445 //try to fill the gap
446 const AliTPCRow &krow=fSectors[s][i];
447 if (accepted>27)
448 if (krow) {
449 for (Int_t i=krow.Find(y-road); i<krow; i++) {
450 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
451 if (c->GetY() > y+road) break;
452 if (c->IsUsed()) continue;
453 if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue;
454 Double_t chi2=seed.GetPredictedChi2(c);
455 if (chi2 > maxchi2) continue;
456 maxchi2=chi2;
457 cl=c;
458 index=krow.GetIndex(i);
459 }
460 }
461 }
462
463 if (cl) {
464 Float_t l=fSectors->GetPadPitchWidth();
024a7fe9 465 Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
466 seed.SetSampledEdx(cl->GetQ()/l*corr,seed.GetNumberOfClusters());
b9de75e1 467 seed.Update(cl,maxchi2,index);
468 }
469
470 }
471
472 seed.SetLabel(nc);
473
474 return 1;
73042f01 475}
476
477//_____________________________________________________________________________
b9de75e1 478void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
73042f01 479 //-----------------------------------------------------------------
480 // This function creates track seeds.
481 //-----------------------------------------------------------------
482 Double_t x[5], c[15];
483
bcc04e2a 484 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
73042f01 485 Double_t cs=cos(alpha), sn=sin(alpha);
486
bcc04e2a 487 Double_t x1 =fSectors->GetX(i1);
488 Double_t xx2=fSectors->GetX(i2);
73042f01 489
bcc04e2a 490 for (Int_t ns=0; ns<fN; ns++) {
491 Int_t nl=fSectors[(ns-1+fN)%fN][i2];
492 Int_t nm=fSectors[ns][i2];
493 Int_t nu=fSectors[(ns+1)%fN][i2];
494 const AliTPCRow& kr1=fSectors[ns][i1];
73042f01 495 for (Int_t is=0; is < kr1; is++) {
496 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
497 for (Int_t js=0; js < nl+nm+nu; js++) {
498 const AliTPCcluster *kcl;
499 Double_t x2, y2, z2;
7f6ddf58 500 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
73042f01 501
502 if (js<nl) {
bcc04e2a 503 const AliTPCRow& kr2=fSectors[(ns-1+fN)%fN][i2];
73042f01 504 kcl=kr2[js];
505 y2=kcl->GetY(); z2=kcl->GetZ();
506 x2= xx2*cs+y2*sn;
507 y2=-xx2*sn+y2*cs;
508 } else
509 if (js<nl+nm) {
bcc04e2a 510 const AliTPCRow& kr2=fSectors[ns][i2];
73042f01 511 kcl=kr2[js-nl];
512 x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
513 } else {
bcc04e2a 514 const AliTPCRow& kr2=fSectors[(ns+1)%fN][i2];
73042f01 515 kcl=kr2[js-nl-nm];
516 y2=kcl->GetY(); z2=kcl->GetZ();
517 x2=xx2*cs-y2*sn;
518 y2=xx2*sn+y2*cs;
519 }
520
7f6ddf58 521 Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2);
73042f01 522 if (TMath::Abs(zz-z2)>5.) continue;
523
524 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
525 if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
526
527 x[0]=y1;
528 x[1]=z1;
b9de75e1 529 x[4]=f1(x1,y1,x2,y2,x3,y3);
530 if (TMath::Abs(x[4]) >= 0.0066) continue;
be5b9287 531 x[2]=f2(x1,y1,x2,y2,x3,y3);
b9de75e1 532 //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
533 x[3]=f3(x1,y1,x2,y2,z1,z2);
534 if (TMath::Abs(x[3]) > 1.2) continue;
be5b9287 535 Double_t a=asin(x[2]);
b9de75e1 536 Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
7f6ddf58 537 if (TMath::Abs(zv-z3)>10.) continue;
73042f01 538
539 Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
540 Double_t sy2=kcl->GetSigmaY2(), sz2=kcl->GetSigmaZ2();
024a7fe9 541 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
542 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
73042f01 543
b9de75e1 544 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
545 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
546 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
be5b9287 547 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
548 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
b9de75e1 549 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
550 Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
551 Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
552 Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
553 Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
73042f01 554
555 c[0]=sy1;
556 c[1]=0.; c[2]=sz1;
b9de75e1 557 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
558 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
559 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
560 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
561 c[13]=f30*sy1*f40+f32*sy2*f42;
562 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
73042f01 563
564 UInt_t index=kr1.GetIndex(is);
565 AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
bcc04e2a 566 Float_t l=fSectors->GetPadPitchWidth();
73042f01 567 track->SetSampledEdx(kr1[is]->GetQ()/l,0);
568
b9de75e1 569 Int_t rc=FollowProlongation(*track, i2);
be9c5115 570 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track;
b9de75e1 571 else fSeeds->AddLast(track);
73042f01 572 }
573 }
574 }
575}
576
577//_____________________________________________________________________________
b9de75e1 578Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
73042f01 579 //-----------------------------------------------------------------
b9de75e1 580 // This function reades track seeds.
73042f01 581 //-----------------------------------------------------------------
582 TDirectory *savedir=gDirectory;
583
b9de75e1 584 TFile *in=(TFile*)inp;
585 if (!in->IsOpen()) {
586 cerr<<"AliTPCtracker::ReadSeeds(): input file is not open !\n";
be5b9287 587 return 1;
73042f01 588 }
589
b9de75e1 590 in->cd();
591 TTree *seedTree=(TTree*)in->Get("Seeds");
592 if (!seedTree) {
593 cerr<<"AliTPCtracker::ReadSeeds(): ";
594 cerr<<"can't get a tree with track seeds !\n";
595 return 2;
596 }
597 AliTPCtrack *seed=new AliTPCtrack;
598 seedTree->SetBranchAddress("tracks",&seed);
599
600 if (fSeeds==0) fSeeds=new TObjArray(15000);
73042f01 601
b9de75e1 602 Int_t n=(Int_t)seedTree->GetEntries();
603 for (Int_t i=0; i<n; i++) {
604 seedTree->GetEvent(i);
605 fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
606 }
607
608 delete seed;
f38c8ae5 609
610 delete seedTree; //Thanks to Mariana Bondila
611
b9de75e1 612 savedir->cd();
73042f01 613
b9de75e1 614 return 0;
615}
73042f01 616
b9de75e1 617//_____________________________________________________________________________
618Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
619 //-----------------------------------------------------------------
620 // This is a track finder.
621 //-----------------------------------------------------------------
622 TDirectory *savedir=gDirectory;
73042f01 623
b9de75e1 624 if (inp) {
625 TFile *in=(TFile*)inp;
626 if (!in->IsOpen()) {
627 cerr<<"AliTPCtracker::Clusters2Tracks(): input file is not open !\n";
628 return 1;
629 }
630 }
73042f01 631
b9de75e1 632 if (!out->IsOpen()) {
633 cerr<<"AliTPCtracker::Clusters2Tracks(): output file is not open !\n";
634 return 2;
73042f01 635 }
636
bcc04e2a 637 LoadClusters();
638
b9de75e1 639 out->cd();
7f6ddf58 640
641 char tname[100];
bcc04e2a 642 sprintf(tname,"TreeT_TPC_%d",GetEventNumber());
7f6ddf58 643 TTree tracktree(tname,"Tree with TPC tracks");
b9de75e1 644 AliTPCtrack *iotrack=0;
645 tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
646
73042f01 647 //find track seeds
b9de75e1 648 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
73042f01 649 Int_t nrows=nlow+nup;
b9de75e1 650 if (fSeeds==0) {
651 Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
bcc04e2a 652 fSectors=fOuterSec; fN=fkNOS;
653 fSeeds=new TObjArray(15000);
b9de75e1 654 MakeSeeds(nup-1, nup-1-gap);
655 MakeSeeds(nup-1-shift, nup-1-shift-gap);
bcc04e2a 656 }
b9de75e1 657 fSeeds->Sort();
73042f01 658
bcc04e2a 659 Int_t found=0;
b9de75e1 660 Int_t nseed=fSeeds->GetEntriesFast();
bcc04e2a 661 for (Int_t i=0; i<nseed; i++) {
662 //tracking in the outer sectors
663 fSectors=fOuterSec; fN=fkNOS;
664
b9de75e1 665 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
bcc04e2a 666 if (!FollowProlongation(t)) {
667 delete fSeeds->RemoveAt(i);
73042f01 668 continue;
669 }
73042f01 670
bcc04e2a 671 //tracking in the inner sectors
672 fSectors=fInnerSec; fN=fkNIS;
73042f01 673
b9de75e1 674 Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
73042f01 675 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
676 if (alpha < 0. ) alpha += 2.*TMath::Pi();
b9de75e1 677 Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
73042f01 678
bcc04e2a 679 alpha=ns*fInnerSec->GetAlpha()+fInnerSec->GetAlphaShift()-t.GetAlpha();
73042f01 680
681 if (t.Rotate(alpha)) {
e3f7105e 682 if (FollowProlongation(t)) {
683 if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
684 t.CookdEdx();
685 CookLabel(pt,0.1); //For comparison only
686 iotrack=pt;
687 tracktree.Fill();
688 UseClusters(&t);
689 cerr<<found++<<'\r';
690 }
691 }
73042f01 692 }
b9de75e1 693 delete fSeeds->RemoveAt(i);
bcc04e2a 694 }
e3f7105e 695
7f6ddf58 696 tracktree.Write();
73042f01 697
b9de75e1 698 cerr<<"Number of found tracks : "<<found<<endl;
699
700 savedir->cd();
701
bcc04e2a 702 UnloadClusters();
703 fSeeds->Clear(); delete fSeeds; fSeeds=0;
704
b9de75e1 705 return 0;
706}
707
708//_____________________________________________________________________________
709Int_t AliTPCtracker::PropagateBack(const TFile *inp, TFile *out) {
710 //-----------------------------------------------------------------
711 // This function propagates tracks back through the TPC.
712 //-----------------------------------------------------------------
713 fSeeds=new TObjArray(15000);
714
715 TFile *in=(TFile*)inp;
716 TDirectory *savedir=gDirectory;
717
718 if (!in->IsOpen()) {
719 cerr<<"AliTPCtracker::PropagateBack(): ";
720 cerr<<"file with back propagated ITS tracks is not open !\n";
721 return 1;
73042f01 722 }
723
b9de75e1 724 if (!out->IsOpen()) {
725 cerr<<"AliTPCtracker::PropagateBack(): ";
726 cerr<<"file for back propagated TPC tracks is not open !\n";
727 return 2;
728 }
729
bcc04e2a 730 LoadClusters();
731
b9de75e1 732 in->cd();
f38c8ae5 733 TTree *bckTree=(TTree*)in->Get("TreeT_ITSb_0");
b9de75e1 734 if (!bckTree) {
735 cerr<<"AliTPCtracker::PropagateBack() ";
736 cerr<<"can't get a tree with back propagated ITS tracks !\n";
737 return 3;
738 }
739 AliTPCtrack *bckTrack=new AliTPCtrack;
740 bckTree->SetBranchAddress("tracks",&bckTrack);
741
f38c8ae5 742 TTree *tpcTree=(TTree*)in->Get("TreeT_TPC_0");
b9de75e1 743 if (!tpcTree) {
744 cerr<<"AliTPCtracker::PropagateBack() ";
745 cerr<<"can't get a tree with TPC tracks !\n";
746 return 4;
747 }
748 AliTPCtrack *tpcTrack=new AliTPCtrack;
749 tpcTree->SetBranchAddress("tracks",&tpcTrack);
750
e3f7105e 751 //*** Prepare an array of tracks to be back propagated
b9de75e1 752 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
753 Int_t nrows=nlow+nup;
754
e3f7105e 755 //
756 // Match ITS tracks with old TPC tracks
757 // the tracks do not have to be sorted [SR, GSI, 18.02.2003]
758 //
759 // the algorithm is linear and uses LUT for sorting
760 // cost of the algorithm: 2 * number of tracks
761 //
762
763 TObjArray tracks(15000);
764 Int_t i=0;
765 Int_t tpcN= (Int_t)tpcTree->GetEntries();
766 Int_t bckN= (bckTree)? (Int_t)bckTree->GetEntries() : 0;
767
768 // look up table
769 const Int_t nLab = 20000; // limit on number of primaries (arbitrary)
770 Int_t lut[nLab];
771 for(Int_t i=0; i<nLab; i++) lut[i] = -1;
772
773 if (bckTree) {
774 for(Int_t i=0; i<bckN; i++) {
775 bckTree->GetEvent(i);
776 Int_t lab = TMath::Abs(bckTrack->GetLabel());
777 if (lab < nLab) lut[lab] = i;
778 }
779 }
780
781 for (Int_t i=0; i<tpcN; i++) {
782 tpcTree->GetEvent(i);
783 Double_t alpha=tpcTrack->GetAlpha();
784
785 if (!bckTree) {
786
787 // No ITS - use TPC track only
788
789 fSeeds->AddLast(new AliTPCseed(*tpcTrack, tpcTrack->GetAlpha()));
790 tracks.AddLast(new AliTPCtrack(*tpcTrack));
791
792 } else {
793
794 // with ITS
795 // discard not prolongated tracks (to be discussed)
796
797 Int_t lab = TMath::Abs(tpcTrack->GetLabel());
798 if (lab > nLab || lut[lab] < 0) continue;
799
800 bckTree->GetEvent(lut[lab]);
801 bckTrack->Rotate(alpha-bckTrack->GetAlpha());
802
803 fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
804 tracks.AddLast(new AliTPCtrack(*tpcTrack));
805 }
806 }
807
808 // end of matching
809
810 /*
b9de75e1 811 TObjArray tracks(15000);
e3f7105e 812 Int_t i=0, j=0;
b9de75e1 813 Int_t tpcN=(Int_t)tpcTree->GetEntries();
814 Int_t bckN=(Int_t)bckTree->GetEntries();
815 if (j<bckN) bckTree->GetEvent(j++);
816 for (i=0; i<tpcN; i++) {
817 tpcTree->GetEvent(i);
818 Double_t alpha=tpcTrack->GetAlpha();
819 if (j<bckN &&
820 TMath::Abs(tpcTrack->GetLabel())==TMath::Abs(bckTrack->GetLabel())) {
821 if (!bckTrack->Rotate(alpha-bckTrack->GetAlpha())) continue;
822 fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
823 bckTree->GetEvent(j++);
824 } else {
825 tpcTrack->ResetCovariance();
826 fSeeds->AddLast(new AliTPCseed(*tpcTrack,alpha));
827 }
828 tracks.AddLast(new AliTPCtrack(*tpcTrack));
829 }
e3f7105e 830 */
831
b9de75e1 832 out->cd();
e3f7105e 833
834 // tree name seedsTPCtoTRD as expected by TRD and as
835 // discussed and decided in Strasbourg (May 2002)
836 // [SR, GSI, 18.02.2003]
837
838 //TTree backTree("TreeT_TPCb_0","Tree with back propagated TPC tracks");
839 TTree backTree("seedsTPCtoTRD_0","Tree with back propagated TPC tracks");
b9de75e1 840 AliTPCtrack *otrack=0;
841 backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
842
843//*** Back propagation through inner sectors
bcc04e2a 844 fSectors=fInnerSec; fN=fkNIS;
b9de75e1 845 Int_t nseed=fSeeds->GetEntriesFast();
846 for (i=0; i<nseed; i++) {
847 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
848 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
73042f01 849
b9de75e1 850 Int_t nc=t.GetNumberOfClusters();
851 s.SetLabel(nc-1); //set number of the cluster to start with
73042f01 852
b9de75e1 853 if (FollowBackProlongation(s,t)) {
854 UseClusters(&s);
855 continue;
856 }
857 delete fSeeds->RemoveAt(i);
858 }
b9de75e1 859
860//*** Back propagation through outer sectors
bcc04e2a 861 fSectors=fOuterSec; fN=fkNOS;
b9de75e1 862 Int_t found=0;
863 for (i=0; i<nseed; i++) {
864 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
865 if (!ps) continue;
866 Int_t nc=s.GetNumberOfClusters();
867 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
868
869 Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
870 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
871 if (alpha < 0. ) alpha += 2.*TMath::Pi();
872 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
873
874 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
875 alpha-=s.GetAlpha();
876
877 if (s.Rotate(alpha)) {
e3f7105e 878 if (FollowBackProlongation(s,t)) {
879 if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
880 s.CookdEdx();
881 s.SetLabel(t.GetLabel());
882 UseClusters(&s,nc);
883
884 // Propagate to outer reference plane for comparison reasons
885 // reason for keeping fParam object [SR, GSI, 18.02.2003]
886
887 ps->PropagateTo(fParam->GetOuterRadiusUp());
888 otrack=ps;
889 backTree.Fill();
890 cerr<<found++<<'\r';
891 continue;
892 }
893 }
b9de75e1 894 }
895 delete fSeeds->RemoveAt(i);
896 }
b9de75e1 897
898 backTree.Write();
73042f01 899 savedir->cd();
b9de75e1 900 cerr<<"Number of seeds: "<<nseed<<endl;
901 cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
902 cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
903
904 delete bckTrack;
905 delete tpcTrack;
be5b9287 906
f38c8ae5 907 delete bckTree; //Thanks to Mariana Bondila
908 delete tpcTree; //Thanks to Mariana Bondila
909
bcc04e2a 910 UnloadClusters();
911
be5b9287 912 return 0;
73042f01 913}
914
b9de75e1 915//_________________________________________________________________________
916AliCluster *AliTPCtracker::GetCluster(Int_t index) const {
917 //--------------------------------------------------------------------
918 // Return pointer to a given cluster
919 //--------------------------------------------------------------------
920 Int_t sec=(index&0xff000000)>>24;
921 Int_t row=(index&0x00ff0000)>>16;
922 Int_t ncl=(index&0x0000ffff)>>00;
923
bcc04e2a 924 const AliTPCcluster *cl=0;
925
926 if (sec<fkNIS) {
927 cl=fInnerSec[sec][row].GetUnsortedCluster(ncl);
928 } else {
929 sec-=fkNIS;
930 cl=fOuterSec[sec][row].GetUnsortedCluster(ncl);
931 }
932
933 return (AliCluster*)cl;
b9de75e1 934}
935
936//__________________________________________________________________________
937void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
938 //--------------------------------------------------------------------
939 //This function "cooks" a track label. If label<0, this track is fake.
940 //--------------------------------------------------------------------
941 Int_t noc=t->GetNumberOfClusters();
942 Int_t *lb=new Int_t[noc];
943 Int_t *mx=new Int_t[noc];
944 AliCluster **clusters=new AliCluster*[noc];
945
946 Int_t i;
947 for (i=0; i<noc; i++) {
948 lb[i]=mx[i]=0;
949 Int_t index=t->GetClusterIndex(i);
950 clusters[i]=GetCluster(index);
951 }
952
953 Int_t lab=123456789;
954 for (i=0; i<noc; i++) {
955 AliCluster *c=clusters[i];
956 lab=TMath::Abs(c->GetLabel(0));
957 Int_t j;
958 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
959 lb[j]=lab;
960 (mx[j])++;
961 }
962
963 Int_t max=0;
964 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
965
966 for (i=0; i<noc; i++) {
967 AliCluster *c=clusters[i];
968 if (TMath::Abs(c->GetLabel(1)) == lab ||
969 TMath::Abs(c->GetLabel(2)) == lab ) max++;
970 }
971
972 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
973
974 else {
975 Int_t tail=Int_t(0.10*noc);
976 max=0;
977 for (i=1; i<=tail; i++) {
978 AliCluster *c=clusters[noc-i];
979 if (lab == TMath::Abs(c->GetLabel(0)) ||
980 lab == TMath::Abs(c->GetLabel(1)) ||
981 lab == TMath::Abs(c->GetLabel(2))) max++;
982 }
983 if (max < Int_t(0.5*tail)) lab=-lab;
984 }
985
986 t->SetLabel(lab);
987
988 delete[] lb;
989 delete[] mx;
990 delete[] clusters;
991}
992
993//_________________________________________________________________________
994void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
995 //-----------------------------------------------------------------------
996 // Setup inner sector
997 //-----------------------------------------------------------------------
998 if (f==0) {
999 fAlpha=par->GetInnerAngle();
1000 fAlphaShift=par->GetInnerAngleShift();
1001 fPadPitchWidth=par->GetInnerPadPitchWidth();
024a7fe9 1002 f1PadPitchLength=par->GetInnerPadPitchLength();
1003 f2PadPitchLength=f1PadPitchLength;
b9de75e1 1004
1005 fN=par->GetNRowLow();
1006 fRow=new AliTPCRow[fN];
1007 for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
1008 } else {
1009 fAlpha=par->GetOuterAngle();
1010 fAlphaShift=par->GetOuterAngleShift();
1011 fPadPitchWidth=par->GetOuterPadPitchWidth();
f03e3423 1012 f1PadPitchLength=par->GetOuter1PadPitchLength();
1013 f2PadPitchLength=par->GetOuter2PadPitchLength();
b9de75e1 1014
1015 fN=par->GetNRowUp();
1016 fRow=new AliTPCRow[fN];
f03e3423 1017 for (Int_t i=0; i<fN; i++){
1018 fRow[i].SetX(par->GetPadRowRadiiUp(i));
1019 }
b9de75e1 1020 }
1021}
1022
73042f01 1023//_________________________________________________________________________
bcc04e2a 1024void AliTPCtracker::
1025AliTPCRow::InsertCluster(const AliTPCcluster* c, Int_t sec, Int_t row) {
73042f01 1026 //-----------------------------------------------------------------------
1027 // Insert a cluster into this pad row in accordence with its y-coordinate
1028 //-----------------------------------------------------------------------
b9de75e1 1029 if (fN==kMaxClusterPerRow) {
73042f01 1030 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
1031 }
bcc04e2a 1032
1033 Int_t index=(((sec<<8)+row)<<16)+fN;
1034
1035 if (fN==0) {
1036 fIndex[0]=index;
1037 fClusterArray[0]=*c; fClusters[fN++]=fClusterArray;
1038 return;
1039 }
1040
1041 if (fN==fSize) {
1042 Int_t size=fSize*2;
1043 AliTPCcluster *buff=new AliTPCcluster[size];
1044 memcpy(buff,fClusterArray,fSize*sizeof(AliTPCcluster));
1045 for (Int_t i=0; i<fN; i++)
1046 fClusters[i]=buff+(fClusters[i]-fClusterArray);
1047 delete[] fClusterArray;
1048 fClusterArray=buff;
1049 fSize=size;
1050 }
1051
73042f01 1052 Int_t i=Find(c->GetY());
1053 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
1054 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
bcc04e2a 1055 fIndex[i]=index;
1056 fClusters[i]=fClusterArray+fN; fClusterArray[fN++]=*c;
73042f01 1057}
1058
1059//___________________________________________________________________
1060Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
1061 //-----------------------------------------------------------------------
1062 // Return the index of the nearest cluster
1063 //-----------------------------------------------------------------------
11c7b548 1064 if(fN<=0) return 0;
73042f01 1065 if (y <= fClusters[0]->GetY()) return 0;
1066 if (y > fClusters[fN-1]->GetY()) return fN;
1067 Int_t b=0, e=fN-1, m=(b+e)/2;
1068 for (; b<e; m=(b+e)/2) {
1069 if (y > fClusters[m]->GetY()) b=m+1;
1070 else e=m;
1071 }
1072 return m;
1073}
1074
1075//_____________________________________________________________________________
1076void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
1077 //-----------------------------------------------------------------
1078 // This funtion calculates dE/dX within the "low" and "up" cuts.
1079 //-----------------------------------------------------------------
1080 Int_t i;
be9c5115 1081 Int_t nc=GetNumberOfClusters();
bcc04e2a 1082
1083 Int_t swap;//stupid sorting
1084 do {
1085 swap=0;
1086 for (i=0; i<nc-1; i++) {
1087 if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
1088 Float_t tmp=fdEdxSample[i];
1089 fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
1090 swap++;
1091 }
1092 } while (swap);
73042f01 1093
be9c5115 1094 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
73042f01 1095 Float_t dedx=0;
bcc04e2a 1096 for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
73042f01 1097 dedx /= (nu-nl+1);
1098 SetdEdx(dedx);
7f6ddf58 1099
1100 //Very rough PID
1101 Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
1102
024a7fe9 1103 Double_t log1=TMath::Log(p+0.45), log2=TMath::Log(p+0.12);
7f6ddf58 1104 if (p<0.6) {
024a7fe9 1105 if (dedx < 34 + 30/(p+0.45)/(p+0.45) + 24*log1) {SetMass(0.13957); return;}
1106 if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.49368); return;}
7f6ddf58 1107 SetMass(0.93827); return;
1108 }
1109
1110 if (p<1.2) {
024a7fe9 1111 if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.13957); return;}
7f6ddf58 1112 SetMass(0.93827); return;
1113 }
1114
1115 SetMass(0.13957); return;
1116
73042f01 1117}
1118
73042f01 1119
1120