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()-1;
 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
 356                 // Checks if one need to deconvolute the signals
 357                 if(bunchData[i]>prevSignal && signalFalling==kTRUE){
 358                   if(i<fDigitReader->GetBunchSize()-1){ // means there are more than one signal left in the bunch
 359                     moreDataInBunch=kTRUE;
 360                     prevSignal=0;
 361                   }
 362                   break;
 363                 }
 364
 365                 // Checks if the signal is 0, then quit processing the data.
 366                 if(bunchData[i]==0 && i<fDigitReader->GetBunchSize()-1){//means we have 0 data fom the rcu, might happen depending on the configuration settings
 367                   moreDataInBunch=kTRUE;
 368                   prevSignal=0;
 369                   break;
 370                 }
 371
 372                 if(prevSignal>bunchData[i]){//means the peak of the signal has been reached and deconvolution will happen if the signal rise again.
 373                   signalFalling=kTRUE;
 374                 }
 375                 candidate.fTotalCharge+=bunchData[i];
 376                 candidate.fTime += time*bunchData[i];
 377                 candidate.fTime2 += time*time*bunchData[i];
 378                 if(bunchData[i]>candidate.fQMax){
 379                   candidate.fQMax=bunchData[i];
 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               if(indexInBunchData<fDigitReader->GetBunchSize()-1){
 394                 moreDataInBunch=kFALSE;
 395               }
 396             }while(moreDataInBunch);
 397           }
 398           else{
 399             const UInt_t *bunchData= fDigitReader->GetSignals();
 400             do{
 401               AliHLTTPCClusters candidate;
 402               const AliHLTTPCDigitData* digits = fDigitReader->GetBunchDigits();
 403               if(fDoMC) fMCDigits.clear();
 404
 405               for(Int_t i=indexInBunchData;i<fDigitReader->GetBunchSize();i++){
 406                 // Checks if one need to deconvolute the signals
 407                 if(bunchData[i]>prevSignal && signalFalling==kTRUE){
 408                   if(i<fDigitReader->GetBunchSize()-1){ // means there are more than one signal left in the bunch
 409                     moreDataInBunch=kTRUE;
 410                     prevSignal=0;
 411                   }
 412                   break;
 413                 }
 414
 415                 // Checks if the signal is 0, then quit processing the data.
 416                 if(bunchData[i]==0 && i<fDigitReader->GetBunchSize()-1){//means we have 0 data fom the rcu, might happen depending on the configuration settings
 417                   moreDataInBunch=kTRUE;
 418                   prevSignal=0;
 419                   break;
 420                 }
 421
 422                 if(prevSignal>bunchData[i]){//means the peak of the signal has been reached and deconvolution will happen if the signal rise again.
 423                   signalFalling=kTRUE;
 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(AliHLTTPCClusterMCLabel * 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(), CompareWeights );
 588             AliHLTTPCClusterMCLabel tmpClusterMCInfo;
 589
 590             AliHLTTPCClusterMCWeight 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         AliHLTTPCClusterMCWeight 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         if(q2 == 0) {
1131           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1132             <<"zero charge "<< list[j].fTotalCharge <<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1133           continue;
1134         }
1135         sy2/=q2;
1136         if(sy2 < 0) {
1137             LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1138               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1139             continue;
1140         } else {
1141           if(!fRawSP){
1142             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
1143             if(sy2 != 0){
1144               fpad2*=0.108; //constants are from offline studies
1145               if(patch<2)
1146                 fpad2*=2.07;
1147             }
1148           } else fpad2=sy2; //take the width not the error
1149         }
1150         //      Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
1151         Float_t sz2=(Float_t)list[j].fTime2*list[j].fTotalCharge - (Float_t)list[j].fTime*list[j].fTime;
1152         sz2/=q2;
1153         if(sz2 < 0){
1154           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1155             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1156           continue;
1157         } else {
1158           if(!fRawSP){
1159             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
1160             if(sz2 != 0) {
1161               ftime2 *= 0.169; //constants are from offline studies
1162               if(patch<2)
1163                 ftime2 *= 1.77;
1164             }
1165           } else ftime2=sz2; //take the width, not the error
1166         }
1167       }
1168       if(fStdout==kTRUE)
1169         HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
1170
1171       if(!fRawSP){
1172         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
1173
1174         if(fOfflineTransform == NULL){
1175           AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
1176
1177           if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
1178                           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
1179           if(fNClusters >= fMaxNClusters)
1180             {
1181               LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
1182                 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
1183               return;
1184             }
1185
1186           fSpacePointData[counter].fX = xyz[0];
1187           //    fSpacePointData[counter].fY = xyz[1];
1188           if(fCurrentSlice<18){
1189             fSpacePointData[counter].fY = xyz[1];
1190           }
1191           else{
1192             fSpacePointData[counter].fY = -1*xyz[1];
1193           }
1194           fSpacePointData[counter].fZ = xyz[2];
1195         }
1196         else{
1197           Double_t x[3]={thisrow,fpad+.5,ftime};
1198           Int_t iSector[1]={thissector};
1199           fOfflineTransform->Transform(x,iSector,0,1);
1200           double y[3] = {x[0], x[1], x[2] };
1201
1202           if( fOfflineTPCParam && thissector<fOfflineTPCParam->GetNSector() ){
1203             TGeoHMatrix  *alignment = fOfflineTPCParam->GetClusterMatrix( thissector );
1204             if ( alignment ) alignment->LocalToMaster( x, y);
1205           }
1206
1207           fSpacePointData[counter].fX = y[0];
1208           fSpacePointData[counter].fY = y[1];
1209           fSpacePointData[counter].fZ = y[2];
1210         }
1211
1212       }
1213       else {
1214         fSpacePointData[counter].fX = fCurrentRow;
1215         fSpacePointData[counter].fY = fpad;
1216         fSpacePointData[counter].fZ = ftime;
1217       }
1218
1219       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
1220       fSpacePointData[counter].fPadRow = fCurrentRow;
1221       fSpacePointData[counter].fSigmaY2 = fpad2;
1222       fSpacePointData[counter].fSigmaZ2  = ftime2;
1223
1224       fSpacePointData[counter].fQMax = list[j].fQMax;
1225
1226       fSpacePointData[counter].SetUsed(kFALSE);         // only used / set in AliHLTTPCDisplay
1227       fSpacePointData[counter].SetTrackNumber(-1);      // only used / set in AliHLTTPCDisplay
1228
1229       Int_t patch=fCurrentPatch;
1230       if(patch==-1) patch=0; //never store negative patch number
1231       fSpacePointData[counter].SetID( fCurrentSlice, patch, counter );
1232
1233 #ifdef do_mc
1234       Int_t trackID[3];
1235       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
1236
1237       fSpacePointData[counter].fTrackID[0] = trackID[0];
1238       fSpacePointData[counter].fTrackID[1] = trackID[1];
1239       fSpacePointData[counter].fTrackID[2] = trackID[2];
1240
1241 #endif
1242
1243       fNClusters++;
1244       counter++;
1245     }
1246 }
1247
1248 // STILL TO FIX  ----------------------------------------------------------------------------
1249
1250 #ifdef do_mc
1251 void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID) const {
1252   // see header file for class documentation
1253
1254   //get mc id
1255   AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
1256
1257   trackID[0]=trackID[1]=trackID[2]=-2;
1258   for(Int_t i=fFirstRow; i<=fLastRow; i++){
1259     if(rowPt->fRow < (UInt_t)fCurrentRow){
1260       AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
1261       continue;
1262     }
1263     AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
1264     for(UInt_t j=0; j<rowPt->fNDigit; j++){
1265       Int_t cpad = digPt[j].fPad;
1266       Int_t ctime = digPt[j].fTime;
1267       if(cpad != pad) continue;
1268       if(ctime != time) continue;
1269
1270       trackID[0] = digPt[j].fTrackID[0];
1271       trackID[1] = digPt[j].fTrackID[1];
1272       trackID[2] = digPt[j].fTrackID[2];
1273
1274       break;
1275     }
1276     break;
1277   }
1278 }
1279 #endif
1280
1281
1282
1283 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list){//This is used when using the AliHLTTPCClusters class for cluster data
1284   // see header file for class documentation
1285
1286   //write cluster to output pointer
1287   Int_t thisrow,thissector;
1288   UInt_t counter = fNClusters;
1289
1290   for(int j=0; j<nclusters; j++)
1291     {
1292       if(!list[j].fFlags) continue; //discard single pad clusters
1293       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
1294
1295       Float_t xyz[3];
1296       Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
1297       Float_t fpad2=fXYErr*fXYErr; //fixed given error
1298       Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
1299       Float_t ftime2=fZErr*fZErr;  //fixed given error
1300
1301
1302       if(fCalcerr) { //calc the errors, otherwice take the fixed error
1303         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
1304         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
1305         Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
1306         sy2/=q2;
1307         if(sy2 < 0) {
1308             LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1309               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1310             continue;
1311         } else {
1312           if(!fRawSP){
1313             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
1314             if(sy2 != 0){
1315               fpad2*=0.108; //constants are from offline studies
1316               if(patch<2)
1317                 fpad2*=2.07;
1318             }
1319           } else fpad2=sy2; //take the width not the error
1320         }
1321         Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
1322         sz2/=q2;
1323         if(sz2 < 0){
1324           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1325             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1326           continue;
1327         } else {
1328           if(!fRawSP){
1329             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
1330             if(sz2 != 0) {
1331               ftime2 *= 0.169; //constants are from offline studies
1332               if(patch<2)
1333                 ftime2 *= 1.77;
1334             }
1335           } else ftime2=sz2; //take the width, not the error
1336         }
1337       }
1338       if(fStdout==kTRUE)
1339         HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
1340
1341       if(!fRawSP){
1342         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
1343         AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
1344
1345         if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
1346           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
1347         if(fNClusters >= fMaxNClusters)
1348           {
1349             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
1350               <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
1351             return;
1352           }
1353
1354         fSpacePointData[counter].fX = xyz[0];
1355         //      fSpacePointData[counter].fY = xyz[1];
1356         if(fCurrentSlice<18){
1357           fSpacePointData[counter].fY = xyz[1];
1358         }
1359         else{
1360           fSpacePointData[counter].fY = -1*xyz[1];
1361         }
1362         fSpacePointData[counter].fZ = xyz[2];
1363
1364       } else {
1365         fSpacePointData[counter].fX = fCurrentRow;
1366         fSpacePointData[counter].fY = fpad;
1367         fSpacePointData[counter].fZ = ftime;
1368       }
1369
1370       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
1371       fSpacePointData[counter].fPadRow = fCurrentRow;
1372       fSpacePointData[counter].fSigmaY2 = fpad2;
1373       fSpacePointData[counter].fSigmaZ2  = ftime2;
1374
1375       fSpacePointData[counter].fQMax = list[j].fQMax;
1376
1377       fSpacePointData[counter].SetUsed(kFALSE);         // only used / set in AliHLTTPCDisplay
1378       fSpacePointData[counter].SetTrackNumber(-1);      // only used / set in AliHLTTPCDisplay
1379
1380       Int_t patch=fCurrentPatch;
1381       if(patch==-1) patch=0; //never store negative patch number
1382       fSpacePointData[counter].SetID( fCurrentSlice, patch, counter );
1383
1384 #ifdef do_mc
1385       Int_t trackID[3];
1386       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
1387
1388       fSpacePointData[counter].fTrackID[0] = trackID[0];
1389       fSpacePointData[counter].fTrackID[1] = trackID[1];
1390       fSpacePointData[counter].fTrackID[2] = trackID[2];
1391
1392 #endif
1393
1394       fNClusters++;
1395       counter++;
1396     }
1397 }