[u/mrichter/AliRoot.git] / HLT / TPCLib / HWCFemulator / AliHLTTPCHWClusterMergerV1.cxx

// $Id: AliHLTTPCHWClusterMergerV1.cxx 53494 2011-12-09 06:24:43Z sgorbuno $

//**************************************************************************
//* This file is property of and copyright by the ALICE HLT Project        * 
//* ALICE Experiment at CERN, All rights reserved.                         *
//*                                                                        *
//* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
//*                  Sergey Gorbunov                                       *
//*                  for The ALICE HLT Project.                            *
//*                                                                        *
//* Permission to use, copy, modify and distribute this software and its   *
//* documentation strictly for non-commercial purposes is hereby granted   *
//* without fee, provided that the above copyright notice appears in all   *
//* copies and that both the copyright notice and this permission notice   *
//* appear in the supporting documentation. The authors make no claims     *
//* about the suitability of this software for any purpose. It is          *
//* provided "as is" without express or implied warranty.                  *
//**************************************************************************

//  @file   AliHLTTPCHWClusterMergerV1.cxx
//  @author Matthias Richter, Sergey Gorbunov
//  @date   2011-11-25
//  @brief  Merger class for HLT TPC Hardware clusters
//          Handles merging of branch border clusters


#include "AliHLTTPCHWClusterMergerV1.h"
#include "AliHLTTPCHWCFSupport.h"
#include "AliHLTTPCDefinitions.h"
#include "AliHLTTPCTransform.h"
#include "AliHLTTPCRawCluster.h"
#include <algorithm>

/** ROOT macro for the implementation of ROOT specific class methods */

AliHLTTPCHWClusterMergerV1::AliHLTTPCHWClusterMergerV1()
  : AliHLTLogging()
  , fNRows(0)
  , fNRowPads(0)
  , fNBorders(0)
  , fMapping(0)
  , fBorders()
  , fpData(NULL)
  , fRawClusterBlocks(NULL)
  , fMCBlocks(NULL)
{
  // see header file for class documentation
  // or
  // refer to README to build package
  // or
  // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  // constructor
  
  const int kNPatchesTotal = fkNSlices*fkNPatches;

  fRawClusterBlocks = new AliHLTComponentBlockData* [kNPatchesTotal];
  fMCBlocks = new AliHLTComponentBlockData* [kNPatchesTotal];
  for( int i=0; i<kNPatchesTotal; i++){
    fRawClusterBlocks[i] = NULL;
    fMCBlocks[i] = NULL;
  }
}

AliHLTTPCHWClusterMergerV1::~AliHLTTPCHWClusterMergerV1()
{
  // destructor
  delete[] fMapping;
  delete[] fRawClusterBlocks;
  delete[] fMCBlocks;
}

Int_t AliHLTTPCHWClusterMergerV1::Init()
{
  // initialisation
  
  Int_t iResult=0;

  delete[] fMapping;
  fMapping = 0;
  fNRows = AliHLTTPCTransform::GetNRows();
  fNRowPads = 0;  
  fNBorders = 0;
  for( int i=0; i<fNRows; i++ ){
    int nPads = AliHLTTPCTransform::GetNPads(i);	
    if( fNRowPads<nPads ) fNRowPads = nPads;
  }
  int nPadsTotal = fNRows*fNRowPads;
  fMapping = new AliHLTInt16_t [nPadsTotal];

  if( !fMapping ){
    HLTError("Can not allocate memory: %d bytes", nPadsTotal*sizeof(AliHLTInt16_t));
    return -ENOMEM;
  }
  for( int i=0; i<nPadsTotal; i++ ){
    fMapping[i] = -1;
  }

  AliHLTTPCHWCFSupport support; 

  for( int iPart=0; iPart<6; iPart++ ){
    const AliHLTUInt32_t *m = support.GetMapping(0,iPart);
    if( !m ){
      HLTError("Can not get mapping for partition %d", iPart);
      iResult = -1;
      continue;
    }
    int nHWAddress=m[0];
    for( int iHW=0; iHW<nHWAddress; iHW++ ){
      AliHLTUInt32_t configWord = m[iHW+1];
      int row = (configWord>>8) & 0x3F;
      int pad =  configWord & 0xFF;
      bool border = (configWord>>14) & 0x1;
      if( !border ) continue;
      row+=AliHLTTPCTransform::GetFirstRow(iPart);
      if( row>fNRows ) continue;
      if( pad>=AliHLTTPCTransform::GetNPads(iPart) ) continue;
      fMapping[row*fNRowPads + pad] = -2-iPart;
    }
  }

  fBorders.clear();
  for( int row=0; row<fNRows; row++ ){
    for( int pad=0; pad<fNRowPads-1; pad++ ){
      AliHLTInt16_t *m = fMapping + row*fNRowPads;
      if( m[pad]<-1 && m[pad+1]<-1 ){		
	fBorders.push_back( AliBorderParam( pad+1., -m[pad]-2) ); // pad, patch
	fBorders.push_back( AliBorderParam( pad+1., -m[pad+1]-2) ); // pad, patch

	m[pad] = fNBorders;	
	for( int i=1; i<fkMergeWidth; i++ ){	  
	  if( pad-i<0 || m[pad-i]!=-1 ) break;
	  m[pad-i] = fNBorders;
	}
	m[pad+1] = fNBorders+1;
	for( int i=1; i<fkMergeWidth; i++ ){	  
	  if( pad+1+i>=fNRowPads || m[pad+1+i]!=-1 ) break;
	  m[pad+1+i] = fNBorders+1;	
	}	
	fNBorders+=2;
      }
    }
  } 

  Clear();

  return iResult;
}


