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