HLT/TPCLib/AliHLTTPCTrackGeometry.cxx

   1 // $Id$
   2
   3 //**************************************************************************
   4 //* This file is property of and copyright by the ALICE HLT Project        *
   5 //* ALICE Experiment at CERN, All rights reserved.                         *
   6 //*                                                                        *
   7 //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
   8 //*                  for The ALICE HLT Project.                            *
   9 //*                                                                        *
  10 //* Permission to use, copy, modify and distribute this software and its   *
  11 //* documentation strictly for non-commercial purposes is hereby granted   *
  12 //* without fee, provided that the above copyright notice appears in all   *
  13 //* copies and that both the copyright notice and this permission notice   *
  14 //* appear in the supporting documentation. The authors make no claims     *
  15 //* about the suitability of this software for any purpose. It is          *
  16 //* provided "as is" without express or implied warranty.                  *
  17 //**************************************************************************
  18
  19 /// @file   AliHLTTPCTrackGeometry.cxx
  20 /// @author Matthias Richter
  21 /// @date   2011-05-20
  22 /// @brief  Desciption of a track by a sequence of track points
  23 ///
  24
  25 #include "AliHLTTPCTrackGeometry.h"
  26 #include "AliHLTTPCTransform.h"
  27 #include "AliHLTTPCSpacePointData.h"
  28 #include "AliHLTTPCClusterDataFormat.h"
  29 #include "AliHLTTPCSpacePointContainer.h"
  30 #include "AliHLTTPCHWCFSpacePointContainer.h"
  31 #include "AliHLTTPCDefinitions.h"
  32 #include "AliHLTComponent.h"
  33 #include "AliHLTGlobalBarrelTrack.h"
  34 #include "AliHLTDataDeflater.h"
  35 #include "AliHLTErrorGuard.h"
  36 #include "TMath.h"
  37 #include "TH2F.h"
  38 #include <memory>
  39 #include <sstream>
  40 using namespace std;
  41
  42 /** ROOT macro for the implementation of ROOT specific class methods */
  43 ClassImp(AliHLTTPCTrackGeometry)
  44
  45 AliHLTTPCTrackGeometry::AliHLTTPCTrackGeometry()
  46   : AliHLTTrackGeometry()
  47   , fRawTrackPoints()
  48   , fDriftTimeFactorA(0.)
  49   , fDriftTimeOffsetA(0.)
  50   , fDriftTimeFactorC(0.)
  51   , fDriftTimeOffsetC(0.)
  52 {
  53   /// standard constructor
  54 }
  55
  56 AliHLTTPCTrackGeometry::AliHLTTPCTrackGeometry(const AliHLTTPCTrackGeometry& src)
  57   : AliHLTTrackGeometry(src)
  58   , fRawTrackPoints(src.fRawTrackPoints)
  59   , fDriftTimeFactorA(0.)
  60   , fDriftTimeOffsetA(0.)
  61   , fDriftTimeFactorC(0.)
  62   , fDriftTimeOffsetC(0.)
  63 {
  64   /// copy constructor
  65 }
  66
  67 AliHLTTPCTrackGeometry& AliHLTTPCTrackGeometry::operator=(const AliHLTTPCTrackGeometry& src)
  68 {
  69   /// assignment operator
  70   AliHLTTrackGeometry::operator=(src);
  71   fRawTrackPoints.assign(src.fRawTrackPoints.begin(), src.fRawTrackPoints.end());
  72   return *this;
  73 }
  74
  75 AliHLTTPCTrackGeometry::~AliHLTTPCTrackGeometry()
  76 {
  77   /// destructor
  78 }
  79
  80 float AliHLTTPCTrackGeometry::GetPlaneAlpha(AliHLTUInt32_t planeId) const
  81 {
  82   /// alpha of the plane
  83   UInt_t slice=AliHLTTPCSpacePointData::GetSlice(planeId);
  84   float alpha=( slice + 0.5 ) * TMath::Pi() / 9.0;
  85   if (alpha>TMath::TwoPi()) alpha-=TMath::TwoPi();
  86   return alpha;
  87 }
  88
  89 float AliHLTTPCTrackGeometry::GetPlaneR(AliHLTUInt32_t planeId) const
  90 {
  91   /// radial distance from global {0,0,0}
  92   UInt_t partition=AliHLTTPCSpacePointData::GetPatch(planeId);
  93   UInt_t number=AliHLTTPCSpacePointData::GetNumber(planeId);
  94   Int_t row=AliHLTTPCTransform::GetFirstRow(partition)+number;
  95   return AliHLTTPCTransform::Row2X(row);
  96 }
  97
  98 float AliHLTTPCTrackGeometry::GetPlaneTheta(AliHLTUInt32_t /*planeId*/) const
  99 {
 100   /// theta of the plane
 101   return 0.0;
 102 }
 103
 104 bool AliHLTTPCTrackGeometry::CheckBounds(AliHLTUInt32_t planeId, float u, float /*v*/) const
 105 {
 106   /// check bounds in u and v coordinate
 107   float r=GetPlaneR(planeId);
 108   if (r<AliHLTTPCTransform::GetFirstRow(0)) return false;
 109
 110   // TODO: check if the pad width needs to be considered here
 111   return TMath::Abs(TMath::ASin(u/r))<=TMath::Pi()/18;
 112 }
 113
 114 int AliHLTTPCTrackGeometry::CalculateTrackPoints(const AliHLTExternalTrackParam& track)
 115 {
 116   /// calculate the track points, expects the global magnetic field to be initialized
 117   AliHLTGlobalBarrelTrack bt(track);
 118   return CalculateTrackPoints(bt);
 119 }
 120
 121 int AliHLTTPCTrackGeometry::CalculateTrackPoints(AliHLTGlobalBarrelTrack& track)
 122 {
 123   /// calculate the track points, expects the global magnetic field to be initialized
 124   int iResult=0;
 125   int firstpadrow=0;
 126   for (;
 127        firstpadrow<AliHLTTPCTransform::GetNRows() &&
 128          AliHLTTPCTransform::Row2X(firstpadrow)+AliHLTTPCTransform::GetPadLength(firstpadrow)<track.GetX();
 129        firstpadrow++);
 130   if (firstpadrow>=AliHLTTPCTransform::GetNRows()) return 0;
 131   iResult=CalculateTrackPoints(track, firstpadrow, 1);
 132   if (iResult>=0 && firstpadrow>0)
 133     iResult=CalculateTrackPoints(track, firstpadrow-1, -1);
 134   return iResult;
 135 }
 136
 137 int 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();
 163     if (planealpha>TMath::TwoPi()) planealpha-=TMath::TwoPi();
 164     int slice=int(9*planealpha/TMath::Pi());
 165     if (z<0) slice+=18;
 166     if (slice>=36) {
 167       HLTError("invalid slice %d calculated from alpha %f", slice, track.GetAlpha());
 168     }
 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) {
 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());
 174     }
 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);
 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;
 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
 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;
 203         // cout << sout.str();
 204       } else {
 205         ALIHLTERRORGUARD(1, "drift time correction not initialized, can not add track points in raw coordinates");
 206       }
 207     }
 208   }
 209   return 0;
 210 }
 211
 212 int 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]);
 218   bool bSpecialRow=row==30 || row==90 || row==139;
 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   }
 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)
 227   row-=AliHLTTPCTransform::GetFirstRow(partition);
 228   UInt_t id=AliHLTTPCSpacePointData::GetID(slice, partition, row);
 229   const AliHLTTrackPoint* point=GetTrackPoint(id);
 230   // track might be outside the partition and cross the central membrane
 231   // search in the other half of the TPC
 232   if (!point && slice<18) {
 233     // search in the neighboring partition on the C side
 234     id=AliHLTTPCSpacePointData::GetID(slice+18, partition, row);
 235     point=GetTrackPoint(id);
 236   } else if (!point && slice>=18) {
 237     // search in the neighboring partition on the A side
 238     id=AliHLTTPCSpacePointData::GetID(slice-18, partition, row);
 239     point=GetTrackPoint(id);
 240   }
 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
 260   if (point) {
 261     planeId=id;
 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     }
 276     return 1;
 277   }
 278   return -ENOENT;
 279 }
 280
 281
 282 int 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
 298 int 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
 314 AliHLTSpacePointContainer* AliHLTTPCTrackGeometry::ConvertToSpacePoints(bool bAssociated) const
 315 {
 316   /// create a collection of all points
 317   std::auto_ptr<AliHLTTPCSpacePointContainer> spacepoints(new AliHLTTPCSpacePointContainer);
 318   if (!spacepoints.get()) return NULL;
 319
 320   const vector<AliHLTTrackPoint>& trackPoints=TrackPoints();
 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++) {
 335       if (bAssociated && !trackPoints[i].HaveAssociatedSpacePoint()) continue;
 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;
 366     bd.fSize=sizeof(AliHLTTPCClusterData)+pClusterData->fSpacePointCnt*sizeof(AliHLTTPCSpacePointData);
 367     AliHLTComponent::SetDataType(bd.fDataType, "CLUSTERS", "TPC ");
 368     bd.fSpecification=kAliHLTVoidDataSpec;//AliHLTTPCDefinitions::EncodeDataSpecification(currentSlice, currentSlice, currentPartition, currentPartition);
 369     spacepoints->AddInputBlock(&bd);
 370   }
 371
 372   return spacepoints.release();
 373 }
 374
 375 const 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
 383 AliHLTTrackGeometry::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
 391 int 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++) {
 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;
 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     }
 420     }
 421   }
 422   return 0;
 423 }
 424
 425 int AliHLTTPCTrackGeometry::Write(const AliHLTGlobalBarrelTrack& track,
 426                                   AliHLTSpacePointContainer* pSpacePoints,
 427                                   AliHLTDataDeflater* pDeflater,
 428                                   AliHLTUInt8_t* outputPtr,
 429                                   AliHLTUInt32_t size,
 430                                   vector<AliHLTUInt32_t>* writtenClusterIds,
 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
 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   }
 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));
 454   int result=WriteAssociatedClusters(pSpacePoints, pDeflater, writtenClusterIds, option);
 455   if (result<0) return result;
 456   pTrackBlock->fSize+=result;
 457   return pTrackBlock->fSize;
 458 }
 459
 460 int AliHLTTPCTrackGeometry::WriteAssociatedClusters(AliHLTSpacePointContainer* pSpacePoints,
 461                                                     AliHLTDataDeflater* pDeflater,
 462                                                     vector<AliHLTUInt32_t>* writtenClusterIds,
 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   }
 491   const AliHLTUInt32_t clusterCountBitLength=AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kClusterCount].fBitLength;
 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) {
 505           bWriteSuccess=bWriteSuccess && pDeflater->OutputBits(clusters.size(), clusterCountBitLength);
 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             }
 512             if (writtenClusterIds) {
 513               writtenClusterIds->push_back(clid->fId);
 514             }
 515
 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;
 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             }
 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);
 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;
 538
 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
 543             AliHLTUInt64_t deltapad64=0;
 544             AliHLTUInt32_t signDeltaPad=0;
 545             if (!isnan(deltapad)) {
 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               }
 563             }
 564             AliHLTUInt64_t deltatime64=0;
 565             AliHLTUInt32_t signDeltaTime=0;
 566             if (!isnan(deltatime)) {
 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               }
 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);
 589             bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kResidualPad    , deltapad64);
 590             bWriteSuccess=bWriteSuccess && pDeflater->OutputBit(signDeltaPad);
 591             bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kResidualTime   , deltatime64);
 592             bWriteSuccess=bWriteSuccess && pDeflater->OutputBit(signDeltaTime);
 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
 653 int 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 }