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

// $Id$

//**************************************************************************
//* 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   AliHLTTPCHWClusterMerger.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 <algorithm>

#include "AliHLTTPCHWClusterMerger.h"
#include "AliHLTTPCTransform.h"
#include "AliHLTTPCSpacePointData.h"
#include "AliHLTTPCHWCFSupport.h"

/** ROOT macro for the implementation of ROOT specific class methods */
ClassImp(AliHLTTPCHWClusterMerger)

AliHLTTPCHWClusterMerger::AliHLTTPCHWClusterMerger()
  : AliHLTLogging()
  , fMapping(0)
  , fNRows(0)
  , fNRowPads(0)
  , fNBorders(0)
  , fBorders(0)
  , fBorderNClusters(0)
  , fBorderFirstCluster(0)
  , fBorderClusters(0)
  , fBorderNClustersTotal(0)
  , fClusters()
  , fRemovedClusterIds()
  , fIter()
  , fEnd()
{
  // 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
}

AliHLTTPCHWClusterMerger::~AliHLTTPCHWClusterMerger()
{
  // destructor
  delete[] fMapping;
  delete [] fBorders;
  delete[] fBorderNClusters;
  delete[] fBorderFirstCluster;
  delete[] fBorderClusters;
}

void AliHLTTPCHWClusterMerger::Init()
{
  // initialisation
  
  delete[] fMapping;
  delete[] fBorders;
  delete[] fBorderNClusters;
  delete[] fBorderFirstCluster;
  delete[] fBorderClusters;
  fMapping = 0;
  fBorders = 0;
  fBorderNClusters = 0;
  fBorderFirstCluster = 0;
  fBorderClusters = 0;
  fBorderNClustersTotal = 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;
  }
  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);
    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;
    }
  }

  std::vector<AliHLTFloat32_t> vBorders;
  for( int row=0; row<fNRows; row++ ){
    for( int pad=0; pad<fNRowPads-1; pad++ ){
      AliHLTInt16_t *m = fMapping + row*fNRowPads;
      if( m[pad]==-2 && m[pad+1]==-2 ){	
	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;
	vBorders.push_back(pad+1.);
	vBorders.push_back(pad+1.);
      }
    }
  } 

  //cout<<" NBorders = "<<fNBorders/2<<endl;

  /*  
  cout<<"Borders:"<<endl;
  for( int row=0; row<fNRows; row++ ){
    for( int pad=0; pad<fNRowPads; pad++ ){
      AliHLTInt16_t *m = fMapping + row*fNRowPads;
      if( m[pad]>=0 ) {
	cout<<row<<" "<<pad<<" "<<m[pad]%2<<endl;
      }
    }
  }
  */

  fBorders = new AliHLTFloat32_t [fNBorders];
  if( !fBorders ){
    HLTError("Can not allocate memory: %d bytes", fNBorders*sizeof(AliHLTFloat32_t));
    fNBorders = 0;
    return;
  }
  for( int i=0; i<fNBorders; i++ ){
    fBorders[i] = vBorders[i];
  }
  fBorderNClusters = new int[fkNSlices*fNBorders];
  if( !fBorderNClusters ){
    HLTError("Can not allocate memory: %d bytes", fkNSlices*fNBorders*sizeof(int));    
    return;
  }

  fBorderFirstCluster = new int[fkNSlices*fNBorders];
  if( !fBorderFirstCluster ){
    HLTError("Can not allocate memory: %d bytes", fkNSlices*fNBorders*sizeof(int));    
    return;
  }

  for( int i=0; i<fkNSlices*fNBorders; i++){
    fBorderNClusters[i] = 0;  
    fBorderFirstCluster[i] = 0;  
  }
}

bool AliHLTTPCHWClusterMerger::CheckCandidate(int /*slice*/, int partition, int partitionrow, float pad, float /*time*/, float sigmaPad2
) 
{
  /// check cluster if it is a candidate for merging
  int slicerow=partitionrow+AliHLTTPCTransform::GetFirstRow(partition);
  
  if( !fMapping ) Init();
  if( !fMapping ) return 0;
  if( slicerow <0 || slicerow>= fNRows ) return 0;
  int iPad = (int) pad;
  if( iPad<0 || iPad>=fNRowPads ) return 0;
  int atBorder = fMapping[slicerow*fNRowPads+iPad];    
  if( atBorder<0 ) return 0;
  float dPad = pad - fBorders[atBorder];
  if( sigmaPad2>1.e-4 ){
    if( dPad*dPad > 12.*sigmaPad2 ) return 0;
  } else {
    if( fabs(dPad)>1. ) return 0;
  }
  return 1;
}

