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