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