correcting the drift time transformation; optional output of cluster id array in...
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCTrackGeometry.cxx
1 // $Id$
2
3 //**************************************************************************
4 //* This file is property of and copyright by the ALICE HLT Project        * 
5 //* ALICE Experiment at CERN, All rights reserved.                         *
6 //*                                                                        *
7 //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
8 //*                  for The ALICE HLT Project.                            *
9 //*                                                                        *
10 //* Permission to use, copy, modify and distribute this software and its   *
11 //* documentation strictly for non-commercial purposes is hereby granted   *
12 //* without fee, provided that the above copyright notice appears in all   *
13 //* copies and that both the copyright notice and this permission notice   *
14 //* appear in the supporting documentation. The authors make no claims     *
15 //* about the suitability of this software for any purpose. It is          *
16 //* provided "as is" without express or implied warranty.                  *
17 //**************************************************************************
18
19 /// @file   AliHLTTPCTrackGeometry.cxx
20 /// @author Matthias Richter
21 /// @date   2011-05-20
22 /// @brief  Desciption of a track by a sequence of track points
23 ///
24
25 #include "AliHLTTPCTrackGeometry.h"
26 #include "AliHLTTPCTransform.h"
27 #include "AliHLTTPCSpacePointData.h"
28 #include "AliHLTTPCClusterDataFormat.h"
29 #include "AliHLTTPCSpacePointContainer.h"
30 #include "AliHLTTPCHWCFSpacePointContainer.h"
31 #include "AliHLTTPCDefinitions.h"
32 #include "AliHLTComponent.h"
33 #include "AliHLTGlobalBarrelTrack.h"
34 #include "AliHLTDataDeflater.h"
35 #include "AliHLTErrorGuard.h"
36 #include "TMath.h"
37 #include "TH2F.h"
38 #include <memory>
39
40 /** ROOT macro for the implementation of ROOT specific class methods */
41 ClassImp(AliHLTTPCTrackGeometry)
42
43 AliHLTTPCTrackGeometry::AliHLTTPCTrackGeometry()
44   : AliHLTTrackGeometry()
45   , fRawTrackPoints()
46   , fDriftTimeFactorA(0.)
47   , fDriftTimeOffsetA(0.)
48   , fDriftTimeFactorC(0.)
49   , fDriftTimeOffsetC(0.)
50 {
51   /// standard constructor
52 }
53
54 AliHLTTPCTrackGeometry::AliHLTTPCTrackGeometry(const AliHLTTPCTrackGeometry& src)
55   : AliHLTTrackGeometry(src)
56   , fRawTrackPoints(src.fRawTrackPoints)
57   , fDriftTimeFactorA(0.)
58   , fDriftTimeOffsetA(0.)
59   , fDriftTimeFactorC(0.)
60   , fDriftTimeOffsetC(0.)
61 {
62   /// copy constructor
63 }
64
65 AliHLTTPCTrackGeometry& AliHLTTPCTrackGeometry::operator=(const AliHLTTPCTrackGeometry& src)
66 {
67   /// assignment operator
68   AliHLTTrackGeometry::operator=(src);
69   fRawTrackPoints.assign(src.fRawTrackPoints.begin(), src.fRawTrackPoints.end());
70   return *this;
71 }
72
73 AliHLTTPCTrackGeometry::~AliHLTTPCTrackGeometry()
74 {
75   /// destructor
76 }
77
78 float AliHLTTPCTrackGeometry::GetPlaneAlpha(AliHLTUInt32_t planeId) const
79 {
80   /// alpha of the plane
81   UInt_t slice=AliHLTTPCSpacePointData::GetSlice(planeId);
82   float alpha=( slice + 0.5 ) * TMath::Pi() / 9.0;
83   if (alpha>TMath::TwoPi()) alpha-=TMath::TwoPi();
84   return alpha;
85 }
86
87 float AliHLTTPCTrackGeometry::GetPlaneR(AliHLTUInt32_t planeId) const
88 {
89   /// radial distance from global {0,0,0}
90   UInt_t partition=AliHLTTPCSpacePointData::GetPatch(planeId);
91   UInt_t number=AliHLTTPCSpacePointData::GetNumber(planeId);
92   Int_t row=AliHLTTPCTransform::GetFirstRow(partition)+number;
93   return AliHLTTPCTransform::Row2X(row);
94 }
95
96 float AliHLTTPCTrackGeometry::GetPlaneTheta(AliHLTUInt32_t /*planeId*/) const
97 {
98   /// theta of the plane
99   return 0.0;
100 }
101
102 bool AliHLTTPCTrackGeometry::CheckBounds(AliHLTUInt32_t planeId, float u, float /*v*/) const
103 {
104   /// check bounds in u and v coordinate
105   float r=GetPlaneR(planeId);
106   if (r<AliHLTTPCTransform::GetFirstRow(0)) return false;
107
108   // TODO: check if the pad width needs to be considered here
109   return TMath::Abs(TMath::ASin(u/r))<=TMath::Pi()/18;
110 }
111
112 int AliHLTTPCTrackGeometry::CalculateTrackPoints(const AliHLTExternalTrackParam& track)
113 {
114   /// calculate the track points, expects the global magnetic field to be initialized
115   AliHLTGlobalBarrelTrack bt(track);
116   return CalculateTrackPoints(bt);
117 }
118
119 int AliHLTTPCTrackGeometry::CalculateTrackPoints(AliHLTGlobalBarrelTrack& track)
120 {
121   /// calculate the track points, expects the global magnetic field to be initialized
122   int iResult=0;
123   int firstpadrow=0;
124   for (;
125        firstpadrow<AliHLTTPCTransform::GetNRows() && 
126          AliHLTTPCTransform::Row2X(firstpadrow)+AliHLTTPCTransform::GetPadLength(firstpadrow)<track.GetX();
127        firstpadrow++);
128   if (firstpadrow>=AliHLTTPCTransform::GetNRows()) return 0;
129   iResult=CalculateTrackPoints(track, firstpadrow, 1);
130   if (iResult>=0 && firstpadrow>0)
131     iResult=CalculateTrackPoints(track, firstpadrow-1, -1);
132   return iResult;
133 }
134
135 int AliHLTTPCTrackGeometry::CalculateTrackPoints(AliHLTGlobalBarrelTrack& track, int firstpadrow, int step)
136 {
137   /// calculate the track points, expects the global magnetic field to be initialized
138   float offsetAlpha=0.0;
139   for (int padrow=firstpadrow; padrow>=0 && padrow<AliHLTTPCTransform::GetNRows(); padrow+=step) {
140     float x=AliHLTTPCTransform::Row2X(padrow);
141     float y=0.0;
142     float z=0.0;
143
144     int maxshift=9;
145     int shift=0;
146     int result=0;
147     do {
148       // start calculation of crossing points with padrow planes in the slice of the first point
149       // plane alpha corresponds to alpha of the track, switch to neighboring slice if the result
150       // is out of bounds
151       if ((result=track.CalculateCrossingPoint(x, track.GetAlpha()-offsetAlpha, y, z))<1) break;
152       float pointAlpha=TMath::ATan(y/x);
153       if (TMath::Abs(pointAlpha)>TMath::Pi()/18) {
154         offsetAlpha+=(pointAlpha>0?-1:1)*TMath::Pi()/9;
155         result=0;
156       }
157     } while (result==0 && shift++<maxshift);
158     if (result<1) continue;
159     float planealpha=track.GetAlpha()-offsetAlpha;
160     if (planealpha<0) planealpha+=TMath::TwoPi();
161     int slice=int(9*planealpha/TMath::Pi());
162     if (z<0) slice+=18;
163     int partition=AliHLTTPCTransform::GetPatch(padrow);
164     int row=padrow-AliHLTTPCTransform::GetFirstRow(partition);
165     UInt_t id=AliHLTTPCSpacePointData::GetID(slice, partition, row);
166     if (TMath::Abs(planealpha-GetPlaneAlpha(id))>0.0001) {
167       HLTError("alpha missmatch for plane %08x (slice %d): alpha from id %f (%.0f), expected %f (%.0f)", id, slice, GetPlaneAlpha(id), 180*GetPlaneAlpha(id)/TMath::Pi(), planealpha, 180*planealpha/TMath::Pi());
168     }
169     if (AddTrackPoint(AliHLTTrackPoint(id, y, z), AliHLTTPCSpacePointData::GetID(slice, partition, 0))>=0) {
170       Float_t rpt[3]={0.,y,z}; // row pad time
171       AliHLTTPCTransform::LocHLT2Raw(rpt, slice, padrow);
172       float m=fDriftTimeFactorA;
173       float n=fDriftTimeOffsetA;
174       if (slice>=18) {
175         m=fDriftTimeFactorC;
176         n=fDriftTimeOffsetC;
177       }
178       if (TMath::Abs(m)>0.) {
179         rpt[2]=(z-n)/m;
180         fRawTrackPoints.push_back(AliHLTTrackPoint(id, rpt[1], rpt[2]));
181       } else {
182         ALIHLTERRORGUARD(1, "drift time correction not initialized, can not add track points in raw coordinates");
183       }
184     }
185   }
186   return 0;
187 }
188
189 int AliHLTTPCTrackGeometry::FindMatchingTrackPoint(AliHLTUInt32_t spacepointId, float spacepoint[3], AliHLTUInt32_t& planeId)
190 {
191   /// find the track point which can be associated to a spacepoint with coordinates and id
192   UInt_t slice=AliHLTTPCSpacePointData::GetSlice(spacepointId);
193   UInt_t partition=AliHLTTPCSpacePointData::GetPatch(spacepointId);
194   int row=AliHLTTPCTransform::GetPadRow(spacepoint[0]);
195   bool bSpecialRow=row==30 || row==90 || row==139;
196   if (row<AliHLTTPCTransform::GetFirstRow(partition) || row>AliHLTTPCTransform::GetLastRow(partition)) {
197     HLTError("row number %d calculated from x value %f is outside slice %d partition %d", row, spacepoint[0], slice, partition);
198     return -EINVAL;
199   }
200
201   // find the crossing point of the track with the padrow plane where
202   // the spacepoint is
203   // 1) calculate plane id from slice, partition and row (within partition)
204   row-=AliHLTTPCTransform::GetFirstRow(partition);
205   UInt_t id=AliHLTTPCSpacePointData::GetID(slice, partition, row);
206   const AliHLTTrackPoint* point=GetTrackPoint(id);
207   // track might be outside the partition and cross the central membrane
208   // search in the other half of the TPC
209   if (!point && slice<18) {
210     // search in the neighboring partition on the C side
211     id=AliHLTTPCSpacePointData::GetID(slice+18, partition, row);
212     point=GetTrackPoint(id);
213   } else if (!point && slice>=18) {
214     // search in the neighboring partition on the A side
215     id=AliHLTTPCSpacePointData::GetID(slice-18, partition, row);
216     point=GetTrackPoint(id);
217   }
218   
219   // search in the neighboring partition, this takes account for rows
220   // 30, 90, and 139 which are partly in one and the other partition
221   if (!point && bSpecialRow) {
222     row+=AliHLTTPCTransform::GetFirstRow(partition);
223     row-=AliHLTTPCTransform::GetFirstRow(partition-1);
224     id=AliHLTTPCSpacePointData::GetID(slice, partition-1, row);
225     point=GetTrackPoint(id);
226     if (!point && slice<18) {
227       // search in the neighboring partition on the C side
228       id=AliHLTTPCSpacePointData::GetID(slice+18, partition-1, row);
229       point=GetTrackPoint(id);
230     } else if (!point && slice>=18) {
231       // search in the neighboring partition on the A side
232       id=AliHLTTPCSpacePointData::GetID(slice-18, partition-1, row);
233       point=GetTrackPoint(id);
234     }
235   }
236
237   if (point) {
238     planeId=id;
239     if (point->HaveAssociatedSpacePoint()) {
240       if (GetVerbosity()>2) HLTInfo("descarding spacepoint 0x%08x z=%f y=%f z=%f: track point 0x%08x already occupied", spacepoint[0], spacepoint[1], spacepoint[2], planeId);
241       return 0; // already occupied
242     }
243     float maxdy=2.;
244     float maxdz=2.;
245     if (TMath::Abs(point->GetU()-spacepoint[1])>maxdy) {
246       if (GetVerbosity()>0) HLTInfo("descarding spacepoint 0x%08x y=%f z=%f: track point 0x%08x y %f outside tolerance %f", spacepoint[1], spacepoint[2], planeId, point->GetU(), maxdy);
247       return -ENOENT;
248     }
249     if (TMath::Abs(point->GetV()-spacepoint[2])>maxdz) {
250       if (GetVerbosity()>0) HLTInfo("descarding spacepoint 0x%08x y=%f z=%f: track point 0x%08x z %f outside tolerance %f", spacepoint[1], spacepoint[2], planeId, point->GetV(), maxdz);
251       return -ENOENT;
252     }
253     return 1;
254   }
255   return -ENOENT;
256 }
257
258
259 int AliHLTTPCTrackGeometry::RegisterTrackPoints(AliHLTTrackGrid* pGrid) const
260 {
261   /// register track points in the index grid, at this step the number
262   /// of tracks in each cell is counted
263   if (!pGrid) return -EINVAL;
264   int iResult=0;
265   for (vector<AliHLTTrackPoint>::const_iterator tp=TrackPoints().begin();
266        tp!=TrackPoints().end() && iResult>=0; tp++) {
267     AliHLTUInt32_t id=tp->GetId();
268     iResult=pGrid->CountSpacePoint(AliHLTTPCSpacePointData::GetSlice(id),
269                                    AliHLTTPCSpacePointData::GetPatch(id),
270                                    AliHLTTPCSpacePointData::GetNumber(id));
271   }
272   return iResult;
273 }
274
275 int AliHLTTPCTrackGeometry::FillTrackPoints(AliHLTTrackGrid* pGrid) const
276 {
277   /// fill track points to index grid
278   if (!pGrid) return -EINVAL;
279   int iResult=0;
280   for (vector<AliHLTTrackPoint>::const_iterator tp=TrackPoints().begin();
281        tp!=TrackPoints().end() && iResult>=0; tp++) {
282     AliHLTUInt32_t id=tp->GetId();
283     iResult=pGrid->AddSpacePoint(GetTrackId(),
284                                  AliHLTTPCSpacePointData::GetSlice(id),
285                                  AliHLTTPCSpacePointData::GetPatch(id),
286                                  AliHLTTPCSpacePointData::GetNumber(id));
287   }
288   return iResult;
289 }
290
291 AliHLTSpacePointContainer* AliHLTTPCTrackGeometry::ConvertToSpacePoints(bool bAssociated) const
292 {
293   /// create a collection of all points
294   std::auto_ptr<AliHLTTPCSpacePointContainer> spacepoints(new AliHLTTPCSpacePointContainer);
295   if (!spacepoints.get()) return NULL;
296
297   const vector<AliHLTTrackPoint>& trackPoints=TrackPoints();
298   unsigned i=0;
299   while (i<trackPoints.size()) {
300     // allocate buffer for all points, even though the buffer might not be filled
301     // completely because of a partition change
302     int nofPoints=trackPoints.size()-i;
303     int blocksize=sizeof(AliHLTTPCClusterData)+nofPoints*sizeof(AliHLTTPCSpacePointData);
304     AliHLTUInt8_t* pBuffer=spacepoints->Alloc(blocksize);
305     if (!pBuffer) return NULL;
306     AliHLTTPCClusterData* pClusterData=reinterpret_cast<AliHLTTPCClusterData*>(pBuffer);
307     pClusterData->fSpacePointCnt=0;
308     AliHLTTPCSpacePointData* pClusters=pClusterData->fSpacePoints;
309     int currentSlice=-1;
310     int currentPartition=-1;
311     for (; i<trackPoints.size(); i++) {
312       if (bAssociated && !trackPoints[i].HaveAssociatedSpacePoint()) continue;
313       AliHLTUInt32_t planeId=trackPoints[i].GetId();
314       int slice=AliHLTTPCSpacePointData::GetSlice(planeId);
315       int partition=AliHLTTPCSpacePointData::GetPatch(planeId);
316       int number=AliHLTTPCSpacePointData::GetNumber(planeId);
317       if ((currentSlice>=0 && currentSlice!=slice) || (currentPartition>=0 && currentPartition!=partition)) {
318         // change of partition or slice, need to go to next block
319         // 2011-07-26 currently all spacepoints go into one block, if separated
320         // blocks per partition are needed one has to leave the inner loop here
321         // and set the data block specification below
322         // Caution: not tested, only the last block seems to make it through
323         //break;
324       }
325       currentSlice=slice;
326       currentPartition=partition;
327       pClusters[pClusterData->fSpacePointCnt].fX=GetPlaneR(planeId);
328       pClusters[pClusterData->fSpacePointCnt].fY=trackPoints[i].GetU();
329       pClusters[pClusterData->fSpacePointCnt].fZ=trackPoints[i].GetV();
330       pClusters[pClusterData->fSpacePointCnt].fID=planeId;
331       pClusters[pClusterData->fSpacePointCnt].fPadRow=AliHLTTPCTransform::GetFirstRow(partition)+number;
332       pClusters[pClusterData->fSpacePointCnt].fSigmaY2=0.;
333       pClusters[pClusterData->fSpacePointCnt].fSigmaZ2=0.;
334       pClusters[pClusterData->fSpacePointCnt].fCharge=0;
335       pClusters[pClusterData->fSpacePointCnt].fQMax=0;
336       pClusters[pClusterData->fSpacePointCnt].fUsed=0;
337       pClusters[pClusterData->fSpacePointCnt].fTrackN=0;
338       pClusterData->fSpacePointCnt++;
339     }
340     AliHLTComponentBlockData bd;
341     AliHLTComponent::FillBlockData(bd);
342     bd.fPtr=pBuffer;
343     bd.fSize=sizeof(AliHLTTPCClusterData)+pClusterData->fSpacePointCnt*sizeof(AliHLTTPCSpacePointData);
344     AliHLTComponent::SetDataType(bd.fDataType, "CLUSTERS", "TPC ");
345     bd.fSpecification=kAliHLTVoidDataSpec;//AliHLTTPCDefinitions::EncodeDataSpecification(currentSlice, currentSlice, currentPartition, currentPartition);
346     spacepoints->AddInputBlock(&bd);
347   }
348
349   return spacepoints.release();
350 }
351
352 const AliHLTTrackGeometry::AliHLTTrackPoint* AliHLTTPCTrackGeometry::GetRawTrackPoint(AliHLTUInt32_t id) const
353 {
354   /// get raw track point of id
355   const AliHLTTrackGeometry::AliHLTTrackPoint* p=find(&fRawTrackPoints[0], &fRawTrackPoints[fRawTrackPoints.size()], id);
356   if (p==&fRawTrackPoints[fRawTrackPoints.size()]) return 0;
357   return p;
358 }
359
360 AliHLTTrackGeometry::AliHLTTrackPoint* AliHLTTPCTrackGeometry::GetRawTrackPoint(AliHLTUInt32_t id)
361 {
362   /// get raw track point of id
363   AliHLTTrackGeometry::AliHLTTrackPoint* p=find(&fRawTrackPoints[0], &fRawTrackPoints[fRawTrackPoints.size()], id);
364   if (p==&fRawTrackPoints[fRawTrackPoints.size()]) return 0;
365   return p;
366 }
367
368 int AliHLTTPCTrackGeometry::FillRawResidual(int coordinate, TH2* histo, AliHLTSpacePointContainer* points) const
369 {
370   // fill residual histogram
371   if (!histo || !points) return -EINVAL;
372   const vector<AliHLTTrackPoint>& trackPoints=TrackPoints();
373   for (vector<AliHLTTrackPoint>::const_iterator trackpoint=trackPoints.begin();
374        trackpoint!=trackPoints.end(); trackpoint++) {
375     if (!trackpoint->HaveAssociatedSpacePoint()) continue;
376     for (vector<AliHLTTrackSpacepoint>::const_iterator sp=(trackpoint->GetSpacepoints()).begin();
377          sp!=(trackpoint->GetSpacepoints()).end(); sp++) {
378     AliHLTUInt32_t spacepointId=sp->fId;
379     vector<AliHLTTrackPoint>::const_iterator rawpoint=find(fRawTrackPoints.begin(), fRawTrackPoints.end(), trackpoint->GetId());
380     if (rawpoint==fRawTrackPoints.end()) {
381       HLTError("can not find track raw coordinates of track point 0x%08x", trackpoint->GetId());
382       continue;
383     }
384     if (!points->Check(spacepointId)) {
385       //HLTError("can not find associated space point 0x%08x of track point 0x%08x", spacepointId, trackpoint->GetId());
386       continue;
387     }
388     float value=0.;
389     if (coordinate==0) {
390       value=rawpoint->GetU()-points->GetY(spacepointId);
391       histo->Fill(GetPlaneR(trackpoint->GetId()), value);
392     } else {
393       value=rawpoint->GetV()-points->GetZ(spacepointId);
394       //histo->Fill(GetPlaneR(trackpoint->GetId()), value);
395       histo->Fill(rawpoint->GetV(), value);
396     }
397     }
398   }
399   return 0;
400 }
401
402 int AliHLTTPCTrackGeometry::Write(const AliHLTGlobalBarrelTrack& track,
403                                   AliHLTSpacePointContainer* pSpacePoints,
404                                   AliHLTDataDeflater* pDeflater,
405                                   AliHLTUInt8_t* outputPtr,
406                                   AliHLTUInt32_t size,
407                                   vector<AliHLTUInt32_t>* writtenClusterIds,
408                                   const char* option) const
409 {
410   // write track block to buffer
411   if (size<=sizeof(AliHLTTPCTrackBlock)) return -ENOSPC;
412   AliHLTTPCTrackBlock* pTrackBlock=reinterpret_cast<AliHLTTPCTrackBlock*>(outputPtr);
413   pTrackBlock->fSize=sizeof(AliHLTTPCTrackBlock); // size of cluster block added later
414   pTrackBlock->fSlice=AliHLTUInt8_t(9*track.GetAlpha()/TMath::Pi());
415   pTrackBlock->fReserved=0;
416   pTrackBlock->fX      = track.GetX();
417   pTrackBlock->fY      = track.GetY();
418   pTrackBlock->fZ      = track.GetZ();
419   pTrackBlock->fSinPsi = track.GetSnp();
420   pTrackBlock->fTgl    = track.GetTgl();
421   pTrackBlock->fq1Pt   = track.GetSigned1Pt();
422
423   pDeflater->Clear();
424   pDeflater->InitBitDataOutput(reinterpret_cast<AliHLTUInt8_t*>(outputPtr+sizeof(AliHLTTPCTrackBlock)), size-sizeof(AliHLTTPCTrackBlock));
425   int result=WriteAssociatedClusters(pSpacePoints, pDeflater, writtenClusterIds, option);
426   if (result<0) return result;
427   pTrackBlock->fSize+=result;
428   return pTrackBlock->fSize;
429 }
430
431 int AliHLTTPCTrackGeometry::WriteAssociatedClusters(AliHLTSpacePointContainer* pSpacePoints,
432                                                     AliHLTDataDeflater* pDeflater,
433                                                     vector<AliHLTUInt32_t>* writtenClusterIds,
434                                                     const char* /*option*/) const
435 {
436   // write associated clusters to buffer via deflater
437   if (!pDeflater || !pSpacePoints) return -EINVAL;
438   AliHLTTPCHWCFSpacePointContainer* pTPCRawSpacePoints=dynamic_cast<AliHLTTPCHWCFSpacePointContainer*>(pSpacePoints);
439   if (!pTPCRawSpacePoints) return -EINVAL;
440   bool bReverse=true;
441   bool bWriteSuccess=true;
442   int writtenClusters=0;
443   // filling of track points starts from first point on track outwards, and
444   // then from that point inwards. That's why the lower padrows might be in
445   // reverse order at the end of the track point array. If the last element
446   // is bigger than the first element, only trackpoints in ascending order
447   // are in the array
448   vector<AliHLTTrackPoint>::const_iterator clrow=fRawTrackPoints.end();
449   if (clrow!=fRawTrackPoints.begin()) {
450     clrow--;
451     AliHLTUInt32_t partition=AliHLTTPCSpacePointData::GetPatch(clrow->GetId());
452     AliHLTUInt32_t partitionrow=AliHLTTPCSpacePointData::GetNumber(clrow->GetId());
453     partitionrow+=AliHLTTPCTransform::GetFirstRow(partition);
454     AliHLTUInt32_t firstpartition=AliHLTTPCSpacePointData::GetPatch(fRawTrackPoints.begin()->GetId());
455     AliHLTUInt32_t firstpartitionrow=AliHLTTPCSpacePointData::GetNumber(fRawTrackPoints.begin()->GetId());
456     firstpartitionrow+=AliHLTTPCTransform::GetFirstRow(firstpartition);
457     if (partitionrow>=firstpartitionrow) {
458       bReverse=false;
459       clrow=fRawTrackPoints.begin();
460     }
461   }
462   unsigned long dataPosition=pDeflater->GetCurrentByteOutputPosition();
463   for (unsigned row=0; row<159 && bWriteSuccess; row++) {
464     if (clrow!=fRawTrackPoints.end()) {
465       AliHLTUInt32_t thisPartition=AliHLTTPCSpacePointData::GetPatch(clrow->GetId());
466       AliHLTUInt32_t thisTrackRow=AliHLTTPCSpacePointData::GetNumber(clrow->GetId());
467       thisTrackRow+=AliHLTTPCTransform::GetFirstRow(thisPartition);
468       if (thisTrackRow==row) {
469         // write clusters
470         const vector<AliHLTTrackSpacepoint>&  clusters=clrow->GetSpacepoints();
471         AliHLTUInt32_t haveClusters=clusters.size()>0;
472         // 1 bit for clusters on that padrow
473         bWriteSuccess=bWriteSuccess && pDeflater->OutputBit(haveClusters);
474         if (haveClusters) {
475           bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kClusterCount, clusters.size());
476           for (vector<AliHLTTrackSpacepoint>::const_iterator clid=clusters.begin();
477                clid!=clusters.end() && bWriteSuccess; clid++) {
478             if (!pSpacePoints->Check(clid->fId)) {
479               HLTError("can not find spacepoint 0x%08x", clid->fId);
480               continue;
481             }
482             if (writtenClusterIds) {
483               writtenClusterIds->push_back(clid->fId);
484             }
485
486             float deltapad =clid->fdU;
487             float deltatime =clid->fdV;
488             float sigmaY2=pSpacePoints->GetYWidth(clid->fId);
489             float sigmaZ2=pSpacePoints->GetZWidth(clid->fId);
490             AliHLTUInt64_t charge=(AliHLTUInt64_t)pSpacePoints->GetCharge(clid->fId);
491             AliHLTUInt64_t qmax=(AliHLTUInt64_t)pTPCRawSpacePoints->GetQMax(clid->fId);
492
493             AliHLTUInt64_t deltapad64=0;
494             AliHLTUInt32_t signDeltaPad=0;
495             if (!isnan(deltapad)) {
496               if (deltapad<0.) {deltapad*=-1; signDeltaPad=1;}
497               deltapad*=AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kResidualPad].fScale;
498               deltapad64=(AliHLTUInt64_t)round(deltapad);
499             }
500             AliHLTUInt64_t deltatime64=0;
501             AliHLTUInt32_t signDeltaTime=0;
502             if (!isnan(deltatime)) {
503               if (deltatime<0.) {deltatime*=-1; signDeltaTime=1;}
504               deltatime*=AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kResidualTime].fScale;
505               deltatime64=(AliHLTUInt64_t)round(deltatime);
506             }
507             AliHLTUInt64_t sigmaY264=0;
508             if (!isnan(sigmaY2)) sigmaY264=(AliHLTUInt64_t)round(sigmaY2*AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kSigmaY2].fScale);
509             AliHLTUInt64_t sigmaZ264=0;
510             if (!isnan(sigmaZ2)) sigmaZ264=(AliHLTUInt64_t)round(sigmaZ2*AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kSigmaZ2].fScale);
511             bWriteSuccess=bWriteSuccess && pDeflater->OutputBit(signDeltaPad);
512             bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kResidualPad    , deltapad64);  
513             bWriteSuccess=bWriteSuccess && pDeflater->OutputBit(signDeltaTime);
514             bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kResidualTime   , deltatime64);
515             bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kSigmaY2        , sigmaY264);
516             bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kSigmaZ2        , sigmaZ264);
517             bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kCharge         , charge);
518             bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kQMax           , qmax);
519             if (bWriteSuccess) writtenClusters++;
520           }
521         }
522
523         // set to next trackpoint
524         if (bReverse) {
525           if (clrow!=fRawTrackPoints.begin()) {
526             AliHLTUInt32_t nextPartition=AliHLTTPCSpacePointData::GetPatch((clrow-1)->GetId());
527             AliHLTUInt32_t nextTrackRow=AliHLTTPCSpacePointData::GetNumber((clrow-1)->GetId());
528             nextTrackRow+=AliHLTTPCTransform::GetFirstRow(nextPartition);
529             if (thisTrackRow+1==nextTrackRow) {
530               clrow--;
531             } else {
532               // switch direction start from beginning
533               clrow=fRawTrackPoints.begin();
534               bReverse=false;
535             }
536           } else {
537             // all trackpoints processed
538             clrow=fRawTrackPoints.end();
539           }
540         } else {
541           clrow++;
542         }
543         continue;
544       } else {
545         // sequence not ordered, search
546         // this has been fixed and the search is no longer necessary
547         // for (clrow=fRawTrackPoints.begin(); clrow!=fRawTrackPoints.end(); clrow++) {
548         //   if ((AliHLTTPCSpacePointData::GetNumber(clrow->GetId())+AliHLTTPCTransform::GetFirstRow(AliHLTTPCSpacePointData::GetPatch(clrow->GetId())))==row) break;
549         // }
550         // if (clrow==fRawTrackPoints.end()) {
551         //   clrow=fRawTrackPoints.begin();
552         //   HLTWarning("no trackpoint on row %d, current point %d", row, thisTrackRow);
553         // }
554       }
555     }
556     // no cluster on that padrow
557     AliHLTUInt32_t haveClusters=0;
558     bWriteSuccess=bWriteSuccess && pDeflater->OutputBit(haveClusters);
559   }
560
561   if (!bWriteSuccess) return -ENOSPC;
562
563   int allClusters=0;
564   for (clrow=fRawTrackPoints.begin(); clrow!=fRawTrackPoints.end(); clrow++) {
565     allClusters+=clrow->GetSpacepoints().size();
566   }
567   if (allClusters!=writtenClusters) {
568     HLTError("track %d mismatch in written clusters: %d but expected %d", GetTrackId(), writtenClusters, allClusters);
569   }
570
571   pDeflater->Pad8Bits();
572   return pDeflater->GetCurrentByteOutputPosition()-dataPosition;
573 }
574
575 // int AliHLTTPCTrackGeometry::Read(AliHLTUInt8_t* buffer,
576 //                               AliHLTUInt32_t size,
577 //                               const char* /*option*/) const
578 // {
579 //   // read track block from buffer
580 //   if (!buffer) return -EINVAL;
581 //   if (size<sizeof(AliHLTTPCTrackBlock)) {
582 //     HLTError("buffer does not contain valid data of track model clusters");
583 //     return -ENODATA;
584 //   }
585 //   AliHLTTPCTrackBlock* pTrackBlock=reinterpret_cast<AliHLTTPCTrackBlock*>(buffer);
586 //   if (pTrackBlock->fSize>size) {
587 //     HLTError("inconsistent track data block of size %d exceeds available buffer of size %d", pTrackBlock->fSize, size);
588 //     return -ENODATA;
589 //   }
590 //   pTrackBlock->fSlice=AliHLTUInt8_t(9*track.GetAlpha()/TMath::Pi());
591 //   pTrackBlock->fReserved=0;
592 //   pTrackBlock->fX      = track.GetX();
593 //   pTrackBlock->fY      = track.GetY();
594 //   pTrackBlock->fZ      = track.GetZ();
595 //   pTrackBlock->fSinPsi = track.GetSnp();
596 //   pTrackBlock->fTgl    = track.GetTgl();
597 //   pTrackBlock->fq1Pt   = track.GetSigned1Pt();
598
599 // }