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