int AliHLTTPCHWClusterMerger::AddCandidate(int slice,
					   int partition,
					   short partitionrow,
					   float pad,
					   float time,
					   float sigmaY2,
					   float sigmaZ2,
					   unsigned short charge,
					   unsigned short qmax,
					   AliHLTUInt32_t id,
					   const AliHLTTPCClusterMCLabel &mc 
					   )
{
  /// add a candidate for merging and register in the index grid
    
  int iBorder = -1;
  int iPad = (int) pad;
  if( !fMapping ) Init();
 
  int slicerow=partitionrow+AliHLTTPCTransform::GetFirstRow(partition);
 
  if (id!=~AliHLTUInt32_t(0)) {
    if (slice<0) {
      slice=AliHLTTPCSpacePointData::GetSlice(id);
    } else if ((unsigned)slice!=AliHLTTPCSpacePointData::GetSlice(id)) {
      HLTError("cluster id 0x%08x is not consistent with specified slice %d", id, slice);
    }
    if (partition<0) {
      partition=AliHLTTPCSpacePointData::GetPatch(id);
    } else if ((unsigned)partition!=AliHLTTPCSpacePointData::GetPatch(id)) {
      HLTError("cluster id 0x%08x is not consistent with specified partition %d", id, partition);
    }
  }

  if( slice<0 || slice>=fkNSlices ){
    HLTError("cluster id 0x%08x has wrong slice number %d", id, slice);    
  }
  
  if( partition<0 || partition>=6 ){
    HLTError("cluster id 0x%08x has wrong partition number %d", id, partition);
  }
  
  if( slice>=0 && slice<fkNSlices && slicerow>=0 && slicerow <fNRows && fMapping && iPad>=0 && iPad < fNRowPads ){
    iBorder = fMapping[slicerow*fNRowPads+iPad];
    
    if( iBorder>=0 ){
      float dPad = pad - fBorders[iBorder];
      if( sigmaY2>1.e-4 ){
	if( dPad*dPad > 12.*sigmaY2 ) iBorder = -1;
      } else {
	if( fabs(dPad)>1. ) iBorder = -1;
      }    
    }
    if( iBorder>=0 ) iBorder = slice*fNBorders + iBorder;
  }

  fClusters.push_back(AliClusterRecord(slice, partition, iBorder, -1, id,
				       AliHLTTPCRawCluster(partitionrow, pad, time, sigmaY2, sigmaZ2, charge, qmax), mc ));

  if( iBorder>=0 ){
    fBorderNClusters[iBorder]++;
    fBorderNClustersTotal++;
  }
  return fClusters.size()-1;
}

int AliHLTTPCHWClusterMerger::FillIndex()
{
  /// loop over cached raw clusters and fill index

  delete[] fBorderClusters;
  fBorderClusters = 0;
  fBorderClusters = new AliBorderRecord[fBorderNClustersTotal];
  
  fBorderFirstCluster[0] = 0;
  for(int i=1; i<fkNSlices*fNBorders; i++ ){
    fBorderFirstCluster[i] = fBorderFirstCluster[i-1] + fBorderNClusters[i-1];
    fBorderNClusters[i-1] = 0;
  }
  fBorderNClusters[fkNSlices*fNBorders-1]=0;
  for (AliHLTUInt32_t pos=0; pos<fClusters.size(); pos++) {
    int iBorder = fClusters[pos].GetBorder();
    if( iBorder<0 ) continue;
    AliBorderRecord &b = fBorderClusters[fBorderFirstCluster[iBorder]+fBorderNClusters[iBorder]];
    b.fClusterRecordID = pos;
    b.fTimeBin = (int)fClusters[pos].GetCluster().GetTime();
    fBorderNClusters[iBorder]++;
  }
  /*
  cout<<"Border clusters: "<<endl;
  for( int iSlice=0; iSlice<fkNSlices; iSlice++){
    for( int ib = 0; ib<fNBorders; ib++ ){
      int ind = iSlice*fNBorders+ib;
      cout<<iSlice<<" "<<ib<<" :"<<fBorderFirstCluster[ind]<<", "<<fBorderNClusters[ind]<<endl;    
    }
  }
  */
  /*
  for( int ib=1; ib<fkNSlices*fNBorders; ib++ ){
    if( fBorderFirstCluster[ib] != fBorderFirstCluster[ib-1]+fBorderNClusters[ib-1] ){
      cout<<"Something wrong with cluster borders !!! "<<endl;
    }
  }
  */
  return 0;
}

