3 //**************************************************************************
4 //* This file is property of and copyright by the ALICE HLT Project *
5 //* ALICE Experiment at CERN, All rights reserved. *
7 //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no> *
9 //* for The ALICE HLT Project. *
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 //**************************************************************************
20 // @file AliHLTTPCHWClusterMerger.cxx
21 // @author Matthias Richter, Sergey Gorbunov
23 // @brief Merger class for HLT TPC Hardware clusters
24 // Handles merging of branch border clusters
28 #include "AliHLTTPCHWClusterMerger.h"
29 #include "AliHLTTPCTransform.h"
30 #include "AliHLTTPCSpacePointData.h"
31 #include "AliHLTTPCHWCFSupport.h"
33 /** ROOT macro for the implementation of ROOT specific class methods */
34 ClassImp(AliHLTTPCHWClusterMerger)
36 AliHLTTPCHWClusterMerger::AliHLTTPCHWClusterMerger()
44 , fBorderFirstCluster(0)
46 , fBorderNClustersTotal(0)
48 , fRemovedClusterIds()
52 // see header file for class documentation
54 // refer to README to build package
56 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
60 AliHLTTPCHWClusterMerger::~AliHLTTPCHWClusterMerger()
65 delete[] fBorderNClusters;
66 delete[] fBorderFirstCluster;
67 delete[] fBorderClusters;
70 void AliHLTTPCHWClusterMerger::Init()
76 delete[] fBorderNClusters;
77 delete[] fBorderFirstCluster;
78 delete[] fBorderClusters;
82 fBorderFirstCluster = 0;
84 fBorderNClustersTotal = 0;
85 fNRows = AliHLTTPCTransform::GetNRows();
88 for( int i=0; i<fNRows; i++ ){
89 int nPads = AliHLTTPCTransform::GetNPads(i);
90 if( fNRowPads<nPads ) fNRowPads = nPads;
92 int nPadsTotal = fNRows*fNRowPads;
93 fMapping = new AliHLTInt16_t [nPadsTotal];
96 HLTError("Can not allocate memory: %d bytes", nPadsTotal*sizeof(AliHLTInt16_t));
99 for( int i=0; i<nPadsTotal; i++ ){
103 AliHLTTPCHWCFSupport support;
105 for( int iPart=0; iPart<6; iPart++ ){
106 const AliHLTUInt32_t *m = support.GetMapping(0,iPart);
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;
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 ){
127 for( int i=1; i<fkMergeWidth; i++ ){
128 if( pad-i<0 || m[pad-i]!=-1 ) break;
129 m[pad-i] = fNBorders;
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;
137 vBorders.push_back(pad+1.);
138 vBorders.push_back(pad+1.);
143 //cout<<" NBorders = "<<fNBorders/2<<endl;
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 + row*fNRowPads;
151 cout<<row<<" "<<pad<<" "<<m[pad]%2<<endl;
157 fBorders = new AliHLTFloat32_t [fNBorders];
159 HLTError("Can not allocate memory: %d bytes", fNBorders*sizeof(AliHLTFloat32_t));
163 for( int i=0; i<fNBorders; i++ ){
164 fBorders[i] = vBorders[i];
166 fBorderNClusters = new int[fkNSlices*fNBorders];
167 if( !fBorderNClusters ){
168 HLTError("Can not allocate memory: %d bytes", fkNSlices*fNBorders*sizeof(int));
172 fBorderFirstCluster = new int[fkNSlices*fNBorders];
173 if( !fBorderFirstCluster ){
174 HLTError("Can not allocate memory: %d bytes", fkNSlices*fNBorders*sizeof(int));
178 for( int i=0; i<fkNSlices*fNBorders; i++){
179 fBorderNClusters[i] = 0;
180 fBorderFirstCluster[i] = 0;
184 bool AliHLTTPCHWClusterMerger::CheckCandidate(int /*slice*/, int partition, int partitionrow, float pad, float /*time*/, float sigmaPad2
187 /// check cluster if it is a candidate for merging
188 int slicerow=partitionrow+AliHLTTPCTransform::GetFirstRow(partition);
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;
195 int atBorder = fMapping[slicerow*fNRowPads+iPad];
196 if( atBorder<0 ) return 0;
197 float dPad = pad - fBorders[atBorder];
198 if( sigmaPad2>1.e-4 ){
199 if( dPad*dPad > 12.*sigmaPad2 ) return 0;
201 if( fabs(dPad)>1. ) return 0;
206 int AliHLTTPCHWClusterMerger::AddCandidate(int slice,
213 unsigned short charge,
216 const AliHLTTPCClusterMCLabel *mc
219 /// add a candidate for merging and register in the index grid
222 int iPad = (int) pad;
223 if( !fMapping ) Init();
225 int slicerow=partitionrow+AliHLTTPCTransform::GetFirstRow(partition);
227 if (id!=~AliHLTUInt32_t(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);
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);
240 if( slice<0 || slice>=fkNSlices ){
241 HLTError("cluster id 0x%08x has wrong slice number %d", id, slice);
244 if( partition<0 || partition>=6 ){
245 HLTError("cluster id 0x%08x has wrong partition number %d", id, partition);
248 if( slice>=0 && slice<fkNSlices && slicerow>=0 && slicerow <fNRows && fMapping && iPad>=0 && iPad < fNRowPads ){
249 iBorder = fMapping[slicerow*fNRowPads+iPad];
252 float dPad = pad - fBorders[iBorder];
254 if( dPad*dPad > 12.*sigmaY2 ) iBorder = -1;
256 if( fabs(dPad)>1. ) iBorder = -1;
259 if( iBorder>=0 ) iBorder = slice*fNBorders + iBorder;
262 fClusters.push_back(AliClusterRecord(slice, partition, iBorder, -1, id,
263 AliHLTTPCRawCluster(partitionrow, pad, time, sigmaY2, sigmaZ2, charge, qmax),
264 mc!=NULL?*mc:AliHLTTPCClusterMCLabel() ));
267 fBorderNClusters[iBorder]++;
268 fBorderNClustersTotal++;
270 return fClusters.size()-1;
273 int AliHLTTPCHWClusterMerger::FillIndex()
275 /// loop over cached raw clusters and fill index
277 delete[] fBorderClusters;
279 fBorderClusters = new AliBorderRecord[fBorderNClustersTotal];
281 fBorderFirstCluster[0] = 0;
282 for(int i=1; i<fkNSlices*fNBorders; i++ ){
283 fBorderFirstCluster[i] = fBorderFirstCluster[i-1] + fBorderNClusters[i-1];
284 fBorderNClusters[i-1] = 0;
286 fBorderNClusters[fkNSlices*fNBorders-1]=0;
287 for (AliHLTUInt32_t pos=0; pos<fClusters.size(); pos++) {
288 int iBorder = fClusters[pos].GetBorder();
289 if( iBorder<0 ) continue;
290 AliBorderRecord &b = fBorderClusters[fBorderFirstCluster[iBorder]+fBorderNClusters[iBorder]];
291 b.fClusterRecordID = pos;
292 b.fTimeBin = (int)fClusters[pos].GetCluster().GetTime();
293 fBorderNClusters[iBorder]++;
296 cout<<"Border clusters: "<<endl;
297 for( int iSlice=0; iSlice<fkNSlices; iSlice++){
298 for( int ib = 0; ib<fNBorders; ib++ ){
299 int ind = iSlice*fNBorders+ib;
300 cout<<iSlice<<" "<<ib<<" :"<<fBorderFirstCluster[ind]<<", "<<fBorderNClusters[ind]<<endl;
305 for( int ib=1; ib<fkNSlices*fNBorders; ib++ ){
306 if( fBorderFirstCluster[ib] != fBorderFirstCluster[ib-1]+fBorderNClusters[ib-1] ){
307 cout<<"Something wrong with cluster borders !!! "<<endl;
314 void AliHLTTPCHWClusterMerger::Clear()
318 fRemovedClusterIds.clear();
320 delete[] fBorderClusters;
322 for( int i=0; i<fkNSlices*fNBorders; i++){
323 fBorderNClusters[i] = 0;
324 fBorderFirstCluster[i] = 0;
326 fBorderNClustersTotal = 0;
330 int AliHLTTPCHWClusterMerger::Merge()
333 /// first all candidates are filled into the index grid, looping over all clusters
334 /// in the grid automatically results in time ordered clusters per row (ordering with
335 /// respect to cells, within a cell two clusters can be reversed)
337 if( !fMapping ) Init();
338 if( !fMapping ) return 0;
339 //static int sLeft = 0, sRight = 0, sMerged = 0;
342 if ((iResult=FillIndex())<0) {
346 for( int i=0; i<fkNSlices*fNBorders; i++){
347 std::sort(fBorderClusters+fBorderFirstCluster[i],fBorderClusters+fBorderFirstCluster[i]+fBorderNClusters[i], CompareTime);
350 for( int iBorder=0; iBorder<fkNSlices*fNBorders; iBorder+=2){
351 AliBorderRecord *border1 = fBorderClusters+fBorderFirstCluster[iBorder];
352 AliBorderRecord *border2 = fBorderClusters+fBorderFirstCluster[iBorder+1];
353 int n1 = fBorderNClusters[iBorder];
354 int n2 = fBorderNClusters[iBorder+1];
358 bool prn = (n1>0) && (n2>0);
360 cout<<"Border "<<iBorder/2<<", left: "<<n1<<" clusters :"<<endl;
361 for( int ib=0; ib<n1; ib++ ){
362 cout<<ib<<": "<<border1[ib].fTimeBin<<endl;
364 cout<<"Border "<<iBorder/2<<", right: "<<n2<<" clusters :"<<endl;
365 for( int ib=0; ib<n2; ib++ ){
366 cout<<ib<<": "<<border2[ib].fTimeBin<<endl;
371 for( int ib2 = 0; (ib1<n1) && (ib2<n2); ib2++ ){
372 // search first cluster at border1 to merge
374 while( ib1<n1 && border1[ib1].fTimeBin>border2[ib2].fTimeBin+fkMergeTimeWindow ){
378 if( ib1>= n1 ) break;
379 if( border1[ib1].fTimeBin < border2[ib2].fTimeBin-fkMergeTimeWindow) continue;
382 //if( prn ) cout<<"merge "<<ib2<<"->"<<ib1<<endl;
383 AliClusterRecord &c1 = fClusters[border1[ib1].fClusterRecordID];
384 AliClusterRecord &c2 = fClusters[border2[ib2].fClusterRecordID];
386 float w1 = c1.Cluster().fCharge;
387 float w2 = c2.Cluster().fCharge;
390 float w = 1./(w1+w2);
394 c1.Cluster().fCharge += c2.Cluster().fCharge;
395 if( c1.Cluster().fQMax < c2.Cluster().fQMax ) c1.Cluster().fQMax = c2.Cluster().fQMax;
397 c1.Cluster().fSigmaY2 =
398 w1*c1.Cluster().fSigmaY2 + w2*c2.Cluster().fSigmaY2
399 + (c1.Cluster().fPad - c2.Cluster().fPad)*(c1.Cluster().fPad - c2.Cluster().fPad)*w1*w2;
401 c1.Cluster().fSigmaZ2 =
402 w1*c1.Cluster().fSigmaZ2 + w2*c2.Cluster().fSigmaZ2
403 + (c1.Cluster().fTime - c2.Cluster().fTime)*(c1.Cluster().fTime - c2.Cluster().fTime)*w1*w2;
405 c1.Cluster().fPad = w1*c1.Cluster().fPad + w2*c2.Cluster().fPad;
406 c1.Cluster().fTime = w1*c1.Cluster().fTime + w2*c2.Cluster().fTime;
410 AliHLTTPCClusterMCWeight labels[6] = {
411 c1.GetMCLabel().fClusterID[0],
412 c1.GetMCLabel().fClusterID[1],
413 c1.GetMCLabel().fClusterID[2],
414 c2.GetMCLabel().fClusterID[0],
415 c2.GetMCLabel().fClusterID[1],
416 c2.GetMCLabel().fClusterID[2]
419 sort(labels, labels+6, CompareMCLabels);
420 for( unsigned int i=1; i<6; i++ ){
421 if(labels[i-1].fMCID==labels[i].fMCID ){
422 labels[i].fWeight+=labels[i-1].fWeight;
423 labels[i-1].fWeight = 0;
427 sort(labels, labels+6, CompareMCWeights );
429 for( unsigned int i=0; i<3; i++ ){
430 c1.MCLabel().fClusterID[i] = labels[i];
434 fRemovedClusterIds.push_back(c2.GetId());
435 HLTDebug("merging %d into %d", border2[ib2].fClusterRecordID, border1[ib1].fClusterRecordID);
436 //cout<<"Set merged flag at position "<<border2[ib2].fClusterRecordID<<endl;
437 c2.SetMergedTo(border1[ib1].fClusterRecordID);
439 // do not merge c1 anymore
444 //cout<<"Merged "<<count<<" clusters "<<" out of "<<fClusters.size()<<endl;
446 // cout<<"L "<<sLeft<<" r "<<sRight<<" m "<<sMerged<<endl;
447 if (iResult<0) return iResult;