Configuration of pt-kick added.
[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$
19364939 18Revision 1.21 2002/10/21 16:49:46 hristov
19Replacing the old sorting by the Root one (M.Ivanov)
20
986921a6 21Revision 1.20 2002/10/14 14:57:43 hristov
22Merging the VirtualMC branch to the main development branch (HEAD)
23
b9d0a01d 24Revision 1.17.4.3 2002/10/11 08:34:48 hristov
25Updating VirtualMC to v3-09-02
26
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//-------------------------------------------------------
84
73042f01 85#include <TObjArray.h>
86#include <TFile.h>
be5b9287 87#include <TTree.h>
19364939 88#include <Riostream.h>
be5b9287 89
90#include "AliTPCtracker.h"
91#include "AliTPCcluster.h"
b9de75e1 92#include "AliTPCParam.h"
73042f01 93#include "AliTPCClustersRow.h"
94
b9de75e1 95//_____________________________________________________________________________
afc42102 96AliTPCtracker::AliTPCtracker(const AliTPCParam *par, Int_t eventn):
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
109 fN=0; fSectors=0;
110
111 fClustersArray.Setup(par);
112 fClustersArray.SetClusterType("AliTPCcluster");
b9de75e1 113
afc42102 114 char cname[100];
115 if (eventn==-1) {
116 sprintf(cname,"TreeC_TPC");
117 }
118 else {
119 sprintf(cname,"TreeC_TPC_%d",eventn);
120 }
121
122 fClustersArray.ConnectTree(cname);
123
124 fEventN = eventn;
b9de75e1 125 fSeeds=0;
126}
127
128//_____________________________________________________________________________
129AliTPCtracker::~AliTPCtracker() {
130 //------------------------------------------------------------------
131 // TPC tracker destructor
132 //------------------------------------------------------------------
133 delete[] fInnerSec;
134 delete[] fOuterSec;
1c8d9aad 135 if (fSeeds) {
136 fSeeds->Delete();
137 delete fSeeds;
138 }
b9de75e1 139}
73042f01 140
141//_____________________________________________________________________________
142Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
143{
144 //
145 // Parametrised error of the cluster reconstruction (pad direction)
146 //
147 // Sigma rphi
148 const Float_t kArphi=0.41818e-2;
149 const Float_t kBrphi=0.17460e-4;
150 const Float_t kCrphi=0.30993e-2;
151 const Float_t kDrphi=0.41061e-3;
152
153 pt=TMath::Abs(pt)*1000.;
154 Double_t x=r/pt;
155 tgl=TMath::Abs(tgl);
156 Double_t s=kArphi - kBrphi*r*tgl + kCrphi*x*x + kDrphi*x;
157 if (s<0.4e-3) s=0.4e-3;
158 s*=1.3; //Iouri Belikov
159
160 return s;
161}
162
163//_____________________________________________________________________________
164Double_t SigmaZ2(Double_t r, Double_t tgl)
165{
166 //
167 // Parametrised error of the cluster reconstruction (drift direction)
168 //
169 // Sigma z
170 const Float_t kAz=0.39614e-2;
171 const Float_t kBz=0.22443e-4;
172 const Float_t kCz=0.51504e-1;
173
174
175 tgl=TMath::Abs(tgl);
176 Double_t s=kAz - kBz*r*tgl + kCz*tgl*tgl;
177 if (s<0.4e-3) s=0.4e-3;
178 s*=1.3; //Iouri Belikov
179
180 return s;
181}
182
183//_____________________________________________________________________________
bdbd0f7a 184Double_t f1(Double_t x1,Double_t y1,
73042f01 185 Double_t x2,Double_t y2,
186 Double_t x3,Double_t y3)
187{
188 //-----------------------------------------------------------------
189 // Initial approximation of the track curvature
190 //-----------------------------------------------------------------
191 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
192 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
193 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
194 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
195 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
196
197 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
198
199 return -xr*yr/sqrt(xr*xr+yr*yr);
200}
201
202
203//_____________________________________________________________________________
bdbd0f7a 204Double_t f2(Double_t x1,Double_t y1,
73042f01 205 Double_t x2,Double_t y2,
206 Double_t x3,Double_t y3)
207{
208 //-----------------------------------------------------------------
209 // Initial approximation of the track curvature times center of curvature
210 //-----------------------------------------------------------------
211 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
212 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
213 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
214 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
215 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
216
217 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
218
219 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
220}
221
222//_____________________________________________________________________________
bdbd0f7a 223Double_t f3(Double_t x1,Double_t y1,
73042f01 224 Double_t x2,Double_t y2,
225 Double_t z1,Double_t z2)
226{
227 //-----------------------------------------------------------------
228 // Initial approximation of the tangent of the track dip angle
229 //-----------------------------------------------------------------
230 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
231}
232
233//_____________________________________________________________________________
b9de75e1 234void AliTPCtracker::LoadOuterSectors() {
235 //-----------------------------------------------------------------
236 // This function fills outer TPC sectors with clusters.
237 //-----------------------------------------------------------------
238 UInt_t index;
239 Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
240 for (Int_t i=0; i<j; i++) {
241 AliSegmentID *s=fClustersArray.LoadEntry(i);
242 Int_t sec,row;
243 AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
244 par->AdjustSectorRow(s->GetID(),sec,row);
245 if (sec<fkNIS*2) continue;
246 AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
247 Int_t ncl=clrow->GetArray()->GetEntriesFast();
248 while (ncl--) {
249 AliTPCcluster *c=(AliTPCcluster*)(*clrow)[ncl];
250 index=(((sec<<8)+row)<<16)+ncl;
251 fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index);
252 }
253 }
254
255 fN=fkNOS;
256 fSectors=fOuterSec;
257}
258
259//_____________________________________________________________________________
260void AliTPCtracker::UnloadOuterSectors() {
261 //-----------------------------------------------------------------
262 // This function clears outer TPC sectors.
263 //-----------------------------------------------------------------
264 Int_t nup=fOuterSec->GetNRows();
265 for (Int_t i=0; i<fkNOS; i++) {
266 for (Int_t j=0; j<nup; j++) {
267 if (fClustersArray.GetRow(i+fkNIS*2,j))
268 fClustersArray.ClearRow(i+fkNIS*2,j);
269 if (fClustersArray.GetRow(i+fkNIS*2+fkNOS,j))
270 fClustersArray.ClearRow(i+fkNIS*2+fkNOS,j);
271 }
272 }
273
274 fN=0;
275 fSectors=0;
276}
277
278//_____________________________________________________________________________
279void AliTPCtracker::LoadInnerSectors() {
280 //-----------------------------------------------------------------
281 // This function fills inner TPC sectors with clusters.
282 //-----------------------------------------------------------------
283 UInt_t index;
284 Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
285 for (Int_t i=0; i<j; i++) {
286 AliSegmentID *s=fClustersArray.LoadEntry(i);
287 Int_t sec,row;
288 AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
289 par->AdjustSectorRow(s->GetID(),sec,row);
290 if (sec>=fkNIS*2) continue;
291 AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
292 Int_t ncl=clrow->GetArray()->GetEntriesFast();
293 while (ncl--) {
294 AliTPCcluster *c=(AliTPCcluster*)(*clrow)[ncl];
295 index=(((sec<<8)+row)<<16)+ncl;
296 fInnerSec[sec%fkNIS][row].InsertCluster(c,index);
297 }
298 }
299
300 fN=fkNIS;
301 fSectors=fInnerSec;
302}
303
304//_____________________________________________________________________________
305void AliTPCtracker::UnloadInnerSectors() {
306 //-----------------------------------------------------------------
307 // This function clears inner TPC sectors.
308 //-----------------------------------------------------------------
309 Int_t nlow=fInnerSec->GetNRows();
310 for (Int_t i=0; i<fkNIS; i++) {
311 for (Int_t j=0; j<nlow; j++) {
312 if (fClustersArray.GetRow(i,j)) fClustersArray.ClearRow(i,j);
313 if (fClustersArray.GetRow(i+fkNIS,j)) fClustersArray.ClearRow(i+fkNIS,j);
314 }
315 }
316
317 fN=0;
318 fSectors=0;
319}
320
321//_____________________________________________________________________________
322Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
73042f01 323 //-----------------------------------------------------------------
324 // This function tries to find a track prolongation.
325 //-----------------------------------------------------------------
b9de75e1 326 Double_t xt=t.GetX();
327 const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip :
328 Int_t(0.5*fSectors->GetNRows());
73042f01 329 Int_t tryAgain=kSKIP;
73042f01 330
b9de75e1 331 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
332 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
333 if (alpha < 0. ) alpha += 2.*TMath::Pi();
334 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
335
f03e3423 336 Int_t nrows=fSectors->GetRowNumber(xt)-1;
337 for (Int_t nr=nrows; nr>=rf; nr--) {
b9de75e1 338 Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
73042f01 339 if (!t.PropagateTo(x)) return 0;
340
341 AliTPCcluster *cl=0;
342 UInt_t index=0;
b9de75e1 343 Double_t maxchi2=kMaxCHI2;
344 const AliTPCRow &krow=fSectors[s][nr];
345 Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
346 Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),pt);
73042f01 347 Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
348 Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
349
b9de75e1 350 if (road>kMaxROAD) {
73042f01 351 if (t.GetNumberOfClusters()>4)
352 cerr<<t.GetNumberOfClusters()
353 <<"FindProlongation warning: Too broad road !\n";
354 return 0;
355 }
356
357 if (krow) {
358 for (Int_t i=krow.Find(y-road); i<krow; i++) {
be5b9287 359 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
73042f01 360 if (c->GetY() > y+road) break;
361 if (c->IsUsed()) continue;
362 if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue;
363 Double_t chi2=t.GetPredictedChi2(c);
364 if (chi2 > maxchi2) continue;
365 maxchi2=chi2;
366 cl=c;
367 index=krow.GetIndex(i);
368 }
369 }
370 if (cl) {
b9de75e1 371 Float_t l=fSectors->GetPadPitchWidth();
024a7fe9 372 Float_t corr=1.; if (nr>63) corr=0.67; // new (third) pad response !
373 t.SetSampledEdx(cl->GetQ()/l*corr,t.GetNumberOfClusters());
be9c5115 374 if (!t.Update(cl,maxchi2,index)) {
375 if (!tryAgain--) return 0;
376 } else tryAgain=kSKIP;
73042f01 377 } else {
378 if (tryAgain==0) break;
379 if (y > ymax) {
b9de75e1 380 s = (s+1) % fN;
381 if (!t.Rotate(fSectors->GetAlpha())) return 0;
73042f01 382 } else if (y <-ymax) {
b9de75e1 383 s = (s-1+fN) % fN;
384 if (!t.Rotate(-fSectors->GetAlpha())) return 0;
73042f01 385 }
386 tryAgain--;
387 }
388 }
389
390 return 1;
b9de75e1 391}
392
393//_____________________________________________________________________________
394Int_t AliTPCtracker::FollowBackProlongation
395(AliTPCseed& seed, const AliTPCtrack &track) {
396 //-----------------------------------------------------------------
397 // This function propagates tracks back through the TPC
398 //-----------------------------------------------------------------
399 Double_t alpha=seed.GetAlpha() - fSectors->GetAlphaShift();
400 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
401 if (alpha < 0. ) alpha += 2.*TMath::Pi();
402 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
403
404 Int_t idx=-1, sec=-1, row=-1;
405 Int_t nc=seed.GetLabel(); //index of the cluster to start with
406 if (nc--) {
407 idx=track.GetClusterIndex(nc);
408 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
409 }
410 if (fSectors==fInnerSec) { if (sec >= 2*fkNIS) row=-1; }
411 else { if (sec < 2*fkNIS) row=-1; }
412
413 Int_t nr=fSectors->GetNRows();
414 for (Int_t i=0; i<nr; i++) {
415 Double_t x=fSectors->GetX(i), ymax=fSectors->GetMaxY(i);
416
417 if (!seed.PropagateTo(x)) return 0;
418
419 Double_t y=seed.GetY();
420 if (y > ymax) {
421 s = (s+1) % fN;
422 if (!seed.Rotate(fSectors->GetAlpha())) return 0;
423 } else if (y <-ymax) {
424 s = (s-1+fN) % fN;
425 if (!seed.Rotate(-fSectors->GetAlpha())) return 0;
426 }
427
428 AliTPCcluster *cl=0;
429 Int_t index=0;
430 Double_t maxchi2=kMaxCHI2;
431 Double_t pt=seed.GetConvConst()/(100/0.299792458/0.2)/seed.Get1Pt();
432 Double_t sy2=SigmaY2(seed.GetX(),seed.GetTgl(),pt);
433 Double_t sz2=SigmaZ2(seed.GetX(),seed.GetTgl());
434 Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ();
435 if (road>kMaxROAD) {
436 cerr<<seed.GetNumberOfClusters()
437 <<"AliTPCtracker::FollowBackProlongation: Too broad road !\n";
438 return 0;
439 }
440
441
442 Int_t accepted=seed.GetNumberOfClusters();
443 if (row==i) {
444 //try to accept already found cluster
445 AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
446 Double_t chi2;
447 if ((chi2=seed.GetPredictedChi2(c))<maxchi2 || accepted<27) {
448 index=idx; cl=c; maxchi2=chi2;
449 } //else cerr<<"AliTPCtracker::FollowBackProlongation: oulier !\n";
450
451 if (nc--) {
452 idx=track.GetClusterIndex(nc);
453 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
454 }
455 if (fSectors==fInnerSec) { if (sec >= 2*fkNIS) row=-1; }
456 else { if (sec < 2*fkNIS) row=-1; }
73042f01 457
b9de75e1 458 }
459 if (!cl) {
460 //try to fill the gap
461 const AliTPCRow &krow=fSectors[s][i];
462 if (accepted>27)
463 if (krow) {
464 for (Int_t i=krow.Find(y-road); i<krow; i++) {
465 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
466 if (c->GetY() > y+road) break;
467 if (c->IsUsed()) continue;
468 if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue;
469 Double_t chi2=seed.GetPredictedChi2(c);
470 if (chi2 > maxchi2) continue;
471 maxchi2=chi2;
472 cl=c;
473 index=krow.GetIndex(i);
474 }
475 }
476 }
477
478 if (cl) {
479 Float_t l=fSectors->GetPadPitchWidth();
024a7fe9 480 Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
481 seed.SetSampledEdx(cl->GetQ()/l*corr,seed.GetNumberOfClusters());
b9de75e1 482 seed.Update(cl,maxchi2,index);
483 }
484
485 }
486
487 seed.SetLabel(nc);
488
489 return 1;
73042f01 490}
491
492//_____________________________________________________________________________
b9de75e1 493void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
73042f01 494 //-----------------------------------------------------------------
495 // This function creates track seeds.
496 //-----------------------------------------------------------------
b9de75e1 497 if (fSeeds==0) fSeeds=new TObjArray(15000);
498
73042f01 499 Double_t x[5], c[15];
500
b9de75e1 501 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
73042f01 502 Double_t cs=cos(alpha), sn=sin(alpha);
503
b9de75e1 504 Double_t x1 =fOuterSec->GetX(i1);
505 Double_t xx2=fOuterSec->GetX(i2);
73042f01 506
b9de75e1 507 for (Int_t ns=0; ns<fkNOS; ns++) {
508 Int_t nl=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
509 Int_t nm=fOuterSec[ns][i2];
510 Int_t nu=fOuterSec[(ns+1)%fkNOS][i2];
511 const AliTPCRow& kr1=fOuterSec[ns][i1];
73042f01 512 for (Int_t is=0; is < kr1; is++) {
513 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
514 for (Int_t js=0; js < nl+nm+nu; js++) {
515 const AliTPCcluster *kcl;
516 Double_t x2, y2, z2;
7f6ddf58 517 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
73042f01 518
519 if (js<nl) {
b9de75e1 520 const AliTPCRow& kr2=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
73042f01 521 kcl=kr2[js];
522 y2=kcl->GetY(); z2=kcl->GetZ();
523 x2= xx2*cs+y2*sn;
524 y2=-xx2*sn+y2*cs;
525 } else
526 if (js<nl+nm) {
b9de75e1 527 const AliTPCRow& kr2=fOuterSec[ns][i2];
73042f01 528 kcl=kr2[js-nl];
529 x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
530 } else {
b9de75e1 531 const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][i2];
73042f01 532 kcl=kr2[js-nl-nm];
533 y2=kcl->GetY(); z2=kcl->GetZ();
534 x2=xx2*cs-y2*sn;
535 y2=xx2*sn+y2*cs;
536 }
537
7f6ddf58 538 Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2);
73042f01 539 if (TMath::Abs(zz-z2)>5.) continue;
540
541 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
542 if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
543
544 x[0]=y1;
545 x[1]=z1;
b9de75e1 546 x[4]=f1(x1,y1,x2,y2,x3,y3);
547 if (TMath::Abs(x[4]) >= 0.0066) continue;
be5b9287 548 x[2]=f2(x1,y1,x2,y2,x3,y3);
b9de75e1 549 //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
550 x[3]=f3(x1,y1,x2,y2,z1,z2);
551 if (TMath::Abs(x[3]) > 1.2) continue;
be5b9287 552 Double_t a=asin(x[2]);
b9de75e1 553 Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
7f6ddf58 554 if (TMath::Abs(zv-z3)>10.) continue;
73042f01 555
556 Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
557 Double_t sy2=kcl->GetSigmaY2(), sz2=kcl->GetSigmaZ2();
024a7fe9 558 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
559 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
73042f01 560
b9de75e1 561 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
562 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
563 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
be5b9287 564 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
565 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
b9de75e1 566 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
567 Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
568 Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
569 Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
570 Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
73042f01 571
572 c[0]=sy1;
573 c[1]=0.; c[2]=sz1;
b9de75e1 574 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
575 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
576 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
577 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
578 c[13]=f30*sy1*f40+f32*sy2*f42;
579 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
73042f01 580
581 UInt_t index=kr1.GetIndex(is);
582 AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
b9de75e1 583 Float_t l=fOuterSec->GetPadPitchWidth();
73042f01 584 track->SetSampledEdx(kr1[is]->GetQ()/l,0);
585
b9de75e1 586 Int_t rc=FollowProlongation(*track, i2);
be9c5115 587 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track;
b9de75e1 588 else fSeeds->AddLast(track);
73042f01 589 }
590 }
591 }
592}
593
594//_____________________________________________________________________________
b9de75e1 595Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
73042f01 596 //-----------------------------------------------------------------
b9de75e1 597 // This function reades track seeds.
73042f01 598 //-----------------------------------------------------------------
599 TDirectory *savedir=gDirectory;
600
b9de75e1 601 TFile *in=(TFile*)inp;
602 if (!in->IsOpen()) {
603 cerr<<"AliTPCtracker::ReadSeeds(): input file is not open !\n";
be5b9287 604 return 1;
73042f01 605 }
606
b9de75e1 607 in->cd();
608 TTree *seedTree=(TTree*)in->Get("Seeds");
609 if (!seedTree) {
610 cerr<<"AliTPCtracker::ReadSeeds(): ";
611 cerr<<"can't get a tree with track seeds !\n";
612 return 2;
613 }
614 AliTPCtrack *seed=new AliTPCtrack;
615 seedTree->SetBranchAddress("tracks",&seed);
616
617 if (fSeeds==0) fSeeds=new TObjArray(15000);
73042f01 618
b9de75e1 619 Int_t n=(Int_t)seedTree->GetEntries();
620 for (Int_t i=0; i<n; i++) {
621 seedTree->GetEvent(i);
622 fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
623 }
624
625 delete seed;
f38c8ae5 626
627 delete seedTree; //Thanks to Mariana Bondila
628
b9de75e1 629 savedir->cd();
73042f01 630
b9de75e1 631 return 0;
632}
73042f01 633
b9de75e1 634//_____________________________________________________________________________
635Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
636 //-----------------------------------------------------------------
637 // This is a track finder.
638 //-----------------------------------------------------------------
639 TDirectory *savedir=gDirectory;
73042f01 640
b9de75e1 641 if (inp) {
642 TFile *in=(TFile*)inp;
643 if (!in->IsOpen()) {
644 cerr<<"AliTPCtracker::Clusters2Tracks(): input file is not open !\n";
645 return 1;
646 }
647 }
73042f01 648
b9de75e1 649 if (!out->IsOpen()) {
650 cerr<<"AliTPCtracker::Clusters2Tracks(): output file is not open !\n";
651 return 2;
73042f01 652 }
653
b9de75e1 654 out->cd();
7f6ddf58 655
656 char tname[100];
657 sprintf(tname,"TreeT_TPC_%d",fEventN);
658 TTree tracktree(tname,"Tree with TPC tracks");
b9de75e1 659 AliTPCtrack *iotrack=0;
660 tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
661
662 LoadOuterSectors();
663
73042f01 664 //find track seeds
b9de75e1 665 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
73042f01 666 Int_t nrows=nlow+nup;
b9de75e1 667 if (fSeeds==0) {
668 Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
669 MakeSeeds(nup-1, nup-1-gap);
670 MakeSeeds(nup-1-shift, nup-1-shift-gap);
671 }
672 fSeeds->Sort();
73042f01 673
674 //tracking in outer sectors
b9de75e1 675 Int_t nseed=fSeeds->GetEntriesFast();
676 Int_t i;
73042f01 677 for (i=0; i<nseed; i++) {
b9de75e1 678 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
679 if (FollowProlongation(t)) {
680 UseClusters(&t);
73042f01 681 continue;
682 }
b9de75e1 683 delete fSeeds->RemoveAt(i);
73042f01 684 }
7f6ddf58 685 //UnloadOuterSectors();
73042f01 686
687 //tracking in inner sectors
b9de75e1 688 LoadInnerSectors();
73042f01 689 Int_t found=0;
690 for (i=0; i<nseed; i++) {
b9de75e1 691 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
73042f01 692 if (!pt) continue;
693 Int_t nc=t.GetNumberOfClusters();
694
b9de75e1 695 Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
73042f01 696 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
697 if (alpha < 0. ) alpha += 2.*TMath::Pi();
b9de75e1 698 Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
73042f01 699
b9de75e1 700 alpha=ns*fInnerSec->GetAlpha() + fInnerSec->GetAlphaShift() - t.GetAlpha();
73042f01 701
702 if (t.Rotate(alpha)) {
b9de75e1 703 if (FollowProlongation(t)) {
73042f01 704 if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
705 t.CookdEdx();
7f6ddf58 706 CookLabel(pt,0.1); //For comparison only
73042f01 707 iotrack=pt;
708 tracktree.Fill();
b9de75e1 709 UseClusters(&t,nc);
710 cerr<<found++<<'\r';
73042f01 711 }
712 }
713 }
b9de75e1 714 delete fSeeds->RemoveAt(i);
73042f01 715 }
b9de75e1 716 UnloadInnerSectors();
7f6ddf58 717 UnloadOuterSectors();
b9de75e1 718
7f6ddf58 719 tracktree.Write();
73042f01 720
b9de75e1 721 cerr<<"Number of found tracks : "<<found<<endl;
722
723 savedir->cd();
724
725 return 0;
726}
727
728//_____________________________________________________________________________
729Int_t AliTPCtracker::PropagateBack(const TFile *inp, TFile *out) {
730 //-----------------------------------------------------------------
731 // This function propagates tracks back through the TPC.
732 //-----------------------------------------------------------------
733 fSeeds=new TObjArray(15000);
734
735 TFile *in=(TFile*)inp;
736 TDirectory *savedir=gDirectory;
737
738 if (!in->IsOpen()) {
739 cerr<<"AliTPCtracker::PropagateBack(): ";
740 cerr<<"file with back propagated ITS tracks is not open !\n";
741 return 1;
73042f01 742 }
743
b9de75e1 744 if (!out->IsOpen()) {
745 cerr<<"AliTPCtracker::PropagateBack(): ";
746 cerr<<"file for back propagated TPC tracks is not open !\n";
747 return 2;
748 }
749
750 in->cd();
f38c8ae5 751 TTree *bckTree=(TTree*)in->Get("TreeT_ITSb_0");
b9de75e1 752 if (!bckTree) {
753 cerr<<"AliTPCtracker::PropagateBack() ";
754 cerr<<"can't get a tree with back propagated ITS tracks !\n";
755 return 3;
756 }
757 AliTPCtrack *bckTrack=new AliTPCtrack;
758 bckTree->SetBranchAddress("tracks",&bckTrack);
759
f38c8ae5 760 TTree *tpcTree=(TTree*)in->Get("TreeT_TPC_0");
b9de75e1 761 if (!tpcTree) {
762 cerr<<"AliTPCtracker::PropagateBack() ";
763 cerr<<"can't get a tree with TPC tracks !\n";
764 return 4;
765 }
766 AliTPCtrack *tpcTrack=new AliTPCtrack;
767 tpcTree->SetBranchAddress("tracks",&tpcTrack);
768
769//*** Prepare an array of tracks to be back propagated
770 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
771 Int_t nrows=nlow+nup;
772
773 TObjArray tracks(15000);
774 Int_t i=0,j=0;
775 Int_t tpcN=(Int_t)tpcTree->GetEntries();
776 Int_t bckN=(Int_t)bckTree->GetEntries();
777 if (j<bckN) bckTree->GetEvent(j++);
778 for (i=0; i<tpcN; i++) {
779 tpcTree->GetEvent(i);
780 Double_t alpha=tpcTrack->GetAlpha();
781 if (j<bckN &&
782 TMath::Abs(tpcTrack->GetLabel())==TMath::Abs(bckTrack->GetLabel())) {
783 if (!bckTrack->Rotate(alpha-bckTrack->GetAlpha())) continue;
784 fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
785 bckTree->GetEvent(j++);
786 } else {
787 tpcTrack->ResetCovariance();
788 fSeeds->AddLast(new AliTPCseed(*tpcTrack,alpha));
789 }
790 tracks.AddLast(new AliTPCtrack(*tpcTrack));
791 }
792
793 out->cd();
f38c8ae5 794 TTree backTree("TreeT_TPCb_0","Tree with back propagated TPC tracks");
b9de75e1 795 AliTPCtrack *otrack=0;
796 backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
797
798//*** Back propagation through inner sectors
799 LoadInnerSectors();
800 Int_t nseed=fSeeds->GetEntriesFast();
801 for (i=0; i<nseed; i++) {
802 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
803 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
73042f01 804
b9de75e1 805 Int_t nc=t.GetNumberOfClusters();
806 s.SetLabel(nc-1); //set number of the cluster to start with
73042f01 807
b9de75e1 808 if (FollowBackProlongation(s,t)) {
809 UseClusters(&s);
810 continue;
811 }
812 delete fSeeds->RemoveAt(i);
813 }
814 UnloadInnerSectors();
815
816//*** Back propagation through outer sectors
817 LoadOuterSectors();
818 Int_t found=0;
819 for (i=0; i<nseed; i++) {
820 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
821 if (!ps) continue;
822 Int_t nc=s.GetNumberOfClusters();
823 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
824
825 Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
826 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
827 if (alpha < 0. ) alpha += 2.*TMath::Pi();
828 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
829
830 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
831 alpha-=s.GetAlpha();
832
833 if (s.Rotate(alpha)) {
834 if (FollowBackProlongation(s,t)) {
835 if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
836 s.CookdEdx();
837 s.SetLabel(t.GetLabel());
838 UseClusters(&s,nc);
839 otrack=ps;
840 backTree.Fill();
841 cerr<<found++<<'\r';
842 continue;
843 }
844 }
845 }
846 delete fSeeds->RemoveAt(i);
847 }
848 UnloadOuterSectors();
849
850 backTree.Write();
73042f01 851 savedir->cd();
b9de75e1 852 cerr<<"Number of seeds: "<<nseed<<endl;
853 cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
854 cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
855
856 delete bckTrack;
857 delete tpcTrack;
be5b9287 858
f38c8ae5 859 delete bckTree; //Thanks to Mariana Bondila
860 delete tpcTree; //Thanks to Mariana Bondila
861
be5b9287 862 return 0;
73042f01 863}
864
865//_________________________________________________________________________
b9de75e1 866AliCluster *AliTPCtracker::GetCluster(Int_t index) const {
867 //--------------------------------------------------------------------
868 // Return pointer to a given cluster
869 //--------------------------------------------------------------------
870 Int_t sec=(index&0xff000000)>>24;
871 Int_t row=(index&0x00ff0000)>>16;
872 Int_t ncl=(index&0x0000ffff)>>00;
873
9e1a0ddb 874 AliTPCClustersRow *clrow=((AliTPCtracker *) this)->fClustersArray.GetRow(sec,row);
b9de75e1 875 return (AliCluster*)(*clrow)[ncl];
876}
877
878//__________________________________________________________________________
879void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
880 //--------------------------------------------------------------------
881 //This function "cooks" a track label. If label<0, this track is fake.
882 //--------------------------------------------------------------------
883 Int_t noc=t->GetNumberOfClusters();
884 Int_t *lb=new Int_t[noc];
885 Int_t *mx=new Int_t[noc];
886 AliCluster **clusters=new AliCluster*[noc];
887
888 Int_t i;
889 for (i=0; i<noc; i++) {
890 lb[i]=mx[i]=0;
891 Int_t index=t->GetClusterIndex(i);
892 clusters[i]=GetCluster(index);
893 }
894
895 Int_t lab=123456789;
896 for (i=0; i<noc; i++) {
897 AliCluster *c=clusters[i];
898 lab=TMath::Abs(c->GetLabel(0));
899 Int_t j;
900 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
901 lb[j]=lab;
902 (mx[j])++;
903 }
904
905 Int_t max=0;
906 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
907
908 for (i=0; i<noc; i++) {
909 AliCluster *c=clusters[i];
910 if (TMath::Abs(c->GetLabel(1)) == lab ||
911 TMath::Abs(c->GetLabel(2)) == lab ) max++;
912 }
913
914 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
915
916 else {
917 Int_t tail=Int_t(0.10*noc);
918 max=0;
919 for (i=1; i<=tail; i++) {
920 AliCluster *c=clusters[noc-i];
921 if (lab == TMath::Abs(c->GetLabel(0)) ||
922 lab == TMath::Abs(c->GetLabel(1)) ||
923 lab == TMath::Abs(c->GetLabel(2))) max++;
924 }
925 if (max < Int_t(0.5*tail)) lab=-lab;
926 }
927
928 t->SetLabel(lab);
929
930 delete[] lb;
931 delete[] mx;
932 delete[] clusters;
933}
934
935//_________________________________________________________________________
936void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
937 //-----------------------------------------------------------------------
938 // Setup inner sector
939 //-----------------------------------------------------------------------
940 if (f==0) {
941 fAlpha=par->GetInnerAngle();
942 fAlphaShift=par->GetInnerAngleShift();
943 fPadPitchWidth=par->GetInnerPadPitchWidth();
024a7fe9 944 f1PadPitchLength=par->GetInnerPadPitchLength();
945 f2PadPitchLength=f1PadPitchLength;
b9de75e1 946
947 fN=par->GetNRowLow();
948 fRow=new AliTPCRow[fN];
949 for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
950 } else {
951 fAlpha=par->GetOuterAngle();
952 fAlphaShift=par->GetOuterAngleShift();
953 fPadPitchWidth=par->GetOuterPadPitchWidth();
f03e3423 954 f1PadPitchLength=par->GetOuter1PadPitchLength();
955 f2PadPitchLength=par->GetOuter2PadPitchLength();
b9de75e1 956
957 fN=par->GetNRowUp();
958 fRow=new AliTPCRow[fN];
f03e3423 959 for (Int_t i=0; i<fN; i++){
960 fRow[i].SetX(par->GetPadRowRadiiUp(i));
961 }
b9de75e1 962 }
963}
964
965//_________________________________________________________________________
73042f01 966void
967AliTPCtracker::AliTPCRow::InsertCluster(const AliTPCcluster* c, UInt_t index) {
968 //-----------------------------------------------------------------------
969 // Insert a cluster into this pad row in accordence with its y-coordinate
970 //-----------------------------------------------------------------------
b9de75e1 971 if (fN==kMaxClusterPerRow) {
73042f01 972 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
973 }
974 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
975 Int_t i=Find(c->GetY());
976 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
977 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
978 fIndex[i]=index; fClusters[i]=c; fN++;
979}
980
981//___________________________________________________________________
982Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
983 //-----------------------------------------------------------------------
984 // Return the index of the nearest cluster
985 //-----------------------------------------------------------------------
11c7b548 986 if(fN<=0) return 0;
73042f01 987 if (y <= fClusters[0]->GetY()) return 0;
988 if (y > fClusters[fN-1]->GetY()) return fN;
989 Int_t b=0, e=fN-1, m=(b+e)/2;
990 for (; b<e; m=(b+e)/2) {
991 if (y > fClusters[m]->GetY()) b=m+1;
992 else e=m;
993 }
994 return m;
995}
996
997//_____________________________________________________________________________
998void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
999 //-----------------------------------------------------------------
1000 // This funtion calculates dE/dX within the "low" and "up" cuts.
1001 //-----------------------------------------------------------------
1002 Int_t i;
be9c5115 1003 Int_t nc=GetNumberOfClusters();
986921a6 1004 Int_t * index = new Int_t[nc];
1005 TMath::Sort(nc, fdEdxSample,index,kFALSE);
73042f01 1006
be9c5115 1007 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
73042f01 1008 Float_t dedx=0;
986921a6 1009 for (i=nl; i<=nu; i++) dedx += fdEdxSample[index[i]];
73042f01 1010 dedx /= (nu-nl+1);
1011 SetdEdx(dedx);
7f6ddf58 1012
986921a6 1013 delete [] index;
1014
7f6ddf58 1015 //Very rough PID
1016 Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
1017
024a7fe9 1018 Double_t log1=TMath::Log(p+0.45), log2=TMath::Log(p+0.12);
7f6ddf58 1019 if (p<0.6) {
024a7fe9 1020 if (dedx < 34 + 30/(p+0.45)/(p+0.45) + 24*log1) {SetMass(0.13957); return;}
1021 if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.49368); return;}
7f6ddf58 1022 SetMass(0.93827); return;
1023 }
1024
1025 if (p<1.2) {
024a7fe9 1026 if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.13957); return;}
7f6ddf58 1027 SetMass(0.93827); return;
1028 }
1029
1030 SetMass(0.13957); return;
1031
73042f01 1032}
1033
73042f01 1034
1035