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