HLT/TPCLib/Ref/AliHLTTPCClustFinderNew.cxx

   1 // @(#) $Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>, Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
   4 //*-- Copyright &copy ALICE HLT Group
   5
   6 #include "AliHLTTPCStandardIncludes.h"
   7
   8 #include "AliHLTTPCLogging.h"
   9 #include "AliHLTTPCClustFinderNew.h"
  10 #include "AliHLTTPCDigitData.h"
  11 #include "AliHLTTPCTransform.h"
  12 #include "AliHLTTPCSpacePointData.h"
  13 #include "AliHLTTPCMemHandler.h"
  14
  15 #if GCCVERSION == 3
  16 using namespace std;
  17 #endif
  18
  19 /** \class AliHLTTPCClustFinderNew
  20 <pre>
  21 //_____________________________________________________________
  22 // AliHLTTPCClustFinderNew
  23 //
  24 // The current cluster finder for HLT
  25 // (Based on STAR L3)
  26 //
  27 // The cluster finder is initialized with the Init function,
  28 // providing the slice and patch information to work on.
  29 // The input is a AliHLTTPCDigitRowData structure using the
  30 // Read function. The resulting space points will be in the
  31 // array given by the SetOutputArray function.
  32 //
  33 // There are several setters which control the behaviour:
  34 //
  35 // - SetXYError(Float_t):   set fixed error in XY direction
  36 // - SetZError(Float_t):    set fixed error in Z  direction
  37 //                            (used if errors are not calculated)
  38 // - SetDeconv(Bool_t):     switch on/off deconvolution
  39 // - SetThreshold(UInt_t):  set charge threshold for cluster
  40 // - SetMatchWidth(UInt_t): set the match distance in
  41 //                            time for sequences to be merged
  42 // - SetSTDOutput(Bool_t):  switch on/off output about found clusters
  43 // - SetCalcErr(Bool_t):    switch on/off calculation of
  44 //                          space point errors (or widths in raw system)
  45 // - SetRawSP(Bool_t):      switch on/off convertion to raw system
  46 //
  47 //
  48 // Example Usage:
  49 //
  50 // AliHLTTPCFileHandler *file = new AliHLTTPCFileHandler();
  51 // file->SetAliInput(digitfile); //give some input file
  52 // for(int slice=0; slice<=35; slice++){
  53 //   for(int patch=0; pat<6; pat++){
  54 //     file->Init(slice,patch);
  55 //     UInt_t ndigits=0;
  56 //     UInt_t maxclusters=100000;
  57 //     UInt_t pointsize = maxclusters*sizeof(AliHLTTPCSpacePointData);
  58 //     AliHLTTPCSpacePointData *points = (AliHLTTPCSpacePointData*)memory->Allocate(pointsize);
  59 //     AliHLTTPCDigitRowData *digits = (AliHLTTPCDigitRowData*)file->AliAltroDigits2Memory(ndigits,event);
  60 //     AliHLTTPCClustFinderNew *cf = new AliHLTTPCClustFinderNew();
  61 //     cf->SetMatchWidth(2);
  62 //     cf->InitSlice(slice,patch,maxclusters);
  63 //     cf->SetSTDOutput(kTRUE);    //Some output to standard IO
  64 //     cf->SetRawSP(kFALSE);       //Convert space points to local system
  65 //     cf->SetThreshold(5);        //Threshold of cluster charge
  66 //     cf->SetDeconv(kTRUE);       //Deconv in pad and time direction
  67 //     cf->SetCalcErr(kTRUE);      //Calculate the errors of the spacepoints
  68 //     cf->SetOutputArray(points); //Move the spacepoints to the array
  69 //     cf->Read(ndigits,digits);   //give the data to the cf
  70 //     cf->ProcessDigits();        //process the rows given by init
  71 //     Int_t npoints = cf->GetNumberOfClusters();
  72 //     AliHLTTPCMemHandler *out= new AliHLTTPCMemHandler();
  73 //     out->SetBinaryOutput(fname);
  74 //     out->Memory2Binary(npoints,points); //store the spacepoints
  75 //     out->CloseBinaryOutput();
  76 //     delete out;
  77 //     file->free();
  78 //     delete cf;
  79 //   }
  80 // }
  81 </pre>
  82 */
  83
  84 ClassImp(AliHLTTPCClustFinderNew)
  85
  86 AliHLTTPCClustFinderNew::AliHLTTPCClustFinderNew()
  87 {
  88   fMatch = 1;
  89   fThreshold = 10;
  90   fXYErr = 0.2;
  91   fZErr = 0.3;
  92   fDeconvPad = kTRUE;
  93   fDeconvTime = kTRUE;
  94   fStdout = kFALSE;
  95   fCalcerr = kTRUE;
  96   fRawSP = kFALSE;
  97   fFirstRow=0;
  98   fLastRow=0;
  99 }
 100
 101 AliHLTTPCClustFinderNew::~AliHLTTPCClustFinderNew()
 102 {
 103 }
 104
 105 void AliHLTTPCClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
 106 {
 107   fNClusters = 0;
 108   fMaxNClusters = nmaxpoints;
 109   fCurrentSlice = slice;
 110   fCurrentPatch = patch;
 111   fFirstRow = firstrow;
 112   fLastRow = lastrow;
 113 }
 114
 115 void AliHLTTPCClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
 116 {
 117   fNClusters = 0;
 118   fMaxNClusters = nmaxpoints;
 119   fCurrentSlice = slice;
 120   fCurrentPatch = patch;
 121   fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
 122   fLastRow=AliHLTTPCTransform::GetLastRow(patch);
 123 }
 124
 125 void AliHLTTPCClustFinderNew::SetOutputArray(AliHLTTPCSpacePointData *pt)
 126 {
 127   fSpacePointData = pt;
 128 }
 129
 130 void AliHLTTPCClustFinderNew::Read(UInt_t ndigits,AliHLTTPCDigitRowData *ptr)
 131 {
 132   fNDigitRowData = ndigits;
 133   fDigitRowData = ptr;
 134 }
 135
 136 void AliHLTTPCClustFinderNew::ProcessDigits()
 137 {
 138   //Loop over rows, and call processrow
 139
 140   AliHLTTPCDigitRowData *tempPt = (AliHLTTPCDigitRowData*)fDigitRowData;
 141
 142   for(Int_t i=fFirstRow; i<=fLastRow; i++)
 143     {
 144       fCurrentRow = i;
 145       if((Int_t)tempPt->fRow!=fCurrentRow){
 146         LOG(AliHLTTPCLog::kWarning,"AliHLTTPCClustFinderNew::ProcessDigits","Digits")
 147           <<"Row number should match! "<<tempPt->fRow<<" "<<fCurrentRow<<ENDLOG;
 148         continue;
 149       }
 150 #if 0
 151       LOG(AliHLTTPCLog::kDebug,"AliHLTTPCClustFinderNew::ProcessDigits","Digits")
 152           << "Row " << AliHLTTPCLog::kDec << tempPt->fRow << " digits: "
 153           << tempPt->fNDigit << " (Offset: " << ((unsigned long)tempPt) - ((unsigned long)fDigitRowData)
 154           << ")." << ENDLOG;
 155 #endif
 156       ProcessRow(tempPt);
 157       Byte_t *tmp = (Byte_t*)tempPt;
 158       Int_t size = sizeof(AliHLTTPCDigitRowData) + tempPt->fNDigit*sizeof(AliHLTTPCDigitData);
 159       tmp += size;
 160       tempPt = (AliHLTTPCDigitRowData*)tmp;
 161     }
 162   LOG(AliHLTTPCLog::kDebug,"AliHLTTPCClustFinderNew::ProcessDigits","Space points")
 163     <<"Cluster finder found "<<fNClusters<<" clusters in slice "<<fCurrentSlice
 164     <<" patch "<<fCurrentPatch<<ENDLOG;
 165 }
 166
 167 void AliHLTTPCClustFinderNew::ProcessRow(AliHLTTPCDigitRowData *tempPt)
 168 {
 169
 170   UInt_t last_pad = 123456789;
 171
 172   ClusterData *pad1[5000]; //2 lists for internal memory=2pads
 173   ClusterData *pad2[5000]; //2 lists for internal memory=2pads
 174   ClusterData clusterlist[10000]; //Clusterlist
 175
 176   ClusterData **currentPt;  //List of pointers to the current pad
 177   ClusterData **previousPt; //List of pointers to the previous pad
 178   currentPt = pad2;
 179   previousPt = pad1;
 180   UInt_t n_previous=0,n_current=0,n_total=0;
 181
 182   //Loop over sequences of this row:
 183   for(UInt_t bin=0; bin<tempPt->fNDigit; bin++)
 184     {
 185       AliHLTTPCDigitData *data = tempPt->fDigitData;
 186 #if 0
 187       LOG( AliHLTTPCLog::kDebug, "AliHLTTPCClustFinderNew::ProcessRow", "Signal" )
 188           << "Pad " << AliHLTTPCLog::kDec << data[bin].fPad
 189           << " (" << bin << ")"
 190           << " time " << data[bin].fTime << ": " << data[bin].fCharge
 191           << " (Offset: " << ((unsigned long)(data+bin)) - ((unsigned long)fDigitRowData)
 192           << ")." << ENDLOG;
 193 #ifdef do_mc
 194       LOG( AliHLTTPCLog::kWarning, "AliHLTTPCClustFinderNew::ProcessRow", "Signal" )
 195           << "Compiled with do_mc" << ENDLOG;
 196 #endif
 197 #endif
 198       if(data[bin].fPad != last_pad)
 199         {
 200           //This is a new pad
 201
 202           //Switch the lists:
 203           if(currentPt == pad2)
 204             {
 205               currentPt = pad1;
 206               previousPt = pad2;
 207             }
 208           else
 209             {
 210               currentPt = pad2;
 211               previousPt = pad1;
 212             }
 213           n_previous = n_current;
 214           n_current = 0;
 215           if(bin[data].fPad != last_pad+1)
 216             {
 217               //this happens if there is a pad with no signal.
 218               n_previous = n_current = 0;
 219             }
 220           last_pad = data[bin].fPad;
 221         }
 222
 223       Bool_t new_cluster = kTRUE;
 224       UInt_t seq_charge=0,seq_average=0,seq_error=0;
 225       UInt_t last_charge=0,last_was_falling=0;
 226       Int_t new_bin=-1;
 227
 228       if(fDeconvTime)
 229         {
 230         redo: //This is a goto.
 231           if(new_bin > -1)
 232             {
 233               bin = new_bin;
 234               new_bin = -1;
 235             }
 236
 237           last_charge=0;
 238           last_was_falling = 0;
 239         }
 240
 241       while(1) //Loop over current sequence
 242         {
 243           if(data[bin].fTime >= AliHLTTPCTransform::GetNTimeBins())
 244             {
 245               LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClustFinderNew::ProcessRow","Digits")
 246                 <<"Timebin out of range "<<(Int_t)data[bin].fTime<<ENDLOG;
 247               break;
 248             }
 249
 250           //Get the current ADC-value
 251           UInt_t charge = data[bin].fCharge;
 252
 253           if(fDeconvTime)
 254             {
 255               //Check if the last pixel in the sequence is smaller than this
 256               if(charge > last_charge)
 257                 {
 258                   if(last_was_falling)
 259                     {
 260                       new_bin = bin;
 261                       break;
 262                     }
 263                 }
 264               else last_was_falling = 1; //last pixel was larger than this
 265               last_charge = charge;
 266             }
 267
 268           //Sum the total charge of this sequence
 269           seq_charge += charge;
 270           seq_average += data[bin].fTime*charge;
 271           seq_error += data[bin].fTime*data[bin].fTime*charge;
 272
 273           //Check where to stop:
 274           if(bin >= tempPt->fNDigit - 1) //out of range
 275             break;
 276           if(data[bin+1].fPad != data[bin].fPad) //new pad
 277             break;
 278           if(data[bin+1].fTime != data[bin].fTime+1) //end of sequence
 279             break;
 280
 281           bin++;
 282         }//end loop over sequence
 283
 284       //Calculate mean of sequence:
 285       Int_t seq_mean=0;
 286       if(seq_charge)
 287         seq_mean = seq_average/seq_charge;
 288       else
 289         {
 290           LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClustFinderNew::ProcessRow","Data")
 291             <<"Error in data given to the cluster finder"<<ENDLOG;
 292           seq_mean = 1;
 293           seq_charge = 1;
 294         }
 295
 296       //Calculate mean in pad direction:
 297       Int_t pad_mean = seq_charge*data[bin].fPad;
 298       Int_t pad_error = data[bin].fPad*pad_mean;
 299
 300       //Compare with results on previous pad:
 301       for(UInt_t p=0; p<n_previous; p++)
 302         {
 303           //dont merge sequences on the same pad twice
 304           if(previousPt[p]->fLastMergedPad==data[bin].fPad) continue;
 305
 306           Int_t difference = seq_mean - previousPt[p]->fMean;
 307           if(difference < -fMatch) break;
 308
 309           if(difference <= fMatch) //There is a match here!!
 310             {
 311               ClusterData *local = previousPt[p];
 312
 313               if(fDeconvPad)
 314                 {
 315                   if(seq_charge > local->fLastCharge)
 316                     {
 317                       if(local->fChargeFalling) //The previous pad was falling
 318                         {
 319                           break; //create a new cluster
 320                         }
 321                     }
 322                   else
 323                     local->fChargeFalling = 1;
 324                   local->fLastCharge = seq_charge;
 325                 }
 326
 327               //Don't create a new cluster, because we found a match
 328               new_cluster = kFALSE;
 329
 330               //Update cluster on current pad with the matching one:
 331               local->fTotalCharge += seq_charge;
 332               local->fPad += pad_mean;
 333               local->fPad2 += pad_error;
 334               local->fTime += seq_average;
 335               local->fTime2 += seq_error;
 336               local->fMean = seq_mean;
 337               local->fFlags++; //means we have more than one pad
 338               local->fLastMergedPad = data[bin].fPad;
 339
 340               currentPt[n_current] = local;
 341               n_current++;
 342
 343               break;
 344             } //Checking for match at previous pad
 345         } //Loop over results on previous pad.
 346
 347       if(new_cluster)
 348         {
 349           //Start a new cluster. Add it to the clusterlist, and update
 350           //the list of pointers to clusters in current pad.
 351           //current pad will be previous pad on next pad.
 352
 353           //Add to the clusterlist:
 354           ClusterData *tmp = &clusterlist[n_total];
 355           tmp->fTotalCharge = seq_charge;
 356           tmp->fPad = pad_mean;
 357           tmp->fPad2 = pad_error;
 358           tmp->fTime = seq_average;
 359           tmp->fTime2 = seq_error;
 360           tmp->fMean = seq_mean;
 361           tmp->fFlags = 0;  //flags for single pad clusters
 362           tmp->fLastMergedPad = data[bin].fPad;
 363
 364           if(fDeconvPad)
 365             {
 366               tmp->fChargeFalling = 0;
 367               tmp->fLastCharge = seq_charge;
 368             }
 369
 370           //Update list of pointers to previous pad:
 371           currentPt[n_current] = &clusterlist[n_total];
 372           n_total++;
 373           n_current++;
 374         }
 375
 376       if(fDeconvTime)
 377         if(new_bin >= 0) goto redo;
 378     }//Loop over digits on this padrow
 379
 380   WriteClusters(n_total,clusterlist);
 381 }
 382
 383 void AliHLTTPCClustFinderNew::WriteClusters(Int_t n_clusters,ClusterData *list)
 384 {
 385   Int_t thisrow,thissector;
 386   UInt_t counter = fNClusters;
 387
 388   for(int j=0; j<n_clusters; j++)
 389     {
 390       if(!list[j].fFlags) continue; //discard single pad clusters
 391       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
 392
 393       Float_t xyz[3];
 394       Float_t fpad =(Float_t)list[j].fPad /(Float_t)list[j].fTotalCharge;
 395       Float_t fpad2=fXYErr*fXYErr; //fixed given error
 396       Float_t ftime =(Float_t)list[j].fTime /(Float_t)list[j].fTotalCharge;
 397       Float_t ftime2=fZErr*fZErr;  //fixed given error
 398
 399       if(fCalcerr) { //calc the errors, otherwice take the fixed error
 400         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
 401
 402         Float_t sy2=(Float_t)list[j].fPad2/(Float_t)list[j].fTotalCharge - fpad*fpad;
 403         if(sy2 < 0) {
 404             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinderNew::WriteClusters","Cluster width")
 405               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
 406             continue;
 407         } else {
 408           if(!fRawSP){
 409             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
 410             if(sy2 != 0){
 411               fpad2*=0.108; //constants are from offline studies
 412               if(patch<2)
 413                 fpad2*=2.07;
 414             }
 415           } else fpad2=sy2; //take the width not the error
 416         }
 417         Float_t sz2=(Float_t)list[j].fTime2/(Float_t)list[j].fTotalCharge - ftime*ftime;
 418         if(sz2 < 0){
 419           LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinderNew::WriteClusters","Cluster width")
 420             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
 421           continue;
 422         } else {
 423           if(!fRawSP){
 424             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
 425             if(sz2 != 0) {
 426               ftime2 *= 0.169; //constants are from offline studies
 427               if(patch<2)
 428                 ftime2 *= 1.77;
 429             }
 430           } else ftime2=sz2; //take the width, not the error
 431         }
 432       }
 433       if(fStdout==kTRUE)
 434           {
 435 #if 0
 436         cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
 437 #endif
 438           }
 439
 440       if(!fRawSP){
 441         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
 442         AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
 443
 444         if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
 445           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
 446         if(fNClusters >= fMaxNClusters)
 447           {
 448             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
 449               <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
 450             return;
 451           }
 452
 453         fSpacePointData[counter].fX = xyz[0];
 454         fSpacePointData[counter].fY = xyz[1];
 455         fSpacePointData[counter].fZ = xyz[2];
 456
 457       } else {
 458         fSpacePointData[counter].fX = fCurrentRow;
 459         fSpacePointData[counter].fY = fpad;
 460         fSpacePointData[counter].fZ = ftime;
 461       }
 462
 463       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
 464       fSpacePointData[counter].fPadRow = fCurrentRow;
 465       fSpacePointData[counter].fSigmaY2 = fpad2;
 466       fSpacePointData[counter].fSigmaZ2  = ftime2;
 467
 468       Int_t patch=fCurrentPatch;
 469       if(patch==-1) patch=0; //never store negative patch number
 470       fSpacePointData[counter].fID = counter
 471         +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
 472 #ifdef do_mc
 473       Int_t trackID[3];
 474       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
 475
 476       fSpacePointData[counter].fTrackID[0] = trackID[0];
 477       fSpacePointData[counter].fTrackID[1] = trackID[1];
 478       fSpacePointData[counter].fTrackID[2] = trackID[2];
 479
 480       //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
 481 #endif
 482
 483       fNClusters++;
 484       counter++;
 485     }
 486 }
 487
 488 #ifdef do_mc
 489 void AliHLTTPCClustFinderNew::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
 490 {
 491   AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
 492
 493   trackID[0]=trackID[1]=trackID[2]=-2;
 494   //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
 495   for(Int_t i=fFirstRow; i<=fLastRow; i++){
 496     if(rowPt->fRow < (UInt_t)fCurrentRow){
 497       AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
 498       continue;
 499     }
 500     AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
 501     for(UInt_t j=0; j<rowPt->fNDigit; j++){
 502       Int_t cpad = digPt[j].fPad;
 503       Int_t ctime = digPt[j].fTime;
 504       if(cpad != pad) continue;
 505       if(ctime != time) continue;
 506
 507       trackID[0] = digPt[j].fTrackID[0];
 508       trackID[1] = digPt[j].fTrackID[1];
 509       trackID[2] = digPt[j].fTrackID[2];
 510
 511       //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
 512       break;
 513     }
 514     break;
 515   }
 516 }
 517 #endif