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