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