HLT/TPCLib/AliHLTTPCClusterFinder.cxx

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