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