Int_t AliHLTTPCHWClusterMergerV1::SetDataBlock(  AliHLTComponentBlockData *block)
{
  //
  // set block of clusters or MC labels
  //

  if( !block ){    
    HLTError("Input NULL pointer to data block");
    return -1;
  }
  
  Int_t blockType=0;
  
  if( block->fDataType == (AliHLTTPCDefinitions::fgkRawClustersDataType | kAliHLTDataOriginTPC) ) blockType=1;
  else if( block->fDataType == (AliHLTTPCDefinitions::fgkAliHLTDataTypeClusterMCInfo | kAliHLTDataOriginTPC) ) blockType=2;
  
  if( blockType==0 ){
    HLTError("Unexpected type of input block:  %d", block->fDataType);
    return -1;
  }
  
  UInt_t minSlice     = AliHLTTPCDefinitions::GetMinSliceNr(*block); 
  UInt_t minPartition = AliHLTTPCDefinitions::GetMinPatchNr(*block);
  UInt_t maxSlice     = AliHLTTPCDefinitions::GetMaxSliceNr(*block); 
  UInt_t maxPartition = AliHLTTPCDefinitions::GetMaxPatchNr(*block);
  if( maxSlice!=minSlice || maxPartition!=minPartition ){
    HLTError("Unexpected input block (slices %d-%d, patches %d-%d): only one TPC partition per block is expected", minSlice, maxSlice, minPartition, maxPartition );
    return -1;
  }

  if( minSlice>=(UInt_t)fkNSlices || minPartition >=(UInt_t)fkNPatches ){
    HLTError("Wrong Slice/Patch number of input block: slice %d, patch %d", minSlice, minPartition );
    return -1;
  }

  Int_t iSlicePatch = minSlice*fkNPatches + minPartition;
  
  if( blockType == 1 ) fRawClusterBlocks[iSlicePatch] = block;
  else if( blockType == 2 )  fMCBlocks[iSlicePatch] = block;

  return 0;
}


void AliHLTTPCHWClusterMergerV1::Clear()
{
  /// cleanup
  
  //  fClusters.~vector<AliClusterRecord>();
  // new(&fClusters) vector<AliClusterRecord>;

  fpData = NULL;
  for( int i=0; i<fkNSlices*fkNPatches; i++){
    fRawClusterBlocks[i] = NULL;
    fMCBlocks[i] = NULL;
  }
}


