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