bugfix(Kenneth): Reset list of cluster candidates after event processing
[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   bool readValue = true;
238   Int_t newRow = 0;    
239   Int_t rowOffset = 0;
240   UShort_t time=0,newTime=0;
241   UInt_t pad=0,newPad=0;
242   AliHLTTPCSignal_t charge=0;
243
244   fNClusters = 0;
245
246   // initialize block for reading packed data
247   fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
248   readValue = fDigitReader->Next();
249
250   // Matthias 08.11.2006 the following return would cause termination without writing the
251   // ClusterData and thus would block the component. I just want to have the commented line
252   // here for information
253   //if (!readValue)return;
254
255   pad = fDigitReader->GetPad();
256   time = fDigitReader->GetTime();
257   fCurrentRow = fDigitReader->GetRow();
258
259   if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
260     rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
261
262   fCurrentRow += rowOffset;
263
264   UInt_t lastpad = 123456789;
265   const UInt_t kPadArraySize=5000;
266   const UInt_t kClusterListSize=10000;
267   AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
268   AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
269   AliClusterData clusterlist[kClusterListSize]; //Clusterlist
270
271   AliClusterData **currentPt;  //List of pointers to the current pad
272   AliClusterData **previousPt; //List of pointers to the previous pad
273   currentPt = pad2;
274   previousPt = pad1;
275   UInt_t nprevious=0,ncurrent=0,ntotal=0;
276
277   /* quick implementation of baseline calculation and zero suppression
278      open a pad object for each pad and delete it after processing.
279      later a list of pad objects with base line history can be used
280      The whole thing only works if we really get unprocessed raw data, if
281      the data is already zero suppressed, there might be gaps in the time
282      bins.
283    */
284   Int_t gatingGridOffset=50;
285   AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
286   // just to make later conversion to a list of objects easier
287   AliHLTTPCPad* pCurrentPad=NULL;
288   if (fSignalThreshold>=0) {
289     pCurrentPad=&baseline;
290     baseline.SetThreshold(fSignalThreshold);
291   }
292
293   while ( readValue ){   // Reads through all digits in block
294
295     if(pad != lastpad){
296       //This is a new pad
297       
298       //Switch the lists:
299       if(currentPt == pad2){
300         currentPt = pad1;
301         previousPt = pad2;
302       }
303       else {
304         currentPt = pad2;
305         previousPt = pad1;
306       }
307       nprevious = ncurrent;
308       ncurrent = 0;
309       if(pad != lastpad+1){
310         //this happens if there is a pad with no signal.
311         nprevious = ncurrent = 0;
312       }
313       lastpad = pad;
314     }
315
316     Bool_t newcluster = kTRUE;
317     UInt_t seqcharge=0,seqaverage=0,seqerror=0;
318     AliHLTTPCSignal_t lastcharge=0;
319     UInt_t bLastWasFalling=0;
320     Int_t newbin=-1;
321
322
323     if(fDeconvTime){
324       redo: //This is a goto.
325       
326       if(newbin > -1){
327         //bin = newbin;
328         newbin = -1;
329       }
330           
331       lastcharge=0;
332       bLastWasFalling = 0;
333     }
334
335     while(1){ //Loop over time bins of current pad
336       // read all the values for one pad at once to calculate the base line
337       if (pCurrentPad) {
338         if (!pCurrentPad->IsStarted()) {
339           //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
340           pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
341           if ((pCurrentPad->StartEvent())>=0) {
342             do {
343               if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
344               if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
345               pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
346               //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
347             } while ((readValue = fDigitReader->Next())!=0);
348           }
349           pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
350           if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
351             //HLTDebug("no data available after zero suppression");
352             pCurrentPad->StopEvent();
353             pCurrentPad->ResetHistory();
354             break;
355           }
356           time=pCurrentPad->GetCurrentPosition();
357           if (time>pCurrentPad->GetSize()) {
358             HLTError("invalid time bin for pad");
359             break;
360           }
361         }
362       }
363
364       if (pCurrentPad) {
365         Float_t occupancy=pCurrentPad->GetOccupancy();
366         //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
367         if ( occupancy < fOccupancyLimit ) {
368           charge = pCurrentPad->GetCorrectedData();
369         } else {
370           charge = 0;
371           //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
372         }
373       } else {
374         charge = fDigitReader->GetSignal();
375       }
376       //HLTDebug("get next charge value: position %d charge %d", time, charge);
377
378
379       // CHARGE DEBUG
380       if (fDigitReader->GetRow() == 90){
381 /////     LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90")  << "PAD=" <<  fDigitReader->GetPad() << "  TIME=" <<  fDigitReader->GetTime() 
382           //                                       << "  SIGNAL=" <<  fDigitReader->GetSignal() << ENDLOG;
383
384       }
385
386       if(time >= AliHLTTPCTransform::GetNTimeBins()){
387         HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
388         break;
389       }
390
391
392       //Get the current ADC-value
393       if(fDeconvTime){
394
395         //Check if the last pixel in the sequence is smaller than this
396         if(charge > lastcharge){
397           if(bLastWasFalling){
398             newbin = 1;
399             break;
400           }
401         }
402         else bLastWasFalling = 1; //last pixel was larger than this
403         lastcharge = charge;
404       }
405           
406       //Sum the total charge of this sequence
407       seqcharge += charge;
408       seqaverage += time*charge;
409       seqerror += time*time*charge;
410       
411       if (pCurrentPad) {
412         
413         if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
414           pCurrentPad->StopEvent();
415           pCurrentPad->ResetHistory();
416           if(readValue) {
417             newPad = fDigitReader->GetPad();
418             newTime = fDigitReader->GetTime();
419             newRow = fDigitReader->GetRow() + rowOffset;
420           }
421           break;
422         }
423
424         newPad=pCurrentPad->GetPadNumber();
425         newTime=pCurrentPad->GetCurrentPosition();
426         newRow=pCurrentPad->GetRowNumber();
427       } else {
428       readValue = fDigitReader->Next();
429       //Check where to stop:
430       if(!readValue) break; //No more value
431
432       newPad = fDigitReader->GetPad();
433       newTime = fDigitReader->GetTime();
434       newRow = fDigitReader->GetRow() + rowOffset;
435       }
436
437       if(newPad != pad)break; //new pad
438       if(newTime != time+1) break; //end of sequence
439
440       // pad = newpad;    is equal
441       time = newTime;
442
443     }//end loop over sequence
444
445     //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
446     //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
447     if (seqcharge<=0) {
448       // with active zero suppression zero values are possible
449       continue;
450     }
451
452     //Calculate mean of sequence:
453     Int_t seqmean=0;
454     if(seqcharge)
455       seqmean = seqaverage/seqcharge;
456     else{
457       LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
458         <<"Error in data given to the cluster finder"<<ENDLOG;
459       seqmean = 1;
460       seqcharge = 1;
461     }
462
463     //Calculate mean in pad direction:
464     Int_t padmean = seqcharge*pad;
465     Int_t paderror = pad*padmean;
466
467     //Compare with results on previous pad:
468     for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
469       
470       //dont merge sequences on the same pad twice
471       if(previousPt[p]->fLastMergedPad==pad) continue;
472
473       Int_t difference = seqmean - previousPt[p]->fMean;
474       if(difference < -fMatch) break;
475
476       if(difference <= fMatch){ //There is a match here!!
477         AliClusterData *local = previousPt[p];
478         
479         if(fDeconvPad){
480           if(seqcharge > local->fLastCharge){
481             if(local->fChargeFalling){ //The previous pad was falling
482               break; //create a new cluster
483             }               
484           }
485           else local->fChargeFalling = 1;
486           local->fLastCharge = seqcharge;
487         }
488               
489         //Don't create a new cluster, because we found a match
490         newcluster = kFALSE;
491               
492         //Update cluster on current pad with the matching one:
493         local->fTotalCharge += seqcharge;
494         local->fPad += padmean;
495         local->fPad2 += paderror;
496         local->fTime += seqaverage;
497         local->fTime2 += seqerror;
498         local->fMean = seqmean;
499         local->fFlags++; //means we have more than one pad 
500         local->fLastMergedPad = pad;
501
502         currentPt[ncurrent] = local;
503         ncurrent++;
504               
505         break;
506       } //Checking for match at previous pad
507     } //Loop over results on previous pad.
508
509     if(newcluster && ncurrent<kPadArraySize){
510       //Start a new cluster. Add it to the clusterlist, and update
511       //the list of pointers to clusters in current pad.
512       //current pad will be previous pad on next pad.
513
514       //Add to the clusterlist:
515       AliClusterData *tmp = &clusterlist[ntotal];
516       tmp->fTotalCharge = seqcharge;
517       tmp->fPad = padmean;
518       tmp->fPad2 = paderror;
519       tmp->fTime = seqaverage;
520       tmp->fTime2 = seqerror;
521       tmp->fMean = seqmean;
522       tmp->fFlags = 0;  //flags for single pad clusters
523       tmp->fLastMergedPad = pad;
524
525       if(fDeconvPad){
526         tmp->fChargeFalling = 0;
527         tmp->fLastCharge = seqcharge;
528       }
529
530       //Update list of pointers to previous pad:
531       currentPt[ncurrent] = &clusterlist[ntotal];
532       ntotal++;
533       ncurrent++;
534     }
535
536     if(fDeconvTime)
537       if(newbin >= 0) goto redo;
538   
539     // to prevent endless loop  
540     if(time >= AliHLTTPCTransform::GetNTimeBins()){
541       HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
542       break;
543     }
544
545
546     if(!readValue) break; //No more value
547     
548     if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
549       HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
550       break;
551     }
552
553     if(fCurrentRow != newRow){
554       WriteClusters(ntotal,clusterlist);
555
556       lastpad = 123456789;
557
558       currentPt = pad2;
559       previousPt = pad1;
560       nprevious=0;
561       ncurrent=0;
562       ntotal=0;
563       
564       fCurrentRow = newRow;
565     }
566
567     pad = newPad;
568     time = newTime;
569
570   } // END while(readValue)
571
572   WriteClusters(ntotal,clusterlist);
573
574   HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
575
576 } // ENDEND
577
578 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list)
579 {
580   //write cluster to output pointer
581   Int_t thisrow,thissector;
582   UInt_t counter = fNClusters;
583   
584   for(int j=0; j<nclusters; j++)
585     {
586
587
588
589       if(!list[j].fFlags) continue; //discard single pad clusters
590       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
591
592       Float_t xyz[3];      
593       Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
594       Float_t fpad2=fXYErr*fXYErr; //fixed given error
595       Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
596       Float_t ftime2=fZErr*fZErr;  //fixed given error
597
598
599 #if UNSORTED
600       fCurrentRow=list[j].fRow;
601 #endif
602    
603       if(fCalcerr) { //calc the errors, otherwice take the fixed error 
604         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
605         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
606         Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
607         sy2/=q2;
608         if(sy2 < 0) {
609             LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
610               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
611             continue;
612         } else {
613           if(!fRawSP){
614             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
615             if(sy2 != 0){
616               fpad2*=0.108; //constants are from offline studies
617               if(patch<2)
618                 fpad2*=2.07;
619             }
620           } else fpad2=sy2; //take the width not the error
621         }
622         Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
623         sz2/=q2;
624         if(sz2 < 0){
625           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
626             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
627           continue;
628         } else {
629           if(!fRawSP){
630             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
631             if(sz2 != 0) {
632               ftime2 *= 0.169; //constants are from offline studies
633               if(patch<2)
634                 ftime2 *= 1.77;
635             }
636           } else ftime2=sz2; //take the width, not the error
637         }
638       }
639       if(fStdout==kTRUE)
640         cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
641       
642       if(!fRawSP){
643         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
644         AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
645         
646         if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
647           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
648         if(fNClusters >= fMaxNClusters)
649           {
650             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
651               <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
652             return;
653           }  
654         
655         fSpacePointData[counter].fX = xyz[0];
656         fSpacePointData[counter].fY = xyz[1];
657         fSpacePointData[counter].fZ = xyz[2];
658         
659       } else {
660         fSpacePointData[counter].fX = fCurrentRow;
661         fSpacePointData[counter].fY = fpad;
662         fSpacePointData[counter].fZ = ftime;
663       }
664       
665       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
666       fSpacePointData[counter].fPadRow = fCurrentRow;
667       fSpacePointData[counter].fSigmaY2 = fpad2;
668       fSpacePointData[counter].fSigmaZ2  = ftime2;
669
670       fSpacePointData[counter].fUsed = kFALSE;         // only used / set in AliHLTTPCDisplay
671       fSpacePointData[counter].fTrackN = -1;           // only used / set in AliHLTTPCDisplay
672
673       Int_t patch=fCurrentPatch;
674       if(patch==-1) patch=0; //never store negative patch number
675       fSpacePointData[counter].fID = counter
676         +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
677
678 #ifdef do_mc
679       Int_t trackID[3];
680       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
681
682       fSpacePointData[counter].fTrackID[0] = trackID[0];
683       fSpacePointData[counter].fTrackID[1] = trackID[1];
684       fSpacePointData[counter].fTrackID[2] = trackID[2];
685
686       //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
687 #endif
688       
689       fNClusters++;
690       counter++;
691     }
692 }
693
694 // STILL TO FIX  ----------------------------------------------------------------------------
695
696 #ifdef do_mc
697 void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
698 {
699   //get mc id
700   AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
701   
702   trackID[0]=trackID[1]=trackID[2]=-2;
703   //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
704   for(Int_t i=fFirstRow; i<=fLastRow; i++){
705     if(rowPt->fRow < (UInt_t)fCurrentRow){
706       AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
707       continue;
708     }
709     AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
710     for(UInt_t j=0; j<rowPt->fNDigit; j++){
711       Int_t cpad = digPt[j].fPad;
712       Int_t ctime = digPt[j].fTime;
713       if(cpad != pad) continue;
714       if(ctime != time) continue;
715
716       trackID[0] = digPt[j].fTrackID[0];
717       trackID[1] = digPt[j].fTrackID[1];
718       trackID[2] = digPt[j].fTrackID[2];
719       
720       //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
721       break;
722     }
723     break;
724   }
725 }
726 #endif
727
728 //----------------------------------Methods for the new unsorted way of reading the data --------------------------------
729
730 void AliHLTTPCClusterFinder::SetPadArray(AliHLTTPCPadArray * padArray){
731   fPadArray=padArray;
732 }
733
734 void AliHLTTPCClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size){
735   //set input pointer
736   fPtr = (UChar_t*)ptr;
737   fSize = size;
738
739   fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
740
741   fPadArray->SetDigitReader(fDigitReader);
742   fPadArray->ReadData();
743 }
744 void AliHLTTPCClusterFinder::FindClusters(){
745   fPadArray->FindClusterCandidates();
746   fPadArray->FindClusters(fMatch);
747
748   AliClusterData * clusterlist = new AliClusterData[fPadArray->fClusters.size()]; //Clusterlist
749   for(int i=0;i<fPadArray->fClusters.size();i++){
750     clusterlist[i].fTotalCharge = fPadArray->fClusters[i].fTotalCharge;
751     clusterlist[i].fPad = fPadArray->fClusters[i].fPad;
752     clusterlist[i].fPad2 = fPadArray->fClusters[i].fPad2;
753     clusterlist[i].fTime = fPadArray->fClusters[i].fTime;
754     clusterlist[i].fTime2 = fPadArray->fClusters[i].fTime2;
755     clusterlist[i].fMean = fPadArray->fClusters[i].fMean;
756     clusterlist[i].fFlags = fPadArray->fClusters[i].fFlags;
757     clusterlist[i].fChargeFalling = fPadArray->fClusters[i].fChargeFalling;
758     clusterlist[i].fLastCharge = fPadArray->fClusters[i].fLastCharge;
759     clusterlist[i].fLastMergedPad = fPadArray->fClusters[i].fLastMergedPad;
760     clusterlist[i].fRow = fPadArray->fClusters[i].fRowNumber;
761   }
762   WriteClusters(fPadArray->fClusters.size(),clusterlist);
763   delete [] clusterlist;
764   fPadArray->DataToDefault();
765 }
766 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list)//This is used when using the AliHLTTPCClusters class for cluster data
767 {
768   //write cluster to output pointer
769   Int_t thisrow,thissector;
770   UInt_t counter = fNClusters;
771   
772   for(int j=0; j<nclusters; j++)
773     {
774       if(!list[j].fFlags) continue; //discard single pad clusters
775       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
776
777       Float_t xyz[3];      
778       Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
779       Float_t fpad2=fXYErr*fXYErr; //fixed given error
780       Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
781       Float_t ftime2=fZErr*fZErr;  //fixed given error
782
783
784       if(fCalcerr) { //calc the errors, otherwice take the fixed error 
785         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
786         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
787         Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
788         sy2/=q2;
789         if(sy2 < 0) {
790             LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
791               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
792             continue;
793         } else {
794           if(!fRawSP){
795             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
796             if(sy2 != 0){
797               fpad2*=0.108; //constants are from offline studies
798               if(patch<2)
799                 fpad2*=2.07;
800             }
801           } else fpad2=sy2; //take the width not the error
802         }
803         Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
804         sz2/=q2;
805         if(sz2 < 0){
806           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
807             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
808           continue;
809         } else {
810           if(!fRawSP){
811             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
812             if(sz2 != 0) {
813               ftime2 *= 0.169; //constants are from offline studies
814               if(patch<2)
815                 ftime2 *= 1.77;
816             }
817           } else ftime2=sz2; //take the width, not the error
818         }
819       }
820       if(fStdout==kTRUE)
821         cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
822       
823       if(!fRawSP){
824         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
825         AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
826         
827         if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
828           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
829         if(fNClusters >= fMaxNClusters)
830           {
831             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
832               <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
833             return;
834           }  
835         
836         fSpacePointData[counter].fX = xyz[0];
837         fSpacePointData[counter].fY = xyz[1];
838         fSpacePointData[counter].fZ = xyz[2];
839         
840       } else {
841         fSpacePointData[counter].fX = fCurrentRow;
842         fSpacePointData[counter].fY = fpad;
843         fSpacePointData[counter].fZ = ftime;
844       }
845       
846       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
847       fSpacePointData[counter].fPadRow = fCurrentRow;
848       fSpacePointData[counter].fSigmaY2 = fpad2;
849       fSpacePointData[counter].fSigmaZ2  = ftime2;
850
851       fSpacePointData[counter].fUsed = kFALSE;         // only used / set in AliHLTTPCDisplay
852       fSpacePointData[counter].fTrackN = -1;           // only used / set in AliHLTTPCDisplay
853
854       Int_t patch=fCurrentPatch;
855       if(patch==-1) patch=0; //never store negative patch number
856       fSpacePointData[counter].fID = counter
857         +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
858
859 #ifdef do_mc
860       Int_t trackID[3];
861       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
862
863       fSpacePointData[counter].fTrackID[0] = trackID[0];
864       fSpacePointData[counter].fTrackID[1] = trackID[1];
865       fSpacePointData[counter].fTrackID[2] = trackID[2];
866
867       //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
868 #endif
869       
870       fNClusters++;
871       counter++;
872     }
873 }