HLT/src/AliL3ClustFinderNew.cxx

   1 //$Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ASV
   5
   6 #include <iostream.h>
   7 #include "AliL3Logging.h"
   8 #include "AliL3ClustFinderNew.h"
   9 #include "AliL3DigitData.h"
  10 #include "AliL3Transform.h"
  11 #include "AliL3SpacePointData.h"
  12
  13 //_____________________________________________________________
  14 // AliL3ClustFinderNew
  15 //
  16 // The current cluster finder for HLT
  17 // Based on STAR L3
  18
  19 ClassImp(AliL3ClustFinderNew)
  20
  21 AliL3ClustFinderNew::AliL3ClustFinderNew()
  22 {
  23   fMatch = 4;
  24   fThreshold =10;
  25   fDeconvPad = kTRUE;
  26   fDeconvTime = kTRUE;
  27 }
  28
  29 AliL3ClustFinderNew::AliL3ClustFinderNew(AliL3Transform *transform)
  30 {
  31   fTransform = transform;
  32   fMatch = 4;
  33   fThreshold =10;
  34   fDeconvTime = kTRUE;
  35   fDeconvPad = kTRUE;
  36 }
  37
  38 AliL3ClustFinderNew::~AliL3ClustFinderNew()
  39 {
  40
  41
  42 }
  43
  44 void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
  45 {
  46   fNClusters = 0;
  47   fMaxNClusters = nmaxpoints;
  48   fCurrentSlice = slice;
  49   fCurrentPatch = patch;
  50   fFirstRow = firstrow;
  51   fLastRow = lastrow;
  52   fDeconvTime = kTRUE;
  53   fDeconvPad = kTRUE;
  54 }
  55
  56 void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
  57 {
  58   fNClusters = 0;
  59   fMaxNClusters = nmaxpoints;
  60   fCurrentSlice = slice;
  61   fCurrentPatch = patch;
  62 }
  63
  64 void AliL3ClustFinderNew::SetOutputArray(AliL3SpacePointData *pt)
  65 {
  66
  67   fSpacePointData = pt;
  68 }
  69
  70
  71 void AliL3ClustFinderNew::Read(UInt_t ndigits,AliL3DigitRowData *ptr)
  72 {
  73   fNDigitRowData = ndigits;
  74   fDigitRowData = ptr;
  75
  76 }
  77
  78 void AliL3ClustFinderNew::ProcessDigits()
  79 {
  80   //Loop over rows, and call processrow
  81
  82
  83   AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fDigitRowData;
  84
  85   for(Int_t i=fFirstRow; i<=fLastRow; i++)
  86     {
  87       fCurrentRow = i;
  88       ProcessRow(tempPt);
  89       Byte_t *tmp = (Byte_t*)tempPt;
  90       Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
  91       tmp += size;
  92       tempPt = (AliL3DigitRowData*)tmp;
  93     }
  94 }
  95
  96
  97
  98 void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
  99 {
 100
 101   UInt_t last_pad = 123456789;
 102
 103   ClusterData *pad1[2500]; //2 lists for internal memory=2pads
 104   ClusterData *pad2[2500]; //2 lists for internal memory=2pads
 105   ClusterData clusterlist[5000]; //Clusterlist
 106
 107   ClusterData **currentPt; //List of pointers to the current pad
 108   ClusterData **previousPt; //List of pointers to the previous pad
 109   currentPt = pad2;
 110   previousPt = pad1;
 111   UInt_t n_previous=0,n_current=0,n_total=0;
 112
 113   //Loop over sequences of this row:
 114   for(UInt_t bin=0; bin<tempPt->fNDigit; bin++)
 115     {
 116
 117       AliL3DigitData *data = tempPt->fDigitData;
 118       if(data[bin].fPad != last_pad)
 119         {
 120           //This is a new pad
 121
 122           //Switch:
 123           if(currentPt == pad2)
 124             {
 125               currentPt = pad1;
 126               previousPt = pad2;
 127
 128             }
 129           else
 130             {
 131               currentPt = pad2;
 132               previousPt = pad1;
 133             }
 134           n_previous = n_current;
 135           n_current = 0;
 136           if(bin[data].fPad != last_pad+1)
 137             {
 138               //this happens if there is a pad with no signal.
 139               n_previous = n_current = 0;
 140             }
 141           last_pad = data[bin].fPad;
 142         }
 143
 144       Bool_t new_cluster = kTRUE;
 145
 146       UInt_t seq_charge=0,seq_average=0;
 147
 148       UInt_t last_charge=0,last_was_falling=0;
 149       Int_t new_bin=-1;
 150       if(fDeconvTime)
 151         {
 152         redo: //This is a goto.
 153           if(new_bin > -1)
 154             {
 155               bin = new_bin;
 156               new_bin = -1;
 157             }
 158
 159           last_charge=0;
 160           last_was_falling = 0;
 161         }
 162
 163
 164       while(1) //Loop over current sequence
 165         {
 166           if(data[bin].fTime >= fTransform->GetNTimeBins())
 167             {
 168               LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Digits")
 169                 <<"Timebin out of range "<<(Int_t)data[bin].fTime<<ENDLOG;
 170               break;
 171             }
 172
 173           //Get the current ADC-value
 174           UInt_t charge = data[bin].fCharge;
 175
 176           if(fDeconvTime)
 177             {
 178               //Check if the last pixel in the sequence is smaller than this
 179               if(charge > last_charge)
 180                 {
 181                   if(last_was_falling)
 182                     {
 183                       new_bin = bin;
 184                       break;
 185                     }
 186                 }
 187               else last_was_falling = 1; //last pixel was larger than this
 188               last_charge = charge;
 189             }
 190
 191           //Sum the total charge of this sequence
 192           seq_charge += charge;
 193           seq_average += data[bin].fTime*charge;
 194
 195           //Check where to stop:
 196           if(data[bin+1].fPad != data[bin].fPad) //new pad
 197             break;
 198           if(data[bin+1].fTime != data[bin].fTime+1) //end of sequence
 199             break;
 200
 201           bin++;
 202         }//end loop over sequence
 203
 204       //Calculate mean of sequence:
 205       Int_t seq_mean=0;
 206       if(seq_charge)
 207         seq_mean = seq_average/seq_charge;
 208       else
 209         {
 210           LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Data")
 211             <<"Error in data given to the cluster finder"<<ENDLOG;
 212           seq_mean = 1;
 213           seq_charge = 1;
 214         }
 215
 216       //Calculate mean in pad direction:
 217       Int_t pad_mean = seq_charge*data[bin].fPad;
 218
 219       //Compare with results on previous pad:
 220       for(UInt_t p=0; p<n_previous; p++)
 221         {
 222           Int_t difference = seq_mean - previousPt[p]->fMean;
 223           if(difference < -fMatch)
 224             break;
 225
 226           if(difference <= fMatch)//There is a match here!!
 227             {
 228
 229               ClusterData *local = previousPt[p];
 230               if(fDeconvPad)
 231                 {
 232                   if(seq_charge > local->fLastCharge)
 233                     {
 234                       if(local->fChargeFalling)//The previous pad was falling
 235                         {
 236                           break;//create a new cluster
 237                         }
 238                     }
 239                   else
 240                     local->fChargeFalling = 1;
 241                   local->fLastCharge = seq_charge;
 242                 }
 243
 244               //Don't create a new cluster, because we found a match
 245               new_cluster = kFALSE;
 246
 247               //Update cluster on current pad with the matching one:
 248               //currentPt[n_current] = previousPt[p];
 249
 250               local->fTotalCharge += seq_charge;
 251               local->fPad += pad_mean;
 252               local->fTime += seq_average;
 253               local->fMean = seq_mean;
 254               local->fFlags++; //means we have more than one pad
 255
 256               currentPt[n_current] = local;
 257               n_current++;
 258
 259               break;
 260             }//Checking for match at previous pad
 261         }//Loop over results on previous pad.
 262
 263       if(new_cluster)
 264         {
 265           //Start a new cluster. Add it to the clusterlist, and update
 266           //the list of pointers to clusters in current pad.
 267           //current pad will be previous pad on next pad.
 268
 269           //Add to the clusterlist:
 270           ClusterData *tmp = &clusterlist[n_total];
 271           tmp->fTotalCharge = seq_charge;
 272           tmp->fPad = pad_mean;
 273           tmp->fTime = seq_average;
 274           tmp->fMean = seq_mean;
 275           tmp->fFlags = 0;  //flags for 1 pad clusters
 276           if(fDeconvPad)
 277             {
 278               tmp->fChargeFalling = 0;
 279               tmp->fLastCharge = seq_charge;
 280             }
 281
 282           //Update list of pointers to previous pad:
 283           currentPt[n_current] = &clusterlist[n_total];
 284           n_total++;
 285           n_current++;
 286         }
 287       if(fDeconvTime)
 288         if(new_bin >= 0) goto redo;
 289
 290     }//Loop over digits on this padrow
 291
 292   WriteClusters(n_total,clusterlist);
 293
 294
 295 }
 296
 297 void AliL3ClustFinderNew::WriteClusters(Int_t n_clusters,ClusterData *list)
 298 {
 299   Int_t thisrow,thissector;
 300   UInt_t counter = fNClusters;
 301
 302   for(int j=0; j<n_clusters; j++)
 303     {
 304       if(!list[j].fFlags) continue; //discard 1 pad clusters
 305       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
 306       Float_t xyz[3];
 307
 308       Float_t fpad=(Float_t)list[j].fPad/(Float_t)list[j].fTotalCharge;
 309       Float_t ftime=(Float_t)list[j].fTime/(Float_t)list[j].fTotalCharge;
 310       //printf("padrow %d number of pads %d totalcharge %d\n",fCurrentRow,list[j].fFlags,list[j].fTotalCharge);
 311       fTransform->Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
 312       fTransform->Raw2Local(xyz,thissector,thisrow,fpad,ftime);
 313       if(xyz[0]==0) LOG(AliL3Log::kError,"AliL3ClustFinder","Cluster Finder")
 314         <<AliL3Log::kDec<<"Zero cluster"<<ENDLOG;
 315       if(fNClusters >= fMaxNClusters)
 316         {
 317           LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder")
 318             <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG;
 319           return;
 320           }
 321       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
 322       fSpacePointData[counter].fX = xyz[0];
 323       fSpacePointData[counter].fY = xyz[1];
 324       fSpacePointData[counter].fZ = xyz[2];
 325       fSpacePointData[counter].fPadRow = fCurrentRow;
 326       fSpacePointData[counter].fXYErr = fXYErr;
 327       fSpacePointData[counter].fZErr = fZErr;
 328       fSpacePointData[counter].fID = counter
 329         +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);//uli
 330
 331       fNClusters++;
 332       counter++;
 333
 334     }
 335
 336 }