]>
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$ | |
ca6c93d6 | 18 | Revision 1.16 2002/06/12 09:54:36 cblume |
19 | Update of tracking code provided by Sergei | |
20 | ||
0a29d0f1 | 21 | Revision 1.14 2001/11/14 10:50:46 cblume |
22 | Changes in digits IO. Add merging of summable digits | |
23 | ||
abaf1f1d | 24 | Revision 1.13 2001/05/30 12:17:47 hristov |
25 | Loop variables declared once | |
26 | ||
3ab6f951 | 27 | Revision 1.12 2001/05/28 17:07:58 hristov |
28 | 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) | |
29 | ||
71d9fa7b | 30 | Revision 1.8 2000/12/20 13:00:44 cblume |
31 | Modifications for the HP-compiler | |
32 | ||
d1b06c24 | 33 | Revision 1.7 2000/12/08 16:07:02 cblume |
34 | Update of the tracking by Sergei | |
35 | ||
bbf92647 | 36 | Revision 1.6 2000/11/30 17:38:08 cblume |
37 | Changes to get in line with new STEER and EVGEN | |
38 | ||
1819f4bb | 39 | Revision 1.5 2000/11/14 14:40:27 cblume |
40 | Correction for the Sun compiler (kTRUE and kFALSE) | |
41 | ||
57527628 | 42 | Revision 1.4 2000/11/10 14:57:52 cblume |
43 | Changes in the geometry constants for the DEC compiler | |
44 | ||
dd56b762 | 45 | Revision 1.3 2000/10/15 23:40:01 cblume |
46 | Remove AliTRDconst | |
47 | ||
0e9c2ad5 | 48 | Revision 1.2 2000/10/06 16:49:46 cblume |
49 | Made Getters const | |
50 | ||
46d29e70 | 51 | Revision 1.1.2.2 2000/10/04 16:34:58 cblume |
52 | Replace include files by forward declarations | |
53 | ||
54 | Revision 1.1.2.1 2000/09/22 14:47:52 cblume | |
55 | Add the tracking code | |
56 | ||
57 | */ | |
bbf92647 | 58 | |
1819f4bb | 59 | #include <iostream.h> |
46d29e70 | 60 | |
61 | #include <TFile.h> | |
46d29e70 | 62 | #include <TBranch.h> |
5443e65e | 63 | #include <TTree.h> |
64 | #include <TObjArray.h> | |
46d29e70 | 65 | |
46d29e70 | 66 | #include "AliTRDgeometry.h" |
5443e65e | 67 | #include "AliTRDparameter.h" |
68 | #include "AliTRDgeometryDetail.h" | |
46d29e70 | 69 | #include "AliTRDcluster.h" |
70 | #include "AliTRDtrack.h" | |
5443e65e | 71 | #include "../TPC/AliTPCtrack.h" |
46d29e70 | 72 | |
73 | #include "AliTRDtracker.h" | |
74 | ||
75 | ClassImp(AliTRDtracker) | |
76 | ||
5443e65e | 77 | const Float_t AliTRDtracker::fSeedDepth = 0.5; |
78 | const Float_t AliTRDtracker::fSeedStep = 0.10; | |
79 | const Float_t AliTRDtracker::fSeedGap = 0.25; | |
80 | ||
81 | const Float_t AliTRDtracker::fMaxSeedDeltaZ12 = 40.; | |
82 | const Float_t AliTRDtracker::fMaxSeedDeltaZ = 25.; | |
83 | const Float_t AliTRDtracker::fMaxSeedC = 0.0052; | |
84 | const Float_t AliTRDtracker::fMaxSeedTan = 1.2; | |
85 | const Float_t AliTRDtracker::fMaxSeedVertexZ = 150.; | |
a819a5f7 | 86 | |
5443e65e | 87 | const Double_t AliTRDtracker::fSeedErrorSY = 0.2; |
88 | const Double_t AliTRDtracker::fSeedErrorSY3 = 2.5; | |
89 | const Double_t AliTRDtracker::fSeedErrorSZ = 0.1; | |
bbf92647 | 90 | |
5443e65e | 91 | const Float_t AliTRDtracker::fMinClustersInSeed = 0.7; |
bbf92647 | 92 | |
5443e65e | 93 | const Float_t AliTRDtracker::fMinClustersInTrack = 0.5; |
94 | const Float_t AliTRDtracker::fMinFractionOfFoundClusters = 0.8; | |
bbf92647 | 95 | |
5443e65e | 96 | const Float_t AliTRDtracker::fSkipDepth = 0.05; |
97 | const Float_t AliTRDtracker::fLabelFraction = 0.8; | |
98 | const Float_t AliTRDtracker::fWideRoad = 20.; | |
99 | ||
100 | const Double_t AliTRDtracker::fMaxChi2 = 24.; | |
a819a5f7 | 101 | |
bbf92647 | 102 | |
46d29e70 | 103 | //____________________________________________________________________ |
5443e65e | 104 | AliTRDtracker::AliTRDtracker(const TFile *geomfile) |
46d29e70 | 105 | { |
5443e65e | 106 | // |
107 | // Main constructor | |
108 | // | |
46d29e70 | 109 | |
5443e65e | 110 | Float_t fTzero = 0; |
111 | Int_t version = 2; | |
112 | ||
113 | fAddTRDseeds = kFALSE; | |
bbf92647 | 114 | |
5443e65e | 115 | fGeom = NULL; |
116 | ||
117 | TDirectory *savedir=gDirectory; | |
118 | TFile *in=(TFile*)geomfile; | |
119 | if (!in->IsOpen()) { | |
120 | printf("AliTRDtracker::AliTRDtracker(): geometry file is not open!\n"); | |
121 | printf(" DETAIL TRD geometry and DEFAULT TRD parameter will be used\n"); | |
122 | } | |
123 | else { | |
124 | in->cd(); | |
125 | in->ls(); | |
126 | fGeom = (AliTRDgeometry*) in->Get("TRDgeometry"); | |
127 | fPar = (AliTRDparameter*) in->Get("TRDparameter"); | |
128 | fGeom->Dump(); | |
129 | } | |
46d29e70 | 130 | |
5443e65e | 131 | if(fGeom) { |
132 | // fTzero = geo->GetT0(); | |
133 | fTzero = 0.; | |
134 | version = fGeom->IsVersion(); | |
135 | printf("Found geometry version %d on file \n", version); | |
136 | } | |
137 | else { | |
138 | printf("AliTRDtracker::AliTRDtracker(): cann't find TRD geometry!\n"); | |
139 | printf(" DETAIL TRD geometry and DEFAULT TRD parameter will be used\n"); | |
140 | fGeom = new AliTRDgeometryDetail(); | |
141 | fPar = new AliTRDparameter(); | |
142 | } | |
143 | ||
144 | savedir->cd(); | |
46d29e70 | 145 | |
5443e65e | 146 | |
147 | // fGeom->SetT0(fTzero); | |
0a29d0f1 | 148 | |
46d29e70 | 149 | fEvent = 0; |
bbf92647 | 150 | |
46d29e70 | 151 | fNclusters = 0; |
152 | fClusters = new TObjArray(2000); | |
153 | fNseeds = 0; | |
5443e65e | 154 | fSeeds = new TObjArray(2000); |
46d29e70 | 155 | fNtracks = 0; |
5443e65e | 156 | fTracks = new TObjArray(1000); |
a819a5f7 | 157 | |
5443e65e | 158 | for(Int_t geom_s = 0; geom_s < kTRACKING_SECTORS; geom_s++) { |
159 | Int_t tr_s = CookSectorIndex(geom_s); | |
160 | fTrSec[tr_s] = new AliTRDtrackingSector(fGeom, geom_s, fPar); | |
161 | } | |
a819a5f7 | 162 | |
5443e65e | 163 | fSY2corr = 0.025; |
164 | fSZ2corr = 12.; | |
bbf92647 | 165 | |
5443e65e | 166 | // calculate max gap on track |
a819a5f7 | 167 | |
5443e65e | 168 | Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region |
169 | Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region | |
a819a5f7 | 170 | |
5443e65e | 171 | Double_t dx = (Double_t) fPar->GetTimeBinSize(); |
172 | Int_t tbAmp = fPar->GetTimeBefore(); | |
173 | Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx); | |
174 | Int_t tbDrift = fPar->GetTimeMax(); | |
175 | Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx); | |
a819a5f7 | 176 | |
5443e65e | 177 | tbDrift = TMath::Min(tbDrift,maxDrift); |
178 | tbAmp = TMath::Min(tbAmp,maxAmp); | |
46d29e70 | 179 | |
5443e65e | 180 | fTimeBinsPerPlane = tbAmp + tbDrift; |
181 | fMaxGap = (Int_t) (fTimeBinsPerPlane * fGeom->Nplan() * fSkipDepth); | |
46d29e70 | 182 | |
5443e65e | 183 | fVocal = kFALSE; |
0a29d0f1 | 184 | |
5443e65e | 185 | } |
46d29e70 | 186 | |
5443e65e | 187 | //___________________________________________________________________ |
188 | AliTRDtracker::~AliTRDtracker() | |
46d29e70 | 189 | { |
5443e65e | 190 | delete fClusters; |
191 | delete fTracks; | |
192 | delete fSeeds; | |
193 | delete fGeom; | |
194 | delete fPar; | |
0a29d0f1 | 195 | |
5443e65e | 196 | for(Int_t geom_s = 0; geom_s < kTRACKING_SECTORS; geom_s++) { |
197 | delete fTrSec[geom_s]; | |
198 | } | |
199 | } | |
46d29e70 | 200 | |
201 | //_____________________________________________________________________ | |
5443e65e | 202 | inline Double_t f1trd(Double_t x1,Double_t y1, |
203 | Double_t x2,Double_t y2, | |
204 | Double_t x3,Double_t y3) | |
46d29e70 | 205 | { |
0a29d0f1 | 206 | // |
46d29e70 | 207 | // Initial approximation of the track curvature |
0a29d0f1 | 208 | // |
46d29e70 | 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 | //_____________________________________________________________________ | |
5443e65e | 221 | inline Double_t f2trd(Double_t x1,Double_t y1, |
222 | Double_t x2,Double_t y2, | |
223 | Double_t x3,Double_t y3) | |
46d29e70 | 224 | { |
0a29d0f1 | 225 | // |
5443e65e | 226 | // Initial approximation of the track curvature times X coordinate |
227 | // of the center of curvature | |
0a29d0f1 | 228 | // |
46d29e70 | 229 | |
230 | Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1); | |
231 | Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)- | |
232 | (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2)); | |
233 | Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)- | |
234 | (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1)); | |
235 | ||
236 | Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b); | |
237 | ||
238 | return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr); | |
239 | } | |
240 | ||
241 | //_____________________________________________________________________ | |
5443e65e | 242 | inline Double_t f3trd(Double_t x1,Double_t y1, |
243 | Double_t x2,Double_t y2, | |
244 | Double_t z1,Double_t z2) | |
46d29e70 | 245 | { |
0a29d0f1 | 246 | // |
46d29e70 | 247 | // Initial approximation of the tangent of the track dip angle |
0a29d0f1 | 248 | // |
46d29e70 | 249 | |
250 | return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); | |
251 | } | |
252 | ||
46d29e70 | 253 | //___________________________________________________________________ |
5443e65e | 254 | Int_t AliTRDtracker::Clusters2Tracks(const TFile *inp, TFile *out) |
46d29e70 | 255 | { |
0a29d0f1 | 256 | // |
5443e65e | 257 | // Finds tracks within the TRD. File <inp> is expected to contain seeds |
258 | // at the outer part of the TRD. If <inp> is NULL, the seeds | |
259 | // are found within the TRD if fAddTRDseeds is TRUE. | |
260 | // The tracks are propagated to the innermost time bin | |
261 | // of the TRD and stored in file <out>. | |
0a29d0f1 | 262 | // |
a819a5f7 | 263 | |
5443e65e | 264 | LoadEvent(); |
265 | ||
266 | TDirectory *savedir=gDirectory; | |
a819a5f7 | 267 | |
5443e65e | 268 | char tname[100]; |
a819a5f7 | 269 | |
5443e65e | 270 | if (!out->IsOpen()) { |
271 | cerr<<"AliTRDtracker::Clusters2Tracks(): output file is not open !\n"; | |
272 | return 1; | |
273 | } | |
46d29e70 | 274 | |
5443e65e | 275 | sprintf(tname,"seedTRDtoTPC_%d",fEvent); |
276 | TTree tpc_tree(tname,"Tree with seeds from TRD at outer TPC pad row"); | |
277 | AliTPCtrack *iotrack=0; | |
278 | tpc_tree.Branch("tracks","AliTPCtrack",&iotrack,32000,0); | |
279 | ||
280 | sprintf(tname,"TreeT%d_TRD",fEvent); | |
281 | TTree trd_tree(tname,"TRD tracks at inner TRD time bin"); | |
282 | AliTRDtrack *iotrack_trd=0; | |
283 | trd_tree.Branch("tracks","AliTRDtrack",&iotrack_trd,32000,0); | |
284 | ||
285 | Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins(); | |
286 | Float_t foundMin = fMinClustersInTrack * timeBins; | |
287 | ||
288 | if (inp) { | |
289 | TFile *in=(TFile*)inp; | |
290 | if (!in->IsOpen()) { | |
291 | cerr<<"AliTRDtracker::Clusters2Tracks(): file with seeds is not open !\n"; | |
292 | cerr<<" ... going for seeds finding inside the TRD\n"; | |
293 | } | |
294 | else { | |
295 | in->cd(); | |
296 | sprintf(tname,"TRDb_%d",fEvent); | |
297 | TTree *seedTree=(TTree*)in->Get(tname); | |
298 | if (!seedTree) { | |
299 | cerr<<"AliTRDtracker::Clusters2Tracks(): "; | |
300 | cerr<<"can't get a tree with track seeds !\n"; | |
301 | return 3; | |
302 | } | |
303 | AliTRDtrack *seed=new AliTRDtrack; | |
304 | seedTree->SetBranchAddress("tracks",&seed); | |
305 | ||
306 | Int_t n=(Int_t)seedTree->GetEntries(); | |
307 | for (Int_t i=0; i<n; i++) { | |
308 | seedTree->GetEvent(i); | |
309 | seed->ResetCovariance(); | |
310 | fSeeds->AddLast(new AliTRDtrack(*seed)); | |
311 | fNseeds++; | |
312 | } | |
313 | delete seed; | |
314 | } | |
315 | } | |
46d29e70 | 316 | |
5443e65e | 317 | out->cd(); |
46d29e70 | 318 | |
5443e65e | 319 | // find tracks from loaded seeds |
a819a5f7 | 320 | |
5443e65e | 321 | Int_t nseed=fSeeds->GetEntriesFast(); |
322 | Int_t i, found = 0; | |
323 | Int_t innerTB = fTrSec[0]->GetInnerTimeBin(); | |
324 | ||
325 | for (i=0; i<nseed; i++) { | |
326 | AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt; | |
327 | FollowProlongation(t, innerTB); | |
328 | if (t.GetNumberOfClusters() >= foundMin) { | |
329 | UseClusters(&t); | |
330 | CookLabel(pt, 1-fLabelFraction); | |
331 | // t.CookdEdx(); | |
332 | } | |
333 | iotrack_trd = pt; | |
334 | trd_tree.Fill(); | |
335 | cerr<<found++<<'\r'; | |
336 | ||
337 | if(PropagateToTPC(t)) { | |
338 | AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha()); | |
339 | iotrack = tpc; | |
340 | tpc_tree.Fill(); | |
341 | delete tpc; | |
342 | } | |
343 | delete fSeeds->RemoveAt(i); | |
344 | fNseeds--; | |
345 | } | |
a819a5f7 | 346 | |
5443e65e | 347 | cout<<"Number of loaded seeds: "<<nseed<<endl; |
348 | cout<<"Number of found tracks from loaded seeds: "<<found<<endl; | |
a819a5f7 | 349 | |
5443e65e | 350 | // after tracks from loaded seeds are found and the corresponding |
351 | // clusters are used, look for additional seeds from TRD | |
46d29e70 | 352 | |
5443e65e | 353 | if(fAddTRDseeds) { |
354 | // Find tracks for the seeds in the TRD | |
355 | Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins(); | |
356 | ||
357 | Int_t nSteps = (Int_t) (fSeedDepth / fSeedStep); | |
358 | Int_t gap = (Int_t) (timeBins * fSeedGap); | |
359 | Int_t step = (Int_t) (timeBins * fSeedStep); | |
360 | ||
361 | // make a first turn with tight cut on initial curvature | |
362 | for(Int_t turn = 1; turn <= 2; turn++) { | |
363 | if(turn == 2) { | |
364 | nSteps = (Int_t) (fSeedDepth / (3*fSeedStep)); | |
365 | step = (Int_t) (timeBins * (3*fSeedStep)); | |
366 | } | |
367 | for(Int_t i=0; i<nSteps; i++) { | |
368 | Int_t outer=timeBins-1-i*step; | |
369 | Int_t inner=outer-gap; | |
46d29e70 | 370 | |
5443e65e | 371 | nseed=fSeeds->GetEntriesFast(); |
372 | printf("\n initial number of seeds %d\n", nseed); | |
373 | ||
374 | MakeSeeds(inner, outer, turn); | |
375 | ||
376 | nseed=fSeeds->GetEntriesFast(); | |
377 | printf("\n number of seeds after MakeSeeds %d\n", nseed); | |
378 | ||
379 | ||
380 | for (Int_t i=0; i<nseed; i++) { | |
381 | AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt; | |
382 | FollowProlongation(t,innerTB); | |
383 | if (t.GetNumberOfClusters() >= foundMin) { | |
384 | UseClusters(&t); | |
385 | CookLabel(pt, 1-fLabelFraction); | |
386 | t.CookdEdx(); | |
387 | cerr<<found++<<'\r'; | |
388 | iotrack_trd = pt; | |
389 | trd_tree.Fill(); | |
390 | if(PropagateToTPC(t)) { | |
391 | AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha()); | |
392 | iotrack = tpc; | |
393 | tpc_tree.Fill(); | |
394 | delete tpc; | |
395 | } | |
396 | } | |
397 | delete fSeeds->RemoveAt(i); | |
398 | fNseeds--; | |
399 | } | |
46d29e70 | 400 | } |
5443e65e | 401 | } |
402 | } | |
403 | tpc_tree.Write(); | |
404 | trd_tree.Write(); | |
405 | ||
406 | cout<<"Total number of found tracks: "<<found<<endl; | |
407 | ||
408 | UnloadEvent(); | |
409 | ||
410 | savedir->cd(); | |
411 | ||
412 | return 0; | |
413 | } | |
414 | ||
415 | ||
416 | ||
417 | //_____________________________________________________________________________ | |
418 | Int_t AliTRDtracker::PropagateBack(const TFile *inp, TFile *out) { | |
419 | // | |
420 | // Reads seeds from file <inp>. The seeds are AliTPCtrack's found and | |
421 | // backpropagated by the TPC tracker. Each seed is first propagated | |
422 | // to the TRD, and then its prolongation is searched in the TRD. | |
423 | // If sufficiently long continuation of the track is found in the TRD | |
424 | // the track is updated, otherwise it's stored as originaly defined | |
425 | // by the TPC tracker. | |
426 | // | |
427 | ||
428 | ||
429 | LoadEvent(); | |
430 | ||
431 | TDirectory *savedir=gDirectory; | |
432 | ||
433 | TFile *in=(TFile*)inp; | |
434 | ||
435 | if (!in->IsOpen()) { | |
436 | cerr<<"AliTRDtracker::PropagateBack(): "; | |
437 | cerr<<"file with back propagated TPC tracks is not open !\n"; | |
438 | return 1; | |
439 | } | |
440 | ||
441 | if (!out->IsOpen()) { | |
442 | cerr<<"AliTRDtracker::PropagateBack(): "; | |
443 | cerr<<"file for back propagated TRD tracks is not open !\n"; | |
444 | return 2; | |
445 | } | |
446 | ||
447 | in->cd(); | |
448 | char tname[100]; | |
449 | sprintf(tname,"seedsTPCtoTRD_%d",fEvent); | |
450 | TTree *seedTree=(TTree*)in->Get(tname); | |
451 | if (!seedTree) { | |
452 | cerr<<"AliTRDtracker::PropagateBack(): "; | |
453 | cerr<<"can't get a tree with seeds from TPC !\n"; | |
454 | cerr<<"check if your version of TPC tracker creates tree "<<tname<<"\n"; | |
455 | return 3; | |
456 | } | |
46d29e70 | 457 | |
5443e65e | 458 | AliTPCtrack *seed=new AliTPCtrack; |
459 | seedTree->SetBranchAddress("tracks",&seed); | |
460 | ||
461 | Int_t n=(Int_t)seedTree->GetEntries(); | |
462 | for (Int_t i=0; i<n; i++) { | |
463 | seedTree->GetEvent(i); | |
464 | Int_t lbl = seed->GetLabel(); | |
465 | AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha()); | |
466 | tr->SetSeedLabel(lbl); | |
467 | fSeeds->AddLast(tr); | |
468 | fNseeds++; | |
469 | } | |
a819a5f7 | 470 | |
5443e65e | 471 | delete seed; |
472 | delete seedTree; | |
a819a5f7 | 473 | |
5443e65e | 474 | out->cd(); |
a819a5f7 | 475 | |
5443e65e | 476 | AliTPCtrack *otrack=0; |
a819a5f7 | 477 | |
5443e65e | 478 | sprintf(tname,"seedsTRDtoTOF1_%d",fEvent); |
479 | TTree tofTree1(tname,"Tracks back propagated through TPC and TRD"); | |
480 | tofTree1.Branch("tracks","AliTPCtrack",&otrack,32000,0); | |
a819a5f7 | 481 | |
5443e65e | 482 | sprintf(tname,"seedsTRDtoTOF2_%d",fEvent); |
483 | TTree tofTree2(tname,"Tracks back propagated through TPC and TRD"); | |
484 | tofTree2.Branch("tracks","AliTPCtrack",&otrack,32000,0); | |
46d29e70 | 485 | |
5443e65e | 486 | sprintf(tname,"seedsTRDtoPHOS_%d",fEvent); |
487 | TTree phosTree(tname,"Tracks back propagated through TPC and TRD"); | |
488 | phosTree.Branch("tracks","AliTPCtrack",&otrack,32000,0); | |
a819a5f7 | 489 | |
5443e65e | 490 | sprintf(tname,"seedsTRDtoRICH_%d",fEvent); |
491 | TTree richTree(tname,"Tracks back propagated through TPC and TRD"); | |
492 | richTree.Branch("tracks","AliTPCtrack",&otrack,32000,0); | |
a819a5f7 | 493 | |
5443e65e | 494 | sprintf(tname,"TRDb_%d",fEvent); |
495 | TTree trdTree(tname,"Back propagated TRD tracks at outer TRD time bin"); | |
496 | AliTRDtrack *otrack_trd=0; | |
497 | trdTree.Branch("tracks","AliTRDtrack",&otrack_trd,32000,0); | |
498 | ||
499 | Int_t found=0; | |
46d29e70 | 500 | |
5443e65e | 501 | Int_t nseed=fSeeds->GetEntriesFast(); |
46d29e70 | 502 | |
5443e65e | 503 | Float_t foundMin = fMinClustersInTrack * fTimeBinsPerPlane * fGeom->Nplan(); |
a819a5f7 | 504 | |
5443e65e | 505 | Int_t outermost_tb = fTrSec[0]->GetOuterTimeBin(); |
a819a5f7 | 506 | |
5443e65e | 507 | for (Int_t i=0; i<nseed; i++) { |
a819a5f7 | 508 | |
5443e65e | 509 | AliTRDtrack *ps=(AliTRDtrack*)fSeeds->UncheckedAt(i), &s=*ps; |
510 | Int_t expectedClr = FollowBackProlongation(s); | |
511 | Int_t foundClr = s.GetNumberOfClusters(); | |
512 | Int_t last_tb = fTrSec[0]->GetLayerNumber(s.GetX()); | |
a819a5f7 | 513 | |
5443e65e | 514 | printf("seed %d: found %d out of %d expected clusters, Min is %f\n", |
515 | i, foundClr, expectedClr, foundMin); | |
a819a5f7 | 516 | |
5443e65e | 517 | if (foundClr >= foundMin) { |
518 | s.CookdEdx(); | |
519 | CookLabel(ps, 1-fLabelFraction); | |
520 | UseClusters(ps); | |
521 | otrack_trd=ps; | |
522 | trdTree.Fill(); | |
523 | cout<<found++<<'\r'; | |
46d29e70 | 524 | } |
a819a5f7 | 525 | |
5443e65e | 526 | if(((expectedClr < 10) && (last_tb == outermost_tb)) || |
527 | ((expectedClr >= 10) && | |
528 | (((Float_t) foundClr) / ((Float_t) expectedClr) >= | |
529 | fMinFractionOfFoundClusters) && (last_tb == outermost_tb))) { | |
530 | ||
531 | Double_t x_tof = 375.5; | |
532 | ||
533 | if(PropagateToOuterPlane(s,x_tof)) { | |
534 | AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha()); | |
535 | otrack = pt; | |
536 | tofTree1.Fill(); | |
537 | delete pt; | |
538 | ||
539 | x_tof = 381.5; | |
540 | ||
541 | if(PropagateToOuterPlane(s,x_tof)) { | |
542 | AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha()); | |
543 | otrack = pt; | |
544 | tofTree2.Fill(); | |
545 | delete pt; | |
546 | ||
547 | Double_t x_phos = 460.; | |
548 | ||
549 | if(PropagateToOuterPlane(s,x_phos)) { | |
550 | AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha()); | |
551 | otrack = pt; | |
552 | phosTree.Fill(); | |
553 | delete pt; | |
554 | ||
555 | Double_t x_rich = 490+1.267; | |
556 | ||
557 | if(PropagateToOuterPlane(s,x_rich)) { | |
558 | AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha()); | |
559 | otrack = pt; | |
560 | richTree.Fill(); | |
561 | delete pt; | |
562 | } | |
563 | } | |
564 | } | |
565 | } | |
46d29e70 | 566 | } |
567 | } | |
5443e65e | 568 | |
569 | tofTree1.Write(); | |
570 | tofTree2.Write(); | |
571 | phosTree.Write(); | |
572 | richTree.Write(); | |
573 | trdTree.Write(); | |
46d29e70 | 574 | |
5443e65e | 575 | savedir->cd(); |
576 | cerr<<"Number of seeds: "<<nseed<<endl; | |
577 | cerr<<"Number of back propagated TRD tracks: "<<found<<endl; | |
46d29e70 | 578 | |
5443e65e | 579 | UnloadEvent(); |
46d29e70 | 580 | |
5443e65e | 581 | return 0; |
46d29e70 | 582 | |
5443e65e | 583 | } |
46d29e70 | 584 | |
bbf92647 | 585 | |
5443e65e | 586 | //--------------------------------------------------------------------------- |
587 | Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf) | |
588 | { | |
589 | // Starting from current position on track=t this function tries | |
590 | // to extrapolate the track up to timeBin=0 and to confirm prolongation | |
591 | // if a close cluster is found. Returns the number of clusters | |
592 | // expected to be found in sensitive layers | |
bbf92647 | 593 | |
5443e65e | 594 | Float_t wIndex, wTB, wChi2; |
595 | Float_t wYrt, wYclosest, wYcorrect, wYwindow; | |
596 | Float_t wZrt, wZclosest, wZcorrect, wZwindow; | |
597 | Float_t wPx, wPy, wPz, wC; | |
598 | Double_t Px, Py, Pz; | |
599 | Float_t wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2; | |
46d29e70 | 600 | |
5443e65e | 601 | Int_t trackIndex = t.GetLabel(); |
46d29e70 | 602 | |
5443e65e | 603 | Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5); |
46d29e70 | 604 | |
5443e65e | 605 | Int_t try_again=fMaxGap; |
46d29e70 | 606 | |
5443e65e | 607 | Double_t alpha=t.GetAlpha(); |
46d29e70 | 608 | |
5443e65e | 609 | if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); |
610 | if (alpha < 0. ) alpha += 2.*TMath::Pi(); | |
46d29e70 | 611 | |
5443e65e | 612 | Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect; |
46d29e70 | 613 | |
5443e65e | 614 | Double_t x0, rho, x, dx, y, ymax, z; |
46d29e70 | 615 | |
5443e65e | 616 | Int_t expectedNumberOfClusters = 0; |
617 | Bool_t lookForCluster; | |
46d29e70 | 618 | |
5443e65e | 619 | alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning |
46d29e70 | 620 | |
5443e65e | 621 | |
622 | for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr>rf; nr--) { | |
623 | ||
624 | y = t.GetY(); z = t.GetZ(); | |
625 | ||
626 | // first propagate to the inner surface of the current time bin | |
627 | fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster); | |
628 | x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ(); | |
629 | if(!t.PropagateTo(x,x0,rho,0.139)) break; | |
630 | y = t.GetY(); | |
631 | ymax = x*TMath::Tan(0.5*alpha); | |
632 | if (y > ymax) { | |
633 | s = (s+1) % ns; | |
634 | if (!t.Rotate(alpha)) break; | |
635 | if(!t.PropagateTo(x,x0,rho,0.139)) break; | |
636 | } else if (y <-ymax) { | |
637 | s = (s-1+ns) % ns; | |
638 | if (!t.Rotate(-alpha)) break; | |
639 | if(!t.PropagateTo(x,x0,rho,0.139)) break; | |
640 | } | |
46d29e70 | 641 | |
5443e65e | 642 | y = t.GetY(); z = t.GetZ(); |
643 | ||
644 | // now propagate to the middle plane of the next time bin | |
645 | fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster); | |
646 | x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ(); | |
647 | if(!t.PropagateTo(x,x0,rho,0.139)) break; | |
648 | y = t.GetY(); | |
649 | ymax = x*TMath::Tan(0.5*alpha); | |
650 | if (y > ymax) { | |
651 | s = (s+1) % ns; | |
652 | if (!t.Rotate(alpha)) break; | |
653 | if(!t.PropagateTo(x,x0,rho,0.139)) break; | |
654 | } else if (y <-ymax) { | |
655 | s = (s-1+ns) % ns; | |
656 | if (!t.Rotate(-alpha)) break; | |
657 | if(!t.PropagateTo(x,x0,rho,0.139)) break; | |
658 | } | |
46d29e70 | 659 | |
46d29e70 | 660 | |
5443e65e | 661 | if(lookForCluster) { |
a819a5f7 | 662 | |
46d29e70 | 663 | |
5443e65e | 664 | expectedNumberOfClusters++; |
46d29e70 | 665 | |
5443e65e | 666 | wIndex = (Float_t) t.GetLabel(); |
667 | wTB = nr; | |
46d29e70 | 668 | |
46d29e70 | 669 | |
5443e65e | 670 | AliTRDpropagationLayer& time_bin=*(fTrSec[s]->GetLayer(nr-1)); |
46d29e70 | 671 | |
46d29e70 | 672 | |
5443e65e | 673 | Double_t sy2=ExpectedSigmaY2(x,t.GetTgl(),t.GetPt()); |
46d29e70 | 674 | |
46d29e70 | 675 | |
5443e65e | 676 | Double_t sz2=ExpectedSigmaZ2(x,t.GetTgl()); |
46d29e70 | 677 | |
5443e65e | 678 | Double_t road; |
679 | if((t.GetSigmaY2() + sy2) > 0) road=10.*sqrt(t.GetSigmaY2() + sy2); | |
680 | else return expectedNumberOfClusters; | |
681 | ||
682 | wYrt = (Float_t) y; | |
683 | wZrt = (Float_t) z; | |
684 | wYwindow = (Float_t) road; | |
685 | t.GetPxPyPz(Px,Py,Pz); | |
686 | wPx = (Float_t) Px; | |
687 | wPy = (Float_t) Py; | |
688 | wPz = (Float_t) Pz; | |
689 | wC = (Float_t) t.GetC(); | |
690 | wSigmaC2 = (Float_t) t.GetSigmaC2(); | |
691 | wSigmaTgl2 = (Float_t) t.GetSigmaTgl2(); | |
692 | wSigmaY2 = (Float_t) t.GetSigmaY2(); | |
693 | wSigmaZ2 = (Float_t) t.GetSigmaZ2(); | |
694 | wChi2 = -1; | |
695 | ||
696 | if (road>fWideRoad) { | |
697 | if (t.GetNumberOfClusters()>4) | |
698 | cerr<<t.GetNumberOfClusters() | |
699 | <<"FindProlongation warning: Too broad road !\n"; | |
700 | return 0; | |
701 | } | |
702 | ||
703 | AliTRDcluster *cl=0; | |
704 | UInt_t index=0; | |
705 | ||
706 | Double_t max_chi2=fMaxChi2; | |
707 | ||
708 | wYclosest = 12345678; | |
709 | wYcorrect = 12345678; | |
710 | wZclosest = 12345678; | |
711 | wZcorrect = 12345678; | |
712 | wZwindow = TMath::Sqrt(2.25 * 12 * sz2); | |
713 | ||
714 | // Find the closest correct cluster for debugging purposes | |
715 | if (time_bin) { | |
716 | Float_t minDY = 1000000; | |
717 | for (Int_t i=0; i<time_bin; i++) { | |
718 | AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]); | |
719 | if((c->GetLabel(0) != trackIndex) && | |
720 | (c->GetLabel(1) != trackIndex) && | |
721 | (c->GetLabel(2) != trackIndex)) continue; | |
722 | if(TMath::Abs(c->GetY() - y) > minDY) continue; | |
723 | minDY = TMath::Abs(c->GetY() - y); | |
724 | wYcorrect = c->GetY(); | |
725 | wZcorrect = c->GetZ(); | |
726 | wChi2 = t.GetPredictedChi2(c); | |
727 | } | |
728 | } | |
46d29e70 | 729 | |
5443e65e | 730 | // Now go for the real cluster search |
a819a5f7 | 731 | |
5443e65e | 732 | if (time_bin) { |
a819a5f7 | 733 | |
5443e65e | 734 | for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) { |
735 | AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]); | |
736 | if (c->GetY() > y+road) break; | |
737 | if (c->IsUsed() > 0) continue; | |
738 | if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue; | |
739 | Double_t chi2=t.GetPredictedChi2(c); | |
740 | ||
741 | if (chi2 > max_chi2) continue; | |
742 | max_chi2=chi2; | |
743 | cl=c; | |
744 | index=time_bin.GetIndex(i); | |
745 | } | |
746 | ||
747 | if(!cl) { | |
748 | ||
749 | for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) { | |
750 | AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]); | |
751 | ||
752 | if (c->GetY() > y+road) break; | |
753 | if (c->IsUsed() > 0) continue; | |
754 | if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue; | |
755 | ||
756 | Double_t chi2=t.GetPredictedChi2(c); | |
757 | ||
758 | if (chi2 > max_chi2) continue; | |
759 | max_chi2=chi2; | |
760 | cl=c; | |
761 | index=time_bin.GetIndex(i); | |
762 | } | |
763 | } | |
764 | ||
a819a5f7 | 765 | |
5443e65e | 766 | if (cl) { |
767 | wYclosest = cl->GetY(); | |
768 | wZclosest = cl->GetZ(); | |
a819a5f7 | 769 | |
5443e65e | 770 | t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters()); |
771 | if(!t.Update(cl,max_chi2,index)) { | |
772 | if(!try_again--) return 0; | |
773 | } | |
774 | else try_again=fMaxGap; | |
775 | } | |
776 | else { | |
777 | if (try_again==0) break; | |
778 | try_again--; | |
779 | } | |
46d29e70 | 780 | |
5443e65e | 781 | /* |
782 | if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) { | |
783 | ||
784 | printf(" %f", wIndex); //1 | |
785 | printf(" %f", wTB); //2 | |
786 | printf(" %f", wYrt); //3 | |
787 | printf(" %f", wYclosest); //4 | |
788 | printf(" %f", wYcorrect); //5 | |
789 | printf(" %f", wYwindow); //6 | |
790 | printf(" %f", wZrt); //7 | |
791 | printf(" %f", wZclosest); //8 | |
792 | printf(" %f", wZcorrect); //9 | |
793 | printf(" %f", wZwindow); //10 | |
794 | printf(" %f", wPx); //11 | |
795 | printf(" %f", wPy); //12 | |
796 | printf(" %f", wPz); //13 | |
797 | printf(" %f", wSigmaC2*1000000); //14 | |
798 | printf(" %f", wSigmaTgl2*1000); //15 | |
799 | printf(" %f", wSigmaY2); //16 | |
800 | // printf(" %f", wSigmaZ2); //17 | |
801 | printf(" %f", wChi2); //17 | |
802 | printf(" %f", wC); //18 | |
803 | printf("\n"); | |
804 | } | |
805 | */ | |
806 | } | |
807 | } | |
808 | } | |
809 | return expectedNumberOfClusters; | |
810 | ||
811 | ||
812 | } | |
a819a5f7 | 813 | |
5443e65e | 814 | //___________________________________________________________________ |
46d29e70 | 815 | |
5443e65e | 816 | Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t) |
817 | { | |
818 | // Starting from current radial position of track <t> this function | |
819 | // extrapolates the track up to outer timebin and in the sensitive | |
820 | // layers confirms prolongation if a close cluster is found. | |
821 | // Returns the number of clusters expected to be found in sensitive layers | |
46d29e70 | 822 | |
5443e65e | 823 | Float_t wIndex, wTB, wChi2; |
824 | Float_t wYrt, wYclosest, wYcorrect, wYwindow; | |
825 | Float_t wZrt, wZclosest, wZcorrect, wZwindow; | |
826 | Float_t wPx, wPy, wPz, wC; | |
827 | Double_t Px, Py, Pz; | |
828 | Float_t wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2; | |
46d29e70 | 829 | |
5443e65e | 830 | Int_t trackIndex = t.GetLabel(); |
46d29e70 | 831 | |
5443e65e | 832 | Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5); |
a819a5f7 | 833 | |
5443e65e | 834 | Int_t try_again=fMaxGap; |
46d29e70 | 835 | |
5443e65e | 836 | Double_t alpha=t.GetAlpha(); |
46d29e70 | 837 | |
5443e65e | 838 | if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); |
839 | if (alpha < 0. ) alpha += 2.*TMath::Pi(); | |
46d29e70 | 840 | |
5443e65e | 841 | Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect; |
46d29e70 | 842 | |
5443e65e | 843 | Int_t outerTB = fTrSec[0]->GetOuterTimeBin(); |
46d29e70 | 844 | |
5443e65e | 845 | Double_t x0, rho, x, dx, y, ymax, z; |
a819a5f7 | 846 | |
5443e65e | 847 | Bool_t lookForCluster; |
a819a5f7 | 848 | |
5443e65e | 849 | Int_t expectedNumberOfClusters = 0; |
850 | x = t.GetX(); | |
46d29e70 | 851 | |
5443e65e | 852 | alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning |
46d29e70 | 853 | |
46d29e70 | 854 | |
5443e65e | 855 | for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr<outerTB; nr++) { |
46d29e70 | 856 | |
5443e65e | 857 | y = t.GetY(); z = t.GetZ(); |
46d29e70 | 858 | |
5443e65e | 859 | // first propagate to the outer surface of the current time bin |
46d29e70 | 860 | |
5443e65e | 861 | fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster); |
862 | x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ(); | |
46d29e70 | 863 | |
5443e65e | 864 | if(!t.PropagateTo(x,x0,rho,0.139)) break; |
46d29e70 | 865 | |
5443e65e | 866 | y = t.GetY(); |
867 | ymax = x*TMath::Tan(0.5*alpha); | |
a819a5f7 | 868 | |
5443e65e | 869 | if (y > ymax) { |
870 | s = (s+1) % ns; | |
871 | if (!t.Rotate(alpha)) break; | |
872 | if(!t.PropagateTo(x,x0,rho,0.139)) break; | |
873 | } else if (y <-ymax) { | |
874 | s = (s-1+ns) % ns; | |
875 | if (!t.Rotate(-alpha)) break; | |
876 | if(!t.PropagateTo(x,x0,rho,0.139)) break; | |
877 | } | |
878 | y = t.GetY(); z = t.GetZ(); | |
46d29e70 | 879 | |
5443e65e | 880 | // now propagate to the middle plane of the next time bin |
881 | fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster); | |
46d29e70 | 882 | |
5443e65e | 883 | x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ(); |
46d29e70 | 884 | |
5443e65e | 885 | if(!t.PropagateTo(x,x0,rho,0.139)) break; |
46d29e70 | 886 | |
5443e65e | 887 | y = t.GetY(); |
a819a5f7 | 888 | |
5443e65e | 889 | ymax = x*TMath::Tan(0.5*alpha); |
a819a5f7 | 890 | |
5443e65e | 891 | if(fVocal) printf("nr+1=%d, x %f, z %f, y %f, ymax %f\n",nr+1,x,z,y,ymax); |
892 | ||
893 | if (y > ymax) { | |
894 | s = (s+1) % ns; | |
895 | if (!t.Rotate(alpha)) break; | |
896 | if(!t.PropagateTo(x,x0,rho,0.139)) break; | |
897 | } else if (y <-ymax) { | |
898 | s = (s-1+ns) % ns; | |
899 | if (!t.Rotate(-alpha)) break; | |
900 | if(!t.PropagateTo(x,x0,rho,0.139)) break; | |
901 | } | |
46d29e70 | 902 | |
5443e65e | 903 | // printf("label %d, pl %d, lookForCluster %d \n", |
904 | // trackIndex, nr+1, lookForCluster); | |
905 | ||
906 | if(lookForCluster) { | |
907 | expectedNumberOfClusters++; | |
908 | ||
909 | wIndex = (Float_t) t.GetLabel(); | |
910 | wTB = fTrSec[s]->GetLayer(nr+1)->GetTimeBinIndex(); | |
911 | ||
912 | AliTRDpropagationLayer& time_bin=*(fTrSec[s]->GetLayer(nr+1)); | |
913 | Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt()); | |
914 | Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl()); | |
915 | if((t.GetSigmaY2() + sy2) < 0) break; | |
916 | Double_t road = 10.*sqrt(t.GetSigmaY2() + sy2); | |
917 | Double_t y=t.GetY(), z=t.GetZ(); | |
918 | ||
919 | wYrt = (Float_t) y; | |
920 | wZrt = (Float_t) z; | |
921 | wYwindow = (Float_t) road; | |
922 | t.GetPxPyPz(Px,Py,Pz); | |
923 | wPx = (Float_t) Px; | |
924 | wPy = (Float_t) Py; | |
925 | wPz = (Float_t) Pz; | |
926 | wC = (Float_t) t.GetC(); | |
927 | wSigmaC2 = (Float_t) t.GetSigmaC2(); | |
928 | wSigmaTgl2 = (Float_t) t.GetSigmaTgl2(); | |
929 | wSigmaY2 = (Float_t) t.GetSigmaY2(); | |
930 | wSigmaZ2 = (Float_t) t.GetSigmaZ2(); | |
931 | wChi2 = -1; | |
932 | ||
933 | if (road>fWideRoad) { | |
934 | if (t.GetNumberOfClusters()>4) | |
935 | cerr<<t.GetNumberOfClusters() | |
936 | <<"FindProlongation warning: Too broad road !\n"; | |
937 | return 0; | |
938 | } | |
939 | ||
940 | AliTRDcluster *cl=0; | |
941 | UInt_t index=0; | |
942 | ||
943 | Double_t max_chi2=fMaxChi2; | |
944 | ||
945 | wYclosest = 12345678; | |
946 | wYcorrect = 12345678; | |
947 | wZclosest = 12345678; | |
948 | wZcorrect = 12345678; | |
949 | wZwindow = TMath::Sqrt(2.25 * 12 * sz2); | |
950 | ||
951 | // Find the closest correct cluster for debugging purposes | |
952 | if (time_bin) { | |
953 | Float_t minDY = 1000000; | |
954 | for (Int_t i=0; i<time_bin; i++) { | |
955 | AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]); | |
956 | if((c->GetLabel(0) != trackIndex) && | |
957 | (c->GetLabel(1) != trackIndex) && | |
958 | (c->GetLabel(2) != trackIndex)) continue; | |
959 | if(TMath::Abs(c->GetY() - y) > minDY) continue; | |
960 | minDY = TMath::Abs(c->GetY() - y); | |
961 | wYcorrect = c->GetY(); | |
962 | wZcorrect = c->GetZ(); | |
963 | wChi2 = t.GetPredictedChi2(c); | |
46d29e70 | 964 | } |
5443e65e | 965 | } |
46d29e70 | 966 | |
5443e65e | 967 | // Now go for the real cluster search |
46d29e70 | 968 | |
5443e65e | 969 | if (time_bin) { |
970 | ||
971 | for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) { | |
972 | AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]); | |
973 | if (c->GetY() > y+road) break; | |
974 | if (c->IsUsed() > 0) continue; | |
975 | if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue; | |
976 | Double_t chi2=t.GetPredictedChi2(c); | |
977 | ||
978 | if (chi2 > max_chi2) continue; | |
979 | max_chi2=chi2; | |
980 | cl=c; | |
981 | index=time_bin.GetIndex(i); | |
982 | } | |
983 | ||
984 | if(!cl) { | |
985 | ||
986 | for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) { | |
987 | AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]); | |
988 | ||
989 | if (c->GetY() > y+road) break; | |
990 | if (c->IsUsed() > 0) continue; | |
991 | if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue; | |
992 | ||
993 | Double_t chi2=t.GetPredictedChi2(c); | |
994 | ||
995 | if (chi2 > max_chi2) continue; | |
996 | max_chi2=chi2; | |
997 | cl=c; | |
998 | index=time_bin.GetIndex(i); | |
999 | } | |
1000 | } | |
1001 | ||
1002 | if (cl) { | |
1003 | wYclosest = cl->GetY(); | |
1004 | wZclosest = cl->GetZ(); | |
1005 | ||
1006 | t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters()); | |
1007 | if(!t.Update(cl,max_chi2,index)) { | |
1008 | if(!try_again--) return 0; | |
1009 | } | |
1010 | else try_again=fMaxGap; | |
1011 | } | |
1012 | else { | |
1013 | if (try_again==0) break; | |
1014 | try_again--; | |
1015 | } | |
1016 | ||
1017 | /* | |
1018 | if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) { | |
1019 | ||
1020 | printf(" %f", wIndex); //1 | |
1021 | printf(" %f", wTB); //2 | |
1022 | printf(" %f", wYrt); //3 | |
1023 | printf(" %f", wYclosest); //4 | |
1024 | printf(" %f", wYcorrect); //5 | |
1025 | printf(" %f", wYwindow); //6 | |
1026 | printf(" %f", wZrt); //7 | |
1027 | printf(" %f", wZclosest); //8 | |
1028 | printf(" %f", wZcorrect); //9 | |
1029 | printf(" %f", wZwindow); //10 | |
1030 | printf(" %f", wPx); //11 | |
1031 | printf(" %f", wPy); //12 | |
1032 | printf(" %f", wPz); //13 | |
1033 | printf(" %f", wSigmaC2*1000000); //14 | |
1034 | printf(" %f", wSigmaTgl2*1000); //15 | |
1035 | printf(" %f", wSigmaY2); //16 | |
1036 | // printf(" %f", wSigmaZ2); //17 | |
1037 | printf(" %f", wChi2); //17 | |
1038 | printf(" %f", wC); //18 | |
1039 | printf("\n"); | |
1040 | } | |
1041 | */ | |
1042 | } | |
1043 | } | |
1044 | } | |
1045 | return expectedNumberOfClusters; | |
1046 | } | |
1047 | ||
1048 | //___________________________________________________________________ | |
1049 | ||
1050 | Int_t AliTRDtracker::PropagateToOuterPlane(AliTRDtrack& t, Double_t xToGo) | |
1051 | { | |
1052 | // Starting from current radial position of track <t> this function | |
1053 | // extrapolates the track up to radial position <xToGo>. | |
1054 | // Returns 1 if track reaches the plane, and 0 otherwise | |
1055 | ||
1056 | Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5); | |
1057 | ||
1058 | Double_t alpha=t.GetAlpha(); | |
1059 | ||
1060 | if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); | |
1061 | if (alpha < 0. ) alpha += 2.*TMath::Pi(); | |
1062 | ||
1063 | Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect; | |
1064 | ||
1065 | Bool_t lookForCluster; | |
1066 | Double_t x0, rho, x, dx, y, ymax, z; | |
1067 | ||
1068 | x = t.GetX(); | |
1069 | ||
1070 | alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning | |
1071 | ||
1072 | Int_t plToGo = fTrSec[0]->GetLayerNumber(xToGo); | |
1073 | ||
1074 | for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr<plToGo; nr++) { | |
1075 | ||
1076 | y = t.GetY(); z = t.GetZ(); | |
1077 | ||
1078 | // first propagate to the outer surface of the current time bin | |
1079 | fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster); | |
1080 | x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ(); | |
1081 | if(!t.PropagateTo(x,x0,rho,0.139)) return 0; | |
1082 | y = t.GetY(); | |
1083 | ymax = x*TMath::Tan(0.5*alpha); | |
1084 | if (y > ymax) { | |
1085 | s = (s+1) % ns; | |
1086 | if (!t.Rotate(alpha)) return 0; | |
1087 | } else if (y <-ymax) { | |
1088 | s = (s-1+ns) % ns; | |
1089 | if (!t.Rotate(-alpha)) return 0; | |
1090 | } | |
1091 | if(!t.PropagateTo(x,x0,rho,0.139)) return 0; | |
1092 | ||
1093 | y = t.GetY(); z = t.GetZ(); | |
1094 | ||
1095 | // now propagate to the middle plane of the next time bin | |
1096 | fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster); | |
1097 | x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ(); | |
1098 | if(!t.PropagateTo(x,x0,rho,0.139)) return 0; | |
1099 | y = t.GetY(); | |
1100 | ymax = x*TMath::Tan(0.5*alpha); | |
1101 | if (y > ymax) { | |
1102 | s = (s+1) % ns; | |
1103 | if (!t.Rotate(alpha)) return 0; | |
1104 | } else if (y <-ymax) { | |
1105 | s = (s-1+ns) % ns; | |
1106 | if (!t.Rotate(-alpha)) return 0; | |
1107 | } | |
1108 | if(!t.PropagateTo(x,x0,rho,0.139)) return 0; | |
1109 | } | |
1110 | return 1; | |
1111 | } | |
1112 | ||
1113 | //___________________________________________________________________ | |
1114 | ||
1115 | Int_t AliTRDtracker::PropagateToTPC(AliTRDtrack& t) | |
1116 | { | |
1117 | // Starting from current radial position of track <t> this function | |
1118 | // extrapolates the track up to radial position of the outermost | |
1119 | // padrow of the TPC. | |
1120 | // Returns 1 if track reaches the TPC, and 0 otherwise | |
1121 | ||
1122 | Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5); | |
1123 | ||
1124 | Double_t alpha=t.GetAlpha(); | |
1125 | ||
1126 | if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); | |
1127 | if (alpha < 0. ) alpha += 2.*TMath::Pi(); | |
1128 | ||
1129 | Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect; | |
1130 | ||
1131 | Bool_t lookForCluster; | |
1132 | Double_t x0, rho, x, dx, y, ymax, z; | |
1133 | ||
1134 | x = t.GetX(); | |
1135 | ||
1136 | alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning | |
1137 | ||
1138 | Int_t plTPC = fTrSec[0]->GetLayerNumber(246.055); | |
1139 | ||
1140 | for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr>plTPC; nr--) { | |
1141 | ||
1142 | y = t.GetY(); z = t.GetZ(); | |
1143 | ||
1144 | // first propagate to the outer surface of the current time bin | |
1145 | fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster); | |
1146 | x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ(); | |
1147 | if(!t.PropagateTo(x,x0,rho,0.139)) return 0; | |
1148 | y = t.GetY(); | |
1149 | ymax = x*TMath::Tan(0.5*alpha); | |
1150 | if (y > ymax) { | |
1151 | s = (s+1) % ns; | |
1152 | if (!t.Rotate(alpha)) return 0; | |
1153 | } else if (y <-ymax) { | |
1154 | s = (s-1+ns) % ns; | |
1155 | if (!t.Rotate(-alpha)) return 0; | |
1156 | } | |
1157 | if(!t.PropagateTo(x,x0,rho,0.139)) return 0; | |
1158 | ||
1159 | y = t.GetY(); z = t.GetZ(); | |
1160 | ||
1161 | // now propagate to the middle plane of the next time bin | |
1162 | fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster); | |
1163 | x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ(); | |
1164 | if(!t.PropagateTo(x,x0,rho,0.139)) return 0; | |
1165 | y = t.GetY(); | |
1166 | ymax = x*TMath::Tan(0.5*alpha); | |
1167 | if (y > ymax) { | |
1168 | s = (s+1) % ns; | |
1169 | if (!t.Rotate(alpha)) return 0; | |
1170 | } else if (y <-ymax) { | |
1171 | s = (s-1+ns) % ns; | |
1172 | if (!t.Rotate(-alpha)) return 0; | |
1173 | } | |
1174 | if(!t.PropagateTo(x,x0,rho,0.139)) return 0; | |
1175 | } | |
1176 | return 1; | |
1177 | } | |
1178 | ||
1179 | ||
1180 | //_____________________________________________________________________________ | |
1181 | void AliTRDtracker::LoadEvent() | |
1182 | { | |
1183 | // Fills clusters into TRD tracking_sectors | |
1184 | // Note that the numbering scheme for the TRD tracking_sectors | |
1185 | // differs from that of TRD sectors | |
1186 | ||
1187 | ReadClusters(fClusters); | |
1188 | Int_t ncl=fClusters->GetEntriesFast(); | |
1189 | cout<<"\n LoadSectors: sorting "<<ncl<<" clusters"<<endl; | |
1190 | ||
1191 | UInt_t index; | |
1192 | while (ncl--) { | |
1193 | printf("\r %d left ",ncl); | |
1194 | AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl); | |
1195 | Int_t detector=c->GetDetector(), local_time_bin=c->GetLocalTimeBin(); | |
1196 | Int_t sector=fGeom->GetSector(detector); | |
1197 | Int_t plane=fGeom->GetPlane(detector); | |
1198 | ||
1199 | Int_t tracking_sector = CookSectorIndex(sector); | |
1200 | ||
1201 | Int_t gtb = fTrSec[tracking_sector]->CookTimeBinIndex(plane,local_time_bin); | |
1202 | Int_t layer = fTrSec[tracking_sector]->GetLayerNumber(gtb); | |
1203 | ||
1204 | index=ncl; | |
1205 | fTrSec[tracking_sector]->GetLayer(layer)->InsertCluster(c,index); | |
1206 | } | |
1207 | printf("\r\n"); | |
1208 | ||
1209 | } | |
1210 | ||
1211 | //_____________________________________________________________________________ | |
1212 | void AliTRDtracker::UnloadEvent() | |
1213 | { | |
1214 | // | |
1215 | // Clears the arrays of clusters and tracks. Resets sectors and timebins | |
1216 | // | |
1217 | ||
1218 | Int_t i, nentr; | |
1219 | ||
1220 | nentr = fClusters->GetEntriesFast(); | |
1221 | for (i = 0; i < nentr; i++) delete fClusters->RemoveAt(i); | |
1222 | ||
1223 | nentr = fSeeds->GetEntriesFast(); | |
1224 | for (i = 0; i < nentr; i++) delete fSeeds->RemoveAt(i); | |
1225 | ||
1226 | nentr = fTracks->GetEntriesFast(); | |
1227 | for (i = 0; i < nentr; i++) delete fTracks->RemoveAt(i); | |
1228 | ||
1229 | Int_t nsec = AliTRDgeometry::kNsect; | |
1230 | ||
1231 | for (i = 0; i < nsec; i++) { | |
1232 | for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) { | |
1233 | fTrSec[i]->GetLayer(pl)->Clear(); | |
1234 | } | |
1235 | } | |
1236 | ||
1237 | } | |
1238 | ||
1239 | //__________________________________________________________________________ | |
1240 | void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, Int_t turn) | |
1241 | { | |
1242 | // Creates track seeds using clusters in timeBins=i1,i2 | |
1243 | ||
1244 | if(turn > 2) { | |
1245 | cerr<<"MakeSeeds: turn "<<turn<<" exceeds the limit of 2"<<endl; | |
1246 | return; | |
1247 | } | |
1248 | ||
1249 | Double_t x[5], c[15]; | |
1250 | Int_t max_sec=AliTRDgeometry::kNsect; | |
1251 | ||
1252 | Double_t alpha=AliTRDgeometry::GetAlpha(); | |
1253 | Double_t shift=AliTRDgeometry::GetAlpha()/2.; | |
1254 | Double_t cs=cos(alpha), sn=sin(alpha); | |
1255 | Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha); | |
1256 | ||
1257 | ||
1258 | Int_t i2 = fTrSec[0]->GetLayerNumber(inner); | |
1259 | Int_t i1 = fTrSec[0]->GetLayerNumber(outer); | |
1260 | ||
1261 | Double_t x1 =fTrSec[0]->GetX(i1); | |
1262 | Double_t xx2=fTrSec[0]->GetX(i2); | |
1263 | ||
1264 | for (Int_t ns=0; ns<max_sec; ns++) { | |
1265 | ||
1266 | Int_t nl2 = *(fTrSec[(ns-2+max_sec)%max_sec]->GetLayer(i2)); | |
1267 | Int_t nl=(*fTrSec[(ns-1+max_sec)%max_sec]->GetLayer(i2)); | |
1268 | Int_t nm=(*fTrSec[ns]->GetLayer(i2)); | |
1269 | Int_t nu=(*fTrSec[(ns+1)%max_sec]->GetLayer(i2)); | |
1270 | Int_t nu2=(*fTrSec[(ns+2)%max_sec]->GetLayer(i2)); | |
1271 | ||
1272 | AliTRDpropagationLayer& r1=*(fTrSec[ns]->GetLayer(i1)); | |
1273 | ||
1274 | for (Int_t is=0; is < r1; is++) { | |
1275 | Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ(); | |
1276 | ||
1277 | for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) { | |
1278 | ||
1279 | const AliTRDcluster *cl; | |
1280 | Double_t x2, y2, z2; | |
1281 | Double_t x3=0., y3=0.; | |
1282 | ||
1283 | if (js<nl2) { | |
1284 | if(turn != 2) continue; | |
1285 | AliTRDpropagationLayer& r2=*(fTrSec[(ns-2+max_sec)%max_sec]->GetLayer(i2)); | |
1286 | cl=r2[js]; | |
1287 | y2=cl->GetY(); z2=cl->GetZ(); | |
1288 | ||
1289 | x2= xx2*cs2+y2*sn2; | |
1290 | y2=-xx2*sn2+y2*cs2; | |
1291 | } | |
1292 | else if (js<nl2+nl) { | |
1293 | if(turn != 1) continue; | |
1294 | AliTRDpropagationLayer& r2=*(fTrSec[(ns-1+max_sec)%max_sec]->GetLayer(i2)); | |
1295 | cl=r2[js-nl2]; | |
1296 | y2=cl->GetY(); z2=cl->GetZ(); | |
1297 | ||
1298 | x2= xx2*cs+y2*sn; | |
1299 | y2=-xx2*sn+y2*cs; | |
1300 | } | |
1301 | else if (js<nl2+nl+nm) { | |
1302 | if(turn != 1) continue; | |
1303 | AliTRDpropagationLayer& r2=*(fTrSec[ns]->GetLayer(i2)); | |
1304 | cl=r2[js-nl2-nl]; | |
1305 | x2=xx2; y2=cl->GetY(); z2=cl->GetZ(); | |
1306 | } | |
1307 | else if (js<nl2+nl+nm+nu) { | |
1308 | if(turn != 1) continue; | |
1309 | AliTRDpropagationLayer& r2=*(fTrSec[(ns+1)%max_sec]->GetLayer(i2)); | |
1310 | cl=r2[js-nl2-nl-nm]; | |
1311 | y2=cl->GetY(); z2=cl->GetZ(); | |
1312 | ||
1313 | x2=xx2*cs-y2*sn; | |
1314 | y2=xx2*sn+y2*cs; | |
1315 | } | |
1316 | else { | |
1317 | if(turn != 2) continue; | |
1318 | AliTRDpropagationLayer& r2=*(fTrSec[(ns+2)%max_sec]->GetLayer(i2)); | |
1319 | cl=r2[js-nl2-nl-nm-nu]; | |
1320 | y2=cl->GetY(); z2=cl->GetZ(); | |
1321 | ||
1322 | x2=xx2*cs2-y2*sn2; | |
1323 | y2=xx2*sn2+y2*cs2; | |
1324 | } | |
1325 | ||
1326 | if(TMath::Abs(z1-z2) > fMaxSeedDeltaZ12) continue; | |
1327 | ||
1328 | Double_t zz=z1 - z1/x1*(x1-x2); | |
1329 | ||
1330 | if (TMath::Abs(zz-z2)>fMaxSeedDeltaZ) continue; | |
1331 | ||
1332 | Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1); | |
1333 | if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;} | |
1334 | ||
1335 | x[0]=y1; | |
1336 | x[1]=z1; | |
1337 | x[2]=f1trd(x1,y1,x2,y2,x3,y3); | |
1338 | ||
1339 | if (TMath::Abs(x[2]) > fMaxSeedC) continue; | |
1340 | ||
1341 | x[3]=f2trd(x1,y1,x2,y2,x3,y3); | |
1342 | ||
1343 | if (TMath::Abs(x[2]*x1-x[3]) >= 0.99999) continue; | |
1344 | ||
1345 | x[4]=f3trd(x1,y1,x2,y2,z1,z2); | |
1346 | ||
1347 | if (TMath::Abs(x[4]) > fMaxSeedTan) continue; | |
1348 | ||
1349 | Double_t a=asin(x[3]); | |
1350 | Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3])); | |
1351 | ||
1352 | if (TMath::Abs(zv)>fMaxSeedVertexZ) continue; | |
1353 | ||
1354 | Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2(); | |
1355 | Double_t sy2=cl->GetSigmaY2(), sz2=cl->GetSigmaZ2(); | |
1356 | Double_t sy3=fSeedErrorSY3, sy=fSeedErrorSY, sz=fSeedErrorSZ; | |
1357 | ||
1358 | Double_t f20=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy; | |
1359 | Double_t f22=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy; | |
1360 | Double_t f24=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy; | |
1361 | Double_t f30=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[3])/sy; | |
1362 | Double_t f32=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[3])/sy; | |
1363 | Double_t f34=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[3])/sy; | |
1364 | Double_t f40=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[4])/sy; | |
1365 | Double_t f41=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[4])/sz; | |
1366 | Double_t f42=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[4])/sy; | |
1367 | Double_t f43=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[4])/sz; | |
1368 | ||
1369 | c[0]=sy1; | |
1370 | c[1]=0.; c[2]=sz1; | |
1371 | c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f24*sy3*f24; | |
1372 | c[6]=f30*sy1; c[7]=0.; c[8]=f30*sy1*f20+f32*sy2*f22+f34*sy3*f24; | |
1373 | c[9]=f30*sy1*f30+f32*sy2*f32+f34*sy3*f34; | |
1374 | c[10]=f40*sy1; c[11]=f41*sz1; c[12]=f40*sy1*f20+f42*sy2*f22; | |
1375 | c[13]=f40*sy1*f30+f42*sy2*f32; | |
1376 | c[14]=f40*sy1*f40+f41*sz1*f41+f42*sy2*f42+f43*sz2*f43; | |
1377 | ||
1378 | UInt_t index=r1.GetIndex(is); | |
1379 | ||
1380 | AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift); | |
1381 | ||
1382 | Int_t rc=FollowProlongation(*track, i2); | |
1383 | ||
1384 | if ((rc < 1) || | |
1385 | (track->GetNumberOfClusters() < | |
1386 | (outer-inner)*fMinClustersInSeed)) delete track; | |
1387 | else { | |
1388 | fSeeds->AddLast(track); fNseeds++; | |
1389 | cerr<<"\r found seed "<<fNseeds; | |
1390 | } | |
1391 | } | |
1392 | } | |
1393 | } | |
1394 | } | |
1395 | ||
1396 | //_____________________________________________________________________________ | |
1397 | void AliTRDtracker::ReadClusters(TObjArray *array, const TFile *inp) | |
1398 | { | |
1399 | // | |
a819a5f7 | 1400 | // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0) |
1401 | // from the file. The names of the cluster tree and branches | |
1402 | // should match the ones used in AliTRDclusterizer::WriteClusters() | |
1403 | // | |
46d29e70 | 1404 | |
a819a5f7 | 1405 | TDirectory *savedir=gDirectory; |
1406 | ||
5443e65e | 1407 | if (inp) { |
1408 | TFile *in=(TFile*)inp; | |
1409 | if (!in->IsOpen()) { | |
1410 | cerr<<"AliTRDtracker::ReadClusters(): input file is not open !\n"; | |
1411 | return; | |
1412 | } | |
1413 | else{ | |
1414 | in->cd(); | |
1415 | } | |
1416 | } | |
a819a5f7 | 1417 | |
abaf1f1d | 1418 | Char_t treeName[12]; |
1419 | sprintf(treeName,"TreeR%d_TRD",fEvent); | |
5443e65e | 1420 | TTree *ClusterTree = (TTree*) gDirectory->Get(treeName); |
1421 | ||
1422 | TObjArray *ClusterArray = new TObjArray(400); | |
1423 | ||
1424 | ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray); | |
1425 | ||
1426 | Int_t nEntries = (Int_t) ClusterTree->GetEntries(); | |
1427 | printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName()); | |
a819a5f7 | 1428 | |
a819a5f7 | 1429 | // Loop through all entries in the tree |
1430 | Int_t nbytes; | |
1431 | AliTRDcluster *c = 0; | |
5443e65e | 1432 | printf("\n"); |
a819a5f7 | 1433 | |
1434 | for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) { | |
1435 | ||
1436 | // Import the tree | |
5443e65e | 1437 | nbytes += ClusterTree->GetEvent(iEntry); |
1438 | ||
a819a5f7 | 1439 | // Get the number of points in the detector |
5443e65e | 1440 | Int_t nCluster = ClusterArray->GetEntriesFast(); |
1441 | printf("\r Read %d clusters from entry %d", nCluster, iEntry); | |
1442 | ||
a819a5f7 | 1443 | // Loop through all TRD digits |
1444 | for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { | |
5443e65e | 1445 | c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster); |
a819a5f7 | 1446 | AliTRDcluster *co = new AliTRDcluster(*c); |
1447 | co->SetSigmaY2(c->GetSigmaY2() * fSY2corr); | |
5443e65e | 1448 | Int_t ltb = co->GetLocalTimeBin(); |
1449 | if(ltb != 0) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr); | |
1450 | ||
a819a5f7 | 1451 | array->AddLast(co); |
5443e65e | 1452 | delete ClusterArray->RemoveAt(iCluster); |
a819a5f7 | 1453 | } |
1454 | } | |
1455 | ||
5443e65e | 1456 | delete ClusterArray; |
1457 | savedir->cd(); | |
1458 | ||
a819a5f7 | 1459 | } |
1460 | ||
5443e65e | 1461 | //______________________________________________________________________ |
1462 | void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename) | |
a819a5f7 | 1463 | { |
1464 | // | |
5443e65e | 1465 | // Reads AliTRDclusters from file <filename>. The names of the cluster |
1466 | // tree and branches should match the ones used in | |
1467 | // AliTRDclusterizer::WriteClusters() | |
1468 | // if <array> == 0, clusters are added into AliTRDtracker fCluster array | |
a819a5f7 | 1469 | // |
1470 | ||
5443e65e | 1471 | TDirectory *savedir=gDirectory; |
46d29e70 | 1472 | |
5443e65e | 1473 | TFile *file = TFile::Open(filename); |
1474 | if (!file->IsOpen()) { | |
1475 | cerr<<"Can't open file with TRD clusters"<<endl; | |
1476 | return; | |
1477 | } | |
46d29e70 | 1478 | |
5443e65e | 1479 | Char_t treeName[12]; |
1480 | sprintf(treeName,"TreeR%d_TRD",fEvent); | |
1481 | TTree *ClusterTree = (TTree*) gDirectory->Get(treeName); | |
46d29e70 | 1482 | |
5443e65e | 1483 | if (!ClusterTree) { |
1484 | cerr<<"AliTRDtracker::ReadClusters(): "; | |
1485 | cerr<<"can't get a tree with clusters !\n"; | |
1486 | return; | |
1487 | } | |
46d29e70 | 1488 | |
5443e65e | 1489 | TObjArray *ClusterArray = new TObjArray(400); |
46d29e70 | 1490 | |
5443e65e | 1491 | ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray); |
46d29e70 | 1492 | |
5443e65e | 1493 | Int_t nEntries = (Int_t) ClusterTree->GetEntries(); |
1494 | cout<<"found "<<nEntries<<" in ClusterTree"<<endl; | |
a819a5f7 | 1495 | |
5443e65e | 1496 | // Loop through all entries in the tree |
1497 | Int_t nbytes; | |
1498 | AliTRDcluster *c = 0; | |
bbf92647 | 1499 | |
5443e65e | 1500 | printf("\n"); |
bbf92647 | 1501 | |
5443e65e | 1502 | for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) { |
46d29e70 | 1503 | |
5443e65e | 1504 | // Import the tree |
1505 | nbytes += ClusterTree->GetEvent(iEntry); | |
0a29d0f1 | 1506 | |
5443e65e | 1507 | // Get the number of points in the detector |
1508 | Int_t nCluster = ClusterArray->GetEntriesFast(); | |
1509 | printf("\r Read %d clusters from entry %d", nCluster, iEntry); | |
1510 | ||
1511 | // Loop through all TRD digits | |
1512 | for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { | |
1513 | c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster); | |
1514 | AliTRDcluster *co = new AliTRDcluster(*c); | |
1515 | co->SetSigmaY2(c->GetSigmaY2() * fSY2corr); | |
1516 | Int_t ltb = co->GetLocalTimeBin(); | |
1517 | if(ltb != 0) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr); | |
1518 | array->AddLast(co); | |
1519 | delete ClusterArray->RemoveAt(iCluster); | |
1520 | } | |
46d29e70 | 1521 | } |
0a29d0f1 | 1522 | |
5443e65e | 1523 | file->Close(); |
1524 | delete ClusterArray; | |
1525 | savedir->cd(); | |
1526 | ||
1527 | } | |
1528 | ||
46d29e70 | 1529 | |
1530 | //__________________________________________________________________ | |
5443e65e | 1531 | void AliTRDtracker::CookLabel(AliKalmanTrack* pt, Float_t wrong) const { |
46d29e70 | 1532 | |
1533 | Int_t label=123456789, index, i, j; | |
5443e65e | 1534 | Int_t ncl=pt->GetNumberOfClusters(); |
1535 | const Int_t range = fTrSec[0]->GetOuterTimeBin()+1; | |
1536 | ||
1537 | Bool_t label_added; | |
46d29e70 | 1538 | |
5443e65e | 1539 | // Int_t s[range][2]; |
1540 | Int_t **s = new Int_t* [range]; | |
1541 | for (i=0; i<range; i++) { | |
d1b06c24 | 1542 | s[i] = new Int_t[2]; |
1543 | } | |
5443e65e | 1544 | for (i=0; i<range; i++) { |
46d29e70 | 1545 | s[i][0]=-1; |
1546 | s[i][1]=0; | |
1547 | } | |
1548 | ||
1549 | Int_t t0,t1,t2; | |
1550 | for (i=0; i<ncl; i++) { | |
5443e65e | 1551 | index=pt->GetClusterIndex(i); |
46d29e70 | 1552 | AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index); |
5443e65e | 1553 | t0=c->GetLabel(0); |
1554 | t1=c->GetLabel(1); | |
1555 | t2=c->GetLabel(2); | |
46d29e70 | 1556 | } |
1557 | ||
1558 | for (i=0; i<ncl; i++) { | |
5443e65e | 1559 | index=pt->GetClusterIndex(i); |
46d29e70 | 1560 | AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index); |
1561 | for (Int_t k=0; k<3; k++) { | |
5443e65e | 1562 | label=c->GetLabel(k); |
1563 | label_added=kFALSE; j=0; | |
46d29e70 | 1564 | if (label >= 0) { |
5443e65e | 1565 | while ( (!label_added) && ( j < range ) ) { |
46d29e70 | 1566 | if (s[j][0]==label || s[j][1]==0) { |
1567 | s[j][0]=label; | |
1568 | s[j][1]=s[j][1]+1; | |
5443e65e | 1569 | label_added=kTRUE; |
46d29e70 | 1570 | } |
1571 | j++; | |
1572 | } | |
1573 | } | |
1574 | } | |
1575 | } | |
1576 | ||
46d29e70 | 1577 | Int_t max=0; |
1578 | label = -123456789; | |
1579 | ||
5443e65e | 1580 | for (i=0; i<range; i++) { |
46d29e70 | 1581 | if (s[i][1]>max) { |
1582 | max=s[i][1]; label=s[i][0]; | |
1583 | } | |
1584 | } | |
5443e65e | 1585 | |
1586 | for (i=0; i<range; i++) { | |
1587 | delete []s[i]; | |
1588 | } | |
1589 | ||
d1b06c24 | 1590 | delete []s; |
5443e65e | 1591 | |
1592 | if ((1.- Float_t(max)/ncl) > wrong) label=-label; | |
1593 | ||
1594 | pt->SetLabel(label); | |
1595 | ||
46d29e70 | 1596 | } |
1597 | ||
5443e65e | 1598 | |
1599 | //__________________________________________________________________ | |
ca6c93d6 | 1600 | void AliTRDtracker::UseClusters(const AliKalmanTrack* t, Int_t from) const { |
5443e65e | 1601 | Int_t ncl=t->GetNumberOfClusters(); |
1602 | for (Int_t i=from; i<ncl; i++) { | |
1603 | Int_t index = t->GetClusterIndex(i); | |
1604 | AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index); | |
1605 | c->Use(); | |
1606 | } | |
1607 | } | |
1608 | ||
1609 | ||
1610 | //_____________________________________________________________________ | |
1611 | Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt) | |
1612 | { | |
1613 | // Parametrised "expected" error of the cluster reconstruction in Y | |
1614 | ||
1615 | Double_t s = 0.08 * 0.08; | |
1616 | return s; | |
1617 | } | |
1618 | ||
1619 | //_____________________________________________________________________ | |
1620 | Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl) | |
0a29d0f1 | 1621 | { |
5443e65e | 1622 | // Parametrised "expected" error of the cluster reconstruction in Z |
1623 | ||
1624 | Double_t s = 6 * 6 /12.; | |
1625 | return s; | |
1626 | } | |
1627 | ||
1628 | ||
1629 | //_____________________________________________________________________ | |
1630 | Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t local_tb) const | |
1631 | { | |
1632 | // | |
1633 | // Returns radial position which corresponds to time bin <local_tb> | |
1634 | // in tracking sector <sector> and plane <plane> | |
1635 | // | |
1636 | ||
1637 | Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, local_tb); | |
1638 | Int_t pl = fTrSec[sector]->GetLayerNumber(index); | |
1639 | return fTrSec[sector]->GetLayer(pl)->GetX(); | |
1640 | ||
1641 | } | |
1642 | ||
1643 | ||
1644 | //_______________________________________________________ | |
1645 | AliTRDtracker::AliTRDpropagationLayer::AliTRDpropagationLayer(Double_t x, | |
1646 | Double_t dx, Double_t rho, Double_t x0, Int_t tb_index) | |
1647 | { | |
0a29d0f1 | 1648 | // |
5443e65e | 1649 | // AliTRDpropagationLayer constructor |
0a29d0f1 | 1650 | // |
46d29e70 | 1651 | |
5443e65e | 1652 | fN = 0; fX = x; fdX = dx; fRho = rho; fX0 = x0; |
1653 | fClusters = NULL; fIndex = NULL; fTimeBinIndex = tb_index; | |
1654 | ||
46d29e70 | 1655 | |
5443e65e | 1656 | for(Int_t i=0; i < (Int_t) kZONES; i++) { |
1657 | fZc[i]=0; fZmax[i] = 0; | |
a819a5f7 | 1658 | } |
5443e65e | 1659 | |
1660 | fYmax = 0; | |
1661 | ||
1662 | if(fTimeBinIndex >= 0) { | |
ca6c93d6 | 1663 | fClusters = new AliTRDcluster*[kMAX_CLUSTER_PER_TIME_BIN]; |
1664 | fIndex = new UInt_t[kMAX_CLUSTER_PER_TIME_BIN]; | |
a819a5f7 | 1665 | } |
46d29e70 | 1666 | |
5443e65e | 1667 | fHole = kFALSE; |
1668 | fHoleZc = 0; | |
1669 | fHoleZmax = 0; | |
1670 | fHoleYc = 0; | |
1671 | fHoleYmax = 0; | |
1672 | fHoleRho = 0; | |
1673 | fHoleX0 = 0; | |
1674 | ||
1675 | } | |
1676 | ||
1677 | //_______________________________________________________ | |
1678 | void AliTRDtracker::AliTRDpropagationLayer::SetHole( | |
1679 | Double_t Zmax, Double_t Ymax, Double_t rho, | |
1680 | Double_t x0, Double_t Yc, Double_t Zc) | |
1681 | { | |
1682 | // | |
1683 | // Sets hole in the layer | |
1684 | // | |
1685 | ||
1686 | fHole = kTRUE; | |
1687 | fHoleZc = Zc; | |
1688 | fHoleZmax = Zmax; | |
1689 | fHoleYc = Yc; | |
1690 | fHoleYmax = Ymax; | |
1691 | fHoleRho = rho; | |
1692 | fHoleX0 = x0; | |
1693 | } | |
1694 | ||
46d29e70 | 1695 | |
5443e65e | 1696 | //_______________________________________________________ |
1697 | AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, Int_t gs, AliTRDparameter* par) | |
1698 | { | |
1699 | // | |
1700 | // AliTRDtrackingSector Constructor | |
1701 | // | |
1702 | ||
1703 | fGeom = geo; | |
1704 | fPar = par; | |
1705 | fGeomSector = gs; | |
1706 | fTzeroShift = 0.13; | |
1707 | fN = 0; | |
1708 | ||
1709 | for(UInt_t i=0; i < kMAX_TIME_BIN_INDEX; i++) fTimeBinIndex[i] = -1; | |
1710 | ||
1711 | ||
1712 | AliTRDpropagationLayer* ppl; | |
1713 | ||
1714 | Double_t x, xin, xout, dx, rho, x0; | |
1715 | Int_t steps; | |
46d29e70 | 1716 | |
5443e65e | 1717 | // set time bins in the gas of the TPC |
46d29e70 | 1718 | |
5443e65e | 1719 | xin = 246.055; xout = 254.055; steps = 20; dx = (xout-xin)/steps; |
1720 | rho = 0.9e-3; x0 = 28.94; | |
1721 | ||
1722 | for(Int_t i=0; i<steps; i++) { | |
1723 | x = xin + i*dx + dx/2; | |
1724 | ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1); | |
1725 | InsertLayer(ppl); | |
46d29e70 | 1726 | } |
1727 | ||
5443e65e | 1728 | // set time bins in the outer field cage vessel |
46d29e70 | 1729 | |
5443e65e | 1730 | dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; x0 = 44.77; // Tedlar |
1731 | ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1); | |
1732 | InsertLayer(ppl); | |
46d29e70 | 1733 | |
5443e65e | 1734 | dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; x0 = 44.86; // prepreg |
1735 | ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1); | |
1736 | InsertLayer(ppl); | |
1737 | ||
1738 | dx = 2.; xin = xout; xout = xin + dx; rho = 1.45*0.02; x0 = 41.28; // Nomex | |
1739 | steps = 5; dx = (xout - xin)/steps; | |
1740 | for(Int_t i=0; i<steps; i++) { | |
1741 | x = xin + i*dx + dx/2; | |
1742 | ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1); | |
1743 | InsertLayer(ppl); | |
1744 | } | |
1745 | ||
1746 | dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; x0 = 44.86; // prepreg | |
1747 | ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1); | |
1748 | InsertLayer(ppl); | |
1749 | ||
1750 | dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; x0 = 44.77; // Tedlar | |
1751 | ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1); | |
1752 | InsertLayer(ppl); | |
0a29d0f1 | 1753 | |
5443e65e | 1754 | |
1755 | // set time bins in CO2 | |
1756 | ||
1757 | xin = xout; xout = 275.0; | |
1758 | steps = 50; dx = (xout - xin)/steps; | |
1759 | rho = 1.977e-3; x0 = 36.2; | |
1760 | ||
1761 | for(Int_t i=0; i<steps; i++) { | |
1762 | x = xin + i*dx + dx/2; | |
1763 | ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1); | |
1764 | InsertLayer(ppl); | |
1765 | } | |
1766 | ||
1767 | // set time bins in the outer containment vessel | |
1768 | ||
1769 | dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; x0 = 24.01; // Al | |
1770 | ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1); | |
1771 | InsertLayer(ppl); | |
1772 | ||
1773 | dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; x0 = 44.77; // Tedlar | |
1774 | ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1); | |
1775 | InsertLayer(ppl); | |
1776 | ||
1777 | dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; x0 = 44.86; // prepreg | |
1778 | ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1); | |
1779 | InsertLayer(ppl); | |
1780 | ||
1781 | dx = 3.; xin = xout; xout = xin + dx; rho = 1.45*0.02; x0 = 41.28; // Nomex | |
1782 | steps = 10; dx = (xout - xin)/steps; | |
1783 | for(Int_t i=0; i<steps; i++) { | |
1784 | x = xin + i*dx + dx/2; | |
1785 | ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1); | |
1786 | InsertLayer(ppl); | |
1787 | } | |
1788 | ||
1789 | dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; x0 = 44.86; // prepreg | |
1790 | ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1); | |
1791 | InsertLayer(ppl); | |
1792 | ||
1793 | dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; x0 = 44.77; // Tedlar | |
1794 | ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1); | |
1795 | InsertLayer(ppl); | |
1796 | ||
1797 | dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; x0 = 24.01; // Al | |
1798 | ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1); | |
1799 | InsertLayer(ppl); | |
1800 | ||
1801 | Double_t xtrd = (Double_t) fGeom->Rmin(); | |
1802 | ||
1803 | // add layers between TPC and TRD (Air temporarily) | |
1804 | xin = xout; xout = xtrd; | |
1805 | steps = 50; dx = (xout - xin)/steps; | |
1806 | rho = 1.2e-3; x0 = 36.66; | |
1807 | ||
1808 | for(Int_t i=0; i<steps; i++) { | |
1809 | x = xin + i*dx + dx/2; | |
1810 | ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1); | |
1811 | InsertLayer(ppl); | |
1812 | } | |
1813 | ||
1814 | ||
1815 | Double_t alpha=AliTRDgeometry::GetAlpha(); | |
1816 | ||
1817 | // add layers for each of the planes | |
1818 | ||
1819 | Double_t dxRo = (Double_t) fGeom->CroHght(); // Rohacell | |
1820 | Double_t dxSpace = (Double_t) fGeom->Cspace(); // Spacing between planes | |
1821 | Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region | |
1822 | Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region | |
1823 | Double_t dxRad = (Double_t) fGeom->CraHght(); // Radiator | |
1824 | Double_t dxTEC = dxRad + dxDrift + dxAmp + dxRo; | |
1825 | Double_t dxPlane = dxTEC + dxSpace; | |
1826 | ||
1827 | Int_t tbBefore = (Int_t) (dxAmp/fPar->GetTimeBinSize()); | |
1828 | if(tbBefore > fPar->GetTimeBefore()) tbBefore = fPar->GetTimeBefore(); | |
1829 | ||
1830 | Int_t tb, tb_index; | |
1831 | const Int_t nChambers = AliTRDgeometry::Ncham(); | |
ca6c93d6 | 1832 | Double_t Ymax = 0, holeYmax = 0; |
1833 | Double_t * Zc = new Double_t[nChambers]; | |
1834 | Double_t * Zmax = new Double_t[nChambers]; | |
5443e65e | 1835 | Double_t holeZmax = 1000.; // the whole sector is missing |
1836 | ||
1837 | ||
1838 | for(Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) { | |
1839 | ||
1840 | // Radiator | |
1841 | xin = xtrd + plane * dxPlane; xout = xin + dxRad; | |
1842 | steps = 12; dx = (xout - xin)/steps; rho = 0.074; x0 = 40.6; | |
1843 | for(Int_t i=0; i<steps; i++) { | |
1844 | x = xin + i*dx + dx/2; | |
1845 | ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1); | |
1846 | if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) { | |
1847 | holeYmax = x*TMath::Tan(0.5*alpha); | |
1848 | ppl->SetHole(holeYmax, holeZmax); | |
1849 | } | |
1850 | if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) { | |
1851 | holeYmax = x*TMath::Tan(0.5*alpha); | |
1852 | ppl->SetHole(holeYmax, holeZmax); | |
1853 | } | |
1854 | InsertLayer(ppl); | |
1855 | } | |
1856 | ||
1857 | Ymax = fGeom->GetChamberWidth(plane)/2; | |
1858 | for(Int_t ch = 0; ch < nChambers; ch++) { | |
1859 | Zmax[ch] = fGeom->GetChamberLength(plane,ch)/2; | |
1860 | Float_t pad = fPar->GetRowPadSize(plane,ch,0); | |
1861 | Float_t row0 = fPar->GetRow0(plane,ch,0); | |
1862 | Int_t nPads = fPar->GetRowMax(plane,ch,0); | |
1863 | Zc[ch] = (pad * nPads)/2 + row0 - pad/2; | |
1864 | } | |
1865 | ||
1866 | dx = fPar->GetTimeBinSize(); | |
1867 | rho = 0.00295 * 0.85; x0 = 11.0; | |
1868 | ||
1869 | Double_t x0 = (Double_t) fPar->GetTime0(plane); | |
1870 | Double_t xbottom = x0 - dxDrift; | |
1871 | Double_t xtop = x0 + dxAmp; | |
1872 | ||
1873 | // Amplification region | |
1874 | ||
1875 | steps = (Int_t) (dxAmp/dx); | |
1876 | ||
1877 | for(tb = 0; tb < steps; tb++) { | |
1878 | x = x0 + tb * dx + dx/2; | |
1879 | tb_index = CookTimeBinIndex(plane, -tb-1); | |
1880 | ppl = new AliTRDpropagationLayer(x,dx,rho,x0,tb_index); | |
1881 | ppl->SetYmax(Ymax); | |
1882 | for(Int_t ch = 0; ch < nChambers; ch++) { | |
1883 | ppl->SetZmax(ch, Zc[ch], Zmax[ch]); | |
1884 | } | |
1885 | if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) { | |
1886 | holeYmax = x*TMath::Tan(0.5*alpha); | |
1887 | ppl->SetHole(holeYmax, holeZmax); | |
1888 | } | |
1889 | if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) { | |
1890 | holeYmax = x*TMath::Tan(0.5*alpha); | |
1891 | ppl->SetHole(holeYmax, holeZmax); | |
1892 | } | |
1893 | InsertLayer(ppl); | |
1894 | } | |
1895 | tb_index = CookTimeBinIndex(plane, -steps); | |
1896 | x = (x + dx/2 + xtop)/2; | |
1897 | dx = 2*(xtop-x); | |
1898 | ppl = new AliTRDpropagationLayer(x,dx,rho,x0,tb_index); | |
1899 | ppl->SetYmax(Ymax); | |
1900 | for(Int_t ch = 0; ch < nChambers; ch++) { | |
1901 | ppl->SetZmax(ch, Zc[ch], Zmax[ch]); | |
1902 | } | |
1903 | if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) { | |
1904 | holeYmax = x*TMath::Tan(0.5*alpha); | |
1905 | ppl->SetHole(holeYmax, holeZmax); | |
1906 | } | |
1907 | if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) { | |
1908 | holeYmax = x*TMath::Tan(0.5*alpha); | |
1909 | ppl->SetHole(holeYmax, holeZmax); | |
1910 | } | |
1911 | InsertLayer(ppl); | |
1912 | ||
1913 | // Drift region | |
1914 | dx = fPar->GetTimeBinSize(); | |
1915 | steps = (Int_t) (dxDrift/dx); | |
1916 | ||
1917 | for(tb = 0; tb < steps; tb++) { | |
1918 | x = x0 - tb * dx - dx/2; | |
1919 | tb_index = CookTimeBinIndex(plane, tb); | |
1920 | ||
1921 | ppl = new AliTRDpropagationLayer(x,dx,rho,x0,tb_index); | |
1922 | ppl->SetYmax(Ymax); | |
1923 | for(Int_t ch = 0; ch < nChambers; ch++) { | |
1924 | ppl->SetZmax(ch, Zc[ch], Zmax[ch]); | |
1925 | } | |
1926 | if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) { | |
1927 | holeYmax = x*TMath::Tan(0.5*alpha); | |
1928 | ppl->SetHole(holeYmax, holeZmax); | |
1929 | } | |
1930 | if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) { | |
1931 | holeYmax = x*TMath::Tan(0.5*alpha); | |
1932 | ppl->SetHole(holeYmax, holeZmax); | |
1933 | } | |
1934 | InsertLayer(ppl); | |
1935 | } | |
1936 | tb_index = CookTimeBinIndex(plane, steps); | |
1937 | x = (x - dx/2 + xbottom)/2; | |
1938 | dx = 2*(x-xbottom); | |
1939 | ppl = new AliTRDpropagationLayer(x,dx,rho,x0,tb_index); | |
1940 | ppl->SetYmax(Ymax); | |
1941 | for(Int_t ch = 0; ch < nChambers; ch++) { | |
1942 | ppl->SetZmax(ch, Zc[ch], Zmax[ch]); | |
1943 | } | |
1944 | if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) { | |
1945 | holeYmax = x*TMath::Tan(0.5*alpha); | |
1946 | ppl->SetHole(holeYmax, holeZmax); | |
1947 | } | |
1948 | if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) { | |
1949 | holeYmax = x*TMath::Tan(0.5*alpha); | |
1950 | ppl->SetHole(holeYmax, holeZmax); | |
1951 | } | |
1952 | InsertLayer(ppl); | |
1953 | ||
1954 | // Pad Plane | |
1955 | xin = xtop; dx = 0.025; xout = xin + dx; rho = 1.7; x0 = 33.0; | |
1956 | ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1); | |
1957 | if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) { | |
1958 | holeYmax = (xin+dx/2)*TMath::Tan(0.5*alpha); | |
1959 | ppl->SetHole(holeYmax, holeZmax); | |
1960 | } | |
1961 | if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) { | |
1962 | holeYmax = (xin+dx/2)*TMath::Tan(0.5*alpha); | |
1963 | ppl->SetHole(holeYmax, holeZmax); | |
1964 | } | |
1965 | InsertLayer(ppl); | |
1966 | ||
1967 | // Rohacell | |
1968 | xin = xout; xout = xtrd + (plane + 1) * dxPlane - dxSpace; | |
1969 | steps = 5; dx = (xout - xin)/steps; rho = 0.074; x0 = 40.6; | |
1970 | for(Int_t i=0; i<steps; i++) { | |
1971 | x = xin + i*dx + dx/2; | |
1972 | ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1); | |
1973 | if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) { | |
1974 | holeYmax = x*TMath::Tan(0.5*alpha); | |
1975 | ppl->SetHole(holeYmax, holeZmax); | |
1976 | } | |
1977 | if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) { | |
1978 | holeYmax = x*TMath::Tan(0.5*alpha); | |
1979 | ppl->SetHole(holeYmax, holeZmax); | |
1980 | } | |
1981 | InsertLayer(ppl); | |
1982 | } | |
1983 | ||
1984 | // Space between the chambers, air | |
1985 | xin = xout; xout = xtrd + (plane + 1) * dxPlane; | |
1986 | steps = 5; dx = (xout - xin)/steps; rho = 1.29e-3; x0 = 36.66; | |
1987 | for(Int_t i=0; i<steps; i++) { | |
1988 | x = xin + i*dx + dx/2; | |
1989 | ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1); | |
1990 | if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) { | |
1991 | holeYmax = x*TMath::Tan(0.5*alpha); | |
1992 | ppl->SetHole(holeYmax, holeZmax); | |
1993 | } | |
1994 | if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) { | |
1995 | holeYmax = x*TMath::Tan(0.5*alpha); | |
1996 | ppl->SetHole(holeYmax, holeZmax); | |
1997 | } | |
1998 | InsertLayer(ppl); | |
1999 | } | |
2000 | } | |
2001 | ||
2002 | // Space between the TRD and RICH | |
2003 | Double_t xRICH = 500.; | |
2004 | xin = xout; xout = xRICH; | |
2005 | steps = 200; dx = (xout - xin)/steps; rho = 1.29e-3; x0 = 36.66; | |
2006 | for(Int_t i=0; i<steps; i++) { | |
2007 | x = xin + i*dx + dx/2; | |
2008 | ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1); | |
2009 | InsertLayer(ppl); | |
2010 | } | |
2011 | ||
2012 | MapTimeBinLayers(); | |
ca6c93d6 | 2013 | delete [] Zc; |
2014 | delete [] Zmax; | |
5443e65e | 2015 | |
2016 | } | |
2017 | ||
2018 | //______________________________________________________ | |
2019 | ||
2020 | Int_t AliTRDtracker::AliTRDtrackingSector::CookTimeBinIndex(Int_t plane, Int_t local_tb) const | |
2021 | { | |
2022 | // | |
2023 | // depending on the digitization parameters calculates "global" | |
2024 | // time bin index for timebin <local_tb> in plane <plane> | |
2025 | // | |
2026 | ||
2027 | Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region | |
2028 | Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region | |
2029 | Double_t dx = (Double_t) fPar->GetTimeBinSize(); | |
2030 | ||
2031 | Int_t tbAmp = fPar->GetTimeBefore(); | |
2032 | Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx); | |
2033 | Int_t tbDrift = fPar->GetTimeMax(); | |
2034 | Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx); | |
2035 | ||
2036 | Int_t tb_per_plane = TMath::Min(tbAmp,maxAmp) + TMath::Min(tbDrift,maxDrift); | |
2037 | ||
2038 | Int_t gtb = (plane+1) * tb_per_plane - local_tb - 1; | |
2039 | ||
2040 | if((local_tb < 0) && | |
2041 | (TMath::Abs(local_tb) > TMath::Min(tbAmp,maxAmp))) return -1; | |
2042 | if(local_tb >= TMath::Min(tbDrift,maxDrift)) return -1; | |
2043 | ||
2044 | return gtb; | |
2045 | ||
2046 | ||
2047 | } | |
2048 | ||
2049 | //______________________________________________________ | |
2050 | ||
2051 | void AliTRDtracker::AliTRDtrackingSector::MapTimeBinLayers() | |
2052 | { | |
2053 | // | |
2054 | // For all sensitive time bins sets corresponding layer index | |
2055 | // in the array fTimeBins | |
2056 | // | |
2057 | ||
2058 | Int_t index; | |
2059 | ||
2060 | for(Int_t i = 0; i < fN; i++) { | |
2061 | index = fLayers[i]->GetTimeBinIndex(); | |
2062 | ||
2063 | // printf("gtb %d -> pl %d -> x %f \n", index, i, fLayers[i]->GetX()); | |
2064 | ||
2065 | if(index < 0) continue; | |
2066 | if(index >= (Int_t) kMAX_TIME_BIN_INDEX) { | |
2067 | printf("*** AliTRDtracker::MapTimeBinLayers: \n"); | |
2068 | printf(" index %d exceeds allowed maximum of %d!\n", | |
2069 | index, kMAX_TIME_BIN_INDEX-1); | |
2070 | continue; | |
2071 | } | |
2072 | fTimeBinIndex[index] = i; | |
2073 | } | |
2074 | ||
2075 | Double_t x1, dx1, x2, dx2, gap; | |
2076 | ||
2077 | for(Int_t i = 0; i < fN-1; i++) { | |
2078 | x1 = fLayers[i]->GetX(); | |
2079 | dx1 = fLayers[i]->GetdX(); | |
2080 | x2 = fLayers[i+1]->GetX(); | |
2081 | dx2 = fLayers[i+1]->GetdX(); | |
2082 | gap = (x2 - dx2/2) - (x1 + dx1/2); | |
2083 | if(gap < -0.01) { | |
2084 | printf("*** warning: layers %d and %d are overlayed:\n",i,i+1); | |
2085 | printf(" %f + %f + %f > %f\n", x1, dx1/2, dx2/2, x2); | |
2086 | } | |
2087 | if(gap > 0.01) { | |
2088 | printf("*** warning: layers %d and %d have a large gap:\n",i,i+1); | |
2089 | printf(" (%f - %f) - (%f + %f) = %f\n", | |
2090 | x2, dx2/2, x1, dx1, gap); | |
2091 | } | |
2092 | } | |
2093 | } | |
2094 | ||
2095 | ||
2096 | //______________________________________________________ | |
2097 | ||
2098 | ||
2099 | Int_t AliTRDtracker::AliTRDtrackingSector::GetLayerNumber(Double_t x) const | |
2100 | { | |
2101 | // | |
2102 | // Returns the number of time bin which in radial position is closest to <x> | |
2103 | // | |
2104 | ||
2105 | if(x >= fLayers[fN-1]->GetX()) return fN-1; | |
2106 | if(x <= fLayers[0]->GetX()) return 0; | |
2107 | ||
2108 | Int_t b=0, e=fN-1, m=(b+e)/2; | |
2109 | for (; b<e; m=(b+e)/2) { | |
2110 | if (x > fLayers[m]->GetX()) b=m+1; | |
2111 | else e=m; | |
2112 | } | |
2113 | if(TMath::Abs(x - fLayers[m]->GetX()) > | |
2114 | TMath::Abs(x - fLayers[m+1]->GetX())) return m+1; | |
2115 | else return m; | |
2116 | ||
2117 | } | |
2118 | ||
2119 | //______________________________________________________ | |
2120 | ||
2121 | Int_t AliTRDtracker::AliTRDtrackingSector::GetInnerTimeBin() const | |
2122 | { | |
2123 | // | |
2124 | // Returns number of the innermost SENSITIVE propagation layer | |
2125 | // | |
2126 | ||
2127 | return GetLayerNumber(0); | |
2128 | } | |
2129 | ||
2130 | //______________________________________________________ | |
2131 | ||
2132 | Int_t AliTRDtracker::AliTRDtrackingSector::GetOuterTimeBin() const | |
2133 | { | |
2134 | // | |
2135 | // Returns number of the outermost SENSITIVE time bin | |
2136 | // | |
2137 | ||
2138 | return GetLayerNumber(GetNumberOfTimeBins() - 1); | |
46d29e70 | 2139 | } |
2140 | ||
5443e65e | 2141 | //______________________________________________________ |
2142 | ||
2143 | Int_t AliTRDtracker::AliTRDtrackingSector::GetNumberOfTimeBins() const | |
2144 | { | |
2145 | // | |
2146 | // Returns number of SENSITIVE time bins | |
2147 | // | |
2148 | ||
2149 | Int_t tb, layer; | |
2150 | for(tb = kMAX_TIME_BIN_INDEX-1; tb >=0; tb--) { | |
2151 | layer = GetLayerNumber(tb); | |
2152 | if(layer>=0) break; | |
2153 | } | |
2154 | return tb+1; | |
2155 | } | |
2156 | ||
2157 | //______________________________________________________ | |
2158 | ||
2159 | void AliTRDtracker::AliTRDtrackingSector::InsertLayer(AliTRDpropagationLayer* pl) | |
2160 | { | |
2161 | // | |
2162 | // Insert layer <pl> in fLayers array. | |
2163 | // Layers are sorted according to X coordinate. | |
2164 | ||
2165 | if ( fN == ((Int_t) kMAX_LAYERS_PER_SECTOR)) { | |
2166 | printf("AliTRDtrackingSector::InsertLayer(): Too many layers !\n"); | |
2167 | return; | |
2168 | } | |
2169 | if (fN==0) {fLayers[fN++] = pl; return;} | |
2170 | Int_t i=Find(pl->GetX()); | |
2171 | ||
2172 | memmove(fLayers+i+1 ,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*)); | |
2173 | fLayers[i]=pl; fN++; | |
2174 | ||
2175 | } | |
2176 | ||
2177 | //______________________________________________________ | |
2178 | ||
2179 | Int_t AliTRDtracker::AliTRDtrackingSector::Find(Double_t x) const | |
2180 | { | |
2181 | // | |
2182 | // Returns index of the propagation layer nearest to X | |
2183 | // | |
2184 | ||
2185 | if (x <= fLayers[0]->GetX()) return 0; | |
2186 | if (x > fLayers[fN-1]->GetX()) return fN; | |
2187 | Int_t b=0, e=fN-1, m=(b+e)/2; | |
2188 | for (; b<e; m=(b+e)/2) { | |
2189 | if (x > fLayers[m]->GetX()) b=m+1; | |
2190 | else e=m; | |
2191 | } | |
2192 | return m; | |
2193 | } | |
2194 | ||
2195 | //______________________________________________________ | |
2196 | ||
2197 | void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters( | |
2198 | Double_t y, Double_t z, Double_t &dx, Double_t &rho, Double_t &x0, | |
2199 | Bool_t &lookForCluster) const | |
2200 | { | |
2201 | // | |
2202 | // Returns radial step <dx>, density <rho>, rad. length <x0>, | |
2203 | // and sensitivity <lookForCluster> in point <y,z> | |
2204 | // | |
2205 | ||
2206 | dx = fdX; | |
2207 | rho = fRho; | |
2208 | x0 = fX0; | |
2209 | lookForCluster = kFALSE; | |
2210 | ||
2211 | // check dead regions | |
2212 | if(fTimeBinIndex >= 0) { | |
2213 | for(Int_t ch = 0; ch < (Int_t) kZONES; ch++) { | |
2214 | if(TMath::Abs(z - fZc[ch]) < fZmax[ch]) | |
2215 | lookForCluster = kTRUE; | |
2216 | } | |
2217 | if(TMath::Abs(y) > fYmax) lookForCluster = kFALSE; | |
2218 | if(!lookForCluster) { | |
2219 | // rho = 1.7; x0 = 33.0; // G10 | |
2220 | } | |
2221 | } | |
2222 | ||
2223 | // check hole | |
2224 | if(fHole && (TMath::Abs(y - fHoleYc) < fHoleYmax) && | |
2225 | (TMath::Abs(z - fHoleZc) < fHoleZmax)) { | |
2226 | lookForCluster = kFALSE; | |
2227 | rho = fHoleRho; | |
2228 | x0 = fHoleX0; | |
2229 | } | |
2230 | ||
2231 | return; | |
2232 | } | |
2233 | ||
2234 | //______________________________________________________ | |
2235 | ||
2236 | void AliTRDtracker::AliTRDpropagationLayer::InsertCluster(AliTRDcluster* c, | |
2237 | UInt_t index) { | |
2238 | ||
2239 | // Insert cluster in cluster array. | |
2240 | // Clusters are sorted according to Y coordinate. | |
2241 | ||
2242 | if(fTimeBinIndex < 0) { | |
2243 | printf("*** attempt to insert cluster into non-sensitive time bin!\n"); | |
2244 | return; | |
2245 | } | |
2246 | ||
2247 | if (fN== (Int_t) kMAX_CLUSTER_PER_TIME_BIN) { | |
2248 | printf("AliTRDpropagationLayer::InsertCluster(): Too many clusters !\n"); | |
2249 | return; | |
2250 | } | |
2251 | if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;} | |
2252 | Int_t i=Find(c->GetY()); | |
2253 | memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTRDcluster*)); | |
2254 | memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t)); | |
2255 | fIndex[i]=index; fClusters[i]=c; fN++; | |
2256 | } | |
2257 | ||
2258 | //______________________________________________________ | |
2259 | ||
2260 | Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Double_t y) const { | |
2261 | ||
2262 | // Returns index of the cluster nearest in Y | |
2263 | ||
2264 | if (y <= fClusters[0]->GetY()) return 0; | |
2265 | if (y > fClusters[fN-1]->GetY()) return fN; | |
2266 | Int_t b=0, e=fN-1, m=(b+e)/2; | |
2267 | for (; b<e; m=(b+e)/2) { | |
2268 | if (y > fClusters[m]->GetY()) b=m+1; | |
2269 | else e=m; | |
2270 | } | |
2271 | return m; | |
2272 | } | |
2273 | ||
2274 | ||
2275 | ||
2276 | ||
2277 | ||
2278 | ||
2279 | ||
2280 | ||
2281 |