modifs for ITSU v1
[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)
17 :fActiveID(-1)
18 ,fNSensors(0)
5009207f 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)
5009207f 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)
5009207f 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)
5009207f 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
5009207f 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
5009207f 77 printf("Lr %-15s %d (act:%+d), NSens: %4d in %d 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;
5009207f 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];
5009207f 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;
5009207f 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);
5009207f 132 if (phiMin>1e8) phiMin=phi;
a11ef2e4 133 else if (!OKforPhiMin(phiMin,phi)) phiMin=phi;
5009207f 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 }
5009207f 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//______________________________________________________
5009207f 209Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[kNNeighbors])
a11ef2e4 210{
211 // find sensors having intersection with track
212 // impPar contains: lab phi of track, dphi, labZ, dz
5785a9d8 213 //
5009207f 214 return 0;
215 /*
a11ef2e4 216 double z = impPar[2];
0091e9f0 217 if (z>fZMax+impPar[3]) return 0; // outside of Z coverage
a11ef2e4 218 z -= fZMin;
0091e9f0 219 if (z<-impPar[3]) return 0; // outside of Z coverage
5009207f 220 int sensInRow = int(z*fSensDZInv);
221 if (sensInRow<0) sensInRow = 0;
222 else if (sensInRow>=fNSensorsPerRow) sensInRow = fNSensorsPerRow-1;
a11ef2e4 223 //
224 double phi = impPar[0] - fPhiOffs;
225 BringTo02Pi(phi);
5009207f 226 int rowID = int(phi*fSensDPhiInv); // stave id
227 //
228 int sensID = rowID*fNSensorsPerRow + sensInRow; // most probable candidate
a11ef2e4 229 int nsens = 0;
230 //
5009207f 231 AliITSURecoSens* sensN,*sens = GetSensor(sesnID);
232 //
19ea9736 233 int res = sens->CheckCoverage(impPar);
234
5009207f 235 // make sure this is best matching sensor
236 if (sens->GetZMin()<impPar[2] && sens->GetZMax()>impPar[2] &&
237 OKforPhiMin(sens->GetPhiMin(),impPar[0]) && OKforPhiMax(sens->GetPhiMax(),impPar[0]) )
238
a11ef2e4 239 sensors[nsens++] = sens;
5785a9d8 240 //
a11ef2e4 241 // check neighbours
242 double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3], phiMn=impPar[0]-impPar[1], phiMx=impPar[0]+impPar[1];
243 BringTo02Pi(phiMn);
244 BringTo02Pi(phiMx);
245 //
173b3073 246 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbR)); // neighbor on the right (smaller phi)
247 if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax())) sensors[nsens++] = sensN;
a11ef2e4 248 //
173b3073 249 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTR)); // neighbor on the top right (smaller phi, larger Z)
250 if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
a11ef2e4 251 //
173b3073 252 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbT)); // neighbor on the top (larger Z)
253 if (sensN && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
a11ef2e4 254 //
173b3073 255 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTL)); // neighbor on the top left (larger Z, larger phi)
256 if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
a11ef2e4 257 //
173b3073 258 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbL)); // neighbor on the left (larger phi)
259 if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin())) sensors[nsens++] = sensN;
a11ef2e4 260 //
173b3073 261 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBL)); // neighbor on the bottom left (smaller Z, larger phi)
262 if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
a11ef2e4 263 //
173b3073 264 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbB)); // neighbor on the bottom (smaller Z)
a11ef2e4 265 if (sensN && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
266 //
173b3073 267 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBR)); // neighbor on the bottom right (smaller Z, smaller phi)
268 if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
a11ef2e4 269 //
270 return nsens;
5009207f 271 */
a11ef2e4 272}
5009207f 273
c03e4f8a 274//*/
5785a9d8 275/*
276Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[AliITSURecoSens::kNNeighbors],int mcLab)
277{
278 // find sensors having intersection with track
279 // impPar contains: lab phi of track, dphi, labZ, dz
280
281 //tmp>>>
282 int nFnd = 0;
283 int fndSens[50];
284 if (mcLab>=0) { // find correct sensors from MC info
285 int ncl = GetNClusters();
286 for (int icl=ncl;icl--;) {
287 AliCluster* cl = GetCluster(icl);
288 for (int ilb=0;ilb<3;ilb++) {
289 if (cl->GetLabel(ilb)<0) break;
290 if (cl->GetLabel(ilb)==mcLab) {fndSens[nFnd++] = cl->GetVolumeId(); break;}
291 }
292 }
293 if (nFnd>0) {
852af72e 294 Int_t layS,staS,sensS;
5785a9d8 295 for (int is=0;is<nFnd;is++) {
852af72e 296 fITSGeom->GetChipId(fndSens[is],layS,staS,sensS);
297 printf("SNMC#%d(%d): %d %d %d | ",is,mcLab,layS,staS,sensS); GetSensorFromID(fndSens[is])->Print();
5785a9d8 298 }
299 }
300 }
301
302 //tmp<<<
303
304 double z = impPar[2];
305 if (z>fZMax+impPar[3]) {
306 if (nFnd>0) printf("MissedSens!!!\n");
307 return 0; // outside of Z coverage
308 }
309 z -= fZMin;
310 if (z<-impPar[3]) {
311 if (nFnd>0) printf("MissedSens!!!\n");
312 return 0; // outside of Z coverage
313 }
852af72e 314 int sensInSta = int(z*fSensDZInv);
315 if (sensInSta<0) sensInSta = 0;
316 else if (sensInSta>=fNSensInStave) sensInSta = fNSensInStave-1;
5785a9d8 317 //
318 double phi = impPar[0] - fPhiOffs;
319 BringTo02Pi(phi);
5009207f 320 int staID = int(phi*fSensDPhiInv); // stave id
5785a9d8 321 int nsens = 0;
322 //
852af72e 323 AliITSURecoSens* sensN,*sens = GetSensor(staID*fNSensInStave+sensInSta);
5785a9d8 324 sensors[nsens++] = sens;
325 //
326 // check neighbours
327 double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3], phiMn=impPar[0]-impPar[1], phiMx=impPar[0]+impPar[1];
328 BringTo02Pi(phiMn);
329 BringTo02Pi(phiMx);
330 //
331 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbR)); // neighbor on the right (smaller phi)
332 if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax())) sensors[nsens++] = sensN;
333 //
334 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTR)); // neighbor on the top right (smaller phi, larger Z)
335 if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
336 //
337 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbT)); // neighbor on the top (larger Z)
338 if (sensN && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
339 //
340 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTL)); // neighbor on the top left (larger Z, larger phi)
341 if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
342 //
343 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbL)); // neighbor on the left (larger phi)
344 if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin())) sensors[nsens++] = sensN;
345 //
346 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBL)); // neighbor on the bottom left (smaller Z, larger phi)
347 if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
348 //
349 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbB)); // neighbor on the bottom (smaller Z)
350 if (sensN && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
351 //
352 sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBR)); // neighbor on the bottom right (smaller Z, smaller phi)
353 if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
354 //
355 if (mcLab>=0) {
852af72e 356 Int_t layS,staS,sensS;
5785a9d8 357 printf("Found %d sensors for phi %.3f : %.3f | Z %.4f %.4f\n", nsens,phiMn,phiMx,zMn,zMx);
358 for (int is=0;is<nsens;is++) {
852af72e 359 fITSGeom->GetChipId(sensors[is]->GetID()+fITSGeom->GetFirstModIndex(fActiveID),layS,staS,sensS);
360 printf("*SNF#%d: %d %d %d | ",is,layS,staS,sensS); sensors[is]->Print();
5785a9d8 361 }
362 for (int ism=0;ism<nFnd;ism++) {
363 AliITSURecoSens* snMC = GetSensorFromID(fndSens[ism]);
364 Bool_t ok=kFALSE;
365 for (int isf=0;isf<nsens;isf++) {
366 if (snMC==sensors[isf]) {ok=kTRUE;break;}
367 }
368 if (!ok) printf("MissedSens %d!!!\n",ism);
369 }
370 }
371 return nsens;
372}
373*/
32d38de2 374//______________________________________________________
375void AliITSURecoLayer::ProcessClusters(Int_t mode)
376{
377 // register in each sensor of the layer its cluster.
378 // the clusters of the layer must be sorted per sensor
379 int ncl = fClusters->GetEntriesFast();
380 int curSensID = -1;
0091e9f0 381 for (int i=fNSensors;i--;) GetSensor(i)->SetNClusters(0);
32d38de2 382 AliITSURecoSens* curSens = 0;
383 for (int icl=0;icl<ncl;icl++) {
384 AliITSUClusterPix* cl = (AliITSUClusterPix*) fClusters->UncheckedAt(icl);
385 cl->GoToFrameTrk();
386 int vID = cl->GetVolumeId();
387 if (vID<curSensID) {AliFatal("Clusters are not sorted in increasing sensorID");}
388 if (vID>curSensID) {
389 if (curSens) curSens->ProcessClusters(mode); // prepare clusters for reconstruction
852af72e 390 curSens = GetSensor(vID - fITSGeom->GetFirstChipIndex(fActiveID));
32d38de2 391 curSensID = vID;
392 curSens->SetFirstClusterId(icl);
393 }
394 curSens->IncNClusters();
395 }
396 if (curSens) curSens->ProcessClusters(mode); // last sensor was not processed yet
397 //
398}
dde91d5d 399
400//______________________________________________________
401Bool_t AliITSURecoLayer::IsEqual(const TObject* obj) const
402{
403 // check if layers are equal in R
44785f3e 404 const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
dde91d5d 405 return Abs(lr->GetR()-GetR())<1e-6 ? kTRUE : kFALSE;
406}
407
408//______________________________________________________
409Int_t AliITSURecoLayer::Compare(const TObject* obj) const
410{
411 // compare two layers
44785f3e 412 const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
dde91d5d 413 double dr = GetR() - lr->GetR();
414 if (Abs(dr)<1e-6) return 0;
415 return dr>0 ? 1:-1;
416 //
417}
8b16dbae 418
419//_________________________________________________________________
420AliITSURecoSens* AliITSURecoLayer::GetSensorFromID(Int_t i) const
421{
422 // get sensor from its global id
852af72e 423 i -= fITSGeom->GetFirstChipIndex(fActiveID);
424 if (i<0||i>=fNSensors) AliFatal(Form("Sensor with id=%d is not in layer %d",i+fITSGeom->GetFirstChipIndex(fActiveID),fActiveID));
5009207f 425 return GetSensor(SensVIDtoMatrixID(i));
8b16dbae 426}