HLT/misc/AliL3VHDLClusterFinder.cxx

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