adding map include file to solve compilation issue (bug https://savannah.cern.ch...
[u/mrichter/AliRoot.git] / HLT / ITS / tracking / AliHLTITSLayer.cxx
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"
20 #include <algorithm>
21
22
23 //------------------------------------------------------------------------
24 AliHLTITSLayer::AliHLTITSLayer():
25 fR(0),
26 fPhiOffset(0),
27 fNladders(0),
28 fZOffset(0),
29 fNdetectors(0),
30 fDetectors(0),
31 fN(0),
32 fDy5(0),
33 fDy10(0),
34 fDy20(0),
35 fClustersCs(0),
36 fClusterIndexCs(0),
37 fYcs(0),
38 fZcs(0),
39 fNcs(0),
40 fCurrentSlice(-1),
41 fZmax(0),
42 fYmin(0),
43 fYmax(0),
44 fI(0),
45 fImax(0),
46 fSkip(0),
47 fAccepted(0),
48 fRoad(0){
49   //--------------------------------------------------------------------
50   //default AliHLTITSLayer constructor
51   //--------------------------------------------------------------------
52 }
53 //------------------------------------------------------------------------
54 AliHLTITSLayer::
55 AliHLTITSLayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
56 fR(r),
57 fPhiOffset(p),
58 fNladders(nl),
59 fZOffset(z),
60 fNdetectors(nd),
61 fDetectors(0),
62 fN(0),
63 fDy5(0),
64 fDy10(0),
65 fDy20(0),
66 fClustersCs(0),
67 fClusterIndexCs(0),
68 fYcs(0),
69 fZcs(0),
70 fNcs(0),
71 fCurrentSlice(-1),
72 fZmax(0),
73 fYmin(0),
74 fYmax(0),
75 fI(0),
76 fImax(0),
77 fSkip(0),
78 fAccepted(0),
79 fRoad(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 //------------------------------------------------------------------------
87 AliHLTITSLayer::AliHLTITSLayer(const AliHLTITSLayer& layer):
88 fR(layer.fR),
89 fPhiOffset(layer.fPhiOffset),
90 fNladders(layer.fNladders),
91 fZOffset(layer.fZOffset),
92 fNdetectors(layer.fNdetectors),
93 fDetectors(layer.fDetectors),
94 fN(layer.fN),
95 fDy5(layer.fDy5),
96 fDy10(layer.fDy10),
97 fDy20(layer.fDy20),
98 fClustersCs(layer.fClustersCs),
99 fClusterIndexCs(layer.fClusterIndexCs),
100 fYcs(layer.fYcs),
101 fZcs(layer.fZcs),
102 fNcs(layer.fNcs),
103 fCurrentSlice(layer.fCurrentSlice),
104 fZmax(layer.fZmax),
105 fYmin(layer.fYmin),
106 fYmax(layer.fYmax),
107 fI(layer.fI),
108 fImax(layer.fImax),
109 fSkip(layer.fSkip),
110 fAccepted(layer.fAccepted),
111 fRoad(layer.fRoad){
112   //Copy constructor
113 }
114 //------------------------------------------------------------------------
115 AliHLTITSLayer::~AliHLTITSLayer() {
116   //--------------------------------------------------------------------
117   // AliHLTITSLayer destructor
118   //--------------------------------------------------------------------
119   delete [] fDetectors;
120 }
121 //------------------------------------------------------------------------
122 void AliHLTITSLayer::ResetClusters() {
123   //--------------------------------------------------------------------
124   // This function removes loaded clusters
125   //--------------------------------------------------------------------
126   
127   fN=0;
128   fI=0;
129 }
130
131
132 //------------------------------------------------------------------------
133 void 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 //------------------------------------------------------------------------
144 Int_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 }
163
164
165 //------------------------------------------------------------------------
166 void  AliHLTITSLayer::SortClusters()
167 {
168   //
169   //sort clusters
170   //
171  
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);
180   
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 //------------------------------------------------------------------------
273 Int_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 //------------------------------------------------------------------------
299 void AliHLTITSLayer::
300 SelectClusters(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 //------------------------------------------------------------------------
379 Int_t AliHLTITSLayer::
380 FindDetectorIndex(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 //------------------------------------------------------------------------
418 const 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 //------------------------------------------------------------------------
450 Double_t AliHLTITSLayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
451 const {
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 //------------------------------------------------------------------------
549 Int_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