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