HLT/src/AliL3ClustFinder.cxx

   1 //*-- Author : Anders Strand Vestbo <mailto:vestbo@fi.uib.no>
   2 // @(#) 16.12.2000
   3 //*-- Copyright &copy ASV
   4
   5 #include <iostream.h>
   6
   7 #include "AliL3Logging.h"
   8 #include "AliL3ClustFinder.h"
   9 #include "AliL3Transform.h"
  10 #include "AliL3DigitData.h"
  11 #include "AliL3SpacePointData.h"
  12
  13 //
  14 // AliL3ClustFinder
  15 //
  16
  17
  18 ClassImp(AliL3ClustFinder)
  19
  20 AliL3ClustFinder::AliL3ClustFinder()
  21 {
  22
  23   fDeconvTime=false;
  24   fDeconvPad=false;
  25   fXYErr = 0;
  26   fZErr = 0;
  27 }
  28
  29 AliL3ClustFinder::AliL3ClustFinder(AliL3Transform *transform)
  30 {
  31
  32   fTransform = transform;
  33
  34   fThreshold = 20;
  35   fDeconvTime = true;
  36   fDeconvPad = true;
  37 }
  38
  39
  40 AliL3ClustFinder::~AliL3ClustFinder()
  41 {
  42   //destructor
  43 }
  44
  45
  46   void AliL3ClustFinder::SetTransformer( AliL3Transform *transform )
  47     {
  48     fTransform = transform;
  49     }
  50
  51
  52 void AliL3ClustFinder::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
  53 {
  54   fNClusters = 0;
  55   fMaxNClusters = nmaxpoints;
  56   fCurrentSlice = slice;
  57   fCurrentPatch = patch;
  58   fFirstRow = firstrow;
  59   fLastRow = lastrow;
  60 }
  61
  62
  63 void AliL3ClustFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
  64 {
  65   fNClusters = 0;
  66   fMaxNClusters = nmaxpoints;
  67   fCurrentSlice = slice;
  68   fCurrentPatch = patch;
  69 }
  70
  71 void AliL3ClustFinder::SetOutputArray(AliL3SpacePointData *pt)
  72 {
  73
  74   fSpacePointData = pt;
  75 }
  76
  77 void AliL3ClustFinder::Read(UInt_t ndigits,AliL3DigitRowData *ptr)
  78 {
  79   fNDigitRowData = ndigits;
  80   fDigitRowData = ptr;
  81
  82 }
  83
  84 void AliL3ClustFinder::ProcessDigits()
  85 {
  86   //Loop over rows, and call processrow
  87
  88
  89   AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fDigitRowData;
  90
  91   for(Int_t i=fFirstRow; i<=fLastRow; i++)
  92     {
  93       fCurrentRow = i;
  94       ProcessRow(tempPt);
  95       Byte_t *tmp = (Byte_t*)tempPt;
  96       Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
  97       tmp += size;
  98       tempPt = (AliL3DigitRowData*)tmp;
  99     }
 100
 101   LOG(AliL3Log::kInformational,"AliL3ClustFinder","Cluster Finder")
 102     <<AliL3Log::kDec<<"ClusterFinder found "<<fNClusters<<" clusters"<<ENDLOG;
 103 }
 104
 105 void AliL3ClustFinder::ProcessRow(AliL3DigitRowData* tempPt)
 106 {
 107
 108   Int_t cl_found,trackID[3];
 109   //UInt_t rows;
 110
 111   //Reserve array for found clusters:
 112   resx aliresx[5000];
 113   bzero((char *)aliresx,5000*sizeof(resx));
 114
 115   //local results
 116   resx *r;
 117 //  UInt_t pres1[2500],pres2[2500];
 118   resx *pres1[2500];
 119   resx *pres2[2500];
 120   resx rr_local;
 121   resx *rr = &rr_local;
 122
 123   //Set variables for each row:
 124   Int_t pres_cou1,pres_cou2;
 125 //  UInt_t *r1,*r2;
 126   resx **r1,**r2;
 127   pres_cou1=pres_cou2=0;
 128   r1=pres2;
 129   r2=pres1;
 130   r=aliresx;
 131
 132   //variables for each pad:
 133   UInt_t start;
 134   Int_t start_new=-1;
 135   Int_t go_on;
 136   Int_t cl_counter=0,rows;
 137   Int_t prev_start;
 138
 139   UInt_t last_pad = 123456789;
 140
 141   //Loop over digits in this row:
 142   for(start=0; start<tempPt->fNDigit; start++)
 143     {
 144       AliL3DigitData *bins = tempPt->fDigitData;
 145
 146       if(bins[start].fPad!=last_pad)
 147         {//new pad
 148
 149           last_pad=bins[start].fPad;
 150
 151           if(r2==pres2)
 152             {
 153               r1=pres2;
 154               r2=pres1;
 155             }
 156           else
 157             {
 158               r1=pres1;
 159               r2=pres2;
 160             }
 161           pres_cou1=pres_cou2;
 162           pres_cou2=0;
 163
 164           start_new=-1;
 165
 166           cl_counter=0;
 167           prev_start=0;
 168
 169
 170         }//Looping over pads
 171
 172       UInt_t av,charge;
 173       UInt_t mean;
 174       UInt_t flags;
 175       Int_t last_falling=0;
 176 //      UInt_t *ri;
 177       resx **ri;
 178       UInt_t tmp_charge;
 179
 180       //Set flag:
 181       flags=0;
 182
 183       cl_counter++;
 184       if(cl_counter>MAX_C) {LOG(AliL3Log::kWarning,"AliL3ClustFinder::ProcessRow","Cluster Finder")
 185       <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG; return;} //too many clusters
 186
 187
 188       if(fDeconvTime)
 189         {
 190
 191           //this is a goto label:
 192         redo: ;
 193         //divide sequenz in two in case there is a change in slope:
 194         if(start_new>-1)
 195           {
 196             start=start_new;
 197             start_new=-1;
 198           }
 199         }
 200       //average and charge
 201       av=charge=0;
 202
 203       //ok, create a new cluster:
 204       go_on=1;
 205
 206       if(fDeconvTime)
 207         {
 208
 209           last_falling=flags=0;
 210         }
 211
 212       if(fDeconvPad)
 213         {
 214           flags=0;
 215         }
 216
 217
 218       //block to introduce new variables:
 219       {
 220         UInt_t last_a;
 221         last_a=0;
 222
 223         //Loop over this sequenz:
 224         while(1)
 225           {
 226             UInt_t aa;
 227             //Int_t start_temp=start;
 228
 229             //get the adc-value:
 230             aa=bins[start].fCharge;
 231
 232
 233             if(fDeconvTime)
 234               {
 235                 //Check if the last pixel in this sequenz is bigger or smaller than this:
 236                 if(aa>last_a)
 237                   {
 238                     if(last_falling)
 239                       {
 240                         start_new=start;
 241                         break;
 242                       }
 243                   }
 244                 else last_falling=1;
 245                 last_a=aa;
 246
 247               }
 248             //sum of total charge of this cluster on this pad
 249             charge +=aa;
 250             //this one is needed to determine the mean over all pads:
 251             //av += start * (aa);
 252             av+=bins[start].fTime*(aa);
 253
 254             if((bins[start+1].fTime)!=(bins[start].fTime+1) ||
 255                ((bins[start+1].fPad)!=(bins[start].fPad))) break;
 256
 257             start++;
 258           }//Loop over sequenz
 259       }//block of new variables
 260
 261       if(fDeconvTime)
 262         {
 263           if(start_new>0) flags = FLAG_DOUBLE_T;
 264         }
 265       //calculate mean in time direction for this sequenz:
 266       if(charge)
 267         {
 268           //mean=cog for this sequenz:
 269           mean=av/charge;
 270         }
 271       else
 272         {
 273           charge=1;
 274           mean=1;
 275         }
 276       //get the pointer to results on previous pad(s)
 277       ri=r1;
 278
 279       //Center of gravity in pad direction
 280       //tmp_charge=charge*bins[start-1].fPad;
 281       tmp_charge=charge*bins[start].fPad;
 282
 283       //compare with prevously stored results (previous pads)
 284       for(int k=0; k<pres_cou1; k++)
 285         {
 286           //variable to store difference:
 287           Int_t v;
 288
 289           //structure that contains means (+more) of last pads(s)
 290           resx *rr_tmp;
 291
 292           //go get it
 293           rr_tmp = (resx *) *ri;
 294
 295           //difference between mean and mean on previous pads
 296           v=mean-rr_tmp->mean;
 297
 298
 299           if(v<-PARAM1) break;
 300
 301           //goto next mean on previous pad
 302           ri++;
 303
 304           //is there a matching cluster on previous pad:
 305           if(v<=PARAM1)
 306             {
 307               //ok, we got a new one, so we will add the new and the old
 308               rr=&rr_local;
 309
 310               //Add this function somewhere
 311               memcpy7((UInt_t *)rr,(UInt_t *)rr_tmp);
 312
 313               /* in case we found another sequenz-section on this pad
 314                  don't loop over all sequenzes on last pad but start
 315                  where we found a matching candidate for this sequenz-
 316                  section, adjust upper limit (=pres_cou1)*/
 317               r1=ri;
 318               pres_cou1-=1+k;
 319
 320               if(fDeconvPad)
 321                 {
 322                   if(charge>rr->scharge)
 323                     {
 324                       if(rr->falling)
 325                         {
 326                           //previous state was falling
 327                           flags |= FLAG_DOUBLE_PAD;
 328                           rr_tmp->flags |= flags;
 329                           //and create a new one
 330                           break;
 331                         }
 332                     }
 333                   else
 334                     {
 335                       rr->falling = 1;
 336                     }
 337
 338                 }
 339               //don't create a new one
 340               go_on=0;
 341
 342               //store the pointer to the matched in "2"
 343 //            r2[pres_cou2++] = (UInt_t) rr_tmp;
 344               r2[pres_cou2++] =  rr_tmp;
 345
 346               //calculate and fill the new means, charge etc. in rr_tmp
 347               //charge on this pad to determine whether we cut in pad direction
 348               rr->scharge = charge;
 349               //unset FLAG_ONEPAD since we have at least 2 pads
 350               rr->flags &= (~FLAG_ONEPAD);
 351               //summ the charges for de/dx
 352               rr->charge += charge;
 353               //calculate mean in pad direction
 354               rr->pad += tmp_charge;
 355               //calculate mean in time direction
 356               rr->t+=av;
 357               //store the mean of this sequenz on this pad
 358               rr->mean = mean;
 359
 360               //Copy it back to storage
 361               memcpy7((UInt_t *)rr_tmp,(UInt_t *)rr);
 362
 363               //we found a matching one so stop looping over pads on previous pad
 364               break;
 365
 366             }//matching cluster on previous pads
 367         }//loop: means on previous pad
 368
 369       //here we create from scratch new clusters
 370       if(go_on)
 371         {
 372           //store this one beacuse it is the first
 373 //        r2[pres_cou2++]= (UInt_t )r;
 374           r2[pres_cou2++]= r;
 375
 376           mstore((UInt_t *)r,av,tmp_charge,charge,flags|FLAG_ONEPAD,mean,trackID);
 377           r++;
 378         }
 379
 380       if(fDeconvTime)
 381         {
 382
 383           //maybe we still have to do something
 384           if(start_new>=0) goto redo;
 385         }
 386
 387     }//end of loop over digits
 388
 389
 390   //number of clusters found
 391   cl_found = r-aliresx;
 392
 393   //tot_cl_found+=cl_found;
 394
 395   //there was something on this padrow
 396   rows++;
 397
 398   WriteClusters(cl_found,aliresx);
 399
 400 }
 401
 402 void AliL3ClustFinder::memcpy7(UInt_t *dst, UInt_t *src)
 403 {
 404   Int_t i;
 405   for(i=0; i<7; i++)
 406     {
 407       *dst++ = *src++;
 408     }
 409   return;
 410 }
 411
 412 void AliL3ClustFinder::mstore(UInt_t *r,UInt_t av,UInt_t pad,UInt_t ch,UInt_t flags,UInt_t mean,Int_t *trackID)
 413 {
 414   resx *rr = (resx *) r;
 415
 416   rr->mean = mean;
 417   rr->charge = ch;
 418   rr->flags =flags;
 419   rr->scharge = ch;
 420   rr->pad = pad;
 421   rr->t = av;
 422   rr->falling=0;
 423
 424   rr->trackID[0]=trackID[0];
 425   rr->trackID[1]=trackID[1];
 426   rr->trackID[2]=trackID[2];
 427   return;
 428 }
 429
 430 void AliL3ClustFinder::WriteClusters(Int_t ncl,resx *r)
 431 {
 432 /*
 433    if(fNClusters >= fMaxNClusters)
 434    {
 435    LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder")
 436    <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG;
 437    return;
 438    }
 439 */
 440   Int_t thisrow,thissector;
 441   UInt_t counter = fNClusters;
 442
 443   if(!fXYErr || !fZErr)
 444     {
 445       LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Errors")
 446         <<"No Errors"<<ENDLOG;
 447       return;
 448     }
 449
 450   for(int j=0; j<ncl; j++)
 451     {
 452       if(r[j].flags & FLAG_ONEPAD) continue; //discard 1 pad clusters
 453       if(r[j].charge < fThreshold) continue; //noise cluster
 454       Float_t xyz[3];
 455
 456       Float_t fpad=(Float_t)r[j].pad/(Float_t)r[j].charge;
 457       Float_t ftime=(Float_t)r[j].t/(Float_t)r[j].charge;
 458
 459       fTransform->Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
 460       //if(fCurrentRow > 54) {thisrow = fCurrentRow-55; thissector = fCurrentSlice+36;}
 461       //else {thisrow = fCurrentRow; thissector = fCurrentSlice;}
 462       fTransform->Raw2Local(xyz,thissector,thisrow,fpad,ftime);
 463       if(xyz[0]==0) LOG(AliL3Log::kError,"AliL3ClustFinder","Cluster Finder")
 464         <<AliL3Log::kDec<<"Zero cluster"<<ENDLOG;
 465       if(fNClusters >= fMaxNClusters)
 466           {
 467           LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder")
 468               <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG;
 469           return;
 470           }
 471       fSpacePointData[counter].fX = xyz[0];
 472       fSpacePointData[counter].fY = xyz[1];
 473       fSpacePointData[counter].fZ = xyz[2];
 474       fSpacePointData[counter].fPadRow = fCurrentRow;
 475       fSpacePointData[counter].fXYErr = fXYErr;
 476       fSpacePointData[counter].fZErr = fZErr;
 477       fSpacePointData[counter].fID = counter
 478                   +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);
 479
 480       fNClusters++;
 481       counter++;
 482
 483     }
 484
 485
 486 }