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