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