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