HLT/TPCLib/AliHLTTPCClusterFinder.cxx

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