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 <mailto:vestbo@fi.uib.no>, *
9 * Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de> *
10 * Jochen Thaeder <mailto:thaeder@kip.uni-heidelberg.de> *
11 * for The ALICE Off-line Project. *
13 * Permission to use, copy, modify and distribute this software and its *
14 * documentation strictly for non-commercial purposes is hereby granted *
15 * without fee, provided that the above copyright notice appears in all *
16 * copies and that both the copyright notice and this permission notice *
17 * appear in the supporting documentation. The authors make no claims *
18 * about the suitability of this software for any purpose. It is *
19 * provided "as is" without express or implied warranty. *
20 **************************************************************************/
22 /** @file AliHLTTPCClusterFinder.cxx
23 @author Anders Vestbo <mailto:vestbo@fi.uib.no>,
24 Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
25 Jochen Thaeder <mailto:thaeder@kip.uni-heidelberg.de>
27 @brief Cluster Finder for the TPC
30 #include "AliHLTTPCDigitReader.h"
31 #include "AliHLTTPCRootTypes.h"
32 #include "AliHLTTPCLogging.h"
33 #include "AliHLTTPCClusterFinder.h"
34 #include "AliHLTTPCDigitData.h"
35 #include "AliHLTTPCTransform.h"
36 #include "AliHLTTPCSpacePointData.h"
37 #include "AliHLTTPCMemHandler.h"
38 #include "AliHLTTPCPad.h"
44 /** \class AliHLTTPCClusterFinder
46 // The current cluster finder for HLT
49 // The cluster finder is initialized with the Init function,
50 // providing the slice and patch information to work on.
52 // The input is a provided by the AliHLTTPCDigitReader class,
53 // using the init() funktion, and the next() funktion in order
54 // to get the next bin. Either packed or unpacked data can be
55 // processed, dependent if one uses AliHLTTPCDigitReaderPacked
56 // class or AliHLTTPCDigitReaderUnpacked class in the
57 // Clusterfinder Component.
58 // The resulting space points will be in the
59 // array given by the SetOutputArray function.
61 // There are several setters which control the behaviour:
63 // - SetXYError(Float_t): set fixed error in XY direction
64 // - SetZError(Float_t): set fixed error in Z direction
65 // (used if errors are not calculated)
66 // - SetDeconv(Bool_t): switch on/off deconvolution
67 // - SetThreshold(UInt_t): set charge threshold for cluster
68 // - SetMatchWidth(UInt_t): set the match distance in
69 // time for sequences to be merged
70 // - SetSTDOutput(Bool_t): switch on/off output about found clusters
71 // - SetCalcErr(Bool_t): switch on/off calculation of
72 // space point errors (or widths in raw system)
73 // - SetRawSP(Bool_t): switch on/off convertion to raw system
78 // AliHLTTPCFileHandler *file = new AliHLTTPCFileHandler();
79 // file->SetAliInput(digitfile); //give some input file
80 // for(int slice=0; slice<=35; slice++){
81 // for(int patch=0; pat<6; pat++){
82 // file->Init(slice,patch);
84 // UInt_t maxclusters=100000;
85 // UInt_t pointsize = maxclusters*sizeof(AliHLTTPCSpacePointData);
86 // AliHLTTPCSpacePointData *points = (AliHLTTPCSpacePointData*)memory->Allocate(pointsize);
87 // AliHLTTPCDigitRowData *digits = (AliHLTTPCDigitRowData*)file->AliAltroDigits2Memory(ndigits,event);
88 // AliHLTTPCClusterFinder *cf = new AliHLTTPCClusterFinder();
89 // cf->SetMatchWidth(2);
90 // cf->InitSlice( slice, patch, row[0], row[1], maxPoints );
91 // cf->SetSTDOutput(kTRUE); //Some output to standard IO
92 // cf->SetRawSP(kFALSE); //Convert space points to local system
93 // cf->SetThreshold(5); //Threshold of cluster charge
94 // cf->SetDeconv(kTRUE); //Deconv in pad and time direction
95 // cf->SetCalcErr(kTRUE); //Calculate the errors of the spacepoints
96 // cf->SetOutputArray(points); //Move the spacepoints to the array
97 // cf->Read(iter->fPtr, iter->fSize ); //give the data to the cf
98 // cf->ProcessDigits(); //process the rows given by init
99 // Int_t npoints = cf->GetNumberOfClusters();
100 // AliHLTTPCMemHandler *out= new AliHLTTPCMemHandler();
101 // out->SetBinaryOutput(fname);
102 // out->Memory2Binary(npoints,points); //store the spacepoints
103 // out->CloseBinaryOutput();
111 ClassImp(AliHLTTPCClusterFinder)
113 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
115 fSpacePointData(NULL),
131 fSignalThreshold(-1),
141 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder(const AliHLTTPCClusterFinder& src)
143 fSpacePointData(NULL),
144 fDigitReader(src.fDigitReader),
147 fDeconvTime(src.fDeconvTime),
148 fDeconvPad(src.fDeconvPad),
149 fStdout(src.fStdout),
150 fCalcerr(src.fCalcerr),
152 fFirstRow(src.fFirstRow),
153 fLastRow(src.fLastRow),
154 fCurrentRow(src.fCurrentRow),
155 fCurrentSlice(src.fCurrentSlice),
156 fCurrentPatch(src.fCurrentPatch),
158 fThreshold(src.fThreshold),
159 fSignalThreshold(src.fSignalThreshold),
160 fNClusters(src.fNClusters),
161 fMaxNClusters(src.fMaxNClusters),
164 fOccupancyLimit(src.fOccupancyLimit)
168 AliHLTTPCClusterFinder& AliHLTTPCClusterFinder::operator=(const AliHLTTPCClusterFinder& src)
171 fThreshold=src.fThreshold;
172 fSignalThreshold=src.fSignalThreshold;
173 fOccupancyLimit=src.fOccupancyLimit;
176 fDeconvPad=src.fDeconvPad;
177 fDeconvTime=src.fDeconvTime;
179 fCalcerr=src.fCalcerr;
181 fFirstRow=src.fFirstRow;
182 fLastRow=src.fLastRow;
183 fDigitReader=src.fDigitReader;
187 AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder()
192 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
196 fMaxNClusters = nmaxpoints;
197 fCurrentSlice = slice;
198 fCurrentPatch = patch;
199 fFirstRow = firstrow;
203 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
207 fMaxNClusters = nmaxpoints;
208 fCurrentSlice = slice;
209 fCurrentPatch = patch;
210 fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
211 fLastRow=AliHLTTPCTransform::GetLastRow(patch);
214 void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt)
216 //set pointer to output
217 fSpacePointData = pt;
220 void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
222 fPtr = (UChar_t*)ptr;
226 void AliHLTTPCClusterFinder::ProcessDigits()
228 bool readValue = true;
231 UShort_t time=0,newTime=0;
232 UInt_t pad=0,newPad=0;
233 AliHLTTPCSignal_t charge=0;
237 // initialize block for reading packed data
238 fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
239 readValue = fDigitReader->Next();
241 // Matthias 08.11.2006 the following return would cause termination without writing the
242 // ClusterData and thus would block the component. I just want to have the commented line
243 // here for information
244 //if (!readValue)return;
246 pad = fDigitReader->GetPad();
247 time = fDigitReader->GetTime();
248 fCurrentRow = fDigitReader->GetRow();
250 if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
251 rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
253 fCurrentRow += rowOffset;
255 UInt_t lastpad = 123456789;
256 const UInt_t kPadArraySize=5000;
257 const UInt_t kClusterListSize=10000;
258 AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
259 AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
260 AliClusterData clusterlist[kClusterListSize]; //Clusterlist
262 AliClusterData **currentPt; //List of pointers to the current pad
263 AliClusterData **previousPt; //List of pointers to the previous pad
266 UInt_t nprevious=0,ncurrent=0,ntotal=0;
268 /* quick implementation of baseline calculation and zero suppression
269 open a pad object for each pad and delete it after processing.
270 later a list of pad objects with base line history can be used
271 The whole thing only works if we really get unprocessed raw data, if
272 the data is already zero suppressed, there might be gaps in the time
275 Int_t gatingGridOffset=50;
276 AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
277 // just to make later conversion to a list of objects easier
278 AliHLTTPCPad* pCurrentPad=NULL;
279 if (fSignalThreshold>=0) {
280 pCurrentPad=&baseline;
281 baseline.SetThreshold(fSignalThreshold);
284 while ( readValue ){ // Reads through all digits in block
290 if(currentPt == pad2){
298 nprevious = ncurrent;
300 if(pad != lastpad+1){
301 //this happens if there is a pad with no signal.
302 nprevious = ncurrent = 0;
307 Bool_t newcluster = kTRUE;
308 UInt_t seqcharge=0,seqaverage=0,seqerror=0;
309 AliHLTTPCSignal_t lastcharge=0;
310 UInt_t bLastWasFalling=0;
315 redo: //This is a goto.
326 while(1){ //Loop over time bins of current pad
327 // read all the values for one pad at once to calculate the base line
329 if (!pCurrentPad->IsStarted()) {
330 //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
331 pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
332 if ((pCurrentPad->StartEvent())>=0) {
334 if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
335 if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
336 pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
337 //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
338 } while ((readValue = fDigitReader->Next())!=0);
340 pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
341 if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
342 //HLTDebug("no data available after zero suppression");
343 pCurrentPad->StopEvent();
344 pCurrentPad->ResetHistory();
347 time=pCurrentPad->GetCurrentPosition();
348 if (time>pCurrentPad->GetSize()) {
349 HLTError("invalid time bin for pad");
356 Float_t occupancy=pCurrentPad->GetOccupancy();
357 //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
358 if ( occupancy < fOccupancyLimit ) {
359 charge = pCurrentPad->GetCorrectedData();
362 //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
365 charge = fDigitReader->GetSignal();
367 //HLTDebug("get next charge value: position %d charge %d", time, charge);
371 if (fDigitReader->GetRow() == 90){
372 ///// LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90") << "PAD=" << fDigitReader->GetPad() << " TIME=" << fDigitReader->GetTime()
373 // << " SIGNAL=" << fDigitReader->GetSignal() << ENDLOG;
377 if(time >= AliHLTTPCTransform::GetNTimeBins()){
378 HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
383 //Get the current ADC-value
386 //Check if the last pixel in the sequence is smaller than this
387 if(charge > lastcharge){
393 else bLastWasFalling = 1; //last pixel was larger than this
397 //Sum the total charge of this sequence
399 seqaverage += time*charge;
400 seqerror += time*time*charge;
404 if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
405 pCurrentPad->StopEvent();
406 pCurrentPad->ResetHistory();
408 newPad = fDigitReader->GetPad();
409 newTime = fDigitReader->GetTime();
410 newRow = fDigitReader->GetRow() + rowOffset;
415 newPad=pCurrentPad->GetPadNumber();
416 newTime=pCurrentPad->GetCurrentPosition();
417 newRow=pCurrentPad->GetRowNumber();
419 readValue = fDigitReader->Next();
420 //Check where to stop:
421 if(!readValue) break; //No more value
423 newPad = fDigitReader->GetPad();
424 newTime = fDigitReader->GetTime();
425 newRow = fDigitReader->GetRow() + rowOffset;
428 if(newPad != pad)break; //new pad
429 if(newTime != time+1) break; //end of sequence
431 // pad = newpad; is equal
434 }//end loop over sequence
436 //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
437 //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
439 // with active zero suppression zero values are possible
443 //Calculate mean of sequence:
446 seqmean = seqaverage/seqcharge;
448 LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
449 <<"Error in data given to the cluster finder"<<ENDLOG;
454 //Calculate mean in pad direction:
455 Int_t padmean = seqcharge*pad;
456 Int_t paderror = pad*padmean;
458 //Compare with results on previous pad:
459 for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
461 //dont merge sequences on the same pad twice
462 if(previousPt[p]->fLastMergedPad==pad) continue;
464 Int_t difference = seqmean - previousPt[p]->fMean;
465 if(difference < -fMatch) break;
467 if(difference <= fMatch){ //There is a match here!!
468 AliClusterData *local = previousPt[p];
471 if(seqcharge > local->fLastCharge){
472 if(local->fChargeFalling){ //The previous pad was falling
473 break; //create a new cluster
476 else local->fChargeFalling = 1;
477 local->fLastCharge = seqcharge;
480 //Don't create a new cluster, because we found a match
483 //Update cluster on current pad with the matching one:
484 local->fTotalCharge += seqcharge;
485 local->fPad += padmean;
486 local->fPad2 += paderror;
487 local->fTime += seqaverage;
488 local->fTime2 += seqerror;
489 local->fMean = seqmean;
490 local->fFlags++; //means we have more than one pad
491 local->fLastMergedPad = pad;
493 currentPt[ncurrent] = local;
497 } //Checking for match at previous pad
498 } //Loop over results on previous pad.
500 if(newcluster && ncurrent<kPadArraySize){
501 //Start a new cluster. Add it to the clusterlist, and update
502 //the list of pointers to clusters in current pad.
503 //current pad will be previous pad on next pad.
505 //Add to the clusterlist:
506 AliClusterData *tmp = &clusterlist[ntotal];
507 tmp->fTotalCharge = seqcharge;
509 tmp->fPad2 = paderror;
510 tmp->fTime = seqaverage;
511 tmp->fTime2 = seqerror;
512 tmp->fMean = seqmean;
513 tmp->fFlags = 0; //flags for single pad clusters
514 tmp->fLastMergedPad = pad;
517 tmp->fChargeFalling = 0;
518 tmp->fLastCharge = seqcharge;
521 //Update list of pointers to previous pad:
522 currentPt[ncurrent] = &clusterlist[ntotal];
528 if(newbin >= 0) goto redo;
530 // to prevent endless loop
531 if(time >= AliHLTTPCTransform::GetNTimeBins()){
532 HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
537 if(!readValue) break; //No more value
539 if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
540 HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
544 if(fCurrentRow != newRow){
545 WriteClusters(ntotal,clusterlist);
555 fCurrentRow = newRow;
561 } // END while(readValue)
563 WriteClusters(ntotal,clusterlist);
565 HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
569 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list)
571 //write cluster to output pointer
572 Int_t thisrow,thissector;
573 UInt_t counter = fNClusters;
575 for(int j=0; j<nclusters; j++)
577 if(!list[j].fFlags) continue; //discard single pad clusters
578 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
581 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
582 Float_t fpad2=fXYErr*fXYErr; //fixed given error
583 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
584 Float_t ftime2=fZErr*fZErr; //fixed given error
591 if(fCalcerr) { //calc the errors, otherwice take the fixed error
592 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
593 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
594 Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
597 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
598 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
602 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
604 fpad2*=0.108; //constants are from offline studies
608 } else fpad2=sy2; //take the width not the error
610 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
613 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
614 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
618 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
620 ftime2 *= 0.169; //constants are from offline studies
624 } else ftime2=sz2; //take the width, not the error
628 cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
631 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
632 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
634 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
635 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
636 if(fNClusters >= fMaxNClusters)
638 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
639 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
643 fSpacePointData[counter].fX = xyz[0];
644 fSpacePointData[counter].fY = xyz[1];
645 fSpacePointData[counter].fZ = xyz[2];
648 fSpacePointData[counter].fX = fCurrentRow;
649 fSpacePointData[counter].fY = fpad;
650 fSpacePointData[counter].fZ = ftime;
653 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
654 fSpacePointData[counter].fPadRow = fCurrentRow;
655 fSpacePointData[counter].fSigmaY2 = fpad2;
656 fSpacePointData[counter].fSigmaZ2 = ftime2;
658 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
659 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
661 Int_t patch=fCurrentPatch;
662 if(patch==-1) patch=0; //never store negative patch number
663 fSpacePointData[counter].fID = counter
664 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
668 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
670 fSpacePointData[counter].fTrackID[0] = trackID[0];
671 fSpacePointData[counter].fTrackID[1] = trackID[1];
672 fSpacePointData[counter].fTrackID[2] = trackID[2];
674 //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
682 // STILL TO FIX ----------------------------------------------------------------------------
685 void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
688 AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
690 trackID[0]=trackID[1]=trackID[2]=-2;
691 //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
692 for(Int_t i=fFirstRow; i<=fLastRow; i++){
693 if(rowPt->fRow < (UInt_t)fCurrentRow){
694 AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
697 AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
698 for(UInt_t j=0; j<rowPt->fNDigit; j++){
699 Int_t cpad = digPt[j].fPad;
700 Int_t ctime = digPt[j].fTime;
701 if(cpad != pad) continue;
702 if(ctime != time) continue;
704 trackID[0] = digPt[j].fTrackID[0];
705 trackID[1] = digPt[j].fTrackID[1];
706 trackID[2] = digPt[j].fTrackID[2];
708 //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;