HLT/TPCLib/comp/AliHLTTPCCompModelConverter.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: Timm Steinbeck <timm@kip.uni-heidelberg.de>           *
   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   AliHLTTPCCompModelConverter.cxx
  20     @author Timm Steinbeck
  21     @author changed by J. Wagner
  22     @date   17-11-2007
  23     @brief  A copy processing component for the HLT. */
  24
  25 #if __GNUC__ >= 3
  26 using namespace std;
  27 #endif
  28
  29 #include "AliHLTTPCCompModelConverter.h"
  30 #include "AliHLTTPCTransform.h"
  31 #include "AliHLTTPCTrack.h"
  32 #include "AliHLTTPCModelTrack.h"
  33 #include "AliHLTTPCCompDataCompressorHelper.h"
  34 #include "AliHLTTPCSpacePointContainer.h"
  35 #include <cerrno>
  36
  37 AliHLTTPCCompModelConverter::AliHLTTPCCompModelConverter()
  38   : AliHLTLogging()
  39   , fInputTrackArray()
  40   , fOutputTrackArray("AliHLTTPCModelTrack")
  41   , fModelAnalysisInstance(NULL)
  42   , fMinHits(0)
  43     {
  44       // see header file for class documentation
  45       for ( UInt_t slice=0; slice<36; slice++ )
  46         for ( UInt_t patch=0; patch<6; patch++ )
  47           {
  48             fClusterUsedSizes[slice][patch] = 0;
  49             fClusterUsed[slice][patch] = NULL;
  50           }
  51       Init();
  52       fMinHits = 5;
  53     }
  54
  55 AliHLTTPCCompModelConverter::AliHLTTPCCompModelConverter(AliHLTTPCCompModelAnalysis* modelanalysis)
  56   : AliHLTLogging()
  57   , fInputTrackArray()
  58   , fOutputTrackArray("AliHLTTPCModelTrack")
  59   , fModelAnalysisInstance(modelanalysis)
  60   , fMinHits(0)
  61     {
  62       // see header file for class documentation
  63       for ( UInt_t slice=0; slice<36; slice++ )
  64         for ( UInt_t patch=0; patch<6; patch++ )
  65           {
  66             fClusterUsedSizes[slice][patch] = 0;
  67             fClusterUsed[slice][patch] = NULL;
  68           }
  69       Init();
  70       fMinHits = 5;
  71     }
  72
  73 AliHLTTPCCompModelConverter::~AliHLTTPCCompModelConverter()
  74     {
  75       // see header file for class documentation
  76       for ( UInt_t slice=0; slice<36; slice++ )
  77         for ( UInt_t patch=0; patch<6; patch++ )
  78           {
  79             if ( fClusterUsed[slice][patch] )
  80               {
  81                 delete [] fClusterUsed[slice][patch];
  82                 fClusterUsed[slice][patch] = NULL;
  83               }
  84           }
  85
  86     }
  87
  88 int AliHLTTPCCompModelConverter::Init()
  89     {
  90       // see header file for class documentation
  91       fInputTrackArray.clear();
  92       fOutputTrackArray.Reset();
  93       for ( UInt_t slice=0; slice<36; slice++ )
  94         for ( UInt_t patch=0; patch<6; patch++ )
  95           fClusters[slice][patch] = NULL;
  96
  97       return 0;
  98     }
  99
 100 int AliHLTTPCCompModelConverter::SetInputTracks( const AliHLTTracksData* pTracks, unsigned sizeInByte )
 101     {
 102       // see header file for class documentation
 103       HLTDebug( "Filling %u tracks", (unsigned)pTracks->fCount );
 104       return AliHLTGlobalBarrelTrack::ConvertTrackDataArray(pTracks, sizeInByte, fInputTrackArray);
 105     }
 106
 107 int AliHLTTPCCompModelConverter::SetInputClusters( AliHLTTPCClusterData* clusters, UInt_t slice, UInt_t patch )
 108     {
 109       // set input clusters for a partition
 110       if ( slice>=36 || patch>=6 )
 111         return -EINVAL;
 112       if ( fClusters[slice][patch] )
 113         return -EBUSY;
 114       fClusters[slice][patch] = clusters;
 115       if ( fClusterUsedSizes[slice][patch]<clusters->fSpacePointCnt ||
 116            fClusterUsedSizes[slice][patch]>clusters->fSpacePointCnt*8 )
 117         {
 118           delete [] fClusterUsed[slice][patch];
 119           fClusterUsed[slice][patch] = NULL;
 120         }
 121       if ( !fClusterUsed[slice][patch] )
 122         {
 123           fClusterUsed[slice][patch] = new bool[clusters->fSpacePointCnt];
 124           if ( !fClusterUsed[slice][patch] )
 125             {
 126               return -ENOMEM;
 127             }
 128         }
 129       for ( unsigned long nn=0; nn<clusters->fSpacePointCnt; nn++ )
 130         fClusterUsed[slice][patch][nn]=false;
 131
 132       return 0;
 133     }
 134
 135 void AliHLTTPCCompModelConverter::Convert()
 136     {
 137       // convert tracks to the track model
 138
 139       // FIXME: the tracks have supposedly been sorted according to momentum to
 140       // process the high momentum tracks first, QSort however builds upon
 141       // AliTPCtrack::Compare, it's implementation has been changed
 142       //fInputTrackArray.QSort();
 143
 144       for(unsigned i=0; i<fInputTrackArray.size(); i++)
 145         {
 146           AliHLTGlobalBarrelTrack *intrack = &(fInputTrackArray[i]);
 147           if(!intrack) continue;
 148
 149           if(intrack->GetNumberOfPoints()<fMinHits)
 150             {
 151               HLTDebug("Track %d with %d clusters is below minimum of %d clusters",i,intrack->GetNHits(),fMinHits);
 152               continue;
 153             };
 154
 155           // LOSS OF TRACKS due to following statement possible!
 156           // FIXME: make this configurable
 157           if(intrack->Pt()<0.1)
 158             {
 159               HLTDebug("Discarding track with low pt.");
 160               if(fModelAnalysisInstance)
 161                 {
 162                   if(fModelAnalysisInstance->GetfModelAnalysis()) // analysis of model
 163                     {
 164                       HLTError("FIXME: correct parameter");
 165                       //fModelAnalysisInstance->MarkTrashTrack(intrack);
 166                     }
 167                 }
 168
 169               continue;
 170             }
 171
 172           // check if this is necessary
 173           //intrack->CalculateHelix();
 174
 175           AliHLTTPCModelTrack *outtrack = (AliHLTTPCModelTrack*)fOutputTrackArray.NextTrack();
 176           outtrack->SetNHits(intrack->GetNumberOfPoints());
 177           //outtrack->SetRowRange(intrack->GetFirstRow(),intrack->GetLastRow());
 178           //outtrack->SetFirstPoint(intrack->GetFirstPointX(),intrack->GetFirstPointY(),intrack->GetFirstPointZ());
 179           //outtrack->SetLastPoint(intrack->GetLastPointX(),intrack->GetLastPointY(),intrack->GetLastPointZ());
 180           outtrack->SetPt(intrack->Pt());
 181           outtrack->SetPsi(intrack->GetSnp());
 182           outtrack->SetTgl(intrack->GetTgl());
 183           // FIXME: charge is not propagated with AliHLTExternalTrackParam, check if it is needed
 184           //outtrack->SetCharge(intrack->GetCharge());
 185           outtrack->CalculateHelix();
 186           Int_t nhits = intrack->GetNumberOfPoints();
 187           const UInt_t *hitids = intrack->GetPoints();
 188           Int_t origslice = AliHLTTPCSpacePointData::GetSlice(hitids[nhits-1]);
 189           outtrack->Init(origslice,-1);
 190
 191           for(Int_t j=nhits-1; j>=0; j--)
 192             {
 193               UInt_t id=hitids[j];
 194               Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
 195               Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
 196               UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id);
 197
 198               //UInt_t size;
 199               if ( !fClusters[slice][patch] )
 200                 {
 201                   //HLTWarning( "No clusters for slice %d, patch %d", slice, patch );
 202                   continue;
 203                 }
 204             if ( !fClusterUsed[slice][patch] )
 205               {
 206                 HLTWarning( "No cluster used data for slice %d, patch %d", slice, patch );
 207                 continue;
 208               }
 209             if ( fClusters[slice][patch]->fSpacePointCnt<=pos )
 210                 {
 211                   HLTWarning( "Clusters position %d too large in slice %d, patch %d (%u max.)", pos,
 212                               slice, patch, fClusters[slice][patch]->fSpacePointCnt );
 213                   continue;
 214                 }
 215
 216             AliHLTTPCSpacePointData *points = fClusters[slice][patch]->fSpacePoints;
 217             bool* clustersUsed = fClusterUsed[slice][patch];
 218             Float_t xyz[3] = {points[pos].fX,points[pos].fY,points[pos].fZ};
 219             Int_t padrow = points[pos].fPadRow;
 220             //HLTInfo("track %d cluster %d: slice %d part %d nr %d:\n\t x=%f y=%f z=%f padrow=%d", i, j, slice, patch, pos, xyz[0], xyz[1], xyz[2], padrow);
 221
 222             //Calculate the crossing point between track and padrow
 223             Float_t angle = 0; //Perpendicular to padrow in local coordinates
 224             AliHLTTPCTransform::Local2GlobalAngle(&angle,slice);
 225             // FIXME: CalculateReferencePoint not available
 226             // if(!intrack->CalculateReferencePoint(angle,AliHLTTPCTransform::Row2X(padrow)))
 227             //   {
 228             //  HLTError( "AliHLTDataCompressor::FillData : Error in crossing point calc on slice %d, padrow %d", slice, padrow );
 229             //  break;
 230             //  //outtrack->Print(kFALSE);
 231             //  //exit(5);
 232             //   }
 233
 234             Float_t xyzCross[3] = {intrack->GetX(),intrack->GetY(),intrack->GetZ()};
 235
 236             Int_t sector,row;
 237             AliHLTTPCTransform::Slice2Sector(slice,padrow,sector,row);
 238             AliHLTTPCTransform::Global2Raw(xyzCross,sector,row);
 239 #if 1
 240             AliHLTTPCTransform::Local2Raw(xyz,sector,row);
 241 #else
 242             AliHLTTPCTransform::Global2Raw(xyz,sector,row);
 243 #endif
 244
 245             outtrack->SetPadHit(padrow,xyzCross[1]);
 246             outtrack->SetTimeHit(padrow,xyzCross[2]);
 247             //HLTInfo("track %d cluster %d: xyzCross={%f,%f,%f}", i, j, xyzCross[0], xyzCross[1], xyzCross[2]);
 248             //HLTInfo("track %d cluster %d:      xyz={%f,%f,%f}", i, j, xyz[0], xyz[1], xyz[2]);
 249
 250             //outtrack->SetCrossingAngleLUT(padrow,intrack->GetCrossingAngle(padrow,slice));
 251             outtrack->CalculateClusterWidths(padrow,kTRUE); // calculates parSigmas (with parametrisation) in raw coordinates
 252             //HLTInfo("angle %f", outtrack->GetCrossingAngleLUT(padrow));
 253             //HLTInfo("parsigma %f",outtrack->GetParSigmaY2(padrow));
 254             patch = AliHLTTPCTransform::GetPatch(padrow);
 255             // sigmay in units of pads (quantisation!)
 256             Float_t sigmaY2 = points[pos].fSigmaY2 / pow(AliHLTTPCTransform::GetPadPitchWidth(patch),2);
 257             //HLTInfo("sigmaY conv.: %f", points[pos].fSigmaY2);
 258
 259             //HLTInfo("parSigmaY2 = %f",AliHLTTPCTransform::GetParSigmaY2(padrow, xyzCross[2],angle));
 260             //Float_t testsigma = 0.0;
 261             //outtrack->GetSigmaY2(padrow, testsigma);
 262             //HLTInfo("DSigmaY2 = %f",testsigma);
 263
 264             //HLTInfo("sigmaY2 float: %f",sigmaY2);
 265             Float_t sigmaZ2 = points[pos].fSigmaZ2 / pow(AliHLTTPCTransform::GetZWidth(),2);
 266             outtrack->SetCluster(padrow,xyz[1],xyz[2],points[pos].fCharge,sigmaY2,sigmaZ2,3);
 267             //AliHLTTPCClusterModel* test1 = outtrack->GetClusterModel(padrow);
 268             //HLTInfo("Dsigma %f",test1->fDSigmaY);
 269
 270             //IMPORTANT: Set the slice in which cluster is, you need it in AliHLTTPCModelTrack::FillTrack!
 271             outtrack->GetClusterModel(padrow)->fSlice=slice;
 272 #ifdef MODELDEBUG
 273             outtrack->GetClusterModel(padrow)->fID=points[pos].fID;
 274             HLTDebug( "Track %d cluster for padrow %d ID: %u (0x%08X) - fSlice: %u", i, padrow,
 275                       outtrack->GetClusterModel(padrow)->fID, outtrack->GetClusterModel(padrow)->fID,
 276                       (unsigned)outtrack->GetClusterModel(padrow)->fSlice );
 277 #endif
 278             //points[pos].fCharge = 0;//Mark this cluster as used.
 279             clustersUsed[pos] = true;//Mark this cluster as used.
 280             //fNusedClusters++;
 281             } //end of clusters for each track
 282
 283         //outtrack->SetNClusters(AliHLTTPCTransform::GetNRows(-1)); // Equivalent call in ExpandTrackData
 284         } // end of track-loop
 285     ExpandTrackData();
 286
 287     // validation test for clusternumbers of tracks:
 288     //for(unsigned long jj = 0; jj < (unsigned long) fOutputTrackArray.GetNTracks(); jj++)
 289     //  {
 290     //  AliHLTTPCModelTrack *track = (AliHLTTPCModelTrack*)fOutputTrackArray.GetCheckedTrack(jj);
 291     //  Int_t nhits = track->GetNHits();
 292     //  HLTInfo("Number of clusters for track %lu is %d",jj, nhits);
 293     //  }
 294
 295     //comp->WriteFile(fOutputTrackArray);
 296
 297     }
 298
 299 void AliHLTTPCCompModelConverter::ExpandTrackData()
 300     {
 301       // see header file for class documentation
 302       //Loop over tracks and try to assign unused clusters.
 303       //Only clusters which are closer than the max. residual are taken.
 304
 305       HLTDebug( "Expanding %lu tracks", (unsigned long)fOutputTrackArray.GetNTracks() );
 306       for(Int_t i=0; i<fOutputTrackArray.GetNTracks(); i++)
 307         {
 308           AliHLTTPCModelTrack *track = (AliHLTTPCModelTrack*)fOutputTrackArray.GetCheckedTrack(i);
 309
 310           if(!track) continue;
 311
 312           // tracks that hit every row already cannot be expanded in the current model!
 313           if(track->GetNHits() == AliHLTTPCTransform::GetNRows()) continue;
 314
 315           Int_t nhits = track->GetNHits();
 316
 317           // validation test
 318           //HLTInfo("Before expansion: track %u with number of clusters %d", i, nhits);
 319
 320           Int_t lastSlice=-1;
 321           for(Int_t padrow=AliHLTTPCTransform::GetNRows()-1; padrow>=0; padrow--)
 322             {
 323               if(track->IsPresent(padrow))
 324                 {
 325                   lastSlice = track->GetClusterModel(padrow)->fSlice;
 326                   continue;
 327                 }
 328
 329               if(lastSlice < 0) //the outer cluster is missing, so skip it - it will be written anyhow.
 330                 continue;
 331
 332               //Check the slice of the next padrow:
 333               Int_t nextPadrow = padrow-1;
 334               Int_t nextSlice = -1;
 335               while(nextPadrow >=0)
 336                 {
 337                   if(track->IsPresent(nextPadrow))
 338                     {
 339                       nextSlice = track->GetClusterModel(nextPadrow)->fSlice;
 340                       break;
 341                     }
 342                   nextPadrow--;
 343                 }
 344               if(nextSlice>=0)
 345                 if(nextSlice != lastSlice)//The track crosses a slice boundary here
 346                   continue;
 347
 348               //UInt_t size;
 349               if ( !fClusters[lastSlice][0] )
 350                 {
 351                   HLTWarning( "No clusters for slice %d, patch %d", lastSlice, 0 );
 352                   continue;
 353                 }
 354               if ( !fClusterUsed[lastSlice][0] )
 355                 {
 356                   HLTWarning( "No cluster used data for slice %d, patch %d", lastSlice, 0 );
 357                   continue;
 358                 }
 359               AliHLTTPCSpacePointData *points = fClusters[lastSlice][0]->fSpacePoints;//->GetDataPointer(size);
 360               bool* clustersUsed = fClusterUsed[lastSlice][0];
 361
 362               Float_t globalangle = 0;
 363               AliHLTTPCTransform::Local2GlobalAngle(&globalangle,lastSlice);
 364               if(!track->CalculateReferencePoint(globalangle,AliHLTTPCTransform::Row2X(padrow)))
 365                 continue;
 366               Float_t xyzCross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
 367               AliHLTTPCTransform::Global2LocHLT(xyzCross,lastSlice);
 368               Float_t mindist = 123456789;
 369               AliHLTTPCSpacePointData *closest=0;
 370               UInt_t closestJ=0;
 371               for(UInt_t j=0; j<fClusters[lastSlice][0]->fSpacePointCnt; j++)
 372                 {
 373                   //if(points[j].fCharge == 0) continue;// || points[j].fPadRow != padrow) continue;
 374                   if (clustersUsed[j]) continue; // Cluster already used
 375                   if(points[j].fPadRow < padrow) continue;
 376                   if(points[j].fPadRow > padrow) break;
 377                   Float_t xyz[3] = {points[j].fX,points[j].fY,points[j].fZ};
 378 #if 1
 379 #else
 380                   AliHLTTPCTransform::Global2LocHLT(xyz,lastSlice);
 381 #endif
 382
 383                   //Check for overflow:
 384                   Int_t temp = (Int_t)rint((xyzCross[1]-xyz[1])/AliHLTTPCCompDataCompressorHelper::GetXYResidualStep(padrow));
 385                   if( abs(temp) > 1<<(AliHLTTPCCompDataCompressorHelper::GetNPadBits()-1))
 386                     continue;
 387
 388                   temp = (Int_t)rint((xyzCross[2]-xyz[2])/AliHLTTPCCompDataCompressorHelper::GetZResidualStep(padrow));
 389                   if( abs(temp) > 1<<(AliHLTTPCCompDataCompressorHelper::GetNTimeBits()-1))
 390                     continue;
 391
 392                   Float_t dist = sqrt( pow(xyzCross[1]-xyz[1],2) + pow(xyzCross[2]-xyz[2],2) );
 393                   if(dist < mindist)
 394                     {
 395                       closest = &points[j];
 396                       closestJ = j;
 397                       mindist = dist;
 398                     }
 399                 }
 400               if(closest) //there was a cluster assigned
 401                 {
 402                   Int_t sector,row;
 403                   Float_t xyz[3] = {closest->fX,closest->fY,closest->fZ};
 404                   AliHLTTPCTransform::Slice2Sector(lastSlice,padrow,sector,row);
 405                   AliHLTTPCTransform::Local2Raw(xyzCross,sector,row);
 406 #if 1
 407                   AliHLTTPCTransform::Local2Raw(xyz,sector,row);
 408 #else
 409                   AliHLTTPCTransform::Global2Raw(xyz,sector,row);
 410 #endif
 411
 412                   track->SetPadHit(padrow,xyzCross[1]);
 413                   track->SetTimeHit(padrow,xyzCross[2]);
 414
 415                   Float_t angle = track->GetCrossingAngle(padrow,lastSlice);
 416                   track->SetCrossingAngleLUT(padrow,angle);
 417                   track->CalculateClusterWidths(padrow,kTRUE);
 418                   Int_t patch = AliHLTTPCTransform::GetPatch(padrow);
 419                   Float_t sigmaY2 = closest->fSigmaY2 / pow(AliHLTTPCTransform::GetPadPitchWidth(patch),2);
 420                   Float_t sigmaZ2 = closest->fSigmaZ2 / pow(AliHLTTPCTransform::GetZWidth(),2);
 421                   track->SetCluster(padrow,xyz[1],xyz[2],closest->fCharge,sigmaY2,sigmaZ2,3);
 422                   //AliHLTTPCClusterModel* test1 = track->GetClusterModel(padrow);
 423                   //HLTInfo("Dsigma %f",test1->fDSigmaY);
 424
 425                   nhits++;
 426
 427                   //IMPORTANT: Set the slice in which cluster is, you need it in AliHLTTPCModelTrack::FillTrack!
 428                   track->GetClusterModel(padrow)->fSlice=lastSlice;
 429 #ifdef MODELDEBUG
 430                   track->GetClusterModel(padrow)->fID=closest->fID;
 431                   HLTDebug( "Track %d cluster for padrow %d ID: %u (0x%08X) - fSlice: %u", i, padrow,
 432                             track->GetClusterModel(padrow)->fID, track->GetClusterModel(padrow)->fID,
 433                             track->GetClusterModel(padrow)->fSlice );
 434 #endif
 435                   //closest->fCharge = 0;//Mark this cluster as used.
 436                   clustersUsed[closestJ] = true;//Mark this cluster as used.
 437                 }
 438             }
 439           track->SetNClusters(AliHLTTPCTransform::GetNRows());
 440           //cout<<"Track was assigned "<<nhits<<" clusters"<<endl;
 441
 442           // validation test
 443           //HLTInfo( "After expansion: track %d with clusters %u", i, nhits);
 444         }
 445
 446     }
 447
 448 unsigned long AliHLTTPCCompModelConverter::GetOutputModelDataSize() const
 449     {
 450       // see header file for class documentation
 451       unsigned long dataSize=0;
 452       Short_t ntracks = fOutputTrackArray.GetNTracks();
 453
 454       dataSize += sizeof(AliHLTUInt32_t);
 455
 456       for(Int_t i=0; i<ntracks; i++)
 457         {
 458           AliHLTTPCModelTrack *track = (AliHLTTPCModelTrack*)fOutputTrackArray.GetCheckedTrack(i);
 459           if ( !track )
 460             continue;
 461
 462           dataSize += sizeof(AliHLTTPCTrackModel)+track->GetNClusters()*sizeof(AliHLTTPCClusterModel);
 463         }
 464       return dataSize;
 465     }
 466
 467 int AliHLTTPCCompModelConverter::OutputModelData( AliHLTUInt8_t* data, unsigned long& /*dataSize*/ ) const
 468     {
 469       // see header file for class documentation
 470       unsigned long dataOffset=0;
 471       Short_t ntracks = fOutputTrackArray.GetNTracks();
 472
 473       AliHLTTPCClusterModel *clusters=0;
 474       AliHLTTPCTrackModel *model=0;
 475
 476       *(AliHLTUInt32_t*)data = 0; // Write format version number
 477       dataOffset += sizeof(AliHLTUInt32_t);
 478
 479       for(Int_t i=0; i<ntracks; i++)
 480         {
 481           AliHLTTPCModelTrack *track = (AliHLTTPCModelTrack*)fOutputTrackArray.GetCheckedTrack(i);
 482           if ( !track )
 483             continue;
 484
 485           track->FillModel();
 486           model = track->GetModel();
 487
 488           clusters = track->GetClusters();
 489
 490           // validation test
 491           //HLTInfo( "Track %d clusters: %u", i, (unsigned)track->GetNPresentClusters() );
 492
 493           for ( Int_t jj=0; jj<track->GetNClusters(); jj++ )
 494             {
 495               //HLTDebug( "  Cluster %d fPresent: %u", jj, (unsigned)clusters[jj].fPresent );
 496             }
 497
 498           memcpy( data+dataOffset, model, sizeof(AliHLTTPCTrackModel) );
 499           dataOffset += sizeof(AliHLTTPCTrackModel);
 500
 501           memcpy( data+dataOffset, clusters, track->GetNClusters()*sizeof(AliHLTTPCClusterModel) );
 502           dataOffset += track->GetNClusters()*sizeof(AliHLTTPCClusterModel);
 503         }
 504       return 0;
 505     }
 506
 507 void AliHLTTPCCompModelConverter::SelectRemainingClusters()
 508     {
 509       // see header file for class documentation
 510       //Select which remaining clusters to write in addition to the compressed data.
 511       //In particular one can here make sure that "important" clusters are not missed:
 512       //The offline track finder perform seed finding in the outer padrows;
 513       //the first seeding is using pair of points on outermost padrow and
 514       //0.125*nrows more rows towards the vertex. The second seeding uses pair
 515       //of points on the outermost padrow-0.5*0.125*nrows and 0.125*nrows + 0.5*0.125*nrows
 516       //more rows towards the vertex. In order to evaluate the seeds, the track offline
 517       //track finder checks whether a certain amount of possible clusters (padrows) is
 518       //attached to the track, and then the kalman filtering starts.
 519       //To ensure a minimal loss off efficiency, all clusters in this region should be
 520       //intact.....
 521
 522       Int_t nrows = AliHLTTPCTransform::GetNRows();
 523       Int_t gap=(Int_t)(0.125*nrows), shift=(Int_t)(0.5*gap);
 524
 525       for(Int_t slice=0; slice<36; slice++)
 526         {
 527           for(Int_t patch=0; patch<6; patch++)
 528             {
 529               if ( !fClusters[slice][patch] )
 530                 continue;
 531               AliHLTTPCSpacePointData *points = fClusters[slice][patch]->fSpacePoints;
 532               bool* clustersUsed = fClusterUsed[slice][patch];
 533               for(UInt_t i=0; i<fClusters[slice][patch]->fSpacePointCnt; i++)
 534                 {
 535                   //if(points[i].fCharge == 0) continue; //Already removed
 536                   if (clustersUsed[i]) continue; //Already removed
 537                   Int_t padrow = (Int_t)points[i].fPadRow;
 538
 539                   //Check the widths (errors) of the cluster, and remove big bastards:
 540                   Float_t padw = sqrt(points[i].fSigmaY2) / AliHLTTPCTransform::GetPadPitchWidth(AliHLTTPCTransform::GetPatch(padrow));
 541                   Float_t timew = sqrt(points[i].fSigmaZ2) / AliHLTTPCTransform::GetZWidth();
 542                   if(padw >= 2.55 || timew >= 2.55)//Because we use 1 byte to store
 543                     {
 544                       //points[i].fCharge = 0;
 545                       clustersUsed[i] = true;
 546                       continue;
 547                     }
 548
 549                   Float_t xyz[3] = {points[i].fX,points[i].fY,points[i].fZ};
 550                   Int_t sector,row;
 551                   AliHLTTPCTransform::Slice2Sector(slice,padrow,sector,row);
 552                   AliHLTTPCTransform::Global2Raw(xyz,sector,row);
 553
 554                   if(padrow >= nrows-1-gap-shift) continue;//save all the clusters in this region
 555
 556                   //if(padrow >= nrows-1-shift) continue;
 557
 558                   //Save the clusters at the borders:
 559                   //if(xyz[1] < 3 || xyz[1] >= AliHLTTPCTransform::GetNPads(padrow)-4)
 560                   // continue;
 561
 562                   //Save clusters on padrows used for offline seeding:
 563                   if(padrow == nrows - 1 || padrow == nrows - 1 - gap ||                 //First seeding
 564                      padrow == nrows - 1 - shift || padrow == nrows - 1 - gap - shift)   //Second seeding
 565                     continue;
 566
 567                   //Cluster did not meet any of the above criteria, so disregard it:
 568                   //points[i].fCharge = 0;
 569                   clustersUsed[i] = true;
 570                 }
 571             }
 572         }
 573
 574     }
 575
 576 unsigned long AliHLTTPCCompModelConverter::GetRemainingClustersOutputDataSize() const
 577     {
 578       // see header file for class documentation
 579       int iResult=0;
 580 #if 0
 581     for ( UInt_t slice=0; slice<36; slice++ )
 582         for ( UInt_t patch=0; patch<6; patch++ )
 583             {
 584             bool* clustersUsed = fClusterUsed[slice][patch];
 585             if ( !clustersUsed || !fClusters[slice][patch] )
 586                 continue;
 587             for ( UInt_t pos=0; pos<fClusters[slice][patch]->fSpacePointCnt; pos++ )
 588                 {
 589                 if ( !clustersUsed[pos] )
 590                     clusterCnt++;
 591                 }
 592             }
 593     return clusterCnt*sizeof(AliHLTTPCClusterModel);
 594 #else
 595     const Int_t nrows = AliHLTTPCTransform::GetNRows();
 596     Int_t * npoints = new Int_t[nrows];
 597     unsigned long dataWritten = 0;
 598
 599     dataWritten += sizeof(AliHLTUInt32_t);
 600
 601     // FIXME: get rid of hardcoded numbers
 602     for(Int_t slice=0; slice<35 && iResult>=0; slice++)
 603         {
 604         for(Int_t patch=0; patch < 6 && iResult>=0; patch++)
 605             {
 606             if ( !fClusters[slice][patch] )
 607                 {
 608                 dataWritten++;
 609                 continue;
 610                 }
 611             AliHLTTPCSpacePointData *points = fClusters[slice][patch]->fSpacePoints;
 612             bool* clustersUsed = fClusterUsed[slice][patch];
 613             if ( !clustersUsed )
 614                 continue;
 615             memset(npoints,0,nrows*sizeof(Int_t));
 616             Int_t nonZeroRows=0;
 617
 618             for(UInt_t j=0; j<fClusters[slice][patch]->fSpacePointCnt; j++)
 619                 {
 620                 //if(points[j].fCharge == 0) continue; //has been used
 621                 if ( clustersUsed[j] ) continue; //has been used
 622                 if ( !npoints[points[j].fPadRow] )
 623                     nonZeroRows++;
 624                 npoints[points[j].fPadRow]++;
 625                 }
 626
 627             dataWritten++;
 628
 629             Int_t size =0;
 630             Byte_t *data = 0;
 631             AliHLTTPCRemainingRow *tempPt=0;
 632
 633             Int_t lastRow = -2;
 634             Int_t localcounter=0;
 635
 636             for(UInt_t j=0; j<fClusters[slice][patch]->fSpacePointCnt; j++)
 637                 {
 638                 //if(points[j].fCharge == 0) continue; //has been used
 639                 if ( clustersUsed[j] ) continue; //has been used
 640
 641                 Int_t padrow = points[j].fPadRow;
 642                 if(padrow != lastRow)
 643                     {
 644                     if(lastRow != -2)
 645                         {
 646                         if(!tempPt)
 647                             {
 648                             HLTError( "Zero row pointer " );
 649                             iResult=-EINVAL;
 650                             break;
 651                             }
 652                         if(localcounter != tempPt->fNClusters)
 653                             {
 654                             HLTError( "Mismatching clustercounter %lu - %d ",
 655                                       (unsigned long)localcounter, (Int_t)tempPt->fNClusters );
 656                             iResult=EINVAL;
 657                             break;
 658                             }
 659                         dataWritten += size;
 660                         }
 661                     if(data)
 662                         delete [] data;
 663                     size = sizeof(AliHLTTPCRemainingRow) + npoints[padrow]*sizeof(AliHLTTPCRemainingCluster);
 664                     data = new Byte_t[size];
 665                     tempPt = reinterpret_cast<AliHLTTPCRemainingRow*>(data);
 666
 667                     localcounter=0;
 668                     tempPt->fPadRow = padrow;
 669                     tempPt->fNClusters = npoints[padrow];
 670                     lastRow = padrow;
 671                     }
 672                 if(localcounter >= npoints[padrow])
 673                     {
 674                     HLTError( "Cluster counter out of range: %lu - %lu",
 675                               (unsigned long)localcounter, (unsigned long)npoints[padrow] );
 676                     iResult=-EINVAL;
 677                     break;
 678                     }
 679
 680                 localcounter++;
 681                 }
 682
 683             //Write the last row:
 684             if ( tempPt )
 685                 {
 686                 dataWritten += size;
 687                 }
 688             if(data)
 689               delete [] data;
 690             }
 691         }
 692     delete [] npoints;
 693     // FIXME check the caller and propagate an error condition
 694     if (iResult<0) return 0;
 695     return dataWritten;
 696 #endif
 697     }
 698
 699 int AliHLTTPCCompModelConverter::GetRemainingClusters( AliHLTUInt8_t* const pTgt, unsigned long& dataSize ) const
 700     {
 701       // see header file for class documentation
 702       int iResult=0;
 703
 704       // FIXME: almost identical code to  GetRemainingClustersOutputDataSize
 705       // try to combine
 706       const Int_t nrows = AliHLTTPCTransform::GetNRows();
 707       Int_t * npoints = new Int_t[nrows];
 708       unsigned long dataWritten = 0;
 709       AliHLTUInt8_t* writePtr = pTgt;
 710
 711       *(AliHLTUInt32_t*)writePtr = 0; // Write format version
 712       dataWritten += sizeof(AliHLTUInt32_t);
 713       writePtr += sizeof(AliHLTUInt32_t);
 714
 715       for(Int_t slice=0; slice<=35 && iResult>=0; slice++)
 716         {
 717           for(Int_t patch=0; patch < 6 && iResult>=0; patch++)
 718             {
 719               if ( !fClusters[slice][patch] )
 720                 {
 721                   *writePtr = (AliHLTUInt8_t)0;
 722                   writePtr++;
 723                   dataWritten++;
 724                   continue;
 725                 }
 726               AliHLTTPCSpacePointData *points = fClusters[slice][patch]->fSpacePoints;
 727               bool* clustersUsed = fClusterUsed[slice][patch];
 728               if ( !clustersUsed )
 729                 continue;
 730               memset(npoints,0,nrows*sizeof(Int_t));
 731               Int_t nonZeroRows=0;
 732
 733               for(UInt_t j=0; j<fClusters[slice][patch]->fSpacePointCnt; j++)
 734                 {
 735                   //if(points[j].fCharge == 0) continue; //has been used
 736                   if ( clustersUsed[j] ) continue; //has been used
 737                   if ( !npoints[points[j].fPadRow] )
 738                     nonZeroRows++;
 739                   npoints[points[j].fPadRow]++;
 740                 }
 741
 742               *writePtr = (AliHLTUInt8_t)nonZeroRows;
 743               writePtr++;
 744               dataWritten++;
 745
 746               Int_t size =0;
 747               Byte_t *data = 0;
 748               AliHLTTPCRemainingRow *tempPt=0;
 749
 750               Int_t lastRow = -2;
 751               Int_t localcounter=0;
 752
 753               for(UInt_t j=0; j<fClusters[slice][patch]->fSpacePointCnt; j++)
 754                 {
 755                   //if(points[j].fCharge == 0) continue; //has been used
 756                   if ( clustersUsed[j] ) continue; //has been used
 757
 758                   Int_t padrow = points[j].fPadRow;
 759                   if(padrow != lastRow)
 760                     {
 761                       if(lastRow != -2)
 762                         {
 763                           if(!tempPt)
 764                             {
 765                               HLTError( "Zero row pointer " );
 766                               iResult=-EINVAL;
 767                               break;
 768                             }
 769                           if(localcounter != tempPt->fNClusters)
 770                             {
 771                               HLTError( "Mismatching clustercounter %lu - %d ",
 772                                         (unsigned long)localcounter, (Int_t)tempPt->fNClusters );
 773                               iResult=-EINVAL;
 774                               break;
 775                             }
 776                           //cout<<"Writing row "<<(int)tempPt->fPadRow<<" with "<<(int)tempPt->fNClusters<<" clusters"<<endl;
 777                           //fwrite(tempPt,size,1,outfile);
 778                           if ( dataWritten+size > dataSize )
 779                             {
 780                               HLTWarning( "Cannot write remaining clusters to output. Data size too large (exceeding %lu bytes)", (unsigned long)dataSize );
 781                               iResult=-ENOBUFS;
 782                               break;
 783                             }
 784                           memcpy( writePtr, tempPt, size );
 785                           dataWritten += size;
 786                           writePtr += size;
 787                         }
 788                       if(data)
 789                         delete [] data;
 790                       size = sizeof(AliHLTTPCRemainingRow) + npoints[padrow]*sizeof(AliHLTTPCRemainingCluster);
 791                       data = new Byte_t[size];
 792                       tempPt = (AliHLTTPCRemainingRow*)data;
 793
 794                       localcounter=0;
 795                       tempPt->fPadRow = padrow;
 796                       tempPt->fNClusters = npoints[padrow];
 797                       lastRow = padrow;
 798                     }
 799                   if(localcounter >= npoints[padrow])
 800                     {
 801                       HLTError( "Cluster counter out of range: %lu - %lu",
 802                                 (unsigned long)localcounter, (unsigned long)npoints[padrow] );
 803                       iResult=EINVAL;
 804                       break;
 805                     }
 806
 807                   Float_t xyz[3] = {points[j].fX,points[j].fY,points[j].fZ};
 808                   Int_t sector,row;
 809                   AliHLTTPCTransform::Slice2Sector(slice,padrow,sector,row);
 810 #if 1
 811                   AliHLTTPCTransform::Local2Raw(xyz,sector,row);
 812 #else
 813                   AliHLTTPCTransform::Global2Raw(xyz,sector,row);
 814 #endif
 815
 816                   Float_t padw = points[j].fSigmaY2 / pow(AliHLTTPCTransform::GetPadPitchWidth(AliHLTTPCTransform::GetPatch(padrow)),2);
 817                   Float_t timew = points[j].fSigmaZ2 / pow(AliHLTTPCTransform::GetZWidth(),2);
 818                   tempPt->fClusters[localcounter].fPad = xyz[1];
 819                   tempPt->fClusters[localcounter].fTime = xyz[2];
 820                   tempPt->fClusters[localcounter].fCharge = points[j].fCharge;
 821                   tempPt->fClusters[localcounter].fSigmaY2 = padw;
 822                   tempPt->fClusters[localcounter].fSigmaZ2 = timew;
 823 #ifdef MODELDEBUG
 824                   tempPt->fClusters[localcounter].fID = points[j].fID;
 825 #endif
 826                   localcounter++;
 827                   if(fModelAnalysisInstance)
 828                     {
 829                       if(fModelAnalysisInstance->GetfModelAnalysis())
 830                         {
 831                           fModelAnalysisInstance->MarkTrashCluster(fClusters[slice][patch], slice, patch);
 832                         }
 833                     }
 834                 }
 835
 836
 837               //Write the last row:
 838               if ( dataWritten+size > dataSize )
 839                 {
 840                 HLTWarning( "Cannot write remaining clusters to output. Data size too large (exceeding %lu bytes)", (unsigned long)dataSize );
 841                 iResult=-ENOBUFS;
 842                 if(data)
 843                   delete [] data;
 844                 break;
 845                 }
 846               if ( tempPt )
 847                 {
 848                   memcpy( writePtr, tempPt, size );
 849                   dataWritten += size;
 850                   writePtr += size;
 851                 }
 852               if(data)
 853                 delete [] data;
 854             }
 855         }
 856       dataSize = dataWritten;
 857
 858       delete [] npoints;
 859       return iResult;
 860     }