HLT/src/AliL3ClustFinderNew.cxx

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