1 // $Id: AliHLTTPCHWClusterMergerV1.cxx 53494 2011-12-09 06:24:43Z sgorbuno $
3 //**************************************************************************
4 //* This file is property of and copyright by the ALICE HLT Project *
5 //* ALICE Experiment at CERN, All rights reserved. *
7 //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no> *
9 //* for The ALICE HLT Project. *
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 //**************************************************************************
20 // @file AliHLTTPCHWClusterMergerV1.cxx
21 // @author Matthias Richter, Sergey Gorbunov
23 // @brief Merger class for HLT TPC Hardware clusters
24 // Handles merging of branch border clusters
27 #include "AliHLTTPCHWClusterMergerV1.h"
28 #include "AliHLTTPCHWCFSupport.h"
29 #include "AliHLTTPCDefinitions.h"
30 #include "AliHLTTPCTransform.h"
31 #include "AliHLTTPCRawCluster.h"
34 /** ROOT macro for the implementation of ROOT specific class methods */
36 AliHLTTPCHWClusterMergerV1::AliHLTTPCHWClusterMergerV1()
44 , fRawClusterBlocks(NULL)
47 // see header file for class documentation
49 // refer to README to build package
51 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
54 const int kNPatchesTotal = fkNSlices*fkNPatches;
56 fRawClusterBlocks = new AliHLTComponentBlockData* [kNPatchesTotal];
57 fMCBlocks = new AliHLTComponentBlockData* [kNPatchesTotal];
58 for( int i=0; i<kNPatchesTotal; i++){
59 fRawClusterBlocks[i] = NULL;
64 AliHLTTPCHWClusterMergerV1::~AliHLTTPCHWClusterMergerV1()
68 delete[] fRawClusterBlocks;
72 Int_t AliHLTTPCHWClusterMergerV1::Init()
80 fNRows = AliHLTTPCTransform::GetNRows();
83 for( int i=0; i<fNRows; i++ ){
84 int nPads = AliHLTTPCTransform::GetNPads(i);
85 if( fNRowPads<nPads ) fNRowPads = nPads;
87 int nPadsTotal = fNRows*fNRowPads;
88 fMapping = new AliHLTInt16_t [nPadsTotal];
91 HLTError("Can not allocate memory: %d bytes", nPadsTotal*sizeof(AliHLTInt16_t));
94 for( int i=0; i<nPadsTotal; i++ ){
98 AliHLTTPCHWCFSupport support;
100 for( int iPart=0; iPart<6; iPart++ ){
101 const AliHLTUInt32_t *m = support.GetMapping(0,iPart);
103 HLTError("Can not get mapping for partition %d", iPart);
108 for( int iHW=0; iHW<nHWAddress; iHW++ ){
109 AliHLTUInt32_t configWord = m[iHW+1];
110 int row = (configWord>>8) & 0x3F;
111 int pad = configWord & 0xFF;
112 bool border = (configWord>>14) & 0x1;
113 if( !border ) continue;
114 row+=AliHLTTPCTransform::GetFirstRow(iPart);
115 if( row>fNRows ) continue;
116 if( pad>=AliHLTTPCTransform::GetNPads(iPart) ) continue;
117 fMapping[row*fNRowPads + pad] = -2-iPart;
122 for( int row=0; row<fNRows; row++ ){
123 for( int pad=0; pad<fNRowPads-1; pad++ ){
124 AliHLTInt16_t *m = fMapping + row*fNRowPads;
125 if( m[pad]<-1 && m[pad+1]<-1 ){
126 fBorders.push_back( AliBorderParam( pad+1., -m[pad]-2) ); // pad, patch
127 fBorders.push_back( AliBorderParam( pad+1., -m[pad+1]-2) ); // pad, patch
130 for( int i=1; i<fkMergeWidth; i++ ){
131 if( pad-i<0 || m[pad-i]!=-1 ) break;
132 m[pad-i] = fNBorders;
134 m[pad+1] = fNBorders+1;
135 for( int i=1; i<fkMergeWidth; i++ ){
136 if( pad+1+i>=fNRowPads || m[pad+1+i]!=-1 ) break;
137 m[pad+1+i] = fNBorders+1;
151 Int_t AliHLTTPCHWClusterMergerV1::SetDataBlock( AliHLTComponentBlockData *block)
154 // set block of clusters or MC labels
158 HLTError("Input NULL pointer to data block");
164 if( block->fDataType == (AliHLTTPCDefinitions::fgkRawClustersDataType | kAliHLTDataOriginTPC) ) blockType=1;
165 else if( block->fDataType == (AliHLTTPCDefinitions::fgkAliHLTDataTypeClusterMCInfo | kAliHLTDataOriginTPC) ) blockType=2;
168 HLTError("Unexpected type of input block: %d", block->fDataType);
172 UInt_t minSlice = AliHLTTPCDefinitions::GetMinSliceNr(*block);
173 UInt_t minPartition = AliHLTTPCDefinitions::GetMinPatchNr(*block);
174 UInt_t maxSlice = AliHLTTPCDefinitions::GetMaxSliceNr(*block);
175 UInt_t maxPartition = AliHLTTPCDefinitions::GetMaxPatchNr(*block);
176 if( maxSlice!=minSlice || maxPartition!=minPartition ){
177 HLTError("Unexpected input block (slices %d-%d, patches %d-%d): only one TPC partition per block is expected", minSlice, maxSlice, minPartition, maxPartition );
181 if( minSlice>=(UInt_t)fkNSlices || minPartition >=(UInt_t)fkNPatches ){
182 HLTError("Wrong Slice/Patch number of input block: slice %d, patch %d", minSlice, minPartition );
186 Int_t iSlicePatch = minSlice*fkNPatches + minPartition;
188 if( blockType == 1 ) fRawClusterBlocks[iSlicePatch] = block;
189 else if( blockType == 2 ) fMCBlocks[iSlicePatch] = block;
195 void AliHLTTPCHWClusterMergerV1::Clear()
199 // fClusters.~vector<AliClusterRecord>();
200 // new(&fClusters) vector<AliClusterRecord>;
203 for( int i=0; i<fkNSlices*fkNPatches; i++){
204 fRawClusterBlocks[i] = NULL;
210 int AliHLTTPCHWClusterMergerV1::Merge()
214 if( !fMapping ) Init();
215 if( !fMapping ) return 0;
218 HLTError("Pointer to input data is not set");
226 vector<AliBorderRecord> *fBorderClusters = new vector<AliBorderRecord> [fNBorders];
228 for( int iSlice=0; iSlice<fkNSlices; iSlice++){
230 for( int iBorder=0; iBorder<fNBorders; iBorder++) fBorderClusters[iBorder].clear();
232 for( int iPatch=0; iPatch<fkNPatches; iPatch++){
233 int iSlicePatch = iSlice*fkNPatches+iPatch;
234 AliHLTComponentBlockData *block = fRawClusterBlocks[iSlicePatch];
235 if( !block ) continue;
236 AliHLTTPCRawClusterData* clusters= (AliHLTTPCRawClusterData*)( fpData + block->fOffset);
237 if(!clusters) continue;
239 AliHLTComponentBlockData *mcblock = fMCBlocks[iSlicePatch];
240 AliHLTTPCClusterMCData* mclabels = 0;
241 if( mcblock ) mclabels = (AliHLTTPCClusterMCData*)(fpData + mcblock->fOffset );
242 if( mclabels && mclabels->fCount != clusters->fCount ){
243 HLTError("Slice %d, patch %d: Number of MC labels (%d) not equal to N clusters (%d)", iSlice, iPatch, mclabels->fCount, clusters->fCount );
244 mclabels->fCount = 0;
245 mcblock->fSize = sizeof(AliHLTTPCClusterMCData);
247 fMCBlocks[iSlicePatch] = 0;
250 for( UInt_t iCluster=0; iCluster<clusters->fCount; iCluster++){
251 AliHLTTPCRawCluster &cluster = clusters->fClusters[iCluster];
253 // check if the cluster is at the branch border
254 Int_t patchRow = cluster.GetPadRow();
255 int sliceRow=patchRow+AliHLTTPCTransform::GetFirstRow(iPatch);
256 float pad = cluster.GetPad();
257 int iPad = (int) pad;
261 if( sliceRow>=0 && sliceRow <fNRows && iPad>=0 && iPad < fNRowPads ){
262 iBorder = fMapping[sliceRow*fNRowPads+iPad];
264 float dPad = pad - fBorders[iBorder].fPadPosition;
265 if( cluster.GetSigmaY2()>1.e-4 ){
266 if( dPad*dPad > 12.*cluster.GetSigmaY2() ) iBorder = -1;
268 if( fabs(dPad)>1. ) iBorder = -1;
273 AliHLTTPCClusterMCLabel *mc=0;
274 if( mclabels ) mc = &( mclabels->fLabels[iCluster] );
275 fBorderClusters[iBorder].push_back( AliBorderRecord( &cluster, mc, (int)cluster.GetTime()) );
281 for( int iBorder=0; iBorder<fNBorders; iBorder+=2){
282 vector<AliBorderRecord> &border1 = fBorderClusters[iBorder];
283 vector<AliBorderRecord> &border2 = fBorderClusters[iBorder+1];
284 UInt_t n1 = border1.size();
285 UInt_t n2 = border2.size();
286 if( n1==0 || n2==0 ) continue;
290 std::sort(border1.begin(),border1.end(), CompareTime);
291 std::sort(border2.begin(),border2.end(), CompareTime);
296 for( UInt_t ib2 = 0; (ib1<n1) && (ib2<n2); ib2++ ){
298 // search first cluster at border1 to merge
300 while( ib1<n1 && border1[ib1].fTimeBin>border2[ib2].fTimeBin+fkMergeTimeWindow ){
304 if( ib1>= n1 ) break;
305 if( border1[ib1].fTimeBin < border2[ib2].fTimeBin-fkMergeTimeWindow) continue;
309 AliHLTTPCRawCluster *c1 = border1[ib1].fCluster;
310 AliHLTTPCRawCluster *c2 = border2[ib2].fCluster;
311 AliHLTTPCClusterMCLabel *mc1 = border1[ib1].fMC;
312 AliHLTTPCClusterMCLabel *mc2 = border2[ib2].fMC;
315 if( c1->GetPad()<-1 || c2->GetPad()<-1 ) continue;
317 float w1 = c1->GetCharge();
318 float w2 = c2->GetCharge();
321 float w = 1./(w1+w2);
325 c1->SetCharge( c1->GetCharge() + c2->GetCharge() );
326 if( c1->GetQMax() < c2->GetQMax() ) c1->SetQMax( c2->GetQMax() );
328 c1->SetSigmaY2( w1*c1->GetSigmaY2() + w2*c2->GetSigmaY2()
329 + (c1->GetPad() - c2->GetPad())*(c1->GetPad() - c2->GetPad())*w1*w2 );
331 c1->SetSigmaZ2( w1*c1->GetSigmaZ2() + w2*c2->GetSigmaZ2()
332 + (c1->GetTime() - c2->GetTime())*(c1->GetTime() - c2->GetTime())*w1*w2 );
334 c1->SetPad( w1*c1->GetPad() + w2*c2->GetPad() );
335 c1->SetTime( w1*c1->GetTime() + w2*c2->GetTime() );
340 AliHLTTPCClusterMCWeight labels[6] = {
349 sort(labels, labels+6, CompareMCLabels);
350 for( unsigned int i=1; i<6; i++ ){
351 if(labels[i-1].fMCID==labels[i].fMCID ){
352 labels[i].fWeight+=labels[i-1].fWeight;
353 labels[i-1].fWeight = 0;
357 sort(labels, labels+6, CompareMCWeights );
359 for( unsigned int i=0; i<3; i++ ){
360 mc1->fClusterID[i] = labels[i];
369 // do not merge to c1 anymore
375 // remove merged clusters from data
377 for( int iPatch=0; iPatch<fkNPatches; iPatch++){
378 int iSlicePatch = iSlice*fkNPatches+iPatch;
379 AliHLTComponentBlockData *block = fRawClusterBlocks[iSlicePatch];
380 if( !block ) continue;
381 AliHLTTPCRawClusterData* clusters= (AliHLTTPCRawClusterData*)(fpData + block->fOffset);
382 if(!clusters) continue;
383 AliHLTUInt32_t nClustersOrig = clusters->fCount;
385 AliHLTComponentBlockData *mcblock = fMCBlocks[iSlicePatch];
386 AliHLTTPCClusterMCData* mclabels = 0;
387 if( mcblock ) mclabels = (AliHLTTPCClusterMCData*)(fpData + mcblock->fOffset);
388 if( mclabels && mclabels->fCount != nClustersOrig ) mclabels = 0;
391 for( UInt_t iCluster=0; iCluster<nClustersOrig; iCluster++){
392 AliHLTTPCRawCluster &cluster = clusters->fClusters[iCluster];
393 if( cluster.GetPad()<-1 ) continue; // the cluster has been merged
394 if( clusters->fCount != iCluster ){
395 clusters->fClusters[clusters->fCount] = cluster;
396 if( mclabels ) mclabels->fLabels[clusters->fCount] = mclabels->fLabels[iCluster];
401 block->fSize = sizeof(AliHLTTPCRawClusterData) + clusters->fCount*sizeof(AliHLTTPCRawCluster);
403 mclabels->fCount = clusters->fCount;
404 mcblock->fSize = sizeof(AliHLTTPCClusterMCData) + mclabels->fCount*sizeof(AliHLTTPCClusterMCLabel);
410 delete[] fBorderClusters;
412 if (iResult<0) return iResult;