HLT/misc/AliL3VHDLClusterFinder.cxx

   1 /* $Id$
   2 // Author: Constantin Loizides <mailto:loizides@ikf.physik.uni-frankfurt.de>
   3 */
   4
   5 #include <math.h>
   6 #include <time.h>
   7 #include <iostream.h>
   8 #include <stdio.h>
   9 #include <string.h>
  10
  11 #include "AliL3VHDLClusterFinder.h"
  12
  13 /** \class AliL3VHDLClusterFinder
  14 //<pre>
  15 //____________________________________________________
  16 // AliL3VHDLClusterFinder
  17 //
  18 // The current VHDL cluster finder for HLT
  19 // Based on STAR L3
  20 //
  21 // Most important parameters:
  22 // fThreshold - threshold for noise clusters
  23 // fMatch - length in time for overlapping sequences
  24 //</pre>
  25 */
  26
  27 ClassImp(AliL3VHDLClusterFinder)
  28
  29 AliL3VHDLClusterFinder::AliL3VHDLClusterFinder()
  30 {
  31   fMatch = 4;
  32   fThreshold = 10;
  33   fMinMerge = 1;
  34   fNClusters=0;
  35
  36   fXYErr = 0.2;
  37   fZErr = 0.3;
  38   fDeconvPad = kTRUE;
  39   fDeconvTime = kTRUE;
  40   fstdout = kFALSE;
  41   fcalcerr = kTRUE;
  42
  43   Clear();
  44 #ifdef DEBUG
  45   fdeb=fopen("vhdlclusterfinder.debug","w");
  46   //fdeb=stderr;
  47 #endif
  48 }
  49
  50 AliL3VHDLClusterFinder::~AliL3VHDLClusterFinder()
  51 {
  52 #ifdef DEBUG
  53   fclose(fdeb);
  54 #endif
  55 }
  56
  57 void AliL3VHDLClusterFinder::ProcessDigits()
  58 {
  59   //Loop over data like the analyzer of the VHDL code
  60   const UChar_t n=255;
  61   UShort_t rrow=0,rtime=0;
  62   UChar_t rpad=0,i=n;
  63   UShort_t *charges=new UShort_t[n];
  64   Int_t tc=0,mp=0,mt=0,sp=0,st=0;
  65
  66   fNClusters=0;
  67   fRow=0;
  68   fPad=0;
  69   Clear();
  70
  71   //loop over input data
  72   while(fAltromem.ReadSequence(rrow,rpad,rtime,i,&charges)){
  73     tc=0;mp=0;mt=0;sp=0;st=0;
  74 #if 0
  75     cout << "Padrow " << (int)rrow << " pad " << (int)rpad << " time " <<(int)rtime << " charges ";
  76     for(UChar_t ii=0;ii<i;ii++) cout << (int)charges[ii] << " ";
  77     cout << endl;
  78 #endif
  79 #ifdef DEBUG
  80     fprintf(fdeb,"ProcessDigits: Input Data: %d %d %d charges:",(int)rrow,(int)rpad,(int)rtime);
  81     for(UChar_t ii=0;ii<i;ii++) fprintf(fdeb," %d",(int)charges[ii]);
  82     fprintf(fdeb,"\n");
  83 #endif
  84
  85     fNRow=rrow;
  86     fNPad=rpad;
  87
  88     //calculate sequence values
  89     //no deconvulution so far
  90     for(UChar_t ii=0;ii<i;ii++){
  91       tc+=charges[ii];
  92       mt+=(rtime-ii)*charges[ii];
  93       st+=(rtime-ii)*(rtime-ii)*charges[ii];
  94     }
  95     mp=rpad*tc;
  96     sp=rpad*rpad*tc;
  97
  98     fSeq.fTotalCharge=tc;
  99     fSeq.fPad=mp;
 100     fSeq.fTime=mt;
 101     fSeq.fPad2=sp;
 102     fSeq.fTime2=st;
 103     fSeq.fMean=0;
 104     //if(tc!=0) fSeq.fMean=mt/tc;
 105     if(tc!=0) fSeq.fMean=rtime-i/2;
 106     fSeq.fMerge=0;
 107     fSeq.fRow=rrow;
 108     fSeq.fLastPad=rpad;
 109
 110     //work on this sequence
 111     ProcessSequence();
 112     //output one cluster
 113     OutputMemory();
 114
 115 #ifdef DEBUG
 116     fflush(fdeb);
 117 #endif
 118     i=n; //store size of charges array
 119   } //loop over data
 120
 121   //flush everything left
 122   FlushMemory();
 123   while(fOP!=fFP) OutputMemory();
 124 }
 125
 126 void AliL3VHDLClusterFinder::ProcessSequence()
 127 {
 128   if(fNRow!=fRow) FlushMemory();
 129   else if(fNPad==fPad+1) PrepareMemory();
 130   else if(fNPad!=fPad) FlushMemory();
 131
 132   //store new row and pad values
 133   fRow=fNRow;
 134   fPad=fNPad;
 135
 136 #ifdef DEBUG
 137   fprintf(fdeb,"ProcessSequence: Mean=%d TC=%d ",fSeq.fMean,fSeq.fTotalCharge);
 138 #endif
 139
 140   CompareSeq(); //merge or insert
 141 }
 142
 143 void AliL3VHDLClusterFinder::PrepareMemory()
 144 {
 145 #ifdef DEBUG
 146   fprintf(fdeb,"PrepareMemory %d %d %d\n",fRP,fEP,fWP);
 147 #endif
 148   fRP=fEP;
 149   fEP=fWP;
 150 }
 151
 152 void AliL3VHDLClusterFinder::FlushMemory()
 153 {
 154 #ifdef DEBUG
 155   fprintf(fdeb,"FlushMemory %d %d %d %d\n",fFP,fRP,fEP,fWP);
 156 #endif
 157   fFP=fWP;
 158   fRP=fWP;
 159   fEP=fWP;
 160 }
 161
 162 void AliL3VHDLClusterFinder::CompareSeq()
 163 {
 164
 165   while(fRP!=fEP){
 166     Int_t diff=fSeqs[fPList[fRP]].fMean-fSeq.fMean;
 167
 168     if(diff>fMatch){ //no match
 169       IncRPointer(); //cluster finished
 170       continue;
 171     } else if(diff<-fMatch){ //no match
 172       break;                 //insert new cluster
 173     }
 174     else { //match found, merge it
 175       MergeSeq();
 176       return;
 177     }
 178   }
 179
 180   InsertSeq(); //start new cluster
 181 }
 182
 183 void AliL3VHDLClusterFinder::MergeSeq()
 184 {
 185 #ifdef DEBUG
 186   fprintf(fdeb,"merged with Mean=%d TC=%d (new Merge=%d) %d %d\n",fSeqs[fPList[fRP]].fMean,fSeqs[fPList[fRP]].fTotalCharge,fSeqs[fPList[fRP]].fMerge+1,fRow,fPad);
 187 #endif
 188   if(fSeqs[fPList[fRP]].fRow==fSeq.fRow){
 189     LOG(AliL3Log::kWarning,"AliL3VHDLClusterFinder::","Memory Check")
 190       <<"Sequences can be merged on the same rows only."<<ENDLOG;
 191   }
 192   if(fSeqs[fPList[fRP]].fLastPad+1!=fSeq.fLastPad){
 193     LOG(AliL3Log::kWarning,"AliL3VHDLClusterFinder::","Memory Check")
 194       <<"Sequences can be merged on consecutive pads only."<<ENDLOG;
 195   }
 196
 197   fSeqs[fPList[fRP]].fMean=fSeq.fMean; //take the new mean
 198   fSeqs[fPList[fRP]].fLastPad=fSeq.fLastPad;
 199   fSeqs[fPList[fRP]].fTotalCharge+=fSeq.fTotalCharge;
 200   fSeqs[fPList[fRP]].fPad+=fSeq.fPad;
 201   fSeqs[fPList[fRP]].fTime+=fSeq.fTime;
 202   fSeqs[fPList[fRP]].fPad2+=fSeq.fPad2;
 203   fSeqs[fPList[fRP]].fTime2+=fSeq.fTime2;
 204   fSeqs[fPList[fRP]].fMerge+=1;
 205
 206   fPList[fWP]=fPList[fRP];
 207   fPList[fRP]=N_clmem; //mark it to be free
 208
 209   IncRPointer();
 210   IncWPointer();
 211 }
 212
 213 void AliL3VHDLClusterFinder::InsertSeq()
 214 {
 215 #ifdef DEBUG
 216   fprintf(fdeb,"inserted %d %d\n",fRow,fPad);
 217 #endif
 218   NextFreeIndex();    //get next index
 219   fSeqs[fFirst]=fSeq; //store data
 220   fPList[fWP]=fFirst; //store index
 221   IncWPointer();
 222 }
 223
 224 void AliL3VHDLClusterFinder::OutputMemory()
 225 {
 226   Float_t mtime=0,mpad=0;
 227   Float_t mtime2=0,mpad2=0;
 228   UInt_t tc,row,mno;
 229
 230   if(fOP!=fFP){
 231   //while(fOP!=fFP){
 232
 233     UInt_t index=fPList[fOP];
 234     fPList[fOP]=N_clmem;
 235     //cout << fOP << " - " << fFP << " - " << index << endl;
 236     IncPointer(fOP);
 237     if(index>=N_clmem) return; //nothing to do
 238     //if(index>=N_clmem) continue; //nothing to do
 239
 240     tc=fSeqs[index].fTotalCharge;
 241     mno=fSeqs[index].fMerge;
 242     row=fSeqs[index].fRow;
 243     if(tc!=0){
 244       mtime=(Float_t)(fSeqs[index].fTime)/tc;
 245       mtime2=sqrt((Float_t)(fSeqs[index].fTime2)/tc-mtime*mtime);
 246     }
 247     if(tc!=0){
 248       mpad=(Float_t)(fSeqs[index].fPad)/tc;
 249       mpad2=sqrt((Float_t)(fSeqs[index].fPad2)/tc-mpad*mpad);
 250     }
 251
 252     if(mno<fMinMerge){ //noise cluster
 253 #ifdef DEBUG
 254     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);
 255 #endif
 256       FreeSeq(index);
 257       return;
 258       //continue;
 259     }
 260
 261     if(tc<fThreshold){ //total charge below threshold
 262 #ifdef DEBUG
 263     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);
 264 #endif
 265       FreeSeq(index);
 266       return;
 267       //continue;
 268     }
 269
 270 #ifdef DEBUG
 271     fprintf(fdeb,"OutputMemory %d: Row %d Pad %.3f Time %.3f TC %d Merge %d\n",fNClusters,row,mpad,mtime,tc,mno);
 272 #endif
 273
 274     if(stdout)
 275       cout<<"WriteCluster: padrow "<<row<<" pad "<<mpad<< " +- "<<mpad2<<" time "<<mtime<<" +- "<<mtime2<<" charge "<<tc<<endl;
 276
 277     fNClusters++;
 278     FreeSeq(index);
 279   }
 280 }
 281
 282 inline void AliL3VHDLClusterFinder::FreeSeq(UShort_t i)
 283 {
 284   ClearSeq(i);
 285   IncPointer(fLast,1,N_clmem);
 286 }
 287
 288 void AliL3VHDLClusterFinder::Clear()
 289 {
 290   fFirst=0; //first list pointer
 291   fLast=0;  //last  list pointer
 292
 293   fWP=0; //write pointer
 294   fRP=0; //read pointer
 295   fEP=0; //end read pointer
 296   fFP=0; //end flush pointer
 297   fOP=0; //output pointer
 298
 299   fSeq.fRow=0;
 300   fSeq.fLastPad=0;
 301   fSeq.fTotalCharge=0;
 302   fSeq.fPad=0;
 303   fSeq.fTime=0;
 304   fSeq.fPad2=0;
 305   fSeq.fTime2=0;
 306   fSeq.fMean=0;
 307   fSeq.fMerge=0;
 308
 309   for(UInt_t i=0;i<N_mem;i++) fPList[i]=N_clmem;
 310   for(UInt_t i=0;i<N_clmem;i++) ClearSeq(i);
 311 }
 312
 313 void AliL3VHDLClusterFinder::ClearSeq(UShort_t i){
 314   fSeqs[i].fRow=0;
 315   fSeqs[i].fLastPad=0;
 316   fSeqs[i].fTotalCharge=0;
 317   fSeqs[i].fPad=0;
 318   fSeqs[i].fTime=0;
 319   fSeqs[i].fPad2=0;
 320   fSeqs[i].fTime2=0;
 321   fSeqs[i].fMean=0;
 322   fSeqs[i].fMerge=0;
 323 }
 324
 325 void AliL3VHDLClusterFinder::IncPointer(UShort_t &p, Short_t add, UShort_t N){
 326   Short_t pp=p;
 327   pp+=add;
 328   if(pp>=N) p=UShort_t(pp-N);
 329   else if(pp<0) p=UShort_t(pp+N);
 330   else p=UShort_t(pp);
 331 }
 332
 333 inline void AliL3VHDLClusterFinder::IncRPointer(){
 334   IncPointer(fRP);
 335 }
 336
 337 inline void AliL3VHDLClusterFinder::IncWPointer(){
 338   IncPointer(fWP);
 339
 340   if(fWP==fOP){
 341     LOG(AliL3Log::kWarning,"AliL3VHDLClusterFinder::IncWPointer","Memory Check")
 342       <<"Write pointer overwrites output pointer."<<ENDLOG;
 343   }
 344 }
 345
 346 inline void AliL3VHDLClusterFinder::NextFreeIndex(){
 347   IncPointer(fFirst,1,N_clmem);
 348   if(fFirst==fLast) {
 349     LOG(AliL3Log::kFatal,"AliL3VHDLClusterFinder::GetFreeIndex","Memory Check")
 350       <<"No space left in sequence list: "<<fFirst<<"=="<<fLast<<ENDLOG;
 351   }
 352 }
 353
 354 #if 0
 355 void AliL3ClustFinderNew::WriteClusters(Int_t n_clusters,ClusterData *list)
 356 {
 357   Int_t thisrow,thissector;
 358   UInt_t counter = fNClusters;
 359
 360   for(int j=0; j<n_clusters; j++)
 361     {
 362       if(!list[j].fFlags) continue; //discard 1 pad clusters
 363       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
 364
 365       Float_t xyz[3];
 366       Float_t fpad =(Float_t)list[j].fPad /(Float_t)list[j].fTotalCharge;
 367       Float_t fpad2=fXYErr;
 368       Float_t ftime =(Float_t)list[j].fTime /(Float_t)list[j].fTotalCharge;
 369       Float_t ftime2=fZErr;
 370
 371       if(fcalcerr) {
 372         fpad2=(Float_t)list[j].fPad2/(Float_t)list[j].fTotalCharge - fpad*fpad;
 373         fpad2 = sqrt(fpad2);
 374         ftime2=(Float_t)list[j].fTime2/(Float_t)list[j].fTotalCharge - ftime*ftime;
 375         ftime2 = sqrt(ftime2);
 376       }
 377
 378       if(fstdout==kTRUE)
 379         cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << "+-"<<fpad2<<" time "<<ftime<<"+-"<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
 380
 381       AliL3Transform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
 382       AliL3Transform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
 383       if(xyz[0]==0) LOG(AliL3Log::kError,"AliL3ClustFinder","Cluster Finder")
 384         <<AliL3Log::kDec<<"Zero cluster"<<ENDLOG;
 385       if(fNClusters >= fMaxNClusters)
 386         {
 387           LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder")
 388             <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG;
 389           return;
 390         }
 391       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
 392       fSpacePointData[counter].fX = xyz[0];
 393       fSpacePointData[counter].fY = xyz[1];
 394       fSpacePointData[counter].fZ = xyz[2];
 395       fSpacePointData[counter].fPadRow = fCurrentRow;
 396       fSpacePointData[counter].fXYErr = fpad2;
 397       fSpacePointData[counter].fZErr  = ftime2;
 398       fSpacePointData[counter].fID = counter
 399         +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);//Uli
 400
 401       fNClusters++;
 402       counter++;
 403     }
 404 }
 405 #endif
 406