HLT/TPCLib/AliHLTTPCClusterFinder.cxx

   1 // @(#) $Id$
   2 // Original: AliHLTClustFinderNew.cxx,v 1.29 2005/06/14 10:55:21 cvetan Exp
   3
   4 //**************************************************************************
   5 //* This file is property of and copyright by the ALICE HLT Project        *
   6 //* ALICE Experiment at CERN, All rights reserved.                         *
   7 //*                                                                        *
   8 //* Primary Authors: Anders Vestbo, Constantin Loizides                    *
   9 //* Developers:      Kenneth Aamodt <kenneth.aamodt@student.uib.no>        *
  10 //*                  Kalliopi Kanaki                                       *
  11 //*                  for The ALICE HLT Project.                            *
  12 //*                                                                        *
  13 //* Permission to use, copy, modify and distribute this software and its   *
  14 //* documentation strictly for non-commercial purposes is hereby granted   *
  15 //* without fee, provided that the above copyright notice appears in all   *
  16 //* copies and that both the copyright notice and this permission notice   *
  17 //* appear in the supporting documentation. The authors make no claims     *
  18 //* about the suitability of this software for any purpose. It is          *
  19 //* provided "as is" without express or implied warranty.                  *
  20 //**************************************************************************
  21
  22 /** @file   AliHLTTPCClusterFinder.cxx
  23     @author Kenneth Aamodt, Kalliopi Kanaki
  24     @date
  25     @brief  Cluster Finder for the TPC
  26 */
  27
  28 #include "AliHLTTPCDigitReader.h"
  29 #include "AliHLTTPCRootTypes.h"
  30 #include "AliHLTTPCLogging.h"
  31 #include "AliHLTTPCClusterFinder.h"
  32 #include "AliHLTTPCSpacePointData.h"
  33 #include "AliHLTTPCMemHandler.h"
  34 #include "AliHLTTPCPad.h"
  35 #include <sys/time.h>
  36 #include <algorithm>
  37
  38 #if __GNUC__ >= 3
  39 using namespace std;
  40 #endif
  41
  42 ClassImp(AliHLTTPCClusterFinder)
  43
  44 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
  45   :
  46   fClustersHWAddressVector(),
  47   fRowPadVector(),
  48   fSpacePointData(NULL),
  49   fDigitReader(NULL),
  50   fPtr(NULL),
  51   fSize(0),
  52   fDeconvTime(kFALSE),
  53   fDeconvPad(kFALSE),
  54   fStdout(kFALSE),
  55   fCalcerr(kTRUE),
  56   fRawSP(kFALSE),
  57   fFirstRow(0),
  58   fLastRow(0),
  59   fCurrentRow(0),
  60   fCurrentSlice(0),
  61   fCurrentPatch(0),
  62   fMatch(1),
  63   fThreshold(10),
  64   fNClusters(0),
  65   fMaxNClusters(0),
  66   fXYErr(0.2),
  67   fZErr(0.3),
  68   fOccupancyLimit(1.0),
  69   fUnsorted(0),
  70   fVectorInitialized(kFALSE),
  71   fClusters(),
  72   fClustersMCInfo(),
  73   fMCDigits(),
  74   fNumberOfPadsInRow(NULL),
  75   fNumberOfRows(0),
  76   fRowOfFirstCandidate(0),
  77   fDoPadSelection(kFALSE),
  78   fFirstTimeBin(0),
  79   fLastTimeBin(AliHLTTPCTransform::GetNTimeBins()),
  80   fTotalChargeOfPreviousClusterCandidate(0),
  81   fChargeOfCandidatesFalling(kFALSE),
  82   f32BitFormat(kFALSE),
  83   fDoMC(kFALSE),
  84   fClusterMCVector()
  85 {
  86   //constructor
  87 }
  88
  89 AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder(){
  90   // see header file for class documentation
  91
  92   //destructor
  93   if(fVectorInitialized){
  94     DeInitializePadArray();
  95   }
  96   if(fNumberOfPadsInRow){
  97     delete [] fNumberOfPadsInRow;
  98     fNumberOfPadsInRow=NULL;
  99   }
 100 }
 101
 102 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints){
 103   // see header file for class documentation
 104
 105   //init slice
 106   fNClusters = 0;
 107   fMaxNClusters = nmaxpoints;
 108   fCurrentSlice = slice;
 109   fCurrentPatch = patch;
 110   fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
 111   fLastRow=AliHLTTPCTransform::GetLastRow(patch);
 112 }
 113
 114 void AliHLTTPCClusterFinder::InitializePadArray(){
 115   // see header file for class documentation
 116
 117   if(fCurrentPatch>5||fCurrentPatch<0){
 118     HLTFatal("Patch is not set");
 119     return;
 120   }
 121
 122   HLTDebug("Patch number=%d",fCurrentPatch);
 123
 124   fFirstRow = AliHLTTPCTransform::GetFirstRow(fCurrentPatch);
 125   fLastRow = AliHLTTPCTransform::GetLastRow(fCurrentPatch);
 126
 127   fNumberOfRows=fLastRow-fFirstRow+1;
 128   fNumberOfPadsInRow= new UInt_t[fNumberOfRows];
 129
 130   memset( fNumberOfPadsInRow, 0, sizeof(Int_t)*(fNumberOfRows));
 131
 132   for(UInt_t i=0;i<fNumberOfRows;i++){
 133     fNumberOfPadsInRow[i]=AliHLTTPCTransform::GetNPads(i+fFirstRow);
 134     AliHLTTPCPadVector tmpRow;
 135     for(UInt_t j=0;j<=fNumberOfPadsInRow[i];j++){
 136       AliHLTTPCPad *tmpPad = new AliHLTTPCPad(2);
 137       tmpPad->SetID(i,j);
 138       tmpRow.push_back(tmpPad);
 139     }
 140     fRowPadVector.push_back(tmpRow);
 141   }
 142   fVectorInitialized=kTRUE;
 143 }
 144
 145 Int_t AliHLTTPCClusterFinder::DeInitializePadArray(){
 146   // see header file for class documentation
 147
 148   for(UInt_t i=0;i<fNumberOfRows;i++){
 149     for(UInt_t j=0;j<=fNumberOfPadsInRow[i];j++){
 150       delete fRowPadVector[i][j];
 151       fRowPadVector[i][j]=NULL;
 152     }
 153     fRowPadVector[i].clear();
 154   }
 155   fRowPadVector.clear();
 156   return 1;
 157 }
 158
 159
 160
 161 void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt){
 162   // see header file for class documentation
 163   //set pointer to output
 164   fSpacePointData = pt;
 165 }
 166
 167
 168 void AliHLTTPCClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size){
 169   // see header file for class documentation
 170   //set input pointer
 171   fPtr = (UChar_t*)ptr;
 172   fSize = size;
 173
 174   if(!fVectorInitialized){
 175     InitializePadArray();
 176   }
 177
 178   if (fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice)<0) {
 179     HLTError("failed setting up digit reader (InitBlock)");
 180     return;
 181   }
 182
 183   while(fDigitReader->NextChannel()){
 184     UInt_t row=fDigitReader->GetRow();
 185     UInt_t pad=fDigitReader->GetPad();
 186
 187     if(row>=fRowPadVector.size()){
 188       HLTError("Row number is to large: %d, max is %d",row,fRowPadVector.size()-1);
 189       continue;
 190     }
 191     if(pad>=fRowPadVector[row].size()){
 192       HLTError("Pad number is to large: %d, max is %d",pad,fRowPadVector[row].size());
 193       continue;
 194     }
 195
 196     while(fDigitReader->NextBunch()){
 197       if(fDigitReader->GetBunchSize()>1){//to remove single timebin values, this will have to change at some point
 198         UInt_t time = fDigitReader->GetTime();
 199         if((Int_t)time>=fFirstTimeBin && (Int_t)time+fDigitReader->GetBunchSize()<=fLastTimeBin){
 200           // Kenneth: 20-04-09. The following if have been added because of inconsistency in the 40 bit decoder and the 32 bit decoder.
 201           // GetSignals() in the 40 bit decoder returns an array of UInt_t while the 32 bit one returns UShort_t
 202           // The same is true for the function ReadDataUnsortedDeconvoluteTime() below.
 203           // In addition the signals are organized in the opposite direction
 204           if(f32BitFormat){
 205             const UShort_t *bunchData= fDigitReader->GetSignalsShort();
 206             AliHLTTPCClusters candidate;
 207             for(Int_t i=fDigitReader->GetBunchSize()-1;i>=0;i--){
 208               candidate.fTotalCharge+=bunchData[i];
 209               candidate.fTime += time*bunchData[i];
 210               candidate.fTime2 += time*time*bunchData[i];
 211               if(bunchData[i]>candidate.fQMax){
 212                 candidate.fQMax=bunchData[i];
 213               }
 214               time++;
 215             }
 216             if(candidate.fTotalCharge>0){
 217               candidate.fMean=candidate.fTime/candidate.fTotalCharge;
 218               candidate.fPad=candidate.fTotalCharge*pad;
 219               candidate.fPad2=candidate.fPad*pad;
 220               candidate.fLastMergedPad=pad;
 221               candidate.fRowNumber=row+fDigitReader->GetRowOffset();
 222             }
 223             if(fRowPadVector[row][pad] != NULL){
 224               fRowPadVector[row][pad]->AddClusterCandidate(candidate);
 225             }
 226           }
 227           else{
 228             const UInt_t *bunchData= fDigitReader->GetSignals();
 229             AliHLTTPCClusters candidate;
 230             if(fDoMC){
 231               const AliHLTTPCDigitData* digits = fDigitReader->GetBunchDigits();
 232               for(Int_t i=0;i<fDigitReader->GetBunchSize();i++){
 233                 candidate.fTotalCharge+=bunchData[i];
 234                 candidate.fTime += time*bunchData[i];
 235                 candidate.fTime2 += time*time*bunchData[i];
 236                 if(bunchData[i]>candidate.fQMax){
 237                   candidate.fQMax=bunchData[i];
 238                 }
 239                 fMCDigits.push_back(digits[i]);
 240                 time++;
 241               }
 242             }
 243             else{
 244               for(Int_t i=0;i<fDigitReader->GetBunchSize();i++){
 245                 candidate.fTotalCharge+=bunchData[i];
 246                 candidate.fTime += time*bunchData[i];
 247                 candidate.fTime2 += time*time*bunchData[i];
 248                 if(bunchData[i]>candidate.fQMax){
 249                   candidate.fQMax=bunchData[i];
 250                 }
 251                 time++;
 252               }
 253             }
 254             if(candidate.fTotalCharge>0){
 255               candidate.fMean=candidate.fTime/candidate.fTotalCharge;
 256               candidate.fPad=candidate.fTotalCharge*pad;
 257               candidate.fPad2=candidate.fPad*pad;
 258               candidate.fLastMergedPad=pad;
 259               candidate.fRowNumber=row+fDigitReader->GetRowOffset();
 260             }
 261             if(fRowPadVector[row][pad] != NULL){
 262               fRowPadVector[row][pad]->AddClusterCandidate(candidate);
 263               if(fDoMC){
 264                 fRowPadVector[row][pad]->AddCandidateDigits(fMCDigits);
 265                 fMCDigits.clear();
 266               }
 267             }
 268           }
 269         }
 270       }
 271     }
 272   }
 273 }
 274
 275 void AliHLTTPCClusterFinder::ReadDataUnsortedDeconvoluteTime(void* ptr,unsigned long size){
 276   // see header file for class documentation
 277
 278   //set input pointer
 279   fPtr = (UChar_t*)ptr;
 280   fSize = size;
 281
 282   if(!fVectorInitialized){
 283     InitializePadArray();
 284   }
 285
 286   if (fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice)<0) {
 287     HLTError("failed setting up digit reader (InitBlock)");
 288     return;
 289   }
 290
 291   while(fDigitReader->NextChannel()){
 292     UInt_t row=fDigitReader->GetRow();
 293     UInt_t pad=fDigitReader->GetPad();
 294
 295     while(fDigitReader->NextBunch()){
 296       if(fDigitReader->GetBunchSize()>1){//to remove single timebin values, this will have to change at some point
 297         UInt_t time = fDigitReader->GetTime();
 298         if((Int_t)time>=fFirstTimeBin && (Int_t)time+fDigitReader->GetBunchSize()<=fLastTimeBin){
 299           Int_t indexInBunchData=0;
 300           Bool_t moreDataInBunch=kFALSE;
 301           UInt_t prevSignal=0;
 302           Bool_t signalFalling=kFALSE;
 303
 304           // Kenneth: 20-04-09. The following if have been added because of inconsistency in the 40 bit decoder and the 32 bit decoder.
 305           // GetSignals() in the 40 bit decoder returns an array of UInt_t while the 32 bit one returns UShort_t
 306           // The same is true for the function ReadDataUnsorted() above.
 307           // In addition the signals are organized in the opposite direction
 308           if(f32BitFormat){
 309             indexInBunchData = fDigitReader->GetBunchSize();
 310             const UShort_t *bunchData= fDigitReader->GetSignalsShort();
 311             do{
 312               AliHLTTPCClusters candidate;
 313               //for(Int_t i=indexInBunchData;i<fDigitReader->GetBunchSize();i++){
 314               for(Int_t i=indexInBunchData;i>=0;i--){
 315                 // Checks if one need to deconvolute the signals
 316                 if(bunchData[i]>prevSignal && signalFalling==kTRUE){
 317                   if(i<fDigitReader->GetBunchSize()-1){ // means there are more than one signal left in the bunch
 318                     moreDataInBunch=kTRUE;
 319                     prevSignal=0;
 320                   }
 321                   break;
 322                 }
 323
 324                 // Checks if the signal is 0, then quit processing the data.
 325                 if(bunchData[i]==0 && i<fDigitReader->GetBunchSize()-1){//means we have 0 data fom the rcu, might happen depending on the configuration settings
 326                   moreDataInBunch=kTRUE;
 327                   prevSignal=0;
 328                   break;
 329                 }
 330
 331                 if(prevSignal>bunchData[i]){//means the peak of the signal has been reached and deconvolution will happen if the signal rise again.
 332                   signalFalling=kTRUE;
 333                 }
 334                 candidate.fTotalCharge+=bunchData[i];
 335                 candidate.fTime += time*bunchData[i];
 336                 candidate.fTime2 += time*time*bunchData[i];
 337                 if(bunchData[i]>candidate.fQMax){
 338                   candidate.fQMax=bunchData[i];
 339                 }
 340
 341                 prevSignal=bunchData[i];
 342                 time++;
 343                 indexInBunchData--;
 344               }
 345               if(candidate.fTotalCharge>0){
 346                 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
 347                 candidate.fPad=candidate.fTotalCharge*pad;
 348                 candidate.fPad2=candidate.fPad*pad;
 349                 candidate.fLastMergedPad=pad;
 350                 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
 351               }
 352               fRowPadVector[row][pad]->AddClusterCandidate(candidate);
 353               fRowPadVector[row][pad]->AddCandidateDigits(fMCDigits);
 354               if(indexInBunchData<fDigitReader->GetBunchSize()-1){
 355                 moreDataInBunch=kFALSE;
 356               }
 357             }while(moreDataInBunch);
 358           }
 359           else{
 360             const UInt_t *bunchData= fDigitReader->GetSignals();
 361             do{
 362               AliHLTTPCClusters candidate;
 363               for(Int_t i=indexInBunchData;i<fDigitReader->GetBunchSize();i++){
 364                 // Checks if one need to deconvolute the signals
 365                 if(bunchData[i]>prevSignal && signalFalling==kTRUE){
 366                   if(i<fDigitReader->GetBunchSize()-1){ // means there are more than one signal left in the bunch
 367                     moreDataInBunch=kTRUE;
 368                     prevSignal=0;
 369                   }
 370                   break;
 371                 }
 372
 373                 // Checks if the signal is 0, then quit processing the data.
 374                 if(bunchData[i]==0 && i<fDigitReader->GetBunchSize()-1){//means we have 0 data fom the rcu, might happen depending on the configuration settings
 375                   moreDataInBunch=kTRUE;
 376                   prevSignal=0;
 377                   break;
 378                 }
 379
 380                 if(prevSignal>bunchData[i]){//means the peak of the signal has been reached and deconvolution will happen if the signal rise again.
 381                   signalFalling=kTRUE;
 382                 }
 383                 candidate.fTotalCharge+=bunchData[i];
 384                 candidate.fTime += time*bunchData[i];
 385                 candidate.fTime2 += time*time*bunchData[i];
 386                 if(bunchData[i]>candidate.fQMax){
 387                   candidate.fQMax=bunchData[i];
 388                 }
 389
 390                 prevSignal=bunchData[i];
 391                 time++;
 392                 indexInBunchData++;
 393               }
 394               if(candidate.fTotalCharge>0){
 395                 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
 396                 candidate.fPad=candidate.fTotalCharge*pad;
 397                 candidate.fPad2=candidate.fPad*pad;
 398                 candidate.fLastMergedPad=pad;
 399                 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
 400               }
 401               fRowPadVector[row][pad]->AddClusterCandidate(candidate);
 402               if(indexInBunchData<fDigitReader->GetBunchSize()-1){
 403                 moreDataInBunch=kFALSE;
 404               }
 405             }while(moreDataInBunch);
 406           }
 407         }
 408       }
 409     }
 410   }
 411 }
 412
 413 Bool_t AliHLTTPCClusterFinder::ComparePads(AliHLTTPCPad *nextPad,AliHLTTPCClusters* cluster,Int_t nextPadToRead){
 414   // see header file for class documentation
 415
 416   //Checking if we have a match on the next pad
 417   for(UInt_t candidateNumber=0;candidateNumber<nextPad->fClusterCandidates.size();candidateNumber++){
 418     AliHLTTPCClusters *candidate =&nextPad->fClusterCandidates[candidateNumber];
 419     if(cluster->fMean-candidate->fMean==1 || candidate->fMean-cluster->fMean==1 || cluster->fMean-candidate->fMean==0){
 420       if(fDeconvPad){
 421         if(candidate->fTotalCharge<fTotalChargeOfPreviousClusterCandidate){//peak is reached
 422           fChargeOfCandidatesFalling=kTRUE;
 423         }
 424         if(candidate->fTotalCharge>fTotalChargeOfPreviousClusterCandidate && fChargeOfCandidatesFalling==kTRUE){//we have deconvolution
 425           return kFALSE;
 426         }
 427       }
 428       cluster->fMean=candidate->fMean;
 429       cluster->fTotalCharge+=candidate->fTotalCharge;
 430       cluster->fTime += candidate->fTime;
 431       cluster->fTime2 += candidate->fTime2;
 432       cluster->fPad+=candidate->fPad;
 433       cluster->fPad2+=candidate->fPad2;
 434       cluster->fLastMergedPad=candidate->fPad;
 435       if(candidate->fQMax>cluster->fQMax){
 436         cluster->fQMax=candidate->fQMax;
 437       }
 438       if(fDoMC){
 439         FillMCClusterVector(nextPad->GetCandidateDigits(candidateNumber));
 440       }
 441
 442       if(fDoPadSelection){
 443         UInt_t rowNo = nextPad->GetRowNumber();
 444         UInt_t padNo = nextPad->GetPadNumber();
 445         if(padNo-1>0){
 446           fRowPadVector[rowNo][padNo-2]->fSelectedPad=kTRUE;
 447           fRowPadVector[rowNo][padNo-2]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo-2);
 448         }
 449         fRowPadVector[rowNo][padNo-1]->fSelectedPad=kTRUE;// quick solution to set the first pad to selected
 450         fRowPadVector[rowNo][padNo-1]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo-1);
 451         fRowPadVector[rowNo][padNo]->fSelectedPad=kTRUE;
 452         fRowPadVector[rowNo][padNo]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo);
 453       }
 454
 455       //setting the matched pad to used
 456       nextPad->fUsedClusterCandidates[candidateNumber]=1;
 457       nextPadToRead++;
 458       if(nextPadToRead<(Int_t)fNumberOfPadsInRow[fRowOfFirstCandidate]){
 459         nextPad=fRowPadVector[fRowOfFirstCandidate][nextPadToRead];
 460         ComparePads(nextPad,cluster,nextPadToRead);
 461       }
 462       else{
 463         return kFALSE;
 464       }
 465     }
 466     else{
 467       return kFALSE;
 468     }
 469   }
 470   return kFALSE;
 471 }
 472
 473 Int_t AliHLTTPCClusterFinder::FillHWAddressList(AliHLTUInt16_t *hwaddlist, Int_t maxHWadd){
 474   // see header file for class documentation
 475
 476   Int_t counter=0;
 477   for(UInt_t row=0;row<fNumberOfRows;row++){
 478     for(UInt_t pad=0;pad<fNumberOfPadsInRow[row]-1;pad++){
 479       if(fRowPadVector[row][pad]->fSelectedPad){
 480        if(counter<maxHWadd){
 481          hwaddlist[counter]=(AliHLTUInt16_t)fRowPadVector[row][pad]->fHWAddress;
 482          counter++;
 483        }
 484        else{
 485          HLTWarning("To many hardwareaddresses, skip adding");
 486        }
 487
 488       }
 489     }
 490   }
 491   return counter;
 492 }
 493
 494
 495 Int_t AliHLTTPCClusterFinder::FillOutputMCInfo(AliHLTTPCClusterFinder::ClusterMCInfo * outputMCInfo, Int_t maxNumberOfClusterMCInfo){
 496   // see header file for class documentation
 497
 498   Int_t counter=0;
 499
 500   for(UInt_t mc=0;mc<fClustersMCInfo.size();mc++){
 501     if(counter<maxNumberOfClusterMCInfo){
 502       outputMCInfo[counter] = fClustersMCInfo[mc];
 503       counter++;
 504     }
 505     else{
 506       HLTWarning("To much MCInfo has been added (no more space), skip adding");
 507     }
 508   }
 509   return counter;
 510 }
 511
 512 void AliHLTTPCClusterFinder::FindClusters(){
 513   // see header file for function documentation
 514
 515   AliHLTTPCClusters* tmpCandidate=NULL;
 516   for(UInt_t row=0;row<fNumberOfRows;row++){
 517     fRowOfFirstCandidate=row;
 518     for(UInt_t pad=0;pad<fNumberOfPadsInRow[row];pad++){
 519       AliHLTTPCPad *tmpPad=fRowPadVector[row][pad];
 520       for(size_t candidate=0;candidate<tmpPad->fClusterCandidates.size();candidate++){
 521         if(tmpPad->fUsedClusterCandidates[candidate]){
 522           continue;
 523         }
 524         tmpCandidate=&tmpPad->fClusterCandidates[candidate];
 525         UInt_t tmpTotalCharge=tmpCandidate->fTotalCharge;
 526
 527         if(fDoMC){
 528           fClusterMCVector.clear();
 529           FillMCClusterVector(tmpPad->GetCandidateDigits(candidate));
 530         }
 531
 532         ComparePads(fRowPadVector[row][pad+1],tmpCandidate,pad+1);
 533         if(tmpCandidate->fTotalCharge>tmpTotalCharge){
 534           //we have a cluster
 535           fClusters.push_back(*tmpCandidate);
 536           if(fDoMC){
 537             //sort the vector (large->small) according to weight and remove elements above 2 (keep 0 1 and 2)
 538             //sort(fClusterMCVector,fClusterMCVector.size(), MCWeight::CompareWeights );
 539             ClusterMCInfo tmpClusterMCInfo;
 540
 541             MCWeight zeroMC;
 542             zeroMC.fMCID=-1;
 543             zeroMC.fWeight=0;
 544
 545             if(fClusterMCVector.size()>0){
 546               tmpClusterMCInfo.fClusterID[0]=fClusterMCVector.at(0);
 547             }
 548             else{
 549               tmpClusterMCInfo.fClusterID[0]=zeroMC;
 550             }
 551
 552             if(fClusterMCVector.size()>1){
 553             tmpClusterMCInfo.fClusterID[1]=fClusterMCVector.at(1);
 554             }
 555             else{
 556               tmpClusterMCInfo.fClusterID[1]=zeroMC;
 557             }
 558
 559             if(fClusterMCVector.size()>2){
 560             tmpClusterMCInfo.fClusterID[2]=fClusterMCVector.at(2);
 561             }
 562             else{
 563               tmpClusterMCInfo.fClusterID[2]=zeroMC;
 564             }
 565
 566             fClustersMCInfo.push_back(tmpClusterMCInfo);
 567           }
 568
 569         }
 570       }
 571       tmpPad->ClearCandidates();
 572     }
 573     fRowPadVector[row][fNumberOfPadsInRow[row]]->ClearCandidates();
 574   }
 575
 576   HLTInfo("Found %d clusters.",fClusters.size());
 577
 578   //TODO:  Change so it stores AliHLTTPCSpacePointData directly, instead of this copying
 579
 580   AliClusterData * clusterlist = new AliClusterData[fClusters.size()]; //Clusterlist
 581   for(unsigned int i=0;i<fClusters.size();i++){
 582     clusterlist[i].fTotalCharge = fClusters[i].fTotalCharge;
 583     clusterlist[i].fPad = fClusters[i].fPad;
 584     clusterlist[i].fPad2 = fClusters[i].fPad2;
 585     clusterlist[i].fTime = fClusters[i].fTime;
 586     clusterlist[i].fTime2 = fClusters[i].fTime2;
 587     clusterlist[i].fMean = fClusters[i].fMean;
 588     clusterlist[i].fFlags = fClusters[i].fFlags;
 589     clusterlist[i].fChargeFalling = fClusters[i].fChargeFalling;
 590     clusterlist[i].fLastCharge = fClusters[i].fLastCharge;
 591     clusterlist[i].fLastMergedPad = fClusters[i].fLastMergedPad;
 592     clusterlist[i].fRow = fClusters[i].fRowNumber;
 593     clusterlist[i].fQMax = fClusters[i].fQMax;
 594   }
 595
 596   WriteClusters(fClusters.size(),clusterlist);
 597   delete [] clusterlist;
 598   fClusters.clear();
 599 }
 600
 601
 602 //---------------------------------- Under this line the old sorted clusterfinder functions can be found --------------------------------
 603
 604
 605 void AliHLTTPCClusterFinder::PrintClusters(){
 606   // see header file for class documentation
 607
 608   for(size_t i=0;i<fClusters.size();i++){
 609     HLTInfo("Cluster number: %d",i);
 610     HLTInfo("Row: %d \t Pad: %d",fClusters[i].fRowNumber,fClusters[i].fPad/fClusters[i].fTotalCharge);
 611     HLTInfo("Total Charge:   %d",fClusters[i].fTotalCharge);
 612     HLTInfo("fPad:           %d",fClusters[i].fPad);
 613     HLTInfo("PadError:       %d",fClusters[i].fPad2);
 614     HLTInfo("TimeMean:       %d",fClusters[i].fTime/fClusters[i].fTotalCharge);
 615     HLTInfo("TimeError:      %d",fClusters[i].fTime2);
 616     HLTInfo("EndOfCluster:");
 617   }
 618 }
 619
 620 void AliHLTTPCClusterFinder::FillMCClusterVector(vector<AliHLTTPCDigitData> digitData){
 621
 622   for(UInt_t d=0;d<digitData.size();d++){
 623     Int_t nIDsInDigit = (digitData.at(d).fTrackID[0]>=0) + (digitData.at(d).fTrackID[1]>=0) + (digitData.at(d).fTrackID[2]>=0);
 624     for(Int_t id=0; id<3; id++){
 625       if(digitData.at(d).fTrackID[id]>=0){
 626         Bool_t matchFound = kFALSE;
 627         MCWeight mc;
 628         mc.fMCID = digitData.at(d).fTrackID[id];
 629         mc.fWeight = ((Float_t)digitData.at(d).fCharge)/nIDsInDigit;
 630         for(UInt_t i=0;i<fClusterMCVector.size();i++){
 631           if(mc.fMCID == fClusterMCVector.at(i).fMCID){
 632             fClusterMCVector.at(i).fWeight += mc.fWeight;
 633             matchFound = kTRUE;
 634           }
 635         }
 636         if(matchFound == kFALSE){
 637           fClusterMCVector.push_back(mc);
 638         }
 639       }
 640     }
 641   }
 642 }
 643
 644
 645 void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
 646   //set input pointer
 647   fPtr = (UChar_t*)ptr;
 648   fSize = size;
 649 }
 650
 651 void AliHLTTPCClusterFinder::ProcessDigits(){
 652   // see header file for class documentation
 653
 654   int iResult=0;
 655   bool readValue = true;
 656   Int_t newRow = 0;
 657   Int_t rowOffset = 0;
 658   UShort_t time=0,newTime=0;
 659   UInt_t pad=0,newPad=0;
 660   AliHLTTPCSignal_t charge=0;
 661
 662   fNClusters = 0;
 663
 664   // initialize block for reading packed data
 665   iResult=fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
 666   if (iResult<0) return;
 667
 668   readValue = fDigitReader->Next();
 669
 670   // Matthias 08.11.2006 the following return would cause termination without writing the
 671   // ClusterData and thus would block the component. I just want to have the commented line
 672   // here for information
 673   //if (!readValue)return;
 674
 675   pad = fDigitReader->GetPad();
 676   time = fDigitReader->GetTime();
 677   fCurrentRow = fDigitReader->GetRow();
 678
 679   if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
 680     rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
 681
 682   fCurrentRow += rowOffset;
 683
 684   UInt_t lastpad = 123456789;
 685   const UInt_t kPadArraySize=5000;
 686   const UInt_t kClusterListSize=10000;
 687   AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
 688   AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
 689   AliClusterData clusterlist[kClusterListSize]; //Clusterlist
 690
 691   AliClusterData **currentPt;  //List of pointers to the current pad
 692   AliClusterData **previousPt; //List of pointers to the previous pad
 693   currentPt = pad2;
 694   previousPt = pad1;
 695   UInt_t nprevious=0,ncurrent=0,ntotal=0;
 696
 697   /* quick implementation of baseline calculation and zero suppression
 698      open a pad object for each pad and delete it after processing.
 699      later a list of pad objects with base line history can be used
 700      The whole thing only works if we really get unprocessed raw data, if
 701      the data is already zero suppressed, there might be gaps in the time
 702      bins.
 703    */
 704   Int_t gatingGridOffset=50;
 705   if(fFirstTimeBin>0){
 706     gatingGridOffset=fFirstTimeBin;
 707   }
 708   AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
 709   // just to make later conversion to a list of objects easier
 710   AliHLTTPCPad* pCurrentPad=NULL;
 711   /*
 712     if (fSignalThreshold>=0) {
 713     pCurrentPad=&baseline;
 714     baseline.SetThreshold(fSignalThreshold);
 715   }
 716   */
 717   while ( readValue!=0 && iResult>=0){   // Reads through all digits in block
 718     iResult=0;
 719
 720     if(pad != lastpad){
 721       //This is a new pad
 722
 723       //Switch the lists:
 724       if(currentPt == pad2){
 725         currentPt = pad1;
 726         previousPt = pad2;
 727       }
 728       else {
 729         currentPt = pad2;
 730         previousPt = pad1;
 731       }
 732       nprevious = ncurrent;
 733       ncurrent = 0;
 734       if(pad != lastpad+1){
 735         //this happens if there is a pad with no signal.
 736         nprevious = ncurrent = 0;
 737       }
 738       lastpad = pad;
 739     }
 740
 741     Bool_t newcluster = kTRUE;
 742     UInt_t seqcharge=0,seqaverage=0,seqerror=0;
 743     AliHLTTPCSignal_t lastcharge=0;
 744     UInt_t bLastWasFalling=0;
 745     Int_t newbin=-1;
 746
 747
 748     if(fDeconvTime){
 749       redo: //This is a goto.
 750
 751       if(newbin > -1){
 752         //bin = newbin;
 753         newbin = -1;
 754       }
 755
 756       lastcharge=0;
 757       bLastWasFalling = 0;
 758     }
 759
 760     while(iResult>=0){ //Loop over time bins of current pad
 761       // read all the values for one pad at once to calculate the base line
 762       if (pCurrentPad) {
 763         if (!pCurrentPad->IsStarted()) {
 764           //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
 765           pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
 766           if ((pCurrentPad->StartEvent())>=0) {
 767             do {
 768               if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
 769               if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
 770               pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
 771               //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
 772             } while ((readValue = fDigitReader->Next())!=0);
 773           }
 774           pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
 775           if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
 776             //HLTDebug("no data available after zero suppression");
 777             pCurrentPad->StopEvent();
 778             pCurrentPad->ResetHistory();
 779             break;
 780           }
 781           time=pCurrentPad->GetCurrentPosition();
 782           if (time>pCurrentPad->GetSize()) {
 783             HLTError("invalid time bin for pad");
 784             break;
 785           }
 786         }
 787       }
 788       if (pCurrentPad) {
 789         Float_t occupancy=pCurrentPad->GetOccupancy();
 790         //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
 791         if ( occupancy < fOccupancyLimit ) {
 792           charge = pCurrentPad->GetCorrectedData();
 793         } else {
 794           charge = 0;
 795           //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
 796         }
 797       } else {
 798         charge = fDigitReader->GetSignal();
 799       }
 800       //HLTDebug("get next charge value: position %d charge %d", time, charge);
 801
 802
 803       // CHARGE DEBUG
 804       if (fDigitReader->GetRow() == 90){
 805 /////     LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90")  << "PAD=" <<  fDigitReader->GetPad() << "  TIME=" <<  fDigitReader->GetTime()
 806           //                                       << "  SIGNAL=" <<  fDigitReader->GetSignal() << ENDLOG;
 807
 808       }
 809
 810       if(time >= AliHLTTPCTransform::GetNTimeBins()){
 811         HLTWarning("Pad %d: Timebin (%d) out of range (%d)", pad, time, AliHLTTPCTransform::GetNTimeBins());
 812         iResult=-ERANGE;
 813       }
 814
 815
 816       //Get the current ADC-value
 817       if(fDeconvTime){
 818
 819         //Check if the last pixel in the sequence is smaller than this
 820         if(charge > lastcharge){
 821           if(bLastWasFalling){
 822             newbin = 1;
 823             break;
 824           }
 825         }
 826         else bLastWasFalling = 1; //last pixel was larger than this
 827         lastcharge = charge;
 828       }
 829
 830       //Sum the total charge of this sequence
 831       seqcharge += charge;
 832       seqaverage += time*charge;
 833       seqerror += time*time*charge;
 834
 835       if (pCurrentPad) {
 836
 837         if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
 838           pCurrentPad->StopEvent();
 839           pCurrentPad->ResetHistory();
 840           if(readValue) {
 841             newPad = fDigitReader->GetPad();
 842             newTime = fDigitReader->GetTime();
 843             newRow = fDigitReader->GetRow() + rowOffset;
 844           }
 845           break;
 846         }
 847
 848         newPad=pCurrentPad->GetPadNumber();
 849         newTime=pCurrentPad->GetCurrentPosition();
 850         newRow=pCurrentPad->GetRowNumber();
 851       } else {
 852       readValue = fDigitReader->Next();
 853       //Check where to stop:
 854       if(!readValue) break; //No more value
 855
 856       newPad = fDigitReader->GetPad();
 857       newTime = fDigitReader->GetTime();
 858       newRow = fDigitReader->GetRow() + rowOffset;
 859       }
 860
 861       if(newPad != pad)break; //new pad
 862       if(newTime != time+1) break; //end of sequence
 863       if(iResult<0) break;
 864
 865       // pad = newpad;    is equal
 866       time = newTime;
 867
 868     }//end loop over sequence
 869
 870     //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
 871     //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
 872     if (seqcharge<=0) {
 873       // with active zero suppression zero values are possible
 874       continue;
 875     }
 876
 877     //Calculate mean of sequence:
 878     Int_t seqmean=0;
 879     if(seqcharge)
 880       seqmean = seqaverage/seqcharge;
 881     else{
 882       LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
 883         <<"Error in data given to the cluster finder"<<ENDLOG;
 884       seqmean = 1;
 885       seqcharge = 1;
 886     }
 887
 888     //Calculate mean in pad direction:
 889     Int_t padmean = seqcharge*pad;
 890     Int_t paderror = pad*padmean;
 891
 892     //Compare with results on previous pad:
 893     for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
 894
 895       //dont merge sequences on the same pad twice
 896       if(previousPt[p]->fLastMergedPad==pad) continue;
 897
 898       Int_t difference = seqmean - previousPt[p]->fMean;
 899       if(difference < -fMatch) break;
 900
 901       if(difference <= fMatch){ //There is a match here!!
 902         AliClusterData *local = previousPt[p];
 903
 904         if(fDeconvPad){
 905           if(seqcharge > local->fLastCharge){
 906             if(local->fChargeFalling){ //The previous pad was falling
 907               break; //create a new cluster
 908             }
 909           }
 910           else local->fChargeFalling = 1;
 911           local->fLastCharge = seqcharge;
 912         }
 913
 914         //Don't create a new cluster, because we found a match
 915         newcluster = kFALSE;
 916
 917         //Update cluster on current pad with the matching one:
 918         local->fTotalCharge += seqcharge;
 919         local->fPad += padmean;
 920         local->fPad2 += paderror;
 921         local->fTime += seqaverage;
 922         local->fTime2 += seqerror;
 923         local->fMean = seqmean;
 924         local->fFlags++; //means we have more than one pad
 925         local->fLastMergedPad = pad;
 926
 927         currentPt[ncurrent] = local;
 928         ncurrent++;
 929
 930         break;
 931       } //Checking for match at previous pad
 932     } //Loop over results on previous pad.
 933
 934     if(newcluster && ncurrent<kPadArraySize){
 935       //Start a new cluster. Add it to the clusterlist, and update
 936       //the list of pointers to clusters in current pad.
 937       //current pad will be previous pad on next pad.
 938
 939       //Add to the clusterlist:
 940       AliClusterData *tmp = &clusterlist[ntotal];
 941       tmp->fTotalCharge = seqcharge;
 942       tmp->fPad = padmean;
 943       tmp->fPad2 = paderror;
 944       tmp->fTime = seqaverage;
 945       tmp->fTime2 = seqerror;
 946       tmp->fMean = seqmean;
 947       tmp->fFlags = 0;  //flags for single pad clusters
 948       tmp->fLastMergedPad = pad;
 949
 950       if(fDeconvPad){
 951         tmp->fChargeFalling = 0;
 952         tmp->fLastCharge = seqcharge;
 953       }
 954
 955       //Update list of pointers to previous pad:
 956       currentPt[ncurrent] = &clusterlist[ntotal];
 957       ntotal++;
 958       ncurrent++;
 959     }
 960
 961     if(fDeconvTime)
 962       if(newbin >= 0) goto redo;
 963
 964     // to prevent endless loop
 965     if(time >= AliHLTTPCTransform::GetNTimeBins()){
 966       HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
 967       break;
 968     }
 969
 970
 971     if(!readValue) break; //No more value
 972
 973     if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
 974       HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
 975       break;
 976     }
 977
 978     if(fCurrentRow != newRow){
 979       WriteClusters(ntotal,clusterlist);
 980
 981       lastpad = 123456789;
 982
 983       currentPt = pad2;
 984       previousPt = pad1;
 985       nprevious=0;
 986       ncurrent=0;
 987       ntotal=0;
 988
 989       fCurrentRow = newRow;
 990     }
 991
 992     pad = newPad;
 993     time = newTime;
 994
 995   } // END while(readValue)
 996
 997   WriteClusters(ntotal,clusterlist);
 998
 999   HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