void AliHLTTPCHWClusterMerger::Clear()
{
  /// cleanup
  fClusters.clear();
  fRemovedClusterIds.clear();
  
  delete[] fBorderClusters;
  fBorderClusters = 0;
  for( int i=0; i<fkNSlices*fNBorders; i++){
    fBorderNClusters[i] = 0;  
    fBorderFirstCluster[i] = 0; 
  }
  fBorderNClustersTotal = 0;
}


int AliHLTTPCHWClusterMerger::Merge()
{
  /// merge clusters
  /// first all candidates are filled into the index grid, looping over all clusters
  /// in the grid automatically results in time ordered clusters per row (ordering with
  /// respect to cells, within a cell two clusters can be reversed)
  
  if( !fMapping ) Init();
  if( !fMapping ) return 0;
  //static int sLeft = 0, sRight = 0, sMerged = 0;
  int iResult=0;
  int count=0;
  if ((iResult=FillIndex())<0) {
    return iResult;
  }
  // sort 
  for( int i=0; i<fkNSlices*fNBorders; i++){
    std::sort(fBorderClusters+fBorderFirstCluster[i],fBorderClusters+fBorderFirstCluster[i]+fBorderNClusters[i], CompareTime);
  }

  for( int iBorder=0; iBorder<fkNSlices*fNBorders; iBorder+=2){
    AliBorderRecord *border1 = fBorderClusters+fBorderFirstCluster[iBorder];
    AliBorderRecord *border2 = fBorderClusters+fBorderFirstCluster[iBorder+1];
    int n1 = fBorderNClusters[iBorder];
    int n2 = fBorderNClusters[iBorder+1];
    /*
sLeft+=n1;
    sRight+=n2;
    bool prn = (n1>0) && (n2>0);
    if( prn ){
      cout<<"Border "<<iBorder/2<<", left: "<<n1<<" clusters :"<<endl;
      for( int ib=0; ib<n1; ib++ ){
	cout<<ib<<": "<<border1[ib].fTimeBin<<endl;
      }
      cout<<"Border "<<iBorder/2<<", right: "<<n2<<" clusters :"<<endl;
      for( int ib=0; ib<n2; ib++ ){
	cout<<ib<<": "<<border2[ib].fTimeBin<<endl;
      }
    }
    */
    int ib1 = 0;
    for( int 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 
      //if( prn ) cout<<"merge "<<ib2<<"->"<<ib1<<endl;
      AliClusterRecord &c1 = fClusters[border1[ib1].fClusterRecordID];
      AliClusterRecord &c2 = fClusters[border2[ib2].fClusterRecordID];
      
      float w1 = c1.Cluster().fCharge;
      float w2 = c2.Cluster().fCharge;
      if( w1<=0 ) w1 = 1;
      if( w2<=0 ) w2 = 1;
      float w = 1./(w1+w2);
      w1*=w;
      w2*=w;
      
      c1.Cluster().fCharge += c2.Cluster().fCharge;
      c1.Cluster().fQMax += c2.Cluster().fQMax;
      
      c1.Cluster().fSigmaY2 = 
	w1*c1.Cluster().fSigmaY2 + w2*c2.Cluster().fSigmaY2
	+ (c1.Cluster().fPad - c2.Cluster().fPad)*(c1.Cluster().fPad - c2.Cluster().fPad)*w1*w2;
      
      c1.Cluster().fSigmaZ2 = 
	w1*c1.Cluster().fSigmaZ2 + w2*c2.Cluster().fSigmaZ2
	+ (c1.Cluster().fTime - c2.Cluster().fTime)*(c1.Cluster().fTime - c2.Cluster().fTime)*w1*w2;
      
      c1.Cluster().fPad  = w1*c1.Cluster().fPad + w2*c2.Cluster().fPad;
      c1.Cluster().fTime = w1*c1.Cluster().fTime + w2*c2.Cluster().fTime;
      
      // merge MC labels
      
      AliHLTTPCClusterMCWeight labels[6] = {
	c1.GetMCLabel().fClusterID[0],
	c1.GetMCLabel().fClusterID[1],
	c1.GetMCLabel().fClusterID[2],
	c2.GetMCLabel().fClusterID[0],
	c2.GetMCLabel().fClusterID[1],
	c2.GetMCLabel().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<6; i++ ){
	c1.MCLabel().fClusterID[i] = labels[i];
      }

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