Line data Source code
1 : // $Id$
2 :
3 : //**************************************************************************
4 : //* This file is property of and copyright by the ALICE HLT Project *
5 : //* ALICE Experiment at CERN, All rights reserved. *
6 : //* *
7 : //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no> *
8 : //* Sergey Gorbunov *
9 : //* for The ALICE HLT 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 : // @file AliHLTTPCHWClusterMerger.cxx
21 : // @author Matthias Richter, Sergey Gorbunov
22 : // @date 2011-11-25
23 : // @brief Merger class for HLT TPC Hardware clusters
24 : // Handles merging of branch border clusters
25 :
26 : #include <algorithm>
27 :
28 : #include "AliHLTTPCHWClusterMerger.h"
29 : #include "AliHLTTPCGeometry.h"
30 : #include "AliHLTTPCSpacePointData.h"
31 : #include "AliHLTTPCHWCFSupport.h"
32 :
33 : /** ROOT macro for the implementation of ROOT specific class methods */
34 6 : ClassImp(AliHLTTPCHWClusterMerger)
35 :
36 : AliHLTTPCHWClusterMerger::AliHLTTPCHWClusterMerger()
37 0 : : AliHLTLogging()
38 0 : , fMapping(0)
39 0 : , fNRows(0)
40 0 : , fNRowPads(0)
41 0 : , fNBorders(0)
42 0 : , fBorders(0)
43 0 : , fBorderNClusters(0)
44 0 : , fBorderFirstCluster(0)
45 0 : , fBorderClusters(0)
46 0 : , fBorderNClustersTotal(0)
47 0 : , fClusters()
48 0 : , fRemovedClusterIds()
49 0 : , fIter()
50 0 : , fEnd()
51 0 : {
52 : // see header file for class documentation
53 : // or
54 : // refer to README to build package
55 : // or
56 : // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
57 : // constructor
58 0 : }
59 :
60 : AliHLTTPCHWClusterMerger::~AliHLTTPCHWClusterMerger()
61 0 : {
62 : // destructor
63 0 : delete[] fMapping;
64 0 : delete [] fBorders;
65 0 : delete[] fBorderNClusters;
66 0 : delete[] fBorderFirstCluster;
67 0 : delete[] fBorderClusters;
68 0 : }
69 :
70 : void AliHLTTPCHWClusterMerger::Init()
71 : {
72 : // initialisation
73 :
74 0 : delete[] fMapping;
75 0 : delete[] fBorders;
76 0 : delete[] fBorderNClusters;
77 0 : delete[] fBorderFirstCluster;
78 0 : delete[] fBorderClusters;
79 0 : fMapping = 0;
80 0 : fBorders = 0;
81 0 : fBorderNClusters = 0;
82 0 : fBorderFirstCluster = 0;
83 0 : fBorderClusters = 0;
84 0 : fBorderNClustersTotal = 0;
85 0 : fNRows = AliHLTTPCGeometry::GetNRows();
86 0 : fNRowPads = 0;
87 0 : fNBorders = 0;
88 0 : for( int i=0; i<fNRows; i++ ){
89 0 : int nPads = AliHLTTPCGeometry::GetNPads(i);
90 0 : if( fNRowPads<nPads ) fNRowPads = nPads;
91 : }
92 0 : int nPadsTotal = fNRows*fNRowPads;
93 0 : fMapping = new AliHLTInt16_t [nPadsTotal];
94 :
95 0 : if( !fMapping ){
96 0 : HLTError("Can not allocate memory: %d bytes", nPadsTotal*sizeof(AliHLTInt16_t));
97 0 : return;
98 : }
99 0 : for( int i=0; i<nPadsTotal; i++ ){
100 0 : fMapping[i] = -1;
101 : }
102 :
103 0 : AliHLTTPCHWCFSupport support;
104 :
105 0 : for( int iPart=0; iPart<6; iPart++ ){
106 0 : const AliHLTUInt32_t *m = support.GetMapping(0,iPart);
107 0 : int nHWAddress=m[0];
108 0 : for( int iHW=0; iHW<nHWAddress; iHW++ ){
109 0 : AliHLTUInt32_t configWord = m[iHW+1];
110 0 : int row = (configWord>>8) & 0x3F;
111 0 : int pad = configWord & 0xFF;
112 0 : bool border = (configWord>>14) & 0x1;
113 0 : if( !border ) continue;
114 0 : row+=AliHLTTPCGeometry::GetFirstRow(iPart);
115 0 : if( row>=fNRows ) continue;
116 0 : if( pad>=AliHLTTPCGeometry::GetNPads(row) ) continue;
117 0 : fMapping[row*fNRowPads + pad] = -2;
118 0 : }
119 : }
120 :
121 0 : std::vector<AliHLTFloat32_t> vBorders;
122 0 : for( int row=0; row<fNRows; row++ ){
123 0 : for( int pad=0; pad<fNRowPads-1; pad++ ){
124 0 : AliHLTInt16_t *m = fMapping + row*fNRowPads;
125 0 : if( m[pad]==-2 && m[pad+1]==-2 ){
126 0 : m[pad] = fNBorders;
127 0 : for( int i=1; i<fkMergeWidth; i++ ){
128 0 : if( pad-i<0 || m[pad-i]!=-1 ) break;
129 0 : m[pad-i] = fNBorders;
130 : }
131 0 : m[pad+1] = fNBorders+1;
132 0 : for( int i=1; i<fkMergeWidth; i++ ){
133 0 : if( pad+1+i>=fNRowPads || m[pad+1+i]!=-1 ) break;
134 0 : m[pad+1+i] = fNBorders+1;
135 : }
136 0 : fNBorders+=2;
137 0 : vBorders.push_back(pad+1.);
138 0 : vBorders.push_back(pad+1.);
139 : }
140 : }
141 : }
142 :
143 : //cout<<" NBorders = "<<fNBorders/2<<endl;
144 :
145 : /*
146 : cout<<"Borders:"<<endl;
147 : for( int row=0; row<fNRows; row++ ){
148 : for( int pad=0; pad<fNRowPads; pad++ ){
149 : AliHLTInt16_t *m = fMapping + row*fNRowPads;
150 : if( m[pad]>=0 ) {
151 : cout<<row<<" "<<pad<<" "<<m[pad]%2<<endl;
152 : }
153 : }
154 : }
155 : */
156 :
157 0 : fBorders = new AliHLTFloat32_t [fNBorders];
158 0 : if( !fBorders ){
159 0 : HLTError("Can not allocate memory: %d bytes", fNBorders*sizeof(AliHLTFloat32_t));
160 0 : fNBorders = 0;
161 0 : return;
162 : }
163 0 : for( int i=0; i<fNBorders; i++ ){
164 0 : fBorders[i] = vBorders[i];
165 : }
166 0 : fBorderNClusters = new int[fkNSlices*fNBorders];
167 0 : if( !fBorderNClusters ){
168 0 : HLTError("Can not allocate memory: %d bytes", fkNSlices*fNBorders*sizeof(int));
169 0 : return;
170 : }
171 :
172 0 : fBorderFirstCluster = new int[fkNSlices*fNBorders];
173 0 : if( !fBorderFirstCluster ){
174 0 : HLTError("Can not allocate memory: %d bytes", fkNSlices*fNBorders*sizeof(int));
175 0 : return;
176 : }
177 :
178 0 : for( int i=0; i<fkNSlices*fNBorders; i++){
179 0 : fBorderNClusters[i] = 0;
180 0 : fBorderFirstCluster[i] = 0;
181 : }
182 0 : }
183 :
184 : bool AliHLTTPCHWClusterMerger::CheckCandidate(int /*slice*/, int partition, int partitionrow, float pad, float /*time*/, float sigmaPad2
185 : )
186 : {
187 : /// check cluster if it is a candidate for merging
188 0 : int slicerow=partitionrow+AliHLTTPCGeometry::GetFirstRow(partition);
189 :
190 0 : if( !fMapping ) Init();
191 0 : if( !fMapping ) return 0;
192 0 : if( slicerow <0 || slicerow>= fNRows ) return 0;
193 0 : int iPad = (int) pad;
194 0 : if( iPad<0 || iPad>=fNRowPads ) return 0;
195 0 : int atBorder = fMapping[slicerow*fNRowPads+iPad];
196 0 : if( atBorder<0 ) return 0;
197 0 : float dPad = pad - fBorders[atBorder];
198 0 : if( sigmaPad2>1.e-4 ){
199 0 : if( dPad*dPad > 12.*sigmaPad2 ) return 0;
200 : } else {
201 0 : if( fabs(dPad)>1. ) return 0;
202 : }
203 0 : return 1;
204 0 : }
205 :
206 : int AliHLTTPCHWClusterMerger::AddCandidate(int slice,
207 : int partition,
208 : short partitionrow,
209 : float pad,
210 : float time,
211 : float sigmaPad2,
212 : float sigmaTime2,
213 : unsigned short charge,
214 : unsigned short qmax,
215 : unsigned short flags,
216 : AliHLTUInt32_t id,
217 : const AliHLTTPCClusterMCLabel *mc
218 : )
219 : {
220 : /// add a candidate for merging and register in the index grid
221 :
222 : int iBorder = -1;
223 0 : int iPad = (int) pad;
224 0 : if( !fMapping ) Init();
225 :
226 0 : int slicerow=partitionrow+AliHLTTPCGeometry::GetFirstRow(partition);
227 :
228 0 : if (id!=~AliHLTUInt32_t(0)) {
229 0 : if (slice<0) {
230 0 : slice=AliHLTTPCSpacePointData::GetSlice(id);
231 0 : } else if ((unsigned)slice!=AliHLTTPCSpacePointData::GetSlice(id)) {
232 0 : HLTError("cluster id 0x%08x is not consistent with specified slice %d", id, slice);
233 : }
234 0 : if (partition<0) {
235 0 : partition=AliHLTTPCSpacePointData::GetPatch(id);
236 0 : } else if ((unsigned)partition!=AliHLTTPCSpacePointData::GetPatch(id)) {
237 0 : HLTError("cluster id 0x%08x is not consistent with specified partition %d", id, partition);
238 : }
239 : }
240 :
241 0 : if( slice<0 || slice>=fkNSlices ){
242 0 : HLTError("cluster id 0x%08x has wrong slice number %d", id, slice);
243 : }
244 :
245 0 : if( partition<0 || partition>=6 ){
246 0 : HLTError("cluster id 0x%08x has wrong partition number %d", id, partition);
247 : }
248 :
249 0 : if( slice>=0 && slice<fkNSlices && slicerow>=0 && slicerow <fNRows && fMapping && iPad>=0 && iPad < fNRowPads ){
250 0 : iBorder = fMapping[slicerow*fNRowPads+iPad];
251 :
252 0 : if( iBorder>=0 ){
253 0 : float dPad = pad - fBorders[iBorder];
254 0 : if( sigmaPad2>1.e-4 ){
255 0 : if( dPad*dPad > 12.*sigmaPad2 ) iBorder = -1;
256 : } else {
257 0 : if( fabs(dPad)>1. ) iBorder = -1;
258 : }
259 0 : }
260 0 : if( iBorder>=0 ) iBorder = slice*fNBorders + iBorder;
261 : }
262 :
263 0 : fClusters.push_back(AliClusterRecord(slice, partition, iBorder, -1, id,
264 0 : AliHLTTPCRawCluster(partitionrow, pad, time, sigmaPad2, sigmaTime2, charge, qmax, flags),
265 0 : mc!=NULL?*mc:AliHLTTPCClusterMCLabel() ));
266 :
267 0 : if( iBorder>=0 ){
268 0 : fBorderNClusters[iBorder]++;
269 0 : fBorderNClustersTotal++;
270 0 : }
271 0 : return fClusters.size()-1;
272 0 : }
273 :
274 : int AliHLTTPCHWClusterMerger::FillIndex()
275 : {
276 : /// loop over cached raw clusters and fill index
277 :
278 0 : delete[] fBorderClusters;
279 0 : fBorderClusters = 0;
280 0 : fBorderClusters = new AliBorderRecord[fBorderNClustersTotal];
281 :
282 0 : fBorderFirstCluster[0] = 0;
283 0 : for(int i=1; i<fkNSlices*fNBorders; i++ ){
284 0 : fBorderFirstCluster[i] = fBorderFirstCluster[i-1] + fBorderNClusters[i-1];
285 0 : fBorderNClusters[i-1] = 0;
286 : }
287 0 : fBorderNClusters[fkNSlices*fNBorders-1]=0;
288 0 : for (AliHLTUInt32_t pos=0; pos<fClusters.size(); pos++) {
289 0 : int iBorder = fClusters[pos].GetBorder();
290 0 : if( iBorder<0 ) continue;
291 0 : AliBorderRecord &b = fBorderClusters[fBorderFirstCluster[iBorder]+fBorderNClusters[iBorder]];
292 0 : b.fClusterRecordID = pos;
293 0 : b.fTimeBin = (int)fClusters[pos].GetCluster().GetTime();
294 0 : fBorderNClusters[iBorder]++;
295 0 : }
296 : /*
297 : cout<<"Border clusters: "<<endl;
298 : for( int iSlice=0; iSlice<fkNSlices; iSlice++){
299 : for( int ib = 0; ib<fNBorders; ib++ ){
300 : int ind = iSlice*fNBorders+ib;
301 : cout<<iSlice<<" "<<ib<<" :"<<fBorderFirstCluster[ind]<<", "<<fBorderNClusters[ind]<<endl;
302 : }
303 : }
304 : */
305 : /*
306 : for( int ib=1; ib<fkNSlices*fNBorders; ib++ ){
307 : if( fBorderFirstCluster[ib] != fBorderFirstCluster[ib-1]+fBorderNClusters[ib-1] ){
308 : cout<<"Something wrong with cluster borders !!! "<<endl;
309 : }
310 : }
311 : */
312 0 : return 0;
313 : }
314 :
315 : void AliHLTTPCHWClusterMerger::Clear()
316 : {
317 : /// cleanup
318 0 : fClusters.clear();
319 0 : fRemovedClusterIds.clear();
320 :
321 0 : delete[] fBorderClusters;
322 0 : fBorderClusters = 0;
323 0 : for( int i=0; i<fkNSlices*fNBorders; i++){
324 0 : fBorderNClusters[i] = 0;
325 0 : fBorderFirstCluster[i] = 0;
326 : }
327 0 : fBorderNClustersTotal = 0;
328 0 : }
329 :
330 :
331 : int AliHLTTPCHWClusterMerger::Merge()
332 : {
333 : /// merge clusters
334 : /// first all candidates are filled into the index grid, looping over all clusters
335 : /// in the grid automatically results in time ordered clusters per row (ordering with
336 : /// respect to cells, within a cell two clusters can be reversed)
337 :
338 0 : if( !fMapping ) Init();
339 0 : if( !fMapping ) return 0;
340 : //static int sLeft = 0, sRight = 0, sMerged = 0;
341 : int iResult=0;
342 : int count=0;
343 0 : if ((iResult=FillIndex())<0) {
344 0 : return iResult;
345 : }
346 : // sort
347 0 : for( int i=0; i<fkNSlices*fNBorders; i++){
348 0 : std::sort(fBorderClusters+fBorderFirstCluster[i],fBorderClusters+fBorderFirstCluster[i]+fBorderNClusters[i], CompareTime);
349 : }
350 :
351 0 : for( int iBorder=0; iBorder<fkNSlices*fNBorders; iBorder+=2){
352 0 : AliBorderRecord *border1 = fBorderClusters+fBorderFirstCluster[iBorder];
353 0 : AliBorderRecord *border2 = fBorderClusters+fBorderFirstCluster[iBorder+1];
354 0 : int n1 = fBorderNClusters[iBorder];
355 0 : int n2 = fBorderNClusters[iBorder+1];
356 : /*
357 : sLeft+=n1;
358 : sRight+=n2;
359 : bool prn = (n1>0) && (n2>0);
360 : if( prn ){
361 : cout<<"Border "<<iBorder/2<<", left: "<<n1<<" clusters :"<<endl;
362 : for( int ib=0; ib<n1; ib++ ){
363 : cout<<ib<<": "<<border1[ib].fTimeBin<<endl;
364 : }
365 : cout<<"Border "<<iBorder/2<<", right: "<<n2<<" clusters :"<<endl;
366 : for( int ib=0; ib<n2; ib++ ){
367 : cout<<ib<<": "<<border2[ib].fTimeBin<<endl;
368 : }
369 : }
370 : */
371 : int ib1 = 0;
372 0 : for( int ib2 = 0; (ib1<n1) && (ib2<n2); ib2++ ){
373 : // search first cluster at border1 to merge
374 :
375 0 : while( ib1<n1 && border1[ib1].fTimeBin>border2[ib2].fTimeBin+fkMergeTimeWindow ){
376 0 : ib1++;
377 : }
378 :
379 0 : if( ib1>= n1 ) break;
380 0 : if( border1[ib1].fTimeBin < border2[ib2].fTimeBin-fkMergeTimeWindow) continue;
381 :
382 : // merge c2->c1
383 : //if( prn ) cout<<"merge "<<ib2<<"->"<<ib1<<endl;
384 0 : AliClusterRecord &c1 = fClusters[border1[ib1].fClusterRecordID];
385 0 : AliClusterRecord &c2 = fClusters[border2[ib2].fClusterRecordID];
386 :
387 0 : float w1 = c1.Cluster().fCharge;
388 0 : float w2 = c2.Cluster().fCharge;
389 0 : if( w1<=0 ) w1 = 1;
390 0 : if( w2<=0 ) w2 = 1;
391 0 : float w = 1./(w1+w2);
392 0 : w1*=w;
393 0 : w2*=w;
394 :
395 0 : c1.Cluster().fCharge += c2.Cluster().fCharge;
396 0 : if( c1.Cluster().fQMax < c2.Cluster().fQMax ) c1.Cluster().fQMax = c2.Cluster().fQMax;
397 :
398 0 : c1.Cluster().fSigmaPad2 =
399 0 : w1*c1.Cluster().fSigmaPad2 + w2*c2.Cluster().fSigmaPad2
400 0 : + (c1.Cluster().fPad - c2.Cluster().fPad)*(c1.Cluster().fPad - c2.Cluster().fPad)*w1*w2;
401 :
402 0 : c1.Cluster().fSigmaTime2 =
403 0 : w1*c1.Cluster().fSigmaTime2 + w2*c2.Cluster().fSigmaTime2
404 0 : + (c1.Cluster().fTime - c2.Cluster().fTime)*(c1.Cluster().fTime - c2.Cluster().fTime)*w1*w2;
405 :
406 0 : c1.Cluster().fPad = w1*c1.Cluster().fPad + w2*c2.Cluster().fPad;
407 0 : c1.Cluster().fTime = w1*c1.Cluster().fTime + w2*c2.Cluster().fTime;
408 :
409 0 : c1.Cluster().fFlags |= c2.Cluster().fFlags;
410 :
411 : // merge MC labels
412 :
413 0 : AliHLTTPCClusterMCWeight labels[6] = {
414 0 : c1.GetMCLabel().fClusterID[0],
415 0 : c1.GetMCLabel().fClusterID[1],
416 0 : c1.GetMCLabel().fClusterID[2],
417 0 : c2.GetMCLabel().fClusterID[0],
418 0 : c2.GetMCLabel().fClusterID[1],
419 0 : c2.GetMCLabel().fClusterID[2]
420 : };
421 :
422 0 : sort(labels, labels+6, CompareMCLabels);
423 0 : for( unsigned int i=1; i<6; i++ ){
424 0 : if(labels[i-1].fMCID==labels[i].fMCID ){
425 0 : labels[i].fWeight+=labels[i-1].fWeight;
426 0 : labels[i-1].fWeight = 0;
427 0 : }
428 : }
429 :
430 0 : sort(labels, labels+6, CompareMCWeights );
431 :
432 0 : for( unsigned int i=0; i<3; i++ ){
433 0 : c1.MCLabel().fClusterID[i] = labels[i];
434 : }
435 :
436 : // wipe c2
437 0 : fRemovedClusterIds.push_back(c2.GetId());
438 : HLTDebug("merging %d into %d", border2[ib2].fClusterRecordID, border1[ib1].fClusterRecordID);
439 : //cout<<"Set merged flag at position "<<border2[ib2].fClusterRecordID<<endl;
440 0 : c2.SetMergedTo(border1[ib1].fClusterRecordID);
441 0 : count++;
442 : // do not merge c1 anymore
443 0 : ib1++;
444 0 : }
445 : } // iBorder
446 :
447 : //cout<<"Merged "<<count<<" clusters "<<" out of "<<fClusters.size()<<endl;
448 : //sMerged+= count;
449 : // cout<<"L "<<sLeft<<" r "<<sRight<<" m "<<sMerged<<endl;
450 0 : if (iResult<0) return iResult;
451 0 : return count;
452 0 : }
453 :
|
