HLT/TPCLib/AliHLTTPCClusterFinder.cxx

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