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
 381                 prevSignal=bunchData[i];
 382                 time++;
 383                 indexInBunchData--;
 384               }
 385               if(candidate.fTotalCharge>0){
 386                 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
 387                 candidate.fPad=candidate.fTotalCharge*pad;
 388                 candidate.fPad2=candidate.fPad*pad;
 389                 candidate.fLastMergedPad=pad;
 390                 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
 391               }
 392               fRowPadVector[row][pad]->AddClusterCandidate(candidate);
 393               fRowPadVector[row][pad]->AddCandidateDigits(fMCDigits);
 394               if(indexInBunchData<fDigitReader->GetBunchSize()-1){
 395                 moreDataInBunch=kFALSE;
 396               }
 397             }while(moreDataInBunch);
 398           }
 399           else{
 400             const UInt_t *bunchData= fDigitReader->GetSignals();
 401             do{
 402               AliHLTTPCClusters candidate;
 403               for(Int_t i=indexInBunchData;i<fDigitReader->GetBunchSize();i++){
 404                 // Checks if one need to deconvolute the signals
 405                 if(bunchData[i]>prevSignal && signalFalling==kTRUE){
 406                   if(i<fDigitReader->GetBunchSize()-1){ // means there are more than one signal left in the bunch
 407                     moreDataInBunch=kTRUE;
 408                     prevSignal=0;
 409                   }
 410                   break;
 411                 }
 412
 413                 // Checks if the signal is 0, then quit processing the data.
 414                 if(bunchData[i]==0 && i<fDigitReader->GetBunchSize()-1){//means we have 0 data fom the rcu, might happen depending on the configuration settings
 415                   moreDataInBunch=kTRUE;
 416                   prevSignal=0;
 417                   break;
 418                 }
 419
 420                 if(prevSignal>bunchData[i]){//means the peak of the signal has been reached and deconvolution will happen if the signal rise again.
 421                   signalFalling=kTRUE;
 422                 }
 423
 424                 candidate.fTotalCharge+=bunchData[i];
 425                 candidate.fTime += time*bunchData[i];
 426                 candidate.fTime2 += time*time*bunchData[i];
 427                 if(bunchData[i]>candidate.fQMax){
 428                   candidate.fQMax=bunchData[i];
 429                 }
 430
 431                 prevSignal=bunchData[i];
 432                 time++;
 433                 indexInBunchData++;
 434               }
 435               if(candidate.fTotalCharge>0){
 436                 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
 437                 candidate.fPad=candidate.fTotalCharge*pad;
 438                 candidate.fPad2=candidate.fPad*pad;
 439                 candidate.fLastMergedPad=pad;
 440                 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
 441               }
 442               fRowPadVector[row][pad]->AddClusterCandidate(candidate);
 443               if(indexInBunchData<fDigitReader->GetBunchSize()-1){
 444                 moreDataInBunch=kFALSE;
 445               }
 446             }while(moreDataInBunch);
 447           }
 448         }
 449       }
 450     }
 451   }
 452 }
 453
 454 Bool_t AliHLTTPCClusterFinder::ComparePads(AliHLTTPCPad *nextPad,AliHLTTPCClusters* cluster,Int_t nextPadToRead){
 455   // see header file for class documentation
 456
 457   //Checking if we have a match on the next pad
 458   for(UInt_t candidateNumber=0;candidateNumber<nextPad->fClusterCandidates.size();candidateNumber++){
 459     if(nextPad->fUsedClusterCandidates[candidateNumber] == 1){
 460       continue;
 461     }
 462     AliHLTTPCClusters *candidate =&nextPad->fClusterCandidates[candidateNumber];
 463     //    if(cluster->fMean-candidate->fMean==1 || candidate->fMean-cluster->fMean==1 || cluster->fMean-candidate->fMean==0){
 464
 465     if( abs((Int_t)(cluster->fMean - candidate->fMean)) <= fTimeMeanDiff ){
 466       if(fDeconvPad){
 467         if(candidate->fTotalCharge<fTotalChargeOfPreviousClusterCandidate){//peak is reached
 468           fChargeOfCandidatesFalling=kTRUE;
 469         }
 470         if(candidate->fTotalCharge>fTotalChargeOfPreviousClusterCandidate && fChargeOfCandidatesFalling==kTRUE){//we have deconvolution
 471           return kFALSE;
 472         }
 473       }
 474       cluster->fMean=candidate->fMean;
 475       cluster->fTotalCharge+=candidate->fTotalCharge;
 476       cluster->fTime += candidate->fTime;
 477       cluster->fTime2 += candidate->fTime2;
 478       cluster->fPad+=candidate->fPad;
 479       cluster->fPad2+=candidate->fPad2;
 480       cluster->fLastMergedPad=candidate->fPad;
 481       if(candidate->fQMax>cluster->fQMax){
 482         cluster->fQMax=candidate->fQMax;
 483       }
 484       if(fDoMC){
 485         FillMCClusterVector(nextPad->GetCandidateDigits(candidateNumber));
 486       }
 487
 488       if(fDoPadSelection){
 489         UInt_t rowNo = nextPad->GetRowNumber();
 490         UInt_t padNo = nextPad->GetPadNumber();
 491         if(padNo-1>0){
 492           fRowPadVector[rowNo][padNo-2]->fSelectedPad=kTRUE;
 493           fRowPadVector[rowNo][padNo-2]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo-2);
 494         }
 495         fRowPadVector[rowNo][padNo-1]->fSelectedPad=kTRUE;// quick solution to set the first pad to selected
 496         fRowPadVector[rowNo][padNo-1]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo-1);
 497         fRowPadVector[rowNo][padNo]->fSelectedPad=kTRUE;
 498         fRowPadVector[rowNo][padNo]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo);
 499       }
 500
 501       //setting the matched pad to used
 502       nextPad->fUsedClusterCandidates[candidateNumber]=1;
 503       nextPadToRead++;
 504       if(nextPadToRead<(Int_t)fNumberOfPadsInRow[fRowOfFirstCandidate]){
 505         nextPad=fRowPadVector[fRowOfFirstCandidate][nextPadToRead];
 506         ComparePads(nextPad,cluster,nextPadToRead);
 507       }
 508       else{
 509         return kFALSE;
 510       }
 511     }
 512   }
 513   return kFALSE;
 514 }
 515
 516 Int_t AliHLTTPCClusterFinder::FillHWAddressList(AliHLTUInt16_t *hwaddlist, Int_t maxHWadd){
 517   // see header file for class documentation
 518
 519   Int_t counter=0;
 520   for(UInt_t row=0;row<fNumberOfRows;row++){
 521     for(UInt_t pad=0;pad<fNumberOfPadsInRow[row]-1;pad++){
 522       if(fRowPadVector[row][pad]->fSelectedPad){
 523        if(counter<maxHWadd){
 524          hwaddlist[counter]=(AliHLTUInt16_t)fRowPadVector[row][pad]->fHWAddress;
 525          counter++;
 526        }
 527        else{
 528          HLTWarning("To many hardwareaddresses, skip adding");
 529        }
 530
 531       }
 532     }
 533   }
 534   return counter;
 535 }
 536
 537
 538 Int_t AliHLTTPCClusterFinder::FillOutputMCInfo(AliHLTTPCClusterFinder::ClusterMCInfo * outputMCInfo, Int_t maxNumberOfClusterMCInfo){
 539   // see header file for class documentation
 540
 541   Int_t counter=0;
 542
 543   for(UInt_t mc=0;mc<fClustersMCInfo.size();mc++){
 544     if(counter<maxNumberOfClusterMCInfo){
 545       outputMCInfo[counter] = fClustersMCInfo[mc];
 546       counter++;
 547     }
 548     else{
 549       HLTWarning("To much MCInfo has been added (no more space), skip adding");
 550     }
 551   }
 552   return counter;
 553 }
 554
 555 void AliHLTTPCClusterFinder::FindClusters(){
 556   // see header file for function documentation
 557
 558   AliHLTTPCClusters* tmpCandidate=NULL;
 559   for(UInt_t row=0;row<fNumberOfRows;row++){
 560     fRowOfFirstCandidate=row;
 561     for(UInt_t pad=0;pad<fNumberOfPadsInRow[row];pad++){
 562       AliHLTTPCPad *tmpPad=fRowPadVector[row][pad];
 563       for(size_t candidate=0;candidate<tmpPad->fClusterCandidates.size();candidate++){
 564         if(tmpPad->fUsedClusterCandidates[candidate]){
 565           continue;
 566         }
 567         tmpCandidate=&tmpPad->fClusterCandidates[candidate];
 568         UInt_t tmpTotalCharge=tmpCandidate->fTotalCharge;
 569
 570         if(fDoMC){
 571           fClusterMCVector.clear();
 572           FillMCClusterVector(tmpPad->GetCandidateDigits(candidate));
 573         }
 574
 575         ComparePads(fRowPadVector[row][pad+1],tmpCandidate,pad+1);
 576         if(tmpCandidate->fTotalCharge>tmpTotalCharge){
 577           //we have a cluster
 578           fClusters.push_back(*tmpCandidate);
 579           if(fDoMC){
 580             //sort the vector (large->small) according to weight and remove elements above 2 (keep 0 1 and 2)
 581             sort(fClusterMCVector.begin(),fClusterMCVector.end(), MCWeight::CompareWeights );
 582             ClusterMCInfo tmpClusterMCInfo;
 583
 584             MCWeight zeroMC;
 585             zeroMC.fMCID=-1;
 586             zeroMC.fWeight=0;
 587
 588             if(fClusterMCVector.size()>0){
 589               tmpClusterMCInfo.fClusterID[0]=fClusterMCVector.at(0);
 590             }
 591             else{
 592               tmpClusterMCInfo.fClusterID[0]=zeroMC;
 593             }
 594
 595             if(fClusterMCVector.size()>1){
 596             tmpClusterMCInfo.fClusterID[1]=fClusterMCVector.at(1);
 597             }
 598             else{
 599               tmpClusterMCInfo.fClusterID[1]=zeroMC;
 600             }
 601
 602             if(fClusterMCVector.size()>2){
 603             tmpClusterMCInfo.fClusterID[2]=fClusterMCVector.at(2);
 604             }
 605             else{
 606               tmpClusterMCInfo.fClusterID[2]=zeroMC;
 607             }
 608
 609             fClustersMCInfo.push_back(tmpClusterMCInfo);
 610           }
 611
 612         }
 613       }
 614       tmpPad->ClearCandidates();
 615     }
 616     fRowPadVector[row][fNumberOfPadsInRow[row]]->ClearCandidates();
 617   }
 618
 619   HLTInfo("Found %d clusters.",fClusters.size());
 620
 621   //TODO:  Change so it stores AliHLTTPCSpacePointData directly, instead of this copying
 622
 623   AliClusterData * clusterlist = new AliClusterData[fClusters.size()]; //Clusterlist
 624   for(unsigned int i=0;i<fClusters.size();i++){
 625     clusterlist[i].fTotalCharge = fClusters[i].fTotalCharge;
 626     clusterlist[i].fPad = fClusters[i].fPad;
 627     clusterlist[i].fPad2 = fClusters[i].fPad2;
 628     clusterlist[i].fTime = fClusters[i].fTime;
 629     clusterlist[i].fTime2 = fClusters[i].fTime2;
 630     clusterlist[i].fMean = fClusters[i].fMean;
 631     clusterlist[i].fFlags = fClusters[i].fFlags;
 632     clusterlist[i].fChargeFalling = fClusters[i].fChargeFalling;
 633     clusterlist[i].fLastCharge = fClusters[i].fLastCharge;
 634     clusterlist[i].fLastMergedPad = fClusters[i].fLastMergedPad;
 635     clusterlist[i].fRow = fClusters[i].fRowNumber;
 636     clusterlist[i].fQMax = fClusters[i].fQMax;
 637   }
 638
 639   WriteClusters(fClusters.size(),clusterlist);
 640   delete [] clusterlist;
 641   fClusters.clear();
 642   if( fReleaseMemory ) DeInitializePadArray();// call this when  the -releaseMemory flag is set
 643 }
 644
 645
 646 Bool_t AliHLTTPCClusterFinder::UpdateCalibDB(){
 647
 648   //update the db
 649   AliTPCcalibDB::Instance()->Update();
 650
 651   Bool_t ret = 1;
 652
 653   //uptate the transform class
 654
 655   fOfflineTransform = AliTPCcalibDB::Instance()->GetTransform();
 656   if(!fOfflineTransform){
 657     HLTError("AliHLTTPCClusterFinder()::UpdateCAlibDB::  Offline transform not in AliTPCcalibDB.");
 658     ret = 0;
 659   }
 660   else{
 661     fOfflineTransform->SetCurrentRecoParam(&fOfflineTPCRecoParam);
 662   }
 663
 664   fOfflineTPCParam = AliTPCcalibDB::Instance()->GetParameters();
 665   if( !fOfflineTPCParam ){
 666     HLTError("AliHLTTPCClusterFinder()::UpdateCAlibDB::  Offline TPC parameters not in AliTPCcalibDB.");
 667     ret = 0;
 668   } else {
 669     fOfflineTPCParam->Update();
 670     fOfflineTPCParam->ReadGeoMatrices();
 671   }
 672
 673   return ret;
 674 }
 675
 676 //---------------------------------- Under this line the old sorted clusterfinder functions can be found --------------------------------
 677
 678
 679 void AliHLTTPCClusterFinder::PrintClusters(){
 680   // see header file for class documentation
 681
 682   for(size_t i=0;i<fClusters.size();i++){
 683     HLTInfo("Cluster number: %d",i);
 684     HLTInfo("Row: %d \t Pad: %d",fClusters[i].fRowNumber,fClusters[i].fPad/fClusters[i].fTotalCharge);
 685     HLTInfo("Total Charge:   %d",fClusters[i].fTotalCharge);
 686     HLTInfo("fPad:           %d",fClusters[i].fPad);
 687     HLTInfo("PadError:       %d",fClusters[i].fPad2);
 688     HLTInfo("TimeMean:       %d",fClusters[i].fTime/fClusters[i].fTotalCharge);
 689     HLTInfo("TimeError:      %d",fClusters[i].fTime2);
 690     HLTInfo("EndOfCluster:");
 691   }
 692 }
 693
 694 void AliHLTTPCClusterFinder::FillMCClusterVector(vector<AliHLTTPCDigitData> digitData){
 695   // see header file for class documentation
 696
 697   for(UInt_t d=0;d<digitData.size();d++){
 698     Int_t nIDsInDigit = (digitData.at(d).fTrackID[0]>=0) + (digitData.at(d).fTrackID[1]>=0) + (digitData.at(d).fTrackID[2]>=0);
 699     for(Int_t id=0; id<3; id++){
 700       if(digitData.at(d).fTrackID[id]>=0){
 701         Bool_t matchFound = kFALSE;
 702         MCWeight mc;
 703         mc.fMCID = digitData.at(d).fTrackID[id];
 704         mc.fWeight = ((Float_t)digitData.at(d).fCharge)/nIDsInDigit;
 705         for(UInt_t i=0;i<fClusterMCVector.size();i++){
 706           if(mc.fMCID == fClusterMCVector.at(i).fMCID){
 707             fClusterMCVector.at(i).fWeight += mc.fWeight;
 708             matchFound = kTRUE;
 709           }
 710         }
 711         if(matchFound == kFALSE){
 712           fClusterMCVector.push_back(mc);
 713         }
 714       }
 715     }
 716   }
 717 }
 718
 719
 720 void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
 721   //set input pointer
 722   fPtr = (UChar_t*)ptr;
 723   fSize = size;
 724 }
 725
 726 void AliHLTTPCClusterFinder::ProcessDigits(){
 727   // see header file for class documentation
 728
 729   int iResult=0;
 730   bool readValue = true;
 731   Int_t newRow = 0;
 732   Int_t rowOffset = 0;
 733   UShort_t time=0,newTime=0;
 734   UInt_t pad=0,newPad=0;
 735   AliHLTTPCSignal_t charge=0;
 736
 737   fNClusters = 0;
 738
 739   // initialize block for reading packed data
 740   iResult=fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
 741   if (iResult<0) return;
 742
 743   readValue = fDigitReader->Next();
 744
 745   // Matthias 08.11.2006 the following return would cause termination without writing the
 746   // ClusterData and thus would block the component. I just want to have the commented line
 747   // here for information
 748   //if (!readValue)return;
 749
 750   pad = fDigitReader->GetPad();
 751   time = fDigitReader->GetTime();
 752   fCurrentRow = fDigitReader->GetRow();
 753
 754   if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
 755     rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
 756
 757   fCurrentRow += rowOffset;
 758
 759   UInt_t lastpad = 123456789;
 760   const UInt_t kPadArraySize=5000;
 761   const UInt_t kClusterListSize=10000;
 762   AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
 763   AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
 764   AliClusterData clusterlist[kClusterListSize]; //Clusterlist
 765
 766   AliClusterData **currentPt;  //List of pointers to the current pad
 767   AliClusterData **previousPt; //List of pointers to the previous pad
 768   currentPt = pad2;
 769   previousPt = pad1;
 770   UInt_t nprevious=0,ncurrent=0,ntotal=0;
 771
 772   /* quick implementation of baseline calculation and zero suppression
 773      open a pad object for each pad and delete it after processing.
 774      later a list of pad objects with base line history can be used
 775      The whole thing only works if we really get unprocessed raw data, if
 776      the data is already zero suppressed, there might be gaps in the time
 777      bins.
 778    */
 779   Int_t gatingGridOffset=50;
 780   if(fFirstTimeBin>0){
 781     gatingGridOffset=fFirstTimeBin;
 782   }
 783   AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
 784   // just to make later conversion to a list of objects easier
 785   AliHLTTPCPad* pCurrentPad=NULL;
 786   /*
 787     if (fSignalThreshold>=0) {
 788     pCurrentPad=&baseline;
 789     baseline.SetThreshold(fSignalThreshold);
 790   }
 791   */
 792   while ( readValue!=0 && iResult>=0){   // Reads through all digits in block
 793     iResult=0;
 794
 795     if(pad != lastpad){
 796       //This is a new pad
 797
 798       //Switch the lists:
 799       if(currentPt == pad2){
 800         currentPt = pad1;
 801         previousPt = pad2;
 802       }
 803       else {
 804         currentPt = pad2;
 805         previousPt = pad1;
 806       }
 807       nprevious = ncurrent;
 808       ncurrent = 0;
 809       if(pad != lastpad+1){
 810         //this happens if there is a pad with no signal.
 811         nprevious = ncurrent = 0;
 812       }
 813       lastpad = pad;
 814     }
 815
 816     Bool_t newcluster = kTRUE;
 817     UInt_t seqcharge=0,seqaverage=0,seqerror=0;
 818     AliHLTTPCSignal_t lastcharge=0;
 819     UInt_t bLastWasFalling=0;
 820     Int_t newbin=-1;
 821
 822
 823     if(fDeconvTime){
 824       redo: //This is a goto.
 825
 826       if(newbin > -1){
 827         //bin = newbin;
 828         newbin = -1;
 829       }
 830
 831       lastcharge=0;
 832       bLastWasFalling = 0;
 833     }
 834
 835     while(iResult>=0){ //Loop over time bins of current pad
 836       // read all the values for one pad at once to calculate the base line
 837       if (pCurrentPad) {
 838         if (!pCurrentPad->IsStarted()) {
 839           //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
 840           pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
 841           if ((pCurrentPad->StartEvent())>=0) {
 842             do {
 843               if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
 844               if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
 845               pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
 846               //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
 847             } while ((readValue = fDigitReader->Next())!=0);
 848           }
 849           pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
 850           if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
 851             //HLTDebug("no data available after zero suppression");
 852             pCurrentPad->StopEvent();
 853             pCurrentPad->ResetHistory();
 854             break;
 855           }
 856           time=pCurrentPad->GetCurrentPosition();
 857           if (time>pCurrentPad->GetSize()) {
 858             HLTError("invalid time bin for pad");
 859             break;
 860           }
 861         }
 862       }
 863       if (pCurrentPad) {
 864         Float_t occupancy=pCurrentPad->GetOccupancy();
 865         //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
 866         if ( occupancy < fOccupancyLimit ) {
 867           charge = pCurrentPad->GetCorrectedData();
 868         } else {
 869           charge = 0;
 870           //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
 871         }
 872       } else {
 873         charge = fDigitReader->GetSignal();
 874       }
 875       //HLTDebug("get next charge value: position %d charge %d", time, charge);
 876
 877
 878       // CHARGE DEBUG
 879       if (fDigitReader->GetRow() == 90){
 880 /////     LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90")  << "PAD=" <<  fDigitReader->GetPad() << "  TIME=" <<  fDigitReader->GetTime()
 881           //                                       << "  SIGNAL=" <<  fDigitReader->GetSignal() << ENDLOG;
 882
 883       }
 884
 885       if(time >= AliHLTTPCTransform::GetNTimeBins()){
 886         HLTWarning("Pad %d: Timebin (%d) out of range (%d)", pad, time, AliHLTTPCTransform::GetNTimeBins());
 887         iResult=-ERANGE;
 888       }
 889
 890
 891       //Get the current ADC-value
 892       if(fDeconvTime){
 893
 894         //Check if the last pixel in the sequence is smaller than this
 895         if(charge > lastcharge){
 896           if(bLastWasFalling){
 897             newbin = 1;
 898             break;
 899           }
 900         }
 901         else bLastWasFalling = 1; //last pixel was larger than this
 902         lastcharge = charge;
 903       }
 904
 905       //Sum the total charge of this sequence
 906       seqcharge += charge;
 907       seqaverage += time*charge;
 908       seqerror += time*time*charge;
 909
 910       if (pCurrentPad) {
 911
 912         if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
 913           pCurrentPad->StopEvent();
 914           pCurrentPad->ResetHistory();
 915           if(readValue) {
 916             newPad = fDigitReader->GetPad();
 917             newTime = fDigitReader->GetTime();
 918             newRow = fDigitReader->GetRow() + rowOffset;
 919           }
 920           break;
 921         }
 922
 923         newPad=pCurrentPad->GetPadNumber();
 924         newTime=pCurrentPad->GetCurrentPosition();
 925         newRow=pCurrentPad->GetRowNumber();
 926       } else {
 927       readValue = fDigitReader->Next();
 928       //Check where to stop:
 929       if(!readValue) break; //No more value
 930
 931       newPad = fDigitReader->GetPad();
 932       newTime = fDigitReader->GetTime();
 933       newRow = fDigitReader->GetRow() + rowOffset;
 934       }
 935
 936       if(newPad != pad)break; //new pad
 937       if(newTime != time+1) break; //end of sequence
 938       if(iResult<0) break;
 939
 940       // pad = newpad;    is equal
 941       time = newTime;
 942
 943     }//end loop over sequence
 944
 945     //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
 946     //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
 947     if (seqcharge<=0) {
 948       // with active zero suppression zero values are possible
 949       continue;
 950     }
 951
 952     //Calculate mean of sequence:
 953     Int_t seqmean=0;
 954     if(seqcharge)
 955       seqmean = seqaverage/seqcharge;
 956     else{
 957       LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
 958         <<"Error in data given to the cluster finder"<<ENDLOG;
 959       seqmean = 1;
 960       seqcharge = 1;
 961     }
 962
 963     //Calculate mean in pad direction:
 964     Int_t padmean = seqcharge*pad;
 965     Int_t paderror = pad*padmean;
 966
 967     //Compare with results on previous pad:
 968     for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
 969
 970       //dont merge sequences on the same pad twice
 971       if(previousPt[p]->fLastMergedPad==pad) continue;
 972
 973       Int_t difference = seqmean - previousPt[p]->fMean;
 974       if(difference < -fMatch) break;
 975
 976       if(difference <= fMatch){ //There is a match here!!
 977         AliClusterData *local = previousPt[p];
 978
 979         if(fDeconvPad){
 980           if(seqcharge > local->fLastCharge){
 981             if(local->fChargeFalling){ //The previous pad was falling
 982               break; //create a new cluster
 983             }
 984           }
 985           else local->fChargeFalling = 1;
 986           local->fLastCharge = seqcharge;
 987         }
 988
 989         //Don't create a new cluster, because we found a match
 990         newcluster = kFALSE;
 991
 992         //Update cluster on current pad with the matching one:
 993         local->fTotalCharge += seqcharge;
 994         local->fPad += padmean;
 995         local->fPad2 += paderror;
 996         local->fTime += seqaverage;
 997         local->fTime2 += seqerror;
 998         local->fMean = seqmean;
 999         local->fFlags++; //means we have more than one pad
1000         local->fLastMergedPad = pad;
1001
1002         currentPt[ncurrent] = local;
1003         ncurrent++;
1004
1005         break;
1006       } //Checking for match at previous pad
1007     } //Loop over results on previous pad.
1008
1009     if(newcluster && ncurrent<kPadArraySize){
1010       //Start a new cluster. Add it to the clusterlist, and update
1011       //the list of pointers to clusters in current pad.
1012       //current pad will be previous pad on next pad.
1013
1014       //Add to the clusterlist:
1015       AliClusterData *tmp = &clusterlist[ntotal];
1016       tmp->fTotalCharge = seqcharge;
1017       tmp->fPad = padmean;
1018       tmp->fPad2 = paderror;
1019       tmp->fTime = seqaverage;
1020       tmp->fTime2 = seqerror;
1021       tmp->fMean = seqmean;
1022       tmp->fFlags = 0;  //flags for single pad clusters
1023       tmp->fLastMergedPad = pad;
1024
1025       if(fDeconvPad){
1026         tmp->fChargeFalling = 0;
1027         tmp->fLastCharge = seqcharge;
1028       }
1029
1030       //Update list of pointers to previous pad:
1031       currentPt[ncurrent] = &clusterlist[ntotal];
1032       ntotal++;
1033       ncurrent++;
1034     }
1035
1036     if(fDeconvTime)
1037       if(newbin >= 0) goto redo;
1038
1039     // to prevent endless loop
1040     if(time >= AliHLTTPCTransform::GetNTimeBins()){
1041       HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
1042       break;
1043     }
1044
1045
1046     if(!readValue) break; //No more value
1047
1048     if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
1049       HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
1050       break;
1051     }
1052
1053     if(fCurrentRow != newRow){
1054       WriteClusters(ntotal,clusterlist);
1055
1056       lastpad = 123456789;
1057
1058       currentPt = pad2;
1059       previousPt = pad1;
1060       nprevious=0;
1061       ncurrent=0;
1062       ntotal=0;
1063
1064       fCurrentRow = newRow;
1065     }
1066
1067     pad = newPad;
1068     time = newTime;
1069
1070   } // END while(readValue)
1071
1072   WriteClusters(ntotal,clusterlist);
1073
1074   HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
1075
1076 } // ENDEND
1077
1078 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list){
1079   // see header file for class documentation
1080
1081   //write cluster to output pointer
1082   Int_t thisrow=-1,thissector=-1;
1083   UInt_t counter = fNClusters;
1084
1085   for(int j=0; j<nclusters; j++)
1086     {
1087
1088
1089
1090       if(!list[j].fFlags){
1091         if(fDoMC){
1092           if(j+(Int_t)fClustersMCInfo.size()-nclusters >=0 && j+fClustersMCInfo.size()-nclusters < fClustersMCInfo.size()){
1093             fClustersMCInfo.erase(fClustersMCInfo.begin()+j+fClustersMCInfo.size()-nclusters); // remove the mc info for this cluster since it is not taken into account
1094           }
1095         }
1096         continue; //discard single pad clusters
1097       }
1098       if(list[j].fTotalCharge < fThreshold){
1099         if(fDoMC){
1100           if(j+(Int_t)fClustersMCInfo.size()-nclusters >=0 && j+fClustersMCInfo.size()-nclusters < fClustersMCInfo.size()){
1101             fClustersMCInfo.erase(fClustersMCInfo.begin()+j+fClustersMCInfo.size()-nclusters); // remove the mc info for this cluster since it is not taken into account
1102           }
1103         }
1104         continue; //noise cluster
1105       }
1106       Float_t xyz[3];
1107       Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
1108       Float_t fpad2=fXYErr*fXYErr; //fixed given error
1109       Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
1110       Float_t ftime2=fZErr*fZErr;  //fixed given error
1111
1112
1113
1114       if(fUnsorted){
1115         fCurrentRow=list[j].fRow;
1116       }
1117
1118
1119       if(fCalcerr) { //calc the errors, otherwice take the fixed error
1120         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
1121         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
1122         //      Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
1123         Float_t sy2=(Float_t)list[j].fPad2 * list[j].fTotalCharge - (Float_t)list[j].fPad * list[j].fPad;
1124         sy2/=q2;
1125         if(sy2 < 0) {
1126             LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1127               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1128             continue;
1129         } else {
1130           if(!fRawSP){
1131             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
1132             if(sy2 != 0){
1133               fpad2*=0.108; //constants are from offline studies
1134               if(patch<2)
1135                 fpad2*=2.07;
1136             }
1137           } else fpad2=sy2; //take the width not the error
1138         }
1139         //      Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
1140         Float_t sz2=(Float_t)list[j].fTime2*list[j].fTotalCharge - (Float_t)list[j].fTime*list[j].fTime;
1141         sz2/=q2;
1142         if(sz2 < 0){
1143           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1144             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1145           continue;
1146         } else {
1147           if(!fRawSP){
1148             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
1149             if(sz2 != 0) {
1150               ftime2 *= 0.169; //constants are from offline studies
1151               if(patch<2)
1152                 ftime2 *= 1.77;
1153             }
1154           } else ftime2=sz2; //take the width, not the error
1155         }
1156       }
1157       if(fStdout==kTRUE)
1158         HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
1159
1160       if(!fRawSP){
1161         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
1162
1163         if(fOfflineTransform == NULL){
1164           AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
1165
1166           if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
1167                           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
1168           if(fNClusters >= fMaxNClusters)
1169             {
1170               LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
1171                 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
1172               return;
1173             }
1174
1175           fSpacePointData[counter].fX = xyz[0];
1176           //    fSpacePointData[counter].fY = xyz[1];
1177           if(fCurrentSlice<18){
1178             fSpacePointData[counter].fY = xyz[1];
1179           }
1180           else{
1181             fSpacePointData[counter].fY = -1*xyz[1];
1182           }
1183           fSpacePointData[counter].fZ = xyz[2];
1184         }
1185         else{
1186           Double_t x[3]={thisrow,fpad+.5,ftime};
1187           Int_t iSector[1]={thissector};
1188           fOfflineTransform->Transform(x,iSector,0,1);
1189           double y[3] = {x[0], x[1], x[2] };
1190
1191           if( fOfflineTPCParam && thissector<fOfflineTPCParam->GetNSector() ){
1192             TGeoHMatrix  *alignment = fOfflineTPCParam->GetClusterMatrix( thissector );
1193             if ( alignment ) alignment->LocalToMaster( x, y);
1194           }
1195
1196           fSpacePointData[counter].fX = y[0];
1197           fSpacePointData[counter].fY = y[1];
1198           fSpacePointData[counter].fZ = y[2];
1199         }
1200
1201       }
1202       else {
1203         fSpacePointData[counter].fX = fCurrentRow;
1204         fSpacePointData[counter].fY = fpad;
1205         fSpacePointData[counter].fZ = ftime;
1206       }
1207
1208       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
1209       fSpacePointData[counter].fPadRow = fCurrentRow;
1210       fSpacePointData[counter].fSigmaY2 = fpad2;
1211       fSpacePointData[counter].fSigmaZ2  = ftime2;
1212
1213       fSpacePointData[counter].fQMax = list[j].fQMax;
1214
1215       fSpacePointData[counter].fUsed = kFALSE;         // only used / set in AliHLTTPCDisplay
1216       fSpacePointData[counter].fTrackN = -1;           // only used / set in AliHLTTPCDisplay
1217
1218       Int_t patch=fCurrentPatch;
1219       if(patch==-1) patch=0; //never store negative patch number
1220       fSpacePointData[counter].fID = counter
1221         +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
1222
1223 #ifdef do_mc
1224       Int_t trackID[3];
1225       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
1226
1227       fSpacePointData[counter].fTrackID[0] = trackID[0];
1228       fSpacePointData[counter].fTrackID[1] = trackID[1];
1229       fSpacePointData[counter].fTrackID[2] = trackID[2];
1230
1231 #endif
1232
1233       fNClusters++;
1234       counter++;
1235     }
1236 }
1237
1238 // STILL TO FIX  ----------------------------------------------------------------------------
1239
1240 #ifdef do_mc
1241 void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID) const {
1242   // see header file for class documentation
1243
1244   //get mc id
1245   AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
1246
1247   trackID[0]=trackID[1]=trackID[2]=-2;
1248   for(Int_t i=fFirstRow; i<=fLastRow; i++){
1249     if(rowPt->fRow < (UInt_t)fCurrentRow){
1250       AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
1251       continue;
1252     }
1253     AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
1254     for(UInt_t j=0; j<rowPt->fNDigit; j++){
1255       Int_t cpad = digPt[j].fPad;
1256       Int_t ctime = digPt[j].fTime;
1257       if(cpad != pad) continue;
1258       if(ctime != time) continue;
1259
1260       trackID[0] = digPt[j].fTrackID[0];
1261       trackID[1] = digPt[j].fTrackID[1];
1262       trackID[2] = digPt[j].fTrackID[2];
1263
1264       break;
1265     }
1266     break;
1267   }
1268 }
1269 #endif
1270
1271
1272
1273 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list){//This is used when using the AliHLTTPCClusters class for cluster data
1274   // see header file for class documentation
1275
1276   //write cluster to output pointer
1277   Int_t thisrow,thissector;
1278   UInt_t counter = fNClusters;
1279
1280   for(int j=0; j<nclusters; j++)
1281     {
1282       if(!list[j].fFlags) continue; //discard single pad clusters
1283       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
1284
1285       Float_t xyz[3];
1286       Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
1287       Float_t fpad2=fXYErr*fXYErr; //fixed given error
1288       Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
1289       Float_t ftime2=fZErr*fZErr;  //fixed given error
1290
1291
1292       if(fCalcerr) { //calc the errors, otherwice take the fixed error
1293         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
1294         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
1295         Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
1296         sy2/=q2;
1297         if(sy2 < 0) {
1298             LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1299               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1300             continue;
1301         } else {
1302           if(!fRawSP){
1303             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
1304             if(sy2 != 0){
1305               fpad2*=0.108; //constants are from offline studies
1306               if(patch<2)
1307                 fpad2*=2.07;
1308             }
1309           } else fpad2=sy2; //take the width not the error
1310         }
1311         Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
1312         sz2/=q2;
1313         if(sz2 < 0){
1314           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1315             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1316           continue;
1317         } else {
1318           if(!fRawSP){
1319             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
1320             if(sz2 != 0) {
1321               ftime2 *= 0.169; //constants are from offline studies
1322               if(patch<2)
1323                 ftime2 *= 1.77;
1324             }
1325           } else ftime2=sz2; //take the width, not the error
1326         }
1327       }
1328       if(fStdout==kTRUE)
1329         HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
1330
1331       if(!fRawSP){
1332         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
1333         AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
1334
1335         if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
1336           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
1337         if(fNClusters >= fMaxNClusters)
1338           {
1339             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
1340               <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
1341             return;
1342           }
1343
1344         fSpacePointData[counter].fX = xyz[0];
1345         //      fSpacePointData[counter].fY = xyz[1];
1346         if(fCurrentSlice<18){
1347           fSpacePointData[counter].fY = xyz[1];
1348         }
1349         else{
1350           fSpacePointData[counter].fY = -1*xyz[1];
1351         }
1352         fSpacePointData[counter].fZ = xyz[2];
1353
1354       } else {
1355         fSpacePointData[counter].fX = fCurrentRow;
1356         fSpacePointData[counter].fY = fpad;
1357         fSpacePointData[counter].fZ = ftime;
1358       }
1359
1360       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
1361       fSpacePointData[counter].fPadRow = fCurrentRow;
1362       fSpacePointData[counter].fSigmaY2 = fpad2;
1363       fSpacePointData[counter].fSigmaZ2  = ftime2;
1364
1365       fSpacePointData[counter].fQMax = list[j].fQMax;
1366
1367       fSpacePointData[counter].fUsed = kFALSE;         // only used / set in AliHLTTPCDisplay
1368       fSpacePointData[counter].fTrackN = -1;           // only used / set in AliHLTTPCDisplay
1369
1370       Int_t patch=fCurrentPatch;
1371       if(patch==-1) patch=0; //never store negative patch number
1372       fSpacePointData[counter].fID = counter
1373         +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
1374
1375 #ifdef do_mc
1376       Int_t trackID[3];
1377       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
1378
1379       fSpacePointData[counter].fTrackID[0] = trackID[0];
1380       fSpacePointData[counter].fTrackID[1] = trackID[1];
1381       fSpacePointData[counter].fTrackID[2] = trackID[2];
1382
1383 #endif
1384
1385       fNClusters++;
1386       counter++;
1387     }
1388 }