Created cluster finder class that simulates the VHDL cluster finder on Altro data.
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 14 Jun 2002 21:40:18 +0000 (21:40 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 14 Jun 2002 21:40:18 +0000 (21:40 +0000)
HLT/misc/AliL3MiscLinkDef.h
HLT/misc/AliL3VHDLClusterFinder.cxx [new file with mode: 0644]
HLT/misc/AliL3VHDLClusterFinder.h [new file with mode: 0644]
HLT/misc/Makefile

index ec6d87d..f28b4e2 100644 (file)
@@ -9,5 +9,6 @@
 #pragma link C++ class AliTransBit_v2; 
 #pragma link C++ class AliL3AltroMemHandler;
 #pragma link C++ class AliL3DataHandler;
+#pragma link C++ class AliL3VHDLClusterFinder;
 #endif
 
diff --git a/HLT/misc/AliL3VHDLClusterFinder.cxx b/HLT/misc/AliL3VHDLClusterFinder.cxx
new file mode 100644 (file)
index 0000000..f451903
--- /dev/null
@@ -0,0 +1,404 @@
+/* $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
+
diff --git a/HLT/misc/AliL3VHDLClusterFinder.h b/HLT/misc/AliL3VHDLClusterFinder.h
new file mode 100644 (file)
index 0000000..e5b7539
--- /dev/null
@@ -0,0 +1,98 @@
+#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
index 4e31ed4..8d07a72 100644 (file)
@@ -52,7 +52,7 @@ else
 INCLUDES += -I/prog/alice/level3/kip/MLUC/include
 endif
 
-SRCS   = AliTransBit.cxx AliL3AltroMemHandler.cxx AliL3DataHandler.cxx
+SRCS   = AliTransBit.cxx AliL3AltroMemHandler.cxx AliL3DataHandler.cxx AliL3VHDLClusterFinder.cxx
 
 DICT = AliL3MiscCint.cxx
 DICTH = AliL3MiscCint.h
@@ -101,13 +101,6 @@ so :
 
 
 
-
-
-
-
-
-
-