Line data Source code
1 : // $Id: AliHLTTPCHWClusterMergerV1.cxx 53494 2011-12-09 06:24:43Z sgorbuno $
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 AliHLTTPCHWClusterMergerV1.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 :
27 : #include "AliHLTTPCHWClusterMergerV1.h"
28 : #include "AliHLTTPCHWCFSupport.h"
29 : #include "AliHLTTPCDefinitions.h"
30 : #include "AliHLTTPCGeometry.h"
31 : #include "AliHLTTPCRawCluster.h"
32 : #include <algorithm>
33 :
34 : /** ROOT macro for the implementation of ROOT specific class methods */
35 :
36 : AliHLTTPCHWClusterMergerV1::AliHLTTPCHWClusterMergerV1()
37 0 : : AliHLTLogging()
38 0 : , fNRows(0)
39 0 : , fNRowPads(0)
40 0 : , fNBorders(0)
41 0 : , fMapping(0)
42 0 : , fBorders()
43 0 : , fpData(NULL)
44 0 : , fRawClusterBlocks(NULL)
45 0 : , fMCBlocks(NULL)
46 0 : {
47 : // see header file for class documentation
48 : // or
49 : // refer to README to build package
50 : // or
51 : // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
52 : // constructor
53 :
54 : const int kNPatchesTotal = fkNSlices*fkNPatches;
55 :
56 0 : fRawClusterBlocks = new AliHLTComponentBlockData* [kNPatchesTotal];
57 0 : fMCBlocks = new AliHLTComponentBlockData* [kNPatchesTotal];
58 0 : for( int i=0; i<kNPatchesTotal; i++){
59 0 : fRawClusterBlocks[i] = NULL;
60 0 : fMCBlocks[i] = NULL;
61 : }
62 0 : }
63 :
64 : AliHLTTPCHWClusterMergerV1::~AliHLTTPCHWClusterMergerV1()
65 0 : {
66 : // destructor
67 0 : delete[] fMapping;
68 0 : delete[] fRawClusterBlocks;
69 0 : delete[] fMCBlocks;
70 0 : }
71 :
72 : Int_t AliHLTTPCHWClusterMergerV1::Init( Bool_t processingRCU2Data )
73 : {
74 : // initialisation
75 :
76 : Int_t iResult=0;
77 :
78 0 : delete[] fMapping;
79 0 : fMapping = 0;
80 0 : fNRows = AliHLTTPCGeometry::GetNRows();
81 0 : fNRowPads = 0;
82 0 : fNBorders = 0;
83 0 : for( int i=0; i<fNRows; i++ ){
84 0 : int nPads = AliHLTTPCGeometry::GetNPads(i);
85 0 : if( fNRowPads<nPads ) fNRowPads = nPads;
86 : }
87 0 : int nPadsTotal = fNRows*fNRowPads;
88 0 : fMapping = new AliHLTInt16_t [nPadsTotal];
89 :
90 0 : if( !fMapping ){
91 0 : HLTError("Can not allocate memory: %d bytes", nPadsTotal*sizeof(AliHLTInt16_t));
92 0 : return -ENOMEM;
93 : }
94 0 : for( int i=0; i<nPadsTotal; i++ ){
95 0 : fMapping[i] = -1;
96 : }
97 :
98 0 : AliHLTTPCHWCFSupport support;
99 0 : support.SetProcessingRCU2Data( processingRCU2Data );
100 :
101 0 : for( int iPart=0; iPart<6; iPart++ ){
102 0 : const AliHLTUInt32_t *m = support.GetMapping(0,iPart);
103 0 : if( !m ){
104 0 : HLTError("Can not get mapping for partition %d", iPart);
105 : iResult = -1;
106 0 : continue;
107 : }
108 0 : int nHWAddress=m[0];
109 0 : for( int iHW=0; iHW<nHWAddress; iHW++ ){
110 0 : AliHLTUInt32_t configWord = m[iHW+1];
111 0 : int row = (configWord>>8) & 0x3F;
112 0 : int pad = configWord & 0xFF;
113 0 : bool border = (configWord>>14) & 0x1;
114 0 : if( !border ) continue;
115 0 : row+=AliHLTTPCGeometry::GetFirstRow(iPart);
116 0 : if( row>=fNRows ) continue;
117 0 : if( pad>=AliHLTTPCGeometry::GetNPads(row) ) continue;
118 0 : fMapping[row*fNRowPads + pad] = -2-iPart;
119 0 : }
120 0 : }
121 :
122 0 : fBorders.clear();
123 0 : for( int row=0; row<fNRows; row++ ){
124 0 : for( int pad=0; pad<fNRowPads-1; pad++ ){
125 0 : AliHLTInt16_t *m = fMapping + row*fNRowPads;
126 0 : if( m[pad]<-1 && m[pad+1]<-1 ){
127 0 : fBorders.push_back( AliBorderParam( pad+1., -m[pad]-2) ); // pad, patch
128 0 : fBorders.push_back( AliBorderParam( pad+1., -m[pad+1]-2) ); // pad, patch
129 :
130 0 : m[pad] = fNBorders;
131 0 : for( int i=1; i<fkMergeWidth; i++ ){
132 0 : if( pad-i<0 || m[pad-i]!=-1 ) break;
133 0 : m[pad-i] = fNBorders;
134 : }
135 0 : m[pad+1] = fNBorders+1;
136 0 : for( int i=1; i<fkMergeWidth; i++ ){
137 0 : if( pad+1+i>=fNRowPads || m[pad+1+i]!=-1 ) break;
138 0 : m[pad+1+i] = fNBorders+1;
139 : }
140 0 : fNBorders+=2;
141 0 : }
142 : }
143 : }
144 :
145 0 : Clear();
146 :
147 : return iResult;
148 0 : }
149 :
150 :
151 :
152 : Int_t AliHLTTPCHWClusterMergerV1::SetDataBlock( AliHLTComponentBlockData *block)
153 : {
154 : //
155 : // set block of clusters or MC labels
156 : //
157 :
158 0 : if( !block ){
159 0 : HLTError("Input NULL pointer to data block");
160 0 : return -1;
161 : }
162 :
163 : Int_t blockType=0;
164 :
165 0 : if( block->fDataType == (AliHLTTPCDefinitions::fgkRawClustersDataType | kAliHLTDataOriginTPC) ) blockType=1;
166 0 : else if( block->fDataType == (AliHLTTPCDefinitions::fgkAliHLTDataTypeClusterMCInfo | kAliHLTDataOriginTPC) ) blockType=2;
167 :
168 0 : if( blockType==0 ){
169 0 : HLTError("Unexpected type of input block: %d", block->fDataType);
170 0 : return -1;
171 : }
172 :
173 0 : UInt_t minSlice = AliHLTTPCDefinitions::GetMinSliceNr(*block);
174 0 : UInt_t minPartition = AliHLTTPCDefinitions::GetMinPatchNr(*block);
175 0 : UInt_t maxSlice = AliHLTTPCDefinitions::GetMaxSliceNr(*block);
176 0 : UInt_t maxPartition = AliHLTTPCDefinitions::GetMaxPatchNr(*block);
177 0 : if( maxSlice!=minSlice || maxPartition!=minPartition ){
178 0 : HLTError("Unexpected input block (slices %d-%d, patches %d-%d): only one TPC partition per block is expected", minSlice, maxSlice, minPartition, maxPartition );
179 0 : return -1;
180 : }
181 :
182 0 : if( minSlice>=(UInt_t)fkNSlices || minPartition >=(UInt_t)fkNPatches ){
183 0 : HLTError("Wrong Slice/Patch number of input block: slice %d, patch %d", minSlice, minPartition );
184 0 : return -1;
185 : }
186 :
187 0 : Int_t iSlicePatch = minSlice*fkNPatches + minPartition;
188 :
189 0 : if( blockType == 1 ) fRawClusterBlocks[iSlicePatch] = block;
190 0 : else if( blockType == 2 ) fMCBlocks[iSlicePatch] = block;
191 :
192 : return 0;
193 0 : }
194 :
195 :
196 : void AliHLTTPCHWClusterMergerV1::Clear()
197 : {
198 : /// cleanup
199 :
200 : // fClusters.~vector<AliClusterRecord>();
201 : // new(&fClusters) vector<AliClusterRecord>;
202 :
203 0 : fpData = NULL;
204 0 : for( int i=0; i<fkNSlices*fkNPatches; i++){
205 0 : fRawClusterBlocks[i] = NULL;
206 0 : fMCBlocks[i] = NULL;
207 : }
208 0 : }
209 :
210 :
211 : int AliHLTTPCHWClusterMergerV1::Merge()
212 : {
213 : /// merge clusters
214 :
215 0 : if( !fMapping ){
216 0 : HLTError("Mapping is not initialised");
217 0 : return 0;
218 : }
219 :
220 0 : if( !fpData ){
221 0 : HLTError("Pointer to input data is not set");
222 0 : return -1;
223 : }
224 :
225 :
226 : int iResult = 0;
227 : int count = 0;
228 :
229 0 : vector<AliBorderRecord> *fBorderClusters = new vector<AliBorderRecord> [fNBorders];
230 :
231 0 : for( int iSlice=0; iSlice<fkNSlices; iSlice++){
232 :
233 0 : for( int iBorder=0; iBorder<fNBorders; iBorder++) fBorderClusters[iBorder].clear();
234 :
235 0 : for( int iPatch=0; iPatch<fkNPatches; iPatch++){
236 0 : int iSlicePatch = iSlice*fkNPatches+iPatch;
237 0 : AliHLTComponentBlockData *block = fRawClusterBlocks[iSlicePatch];
238 0 : if( !block ) continue;
239 0 : AliHLTTPCRawClusterData* clusters= (AliHLTTPCRawClusterData*)( fpData + block->fOffset);
240 0 : if(!clusters) continue;
241 :
242 0 : AliHLTComponentBlockData *mcblock = fMCBlocks[iSlicePatch];
243 : AliHLTTPCClusterMCData* mclabels = 0;
244 0 : if( mcblock ) mclabels = (AliHLTTPCClusterMCData*)(fpData + mcblock->fOffset );
245 0 : if( mclabels && mclabels->fCount != clusters->fCount ){
246 0 : HLTError("Slice %d, patch %d: Number of MC labels (%d) not equal to N clusters (%d)", iSlice, iPatch, mclabels->fCount, clusters->fCount );
247 0 : mclabels->fCount = 0;
248 0 : mcblock->fSize = sizeof(AliHLTTPCClusterMCData);
249 : mclabels = 0;
250 0 : fMCBlocks[iSlicePatch] = 0;
251 0 : }
252 :
253 0 : for( UInt_t iCluster=0; iCluster<clusters->fCount; iCluster++){
254 0 : AliHLTTPCRawCluster &cluster = clusters->fClusters[iCluster];
255 :
256 : // check if the cluster is at the branch border
257 0 : Int_t patchRow = cluster.GetPadRow();
258 0 : int sliceRow=patchRow+AliHLTTPCGeometry::GetFirstRow(iPatch);
259 0 : float pad = cluster.GetPad();
260 0 : int iPad = (int) pad;
261 :
262 : int iBorder = -1;
263 :
264 0 : if( sliceRow>=0 && sliceRow <fNRows && iPad>=0 && iPad < fNRowPads ){
265 0 : iBorder = fMapping[sliceRow*fNRowPads+iPad];
266 0 : if( iBorder>=0 ){
267 0 : float dPad = pad - fBorders[iBorder].fPadPosition;
268 0 : if( cluster.GetSigmaPad2()>1.e-4 ){
269 0 : if( dPad*dPad > 12.*cluster.GetSigmaPad2() ) iBorder = -1;
270 : } else {
271 0 : if( fabs(dPad)>1. ) iBorder = -1;
272 : }
273 0 : }
274 : }
275 0 : if( iBorder>=0 ){
276 : AliHLTTPCClusterMCLabel *mc=0;
277 0 : if( mclabels ) mc = &( mclabels->fLabels[iCluster] );
278 0 : fBorderClusters[iBorder].push_back( AliBorderRecord( &cluster, mc, (int)cluster.GetTime()) );
279 0 : }
280 : }
281 0 : } // patches
282 :
283 :
284 0 : for( int iBorder=0; iBorder<fNBorders; iBorder+=2){
285 0 : vector<AliBorderRecord> &border1 = fBorderClusters[iBorder];
286 0 : vector<AliBorderRecord> &border2 = fBorderClusters[iBorder+1];
287 0 : UInt_t n1 = border1.size();
288 0 : UInt_t n2 = border2.size();
289 0 : if( n1==0 || n2==0 ) continue;
290 :
291 : // sort
292 :
293 0 : std::sort(border1.begin(),border1.end(), CompareTime);
294 0 : std::sort(border2.begin(),border2.end(), CompareTime);
295 :
296 : // merge
297 :
298 : UInt_t ib1 = 0;
299 0 : for( UInt_t ib2 = 0; (ib1<n1) && (ib2<n2); ib2++ ){
300 :
301 : // search first cluster at border1 to merge
302 :
303 0 : while( ib1<n1 && border1[ib1].fTimeBin>border2[ib2].fTimeBin+fkMergeTimeWindow ){
304 0 : ib1++;
305 : }
306 :
307 0 : if( ib1>= n1 ) break;
308 0 : if( border1[ib1].fTimeBin < border2[ib2].fTimeBin-fkMergeTimeWindow) continue;
309 :
310 : // merge c2->c1
311 :
312 0 : AliHLTTPCRawCluster *c1 = border1[ib1].fCluster;
313 0 : AliHLTTPCRawCluster *c2 = border2[ib2].fCluster;
314 0 : AliHLTTPCClusterMCLabel *mc1 = border1[ib1].fMC;
315 0 : AliHLTTPCClusterMCLabel *mc2 = border2[ib2].fMC;
316 :
317 :
318 0 : if( c1->GetPad()<-1 || c2->GetPad()<-1 ) continue;
319 :
320 0 : float w1 = c1->GetCharge();
321 0 : float w2 = c2->GetCharge();
322 0 : if( w1<=0 ) w1 = 1;
323 0 : if( w2<=0 ) w2 = 1;
324 0 : float w = 1./(w1+w2);
325 0 : w1*=w;
326 0 : w2*=w;
327 :
328 0 : c1->SetCharge( c1->GetCharge() + c2->GetCharge() );
329 0 : if( c1->GetQMax() < c2->GetQMax() ) c1->SetQMax( c2->GetQMax() );
330 :
331 0 : c1->SetSigmaPad2( w1*c1->GetSigmaPad2() + w2*c2->GetSigmaPad2()
332 0 : + (c1->GetPad() - c2->GetPad())*(c1->GetPad() - c2->GetPad())*w1*w2 );
333 :
334 0 : c1->SetSigmaTime2( w1*c1->GetSigmaTime2() + w2*c2->GetSigmaTime2()
335 0 : + (c1->GetTime() - c2->GetTime())*(c1->GetTime() - c2->GetTime())*w1*w2 );
336 :
337 0 : c1->SetPad( w1*c1->GetPad() + w2*c2->GetPad() );
338 0 : c1->SetTime( w1*c1->GetTime() + w2*c2->GetTime() );
339 :
340 : // merge MC labels
341 0 : if( mc1 && mc2 ){
342 :
343 0 : AliHLTTPCClusterMCWeight labels[6] = {
344 0 : mc1->fClusterID[0],
345 0 : mc1->fClusterID[1],
346 0 : mc1->fClusterID[2],
347 0 : mc2->fClusterID[0],
348 0 : mc2->fClusterID[1],
349 0 : mc2->fClusterID[2]
350 : };
351 :
352 0 : sort(labels, labels+6, CompareMCLabels);
353 0 : for( unsigned int i=1; i<6; i++ ){
354 0 : if(labels[i-1].fMCID==labels[i].fMCID ){
355 0 : labels[i].fWeight+=labels[i-1].fWeight;
356 0 : labels[i-1].fWeight = 0;
357 0 : }
358 : }
359 :
360 0 : sort(labels, labels+6, CompareMCWeights );
361 :
362 0 : for( unsigned int i=0; i<3; i++ ){
363 0 : mc1->fClusterID[i] = labels[i];
364 : }
365 :
366 0 : } // MC labels
367 :
368 : // wipe c2
369 0 : c2->SetPad(-10);
370 :
371 0 : count++;
372 : // do not merge to c1 anymore
373 0 : ib1++;
374 0 : }
375 0 : } // iBorder
376 :
377 :
378 : // remove merged clusters from data
379 :
380 0 : for( int iPatch=0; iPatch<fkNPatches; iPatch++){
381 0 : int iSlicePatch = iSlice*fkNPatches+iPatch;
382 0 : AliHLTComponentBlockData *block = fRawClusterBlocks[iSlicePatch];
383 0 : if( !block ) continue;
384 0 : AliHLTTPCRawClusterData* clusters= (AliHLTTPCRawClusterData*)(fpData + block->fOffset);
385 0 : if(!clusters) continue;
386 0 : AliHLTUInt32_t nClustersOrig = clusters->fCount;
387 :
388 0 : AliHLTComponentBlockData *mcblock = fMCBlocks[iSlicePatch];
389 : AliHLTTPCClusterMCData* mclabels = 0;
390 0 : if( mcblock ) mclabels = (AliHLTTPCClusterMCData*)(fpData + mcblock->fOffset);
391 0 : if( mclabels && mclabels->fCount != nClustersOrig ) mclabels = 0;
392 :
393 0 : clusters->fCount=0;
394 0 : for( UInt_t iCluster=0; iCluster<nClustersOrig; iCluster++){
395 0 : AliHLTTPCRawCluster &cluster = clusters->fClusters[iCluster];
396 0 : if( cluster.GetPad()<-1 ) continue; // the cluster has been merged
397 0 : if( clusters->fCount != iCluster ){
398 0 : clusters->fClusters[clusters->fCount] = cluster;
399 0 : if( mclabels ) mclabels->fLabels[clusters->fCount] = mclabels->fLabels[iCluster];
400 : }
401 0 : clusters->fCount++;
402 0 : }
403 :
404 0 : block->fSize = sizeof(AliHLTTPCRawClusterData) + clusters->fCount*sizeof(AliHLTTPCRawCluster);
405 0 : if( mclabels ){
406 0 : mclabels->fCount = clusters->fCount;
407 0 : mcblock->fSize = sizeof(AliHLTTPCClusterMCData) + mclabels->fCount*sizeof(AliHLTTPCClusterMCLabel);
408 0 : }
409 0 : }
410 :
411 : } // iSlice
412 :
413 0 : delete[] fBorderClusters;
414 :
415 0 : if (iResult<0) return iResult;
416 0 : return count;
417 0 : }
418 :
|