int AliHLTTPCHWClusterMergerV1::Merge()
{
  /// merge clusters
   
  if( !fMapping ) Init();
  if( !fMapping ) return 0;

  if( !fpData ){    
    HLTError("Pointer to input data is not set");
    return -1;
  }

  
  int iResult = 0;
  int count = 0;

  vector<AliBorderRecord> *fBorderClusters = new vector<AliBorderRecord> [fNBorders];

  for( int iSlice=0; iSlice<fkNSlices; iSlice++){
    
    for( int iBorder=0; iBorder<fNBorders; iBorder++) fBorderClusters[iBorder].clear();
    
    for( int iPatch=0; iPatch<fkNPatches; iPatch++){
      int iSlicePatch =  iSlice*fkNPatches+iPatch;
      AliHLTComponentBlockData *block = fRawClusterBlocks[iSlicePatch];
      if( !block ) continue;
      AliHLTTPCRawClusterData* clusters= (AliHLTTPCRawClusterData*)( fpData + block->fOffset);
      if(!clusters) continue;
  
      AliHLTComponentBlockData *mcblock = fMCBlocks[iSlicePatch];
      AliHLTTPCClusterMCData* mclabels = 0;
      if( mcblock ) mclabels = (AliHLTTPCClusterMCData*)(fpData + mcblock->fOffset );
      if( mclabels && mclabels->fCount != clusters->fCount ){
	HLTError("Slice %d, patch %d: Number of MC labels (%d) not equal to N clusters (%d)", iSlice, iPatch,  mclabels->fCount, clusters->fCount );
	mclabels->fCount = 0;
	mcblock->fSize = sizeof(AliHLTTPCClusterMCData);
	mclabels = 0;
	fMCBlocks[iSlicePatch] = 0;
      }
      
      for( UInt_t iCluster=0; iCluster<clusters->fCount; iCluster++){
	AliHLTTPCRawCluster &cluster = clusters->fClusters[iCluster];

	// check if the cluster is at the branch border    
	Int_t patchRow = cluster.GetPadRow();	
	int sliceRow=patchRow+AliHLTTPCTransform::GetFirstRow(iPatch);
 	float pad = cluster.GetPad();
	int iPad = (int) pad;

	int iBorder = -1;

	if( sliceRow>=0 && sliceRow <fNRows && iPad>=0 && iPad < fNRowPads ){
	  iBorder = fMapping[sliceRow*fNRowPads+iPad];    
	  if( iBorder>=0 ){
	    float dPad = pad - fBorders[iBorder].fPadPosition;
	    if( cluster.GetSigmaY2()>1.e-4 ){
	      if( dPad*dPad > 12.*cluster.GetSigmaY2() ) iBorder = -1;
	    } else {
	      if( fabs(dPad)>1. ) iBorder = -1;
	    }    
	  }
	}
	if( iBorder>=0 ){
	  AliHLTTPCClusterMCLabel *mc=0;
	  if( mclabels ) mc = &( mclabels->fLabels[iCluster] );
	  fBorderClusters[iBorder].push_back( AliBorderRecord( &cluster, mc, (int)cluster.GetTime()) );
	}
      }
    } // patches
  
    
    for( int iBorder=0; iBorder<fNBorders; iBorder+=2){      
      vector<AliBorderRecord> &border1 = fBorderClusters[iBorder];
      vector<AliBorderRecord> &border2 = fBorderClusters[iBorder+1];
      UInt_t n1 = border1.size();
      UInt_t n2 = border2.size();      
      if( n1==0 || n2==0 ) continue;

      // sort 
      
      std::sort(border1.begin(),border1.end(), CompareTime);
      std::sort(border2.begin(),border2.end(), CompareTime);
      
      // merge

      UInt_t ib1 = 0;
      for( UInt_t ib2 = 0; (ib1<n1) && (ib2<n2); ib2++ ){
	
	// search first cluster at border1 to merge
	
	while( ib1<n1 && border1[ib1].fTimeBin>border2[ib2].fTimeBin+fkMergeTimeWindow ){
	  ib1++;
	}
	
	if( ib1>= n1 ) break;
	if( border1[ib1].fTimeBin < border2[ib2].fTimeBin-fkMergeTimeWindow) continue;
      
	// merge c2->c1 
	
	AliHLTTPCRawCluster *c1 = border1[ib1].fCluster;
	AliHLTTPCRawCluster *c2 = border2[ib2].fCluster;
	AliHLTTPCClusterMCLabel *mc1 = border1[ib1].fMC;
	AliHLTTPCClusterMCLabel *mc2 = border2[ib2].fMC;


	if( c1->GetPad()<-1 || c2->GetPad()<-1 ) continue;

	float w1 = c1->GetCharge();
	float w2 = c2->GetCharge();
	if( w1<=0 ) w1 = 1;
	if( w2<=0 ) w2 = 1;
	float w = 1./(w1+w2);
	w1*=w;
	w2*=w;
	
	c1->SetCharge( c1->GetCharge() + c2->GetCharge() );
	if( c1->GetQMax() < c2->GetQMax() ) c1->SetQMax( c2->GetQMax() );
	
	c1->SetSigmaY2(  w1*c1->GetSigmaY2() + w2*c2->GetSigmaY2()
			 + (c1->GetPad() - c2->GetPad())*(c1->GetPad() - c2->GetPad())*w1*w2 );
	
	c1->SetSigmaZ2( w1*c1->GetSigmaZ2() + w2*c2->GetSigmaZ2()
			+ (c1->GetTime() - c2->GetTime())*(c1->GetTime() - c2->GetTime())*w1*w2 );
      
	c1->SetPad( w1*c1->GetPad() + w2*c2->GetPad() );
	c1->SetTime( w1*c1->GetTime() + w2*c2->GetTime() );
      
	// merge MC labels
	if( mc1 && mc2 ){

	  AliHLTTPCClusterMCWeight labels[6] = {
	    mc1->fClusterID[0],
	    mc1->fClusterID[1],
	    mc1->fClusterID[2],
	    mc2->fClusterID[0],
	    mc2->fClusterID[1],
	    mc2->fClusterID[2]
	  };

	  sort(labels, labels+6, CompareMCLabels);
	  for( unsigned int i=1; i<6; i++ ){
	    if(labels[i-1].fMCID==labels[i].fMCID ){
	      labels[i].fWeight+=labels[i-1].fWeight;
	      labels[i-1].fWeight = 0;
	    }
	  }
	
	  sort(labels, labels+6, CompareMCWeights );
	  
	  for( unsigned int i=0; i<3; i++ ){
	    mc1->fClusterID[i] = labels[i];
	  }

	} // MC labels

	// wipe c2
	c2->SetPad(-10);

	count++;
	// do not merge to c1 anymore
	ib1++;
      }    
    } // iBorder
    

    // remove merged clusters from data

    for( int iPatch=0; iPatch<fkNPatches; iPatch++){
      int iSlicePatch =  iSlice*fkNPatches+iPatch;
      AliHLTComponentBlockData *block = fRawClusterBlocks[iSlicePatch];
      if( !block ) continue;
      AliHLTTPCRawClusterData* clusters= (AliHLTTPCRawClusterData*)(fpData + block->fOffset);
      if(!clusters) continue;
      AliHLTUInt32_t nClustersOrig = clusters->fCount;
      
      AliHLTComponentBlockData *mcblock = fMCBlocks[iSlicePatch];
      AliHLTTPCClusterMCData* mclabels = 0;
      if( mcblock ) mclabels = (AliHLTTPCClusterMCData*)(fpData + mcblock->fOffset);
      if( mclabels && mclabels->fCount != nClustersOrig ) mclabels = 0;      

      clusters->fCount=0;
      for( UInt_t  iCluster=0; iCluster<nClustersOrig; iCluster++){
	AliHLTTPCRawCluster &cluster = clusters->fClusters[iCluster];
	if( cluster.GetPad()<-1 ) continue; // the cluster has been merged
	if( clusters->fCount != iCluster ){
	  clusters->fClusters[clusters->fCount] = cluster;
	  if( mclabels ) mclabels->fLabels[clusters->fCount] = mclabels->fLabels[iCluster];
	}
	clusters->fCount++;
      }
      
      block->fSize = sizeof(AliHLTTPCRawClusterData) + clusters->fCount*sizeof(AliHLTTPCRawCluster);
      if( mclabels ){
	mclabels->fCount = clusters->fCount;
	mcblock->fSize = sizeof(AliHLTTPCClusterMCData) + mclabels->fCount*sizeof(AliHLTTPCClusterMCLabel);	
      }
    }

  } // iSlice

  delete[] fBorderClusters;

  if (iResult<0) return iResult;
  return count;
}
Commit	Line	Data
0e3c57d2	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()
e43f6920	43	, fpData(NULL)
0e3c57d2	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 + rowfNRowPads;
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
e43f6920	144	Clear();
e43f6920	145
0e3c57d2	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	}
0e3c57d2	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
e43f6920	202	fpData = NULL;
0e3c57d2	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
e43f6920	217	if( !fpData ){
	218	HLTError("Pointer to input data is not set");
	219	return -1;
	220	}
	221
	222
