--- /dev/null
+/* $Id$
+// Author: Constantin Loizides <mailto:loizides@ikf.physik.uni-frankfurt.de>
+*/
+
+#include <math.h>
+#include <time.h>
+#include <iostream.h>
+#include <stdio.h>
+#include <string.h>
+
+#include "AliL3VHDLClusterFinder.h"
+
+/** \class AliL3VHDLClusterFinder
+//<pre>
+//____________________________________________________
+// AliL3VHDLClusterFinder
+//
+// The current VHDL cluster finder for HLT
+// Based on STAR L3
+//
+// Most important parameters:
+// fThreshold - threshold for noise clusters
+// fMatch - length in time for overlapping sequences
+//</pre> */
+
+ClassImp(AliL3VHDLClusterFinder)
+
+AliL3VHDLClusterFinder::AliL3VHDLClusterFinder()
+{
+ fMatch = 4;
+ fThreshold = 10;
+ fMinMerge = 1;
+ fNClusters=0;
+
+ fXYErr = 0.2;
+ fZErr = 0.3;
+ fDeconvPad = kTRUE;
+ fDeconvTime = kTRUE;
+ fstdout = kFALSE;
+ fcalcerr = kTRUE;
+
+ Clear();
+#ifdef DEBUG
+ fdeb=fopen("vhdlclusterfinder.debug","w");
+ //fdeb=stderr;
+#endif
+}
+
+AliL3VHDLClusterFinder::~AliL3VHDLClusterFinder()
+{
+#ifdef DEBUG
+ fclose(fdeb);
+#endif
+}
+
+void AliL3VHDLClusterFinder::ProcessDigits()
+{
+ //Loop over data like the analyzer of the VHDL code
+ const UChar_t n=255;
+ UShort_t rrow=0,rtime=0;
+ UChar_t rpad=0,i=n;
+ UShort_t *charges=new UShort_t[n];
+ Int_t tc=0,mp=0,mt=0,sp=0,st=0;
+
+ fNClusters=0;
+ fRow=0;
+ fPad=0;
+ Clear();
+
+ //loop over input data
+ while(fAltromem.ReadSequence(rrow,rpad,rtime,i,&charges)){
+ tc=0;mp=0;mt=0;sp=0;st=0;
+#if 0
+ cout << "Padrow " << (int)rrow << " pad " << (int)rpad << " time " <<(int)rtime << " charges ";
+ for(UChar_t ii=0;ii<i;ii++) cout << (int)charges[ii] << " ";
+ cout << endl;
+#endif
+#ifdef DEBUG
+ fprintf(fdeb,"ProcessDigits: Input Data: %d %d %d charges:",(int)rrow,(int)rpad,(int)rtime);
+ for(UChar_t ii=0;ii<i;ii++) fprintf(fdeb," %d",(int)charges[ii]);
+ fprintf(fdeb,"\n");
+#endif
+
+ fNRow=rrow;
+ fNPad=rpad;
+
+ //calculate sequence values
+ //no deconvulution so far
+ for(UChar_t ii=0;ii<i;ii++){
+ tc+=charges[ii];
+ mt+=(rtime-ii)*charges[ii];
+ st+=(rtime-ii)*(rtime-ii)*charges[ii];
+ }
+ mp=rpad*tc;
+ sp=rpad*rpad*tc;
+
+ fSeq.fTotalCharge=tc;
+ fSeq.fPad=mp;
+ fSeq.fTime=mt;
+ fSeq.fPad2=sp;
+ fSeq.fTime2=st;
+ fSeq.fMean=0;
+ if(tc!=0) fSeq.fMean=mt/tc;
+ fSeq.fMerge=0;
+ fSeq.fRow=rrow;
+ fSeq.fLastPad=rpad;
+
+ //work on this sequence
+ ProcessSequence();
+ //output one cluster
+ OutputMemory();
+
+#ifdef DEBUG
+ fflush(fdeb);
+#endif
+ i=n; //store size of charges array
+ } //loop over data
+
+ //flush everything left
+ FlushMemory();
+ while(fOP!=fFP) OutputMemory();
+}
+
+void AliL3VHDLClusterFinder::ProcessSequence()
+{
+ if(fNRow!=fRow) FlushMemory();
+ else if(fNPad==fPad+1) PrepareMemory();
+ else if(fNPad!=fPad) FlushMemory();
+
+ //store new row and pad values
+ fRow=fNRow;
+ fPad=fNPad;
+
+#ifdef DEBUG
+ fprintf(fdeb,"ProcessSequence: Mean=%d TC=%d ",fSeq.fMean,fSeq.fTotalCharge);
+#endif
+
+ CompareSeq(); //merge or insert
+}
+
+void AliL3VHDLClusterFinder::PrepareMemory()
+{
+#ifdef DEBUG
+ fprintf(fdeb,"PrepareMemory %d %d %d\n",fRP,fEP,fWP);
+#endif
+ fRP=fEP;
+ fEP=fWP;
+}
+
+void AliL3VHDLClusterFinder::FlushMemory()
+{
+#ifdef DEBUG
+ fprintf(fdeb,"FlushMemory %d %d %d %d\n",fFP,fRP,fEP,fWP);
+#endif
+ fFP=fWP;
+ fRP=fWP;
+ fEP=fWP;
+}
+
+void AliL3VHDLClusterFinder::CompareSeq()
+{
+
+ while(fRP!=fEP){
+ Int_t diff=fSeqs[fPList[fRP]].fMean-fSeq.fMean;
+
+ if(diff>fMatch){ //no match
+ IncRPointer(); //cluster finished
+ continue;
+ } else if(diff<-fMatch){ //no match
+ break; //insert new cluster
+ }
+ else { //match found, merge it
+ MergeSeq();
+ return;
+ }
+ }
+
+ InsertSeq(); //start new cluster
+}
+
+void AliL3VHDLClusterFinder::MergeSeq()
+{
+#ifdef DEBUG
+ fprintf(fdeb,"merged with Mean=%d TC=%d (new Merge=%d)\n",fSeqs[fPList[fRP]].fMean,fSeqs[fPList[fRP]].fTotalCharge,fSeqs[fPList[fRP]].fMerge+1);
+#endif
+ if(fSeqs[fPList[fRP]].fRow==fSeq.fRow){
+ LOG(AliL3Log::kWarning,"AliL3VHDLClusterFinder::","Memory Check")
+ <<"Sequences can be merged on the same rows only."<<ENDLOG;
+ }
+ if(fSeqs[fPList[fRP]].fLastPad+1!=fSeq.fLastPad){
+ LOG(AliL3Log::kWarning,"AliL3VHDLClusterFinder::","Memory Check")
+ <<"Sequences can be merged on consecutive pads only."<<ENDLOG;
+ }
+
+ fSeqs[fPList[fRP]].fMean=fSeq.fMean; //take the new mean
+ fSeqs[fPList[fRP]].fLastPad=fSeq.fLastPad;
+ fSeqs[fPList[fRP]].fTotalCharge+=fSeq.fTotalCharge;
+ fSeqs[fPList[fRP]].fPad+=fSeq.fPad;
+ fSeqs[fPList[fRP]].fTime+=fSeq.fTime;
+ fSeqs[fPList[fRP]].fPad2+=fSeq.fPad2;
+ fSeqs[fPList[fRP]].fTime2+=fSeq.fTime2;
+ fSeqs[fPList[fRP]].fMerge+=1;
+
+ fPList[fWP]=fPList[fRP];
+ fPList[fRP]=N_clmem; //mark it to be free
+
+ IncRPointer();
+ IncWPointer();
+}
+
+void AliL3VHDLClusterFinder::InsertSeq()
+{
+#ifdef DEBUG
+ fprintf(fdeb,"inserted\n");
+#endif
+ NextFreeIndex(); //get next index
+ fSeqs[fFirst]=fSeq; //store data
+ fPList[fWP]=fFirst; //store index
+ IncWPointer();
+}
+
+void AliL3VHDLClusterFinder::OutputMemory()
+{
+ Float_t mtime=0,mpad=0;
+ Float_t mtime2=0,mpad2=0;
+ UInt_t tc,row,mno;
+
+ if(fOP!=fFP){
+ //while(fOP!=fFP){
+
+ UInt_t index=fPList[fOP];
+ fPList[fOP]=N_clmem;
+ //cout << fOP << " - " << fFP << " - " << index << endl;
+ IncPointer(fOP);
+ if(index>=N_clmem) return; //nothing to do
+ //if(index>=N_clmem) continue; //nothing to do
+
+ tc=fSeqs[index].fTotalCharge;
+ mno=fSeqs[index].fMerge;
+ row=fSeqs[index].fRow;
+ if(tc!=0){
+ mtime=(Float_t)(fSeqs[index].fTime)/tc;
+ mtime2=sqrt((Float_t)(fSeqs[index].fTime2)/tc-mtime*mtime);
+ }
+ if(tc!=0){
+ mpad=(Float_t)(fSeqs[index].fPad)/tc;
+ mpad2=sqrt((Float_t)(fSeqs[index].fPad2)/tc-mpad*mpad);
+ }
+
+ if(mno<fMinMerge){ //noise cluster
+#ifdef DEBUG
+ fprintf(fdeb,"OutputMemory %d %d: noise cluster (merge): Row %d Pad %.3f Time %.3f TC %d Merge %d\n",fOP,fFP,row,mpad,mtime,tc,mno);
+#endif
+ FreeSeq(index);
+ return;
+ //continue;
+ }
+
+ if(tc<fThreshold){ //total charge below threshold
+#ifdef DEBUG
+ fprintf(fdeb,"OutputMemory %d %d: noise cluster (threshold): Row %d Pad %.3f Time %.3f TC %d Merge %d\n",fOP,fFP,row,mpad,mtime,tc,mno);
+#endif
+ FreeSeq(index);
+ return;
+ //continue;
+ }
+
+#ifdef DEBUG
+ fprintf(fdeb,"OutputMemory %d: Row %d Pad %.3f Time %.3f TC %d Merge %d\n",fNClusters,row,mpad,mtime,tc,mno);
+#endif
+
+ if(stdout)
+ cout<<"WriteCluster: padrow "<<row<<" pad "<<mpad<< " +- "<<mpad2<<" time "<<mtime<<" +- "<<mtime2<<" charge "<<tc<<endl;
+
+ fNClusters++;
+ FreeSeq(index);
+ }
+}
+
+inline void AliL3VHDLClusterFinder::FreeSeq(UShort_t i)
+{
+ ClearSeq(i);
+ IncPointer(fLast,1,N_clmem);
+}
+
+void AliL3VHDLClusterFinder::Clear()
+{
+ fFirst=0; //first list pointer
+ fLast=0; //last list pointer
+
+ fWP=0; //write pointer
+ fRP=0; //read pointer
+ fEP=0; //end read pointer
+ fFP=0; //end flush pointer
+ fOP=0; //output pointer
+
+ fSeq.fRow=0;
+ fSeq.fLastPad=0;
+ fSeq.fTotalCharge=0;
+ fSeq.fPad=0;
+ fSeq.fTime=0;
+ fSeq.fPad2=0;
+ fSeq.fTime2=0;
+ fSeq.fMean=0;
+ fSeq.fMerge=0;
+
+ for(UInt_t i=0;i<N_mem;i++) fPList[i]=N_clmem;
+ for(UInt_t i=0;i<N_clmem;i++) ClearSeq(i);
+}
+
+void AliL3VHDLClusterFinder::ClearSeq(UShort_t i){
+ fSeqs[i].fRow=0;
+ fSeqs[i].fLastPad=0;
+ fSeqs[i].fTotalCharge=0;
+ fSeqs[i].fPad=0;
+ fSeqs[i].fTime=0;
+ fSeqs[i].fPad2=0;
+ fSeqs[i].fTime2=0;
+ fSeqs[i].fMean=0;
+ fSeqs[i].fMerge=0;
+}
+
+void AliL3VHDLClusterFinder::IncPointer(UShort_t &p, Short_t add, UShort_t N){
+ Short_t pp=p;
+ pp+=add;
+ if(pp>=N) p=UShort_t(pp-N);
+ else if(pp<0) p=UShort_t(pp+N);
+ else p=UShort_t(pp);
+}
+
+inline void AliL3VHDLClusterFinder::IncRPointer(){
+ IncPointer(fRP);
+}
+
+inline void AliL3VHDLClusterFinder::IncWPointer(){
+ IncPointer(fWP);
+
+ if(fWP==fOP){
+ LOG(AliL3Log::kWarning,"AliL3VHDLClusterFinder::IncWPointer","Memory Check")
+ <<"Write pointer overwrites output pointer."<<ENDLOG;
+ }
+}
+
+inline void AliL3VHDLClusterFinder::NextFreeIndex(){
+ IncPointer(fFirst,1,N_clmem);
+ if(fFirst==fLast) {
+ LOG(AliL3Log::kFatal,"AliL3VHDLClusterFinder::GetFreeIndex","Memory Check")
+ <<"No space left in sequence list: "<<fFirst<<"=="<<fLast<<ENDLOG;
+ }
+}
+
+#if 0
+void AliL3ClustFinderNew::WriteClusters(Int_t n_clusters,ClusterData *list)
+{
+ Int_t thisrow,thissector;
+ UInt_t counter = fNClusters;
+
+ for(int j=0; j<n_clusters; j++)
+ {
+ if(!list[j].fFlags) continue; //discard 1 pad clusters
+ if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
+
+ Float_t xyz[3];
+ Float_t fpad =(Float_t)list[j].fPad /(Float_t)list[j].fTotalCharge;
+ Float_t fpad2=fXYErr;
+ Float_t ftime =(Float_t)list[j].fTime /(Float_t)list[j].fTotalCharge;
+ Float_t ftime2=fZErr;
+
+ if(fcalcerr) {
+ fpad2=(Float_t)list[j].fPad2/(Float_t)list[j].fTotalCharge - fpad*fpad;
+ fpad2 = sqrt(fpad2);
+ ftime2=(Float_t)list[j].fTime2/(Float_t)list[j].fTotalCharge - ftime*ftime;
+ ftime2 = sqrt(ftime2);
+ }
+
+ if(fstdout==kTRUE)
+ cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << "+-"<<fpad2<<" time "<<ftime<<"+-"<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
+
+ AliL3Transform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
+ AliL3Transform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
+ if(xyz[0]==0) LOG(AliL3Log::kError,"AliL3ClustFinder","Cluster Finder")
+ <<AliL3Log::kDec<<"Zero cluster"<<ENDLOG;
+ if(fNClusters >= fMaxNClusters)
+ {
+ LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder")
+ <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG;
+ return;
+ }
+ fSpacePointData[counter].fCharge = list[j].fTotalCharge;
+ fSpacePointData[counter].fX = xyz[0];
+ fSpacePointData[counter].fY = xyz[1];
+ fSpacePointData[counter].fZ = xyz[2];
+ fSpacePointData[counter].fPadRow = fCurrentRow;
+ fSpacePointData[counter].fXYErr = fpad2;
+ fSpacePointData[counter].fZErr = ftime2;
+ fSpacePointData[counter].fID = counter
+ +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);//Uli
+
+ fNClusters++;
+ counter++;
+ }
+}
+#endif
+
--- /dev/null
+#ifndef AliL3_VHDLClusterFinder
+#define AliL3_VHDLClusterFinder
+
+#include "AliL3RootTypes.h"
+#include "AliL3Logging.h"
+#include "AliL3AltroMemHandler.h"
+
+//#define DEBUG
+
+struct VHDLClusterData
+{
+ UInt_t fTotalCharge;
+ UInt_t fPad; //mean in pad
+ UInt_t fTime; //mean in time
+ UInt_t fPad2; //for error in XY direction
+ UInt_t fTime2; //for error in Z direction
+ UInt_t fMean; //mean for comparism
+ UInt_t fMerge; //number of merges
+ UShort_t fRow; //row of cluster
+ UShort_t fLastPad; //last pad on merge
+ //UInt_t fChargeFalling; //for deconvolution
+ //UInt_t fLastCharge; //for deconvolution
+};
+typedef struct VHDLClusterData VCData;
+
+//size of ring buffer
+#define N_mem 2500
+//size of cluster list
+#define N_clmem 5000
+
+
+class AliL3VHDLClusterFinder {
+
+ private:
+ AliL3AltroMemHandler fAltromem; //!
+ VCData fSeq; //!
+ VCData fSeqs[N_clmem]; //!
+ UShort_t fPList[N_mem];
+ UShort_t fRow,fNRow;
+ UChar_t fPad,fNPad;
+ UShort_t fRP,fWP,fOP,fEP,fFP; //pointer in ringbuffer
+ UShort_t fLast,fFirst; //free area in memory
+
+ Bool_t fDeconvTime; //not used for now
+ Bool_t fDeconvPad;
+ Bool_t fstdout;
+ Bool_t fcalcerr;
+ Float_t fXYErr;
+ Float_t fZErr;
+
+ Int_t fMatch; //match distance
+ UInt_t fThreshold; //threshold for cluster
+ UInt_t fMinMerge; //minimum number of merges for cluster
+ Int_t fNClusters; //number of found clusters
+
+#ifdef DEBUG
+ FILE *fdeb; //!
+#endif
+
+ void Clear();
+ void ClearSeq(UShort_t i);
+ void FreeSeq(UShort_t i);
+ void IncPointer(UShort_t &p,Short_t add=1,UShort_t N=N_mem);
+ void IncWPointer();
+ void IncRPointer();
+ void NextFreeIndex();
+ void FlushMemory();
+ void PrepareMemory();
+ void OutputMemory();
+ void CompareSeq();
+ void MergeSeq();
+ void InsertSeq();
+ void ProcessSequence();
+ //void WriteClusters(Int_t n_clusters,ClusterData *list);
+
+ public:
+ AliL3VHDLClusterFinder();
+ virtual ~AliL3VHDLClusterFinder();
+
+ void ProcessDigits();
+
+ void SetXYError(Float_t f) {fXYErr=f;}
+ void SetZError(Float_t f) {fZErr=f;}
+ void SetDeconv(Bool_t f) {fDeconvPad=f; fDeconvTime=f;}
+ void SetThreshold(UInt_t i=10) {fThreshold=i;}
+ void SetMatchWidth(UInt_t i=4) {fMatch=i;}
+ void SetMergeMinimum(UInt_t i=1) {fMinMerge=i;}
+ void SetSTDOutput(Bool_t f=kFALSE) {fstdout=f;}
+ void SetCalcErr(Bool_t f=kTRUE) {fcalcerr=f;}
+ void SetASCIIInput(FILE *f){fAltromem.SetASCIIInput(f);}
+
+ Int_t GetNumberOfClusters() {return fNClusters;}
+
+ ClassDef(AliL3VHDLClusterFinder,1)
+
+};
+
+#endif
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff