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