HLT/TPCLib/AliHLTTPCClusterFinder.cxx

   1 // @(#) $Id$
   2 // Original: AliL3ClustFinderNew.cxx,v 1.29 2005/06/14 10:55:21 cvetan Exp
   3
   4 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>,
   5 //         Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
   6 //         Jochen Thaeder <mailto:thaeder@kip.uni-heidelberg.de>
   7
   8 //*-- Copyright &copy ALICE HLT Group
   9
  10
  11 #include "AliHLTTPCDigitReader.h"
  12 #include "AliHLTTPCRootTypes.h"
  13 #include "AliHLTTPCLogging.h"
  14 #include "AliHLTTPCClusterFinder.h"
  15 #include "AliHLTTPCDigitData.h"
  16 #include "AliHLTTPCTransform.h"
  17 #include "AliHLTTPCSpacePointData.h"
  18 #include "AliHLTTPCMemHandler.h"
  19 #include "AliHLTTPCPad.h"
  20
  21 #if __GNUC__ >= 3
  22 using namespace std;
  23 #endif
  24
  25 /** \class AliHLTTPCClusterFinder
  26 <pre>
  27 //_____________________________________________________________
  28 // AliHLTTPCClusterFinder
  29 //
  30 // The current cluster finder for HLT
  31 // (Based on STAR L3)
  32 //
  33 // The cluster finder is initialized with the Init function,
  34 // providing the slice and patch information to work on.
  35 //
  36 // The input is a provided by the AliHLTTPCDigitReader class,
  37 // using the init() funktion, and the next() funktion in order
  38 // to get the next bin. Either packed or unpacked data can be
  39 // processed, dependent if one uses AliHLTTPCDigitReaderPacked
  40 // class or AliHLTTPCDigitReaderUnpacked class in the
  41 // Clusterfinder Component.
  42 // The resulting space points will be in the
  43 // array given by the SetOutputArray function.
  44 //
  45 // There are several setters which control the behaviour:
  46 //
  47 // - SetXYError(Float_t):   set fixed error in XY direction
  48 // - SetZError(Float_t):    set fixed error in Z  direction
  49 //                            (used if errors are not calculated)
  50 // - SetDeconv(Bool_t):     switch on/off deconvolution
  51 // - SetThreshold(UInt_t):  set charge threshold for cluster
  52 // - SetMatchWidth(UInt_t): set the match distance in
  53 //                            time for sequences to be merged
  54 // - SetSTDOutput(Bool_t):  switch on/off output about found clusters
  55 // - SetCalcErr(Bool_t):    switch on/off calculation of
  56 //                          space point errors (or widths in raw system)
  57 // - SetRawSP(Bool_t):      switch on/off convertion to raw system
  58 //
  59 //
  60 // Example Usage:
  61 //
  62 // AliHLTTPCFileHandler *file = new AliHLTTPCFileHandler();
  63 // file->SetAliInput(digitfile); //give some input file
  64 // for(int slice=0; slice<=35; slice++){
  65 //   for(int patch=0; pat<6; pat++){
  66 //     file->Init(slice,patch);
  67 //     UInt_t ndigits=0;
  68 //     UInt_t maxclusters=100000;
  69 //     UInt_t pointsize = maxclusters*sizeof(AliHLTTPCSpacePointData);
  70 //     AliHLTTPCSpacePointData *points = (AliHLTTPCSpacePointData*)memory->Allocate(pointsize);
  71 //     AliHLTTPCDigitRowData *digits = (AliHLTTPCDigitRowData*)file->AliAltroDigits2Memory(ndigits,event);
  72 //     AliHLTTPCClusterFinder *cf = new AliHLTTPCClusterFinder();
  73 //     cf->SetMatchWidth(2);
  74 //     cf->InitSlice( slice, patch, row[0], row[1], maxPoints );
  75 //     cf->SetSTDOutput(kTRUE);    //Some output to standard IO
  76 //     cf->SetRawSP(kFALSE);       //Convert space points to local system
  77 //     cf->SetThreshold(5);        //Threshold of cluster charge
  78 //     cf->SetDeconv(kTRUE);       //Deconv in pad and time direction
  79 //     cf->SetCalcErr(kTRUE);      //Calculate the errors of the spacepoints
  80 //     cf->SetOutputArray(points); //Move the spacepoints to the array
  81 //     cf->Read(iter->fPtr, iter->fSize ); //give the data to the cf
  82 //     cf->ProcessDigits();        //process the rows given by init
  83 //     Int_t npoints = cf->GetNumberOfClusters();
  84 //     AliHLTTPCMemHandler *out= new AliHLTTPCMemHandler();
  85 //     out->SetBinaryOutput(fname);
  86 //     out->Memory2Binary(npoints,points); //store the spacepoints
  87 //     out->CloseBinaryOutput();
  88 //     delete out;
  89 //     file->free();
  90 //     delete cf;
  91 //   }
  92 // }
  93 </pre>
  94 */
  95
  96 ClassImp(AliHLTTPCClusterFinder)
  97
  98 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
  99   :
 100   fMatch(1),
 101   fThreshold(10),
 102   fSignalThreshold(-1),
 103   fXYErr(0.2),
 104   fZErr(0.3),
 105   fDeconvPad(kTRUE),
 106   fDeconvTime(kTRUE),
 107   fStdout(kFALSE),
 108   fCalcerr(kTRUE),
 109   fRawSP(kFALSE),
 110   fFirstRow(0),
 111   fLastRow(0),
 112   fDigitReader(NULL)
 113 {
 114   //constructor
 115 }
 116
 117 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder(const AliHLTTPCClusterFinder& src)
 118   :
 119   fMatch(src.fMatch),
 120   fThreshold(src.fThreshold),
 121   fSignalThreshold(src.fSignalThreshold),
 122   fXYErr(src.fXYErr),
 123   fZErr(src.fZErr),
 124   fDeconvPad(src.fDeconvPad),
 125   fDeconvTime(src.fDeconvTime),
 126   fStdout(src.fStdout),
 127   fCalcerr(src.fCalcerr),
 128   fRawSP(src.fRawSP),
 129   fFirstRow(src.fFirstRow),
 130   fLastRow(src.fLastRow),
 131   fDigitReader(src.fDigitReader)
 132 {
 133 }
 134
 135 AliHLTTPCClusterFinder& AliHLTTPCClusterFinder::operator=(const AliHLTTPCClusterFinder& src)
 136 {
 137   fMatch=src.fMatch;
 138   fThreshold=src.fThreshold;
 139   fSignalThreshold=src.fSignalThreshold;
 140   fXYErr=src.fXYErr;
 141   fZErr=src.fZErr;
 142   fDeconvPad=src.fDeconvPad;
 143   fDeconvTime=src.fDeconvTime;
 144   fStdout=src.fStdout;
 145   fCalcerr=src.fCalcerr;
 146   fRawSP=src.fRawSP;
 147   fFirstRow=src.fFirstRow;
 148   fLastRow=src.fLastRow;
 149   fDigitReader=src.fDigitReader;
 150   return (*this);
 151 }
 152
 153 AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder()
 154 {
 155   //destructor
 156 }
 157
 158 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
 159 {
 160   //init slice
 161   fNClusters = 0;
 162   fMaxNClusters = nmaxpoints;
 163   fCurrentSlice = slice;
 164   fCurrentPatch = patch;
 165   fFirstRow = firstrow;
 166   fLastRow = lastrow;
 167 }
 168
 169 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
 170 {
 171   //init slice
 172   fNClusters = 0;
 173   fMaxNClusters = nmaxpoints;
 174   fCurrentSlice = slice;
 175   fCurrentPatch = patch;
 176   fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
 177   fLastRow=AliHLTTPCTransform::GetLastRow(patch);
 178 }
 179
 180 void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt)
 181 {
 182   //set pointer to output
 183   fSpacePointData = pt;
 184 }
 185
 186 void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
 187   //set input pointer
 188   fPtr = (UChar_t*)ptr;
 189   fSize = size;
 190 }
 191
 192 void AliHLTTPCClusterFinder::ProcessDigits()
 193 {
 194   bool readValue = true;
 195   Int_t newRow = 0;
 196   Int_t rowOffset = 0;
 197   UShort_t time=0,newTime=0;
 198   UInt_t pad=0,newPad=0;
 199   AliHLTTPCSignal_t charge=0;
 200
 201   fNClusters = 0;
 202
 203   // initialize block for reading packed data
 204   fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
 205   readValue = fDigitReader->Next();
 206
 207   if (!readValue)return;
 208
 209   pad = fDigitReader->GetPad();
 210   time = fDigitReader->GetTime();
 211   fCurrentRow = fDigitReader->GetRow();
 212
 213   if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
 214     rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
 215
 216   fCurrentRow += rowOffset;
 217
 218   UInt_t lastpad = 123456789;
 219   const Int_t kPadArraySize=5000;
 220   const Int_t kClusterListSize=10000;
 221   AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
 222   AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
 223   AliClusterData clusterlist[kClusterListSize]; //Clusterlist
 224
 225   AliClusterData **currentPt;  //List of pointers to the current pad
 226   AliClusterData **previousPt; //List of pointers to the previous pad
 227   currentPt = pad2;
 228   previousPt = pad1;
 229   UInt_t nprevious=0,ncurrent=0,ntotal=0;
 230
 231   /* quick implementation of baseline calculation and zero suppression
 232      open a pad object for each pad and delete it after processing.
 233      later a list of pad objects with base line history can be used
 234      The whole thing only works if we really get unprocessed raw data, if
 235      the data is already zero suppressed, there might be gaps in the time
 236      bins.
 237    */
 238   Int_t gatingGridOffset=50;
 239   AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
 240   // just to make later conversion to a list of objects easier
 241   AliHLTTPCPad* pCurrentPad=NULL;
 242   if (fSignalThreshold>=0) {
 243     pCurrentPad=&baseline;
 244     baseline.SetThreshold(fSignalThreshold);
 245   }
 246
 247   while ( readValue ){   // Reads through all digits in block
 248
 249     if(pad != lastpad){
 250       //This is a new pad
 251
 252       //Switch the lists:
 253       if(currentPt == pad2){
 254         currentPt = pad1;
 255         previousPt = pad2;
 256       }
 257       else {
 258         currentPt = pad2;
 259         previousPt = pad1;
 260       }
 261       nprevious = ncurrent;
 262       ncurrent = 0;
 263       if(pad != lastpad+1){
 264         //this happens if there is a pad with no signal.
 265         nprevious = ncurrent = 0;
 266       }
 267       lastpad = pad;
 268     }
 269
 270     Bool_t newcluster = kTRUE;
 271     UInt_t seqcharge=0,seqaverage=0,seqerror=0;
 272     UInt_t lastcharge=0,lastwas_falling=0;
 273     Int_t newbin=-1;
 274
 275
 276     if(fDeconvTime){
 277       redo: //This is a goto.
 278
 279       if(newbin > -1){
 280         //bin = newbin;
 281         newbin = -1;
 282       }
 283
 284       lastcharge=0;
 285       lastwas_falling = 0;
 286     }
 287
 288     while(1){ //Loop over time bins of current pad
 289       // read all the values for one pad at once to calculate the base line
 290       if (pCurrentPad) {
 291         if (!pCurrentPad->IsStarted()) {
 292           //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
 293           pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
 294           if ((pCurrentPad->StartEvent())>=0) {
 295             do {
 296               if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
 297               if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
 298               pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
 299               //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
 300             } while ((readValue = fDigitReader->Next())!=0);
 301           }
 302           pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
 303           if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
 304             //HLTDebug("no data available after zero suppression");
 305             pCurrentPad->StopEvent();
 306             pCurrentPad->ResetHistory();
 307             break;
 308           }
 309           time=pCurrentPad->GetCurrentPosition();
 310           if (time>pCurrentPad->GetSize()) {
 311             HLTError("invalid time bin for pad");
 312             break;
 313           }
 314         }
 315       }
 316
 317       if (pCurrentPad) {
 318         charge = pCurrentPad->GetCorrectedData();
 319       } else {
 320       charge = fDigitReader->GetSignal();
 321       }
 322       //HLTDebug("get next charge value: position %d charge %d", time, charge);
 323
 324
 325       // CHARGE DEBUG
 326       if (fDigitReader->GetRow() == 90){
 327 /////     LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90")  << "PAD=" <<  fDigitReader->GetPad() << "  TIME=" <<  fDigitReader->GetTime()
 328           //                                       << "  SIGNAL=" <<  fDigitReader->GetSignal() << ENDLOG;
 329
 330       }
 331
 332       if(time >= AliHLTTPCTransform::GetNTimeBins()){
 333         HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
 334         break;
 335       }
 336
 337
 338       //Get the current ADC-value
 339       if(fDeconvTime){
 340
 341         //Check if the last pixel in the sequence is smaller than this
 342         if(charge > lastcharge){
 343           if(lastwas_falling){
 344             newbin = 1;
 345             break;
 346           }
 347         }
 348         else lastwas_falling = 1; //last pixel was larger than this
 349         lastcharge = charge;
 350       }
 351
 352       //Sum the total charge of this sequence
 353       seqcharge += charge;
 354       seqaverage += time*charge;
 355       seqerror += time*time*charge;
 356
 357       if (pCurrentPad) {
 358
 359         if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
 360           pCurrentPad->StopEvent();
 361           pCurrentPad->ResetHistory();
 362           if(readValue) {
 363             newPad = fDigitReader->GetPad();
 364             newTime = fDigitReader->GetTime();
 365             newRow = fDigitReader->GetRow() + rowOffset;
 366           }
 367           break;
 368         }
 369
 370         newPad=pCurrentPad->GetPadNumber();
 371         newTime=pCurrentPad->GetCurrentPosition();
 372         newRow=pCurrentPad->GetRowNumber();
 373       } else {
 374       readValue = fDigitReader->Next();
 375       //Check where to stop:
 376       if(!readValue) break; //No more value
 377
 378       newPad = fDigitReader->GetPad();
 379       newTime = fDigitReader->GetTime();
 380       newRow = fDigitReader->GetRow() + rowOffset;
 381       }
 382
 383       if(newPad != pad)break; //new pad
 384       if(newTime != time+1) break; //end of sequence
 385
 386       // pad = newpad;    is equal
 387       time = newTime;
 388
 389     }//end loop over sequence
 390
 391     //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
 392     //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
 393     if (seqcharge<=0) {
 394       // with active zero suppression zero values are possible
 395       continue;
 396     }
 397
 398     //Calculate mean of sequence:
 399     Int_t seqmean=0;
 400     if(seqcharge)
 401       seqmean = seqaverage/seqcharge;
 402     else{
 403       LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
 404         <<"Error in data given to the cluster finder"<<ENDLOG;
 405       seqmean = 1;
 406       seqcharge = 1;
 407     }
 408
 409     //Calculate mean in pad direction:
 410     Int_t padmean = seqcharge*pad;
 411     Int_t paderror = pad*padmean;
 412
 413
 414     //Compare with results on previous pad:
 415     for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
 416
 417       //dont merge sequences on the same pad twice
 418       if(previousPt[p]->fLastMergedPad==pad) continue;
 419
 420       Int_t difference = seqmean - previousPt[p]->fMean;
 421       if(difference < -fMatch) break;
 422
 423       if(difference <= fMatch){ //There is a match here!!
 424         AliClusterData *local = previousPt[p];
 425
 426         if(fDeconvPad){
 427           if(seqcharge > local->fLastCharge){
 428             if(local->fChargeFalling){ //The previous pad was falling
 429               break; //create a new cluster
 430             }
 431           }
 432           else local->fChargeFalling = 1;
 433           local->fLastCharge = seqcharge;
 434         }
 435
 436         //Don't create a new cluster, because we found a match
 437         newcluster = kFALSE;
 438
 439         //Update cluster on current pad with the matching one:
 440         local->fTotalCharge += seqcharge;
 441         local->fPad += padmean;
 442         local->fPad2 += paderror;
 443         local->fTime += seqaverage;
 444         local->fTime2 += seqerror;
 445         local->fMean = seqmean;
 446         local->fFlags++; //means we have more than one pad
 447         local->fLastMergedPad = pad;
 448
 449         currentPt[ncurrent] = local;
 450         ncurrent++;
 451
 452         break;
 453       } //Checking for match at previous pad
 454     } //Loop over results on previous pad.
 455
 456
 457     if(newcluster && ncurrent<kPadArraySize){
 458       //Start a new cluster. Add it to the clusterlist, and update
 459       //the list of pointers to clusters in current pad.
 460       //current pad will be previous pad on next pad.
 461
 462       //Add to the clusterlist:
 463       AliClusterData *tmp = &clusterlist[ntotal];
 464       tmp->fTotalCharge = seqcharge;
 465       tmp->fPad = padmean;
 466       tmp->fPad2 = paderror;
 467       tmp->fTime = seqaverage;
 468       tmp->fTime2 = seqerror;
 469       tmp->fMean = seqmean;
 470       tmp->fFlags = 0;  //flags for single pad clusters
 471       tmp->fLastMergedPad = pad;
 472
 473       if(fDeconvPad){
 474         tmp->fChargeFalling = 0;
 475         tmp->fLastCharge = seqcharge;
 476       }
 477
 478       //Update list of pointers to previous pad:
 479       currentPt[ncurrent] = &clusterlist[ntotal];
 480       ntotal++;
 481       ncurrent++;
 482     }
 483
 484     if(fDeconvTime)
 485       if(newbin >= 0) goto redo;
 486
 487     // to prevent endless loop
 488     if(time >= AliHLTTPCTransform::GetNTimeBins()){
 489       HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
 490       break;
 491     }
 492
 493
 494     if(!readValue) break; //No more value
 495
 496     if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
 497       HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
 498       break;
 499     }
 500
 501     if(fCurrentRow != newRow){
 502       WriteClusters(ntotal,clusterlist);
 503
 504       lastpad = 123456789;
 505
 506       currentPt = pad2;
 507       previousPt = pad1;
 508       nprevious=0;
 509       ncurrent=0;
 510       ntotal=0;
 511
 512       fCurrentRow = newRow;
 513     }
 514
 515     pad = newPad;
 516     time = newTime;
 517
 518   } // END while(readValue)
 519
 520   WriteClusters(ntotal,clusterlist);
 521
 522   HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
 523
 524 } // ENDEND
 525
 526 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list)
 527 {
 528   //write cluster to output pointer
 529   Int_t thisrow,thissector;
 530   UInt_t counter = fNClusters;
 531
 532   for(int j=0; j<nclusters; j++)
 533     {
 534       if(!list[j].fFlags) continue; //discard single pad clusters
 535       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
 536
 537       Float_t xyz[3];
 538       Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
 539       Float_t fpad2=fXYErr*fXYErr; //fixed given error
 540       Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
 541       Float_t ftime2=fZErr*fZErr;  //fixed given error
 542
 543
 544
 545
 546
 547
 548       if(fCalcerr) { //calc the errors, otherwice take the fixed error
 549         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
 550         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
 551         Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
 552         sy2/=q2;
 553         if(sy2 < 0) {
 554             LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
 555               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
 556             continue;
 557         } else {
 558           if(!fRawSP){
 559             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
 560             if(sy2 != 0){
 561               fpad2*=0.108; //constants are from offline studies
 562               if(patch<2)
 563                 fpad2*=2.07;
 564             }
 565           } else fpad2=sy2; //take the width not the error
 566         }
 567         Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
 568         sz2/=q2;
 569         if(sz2 < 0){
 570           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
 571             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
 572           continue;
 573         } else {
 574           if(!fRawSP){
 575             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
 576             if(sz2 != 0) {
 577               ftime2 *= 0.169; //constants are from offline studies
 578               if(patch<2)
 579                 ftime2 *= 1.77;
 580             }
 581           } else ftime2=sz2; //take the width, not the error
 582         }
 583       }
 584       if(fStdout==kTRUE)
 585         cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
 586
 587       if(!fRawSP){
 588         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
 589         AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
 590
 591         if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
 592           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
 593         if(fNClusters >= fMaxNClusters)
 594           {
 595             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
 596               <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
 597             return;
 598           }
 599
 600         fSpacePointData[counter].fX = xyz[0];
 601         fSpacePointData[counter].fY = xyz[1];
 602         fSpacePointData[counter].fZ = xyz[2];
 603
 604       } else {
 605         fSpacePointData[counter].fX = fCurrentRow;
 606         fSpacePointData[counter].fY = fpad;
 607         fSpacePointData[counter].fZ = ftime;
 608       }
 609
 610       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
 611       fSpacePointData[counter].fPadRow = fCurrentRow;
 612       fSpacePointData[counter].fSigmaY2 = fpad2;
 613       fSpacePointData[counter].fSigmaZ2  = ftime2;
 614
 615       fSpacePointData[counter].fUsed = kFALSE;         // only used / set in AliHLTTPCDisplay
 616       fSpacePointData[counter].fTrackN = -1;           // only used / set in AliHLTTPCDisplay
 617
 618       Int_t patch=fCurrentPatch;
 619       if(patch==-1) patch=0; //never store negative patch number
 620       fSpacePointData[counter].fID = counter
 621         +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
 622
 623 #ifdef do_mc
 624       Int_t trackID[3];
 625       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
 626
 627       fSpacePointData[counter].fTrackID[0] = trackID[0];
 628       fSpacePointData[counter].fTrackID[1] = trackID[1];
 629       fSpacePointData[counter].fTrackID[2] = trackID[2];
 630
 631       //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
 632 #endif
 633
 634       fNClusters++;
 635       counter++;
 636     }
 637 }
 638
 639 // STILL TO FIX  ----------------------------------------------------------------------------
 640
 641 #ifdef do_mc
 642 void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
 643 {
 644   //get mc id
 645   AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
 646
 647   trackID[0]=trackID[1]=trackID[2]=-2;
 648   //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
 649   for(Int_t i=fFirstRow; i<=fLastRow; i++){
 650     if(rowPt->fRow < (UInt_t)fCurrentRow){
 651       AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
 652       continue;
 653     }
 654     AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
 655     for(UInt_t j=0; j<rowPt->fNDigit; j++){
 656       Int_t cpad = digPt[j].fPad;
 657       Int_t ctime = digPt[j].fTime;
 658       if(cpad != pad) continue;
 659       if(ctime != time) continue;
 660
 661       trackID[0] = digPt[j].fTrackID[0];
 662       trackID[1] = digPt[j].fTrackID[1];
 663       trackID[2] = digPt[j].fTrackID[2];
 664
 665       //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
 666       break;
 667     }
 668     break;
 669   }
 670 }
 671 #endif