Code cleanup and bug fix Fixed valgrind warning concerning conditional jump on undefi...
[u/mrichter/AliRoot.git] / ITS / UPGRADE / ITSUpgradeRec / AliITSUCATrackingStation.cxx
1 //-*- Mode: C++ -*-
2 // **************************************************************************
3 // This file is property of and copyright by the ALICE ITSU Project         *
4 // ALICE Experiment at CERN, All rights reserved.                           *
5 //                                                                          *
6 // Primary Author: Ruben Shahoyan                                           *
7 //                                                                          *
8 // Adapted to ITSU: Maximiliano Puccio <maximiliano.puccio@cern.ch>         *
9 //                  for the ITS Upgrade project                             *
10 //                                                                          *
11 // Permission to use, copy, modify and distribute this software and its     *
12 // documentation strictly for non-commercial purposes is hereby granted     *
13 // without fee, provided that the above copyright notice appears in all     *
14 // copies and that both the copyright notice and this permission notice     *
15 // appear in the supporting documentation. The authors make no claims       *
16 // about the suitability of this software for any purpose. It is            *
17 // provided "as is" without express or implied warranty.                    *
18 //                                                                          *
19 //***************************************************************************
20
21 #include "AliITSUCATrackingStation.h"
22 #include <TMath.h>
23 #include "AliITSUGeomTGeo.h"
24 #include "AliVertex.h"
25 #include <TRandom.h>
26 #include <TStopwatch.h>
27 #include <TString.h>
28 #include "AliITSUAux.h"
29 #include "AliITSURecoSens.h"
30 #include <Riostream.h>
31
32 using AliITSUAux::BringTo02Pi;
33
34 //_________________________________________________________________
35 AliITSUCATrackingStation::AliITSUCATrackingStation()
36 :fClusters(0x0)
37 ,fVIDOffset(0)
38 ,fNClusters(0)
39 ,fZMin(0)
40 ,fZMax(0)
41 ,fDZInv(-1)
42 ,fDPhiInv(-1)
43 ,fNZBins(20)
44 ,fNPhiBins(20)
45 ,fQueryZBmin(-1)
46 ,fQueryZBmax(-1)
47 ,fQueryPhiBmin(-1)
48 ,fQueryPhiBmax(-1)
49 ,fBins(0)
50 ,fOccBins(0)
51 ,fNOccBins(0)
52 ,fNFoundClusters(0)
53 ,fFoundClusterIterator(0)
54 ,fFoundBinIterator(0)
55 ,fIndex()
56 ,fFoundBins(0)
57 ,fSortedClInfo(0)
58 ,fDetectors()
59 {
60   // def. c-tor
61 }
62
63 //_________________________________________________________________
64 AliITSUCATrackingStation::AliITSUCATrackingStation(int nzbins,int nphibins,
65                                                    AliITSURecoLayer *lr, AliITSUGeomTGeo *geo)
66 :fClusters(0x0)
67 ,fVIDOffset()
68 ,fNClusters(0)
69 ,fZMin(0.f)
70 ,fZMax(0.f)
71 ,fDZInv(-1)
72 ,fDPhiInv(-1)
73 ,fNZBins(nzbins)
74 ,fNPhiBins(nphibins)
75 ,fQueryZBmin(-1)
76 ,fQueryZBmax(-1)
77 ,fQueryPhiBmin(-1)
78 ,fQueryPhiBmax(-1)
79 ,fBins(0)
80 ,fOccBins(0)
81 ,fNOccBins(0)
82 ,fNFoundClusters(0)
83 ,fFoundClusterIterator(0)
84 ,fFoundBinIterator(0)
85 ,fIndex()
86 ,fFoundBins(0)
87 ,fSortedClInfo(0)
88 ,fDetectors()
89 {
90   // c-tor
91   Init(lr,geo);
92 }
93
94 //_________________________________________________________________
95 AliITSUCATrackingStation::~AliITSUCATrackingStation()
96 {
97   // d-tor
98   delete[] fBins;
99   delete[] fOccBins;
100 //  delete fClusters;
101 }
102
103 //_________________________________________________________________
104 void AliITSUCATrackingStation::Init(AliITSURecoLayer *lr, AliITSUGeomTGeo *geo)
105 {
106   if (fClusters) {
107     printf("Already initialized\n");
108     return;
109   }
110   
111   fClusters = *lr->GetClustersAddress();
112   fZMin = lr->GetZMin();
113   fZMax = lr->GetZMax();
114   if (fNZBins < 1)   fNZBins = 2;
115   if (fNPhiBins < 1) fNPhiBins = 1;
116   fDZInv   = fNZBins / (fZMax - fZMin);
117   fDPhiInv = fNPhiBins / TMath::TwoPi();
118   //
119   fBins = new ClBinInfo_t[fNZBins * fNPhiBins];
120   fOccBins = new int[fNZBins * fNPhiBins];
121   fNClusters = fClusters->GetEntriesFast();
122   if(fNClusters == 0) return;
123   fSortedClInfo.reserve(fClusters->GetEntriesFast());
124   fVIDOffset = ((AliITSUClusterPix*)fClusters->UncheckedAt(0))->GetVolumeId();
125   //
126   // prepare detectors info
127   int detID = -1;
128   fIndex.resize(lr->GetNSensors(),-1);
129   fDetectors.reserve(lr->GetNSensors());
130   for (int iCl = 0; iCl < fClusters->GetEntriesFast(); ++iCl)
131   { //Fill this layer with detectors
132     AliITSUClusterPix* c = (AliITSUClusterPix*)fClusters->UncheckedAt(iCl);
133     if (detID == c->GetVolumeId())
134     {
135       continue;
136     }
137     detID = c->GetVolumeId();
138     ITSDetInfo_t det;
139     det.index = iCl;
140     //
141     TGeoHMatrix m;
142     geo->GetOrigMatrix(detID,m);
143     //
144     fIndex[detID - fVIDOffset] = fDetectors.size();
145     AliITSURecoSens* sens = lr->GetSensorFromID(detID);
146     const TGeoHMatrix *tm = geo->GetMatrixT2L(detID);
147     m.Multiply(tm);
148     double txyz[3] = {0.,0.,0.}, xyz[3] = {0.,0.,0.};
149     m.LocalToMaster(txyz,xyz);
150     det.xTF = sens->GetXTF(); // TMath::Sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]);
151
152     det.phiTF = sens->GetPhiTF();//TMath::ATan2(xyz[1],xyz[0]);
153     det.sinTF = TMath::Sin(det.phiTF);
154     det.cosTF = TMath::Cos(det.phiTF);
155     //
156     // compute the real radius (with misalignment)
157     TGeoHMatrix mmisal(*(geo->GetMatrix(detID)));
158     mmisal.Multiply(tm);
159     xyz[0] = 0.;
160     xyz[1] = 0.;
161     xyz[2] = 0.;
162     mmisal.LocalToMaster(txyz,xyz);
163     det.xTFmisal = TMath::Sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]);
164 //    if (!TMath::AreEqualAbs(det.xTFmisal, sens->GetXTF(), 1e-9)) {
165 //      cout << "Error: " << lr->GetActiveID() << " " << detID << " " << det.xTFmisal - sens->GetXTF()<< endl;
166 //    }
167     
168     fDetectors.push_back(det);
169   } // end loop on detectors
170 }
171
172 //_________________________________________________________________
173 void AliITSUCATrackingStation::SortClusters(const AliVertex* vtx)
174 {
175   // sort clusters and build fast lookup table
176   //
177   ClearSortedInfo();
178   fSortedClInfo.reserve(fNClusters);
179   //
180   ClsInfo_t cl;
181   for (int icl = fNClusters;icl--;) {
182     AliITSUClusterPix* cluster = (AliITSUClusterPix*)fClusters->UncheckedAt(icl);
183     cluster->GetGlobalXYZ( (float*)&cl );
184     //
185     if (vtx) { // phi and r will be computed wrt vertex
186       cl.x -= vtx->GetX();
187       cl.y -= vtx->GetY();
188     }
189     //
190     cl.r = TMath::Sqrt(cl.x*cl.x + cl.y*cl.y);
191     cl.phi = TMath::ATan2(cl.y,cl.x);
192     BringTo02Pi(cl.phi);
193     cl.index = icl;
194     cl.zphibin = GetBinIndex(GetZBin(cl.z),GetPhiBin(cl.phi));
195     cl.detid = cluster->GetVolumeId() - fVIDOffset;
196     //
197     fSortedClInfo.push_back( cl );
198     //
199   }
200   sort(fSortedClInfo.begin(), fSortedClInfo.end()); // sort in phi, z
201   //
202   // fill cells in phi,z
203   int currBin = -1;
204   for (int icl = 0;icl < fNClusters; ++icl)
205   {
206     ClsInfo_t &t = fSortedClInfo[icl];
207     if (t.zphibin > currBin)
208     { // register new occupied bin
209       currBin = t.zphibin;
210       fBins[currBin].first = icl;
211       fBins[currBin].index = fNOccBins;
212       fOccBins[fNOccBins++] = currBin;
213     }
214     fBins[currBin].ncl++;
215   }
216   //  Print("clb"); //RS
217 }
218
219 //_________________________________________________________________
220 void AliITSUCATrackingStation::Clear(Option_t *)
221 {
222   // clear cluster info
223   ClearSortedInfo();
224   fIndex.clear();
225   fNClusters = 0;
226   if (fClusters) fClusters = 0x0;
227   //
228 }
229
230 //_________________________________________________________________
231 void AliITSUCATrackingStation::ClearSortedInfo()
232 {
233   // clear cluster info
234   fSortedClInfo.clear();
235   memset(fBins,0,fNZBins * fNPhiBins * sizeof(ClBinInfo_t));
236   memset(fOccBins,0,fNZBins * fNPhiBins * sizeof(int));
237   fNOccBins = 0;
238 }
239
240 //_________________________________________________________________
241 void AliITSUCATrackingStation::Print(Option_t *opt) const
242 {
243   // dump cluster bins info
244   TString opts = opt;
245   opts.ToLower();
246   printf("Stored %d clusters in %d occupied bins\n",fNClusters,fNOccBins);
247   //
248   if (opts.Contains("c")) {
249     printf("\nCluster info\n");
250     for (int i = 0; i < fNClusters;i++)
251     {
252       const ClsInfo_t &t = fSortedClInfo[i];
253       printf("#%5d Bin(phi/z):%03d/%03d Z:%+8.3f Phi:%+6.3f R:%7.3f Ind:%d ",
254              i,t.zphibin/fNZBins,t.zphibin%fNZBins,t.z,t.phi,t.r,t.index);
255       if (opts.Contains("l")) { // mc labels
256         AliITSUClusterPix* rp = (AliITSUClusterPix*)fClusters->UncheckedAt(t.index);
257         for (int l = 0;l < 3; l++) if (rp->GetLabel(l) >= 0) printf("| %d ",rp->GetLabel(l));
258       }
259       printf("\n");
260     }
261   }
262   //
263   if (opts.Contains("b")) {
264     printf("\nBins info (occupied only)\n");
265     for (int i=0;i<fNOccBins;i++) {
266       printf("%4d %5d(phi/z: %03d/%03d) -> %3d cl from %d\n",i,fOccBins[i],fOccBins[i]/fNZBins,fOccBins[i]%fNZBins,
267              fBins[fOccBins[i]].ncl,fBins[fOccBins[i]].first);
268     }
269   }
270   //
271 }
272
273 //_____________________________________________________________
274 int AliITSUCATrackingStation::SelectClusters(float zmin,float zmax,float phimin,float phimax)
275 {
276   // prepare occupied bins in the requested region
277   //printf("Select: Z %f %f | Phi: %f %f\n",zmin,zmax,phimin,phimax);
278   if (!fNOccBins) return 0;
279   if (zmax < fZMin || zmin > fZMax || zmin > zmax) return 0;
280   fFoundBins.clear();
281   
282   fQueryZBmin = GetZBin(zmin);
283   if (fQueryZBmin < 0) fQueryZBmin = 0;
284   fQueryZBmax = GetZBin(zmax);
285   if (fQueryZBmax >= fNZBins) fQueryZBmax = fNZBins - 1;
286   BringTo02Pi(phimin);
287   BringTo02Pi(phimax);
288   fQueryPhiBmin = GetPhiBin(phimin);
289   fQueryPhiBmax = GetPhiBin(phimax);
290   int dbz = 0;
291   fNFoundClusters = 0;
292   int nbcheck = fQueryPhiBmax - fQueryPhiBmin + 1; //TODO:(MP) check if a circular buffer is feasible
293   if (nbcheck > 0)
294   { // no wrapping around 0-2pi, fast case
295     for (int ip = fQueryPhiBmin;ip <= fQueryPhiBmax;ip++) {
296       int binID = GetBinIndex(fQueryZBmin,ip);
297       if ( !(dbz = (fQueryZBmax-fQueryZBmin)) ) { // just one Z bin in the query range
298         ClBinInfo_t& binInfo = fBins[binID];
299         if (!binInfo.ncl) continue;
300         fNFoundClusters += binInfo.ncl;
301         fFoundBins.push_back(binID);
302         continue;
303       }
304       int binMax = binID + dbz;
305       for ( ; binID <= binMax; binID++) {
306         ClBinInfo_t& binInfo = fBins[binID];
307         if (!binInfo.ncl) continue;
308         fNFoundClusters += binInfo.ncl;
309         fFoundBins.push_back(binID);
310       }
311     }
312   }
313   else
314   {  // wrapping
315     nbcheck += fNPhiBins;
316     for (int ip0 = 0;ip0 <= nbcheck;ip0++) {
317       int ip = fQueryPhiBmin + ip0;
318       if (ip >= fNPhiBins) ip -= fNPhiBins;
319       int binID = GetBinIndex(fQueryZBmin,ip);
320       if ( !(dbz = (fQueryZBmax - fQueryZBmin)) ) { // just one Z bin in the query range
321         ClBinInfo_t& binInfo = fBins[binID];
322         if (!binInfo.ncl) continue;
323         fNFoundClusters += binInfo.ncl;
324         fFoundBins.push_back(binID);
325         continue;
326       }
327       int binMax = binID + dbz;
328       for (;binID <= binMax;binID++) {
329         ClBinInfo_t& binInfo = fBins[binID];
330         if (!binInfo.ncl) continue;
331         fNFoundClusters += binInfo.ncl;
332         fFoundBins.push_back(binID);
333       }
334     }
335   }
336   fFoundClusterIterator = fFoundBinIterator = 0;
337   /*
338    //printf("Selected -> %d cl in %d bins\n",fNFoundClusters,(int)fFoundBins.size());
339    for (int i=0;i<(int)fFoundBins.size();i++) {
340    int bn = fFoundBins[i];
341    ClBinInfo_t& bin=fBins[bn];
342    printf("#%d b:%d 1st: %3d Ncl:%d\n",i,bn,bin.first,bin.ncl);
343    }
344    printf("\n");
345    */
346   return fNFoundClusters;
347 }
348
349 //_____________________________________________________________
350 int AliITSUCATrackingStation::GetNextClusterInfoID()
351 {
352   if (fFoundBinIterator < 0) return 0;
353   int currBin = fFoundBins[fFoundBinIterator];
354   if (fFoundClusterIterator < fBins[currBin].ncl) { // same bin
355     return fBins[currBin].first + fFoundClusterIterator++;
356   }
357   if (++fFoundBinIterator < int(fFoundBins.size())) {  // need to change bin
358     currBin = fFoundBins[fFoundBinIterator];
359     fFoundClusterIterator = 1;
360     return fBins[currBin].first;
361   }
362   fFoundBinIterator = -1;
363   return -1;
364 }
365
366 //_____________________________________________________________
367 void AliITSUCATrackingStation::ResetFoundIterator()
368 {
369   // prepare for a new loop over found clusters
370   if (fNFoundClusters)  fFoundClusterIterator = fFoundBinIterator = 0;
371 }