// @(#) $Id$
-// Original: AliL3ClustFinderNew.cxx,v 1.29 2005/06/14 10:55:21 cvetan Exp
-
-// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>,
-// Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
-// Jochen Thaeder <mailto:thaeder@kip.uni-heidelberg.de>
-
-//*-- Copyright © ALICE HLT Group
-
+// Original: AliHLTClustFinderNew.cxx,v 1.29 2005/06/14 10:55:21 cvetan Exp
+
+/**************************************************************************
+ * This file is property of and copyright by the ALICE HLT Project *
+ * ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Primary Authors: Anders Vestbo, Constantin Loizides, Jochen Thaeder *
+ * Kenneth Aamodt <kenneth.aamodt@student.uib.no> *
+ * for The ALICE HLT Project. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/** @file AliHLTTPCClusterFinder.cxx
+ @author Anders Vestbo, Constantin Loizides, Jochen Thaeder
+ Kenneth Aamodt kenneth.aamodt@student.uib.no
+ @date
+ @brief Cluster Finder for the TPC
+*/
#include "AliHLTTPCDigitReader.h"
#include "AliHLTTPCRootTypes.h"
#include "AliHLTTPCSpacePointData.h"
#include "AliHLTTPCMemHandler.h"
#include "AliHLTTPCPad.h"
+#include "AliHLTTPCPadArray.h"
#if __GNUC__ >= 3
using namespace std;
#endif
/** \class AliHLTTPCClusterFinder
-<pre>
-//_____________________________________________________________
-// AliHLTTPCClusterFinder
//
// The current cluster finder for HLT
// (Based on STAR L3)
-//
+//
+//Basically we have two versions for the cluster finder now.
+//The default version, reads the data pad by pad, and find the
+//clusters as it reads the data. The other version has now been
+//developed to cope with unsorted data. New methods for the unsorted
+//version can be found at the end of the default one i the source file.
+//Currently the new version is only build to manage zero-suppressed data.
+//More functionality will be added later.
+//
// The cluster finder is initialized with the Init function,
// providing the slice and patch information to work on.
//
// delete cf;
// }
// }
-</pre>
*/
ClassImp(AliHLTTPCClusterFinder)
AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
:
- fMatch(1),
- fThreshold(10),
- fXYErr(0.2),
- fZErr(0.3),
- fDeconvPad(kTRUE),
+ fSpacePointData(NULL),
+ fDigitReader(NULL),
+ fPtr(NULL),
+ fSize(0),
fDeconvTime(kTRUE),
+ fDeconvPad(kTRUE),
fStdout(kFALSE),
fCalcerr(kTRUE),
fRawSP(kFALSE),
fFirstRow(0),
fLastRow(0),
- fDigitReader(NULL)
+ fCurrentRow(0),
+ fCurrentSlice(0),
+ fCurrentPatch(0),
+ fMatch(1),
+ fThreshold(10),
+ fSignalThreshold(-1),
+ fNClusters(0),
+ fMaxNClusters(0),
+ fXYErr(0.2),
+ fZErr(0.3),
+ fOccupancyLimit(1.0),
+ fPadArray(NULL)
{
//constructor
}
AliHLTTPCClusterFinder::AliHLTTPCClusterFinder(const AliHLTTPCClusterFinder& src)
:
- fMatch(src.fMatch),
- fThreshold(src.fThreshold),
- fXYErr(src.fXYErr),
- fZErr(src.fZErr),
- fDeconvPad(src.fDeconvPad),
+ fSpacePointData(NULL),
+ fDigitReader(src.fDigitReader),
+ fPtr(NULL),
+ fSize(src.fSize),
fDeconvTime(src.fDeconvTime),
+ fDeconvPad(src.fDeconvPad),
fStdout(src.fStdout),
fCalcerr(src.fCalcerr),
fRawSP(src.fRawSP),
fFirstRow(src.fFirstRow),
fLastRow(src.fLastRow),
- fDigitReader(src.fDigitReader)
+ fCurrentRow(src.fCurrentRow),
+ fCurrentSlice(src.fCurrentSlice),
+ fCurrentPatch(src.fCurrentPatch),
+ fMatch(src.fMatch),
+ fThreshold(src.fThreshold),
+ fSignalThreshold(src.fSignalThreshold),
+ fNClusters(src.fNClusters),
+ fMaxNClusters(src.fMaxNClusters),
+ fXYErr(src.fXYErr),
+ fZErr(src.fZErr),
+ fOccupancyLimit(src.fOccupancyLimit),
+ fPadArray(NULL)
{
}
{
fMatch=src.fMatch;
fThreshold=src.fThreshold;
+ fSignalThreshold=src.fSignalThreshold;
+ fOccupancyLimit=src.fOccupancyLimit;
fXYErr=src.fXYErr;
fZErr=src.fZErr;
fDeconvPad=src.fDeconvPad;
void AliHLTTPCClusterFinder::ProcessDigits()
{
+ int iResult=0;
bool readValue = true;
Int_t newRow = 0;
Int_t rowOffset = 0;
fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
readValue = fDigitReader->Next();
- if (!readValue)return;
+ // Matthias 08.11.2006 the following return would cause termination without writing the
+ // ClusterData and thus would block the component. I just want to have the commented line
+ // here for information
+ //if (!readValue)return;
pad = fDigitReader->GetPad();
time = fDigitReader->GetTime();
fCurrentRow += rowOffset;
UInt_t lastpad = 123456789;
- AliClusterData *pad1[5000]; //2 lists for internal memory=2pads
- AliClusterData *pad2[5000]; //2 lists for internal memory=2pads
- AliClusterData clusterlist[10000]; //Clusterlist
+ const UInt_t kPadArraySize=5000;
+ const UInt_t kClusterListSize=10000;
+ AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
+ AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
+ AliClusterData clusterlist[kClusterListSize]; //Clusterlist
AliClusterData **currentPt; //List of pointers to the current pad
AliClusterData **previousPt; //List of pointers to the previous pad
Int_t gatingGridOffset=50;
AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
// just to make later conversion to a list of objects easier
- AliHLTTPCPad* pCurrentPad=&baseline;
+ AliHLTTPCPad* pCurrentPad=NULL;
+ if (fSignalThreshold>=0) {
+ pCurrentPad=&baseline;
+ baseline.SetThreshold(fSignalThreshold);
+ }
- while ( readValue ){ // Reads through all digits in block
+ while ( readValue!=0 && iResult>=0){ // Reads through all digits in block
+ iResult=0;
if(pad != lastpad){
//This is a new pad
Bool_t newcluster = kTRUE;
UInt_t seqcharge=0,seqaverage=0,seqerror=0;
- UInt_t lastcharge=0,lastwas_falling=0;
+ AliHLTTPCSignal_t lastcharge=0;
+ UInt_t bLastWasFalling=0;
Int_t newbin=-1;
}
lastcharge=0;
- lastwas_falling = 0;
+ bLastWasFalling = 0;
}
- while(1){ //Loop over time bins of current pad
+ while(iResult>=0){ //Loop over time bins of current pad
// read all the values for one pad at once to calculate the base line
if (pCurrentPad) {
if (!pCurrentPad->IsStarted()) {
}
pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
- HLTDebug("no data available after zero suppression");
+ //HLTDebug("no data available after zero suppression");
pCurrentPad->StopEvent();
pCurrentPad->ResetHistory();
break;
}
if (pCurrentPad) {
- charge = pCurrentPad->GetCorrectedData();
+ Float_t occupancy=pCurrentPad->GetOccupancy();
+ //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
+ if ( occupancy < fOccupancyLimit ) {
+ charge = pCurrentPad->GetCorrectedData();
+ } else {
+ charge = 0;
+ //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
+ }
} else {
- charge = fDigitReader->GetSignal();
+ charge = fDigitReader->GetSignal();
}
//HLTDebug("get next charge value: position %d charge %d", time, charge);
}
if(time >= AliHLTTPCTransform::GetNTimeBins()){
- HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
- break;
+ HLTWarning("Pad %d: Timebin (%d) out of range (%d)", pad, time, AliHLTTPCTransform::GetNTimeBins());
+ iResult=-ERANGE;
}
//Check if the last pixel in the sequence is smaller than this
if(charge > lastcharge){
- if(lastwas_falling){
+ if(bLastWasFalling){
newbin = 1;
break;
}
}
- else lastwas_falling = 1; //last pixel was larger than this
+ else bLastWasFalling = 1; //last pixel was larger than this
lastcharge = charge;
}
if(newPad != pad)break; //new pad
if(newTime != time+1) break; //end of sequence
+ if(iResult<0) break;
// pad = newpad; is equal
time = newTime;
Int_t padmean = seqcharge*pad;
Int_t paderror = pad*padmean;
-
//Compare with results on previous pad:
- for(UInt_t p=0; p<nprevious; p++){
+ for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
//dont merge sequences on the same pad twice
if(previousPt[p]->fLastMergedPad==pad) continue;
} //Checking for match at previous pad
} //Loop over results on previous pad.
-
- if(newcluster){
+ if(newcluster && ncurrent<kPadArraySize){
//Start a new cluster. Add it to the clusterlist, and update
//the list of pointers to clusters in current pad.
//current pad will be previous pad on next pad.
if(!readValue) break; //No more value
+
+ if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
+ HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
+ break;
+ }
if(fCurrentRow != newRow){
WriteClusters(ntotal,clusterlist);
} // END while(readValue)
- if (pCurrentPad) pCurrentPad->StopEvent();
-
WriteClusters(ntotal,clusterlist);
HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
for(int j=0; j<nclusters; j++)
{
+
+
+
if(!list[j].fFlags) continue; //discard single pad clusters
if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
Float_t ftime2=fZErr*fZErr; //fixed given error
+#if UNSORTED
+ fCurrentRow=list[j].fRow;
+#endif
-
-
-
if(fCalcerr) { //calc the errors, otherwice take the fixed error
Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
}
}
#endif
+
+//----------------------------------Methods for the new unsorted way of reading the data --------------------------------
+
+void AliHLTTPCClusterFinder::SetPadArray(AliHLTTPCPadArray * padArray){
+ fPadArray=padArray;
+}
+
+void AliHLTTPCClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size){
+ //set input pointer
+ fPtr = (UChar_t*)ptr;
+ fSize = size;
+
+ fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
+
+ fPadArray->SetDigitReader(fDigitReader);
+ fPadArray->ReadData();
+}
+void AliHLTTPCClusterFinder::FindClusters(){
+ fPadArray->FindClusterCandidates();
+ fPadArray->FindClusters(fMatch);
+
+ AliClusterData * clusterlist = new AliClusterData[fPadArray->fClusters.size()]; //Clusterlist
+ for(unsigned int i=0;i<fPadArray->fClusters.size();i++){
+ clusterlist[i].fTotalCharge = fPadArray->fClusters[i].fTotalCharge;
+ clusterlist[i].fPad = fPadArray->fClusters[i].fPad;
+ clusterlist[i].fPad2 = fPadArray->fClusters[i].fPad2;
+ clusterlist[i].fTime = fPadArray->fClusters[i].fTime;
+ clusterlist[i].fTime2 = fPadArray->fClusters[i].fTime2;
+ clusterlist[i].fMean = fPadArray->fClusters[i].fMean;
+ clusterlist[i].fFlags = fPadArray->fClusters[i].fFlags;
+ clusterlist[i].fChargeFalling = fPadArray->fClusters[i].fChargeFalling;
+ clusterlist[i].fLastCharge = fPadArray->fClusters[i].fLastCharge;
+ clusterlist[i].fLastMergedPad = fPadArray->fClusters[i].fLastMergedPad;
+ clusterlist[i].fRow = fPadArray->fClusters[i].fRowNumber;
+ }
+ WriteClusters(fPadArray->fClusters.size(),clusterlist);
+ delete [] clusterlist;
+ fPadArray->DataToDefault();
+}
+void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list)//This is used when using the AliHLTTPCClusters class for cluster data
+{
+ //write cluster to output pointer
+ Int_t thisrow,thissector;
+ UInt_t counter = fNClusters;
+
+ for(int j=0; j<nclusters; j++)
+ {
+ if(!list[j].fFlags) continue; //discard single pad clusters
+ if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
+
+ Float_t xyz[3];
+ Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
+ Float_t fpad2=fXYErr*fXYErr; //fixed given error
+ Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
+ Float_t ftime2=fZErr*fZErr; //fixed given error
+
+
+ if(fCalcerr) { //calc the errors, otherwice take the fixed error
+ Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
+ UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
+ Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
+ sy2/=q2;
+ if(sy2 < 0) {
+ LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
+ <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
+ continue;
+ } else {
+ if(!fRawSP){
+ fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
+ if(sy2 != 0){
+ fpad2*=0.108; //constants are from offline studies
+ if(patch<2)
+ fpad2*=2.07;
+ }
+ } else fpad2=sy2; //take the width not the error
+ }
+ Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
+ sz2/=q2;
+ if(sz2 < 0){
+ LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
+ <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
+ continue;
+ } else {
+ if(!fRawSP){
+ ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
+ if(sz2 != 0) {
+ ftime2 *= 0.169; //constants are from offline studies
+ if(patch<2)
+ ftime2 *= 1.77;
+ }
+ } else ftime2=sz2; //take the width, not the error
+ }
+ }
+ if(fStdout==kTRUE)
+ cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
+
+ if(!fRawSP){
+ AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
+ AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
+
+ if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
+ <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
+ if(fNClusters >= fMaxNClusters)
+ {
+ LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
+ <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
+ return;
+ }
+
+ fSpacePointData[counter].fX = xyz[0];
+ fSpacePointData[counter].fY = xyz[1];
+ fSpacePointData[counter].fZ = xyz[2];
+
+ } else {
+ fSpacePointData[counter].fX = fCurrentRow;
+ fSpacePointData[counter].fY = fpad;
+ fSpacePointData[counter].fZ = ftime;
+ }
+
+ fSpacePointData[counter].fCharge = list[j].fTotalCharge;
+ fSpacePointData[counter].fPadRow = fCurrentRow;
+ fSpacePointData[counter].fSigmaY2 = fpad2;
+ fSpacePointData[counter].fSigmaZ2 = ftime2;
+
+ fSpacePointData[counter].fUsed = kFALSE; // only used / set in AliHLTTPCDisplay
+ fSpacePointData[counter].fTrackN = -1; // only used / set in AliHLTTPCDisplay
+
+ Int_t patch=fCurrentPatch;
+ if(patch==-1) patch=0; //never store negative patch number
+ fSpacePointData[counter].fID = counter
+ +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
+
+#ifdef do_mc
+ Int_t trackID[3];
+ GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
+
+ fSpacePointData[counter].fTrackID[0] = trackID[0];
+ fSpacePointData[counter].fTrackID[1] = trackID[1];
+ fSpacePointData[counter].fTrackID[2] = trackID[2];
+
+ //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
+#endif
+
+ fNClusters++;
+ counter++;
+ }
+}