Line data Source code
1 : #include <TMath.h>
2 : #include "AliITSSAPLayer.h"
3 : #include "AliITSSAPAux.h"
4 : #include "AliITSRecPoint.h"
5 : #include "AliITSgeomTGeo.h"
6 : #include "AliVertex.h"
7 : #include <TRandom.h>
8 : #include <TStopwatch.h>
9 : #include <TString.h>
10 :
11 :
12 : using namespace AliITSSAPAux;
13 :
14 : //_________________________________________________________________
15 : AliITSSAPLayer::AliITSSAPLayer() :
16 0 : fClusters(0)
17 0 : ,fLrID(-1)
18 0 : ,fVIDOffset(0)
19 0 : ,fNClusters(0)
20 0 : ,fZMin(0)
21 0 : ,fZMax(0)
22 0 : ,fDZInv(-1)
23 0 : ,fDPhiInv(-1)
24 0 : ,fNZBins(20)
25 0 : ,fNPhiBins(20)
26 0 : ,fQueryZBmin(-1)
27 0 : ,fQueryZBmax(-1)
28 0 : ,fQueryPhiBmin(-1)
29 0 : ,fQueryPhiBmax(-1)
30 0 : ,fBins(0)
31 0 : ,fOccBins(0)
32 0 : ,fNOccBins(0)
33 0 : ,fNFoundClusters(0)
34 0 : ,fFoundClusterIterator(0)
35 0 : ,fFoundBinIterator(0)
36 0 : ,fFoundBins()
37 0 : ,fSortedClInfo()
38 0 : ,fDetectors()
39 0 : {
40 : // def. c-tor
41 0 : }
42 :
43 : //_________________________________________________________________
44 : AliITSSAPLayer::AliITSSAPLayer(int id, float zspan,int nzbins,int nphibins, int buffer) :
45 0 : fClusters(0)
46 0 : ,fLrID(id)
47 0 : ,fVIDOffset((id+1)*2048)
48 0 : ,fNClusters(0)
49 0 : ,fZMin(-zspan)
50 0 : ,fZMax(zspan)
51 0 : ,fDZInv(-1)
52 0 : ,fDPhiInv(-1)
53 0 : ,fNZBins(nzbins)
54 0 : ,fNPhiBins(nphibins)
55 0 : ,fQueryZBmin(-1)
56 0 : ,fQueryZBmax(-1)
57 0 : ,fQueryPhiBmin(-1)
58 0 : ,fQueryPhiBmax(-1)
59 0 : ,fBins(0)
60 0 : ,fOccBins(0)
61 0 : ,fNOccBins(0)
62 0 : ,fNFoundClusters(0)
63 0 : ,fFoundClusterIterator(0)
64 0 : ,fFoundBinIterator(0)
65 0 : ,fFoundBins()
66 0 : ,fSortedClInfo()
67 0 : ,fDetectors()
68 0 : {
69 : // c-tor
70 0 : Init(buffer);
71 0 : }
72 :
73 : //_________________________________________________________________
74 : AliITSSAPLayer::~AliITSSAPLayer()
75 0 : {
76 : // d-tor
77 0 : delete[] fBins;
78 0 : delete[] fOccBins;
79 0 : delete fClusters;
80 0 : }
81 :
82 : //_________________________________________________________________
83 : void AliITSSAPLayer::Init(int buffer)
84 : {
85 0 : if (fClusters) {
86 0 : printf("Already initialized\n");
87 0 : return;
88 : }
89 0 : if (fNZBins<1) fNZBins = 2;
90 0 : if (fNPhiBins<1) fNPhiBins = 1;
91 0 : fDZInv = fNZBins/(fZMax-fZMin);
92 0 : fDPhiInv = fNPhiBins/TMath::TwoPi();
93 : //
94 0 : fBins = new ClBinInfo_t[fNZBins*fNPhiBins];
95 0 : fOccBins = new int[fNZBins*fNPhiBins];
96 0 : if (buffer<100) buffer = 100;
97 0 : fClusters = new TObjArray(buffer);
98 0 : fSortedClInfo.reserve(buffer);
99 : //
100 : // prepare detectors info
101 0 : int id1 = fLrID+1;
102 0 : Int_t nlad=AliITSgeomTGeo::GetNLadders(id1);
103 0 : Int_t ndet=AliITSgeomTGeo::GetNDetectors(id1);
104 : int detID = 0;
105 0 : for (Int_t j=1; j<nlad+1; j++) {
106 0 : for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
107 0 : ITSDetInfo_t det;
108 0 : det.index = detID++;
109 : //
110 0 : TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(id1,j,k,m);
111 0 : const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(id1,j,k);
112 0 : m.Multiply(tm);
113 0 : Double_t txyz[3] = {0.}, xyz[3] = {0.};
114 0 : m.LocalToMaster(txyz,xyz);
115 0 : det.xTF = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
116 0 : det.phiTF = TMath::ATan2(xyz[1],xyz[0]);
117 : //BringTo02Pi(det.phiTF);
118 0 : det.sinTF = TMath::Sin(det.phiTF);
119 0 : det.cosTF = TMath::Cos(det.phiTF);
120 : //
121 : // compute the real radius (with misalignment)
122 0 : TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(id1,j,k)));
123 0 : mmisal.Multiply(tm);
124 0 : xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
125 0 : mmisal.LocalToMaster(txyz,xyz);
126 0 : det.xTFmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
127 : //
128 0 : fDetectors.push_back(det);
129 0 : } // end loop on detectors
130 : } // end loop on ladders
131 :
132 0 : }
133 :
134 : //_________________________________________________________________
135 : void AliITSSAPLayer::SortClusters(const AliVertex* vtx)
136 : {
137 : // sort clusters and build fast lookup table
138 : //
139 0 : ClearSortedInfo();
140 0 : fSortedClInfo.reserve(fNClusters);
141 : //
142 0 : ClsInfo_t cl;
143 0 : for (int icl=fNClusters;icl--;) {
144 0 : AliITSRecPoint* cluster = (AliITSRecPoint*)fClusters->UncheckedAt(icl);
145 0 : cluster->GetGlobalXYZ( (float*)&cl );
146 : //
147 0 : if (vtx) { // phi and r will be computed wrt vertex
148 0 : cl.x -= vtx->GetX();
149 0 : cl.y -= vtx->GetY();
150 0 : }
151 : //
152 0 : cl.r = TMath::Sqrt(cl.x*cl.x + cl.y*cl.y);
153 0 : cl.phi = TMath::ATan2(cl.y,cl.x);
154 0 : BringTo02Pi(cl.phi);
155 0 : cl.index = icl;
156 0 : cl.zphibin = GetBinIndex(GetZBin(cl.z),GetPhiBin(cl.phi));
157 0 : cl.detid = cluster->GetVolumeId() - fVIDOffset;
158 : //
159 0 : fSortedClInfo.push_back( cl );
160 : //
161 : }
162 0 : sort(fSortedClInfo.begin(), fSortedClInfo.end()); // sort in phi, z
163 : //
164 : // fill cells in phi,z
165 : int currBin = -1;
166 0 : for (int icl=0;icl<fNClusters;icl++) {
167 0 : ClsInfo_t &t = fSortedClInfo[icl];
168 0 : if (t.zphibin>currBin) { // register new occupied bin
169 : currBin = t.zphibin;
170 0 : fBins[currBin].first = icl;
171 0 : fBins[currBin].index = fNOccBins;
172 0 : fOccBins[fNOccBins++] = currBin;
173 0 : }
174 0 : fBins[currBin].ncl++;
175 : }
176 : // Print("clb"); //RS
177 0 : }
178 :
179 : //_________________________________________________________________
180 : void AliITSSAPLayer::Clear(Option_t *)
181 : {
182 : // clear cluster info
183 0 : ClearSortedInfo();
184 0 : fNClusters = 0;
185 0 : if (fClusters) fClusters->Clear();
186 : //
187 0 : }
188 :
189 : //_________________________________________________________________
190 : void AliITSSAPLayer::ClearSortedInfo()
191 : {
192 : // clear cluster info
193 0 : fSortedClInfo.clear();
194 0 : memset(fBins,0,fNZBins*fNPhiBins*sizeof(ClBinInfo_t));
195 0 : memset(fOccBins,0,fNZBins*fNPhiBins*sizeof(int));
196 0 : fNOccBins = 0;
197 0 : }
198 :
199 : //_________________________________________________________________
200 : void AliITSSAPLayer::Print(Option_t *opt) const
201 : {
202 : // dump cluster bins info
203 0 : TString opts = opt;
204 0 : opts.ToLower();
205 0 : printf("Stored %d clusters in %d occupied bins\n",fNClusters,fNOccBins);
206 : //
207 0 : if (opts.Contains("c")) {
208 0 : printf("\nCluster info\n");
209 0 : for (int i=0;i<fNClusters;i++) {
210 0 : const ClsInfo_t &t = fSortedClInfo[i];
211 0 : printf("#%5d Bin(phi/z):%03d/%03d Z:%+8.3f Phi:%+6.3f R:%7.3f Ind:%d ",
212 0 : i,t.zphibin/fNZBins,t.zphibin%fNZBins,t.z,t.phi,t.r,t.index);
213 0 : if (opts.Contains("l")) { // mc labels
214 0 : AliITSRecPoint* rp = (AliITSRecPoint*)fClusters->UncheckedAt(t.index);
215 0 : for (int l=0;l<3;l++) if (rp->GetLabel(l)>=0) printf("| %d ",rp->GetLabel(l));
216 0 : }
217 0 : printf("\n");
218 : }
219 0 : }
220 : //
221 0 : if (opts.Contains("b")) {
222 0 : printf("\nBins info (occupied only)\n");
223 0 : for (int i=0;i<fNOccBins;i++) {
224 0 : printf("%4d %5d(phi/z: %03d/%03d) -> %3d cl from %d\n",i,fOccBins[i],fOccBins[i]/fNZBins,fOccBins[i]%fNZBins,
225 0 : fBins[fOccBins[i]].ncl,fBins[fOccBins[i]].first);
226 : }
227 0 : }
228 : //
229 0 : }
230 :
231 : //_____________________________________________________________
232 : int AliITSSAPLayer::SelectClusters(float zmin,float zmax,float phimin,float phimax)
233 : {
234 : // prepare occupied bins in the requested region
235 : //printf("Select: Z %f %f | Phi: %f %f\n",zmin,zmax,phimin,phimax);
236 0 : if (!fNOccBins) return 0;
237 0 : if (zmax<fZMin || zmin>fZMax || zmin>zmax) return 0;
238 0 : fFoundBins.clear();
239 :
240 0 : fQueryZBmin = GetZBin(zmin);
241 0 : if (fQueryZBmin<0) fQueryZBmin = 0;
242 0 : fQueryZBmax = GetZBin(zmax);
243 0 : if (fQueryZBmax>=fNZBins) fQueryZBmax = fNZBins-1;
244 0 : BringTo02Pi(phimin);
245 0 : BringTo02Pi(phimax);
246 0 : fQueryPhiBmin = GetPhiBin(phimin);
247 0 : fQueryPhiBmax = GetPhiBin(phimax);
248 : int dbz=0;
249 0 : fNFoundClusters = 0;
250 0 : int nbcheck = fQueryPhiBmax - fQueryPhiBmin + 1;
251 0 : if (nbcheck>0) { // no wrapping around 0-2pi, fast case
252 0 : for (int ip=fQueryPhiBmin;ip<=fQueryPhiBmax;ip++) {
253 0 : int binID = GetBinIndex(fQueryZBmin,ip);
254 0 : if ( !(dbz=(fQueryZBmax-fQueryZBmin)) ) { // just one Z bin in the query range
255 0 : ClBinInfo_t& binInfo = fBins[binID];
256 0 : if (!binInfo.ncl) continue;
257 0 : fNFoundClusters += binInfo.ncl;
258 0 : fFoundBins.push_back(binID);
259 0 : continue;
260 : }
261 0 : int binMax = binID+dbz;
262 0 : for (;binID<=binMax;binID++) {
263 0 : ClBinInfo_t& binInfo = fBins[binID];
264 0 : if (!binInfo.ncl) continue;
265 0 : fNFoundClusters += binInfo.ncl;
266 0 : fFoundBins.push_back(binID);
267 0 : }
268 0 : }
269 0 : }
270 : else { // wrapping
271 0 : nbcheck += fNPhiBins;
272 0 : for (int ip0=0;ip0<=nbcheck;ip0++) {
273 0 : int ip = fQueryPhiBmin+ip0;
274 0 : if (ip>=fNPhiBins) ip -= fNPhiBins;
275 0 : int binID = GetBinIndex(fQueryZBmin,ip);
276 0 : if ( !(dbz=(fQueryZBmax-fQueryZBmin)) ) { // just one Z bin in the query range
277 0 : ClBinInfo_t& binInfo = fBins[binID];
278 0 : if (!binInfo.ncl) continue;
279 0 : fNFoundClusters += binInfo.ncl;
280 0 : fFoundBins.push_back(binID);
281 0 : continue;
282 : }
283 0 : int binMax = binID+dbz;
284 0 : for (;binID<=binMax;binID++) {
285 0 : ClBinInfo_t& binInfo = fBins[binID];
286 0 : if (!binInfo.ncl) continue;
287 0 : fNFoundClusters += binInfo.ncl;
288 0 : fFoundBins.push_back(binID);
289 0 : }
290 0 : }
291 : }
292 0 : fFoundClusterIterator = fFoundBinIterator = 0;
293 : /*
294 : //printf("Selected -> %d cl in %d bins\n",fNFoundClusters,(int)fFoundBins.size());
295 : for (int i=0;i<(int)fFoundBins.size();i++) {
296 : int bn = fFoundBins[i];
297 : ClBinInfo_t& bin=fBins[bn];
298 : printf("#%d b:%d 1st: %3d Ncl:%d\n",i,bn,bin.first,bin.ncl);
299 : }
300 : printf("\n");
301 : */
302 0 : return fNFoundClusters;
303 0 : }
304 :
305 : //_____________________________________________________________
306 : int AliITSSAPLayer::GetNextClusterInfoID()
307 : {
308 0 : if (fFoundBinIterator<0) return 0;
309 0 : int currBin = fFoundBins[fFoundBinIterator];
310 0 : if (fFoundClusterIterator<fBins[currBin].ncl) { // same bin
311 0 : return fBins[currBin].first+fFoundClusterIterator++;
312 : }
313 0 : if (++fFoundBinIterator<int(fFoundBins.size())) { // need to change bin
314 0 : currBin = fFoundBins[fFoundBinIterator];
315 0 : fFoundClusterIterator = 1;
316 0 : return fBins[currBin].first;
317 : }
318 0 : fFoundBinIterator = -1;
319 0 : return -1;
320 0 : }
321 :
322 : //_____________________________________________________________
323 : void AliITSSAPLayer::ResetFoundIterator()
324 : {
325 : // prepare for a new loop over found clusters
326 0 : if (fNFoundClusters) fFoundClusterIterator = fFoundBinIterator = 0;
327 0 : }
|
