2 // Original: AliHLTClustFinderNew.cxx,v 1.29 2005/06/14 10:55:21 cvetan Exp
4 /**************************************************************************
5 * This file is property of and copyright by the ALICE HLT Project *
6 * ALICE Experiment at CERN, All rights reserved. *
8 * Primary Authors: Anders Vestbo, Constantin Loizides, Jochen Thaeder *
9 * Kenneth Aamodt <kenneth.aamodt@student.uib.no> *
10 * for The ALICE HLT Project. *
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 **************************************************************************/
21 /** @file AliHLTTPCClusterFinder.cxx
22 @author Anders Vestbo, Constantin Loizides, Jochen Thaeder
23 Kenneth Aamodt kenneth.aamodt@student.uib.no
25 @brief Cluster Finder for the TPC
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"
43 /** \class AliHLTTPCClusterFinder
45 // The current cluster finder for HLT
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.
56 // The cluster finder is initialized with the Init function,
57 // providing the slice and patch information to work on.
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.
68 // There are several setters which control the behaviour:
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
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);
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();
118 ClassImp(AliHLTTPCClusterFinder)
120 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
122 fSpacePointData(NULL),
138 fSignalThreshold(-1),
143 fOccupancyLimit(1.0),
150 AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder()
155 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
159 fMaxNClusters = nmaxpoints;
160 fCurrentSlice = slice;
161 fCurrentPatch = patch;
162 fFirstRow = firstrow;
166 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
170 fMaxNClusters = nmaxpoints;
171 fCurrentSlice = slice;
172 fCurrentPatch = patch;
173 fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
174 fLastRow=AliHLTTPCTransform::GetLastRow(patch);
177 void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt)
179 //set pointer to output
180 fSpacePointData = pt;
183 void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
185 fPtr = (UChar_t*)ptr;
189 void AliHLTTPCClusterFinder::ProcessDigits()
192 bool readValue = true;
195 UShort_t time=0,newTime=0;
196 UInt_t pad=0,newPad=0;
197 AliHLTTPCSignal_t charge=0;
201 // initialize block for reading packed data
202 fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
203 readValue = fDigitReader->Next();
205 // Matthias 08.11.2006 the following return would cause termination without writing the
206 // ClusterData and thus would block the component. I just want to have the commented line
207 // here for information
208 //if (!readValue)return;
210 pad = fDigitReader->GetPad();
211 time = fDigitReader->GetTime();
212 fCurrentRow = fDigitReader->GetRow();
214 if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
215 rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
217 fCurrentRow += rowOffset;
219 UInt_t lastpad = 123456789;
220 const UInt_t kPadArraySize=5000;
221 const UInt_t kClusterListSize=10000;
222 AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
223 AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
224 AliClusterData clusterlist[kClusterListSize]; //Clusterlist
226 AliClusterData **currentPt; //List of pointers to the current pad
227 AliClusterData **previousPt; //List of pointers to the previous pad
230 UInt_t nprevious=0,ncurrent=0,ntotal=0;
232 /* quick implementation of baseline calculation and zero suppression
233 open a pad object for each pad and delete it after processing.
234 later a list of pad objects with base line history can be used
235 The whole thing only works if we really get unprocessed raw data, if
236 the data is already zero suppressed, there might be gaps in the time
239 Int_t gatingGridOffset=50;
240 AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
241 // just to make later conversion to a list of objects easier
242 AliHLTTPCPad* pCurrentPad=NULL;
243 if (fSignalThreshold>=0) {
244 pCurrentPad=&baseline;
245 baseline.SetThreshold(fSignalThreshold);
248 while ( readValue!=0 && iResult>=0){ // Reads through all digits in block
255 if(currentPt == pad2){
263 nprevious = ncurrent;
265 if(pad != lastpad+1){
266 //this happens if there is a pad with no signal.
267 nprevious = ncurrent = 0;
272 Bool_t newcluster = kTRUE;
273 UInt_t seqcharge=0,seqaverage=0,seqerror=0;
274 AliHLTTPCSignal_t lastcharge=0;
275 UInt_t bLastWasFalling=0;
280 redo: //This is a goto.
291 while(iResult>=0){ //Loop over time bins of current pad
292 // read all the values for one pad at once to calculate the base line
294 if (!pCurrentPad->IsStarted()) {
295 //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
296 pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
297 if ((pCurrentPad->StartEvent())>=0) {
299 if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
300 if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
301 pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
302 //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
303 } while ((readValue = fDigitReader->Next())!=0);
305 pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
306 if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
307 //HLTDebug("no data available after zero suppression");
308 pCurrentPad->StopEvent();
309 pCurrentPad->ResetHistory();
312 time=pCurrentPad->GetCurrentPosition();
313 if (time>pCurrentPad->GetSize()) {
314 HLTError("invalid time bin for pad");
321 Float_t occupancy=pCurrentPad->GetOccupancy();
322 //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
323 if ( occupancy < fOccupancyLimit ) {
324 charge = pCurrentPad->GetCorrectedData();
327 //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
330 charge = fDigitReader->GetSignal();
332 //HLTDebug("get next charge value: position %d charge %d", time, charge);
336 if (fDigitReader->GetRow() == 90){
337 ///// LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90") << "PAD=" << fDigitReader->GetPad() << " TIME=" << fDigitReader->GetTime()
338 // << " SIGNAL=" << fDigitReader->GetSignal() << ENDLOG;
342 if(time >= AliHLTTPCTransform::GetNTimeBins()){
343 HLTWarning("Pad %d: Timebin (%d) out of range (%d)", pad, time, AliHLTTPCTransform::GetNTimeBins());
348 //Get the current ADC-value
351 //Check if the last pixel in the sequence is smaller than this
352 if(charge > lastcharge){
358 else bLastWasFalling = 1; //last pixel was larger than this
362 //Sum the total charge of this sequence
364 seqaverage += time*charge;
365 seqerror += time*time*charge;
369 if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
370 pCurrentPad->StopEvent();
371 pCurrentPad->ResetHistory();
373 newPad = fDigitReader->GetPad();
374 newTime = fDigitReader->GetTime();
375 newRow = fDigitReader->GetRow() + rowOffset;
380 newPad=pCurrentPad->GetPadNumber();
381 newTime=pCurrentPad->GetCurrentPosition();
382 newRow=pCurrentPad->GetRowNumber();
384 readValue = fDigitReader->Next();
385 //Check where to stop:
386 if(!readValue) break; //No more value
388 newPad = fDigitReader->GetPad();
389 newTime = fDigitReader->GetTime();
390 newRow = fDigitReader->GetRow() + rowOffset;
393 if(newPad != pad)break; //new pad
394 if(newTime != time+1) break; //end of sequence
397 // pad = newpad; is equal
400 }//end loop over sequence
402 //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
403 //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
405 // with active zero suppression zero values are possible
409 //Calculate mean of sequence:
412 seqmean = seqaverage/seqcharge;
414 LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
415 <<"Error in data given to the cluster finder"<<ENDLOG;
420 //Calculate mean in pad direction:
421 Int_t padmean = seqcharge*pad;
422 Int_t paderror = pad*padmean;
424 //Compare with results on previous pad:
425 for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
427 //dont merge sequences on the same pad twice
428 if(previousPt[p]->fLastMergedPad==pad) continue;
430 Int_t difference = seqmean - previousPt[p]->fMean;
431 if(difference < -fMatch) break;
433 if(difference <= fMatch){ //There is a match here!!
434 AliClusterData *local = previousPt[p];
437 if(seqcharge > local->fLastCharge){
438 if(local->fChargeFalling){ //The previous pad was falling
439 break; //create a new cluster
442 else local->fChargeFalling = 1;
443 local->fLastCharge = seqcharge;
446 //Don't create a new cluster, because we found a match
449 //Update cluster on current pad with the matching one:
450 local->fTotalCharge += seqcharge;
451 local->fPad += padmean;
452 local->fPad2 += paderror;
453 local->fTime += seqaverage;
454 local->fTime2 += seqerror;
455 local->fMean = seqmean;
456 local->fFlags++; //means we have more than one pad
457 local->fLastMergedPad = pad;
459 currentPt[ncurrent] = local;
463 } //Checking for match at previous pad
464 } //Loop over results on previous pad.
466 if(newcluster && ncurrent<kPadArraySize){
467 //Start a new cluster. Add it to the clusterlist, and update
468 //the list of pointers to clusters in current pad.
469 //current pad will be previous pad on next pad.
471 //Add to the clusterlist:
472 AliClusterData *tmp = &clusterlist[ntotal];
473 tmp->fTotalCharge = seqcharge;
475 tmp->fPad2 = paderror;
476 tmp->fTime = seqaverage;
477 tmp->fTime2 = seqerror;
478 tmp->fMean = seqmean;
479 tmp->fFlags = 0; //flags for single pad clusters
480 tmp->fLastMergedPad = pad;
483 tmp->fChargeFalling = 0;
484 tmp->fLastCharge = seqcharge;
487 //Update list of pointers to previous pad:
488 currentPt[ncurrent] = &clusterlist[ntotal];
494 if(newbin >= 0) goto redo;
496 // to prevent endless loop
497 if(time >= AliHLTTPCTransform::GetNTimeBins()){
498 HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
503 if(!readValue) break; //No more value
505 if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
506 HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
510 if(fCurrentRow != newRow){
511 WriteClusters(ntotal,clusterlist);
521 fCurrentRow = newRow;
527 } // END while(readValue)
529 WriteClusters(ntotal,clusterlist);
531 HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
535 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list)
537 //write cluster to output pointer
538 Int_t thisrow=-1,thissector=-1;
539 UInt_t counter = fNClusters;
541 for(int j=0; j<nclusters; j++)
546 if(!list[j].fFlags) continue; //discard single pad clusters
547 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
550 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
551 Float_t fpad2=fXYErr*fXYErr; //fixed given error
552 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
553 Float_t ftime2=fZErr*fZErr; //fixed given error
558 fCurrentRow=list[j].fRow;
562 if(fCalcerr) { //calc the errors, otherwice take the fixed error
563 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
564 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
565 Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
568 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
569 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
573 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
575 fpad2*=0.108; //constants are from offline studies
579 } else fpad2=sy2; //take the width not the error
581 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
584 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
585 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
589 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
591 ftime2 *= 0.169; //constants are from offline studies
595 } else ftime2=sz2; //take the width, not the error
599 cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
602 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
603 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
605 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
606 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
607 if(fNClusters >= fMaxNClusters)
609 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
610 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
614 fSpacePointData[counter].fX = xyz[0];
615 fSpacePointData[counter].fY = xyz[1];
616 fSpacePointData[counter].fZ = xyz[2];
619 fSpacePointData[counter].fX = fCurrentRow;
620 fSpacePointData[counter].fY = fpad;
621 fSpacePointData[counter].fZ = ftime;
624 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
625 fSpacePointData[counter].fPadRow = fCurrentRow;
626 fSpacePointData[counter].fSigmaY2 = fpad2;
627 fSpacePointData[counter].fSigmaZ2 = ftime2;
629 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
630 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
632 Int_t patch=fCurrentPatch;
633 if(patch==-1) patch=0; //never store negative patch number
634 fSpacePointData[counter].fID = counter
635 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
639 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
641 fSpacePointData[counter].fTrackID[0] = trackID[0];
642 fSpacePointData[counter].fTrackID[1] = trackID[1];
643 fSpacePointData[counter].fTrackID[2] = trackID[2];
645 //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
653 // STILL TO FIX ----------------------------------------------------------------------------
656 void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
659 AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
661 trackID[0]=trackID[1]=trackID[2]=-2;
662 //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
663 for(Int_t i=fFirstRow; i<=fLastRow; i++){
664 if(rowPt->fRow < (UInt_t)fCurrentRow){
665 AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
668 AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
669 for(UInt_t j=0; j<rowPt->fNDigit; j++){
670 Int_t cpad = digPt[j].fPad;
671 Int_t ctime = digPt[j].fTime;
672 if(cpad != pad) continue;
673 if(ctime != time) continue;
675 trackID[0] = digPt[j].fTrackID[0];
676 trackID[1] = digPt[j].fTrackID[1];
677 trackID[2] = digPt[j].fTrackID[2];
679 //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
687 //----------------------------------Methods for the new unsorted way of reading the data --------------------------------
689 void AliHLTTPCClusterFinder::SetPadArray(AliHLTTPCPadArray * padArray){
690 // see header file for function documentation
694 void AliHLTTPCClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size){
696 fPtr = (UChar_t*)ptr;
699 fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
701 fPadArray->SetDigitReader(fDigitReader);
703 if(fSignalThreshold>0){
704 fPadArray->SetSignalThreshold(fSignalThreshold);
706 if(fNSigmaThreshold>0){
707 fPadArray->SetNSigmaThreshold(fNSigmaThreshold);
709 fPadArray->ReadData();
711 void AliHLTTPCClusterFinder::FindClusters(){
712 // see header file for function documentation
713 fPadArray->FindClusterCandidates();
714 fPadArray->FindClusters(fMatch);
716 AliClusterData * clusterlist = new AliClusterData[fPadArray->fClusters.size()]; //Clusterlist
717 for(unsigned int i=0;i<fPadArray->fClusters.size();i++){
718 clusterlist[i].fTotalCharge = fPadArray->fClusters[i].fTotalCharge;
719 clusterlist[i].fPad = fPadArray->fClusters[i].fPad;
720 clusterlist[i].fPad2 = fPadArray->fClusters[i].fPad2;
721 clusterlist[i].fTime = fPadArray->fClusters[i].fTime;
722 clusterlist[i].fTime2 = fPadArray->fClusters[i].fTime2;
723 clusterlist[i].fMean = fPadArray->fClusters[i].fMean;
724 clusterlist[i].fFlags = fPadArray->fClusters[i].fFlags;
725 clusterlist[i].fChargeFalling = fPadArray->fClusters[i].fChargeFalling;
726 clusterlist[i].fLastCharge = fPadArray->fClusters[i].fLastCharge;
727 clusterlist[i].fLastMergedPad = fPadArray->fClusters[i].fLastMergedPad;
728 clusterlist[i].fRow = fPadArray->fClusters[i].fRowNumber;
730 WriteClusters(fPadArray->fClusters.size(),clusterlist);
731 delete [] clusterlist;
732 fPadArray->DataToDefault();
734 Int_t AliHLTTPCClusterFinder::GetActivePads(AliHLTTPCPadArray::AliHLTTPCActivePads* activePads,Int_t maxActivePads){
735 // see header file for function documentation
736 return fPadArray->GetActivePads((AliHLTTPCPadArray::AliHLTTPCActivePads*)activePads , maxActivePads);
738 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list)//This is used when using the AliHLTTPCClusters class for cluster data
740 //write cluster to output pointer
741 Int_t thisrow,thissector;
742 UInt_t counter = fNClusters;
744 for(int j=0; j<nclusters; j++)
746 if(!list[j].fFlags) continue; //discard single pad clusters
747 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
750 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
751 Float_t fpad2=fXYErr*fXYErr; //fixed given error
752 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
753 Float_t ftime2=fZErr*fZErr; //fixed given error
756 if(fCalcerr) { //calc the errors, otherwice take the fixed error
757 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
758 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
759 Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
762 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
763 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
767 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
769 fpad2*=0.108; //constants are from offline studies
773 } else fpad2=sy2; //take the width not the error
775 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
778 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
779 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
783 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
785 ftime2 *= 0.169; //constants are from offline studies
789 } else ftime2=sz2; //take the width, not the error
793 cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
796 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
797 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
799 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
800 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
801 if(fNClusters >= fMaxNClusters)
803 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
804 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
808 fSpacePointData[counter].fX = xyz[0];
809 fSpacePointData[counter].fY = xyz[1];
810 fSpacePointData[counter].fZ = xyz[2];
813 fSpacePointData[counter].fX = fCurrentRow;
814 fSpacePointData[counter].fY = fpad;
815 fSpacePointData[counter].fZ = ftime;
818 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
819 fSpacePointData[counter].fPadRow = fCurrentRow;
820 fSpacePointData[counter].fSigmaY2 = fpad2;
821 fSpacePointData[counter].fSigmaZ2 = ftime2;
823 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
824 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
826 Int_t patch=fCurrentPatch;
827 if(patch==-1) patch=0; //never store negative patch number
828 fSpacePointData[counter].fID = counter
829 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
833 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
835 fSpacePointData[counter].fTrackID[0] = trackID[0];
836 fSpacePointData[counter].fTrackID[1] = trackID[1];
837 fSpacePointData[counter].fTrackID[2] = trackID[2];
839 //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;