Fix in layers rmin/rmax calculation and their accounting in track navigation.
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSURecoLayer.cxx
1 #include <TClonesArray.h>
2 #include "AliITSURecoLayer.h"
3 #include "AliITSsegmentation.h"
4 #include "AliITSUAux.h"
5 #include "AliITSUClusterPix.h"
6 #include "AliITSUGeomTGeo.h"
7 #include "AliLog.h"
8
9 using namespace AliITSUAux;
10 using namespace TMath;
11
12 ClassImp(AliITSURecoLayer)
13
14
15 //______________________________________________________
16 AliITSURecoLayer::AliITSURecoLayer(const char* name)
17   :fActiveID(-1)
18   ,fNSensors(0)
19   ,fNSensInLadder(0)
20   ,fNLadders(0)
21   ,fR(0)
22   ,fRMax(0)
23   ,fRMin(0)
24   ,fZMax(0)
25   ,fZMin(0)
26   ,fPhiLadMax(0)
27   ,fPhiLadMin(0)
28   ,fPhiOffs(0)
29   ,fSensDZInv(0)
30   ,fDPhiLadInv(0)
31   ,fMaxStep(0.5)
32   ,fSensors(0)
33   ,fITSGeom(0)
34   ,fClusters(0)
35 {
36   // def. c-tor
37   SetNameTitle(name,name);
38 }
39
40 //______________________________________________________
41 AliITSURecoLayer::AliITSURecoLayer(const char* name, Int_t activeID, AliITSUGeomTGeo* gm)
42   :fActiveID(activeID)
43   ,fNSensors(0)
44   ,fNSensInLadder(0)
45   ,fNLadders(0)
46   ,fR(0)
47   ,fRMax(0)
48   ,fRMin(0)
49   ,fZMax(0)
50   ,fZMin(0)
51   ,fPhiLadMax(0)
52   ,fPhiLadMin(0)
53   ,fPhiOffs(0)
54   ,fSensDZInv(0)
55   ,fDPhiLadInv(0)
56   ,fMaxStep(0.5)
57   ,fSensors(0)
58   ,fITSGeom(gm)
59   ,fClusters(0)
60 {
61   // def. c-tor
62   SetNameTitle(name,name);
63   Build();
64 }
65
66 //______________________________________________________
67 AliITSURecoLayer::~AliITSURecoLayer()
68 {
69   // def. d-tor
70   delete[] fSensors;
71   delete[] fPhiLadMax;
72   delete[] fPhiLadMin;
73   if (GetOwnsClusterArray()) delete fClusters;
74 }
75
76 //______________________________________________________
77 void AliITSURecoLayer::Print(Option_t* opt) const                             
78 {
79   //print 
80   printf("Lr %-15s %d (act:%+d), NSens: %4d | MaxStep:%.2f ",GetName(),GetID(),GetActiveID(),GetNSensors(),fMaxStep);
81   printf("%6.3f<R<%6.3f | %+8.3f<Z<%+8.3f dZ:%6.3f\n",fRMin,fRMax,fZMin,fZMax,fSensDZInv>0 ? 1/fSensDZInv : 0);
82   TString opts = opt; opts.ToLower();
83   if (opts.Contains("sn")) for (int i=0;i<GetNSensors();i++) GetSensor(i)->Print(opt);
84 }
85
86 //______________________________________________________
87 void AliITSURecoLayer::Build()
88 {
89   // build internal structures
90   const double kSafeR = 0.05; // safety margin for Rmin,Rmax of the layer
91   if (fActiveID<0) return;
92   fNLadders = fITSGeom->GetNLadders(fActiveID);
93   fNSensInLadder = fITSGeom->GetNDetectors(fActiveID);
94   fNSensors = fNLadders*fNSensInLadder;
95   fSensors = new AliITSURecoSens*[fNSensors];
96   const AliITSsegmentation* kSegm = fITSGeom->GetSegmentation(fActiveID);
97   //
98   // name layer according its active id, detector type and segmentation tyoe
99   TGeoHMatrix mmod;
100   const TGeoHMatrix* mt2l;
101   fRMin=fZMin=1e9;
102   fRMax=fZMax=-1e9;
103   double phiTF,rTF, loc[3]={0,0,0},glo[3];
104   fNSensors = 0;
105   fPhiLadMin = new Double_t[fNLadders];
106   fPhiLadMax = new Double_t[fNLadders];
107   fSensDZInv = 0;
108   fDPhiLadInv = fNLadders/TwoPi();
109   //
110   for (int ild=0;ild<fNLadders;ild++) {
111     fPhiLadMin[ild] = 1e9;
112     fPhiLadMax[ild] = -1e9;
113     //
114     for (int idt=0;idt<fNSensInLadder;idt++) {
115       AliITSURecoSens* sens = new AliITSURecoSens(fNSensors++);
116       fSensors[ild*fNSensInLadder+idt] = sens;
117       //
118       double phiMin=1e9,phiMax=-1e9,zMin=1e9,zMax=-1e9;
119       mmod = *fITSGeom->GetMatrixSens(fActiveID,ild,idt);
120       for (int ix=0;ix<2;ix++) {
121         loc[0] = (ix-0.5)*kSegm->Dx(); // +-DX/2
122         for (int iy=0;iy<2;iy++) {
123           loc[1] = (iy-0.5)*kSegm->Dy(); // +-DY/2
124           for (int iz=0;iz<2;iz++) {
125             loc[2] = (iz-0.5)*kSegm->Dz(); // +-DZ/2
126             //
127             mmod.LocalToMaster(loc,glo);
128             double phi = ATan2(glo[1],glo[0]);
129             double r   = glo[0]*glo[0] + glo[1]*glo[1];
130             if (fRMin>r) fRMin = r;
131             if (fRMax<r) fRMax = r;
132             BringTo02Pi(phi);
133             if      (phiMin>1e8) phiMin=phi; 
134             else if (!OKforPhiMin(phiMin,phi)) phiMin=phi;
135             if      (phiMax<-1e8) phiMax=phi;
136             else if (!OKforPhiMax(phiMax,phi)) phiMax=phi;            
137             if (glo[2]>zMax) zMax=glo[2];
138             if (glo[2]<zMin) zMin=glo[2];
139           }
140         }
141       }
142       sens->SetBoundaries(phiMin,phiMax,zMin,zMax);
143       mt2l = fITSGeom->GetMatrixT2L(fActiveID,ild,idt);
144       mmod.Multiply(mt2l);      
145       loc[0]=loc[1]=loc[2]=0;
146       mmod.LocalToMaster(loc,glo);
147       rTF   = Sqrt(glo[0]*glo[0] + glo[1]*glo[1]);  //  tracking params (misaligned)
148       phiTF = ATan2(glo[1],glo[0]);
149       BringTo02Pi(phiTF);
150       //
151       sens->SetXTF(rTF);
152       sens->SetPhiTF(phiTF);
153       //
154       if      (fPhiLadMin[ild]>1e8)  fPhiLadMin[ild] = phiMin;
155       else if (!OKforPhiMin(fPhiLadMin[ild],phiMin)) fPhiLadMin[ild] = phiMin;
156       if      (fPhiLadMax[ild]<-1e8) fPhiLadMax[ild] = phiMax;
157       else if (!OKforPhiMax(fPhiLadMax[ild],phiMax)) fPhiLadMax[ild] = phiMax;
158       if (fZMin>zMin) fZMin = zMin;
159       if (fZMax<zMax) fZMax = zMax;
160       //
161       if (idt>0) fSensDZInv += zMax - GetSensor(ild,idt-1)->GetZMax(); // z interval to previous
162     }
163   }
164   //
165   fRMin = Sqrt(fRMin);
166   fRMax = Sqrt(fRMax);
167   fR = 0.5*(fRMin+fRMax);
168   fRMin -= kSafeR;
169   fRMax += kSafeR;
170   double dz = fNSensInLadder>0 ? fSensDZInv/(fNSensInLadder-1)/fNLadders : fZMax-fZMin;
171   fSensDZInv = 1./dz;
172
173   const int kNBId[3][3] = { 
174     {AliITSURecoSens::kNghbBL,AliITSURecoSens::kNghbB,AliITSURecoSens::kNghbBR},
175     {AliITSURecoSens::kNghbL,          -1            ,AliITSURecoSens::kNghbR },
176     {AliITSURecoSens::kNghbTL,AliITSURecoSens::kNghbT,AliITSURecoSens::kNghbTR}
177   };
178
179   // add neighbours info
180   double zTol = 0.45*dz, phiTol = 0.45*TwoPi()/fNLadders;
181   for (int ild=0;ild<fNLadders;ild++) {
182     for (int idt=0;idt<fNSensInLadder;idt++) {
183       AliITSURecoSens* sens = GetSensor(ild,idt);
184       //
185       for (int ils=-1;ils<=1;ils++) {
186         int ildN = ild+ils;  // ladders of neighbouring sensors
187         if (ildN<0) ildN = fNLadders-1; else if (ildN==fNLadders) ildN = 0;
188         for (int ids=-1;ids<=1;ids++) {
189           int idtN = idt+ids;
190           if (idtN<0 || idtN==fNSensInLadder || (ids==0&&ils==0)) continue;
191           AliITSURecoSens* sensN = GetSensor(ildN,idtN); // potential neighbor
192           int neighbID = ildN*fNSensInLadder+idtN;
193           //      
194           int zType = 1;  // side
195           if (sens->GetZMin()-zTol  > sensN->GetZMax()) continue; // too large distance
196           if (sensN->GetZMin()-zTol > sens->GetZMax() ) continue; // too large distance
197           if      (sens->GetZMin()-zTol>sensN->GetZMin()) zType =  0;     // bottom
198           else if (sensN->GetZMin()-zTol>sens->GetZMin()) zType =  2;     // top
199           //
200           int phiType = 1;
201
202           double phiTstMn = sensN->GetPhiMin()-phiTol;
203           BringTo02Pi(phiTstMn);
204           if (!OKforPhiMax(sens->GetPhiMax(),phiTstMn)) continue; // too large angle      
205           double phiTstMx = sensN->GetPhiMax()+phiTol;    
206           BringTo02Pi(phiTstMx);
207           if (!OKforPhiMin(sens->GetPhiMin(),phiTstMx)) continue; // too large angle
208           //
209           phiTstMn = sensN->GetPhiMin()+phiTol;
210           BringTo02Pi(phiTstMn);
211           phiTstMx = sensN->GetPhiMax()-phiTol;   
212           BringTo02Pi(phiTstMx);
213           if      (!OKforPhiMax(sens->GetPhiMax(),phiTstMx)) phiType = 0; // left
214           else if (!OKforPhiMin(sens->GetPhiMin(),phiTstMn)) phiType = 2; // right
215           //
216           sens->SetNeighborID(kNBId[zType][phiType], neighbID);
217         } // phi scan
218       } // z scan
219     } // sensors
220   } // ladders
221   //
222 }
223
224 //______________________________________________________
225 Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[AliITSURecoSens::kNNeighbors])
226 {
227   // find sensors having intersection with track
228   // impPar contains: lab phi of track, dphi, labZ, dz
229   double z = impPar[2];
230   if (z>fZMax+impPar[3]) return 0; // outside of Z coverage
231   z -= fZMin;
232   if (z<-impPar[3]) return 0; // outside of Z coverage
233   int sensInLad = int(z*fSensDZInv);
234   if      (sensInLad<0) sensInLad = 0;
235   else if (sensInLad>=fNSensInLadder) sensInLad = fNSensInLadder-1;
236   //
237   double phi = impPar[0] - fPhiOffs;
238   BringTo02Pi(phi);
239   int ladID = int(phi*fDPhiLadInv);  // ladder id
240   int nsens = 0;
241   //
242   AliITSURecoSens* sensN,*sens = GetSensor(ladID*fNSensInLadder+sensInLad);
243   sensors[nsens++] = sens;
244   // check neighbours
245   double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3], phiMn=impPar[0]-impPar[1], phiMx=impPar[0]+impPar[1];
246   BringTo02Pi(phiMn);
247   BringTo02Pi(phiMx);
248   //
249   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbR)); // neighbor on the right (smaller phi)
250   if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax())) sensors[nsens++] = sensN;
251   //
252   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTR)); // neighbor on the top right (smaller phi, larger Z)
253   if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
254   //
255   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbT)); // neighbor on the top (larger Z)
256   if (sensN && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
257   //
258   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTL)); // neighbor on the top left (larger Z, larger phi)
259   if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
260   //
261   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbL)); // neighbor on the left (larger phi)
262   if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin())) sensors[nsens++] = sensN;
263   //
264   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBL)); // neighbor on the bottom left (smaller Z, larger phi)
265   if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
266   //
267   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbB));  // neighbor on the bottom (smaller Z)
268   if (sensN && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
269   //
270   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBR)); // neighbor on the bottom right (smaller Z, smaller phi)
271   if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
272   //
273   return nsens;
274 }
275
276 //______________________________________________________
277 void AliITSURecoLayer::ProcessClusters(Int_t mode)
278 {
279   // register in each sensor of the layer its cluster.
280   // the clusters of the layer must be sorted per sensor
281   int ncl = fClusters->GetEntriesFast();
282   int curSensID = -1;
283   for (int i=fNSensors;i--;) GetSensor(i)->SetNClusters(0);
284   AliITSURecoSens* curSens = 0;
285   for (int icl=0;icl<ncl;icl++) {
286     AliITSUClusterPix* cl = (AliITSUClusterPix*) fClusters->UncheckedAt(icl);
287     cl->GoToFrameTrk();
288     int vID = cl->GetVolumeId();
289     if (vID<curSensID) {AliFatal("Clusters are not sorted in increasing sensorID");}
290     if (vID>curSensID) {
291       if (curSens) curSens->ProcessClusters(mode);    // prepare clusters for reconstruction
292       curSens   = GetSensor(vID - fITSGeom->GetFirstModIndex(fActiveID));
293       curSensID = vID;
294       curSens->SetFirstClusterId(icl);
295     }
296     curSens->IncNClusters();
297   }
298   if (curSens) curSens->ProcessClusters(mode); // last sensor was not processed yet
299   //
300 }
301
302 //______________________________________________________
303 Bool_t AliITSURecoLayer::IsEqual(const TObject* obj) const
304 {
305   // check if layers are equal in R
306   const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
307   return Abs(lr->GetR()-GetR())<1e-6 ? kTRUE : kFALSE;
308 }
309
310 //______________________________________________________
311 Int_t  AliITSURecoLayer::Compare(const TObject* obj) const
312 {
313   // compare two layers
314   const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
315   double dr = GetR() - lr->GetR();
316   if (Abs(dr)<1e-6) return 0;
317   return dr>0 ? 1:-1;
318   //      
319 }
320
321 //_________________________________________________________________
322 AliITSURecoSens* AliITSURecoLayer::GetSensorFromID(Int_t i) const 
323 {
324   // get sensor from its global id
325   i -= fITSGeom->GetFirstModIndex(fActiveID);
326   if (i<0||i>=fNSensors) AliFatal(Form("Sensor with id=%d is not in layer %d",i+fITSGeom->GetFirstModIndex(fActiveID),fActiveID));
327   return (AliITSURecoSens*)fSensors[i];
328 }