HLT/TPCLib/AliHLTTPCClusterFinder.cxx

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