HLT/misc/AliL3VHDLClusterFinder.cxx

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