1 #include <TClonesArray.h>
2 #include "AliITSURecoLayer.h"
3 #include "AliITSsegmentation.h"
4 #include "AliITSUAux.h"
5 #include "AliITSUClusterPix.h"
6 #include "AliITSUGeomTGeo.h"
9 using namespace AliITSUAux;
10 using namespace TMath;
12 ClassImp(AliITSURecoLayer)
15 //______________________________________________________
16 AliITSURecoLayer::AliITSURecoLayer(const char* name)
21 ,fSensVIDtoMatrixID(0)
36 SetNameTitle(name,name);
39 //______________________________________________________
40 AliITSURecoLayer::AliITSURecoLayer(const char* name, Int_t activeID, AliITSUGeomTGeo* gm)
45 ,fSensVIDtoMatrixID(0)
60 SetNameTitle(name,name);
64 //______________________________________________________
65 AliITSURecoLayer::~AliITSURecoLayer()
69 delete[] fSensVIDtoMatrixID;
70 if (GetOwnsClusterArray()) delete fClusters;
73 //______________________________________________________
74 void AliITSURecoLayer::Print(Option_t* opt) const
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);
85 //______________________________________________________
86 void AliITSURecoLayer::Build()
88 // build internal structures
89 const double kSafeR = 0.05; // safety margin for Rmin,Rmax of the layer
90 if (fActiveID<0) return;
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;
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);
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);
102 fNSensors = fITSGeom->GetNChipsPerLayer(fActiveID);
103 fNSensorsPerRow = fNSensors/fNSensorRows;
105 fSensors = new TObjArray(fNSensors);
106 fSensVIDtoMatrixID = new Int_t[fNSensors];
107 const AliITSsegmentation* kSegm = fITSGeom->GetSegmentation(fActiveID);
110 const TGeoHMatrix* mt2l;
111 double phiTF,rTF, loc[3]={0,0,0},glo[3];
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);
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]);
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];
141 sens->SetBoundaries(phiMin,phiMax,zMin,zMax);
144 fSensors->Sort(); // sort sensors to get the neighborhood correct
146 // now fill real sensor angles, Z's, accounting for misalignment
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
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;
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];
179 mt2l = fITSGeom->GetMatrixT2L( sens->GetID() );
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]);
188 sens->SetPhiTF(phiTF);
189 sens->SetBoundaries(phiMin,phiMax,zMin,zMax);
190 if (fZMin>zMin) fZMin = zMin;
191 if (fZMax<zMax) fZMax = zMax;
193 if (sensI<fNSensorsPerRow) fPhiOffs += MeanPhiSmall(phiMax,phiMin);
196 fPhiOffs /= fNSensorsPerRow; // average phi of the 1st row
197 fSensDZInv = fNSensorsPerRow/(fZMax-fZMin);
198 fSensDPhiInv = fNSensorRows/(2*Pi());
202 fR = 0.5*(fRMin+fRMax);
208 //______________________________________________________
209 Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[kMaxSensMatching])
211 // find sensors having intersection with track
212 // impPar contains: lab phi of track, dphi, labZ, dz
214 double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3];
215 if (zMn>fZMax) return 0;
216 if (zMx<fZMin) return 0;
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);
225 // due to the misalignments the actual sensorID's might be shifted
227 AliITSURecoSens* sensPrev=0, *sens = GetSensor(rowCenID,zCenID);
229 // printf("Guess: Primary Sensor: phiID: %d zID: %d ->",rowCenID,zCenID); sens->Print();
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)
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
243 // printf("Found: Primary Sensor: phiID: %d zID: %d ->",rowCenID,zCenID); sens->Print();
246 sensors[nFnd++] = sens;
248 double phiMn = impPar[0]-impPar[1], phiMx = impPar[0]+impPar[1];
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}
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);
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);
271 if (idz>0) {if (zMx<sens->GetZMin()) continue;}
272 else if (idz<0) {if (zMn>sens->GetZMax()) continue;}
275 if (idphi>0) {if (!OKforPhiMin(sens->GetPhiMin(),phiMx)) continue;}
276 else if (idphi<0) {if (!OKforPhiMax(sens->GetPhiMax(),phiMn)) continue;}
278 // printf("Add %d\n",nFnd);
279 sensors[nFnd++] = sens;
280 if (nFnd==kMaxSensMatching) break;
287 Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[AliITSURecoSens::kNNeighbors],int mcLab)
289 // find sensors having intersection with track
290 // impPar contains: lab phi of track, dphi, labZ, dz
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;}
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();
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
322 if (nFnd>0) printf("MissedSens!!!\n");
323 return 0; // outside of Z coverage
325 int sensInSta = int(z*fSensDZInv);
326 if (sensInSta<0) sensInSta = 0;
327 else if (sensInSta>=fNSensInStave) sensInSta = fNSensInStave-1;
329 double phi = impPar[0] - fPhiOffs;
331 int staID = int(phi*fSensDPhiInv); // stave id
334 AliITSURecoSens* sensN,*sens = GetSensor(staID*fNSensInStave+sensInSta);
335 sensors[nsens++] = sens;
338 double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3], phiMn=impPar[0]-impPar[1], phiMx=impPar[0]+impPar[1];
342 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbR)); // neighbor on the right (smaller phi)
343 if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax())) sensors[nsens++] = sensN;
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;
348 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbT)); // neighbor on the top (larger Z)
349 if (sensN && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
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;
354 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbL)); // neighbor on the left (larger phi)
355 if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin())) sensors[nsens++] = sensN;
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;
360 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbB)); // neighbor on the bottom (smaller Z)
361 if (sensN && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
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;
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();
373 for (int ism=0;ism<nFnd;ism++) {
374 AliITSURecoSens* snMC = GetSensorFromID(fndSens[ism]);
376 for (int isf=0;isf<nsens;isf++) {
377 if (snMC==sensors[isf]) {ok=kTRUE;break;}
379 if (!ok) printf("MissedSens %d!!!\n",ism);
385 //______________________________________________________
386 void AliITSURecoLayer::ProcessClusters(Int_t mode)
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();
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);
397 int vID = cl->GetVolumeId();
398 if (vID<curSensID) {AliFatal("Clusters are not sorted in increasing sensorID");}
400 if (curSens) curSens->ProcessClusters(mode); // prepare clusters for reconstruction
401 curSens = GetSensorFromID(vID); //GetSensor(vID - fITSGeom->GetFirstChipIndex(fActiveID));
403 curSens->SetFirstClusterId(icl);
405 curSens->IncNClusters();
407 if (curSens) curSens->ProcessClusters(mode); // last sensor was not processed yet
411 //______________________________________________________
412 Bool_t AliITSURecoLayer::IsEqual(const TObject* obj) const
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;
419 //______________________________________________________
420 Int_t AliITSURecoLayer::Compare(const TObject* obj) const
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;
430 //_________________________________________________________________
431 AliITSURecoSens* AliITSURecoLayer::GetSensorFromID(Int_t i) const
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));