adding map include file to solve compilation issue (bug https://savannah.cern.ch...
[u/mrichter/AliRoot.git] / HLT / ITS / tracking / AliHLTITSLayer.cxx
CommitLineData
2f399afc 1// **************************************************************************
2// This file is property of and copyright by the ALICE HLT Project *
3// ALICE Experiment at CERN, All rights reserved. *
4// *
5// Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
6// for The ALICE HLT Project. *
7// *
8// Permission to use, copy, modify and distribute this software and its *
9// documentation strictly for non-commercial purposes is hereby granted *
10// without fee, provided that the above copyright notice appears in all *
11// copies and that both the copyright notice and this permission notice *
12// appear in the supporting documentation. The authors make no claims *
13// about the suitability of this software for any purpose. It is *
14// provided "as is" without express or implied warranty. *
15// *
16//***************************************************************************
17
18
19#include "AliHLTITSLayer.h"
9eee44f6 20#include <algorithm>
21
2f399afc 22
23//------------------------------------------------------------------------
24AliHLTITSLayer::AliHLTITSLayer():
25fR(0),
26fPhiOffset(0),
27fNladders(0),
28fZOffset(0),
29fNdetectors(0),
30fDetectors(0),
31fN(0),
32fDy5(0),
33fDy10(0),
34fDy20(0),
35fClustersCs(0),
36fClusterIndexCs(0),
37fYcs(0),
38fZcs(0),
39fNcs(0),
40fCurrentSlice(-1),
41fZmax(0),
42fYmin(0),
43fYmax(0),
44fI(0),
45fImax(0),
46fSkip(0),
47fAccepted(0),
48fRoad(0){
49 //--------------------------------------------------------------------
50 //default AliHLTITSLayer constructor
51 //--------------------------------------------------------------------
2f399afc 52}
53//------------------------------------------------------------------------
54AliHLTITSLayer::
55AliHLTITSLayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
56fR(r),
57fPhiOffset(p),
58fNladders(nl),
59fZOffset(z),
60fNdetectors(nd),
61fDetectors(0),
62fN(0),
63fDy5(0),
64fDy10(0),
65fDy20(0),
66fClustersCs(0),
67fClusterIndexCs(0),
68fYcs(0),
69fZcs(0),
70fNcs(0),
71fCurrentSlice(-1),
72fZmax(0),
73fYmin(0),
74fYmax(0),
75fI(0),
76fImax(0),
77fSkip(0),
78fAccepted(0),
79fRoad(0) {
80 //--------------------------------------------------------------------
81 //main AliHLTITSLayer constructor
82 //--------------------------------------------------------------------
83 fDetectors=new AliHLTITSDetector[fNladders*fNdetectors];
84 fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
85}
86//------------------------------------------------------------------------
87AliHLTITSLayer::AliHLTITSLayer(const AliHLTITSLayer& layer):
88fR(layer.fR),
89fPhiOffset(layer.fPhiOffset),
90fNladders(layer.fNladders),
91fZOffset(layer.fZOffset),
92fNdetectors(layer.fNdetectors),
93fDetectors(layer.fDetectors),
94fN(layer.fN),
95fDy5(layer.fDy5),
96fDy10(layer.fDy10),
97fDy20(layer.fDy20),
98fClustersCs(layer.fClustersCs),
99fClusterIndexCs(layer.fClusterIndexCs),
100fYcs(layer.fYcs),
101fZcs(layer.fZcs),
102fNcs(layer.fNcs),
103fCurrentSlice(layer.fCurrentSlice),
104fZmax(layer.fZmax),
105fYmin(layer.fYmin),
106fYmax(layer.fYmax),
107fI(layer.fI),
108fImax(layer.fImax),
109fSkip(layer.fSkip),
110fAccepted(layer.fAccepted),
111fRoad(layer.fRoad){
112 //Copy constructor
113}
114//------------------------------------------------------------------------
115AliHLTITSLayer::~AliHLTITSLayer() {
116 //--------------------------------------------------------------------
117 // AliHLTITSLayer destructor
118 //--------------------------------------------------------------------
119 delete [] fDetectors;
2f399afc 120}
121//------------------------------------------------------------------------
122void AliHLTITSLayer::ResetClusters() {
123 //--------------------------------------------------------------------
124 // This function removes loaded clusters
125 //--------------------------------------------------------------------
2f399afc 126
127 fN=0;
128 fI=0;
129}
2f399afc 130
ef6e2aa2 131
2f399afc 132//------------------------------------------------------------------------
133void AliHLTITSLayer::ResetRoad() {
134 //--------------------------------------------------------------------
135 // This function calculates the road defined by the cluster density
136 //--------------------------------------------------------------------
137 Int_t n=0;
138 for (Int_t i=0; i<fN; i++) {
139 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
140 }
141 if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
142}
143//------------------------------------------------------------------------
144Int_t AliHLTITSLayer::InsertCluster(AliITSRecPoint *cl) {
145 //--------------------------------------------------------------------
146 //This function adds a cluster to this layer
147 //--------------------------------------------------------------------
148 if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
149 ::Error("InsertCluster","Too many clusters !\n");
150 return 1;
151 }
152 fCurrentSlice=-1;
153 fClusters[fN]=cl;
154 fN++;
155 AliHLTITSDetector &det=GetDetector(cl->GetDetectorIndex());
156 if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
157 if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
158 if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
159 if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
160
161 return 0;
162}
9eee44f6 163
164
2f399afc 165//------------------------------------------------------------------------
166void AliHLTITSLayer::SortClusters()
167{
168 //
169 //sort clusters
170 //
9eee44f6 171
2f399afc 172 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
173 Float_t *z = new Float_t[fN];
174 Int_t * index = new Int_t[fN];
175 //
176 for (Int_t i=0;i<fN;i++){
177 z[i] = fClusters[i]->GetZ();
178 }
179 TMath::Sort(fN,z,index,kFALSE);
9eee44f6 180
2f399afc 181 for (Int_t i=0;i<fN;i++){
182 clusters[i] = fClusters[index[i]];
183 }
184 //
185 for (Int_t i=0;i<fN;i++){
186 fClusters[i] = clusters[i];
187 fZ[i] = fClusters[i]->GetZ();
188 AliHLTITSDetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
189 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
190 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
191 fY[i] = y;
192 }
193 delete[] index;
194 delete[] z;
195 delete[] clusters;
196 //
197
198 fYB[0]=10000000;
199 fYB[1]=-10000000;
200 for (Int_t i=0;i<fN;i++){
201 if (fY[i]<fYB[0]) fYB[0]=fY[i];
202 if (fY[i]>fYB[1]) fYB[1]=fY[i];
203 fClusterIndex[i] = i;
204 }
205 //
206 // fill slices
207 fDy5 = (fYB[1]-fYB[0])/5.;
208 fDy10 = (fYB[1]-fYB[0])/10.;
209 fDy20 = (fYB[1]-fYB[0])/20.;
210 for (Int_t i=0;i<6;i++) fN5[i] =0;
211 for (Int_t i=0;i<11;i++) fN10[i]=0;
212 for (Int_t i=0;i<21;i++) fN20[i]=0;
213 //
214 for (Int_t i=0;i<6;i++) {fBy5[i][0] = fYB[0]+(i-0.75)*fDy5; fBy5[i][1] = fYB[0]+(i+0.75)*fDy5;}
215 for (Int_t i=0;i<11;i++) {fBy10[i][0] = fYB[0]+(i-0.75)*fDy10; fBy10[i][1] = fYB[0]+(i+0.75)*fDy10;}
216 for (Int_t i=0;i<21;i++) {fBy20[i][0] = fYB[0]+(i-0.75)*fDy20; fBy20[i][1] = fYB[0]+(i+0.75)*fDy20;}
217 //
218 //
219 for (Int_t i=0;i<fN;i++)
220 for (Int_t irot=-1;irot<=1;irot++){
221 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
222 // slice 5
223 for (Int_t slice=0; slice<6;slice++){
224 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
225 fClusters5[slice][fN5[slice]] = fClusters[i];
226 fY5[slice][fN5[slice]] = curY;
227 fZ5[slice][fN5[slice]] = fZ[i];
228 fClusterIndex5[slice][fN5[slice]]=i;
229 fN5[slice]++;
230 }
231 }
232 // slice 10
233 for (Int_t slice=0; slice<11;slice++){
234 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
235 fClusters10[slice][fN10[slice]] = fClusters[i];
236 fY10[slice][fN10[slice]] = curY;
237 fZ10[slice][fN10[slice]] = fZ[i];
238 fClusterIndex10[slice][fN10[slice]]=i;
239 fN10[slice]++;
240 }
241 }
242 // slice 20
243 for (Int_t slice=0; slice<21;slice++){
244 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
245 fClusters20[slice][fN20[slice]] = fClusters[i];
246 fY20[slice][fN20[slice]] = curY;
247 fZ20[slice][fN20[slice]] = fZ[i];
248 fClusterIndex20[slice][fN20[slice]]=i;
249 fN20[slice]++;
250 }
251 }
252 }
253
254 //
255 // consistency check
256 //
257 for (Int_t i=0;i<fN-1;i++){
258 if (fZ[i]>fZ[i+1]){
259 printf("Bug\n");
260 }
261 }
262 //
263 for (Int_t slice=0;slice<21;slice++)
264 for (Int_t i=0;i<fN20[slice]-1;i++){
265 if (fZ20[slice][i]>fZ20[slice][i+1]){
266 printf("Bug\n");
267 }
268 }
269
270
271}
272//------------------------------------------------------------------------
273Int_t AliHLTITSLayer::FindClusterIndex(Float_t z) const {
274 //--------------------------------------------------------------------
275 // This function returns the index of the nearest cluster
276 //--------------------------------------------------------------------
277 Int_t ncl=0;
278 const Float_t *zcl;
279 if (fCurrentSlice<0) {
280 ncl = fN;
281 zcl = fZ;
282 }
283 else{
284 ncl = fNcs;
285 zcl = fZcs;;
286 }
287
288 if (ncl==0) return 0;
289 Int_t b=0, e=ncl-1, m=(b+e)/2;
290 for (; b<e; m=(b+e)/2) {
291 // if (z > fClusters[m]->GetZ()) b=m+1;
292 if (z > zcl[m]) b=m+1;
293 else e=m;
294 }
295 return m;
296}
297
298//------------------------------------------------------------------------
299void AliHLTITSLayer::
300SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
301 //--------------------------------------------------------------------
302 // This function sets the "window"
303 //--------------------------------------------------------------------
304
305 Double_t circle=2*TMath::Pi()*fR;
306 fYmin = ymin; fYmax =ymax;
307 Float_t ymiddle = (fYmin+fYmax)*0.5;
308 if (ymiddle<fYB[0]) {
309 fYmin+=circle; fYmax+=circle; ymiddle+=circle;
310 } else if (ymiddle>fYB[1]) {
311 fYmin-=circle; fYmax-=circle; ymiddle-=circle;
312 }
313
314 //
315 fCurrentSlice =-1;
316 // defualt take all
317 fClustersCs = fClusters;
318 fClusterIndexCs = fClusterIndex;
319 fYcs = fY;
320 fZcs = fZ;
321 fNcs = fN;
322 //
323 //is in 20 slice?
324 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
325 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
326 if (slice<0) slice=0;
327 if (slice>20) slice=20;
328 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
329 if (isOK) {
330 fCurrentSlice=slice;
331 fClustersCs = fClusters20[fCurrentSlice];
332 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
333 fYcs = fY20[fCurrentSlice];
334 fZcs = fZ20[fCurrentSlice];
335 fNcs = fN20[fCurrentSlice];
336 }
337 }
338 //
339 //is in 10 slice?
340 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
341 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
342 if (slice<0) slice=0;
343 if (slice>10) slice=10;
344 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
345 if (isOK) {
346 fCurrentSlice=slice;
347 fClustersCs = fClusters10[fCurrentSlice];
348 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
349 fYcs = fY10[fCurrentSlice];
350 fZcs = fZ10[fCurrentSlice];
351 fNcs = fN10[fCurrentSlice];
352 }
353 }
354 //
355 //is in 5 slice?
356 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
357 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
358 if (slice<0) slice=0;
359 if (slice>5) slice=5;
360 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
361 if (isOK) {
362 fCurrentSlice=slice;
363 fClustersCs = fClusters5[fCurrentSlice];
364 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
365 fYcs = fY5[fCurrentSlice];
366 fZcs = fZ5[fCurrentSlice];
367 fNcs = fN5[fCurrentSlice];
368 }
369 }
370 //
371 fI=FindClusterIndex(zmin); fZmax=zmax;
372 fImax = TMath::Min(FindClusterIndex(zmax)+1,fNcs);
373 fSkip = 0;
374 fAccepted =0;
375
376 return;
377}
378//------------------------------------------------------------------------
379Int_t AliHLTITSLayer::
380FindDetectorIndex(Double_t phi, Double_t z) const {
381 //--------------------------------------------------------------------
382 //This function finds the detector crossed by the track
383 //--------------------------------------------------------------------
384 Double_t dphi;
385 if (fZOffset<0) // old geometry
386 dphi = -(phi-fPhiOffset);
387 else // new geometry
388 dphi = phi-fPhiOffset;
389
390
391 if (dphi < 0) dphi += 2*TMath::Pi();
392 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
393 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
394 if (np>=fNladders) np-=fNladders;
395 if (np<0) np+=fNladders;
396
397
398 Double_t dz=fZOffset-z;
399 Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
400 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
401 if (nz>=fNdetectors) return -1;
402 if (nz<0) return -1;
403
404 // ad hoc correction for 3rd ladder of SDD inner layer,
405 // which is reversed (rotated by pi around local y)
406 // this correction is OK only from AliITSv11Hybrid onwards
407 if (GetR()>12. && GetR()<20.) { // SDD inner
408 if(np==2) { // 3rd ladder
409 nz = (fNdetectors-1) - nz;
410 }
411 }
412 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
413
414
415 return np*fNdetectors + nz;
416}
417//------------------------------------------------------------------------
418const AliITSRecPoint *AliHLTITSLayer::GetNextCluster(Int_t &ci,Bool_t test)
419{
420 //--------------------------------------------------------------------
421 // This function returns clusters within the "window"
422 //--------------------------------------------------------------------
423
424 if (fCurrentSlice<0) {
425 Double_t rpi2 = 2.*fR*TMath::Pi();
426 for (Int_t i=fI; i<fImax; i++) {
427 Double_t y = fY[i];
428 if (fYmax<y) y -= rpi2;
429 if (fYmin>y) y += rpi2;
430 if (y<fYmin) continue;
431 if (y>fYmax) continue;
432 if (fClusters[i]->GetQ()==0&&fSkip==2) continue;
433 ci=i;
434 if (!test) fI=i+1;
435 return fClusters[i];
436 }
437 } else {
438 for (Int_t i=fI; i<fImax; i++) {
439 if (fYcs[i]<fYmin) continue;
440 if (fYcs[i]>fYmax) continue;
441 if (fClustersCs[i]->GetQ()==0&&fSkip==2) continue;
442 ci=fClusterIndexCs[i];
443 if (!test) fI=i+1;
444 return fClustersCs[i];
445 }
446 }
447 return 0;
448}
449//------------------------------------------------------------------------
450Double_t AliHLTITSLayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
451const {
452 //--------------------------------------------------------------------
453 // This function returns the layer thickness at this point (units X0)
454 //--------------------------------------------------------------------
455 Double_t d=0.0085;
456 x0=AliITSRecoParam::GetX0Air();
457 if (43<fR&&fR<45) { //SSD2
458 Double_t dd=0.0034;
459 d=dd;
460 if (TMath::Abs(y-0.00)>3.40) d+=dd;
461 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
462 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
463 for (Int_t i=0; i<12; i++) {
464 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
465 if (TMath::Abs(y-0.00)>3.40) d+=dd;
466 d+=0.0034;
467 break;
468 }
469 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
470 if (TMath::Abs(y-0.00)>3.40) d+=dd;
471 d+=0.0034;
472 break;
473 }
474 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
475 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
476 }
477 } else
478 if (37<fR&&fR<41) { //SSD1
479 Double_t dd=0.0034;
480 d=dd;
481 if (TMath::Abs(y-0.00)>3.40) d+=dd;
482 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
483 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
484 for (Int_t i=0; i<11; i++) {
485 if (TMath::Abs(z-3.9*i)<0.15) {
486 if (TMath::Abs(y-0.00)>3.40) d+=dd;
487 d+=dd;
488 break;
489 }
490 if (TMath::Abs(z+3.9*i)<0.15) {
491 if (TMath::Abs(y-0.00)>3.40) d+=dd;
492 d+=dd;
493 break;
494 }
495 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
496 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
497 }
498 } else
499 if (13<fR&&fR<26) { //SDD
500 Double_t dd=0.0033;
501 d=dd;
502 if (TMath::Abs(y-0.00)>3.30) d+=dd;
503
504 if (TMath::Abs(y-1.80)<0.55) {
505 d+=0.016;
506 for (Int_t j=0; j<20; j++) {
507 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
508 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
509 }
510 }
511 if (TMath::Abs(y+1.80)<0.55) {
512 d+=0.016;
513 for (Int_t j=0; j<20; j++) {
514 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
515 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
516 }
517 }
518
519 for (Int_t i=0; i<4; i++) {
520 if (TMath::Abs(z-7.3*i)<0.60) {
521 d+=dd;
522 if (TMath::Abs(y-0.00)>3.30) d+=dd;
523 break;
524 }
525 if (TMath::Abs(z+7.3*i)<0.60) {
526 d+=dd;
527 if (TMath::Abs(y-0.00)>3.30) d+=dd;
528 break;
529 }
530 }
531 } else
532 if (6<fR&&fR<8) { //SPD2
533 Double_t dd=0.0063; x0=21.5;
534 d=dd;
535 if (TMath::Abs(y-3.08)>0.5) d+=dd;
536 if (TMath::Abs(y-3.03)<0.10) d+=0.014;
537 } else
538 if (3<fR&&fR<5) { //SPD1
539 Double_t dd=0.0063; x0=21.5;
540 d=dd;
541 if (TMath::Abs(y+0.21)>0.6) d+=dd;
542 if (TMath::Abs(y+0.10)<0.10) d+=0.014;
543 }
544
545 return d;
546}
547
548//------------------------------------------------------------------------
549Int_t AliHLTITSLayer::InRoad() const {
550 //-------------------------------------------------------------------
551 // This function returns number of clusters within the "window"
552 //--------------------------------------------------------------------
553 Int_t ncl=0;
554 for (Int_t i=fI; i<fN; i++) {
555 const AliITSRecPoint *c=fClusters[i];
556 if (c->GetZ() > fZmax) break;
557 if (c->IsUsed()) continue;
558 const AliHLTITSDetector &det=GetDetector(c->GetDetectorIndex());
559 Double_t y=fR*det.GetPhi() + c->GetY();
560
561 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
562 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
563
564 if (y<fYmin) continue;
565 if (y>fYmax) continue;
566 ncl++;
567 }
568 return ncl;
569}
570