HLT/TPCLib/AliHLTTPCClusterFinder.cxx

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