HLT/TPCLib/HWCFemulator/AliHLTTPCHWClusterMerger.cxx

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