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