HLT/TPCLib/AliHLTTPCClusterFinder.cxx

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