HLT/TPCLib/AliHLTTPCTrackGeometry.cxx

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