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