]> 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
Add support for 64 bits
[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$
18Revision 1.35.2.3 2003/07/15 09:58:03 hristov
19Corrected back-propagation (Yu.Belikov)
20
21Revision 1.35.2.2 2003/07/14 09:19:33 hristov
22TOF included in the combined PID (Yu.Belikov)
23
24Revision 1.35.2.1 2003/07/11 10:53:01 hristov
25Inward refit for TPC and TRD in the ESD schema (T.Kuhr)
26
27Revision 1.35 2003/05/23 10:08:51 hristov
28SetLabel replaced by SetNumber (Yu.Belikov)
29
30Revision 1.34 2003/05/22 13:57:48 hristov
31First implementation of ESD classes (Yu.Belikov)
32
33Revision 1.32 2003/04/10 10:36:54 hristov
34Code for unified TPC/TRD tracking (S.Radomski)
35
36Revision 1.31 2003/03/19 17:14:11 hristov
37Load/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)
38
39Revision 1.30 2003/02/28 16:13:32 hristov
40Typos corrected
41
42Revision 1.29 2003/02/28 15:18:16 hristov
43Corrections suggested by J.Chudoba
44
45Revision 1.28 2003/02/27 16:15:52 hristov
46Code for inward refitting (S.Radomski)
47
48Revision 1.27 2003/02/25 16:47:58 hristov
49allow back propagation for more than 1 event (J.Chudoba)
50
51Revision 1.26 2003/02/19 08:49:46 hristov
52Track time measurement (S.Radomski)
53
54Revision 1.25 2003/01/28 16:43:35 hristov
55Additional protection: to be discussed with the Root team (M.Ivanov)
56
57Revision 1.24 2002/11/19 16:13:24 hristov
58stdlib.h included to declare exit() on HP
59
60Revision 1.23 2002/11/19 11:50:08 hristov
61Removing CONTAINERS (Yu.Belikov)
62
63Revision 1.19 2002/07/19 07:31:40 kowal2
64Improvement in tracking by J. Belikov
65
66Revision 1.18 2002/05/13 07:33:52 kowal2
67Added protection in Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const
68in the case of defined region of interests.
69
70Revision 1.17 2002/03/18 17:59:13 kowal2
71Chnges in the pad geometry - 3 pad lengths introduced.
72
73Revision 1.16 2001/11/08 16:39:03 hristov
74Additional protection (M.Masera)
75
76Revision 1.15 2001/11/08 16:36:33 hristov
77Updated 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)
78
79Revision 1.14 2001/10/21 19:04:55 hristov
80Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
81
82Revision 1.13 2001/07/23 12:01:30 hristov
83Initialisation added
84
85Revision 1.12 2001/07/20 14:32:44 kowal2
86Processing of many events possible now
87
88Revision 1.11 2001/05/23 08:50:10 hristov
89Weird inline removed
90
91Revision 1.10 2001/05/16 14:57:25 alibrary
92New files for folders and Stack
93
94Revision 1.9 2001/05/11 07:16:56 hristov
95Fix needed on Sun and Alpha
96
97Revision 1.8 2001/05/08 15:00:15 hristov
98Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
99
100Revision 1.5 2000/12/20 07:51:59 kowal2
101Changes suggested by Alessandra and Paolo to avoid overlapped
102data fields in encapsulated classes.
103
104Revision 1.4 2000/11/02 07:27:16 kowal2
105code corrections
106
107Revision 1.2 2000/06/30 12:07:50 kowal2
108Updated from the TPC-PreRelease branch
109
110Revision 1.1.2.1 2000/06/25 08:53:55 kowal2
111Splitted from AliTPCtracking
112
113*/
73042f01 114
115//-------------------------------------------------------
116// Implementation of the TPC tracker
117//
118// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
119//-------------------------------------------------------
73042f01 120#include <TObjArray.h>
c630aafd 121#include <TError.h>
73042f01 122#include <TFile.h>
be5b9287 123#include <TTree.h>
c630aafd 124#include <stdlib.h>
125
126#include "AliESD.h"
be5b9287 127
128#include "AliTPCtracker.h"
129#include "AliTPCcluster.h"
b9de75e1 130#include "AliTPCParam.h"
c630aafd 131#include "AliClusters.h"
132
133ClassImp(AliTPCtracker)
83d73500 134
b9de75e1 135//_____________________________________________________________________________
c630aafd 136AliTPCtracker::AliTPCtracker(const AliTPCParam *par):
137AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
b9de75e1 138{
139 //---------------------------------------------------------------------
140 // The main TPC tracker constructor
141 //---------------------------------------------------------------------
142 fInnerSec=new AliTPCSector[fkNIS];
143 fOuterSec=new AliTPCSector[fkNOS];
144
145 Int_t i;
146 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
147 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
148
c630aafd 149 fParam = (AliTPCParam*) par;
150 fSeeds=0;
88cb7938 151
c630aafd 152 // [SR 17.03.2003]
153
154 fBarrelFile = 0;
155 fBarrelTree = 0;
156 fBarrelArray = 0;
157 fBarrelTrack = 0;
9c9d2487 158
b9de75e1 159}
160
161//_____________________________________________________________________________
162AliTPCtracker::~AliTPCtracker() {
163 //------------------------------------------------------------------
164 // TPC tracker destructor
165 //------------------------------------------------------------------
166 delete[] fInnerSec;
167 delete[] fOuterSec;
1c8d9aad 168 if (fSeeds) {
169 fSeeds->Delete();
170 delete fSeeds;
171 }
c630aafd 172
173 // [SR, 01.04.2003]
174 if (fBarrelFile) {
175 fBarrelFile->Close();
176 delete fBarrelFile;
177 }
178}
179
180//_____________________________________________________________________________
181
182void AliTPCtracker::SetBarrelTree(const char *mode) {
183 //
184 // Creates a tree for BarrelTracks
185 // mode = "back" or "refit"
186 //
187 // [SR, 01.04.2003]
188 //
189
190 if (!IsStoringBarrel()) return;
191
192 TDirectory *sav = gDirectory;
193 if (!fBarrelFile) fBarrelFile = new TFile("AliBarrelTracks.root", "UPDATE");
194
195 char buff[40];
196 sprintf(buff, "BarrelTPC_%d_%s", GetEventNumber(), mode);
197
198 fBarrelFile->cd();
199 fBarrelTree = new TTree(buff, "Barrel TPC tracks");
200
201 if (!fBarrelArray) fBarrelArray = new TClonesArray("AliBarrelTrack", 4);
202 for(Int_t i=0; i<4; i++) new((*fBarrelArray)[i]) AliBarrelTrack();
203
204 fBarrelTree->Branch("tracks", &fBarrelArray);
205
206 sav->cd();
207}
208//_____________________________________________________________________________
209
210void AliTPCtracker::StoreBarrelTrack(AliTPCtrack *ps, Int_t refPlane, Int_t isIn) {
211 //
212 // Stores Track at a given reference plane
213 //
214 // refPlane: 1-4
215 // isIn: 1 - backward, 2 - refit
216 //
217
218 if (!IsStoringBarrel()) return;
219 if (refPlane < 0 || refPlane > 4) return;
220 if (isIn > 2) return;
221
222 static Int_t nClusters;
223 static Int_t nWrong;
224 static Double_t chi2;
225 static Int_t index;
226
227 Int_t newClusters, newWrong;
228 Double_t newChi2;
229
230 if ( (refPlane == 1 && isIn == kTrackBack) ||
231 (refPlane == 4 && isIn == kTrackRefit) ) {
232
233 fBarrelArray->Clear();
234 nClusters = nWrong = 0;
235 chi2 = 0.0;
236 index = 0;
237 }
238
239 // propagate
240 Double_t refX = 0;
241 if (refPlane == 1) refX = fParam->GetInnerRadiusLow();
242 if (refPlane == 2) refX = fParam->GetInnerRadiusUp();
243 if (refPlane == 3) refX = fParam->GetOuterRadiusLow();
244 if (refPlane == 4) refX = fParam->GetOuterRadiusUp();
245
246 ps->PropagateTo(refX);
247
248 fBarrelTrack = (AliBarrelTrack*)(*fBarrelArray)[index++];
249 ps->GetBarrelTrack(fBarrelTrack);
250
251 newClusters = ps->GetNumberOfClusters() - nClusters;
252 newWrong = ps->GetNWrong() - nWrong;
253 newChi2 = ps->GetChi2() - chi2;
254
255 nClusters = ps->GetNumberOfClusters();
256 nWrong = ps->GetNWrong();
257 chi2 = ps->GetChi2();
258
259 fBarrelTrack->SetNClusters(newClusters, newChi2);
260 fBarrelTrack->SetNWrongClusters(newWrong);
261 fBarrelTrack->SetRefPlane(refPlane, isIn);
b9de75e1 262}
73042f01 263
264//_____________________________________________________________________________
265Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
266{
267 //
268 // Parametrised error of the cluster reconstruction (pad direction)
269 //
270 // Sigma rphi
271 const Float_t kArphi=0.41818e-2;
272 const Float_t kBrphi=0.17460e-4;
273 const Float_t kCrphi=0.30993e-2;
274 const Float_t kDrphi=0.41061e-3;
275
276 pt=TMath::Abs(pt)*1000.;
277 Double_t x=r/pt;
278 tgl=TMath::Abs(tgl);
279 Double_t s=kArphi - kBrphi*r*tgl + kCrphi*x*x + kDrphi*x;
280 if (s<0.4e-3) s=0.4e-3;
281 s*=1.3; //Iouri Belikov
282
283 return s;
284}
285
286//_____________________________________________________________________________
287Double_t SigmaZ2(Double_t r, Double_t tgl)
288{
289 //
290 // Parametrised error of the cluster reconstruction (drift direction)
291 //
292 // Sigma z
293 const Float_t kAz=0.39614e-2;
294 const Float_t kBz=0.22443e-4;
295 const Float_t kCz=0.51504e-1;
296
297
298 tgl=TMath::Abs(tgl);
299 Double_t s=kAz - kBz*r*tgl + kCz*tgl*tgl;
300 if (s<0.4e-3) s=0.4e-3;
301 s*=1.3; //Iouri Belikov
302
303 return s;
304}
305
306//_____________________________________________________________________________
bdbd0f7a 307Double_t f1(Double_t x1,Double_t y1,
73042f01 308 Double_t x2,Double_t y2,
309 Double_t x3,Double_t y3)
310{
311 //-----------------------------------------------------------------
312 // Initial approximation of the track curvature
313 //-----------------------------------------------------------------
314 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
315 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
316 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
317 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
318 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
319
320 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
321
322 return -xr*yr/sqrt(xr*xr+yr*yr);
323}
324
325
326//_____________________________________________________________________________
bdbd0f7a 327Double_t f2(Double_t x1,Double_t y1,
73042f01 328 Double_t x2,Double_t y2,
329 Double_t x3,Double_t y3)
330{
331 //-----------------------------------------------------------------
332 // Initial approximation of the track curvature times center of curvature
333 //-----------------------------------------------------------------
334 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
335 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
336 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
337 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
338 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
339
340 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
341
342 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
343}
344
345//_____________________________________________________________________________
bdbd0f7a 346Double_t f3(Double_t x1,Double_t y1,
73042f01 347 Double_t x2,Double_t y2,
348 Double_t z1,Double_t z2)
349{
350 //-----------------------------------------------------------------
351 // Initial approximation of the tangent of the track dip angle
352 //-----------------------------------------------------------------
353 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
354}
355
83d73500 356//_____________________________________________________________________________
c630aafd 357Int_t AliTPCtracker::LoadClusters(TTree *cTree) {
b9de75e1 358 //-----------------------------------------------------------------
c630aafd 359 // This function loads TPC clusters.
b9de75e1 360 //-----------------------------------------------------------------
c630aafd 361 TBranch *branch=cTree->GetBranch("Segment");
362 if (!branch) {
363 Error("LoadClusters","Can't get the branch !");
364 return 1;
b9de75e1 365 }
366
c630aafd 367 AliClusters carray, *addr=&carray;
368 carray.SetClass("AliTPCcluster");
369 carray.SetArray(0);
370 branch->SetAddress(&addr);
b9de75e1 371
c630aafd 372 Int_t nentr=(Int_t)cTree->GetEntries();
83d73500 373
c630aafd 374 for (Int_t i=0; i<nentr; i++) {
375 cTree->GetEvent(i);
bcc04e2a 376
c630aafd 377 Int_t ncl=carray.GetArray()->GetEntriesFast();
378
379 Int_t nir=fInnerSec->GetNRows(), nor=fOuterSec->GetNRows();
380 Int_t id=carray.GetID();
381 if ((id<0) || (id>2*(fkNIS*nir + fkNOS*nor))) {
382 Error("LoadClusters","Wrong index !");
383 exit(1);
bcc04e2a 384 }
c630aafd 385 Int_t outindex = 2*fkNIS*nir;
386 if (id<outindex) {
387 Int_t sec = id/nir;
388 Int_t row = id - sec*nir;
389 sec %= fkNIS;
390 AliTPCRow &padrow=fInnerSec[sec][row];
391 while (ncl--) {
392 AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
393 padrow.InsertCluster(c,sec,row);
394 }
395 } else {
396 id -= outindex;
397 Int_t sec = id/nor;
398 Int_t row = id - sec*nor;
399 sec %= fkNOS;
400 AliTPCRow &padrow=fOuterSec[sec][row];
401 while (ncl--) {
402 AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
403 padrow.InsertCluster(c,sec+fkNIS,row);
404 }
405 }
406 carray.GetArray()->Clear();
bcc04e2a 407 }
88cb7938 408
c630aafd 409 return 0;
b9de75e1 410}
411
412//_____________________________________________________________________________
c630aafd 413void AliTPCtracker::UnloadClusters() {
b9de75e1 414 //-----------------------------------------------------------------
c630aafd 415 // This function unloads TPC clusters.
b9de75e1 416 //-----------------------------------------------------------------
c630aafd 417 Int_t i;
418 for (i=0; i<fkNIS; i++) {
419 Int_t nr=fInnerSec->GetNRows();
420 for (Int_t n=0; n<nr; n++) fInnerSec[i][n].ResetClusters();
421 }
422 for (i=0; i<fkNOS; i++) {
423 Int_t nr=fOuterSec->GetNRows();
424 for (Int_t n=0; n<nr; n++) fOuterSec[i][n].ResetClusters();
b9de75e1 425 }
b9de75e1 426}
427
428//_____________________________________________________________________________
429Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
73042f01 430 //-----------------------------------------------------------------
431 // This function tries to find a track prolongation.
432 //-----------------------------------------------------------------
b9de75e1 433 Double_t xt=t.GetX();
434 const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip :
435 Int_t(0.5*fSectors->GetNRows());
73042f01 436 Int_t tryAgain=kSKIP;
73042f01 437
b9de75e1 438 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
439 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
440 if (alpha < 0. ) alpha += 2.*TMath::Pi();
441 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
442
f03e3423 443 Int_t nrows=fSectors->GetRowNumber(xt)-1;
444 for (Int_t nr=nrows; nr>=rf; nr--) {
b9de75e1 445 Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
73042f01 446 if (!t.PropagateTo(x)) return 0;
447
448 AliTPCcluster *cl=0;
449 UInt_t index=0;
b9de75e1 450 Double_t maxchi2=kMaxCHI2;
451 const AliTPCRow &krow=fSectors[s][nr];
452 Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
453 Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),pt);
73042f01 454 Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
455 Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
456
b9de75e1 457 if (road>kMaxROAD) {
73042f01 458 if (t.GetNumberOfClusters()>4)
c630aafd 459 Warning("FindProlongation","Too broad road !");
73042f01 460 return 0;
461 }
462
463 if (krow) {
464 for (Int_t i=krow.Find(y-road); i<krow; i++) {
c630aafd 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.*(t.GetSigmaZ2()+sz2)) continue;
469 Double_t chi2=t.GetPredictedChi2(c);
470 if (chi2 > maxchi2) continue;
471 maxchi2=chi2;
472 cl=c;
73042f01 473 index=krow.GetIndex(i);
474 }
475 }
476 if (cl) {
b9de75e1 477 Float_t l=fSectors->GetPadPitchWidth();
024a7fe9 478 Float_t corr=1.; if (nr>63) corr=0.67; // new (third) pad response !
479 t.SetSampledEdx(cl->GetQ()/l*corr,t.GetNumberOfClusters());
be9c5115 480 if (!t.Update(cl,maxchi2,index)) {
481 if (!tryAgain--) return 0;
482 } else tryAgain=kSKIP;
73042f01 483 } else {
484 if (tryAgain==0) break;
485 if (y > ymax) {
b9de75e1 486 s = (s+1) % fN;
487 if (!t.Rotate(fSectors->GetAlpha())) return 0;
73042f01 488 } else if (y <-ymax) {
b9de75e1 489 s = (s-1+fN) % fN;
490 if (!t.Rotate(-fSectors->GetAlpha())) return 0;
73042f01 491 }
492 tryAgain--;
493 }
494 }
495
496 return 1;
b9de75e1 497}
c630aafd 498//_____________________________________________________________________________
499
500Int_t AliTPCtracker::FollowRefitInward(AliTPCseed *seed, AliTPCtrack *track) {
501 //
502 // This function propagates seed inward TPC using old clusters
503 // from the track.
504 //
505 // Sylwester Radomski, GSI
506 // 26.02.2003
507 //
508
509 // loop over rows
510
511 Int_t nRows = fSectors->GetNRows();
512 for (Int_t iRow = nRows; iRow >= 0; iRow--) {
513
514 Double_t x = fSectors->GetX(iRow);
515 if (!seed->PropagateTo(x)) return 0;
516
517 // try to find an assigned cluster in this row
518
519 AliTPCcluster* cluster = NULL;
520 Int_t idx = -1;
521 Int_t sec = -1;
522 for (Int_t iCluster = 0; iCluster < track->GetNumberOfClusters(); iCluster++) {
523 idx = track->GetClusterIndex(iCluster);
524 sec = (idx&0xff000000)>>24;
525 Int_t row = (idx&0x00ff0000)>>16;
526 if (((fSectors == fInnerSec) && (sec >= fkNIS)) ||
527 ((fSectors == fOuterSec) && (sec < fkNIS))) continue;
528 if (row == iRow) {
529 cluster = (AliTPCcluster*) GetCluster(idx);
530 break;
531 }
532 }
533
534 // update the track seed with the found cluster
535
536 if (cluster) {
537 Double_t dAlpha = fParam->GetAngle(sec) - seed->GetAlpha();
538 if (TMath::Abs(dAlpha) > 0.0001) {
539 if (!seed->Rotate(dAlpha)) return 0;
540 if (!seed->PropagateTo(x)) return 0;
541 }
542
543 seed->Update(cluster, seed->GetPredictedChi2(cluster), idx);
544 }
545 }
546
547 return 1;
548}
b9de75e1 549
550//_____________________________________________________________________________
551Int_t AliTPCtracker::FollowBackProlongation
552(AliTPCseed& seed, const AliTPCtrack &track) {
553 //-----------------------------------------------------------------
554 // This function propagates tracks back through the TPC
555 //-----------------------------------------------------------------
556 Double_t alpha=seed.GetAlpha() - fSectors->GetAlphaShift();
557 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
558 if (alpha < 0. ) alpha += 2.*TMath::Pi();
559 Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
560
561 Int_t idx=-1, sec=-1, row=-1;
c630aafd 562 Int_t nc=seed.GetNumber();
563
b9de75e1 564 if (nc--) {
565 idx=track.GetClusterIndex(nc);
566 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
567 }
c630aafd 568 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
569 else { if (sec < fkNIS) row=-1; }
b9de75e1 570
571 Int_t nr=fSectors->GetNRows();
572 for (Int_t i=0; i<nr; i++) {
573 Double_t x=fSectors->GetX(i), ymax=fSectors->GetMaxY(i);
c630aafd 574 Double_t y=seed.GetYat(x);
575
b9de75e1 576 if (y > ymax) {
577 s = (s+1) % fN;
578 if (!seed.Rotate(fSectors->GetAlpha())) return 0;
579 } else if (y <-ymax) {
580 s = (s-1+fN) % fN;
581 if (!seed.Rotate(-fSectors->GetAlpha())) return 0;
582 }
583
c630aafd 584 if (!seed.PropagateTo(x)) return 0;
585
b9de75e1 586 AliTPCcluster *cl=0;
587 Int_t index=0;
588 Double_t maxchi2=kMaxCHI2;
589 Double_t pt=seed.GetConvConst()/(100/0.299792458/0.2)/seed.Get1Pt();
590 Double_t sy2=SigmaY2(seed.GetX(),seed.GetTgl(),pt);
591 Double_t sz2=SigmaZ2(seed.GetX(),seed.GetTgl());
592 Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ();
593 if (road>kMaxROAD) {
c630aafd 594 Warning("FollowBackProlongation","Too broad road !");
b9de75e1 595 return 0;
596 }
597
b9de75e1 598 Int_t accepted=seed.GetNumberOfClusters();
599 if (row==i) {
600 //try to accept already found cluster
601 AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
602 Double_t chi2;
603 if ((chi2=seed.GetPredictedChi2(c))<maxchi2 || accepted<27) {
604 index=idx; cl=c; maxchi2=chi2;
605 } //else cerr<<"AliTPCtracker::FollowBackProlongation: oulier !\n";
606
607 if (nc--) {
608 idx=track.GetClusterIndex(nc);
609 sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
610 }
c630aafd 611 if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
612 else { if (sec < fkNIS) row=-1; }
73042f01 613
b9de75e1 614 }
615 if (!cl) {
616 //try to fill the gap
617 const AliTPCRow &krow=fSectors[s][i];
618 if (accepted>27)
619 if (krow) {
620 for (Int_t i=krow.Find(y-road); i<krow; i++) {
c630aafd 621 AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
622 if (c->GetY() > y+road) break;
623 if (c->IsUsed()) continue;
624 if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue;
625 Double_t chi2=seed.GetPredictedChi2(c);
626 if (chi2 > maxchi2) continue;
627 maxchi2=chi2;
628 cl=c;
b9de75e1 629 index=krow.GetIndex(i);
630 }
631 }
632 }
633
634 if (cl) {
635 Float_t l=fSectors->GetPadPitchWidth();
024a7fe9 636 Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
637 seed.SetSampledEdx(cl->GetQ()/l*corr,seed.GetNumberOfClusters());
b9de75e1 638 seed.Update(cl,maxchi2,index);
639 }
640
641 }
642
c630aafd 643 seed.SetNumber(nc);
b9de75e1 644
645 return 1;
73042f01 646}
647
648//_____________________________________________________________________________
b9de75e1 649void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
73042f01 650 //-----------------------------------------------------------------
651 // This function creates track seeds.
652 //-----------------------------------------------------------------
653 Double_t x[5], c[15];
654
c630aafd 655 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
73042f01 656 Double_t cs=cos(alpha), sn=sin(alpha);
657
c630aafd 658 Double_t x1 =fSectors->GetX(i1);
659 Double_t xx2=fSectors->GetX(i2);
73042f01 660
c630aafd 661 for (Int_t ns=0; ns<fN; ns++) {
662 Int_t nl=fSectors[(ns-1+fN)%fN][i2];
663 Int_t nm=fSectors[ns][i2];
664 Int_t nu=fSectors[(ns+1)%fN][i2];
665 const AliTPCRow& kr1=fSectors[ns][i1];
73042f01 666 for (Int_t is=0; is < kr1; is++) {
667 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
668 for (Int_t js=0; js < nl+nm+nu; js++) {
c630aafd 669 const AliTPCcluster *kcl;
73042f01 670 Double_t x2, y2, z2;
7f6ddf58 671 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
73042f01 672
c630aafd 673 if (js<nl) {
674 const AliTPCRow& kr2=fSectors[(ns-1+fN)%fN][i2];
675 kcl=kr2[js];
73042f01 676 y2=kcl->GetY(); z2=kcl->GetZ();
677 x2= xx2*cs+y2*sn;
678 y2=-xx2*sn+y2*cs;
c630aafd 679 } else
680 if (js<nl+nm) {
681 const AliTPCRow& kr2=fSectors[ns][i2];
682 kcl=kr2[js-nl];
73042f01 683 x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
c630aafd 684 } else {
685 const AliTPCRow& kr2=fSectors[(ns+1)%fN][i2];
686 kcl=kr2[js-nl-nm];
73042f01 687 y2=kcl->GetY(); z2=kcl->GetZ();
688 x2=xx2*cs-y2*sn;
689 y2=xx2*sn+y2*cs;
c630aafd 690 }
73042f01 691
7f6ddf58 692 Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2);
73042f01 693 if (TMath::Abs(zz-z2)>5.) continue;
694
695 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
c630aafd 696 if (d==0.) {
697 Warning("MakeSeeds","Straight seed !");
698 continue;
699 }
700 x[0]=y1;
701 x[1]=z1;
702 x[4]=f1(x1,y1,x2,y2,x3,y3);
703 if (TMath::Abs(x[4]) >= 0.0066) continue;
704 x[2]=f2(x1,y1,x2,y2,x3,y3);
705 //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
706 x[3]=f3(x1,y1,x2,y2,z1,z2);
707 if (TMath::Abs(x[3]) > 1.2) continue;
708 Double_t a=asin(x[2]);
709 Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
710 if (TMath::Abs(zv-z3)>10.) continue;
73042f01 711
712 Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
713 Double_t sy2=kcl->GetSigmaY2(), sz2=kcl->GetSigmaZ2();
024a7fe9 714 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
715 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
73042f01 716
c630aafd 717 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
718 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
719 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
720 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
721 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
722 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
723 Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
724 Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
725 Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
726 Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
73042f01 727
728 c[0]=sy1;
729 c[1]=0.; c[2]=sz1;
b9de75e1 730 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
731 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
732 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
733 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
734 c[13]=f30*sy1*f40+f32*sy2*f42;
735 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
73042f01 736
737 UInt_t index=kr1.GetIndex(is);
c630aafd 738 AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
739 Float_t l=fSectors->GetPadPitchWidth();
73042f01 740 track->SetSampledEdx(kr1[is]->GetQ()/l,0);
741
b9de75e1 742 Int_t rc=FollowProlongation(*track, i2);
be9c5115 743 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track;
c630aafd 744 else fSeeds->AddLast(track);
73042f01 745 }
746 }
747 }
748}
749
750//_____________________________________________________________________________
b9de75e1 751Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
73042f01 752 //-----------------------------------------------------------------
b9de75e1 753 // This function reades track seeds.
73042f01 754 //-----------------------------------------------------------------
755 TDirectory *savedir=gDirectory;
756
b9de75e1 757 TFile *in=(TFile*)inp;
758 if (!in->IsOpen()) {
c630aafd 759 Error("ReadSeeds","Input file has not been open !");
be5b9287 760 return 1;
73042f01 761 }
762
b9de75e1 763 in->cd();
764 TTree *seedTree=(TTree*)in->Get("Seeds");
765 if (!seedTree) {
c630aafd 766 Error("ReadSeeds","Can't get a tree with track seeds !");
b9de75e1 767 return 2;
768 }
769 AliTPCtrack *seed=new AliTPCtrack;
770 seedTree->SetBranchAddress("tracks",&seed);
771
772 if (fSeeds==0) fSeeds=new TObjArray(15000);
73042f01 773
b9de75e1 774 Int_t n=(Int_t)seedTree->GetEntries();
775 for (Int_t i=0; i<n; i++) {
776 seedTree->GetEvent(i);
777 fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
778 }
779
780 delete seed;
f38c8ae5 781
782 delete seedTree; //Thanks to Mariana Bondila
783
b9de75e1 784 savedir->cd();
73042f01 785
b9de75e1 786 return 0;
787}
73042f01 788
cab3f0a3 789//_____________________________________________________________________________
c630aafd 790Int_t AliTPCtracker::Clusters2Tracks(AliESD *event) {
b9de75e1 791 //-----------------------------------------------------------------
792 // This is a track finder.
c630aafd 793 // The clusters must be already loaded !
b9de75e1 794 //-----------------------------------------------------------------
b9de75e1 795
73042f01 796 //find track seeds
b9de75e1 797 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
73042f01 798 Int_t nrows=nlow+nup;
b9de75e1 799 if (fSeeds==0) {
800 Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
c630aafd 801 fSectors=fOuterSec; fN=fkNOS;
802 fSeeds=new TObjArray(15000);
b9de75e1 803 MakeSeeds(nup-1, nup-1-gap);
804 MakeSeeds(nup-1-shift, nup-1-shift-gap);
c630aafd 805 }
b9de75e1 806 fSeeds->Sort();
73042f01 807
b9de75e1 808 Int_t nseed=fSeeds->GetEntriesFast();
c630aafd 809 for (Int_t i=0; i<nseed; i++) {
810 //tracking in the outer sectors
811 fSectors=fOuterSec; fN=fkNOS;
812
b9de75e1 813 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
c630aafd 814 if (!FollowProlongation(t)) {
815 delete fSeeds->RemoveAt(i);
73042f01 816 continue;
817 }
c630aafd 818
819 //tracking in the inner sectors
820 fSectors=fInnerSec; fN=fkNIS;
821
822 Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
823 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
824 if (alpha < 0. ) alpha += 2.*TMath::Pi();
825 Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
826
827 alpha=ns*fInnerSec->GetAlpha()+fInnerSec->GetAlphaShift()-t.GetAlpha();
828
829 if (t.Rotate(alpha)) {
830 if (FollowProlongation(t)) {
831 if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
832 t.CookdEdx();
833 CookLabel(pt,0.1); //For comparison only
834 pt->PropagateTo(fParam->GetInnerRadiusLow());
835 AliESDtrack iotrack;
836 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
837
838 event->AddTrack(&iotrack);
839
840 UseClusters(&t);
841 }
842 }
843 }
88cb7938 844 delete fSeeds->RemoveAt(i);
c630aafd 845 }
846
847 Info("Clusters2Tracks","Number of found tracks : %d",
848 event->GetNumberOfTracks());
849
850 fSeeds->Clear(); delete fSeeds; fSeeds=0;
851
852 return 0;
853}
854
855//_____________________________________________________________________________
856Int_t AliTPCtracker::Clusters2Tracks(TTree *cTree, TTree *tTree) {
857 //-----------------------------------------------------------------
858 // This is a track finder.
859 //-----------------------------------------------------------------
860 if (LoadClusters(cTree)) {
861 Error("Clusters2Tracks","Problem with loading the clusters...");
862 return 1;
863 }
864
865 AliTPCtrack *iotrack=0;
866 TBranch *branch=tTree->GetBranch("tracks");
867 if (!branch) tTree->Branch("tracks","AliTPCtrack",&iotrack,32000,3);
868 else branch->SetAddress(&iotrack);
869
870 //find track seeds
871 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
872 Int_t nrows=nlow+nup;
873 if (fSeeds==0) {
874 Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
875 fSectors=fOuterSec; fN=fkNOS;
876 fSeeds=new TObjArray(15000);
877 MakeSeeds(nup-1, nup-1-gap);
878 MakeSeeds(nup-1-shift, nup-1-shift-gap);
879 }
880 fSeeds->Sort();
73042f01 881
88cb7938 882 Int_t found=0;
c630aafd 883 Int_t nseed=fSeeds->GetEntriesFast();
884
885 for (Int_t i=0; i<nseed; i++) {
886 //tracking in the outer sectors
887 fSectors=fOuterSec; fN=fkNOS;
888
88cb7938 889 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
c630aafd 890 if (!FollowProlongation(t)) {
891 delete fSeeds->RemoveAt(i);
892 continue;
893 }
894
895 //tracking in the inner sectors
896 fSectors=fInnerSec; fN=fkNIS;
73042f01 897
b9de75e1 898 Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
73042f01 899 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
900 if (alpha < 0. ) alpha += 2.*TMath::Pi();
b9de75e1 901 Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
73042f01 902
c630aafd 903 alpha=ns*fInnerSec->GetAlpha()+fInnerSec->GetAlphaShift()-t.GetAlpha();
73042f01 904
905 if (t.Rotate(alpha)) {
c630aafd 906 if (FollowProlongation(t)) {
907 if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
908 t.CookdEdx();
909 CookLabel(pt,0.1); //For comparison only
910 pt->PropagateTo(fParam->GetInnerRadiusLow());
911 iotrack=pt;
912 tTree->Fill();
913 UseClusters(&t);
914 found++;
915 }
916 }
73042f01 917 }
b9de75e1 918 delete fSeeds->RemoveAt(i);
c630aafd 919 }
920
921 Info("Clusters2Tracks","Number of found tracks : %d",found);
922
923 UnloadClusters();
924
925 fSeeds->Clear(); delete fSeeds; fSeeds=0;
926
927 return 0;
928}
929
930//_____________________________________________________________________________
931Int_t AliTPCtracker::RefitInward(AliESD* event) {
932 //
933 // The function propagates tracks throught TPC inward
934 // using already associated clusters.
935 // The clusters must be already loaded !
936 //
937
938 Int_t nTracks = event->GetNumberOfTracks();
939 Int_t nRefited = 0;
940
941 for (Int_t i = 0; i < nTracks; i++) {
942 AliESDtrack* track = event->GetTrack(i);
943 ULong_t status = track->GetStatus();
944
945 if ( (status & AliESDtrack::kTPCout ) == 0 ) continue;
946 if ( (status & AliESDtrack::kTPCrefit) != 0 ) continue;
947
948 AliTPCtrack* tpcTrack = new AliTPCtrack(*track);
949 AliTPCseed* seed = new AliTPCseed(*tpcTrack, tpcTrack->GetAlpha());
08dba16b 950
c630aafd 951 fSectors = fOuterSec;
952
953 Int_t res = FollowRefitInward(seed, tpcTrack);
954 UseClusters(seed);
955 Int_t nc = seed->GetNumberOfClusters();
956
957 fSectors = fInnerSec;
958
959 res = FollowRefitInward(seed, tpcTrack);
960 UseClusters(seed, nc);
961
962 if (res) {
963 seed->PropagateTo(fParam->GetInnerRadiusLow());
964 seed->SetLabel(tpcTrack->GetLabel());
965 seed->SetdEdx(tpcTrack->GetdEdx());
966 track->UpdateTrackParams(seed, AliESDtrack::kTPCrefit);
967 nRefited++;
968 }
969
970 delete seed;
971 delete tpcTrack;
972 }
973
974 Info("RefitInward","Number of refitted tracks : %d",nRefited);
975
976 return 0;
977}
978
979//_____________________________________________________________________________
980Int_t AliTPCtracker::RefitInward(TTree *in, TTree *out) {
981 //
982 // The function propagates tracks throught TPC inward
983 // using already associated clusters.
984 //
985 Error("RefitInward","This method is not converted to NewIO yet\n");
986 return 1;
987 /*
988 if (!inSeeds->IsOpen()) {
989 cout << "Input File with seeds not open !\n" << endl;
990 return 1;
991 }
992
993 if (!inTracks->IsOpen()) {
994 cout << "Input File not open !\n" << endl;
995 return 2;
996 }
9c9d2487 997
c630aafd 998 if (!outTracks->IsOpen()) {
999 cout << "Output File not open !\n" << endl;
1000 return 3;
1001 }
1002
1003 TDirectory *savedir = gDirectory;
1004 LoadClusters();
1005
1006 // Connect to input seeds tree
1007 inSeeds->cd();
1008 char tname[100];
1009 sprintf(tname, "seedTRDtoTPC_%d", GetEventNumber());
1010 TTree *seedsTree = (TTree*)inSeeds->Get(tname);
1011 TBranch *inSeedBranch = seedsTree->GetBranch("tracks");
1012 AliTPCtrack *inSeedTrack = 0;
1013 inSeedBranch->SetAddress(&inSeedTrack);
1014
1015 Int_t nSeeds = (Int_t)seedsTree->GetEntries();
1016
1017 // Connect to input tree
1018 inTracks->cd();
1019 sprintf(tname,"TreeT_TPC_%d",GetEventNumber());
1020// sprintf(tname,"seedsTPCtoTRD_%d",GetEventNumber());
1021 TTree *inputTree = (TTree*)inTracks->Get(tname);
1022 TBranch *inBranch = inputTree->GetBranch("tracks");
1023 AliTPCtrack *inTrack = 0;
1024 inBranch->SetAddress(&inTrack);
1025
1026 Int_t nTracks = (Int_t)inputTree->GetEntries();
1027
1028 // Create output tree
1029 outTracks->cd();
1030 AliTPCtrack *outTrack = new AliTPCtrack();
1031 sprintf(tname, "tracksTPC_%d", GetEventNumber());
1032 TTree *outputTree = new TTree(tname,"Refited TPC tracks");
1033 outputTree->Branch("tracks", "AliTPCtrack", &outTrack, 32000, 3);
1034
1035 Int_t nRefited = 0;
1036
1037 // create table of track labels
1038 Int_t* inLab = new Int_t[nTracks];
1039 for (Int_t i = 0; i < nTracks; i++) inLab[i] = -1;
1040
1041 // [SR, 01.04.2003] - Barrel tracks
1042 if (IsStoringBarrel()) SetBarrelTree("refit");
1043
1044 for(Int_t t=0; t < nSeeds; t++) {
1045
1046 seedsTree->GetEntry(t);
1047 // find TPC track for seed
1048 Int_t lab = TMath::Abs(inSeedTrack->GetLabel());
1049 for(Int_t i=0; i < nTracks; i++) {
1050 if (inLab[i] == lab) {
1051 inputTree->GetEntry(i);
1052 break;
1053 } else if (inLab[i] < 0) {
1054 inputTree->GetEntry(i);
1055 inLab[i] = TMath::Abs(inTrack->GetLabel());
1056 if (inLab[i] == lab) break;
1057 }
1058 }
1059 if (TMath::Abs(inTrack->GetLabel()) != lab) continue;
1060 AliTPCseed *seed = new AliTPCseed(*inSeedTrack, inTrack->GetAlpha());
1061 if (IsStoringBarrel()) StoreBarrelTrack(seed, 4, 2);
1062
1063 fSectors = fOuterSec;
1064
1065 Int_t res = FollowRefitInward(seed, inTrack);
1066 UseClusters(seed);
1067 Int_t nc = seed->GetNumberOfClusters();
1068
1069 if (IsStoringBarrel()) StoreBarrelTrack(seed, 3, 2);
1070 if (IsStoringBarrel()) StoreBarrelTrack(seed, 2, 2);
1071
1072 fSectors = fInnerSec;
1073
1074 res = FollowRefitInward(seed, inTrack);
1075 UseClusters(seed, nc);
1076
1077 if (IsStoringBarrel()) StoreBarrelTrack(seed, 1, 2);
1078
1079 if (res) {
1080 seed->PropagateTo(fParam->GetInnerRadiusLow());
1081 outTrack = (AliTPCtrack*)seed;
1082 outTrack->SetLabel(inTrack->GetLabel());
1083 outTrack->SetdEdx(inTrack->GetdEdx());
1084 outputTree->Fill();
1085 nRefited++;
1086 }
1087
1088 if (IsStoringBarrel()) fBarrelTree->Fill();
1089 delete seed;
1090 }
1091
1092 cout << "Refitted " << nRefited << " tracks." << endl;
1093
1094 outTracks->cd();
1095 outputTree->Write();
1096
1097 delete[] inLab;
1098
1099 // [SR, 01.04.2003]
1100 if (IsStoringBarrel()) {
1101 fBarrelFile->cd();
1102 fBarrelTree->Write();
1103 fBarrelFile->Flush();
1104
1105 delete fBarrelArray;
1106 delete fBarrelTree;
1107 }
1108
1109 if (inputTree) delete inputTree;
1110 if (outputTree) delete outputTree;
1111
1112 savedir->cd();
1113 UnloadClusters();
1114
1115 return 0;
1116 */
1117}
1118
1119Int_t AliTPCtracker::PropagateBack(AliESD *event) {
1120 //-----------------------------------------------------------------
1121 // This function propagates tracks back through the TPC.
1122 // The clusters must be already loaded !
1123 //-----------------------------------------------------------------
1124 Int_t nentr=event->GetNumberOfTracks();
1125 Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
1126
1127 Int_t ntrk=0;
1128 for (Int_t i=0; i<nentr; i++) {
1129 AliESDtrack *esd=event->GetTrack(i);
1130 ULong_t status=esd->GetStatus();
1131
1132 if ( (status & AliESDtrack::kTPCin ) == 0 ) continue;
1133 if ( (status & AliESDtrack::kTPCout) != 0 ) continue;
1134
1135 const AliTPCtrack t(*esd);
1136 AliTPCseed s(t,t.GetAlpha());
1137
1138 if ( (status & AliESDtrack::kITSout) == 0 ) s.ResetCovariance();
1139
1140 s.ResetNWrong();
1141 s.ResetNRotation();
1142
1143 Int_t nc=t.GetNumberOfClusters();
1144 s.SetNumber(nc); //set number of the cluster to start with
1145
1146 //inner sectors
1147 fSectors=fInnerSec; fN=fkNIS;
1148
1149 Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
1150 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1151 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1152 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
1153 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
1154 alpha-=s.GetAlpha();
1155
1156 if (!s.Rotate(alpha)) continue;
1157 if (!FollowBackProlongation(s,t)) continue;
1158
1159 UseClusters(&s);
1160
1161 //outer sectors
1162 fSectors=fOuterSec; fN=fkNOS;
1163
1164 nc=s.GetNumberOfClusters();
1165
1166 alpha=s.GetAlpha() - fSectors->GetAlphaShift();
1167 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1168 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1169 ns=Int_t(alpha/fSectors->GetAlpha())%fN;
1170
1171 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
1172 alpha-=s.GetAlpha();
1173
1174 if (!s.Rotate(alpha)) continue;
1175 if (!FollowBackProlongation(s,t)) continue;
1176 {
1177 Int_t nrows=fOuterSec->GetNRows()+fInnerSec->GetNRows();
1178 if (s.GetNumberOfClusters() < Int_t(0.4*nrows)) continue;
1179 }
1180 s.PropagateTo(fParam->GetOuterRadiusUp());
1181 s.CookdEdx();
1182 CookLabel(&s,0.1); //For comparison only
1183 UseClusters(&s,nc);
1184 esd->UpdateTrackParams(&s,AliESDtrack::kTPCout);
1185 ntrk++;
1186 }
1187 Info("PropagateBack","Number of back propagated tracks: %d",ntrk);
08dba16b 1188
1189 return 0;
1190}
b9de75e1 1191
1192//_____________________________________________________________________________
c630aafd 1193Int_t AliTPCtracker::PropagateBack(TTree *in, TTree *out) {
cab3f0a3 1194 //-----------------------------------------------------------------
1195 // This function propagates tracks back through the TPC.
cab3f0a3 1196 //-----------------------------------------------------------------
c630aafd 1197 Error("PropagateBack","This method is not converted to NewIO yet\n");
88cb7938 1198 return 1;
c630aafd 1199 /*
b9de75e1 1200 fSeeds=new TObjArray(15000);
1201
c630aafd 1202 TFile *in=(TFile*)inp;
1203 TFile *in2=(TFile*)inp2;
1204 TDirectory *savedir=gDirectory;
b9de75e1 1205
c630aafd 1206 if (!in->IsOpen()) {
1207 cerr<<"AliTPCtracker::PropagateBack(): ";
1208 cerr<<"file with TPC (or back propagated ITS) tracks is not open !\n";
e24ea474 1209 return 1;
c630aafd 1210 }
73042f01 1211
c630aafd 1212 if (!out->IsOpen()) {
1213 cerr<<"AliTPCtracker::PropagateBack(): ";
1214 cerr<<"file for back propagated TPC tracks is not open !\n";
1215 return 2;
b9de75e1 1216 }
1217
c630aafd 1218 LoadClusters();
bcc04e2a 1219
c630aafd 1220 in->cd();
1221 char tName[100];
1222 sprintf(tName,"TreeT_ITSb_%d",GetEventNumber());
1223 TTree *bckTree=(TTree*)in->Get(tName);
1224 if (!bckTree && inp2) bckTree=(TTree*)in2->Get(tName);
b9de75e1 1225 if (!bckTree) {
1226 cerr<<"AliTPCtracker::PropagateBack() ";
1227 cerr<<"can't get a tree with back propagated ITS tracks !\n";
c630aafd 1228 //return 3;
b9de75e1 1229 }
c630aafd 1230 AliTPCtrack *bckTrack=new AliTPCtrack;
1231 if (bckTree) bckTree->SetBranchAddress("tracks",&bckTrack);
b9de75e1 1232
c630aafd 1233 sprintf(tName,"TreeT_TPC_%d",GetEventNumber());
1234 TTree *tpcTree=(TTree*)in->Get(tName);
b9de75e1 1235 if (!tpcTree) {
1236 cerr<<"AliTPCtracker::PropagateBack() ";
1237 cerr<<"can't get a tree with TPC tracks !\n";
1238 return 4;
1239 }
1240 AliTPCtrack *tpcTrack=new AliTPCtrack;
1241 tpcTree->SetBranchAddress("tracks",&tpcTrack);
1242
c630aafd 1243 // [SR, 01.04.2003] - Barrel tracks
1244 if (IsStoringBarrel()) SetBarrelTree("back");
1245
1246 // Prepare an array of tracks to be back propagated
b9de75e1 1247 Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
1248 Int_t nrows=nlow+nup;
1249
c630aafd 1250 //
1251 // Match ITS tracks with old TPC tracks
1252 // the tracks do not have to be sorted [SR, GSI, 18.02.2003]
1253 //
1254 // the algorithm is linear and uses LUT for sorting
1255 // cost of the algorithm: 2 * number of tracks
1256 //
1257
b9de75e1 1258 TObjArray tracks(15000);
c630aafd 1259 Int_t i=0;
1260 Int_t tpcN= (Int_t)tpcTree->GetEntries();
1261 Int_t bckN= (bckTree)? (Int_t)bckTree->GetEntries() : 0;
1262
1263 // look up table
1264 const Int_t nLab = 200000; // limit on number of primaries (arbitrary)
1265 Int_t lut[nLab];
1266 for(Int_t i=0; i<nLab; i++) lut[i] = -1;
1267
1268 if (bckTree) {
1269 for(Int_t i=0; i<bckN; i++) {
1270 bckTree->GetEvent(i);
1271 Int_t lab = TMath::Abs(bckTrack->GetLabel());
1272 if (lab < nLab) lut[lab] = i;
1273 else {
1274 cerr << "AliTPCtracker: The size of the LUT is too small\n";
1275 }
1276 }
1277 }
1278
1279 for (Int_t i=0; i<tpcN; i++) {
1280 tpcTree->GetEvent(i);
1281 Double_t alpha=tpcTrack->GetAlpha();
1282
1283 if (!bckTree) {
1284
1285 // No ITS - use TPC track only
1286
1287 tpcTrack->ResetCovariance();
1288 AliTPCseed * seed = new AliTPCseed(*tpcTrack, tpcTrack->GetAlpha());
1289
1290 fSeeds->AddLast(seed);
1291 tracks.AddLast(new AliTPCtrack(*tpcTrack));
1292
1293 } else {
1294
1295 // with ITS
1296 // discard not prolongated tracks (to be discussed)
1297
1298 Int_t lab = TMath::Abs(tpcTrack->GetLabel());
1299 if (lab > nLab || lut[lab] < 0) continue;
1300
1301 bckTree->GetEvent(lut[lab]);
1302 bckTrack->Rotate(alpha-bckTrack->GetAlpha());
1303
1304 fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
1305 tracks.AddLast(new AliTPCtrack(*tpcTrack));
1306 }
b9de75e1 1307 }
08dba16b 1308
c630aafd 1309 // end of matching
1310
88cb7938 1311 out->cd();
c630aafd 1312
1313 // tree name seedsTPCtoTRD as expected by TRD and as
1314 // discussed and decided in Strasbourg (May 2002)
1315 // [SR, GSI, 18.02.2003]
1316
1317 sprintf(tName,"seedsTPCtoTRD_%d",GetEventNumber());
1318 TTree backTree(tName,"Tree with back propagated TPC tracks");
1319
b9de75e1 1320 AliTPCtrack *otrack=0;
1321 backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
1322
c630aafd 1323 //
1324 Int_t nRefPlane;
1325 Int_t found=0;
b9de75e1 1326 Int_t nseed=fSeeds->GetEntriesFast();
c630aafd 1327
1328 // loop changed [SR, 01.04.2003]
1329
b9de75e1 1330 for (i=0; i<nseed; i++) {
c630aafd 1331
1332 nRefPlane = 1;
1333 if (IsStoringBarrel()) fBarrelArray->Clear();
1334
b9de75e1 1335 AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
1336 const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
73042f01 1337
c630aafd 1338 ps->ResetNWrong();
1339 ps->ResetNRotation();
1340
1341 if (IsStoringBarrel()) StoreBarrelTrack(ps, 1, 1);
1342
1343 // Load outer sectors
1344 fSectors=fInnerSec;
1345 fN=fkNIS;
1346
b9de75e1 1347 Int_t nc=t.GetNumberOfClusters();
c630aafd 1348 //s.SetLabel(nc-1); //set number of the cluster to start with
1349 s.SetNumber(nc);
73042f01 1350
c630aafd 1351 if (FollowBackProlongation(s,t)) UseClusters(&s);
1352 else {
1353 fSeeds->RemoveAt(i);
b9de75e1 1354 continue;
1355 }
b9de75e1 1356
c630aafd 1357 if (IsStoringBarrel()) StoreBarrelTrack(ps, 2, 1);
1358 if (IsStoringBarrel()) StoreBarrelTrack(ps, 3, 1);
1359
1360 // Load outer sectors
1361 fSectors=fOuterSec;
1362 fN=fkNOS;
1363
1364 nc=s.GetNumberOfClusters();
b9de75e1 1365
1366 Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
1367 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1368 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1369 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
1370
1371 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
1372 alpha-=s.GetAlpha();
1373
1374 if (s.Rotate(alpha)) {
c630aafd 1375 if (FollowBackProlongation(s,t)) {
1376 if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
1377 s.CookdEdx();
1378 s.SetLabel(t.GetLabel());
1379 UseClusters(&s,nc);
1380
1381 // Propagate to outer reference plane for comparison reasons
1382 // reason for keeping fParam object [SR, GSI, 18.02.2003]
1383
1384 ps->PropagateTo(fParam->GetOuterRadiusUp());
1385 otrack=ps;
1386 backTree.Fill();
1387
1388 if (IsStoringBarrel()) StoreBarrelTrack(ps, 4, 1);
1389 cerr<<found++<<'\r';
1390 }
e3f7105e 1391 }
b9de75e1 1392 }
c630aafd 1393
1394 if (IsStoringBarrel()) fBarrelTree->Fill();
b9de75e1 1395 delete fSeeds->RemoveAt(i);
1396 }
b9de75e1 1397
c630aafd 1398
1399 out->cd();
b9de75e1 1400 backTree.Write();
9c9d2487 1401
c630aafd 1402 // [SR, 01.04.2003]
1403 if (IsStoringBarrel()) {
1404 fBarrelFile->cd();
1405 fBarrelTree->Write();
1406 fBarrelFile->Flush();
1407
1408 delete fBarrelArray;
1409 delete fBarrelTree;
1410 }
1411
1412 savedir->cd();
b9de75e1 1413 cerr<<"Number of seeds: "<<nseed<<endl;
1414 cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
1415 cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
1416
1417 delete bckTrack;
1418 delete tpcTrack;
be5b9287 1419
c630aafd 1420 if (bckTree) delete bckTree; //Thanks to Mariana Bondila
f38c8ae5 1421 delete tpcTree; //Thanks to Mariana Bondila
1422
c630aafd 1423 UnloadClusters();
1424
be5b9287 1425 return 0;
c630aafd 1426 */
73042f01 1427}
1428
b9de75e1 1429//_________________________________________________________________________
1430AliCluster *AliTPCtracker::GetCluster(Int_t index) const {
1431 //--------------------------------------------------------------------
1432 // Return pointer to a given cluster
1433 //--------------------------------------------------------------------
1434 Int_t sec=(index&0xff000000)>>24;
1435 Int_t row=(index&0x00ff0000)>>16;
1436 Int_t ncl=(index&0x0000ffff)>>00;
1437
c630aafd 1438 const AliTPCcluster *cl=0;
1439
1440 if (sec<fkNIS) {
1441 cl=fInnerSec[sec][row].GetUnsortedCluster(ncl);
1442 } else {
1443 sec-=fkNIS;
1444 cl=fOuterSec[sec][row].GetUnsortedCluster(ncl);
1445 }
1446
1447 return (AliCluster*)cl;
b9de75e1 1448}
1449
1450//__________________________________________________________________________
1451void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
1452 //--------------------------------------------------------------------
1453 //This function "cooks" a track label. If label<0, this track is fake.
1454 //--------------------------------------------------------------------
1455 Int_t noc=t->GetNumberOfClusters();
1456 Int_t *lb=new Int_t[noc];
1457 Int_t *mx=new Int_t[noc];
1458 AliCluster **clusters=new AliCluster*[noc];
1459
1460 Int_t i;
1461 for (i=0; i<noc; i++) {
1462 lb[i]=mx[i]=0;
1463 Int_t index=t->GetClusterIndex(i);
1464 clusters[i]=GetCluster(index);
1465 }
1466
1467 Int_t lab=123456789;
1468 for (i=0; i<noc; i++) {
1469 AliCluster *c=clusters[i];
1470 lab=TMath::Abs(c->GetLabel(0));
1471 Int_t j;
1472 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
1473 lb[j]=lab;
1474 (mx[j])++;
1475 }
1476
1477 Int_t max=0;
1478 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
1479
1480 for (i=0; i<noc; i++) {
1481 AliCluster *c=clusters[i];
1482 if (TMath::Abs(c->GetLabel(1)) == lab ||
1483 TMath::Abs(c->GetLabel(2)) == lab ) max++;
1484 }
1485
1486 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
1487
1488 else {
1489 Int_t tail=Int_t(0.10*noc);
1490 max=0;
1491 for (i=1; i<=tail; i++) {
1492 AliCluster *c=clusters[noc-i];
1493 if (lab == TMath::Abs(c->GetLabel(0)) ||
1494 lab == TMath::Abs(c->GetLabel(1)) ||
1495 lab == TMath::Abs(c->GetLabel(2))) max++;
1496 }
1497 if (max < Int_t(0.5*tail)) lab=-lab;
1498 }
1499
1500 t->SetLabel(lab);
1501
1502 delete[] lb;
1503 delete[] mx;
1504 delete[] clusters;
1505}
1506
1507//_________________________________________________________________________
1508void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
1509 //-----------------------------------------------------------------------
1510 // Setup inner sector
1511 //-----------------------------------------------------------------------
1512 if (f==0) {
1513 fAlpha=par->GetInnerAngle();
1514 fAlphaShift=par->GetInnerAngleShift();
1515 fPadPitchWidth=par->GetInnerPadPitchWidth();
024a7fe9 1516 f1PadPitchLength=par->GetInnerPadPitchLength();
1517 f2PadPitchLength=f1PadPitchLength;
b9de75e1 1518
1519 fN=par->GetNRowLow();
1520 fRow=new AliTPCRow[fN];
1521 for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
1522 } else {
1523 fAlpha=par->GetOuterAngle();
1524 fAlphaShift=par->GetOuterAngleShift();
1525 fPadPitchWidth=par->GetOuterPadPitchWidth();
f03e3423 1526 f1PadPitchLength=par->GetOuter1PadPitchLength();
1527 f2PadPitchLength=par->GetOuter2PadPitchLength();
b9de75e1 1528
1529 fN=par->GetNRowUp();
1530 fRow=new AliTPCRow[fN];
f03e3423 1531 for (Int_t i=0; i<fN; i++){
1532 fRow[i].SetX(par->GetPadRowRadiiUp(i));
1533 }
b9de75e1 1534 }
1535}
1536
73042f01 1537//_________________________________________________________________________
c630aafd 1538void AliTPCtracker::
1539AliTPCRow::InsertCluster(const AliTPCcluster* c, Int_t sec, Int_t row) {
73042f01 1540 //-----------------------------------------------------------------------
1541 // Insert a cluster into this pad row in accordence with its y-coordinate
1542 //-----------------------------------------------------------------------
b9de75e1 1543 if (fN==kMaxClusterPerRow) {
c630aafd 1544 ::Error("InsertCluster","Too many clusters !");
1545 return;
1546 }
1547
1548 Int_t index=(((sec<<8)+row)<<16)+fN;
1549
1550 if (fN==0) {
1551 fSize=kMaxClusterPerRow/8;
1552 fClusterArray=new AliTPCcluster[fSize];
1553 fIndex[0]=index;
1554 fClusterArray[0]=*c; fClusters[fN++]=fClusterArray;
1555 return;
73042f01 1556 }
c630aafd 1557
1558 if (fN==fSize) {
1559 Int_t size=fSize*2;
1560 AliTPCcluster *buff=new AliTPCcluster[size];
1561 memcpy(buff,fClusterArray,fSize*sizeof(AliTPCcluster));
1562 for (Int_t i=0; i<fN; i++)
1563 fClusters[i]=buff+(fClusters[i]-fClusterArray);
1564 delete[] fClusterArray;
1565 fClusterArray=buff;
1566 fSize=size;
1567 }
1568
73042f01 1569 Int_t i=Find(c->GetY());
1570 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
1571 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
c630aafd 1572 fIndex[i]=index;
1573 fClusters[i]=fClusterArray+fN; fClusterArray[fN++]=*c;
73042f01 1574}
1575
1576//___________________________________________________________________
1577Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
1578 //-----------------------------------------------------------------------
1579 // Return the index of the nearest cluster
1580 //-----------------------------------------------------------------------
11c7b548 1581 if(fN<=0) return 0;
73042f01 1582 if (y <= fClusters[0]->GetY()) return 0;
1583 if (y > fClusters[fN-1]->GetY()) return fN;
1584 Int_t b=0, e=fN-1, m=(b+e)/2;
1585 for (; b<e; m=(b+e)/2) {
1586 if (y > fClusters[m]->GetY()) b=m+1;
1587 else e=m;
1588 }
1589 return m;
1590}
1591
1592//_____________________________________________________________________________
1593void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
1594 //-----------------------------------------------------------------
1595 // This funtion calculates dE/dX within the "low" and "up" cuts.
1596 //-----------------------------------------------------------------
1597 Int_t i;
be9c5115 1598 Int_t nc=GetNumberOfClusters();
c630aafd 1599
1600 Int_t swap;//stupid sorting
1601 do {
1602 swap=0;
1603 for (i=0; i<nc-1; i++) {
1604 if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
1605 Float_t tmp=fdEdxSample[i];
1606 fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
1607 swap++;
1608 }
1609 } while (swap);
73042f01 1610
be9c5115 1611 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
73042f01 1612 Float_t dedx=0;
c630aafd 1613 for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
73042f01 1614 dedx /= (nu-nl+1);
1615 SetdEdx(dedx);
7f6ddf58 1616
1617 //Very rough PID
1618 Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
1619
024a7fe9 1620 Double_t log1=TMath::Log(p+0.45), log2=TMath::Log(p+0.12);
7f6ddf58 1621 if (p<0.6) {
024a7fe9 1622 if (dedx < 34 + 30/(p+0.45)/(p+0.45) + 24*log1) {SetMass(0.13957); return;}
1623 if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.49368); return;}
7f6ddf58 1624 SetMass(0.93827); return;
1625 }
1626
1627 if (p<1.2) {
024a7fe9 1628 if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.13957); return;}
7f6ddf58 1629 SetMass(0.93827); return;
1630 }
1631
1632 SetMass(0.13957); return;
1633
73042f01 1634}
1635
73042f01 1636
1637