- configure adapted to the new directory structure of the HOMER module in PubSub
[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   fOccupancyLimit(1.0),
104   fXYErr(0.2),
105   fZErr(0.3),
106   fDeconvPad(kTRUE),
107   fDeconvTime(kTRUE),
108   fStdout(kFALSE),
109   fCalcerr(kTRUE),
110   fRawSP(kFALSE),
111   fFirstRow(0),
112   fLastRow(0),
113   fDigitReader(NULL)
114 {
115   //constructor
116 }
117
118 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder(const AliHLTTPCClusterFinder& src)
119   :
120   fMatch(src.fMatch),
121   fThreshold(src.fThreshold),
122   fSignalThreshold(src.fSignalThreshold),
123   fOccupancyLimit(src.fOccupancyLimit),
124   fXYErr(src.fXYErr),
125   fZErr(src.fZErr),
126   fDeconvPad(src.fDeconvPad),
127   fDeconvTime(src.fDeconvTime),
128   fStdout(src.fStdout),
129   fCalcerr(src.fCalcerr),
130   fRawSP(src.fRawSP),
131   fFirstRow(src.fFirstRow),
132   fLastRow(src.fLastRow),
133   fDigitReader(src.fDigitReader)
134 {
135 }
136
137 AliHLTTPCClusterFinder& AliHLTTPCClusterFinder::operator=(const AliHLTTPCClusterFinder& src)
138 {
139   fMatch=src.fMatch;
140   fThreshold=src.fThreshold;
141   fSignalThreshold=src.fSignalThreshold;
142   fOccupancyLimit=src.fOccupancyLimit;
143   fXYErr=src.fXYErr;
144   fZErr=src.fZErr;
145   fDeconvPad=src.fDeconvPad;
146   fDeconvTime=src.fDeconvTime;
147   fStdout=src.fStdout;
148   fCalcerr=src.fCalcerr;
149   fRawSP=src.fRawSP;
150   fFirstRow=src.fFirstRow;
151   fLastRow=src.fLastRow;
152   fDigitReader=src.fDigitReader;
153   return (*this);
154 }
155
156 AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder()
157 {
158   //destructor
159 }
160  
161 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
162 {
163   //init slice
164   fNClusters = 0;
165   fMaxNClusters = nmaxpoints;
166   fCurrentSlice = slice;
167   fCurrentPatch = patch;
168   fFirstRow = firstrow;
169   fLastRow = lastrow;
170 }
171
172 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
173 {
174   //init slice
175   fNClusters = 0;
176   fMaxNClusters = nmaxpoints;
177   fCurrentSlice = slice;
178   fCurrentPatch = patch;
179   fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
180   fLastRow=AliHLTTPCTransform::GetLastRow(patch);
181 }
182
183 void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt)
184 {
185   //set pointer to output
186   fSpacePointData = pt;
187 }
188
189 void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
190   //set input pointer
191   fPtr = (UChar_t*)ptr;
192   fSize = size;
193 }
194
195 void AliHLTTPCClusterFinder::ProcessDigits()
196 {
197   bool readValue = true;
198   Int_t newRow = 0;    
199   Int_t rowOffset = 0;
200   UShort_t time=0,newTime=0;
201   UInt_t pad=0,newPad=0;
202   AliHLTTPCSignal_t charge=0;
203
204   fNClusters = 0;
205
206   // initialize block for reading packed data
207   fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
208   readValue = fDigitReader->Next();
209
210   // Matthias 08.11.2006 the following return would cause termination without writing the
211   // ClusterData and thus would block the component. I just wnt to have the commented line
212   // here for information
213   //if (!readValue)return;
214
215   pad = fDigitReader->GetPad();
216   time = fDigitReader->GetTime();
217   fCurrentRow = fDigitReader->GetRow();
218
219   if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
220     rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
221
222   fCurrentRow += rowOffset;
223
224   UInt_t lastpad = 123456789;
225   const Int_t kPadArraySize=5000;
226   const Int_t kClusterListSize=10000;
227   AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
228   AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
229   AliClusterData clusterlist[kClusterListSize]; //Clusterlist
230
231   AliClusterData **currentPt;  //List of pointers to the current pad
232   AliClusterData **previousPt; //List of pointers to the previous pad
233   currentPt = pad2;
234   previousPt = pad1;
235   UInt_t nprevious=0,ncurrent=0,ntotal=0;
236
237   /* quick implementation of baseline calculation and zero suppression
238      open a pad object for each pad and delete it after processing.
239      later a list of pad objects with base line history can be used
240      The whole thing only works if we really get unprocessed raw data, if
241      the data is already zero suppressed, there might be gaps in the time
242      bins.
243    */
244   Int_t gatingGridOffset=50;
245   AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
246   // just to make later conversion to a list of objects easier
247   AliHLTTPCPad* pCurrentPad=NULL;
248   if (fSignalThreshold>=0) {
249     pCurrentPad=&baseline;
250     baseline.SetThreshold(fSignalThreshold);
251   }
252
253   while ( readValue ){   // Reads through all digits in block
254
255     if(pad != lastpad){
256       //This is a new pad
257       
258       //Switch the lists:
259       if(currentPt == pad2){
260         currentPt = pad1;
261         previousPt = pad2;
262       }
263       else {
264         currentPt = pad2;
265         previousPt = pad1;
266       }
267       nprevious = ncurrent;
268       ncurrent = 0;
269       if(pad != lastpad+1){
270         //this happens if there is a pad with no signal.
271         nprevious = ncurrent = 0;
272       }
273       lastpad = pad;
274     }
275
276     Bool_t newcluster = kTRUE;
277     UInt_t seqcharge=0,seqaverage=0,seqerror=0;
278     UInt_t lastcharge=0,lastwas_falling=0;
279     Int_t newbin=-1;
280
281
282     if(fDeconvTime){
283       redo: //This is a goto.
284       
285       if(newbin > -1){
286         //bin = newbin;
287         newbin = -1;
288       }
289           
290       lastcharge=0;
291       lastwas_falling = 0;
292     }
293
294     while(1){ //Loop over time bins of current pad
295       // read all the values for one pad at once to calculate the base line
296       if (pCurrentPad) {
297         if (!pCurrentPad->IsStarted()) {
298           //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
299           pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
300           if ((pCurrentPad->StartEvent())>=0) {
301             do {
302               if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
303               if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
304               pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
305               //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
306             } while ((readValue = fDigitReader->Next())!=0);
307           }
308           pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
309           if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
310             //HLTDebug("no data available after zero suppression");
311             pCurrentPad->StopEvent();
312             pCurrentPad->ResetHistory();
313             break;
314           }
315           time=pCurrentPad->GetCurrentPosition();
316           if (time>pCurrentPad->GetSize()) {
317             HLTError("invalid time bin for pad");
318             break;
319           }
320         }
321       }
322
323       if (pCurrentPad) {
324         Float_t occupancy=pCurrentPad->GetOccupancy();
325         //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
326         if ( occupancy < fOccupancyLimit ) {
327           charge = pCurrentPad->GetCorrectedData();
328         } else {
329           charge = 0;
330           //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
331         }
332       } else {
333         charge = fDigitReader->GetSignal();
334       }
335       //HLTDebug("get next charge value: position %d charge %d", time, charge);
336
337
338       // CHARGE DEBUG
339       if (fDigitReader->GetRow() == 90){
340 /////     LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90")  << "PAD=" <<  fDigitReader->GetPad() << "  TIME=" <<  fDigitReader->GetTime() 
341           //                                       << "  SIGNAL=" <<  fDigitReader->GetSignal() << ENDLOG;
342
343       }
344
345       if(time >= AliHLTTPCTransform::GetNTimeBins()){
346         HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
347         break;
348       }
349
350
351       //Get the current ADC-value
352       if(fDeconvTime){
353
354         //Check if the last pixel in the sequence is smaller than this
355         if(charge > lastcharge){
356           if(lastwas_falling){
357             newbin = 1;
358             break;
359           }
360         }
361         else lastwas_falling = 1; //last pixel was larger than this
362         lastcharge = charge;
363       }
364           
365       //Sum the total charge of this sequence
366       seqcharge += charge;
367       seqaverage += time*charge;
368       seqerror += time*time*charge;
369       
370       if (pCurrentPad) {
371         
372         if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
373           pCurrentPad->StopEvent();
374           pCurrentPad->ResetHistory();
375           if(readValue) {
376             newPad = fDigitReader->GetPad();
377             newTime = fDigitReader->GetTime();
378             newRow = fDigitReader->GetRow() + rowOffset;
379           }
380           break;
381         }
382
383         newPad=pCurrentPad->GetPadNumber();
384         newTime=pCurrentPad->GetCurrentPosition();
385         newRow=pCurrentPad->GetRowNumber();
386       } else {
387       readValue = fDigitReader->Next();
388       //Check where to stop:
389       if(!readValue) break; //No more value
390
391       newPad = fDigitReader->GetPad();
392       newTime = fDigitReader->GetTime();
393       newRow = fDigitReader->GetRow() + rowOffset;
394       }
395
396       if(newPad != pad)break; //new pad
397       if(newTime != time+1) break; //end of sequence
398
399       // pad = newpad;    is equal
400       time = newTime;
401
402     }//end loop over sequence
403
404     //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
405     //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
406     if (seqcharge<=0) {
407       // with active zero suppression zero values are possible
408       continue;
409     }
410
411     //Calculate mean of sequence:
412     Int_t seqmean=0;
413     if(seqcharge)
414       seqmean = seqaverage/seqcharge;
415     else{
416       LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
417         <<"Error in data given to the cluster finder"<<ENDLOG;
418       seqmean = 1;
419       seqcharge = 1;
420     }
421
422     //Calculate mean in pad direction:
423     Int_t padmean = seqcharge*pad;
424     Int_t paderror = pad*padmean;
425
426
427     //Compare with results on previous pad:
428     for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
429       
430       //dont merge sequences on the same pad twice
431       if(previousPt[p]->fLastMergedPad==pad) continue;
432
433       Int_t difference = seqmean - previousPt[p]->fMean;
434       if(difference < -fMatch) break;
435
436       if(difference <= fMatch){ //There is a match here!!
437         AliClusterData *local = previousPt[p];
438         
439         if(fDeconvPad){
440           if(seqcharge > local->fLastCharge){
441             if(local->fChargeFalling){ //The previous pad was falling
442               break; //create a new cluster
443             }               
444           }
445           else local->fChargeFalling = 1;
446           local->fLastCharge = seqcharge;
447         }
448               
449         //Don't create a new cluster, because we found a match
450         newcluster = kFALSE;
451               
452         //Update cluster on current pad with the matching one:
453         local->fTotalCharge += seqcharge;
454         local->fPad += padmean;
455         local->fPad2 += paderror;
456         local->fTime += seqaverage;
457         local->fTime2 += seqerror;
458         local->fMean = seqmean;
459         local->fFlags++; //means we have more than one pad 
460         local->fLastMergedPad = pad;
461
462         currentPt[ncurrent] = local;
463         ncurrent++;
464               
465         break;
466       } //Checking for match at previous pad
467     } //Loop over results on previous pad.
468
469
470     if(newcluster && ncurrent<kPadArraySize){
471       //Start a new cluster. Add it to the clusterlist, and update
472       //the list of pointers to clusters in current pad.
473       //current pad will be previous pad on next pad.
474
475       //Add to the clusterlist:
476       AliClusterData *tmp = &clusterlist[ntotal];
477       tmp->fTotalCharge = seqcharge;
478       tmp->fPad = padmean;
479       tmp->fPad2 = paderror;
480       tmp->fTime = seqaverage;
481       tmp->fTime2 = seqerror;
482       tmp->fMean = seqmean;
483       tmp->fFlags = 0;  //flags for single pad clusters
484       tmp->fLastMergedPad = pad;
485
486       if(fDeconvPad){
487         tmp->fChargeFalling = 0;
488         tmp->fLastCharge = seqcharge;
489       }
490
491       //Update list of pointers to previous pad:
492       currentPt[ncurrent] = &clusterlist[ntotal];
493       ntotal++;
494       ncurrent++;
495     }
496
497     if(fDeconvTime)
498       if(newbin >= 0) goto redo;
499   
500     // to prevent endless loop  
501     if(time >= AliHLTTPCTransform::GetNTimeBins()){
502       HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
503       break;
504     }
505
506
507     if(!readValue) break; //No more value
508     
509     if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
510       HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
511       break;
512     }
513
514     if(fCurrentRow != newRow){
515       WriteClusters(ntotal,clusterlist);
516
517       lastpad = 123456789;
518
519       currentPt = pad2;
520       previousPt = pad1;
521       nprevious=0;
522       ncurrent=0;
523       ntotal=0;
524       
525       fCurrentRow = newRow;
526     }
527
528     pad = newPad;
529     time = newTime;
530
531   } // END while(readValue)
532
533   WriteClusters(ntotal,clusterlist);
534
535   HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
536
537 } // ENDEND
538
539 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list)
540 {
541   //write cluster to output pointer
542   Int_t thisrow,thissector;
543   UInt_t counter = fNClusters;
544   
545   for(int j=0; j<nclusters; j++)
546     {
547       if(!list[j].fFlags) continue; //discard single pad clusters
548       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
549
550       Float_t xyz[3];      
551       Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
552       Float_t fpad2=fXYErr*fXYErr; //fixed given error
553       Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
554       Float_t ftime2=fZErr*fZErr;  //fixed given error
555
556
557    
558      
559
560
561       if(fCalcerr) { //calc the errors, otherwice take the fixed error 
562         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
563         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
564         Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
565         sy2/=q2;
566         if(sy2 < 0) {
567             LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
568               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
569             continue;
570         } else {
571           if(!fRawSP){
572             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
573             if(sy2 != 0){
574               fpad2*=0.108; //constants are from offline studies
575               if(patch<2)
576                 fpad2*=2.07;
577             }
578           } else fpad2=sy2; //take the width not the error
579         }
580         Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
581         sz2/=q2;
582         if(sz2 < 0){
583           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
584             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
585           continue;
586         } else {
587           if(!fRawSP){
588             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
589             if(sz2 != 0) {
590               ftime2 *= 0.169; //constants are from offline studies
591               if(patch<2)
592                 ftime2 *= 1.77;
593             }
594           } else ftime2=sz2; //take the width, not the error
595         }
596       }
597       if(fStdout==kTRUE)
598         cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
599       
600       if(!fRawSP){
601         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
602         AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
603         
604         if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
605           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
606         if(fNClusters >= fMaxNClusters)
607           {
608             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
609               <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
610             return;
611           }  
612         
613         fSpacePointData[counter].fX = xyz[0];
614         fSpacePointData[counter].fY = xyz[1];
615         fSpacePointData[counter].fZ = xyz[2];
616         
617       } else {
618         fSpacePointData[counter].fX = fCurrentRow;
619         fSpacePointData[counter].fY = fpad;
620         fSpacePointData[counter].fZ = ftime;
621       }
622       
623       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
624       fSpacePointData[counter].fPadRow = fCurrentRow;
625       fSpacePointData[counter].fSigmaY2 = fpad2;
626       fSpacePointData[counter].fSigmaZ2  = ftime2;
627
628       fSpacePointData[counter].fUsed = kFALSE;         // only used / set in AliHLTTPCDisplay
629       fSpacePointData[counter].fTrackN = -1;           // only used / set in AliHLTTPCDisplay
630
631       Int_t patch=fCurrentPatch;
632       if(patch==-1) patch=0; //never store negative patch number
633       fSpacePointData[counter].fID = counter
634         +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
635
636 #ifdef do_mc
637       Int_t trackID[3];
638       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
639
640       fSpacePointData[counter].fTrackID[0] = trackID[0];
641       fSpacePointData[counter].fTrackID[1] = trackID[1];
642       fSpacePointData[counter].fTrackID[2] = trackID[2];
643
644       //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
645 #endif
646       
647       fNClusters++;
648       counter++;
649     }
650 }
651
652 // STILL TO FIX  ----------------------------------------------------------------------------
653
654 #ifdef do_mc
655 void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
656 {
657   //get mc id
658   AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
659   
660   trackID[0]=trackID[1]=trackID[2]=-2;
661   //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
662   for(Int_t i=fFirstRow; i<=fLastRow; i++){
663     if(rowPt->fRow < (UInt_t)fCurrentRow){
664       AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
665       continue;
666     }
667     AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
668     for(UInt_t j=0; j<rowPt->fNDigit; j++){
669       Int_t cpad = digPt[j].fPad;
670       Int_t ctime = digPt[j].fTime;
671       if(cpad != pad) continue;
672       if(ctime != time) continue;
673
674       trackID[0] = digPt[j].fTrackID[0];
675       trackID[1] = digPt[j].fTrackID[1];
676       trackID[2] = digPt[j].fTrackID[2];
677       
678       //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
679       break;
680     }
681     break;
682   }
683 }
684 #endif