HLT/TPCLib/AliHLTTPCClustFinderNew.cxx

   1 // @(#) $Id$
   2 //Original: AliHLTClustFinderNew.cxx,v 1.29 2005/06/14 10:55:21 cvetan
   3
   4 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>, Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
   5 //*-- Copyright &copy ALICE HLT Group
   6
   7 #include "AliHLTTPCRootTypes.h"
   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 __GNUC__ >= 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   //constructor
  89   fMatch = 1;
  90   fThreshold = 10;
  91   fXYErr = 0.2;
  92   fZErr = 0.3;
  93   fDeconvPad = kTRUE;
  94   fDeconvTime = kTRUE;
  95   fStdout = kFALSE;
  96   fCalcerr = kTRUE;
  97   fRawSP = kFALSE;
  98   fFirstRow=0;
  99   fLastRow=0;
 100 }
 101
 102 AliHLTTPCClustFinderNew::~AliHLTTPCClustFinderNew()
 103 {
 104   //destructor
 105   ;
 106 }
 107
 108 void AliHLTTPCClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
 109 {
 110   //init slice
 111   fNClusters = 0;
 112   fMaxNClusters = nmaxpoints;
 113   fCurrentSlice = slice;
 114   fCurrentPatch = patch;
 115   fFirstRow = firstrow;
 116   fLastRow = lastrow;
 117 }
 118
 119 void AliHLTTPCClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
 120 {
 121   //init slice
 122   fNClusters = 0;
 123   fMaxNClusters = nmaxpoints;
 124   fCurrentSlice = slice;
 125   fCurrentPatch = patch;
 126   fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
 127   fLastRow=AliHLTTPCTransform::GetLastRow(patch);
 128 }
 129
 130 void AliHLTTPCClustFinderNew::SetOutputArray(AliHLTTPCSpacePointData *pt)
 131 {
 132   //set pointer to output
 133   fSpacePointData = pt;
 134 }
 135
 136 void AliHLTTPCClustFinderNew::Read(UInt_t ndigits,AliHLTTPCDigitRowData *ptr)
 137 {
 138   //set input pointer
 139   fNDigitRowData = ndigits;
 140   fDigitRowData = ptr;
 141 }
 142
 143 void AliHLTTPCClustFinderNew::ProcessDigits()
 144 {
 145   //Loop over rows, and call processrow
 146   AliHLTTPCDigitRowData *tempPt = (AliHLTTPCDigitRowData*)fDigitRowData;
 147   fNClusters = 0;
 148   for(Int_t i=fFirstRow; i<=fLastRow; i++)
 149     {
 150       fCurrentRow = i;
 151       if((Int_t)tempPt->fRow!=fCurrentRow){
 152         LOG(AliHLTTPCLog::kWarning,"AliHLTTPCClustFinderNew::ProcessDigits","Digits")
 153           <<"Row number should match! "<<tempPt->fRow<<" "<<fCurrentRow<<ENDLOG;
 154         continue;
 155       }
 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::kInformational,"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   //process row
 170   UInt_t lastpad = 123456789;
 171
 172   AliClusterData *pad1[5000]; //2 lists for internal memory=2pads
 173   AliClusterData *pad2[5000]; //2 lists for internal memory=2pads
 174   AliClusterData clusterlist[10000]; //Clusterlist
 175
 176   AliClusterData **currentPt;  //List of pointers to the current pad
 177   AliClusterData **previousPt; //List of pointers to the previous pad
 178   currentPt = pad2;
 179   previousPt = pad1;
 180   UInt_t nprevious=0,ncurrent=0,ntotal=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(data[bin].fPad != lastpad)
 187         {
 188           //This is a new pad
 189
 190           //Switch the lists:
 191           if(currentPt == pad2)
 192             {
 193               currentPt = pad1;
 194               previousPt = pad2;
 195             }
 196           else
 197             {
 198               currentPt = pad2;
 199               previousPt = pad1;
 200             }
 201           nprevious = ncurrent;
 202           ncurrent = 0;
 203           if(bin[data].fPad != lastpad+1)
 204             {
 205               //this happens if there is a pad with no signal.
 206               nprevious = ncurrent = 0;
 207             }
 208           lastpad = data[bin].fPad;
 209         }
 210
 211       Bool_t newcluster = kTRUE;
 212       UInt_t seqcharge=0,seqaverage=0,seqerror=0;
 213       UInt_t lastcharge=0,lastwas_falling=0;
 214       Int_t newbin=-1;
 215
 216       if(fDeconvTime)
 217         {
 218         redo: //This is a goto.
 219           if(newbin > -1)
 220             {
 221               bin = newbin;
 222               newbin = -1;
 223             }
 224
 225           lastcharge=0;
 226           lastwas_falling = 0;
 227         }
 228
 229       while(1) //Loop over current sequence
 230         {
 231           if(data[bin].fTime >= AliHLTTPCTransform::GetNTimeBins())
 232             {
 233               LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClustFinderNew::ProcessRow","Digits")
 234                 <<"Timebin out of range "<<(Int_t)data[bin].fTime<<ENDLOG;
 235               break;
 236             }
 237
 238           //Get the current ADC-value
 239           UInt_t charge = data[bin].fCharge;
 240
 241           if(fDeconvTime)
 242             {
 243               //Check if the last pixel in the sequence is smaller than this
 244               if(charge > lastcharge)
 245                 {
 246                   if(lastwas_falling)
 247                     {
 248                       newbin = bin;
 249                       break;
 250                     }
 251                 }
 252               else lastwas_falling = 1; //last pixel was larger than this
 253               lastcharge = charge;
 254             }
 255
 256           //Sum the total charge of this sequence
 257           seqcharge += charge;
 258           seqaverage += data[bin].fTime*charge;
 259           seqerror += data[bin].fTime*data[bin].fTime*charge;
 260
 261           //Check where to stop:
 262           if(bin >= tempPt->fNDigit - 1) //out of range
 263             break;
 264           if(data[bin+1].fPad != data[bin].fPad) //new pad
 265             break;
 266           if(data[bin+1].fTime != data[bin].fTime+1) //end of sequence
 267             break;
 268
 269           bin++;
 270         }//end loop over sequence
 271
 272       //Calculate mean of sequence:
 273       Int_t seqmean=0;
 274       if(seqcharge)
 275         seqmean = seqaverage/seqcharge;
 276       else
 277         {
 278           LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClustFinderNew::ProcessRow","Data")
 279             <<"Error in data given to the cluster finder"<<ENDLOG;
 280           seqmean = 1;
 281           seqcharge = 1;
 282         }
 283
 284       //Calculate mean in pad direction:
 285       Int_t padmean = seqcharge*data[bin].fPad;
 286       Int_t paderror = data[bin].fPad*padmean;
 287
 288       //Compare with results on previous pad:
 289       for(UInt_t p=0; p<nprevious; p++)
 290         {
 291           //dont merge sequences on the same pad twice
 292           if(previousPt[p]->fLastMergedPad==data[bin].fPad) continue;
 293
 294           Int_t difference = seqmean - previousPt[p]->fMean;
 295           if(difference < -fMatch) break;
 296
 297           if(difference <= fMatch) //There is a match here!!
 298             {
 299               AliClusterData *local = previousPt[p];
 300
 301               if(fDeconvPad)
 302                 {
 303                   if(seqcharge > local->fLastCharge)
 304                     {
 305                       if(local->fChargeFalling) //The previous pad was falling
 306                         {
 307                           break; //create a new cluster
 308                         }
 309                     }
 310                   else
 311                     local->fChargeFalling = 1;
 312                   local->fLastCharge = seqcharge;
 313                 }
 314
 315               //Don't create a new cluster, because we found a match
 316               newcluster = kFALSE;
 317
 318               //Update cluster on current pad with the matching one:
 319               local->fTotalCharge += seqcharge;
 320               local->fPad += padmean;
 321               local->fPad2 += paderror;
 322               local->fTime += seqaverage;
 323               local->fTime2 += seqerror;
 324               local->fMean = seqmean;
 325               local->fFlags++; //means we have more than one pad
 326               local->fLastMergedPad = data[bin].fPad;
 327
 328               currentPt[ncurrent] = local;
 329               ncurrent++;
 330
 331               break;
 332             } //Checking for match at previous pad
 333         } //Loop over results on previous pad.
 334
 335       if(newcluster)
 336         {
 337           //Start a new cluster. Add it to the clusterlist, and update
 338           //the list of pointers to clusters in current pad.
 339           //current pad will be previous pad on next pad.
 340
 341           //Add to the clusterlist:
 342           AliClusterData *tmp = &clusterlist[ntotal];
 343           tmp->fTotalCharge = seqcharge;
 344           tmp->fPad = padmean;
 345           tmp->fPad2 = paderror;
 346           tmp->fTime = seqaverage;
 347           tmp->fTime2 = seqerror;
 348           tmp->fMean = seqmean;
 349           tmp->fFlags = 0;  //flags for single pad clusters
 350           tmp->fLastMergedPad = data[bin].fPad;
 351
 352           if(fDeconvPad)
 353             {
 354               tmp->fChargeFalling = 0;
 355               tmp->fLastCharge = seqcharge;
 356             }
 357
 358           //Update list of pointers to previous pad:
 359           currentPt[ncurrent] = &clusterlist[ntotal];
 360           ntotal++;
 361           ncurrent++;
 362         }
 363
 364       if(fDeconvTime)
 365         if(newbin >= 0) goto redo;
 366     }//Loop over digits on this padrow
 367
 368   WriteClusters(ntotal,clusterlist);
 369 }
 370
 371 void AliHLTTPCClustFinderNew::WriteClusters(Int_t nclusters,AliClusterData *list)
 372 {
 373   //write cluster to output pointer
 374   Int_t thisrow,thissector;
 375   UInt_t counter = fNClusters;
 376
 377   for(int j=0; j<nclusters; j++)
 378     {
 379       if(!list[j].fFlags) continue; //discard single pad clusters
 380       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
 381
 382       Float_t xyz[3];
 383       Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
 384       Float_t fpad2=fXYErr*fXYErr; //fixed given error
 385       Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
 386       Float_t ftime2=fZErr*fZErr;  //fixed given error
 387
 388       if(fCalcerr) { //calc the errors, otherwice take the fixed error
 389         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
 390         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
 391         Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
 392         sy2/=q2;
 393         if(sy2 < 0) {
 394             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinderNew::WriteClusters","Cluster width")
 395               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
 396             continue;
 397         } else {
 398           if(!fRawSP){
 399             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
 400             if(sy2 != 0){
 401               fpad2*=0.108; //constants are from offline studies
 402               if(patch<2)
 403                 fpad2*=2.07;
 404             }
 405           } else fpad2=sy2; //take the width not the error
 406         }
 407         Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
 408         sz2/=q2;
 409         if(sz2 < 0){
 410           LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinderNew::WriteClusters","Cluster width")
 411             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
 412           continue;
 413         } else {
 414           if(!fRawSP){
 415             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
 416             if(sz2 != 0) {
 417               ftime2 *= 0.169; //constants are from offline studies
 418               if(patch<2)
 419                 ftime2 *= 1.77;
 420             }
 421           } else ftime2=sz2; //take the width, not the error
 422         }
 423       }
 424       if(fStdout==kTRUE)
 425         cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
 426
 427       if(!fRawSP){
 428         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
 429         AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
 430
 431         if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
 432           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
 433         if(fNClusters >= fMaxNClusters)
 434           {
 435             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
 436               <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
 437             return;
 438           }
 439
 440         fSpacePointData[counter].fX = xyz[0];
 441         fSpacePointData[counter].fY = xyz[1];
 442         fSpacePointData[counter].fZ = xyz[2];
 443
 444       } else {
 445         fSpacePointData[counter].fX = fCurrentRow;
 446         fSpacePointData[counter].fY = fpad;
 447         fSpacePointData[counter].fZ = ftime;
 448       }
 449
 450       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
 451       fSpacePointData[counter].fPadRow = fCurrentRow;
 452       fSpacePointData[counter].fSigmaY2 = fpad2;
 453       fSpacePointData[counter].fSigmaZ2  = ftime2;
 454
 455       Int_t patch=fCurrentPatch;
 456       if(patch==-1) patch=0; //never store negative patch number
 457       fSpacePointData[counter].fID = counter
 458         +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
 459 #ifdef do_mc
 460       Int_t trackID[3];
 461       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
 462
 463       fSpacePointData[counter].fTrackID[0] = trackID[0];
 464       fSpacePointData[counter].fTrackID[1] = trackID[1];
 465       fSpacePointData[counter].fTrackID[2] = trackID[2];
 466
 467       //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
 468 #endif
 469
 470       fNClusters++;
 471       counter++;
 472     }
 473 }
 474
 475 #ifdef do_mc
 476 void AliHLTTPCClustFinderNew::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
 477 {
 478   //get mc id
 479   AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
 480
 481   trackID[0]=trackID[1]=trackID[2]=-2;
 482   //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
 483   for(Int_t i=fFirstRow; i<=fLastRow; i++){
 484     if(rowPt->fRow < (UInt_t)fCurrentRow){
 485       AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
 486       continue;
 487     }
 488     AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
 489     for(UInt_t j=0; j<rowPt->fNDigit; j++){
 490       Int_t cpad = digPt[j].fPad;
 491       Int_t ctime = digPt[j].fTime;
 492       if(cpad != pad) continue;
 493       if(ctime != time) continue;
 494
 495       trackID[0] = digPt[j].fTrackID[0];
 496       trackID[1] = digPt[j].fTrackID[1];
 497       trackID[2] = digPt[j].fTrackID[2];
 498
 499       //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
 500       break;
 501     }
 502     break;
 503   }
 504 }
 505 #endif