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