Update master to aliroot
[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 #include <sstream>
40 using namespace std;
41
42 /** ROOT macro for the implementation of ROOT specific class methods */
43 ClassImp(AliHLTTPCTrackGeometry)
44
45 AliHLTTPCTrackGeometry::AliHLTTPCTrackGeometry()
46   : AliHLTTrackGeometry()
47   , fRawTrackPoints()
48   , fDriftTimeFactorA(0.)
49   , fDriftTimeOffsetA(0.)
50   , fDriftTimeFactorC(0.)
51   , fDriftTimeOffsetC(0.)
52 {
53   /// standard constructor
54 }
55
56 AliHLTTPCTrackGeometry::AliHLTTPCTrackGeometry(const AliHLTTPCTrackGeometry& src)
57   : AliHLTTrackGeometry(src)
58   , fRawTrackPoints(src.fRawTrackPoints)
59   , fDriftTimeFactorA(0.)
60   , fDriftTimeOffsetA(0.)
61   , fDriftTimeFactorC(0.)
62   , fDriftTimeOffsetC(0.)
63 {
64   /// copy constructor
65 }
66
67 AliHLTTPCTrackGeometry& AliHLTTPCTrackGeometry::operator=(const AliHLTTPCTrackGeometry& src)
68 {
69   /// assignment operator
70   AliHLTTrackGeometry::operator=(src);
71   fRawTrackPoints.assign(src.fRawTrackPoints.begin(), src.fRawTrackPoints.end());
72   return *this;
73 }
74
75 AliHLTTPCTrackGeometry::~AliHLTTPCTrackGeometry()
76 {
77   /// destructor
78 }
79
80 float AliHLTTPCTrackGeometry::GetPlaneAlpha(AliHLTUInt32_t planeId) const
81 {
82   /// alpha of the plane
83   UInt_t slice=AliHLTTPCSpacePointData::GetSlice(planeId);
84   float alpha=( slice + 0.5 ) * TMath::Pi() / 9.0;
85   if (alpha>TMath::TwoPi()) alpha-=TMath::TwoPi();
86   return alpha;
87 }
88
89 float AliHLTTPCTrackGeometry::GetPlaneR(AliHLTUInt32_t planeId) const
90 {
91   /// radial distance from global {0,0,0}
92   UInt_t partition=AliHLTTPCSpacePointData::GetPatch(planeId);
93   UInt_t number=AliHLTTPCSpacePointData::GetNumber(planeId);
94   Int_t row=AliHLTTPCTransform::GetFirstRow(partition)+number;
95   return AliHLTTPCTransform::Row2X(row);
96 }
97
98 float AliHLTTPCTrackGeometry::GetPlaneTheta(AliHLTUInt32_t /*planeId*/) const
99 {
100   /// theta of the plane
101   return 0.0;
102 }
103
104 bool AliHLTTPCTrackGeometry::CheckBounds(AliHLTUInt32_t planeId, float u, float /*v*/) const
105 {
106   /// check bounds in u and v coordinate
107   float r=GetPlaneR(planeId);
108   if (r<AliHLTTPCTransform::GetFirstRow(0)) return false;
109
110   // TODO: check if the pad width needs to be considered here
111   return TMath::Abs(TMath::ASin(u/r))<=TMath::Pi()/18;
112 }
113
114 int AliHLTTPCTrackGeometry::CalculateTrackPoints(const AliHLTExternalTrackParam& track)
115 {
116   /// calculate the track points, expects the global magnetic field to be initialized
117   AliHLTGlobalBarrelTrack bt(track);
118   return CalculateTrackPoints(bt);
119 }
120
121 int AliHLTTPCTrackGeometry::CalculateTrackPoints(AliHLTGlobalBarrelTrack& track)
122 {
123   /// calculate the track points, expects the global magnetic field to be initialized
124   /// depending on the x coordinate of the first track point the corresponding padrow
125   /// is searched and the points are calculated outwards. Eventually the points are
126   /// calculated inwards as a second step.
127   int iResult=0;
128   int firstpadrow=0;
129   for (;
130        firstpadrow<AliHLTTPCTransform::GetNRows() && 
131          AliHLTTPCTransform::Row2X(firstpadrow)+AliHLTTPCTransform::GetPadLength(firstpadrow)<track.GetX();
132        firstpadrow++);
133   if (firstpadrow>=AliHLTTPCTransform::GetNRows()) return 0;
134   // first calculated outwards
135   iResult=CalculateTrackPoints(track, firstpadrow, 1);
136   if (iResult>=0 && firstpadrow>0)
137     // now calculate inwards
138     iResult=CalculateTrackPoints(track, firstpadrow-1, -1);
139   return iResult;
140 }
141
142 int AliHLTTPCTrackGeometry::CalculateTrackPoints(AliHLTGlobalBarrelTrack& track, int firstpadrow, int step)
143 {
144   /// calculate the track points, expects the global magnetic field to be initialized
145   float offsetAlpha=0.0;
146   for (int padrow=firstpadrow; padrow>=0 && padrow<AliHLTTPCTransform::GetNRows(); padrow+=step) {
147     float x=AliHLTTPCTransform::Row2X(padrow);
148     float y=0.0;
149     float z=0.0;
150
151     int maxshift=9;
152     int shift=0;
153     int result=0;
154     do {
155       // start calculation of crossing points with padrow planes in the slice of the first point
156       // plane alpha corresponds to alpha of the track, switch to neighboring slice if the result
157       // is out of bounds
158       if ((result=track.CalculateCrossingPoint(x, track.GetAlpha()-offsetAlpha, y, z))<1) break;
159       float pointAlpha=TMath::ATan(y/x);
160       if (TMath::Abs(pointAlpha)>TMath::Pi()/18) {
161         offsetAlpha+=(pointAlpha>0?-1:1)*TMath::Pi()/9;
162         result=0;
163       }
164     } while (result==0 && shift++<maxshift);
165     if (result<1) continue;
166     float planealpha=track.GetAlpha()-offsetAlpha;
167     if (planealpha<0) planealpha+=TMath::TwoPi();
168     if (planealpha>TMath::TwoPi()) planealpha-=TMath::TwoPi();
169     int slice=int(9*planealpha/TMath::Pi());
170     if (z<0) slice+=18;
171     if (slice>=36) {
172       HLTError("invalid slice %d calculated from alpha %f", slice, track.GetAlpha());
173     }
174     int partition=AliHLTTPCTransform::GetPatch(padrow);
175     int row=padrow-AliHLTTPCTransform::GetFirstRow(partition);
176     UInt_t id=AliHLTTPCSpacePointData::GetID(slice, partition, row);
177     if (TMath::Abs(planealpha-GetPlaneAlpha(id))>0.0001) {
178       HLTError("alpha missmatch for plane %08x (slice %d): alpha from id %f (%.0f deg), expected %f (%.0f deg)", id, slice, GetPlaneAlpha(id), 180*GetPlaneAlpha(id)/TMath::Pi(), planealpha, 180*planealpha/TMath::Pi());
179     }
180     if (AddTrackPoint(AliHLTTrackPoint(id, y, z), AliHLTTPCSpacePointData::GetID(slice, partition, 0))>=0) {
181       Float_t rpt[3]={0.,y,z}; // row pad time
182       AliHLTTPCTransform::LocHLT2Raw(rpt, slice, padrow);
183       float m=fDriftTimeFactorA;
184       float n=fDriftTimeOffsetA;
185       if (slice>=18) {
186         m=fDriftTimeFactorC;
187         n=fDriftTimeOffsetC;
188       }
189       if (TMath::Abs(m)>0.) {
190         rpt[2]=(z-n)/m;
191         if (step>0) {
192           // ascending padrows added at end
193           fRawTrackPoints.push_back(AliHLTTrackPoint(id, rpt[1], rpt[2]));
194         } else {
195           // descending padrows added at begin
196           fRawTrackPoints.insert(fRawTrackPoints.begin(), AliHLTTrackPoint(id, rpt[1], rpt[2]));
197         }
198         // FIXME: implement a Print function for the raw track points
199         // stringstream sout;
200         // sout << "  slice "  << setfill(' ') << setw(3) << fixed << right << slice
201         //      << "  row "    << setfill(' ') << setw(3) << fixed << right << padrow
202         //      << "  id 0x"     << hex << setw(8) << id << dec
203         //      << "  y "      << setfill(' ') << setw(5) << fixed << right << setprecision (2) << y
204         //      << "  z "      << setfill(' ') << setw(5) << fixed << right << setprecision (2) << z
205         //      << "  pad "    << setfill(' ') << setw(5) << fixed << right << setprecision (2) << rpt[1]
206         //      << "  time "   << setfill(' ') << setw(5) << fixed << right << setprecision (2) << rpt[2]
207         //      << endl;
208         // cout << sout.str();
209       } else {
210         ALIHLTERRORGUARD(1, "drift time correction not initialized, can not add track points in raw coordinates");
211       }
212     }
213   }
214   return 0;
215 }
216
217 int AliHLTTPCTrackGeometry::FindMatchingTrackPoint(AliHLTUInt32_t spacepointId, float spacepoint[3], AliHLTUInt32_t& planeId)
218 {
219   /// find the track point which can be associated to a spacepoint with coordinates and id
220   UInt_t slice=AliHLTTPCSpacePointData::GetSlice(spacepointId);
221   UInt_t partition=AliHLTTPCSpacePointData::GetPatch(spacepointId);
222   int row=AliHLTTPCTransform::GetPadRow(spacepoint[0]);
223   bool bSpecialRow=row==30 || row==90 || row==139;
224   if (row<AliHLTTPCTransform::GetFirstRow(partition) || row>AliHLTTPCTransform::GetLastRow(partition)) {
225     HLTError("row number %d calculated from x value %f is outside slice %d partition %d", row, spacepoint[0], slice, partition);
226     return -EINVAL;
227   }
228
229   // find the crossing point of the track with the padrow plane where
230   // the spacepoint is
231   // 1) calculate plane id from slice, partition and row (within partition)
232   row-=AliHLTTPCTransform::GetFirstRow(partition);
233   UInt_t id=AliHLTTPCSpacePointData::GetID(slice, partition, row);
234   const AliHLTTrackPoint* point=GetTrackPoint(id);
235   // track might be outside the partition and cross the central membrane
236   // search in the other half of the TPC
237   if (!point && slice<18) {
238     // search in the neighboring partition on the C side
239     id=AliHLTTPCSpacePointData::GetID(slice+18, partition, row);
240     point=GetTrackPoint(id);
241   } else if (!point && slice>=18) {
242     // search in the neighboring partition on the A side
243     id=AliHLTTPCSpacePointData::GetID(slice-18, partition, row);
244     point=GetTrackPoint(id);
245   }
246   
247   // search in the neighboring partition, this takes account for rows
248   // 30, 90, and 139 which are partly in one and the other partition
249   if (!point && bSpecialRow) {
250     row+=AliHLTTPCTransform::GetFirstRow(partition);
251     row-=AliHLTTPCTransform::GetFirstRow(partition-1);
252     id=AliHLTTPCSpacePointData::GetID(slice, partition-1, row);
253     point=GetTrackPoint(id);
254     if (!point && slice<18) {
255       // search in the neighboring partition on the C side
256       id=AliHLTTPCSpacePointData::GetID(slice+18, partition-1, row);
257       point=GetTrackPoint(id);
258     } else if (!point && slice>=18) {
259       // search in the neighboring partition on the A side
260       id=AliHLTTPCSpacePointData::GetID(slice-18, partition-1, row);
261       point=GetTrackPoint(id);
262     }
263   }
264
265   if (point) {
266     planeId=id;
267     if (point->HaveAssociatedSpacePoint()) {
268       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);
269       return 0; // already occupied
270     }
271     float maxdy=2.;
272     float maxdz=2.;
273     if (TMath::Abs(point->GetU()-spacepoint[1])>maxdy) {
274       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);
275       return -ENOENT;
276     }
277     if (TMath::Abs(point->GetV()-spacepoint[2])>maxdz) {
278       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);
279       return -ENOENT;
280     }
281     return 1;
282   }
283   return -ENOENT;
284 }
285
286
287 int AliHLTTPCTrackGeometry::RegisterTrackPoints(AliHLTTrackGrid* pGrid) const
288 {
289   /// register track points in the index grid, at this step the number
290   /// of tracks in each cell is counted
291   if (!pGrid) return -EINVAL;
292   int iResult=0;
293   for (vector<AliHLTTrackPoint>::const_iterator tp=TrackPoints().begin();
294        tp!=TrackPoints().end() && iResult>=0; tp++) {
295     AliHLTUInt32_t id=tp->GetId();
296     iResult=pGrid->CountSpacePoint(AliHLTTPCSpacePointData::GetSlice(id),
297                                    AliHLTTPCSpacePointData::GetPatch(id),
298                                    AliHLTTPCSpacePointData::GetNumber(id));
299   }
300   return iResult;
301 }
302
303 int AliHLTTPCTrackGeometry::FillTrackPoints(AliHLTTrackGrid* pGrid) const
304 {
305   /// fill track points to index grid
306   if (!pGrid) return -EINVAL;
307   int iResult=0;
308   for (vector<AliHLTTrackPoint>::const_iterator tp=TrackPoints().begin();
309        tp!=TrackPoints().end() && iResult>=0; tp++) {
310     AliHLTUInt32_t id=tp->GetId();
311     iResult=pGrid->AddSpacePoint(GetTrackId(),
312                                  AliHLTTPCSpacePointData::GetSlice(id),
313                                  AliHLTTPCSpacePointData::GetPatch(id),
314                                  AliHLTTPCSpacePointData::GetNumber(id));
315   }
316   return iResult;
317 }
318
319 AliHLTSpacePointContainer* AliHLTTPCTrackGeometry::ConvertToSpacePoints(bool bAssociated) const
320 {
321   /// create a collection of all points
322   std::auto_ptr<AliHLTTPCSpacePointContainer> spacepoints(new AliHLTTPCSpacePointContainer);
323   if (!spacepoints.get()) return NULL;
324
325   const vector<AliHLTTrackPoint>& trackPoints=TrackPoints();
326   unsigned i=0;
327   while (i<trackPoints.size()) {
328     // allocate buffer for all points, even though the buffer might not be filled
329     // completely because of a partition change
330     int nofPoints=trackPoints.size()-i;
331     int blocksize=sizeof(AliHLTTPCClusterData)+nofPoints*sizeof(AliHLTTPCSpacePointData);
332     AliHLTUInt8_t* pBuffer=spacepoints->Alloc(blocksize);
333     if (!pBuffer) return NULL;
334     AliHLTTPCClusterData* pClusterData=reinterpret_cast<AliHLTTPCClusterData*>(pBuffer);
335     pClusterData->fSpacePointCnt=0;
336     AliHLTTPCSpacePointData* pClusters=pClusterData->fSpacePoints;
337     int currentSlice=-1;
338     int currentPartition=-1;
339     for (; i<trackPoints.size(); i++) {
340       if (bAssociated && !trackPoints[i].HaveAssociatedSpacePoint()) continue;
341       AliHLTUInt32_t planeId=trackPoints[i].GetId();
342       int slice=AliHLTTPCSpacePointData::GetSlice(planeId);
343       int partition=AliHLTTPCSpacePointData::GetPatch(planeId);
344       int number=AliHLTTPCSpacePointData::GetNumber(planeId);
345       if ((currentSlice>=0 && currentSlice!=slice) || (currentPartition>=0 && currentPartition!=partition)) {
346         // change of partition or slice, need to go to next block
347         // 2011-07-26 currently all spacepoints go into one block, if separated
348         // blocks per partition are needed one has to leave the inner loop here
349         // and set the data block specification below
350         // Caution: not tested, only the last block seems to make it through
351         //break;
352       }
353       currentSlice=slice;
354       currentPartition=partition;
355       pClusters[pClusterData->fSpacePointCnt].fX=GetPlaneR(planeId);
356       pClusters[pClusterData->fSpacePointCnt].fY=trackPoints[i].GetU();
357       pClusters[pClusterData->fSpacePointCnt].fZ=trackPoints[i].GetV();
358       pClusters[pClusterData->fSpacePointCnt].fID=planeId;
359       pClusters[pClusterData->fSpacePointCnt].fPadRow=AliHLTTPCTransform::GetFirstRow(partition)+number;
360       pClusters[pClusterData->fSpacePointCnt].fSigmaY2=0.;
361       pClusters[pClusterData->fSpacePointCnt].fSigmaZ2=0.;
362       pClusters[pClusterData->fSpacePointCnt].fCharge=0;
363       pClusters[pClusterData->fSpacePointCnt].fQMax=0;
364       pClusters[pClusterData->fSpacePointCnt].fUsed=0;
365       pClusters[pClusterData->fSpacePointCnt].fTrackN=0;
366       pClusterData->fSpacePointCnt++;
367     }
368     AliHLTComponentBlockData bd;
369     AliHLTComponent::FillBlockData(bd);
370     bd.fPtr=pBuffer;
371     bd.fSize=sizeof(AliHLTTPCClusterData)+pClusterData->fSpacePointCnt*sizeof(AliHLTTPCSpacePointData);
372     AliHLTComponent::SetDataType(bd.fDataType, "CLUSTERS", "TPC ");
373     bd.fSpecification=kAliHLTVoidDataSpec;//AliHLTTPCDefinitions::EncodeDataSpecification(currentSlice, currentSlice, currentPartition, currentPartition);
374     spacepoints->AddInputBlock(&bd);
375   }
376
377   return spacepoints.release();
378 }
379
380 const AliHLTTrackGeometry::AliHLTTrackPoint* AliHLTTPCTrackGeometry::GetRawTrackPoint(AliHLTUInt32_t id) const
381 {
382   /// get raw track point of id
383   const AliHLTTrackGeometry::AliHLTTrackPoint* p=find(&fRawTrackPoints[0], &fRawTrackPoints[fRawTrackPoints.size()], id);
384   if (p==&fRawTrackPoints[fRawTrackPoints.size()]) return 0;
385   return p;
386 }
387
388 AliHLTTrackGeometry::AliHLTTrackPoint* AliHLTTPCTrackGeometry::GetRawTrackPoint(AliHLTUInt32_t id)
389 {
390   /// get raw track point of id
391   AliHLTTrackGeometry::AliHLTTrackPoint* p=find(&fRawTrackPoints[0], &fRawTrackPoints[fRawTrackPoints.size()], id);
392   if (p==&fRawTrackPoints[fRawTrackPoints.size()]) return 0;
393   return p;
394 }
395
396 int AliHLTTPCTrackGeometry::FillRawResidual(int coordinate, TH2* histo, AliHLTSpacePointContainer* points) const
397 {
398   // fill residual histogram
399   if (!histo || !points) return -EINVAL;
400   const vector<AliHLTTrackPoint>& trackPoints=TrackPoints();
401   for (vector<AliHLTTrackPoint>::const_iterator trackpoint=trackPoints.begin();
402        trackpoint!=trackPoints.end(); trackpoint++) {
403     if (!trackpoint->HaveAssociatedSpacePoint()) continue;
404     for (vector<AliHLTTrackSpacepoint>::const_iterator sp=(trackpoint->GetSpacepoints()).begin();
405          sp!=(trackpoint->GetSpacepoints()).end(); sp++) {
406     AliHLTUInt32_t spacepointId=sp->fId;
407     vector<AliHLTTrackPoint>::const_iterator rawpoint=find(fRawTrackPoints.begin(), fRawTrackPoints.end(), trackpoint->GetId());
408     if (rawpoint==fRawTrackPoints.end()) {
409       HLTError("can not find track raw coordinates of track point 0x%08x", trackpoint->GetId());
410       continue;
411     }
412     if (!points->Check(spacepointId)) {
413       //HLTError("can not find associated space point 0x%08x of track point 0x%08x", spacepointId, trackpoint->GetId());
414       continue;
415     }
416     float value=0.;
417     if (coordinate==0) {
418       value=rawpoint->GetU()-points->GetY(spacepointId);
419       histo->Fill(GetPlaneR(trackpoint->GetId()), value);
420     } else {
421       value=rawpoint->GetV()-points->GetZ(spacepointId);
422       //histo->Fill(GetPlaneR(trackpoint->GetId()), value);
423       histo->Fill(rawpoint->GetV(), value);
424     }
425     }
426   }
427   return 0;
428 }
429
430 int AliHLTTPCTrackGeometry::Write(const AliHLTGlobalBarrelTrack& track,
431                                   AliHLTSpacePointContainer* pSpacePoints,
432                                   AliHLTDataDeflater* pDeflater,
433                                   AliHLTUInt8_t* outputPtr,
434                                   AliHLTUInt32_t size,
435                                   vector<AliHLTUInt32_t>* writtenClusterIds,
436                                   const char* option) const
437 {
438   // write track block to buffer
439   if (size<=sizeof(AliHLTTPCTrackBlock)) return -ENOSPC;
440   AliHLTTPCTrackBlock* pTrackBlock=reinterpret_cast<AliHLTTPCTrackBlock*>(outputPtr);
441   pTrackBlock->fSize=sizeof(AliHLTTPCTrackBlock); // size of cluster block added later
442   float alpha=track.GetAlpha();
443   while (alpha<0.) alpha+=TMath::TwoPi();
444   while (alpha>TMath::TwoPi()) alpha-=TMath::TwoPi();
445   pTrackBlock->fSlice=AliHLTUInt8_t(9*alpha/TMath::Pi());
446   if (pTrackBlock->fSlice>=36) {
447     HLTError("invalid slice %d calculated from alpha %f", pTrackBlock->fSlice, track.GetAlpha());
448   }
449   pTrackBlock->fReserved=0;
450   pTrackBlock->fX      = track.GetX();
451   pTrackBlock->fY      = track.GetY();
452   pTrackBlock->fZ      = track.GetZ();
453   pTrackBlock->fSinPsi = track.GetSnp();
454   pTrackBlock->fTgl    = track.GetTgl();
455   pTrackBlock->fq1Pt   = track.GetSigned1Pt();
456
457   pDeflater->Clear();
458   pDeflater->InitBitDataOutput(reinterpret_cast<AliHLTUInt8_t*>(outputPtr+sizeof(AliHLTTPCTrackBlock)), size-sizeof(AliHLTTPCTrackBlock));
459   int result=WriteAssociatedClusters(pSpacePoints, pDeflater, writtenClusterIds, option);
460   if (result<0) return result;
461   pTrackBlock->fSize+=result;
462   return pTrackBlock->fSize;
463 }
464
465 int AliHLTTPCTrackGeometry::WriteAssociatedClusters(AliHLTSpacePointContainer* pSpacePoints,
466                                                     AliHLTDataDeflater* pDeflater,
467                                                     vector<AliHLTUInt32_t>* writtenClusterIds,
468                                                     const char* /*option*/) const
469 {
470   // write associated clusters to buffer via deflater
471   if (!pDeflater || !pSpacePoints) return -EINVAL;
472   AliHLTTPCHWCFSpacePointContainer* pTPCRawSpacePoints=dynamic_cast<AliHLTTPCHWCFSpacePointContainer*>(pSpacePoints);
473   if (!pTPCRawSpacePoints) return -EINVAL;
474   bool bReverse=true;
475   bool bWriteSuccess=true;
476   int writtenClusters=0;
477   // filling of track points starts from first point on track outwards, and
478   // then from that point inwards. That's why the lower padrows might be in
479   // reverse order at the end of the track point array. If the last element
480   // is bigger than the first element, only trackpoints in ascending order
481   // are in the array
482   vector<AliHLTTrackPoint>::const_iterator clrow=fRawTrackPoints.end();
483   if (clrow!=fRawTrackPoints.begin()) {
484     clrow--;
485     AliHLTUInt32_t partition=AliHLTTPCSpacePointData::GetPatch(clrow->GetId());
486     AliHLTUInt32_t partitionrow=AliHLTTPCSpacePointData::GetNumber(clrow->GetId());
487     partitionrow+=AliHLTTPCTransform::GetFirstRow(partition);
488     AliHLTUInt32_t firstpartition=AliHLTTPCSpacePointData::GetPatch(fRawTrackPoints.begin()->GetId());
489     AliHLTUInt32_t firstpartitionrow=AliHLTTPCSpacePointData::GetNumber(fRawTrackPoints.begin()->GetId());
490     firstpartitionrow+=AliHLTTPCTransform::GetFirstRow(firstpartition);
491     if (partitionrow>=firstpartitionrow) {
492       bReverse=false;
493       clrow=fRawTrackPoints.begin();
494     }
495   }
496   const AliHLTUInt32_t clusterCountBitLength=AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kClusterCount].fBitLength;
497   unsigned long dataPosition=pDeflater->GetCurrentByteOutputPosition();
498   for (unsigned row=0; row<159 && bWriteSuccess; row++) {
499     if (clrow!=fRawTrackPoints.end()) {
500       AliHLTUInt32_t thisPartition=AliHLTTPCSpacePointData::GetPatch(clrow->GetId());
501       AliHLTUInt32_t thisTrackRow=AliHLTTPCSpacePointData::GetNumber(clrow->GetId());
502       thisTrackRow+=AliHLTTPCTransform::GetFirstRow(thisPartition);
503       if (thisTrackRow==row) {
504         // write clusters
505         const vector<AliHLTTrackSpacepoint>&  clusters=clrow->GetSpacepoints();
506         AliHLTUInt32_t haveClusters=clusters.size()>0;
507         // 1 bit for clusters on that padrow
508         bWriteSuccess=bWriteSuccess && pDeflater->OutputBit(haveClusters);
509         if (haveClusters) {
510           bWriteSuccess=bWriteSuccess && pDeflater->OutputBits(clusters.size(), clusterCountBitLength);
511           for (vector<AliHLTTrackSpacepoint>::const_iterator clid=clusters.begin();
512                clid!=clusters.end() && bWriteSuccess; clid++) {
513             if (!pSpacePoints->Check(clid->fId)) {
514               HLTError("can not find spacepoint 0x%08x", clid->fId);
515               continue;
516             }
517             if (writtenClusterIds) {
518               writtenClusterIds->push_back(clid->fId);
519             }
520
521             // FIXME: there is a bug in the calculation of the residuals stored with the
522             // assiciated space point, calculate again, but needs to be fixed
523             float deltapad =pSpacePoints->GetY(clid->fId)-clrow->GetU();//clid->fdU;
524             float deltatime =pSpacePoints->GetZ(clid->fId)-clrow->GetV();//clid->fdV;
525             if (TMath::Abs(deltapad)>=AliHLTTPCDefinitions::fgkMaxClusterDeltaPad ||
526                 TMath::Abs(deltatime)>=AliHLTTPCDefinitions::fgkMaxClusterDeltaTime) {
527               AliHLTUInt8_t slice = AliHLTTPCSpacePointData::GetSlice(clid->fId);
528               AliHLTUInt8_t partition = AliHLTTPCSpacePointData::GetPatch(clid->fId);
529               HLTFatal("cluster 0x%08x slice %d partition %d: residual out of range - pad %f (max %d), time %f (max %d)", clid->fId, slice, partition, deltapad, AliHLTTPCDefinitions::fgkMaxClusterDeltaPad, deltatime, AliHLTTPCDefinitions::fgkMaxClusterDeltaTime);
530             }
531             float sigmaY2=pSpacePoints->GetYWidth(clid->fId);
532             float sigmaZ2=pSpacePoints->GetZWidth(clid->fId);
533             AliHLTUInt64_t charge=(AliHLTUInt64_t)pSpacePoints->GetCharge(clid->fId);
534             AliHLTUInt64_t qmax=(AliHLTUInt64_t)pTPCRawSpacePoints->GetQMax(clid->fId);
535             // cout << "  row "    << setfill(' ') << setw(3) << fixed << right << row
536             //   << "  pad "    << setfill(' ') << setw(7) << fixed << right << setprecision (4) << pSpacePoints->GetY(clid->fId)
537             //   << "  dpad "   << setfill(' ') << setw(7) << fixed << right << setprecision (4) << deltapad
538             //   << "  time "   << setfill(' ') << setw(7) << fixed << right << setprecision (4) << pSpacePoints->GetZ(clid->fId)
539             //   << "  dtime "  << setfill(' ') << setw(7) << fixed << right << setprecision (4) << deltatime
540             //   << "  charge " << setfill(' ') << setw(5) << fixed << right << charge
541             //   << "  qmax "   << setfill(' ') << setw(4) << fixed << right << qmax
542             //   << endl;
543
544             // time and pad coordinates are scaled and transformed to integer values for
545             // both cluster and track point before calculating the residual. this makes
546             // the compression lossless with respect to the format without track model
547             // compression
548             AliHLTUInt64_t deltapad64=0;
549             AliHLTUInt32_t signDeltaPad=0;
550             if (!isnan(deltapad)) {
551               double clusterpad=pSpacePoints->GetY(clid->fId);
552               double trackpad=clrow->GetU();
553               if (clusterpad<0.) {
554                 HLTError("cluster 0x%08x has negative pad position", clid->fId, clusterpad);
555               }
556               clusterpad*=AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kResidualPad].fScale;
557               trackpad*=AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kResidualPad].fScale;
558               AliHLTUInt64_t clusterpad64=(AliHLTUInt64_t)round(clusterpad);
559               AliHLTUInt64_t trackpad64=0;
560               if (trackpad>0.) trackpad64=(AliHLTUInt64_t)round(trackpad);
561               if (clusterpad64<trackpad64) {
562                 deltapad64=trackpad64-clusterpad64;
563                 signDeltaPad=1;
564               } else {
565                 deltapad64=clusterpad64-trackpad64;
566                 signDeltaPad=0;
567               }
568             }
569             AliHLTUInt64_t deltatime64=0;
570             AliHLTUInt32_t signDeltaTime=0;
571             if (!isnan(deltatime)) {
572               double clustertime=pSpacePoints->GetZ(clid->fId);
573               double tracktime=clrow->GetV();
574               if (clustertime<0.) {
575                 HLTError("cluster 0x%08x has negative time position", clid->fId, clustertime);
576               }
577               clustertime*=AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kResidualTime].fScale;
578               tracktime*=AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kResidualTime].fScale;
579               AliHLTUInt64_t clustertime64=(AliHLTUInt64_t)round(clustertime);
580               AliHLTUInt64_t tracktime64=0;
581               if (tracktime>0.) tracktime64=(AliHLTUInt64_t)round(tracktime);
582               if (clustertime64<tracktime64) {
583                 deltatime64=tracktime64-clustertime64;
584                 signDeltaTime=1;
585               } else {
586                 deltatime64=clustertime64-tracktime64;
587                 signDeltaTime=0;
588               }
589             }
590             AliHLTUInt64_t sigmaY264=0;
591             if (!isnan(sigmaY2)) sigmaY264=(AliHLTUInt64_t)round(sigmaY2*AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kSigmaY2].fScale);
592             AliHLTUInt64_t sigmaZ264=0;
593             if (!isnan(sigmaZ2)) sigmaZ264=(AliHLTUInt64_t)round(sigmaZ2*AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kSigmaZ2].fScale);
594             bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kResidualPad    , deltapad64);  
595             bWriteSuccess=bWriteSuccess && pDeflater->OutputBit(signDeltaPad);
596             bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kResidualTime   , deltatime64);
597             bWriteSuccess=bWriteSuccess && pDeflater->OutputBit(signDeltaTime);
598             bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kSigmaY2        , sigmaY264);
599             bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kSigmaZ2        , sigmaZ264);
600             bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kCharge         , charge);
601             bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kQMax           , qmax);
602             if (bWriteSuccess) writtenClusters++;
603           }
604         }
605
606         // set to next trackpoint
607         if (bReverse) {
608           if (clrow!=fRawTrackPoints.begin()) {
609             AliHLTUInt32_t nextPartition=AliHLTTPCSpacePointData::GetPatch((clrow-1)->GetId());
610             AliHLTUInt32_t nextTrackRow=AliHLTTPCSpacePointData::GetNumber((clrow-1)->GetId());
611             nextTrackRow+=AliHLTTPCTransform::GetFirstRow(nextPartition);
612             if (thisTrackRow+1==nextTrackRow) {
613               clrow--;
614             } else {
615               // switch direction start from beginning
616               clrow=fRawTrackPoints.begin();
617               bReverse=false;
618             }
619           } else {
620             // all trackpoints processed
621             clrow=fRawTrackPoints.end();
622           }
623         } else {
624           clrow++;
625         }
626         continue;
627       } else {
628         // sequence not ordered, search
629         // this has been fixed and the search is no longer necessary
630         // for (clrow=fRawTrackPoints.begin(); clrow!=fRawTrackPoints.end(); clrow++) {
631         //   if ((AliHLTTPCSpacePointData::GetNumber(clrow->GetId())+AliHLTTPCTransform::GetFirstRow(AliHLTTPCSpacePointData::GetPatch(clrow->GetId())))==row) break;
632         // }
633         // if (clrow==fRawTrackPoints.end()) {
634         //   clrow=fRawTrackPoints.begin();
635         //   HLTWarning("no trackpoint on row %d, current point %d", row, thisTrackRow);
636         // }
637       }
638     }
639     // no cluster on that padrow
640     AliHLTUInt32_t haveClusters=0;
641     bWriteSuccess=bWriteSuccess && pDeflater->OutputBit(haveClusters);
642   }
643
644   if (!bWriteSuccess) return -ENOSPC;
645
646   int allClusters=0;
647   for (clrow=fRawTrackPoints.begin(); clrow!=fRawTrackPoints.end(); clrow++) {
648     allClusters+=clrow->GetSpacepoints().size();
649   }
650   if (allClusters!=writtenClusters) {
651     HLTError("track %d mismatch in written clusters: %d but expected %d", GetTrackId(), writtenClusters, allClusters);
652   }
653
654   pDeflater->Pad8Bits();
655   return pDeflater->GetCurrentByteOutputPosition()-dataPosition;
656 }
657
658 int AliHLTTPCTrackGeometry::Read(const AliHLTUInt8_t* buffer,
659                                  AliHLTUInt32_t size,
660                                  float bz,
661                                  AliHLTUInt32_t& clusterBlockSize,
662                                  const char* /*option*/)
663 {
664   // read track block from buffer
665   int iResult=0;
666   if (!buffer) return -EINVAL;
667   if (size<sizeof(AliHLTTPCTrackBlock)) {
668     HLTError("buffer does not contain valid data of track model clusters");
669     return -ENODATA;
670   }
671   const AliHLTTPCTrackBlock* pTrackBlock=reinterpret_cast<const AliHLTTPCTrackBlock*>(buffer);
672   if (pTrackBlock->fSize>size) {
673     HLTError("inconsistent track data block of size %d exceeds available buffer of size %d", pTrackBlock->fSize, size);
674     return -ENODATA;
675   }
676   if (pTrackBlock->fSize<sizeof(AliHLTTPCTrackBlock)) {
677     HLTError("inconsistent size of track data block specified in the header: %d", pTrackBlock->fSize);
678     return -ENODATA;
679   }
680   AliHLTExternalTrackParam param;
681   memset(&param, 0, sizeof(param));
682   param.fAlpha   =( pTrackBlock->fSlice + 0.5 ) * TMath::Pi() / 9.0;
683   if (param.fAlpha>TMath::TwoPi()) param.fAlpha-=TMath::TwoPi();
684   param.fX       = pTrackBlock->fX;
685   param.fY       = pTrackBlock->fY;
686   param.fZ       = pTrackBlock->fZ;
687   param.fSinPsi  = pTrackBlock->fSinPsi;
688   param.fTgl     = pTrackBlock->fTgl;
689   param.fq1Pt    = pTrackBlock->fq1Pt;
690   AliHLTGlobalBarrelTrack track(param);
691   if ((iResult=track.CalculateHelixParams(bz))<0) {
692     HLTError("failed to calculate helix params: %d", iResult);
693     return iResult;
694   }
695   if ((iResult=CalculateTrackPoints(track))<0) {
696     HLTError("failed to calculate track points: %d", iResult);
697     return iResult;
698   }
699   clusterBlockSize=pTrackBlock->fSize-sizeof(AliHLTTPCTrackBlock);
700   return sizeof(AliHLTTPCTrackBlock);
701 }