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