In v1 geometry cluster should be attached to sensor according to sensor VIDs, rather...
[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   ,fNSensorRows(0)
20   ,fNSensorsPerRow(0)
21   ,fSensVIDtoMatrixID(0)
22   ,fR(0)
23   ,fRMax(0)
24   ,fRMin(0)
25   ,fZMax(0)
26   ,fZMin(0)
27   ,fPhiOffs(0)
28   ,fSensDZInv(0)
29   ,fSensDPhiInv(0)
30   ,fMaxStep(0.5)
31   ,fSensors(0)
32   ,fITSGeom(0)
33   ,fClusters(0)
34 {
35   // def. c-tor
36   SetNameTitle(name,name);
37 }
38
39 //______________________________________________________
40 AliITSURecoLayer::AliITSURecoLayer(const char* name, Int_t activeID, AliITSUGeomTGeo* gm)
41   :fActiveID(activeID)
42   ,fNSensors(0)
43   ,fNSensorRows(0)
44   ,fNSensorsPerRow(0)
45   ,fSensVIDtoMatrixID(0)
46   ,fR(0)
47   ,fRMax(0)
48   ,fRMin(0)
49   ,fZMax(0)
50   ,fZMin(0)
51   ,fPhiOffs(0)
52   ,fSensDZInv(0)
53   ,fSensDPhiInv(0)
54   ,fMaxStep(0.5)
55   ,fSensors(0)
56   ,fITSGeom(gm)
57   ,fClusters(0)
58 {
59   // def. c-tor
60   SetNameTitle(name,name);
61   Build();
62 }
63
64 //______________________________________________________
65 AliITSURecoLayer::~AliITSURecoLayer()
66 {
67   // def. d-tor
68   delete fSensors;
69   delete[] fSensVIDtoMatrixID;
70   if (GetOwnsClusterArray()) delete fClusters;
71 }
72
73 //______________________________________________________
74 void AliITSURecoLayer::Print(Option_t* opt) const                             
75 {
76   //print 
77   printf("Lr %-15s %d (act:%+d), NSens: %4d in %3d rows| MaxStep:%.2f ",
78          GetName(),GetID(),GetActiveID(),GetNSensors(),GetNSensorRows(),fMaxStep);
79   printf("%6.3f<R<%6.3f | %+8.3f<Z<%+8.3f dZ:%6.3f dPhi:%6.3f\n",fRMin,fRMax,fZMin,fZMax,
80          fSensDZInv>0 ? 1/fSensDZInv : 0, fSensDPhiInv>0 ? 1/fSensDPhiInv : 0);
81   TString opts = opt; opts.ToLower();
82   if (opts.Contains("sn")) for (int i=0;i<GetNSensors();i++) GetSensor(i)->Print(opt);
83 }
84
85 //______________________________________________________
86 void AliITSURecoLayer::Build()
87 {
88   // build internal structures
89   const double kSafeR = 0.05; // safety margin for Rmin,Rmax of the layer
90   if (fActiveID<0) return;
91   //
92   int nStaves=fITSGeom->GetNStaves(fActiveID);
93   // determine number of sensor rows (sensors aligned at same phi and spanning the Z range of the layer)
94   fNSensorRows = nStaves;
95   //
96   // if the stave has susbtaves, each substave can have multiple rows of sensors (but just 1 row of modules)
97   if (fITSGeom->GetNHalfStaves(fActiveID)>0) fNSensorRows *= fITSGeom->GetNHalfStaves(fActiveID);
98   //
99   // if there are modules defined, the module may have multiple rows of sensors (though not spanning full Z)
100   if (fITSGeom->GetNModules(fActiveID)>0) fNSensorRows *= fITSGeom->GetNChipRowsPerModule(fActiveID);
101   //
102   fNSensors = fITSGeom->GetNChipsPerLayer(fActiveID);
103   fNSensorsPerRow = fNSensors/fNSensorRows;
104   //
105   fSensors = new TObjArray(fNSensors);
106   fSensVIDtoMatrixID = new Int_t[fNSensors];
107   const AliITSsegmentation* kSegm = fITSGeom->GetSegmentation(fActiveID);
108   //
109   TGeoHMatrix mmod;
110   const TGeoHMatrix* mt2l;
111   double phiTF,rTF, loc[3]={0,0,0},glo[3];
112   //
113   int nSensPerStave = fITSGeom->GetNChipsPerStave(fActiveID);
114   for (int staveI=0;staveI<nStaves;staveI++) {
115     for (int sensI=0;sensI<nSensPerStave;sensI++) {
116       int sID = fITSGeom->GetChipIndex(fActiveID,staveI,sensI);
117       AliITSURecoSens* sens = new AliITSURecoSens( sID );
118       fSensors->AddLast(sens);
119       double phiMin=1e9,phiMax=-1e9,zMin=1e9,zMax=-1e9;
120       // this is NOT the sensor matrix, just the ideal chip matrix to get neighbors correct
121       fITSGeom->GetOrigMatrix(sID,mmod); 
122       //
123       for (int ix=0;ix<2;ix++) {       // determine sensor boundaries (ideal)
124         loc[0] = (ix-0.5)*kSegm->Dx(); // +-DX/2
125         for (int iy=0;iy<2;iy++) {
126           loc[1] = (iy-0.5)*kSegm->Dy(); // +-DY/2
127           for (int iz=0;iz<2;iz++) {
128             loc[2] = (iz-0.5)*kSegm->Dz(); // +-DZ/2
129             mmod.LocalToMaster(loc,glo);
130             double phi = ATan2(glo[1],glo[0]);
131             BringTo02Pi(phi);
132             if      (phiMin>1e8)  phiMin=phi;
133             else if (!OKforPhiMin(phiMin,phi)) phiMin=phi;
134             if      (phiMax<-1e8) phiMax=phi; 
135             else if (!OKforPhiMax(phiMax,phi)) phiMax=phi;            
136             if (glo[2]>zMax) zMax=glo[2];
137             if (glo[2]<zMin) zMin=glo[2];
138           }
139         }
140       }
141       sens->SetBoundaries(phiMin,phiMax,zMin,zMax);
142     }
143   }
144   fSensors->Sort(); // sort sensors to get the neighborhood correct
145   //
146   // now fill real sensor angles, Z's, accounting for misalignment
147   fRMin=fZMin=1e9;
148   fRMax=fZMax=-1e9;
149   //
150   fPhiOffs = 0;
151   int firstSensID = fITSGeom->GetFirstChipIndex(fActiveID);
152   for (int sensI=0;sensI<fNSensors;sensI++) {
153     AliITSURecoSens* sens = GetSensor(sensI);
154     mmod = *fITSGeom->GetMatrixSens(sens->GetID());
155     fSensVIDtoMatrixID[sens->GetID() - firstSensID] = sensI;
156     double phiMin=1e9,phiMax=-1e9,zMin=1e9,zMax=-1e9;
157     for (int ix=0;ix<2;ix++) {
158       loc[0] = (ix-0.5)*kSegm->Dx(); // +-DX/2
159       for (int iy=0;iy<2;iy++) {
160         loc[1] = (iy-0.5)*kSegm->Dy(); // +-DY/2
161         for (int iz=0;iz<2;iz++) {
162           loc[2] = (iz-0.5)*kSegm->Dz(); // +-DZ/2
163           //
164           mmod.LocalToMaster(loc,glo);
165           double phi = ATan2(glo[1],glo[0]);
166           double r   = glo[0]*glo[0] + glo[1]*glo[1];
167           if (fRMin>r) fRMin = r;
168           if (fRMax<r) fRMax = r;
169           BringTo02Pi(phi);
170           if      (phiMin>1e8) phiMin=phi; 
171           else if (!OKforPhiMin(phiMin,phi)) phiMin=phi;
172           if      (phiMax<-1e8) phiMax=phi;
173           else if (!OKforPhiMax(phiMax,phi)) phiMax=phi;              
174           if (glo[2]>zMax) zMax=glo[2];
175           if (glo[2]<zMin) zMin=glo[2];
176         }
177       }
178     }
179     mt2l = fITSGeom->GetMatrixT2L( sens->GetID() );
180     mmod.Multiply(mt2l);        
181     loc[0]=loc[1]=loc[2]=0;
182     mmod.LocalToMaster(loc,glo);
183     rTF   = Sqrt(glo[0]*glo[0] + glo[1]*glo[1]);  //  tracking params (misaligned)
184     phiTF = ATan2(glo[1],glo[0]);
185     BringTo02Pi(phiTF);
186     //
187     sens->SetXTF(rTF);
188     sens->SetPhiTF(phiTF);
189     sens->SetBoundaries(phiMin,phiMax,zMin,zMax);
190     if (fZMin>zMin) fZMin = zMin;
191     if (fZMax<zMax) fZMax = zMax;
192     //
193     if (sensI<fNSensorsPerRow) fPhiOffs += MeanPhiSmall(phiMax,phiMin);
194   }
195   //
196   fPhiOffs /= fNSensorsPerRow; // average phi of the 1st row
197   fSensDZInv = fNSensorsPerRow/(fZMax-fZMin);
198   fSensDPhiInv = fNSensorRows/(2*Pi());
199   //
200   fRMin = Sqrt(fRMin);
201   fRMax = Sqrt(fRMax);
202   fR = 0.5*(fRMin+fRMax);
203   fRMin -= kSafeR;
204   fRMax += kSafeR;
205   //
206 }
207
208 //______________________________________________________
209 Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[kMaxSensMatching])
210 {
211   // find sensors having intersection with track
212   // impPar contains: lab phi of track, dphi, labZ, dz
213   //
214   double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3]; 
215   if (zMn>fZMax) return 0;
216   if (zMx<fZMin) return 0;
217   //
218   int zCenID = int((impPar[2]-fZMin)*fSensDZInv);
219   if      (zCenID<0) zCenID = 0;
220   else if (zCenID>=fNSensorsPerRow) zCenID = fNSensorsPerRow-1;
221   double phiCn = impPar[0] - fPhiOffs;
222   //  BringTo02Pi(phiCn); 
223   int rowCenID = int(phiCn*fSensDPhiInv);
224   //
225   // due to the misalignments the actual sensorID's might be shifted
226   int res = 0;
227   AliITSURecoSens* sensPrev=0, *sens = GetSensor(rowCenID,zCenID);
228   //
229   //  printf("Guess: Primary Sensor: phiID: %d zID: %d ->",rowCenID,zCenID); sens->Print();
230   //
231   while ( (res=sens->CheckCoverage(impPar[0], impPar[2])) ) {
232     if      (res&AliITSURecoSens::kRight) {if (++rowCenID==fNSensorRows) rowCenID=0;} // neighbor on the right (larger phi)
233     else if (res&AliITSURecoSens::kLeft)  {if (--rowCenID<0) rowCenID = fNSensorRows-1;}      // neighbor on the left (smaller phi)
234     if      (res&AliITSURecoSens::kUp)    {if (++zCenID==fNSensorsPerRow) zCenID = fNSensorsPerRow-1;}    // neighbor at larger Z (if any)
235     else if (res&AliITSURecoSens::kDown)  {if (--zCenID<0) zCenID = 0;}               // neighbor at smaller Z (if any)
236     //
237     AliITSURecoSens* sensAlt = GetSensor(rowCenID,zCenID);
238     if (sensAlt==sens || sensAlt==sensPrev) break;  // there is no better neighbor (z edge) or the point falls in dead area
239     //
240     sensPrev = sens;
241     sens = sensAlt;
242   }
243   //  printf("Found: Primary Sensor: phiID: %d zID: %d ->",rowCenID,zCenID); sens->Print();
244   //
245   int nFnd = 0;
246   sensors[nFnd++] = sens;
247   //
248   double phiMn = impPar[0]-impPar[1], phiMx = impPar[0]+impPar[1];
249   BringTo02Pi(phiMn);
250   BringTo02Pi(phiMx);
251   //
252   const int kNNeighb = 8;
253   const int kCheckNeighb[2][kNNeighb] = { // phi and Z neighbours to check
254     { 1, 1, 0,-1,-1,-1, 0, 1},
255     { 0, 1, 1, 1, 0,-1,-1,-1}
256   };
257   //  
258   //  printf("Search: %+.4f %+.4f | %+.4f %+.4f\n",phiMn,phiMx, zMn,zMx);
259   for (int inb=kNNeighb;inb--;) {
260     int idz = kCheckNeighb[1][inb];
261     int iz   = zCenID   + idz;
262     //    printf("#%d  dp:%+d dz:%+d IZ: %d\n",inb, kCheckNeighb[0][inb], kCheckNeighb[1][inb], iz);
263     //
264     if (iz<0 || iz>=fNSensorsPerRow) continue;
265     int idphi = kCheckNeighb[0][inb];
266     int iphi = rowCenID + idphi;
267     if      (iphi<0) iphi += fNSensorRows;
268     else if (iphi>=fNSensorRows) iphi -= fNSensorRows;
269     sens = GetSensor(iphi,iz);
270     //
271     if      (idz>0) {if (zMx<sens->GetZMin()) continue;}
272     else if (idz<0) {if (zMn>sens->GetZMax()) continue;}
273     //
274     // Z range matches
275     if      (idphi>0) {if (!OKforPhiMin(sens->GetPhiMin(),phiMx)) continue;}
276     else if (idphi<0) {if (!OKforPhiMax(sens->GetPhiMax(),phiMn)) continue;}
277     //
278     //    printf("Add %d\n",nFnd);
279     sensors[nFnd++] = sens;
280     if (nFnd==kMaxSensMatching) break;
281   }
282   return nFnd;
283 }
284
285 //*/
286 /*
287 Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[AliITSURecoSens::kNNeighbors],int mcLab)
288 {
289   // find sensors having intersection with track
290   // impPar contains: lab phi of track, dphi, labZ, dz
291
292   //tmp>>>
293   int nFnd = 0;
294   int fndSens[50];
295   if (mcLab>=0) { // find correct sensors from MC info
296     int ncl = GetNClusters();
297     for (int icl=ncl;icl--;) {
298       AliCluster* cl = GetCluster(icl);
299       for (int ilb=0;ilb<3;ilb++) {
300         if (cl->GetLabel(ilb)<0) break;
301         if (cl->GetLabel(ilb)==mcLab) {fndSens[nFnd++] = cl->GetVolumeId(); break;}
302       }
303     }
304     if (nFnd>0) {
305       Int_t layS,staS,sensS;
306       for (int is=0;is<nFnd;is++) {
307         fITSGeom->GetChipId(fndSens[is],layS,staS,sensS);
308         printf("SNMC#%d(%d): %d %d %d | ",is,mcLab,layS,staS,sensS); GetSensorFromID(fndSens[is])->Print();
309       }
310     }
311   }
312
313   //tmp<<<
314
315   double z = impPar[2];
316   if (z>fZMax+impPar[3]) {
317     if (nFnd>0) printf("MissedSens!!!\n");
318     return 0; // outside of Z coverage
319   }
320   z -= fZMin;
321   if (z<-impPar[3]) {
322     if (nFnd>0) printf("MissedSens!!!\n");
323     return 0; // outside of Z coverage
324   }
325   int sensInSta = int(z*fSensDZInv);
326   if      (sensInSta<0) sensInSta = 0;
327   else if (sensInSta>=fNSensInStave) sensInSta = fNSensInStave-1;
328   //
329   double phi = impPar[0] - fPhiOffs;
330   BringTo02Pi(phi);
331   int staID = int(phi*fSensDPhiInv);  // stave id
332   int nsens = 0;
333   //
334   AliITSURecoSens* sensN,*sens = GetSensor(staID*fNSensInStave+sensInSta);
335   sensors[nsens++] = sens;
336   //
337   // check neighbours
338   double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3], phiMn=impPar[0]-impPar[1], phiMx=impPar[0]+impPar[1];
339   BringTo02Pi(phiMn);
340   BringTo02Pi(phiMx);
341   //
342   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbR)); // neighbor on the right (smaller phi)
343   if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax())) sensors[nsens++] = sensN;
344   //
345   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTR)); // neighbor on the top right (smaller phi, larger Z)
346   if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
347   //
348   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbT)); // neighbor on the top (larger Z)
349   if (sensN && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
350   //
351   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTL)); // neighbor on the top left (larger Z, larger phi)
352   if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
353   //
354   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbL)); // neighbor on the left (larger phi)
355   if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin())) sensors[nsens++] = sensN;
356   //
357   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBL)); // neighbor on the bottom left (smaller Z, larger phi)
358   if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
359   //
360   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbB));  // neighbor on the bottom (smaller Z)
361   if (sensN && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
362   //
363   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBR)); // neighbor on the bottom right (smaller Z, smaller phi)
364   if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
365   //
366   if (mcLab>=0) {
367     Int_t layS,staS,sensS;
368     printf("Found %d sensors for phi %.3f : %.3f | Z %.4f %.4f\n", nsens,phiMn,phiMx,zMn,zMx); 
369     for (int is=0;is<nsens;is++) {
370       fITSGeom->GetChipId(sensors[is]->GetID()+fITSGeom->GetFirstModIndex(fActiveID),layS,staS,sensS);
371       printf("*SNF#%d: %d %d %d | ",is,layS,staS,sensS); sensors[is]->Print();
372     }
373     for (int ism=0;ism<nFnd;ism++) {
374       AliITSURecoSens* snMC = GetSensorFromID(fndSens[ism]);
375       Bool_t ok=kFALSE;
376       for (int isf=0;isf<nsens;isf++) {
377         if (snMC==sensors[isf]) {ok=kTRUE;break;}
378       }
379       if (!ok) printf("MissedSens %d!!!\n",ism);
380     }
381   }
382   return nsens;
383 }
384 */
385 //______________________________________________________
386 void AliITSURecoLayer::ProcessClusters(Int_t mode)
387 {
388   // register in each sensor of the layer its cluster.
389   // the clusters of the layer must be sorted per sensor
390   int ncl = fClusters->GetEntriesFast();
391   int curSensID = -1;
392   for (int i=fNSensors;i--;) GetSensor(i)->SetNClusters(0);
393   AliITSURecoSens* curSens = 0;
394   for (int icl=0;icl<ncl;icl++) {
395     AliITSUClusterPix* cl = (AliITSUClusterPix*) fClusters->UncheckedAt(icl);
396     cl->GoToFrameTrk();
397     int vID = cl->GetVolumeId();
398     if (vID<curSensID) {AliFatal("Clusters are not sorted in increasing sensorID");}
399     if (vID>curSensID) {
400       if (curSens) curSens->ProcessClusters(mode);    // prepare clusters for reconstruction
401       curSens   = GetSensorFromID(vID); //GetSensor(vID - fITSGeom->GetFirstChipIndex(fActiveID));
402       curSensID = vID;
403       curSens->SetFirstClusterId(icl);
404     }
405     curSens->IncNClusters();
406   }
407   if (curSens) curSens->ProcessClusters(mode); // last sensor was not processed yet
408   //
409 }
410
411 //______________________________________________________
412 Bool_t AliITSURecoLayer::IsEqual(const TObject* obj) const
413 {
414   // check if layers are equal in R
415   const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
416   return Abs(lr->GetR()-GetR())<1e-6 ? kTRUE : kFALSE;
417 }
418
419 //______________________________________________________
420 Int_t  AliITSURecoLayer::Compare(const TObject* obj) const
421 {
422   // compare two layers
423   const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
424   double dr = GetR() - lr->GetR();
425   if (Abs(dr)<1e-6) return 0;
426   return dr>0 ? 1:-1;
427   //      
428 }
429
430 //_________________________________________________________________
431 AliITSURecoSens* AliITSURecoLayer::GetSensorFromID(Int_t i) const 
432 {
433   // get sensor from its global id
434   i -= fITSGeom->GetFirstChipIndex(fActiveID);
435   if (i<0||i>=fNSensors) AliFatal(Form("Sensor with id=%d is not in layer %d",i+fITSGeom->GetFirstChipIndex(fActiveID),fActiveID));
436   return GetSensor(SensVIDtoMatrixID(i));
437 }