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