]>
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$ | |
abaf1f1d | 18 | Revision 1.13 2001/05/30 12:17:47 hristov |
19 | Loop variables declared once | |
20 | ||
3ab6f951 | 21 | Revision 1.12 2001/05/28 17:07:58 hristov |
22 | Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh) | |
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 | ||
a819a5f7 | 73 | const Float_t AliTRDtracker::fSeedDepth = 0.5; |
74 | const Float_t AliTRDtracker::fSeedStep = 0.05; | |
75 | const Float_t AliTRDtracker::fSeedGap = 0.25; | |
76 | ||
77 | const Float_t AliTRDtracker::fMaxSeedDeltaZ12 = 40.; | |
78 | const Float_t AliTRDtracker::fMaxSeedDeltaZ = 25.; | |
79 | const Float_t AliTRDtracker::fMaxSeedC = 0.0052; | |
80 | const Float_t AliTRDtracker::fMaxSeedTan = 1.2; | |
81 | const Float_t AliTRDtracker::fMaxSeedVertexZ = 150.; | |
bbf92647 | 82 | |
a819a5f7 | 83 | const Double_t AliTRDtracker::fSeedErrorSY = 0.2; |
84 | const Double_t AliTRDtracker::fSeedErrorSY3 = 2.5; | |
85 | const Double_t AliTRDtracker::fSeedErrorSZ = 0.1; | |
bbf92647 | 86 | |
a819a5f7 | 87 | const Float_t AliTRDtracker::fMinClustersInSeed = 0.7; |
bbf92647 | 88 | |
a819a5f7 | 89 | const Float_t AliTRDtracker::fMinClustersInTrack = 0.5; |
90 | const Float_t AliTRDtracker::fSkipDepth = 0.05; | |
91 | const Float_t AliTRDtracker::fLabelFraction = 0.5; | |
92 | const Float_t AliTRDtracker::fWideRoad = 20.; | |
93 | ||
94 | const Double_t AliTRDtracker::fMaxChi2 = 24.; | |
bbf92647 | 95 | |
46d29e70 | 96 | //____________________________________________________________________ |
97 | AliTRDtracker::AliTRDtracker() | |
98 | { | |
99 | // | |
100 | // Default constructor | |
101 | // | |
102 | ||
46d29e70 | 103 | fEvent = 0; |
46d29e70 | 104 | fGeom = NULL; |
bbf92647 | 105 | |
46d29e70 | 106 | fNclusters = 0; |
107 | fClusters = NULL; | |
108 | fNseeds = 0; | |
109 | fSeeds = NULL; | |
110 | fNtracks = 0; | |
111 | fTracks = NULL; | |
112 | ||
a819a5f7 | 113 | fSY2corr = 0.025; |
46d29e70 | 114 | } |
115 | ||
116 | //____________________________________________________________________ | |
117 | AliTRDtracker::AliTRDtracker(const Text_t* name, const Text_t* title) | |
118 | :TNamed(name, title) | |
119 | { | |
46d29e70 | 120 | fEvent = 0; |
46d29e70 | 121 | fGeom = NULL; |
bbf92647 | 122 | |
46d29e70 | 123 | fNclusters = 0; |
124 | fClusters = new TObjArray(2000); | |
125 | fNseeds = 0; | |
126 | fSeeds = new TObjArray(20000); | |
127 | fNtracks = 0; | |
128 | fTracks = new TObjArray(10000); | |
129 | ||
a819a5f7 | 130 | fSY2corr = 0.025; |
46d29e70 | 131 | } |
132 | ||
46d29e70 | 133 | //___________________________________________________________________ |
134 | AliTRDtracker::~AliTRDtracker() | |
135 | { | |
46d29e70 | 136 | delete fClusters; |
137 | delete fTracks; | |
138 | delete fSeeds; | |
139 | delete fGeom; | |
140 | } | |
141 | ||
bbf92647 | 142 | //___________________________________________________________________ |
a819a5f7 | 143 | void AliTRDtracker::Clusters2Tracks(TH1F *hs, TH1F *hd) |
bbf92647 | 144 | { |
3ab6f951 | 145 | Int_t i; |
146 | ||
bbf92647 | 147 | Int_t inner, outer; |
148 | Int_t fTotalNofTimeBins = fGeom->GetTimeMax() * AliTRDgeometry::Nplan(); | |
a819a5f7 | 149 | Int_t nSteps = (Int_t) (fSeedDepth / fSeedStep); |
150 | Int_t gap = (Int_t) (fTotalNofTimeBins * fSeedGap); | |
151 | Int_t step = (Int_t) (fTotalNofTimeBins * fSeedStep); | |
152 | ||
153 | ||
154 | // nSteps = 1; | |
155 | ||
156 | if (!fClusters) return; | |
157 | ||
158 | AliTRDtrackingSector fTrSec[AliTRDgeometry::kNsect]; | |
159 | SetUpSectors(fTrSec); | |
160 | ||
161 | // make a first turn looking for seed ends in the same (n,n) | |
162 | // and in the adjacent sectors (n,n+1) | |
bbf92647 | 163 | |
3ab6f951 | 164 | for(i=0; i<nSteps; i++) { |
bbf92647 | 165 | printf("step %d out of %d \n", i+1, nSteps); |
a819a5f7 | 166 | outer=fTotalNofTimeBins-1-i*step; inner=outer-gap; |
167 | MakeSeeds(inner,outer, fTrSec, 1, hs, hd); | |
168 | FindTracks(fTrSec, hs, hd); | |
bbf92647 | 169 | } |
a819a5f7 | 170 | |
171 | // make a second turn looking for seed ends in next-to-adjacent | |
172 | // sectors (n,n+2) | |
173 | ||
3ab6f951 | 174 | for(i=0; i<nSteps; i++) { |
a819a5f7 | 175 | printf("step %d out of %d \n", i+1, nSteps); |
176 | outer=fTotalNofTimeBins-1-i*step; inner=outer-gap; | |
177 | MakeSeeds(inner, outer, fTrSec, 2, hs, hd); | |
178 | FindTracks(fTrSec,hs,hd); | |
179 | } | |
180 | ||
bbf92647 | 181 | } |
46d29e70 | 182 | |
183 | //_____________________________________________________________________ | |
bbf92647 | 184 | Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt) |
46d29e70 | 185 | { |
186 | // Parametrised "expected" error of the cluster reconstruction in Y | |
187 | ||
a819a5f7 | 188 | Double_t s = 0.08 * 0.08; |
46d29e70 | 189 | return s; |
190 | } | |
191 | ||
192 | //_____________________________________________________________________ | |
bbf92647 | 193 | Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl) |
46d29e70 | 194 | { |
195 | // Parametrised "expected" error of the cluster reconstruction in Z | |
196 | ||
a819a5f7 | 197 | Double_t s = 6 * 6 /12.; |
46d29e70 | 198 | return s; |
199 | } | |
200 | ||
201 | //_____________________________________________________________________ | |
bdbd0f7a | 202 | Double_t f1trd(Double_t x1,Double_t y1, |
a819a5f7 | 203 | Double_t x2,Double_t y2, |
204 | Double_t x3,Double_t y3) | |
46d29e70 | 205 | { |
206 | // Initial approximation of the track curvature | |
207 | // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch | |
208 | ||
209 | Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1); | |
210 | Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)- | |
211 | (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2)); | |
212 | Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)- | |
213 | (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1)); | |
214 | ||
215 | Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b); | |
216 | ||
217 | return -xr*yr/sqrt(xr*xr+yr*yr); | |
218 | } | |
219 | ||
220 | //_____________________________________________________________________ | |
bdbd0f7a | 221 | Double_t f2trd(Double_t x1,Double_t y1, |
a819a5f7 | 222 | Double_t x2,Double_t y2, |
223 | Double_t x3,Double_t y3) | |
46d29e70 | 224 | { |
225 | // Initial approximation of the track curvature times center of curvature | |
226 | // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch | |
227 | ||
228 | Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1); | |
229 | Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)- | |
230 | (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2)); | |
231 | Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)- | |
232 | (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1)); | |
233 | ||
234 | Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b); | |
235 | ||
236 | return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr); | |
237 | } | |
238 | ||
239 | //_____________________________________________________________________ | |
bdbd0f7a | 240 | Double_t f3trd(Double_t x1,Double_t y1, |
a819a5f7 | 241 | Double_t x2,Double_t y2, |
242 | Double_t z1,Double_t z2) | |
46d29e70 | 243 | { |
244 | // Initial approximation of the tangent of the track dip angle | |
245 | // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch | |
246 | ||
247 | return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); | |
248 | } | |
249 | ||
250 | ||
251 | //___________________________________________________________________ | |
252 | ||
bbf92647 | 253 | Int_t AliTRDtracker::FindProlongation(AliTRDtrack& t, AliTRDtrackingSector *sec, |
a819a5f7 | 254 | Int_t s, Int_t rf, Int_t matched_index, |
255 | TH1F *hs, TH1F *hd) | |
46d29e70 | 256 | { |
257 | // Starting from current position on track=t this function tries | |
258 | // to extrapolate the track up to timeBin=rf and to confirm prolongation | |
259 | // if a close cluster is found. *sec is a pointer to allocated | |
260 | // array of sectors, in which the initial sector has index=s. | |
261 | ||
a819a5f7 | 262 | |
263 | // TH1F *hsame = hs; | |
264 | // TH1F *hdiff = hd; | |
265 | ||
266 | // Bool_t good_match; | |
267 | ||
bbf92647 | 268 | const Int_t TIME_BINS_TO_SKIP=Int_t(fSkipDepth*sec->GetNtimeBins()); |
46d29e70 | 269 | Int_t try_again=TIME_BINS_TO_SKIP; |
270 | ||
0e9c2ad5 | 271 | Double_t alpha=AliTRDgeometry::GetAlpha(); |
46d29e70 | 272 | |
273 | Int_t ns=Int_t(2*TMath::Pi()/alpha+0.5); | |
274 | ||
a819a5f7 | 275 | Double_t x0, rho; |
276 | ||
46d29e70 | 277 | for (Int_t nr=sec->GetTimeBinNumber(t.GetX())-1; nr>=rf; nr--) { |
a819a5f7 | 278 | |
279 | Double_t x=sec->GetX(nr); | |
280 | Double_t ymax=x*TMath::Tan(0.5*alpha); | |
281 | ||
282 | rho = 0.00295; x0 = 11.09; // TEC | |
283 | if(sec->TECframe(nr,t.GetY(),t.GetZ())) { | |
284 | rho = 1.7; x0 = 33.0; // G10 frame | |
285 | } | |
286 | if(TMath::Abs(x - t.GetX()) > 3) { | |
287 | rho = 0.0559; x0 = 55.6; // radiator | |
288 | } | |
289 | if (!t.PropagateTo(x,x0,rho,0.139)) { | |
46d29e70 | 290 | cerr<<"Can't propagate to x = "<<x<<endl; |
291 | return 0; | |
292 | } | |
293 | ||
46d29e70 | 294 | AliTRDtimeBin& time_bin=sec[s][nr]; |
bbf92647 | 295 | Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt()); |
296 | Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl()); | |
a819a5f7 | 297 | Double_t road=25.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ(); |
46d29e70 | 298 | |
bbf92647 | 299 | if (road>fWideRoad) { |
46d29e70 | 300 | if (t.GetNclusters() > 4) { |
301 | cerr<<t.GetNclusters()<<" FindProlongation: Road is too wide !\n"; | |
302 | } | |
303 | return 0; | |
304 | } | |
305 | ||
a819a5f7 | 306 | AliTRDcluster *cl=0; |
307 | UInt_t index=0; | |
308 | // Int_t ncl = 0; | |
309 | ||
310 | Double_t max_chi2=fMaxChi2; | |
311 | ||
46d29e70 | 312 | if (time_bin) { |
a819a5f7 | 313 | |
46d29e70 | 314 | for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) { |
315 | AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]); | |
a819a5f7 | 316 | |
317 | // good_match = kFALSE; | |
318 | // if((c->GetTrackIndex(0) == matched_index) || | |
319 | // (c->GetTrackIndex(1) == matched_index) || | |
320 | // (c->GetTrackIndex(2) == matched_index)) good_match = kTRUE; | |
321 | // if(hsame) hsame->Fill(TMath::Abs(c->GetY()-y)/road); | |
322 | // if(hdiff) hdiff->Fill(road); | |
323 | ||
46d29e70 | 324 | if (c->GetY() > y+road) break; |
325 | if (c->IsUsed() > 0) continue; | |
326 | ||
a819a5f7 | 327 | // if(good_match) hsame->Fill(TMath::Abs(c->GetZ()-z)); |
328 | // else hdiff->Fill(TMath::Abs(c->GetZ()-z)); | |
329 | ||
330 | // if(!good_match) continue; | |
331 | ||
332 | if((c->GetZ()-z)*(c->GetZ()-z) > 3 * 12 * sz2) continue; | |
46d29e70 | 333 | |
334 | Double_t chi2=t.GetPredictedChi2(c); | |
335 | ||
a819a5f7 | 336 | // if((c->GetTrackIndex(0) == matched_index) || |
337 | // (c->GetTrackIndex(1) == matched_index) || | |
338 | // (c->GetTrackIndex(2) == matched_index)) | |
339 | // hdiff->Fill(chi2); | |
340 | ||
341 | // ncl++; | |
342 | ||
46d29e70 | 343 | if (chi2 > max_chi2) continue; |
344 | max_chi2=chi2; | |
345 | cl=c; | |
346 | index=time_bin.GetIndex(i); | |
347 | } | |
a819a5f7 | 348 | |
349 | if(!cl) { | |
350 | ||
351 | for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) { | |
352 | AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]); | |
353 | ||
354 | if (c->GetY() > y+road) break; | |
355 | if (c->IsUsed() > 0) continue; | |
356 | if((c->GetZ()-z)*(c->GetZ()-z) > 3.25 * 12 * sz2) continue; | |
357 | ||
358 | Double_t chi2=t.GetPredictedChi2(c); | |
359 | ||
360 | // ncl++; | |
361 | ||
362 | if (chi2 > max_chi2) continue; | |
363 | max_chi2=chi2; | |
364 | cl=c; | |
365 | index=time_bin.GetIndex(i); | |
366 | } | |
367 | } | |
368 | ||
46d29e70 | 369 | } |
a819a5f7 | 370 | |
46d29e70 | 371 | if (cl) { |
372 | ||
46d29e70 | 373 | t.Update(cl,max_chi2,index); |
374 | ||
375 | try_again=TIME_BINS_TO_SKIP; | |
376 | } else { | |
377 | if (try_again==0) break; | |
378 | if (y > ymax) { | |
a819a5f7 | 379 | // cerr<<"y > ymax: "<<y<<" > "<<ymax<<endl; |
46d29e70 | 380 | s = (s+1) % ns; |
381 | if (!t.Rotate(alpha)) { | |
382 | cerr<<"Failed to rotate, alpha = "<<alpha<<endl; | |
383 | return 0; | |
384 | } | |
385 | } else if (y <-ymax) { | |
a819a5f7 | 386 | // cerr<<"y < -ymax: "<<y<<" < "<<-ymax<<endl; |
46d29e70 | 387 | s = (s-1+ns) % ns; |
388 | if (!t.Rotate(-alpha)) { | |
389 | cerr<<"Failed to rotate, alpha = "<<alpha<<endl; | |
390 | return 0; | |
391 | } | |
392 | } | |
a819a5f7 | 393 | if(!sec->TECframe(nr,y,z)) try_again--; |
46d29e70 | 394 | } |
395 | } | |
396 | ||
397 | return 1; | |
398 | } | |
399 | ||
400 | ||
401 | ||
402 | //_____________________________________________________________________________ | |
bbf92647 | 403 | void AliTRDtracker::GetEvent(const Char_t *hitfile, const Char_t *clusterfile) |
46d29e70 | 404 | { |
a819a5f7 | 405 | // Opens a ROOT-file with TRD-clusters and reads the cluster-tree in |
46d29e70 | 406 | |
bbf92647 | 407 | ReadClusters(fClusters, clusterfile); |
408 | ||
409 | // get geometry from the file with hits | |
410 | ||
411 | TFile *fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(hitfile); | |
46d29e70 | 412 | if (!fInputFile) { |
413 | printf("AliTRDtracker::Open -- "); | |
bbf92647 | 414 | printf("Open the ALIROOT-file %s.\n",hitfile); |
a819a5f7 | 415 | fInputFile = new TFile(hitfile,"UPDATE"); |
46d29e70 | 416 | } |
417 | else { | |
418 | printf("AliTRDtracker::Open -- "); | |
bbf92647 | 419 | printf("%s is already open.\n",hitfile); |
46d29e70 | 420 | } |
421 | ||
422 | // Get AliRun object from file or create it if not on file | |
46d29e70 | 423 | |
bbf92647 | 424 | gAlice = (AliRun*) fInputFile->Get("gAlice"); |
425 | if (gAlice) { | |
426 | printf("AliTRDtracker::GetEvent -- "); | |
427 | printf("AliRun object found on file.\n"); | |
428 | } | |
429 | else { | |
430 | printf("AliTRDtracker::GetEvent -- "); | |
431 | printf("Could not find AliRun object.\n"); | |
432 | } | |
46d29e70 | 433 | |
a819a5f7 | 434 | /* |
46d29e70 | 435 | // Import the Trees for the event nEvent in the file |
436 | Int_t nparticles = gAlice->GetEvent(fEvent); | |
437 | cerr<<"nparticles = "<<nparticles<<endl; | |
438 | ||
439 | if (nparticles <= 0) { | |
440 | printf("AliTRDtracker::GetEvent -- "); | |
441 | printf("No entries in the trees for event %d.\n",fEvent); | |
442 | } | |
a819a5f7 | 443 | */ |
46d29e70 | 444 | |
445 | AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD"); | |
446 | fGeom = TRD->GetGeometry(); | |
447 | ||
46d29e70 | 448 | } |
449 | ||
450 | ||
451 | //_____________________________________________________________________________ | |
452 | void AliTRDtracker::SetUpSectors(AliTRDtrackingSector *sec) | |
453 | { | |
454 | // Fills clusters into TRD tracking_sectors | |
455 | // Note that the numbering scheme for the TRD tracking_sectors | |
456 | // differs from that of TRD sectors | |
457 | ||
bbf92647 | 458 | for (Int_t i=0; i<AliTRDgeometry::Nsect(); i++) sec[i].SetUp(); |
46d29e70 | 459 | |
460 | // Sort clusters into AliTRDtimeBin's within AliTRDtrackSector's | |
461 | ||
a819a5f7 | 462 | cerr<<"SetUpSectors: sorting clusters"<<endl; |
46d29e70 | 463 | |
464 | Int_t ncl=fClusters->GetEntriesFast(); | |
465 | UInt_t index; | |
466 | while (ncl--) { | |
a819a5f7 | 467 | printf("\r %d left ",ncl); |
46d29e70 | 468 | AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl); |
469 | Int_t detector=c->GetDetector(), local_time_bin=c->GetLocalTimeBin(); | |
470 | Int_t sector=fGeom->GetSector(detector); | |
471 | ||
a819a5f7 | 472 | Int_t tracking_sector = AliTRDgeometry::kNsect - sector - 1; |
46d29e70 | 473 | |
474 | Int_t tb=sec[sector].GetTimeBin(detector,local_time_bin); | |
475 | index=ncl; | |
476 | sec[tracking_sector][tb].InsertCluster(c,index); | |
a819a5f7 | 477 | |
46d29e70 | 478 | } |
479 | printf("\r\n"); | |
480 | } | |
481 | ||
482 | ||
483 | //_____________________________________________________________________________ | |
a819a5f7 | 484 | void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, |
485 | AliTRDtrackingSector* fTrSec, Int_t turn, | |
486 | TH1F *hs, TH1F *hd) | |
46d29e70 | 487 | { |
488 | // Creates track seeds using clusters in timeBins=i1,i2 | |
489 | ||
bbf92647 | 490 | Int_t i2 = inner, i1 = outer; |
a819a5f7 | 491 | Int_t ti[3], to[3]; |
492 | Int_t nprim = 85600/2; | |
46d29e70 | 493 | |
46d29e70 | 494 | |
a819a5f7 | 495 | TH1F *hsame = hs; |
496 | TH1F *hdiff = hd; | |
497 | Bool_t match = false; | |
498 | Int_t matched_index; | |
46d29e70 | 499 | |
500 | // find seeds | |
501 | ||
502 | Double_t x[5], c[15]; | |
dd56b762 | 503 | Int_t max_sec=AliTRDgeometry::kNsect; |
46d29e70 | 504 | |
0e9c2ad5 | 505 | Double_t alpha=AliTRDgeometry::GetAlpha(); |
506 | Double_t shift=AliTRDgeometry::GetAlpha()/2.; | |
46d29e70 | 507 | Double_t cs=cos(alpha), sn=sin(alpha); |
a819a5f7 | 508 | Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha); |
46d29e70 | 509 | |
510 | Double_t x1 =fTrSec[0].GetX(i1); | |
511 | Double_t xx2=fTrSec[0].GetX(i2); | |
512 | ||
a819a5f7 | 513 | |
514 | printf("\n"); | |
515 | ||
516 | if((turn != 1)&&(turn != 2)) { | |
517 | printf("*** Error in MakeSeeds: unexpected turn = %d \n", turn); | |
518 | return; | |
519 | } | |
520 | ||
521 | ||
46d29e70 | 522 | for (Int_t ns=0; ns<max_sec; ns++) { |
523 | ||
a819a5f7 | 524 | printf("\n MakeSeeds: sector %d \n", ns); |
525 | ||
526 | Int_t nl2=fTrSec[(ns-2+max_sec)%max_sec][i2]; | |
46d29e70 | 527 | Int_t nl=fTrSec[(ns-1+max_sec)%max_sec][i2]; |
528 | Int_t nm=fTrSec[ns][i2]; | |
529 | Int_t nu=fTrSec[(ns+1)%max_sec][i2]; | |
a819a5f7 | 530 | Int_t nu2=fTrSec[(ns+2)%max_sec][i2]; |
46d29e70 | 531 | |
532 | AliTRDtimeBin& r1=fTrSec[ns][i1]; | |
533 | ||
534 | for (Int_t is=0; is < r1; is++) { | |
535 | Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ(); | |
a819a5f7 | 536 | for(Int_t ii=0; ii<3; ii++) to[ii] = r1[is]->GetTrackIndex(ii); |
46d29e70 | 537 | |
a819a5f7 | 538 | for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) { |
539 | ||
46d29e70 | 540 | const AliTRDcluster *cl; |
a819a5f7 | 541 | Double_t x2, y2, z2; |
542 | Double_t x3=0., y3=0.; | |
46d29e70 | 543 | |
a819a5f7 | 544 | if (js<nl2) { |
545 | if(turn != 2) continue; | |
546 | AliTRDtimeBin& r2=fTrSec[(ns-2+max_sec)%max_sec][i2]; | |
547 | cl=r2[js]; | |
548 | y2=cl->GetY(); z2=cl->GetZ(); | |
549 | for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii); | |
550 | ||
551 | x2= xx2*cs2+y2*sn2; | |
552 | y2=-xx2*sn2+y2*cs2; | |
553 | } | |
554 | else if (js<nl2+nl) { | |
555 | if(turn != 1) continue; | |
46d29e70 | 556 | AliTRDtimeBin& r2=fTrSec[(ns-1+max_sec)%max_sec][i2]; |
a819a5f7 | 557 | cl=r2[js-nl2]; |
46d29e70 | 558 | y2=cl->GetY(); z2=cl->GetZ(); |
a819a5f7 | 559 | for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii); |
46d29e70 | 560 | |
a819a5f7 | 561 | x2= xx2*cs+y2*sn; |
562 | y2=-xx2*sn+y2*cs; | |
46d29e70 | 563 | |
a819a5f7 | 564 | } |
565 | else if (js<nl2+nl+nm) { | |
566 | if(turn != 1) continue; | |
567 | AliTRDtimeBin& r2=fTrSec[ns][i2]; | |
568 | cl=r2[js-nl2-nl]; | |
569 | x2=xx2; y2=cl->GetY(); z2=cl->GetZ(); | |
570 | for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii); | |
571 | } | |
572 | else if (js<nl2+nl+nm+nu) { | |
573 | if(turn != 1) continue; | |
574 | AliTRDtimeBin& r2=fTrSec[(ns+1)%max_sec][i2]; | |
575 | cl=r2[js-nl2-nl-nm]; | |
576 | y2=cl->GetY(); z2=cl->GetZ(); | |
577 | for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii); | |
46d29e70 | 578 | |
a819a5f7 | 579 | x2=xx2*cs-y2*sn; |
580 | y2=xx2*sn+y2*cs; | |
46d29e70 | 581 | |
a819a5f7 | 582 | } |
583 | else { | |
584 | if(turn != 2) continue; | |
585 | AliTRDtimeBin& r2=fTrSec[(ns+2)%max_sec][i2]; | |
586 | cl=r2[js-nl2-nl-nm-nu]; | |
587 | y2=cl->GetY(); z2=cl->GetZ(); | |
588 | for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii); | |
589 | ||
590 | x2=xx2*cs2-y2*sn2; | |
591 | y2=xx2*sn2+y2*cs2; | |
592 | } | |
593 | ||
46d29e70 | 594 | |
a819a5f7 | 595 | match = false; |
596 | matched_index = -1; | |
597 | for (Int_t ii=0; ii<3; ii++) { | |
598 | // cerr<<"ti["<<ii<<"] = "<<ti[ii]<<"; to["<<ii<<"] = "<<to[ii]<<endl; | |
599 | if(ti[ii] < 0) continue; | |
600 | if(ti[ii] >= nprim) continue; | |
601 | for (Int_t kk=0; kk<3; kk++) { | |
602 | if(to[kk] < 0) continue; | |
603 | if(to[kk] >= nprim) continue; | |
604 | if(ti[ii] == to[kk]) { | |
605 | //cerr<<"ti["<<ii<<"] = "<<ti[ii]<<" = "<<to[kk]<<" = to["<<kk<<"]"<<endl; | |
606 | matched_index = ti[ii]; | |
607 | match = true; | |
608 | } | |
609 | } | |
610 | } | |
46d29e70 | 611 | |
a819a5f7 | 612 | if(TMath::Abs(z1-z2) > fMaxSeedDeltaZ12) continue; |
613 | ||
614 | Double_t zz=z1 - z1/x1*(x1-x2); | |
615 | ||
bbf92647 | 616 | if (TMath::Abs(zz-z2)>fMaxSeedDeltaZ) continue; |
46d29e70 | 617 | |
618 | Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1); | |
619 | if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;} | |
620 | ||
621 | x[0]=y1; | |
622 | x[1]=z1; | |
623 | x[2]=f1trd(x1,y1,x2,y2,x3,y3); | |
624 | ||
a819a5f7 | 625 | if (TMath::Abs(x[2]) > fMaxSeedC) continue; |
46d29e70 | 626 | |
627 | x[3]=f2trd(x1,y1,x2,y2,x3,y3); | |
628 | ||
629 | if (TMath::Abs(x[2]*x1-x[3]) >= 0.99999) continue; | |
630 | ||
631 | x[4]=f3trd(x1,y1,x2,y2,z1,z2); | |
632 | ||
bbf92647 | 633 | if (TMath::Abs(x[4]) > fMaxSeedTan) continue; |
46d29e70 | 634 | |
635 | Double_t a=asin(x[3]); | |
636 | Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3])); | |
a819a5f7 | 637 | |
bbf92647 | 638 | if (TMath::Abs(zv)>fMaxSeedVertexZ) continue; |
46d29e70 | 639 | |
bbf92647 | 640 | Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2(); |
641 | Double_t sy2=cl->GetSigmaY2(), sz2=cl->GetSigmaZ2(); | |
642 | Double_t sy3=fSeedErrorSY3, sy=fSeedErrorSY, sz=fSeedErrorSZ; | |
46d29e70 | 643 | |
644 | Double_t f20=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy; | |
645 | Double_t f22=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy; | |
646 | Double_t f24=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy; | |
647 | Double_t f30=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[3])/sy; | |
648 | Double_t f32=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[3])/sy; | |
649 | Double_t f34=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[3])/sy; | |
650 | Double_t f40=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[4])/sy; | |
651 | Double_t f41=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[4])/sz; | |
652 | Double_t f42=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[4])/sy; | |
653 | Double_t f43=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[4])/sz; | |
654 | ||
655 | c[0]=sy1; | |
656 | c[1]=0.; c[2]=sz1; | |
657 | c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f24*sy3*f24; | |
658 | c[6]=f30*sy1; c[7]=0.; c[8]=f30*sy1*f20+f32*sy2*f22+f34*sy3*f24; | |
659 | c[9]=f30*sy1*f30+f32*sy2*f32+f34*sy3*f34; | |
660 | c[10]=f40*sy1; c[11]=f41*sz1; c[12]=f40*sy1*f20+f42*sy2*f22; | |
661 | c[13]=f40*sy1*f30+f42*sy2*f32; | |
662 | c[14]=f40*sy1*f40+f41*sz1*f41+f42*sy2*f42+f43*sz2*f43; | |
663 | ||
664 | UInt_t index=r1.GetIndex(is); | |
a819a5f7 | 665 | |
666 | AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift); | |
667 | ||
668 | Int_t rc=FindProlongation(*track,fTrSec,ns,i2,matched_index,hsame,hdiff); | |
669 | ||
670 | // if (match) hsame->Fill((Float_t) track->GetNclusters()); | |
671 | // else hdiff->Fill((Float_t) track->GetNclusters()); | |
672 | // delete track; | |
673 | // continue; | |
46d29e70 | 674 | |
bbf92647 | 675 | if ((rc < 0) || |
676 | (track->GetNclusters() < (i1-i2)*fMinClustersInSeed)) delete track; | |
46d29e70 | 677 | else { |
678 | fSeeds->AddLast(track); fNseeds++; | |
a819a5f7 | 679 | printf("\r found seed %d ", fNseeds); |
46d29e70 | 680 | } |
681 | } | |
682 | } | |
683 | } | |
684 | ||
46d29e70 | 685 | fSeeds->Sort(); |
686 | } | |
687 | ||
a819a5f7 | 688 | //_____________________________________________________________________________ |
689 | void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename) | |
46d29e70 | 690 | { |
a819a5f7 | 691 | // |
692 | // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0) | |
693 | // from the file. The names of the cluster tree and branches | |
694 | // should match the ones used in AliTRDclusterizer::WriteClusters() | |
695 | // | |
46d29e70 | 696 | |
a819a5f7 | 697 | TDirectory *savedir=gDirectory; |
698 | ||
699 | TFile *file = TFile::Open(filename); | |
700 | if (!file->IsOpen()) {printf("Can't open file %s !\n",filename); return;} | |
701 | ||
abaf1f1d | 702 | Char_t treeName[12]; |
703 | sprintf(treeName,"TreeR%d_TRD",fEvent); | |
704 | TTree *ClusterTree = (TTree*) file->Get(treeName); | |
46d29e70 | 705 | |
a819a5f7 | 706 | TObjArray *ClusterArray = new TObjArray(400); |
707 | ||
708 | ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray); | |
709 | ||
710 | Int_t nEntries = (Int_t) ClusterTree->GetEntries(); | |
711 | printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName()); | |
712 | ||
713 | // Loop through all entries in the tree | |
714 | Int_t nbytes; | |
715 | AliTRDcluster *c = 0; | |
716 | ||
717 | for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) { | |
718 | ||
719 | // Import the tree | |
720 | nbytes += ClusterTree->GetEvent(iEntry); | |
721 | ||
722 | // Get the number of points in the detector | |
723 | Int_t nCluster = ClusterArray->GetEntriesFast(); | |
724 | printf("Read %d clusters from entry %d \n", nCluster, iEntry); | |
725 | ||
726 | // Loop through all TRD digits | |
727 | for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { | |
728 | c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster); | |
729 | AliTRDcluster *co = new AliTRDcluster(*c); | |
730 | co->SetSigmaY2(c->GetSigmaY2() * fSY2corr); | |
731 | array->AddLast(co); | |
732 | delete ClusterArray->RemoveAt(iCluster); | |
733 | } | |
734 | } | |
735 | ||
736 | file->Close(); | |
737 | delete ClusterArray; | |
738 | savedir->cd(); | |
739 | ||
740 | } | |
741 | ||
742 | //___________________________________________________________________ | |
743 | void AliTRDtracker::FindTracks(AliTRDtrackingSector* fTrSec, TH1F *hs, TH1F *hd) | |
744 | { | |
745 | // | |
746 | // Finds tracks in TRD | |
747 | // | |
748 | ||
749 | TH1F *hsame = hs; | |
750 | TH1F *hdiff = hd; | |
46d29e70 | 751 | |
752 | Int_t num_of_time_bins = fTrSec[0].GetNtimeBins(); | |
753 | Int_t nseed=fSeeds->GetEntriesFast(); | |
754 | ||
755 | Int_t nSeedClusters; | |
756 | for (Int_t i=0; i<nseed; i++) { | |
757 | cerr<<"FindTracks: seed "<<i+1<<" out of "<<nseed<<endl; | |
758 | ||
bbf92647 | 759 | AliTRDtrack& t=*((AliTRDtrack*)fSeeds->UncheckedAt(i)); |
46d29e70 | 760 | |
761 | nSeedClusters = t.GetNclusters(); | |
bbf92647 | 762 | Double_t alpha=t.GetAlpha(); |
46d29e70 | 763 | |
764 | if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); | |
765 | if (alpha < 0. ) alpha += 2.*TMath::Pi(); | |
dd56b762 | 766 | Int_t ns=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect; |
46d29e70 | 767 | |
a819a5f7 | 768 | Int_t label = GetTrackLabel(t); |
769 | ||
770 | if (FindProlongation(t,fTrSec,ns,0,label,hsame,hdiff)) { | |
46d29e70 | 771 | cerr<<"No of clusters in the track = "<<t.GetNclusters()<<endl; |
bbf92647 | 772 | if (t.GetNclusters() >= Int_t(fMinClustersInTrack*num_of_time_bins)) { |
46d29e70 | 773 | Int_t label = GetTrackLabel(t); |
774 | t.SetLabel(label); | |
a819a5f7 | 775 | t.CookdEdx(); |
46d29e70 | 776 | UseClusters(t); |
bbf92647 | 777 | |
778 | AliTRDtrack *pt = new AliTRDtrack(t); | |
779 | fTracks->AddLast(pt); fNtracks++; | |
780 | ||
46d29e70 | 781 | cerr<<"found track "<<fNtracks<<endl; |
782 | } | |
783 | } | |
bbf92647 | 784 | delete fSeeds->RemoveAt(i); |
a819a5f7 | 785 | fNseeds--; |
46d29e70 | 786 | } |
787 | } | |
788 | ||
789 | //__________________________________________________________________ | |
bbf92647 | 790 | void AliTRDtracker::UseClusters(AliTRDtrack t) { |
46d29e70 | 791 | Int_t ncl=t.GetNclusters(); |
792 | for (Int_t i=0; i<ncl; i++) { | |
793 | Int_t index = t.GetClusterIndex(i); | |
794 | AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index); | |
795 | c->Use(); | |
796 | } | |
797 | } | |
798 | ||
799 | //__________________________________________________________________ | |
bbf92647 | 800 | Int_t AliTRDtracker::GetTrackLabel(AliTRDtrack t) { |
46d29e70 | 801 | |
802 | Int_t label=123456789, index, i, j; | |
803 | Int_t ncl=t.GetNclusters(); | |
bbf92647 | 804 | const Int_t range = AliTRDgeometry::kNplan * fGeom->GetTimeMax(); |
46d29e70 | 805 | Bool_t label_added; |
806 | ||
d1b06c24 | 807 | // Int_t s[range][2]; |
808 | Int_t **s = new Int_t* [range]; | |
809 | for (i=0; i<range; i++) { | |
810 | s[i] = new Int_t[2]; | |
811 | } | |
46d29e70 | 812 | for (i=0; i<range; i++) { |
813 | s[i][0]=-1; | |
814 | s[i][1]=0; | |
815 | } | |
816 | ||
817 | Int_t t0,t1,t2; | |
818 | for (i=0; i<ncl; i++) { | |
819 | index=t.GetClusterIndex(i); | |
820 | AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index); | |
821 | t0=c->GetTrackIndex(0); | |
822 | t1=c->GetTrackIndex(1); | |
823 | t2=c->GetTrackIndex(2); | |
824 | } | |
825 | ||
826 | for (i=0; i<ncl; i++) { | |
827 | index=t.GetClusterIndex(i); | |
828 | AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index); | |
829 | for (Int_t k=0; k<3; k++) { | |
830 | label=c->GetTrackIndex(k); | |
57527628 | 831 | label_added=kFALSE; j=0; |
46d29e70 | 832 | if (label >= 0) { |
833 | while ( (!label_added) && ( j < range ) ) { | |
834 | if (s[j][0]==label || s[j][1]==0) { | |
835 | s[j][0]=label; | |
836 | s[j][1]=s[j][1]+1; | |
57527628 | 837 | label_added=kTRUE; |
46d29e70 | 838 | } |
839 | j++; | |
840 | } | |
841 | } | |
842 | } | |
843 | } | |
844 | ||
46d29e70 | 845 | Int_t max=0; |
846 | label = -123456789; | |
847 | ||
848 | for (i=0; i<range; i++) { | |
849 | if (s[i][1]>max) { | |
850 | max=s[i][1]; label=s[i][0]; | |
851 | } | |
852 | } | |
d1b06c24 | 853 | delete []s; |
bbf92647 | 854 | if(max > ncl*fLabelFraction) return label; |
46d29e70 | 855 | else return -1; |
856 | } | |
857 | ||
858 | //___________________________________________________________________ | |
859 | ||
bbf92647 | 860 | Int_t AliTRDtracker::WriteTracks(const Char_t *filename) { |
46d29e70 | 861 | |
862 | TDirectory *savedir=gDirectory; | |
863 | ||
a819a5f7 | 864 | //TFile *out=TFile::Open(filename,"RECREATE"); |
865 | TFile *out = (TFile*) gROOT->GetListOfFiles()->FindObject(filename); | |
866 | if (!out) { | |
867 | printf("AliTRDtracker::Open -- "); | |
868 | printf("Open the ALIROOT-file %s.\n",filename); | |
869 | out = new TFile(filename,"RECREATE"); | |
870 | } | |
871 | else { | |
872 | printf("AliTRDtracker::Open -- "); | |
873 | printf("%s is already open.\n",filename); | |
874 | } | |
46d29e70 | 875 | |
abaf1f1d | 876 | Char_t treeName[12]; |
877 | sprintf(treeName,"TreeT%d_TRD",fEvent); | |
878 | TTree tracktree(treeName,"Tree with TRD tracks"); | |
46d29e70 | 879 | |
880 | AliTRDtrack *iotrack=0; | |
881 | tracktree.Branch("tracks","AliTRDtrack",&iotrack,32000,0); | |
882 | ||
883 | Int_t ntracks=fTracks->GetEntriesFast(); | |
884 | ||
885 | for (Int_t i=0; i<ntracks; i++) { | |
886 | AliTRDtrack *pt=(AliTRDtrack*)fTracks->UncheckedAt(i); | |
887 | iotrack=pt; | |
888 | tracktree.Fill(); | |
889 | cerr<<"WriteTracks: put track "<<i<<" in the tree"<<endl; | |
890 | } | |
891 | ||
892 | tracktree.Write(); | |
893 | out->Close(); | |
894 | ||
895 | savedir->cd(); | |
896 | ||
897 | cerr<<"WriteTracks: done"<<endl; | |
a819a5f7 | 898 | return 1; |
46d29e70 | 899 | } |
900 |