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           //Get the current ADC-value
 167           UInt_t charge = data[bin].fCharge;
 168
 169           if(fDeconvTime)
 170             {
 171               //Check if the last pixel in the sequence is smaller than this
 172               if(charge > last_charge)
 173                 {
 174                   if(last_was_falling)
 175                     {
 176                       new_bin = bin;
 177                       break;
 178                     }
 179                 }
 180               else last_was_falling = 1; //last pixel was larger than this
 181               last_charge = charge;
 182             }
 183
 184           //Sum the total charge of this sequence
 185           seq_charge += charge;
 186           seq_average += data[bin].fTime*charge;
 187
 188           //Check where to stop:
 189           if(data[bin+1].fPad != data[bin].fPad) //new pad
 190             break;
 191           if(data[bin+1].fTime != data[bin].fTime+1) //end of sequence
 192             break;
 193
 194           bin++;
 195         }//end loop over sequence
 196
 197       //Calculate mean of sequence:
 198       Int_t seq_mean=0;
 199       if(seq_charge)
 200         seq_mean = seq_average/seq_charge;
 201       else
 202         {
 203           printf("\nCF: Should not happen\n");
 204           seq_mean = 1;
 205           seq_charge = 1;
 206         }
 207
 208       //Calculate mean in pad direction:
 209       Int_t pad_mean = seq_charge*data[bin].fPad;
 210
 211       //Compare with results on previous pad:
 212       for(UInt_t p=0; p<n_previous; p++)
 213         {
 214           Int_t difference = seq_mean - previousPt[p]->fMean;
 215           if(difference < -fMatch)
 216             break;
 217
 218           if(difference <= fMatch)//There is a match here!!
 219             {
 220
 221               ClusterData *local = previousPt[p];
 222               if(fDeconvPad)
 223                 {
 224                   if(seq_charge > local->fLastCharge)
 225                     {
 226                       if(local->fChargeFalling)//The previous pad was falling
 227                         {
 228                           break;//create a new cluster
 229                         }
 230                     }
 231                   else
 232                     local->fChargeFalling = 1;
 233                   local->fLastCharge = seq_charge;
 234                 }
 235
 236               //Don't create a new cluster, because we found a match
 237               new_cluster = kFALSE;
 238
 239               //Update cluster on current pad with the matching one:
 240               //currentPt[n_current] = previousPt[p];
 241
 242               local->fTotalCharge += seq_charge;
 243               local->fPad += pad_mean;
 244               local->fTime += seq_average;
 245               local->fMean = seq_mean;
 246               local->fFlags++; //means we have more than one pad
 247
 248               currentPt[n_current] = local;
 249               n_current++;
 250
 251               break;
 252             }//Checking for match at previous pad
 253         }//Loop over results on previous pad.
 254
 255       if(new_cluster)
 256         {
 257           //Start a new cluster. Add it to the clusterlist, and update
 258           //the list of pointers to clusters in current pad.
 259           //current pad will be previous pad on next pad.
 260
 261           //Add to the clusterlist:
 262           ClusterData *tmp = &clusterlist[n_total];
 263           tmp->fTotalCharge = seq_charge;
 264           tmp->fPad = pad_mean;
 265           tmp->fTime = seq_average;
 266           tmp->fMean = seq_mean;
 267           tmp->fFlags = 0;  //flags for 1 pad clusters
 268           if(fDeconvPad)
 269             {
 270               tmp->fChargeFalling = 0;
 271               tmp->fLastCharge = seq_charge;
 272             }
 273
 274           //Update list of pointers to previous pad:
 275           currentPt[n_current] = &clusterlist[n_total];
 276           n_total++;
 277           n_current++;
 278         }
 279       if(fDeconvTime)
 280         if(new_bin >= 0) goto redo;
 281
 282     }//Loop over digits on this padrow
 283
 284   WriteClusters(n_total,clusterlist);
 285
 286
 287 }
 288
 289 void AliL3ClustFinderNew::WriteClusters(Int_t n_clusters,ClusterData *list)
 290 {
 291   Int_t thisrow,thissector;
 292   UInt_t counter = fNClusters;
 293
 294   for(int j=0; j<n_clusters; j++)
 295     {
 296       if(!list[j].fFlags) continue; //discard 1 pad clusters
 297       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
 298       Float_t xyz[3];
 299
 300       Float_t fpad=(Float_t)list[j].fPad/(Float_t)list[j].fTotalCharge;
 301       Float_t ftime=(Float_t)list[j].fTime/(Float_t)list[j].fTotalCharge;
 302       //printf("padrow %d number of pads %d totalcharge %d\n",fCurrentRow,list[j].fFlags,list[j].fTotalCharge);
 303       fTransform->Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
 304       fTransform->Raw2Local(xyz,thissector,thisrow,fpad,ftime);
 305       if(xyz[0]==0) LOG(AliL3Log::kError,"AliL3ClustFinder","Cluster Finder")
 306         <<AliL3Log::kDec<<"Zero cluster"<<ENDLOG;
 307       if(fNClusters >= fMaxNClusters)
 308         {
 309           LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder")
 310             <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG;
 311           return;
 312           }
 313       fSpacePointData[counter].fX = xyz[0];
 314       fSpacePointData[counter].fY = xyz[1];
 315       fSpacePointData[counter].fZ = xyz[2];
 316       fSpacePointData[counter].fPadRow = fCurrentRow;
 317       fSpacePointData[counter].fXYErr = fXYErr;
 318       fSpacePointData[counter].fZErr = fZErr;
 319       fSpacePointData[counter].fID = counter
 320         +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);//uli
 321
 322       fNClusters++;
 323       counter++;
 324
 325     }
 326
 327 }