Integrating the Cooked Matrix tracker into the commom reconstruction framework
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSURecoLayer.cxx
CommitLineData
32d38de2 1#include <TClonesArray.h>
a11ef2e4 2#include "AliITSURecoLayer.h"
a11ef2e4 3#include "AliITSsegmentation.h"
a11ef2e4 4#include "AliITSUAux.h"
32d38de2 5#include "AliITSUClusterPix.h"
8b16dbae 6#include "AliITSUGeomTGeo.h"
7#include "AliLog.h"
a11ef2e4 8
9using namespace AliITSUAux;
10using namespace TMath;
11
12ClassImp(AliITSURecoLayer)
13
14
15//______________________________________________________
16AliITSURecoLayer::AliITSURecoLayer(const char* name)
68b7631d 17 :fActiveID(-1)
a11ef2e4 18 ,fNSensors(0)
68b7631d 19 ,fNSensorRows(0)
20 ,fNSensorsPerRow(0)
21 ,fSensVIDtoMatrixID(0)
a11ef2e4 22 ,fR(0)
23 ,fRMax(0)
24 ,fRMin(0)
25 ,fZMax(0)
26 ,fZMin(0)
a11ef2e4 27 ,fPhiOffs(0)
28 ,fSensDZInv(0)
68b7631d 29 ,fSensDPhiInv(0)
32d38de2 30 ,fMaxStep(0.5)
a11ef2e4 31 ,fSensors(0)
32 ,fITSGeom(0)
32d38de2 33 ,fClusters(0)
a11ef2e4 34{
35 // def. c-tor
36 SetNameTitle(name,name);
a11ef2e4 37}
38
39//______________________________________________________
32d38de2 40AliITSURecoLayer::AliITSURecoLayer(const char* name, Int_t activeID, AliITSUGeomTGeo* gm)
a11ef2e4 41 :fActiveID(activeID)
42 ,fNSensors(0)
68b7631d 43 ,fNSensorRows(0)
44 ,fNSensorsPerRow(0)
45 ,fSensVIDtoMatrixID(0)
a11ef2e4 46 ,fR(0)
47 ,fRMax(0)
48 ,fRMin(0)
49 ,fZMax(0)
50 ,fZMin(0)
a11ef2e4 51 ,fPhiOffs(0)
52 ,fSensDZInv(0)
68b7631d 53 ,fSensDPhiInv(0)
32d38de2 54 ,fMaxStep(0.5)
55 ,fSensors(0)
a11ef2e4 56 ,fITSGeom(gm)
32d38de2 57 ,fClusters(0)
a11ef2e4 58{
59 // def. c-tor
60 SetNameTitle(name,name);
32d38de2 61 Build();
a11ef2e4 62}
63
64//______________________________________________________
65AliITSURecoLayer::~AliITSURecoLayer()
66{
67 // def. d-tor
68b7631d 68 delete fSensors;
69 delete[] fSensVIDtoMatrixID;
ee52e7b5 70 if (GetOwnsClusterArray()) delete fClusters;
a11ef2e4 71}
72
73//______________________________________________________
74void AliITSURecoLayer::Print(Option_t* opt) const
75{
76 //print
68b7631d 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);
a11ef2e4 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//______________________________________________________
32d38de2 86void AliITSURecoLayer::Build()
a11ef2e4 87{
88 // build internal structures
b515ad7e 89 const double kSafeR = 0.05; // safety margin for Rmin,Rmax of the layer
32d38de2 90 if (fActiveID<0) return;
68b7631d 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];
546d00d8 107 const AliITSsegmentation* kSegm = fITSGeom->GetSegmentation(fActiveID);
a11ef2e4 108 //
a11ef2e4 109 TGeoHMatrix mmod;
110 const TGeoHMatrix* mt2l;
a11ef2e4 111 double phiTF,rTF, loc[3]={0,0,0},glo[3];
68b7631d 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);
a11ef2e4 119 double phiMin=1e9,phiMax=-1e9,zMin=1e9,zMax=-1e9;
68b7631d 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)
546d00d8 124 loc[0] = (ix-0.5)*kSegm->Dx(); // +-DX/2
a11ef2e4 125 for (int iy=0;iy<2;iy++) {
546d00d8 126 loc[1] = (iy-0.5)*kSegm->Dy(); // +-DY/2
a11ef2e4 127 for (int iz=0;iz<2;iz++) {
546d00d8 128 loc[2] = (iz-0.5)*kSegm->Dz(); // +-DZ/2
a11ef2e4 129 mmod.LocalToMaster(loc,glo);
130 double phi = ATan2(glo[1],glo[0]);
a11ef2e4 131 BringTo02Pi(phi);
68b7631d 132 if (phiMin>1e8) phiMin=phi;
a11ef2e4 133 else if (!OKforPhiMin(phiMin,phi)) phiMin=phi;
68b7631d 134 if (phiMax<-1e8) phiMax=phi;
a11ef2e4 135 else if (!OKforPhiMax(phiMax,phi)) phiMax=phi;
dde91d5d 136 if (glo[2]>zMax) zMax=glo[2];
137 if (glo[2]<zMin) zMin=glo[2];
a11ef2e4 138 }
139 }
140 }
a11ef2e4 141 sens->SetBoundaries(phiMin,phiMax,zMin,zMax);
a11ef2e4 142 }
143 }
68b7631d 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());
a11ef2e4 199 //
200 fRMin = Sqrt(fRMin);
201 fRMax = Sqrt(fRMax);
202 fR = 0.5*(fRMin+fRMax);
b515ad7e 203 fRMin -= kSafeR;
204 fRMax += kSafeR;
a11ef2e4 205 //
a11ef2e4 206}
08419930 207
a11ef2e4 208//______________________________________________________
68b7631d 209Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[kMaxSensMatching])
a11ef2e4 210{
211 // find sensors having intersection with track
212 // impPar contains: lab phi of track, dphi, labZ, dz
5785a9d8 213 //
68b7631d 214 double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3];
215 if (zMn>fZMax) return 0;
216 if (zMx<fZMin) return 0;
a11ef2e4 217 //
68b7631d 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);
a11ef2e4 224 //
68b7631d 225 // due to the misalignments the actual sensorID's might be shifted
226 int res = 0;
227 AliITSURecoSens* sensPrev=0, *sens = GetSensor(rowCenID,zCenID);
5785a9d8 228 //
68b7631d 229 // printf("Guess: Primary Sensor: phiID: %d zID: %d ->",rowCenID,zCenID); sens->Print();
a11ef2e4 230 //
68b7631d 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();
a11ef2e4 244 //
68b7631d 245 int nFnd = 0;
246 sensors[nFnd++] = sens;
a11ef2e4 247 //
68b7631d 248 double phiMn = impPar[0]-impPar[1], phiMx = impPar[0]+impPar[1];
249 BringTo02Pi(phiMn);
250 BringTo02Pi(phiMx);
a11ef2e4 251 //
68b7631d 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;
a11ef2e4 283}
68b7631d 284
c03e4f8a 285//*/
5785a9d8 286/*
287Int_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) {
852af72e 305 Int_t layS,staS,sensS;
5785a9d8 306 for (int is=0;is<nFnd;is++) {
852af72e 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();
5785a9d8 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 }
852af72e 325 int sensInSta = int(z*fSensDZInv);
326 if (sensInSta<0) sensInSta = 0;
327 else if (sensInSta>=fNSensInStave) sensInSta = fNSensInStave-1;
5785a9d8 328 //
329 double phi = impPar[0] - fPhiOffs;
330 BringTo02Pi(phi);
68b7631d 331 int staID = int(phi*fSensDPhiInv); // stave id
5785a9d8 332 int nsens = 0;
333 //
852af72e 334 AliITSURecoSens* sensN,*sens = GetSensor(staID*fNSensInStave+sensInSta);
5785a9d8 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) {
852af72e 367 Int_t layS,staS,sensS;
5785a9d8 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++) {
852af72e 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();
5785a9d8 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*/
32d38de2 385//______________________________________________________
386void 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;
0091e9f0 392 for (int i=fNSensors;i--;) GetSensor(i)->SetNClusters(0);
32d38de2 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
6a6a65db 401 curSens = GetSensorFromID(vID); //GetSensor(vID - fITSGeom->GetFirstChipIndex(fActiveID));
32d38de2 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}
dde91d5d 410
411//______________________________________________________
412Bool_t AliITSURecoLayer::IsEqual(const TObject* obj) const
413{
414 // check if layers are equal in R
44785f3e 415 const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
dde91d5d 416 return Abs(lr->GetR()-GetR())<1e-6 ? kTRUE : kFALSE;
417}
418
419//______________________________________________________
420Int_t AliITSURecoLayer::Compare(const TObject* obj) const
421{
422 // compare two layers
44785f3e 423 const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
dde91d5d 424 double dr = GetR() - lr->GetR();
425 if (Abs(dr)<1e-6) return 0;
426 return dr>0 ? 1:-1;
427 //
428}
8b16dbae 429
430//_________________________________________________________________
431AliITSURecoSens* AliITSURecoLayer::GetSensorFromID(Int_t i) const
432{
433 // get sensor from its global id
852af72e 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));
68b7631d 436 return GetSensor(SensVIDtoMatrixID(i));
8b16dbae 437}