1000
1001 } // ENDEND
1002
1003 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list){
1004   // see header file for class documentation
1005
1006   //write cluster to output pointer
1007   Int_t thisrow=-1,thissector=-1;
1008   UInt_t counter = fNClusters;
1009
1010   for(int j=0; j<nclusters; j++)
1011     {
1012
1013
1014
1015       if(!list[j].fFlags){
1016         if(fDoMC){
1017           fClustersMCInfo.erase(fClustersMCInfo.begin()+j); // remove the mc ifo for this cluster since it is not taken into account
1018         }
1019         continue; //discard single pad clusters
1020       }
1021       if(list[j].fTotalCharge < fThreshold){
1022         if(fDoMC){
1023           fClustersMCInfo.erase(fClustersMCInfo.begin()+j); // remove the mc ifo for this cluster since it is not taken into account
1024         }
1025         continue; //noise cluster
1026       }
1027       Float_t xyz[3];
1028       Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
1029       Float_t fpad2=fXYErr*fXYErr; //fixed given error
1030       Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
1031       Float_t ftime2=fZErr*fZErr;  //fixed given error
1032
1033
1034
1035       if(fUnsorted){
1036         fCurrentRow=list[j].fRow;
1037       }
1038
1039
1040       if(fCalcerr) { //calc the errors, otherwice take the fixed error
1041         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
1042         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
1043         //      Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
1044         Float_t sy2=(Float_t)list[j].fPad2 * list[j].fTotalCharge - (Float_t)list[j].fPad * list[j].fPad;
1045         sy2/=q2;
1046         if(sy2 < 0) {
1047             LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1048               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1049             continue;
1050         } else {
1051           if(!fRawSP){
1052             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
1053             if(sy2 != 0){
1054               fpad2*=0.108; //constants are from offline studies
1055               if(patch<2)
1056                 fpad2*=2.07;
1057             }
1058           } else fpad2=sy2; //take the width not the error
1059         }
1060         //      Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
1061         Float_t sz2=(Float_t)list[j].fTime2*list[j].fTotalCharge - (Float_t)list[j].fTime*list[j].fTime;
1062         sz2/=q2;
1063         if(sz2 < 0){
1064           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1065             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1066           continue;
1067         } else {
1068           if(!fRawSP){
1069             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
1070             if(sz2 != 0) {
1071               ftime2 *= 0.169; //constants are from offline studies
1072               if(patch<2)
1073                 ftime2 *= 1.77;
1074             }
1075           } else ftime2=sz2; //take the width, not the error
1076         }
1077       }
1078       if(fStdout==kTRUE)
1079         HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
1080
1081       if(!fRawSP){
1082         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
1083         AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
1084
1085         if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
1086           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
1087         if(fNClusters >= fMaxNClusters)
1088           {
1089             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
1090               <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
1091             return;
1092           }
1093
1094         fSpacePointData[counter].fX = xyz[0];
1095         //      fSpacePointData[counter].fY = xyz[1];
1096         if(fCurrentSlice<18){
1097           fSpacePointData[counter].fY = xyz[1];
1098         }
1099         else{
1100           fSpacePointData[counter].fY = -1*xyz[1];
1101         }
1102         fSpacePointData[counter].fZ = xyz[2];
1103
1104       } else {
1105         fSpacePointData[counter].fX = fCurrentRow;
1106         fSpacePointData[counter].fY = fpad;
1107         fSpacePointData[counter].fZ = ftime;
1108       }
1109
1110       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
1111       fSpacePointData[counter].fPadRow = fCurrentRow;
1112       fSpacePointData[counter].fSigmaY2 = fpad2;
1113       fSpacePointData[counter].fSigmaZ2  = ftime2;
1114
1115       fSpacePointData[counter].fQMax = list[j].fQMax;
1116
1117       fSpacePointData[counter].fUsed = kFALSE;         // only used / set in AliHLTTPCDisplay
1118       fSpacePointData[counter].fTrackN = -1;           // only used / set in AliHLTTPCDisplay
1119
1120       Int_t patch=fCurrentPatch;
1121       if(patch==-1) patch=0; //never store negative patch number
1122       fSpacePointData[counter].fID = counter
1123         +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
1124
1125 #ifdef do_mc
1126       Int_t trackID[3];
1127       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
1128
1129       fSpacePointData[counter].fTrackID[0] = trackID[0];
1130       fSpacePointData[counter].fTrackID[1] = trackID[1];
1131       fSpacePointData[counter].fTrackID[2] = trackID[2];
1132
1133 #endif
1134
1135       fNClusters++;
1136       counter++;
1137     }
1138 }
1139
1140 // STILL TO FIX  ----------------------------------------------------------------------------
1141
1142 #ifdef do_mc
1143 void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID){
1144   // see header file for class documentation
1145
1146   //get mc id
1147   AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
1148
1149   trackID[0]=trackID[1]=trackID[2]=-2;
1150   for(Int_t i=fFirstRow; i<=fLastRow; i++){
1151     if(rowPt->fRow < (UInt_t)fCurrentRow){
1152       AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
1153       continue;
1154     }
1155     AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
1156     for(UInt_t j=0; j<rowPt->fNDigit; j++){
1157       Int_t cpad = digPt[j].fPad;
1158       Int_t ctime = digPt[j].fTime;
1159       if(cpad != pad) continue;
1160       if(ctime != time) continue;
1161
1162       trackID[0] = digPt[j].fTrackID[0];
1163       trackID[1] = digPt[j].fTrackID[1];
1164       trackID[2] = digPt[j].fTrackID[2];
1165
1166       break;
1167     }
1168     break;
1169   }
1170 }
1171 #endif
1172
1173
1174
1175 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list){//This is used when using the AliHLTTPCClusters class for cluster data
1176   // see header file for class documentation
1177
1178   //write cluster to output pointer
1179   Int_t thisrow,thissector;
1180   UInt_t counter = fNClusters;
1181
1182   for(int j=0; j<nclusters; j++)
1183     {
1184       if(!list[j].fFlags) continue; //discard single pad clusters
1185       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
1186
1187       Float_t xyz[3];
1188       Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
1189       Float_t fpad2=fXYErr*fXYErr; //fixed given error
1190       Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
1191       Float_t ftime2=fZErr*fZErr;  //fixed given error
1192
1193
1194       if(fCalcerr) { //calc the errors, otherwice take the fixed error
1195         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
1196         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
1197         Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
1198         sy2/=q2;
1199         if(sy2 < 0) {
1200             LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1201               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1202             continue;
1203         } else {
1204           if(!fRawSP){
1205             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
1206             if(sy2 != 0){
1207               fpad2*=0.108; //constants are from offline studies
1208               if(patch<2)
1209                 fpad2*=2.07;
1210             }
1211           } else fpad2=sy2; //take the width not the error
1212         }
1213         Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
1214         sz2/=q2;
1215         if(sz2 < 0){
1216           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1217             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1218           continue;
1219         } else {
1220           if(!fRawSP){
1221             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
1222             if(sz2 != 0) {
1223               ftime2 *= 0.169; //constants are from offline studies
1224               if(patch<2)
1225                 ftime2 *= 1.77;
1226             }
1227           } else ftime2=sz2; //take the width, not the error
1228         }
1229       }
1230       if(fStdout==kTRUE)
1231         HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
1232
1233       if(!fRawSP){
1234         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
1235         AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
1236
1237         if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
1238           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
1239         if(fNClusters >= fMaxNClusters)
1240           {
1241             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
1242               <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
1243             return;
1244           }
1245
1246         fSpacePointData[counter].fX = xyz[0];
1247         //      fSpacePointData[counter].fY = xyz[1];
1248         if(fCurrentSlice<18){
1249           fSpacePointData[counter].fY = xyz[1];
1250         }
1251         else{
1252           fSpacePointData[counter].fY = -1*xyz[1];
1253         }
1254         fSpacePointData[counter].fZ = xyz[2];
1255
1256       } else {
1257         fSpacePointData[counter].fX = fCurrentRow;
1258         fSpacePointData[counter].fY = fpad;
1259         fSpacePointData[counter].fZ = ftime;
1260       }
1261
1262       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
1263       fSpacePointData[counter].fPadRow = fCurrentRow;
1264       fSpacePointData[counter].fSigmaY2 = fpad2;
1265       fSpacePointData[counter].fSigmaZ2  = ftime2;
1266
1267       fSpacePointData[counter].fQMax = list[j].fQMax;
1268
1269       fSpacePointData[counter].fUsed = kFALSE;         // only used / set in AliHLTTPCDisplay
1270       fSpacePointData[counter].fTrackN = -1;           // only used / set in AliHLTTPCDisplay
1271
1272       Int_t patch=fCurrentPatch;
1273       if(patch==-1) patch=0; //never store negative patch number
1274       fSpacePointData[counter].fID = counter
1275         +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
1276
1277 #ifdef do_mc
1278       Int_t trackID[3];
1279       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
1280
1281       fSpacePointData[counter].fTrackID[0] = trackID[0];
1282       fSpacePointData[counter].fTrackID[1] = trackID[1];
1283       fSpacePointData[counter].fTrackID[2] = trackID[2];
1284
1285 #endif
1286
1287       fNClusters++;
1288       counter++;
1289     }
1290 }