]>
Commit | Line | Data |
---|---|---|
46d29e70 | 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$ | |
bdbd0f7a | 18 | Revision 1.10 2001/05/07 08:08:05 cblume |
19 | Update of TRD code | |
20 | ||
a2b90f83 | 21 | Revision 1.9 2001/02/14 18:22:26 cblume |
22 | Change in the geometry of the padplane | |
23 | ||
71d9fa7b | 24 | Revision 1.8 2000/12/20 13:00:44 cblume |
25 | Modifications for the HP-compiler | |
26 | ||
d1b06c24 | 27 | Revision 1.7 2000/12/08 16:07:02 cblume |
28 | Update of the tracking by Sergei | |
29 | ||
bbf92647 | 30 | Revision 1.6 2000/11/30 17:38:08 cblume |
31 | Changes to get in line with new STEER and EVGEN | |
32 | ||
1819f4bb | 33 | Revision 1.5 2000/11/14 14:40:27 cblume |
34 | Correction for the Sun compiler (kTRUE and kFALSE) | |
35 | ||
57527628 | 36 | Revision 1.4 2000/11/10 14:57:52 cblume |
37 | Changes in the geometry constants for the DEC compiler | |
38 | ||
dd56b762 | 39 | Revision 1.3 2000/10/15 23:40:01 cblume |
40 | Remove AliTRDconst | |
41 | ||
0e9c2ad5 | 42 | Revision 1.2 2000/10/06 16:49:46 cblume |
43 | Made Getters const | |
44 | ||
46d29e70 | 45 | Revision 1.1.2.2 2000/10/04 16:34:58 cblume |
46 | Replace include files by forward declarations | |
47 | ||
48 | Revision 1.1.2.1 2000/09/22 14:47:52 cblume | |
49 | Add the tracking code | |
50 | ||
51 | */ | |
bbf92647 | 52 | |
1819f4bb | 53 | #include <iostream.h> |
46d29e70 | 54 | |
55 | #include <TFile.h> | |
56 | #include <TROOT.h> | |
57 | #include <TBranch.h> | |
58 | #include <TTree.h> | |
59 | ||
60 | #include "AliRun.h" | |
46d29e70 | 61 | #include "AliTRD.h" |
46d29e70 | 62 | #include "AliTRDgeometry.h" |
63 | #include "AliTRDrecPoint.h" | |
64 | #include "AliTRDcluster.h" | |
65 | #include "AliTRDtrack.h" | |
66 | #include "AliTRDtrackingSector.h" | |
67 | #include "AliTRDtimeBin.h" | |
68 | ||
69 | #include "AliTRDtracker.h" | |
70 | ||
71 | ClassImp(AliTRDtracker) | |
72 | ||
bbf92647 | 73 | const Int_t AliTRDtracker::fSeedGap = 35; |
74 | const Int_t AliTRDtracker::fSeedStep = 5; | |
75 | ||
76 | ||
77 | const Float_t AliTRDtracker::fMinClustersInTrack = 0.5; | |
78 | const Float_t AliTRDtracker::fMinClustersInSeed = 0.5; | |
79 | const Float_t AliTRDtracker::fSeedDepth = 0.5; | |
80 | const Float_t AliTRDtracker::fSkipDepth = 0.2; | |
81 | const Float_t AliTRDtracker::fMaxSeedDeltaZ = 30.; | |
82 | const Float_t AliTRDtracker::fMaxSeedC = 0.01; | |
83 | const Float_t AliTRDtracker::fMaxSeedTan = 1.2; | |
84 | const Float_t AliTRDtracker::fMaxSeedVertexZ = 200.; | |
85 | const Float_t AliTRDtracker::fLabelFraction = 0.5; | |
86 | const Float_t AliTRDtracker::fWideRoad = 30.; | |
87 | ||
88 | const Double_t AliTRDtracker::fMaxChi2 = 12.; | |
89 | const Double_t AliTRDtracker::fSeedErrorSY = 0.1; | |
90 | const Double_t AliTRDtracker::fSeedErrorSY3 = 2.5; | |
91 | const Double_t AliTRDtracker::fSeedErrorSZ = 0.1; | |
92 | ||
46d29e70 | 93 | //____________________________________________________________________ |
94 | AliTRDtracker::AliTRDtracker() | |
95 | { | |
96 | // | |
97 | // Default constructor | |
98 | // | |
99 | ||
46d29e70 | 100 | fEvent = 0; |
46d29e70 | 101 | fGeom = NULL; |
bbf92647 | 102 | |
46d29e70 | 103 | fNclusters = 0; |
104 | fClusters = NULL; | |
105 | fNseeds = 0; | |
106 | fSeeds = NULL; | |
107 | fNtracks = 0; | |
108 | fTracks = NULL; | |
109 | ||
110 | } | |
111 | ||
112 | //____________________________________________________________________ | |
113 | AliTRDtracker::AliTRDtracker(const Text_t* name, const Text_t* title) | |
114 | :TNamed(name, title) | |
115 | { | |
46d29e70 | 116 | fEvent = 0; |
46d29e70 | 117 | fGeom = NULL; |
bbf92647 | 118 | |
46d29e70 | 119 | fNclusters = 0; |
120 | fClusters = new TObjArray(2000); | |
121 | fNseeds = 0; | |
122 | fSeeds = new TObjArray(20000); | |
123 | fNtracks = 0; | |
124 | fTracks = new TObjArray(10000); | |
125 | ||
126 | } | |
127 | ||
46d29e70 | 128 | //___________________________________________________________________ |
129 | AliTRDtracker::~AliTRDtracker() | |
130 | { | |
46d29e70 | 131 | delete fClusters; |
132 | delete fTracks; | |
133 | delete fSeeds; | |
134 | delete fGeom; | |
135 | } | |
136 | ||
bbf92647 | 137 | //___________________________________________________________________ |
138 | void AliTRDtracker::Clusters2Tracks() | |
139 | { | |
140 | Int_t inner, outer; | |
141 | Int_t fTotalNofTimeBins = fGeom->GetTimeMax() * AliTRDgeometry::Nplan(); | |
142 | Int_t nSteps = (Int_t) (fTotalNofTimeBins * fSeedDepth)/fSeedStep; | |
143 | ||
144 | for(Int_t i=0; i<nSteps; i++) { | |
145 | printf("step %d out of %d \n", i+1, nSteps); | |
146 | outer=fTotalNofTimeBins-1-i*fSeedStep; inner=outer-fSeedGap; | |
147 | MakeSeeds(inner,outer); | |
148 | FindTracks(); | |
149 | } | |
150 | } | |
46d29e70 | 151 | |
152 | //_____________________________________________________________________ | |
bbf92647 | 153 | Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt) |
46d29e70 | 154 | { |
155 | // Parametrised "expected" error of the cluster reconstruction in Y | |
156 | ||
bbf92647 | 157 | Double_t s = 0.2; |
46d29e70 | 158 | return s; |
159 | } | |
160 | ||
161 | //_____________________________________________________________________ | |
bbf92647 | 162 | Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl) |
46d29e70 | 163 | { |
164 | // Parametrised "expected" error of the cluster reconstruction in Z | |
165 | ||
71d9fa7b | 166 | // Take an example pad size for the time being, check back |
167 | // again with Sergei | |
168 | Double_t s, pad = fGeom->GetRowPadSize(0,2,0); | |
46d29e70 | 169 | s = pad * pad /12.; |
170 | return s; | |
171 | } | |
172 | ||
173 | //_____________________________________________________________________ | |
bdbd0f7a | 174 | Double_t f1trd(Double_t x1,Double_t y1, |
46d29e70 | 175 | Double_t x2,Double_t y2, |
176 | Double_t x3,Double_t y3) | |
177 | { | |
178 | // Initial approximation of the track curvature | |
179 | // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch | |
180 | ||
181 | Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1); | |
182 | Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)- | |
183 | (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2)); | |
184 | Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)- | |
185 | (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1)); | |
186 | ||
187 | Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b); | |
188 | ||
189 | return -xr*yr/sqrt(xr*xr+yr*yr); | |
190 | } | |
191 | ||
192 | //_____________________________________________________________________ | |
bdbd0f7a | 193 | Double_t f2trd(Double_t x1,Double_t y1, |
46d29e70 | 194 | Double_t x2,Double_t y2, |
195 | Double_t x3,Double_t y3) | |
196 | { | |
197 | // Initial approximation of the track curvature times center of curvature | |
198 | // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch | |
199 | ||
200 | Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1); | |
201 | Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)- | |
202 | (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2)); | |
203 | Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)- | |
204 | (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1)); | |
205 | ||
206 | Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b); | |
207 | ||
208 | return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr); | |
209 | } | |
210 | ||
211 | //_____________________________________________________________________ | |
bdbd0f7a | 212 | Double_t f3trd(Double_t x1,Double_t y1, |
46d29e70 | 213 | Double_t x2,Double_t y2, |
214 | Double_t z1,Double_t z2) | |
215 | { | |
216 | // Initial approximation of the tangent of the track dip angle | |
217 | // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch | |
218 | ||
219 | return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); | |
220 | } | |
221 | ||
222 | ||
223 | //___________________________________________________________________ | |
224 | ||
bbf92647 | 225 | Int_t AliTRDtracker::FindProlongation(AliTRDtrack& t, AliTRDtrackingSector *sec, |
d1b06c24 | 226 | Int_t s, Int_t rf) |
46d29e70 | 227 | { |
228 | // Starting from current position on track=t this function tries | |
229 | // to extrapolate the track up to timeBin=rf and to confirm prolongation | |
230 | // if a close cluster is found. *sec is a pointer to allocated | |
231 | // array of sectors, in which the initial sector has index=s. | |
232 | ||
bbf92647 | 233 | const Int_t TIME_BINS_TO_SKIP=Int_t(fSkipDepth*sec->GetNtimeBins()); |
46d29e70 | 234 | Int_t try_again=TIME_BINS_TO_SKIP; |
235 | ||
0e9c2ad5 | 236 | Double_t alpha=AliTRDgeometry::GetAlpha(); |
46d29e70 | 237 | |
238 | Int_t ns=Int_t(2*TMath::Pi()/alpha+0.5); | |
239 | ||
240 | for (Int_t nr=sec->GetTimeBinNumber(t.GetX())-1; nr>=rf; nr--) { | |
241 | Double_t x=sec->GetX(nr), ymax=sec->GetMaxY(nr); | |
242 | if (!t.PropagateTo(x)) { | |
243 | cerr<<"Can't propagate to x = "<<x<<endl; | |
244 | return 0; | |
245 | } | |
246 | ||
247 | AliTRDcluster *cl=0; | |
248 | UInt_t index=0; | |
249 | ||
bbf92647 | 250 | Double_t max_chi2=fMaxChi2; |
251 | ||
46d29e70 | 252 | AliTRDtimeBin& time_bin=sec[s][nr]; |
bbf92647 | 253 | Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt()); |
254 | Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl()); | |
46d29e70 | 255 | Double_t road=5.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ(); |
256 | ||
bbf92647 | 257 | if (road>fWideRoad) { |
46d29e70 | 258 | if (t.GetNclusters() > 4) { |
259 | cerr<<t.GetNclusters()<<" FindProlongation: Road is too wide !\n"; | |
260 | } | |
261 | return 0; | |
262 | } | |
263 | ||
264 | if (time_bin) { | |
265 | for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) { | |
266 | AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]); | |
267 | if (c->GetY() > y+road) break; | |
268 | if (c->IsUsed() > 0) continue; | |
269 | ||
270 | if((c->GetZ()-z)*(c->GetZ()-z) > 25. + sz2) continue; | |
271 | ||
272 | Double_t chi2=t.GetPredictedChi2(c); | |
273 | ||
274 | if (chi2 > max_chi2) continue; | |
275 | max_chi2=chi2; | |
276 | cl=c; | |
277 | index=time_bin.GetIndex(i); | |
278 | } | |
279 | } | |
280 | if (cl) { | |
281 | ||
282 | // Float_t l=sec->GetPitch(); | |
283 | // t.SetSampledEdx(cl->fQ/l,Int_t(t)); | |
284 | ||
285 | t.Update(cl,max_chi2,index); | |
286 | ||
287 | try_again=TIME_BINS_TO_SKIP; | |
288 | } else { | |
289 | if (try_again==0) break; | |
290 | if (y > ymax) { | |
291 | cerr<<"y > ymax: "<<y<<" > "<<ymax<<endl; | |
292 | s = (s+1) % ns; | |
293 | if (!t.Rotate(alpha)) { | |
294 | cerr<<"Failed to rotate, alpha = "<<alpha<<endl; | |
295 | return 0; | |
296 | } | |
297 | } else if (y <-ymax) { | |
298 | cerr<<"y < -ymax: "<<y<<" < "<<-ymax<<endl; | |
299 | s = (s-1+ns) % ns; | |
300 | if (!t.Rotate(-alpha)) { | |
301 | cerr<<"Failed to rotate, alpha = "<<alpha<<endl; | |
302 | return 0; | |
303 | } | |
304 | } | |
305 | try_again--; | |
306 | } | |
307 | } | |
308 | ||
309 | return 1; | |
310 | } | |
311 | ||
312 | ||
313 | ||
314 | //_____________________________________________________________________________ | |
bbf92647 | 315 | void AliTRDtracker::GetEvent(const Char_t *hitfile, const Char_t *clusterfile) |
46d29e70 | 316 | { |
317 | // Opens a ROOT-file with TRD-clusters and reads in the cluster-tree | |
46d29e70 | 318 | |
bbf92647 | 319 | ReadClusters(fClusters, clusterfile); |
320 | ||
321 | // get geometry from the file with hits | |
322 | ||
323 | TFile *fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(hitfile); | |
46d29e70 | 324 | if (!fInputFile) { |
325 | printf("AliTRDtracker::Open -- "); | |
bbf92647 | 326 | printf("Open the ALIROOT-file %s.\n",hitfile); |
327 | fInputFile = new TFile(hitfile); | |
46d29e70 | 328 | } |
329 | else { | |
330 | printf("AliTRDtracker::Open -- "); | |
bbf92647 | 331 | printf("%s is already open.\n",hitfile); |
46d29e70 | 332 | } |
333 | ||
334 | // Get AliRun object from file or create it if not on file | |
46d29e70 | 335 | |
bbf92647 | 336 | gAlice = (AliRun*) fInputFile->Get("gAlice"); |
337 | if (gAlice) { | |
338 | printf("AliTRDtracker::GetEvent -- "); | |
339 | printf("AliRun object found on file.\n"); | |
340 | } | |
341 | else { | |
342 | printf("AliTRDtracker::GetEvent -- "); | |
343 | printf("Could not find AliRun object.\n"); | |
344 | } | |
46d29e70 | 345 | |
346 | // Import the Trees for the event nEvent in the file | |
347 | Int_t nparticles = gAlice->GetEvent(fEvent); | |
348 | cerr<<"nparticles = "<<nparticles<<endl; | |
349 | ||
350 | if (nparticles <= 0) { | |
351 | printf("AliTRDtracker::GetEvent -- "); | |
352 | printf("No entries in the trees for event %d.\n",fEvent); | |
353 | } | |
354 | ||
355 | AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD"); | |
356 | fGeom = TRD->GetGeometry(); | |
357 | ||
46d29e70 | 358 | } |
359 | ||
360 | ||
361 | //_____________________________________________________________________________ | |
362 | void AliTRDtracker::SetUpSectors(AliTRDtrackingSector *sec) | |
363 | { | |
364 | // Fills clusters into TRD tracking_sectors | |
365 | // Note that the numbering scheme for the TRD tracking_sectors | |
366 | // differs from that of TRD sectors | |
367 | ||
bbf92647 | 368 | for (Int_t i=0; i<AliTRDgeometry::Nsect(); i++) sec[i].SetUp(); |
46d29e70 | 369 | |
370 | // Sort clusters into AliTRDtimeBin's within AliTRDtrackSector's | |
371 | ||
372 | cerr<<"MakeSeeds: sorting clusters"<<endl; | |
373 | ||
374 | Int_t ncl=fClusters->GetEntriesFast(); | |
375 | UInt_t index; | |
376 | while (ncl--) { | |
377 | printf("\r %d left",ncl); | |
378 | AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl); | |
379 | Int_t detector=c->GetDetector(), local_time_bin=c->GetLocalTimeBin(); | |
380 | Int_t sector=fGeom->GetSector(detector); | |
381 | ||
382 | Int_t tracking_sector=sector; | |
dd56b762 | 383 | if(sector > 0) tracking_sector = AliTRDgeometry::kNsect - sector; |
46d29e70 | 384 | |
385 | Int_t tb=sec[sector].GetTimeBin(detector,local_time_bin); | |
386 | index=ncl; | |
387 | sec[tracking_sector][tb].InsertCluster(c,index); | |
388 | } | |
389 | printf("\r\n"); | |
390 | } | |
391 | ||
392 | ||
393 | //_____________________________________________________________________________ | |
394 | void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer) | |
395 | { | |
396 | // Creates track seeds using clusters in timeBins=i1,i2 | |
397 | ||
bbf92647 | 398 | Int_t i2 = inner, i1 = outer; |
46d29e70 | 399 | |
400 | if (!fClusters) return; | |
401 | ||
dd56b762 | 402 | AliTRDtrackingSector fTrSec[AliTRDgeometry::kNsect]; |
46d29e70 | 403 | SetUpSectors(fTrSec); |
404 | ||
405 | // find seeds | |
406 | ||
407 | Double_t x[5], c[15]; | |
dd56b762 | 408 | Int_t max_sec=AliTRDgeometry::kNsect; |
46d29e70 | 409 | |
0e9c2ad5 | 410 | Double_t alpha=AliTRDgeometry::GetAlpha(); |
411 | Double_t shift=AliTRDgeometry::GetAlpha()/2.; | |
46d29e70 | 412 | Double_t cs=cos(alpha), sn=sin(alpha); |
413 | ||
414 | Double_t x1 =fTrSec[0].GetX(i1); | |
415 | Double_t xx2=fTrSec[0].GetX(i2); | |
416 | ||
417 | for (Int_t ns=0; ns<max_sec; ns++) { | |
418 | ||
419 | Int_t nl=fTrSec[(ns-1+max_sec)%max_sec][i2]; | |
420 | Int_t nm=fTrSec[ns][i2]; | |
421 | Int_t nu=fTrSec[(ns+1)%max_sec][i2]; | |
422 | ||
423 | AliTRDtimeBin& r1=fTrSec[ns][i1]; | |
424 | ||
425 | for (Int_t is=0; is < r1; is++) { | |
426 | Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ(); | |
427 | ||
428 | for (Int_t js=0; js < nl+nm+nu; js++) { | |
429 | const AliTRDcluster *cl; | |
430 | Double_t x2, y2, z2; | |
431 | Double_t x3=0.,y3=0.; | |
432 | ||
433 | if (js<nl) { | |
434 | AliTRDtimeBin& r2=fTrSec[(ns-1+max_sec)%max_sec][i2]; | |
435 | cl=r2[js]; | |
436 | y2=cl->GetY(); z2=cl->GetZ(); | |
437 | ||
438 | x2= xx2*cs+y2*sn; | |
439 | y2=-xx2*sn+y2*cs; | |
440 | ||
441 | } else | |
442 | if (js<nl+nm) { | |
443 | AliTRDtimeBin& r2=fTrSec[ns][i2]; | |
444 | cl=r2[js-nl]; | |
445 | x2=xx2; y2=cl->GetY(); z2=cl->GetZ(); | |
446 | } else { | |
447 | AliTRDtimeBin& r2=fTrSec[(ns+1)%max_sec][i2]; | |
448 | cl=r2[js-nl-nm]; | |
449 | y2=cl->GetY(); z2=cl->GetZ(); | |
450 | ||
451 | x2=xx2*cs-y2*sn; | |
452 | y2=xx2*sn+y2*cs; | |
453 | ||
454 | } | |
455 | ||
456 | Double_t zz=z1 - z1/x1*(x1-x2); | |
457 | ||
bbf92647 | 458 | if (TMath::Abs(zz-z2)>fMaxSeedDeltaZ) continue; |
46d29e70 | 459 | |
460 | Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1); | |
461 | if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;} | |
462 | ||
463 | x[0]=y1; | |
464 | x[1]=z1; | |
465 | x[2]=f1trd(x1,y1,x2,y2,x3,y3); | |
466 | ||
bbf92647 | 467 | if (TMath::Abs(x[2]) >= fMaxSeedC) continue; |
46d29e70 | 468 | |
469 | x[3]=f2trd(x1,y1,x2,y2,x3,y3); | |
470 | ||
471 | if (TMath::Abs(x[2]*x1-x[3]) >= 0.99999) continue; | |
472 | ||
473 | x[4]=f3trd(x1,y1,x2,y2,z1,z2); | |
474 | ||
bbf92647 | 475 | if (TMath::Abs(x[4]) > fMaxSeedTan) continue; |
46d29e70 | 476 | |
477 | Double_t a=asin(x[3]); | |
478 | Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3])); | |
bbf92647 | 479 | if (TMath::Abs(zv)>fMaxSeedVertexZ) continue; |
46d29e70 | 480 | |
bbf92647 | 481 | Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2(); |
482 | Double_t sy2=cl->GetSigmaY2(), sz2=cl->GetSigmaZ2(); | |
483 | Double_t sy3=fSeedErrorSY3, sy=fSeedErrorSY, sz=fSeedErrorSZ; | |
46d29e70 | 484 | |
485 | Double_t f20=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy; | |
486 | Double_t f22=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy; | |
487 | Double_t f24=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy; | |
488 | Double_t f30=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[3])/sy; | |
489 | Double_t f32=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[3])/sy; | |
490 | Double_t f34=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[3])/sy; | |
491 | Double_t f40=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[4])/sy; | |
492 | Double_t f41=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[4])/sz; | |
493 | Double_t f42=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[4])/sy; | |
494 | Double_t f43=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[4])/sz; | |
495 | ||
496 | c[0]=sy1; | |
497 | c[1]=0.; c[2]=sz1; | |
498 | c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f24*sy3*f24; | |
499 | c[6]=f30*sy1; c[7]=0.; c[8]=f30*sy1*f20+f32*sy2*f22+f34*sy3*f24; | |
500 | c[9]=f30*sy1*f30+f32*sy2*f32+f34*sy3*f34; | |
501 | c[10]=f40*sy1; c[11]=f41*sz1; c[12]=f40*sy1*f20+f42*sy2*f22; | |
502 | c[13]=f40*sy1*f30+f42*sy2*f32; | |
503 | c[14]=f40*sy1*f40+f41*sz1*f41+f42*sy2*f42+f43*sz2*f43; | |
504 | ||
505 | UInt_t index=r1.GetIndex(is); | |
bbf92647 | 506 | AliTRDtrack *track=new AliTRDtrack(index, x, c, x1, ns*alpha+shift); |
46d29e70 | 507 | |
508 | Int_t rc=FindProlongation(*track,fTrSec,ns,i2); | |
509 | ||
bbf92647 | 510 | if ((rc < 0) || |
511 | (track->GetNclusters() < (i1-i2)*fMinClustersInSeed)) delete track; | |
46d29e70 | 512 | else { |
513 | fSeeds->AddLast(track); fNseeds++; | |
514 | cerr<<"found seed "<<fNseeds<<endl; | |
515 | } | |
516 | } | |
517 | } | |
518 | } | |
519 | ||
46d29e70 | 520 | fSeeds->Sort(); |
521 | } | |
522 | ||
523 | //___________________________________________________________________ | |
524 | void AliTRDtracker::FindTracks() | |
525 | { | |
526 | if (!fClusters) return; | |
527 | ||
dd56b762 | 528 | AliTRDtrackingSector fTrSec[AliTRDgeometry::kNsect]; |
46d29e70 | 529 | SetUpSectors(fTrSec); |
530 | ||
531 | // find tracks | |
532 | ||
533 | Int_t num_of_time_bins = fTrSec[0].GetNtimeBins(); | |
534 | Int_t nseed=fSeeds->GetEntriesFast(); | |
535 | ||
536 | Int_t nSeedClusters; | |
537 | for (Int_t i=0; i<nseed; i++) { | |
538 | cerr<<"FindTracks: seed "<<i+1<<" out of "<<nseed<<endl; | |
539 | ||
bbf92647 | 540 | AliTRDtrack& t=*((AliTRDtrack*)fSeeds->UncheckedAt(i)); |
46d29e70 | 541 | |
542 | nSeedClusters = t.GetNclusters(); | |
bbf92647 | 543 | Double_t alpha=t.GetAlpha(); |
46d29e70 | 544 | |
545 | if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); | |
546 | if (alpha < 0. ) alpha += 2.*TMath::Pi(); | |
dd56b762 | 547 | Int_t ns=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect; |
46d29e70 | 548 | |
549 | if (FindProlongation(t,fTrSec,ns)) { | |
550 | cerr<<"No of clusters in the track = "<<t.GetNclusters()<<endl; | |
bbf92647 | 551 | if (t.GetNclusters() >= Int_t(fMinClustersInTrack*num_of_time_bins)) { |
46d29e70 | 552 | Int_t label = GetTrackLabel(t); |
553 | t.SetLabel(label); | |
46d29e70 | 554 | UseClusters(t); |
bbf92647 | 555 | |
556 | AliTRDtrack *pt = new AliTRDtrack(t); | |
557 | fTracks->AddLast(pt); fNtracks++; | |
558 | ||
46d29e70 | 559 | cerr<<"found track "<<fNtracks<<endl; |
560 | } | |
561 | } | |
bbf92647 | 562 | delete fSeeds->RemoveAt(i); |
46d29e70 | 563 | } |
564 | } | |
565 | ||
566 | //__________________________________________________________________ | |
bbf92647 | 567 | void AliTRDtracker::UseClusters(AliTRDtrack t) { |
46d29e70 | 568 | Int_t ncl=t.GetNclusters(); |
569 | for (Int_t i=0; i<ncl; i++) { | |
570 | Int_t index = t.GetClusterIndex(i); | |
571 | AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index); | |
572 | c->Use(); | |
573 | } | |
574 | } | |
575 | ||
576 | //__________________________________________________________________ | |
bbf92647 | 577 | Int_t AliTRDtracker::GetTrackLabel(AliTRDtrack t) { |
46d29e70 | 578 | |
579 | Int_t label=123456789, index, i, j; | |
580 | Int_t ncl=t.GetNclusters(); | |
bbf92647 | 581 | const Int_t range = AliTRDgeometry::kNplan * fGeom->GetTimeMax(); |
46d29e70 | 582 | Bool_t label_added; |
583 | ||
d1b06c24 | 584 | // Int_t s[range][2]; |
585 | Int_t **s = new Int_t* [range]; | |
586 | for (i=0; i<range; i++) { | |
587 | s[i] = new Int_t[2]; | |
588 | } | |
46d29e70 | 589 | for (i=0; i<range; i++) { |
590 | s[i][0]=-1; | |
591 | s[i][1]=0; | |
592 | } | |
593 | ||
594 | Int_t t0,t1,t2; | |
595 | for (i=0; i<ncl; i++) { | |
596 | index=t.GetClusterIndex(i); | |
597 | AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index); | |
598 | t0=c->GetTrackIndex(0); | |
599 | t1=c->GetTrackIndex(1); | |
600 | t2=c->GetTrackIndex(2); | |
601 | } | |
602 | ||
603 | for (i=0; i<ncl; i++) { | |
604 | index=t.GetClusterIndex(i); | |
605 | AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index); | |
606 | for (Int_t k=0; k<3; k++) { | |
607 | label=c->GetTrackIndex(k); | |
57527628 | 608 | label_added=kFALSE; j=0; |
46d29e70 | 609 | if (label >= 0) { |
610 | while ( (!label_added) && ( j < range ) ) { | |
611 | if (s[j][0]==label || s[j][1]==0) { | |
612 | s[j][0]=label; | |
613 | s[j][1]=s[j][1]+1; | |
57527628 | 614 | label_added=kTRUE; |
46d29e70 | 615 | } |
616 | j++; | |
617 | } | |
618 | } | |
619 | } | |
620 | } | |
621 | ||
46d29e70 | 622 | Int_t max=0; |
623 | label = -123456789; | |
624 | ||
625 | for (i=0; i<range; i++) { | |
626 | if (s[i][1]>max) { | |
627 | max=s[i][1]; label=s[i][0]; | |
628 | } | |
629 | } | |
d1b06c24 | 630 | delete []s; |
bbf92647 | 631 | if(max > ncl*fLabelFraction) return label; |
46d29e70 | 632 | else return -1; |
633 | } | |
634 | ||
635 | //___________________________________________________________________ | |
636 | ||
bbf92647 | 637 | Int_t AliTRDtracker::WriteTracks(const Char_t *filename) { |
46d29e70 | 638 | |
639 | TDirectory *savedir=gDirectory; | |
640 | ||
bbf92647 | 641 | TFile *out=TFile::Open(filename,"RECREATE"); |
46d29e70 | 642 | |
643 | TTree tracktree("TreeT","Tree with TRD tracks"); | |
644 | ||
645 | AliTRDtrack *iotrack=0; | |
646 | tracktree.Branch("tracks","AliTRDtrack",&iotrack,32000,0); | |
647 | ||
648 | Int_t ntracks=fTracks->GetEntriesFast(); | |
649 | ||
650 | for (Int_t i=0; i<ntracks; i++) { | |
651 | AliTRDtrack *pt=(AliTRDtrack*)fTracks->UncheckedAt(i); | |
652 | iotrack=pt; | |
653 | tracktree.Fill(); | |
654 | cerr<<"WriteTracks: put track "<<i<<" in the tree"<<endl; | |
655 | } | |
656 | ||
657 | tracktree.Write(); | |
658 | out->Close(); | |
659 | ||
660 | savedir->cd(); | |
661 | ||
662 | cerr<<"WriteTracks: done"<<endl; | |
663 | return 0; | |
664 | } | |
665 | ||
46d29e70 | 666 | //_____________________________________________________________________________ |
bbf92647 | 667 | void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename, |
d1b06c24 | 668 | Int_t option) |
46d29e70 | 669 | { |
670 | // | |
671 | // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0) | |
672 | // from the file. The names of the cluster tree and branches | |
673 | // should match the ones used in AliTRDclusterizer::WriteClusters() | |
674 | // | |
675 | ||
676 | TDirectory *savedir=gDirectory; | |
677 | ||
678 | TFile *file = TFile::Open(filename); | |
679 | if (!file->IsOpen()) {printf("Can't open file %s !\n",filename); return;} | |
680 | ||
bbf92647 | 681 | TTree *tree = (TTree*)file->Get("ClusterTree"); |
46d29e70 | 682 | Int_t nentr = (Int_t) tree->GetEntries(); |
bbf92647 | 683 | printf("found %d entries in %s.\n",nentr,tree->GetName()); |
684 | ||
685 | TBranch *branch; | |
686 | TObjArray *ioArray = new TObjArray(400); | |
687 | ||
688 | if( option < 0 ) { | |
a2b90f83 | 689 | //branch = tree->GetBranch("RecPoints"); |
690 | // changed CBL | |
691 | branch = tree->GetBranch("TRDrecPoints"); | |
bbf92647 | 692 | |
693 | for (Int_t i=0; i<nentr; i++) { | |
694 | branch->SetAddress(&ioArray); | |
695 | tree->GetEvent(i); | |
696 | Int_t npoints = ioArray->GetEntriesFast(); | |
697 | printf("Read %d rec. points from entry %d \n", npoints, i); | |
698 | ||
699 | for(Int_t j=0; j<npoints; j++) { | |
700 | AliTRDrecPoint *p=(AliTRDrecPoint*)ioArray->UncheckedAt(j); | |
701 | array->AddLast(p); | |
702 | ioArray->RemoveAt(j); | |
703 | } | |
704 | } | |
705 | } | |
706 | else { | |
a2b90f83 | 707 | branch = tree->GetBranch("TRDcluster"); |
bbf92647 | 708 | |
709 | for (Int_t i=0; i<nentr; i++) { | |
710 | branch->SetAddress(&ioArray); | |
711 | tree->GetEvent(i); | |
712 | Int_t npoints = ioArray->GetEntriesFast(); | |
713 | printf("Read %d clusters from entry %d \n", npoints, i); | |
714 | ||
715 | for(Int_t j=0; j<npoints; j++) { | |
716 | AliTRDcluster *c=(AliTRDcluster*)ioArray->UncheckedAt(j); | |
717 | array->AddLast(c); | |
718 | ioArray->RemoveAt(j); | |
719 | } | |
46d29e70 | 720 | } |
721 | } | |
722 | ||
723 | file->Close(); | |
724 | savedir->cd(); | |
bbf92647 | 725 | |
46d29e70 | 726 | } |
727 |