0e3c57d2	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++){
e43f6920	233	int iSlicePatch = iSlice*fkNPatches+iPatch;
0e3c57d2	234	AliHLTComponentBlockData *block = fRawClusterBlocks[iSlicePatch];
0e3c57d2	235	if( !block ) continue;
e43f6920	236	AliHLTTPCRawClusterData* clusters= (AliHLTTPCRawClusterData*)( fpData + block->fOffset);
0e3c57d2	237	if(!clusters) continue;
	238
	239	AliHLTComponentBlockData *mcblock = fMCBlocks[iSlicePatch];
	240	AliHLTTPCClusterMCData* mclabels = 0;
e43f6920	241	if( mcblock ) mclabels = (AliHLTTPCClusterMCData*)(fpData + mcblock->fOffset );
0e3c57d2	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	}
e43f6920	249
0e3c57d2	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( dPaddPad > 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
e43f6920	281	for( int iBorder=0; iBorder<fNBorders; iBorder+=2){
0e3c57d2	282	vector<AliBorderRecord> &border1 = fBorderClusters[iBorder];
	283	vector<AliBorderRecord> &border2 = fBorderClusters[iBorder+1];
	284	UInt_t n1 = border1.size();
e43f6920	285	UInt_t n2 = border2.size();
0e3c57d2	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( w1c1->GetSigmaY2() + w2c2->GetSigmaY2()
	329	+ (c1->GetPad() - c2->GetPad())(c1->GetPad() - c2->GetPad())w1*w2 );
	330
	331	c1->SetSigmaZ2( w1c1->GetSigmaZ2() + w2c2->GetSigmaZ2()
	332	+ (c1->GetTime() - c2->GetTime())(c1->GetTime() - c2->GetTime())w1*w2 );
	333
	334	c1->SetPad( w1c1->GetPad() + w2c2->GetPad() );
	335	c1->SetTime( w1c1->GetTime() + w2c2->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
e43f6920	374
0e3c57d2	375	// remove merged clusters from data
	376
	377	for( int iPatch=0; iPatch<fkNPatches; iPatch++){
e43f6920	378	int iSlicePatch = iSlice*fkNPatches+iPatch;
0e3c57d2	379	AliHLTComponentBlockData *block = fRawClusterBlocks[iSlicePatch];
0e3c57d2	380	if( !block ) continue;
e43f6920	381	AliHLTTPCRawClusterData* clusters= (AliHLTTPCRawClusterData*)(fpData + block->fOffset);
0e3c57d2	382	if(!clusters) continue;
	383	AliHLTUInt32_t nClustersOrig = clusters->fCount;
	384
	385	AliHLTComponentBlockData *mcblock = fMCBlocks[iSlicePatch];
	386	AliHLTTPCClusterMCData* mclabels = 0;
e43f6920	387	if( mcblock ) mclabels = (AliHLTTPCClusterMCData*)(fpData + mcblock->fOffset);
0e3c57d2	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
b332ccfe	394	if( clusters->fCount != iCluster ){
	395	clusters->fClusters[clusters->fCount] = cluster;
	396	if( mclabels ) mclabels->fLabels[clusters->fCount] = mclabels->fLabels[iCluster];
	397	}
0e3c57d2	398	clusters->fCount++;
	399	}
	400
	401	block->fSize = sizeof(AliHLTTPCRawClusterData) + clusters->fCount*sizeof(AliHLTTPCRawCluster);
	402	if( mclabels ){
	403	mclabels->fCount = clusters->fCount;
	404	mcblock->fSize = sizeof(AliHLTTPCClusterMCData) + mclabels->fCount*sizeof(AliHLTTPCClusterMCLabel);
	405	}
	406	}
	407
	408	} // iSlice
	409
	410	delete[] fBorderClusters;
	411
	412	if (iResult<0) return iResult;
	413	return count;
	414	}
	415