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 *
9 //* Developers: Kenneth Aamodt <kenneth.aamodt@student.uib.no> *
11 //* for The ALICE HLT 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 Kenneth Aamodt, Kalliopi Kanaki
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 "AliHLTTPCSpacePointData.h"
34 #include "AliHLTTPCMemHandler.h"
35 #include "AliHLTTPCPad.h"
42 ClassImp(AliHLTTPCClusterFinder)
44 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
46 fClustersHWAddressVector(),
48 fSpacePointData(NULL),
70 fVectorInitialized(kFALSE),
72 fNumberOfPadsInRow(NULL),
74 fRowOfFirstCandidate(0),
75 fDoPadSelection(kFALSE),
77 fLastTimeBin(AliHLTTPCTransform::GetNTimeBins()),
78 fTotalChargeOfPreviousClusterCandidate(0),
79 fChargeOfCandidatesFalling(kFALSE)
84 AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder()
87 if(fVectorInitialized){
88 DeInitializePadArray();
90 if(fNumberOfPadsInRow){
91 delete [] fNumberOfPadsInRow;
92 fNumberOfPadsInRow=NULL;
96 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
100 fMaxNClusters = nmaxpoints;
101 fCurrentSlice = slice;
102 fCurrentPatch = patch;
103 fFirstRow = firstrow;
107 void AliHLTTPCClusterFinder::InitializePadArray(){
108 // see header file for class documentation
110 if(fCurrentPatch>5||fCurrentPatch<0){
111 HLTFatal("Patch is not set");
115 HLTDebug("Patch number=%d",fCurrentPatch);
117 fFirstRow = AliHLTTPCTransform::GetFirstRow(fCurrentPatch);
118 fLastRow = AliHLTTPCTransform::GetLastRow(fCurrentPatch);
120 fNumberOfRows=fLastRow-fFirstRow+1;
121 fNumberOfPadsInRow= new UInt_t[fNumberOfRows];
123 memset( fNumberOfPadsInRow, 0, sizeof(Int_t)*(fNumberOfRows));
125 for(UInt_t i=0;i<fNumberOfRows;i++){
126 fNumberOfPadsInRow[i]=AliHLTTPCTransform::GetNPads(i+fFirstRow);
127 AliHLTTPCPadVector tmpRow;
128 for(UInt_t j=0;j<=fNumberOfPadsInRow[i];j++){
129 AliHLTTPCPad *tmpPad = new AliHLTTPCPad(2);
131 tmpRow.push_back(tmpPad);
133 fRowPadVector.push_back(tmpRow);
135 fVectorInitialized=kTRUE;
138 Int_t AliHLTTPCClusterFinder::DeInitializePadArray()
140 // see header file for class documentation
141 for(UInt_t i=0;i<fNumberOfRows;i++){
142 for(UInt_t j=0;j<=fNumberOfPadsInRow[i];j++){
143 delete fRowPadVector[i][j];
144 fRowPadVector[i][j]=NULL;
146 fRowPadVector[i].clear();
148 fRowPadVector.clear();
153 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
157 fMaxNClusters = nmaxpoints;
158 fCurrentSlice = slice;
159 fCurrentPatch = patch;
160 fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
161 fLastRow=AliHLTTPCTransform::GetLastRow(patch);
164 void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt)
166 //set pointer to output
167 fSpacePointData = pt;
170 void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
172 fPtr = (UChar_t*)ptr;
176 void AliHLTTPCClusterFinder::ProcessDigits()
179 bool readValue = true;
182 UShort_t time=0,newTime=0;
183 UInt_t pad=0,newPad=0;
184 AliHLTTPCSignal_t charge=0;
188 // initialize block for reading packed data
189 iResult=fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
190 if (iResult<0) return;
192 readValue = fDigitReader->Next();
194 // Matthias 08.11.2006 the following return would cause termination without writing the
195 // ClusterData and thus would block the component. I just want to have the commented line
196 // here for information
197 //if (!readValue)return;
199 pad = fDigitReader->GetPad();
200 time = fDigitReader->GetTime();
201 fCurrentRow = fDigitReader->GetRow();
203 if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
204 rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
206 fCurrentRow += rowOffset;
208 UInt_t lastpad = 123456789;
209 const UInt_t kPadArraySize=5000;
210 const UInt_t kClusterListSize=10000;
211 AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
212 AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
213 AliClusterData clusterlist[kClusterListSize]; //Clusterlist
215 AliClusterData **currentPt; //List of pointers to the current pad
216 AliClusterData **previousPt; //List of pointers to the previous pad
219 UInt_t nprevious=0,ncurrent=0,ntotal=0;
221 /* quick implementation of baseline calculation and zero suppression
222 open a pad object for each pad and delete it after processing.
223 later a list of pad objects with base line history can be used
224 The whole thing only works if we really get unprocessed raw data, if
225 the data is already zero suppressed, there might be gaps in the time
228 Int_t gatingGridOffset=50;
230 gatingGridOffset=fFirstTimeBin;
232 AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
233 // just to make later conversion to a list of objects easier
234 AliHLTTPCPad* pCurrentPad=NULL;
236 if (fSignalThreshold>=0) {
237 pCurrentPad=&baseline;
238 baseline.SetThreshold(fSignalThreshold);
241 while ( readValue!=0 && iResult>=0){ // Reads through all digits in block
248 if(currentPt == pad2){
256 nprevious = ncurrent;
258 if(pad != lastpad+1){
259 //this happens if there is a pad with no signal.
260 nprevious = ncurrent = 0;
265 Bool_t newcluster = kTRUE;
266 UInt_t seqcharge=0,seqaverage=0,seqerror=0;
267 AliHLTTPCSignal_t lastcharge=0;
268 UInt_t bLastWasFalling=0;
273 redo: //This is a goto.
284 while(iResult>=0){ //Loop over time bins of current pad
285 // read all the values for one pad at once to calculate the base line
287 if (!pCurrentPad->IsStarted()) {
288 //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
289 pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
290 if ((pCurrentPad->StartEvent())>=0) {
292 if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
293 if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
294 pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
295 //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
296 } while ((readValue = fDigitReader->Next())!=0);
298 pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
299 if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
300 //HLTDebug("no data available after zero suppression");
301 pCurrentPad->StopEvent();
302 pCurrentPad->ResetHistory();
305 time=pCurrentPad->GetCurrentPosition();
306 if (time>pCurrentPad->GetSize()) {
307 HLTError("invalid time bin for pad");
313 Float_t occupancy=pCurrentPad->GetOccupancy();
314 //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
315 if ( occupancy < fOccupancyLimit ) {
316 charge = pCurrentPad->GetCorrectedData();
319 //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
322 charge = fDigitReader->GetSignal();
324 //HLTDebug("get next charge value: position %d charge %d", time, charge);
328 if (fDigitReader->GetRow() == 90){
329 ///// LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90") << "PAD=" << fDigitReader->GetPad() << " TIME=" << fDigitReader->GetTime()
330 // << " SIGNAL=" << fDigitReader->GetSignal() << ENDLOG;
334 if(time >= AliHLTTPCTransform::GetNTimeBins()){
335 HLTWarning("Pad %d: Timebin (%d) out of range (%d)", pad, time, AliHLTTPCTransform::GetNTimeBins());
340 //Get the current ADC-value
343 //Check if the last pixel in the sequence is smaller than this
344 if(charge > lastcharge){
350 else bLastWasFalling = 1; //last pixel was larger than this
354 //Sum the total charge of this sequence
356 seqaverage += time*charge;
357 seqerror += time*time*charge;
361 if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
362 pCurrentPad->StopEvent();
363 pCurrentPad->ResetHistory();
365 newPad = fDigitReader->GetPad();
366 newTime = fDigitReader->GetTime();
367 newRow = fDigitReader->GetRow() + rowOffset;
372 newPad=pCurrentPad->GetPadNumber();
373 newTime=pCurrentPad->GetCurrentPosition();
374 newRow=pCurrentPad->GetRowNumber();
376 readValue = fDigitReader->Next();
377 //Check where to stop:
378 if(!readValue) break; //No more value
380 newPad = fDigitReader->GetPad();
381 newTime = fDigitReader->GetTime();
382 newRow = fDigitReader->GetRow() + rowOffset;
385 if(newPad != pad)break; //new pad
386 if(newTime != time+1) break; //end of sequence
389 // pad = newpad; is equal
392 }//end loop over sequence
394 //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
395 //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
397 // with active zero suppression zero values are possible
401 //Calculate mean of sequence:
404 seqmean = seqaverage/seqcharge;
406 LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
407 <<"Error in data given to the cluster finder"<<ENDLOG;
412 //Calculate mean in pad direction:
413 Int_t padmean = seqcharge*pad;
414 Int_t paderror = pad*padmean;
416 //Compare with results on previous pad:
417 for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
419 //dont merge sequences on the same pad twice
420 if(previousPt[p]->fLastMergedPad==pad) continue;
422 Int_t difference = seqmean - previousPt[p]->fMean;
423 if(difference < -fMatch) break;
425 if(difference <= fMatch){ //There is a match here!!
426 AliClusterData *local = previousPt[p];
429 if(seqcharge > local->fLastCharge){
430 if(local->fChargeFalling){ //The previous pad was falling
431 break; //create a new cluster
434 else local->fChargeFalling = 1;
435 local->fLastCharge = seqcharge;
438 //Don't create a new cluster, because we found a match
441 //Update cluster on current pad with the matching one:
442 local->fTotalCharge += seqcharge;
443 local->fPad += padmean;
444 local->fPad2 += paderror;
445 local->fTime += seqaverage;
446 local->fTime2 += seqerror;
447 local->fMean = seqmean;
448 local->fFlags++; //means we have more than one pad
449 local->fLastMergedPad = pad;
451 currentPt[ncurrent] = local;
455 } //Checking for match at previous pad
456 } //Loop over results on previous pad.
458 if(newcluster && ncurrent<kPadArraySize){
459 //Start a new cluster. Add it to the clusterlist, and update
460 //the list of pointers to clusters in current pad.
461 //current pad will be previous pad on next pad.
463 //Add to the clusterlist:
464 AliClusterData *tmp = &clusterlist[ntotal];
465 tmp->fTotalCharge = seqcharge;
467 tmp->fPad2 = paderror;
468 tmp->fTime = seqaverage;
469 tmp->fTime2 = seqerror;
470 tmp->fMean = seqmean;
471 tmp->fFlags = 0; //flags for single pad clusters
472 tmp->fLastMergedPad = pad;
475 tmp->fChargeFalling = 0;
476 tmp->fLastCharge = seqcharge;
479 //Update list of pointers to previous pad:
480 currentPt[ncurrent] = &clusterlist[ntotal];
486 if(newbin >= 0) goto redo;
488 // to prevent endless loop
489 if(time >= AliHLTTPCTransform::GetNTimeBins()){
490 HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
495 if(!readValue) break; //No more value
497 if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
498 HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
502 if(fCurrentRow != newRow){
503 WriteClusters(ntotal,clusterlist);
513 fCurrentRow = newRow;
519 } // END while(readValue)
521 WriteClusters(ntotal,clusterlist);
523 HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
527 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list)
529 //write cluster to output pointer
530 Int_t thisrow=-1,thissector=-1;
531 UInt_t counter = fNClusters;
533 for(int j=0; j<nclusters; j++)
538 if(!list[j].fFlags) continue; //discard single pad clusters
539 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
542 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
543 Float_t fpad2=fXYErr*fXYErr; //fixed given error
544 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
545 Float_t ftime2=fZErr*fZErr; //fixed given error
550 fCurrentRow=list[j].fRow;
554 if(fCalcerr) { //calc the errors, otherwice take the fixed error
555 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
556 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
557 // Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
558 Float_t sy2=(Float_t)list[j].fPad2 * list[j].fTotalCharge - (Float_t)list[j].fPad * list[j].fPad;
561 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
562 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
566 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
568 fpad2*=0.108; //constants are from offline studies
572 } else fpad2=sy2; //take the width not the error
574 // Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
575 Float_t sz2=(Float_t)list[j].fTime2*list[j].fTotalCharge - (Float_t)list[j].fTime*list[j].fTime;
578 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
579 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
583 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
585 ftime2 *= 0.169; //constants are from offline studies
589 } else ftime2=sz2; //take the width, not the error
593 HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
596 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
597 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
599 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
600 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
601 if(fNClusters >= fMaxNClusters)
603 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
604 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
608 fSpacePointData[counter].fX = xyz[0];
609 // fSpacePointData[counter].fY = xyz[1];
610 if(fCurrentSlice<18){
611 fSpacePointData[counter].fY = xyz[1];
614 fSpacePointData[counter].fY = -1*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].fMaxQ = list[j].fQMax;
631 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
632 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
634 Int_t patch=fCurrentPatch;
635 if(patch==-1) patch=0; //never store negative patch number
636 fSpacePointData[counter].fID = counter
637 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
641 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
643 fSpacePointData[counter].fTrackID[0] = trackID[0];
644 fSpacePointData[counter].fTrackID[1] = trackID[1];
645 fSpacePointData[counter].fTrackID[2] = trackID[2];
654 // STILL TO FIX ----------------------------------------------------------------------------
657 void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
660 AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
662 trackID[0]=trackID[1]=trackID[2]=-2;
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];
686 //----------------------------------Methods for the new unsorted way of reading the data --------------------------------
688 void AliHLTTPCClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size)
691 fPtr = (UChar_t*)ptr;
694 if(!fVectorInitialized){
695 InitializePadArray();
698 if (fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice)<0) {
699 HLTError("failed setting up digit reader (InitBlock)");
703 while(fDigitReader->NextChannel()){
704 UInt_t row=fDigitReader->GetRow();
705 UInt_t pad=fDigitReader->GetPad();
707 if(row>=fRowPadVector.size()){
708 HLTError("Row number is to large: %d, max is %d",row,fRowPadVector.size()-1);
711 if(pad>=fRowPadVector[row].size()){
712 HLTError("Pad number is to large: %d, max is %d",pad,fRowPadVector[row].size());
716 while(fDigitReader->NextBunch()){
717 if(fDigitReader->GetBunchSize()>1){//to remove single timebin values, this will have to change at some point
718 UInt_t time = fDigitReader->GetTime();
719 if((Int_t)time>=fFirstTimeBin && (Int_t)time+fDigitReader->GetBunchSize()<=fLastTimeBin){
720 const UInt_t *bunchData= fDigitReader->GetSignals();
721 AliHLTTPCClusters candidate;
722 for(Int_t i=0;i<fDigitReader->GetBunchSize();i++){
723 candidate.fTotalCharge+=bunchData[i];
724 candidate.fTime += time*bunchData[i];
725 candidate.fTime2 += time*time*bunchData[i];
726 if(bunchData[i]>candidate.fQMax){
727 candidate.fQMax=bunchData[i];
731 if(candidate.fTotalCharge>0){
732 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
733 candidate.fPad=candidate.fTotalCharge*pad;
734 candidate.fPad2=candidate.fPad*pad;
735 candidate.fLastMergedPad=pad;
736 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
738 if(fRowPadVector[row][pad] != NULL){
739 fRowPadVector[row][pad]->AddClusterCandidate(candidate);
747 void AliHLTTPCClusterFinder::ReadDataUnsortedDeconvoluteTime(void* ptr,unsigned long size)
750 fPtr = (UChar_t*)ptr;
753 if(!fVectorInitialized){
754 InitializePadArray();
757 if (fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice)<0) {
758 HLTError("failed setting up digit reader (InitBlock)");
762 while(fDigitReader->NextChannel()){
763 UInt_t row=fDigitReader->GetRow();
764 UInt_t pad=fDigitReader->GetPad();
766 while(fDigitReader->NextBunch()){
767 if(fDigitReader->GetBunchSize()>1){//to remove single timebin values, this will have to change at some point
768 UInt_t time = fDigitReader->GetTime();
769 if((Int_t)time>=fFirstTimeBin && (Int_t)time+fDigitReader->GetBunchSize()<=fLastTimeBin){
770 Int_t indexInBunchData=0;
771 Bool_t moreDataInBunch=kFALSE;
773 Bool_t signalFalling=kFALSE;
774 const UInt_t *bunchData= fDigitReader->GetSignals();
776 AliHLTTPCClusters candidate;
777 for(Int_t i=indexInBunchData;i<fDigitReader->GetBunchSize();i++){
778 // Checks if one need to deconvolute the signals
779 if(bunchData[i]>prevSignal && signalFalling==kTRUE){
780 if(i<fDigitReader->GetBunchSize()-1){ // means there are more than one signal left in the bunch
781 moreDataInBunch=kTRUE;
787 // Checks if the signal is 0, then quit processing the data.
788 if(bunchData[i]==0 && i<fDigitReader->GetBunchSize()-1){//means we have 0 data fom the rcu, might happen depending on the configuration settings
789 moreDataInBunch=kTRUE;
794 if(prevSignal>bunchData[i]){//means the peak of the signal has been reached and deconvolution will happen if the signal rise again.
797 candidate.fTotalCharge+=bunchData[i];
798 candidate.fTime += time*bunchData[i];
799 candidate.fTime2 += time*time*bunchData[i];
800 if(bunchData[i]>candidate.fQMax){
801 candidate.fQMax=bunchData[i];
804 prevSignal=bunchData[i];
808 if(candidate.fTotalCharge>0){
809 candidate.fMean=candidate.fTime/candidate.fTotalCharge;
810 candidate.fPad=candidate.fTotalCharge*pad;
811 candidate.fPad2=candidate.fPad*pad;
812 candidate.fLastMergedPad=pad;
813 candidate.fRowNumber=row+fDigitReader->GetRowOffset();
815 fRowPadVector[row][pad]->AddClusterCandidate(candidate);
816 if(indexInBunchData<fDigitReader->GetBunchSize()-1){
817 moreDataInBunch=kFALSE;
819 }while(moreDataInBunch);
826 Bool_t AliHLTTPCClusterFinder::ComparePads(AliHLTTPCPad *nextPad,AliHLTTPCClusters* cluster,Int_t nextPadToRead){
827 //Checking if we have a match on the next pad
828 for(UInt_t candidateNumber=0;candidateNumber<nextPad->fClusterCandidates.size();candidateNumber++){
829 AliHLTTPCClusters *candidate =&nextPad->fClusterCandidates[candidateNumber];
830 if(cluster->fMean-candidate->fMean==1 || candidate->fMean-cluster->fMean==1 || cluster->fMean-candidate->fMean==0){
832 if(candidate->fTotalCharge<fTotalChargeOfPreviousClusterCandidate){//peak is reached
833 fChargeOfCandidatesFalling=kTRUE;
835 if(candidate->fTotalCharge>fTotalChargeOfPreviousClusterCandidate && fChargeOfCandidatesFalling==kTRUE){//we have deconvolution
839 cluster->fMean=candidate->fMean;
840 cluster->fTotalCharge+=candidate->fTotalCharge;
841 cluster->fTime += candidate->fTime;
842 cluster->fTime2 += candidate->fTime2;
843 cluster->fPad+=candidate->fPad;
844 cluster->fPad2+=candidate->fPad2;
845 cluster->fLastMergedPad=candidate->fPad;
846 if(candidate->fQMax>cluster->fQMax){
847 cluster->fQMax=candidate->fQMax;
851 UInt_t rowNo = nextPad->GetRowNumber();
852 UInt_t padNo = nextPad->GetPadNumber();
854 fRowPadVector[rowNo][padNo-2]->fSelectedPad=kTRUE;
855 fRowPadVector[rowNo][padNo-2]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo-2);
857 fRowPadVector[rowNo][padNo-1]->fSelectedPad=kTRUE;// quick solution to set the first pad to selected
858 fRowPadVector[rowNo][padNo-1]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo-1);
859 fRowPadVector[rowNo][padNo]->fSelectedPad=kTRUE;
860 fRowPadVector[rowNo][padNo]->fHWAddress=(AliHLTUInt16_t)fDigitReader->GetAltroBlockHWaddr(rowNo,padNo);
863 //setting the matched pad to used
864 nextPad->fUsedClusterCandidates[candidateNumber]=1;
866 if(nextPadToRead<(Int_t)fNumberOfPadsInRow[fRowOfFirstCandidate]){
867 nextPad=fRowPadVector[fRowOfFirstCandidate][nextPadToRead];
868 ComparePads(nextPad,cluster,nextPadToRead);
881 Int_t AliHLTTPCClusterFinder::FillHWAddressList(AliHLTUInt16_t *hwaddlist, Int_t maxHWadd){
883 for(UInt_t row=0;row<fNumberOfRows;row++){
884 for(UInt_t pad=0;pad<fNumberOfPadsInRow[row]-1;pad++){
885 if(fRowPadVector[row][pad]->fSelectedPad){
886 if(counter<maxHWadd){
887 hwaddlist[counter]=(AliHLTUInt16_t)fRowPadVector[row][pad]->fHWAddress;
891 HLTWarning("To many hardwareaddresses, skip adding");
900 void AliHLTTPCClusterFinder::FindClusters()
902 // see header file for function documentation
904 AliHLTTPCClusters* tmpCandidate=NULL;
905 for(UInt_t row=0;row<fNumberOfRows;row++){
906 fRowOfFirstCandidate=row;
907 for(UInt_t pad=0;pad<fNumberOfPadsInRow[row];pad++){
908 AliHLTTPCPad *tmpPad=fRowPadVector[row][pad];
909 for(size_t candidate=0;candidate<tmpPad->fClusterCandidates.size();candidate++){
910 if(tmpPad->fUsedClusterCandidates[candidate]){
913 tmpCandidate=&tmpPad->fClusterCandidates[candidate];
914 UInt_t tmpTotalCharge=tmpCandidate->fTotalCharge;
916 ComparePads(fRowPadVector[row][pad+1],tmpCandidate,pad+1);
917 if(tmpCandidate->fTotalCharge>tmpTotalCharge){
919 fClusters.push_back(*tmpCandidate);
922 tmpPad->ClearCandidates();
924 fRowPadVector[row][fNumberOfPadsInRow[row]]->ClearCandidates();
927 HLTInfo("Found %d clusters.",fClusters.size());
929 //TODO: Change so it stores AliHLTTPCSpacePointData directly, instead of this copying
931 AliClusterData * clusterlist = new AliClusterData[fClusters.size()]; //Clusterlist
932 for(unsigned int i=0;i<fClusters.size();i++){
933 clusterlist[i].fTotalCharge = fClusters[i].fTotalCharge;
934 clusterlist[i].fPad = fClusters[i].fPad;
935 clusterlist[i].fPad2 = fClusters[i].fPad2;
936 clusterlist[i].fTime = fClusters[i].fTime;
937 clusterlist[i].fTime2 = fClusters[i].fTime2;
938 clusterlist[i].fMean = fClusters[i].fMean;
939 clusterlist[i].fFlags = fClusters[i].fFlags;
940 clusterlist[i].fChargeFalling = fClusters[i].fChargeFalling;
941 clusterlist[i].fLastCharge = fClusters[i].fLastCharge;
942 clusterlist[i].fLastMergedPad = fClusters[i].fLastMergedPad;
943 clusterlist[i].fRow = fClusters[i].fRowNumber;
946 WriteClusters(fClusters.size(),clusterlist);
947 delete [] clusterlist;
951 void AliHLTTPCClusterFinder::PrintClusters()
953 // see header file for class documentation
954 for(size_t i=0;i<fClusters.size();i++){
955 HLTInfo("Cluster number: %d",i);
956 HLTInfo("Row: %d \t Pad: %d",fClusters[i].fRowNumber,fClusters[i].fPad/fClusters[i].fTotalCharge);
957 HLTInfo("Total Charge: %d",fClusters[i].fTotalCharge);
958 HLTInfo("fPad: %d",fClusters[i].fPad);
959 HLTInfo("PadError: %d",fClusters[i].fPad2);
960 HLTInfo("TimeMean: %d",fClusters[i].fTime/fClusters[i].fTotalCharge);
961 HLTInfo("TimeError: %d",fClusters[i].fTime2);
962 HLTInfo("EndOfCluster:");
966 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list)//This is used when using the AliHLTTPCClusters class for cluster data
968 //write cluster to output pointer
969 Int_t thisrow,thissector;
970 UInt_t counter = fNClusters;
972 for(int j=0; j<nclusters; j++)
974 if(!list[j].fFlags) continue; //discard single pad clusters
975 if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
978 Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
979 Float_t fpad2=fXYErr*fXYErr; //fixed given error
980 Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
981 Float_t ftime2=fZErr*fZErr; //fixed given error
984 if(fCalcerr) { //calc the errors, otherwice take the fixed error
985 Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
986 UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
987 Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
990 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
991 <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
995 fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
997 fpad2*=0.108; //constants are from offline studies
1001 } else fpad2=sy2; //take the width not the error
1003 Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
1006 LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
1007 <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
1011 ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
1013 ftime2 *= 0.169; //constants are from offline studies
1017 } else ftime2=sz2; //take the width, not the error
1021 HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
1024 AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
1025 AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
1027 if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
1028 <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
1029 if(fNClusters >= fMaxNClusters)
1031 LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
1032 <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
1036 fSpacePointData[counter].fX = xyz[0];
1037 // fSpacePointData[counter].fY = xyz[1];
1038 if(fCurrentSlice<18){
1039 fSpacePointData[counter].fY = xyz[1];
1042 fSpacePointData[counter].fY = -1*xyz[1];
1044 fSpacePointData[counter].fZ = xyz[2];
1047 fSpacePointData[counter].fX = fCurrentRow;
1048 fSpacePointData[counter].fY = fpad;
1049 fSpacePointData[counter].fZ = ftime;
1052 fSpacePointData[counter].fCharge = list[j].fTotalCharge;
1053 fSpacePointData[counter].fPadRow = fCurrentRow;
1054 fSpacePointData[counter].fSigmaY2 = fpad2;
1055 fSpacePointData[counter].fSigmaZ2 = ftime2;
1057 fSpacePointData[counter].fMaxQ = list[j].fQMax;
1059 fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
1060 fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
1062 Int_t patch=fCurrentPatch;
1063 if(patch==-1) patch=0; //never store negative patch number
1064 fSpacePointData[counter].fID = counter
1065 +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
1069 GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
1071 fSpacePointData[counter].fTrackID[0] = trackID[0];
1072 fSpacePointData[counter].fTrackID[1] = trackID[1];
1073 fSpacePointData[counter].fTrackID[2] = trackID[2];