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