HLT/TPCLib/AliHLTTPCTrackGeometry.cxx

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