HLT/TPCLib/HWCFemulator/AliHLTTPCHWClusterMergerV1.cxx

   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 "AliHLTTPCTransform.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   : AliHLTLogging()
  38   , fNRows(0)
  39   , fNRowPads(0)
  40   , fNBorders(0)
  41   , fMapping(0)
  42   , fBorders()
  43   , fpData(NULL)
  44   , fRawClusterBlocks(NULL)
  45   , fMCBlocks(NULL)
  46 {
  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   fRawClusterBlocks = new AliHLTComponentBlockData* [kNPatchesTotal];
  57   fMCBlocks = new AliHLTComponentBlockData* [kNPatchesTotal];
  58   for( int i=0; i<kNPatchesTotal; i++){
  59     fRawClusterBlocks[i] = NULL;
  60     fMCBlocks[i] = NULL;
  61   }
  62 }
  63
  64 AliHLTTPCHWClusterMergerV1::~AliHLTTPCHWClusterMergerV1()
  65 {
  66   // destructor
  67   delete[] fMapping;
  68   delete[] fRawClusterBlocks;
  69   delete[] fMCBlocks;
  70 }
  71
  72 Int_t AliHLTTPCHWClusterMergerV1::Init()
  73 {
  74   // initialisation
  75
  76   Int_t iResult=0;
  77
  78   delete[] fMapping;
  79   fMapping = 0;
  80   fNRows = AliHLTTPCTransform::GetNRows();
  81   fNRowPads = 0;
  82   fNBorders = 0;
  83   for( int i=0; i<fNRows; i++ ){
  84     int nPads = AliHLTTPCTransform::GetNPads(i);
  85     if( fNRowPads<nPads ) fNRowPads = nPads;
  86   }
  87   int nPadsTotal = fNRows*fNRowPads;
  88   fMapping = new AliHLTInt16_t [nPadsTotal];
  89
  90   if( !fMapping ){
  91     HLTError("Can not allocate memory: %d bytes", nPadsTotal*sizeof(AliHLTInt16_t));
  92     return -ENOMEM;
  93   }
  94   for( int i=0; i<nPadsTotal; i++ ){
  95     fMapping[i] = -1;
  96   }
  97
  98   AliHLTTPCHWCFSupport support;
  99
 100   for( int iPart=0; iPart<6; iPart++ ){
 101     const AliHLTUInt32_t *m = support.GetMapping(0,iPart);
 102     if( !m ){
 103       HLTError("Can not get mapping for partition %d", iPart);
 104       iResult = -1;
 105       continue;
 106     }
 107     int nHWAddress=m[0];
 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;
 118     }
 119   }
 120
 121   fBorders.clear();
 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
 128
 129         m[pad] = fNBorders;
 130         for( int i=1; i<fkMergeWidth; i++ ){
 131           if( pad-i<0 || m[pad-i]!=-1 ) break;
 132           m[pad-i] = fNBorders;
 133         }
 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;
 138         }
 139         fNBorders+=2;
 140       }
 141     }
 142   }
 143
 144   Clear();
 145
 146   return iResult;
 147 }
 148
 149
 150
 151 Int_t AliHLTTPCHWClusterMergerV1::SetDataBlock(  AliHLTComponentBlockData *block)
 152 {
 153   //
 154   // set block of clusters or MC labels
 155   //
 156
 157   if( !block ){
 158     HLTError("Input NULL pointer to data block");
 159     return -1;
 160   }
 161
 162   Int_t blockType=0;
 163
 164   if( block->fDataType == (AliHLTTPCDefinitions::fgkRawClustersDataType | kAliHLTDataOriginTPC) ) blockType=1;
 165   else if( block->fDataType == (AliHLTTPCDefinitions::fgkAliHLTDataTypeClusterMCInfo | kAliHLTDataOriginTPC) ) blockType=2;
 166
 167   if( blockType==0 ){
 168     HLTError("Unexpected type of input block:  %d", block->fDataType);
 169     return -1;
 170   }
 171
 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 );
 178     return -1;
 179   }
 180
 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 );
 183     return -1;
 184   }
 185
 186   Int_t iSlicePatch = minSlice*fkNPatches + minPartition;
 187
 188   if( blockType == 1 ) fRawClusterBlocks[iSlicePatch] = block;
 189   else if( blockType == 2 )  fMCBlocks[iSlicePatch] = block;
 190
 191   return 0;
 192 }
 193
 194
 195 void AliHLTTPCHWClusterMergerV1::Clear()
 196 {
 197   /// cleanup
 198
 199   //  fClusters.~vector<AliClusterRecord>();
 200   // new(&fClusters) vector<AliClusterRecord>;
 201
 202   fpData = NULL;
 203   for( int i=0; i<fkNSlices*fkNPatches; i++){
 204     fRawClusterBlocks[i] = NULL;
 205     fMCBlocks[i] = NULL;
 206   }
 207 }
 208
 209
 210 int AliHLTTPCHWClusterMergerV1::Merge()
 211 {
 212   /// merge clusters
 213
 214   if( !fMapping ) Init();
 215   if( !fMapping ) return 0;
 216
 217   if( !fpData ){
 218     HLTError("Pointer to input data is not set");
 219     return -1;
 220   }
 221
 222
 223   int iResult = 0;
 224   int count = 0;
 225
 226   vector<AliBorderRecord> *fBorderClusters = new vector<AliBorderRecord> [fNBorders];
 227
 228   for( int iSlice=0; iSlice<fkNSlices; iSlice++){
 229
 230     for( int iBorder=0; iBorder<fNBorders; iBorder++) fBorderClusters[iBorder].clear();
 231
 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;
 238
 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);
 246         mclabels = 0;
 247         fMCBlocks[iSlicePatch] = 0;
 248       }
 249
 250       for( UInt_t iCluster=0; iCluster<clusters->fCount; iCluster++){
 251         AliHLTTPCRawCluster &cluster = clusters->fClusters[iCluster];
 252
 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;
 258
 259         int iBorder = -1;
 260
 261         if( sliceRow>=0 && sliceRow <fNRows && iPad>=0 && iPad < fNRowPads ){
 262           iBorder = fMapping[sliceRow*fNRowPads+iPad];
 263           if( iBorder>=0 ){
 264             float dPad = pad - fBorders[iBorder].fPadPosition;
 265             if( cluster.GetSigmaY2()>1.e-4 ){
 266               if( dPad*dPad > 12.*cluster.GetSigmaY2() ) iBorder = -1;
 267             } else {
 268               if( fabs(dPad)>1. ) iBorder = -1;
 269             }
 270           }
 271         }
 272         if( iBorder>=0 ){
 273           AliHLTTPCClusterMCLabel *mc=0;
 274           if( mclabels ) mc = &( mclabels->fLabels[iCluster] );
 275           fBorderClusters[iBorder].push_back( AliBorderRecord( &cluster, mc, (int)cluster.GetTime()) );
 276         }
 277       }
 278     } // patches
 279
 280
 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;
 287
 288       // sort
 289
 290       std::sort(border1.begin(),border1.end(), CompareTime);
 291       std::sort(border2.begin(),border2.end(), CompareTime);
 292
 293       // merge
 294
 295       UInt_t ib1 = 0;
 296       for( UInt_t ib2 = 0; (ib1<n1) && (ib2<n2); ib2++ ){
 297
 298         // search first cluster at border1 to merge
 299
 300         while( ib1<n1 && border1[ib1].fTimeBin>border2[ib2].fTimeBin+fkMergeTimeWindow ){
 301           ib1++;
 302         }
 303
 304         if( ib1>= n1 ) break;
 305         if( border1[ib1].fTimeBin < border2[ib2].fTimeBin-fkMergeTimeWindow) continue;
 306
 307         // merge c2->c1
 308
 309         AliHLTTPCRawCluster *c1 = border1[ib1].fCluster;
 310         AliHLTTPCRawCluster *c2 = border2[ib2].fCluster;
 311         AliHLTTPCClusterMCLabel *mc1 = border1[ib1].fMC;
 312         AliHLTTPCClusterMCLabel *mc2 = border2[ib2].fMC;
 313
 314
 315         if( c1->GetPad()<-1 || c2->GetPad()<-1 ) continue;
 316
 317         float w1 = c1->GetCharge();
 318         float w2 = c2->GetCharge();
 319         if( w1<=0 ) w1 = 1;
 320         if( w2<=0 ) w2 = 1;
 321         float w = 1./(w1+w2);
 322         w1*=w;
 323         w2*=w;
 324
 325         c1->SetCharge( c1->GetCharge() + c2->GetCharge() );
 326         if( c1->GetQMax() < c2->GetQMax() ) c1->SetQMax( c2->GetQMax() );
 327
 328         c1->SetSigmaY2(  w1*c1->GetSigmaY2() + w2*c2->GetSigmaY2()
 329                          + (c1->GetPad() - c2->GetPad())*(c1->GetPad() - c2->GetPad())*w1*w2 );
 330
 331         c1->SetSigmaZ2( w1*c1->GetSigmaZ2() + w2*c2->GetSigmaZ2()
 332                         + (c1->GetTime() - c2->GetTime())*(c1->GetTime() - c2->GetTime())*w1*w2 );
 333
 334         c1->SetPad( w1*c1->GetPad() + w2*c2->GetPad() );
 335         c1->SetTime( w1*c1->GetTime() + w2*c2->GetTime() );
 336
 337         // merge MC labels
 338         if( mc1 && mc2 ){
 339
 340           AliHLTTPCClusterMCWeight labels[6] = {
 341             mc1->fClusterID[0],
 342             mc1->fClusterID[1],
 343             mc1->fClusterID[2],
 344             mc2->fClusterID[0],
 345             mc2->fClusterID[1],
 346             mc2->fClusterID[2]
 347           };
 348
 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;
 354             }
 355           }
 356
 357           sort(labels, labels+6, CompareMCWeights );
 358
 359           for( unsigned int i=0; i<3; i++ ){
 360             mc1->fClusterID[i] = labels[i];
 361           }
 362
 363         } // MC labels
 364
 365         // wipe c2
 366         c2->SetPad(-10);
 367
 368         count++;
 369         // do not merge to c1 anymore
 370         ib1++;
 371       }
 372     } // iBorder
 373
 374
 375     // remove merged clusters from data
 376
 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;
 384
 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;
 389
 390       clusters->fCount=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 ) continue;
 395         clusters->fClusters[clusters->fCount] = cluster;
 396         if( mclabels ) mclabels->fLabels[clusters->fCount] = mclabels->fLabels[iCluster];
 397         clusters->fCount++;
 398       }
 399
 400       block->fSize = sizeof(AliHLTTPCRawClusterData) + clusters->fCount*sizeof(AliHLTTPCRawCluster);
 401       if( mclabels ){
 402         mclabels->fCount = clusters->fCount;
 403         mcblock->fSize = sizeof(AliHLTTPCClusterMCData) + mclabels->fCount*sizeof(AliHLTTPCClusterMCLabel);
 404       }
 405     }
 406
 407   } // iSlice
 408
 409   delete[] fBorderClusters;
 410
 411   if (iResult<0) return iResult;
 412   return count;
 413 }
 414