3 //**************************************************************************
4 //* This file is property of and copyright by the ALICE HLT Project *
5 //* ALICE Experiment at CERN, All rights reserved. *
7 //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no> *
8 //* for The ALICE HLT Project. *
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 //**************************************************************************
19 /// @file AliHLTTPCTrackGeometry.cxx
20 /// @author Matthias Richter
22 /// @brief Desciption of a track by a sequence of track points
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"
42 /** ROOT macro for the implementation of ROOT specific class methods */
43 ClassImp(AliHLTTPCTrackGeometry)
45 AliHLTTPCTrackGeometry::AliHLTTPCTrackGeometry()
46 : AliHLTTrackGeometry()
48 , fDriftTimeFactorA(0.)
49 , fDriftTimeOffsetA(0.)
50 , fDriftTimeFactorC(0.)
51 , fDriftTimeOffsetC(0.)
53 /// standard constructor
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.)
67 AliHLTTPCTrackGeometry& AliHLTTPCTrackGeometry::operator=(const AliHLTTPCTrackGeometry& src)
69 /// assignment operator
70 AliHLTTrackGeometry::operator=(src);
71 fRawTrackPoints.assign(src.fRawTrackPoints.begin(), src.fRawTrackPoints.end());
75 AliHLTTPCTrackGeometry::~AliHLTTPCTrackGeometry()
80 float AliHLTTPCTrackGeometry::GetPlaneAlpha(AliHLTUInt32_t planeId) const
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();
89 float AliHLTTPCTrackGeometry::GetPlaneR(AliHLTUInt32_t planeId) const
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);
98 float AliHLTTPCTrackGeometry::GetPlaneTheta(AliHLTUInt32_t /*planeId*/) const
100 /// theta of the plane
104 bool AliHLTTPCTrackGeometry::CheckBounds(AliHLTUInt32_t planeId, float u, float /*v*/) const
106 /// check bounds in u and v coordinate
107 float r=GetPlaneR(planeId);
108 if (r<AliHLTTPCTransform::GetFirstRow(0)) return false;
110 // TODO: check if the pad width needs to be considered here
111 return TMath::Abs(TMath::ASin(u/r))<=TMath::Pi()/18;
114 int AliHLTTPCTrackGeometry::CalculateTrackPoints(const AliHLTExternalTrackParam& track)
116 /// calculate the track points, expects the global magnetic field to be initialized
117 AliHLTGlobalBarrelTrack bt(track);
118 return CalculateTrackPoints(bt);
121 int AliHLTTPCTrackGeometry::CalculateTrackPoints(AliHLTGlobalBarrelTrack& track)
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.
130 firstpadrow<AliHLTTPCTransform::GetNRows() &&
131 AliHLTTPCTransform::Row2X(firstpadrow)+AliHLTTPCTransform::GetPadLength(firstpadrow)<track.GetX();
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);
142 int AliHLTTPCTrackGeometry::CalculateTrackPoints(AliHLTGlobalBarrelTrack& track, int firstpadrow, int step)
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);
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
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;
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());
172 HLTError("invalid slice %d calculated from alpha %f", slice, track.GetAlpha());
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());
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;
189 if (TMath::Abs(m)>0.) {
192 // ascending padrows added at end
193 fRawTrackPoints.push_back(AliHLTTrackPoint(id, rpt[1], rpt[2]));
195 // descending padrows added at begin
196 fRawTrackPoints.insert(fRawTrackPoints.begin(), AliHLTTrackPoint(id, rpt[1], rpt[2]));
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]
208 // cout << sout.str();
210 ALIHLTERRORGUARD(1, "drift time correction not initialized, can not add track points in raw coordinates");
217 int AliHLTTPCTrackGeometry::FindMatchingTrackPoint(AliHLTUInt32_t spacepointId, float spacepoint[3], AliHLTUInt32_t& planeId)
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);
229 // find the crossing point of the track with the padrow plane where
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);
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);
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
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);
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);
287 int AliHLTTPCTrackGeometry::RegisterTrackPoints(AliHLTTrackGrid* pGrid) const
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;
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));
303 int AliHLTTPCTrackGeometry::FillTrackPoints(AliHLTTrackGrid* pGrid) const
305 /// fill track points to index grid
306 if (!pGrid) return -EINVAL;
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));
319 AliHLTSpacePointContainer* AliHLTTPCTrackGeometry::ConvertToSpacePoints(bool bAssociated) const
321 /// create a collection of all points
322 std::auto_ptr<AliHLTTPCSpacePointContainer> spacepoints(new AliHLTTPCSpacePointContainer);
323 if (!spacepoints.get()) return NULL;
325 const vector<AliHLTTrackPoint>& trackPoints=TrackPoints();
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;
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
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++;
368 AliHLTComponentBlockData bd;
369 AliHLTComponent::FillBlockData(bd);
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);
377 return spacepoints.release();
380 const AliHLTTrackGeometry::AliHLTTrackPoint* AliHLTTPCTrackGeometry::GetRawTrackPoint(AliHLTUInt32_t id) const
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;
388 AliHLTTrackGeometry::AliHLTTrackPoint* AliHLTTPCTrackGeometry::GetRawTrackPoint(AliHLTUInt32_t id)
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;
396 int AliHLTTPCTrackGeometry::FillRawResidual(int coordinate, TH2* histo, AliHLTSpacePointContainer* points) const
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());
412 if (!points->Check(spacepointId)) {
413 //HLTError("can not find associated space point 0x%08x of track point 0x%08x", spacepointId, trackpoint->GetId());
418 value=rawpoint->GetU()-points->GetY(spacepointId);
419 histo->Fill(GetPlaneR(trackpoint->GetId()), value);
421 value=rawpoint->GetV()-points->GetZ(spacepointId);
422 //histo->Fill(GetPlaneR(trackpoint->GetId()), value);
423 histo->Fill(rawpoint->GetV(), value);
430 int AliHLTTPCTrackGeometry::Write(const AliHLTGlobalBarrelTrack& track,
431 AliHLTSpacePointContainer* pSpacePoints,
432 AliHLTDataDeflater* pDeflater,
433 AliHLTUInt8_t* outputPtr,
435 vector<AliHLTUInt32_t>* writtenClusterIds,
436 const char* option) const
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());
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();
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;
465 int AliHLTTPCTrackGeometry::WriteAssociatedClusters(AliHLTSpacePointContainer* pSpacePoints,
466 AliHLTDataDeflater* pDeflater,
467 vector<AliHLTUInt32_t>* writtenClusterIds,
468 const char* /*option*/) const
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;
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
482 vector<AliHLTTrackPoint>::const_iterator clrow=fRawTrackPoints.end();
483 if (clrow!=fRawTrackPoints.begin()) {
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) {
493 clrow=fRawTrackPoints.begin();
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) {
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);
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);
517 if (writtenClusterIds) {
518 writtenClusterIds->push_back(clid->fId);
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);
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
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
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();
554 HLTError("cluster 0x%08x has negative pad position", clid->fId, clusterpad);
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;
565 deltapad64=clusterpad64-trackpad64;
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);
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;
586 deltatime64=clustertime64-tracktime64;
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++;
606 // set to next trackpoint
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) {
615 // switch direction start from beginning
616 clrow=fRawTrackPoints.begin();
620 // all trackpoints processed
621 clrow=fRawTrackPoints.end();
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;
633 // if (clrow==fRawTrackPoints.end()) {
634 // clrow=fRawTrackPoints.begin();
635 // HLTWarning("no trackpoint on row %d, current point %d", row, thisTrackRow);
639 // no cluster on that padrow
640 AliHLTUInt32_t haveClusters=0;
641 bWriteSuccess=bWriteSuccess && pDeflater->OutputBit(haveClusters);
644 if (!bWriteSuccess) return -ENOSPC;
647 for (clrow=fRawTrackPoints.begin(); clrow!=fRawTrackPoints.end(); clrow++) {
648 allClusters+=clrow->GetSpacepoints().size();
650 if (allClusters!=writtenClusters) {
651 HLTError("track %d mismatch in written clusters: %d but expected %d", GetTrackId(), writtenClusters, allClusters);
654 pDeflater->Pad8Bits();
655 return pDeflater->GetCurrentByteOutputPosition()-dataPosition;
658 int AliHLTTPCTrackGeometry::Read(const AliHLTUInt8_t* buffer,
661 AliHLTUInt32_t& clusterBlockSize,
662 const char* /*option*/)
664 // read track block from buffer
666 if (!buffer) return -EINVAL;
667 if (size<sizeof(AliHLTTPCTrackBlock)) {
668 HLTError("buffer does not contain valid data of track model clusters");
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);
676 if (pTrackBlock->fSize<sizeof(AliHLTTPCTrackBlock)) {
677 HLTError("inconsistent size of track data block specified in the header: %d", pTrackBlock->fSize);
680 AliHLTExternalTrackParam param;
681 memset(¶m, 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);
695 if ((iResult=CalculateTrackPoints(track))<0) {
696 HLTError("failed to calculate track points: %d", iResult);
699 clusterBlockSize=pTrackBlock->fSize-sizeof(AliHLTTPCTrackBlock);
700 return sizeof(AliHLTTPCTrackBlock);