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