- use_newio switch removed from libAliHLTTPC, always on; not removed in
[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  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
6  *                                                                        *
7  * Authors: Anders Vestbo <mailto:vestbo@fi.uib.no>,                      *
8  *          Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>    *
9  *          Jochen Thaeder <mailto:thaeder@kip.uni-heidelberg.de>         *
10  *          for The ALICE Off-line Project.                               *
11  *                                                                        *
12  * Permission to use, copy, modify and distribute this software and its   *
13  * documentation strictly for non-commercial purposes is hereby granted   *
14  * without fee, provided that the above copyright notice appears in all   *
15  * copies and that both the copyright notice and this permission notice   *
16  * appear in the supporting documentation. The authors make no claims     *
17  * about the suitability of this software for any purpose. It is          *
18  * provided "as is" without express or implied warranty.                  *
19  **************************************************************************/
20
21 /** @file   AliHLTTPCClusterFinder.cxx
22     @author Anders Vestbo <mailto:vestbo@fi.uib.no>,                 
23             Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
24             Jochen Thaeder <mailto:thaeder@kip.uni-heidelberg.de>     
25     @date   
26     @brief  Cluster Finder for the TPC
27 */
28
29 #include "AliHLTTPCDigitReader.h"
30 #include "AliHLTTPCRootTypes.h"
31 #include "AliHLTTPCLogging.h"
32 #include "AliHLTTPCClusterFinder.h"
33 #include "AliHLTTPCDigitData.h"
34 #include "AliHLTTPCTransform.h"
35 #include "AliHLTTPCSpacePointData.h"
36 #include "AliHLTTPCMemHandler.h"
37 #include "AliHLTTPCPad.h"
38
39 #if __GNUC__ >= 3
40 using namespace std;
41 #endif
42
43 /** \class AliHLTTPCClusterFinder
44 //
45 // The current cluster finder for HLT
46 // (Based on STAR L3)
47 // 
48 // The cluster finder is initialized with the Init function, 
49 // providing the slice and patch information to work on. 
50 //
51 // The input is a provided by the AliHLTTPCDigitReader class,
52 // using the init() funktion, and the next() funktion in order 
53 // to get the next bin. Either packed or unpacked data can be
54 // processed, dependent if one uses AliHLTTPCDigitReaderPacked 
55 // class or AliHLTTPCDigitReaderUnpacked class in the 
56 // Clusterfinder Component.
57 // The resulting space points will be in the
58 // array given by the SetOutputArray function.
59 // 
60 // There are several setters which control the behaviour:
61 //
62 // - SetXYError(Float_t):   set fixed error in XY direction
63 // - SetZError(Float_t):    set fixed error in Z  direction
64 //                            (used if errors are not calculated) 
65 // - SetDeconv(Bool_t):     switch on/off deconvolution
66 // - SetThreshold(UInt_t):  set charge threshold for cluster
67 // - SetMatchWidth(UInt_t): set the match distance in 
68 //                            time for sequences to be merged 
69 // - SetSTDOutput(Bool_t):  switch on/off output about found clusters   
70 // - SetCalcErr(Bool_t):    switch on/off calculation of 
71 //                          space point errors (or widths in raw system)
72 // - SetRawSP(Bool_t):      switch on/off convertion to raw system
73 //
74 //
75 // Example Usage:
76 //
77 // AliHLTTPCFileHandler *file = new AliHLTTPCFileHandler();
78 // file->SetAliInput(digitfile); //give some input file
79 // for(int slice=0; slice<=35; slice++){
80 //   for(int patch=0; pat<6; pat++){
81 //     file->Init(slice,patch);
82 //     UInt_t ndigits=0;
83 //     UInt_t maxclusters=100000;
84 //     UInt_t pointsize = maxclusters*sizeof(AliHLTTPCSpacePointData);
85 //     AliHLTTPCSpacePointData *points = (AliHLTTPCSpacePointData*)memory->Allocate(pointsize);
86 //     AliHLTTPCDigitRowData *digits = (AliHLTTPCDigitRowData*)file->AliAltroDigits2Memory(ndigits,event);
87 //     AliHLTTPCClusterFinder *cf = new AliHLTTPCClusterFinder();
88 //     cf->SetMatchWidth(2);
89 //     cf->InitSlice( slice, patch, row[0], row[1], maxPoints );
90 //     cf->SetSTDOutput(kTRUE);    //Some output to standard IO
91 //     cf->SetRawSP(kFALSE);       //Convert space points to local system
92 //     cf->SetThreshold(5);        //Threshold of cluster charge
93 //     cf->SetDeconv(kTRUE);       //Deconv in pad and time direction
94 //     cf->SetCalcErr(kTRUE);      //Calculate the errors of the spacepoints
95 //     cf->SetOutputArray(points); //Move the spacepoints to the array
96 //     cf->Read(iter->fPtr, iter->fSize ); //give the data to the cf
97 //     cf->ProcessDigits();        //process the rows given by init
98 //     Int_t npoints = cf->GetNumberOfClusters();
99 //     AliHLTTPCMemHandler *out= new AliHLTTPCMemHandler();
100 //     out->SetBinaryOutput(fname);
101 //     out->Memory2Binary(npoints,points); //store the spacepoints
102 //     out->CloseBinaryOutput();
103 //     delete out;
104 //     file->free();
105 //     delete cf;
106 //   }
107 // }
108 */
109
110 ClassImp(AliHLTTPCClusterFinder)
111
112 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
113   :
114   fMatch(1),
115   fThreshold(10),
116   fSignalThreshold(-1),
117   fOccupancyLimit(1.0),
118   fXYErr(0.2),
119   fZErr(0.3),
120   fDeconvPad(kTRUE),
121   fDeconvTime(kTRUE),
122   fStdout(kFALSE),
123   fCalcerr(kTRUE),
124   fRawSP(kFALSE),
125   fFirstRow(0),
126   fLastRow(0),
127   fDigitReader(NULL)
128 {
129   //constructor
130 }
131
132 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder(const AliHLTTPCClusterFinder& src)
133   :
134   fMatch(src.fMatch),
135   fThreshold(src.fThreshold),
136   fSignalThreshold(src.fSignalThreshold),
137   fOccupancyLimit(src.fOccupancyLimit),
138   fXYErr(src.fXYErr),
139   fZErr(src.fZErr),
140   fDeconvPad(src.fDeconvPad),
141   fDeconvTime(src.fDeconvTime),
142   fStdout(src.fStdout),
143   fCalcerr(src.fCalcerr),
144   fRawSP(src.fRawSP),
145   fFirstRow(src.fFirstRow),
146   fLastRow(src.fLastRow),
147   fDigitReader(src.fDigitReader)
148 {
149 }
150
151 AliHLTTPCClusterFinder& AliHLTTPCClusterFinder::operator=(const AliHLTTPCClusterFinder& src)
152 {
153   fMatch=src.fMatch;
154   fThreshold=src.fThreshold;
155   fSignalThreshold=src.fSignalThreshold;
156   fOccupancyLimit=src.fOccupancyLimit;
157   fXYErr=src.fXYErr;
158   fZErr=src.fZErr;
159   fDeconvPad=src.fDeconvPad;
160   fDeconvTime=src.fDeconvTime;
161   fStdout=src.fStdout;
162   fCalcerr=src.fCalcerr;
163   fRawSP=src.fRawSP;
164   fFirstRow=src.fFirstRow;
165   fLastRow=src.fLastRow;
166   fDigitReader=src.fDigitReader;
167   return (*this);
168 }
169
170 AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder()
171 {
172   //destructor
173 }
174  
175 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
176 {
177   //init slice
178   fNClusters = 0;
179   fMaxNClusters = nmaxpoints;
180   fCurrentSlice = slice;
181   fCurrentPatch = patch;
182   fFirstRow = firstrow;
183   fLastRow = lastrow;
184 }
185
186 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
187 {
188   //init slice
189   fNClusters = 0;
190   fMaxNClusters = nmaxpoints;
191   fCurrentSlice = slice;
192   fCurrentPatch = patch;
193   fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
194   fLastRow=AliHLTTPCTransform::GetLastRow(patch);
195 }
196
197 void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt)
198 {
199   //set pointer to output
200   fSpacePointData = pt;
201 }
202
203 void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
204   //set input pointer
205   fPtr = (UChar_t*)ptr;
206   fSize = size;
207 }
208
209 void AliHLTTPCClusterFinder::ProcessDigits()
210 {
211   bool readValue = true;
212   Int_t newRow = 0;    
213   Int_t rowOffset = 0;
214   UShort_t time=0,newTime=0;
215   UInt_t pad=0,newPad=0;
216   AliHLTTPCSignal_t charge=0;
217
218   fNClusters = 0;
219
220   // initialize block for reading packed data
221   fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
222   readValue = fDigitReader->Next();
223
224   // Matthias 08.11.2006 the following return would cause termination without writing the
225   // ClusterData and thus would block the component. I just want to have the commented line
226   // here for information
227   //if (!readValue)return;
228
229   pad = fDigitReader->GetPad();
230   time = fDigitReader->GetTime();
231   fCurrentRow = fDigitReader->GetRow();
232
233   if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
234     rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
235
236   fCurrentRow += rowOffset;
237
238   UInt_t lastpad = 123456789;
239   const Int_t kPadArraySize=5000;
240   const Int_t kClusterListSize=10000;
241   AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
242   AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
243   AliClusterData clusterlist[kClusterListSize]; //Clusterlist
244
245   AliClusterData **currentPt;  //List of pointers to the current pad
246   AliClusterData **previousPt; //List of pointers to the previous pad
247   currentPt = pad2;
248   previousPt = pad1;
249   UInt_t nprevious=0,ncurrent=0,ntotal=0;
250
251   /* quick implementation of baseline calculation and zero suppression
252      open a pad object for each pad and delete it after processing.
253      later a list of pad objects with base line history can be used
254      The whole thing only works if we really get unprocessed raw data, if
255      the data is already zero suppressed, there might be gaps in the time
256      bins.
257    */
258   Int_t gatingGridOffset=50;
259   AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
260   // just to make later conversion to a list of objects easier
261   AliHLTTPCPad* pCurrentPad=NULL;
262   if (fSignalThreshold>=0) {
263     pCurrentPad=&baseline;
264     baseline.SetThreshold(fSignalThreshold);
265   }
266
267   while ( readValue ){   // Reads through all digits in block
268
269     if(pad != lastpad){
270       //This is a new pad
271       
272       //Switch the lists:
273       if(currentPt == pad2){
274         currentPt = pad1;
275         previousPt = pad2;
276       }
277       else {
278         currentPt = pad2;
279         previousPt = pad1;
280       }
281       nprevious = ncurrent;
282       ncurrent = 0;
283       if(pad != lastpad+1){
284         //this happens if there is a pad with no signal.
285         nprevious = ncurrent = 0;
286       }
287       lastpad = pad;
288     }
289
290     Bool_t newcluster = kTRUE;
291     UInt_t seqcharge=0,seqaverage=0,seqerror=0;
292     UInt_t lastcharge=0,lastwas_falling=0;
293     Int_t newbin=-1;
294
295
296     if(fDeconvTime){
297       redo: //This is a goto.
298       
299       if(newbin > -1){
300         //bin = newbin;
301         newbin = -1;
302       }
303           
304       lastcharge=0;
305       lastwas_falling = 0;
306     }
307
308     while(1){ //Loop over time bins of current pad
309       // read all the values for one pad at once to calculate the base line
310       if (pCurrentPad) {
311         if (!pCurrentPad->IsStarted()) {
312           //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
313           pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
314           if ((pCurrentPad->StartEvent())>=0) {
315             do {
316               if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
317               if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
318               pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
319               //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
320             } while ((readValue = fDigitReader->Next())!=0);
321           }
322           pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
323           if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
324             //HLTDebug("no data available after zero suppression");
325             pCurrentPad->StopEvent();
326             pCurrentPad->ResetHistory();
327             break;
328           }
329           time=pCurrentPad->GetCurrentPosition();
330           if (time>pCurrentPad->GetSize()) {
331             HLTError("invalid time bin for pad");
332             break;
333           }
334         }
335       }
336
337       if (pCurrentPad) {
338         Float_t occupancy=pCurrentPad->GetOccupancy();
339         //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
340         if ( occupancy < fOccupancyLimit ) {
341           charge = pCurrentPad->GetCorrectedData();
342         } else {
343           charge = 0;
344           //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
345         }
346       } else {
347         charge = fDigitReader->GetSignal();
348       }
349       //HLTDebug("get next charge value: position %d charge %d", time, charge);
350
351
352       // CHARGE DEBUG
353       if (fDigitReader->GetRow() == 90){
354 /////     LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90")  << "PAD=" <<  fDigitReader->GetPad() << "  TIME=" <<  fDigitReader->GetTime() 
355           //                                       << "  SIGNAL=" <<  fDigitReader->GetSignal() << ENDLOG;
356
357       }
358
359       if(time >= AliHLTTPCTransform::GetNTimeBins()){
360         HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
361         break;
362       }
363
364
365       //Get the current ADC-value
366       if(fDeconvTime){
367
368         //Check if the last pixel in the sequence is smaller than this
369         if(charge > lastcharge){
370           if(lastwas_falling){
371             newbin = 1;
372             break;
373           }
374         }
375         else lastwas_falling = 1; //last pixel was larger than this
376         lastcharge = charge;
377       }
378           
379       //Sum the total charge of this sequence
380       seqcharge += charge;
381       seqaverage += time*charge;
382       seqerror += time*time*charge;
383       
384       if (pCurrentPad) {
385         
386         if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
387           pCurrentPad->StopEvent();
388           pCurrentPad->ResetHistory();
389           if(readValue) {
390             newPad = fDigitReader->GetPad();
391             newTime = fDigitReader->GetTime();
392             newRow = fDigitReader->GetRow() + rowOffset;
393           }
394           break;
395         }
396
397         newPad=pCurrentPad->GetPadNumber();
398         newTime=pCurrentPad->GetCurrentPosition();
399         newRow=pCurrentPad->GetRowNumber();
400       } else {
401       readValue = fDigitReader->Next();
402       //Check where to stop:
403       if(!readValue) break; //No more value
404
405       newPad = fDigitReader->GetPad();
406       newTime = fDigitReader->GetTime();
407       newRow = fDigitReader->GetRow() + rowOffset;
408       }
409
410       if(newPad != pad)break; //new pad
411       if(newTime != time+1) break; //end of sequence
412
413       // pad = newpad;    is equal
414       time = newTime;
415
416     }//end loop over sequence
417
418     //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
419     //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
420     if (seqcharge<=0) {
421       // with active zero suppression zero values are possible
422       continue;
423     }
424
425     //Calculate mean of sequence:
426     Int_t seqmean=0;
427     if(seqcharge)
428       seqmean = seqaverage/seqcharge;
429     else{
430       LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
431         <<"Error in data given to the cluster finder"<<ENDLOG;
432       seqmean = 1;
433       seqcharge = 1;
434     }
435
436     //Calculate mean in pad direction:
437     Int_t padmean = seqcharge*pad;
438     Int_t paderror = pad*padmean;
439
440
441     //Compare with results on previous pad:
442     for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
443       
444       //dont merge sequences on the same pad twice
445       if(previousPt[p]->fLastMergedPad==pad) continue;
446
447       Int_t difference = seqmean - previousPt[p]->fMean;
448       if(difference < -fMatch) break;
449
450       if(difference <= fMatch){ //There is a match here!!
451         AliClusterData *local = previousPt[p];
452         
453         if(fDeconvPad){
454           if(seqcharge > local->fLastCharge){
455             if(local->fChargeFalling){ //The previous pad was falling
456               break; //create a new cluster
457             }               
458           }
459           else local->fChargeFalling = 1;
460           local->fLastCharge = seqcharge;
461         }
462               
463         //Don't create a new cluster, because we found a match
464         newcluster = kFALSE;
465               
466         //Update cluster on current pad with the matching one:
467         local->fTotalCharge += seqcharge;
468         local->fPad += padmean;
469         local->fPad2 += paderror;
470         local->fTime += seqaverage;
471         local->fTime2 += seqerror;
472         local->fMean = seqmean;
473         local->fFlags++; //means we have more than one pad 
474         local->fLastMergedPad = pad;
475
476         currentPt[ncurrent] = local;
477         ncurrent++;
478               
479         break;
480       } //Checking for match at previous pad
481     } //Loop over results on previous pad.
482
483
484     if(newcluster && ncurrent<kPadArraySize){
485       //Start a new cluster. Add it to the clusterlist, and update
486       //the list of pointers to clusters in current pad.
487       //current pad will be previous pad on next pad.
488
489       //Add to the clusterlist:
490       AliClusterData *tmp = &clusterlist[ntotal];
491       tmp->fTotalCharge = seqcharge;
492       tmp->fPad = padmean;
493       tmp->fPad2 = paderror;
494       tmp->fTime = seqaverage;
495       tmp->fTime2 = seqerror;
496       tmp->fMean = seqmean;
497       tmp->fFlags = 0;  //flags for single pad clusters
498       tmp->fLastMergedPad = pad;
499
500       if(fDeconvPad){
501         tmp->fChargeFalling = 0;
502         tmp->fLastCharge = seqcharge;
503       }
504
505       //Update list of pointers to previous pad:
506       currentPt[ncurrent] = &clusterlist[ntotal];
507       ntotal++;
508       ncurrent++;
509     }
510
511     if(fDeconvTime)
512       if(newbin >= 0) goto redo;
513   
514     // to prevent endless loop  
515     if(time >= AliHLTTPCTransform::GetNTimeBins()){
516       HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
517       break;
518     }
519
520
521     if(!readValue) break; //No more value
522     
523     if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
524       HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
525       break;
526     }
527
528     if(fCurrentRow != newRow){
529       WriteClusters(ntotal,clusterlist);
530
531       lastpad = 123456789;
532
533       currentPt = pad2;
534       previousPt = pad1;
535       nprevious=0;
536       ncurrent=0;
537       ntotal=0;
538       
539       fCurrentRow = newRow;
540     }
541
542     pad = newPad;
543     time = newTime;
544
545   } // END while(readValue)
546
547   WriteClusters(ntotal,clusterlist);
548
549   HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
550
551 } // ENDEND
552
553 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list)
554 {
555   //write cluster to output pointer
556   Int_t thisrow,thissector;
557   UInt_t counter = fNClusters;
558   
559   for(int j=0; j<nclusters; j++)
560     {
561       if(!list[j].fFlags) continue; //discard single pad clusters
562       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
563
564       Float_t xyz[3];      
565       Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
566       Float_t fpad2=fXYErr*fXYErr; //fixed given error
567       Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
568       Float_t ftime2=fZErr*fZErr;  //fixed given error
569
570
571    
572      
573
574
575       if(fCalcerr) { //calc the errors, otherwice take the fixed error 
576         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
577         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
578         Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
579         sy2/=q2;
580         if(sy2 < 0) {
581             LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
582               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
583             continue;
584         } else {
585           if(!fRawSP){
586             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
587             if(sy2 != 0){
588               fpad2*=0.108; //constants are from offline studies
589               if(patch<2)
590                 fpad2*=2.07;
591             }
592           } else fpad2=sy2; //take the width not the error
593         }
594         Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
595         sz2/=q2;
596         if(sz2 < 0){
597           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
598             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
599           continue;
600         } else {
601           if(!fRawSP){
602             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
603             if(sz2 != 0) {
604               ftime2 *= 0.169; //constants are from offline studies
605               if(patch<2)
606                 ftime2 *= 1.77;
607             }
608           } else ftime2=sz2; //take the width, not the error
609         }
610       }
611       if(fStdout==kTRUE)
612         cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
613       
614       if(!fRawSP){
615         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
616         AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
617         
618         if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
619           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
620         if(fNClusters >= fMaxNClusters)
621           {
622             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
623               <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
624             return;
625           }  
626         
627         fSpacePointData[counter].fX = xyz[0];
628         fSpacePointData[counter].fY = xyz[1];
629         fSpacePointData[counter].fZ = xyz[2];
630         
631       } else {
632         fSpacePointData[counter].fX = fCurrentRow;
633         fSpacePointData[counter].fY = fpad;
634         fSpacePointData[counter].fZ = ftime;
635       }
636       
637       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
638       fSpacePointData[counter].fPadRow = fCurrentRow;
639       fSpacePointData[counter].fSigmaY2 = fpad2;
640       fSpacePointData[counter].fSigmaZ2  = ftime2;
641
642       fSpacePointData[counter].fUsed = kFALSE;         // only used / set in AliHLTTPCDisplay
643       fSpacePointData[counter].fTrackN = -1;           // only used / set in AliHLTTPCDisplay
644
645       Int_t patch=fCurrentPatch;
646       if(patch==-1) patch=0; //never store negative patch number
647       fSpacePointData[counter].fID = counter
648         +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
649
650 #ifdef do_mc
651       Int_t trackID[3];
652       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
653
654       fSpacePointData[counter].fTrackID[0] = trackID[0];
655       fSpacePointData[counter].fTrackID[1] = trackID[1];
656       fSpacePointData[counter].fTrackID[2] = trackID[2];
657
658       //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
659 #endif
660       
661       fNClusters++;
662       counter++;
663     }
664 }
665
666 // STILL TO FIX  ----------------------------------------------------------------------------
667
668 #ifdef do_mc
669 void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
670 {
671   //get mc id
672   AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
673   
674   trackID[0]=trackID[1]=trackID[2]=-2;
675   //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
676   for(Int_t i=fFirstRow; i<=fLastRow; i++){
677     if(rowPt->fRow < (UInt_t)fCurrentRow){
678       AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
679       continue;
680     }
681     AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
682     for(UInt_t j=0; j<rowPt->fNDigit; j++){
683       Int_t cpad = digPt[j].fPad;
684       Int_t ctime = digPt[j].fTime;
685       if(cpad != pad) continue;
686       if(ctime != time) continue;
687
688       trackID[0] = digPt[j].fTrackID[0];
689       trackID[1] = digPt[j].fTrackID[1];
690       trackID[2] = digPt[j].fTrackID[2];
691       
692       //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
693       break;
694     }
695     break;
696   }
697 }
698 #endif