f4f5ee257544b30e01b8a5fe80ffdaae299f04b0
[u/mrichter/AliRoot.git] / HLT / TPCLib / comp / AliHLTTPCDataCompressionComponent.cxx
1 // $Id$
2 //**************************************************************************
3 //* This file is property of and copyright by the ALICE HLT Project        * 
4 //* ALICE Experiment at CERN, All rights reserved.                         *
5 //*                                                                        *
6 //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
7 //*                  for The ALICE HLT Project.                            *
8 //*                                                                        *
9 //* Permission to use, copy, modify and distribute this software and its   *
10 //* documentation strictly for non-commercial purposes is hereby granted   *
11 //* without fee, provided that the above copyright notice appears in all   *
12 //* copies and that both the copyright notice and this permission notice   *
13 //* appear in the supporting documentation. The authors make no claims     *
14 //* about the suitability of this software for any purpose. It is          *
15 //* provided "as is" without express or implied warranty.                  *
16 //**************************************************************************
17
18 /// @file   AliHLTTPCDataCompressionComponent.cxx
19 /// @author Matthias Richter
20 /// @date   2011-08-08
21 /// @brief  TPC component for data compression
22 ///
23
24 #include "AliHLTTPCDataCompressionComponent.h"
25 #include "AliHLTTPCDefinitions.h"
26 #include "AliHLTTPCTrackGeometry.h"
27 #include "AliHLTTPCSpacePointContainer.h"
28 #include "AliHLTTPCHWCFSpacePointContainer.h"
29 #include "AliHLTGlobalBarrelTrack.h"
30 #include "AliHLTComponentBenchmark.h"
31 #include "AliHLTDataDeflaterSimple.h"
32 #include "AliHLTDataDeflaterHuffman.h"
33 #include "AliHLTTPCTransform.h"
34 #include "AliHLTTPCClusterMCData.h"
35 #include "AliRawDataHeader.h"
36 #include "AliCDBManager.h"
37 #include "AliCDBPath.h"
38 #include "AliCDBId.h"
39 #include "AliCDBMetaData.h"
40 #include "AliCDBEntry.h"
41 #include "TH1F.h"
42 #include "TFile.h"
43 #include <memory>
44
45 ClassImp(AliHLTTPCDataCompressionComponent)
46
47 AliHLTTPCDataCompressionComponent::AliHLTTPCDataCompressionComponent()
48   : AliHLTProcessor()
49   , fMode(0)
50   , fDeflaterMode(0)
51   , fMaxDeltaPad(AliHLTTPCDefinitions::GetMaxClusterDeltaPad())
52   , fMaxDeltaTime(AliHLTTPCDefinitions::GetMaxClusterDeltaTime())
53   , fRawInputClusters(NULL)
54   , fInputClusters(NULL)
55   , fTrackGrid(NULL)
56   , fSpacePointGrid(NULL)
57   , fpDataDeflater(NULL)
58   , fHistoCompFactor(NULL)
59   , fHistoResidualPad(NULL)
60   , fHistoResidualTime(NULL)
61   , fHistoClustersOnTracks(NULL)
62   , fHistoClusterRatio(NULL)
63   , fHistoTrackClusterRatio(NULL)
64   , fHistogramFile()
65   , fTrainingTableOutput()
66   , fpBenchmark(NULL)
67   , fVerbosity(0)
68 {
69 }
70
71 AliHLTTPCDataCompressionComponent::~AliHLTTPCDataCompressionComponent()
72 {
73   /// destructor
74 }
75
76
77 const char* AliHLTTPCDataCompressionComponent::GetComponentID()
78 {
79   /// inherited from AliHLTComponent: id of the component
80   return "TPCDataCompressor";
81 }
82
83
84 void AliHLTTPCDataCompressionComponent::GetInputDataTypes( AliHLTComponentDataTypeList& tgtList)
85 {
86   /// inherited from AliHLTComponent: list of data types in the vector reference
87   tgtList.clear();
88   tgtList.push_back(AliHLTTPCDefinitions::fgkHWClustersDataType);
89   tgtList.push_back(AliHLTTPCDefinitions::fgkClustersDataType);
90   tgtList.push_back(kAliHLTDataTypeTrack|kAliHLTDataOriginTPC);
91 }
92
93 AliHLTComponentDataType AliHLTTPCDataCompressionComponent::GetOutputDataType()
94 {
95   /// inherited from AliHLTComponent: output data type of the component.
96   return kAliHLTMultipleDataType;
97 }
98
99 int AliHLTTPCDataCompressionComponent::GetOutputDataTypes(AliHLTComponentDataTypeList& tgtList)
100 {
101   /// inherited from AliHLTComponent: multiple output data types of the component.
102   tgtList.clear();
103   tgtList.push_back(AliHLTTPCDefinitions::fgkRawClustersDataType);
104   return tgtList.size();
105 }
106
107 void AliHLTTPCDataCompressionComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
108 {
109   /// inherited from AliHLTComponent: output data size estimator
110   constBase=0;
111   inputMultiplier=1.3;
112 }
113
114 AliHLTComponent* AliHLTTPCDataCompressionComponent::Spawn()
115 {
116   /// inherited from AliHLTComponent: spawn function.
117   return new AliHLTTPCDataCompressionComponent;
118 }
119
120 int AliHLTTPCDataCompressionComponent::DoEvent( const AliHLTComponentEventData& evtData, 
121                                                 const AliHLTComponentBlockData* inputBlocks, 
122                                                 AliHLTComponentTriggerData& /*trigData*/,
123                                                 AliHLTUInt8_t* outputPtr,
124                                                 AliHLTUInt32_t& size,
125                                                 AliHLTComponentBlockDataList& outputBlocks )
126 {
127   /// inherited from AliHLTProcessor: data processing
128   int iResult=0;
129   AliHLTUInt32_t capacity=size;
130   size=0;
131
132   if (!IsDataEvent()) return 0;
133
134   if (!fRawInputClusters) {
135     return -ENODEV;
136   }
137
138   if (GetBenchmarkInstance()) {
139     GetBenchmarkInstance()->StartNewEvent();
140     GetBenchmarkInstance()->Start(0);
141   }
142
143   // Process an event
144   // Loop over all input blocks in the event
145   bool bHaveMC=(GetFirstInputBlock(AliHLTTPCDefinitions::fgkAliHLTDataTypeClusterMCInfo | kAliHLTDataOriginTPC))!=NULL;
146
147   const AliHLTComponentBlockData* pDesc=NULL;
148
149   AliHLTUInt8_t minSlice=0xFF, maxSlice=0xFF, minPatch=0xFF, maxPatch=0xFF;
150   AliHLTUInt32_t inputRawClusterSize=0;
151   AliHLTUInt32_t outputDataSize=0;
152   int allClusters=0;
153   int associatedClusters=0;
154
155   /// input track array
156   vector<AliHLTGlobalBarrelTrack> inputTrackArray;
157
158   if (GetBenchmarkInstance()) {
159     GetBenchmarkInstance()->Start(2);
160   }
161
162   // transformed clusters
163   if (fMode==10) { // FIXME: condition to be adjusted
164     for (pDesc=GetFirstInputBlock(AliHLTTPCDefinitions::fgkClustersDataType);
165          pDesc!=NULL; pDesc=GetNextInputBlock()) {
166       if (GetBenchmarkInstance()) {
167         GetBenchmarkInstance()->AddInput(pDesc->fSize);
168       }
169       AliHLTUInt8_t slice = 0;
170       AliHLTUInt8_t patch = 0;
171       slice = AliHLTTPCDefinitions::GetMinSliceNr( pDesc->fSpecification );
172       patch = AliHLTTPCDefinitions::GetMinPatchNr( pDesc->fSpecification );
173       if ( minSlice==0xFF || slice<minSlice )   minSlice = slice;
174       if ( maxSlice==0xFF || slice>maxSlice )   maxSlice = slice;
175       if ( minPatch==0xFF || patch<minPatch )   minPatch = patch;
176       if ( maxPatch==0xFF || patch>maxPatch )   maxPatch = patch;
177       if (fInputClusters) {
178         fInputClusters->AddInputBlock(pDesc);
179       }
180     }
181     if (GetBenchmarkInstance()) {
182       GetBenchmarkInstance()->Stop(2);
183       GetBenchmarkInstance()->Start(3);
184     }
185   }
186
187   // track data input
188   if (fMode==2) { // FIXME: condition to be adjusted
189     for (pDesc=GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginTPC);
190          pDesc!=NULL; pDesc=GetNextInputBlock()) {
191       if (GetBenchmarkInstance()) {
192         GetBenchmarkInstance()->AddInput(pDesc->fSize);
193       }
194       AliHLTUInt8_t slice = 0;
195       AliHLTUInt8_t patch = 0;
196       slice = AliHLTTPCDefinitions::GetMinSliceNr( pDesc->fSpecification );
197       patch = AliHLTTPCDefinitions::GetMinPatchNr( pDesc->fSpecification );
198       if ( minSlice==0xFF || slice<minSlice )   minSlice = slice;
199       if ( maxSlice==0xFF || slice>maxSlice )   maxSlice = slice;
200       if ( minPatch==0xFF || patch<minPatch )   minPatch = patch;
201       if ( maxPatch==0xFF || patch>maxPatch )   maxPatch = patch;
202       const AliHLTTracksData* pTracks=reinterpret_cast<const AliHLTTracksData*>(pDesc->fPtr);
203       if ((iResult=AliHLTGlobalBarrelTrack::ConvertTrackDataArray(pTracks, pDesc->fSize, inputTrackArray))<0) {
204         return iResult;
205       }
206     }
207   }
208
209   if (GetBenchmarkInstance()) {
210     GetBenchmarkInstance()->Stop(3);
211     GetBenchmarkInstance()->Start(4);
212   }
213
214   // processing
215   for (vector<AliHLTGlobalBarrelTrack>::iterator track=inputTrackArray.begin();
216        track!=inputTrackArray.end();
217        track++) {
218     int trackID=track->GetID();
219     if (trackID<0) {
220       // FIXME: error guard
221       HLTError("invalid track ID");
222       continue;
223     }
224
225     if (fVerbosity>0) {
226       UInt_t nofPoints=track->GetNumberOfPoints();
227       const UInt_t* points=track->GetPoints();
228       for (unsigned i=0; i<nofPoints; i++) {
229         int slice=AliHLTTPCSpacePointData::GetSlice(points[i]);
230         int partition=AliHLTTPCSpacePointData::GetPatch(points[i]);
231         int number=AliHLTTPCSpacePointData::GetNumber(points[i]);
232         HLTInfo("track %d point %d id 0x%08x slice %d partition %d number %d", track->GetID(), i, points[i], slice, partition, number);
233       }
234     }
235
236     AliHLTTrackGeometry* trackpoints=new AliHLTTPCTrackGeometry;
237     if (!trackpoints) continue;
238     trackpoints->SetTrackId(trackID);
239     trackpoints->CalculateTrackPoints(*track);
240     trackpoints->RegisterTrackPoints(fTrackGrid);
241     track->SetTrackGeometry(trackpoints);
242   }
243
244   for (vector<AliHLTGlobalBarrelTrack>::const_iterator track=inputTrackArray.begin();
245        track!=inputTrackArray.end();
246        track++) {
247     AliHLTTrackGeometry* trackpoints=track->GetTrackGeometry();
248     if (!trackpoints) continue;
249     trackpoints->FillTrackPoints(fTrackGrid);
250   }
251   if (fVerbosity>0) {
252     fTrackGrid->Print();
253   }
254
255   if (GetBenchmarkInstance()) {
256     GetBenchmarkInstance()->Stop(4);
257     GetBenchmarkInstance()->Start(5);
258   }
259
260   // loop over raw cluster blocks, assign to tracks and write
261   // unassigned clusters
262   for (pDesc=GetFirstInputBlock(AliHLTTPCDefinitions::fgkHWClustersDataType);
263        pDesc!=NULL; pDesc=GetNextInputBlock()) {
264     if (pDesc->fSize<=sizeof(AliRawDataHeader)) continue;
265     if (GetBenchmarkInstance()) {
266       GetBenchmarkInstance()->Start(1);
267       GetBenchmarkInstance()->AddInput(pDesc->fSize);
268     }
269     AliHLTUInt8_t slice = 0;
270     AliHLTUInt8_t patch = 0;
271     slice = AliHLTTPCDefinitions::GetMinSliceNr( pDesc->fSpecification );
272     patch = AliHLTTPCDefinitions::GetMinPatchNr( pDesc->fSpecification );
273     if ( minSlice==0xFF || slice<minSlice )     minSlice = slice;
274     if ( maxSlice==0xFF || slice>maxSlice )     maxSlice = slice;
275     if ( minPatch==0xFF || patch<minPatch )     minPatch = patch;
276     if ( maxPatch==0xFF || patch>maxPatch )     maxPatch = patch;
277     inputRawClusterSize+=pDesc->fSize;
278
279     // add the data and populate the index grid
280     fRawInputClusters->AddInputBlock(pDesc);
281     fRawInputClusters->PopulateAccessGrid(fSpacePointGrid, pDesc->fSpecification);
282     if (fVerbosity>0 && fSpacePointGrid->GetNumberOfSpacePoints()>0) {
283       HLTInfo("index grid slice %d partition %d", slice, patch);
284       fSpacePointGrid->Print();
285       for (AliHLTSpacePointContainer::AliHLTSpacePointPropertyGrid::iterator& cl=fSpacePointGrid->begin();
286            cl!=fSpacePointGrid->end(); cl++) {
287         AliHLTUInt32_t id=cl.Data().fId;
288         float row=fRawInputClusters->GetX(id);
289         float pad=fRawInputClusters->GetY(id);
290         float time=fRawInputClusters->GetZ(id);
291         HLTInfo("    cluster id 0x%08x: row %f  pad %f  time %f", id, row, pad, time);
292       }
293     }
294     if (GetBenchmarkInstance()) {
295       GetBenchmarkInstance()->Stop(1);
296       GetBenchmarkInstance()->Start(4);
297     }
298
299     // process the clusters per padrow and check the track grid
300     // for tracks crossing that particular padrow
301     if (GetBenchmarkInstance()) {
302       GetBenchmarkInstance()->Stop(4);
303       GetBenchmarkInstance()->Start(5);
304     }
305     allClusters+=fSpacePointGrid->GetNumberOfSpacePoints();
306     iResult=ProcessTrackClusters(&inputTrackArray[0], inputTrackArray.size(), fTrackGrid, fSpacePointGrid, fRawInputClusters, slice, patch);
307     int assignedInThisPartition=0;
308     if (iResult>=0) {
309       assignedInThisPartition=iResult;
310       associatedClusters+=iResult;
311     }
312     iResult=ProcessRemainingClusters(&inputTrackArray[0], inputTrackArray.size(), fTrackGrid, fSpacePointGrid, fRawInputClusters, slice, patch);
313     if (iResult>=0) {
314       if (fSpacePointGrid->GetNumberOfSpacePoints()>0) {
315         if (fVerbosity>0) HLTInfo("associated %d (%d) of %d clusters in slice %d partition %d", iResult+assignedInThisPartition, assignedInThisPartition, fSpacePointGrid->GetNumberOfSpacePoints(), slice, patch);
316       }
317       associatedClusters+=iResult;
318     }
319
320     // write all remaining clusters not yet assigned to tracks
321     // the index grid is used to write sorted in padrow
322     // FIXME: decoder index instead of data specification to be used
323     // use an external access grid to reduce allocated memory
324     // set to NULL after writing the clusters
325     fRawInputClusters->SetSpacePointPropertyGrid(pDesc->fSpecification, fSpacePointGrid);
326     iResult=fRawInputClusters->Write(outputPtr+size, capacity-size, outputBlocks, fpDataDeflater);
327     fRawInputClusters->SetSpacePointPropertyGrid(pDesc->fSpecification, NULL);
328     if (iResult>=0) {
329       size+=iResult;
330       outputDataSize+=iResult;
331       if (GetBenchmarkInstance()) GetBenchmarkInstance()->AddOutput(iResult);
332     }
333     if (GetBenchmarkInstance()) {
334       GetBenchmarkInstance()->Stop(5);
335     }
336
337     // forward MC labels
338     if (bHaveMC) {
339       // loop over input blocks and find MC block of current specification
340       unsigned mcBlock=0;
341       for (; mcBlock<evtData.fBlockCnt; mcBlock++) {
342         if (inputBlocks[mcBlock].fDataType!=(AliHLTTPCDefinitions::fgkAliHLTDataTypeClusterMCInfo | kAliHLTDataOriginTPC) ||
343             inputBlocks[mcBlock].fSpecification!=pDesc->fSpecification) {
344           continue;
345         }
346         if (size+inputBlocks[mcBlock].fSize>capacity) {
347           iResult=-ENOSPC;
348           break;
349         }
350         iResult=ForwardMCLabels(inputBlocks[mcBlock], fSpacePointGrid, outputPtr+size, capacity-size, size, outputBlocks);
351         if (iResult>0) {
352           size+=iResult;
353         }
354
355         HLTDebug("forwarded MC data block of slice %d partition %d", slice, patch);
356         break;
357       }
358       if (mcBlock==evtData.fBlockCnt) {
359         HLTWarning("no mc data found for slice %d partition %d", slice, patch);
360       }
361     }
362
363     fSpacePointGrid->Clear();
364   }
365   if (fHistoClusterRatio && allClusters>0) {
366     if (fVerbosity>0) HLTInfo("associated %d of %d clusters to tracks", associatedClusters, allClusters);
367     float ratio=associatedClusters; ratio/=allClusters;
368     fHistoClusterRatio->Fill(ratio);
369   }
370
371   // output of track model clusters
372   if (iResult>=0) {
373     iResult=WriteTrackClusters(inputTrackArray, fRawInputClusters, fpDataDeflater, outputPtr+size, capacity-size);
374     if (iResult>=0) {
375       AliHLTComponent_BlockData bd;
376       FillBlockData(bd);
377       bd.fOffset        = size;
378       bd.fSize          = iResult;
379       bd.fDataType      = AliHLTTPCDefinitions::ClusterTracksCompressedDataType();
380       bd.fSpecification = AliHLTTPCDefinitions::EncodeDataSpecification(minSlice, maxSlice, minPatch, maxPatch);
381       outputBlocks.push_back(bd);
382       size += bd.fSize;
383       outputDataSize+=bd.fSize;
384       HLTBenchmark("track data block of %d tracks: size %d", inputTrackArray.size(), bd.fSize);
385     }
386   }
387
388   fRawInputClusters->Clear();
389
390   float compressionFactor=(float)inputRawClusterSize;
391   if ((outputDataSize)>0) compressionFactor/=outputDataSize;
392   else compressionFactor=0.;
393   if (fHistoCompFactor) fHistoCompFactor->Fill(compressionFactor);
394
395   if (GetBenchmarkInstance()) {
396     GetBenchmarkInstance()->Stop(0);
397     if (fDeflaterMode!=3) {
398       HLTBenchmark("%s - compression factor %.2f", GetBenchmarkInstance()->GetStatistics(), compressionFactor);
399     } else {
400       HLTBenchmark("%s", GetBenchmarkInstance()->GetStatistics());
401     }
402   }
403
404   if (fInputClusters) {
405     fInputClusters->Clear();
406   }
407   if (fRawInputClusters) {
408     fRawInputClusters->Clear();
409   }
410   if (fTrackGrid) {
411     fTrackGrid->Clear();
412   }
413
414   return iResult;
415 }
416
417 int AliHLTTPCDataCompressionComponent::ProcessTrackClusters(AliHLTGlobalBarrelTrack* pTracks, unsigned nofTracks,
418                                                             AliHLTTrackGeometry::AliHLTTrackGrid* pTrackIndex,
419                                                             AliHLTSpacePointContainer::AliHLTSpacePointPropertyGrid* pClusterIndex,
420                                                             AliHLTSpacePointContainer* pClusters,
421                                                             int slice, int partition) const
422 {
423   // process to assigned track clusters
424   int assignedClusters=0;
425   if (!pTracks || nofTracks==0) return 0;
426
427   vector<AliHLTUInt32_t> processedTracks;
428   for (AliHLTTrackGeometry::AliHLTTrackGrid::iterator& trackId=pTrackIndex->begin(slice, partition, -1);
429        trackId!=pTrackIndex->end(); trackId++) {
430     if (find(processedTracks.begin(), processedTracks.end(), trackId.Data())!=processedTracks.end()) {
431       continue;
432     }
433     unsigned trackindex=0;
434     for (; trackindex<nofTracks; trackindex++) {
435       if ((unsigned)pTracks[trackindex].GetID()==trackId.Data()) break;
436     }
437     if (trackindex>=nofTracks) {
438       HLTError("can not find track of id %d", trackId.Data());
439       continue;
440     }
441     processedTracks.push_back(trackId.Data());
442     AliHLTGlobalBarrelTrack& track=pTracks[trackindex];
443     if (!track.GetTrackGeometry()) {
444       HLTError("can not find track geometry for track %d", trackId.Data());
445       continue;
446     }
447     AliHLTTPCTrackGeometry* pTrackPoints=dynamic_cast<AliHLTTPCTrackGeometry*>(track.GetTrackGeometry());
448     if (!pTrackPoints) {
449       HLTError("invalid track geometry type for track %d, expecting AliHLTTPCTrackGeometry", trackId.Data());
450       continue; 
451     }
452
453     UInt_t nofTrackPoints=track.GetNumberOfPoints();
454     const UInt_t* trackPoints=track.GetPoints();
455     for (unsigned i=0; i<nofTrackPoints; i++) {
456       const AliHLTUInt32_t& clusterId=trackPoints[i];
457       if (AliHLTTPCSpacePointData::GetSlice(clusterId)!=(unsigned)slice ||
458           AliHLTTPCSpacePointData::GetPatch(clusterId)!=(unsigned)partition) {
459         // not in the current partition;
460         continue;
461       }
462           
463       AliHLTSpacePointContainer::AliHLTSpacePointPropertyGrid::iterator& cl=pClusterIndex->find(AliHLTSpacePointContainer::AliHLTSpacePointProperties(clusterId));
464       if (cl==pClusterIndex->end()) {
465         HLTError("can not find cluster no 0x%08x of track %d in index grid", clusterId, track.GetID());
466         continue;
467       }
468
469       int clusterrow=(int)pClusters->GetX(clusterId);
470       float clusterpad=pClusters->GetY(clusterId);
471       float clustertime=pClusters->GetZ(clusterId);
472
473       AliHLTUInt32_t pointId=AliHLTTPCSpacePointData::GetID(slice, partition, clusterrow);
474       AliHLTTrackGeometry::AliHLTTrackPoint* point=pTrackPoints->GetRawTrackPoint(pointId);
475       if (!point) {
476         //HLTError("can not find track point slice %d partition %d padrow %d (0x%08x) of track %d", slice, partition, clusterrow, pointId, trackId.Data());
477         continue;
478       }
479       float pad=point->GetU();
480       float time=point->GetV();
481       if (TMath::Abs(clusterpad-pad)<fMaxDeltaPad &&
482           TMath::Abs(clustertime-time)<fMaxDeltaTime) {
483         // add this cluster to the track point and mark in the index grid
484         cl.Data().fTrackId=track.GetID();
485         point->AddAssociatedSpacePoint(clusterId, clusterpad-pad, clustertime-time);
486         assignedClusters++;
487       }
488     }
489   }
490   return assignedClusters;
491 }
492
493 int AliHLTTPCDataCompressionComponent::ProcessRemainingClusters(AliHLTGlobalBarrelTrack* pTracks, unsigned nofTracks,
494                                                                 AliHLTTrackGeometry::AliHLTTrackGrid* pTrackIndex,
495                                                                 AliHLTSpacePointContainer::AliHLTSpacePointPropertyGrid* pClusterIndex,
496                                                                 AliHLTSpacePointContainer* pClusters,
497                                                                 int slice, int partition) const
498 {
499   // assign remaining clusters to tracks
500   int iResult=0;
501   int associatedClusters=0;
502   if (!pTracks || nofTracks==0) return 0;
503
504   for (int padrow=0; padrow<AliHLTTPCTransform::GetNRows(partition); padrow++) {
505     for (AliHLTTrackGeometry::AliHLTTrackGrid::iterator& trackId=pTrackIndex->begin(slice, partition, padrow);
506          trackId!=pTrackIndex->end(); trackId++) {
507       unsigned i=0;
508       for (; i<nofTracks; i++) {
509         if ((unsigned)pTracks[i].GetID()==trackId.Data()) break;
510       }
511       if (i>=nofTracks) {
512         HLTError("can not find track of id %d", trackId.Data());
513         continue;
514       }
515       AliHLTGlobalBarrelTrack& track=pTracks[i];
516       if (!track.GetTrackGeometry()) {
517         HLTError("can not find track geometry for track %d", trackId.Data());
518         continue;
519       }
520       AliHLTTPCTrackGeometry* pTrackPoints=dynamic_cast<AliHLTTPCTrackGeometry*>(track.GetTrackGeometry());
521       if (!pTrackPoints) {
522         HLTError("invalid track geometry type for track %d, expecting AliHLTTPCTrackGeometry", trackId.Data());
523         continue;       
524       }
525       AliHLTUInt32_t pointId=AliHLTTPCSpacePointData::GetID(slice, partition, padrow);
526       AliHLTTrackGeometry::AliHLTTrackPoint* point=pTrackPoints->GetRawTrackPoint(pointId);
527       if (!point) {
528         //HLTError("can not find track point slice %d partition %d padrow %d (0x%08x) of track %d", slice, partition, padrow, pointId, trackId.Data());
529         continue;
530       }
531       float pad=point->GetU();
532       float time=point->GetV();
533
534       iResult=FindCellClusters(trackId.Data(), padrow, pad, time, pClusterIndex, pClusters, point);
535       if (iResult>0) associatedClusters+=iResult;
536       if (fVerbosity>0) {
537         HLTInfo("trackpoint track %d slice %d partition %d padrow %d: %.3f \t%.3f - associated %d", track.GetID(), slice, partition, padrow, pad, time, iResult);
538       }
539     }
540   }
541   if (iResult<0) return iResult;
542   return associatedClusters;
543 }
544
545 int AliHLTTPCDataCompressionComponent::FindCellClusters(int trackId, int padrow, float pad, float time,
546                                                         AliHLTSpacePointContainer::AliHLTSpacePointPropertyGrid* pClusterIndex,
547                                                         AliHLTSpacePointContainer* pClusters,
548                                                         AliHLTTrackGeometry::AliHLTTrackPoint* pTrackPoint) const
549 {
550   // check index cell for entries and assign to track
551   int count=0;
552   // search a 4x4 matrix out of the 9x9 matrix around the cell addressed by
553   // pad and time
554   int rowindex=pClusterIndex->GetYIndex((float)padrow);
555   int padindex=pClusterIndex->GetYIndex(pad);
556   int timeindex=pClusterIndex->GetZIndex(time);
557   if (pClusterIndex->GetCenterY(padindex)>pad) padindex--;
558   if (pClusterIndex->GetCenterZ(timeindex)>pad) timeindex--;
559   for (int padcount=0; padcount<2; padcount++, padindex++) {
560     if (padindex<0) continue;
561     if (padindex>=pClusterIndex->GetDimensionY()) break;
562     for (int timecount=0; timecount<2; timecount++, timeindex++) {
563       if (timeindex<0) continue;
564       if (timeindex>=pClusterIndex->GetDimensionZ()) break;
565       int cellindex=pClusterIndex->Index(rowindex, padindex, timeindex);
566       pad=pClusterIndex->GetCenterY(cellindex);
567       time=pClusterIndex->GetCenterZ(cellindex);
568       for (AliHLTSpacePointContainer::AliHLTSpacePointPropertyGrid::iterator& cl=pClusterIndex->begin((float)padrow, pad, time);
569            cl!=pClusterIndex->end(); cl++) {
570         if (cl.Data().fTrackId>=0) continue;
571         float clusterpad=pClusters->GetY(cl.Data().fId);
572         float clustertime=pClusters->GetZ(cl.Data().fId);
573         if (TMath::Abs(clusterpad-pad)<fMaxDeltaPad &&
574             TMath::Abs(clustertime-time)<fMaxDeltaTime) {
575           // add this cluster to the track point and mark in the index grid
576           cl.Data().fTrackId=trackId;
577           pTrackPoint->AddAssociatedSpacePoint(cl.Data().fId, clusterpad-pad, clustertime-time);
578           count++;
579         }
580       }
581     }
582   }
583   return count;
584 }
585
586 int AliHLTTPCDataCompressionComponent::WriteTrackClusters(const vector<AliHLTGlobalBarrelTrack>& tracks,
587                                                           AliHLTSpacePointContainer* pSpacePoints,
588                                                           AliHLTDataDeflater* pDeflater,
589                                                           AliHLTUInt8_t* outputPtr,
590                                                           AliHLTUInt32_t capacity) const
591 {
592   // write the track data block including all associated clusters
593   AliHLTUInt32_t size=0;
594   for (vector<AliHLTGlobalBarrelTrack>::const_iterator track=tracks.begin();
595        track!=tracks.end();
596        track++) {
597     if (!track->GetTrackGeometry()) {
598       HLTError("can not find track geometry for track %d", track->GetID());
599       return -EBADF;
600     }
601     AliHLTTPCTrackGeometry* pTrackPoints=dynamic_cast<AliHLTTPCTrackGeometry*>(track->GetTrackGeometry());
602     if (!pTrackPoints) {
603       HLTError("invalid track geometry type for track %d, expecting AliHLTTPCTrackGeometry", track->GetID());
604       return -EBADF;
605     }
606
607     int result=pTrackPoints->Write(*track, pSpacePoints, pDeflater, outputPtr+size, capacity-size);
608     if (result<0) return result;
609     size+=result;
610
611     UInt_t nofTrackPoints=track->GetNumberOfPoints();
612     const UInt_t* trackPoints=track->GetPoints();
613
614     int assignedPoints=0;
615     int assignedTrackPoints=0;
616     const vector<AliHLTTrackGeometry::AliHLTTrackPoint>& rawPoints=pTrackPoints->GetRawPoints();
617     for (vector<AliHLTTrackGeometry::AliHLTTrackPoint>::const_iterator point=rawPoints.begin();
618          point!=rawPoints.end(); point++) {
619       const vector<AliHLTTrackGeometry::AliHLTTrackSpacepoint>& spacePoints=point->GetSpacepoints();
620       for (vector<AliHLTTrackGeometry::AliHLTTrackSpacepoint>::const_iterator spacePoint=spacePoints.begin();
621            spacePoint!=spacePoints.end(); spacePoint++) {
622         float dpad=spacePoint->GetResidual(0);
623         float dtime=spacePoint->GetResidual(1);
624         if (dpad>-1000 && dtime>-1000 && fHistoResidualPad && fHistoResidualTime) {
625           fHistoResidualPad->Fill(dpad);
626           fHistoResidualTime->Fill(dtime);
627         }
628         assignedPoints++;
629         for (unsigned i=0; i<nofTrackPoints; i++) {
630           if (trackPoints[i]==spacePoint->fId) {
631             assignedTrackPoints++;
632             break;
633           }
634         }
635       }
636     }
637     if (fHistoClustersOnTracks) {
638       fHistoClustersOnTracks->Fill(assignedPoints);
639     }
640     if (fHistoTrackClusterRatio && nofTrackPoints>0) {
641       float ratio=assignedTrackPoints; ratio/=nofTrackPoints;
642       fHistoTrackClusterRatio->Fill(ratio);
643     }
644   }
645   return size;
646 }
647
648 int AliHLTTPCDataCompressionComponent::DoInit( int argc, const char** argv )
649 {
650   /// inherited from AliHLTComponent: component initialisation and argument scan.
651   int iResult=0;
652
653   // component configuration
654   //Stage 1: default initialization.
655   //Default values.
656
657   //Stage 2: OCDB.
658   TString cdbPath("HLT/ConfigTPC/");
659   cdbPath += GetComponentID();
660   //
661   iResult = ConfigureFromCDBTObjString(cdbPath);
662   if (iResult < 0) 
663     return iResult;
664
665   //Stage 3: command line arguments.
666   if (argc && (iResult = ConfigureFromArgumentString(argc, argv)) < 0)
667     return iResult;
668
669   std::auto_ptr<AliHLTComponentBenchmark> benchmark(new AliHLTComponentBenchmark);
670   if (benchmark.get()) {
671     benchmark->SetTimer(0,"total");
672     benchmark->SetTimer(1,"rawclusterinput");
673     benchmark->SetTimer(2,"clusterinput");
674     benchmark->SetTimer(3,"trackinput");
675     benchmark->SetTimer(4,"processing");
676     benchmark->SetTimer(5,"output");
677   } else {
678     return -ENOMEM;
679   }
680
681   unsigned spacePointContainerMode=(fMode==2)?AliHLTTPCHWCFSpacePointContainer::kModeCreateMap:0;
682   std::auto_ptr<AliHLTTPCHWCFSpacePointContainer> rawInputClusters(new AliHLTTPCHWCFSpacePointContainer(spacePointContainerMode));
683   std::auto_ptr<AliHLTTPCSpacePointContainer> inputClusters(new AliHLTTPCSpacePointContainer);
684
685   std::auto_ptr<TH1F> histoCompFactor(new TH1F("CompressionFactor",
686                                                "HLT TPC data compression factor",
687                                                100, 0., 10.));
688   std::auto_ptr<TH1F> histoResidualPad(new TH1F("PadResidual",
689                                                 "HLT TPC pad residual",
690                                                 100, -fMaxDeltaPad, fMaxDeltaPad));
691   std::auto_ptr<TH1F> histoResidualTime(new TH1F("TimeResidual",
692                                                  "HLT TPC time residual",
693                                                  100, -fMaxDeltaTime, fMaxDeltaTime));
694   std::auto_ptr<TH1F> histoClustersOnTracks(new TH1F("ClustersOnTracks",
695                                                  "Clusters in track model compression",
696                                                  200, 0., 600));
697   std::auto_ptr<TH1F> histoClusterRatio(new TH1F("ClusterRatio",
698                                                  "Fraction of clusters in track model compression",
699                                                  100, 0., 1.));
700   std::auto_ptr<TH1F> histoTrackClusterRatio(new TH1F("UsedTrackClusters",
701                                                  "Fraction of track clusters in track model compression",
702                                                  100, 0., 1.));
703
704   // track grid: 36 slices, each 6 partitions with max 33 rows
705   fTrackGrid=new AliHLTTrackGeometry::AliHLTTrackGrid(36, 1, 6, 1, 33, 1, 20000);
706   fSpacePointGrid=AliHLTTPCHWCFSpacePointContainer::AllocateIndexGrid();
707
708   if (!rawInputClusters.get() ||
709       !inputClusters.get() ||
710       !fTrackGrid ||
711       !fSpacePointGrid) {
712     if (fTrackGrid) delete fTrackGrid; fTrackGrid=NULL;
713     if (fSpacePointGrid) delete fSpacePointGrid; fSpacePointGrid=NULL;
714     return -ENOMEM;
715   }
716
717   if (fDeflaterMode>0 && (iResult=InitDeflater(fDeflaterMode))<0)
718     return iResult;
719
720   fpBenchmark=benchmark.release();
721   fRawInputClusters=rawInputClusters.release();
722   fInputClusters=inputClusters.release();
723
724   // initialize the histograms if stored at the end
725   // condition might be extended
726   if (!fHistogramFile.IsNull()) {
727     fHistoCompFactor=histoCompFactor.release();
728     fHistoResidualPad=histoResidualPad.release();
729     fHistoResidualTime=histoResidualTime.release();
730     fHistoClustersOnTracks=histoClustersOnTracks.release();
731     fHistoClusterRatio=histoClusterRatio.release();
732     fHistoTrackClusterRatio=histoTrackClusterRatio.release();
733   }
734
735   return iResult;
736 }
737
738 int AliHLTTPCDataCompressionComponent::InitDeflater(int mode)
739 {
740   /// init the data deflater
741   int iResult=0;
742   if (mode==2 || mode==3) {
743     // huffman deflater
744     std::auto_ptr<AliHLTDataDeflaterHuffman> deflater(new AliHLTDataDeflaterHuffman(mode==3));
745     if (!deflater.get()) return -ENOMEM;
746
747     if (!deflater->IsTrainingMode()) {
748       TString cdbPath("HLT/ConfigTPC/");
749       cdbPath += GetComponentID();
750       cdbPath += "HuffmanTables";
751       TObject* pConf=LoadAndExtractOCDBObject(cdbPath);
752       if (!pConf) return -ENOENT;
753       if (dynamic_cast<TList*>(pConf)==NULL) {
754         HLTError("huffman table configuration object of inconsistent type");
755         return -EINVAL;
756       }
757       iResult=deflater->InitDecoders(dynamic_cast<TList*>(pConf));
758       if (iResult<0) return iResult;
759     }
760     
761     unsigned nofParameters=AliHLTTPCDefinitions::GetNumberOfClusterParameterDefinitions();
762     unsigned p=0;
763     for (; p<nofParameters; p++) {
764       const AliHLTTPCDefinitions::AliClusterParameter& parameter=AliHLTTPCDefinitions::fgkClusterParameterDefinitions[p];
765       if (deflater->AddParameterDefinition(parameter.fName,
766                                            parameter.fBitLength)!=(int)parameter.fId) {
767         // for performance reason the parameter id is simply used as index in the array of
768         // definitions, the position must match the id
769         HLTFatal("mismatch between parameter id and position in array for parameter %s, rearrange definitions!", parameter.fName);
770         return -EFAULT;
771       }
772     }
773     fpDataDeflater=deflater.release();
774     return 0;
775   }
776   if (mode==1) {
777     std::auto_ptr<AliHLTDataDeflaterSimple> deflater(new AliHLTDataDeflaterSimple);
778     if (!deflater.get()) return -ENOMEM;
779
780     unsigned nofParameters=AliHLTTPCDefinitions::GetNumberOfClusterParameterDefinitions();
781     unsigned p=0;
782     for (; p<nofParameters; p++) {
783       const AliHLTTPCDefinitions::AliClusterParameter& parameter=AliHLTTPCDefinitions::fgkClusterParameterDefinitions[p];
784       if (deflater->AddParameterDefinition(parameter.fName,
785                                            parameter.fBitLength,
786                                            parameter.fOptional)!=(int)parameter.fId) {
787         // for performance reason the parameter id is simply used as index in the array of
788         // definitions, the position must match the id
789         HLTFatal("mismatch between parameter id and position in array for parameter %s, rearrange definitions!", parameter.fName);
790         return -EFAULT;
791       }
792     }
793     fpDataDeflater=deflater.release();
794     return 0;
795   }
796   HLTError("invalid deflater mode %d, allowed 1=simple 2=huffman", mode);
797   return -EINVAL;
798 }
799
800 int AliHLTTPCDataCompressionComponent::DoDeinit()
801 {
802   /// inherited from AliHLTComponent: component cleanup
803   int iResult=0;
804   if (fpBenchmark) delete fpBenchmark; fpBenchmark=NULL;
805   if (fRawInputClusters) delete fRawInputClusters; fRawInputClusters=NULL;
806   if (fInputClusters) delete fInputClusters; fInputClusters=NULL;
807   if (!fHistogramFile.IsNull()) {
808     TFile out(fHistogramFile, "RECREATE");
809     if (!out.IsZombie()) {
810       out.cd();
811       if (fHistoCompFactor) fHistoCompFactor->Write();
812       if (fHistoResidualPad) fHistoResidualPad->Write();
813       if (fHistoResidualTime) fHistoResidualTime->Write();
814       if (fHistoClusterRatio) fHistoClusterRatio->Write();
815       if (fHistoClustersOnTracks) fHistoClustersOnTracks->Write();
816       if (fHistoTrackClusterRatio) fHistoTrackClusterRatio->Write();
817       out.Close();
818     }
819   }
820   if (fHistoCompFactor) delete fHistoCompFactor;
821   fHistoCompFactor=NULL;
822   if (fHistoResidualPad) delete fHistoResidualPad;
823   fHistoResidualPad=NULL;
824   if (fHistoResidualTime) delete fHistoResidualTime;
825   fHistoResidualTime=NULL;
826   if (fHistoClustersOnTracks) delete fHistoClustersOnTracks;
827   fHistoClustersOnTracks=NULL;
828   if (fHistoClusterRatio) delete fHistoClusterRatio;
829   fHistoClusterRatio=NULL;
830   if (fHistoTrackClusterRatio) delete fHistoTrackClusterRatio;
831   fHistoTrackClusterRatio=NULL;
832
833   if (fpDataDeflater) {
834     if (fDeflaterMode==3) {
835       if (fTrainingTableOutput.IsNull()) {
836         fTrainingTableOutput=GetComponentID();
837         fTrainingTableOutput+="-huffman.root";
838       }
839       // TODO: currently, the code tables are also calculated in FindObject
840       // check if a different function is more appropriate
841       TObject* pConf=fpDataDeflater->FindObject("DeflaterConfiguration");
842       if (pConf) {
843         TString cdbEntryPath("HLT/ConfigTPC/");
844         cdbEntryPath += GetComponentID();
845         cdbEntryPath += "HuffmanTables";
846         AliCDBPath cdbPath(cdbEntryPath);
847         AliCDBId cdbId(cdbPath, AliCDBManager::Instance()->GetRun(), AliCDBRunRange::Infinity(), 0, 0);
848         AliCDBMetaData* cdbMetaData=new AliCDBMetaData;
849         cdbMetaData->SetResponsible("ALICE HLT Matthias.Richter@cern.ch");
850         cdbMetaData->SetComment("Huffman encoder configuration");
851         AliCDBEntry* entry=new AliCDBEntry(pConf, cdbId, cdbMetaData, kTRUE);
852
853         entry->SaveAs(fTrainingTableOutput);
854         // this is a small memory leak
855         // seg fault in ROOT object handling if the two objects are deleted
856         // investigate later
857         //delete entry;
858         //delete cdbMetaData;
859       }
860     }
861     delete fpDataDeflater;
862   }
863   fpDataDeflater=NULL;
864
865
866   if (fTrackGrid) delete fTrackGrid; fTrackGrid=NULL;
867   if (fSpacePointGrid) delete fSpacePointGrid; fSpacePointGrid=NULL;
868
869   return iResult;
870 }
871
872 int AliHLTTPCDataCompressionComponent::ScanConfigurationArgument(int argc, const char** argv)
873 {
874   /// inherited from AliHLTComponent: argument scan
875   int iResult=0;
876   if (argc<1) return 0;
877   int bMissingParam=0;
878   int i=0;
879   TString argument=argv[i];
880
881   do {
882     // -mode
883     if (argument.CompareTo("-mode")==0) {
884       if ((bMissingParam=(++i>=argc))) break;
885       TString parameter=argv[i];
886       if (parameter.IsDigit()) {
887         fMode=parameter.Atoi();
888         return 2;
889       } else {
890         HLTError("invalid parameter for argument %s, expecting number instead of %s", argument.Data(), parameter.Data());
891         return -EPROTO;
892       }
893     }
894
895     // -deflater-mode
896     if (argument.CompareTo("-deflater-mode")==0) {
897       if ((bMissingParam=(++i>=argc))) break;
898       TString parameter=argv[i];
899       if (parameter.IsDigit()) {
900         fDeflaterMode=parameter.Atoi();
901         return 2;
902       } else {
903         HLTError("invalid parameter for argument %s, expecting number instead of %s", argument.Data(), parameter.Data());
904         return -EPROTO;
905       }
906     }
907
908     // -histogram-file
909     if (argument.CompareTo("-histogram-file")==0) {
910       if ((bMissingParam=(++i>=argc))) break;
911       fHistogramFile=argv[i++];
912       return 2;
913     }
914     // -save-histogram-table
915     if (argument.CompareTo("-save-huffman-table")==0) {
916       if ((bMissingParam=(++i>=argc))) break;
917       fTrainingTableOutput=argv[i++];
918       return 2;
919     }
920   } while (0); // using do-while only to have break available
921
922   if (bMissingParam) {
923     HLTError("missing parameter for argument %s", argument.Data());
924     iResult=-EPROTO;
925   }
926
927   return iResult;
928 }
929
930 int AliHLTTPCDataCompressionComponent::ForwardMCLabels(const AliHLTComponentBlockData& desc,
931                                                        AliHLTSpacePointContainer::AliHLTSpacePointPropertyGrid* pIndex,
932                                                        AliHLTUInt8_t* outputPtr,
933                                                        AliHLTUInt32_t size, AliHLTUInt32_t offset,
934                                                        vector<AliHLTComponentBlockData>& outputBlocks) const
935 {
936   // forward the mc labels in the same order as the clusters
937   // sorted according to the index grid
938   if (!outputPtr || !pIndex) return -EINVAL;
939   if (!desc.fPtr) return -ENODATA;
940   if (size<desc.fSize) return -ENOSPC;
941
942   int slice=AliHLTTPCDefinitions::GetMinSliceNr(desc.fSpecification);
943   int part=AliHLTTPCDefinitions::GetMinPatchNr(desc.fSpecification);
944
945   const AliHLTTPCClusterMCData* pInput = reinterpret_cast<const AliHLTTPCClusterMCData*>(desc.fPtr);
946   Int_t nLabels = (Int_t) pInput->fCount;
947   if (nLabels*sizeof(AliHLTTPCClusterMCLabel) + sizeof(AliHLTTPCClusterMCData) != desc.fSize) {
948     HLTError("inconsistent cluster mc data block size, skipping block");
949     return -EBADF;
950   }
951   const AliHLTTPCClusterMCLabel *pInputLabels = pInput->fLabels;
952
953   AliHLTTPCClusterMCData* pOutput = reinterpret_cast<AliHLTTPCClusterMCData*>(outputPtr);
954   AliHLTTPCClusterMCLabel *pOutputLabels = pOutput->fLabels;
955
956   unsigned outIndex=0;
957   for (AliHLTSpacePointContainer::AliHLTSpacePointPropertyGrid::iterator clusterID=pIndex->begin();
958        clusterID!=pIndex->end(); clusterID++, outIndex++) {
959       if (clusterID.Data().fTrackId>-1) {
960         // this is an assigned cluster, skip
961         // TODO: introduce selectors into AliHLTIndexGrid::begin to loop
962         // consistently over entries, e.g. this check has to be done also
963         // in the forwarding of MC labels in
964         // AliHLTTPCHWCFSpacePointContainer::WriteSorted
965         continue;
966       }
967       if ((unsigned)slice!=AliHLTTPCSpacePointData::GetSlice(clusterID.Data().fId) ||
968           (unsigned)part!=AliHLTTPCSpacePointData::GetPatch(clusterID.Data().fId)) {
969         HLTError("spacepoint index 0x%08x out of slice %d partition %d", clusterID.Data().fId, slice, part);
970       }
971       int index=AliHLTTPCSpacePointData::GetNumber(clusterID.Data().fId);
972       pOutputLabels[outIndex]=pInputLabels[index];
973   }
974   if (outIndex==pInput->fCount) {
975     pOutput->fCount=outIndex;
976   } else {
977     HLTError("failed to copy MC label data block 0x%08x: expecting %d, copied %d entries", desc.fSpecification, pInput->fCount, outIndex);
978     return -EBADMSG;
979   }
980
981   AliHLTComponent_BlockData bd;
982   AliHLTComponent::FillBlockData(bd);
983   bd.fSize          = desc.fSize;
984   bd.fOffset        = offset;
985   bd.fSpecification = desc.fSpecification;
986   bd.fDataType      = desc.fDataType;
987   outputBlocks.push_back(bd);
988
989   return bd.fSize;
990 }