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