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