HLT/TPCLib/AliHLTTPCClusterFinder.cxx

   1 // @(#) $Id$
   2 // Original: AliHLTClustFinderNew.cxx,v 1.29 2005/06/14 10:55:21 cvetan Exp
   3
   4 /**************************************************************************
   5  * This file is property of and copyright by the ALICE HLT Project        *
   6  * ALICE Experiment at CERN, All rights reserved.                         *
   7  *                                                                        *
   8  * Primary Authors: Anders Vestbo, Constantin Loizides, Jochen Thaeder    *
   9  *                  Kenneth Aamodt <kenneth.aamodt@student.uib.no>        *
  10  *                  for The ALICE HLT Project.                            *
  11  *                                                                        *
  12  * Permission to use, copy, modify and distribute this software and its   *
  13  * documentation strictly for non-commercial purposes is hereby granted   *
  14  * without fee, provided that the above copyright notice appears in all   *
  15  * copies and that both the copyright notice and this permission notice   *
  16  * appear in the supporting documentation. The authors make no claims     *
  17  * about the suitability of this software for any purpose. It is          *
  18  * provided "as is" without express or implied warranty.                  *
  19  **************************************************************************/
  20
  21 /** @file   AliHLTTPCClusterFinder.cxx
  22     @author Anders Vestbo, Constantin Loizides, Jochen Thaeder
  23             Kenneth Aamodt kenneth.aamodt@student.uib.no
  24     @date
  25     @brief  Cluster Finder for the TPC
  26 */
  27
  28 #include "AliHLTTPCDigitReader.h"
  29 #include "AliHLTTPCRootTypes.h"
  30 #include "AliHLTTPCLogging.h"
  31 #include "AliHLTTPCClusterFinder.h"
  32 #include "AliHLTTPCDigitData.h"
  33 #include "AliHLTTPCTransform.h"
  34 #include "AliHLTTPCSpacePointData.h"
  35 #include "AliHLTTPCMemHandler.h"
  36 #include "AliHLTTPCPad.h"
  37 #include <sys/time.h>
  38
  39 #if __GNUC__ >= 3
  40 using namespace std;
  41 #endif
  42
  43 /** \class AliHLTTPCClusterFinder
  44 //
  45 // The current cluster finder for HLT
  46 // (Based on STAR L3)
  47 //
  48 //Basically we have two versions for the cluster finder now.
  49 //The default version, reads the data pad by pad, and find the
  50 //clusters as it reads the data. The other version has now been
  51 //developed to cope with unsorted data. New methods for the unsorted
  52 //version can  be found at the end of the default one i the source file.
  53 //Currently the new version is only build to manage zero-suppressed data.
  54 //More functionality will be added later.
  55 //
  56 // The cluster finder is initialized with the Init function,
  57 // providing the slice and patch information to work on.
  58 //
  59 // The input is a provided by the AliHLTTPCDigitReader class,
  60 // using the init() funktion, and the next() funktion in order
  61 // to get the next bin. Either packed or unpacked data can be
  62 // processed, dependent if one uses AliHLTTPCDigitReaderPacked
  63 // class or AliHLTTPCDigitReaderUnpacked class in the
  64 // Clusterfinder Component.
  65 // The resulting space points will be in the
  66 // array given by the SetOutputArray function.
  67 //
  68 // There are several setters which control the behaviour:
  69 //
  70 // - SetXYError(Float_t):   set fixed error in XY direction
  71 // - SetZError(Float_t):    set fixed error in Z  direction
  72 //                            (used if errors are not calculated)
  73 // - SetDeconv(Bool_t):     switch on/off deconvolution
  74 // - SetThreshold(UInt_t):  set charge threshold for cluster
  75 // - SetMatchWidth(UInt_t): set the match distance in
  76 //                            time for sequences to be merged
  77 // - SetSTDOutput(Bool_t):  switch on/off output about found clusters
  78 // - SetCalcErr(Bool_t):    switch on/off calculation of
  79 //                          space point errors (or widths in raw system)
  80 // - SetRawSP(Bool_t):      switch on/off convertion to raw system
  81 //
  82 //
  83 // Example Usage:
  84 //
  85 // AliHLTTPCFileHandler *file = new AliHLTTPCFileHandler();
  86 // file->SetAliInput(digitfile); //give some input file
  87 // for(int slice=0; slice<=35; slice++){
  88 //   for(int patch=0; pat<6; pat++){
  89 //     file->Init(slice,patch);
  90 //     UInt_t ndigits=0;
  91 //     UInt_t maxclusters=100000;
  92 //     UInt_t pointsize = maxclusters*sizeof(AliHLTTPCSpacePointData);
  93 //     AliHLTTPCSpacePointData *points = (AliHLTTPCSpacePointData*)memory->Allocate(pointsize);
  94 //     AliHLTTPCDigitRowData *digits = (AliHLTTPCDigitRowData*)file->AliAltroDigits2Memory(ndigits,event);
  95 //     AliHLTTPCClusterFinder *cf = new AliHLTTPCClusterFinder();
  96 //     cf->SetMatchWidth(2);
  97 //     cf->InitSlice( slice, patch, row[0], row[1], maxPoints );
  98 //     cf->SetSTDOutput(kTRUE);    //Some output to standard IO
  99 //     cf->SetRawSP(kFALSE);       //Convert space points to local system
 100 //     cf->SetThreshold(5);        //Threshold of cluster charge
 101 //     cf->SetDeconv(kTRUE);       //Deconv in pad and time direction
 102 //     cf->SetCalcErr(kTRUE);      //Calculate the errors of the spacepoints
 103 //     cf->SetOutputArray(points); //Move the spacepoints to the array
 104 //     cf->Read(iter->fPtr, iter->fSize ); //give the data to the cf
 105 //     cf->ProcessDigits();        //process the rows given by init
 106 //     Int_t npoints = cf->GetNumberOfClusters();
 107 //     AliHLTTPCMemHandler *out= new AliHLTTPCMemHandler();
 108 //     out->SetBinaryOutput(fname);
 109 //     out->Memory2Binary(npoints,points); //store the spacepoints
 110 //     out->CloseBinaryOutput();
 111 //     delete out;
 112 //     file->free();
 113 //     delete cf;
 114 //   }
 115 // }
 116 */
 117
 118 ClassImp(AliHLTTPCClusterFinder)
 119
 120 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
 121   :
 122   fSpacePointData(NULL),
 123   fDigitReader(NULL),
 124   fPtr(NULL),
 125   fSize(0),
 126   fDeconvTime(kTRUE),
 127   fDeconvPad(kTRUE),
 128   fStdout(kFALSE),
 129   fCalcerr(kTRUE),
 130   fRawSP(kFALSE),
 131   fFirstRow(0),
 132   fLastRow(0),
 133   fCurrentRow(0),
 134   fCurrentSlice(0),
 135   fCurrentPatch(0),
 136   fMatch(1),
 137   fThreshold(10),
 138   fSignalThreshold(-1),
 139   fNSigmaThreshold(0),
 140   fNClusters(0),
 141   fMaxNClusters(0),
 142   fXYErr(0.2),
 143   fZErr(0.3),
 144   fOccupancyLimit(1.0),
 145   fUnsorted(0),
 146   fVectorInitialized(kFALSE),
 147   fRowPadVector(),
 148   fClusters(),
 149   fNumberOfPadsInRow(NULL),
 150   fNumberOfRows(0),
 151   fRowOfFirstCandidate(0)
 152 {
 153   //constructor
 154 }
 155
 156 AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder()
 157 {
 158   //destructor
 159   if(fVectorInitialized){
 160     DeInitializePadArray();
 161   }
 162   if(fNumberOfPadsInRow){
 163     delete [] fNumberOfPadsInRow;
 164     fNumberOfPadsInRow=NULL;
 165   }
 166 }
 167
 168 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
 169 {
 170   //init slice
 171   fNClusters = 0;
 172   fMaxNClusters = nmaxpoints;
 173   fCurrentSlice = slice;
 174   fCurrentPatch = patch;
 175   fFirstRow = firstrow;
 176   fLastRow = lastrow;
 177 }
 178
 179 void AliHLTTPCClusterFinder::InitializePadArray(){
 180   // see header file for class documentation
 181
 182   if(fCurrentPatch>5||fCurrentPatch<0){
 183     HLTFatal("Patch is not set");
 184     return;
 185   }
 186
 187   HLTDebug("Patch number=%d",fCurrentPatch);
 188
 189   fFirstRow = AliHLTTPCTransform::GetFirstRow(fCurrentPatch);
 190   fLastRow = AliHLTTPCTransform::GetLastRow(fCurrentPatch);
 191
 192   fNumberOfRows=fLastRow-fFirstRow+1;
 193   fNumberOfPadsInRow= new UInt_t[fNumberOfRows];
 194
 195   memset( fNumberOfPadsInRow, 0, sizeof(Int_t)*(fNumberOfRows));
 196
 197   for(UInt_t i=0;i<fNumberOfRows;i++){
 198     fNumberOfPadsInRow[i]=AliHLTTPCTransform::GetNPads(i+fFirstRow);
 199     AliHLTTPCPadVector tmpRow;
 200     for(UInt_t j=0;j<fNumberOfPadsInRow[i];j++){
 201       AliHLTTPCPad *tmpPad = new AliHLTTPCPad(2);
 202       tmpPad->SetID(i,j);
 203       tmpRow.push_back(tmpPad);
 204     }
 205     fRowPadVector.push_back(tmpRow);
 206   }
 207   fVectorInitialized=kTRUE;
 208 }
 209
 210 Int_t AliHLTTPCClusterFinder::DeInitializePadArray()
 211 {
 212   // see header file for class documentation
 213   for(UInt_t i=0;i<fNumberOfRows;i++){
 214     for(UInt_t j=0;j<fNumberOfPadsInRow[i];j++){
 215       delete fRowPadVector[i][j];
 216       fRowPadVector[i][j]=NULL;
 217     }
 218     fRowPadVector[i].clear();
 219   }
 220   fRowPadVector.clear();
 221   return 1;
 222 }
 223
 224
 225 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
 226 {
 227   //init slice
 228   fNClusters = 0;
 229   fMaxNClusters = nmaxpoints;
 230   fCurrentSlice = slice;
 231   fCurrentPatch = patch;
 232   fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
 233   fLastRow=AliHLTTPCTransform::GetLastRow(patch);
 234 }
 235
 236 void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt)
 237 {
 238   //set pointer to output
 239   fSpacePointData = pt;
 240 }
 241
 242 void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
 243   //set input pointer
 244   fPtr = (UChar_t*)ptr;
 245   fSize = size;
 246 }
 247
 248 void AliHLTTPCClusterFinder::ProcessDigits()
 249 {
 250   int iResult=0;
 251   bool readValue = true;
 252   Int_t newRow = 0;
 253   Int_t rowOffset = 0;
 254   UShort_t time=0,newTime=0;
 255   UInt_t pad=0,newPad=0;
 256   AliHLTTPCSignal_t charge=0;
 257
 258   fNClusters = 0;
 259
 260   // initialize block for reading packed data
 261   iResult=fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
 262   if (iResult<0) return;
 263
 264   readValue = fDigitReader->Next();
 265
 266   // Matthias 08.11.2006 the following return would cause termination without writing the
 267   // ClusterData and thus would block the component. I just want to have the commented line
 268   // here for information
 269   //if (!readValue)return;
 270
 271   pad = fDigitReader->GetPad();
 272   time = fDigitReader->GetTime();
 273   fCurrentRow = fDigitReader->GetRow();
 274
 275   if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
 276     rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
 277
 278   fCurrentRow += rowOffset;
 279
 280   UInt_t lastpad = 123456789;
 281   const UInt_t kPadArraySize=5000;
 282   const UInt_t kClusterListSize=10000;
 283   AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
 284   AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
 285   AliClusterData clusterlist[kClusterListSize]; //Clusterlist
 286
 287   AliClusterData **currentPt;  //List of pointers to the current pad
 288   AliClusterData **previousPt; //List of pointers to the previous pad
 289   currentPt = pad2;
 290   previousPt = pad1;
 291   UInt_t nprevious=0,ncurrent=0,ntotal=0;
 292
 293   /* quick implementation of baseline calculation and zero suppression
 294      open a pad object for each pad and delete it after processing.
 295      later a list of pad objects with base line history can be used
 296      The whole thing only works if we really get unprocessed raw data, if
 297      the data is already zero suppressed, there might be gaps in the time
 298      bins.
 299    */
 300   Int_t gatingGridOffset=50;
 301   AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
 302   // just to make later conversion to a list of objects easier
 303   AliHLTTPCPad* pCurrentPad=NULL;
 304   if (fSignalThreshold>=0) {
 305     pCurrentPad=&baseline;
 306     baseline.SetThreshold(fSignalThreshold);
 307   }
 308
 309   while ( readValue!=0 && iResult>=0){   // Reads through all digits in block
 310     iResult=0;
 311
 312     if(pad != lastpad){
 313       //This is a new pad
 314
 315       //Switch the lists:
 316       if(currentPt == pad2){
 317         currentPt = pad1;
 318         previousPt = pad2;
 319       }
 320       else {
 321         currentPt = pad2;
 322         previousPt = pad1;
 323       }
 324       nprevious = ncurrent;
 325       ncurrent = 0;
 326       if(pad != lastpad+1){
 327         //this happens if there is a pad with no signal.
 328         nprevious = ncurrent = 0;
 329       }
 330       lastpad = pad;
 331     }
 332
 333     Bool_t newcluster = kTRUE;
 334     UInt_t seqcharge=0,seqaverage=0,seqerror=0;
 335     AliHLTTPCSignal_t lastcharge=0;
 336     UInt_t bLastWasFalling=0;
 337     Int_t newbin=-1;
 338
 339
 340     if(fDeconvTime){
 341       redo: //This is a goto.
 342
 343       if(newbin > -1){
 344         //bin = newbin;
 345         newbin = -1;
 346       }
 347
 348       lastcharge=0;
 349       bLastWasFalling = 0;
 350     }
 351
 352     while(iResult>=0){ //Loop over time bins of current pad
 353       // read all the values for one pad at once to calculate the base line
 354       if (pCurrentPad) {
 355         if (!pCurrentPad->IsStarted()) {
 356           //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
 357           pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
 358           if ((pCurrentPad->StartEvent())>=0) {
 359             do {
 360               if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
 361               if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
 362               pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
 363               //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
 364             } while ((readValue = fDigitReader->Next())!=0);
 365           }
 366           pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
 367           if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
 368             //HLTDebug("no data available after zero suppression");
 369             pCurrentPad->StopEvent();
 370             pCurrentPad->ResetHistory();
 371             break;
 372           }
 373           time=pCurrentPad->GetCurrentPosition();
 374           if (time>pCurrentPad->GetSize()) {
 375             HLTError("invalid time bin for pad");
 376             break;
 377           }
 378         }
 379       }
 380       if (pCurrentPad) {
 381         Float_t occupancy=pCurrentPad->GetOccupancy();
 382         //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
 383         if ( occupancy < fOccupancyLimit ) {
 384           charge = pCurrentPad->GetCorrectedData();
 385         } else {
 386           charge = 0;
 387           //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
 388         }
 389       } else {
 390         charge = fDigitReader->GetSignal();
 391       }
 392       //HLTDebug("get next charge value: position %d charge %d", time, charge);
 393
 394
 395       // CHARGE DEBUG
 396       if (fDigitReader->GetRow() == 90){
 397 /////     LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90")  << "PAD=" <<  fDigitReader->GetPad() << "  TIME=" <<  fDigitReader->GetTime()
 398           //                                       << "  SIGNAL=" <<  fDigitReader->GetSignal() << ENDLOG;
 399
 400       }
 401
 402       if(time >= AliHLTTPCTransform::GetNTimeBins()){
 403         HLTWarning("Pad %d: Timebin (%d) out of range (%d)", pad, time, AliHLTTPCTransform::GetNTimeBins());
 404         iResult=-ERANGE;
 405       }
 406
 407
 408       //Get the current ADC-value
 409       if(fDeconvTime){
 410
 411         //Check if the last pixel in the sequence is smaller than this
 412         if(charge > lastcharge){
 413           if(bLastWasFalling){
 414             newbin = 1;
 415             break;
 416           }
 417         }
 418         else bLastWasFalling = 1; //last pixel was larger than this
 419         lastcharge = charge;
 420       }
 421
 422       //Sum the total charge of this sequence
 423       seqcharge += charge;
 424       seqaverage += time*charge;
 425       seqerror += time*time*charge;
 426
 427       if (pCurrentPad) {
 428
 429         if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
 430           pCurrentPad->StopEvent();
 431           pCurrentPad->ResetHistory();
 432           if(readValue) {
 433             newPad = fDigitReader->GetPad();
 434             newTime = fDigitReader->GetTime();
 435             newRow = fDigitReader->GetRow() + rowOffset;
 436           }
 437           break;
 438         }
 439
 440         newPad=pCurrentPad->GetPadNumber();
 441         newTime=pCurrentPad->GetCurrentPosition();
 442         newRow=pCurrentPad->GetRowNumber();
 443       } else {
 444       readValue = fDigitReader->Next();
 445       //Check where to stop:
 446       if(!readValue) break; //No more value
 447
 448       newPad = fDigitReader->GetPad();
 449       newTime = fDigitReader->GetTime();
 450       newRow = fDigitReader->GetRow() + rowOffset;
 451       }
 452
 453       if(newPad != pad)break; //new pad
 454       if(newTime != time+1) break; //end of sequence
 455       if(iResult<0) break;
 456
 457       // pad = newpad;    is equal
 458       time = newTime;
 459
 460     }//end loop over sequence
 461
 462     //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
 463     //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
 464     if (seqcharge<=0) {
 465       // with active zero suppression zero values are possible
 466       continue;
 467     }
 468
 469     //Calculate mean of sequence:
 470     Int_t seqmean=0;
 471     if(seqcharge)
 472       seqmean = seqaverage/seqcharge;
 473     else{
 474       LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
 475         <<"Error in data given to the cluster finder"<<ENDLOG;
 476       seqmean = 1;
 477       seqcharge = 1;
 478     }
 479
 480     //Calculate mean in pad direction:
 481     Int_t padmean = seqcharge*pad;
 482     Int_t paderror = pad*padmean;
 483
 484     //Compare with results on previous pad:
 485     for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
 486
 487       //dont merge sequences on the same pad twice
 488       if(previousPt[p]->fLastMergedPad==pad) continue;
 489
 490       Int_t difference = seqmean - previousPt[p]->fMean;
 491       if(difference < -fMatch) break;
 492
 493       if(difference <= fMatch){ //There is a match here!!
 494         AliClusterData *local = previousPt[p];
 495
 496         if(fDeconvPad){
 497           if(seqcharge > local->fLastCharge){
 498             if(local->fChargeFalling){ //The previous pad was falling
 499               break; //create a new cluster
 500             }
 501           }
 502           else local->fChargeFalling = 1;
 503           local->fLastCharge = seqcharge;
 504         }
 505
 506         //Don't create a new cluster, because we found a match
 507         newcluster = kFALSE;
 508
 509         //Update cluster on current pad with the matching one:
 510         local->fTotalCharge += seqcharge;
 511         local->fPad += padmean;
 512         local->fPad2 += paderror;
 513         local->fTime += seqaverage;
 514         local->fTime2 += seqerror;
 515         local->fMean = seqmean;
 516         local->fFlags++; //means we have more than one pad
 517         local->fLastMergedPad = pad;
 518
 519         currentPt[ncurrent] = local;
 520         ncurrent++;
 521
 522         break;
 523       } //Checking for match at previous pad
 524     } //Loop over results on previous pad.
 525
 526     if(newcluster && ncurrent<kPadArraySize){
 527       //Start a new cluster. Add it to the clusterlist, and update
 528       //the list of pointers to clusters in current pad.
 529       //current pad will be previous pad on next pad.
 530
 531       //Add to the clusterlist:
 532       AliClusterData *tmp = &clusterlist[ntotal];
 533       tmp->fTotalCharge = seqcharge;
 534       tmp->fPad = padmean;
 535       tmp->fPad2 = paderror;
 536       tmp->fTime = seqaverage;
 537       tmp->fTime2 = seqerror;
 538       tmp->fMean = seqmean;
 539       tmp->fFlags = 0;  //flags for single pad clusters
 540       tmp->fLastMergedPad = pad;
 541
 542       if(fDeconvPad){
 543         tmp->fChargeFalling = 0;
 544         tmp->fLastCharge = seqcharge;
 545       }
 546
 547       //Update list of pointers to previous pad:
 548       currentPt[ncurrent] = &clusterlist[ntotal];
 549       ntotal++;
 550       ncurrent++;
 551     }
 552
 553     if(fDeconvTime)
 554       if(newbin >= 0) goto redo;
 555
 556     // to prevent endless loop
 557     if(time >= AliHLTTPCTransform::GetNTimeBins()){
 558       HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
 559       break;
 560     }
 561
 562
 563     if(!readValue) break; //No more value
 564
 565     if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
 566       HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
 567       break;
 568     }
 569
 570     if(fCurrentRow != newRow){
 571       WriteClusters(ntotal,clusterlist);
 572
 573       lastpad = 123456789;
 574
 575       currentPt = pad2;
 576       previousPt = pad1;
 577       nprevious=0;
 578       ncurrent=0;
 579       ntotal=0;
 580
 581       fCurrentRow = newRow;
 582     }
 583
 584     pad = newPad;
 585     time = newTime;
 586
 587   } // END while(readValue)
 588
 589   WriteClusters(ntotal,clusterlist);
 590
 591   HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
 592
 593 } // ENDEND
 594
 595 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list)
 596 {
 597   //write cluster to output pointer
 598   Int_t thisrow=-1,thissector=-1;
 599   UInt_t counter = fNClusters;
 600
 601   for(int j=0; j<nclusters; j++)
 602     {
 603
 604
 605
 606       if(!list[j].fFlags) continue; //discard single pad clusters
 607       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
 608
 609       Float_t xyz[3];
 610       Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
 611       Float_t fpad2=fXYErr*fXYErr; //fixed given error
 612       Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
 613       Float_t ftime2=fZErr*fZErr;  //fixed given error
 614
 615
 616
 617       if(fUnsorted){
 618         fCurrentRow=list[j].fRow;
 619       }
 620
 621
 622       if(fCalcerr) { //calc the errors, otherwice take the fixed error
 623         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
 624         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
 625         Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
 626         sy2/=q2;
 627         if(sy2 < 0) {
 628             LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
 629               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
 630             continue;
 631         } else {
 632           if(!fRawSP){
 633             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
 634             if(sy2 != 0){
 635               fpad2*=0.108; //constants are from offline studies
 636               if(patch<2)
 637                 fpad2*=2.07;
 638             }
 639           } else fpad2=sy2; //take the width not the error
 640         }
 641         Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
 642         sz2/=q2;
 643         if(sz2 < 0){
 644           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
 645             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
 646           continue;
 647         } else {
 648           if(!fRawSP){
 649             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
 650             if(sz2 != 0) {
 651               ftime2 *= 0.169; //constants are from offline studies
 652               if(patch<2)
 653                 ftime2 *= 1.77;
 654             }
 655           } else ftime2=sz2; //take the width, not the error
 656         }
 657       }
 658       if(fStdout==kTRUE)
 659         HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
 660
 661       if(!fRawSP){
 662         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
 663         AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
 664
 665         if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
 666           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
 667         if(fNClusters >= fMaxNClusters)
 668           {
 669             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
 670               <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
 671             return;
 672           }
 673
 674         fSpacePointData[counter].fX = xyz[0];
 675         fSpacePointData[counter].fY = xyz[1];
 676         fSpacePointData[counter].fZ = xyz[2];
 677
 678       } else {
 679         fSpacePointData[counter].fX = fCurrentRow;
 680         fSpacePointData[counter].fY = fpad;
 681         fSpacePointData[counter].fZ = ftime;
 682       }
 683
 684       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
 685       fSpacePointData[counter].fPadRow = fCurrentRow;
 686       fSpacePointData[counter].fSigmaY2 = fpad2;
 687       fSpacePointData[counter].fSigmaZ2  = ftime2;
 688
 689       fSpacePointData[counter].fUsed = kFALSE;         // only used / set in AliHLTTPCDisplay
 690       fSpacePointData[counter].fTrackN = -1;           // only used / set in AliHLTTPCDisplay
 691
 692       Int_t patch=fCurrentPatch;
 693       if(patch==-1) patch=0; //never store negative patch number
 694       fSpacePointData[counter].fID = counter
 695         +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
 696
 697 #ifdef do_mc
 698       Int_t trackID[3];
 699       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
 700
 701       fSpacePointData[counter].fTrackID[0] = trackID[0];
 702       fSpacePointData[counter].fTrackID[1] = trackID[1];
 703       fSpacePointData[counter].fTrackID[2] = trackID[2];
 704
 705 #endif
 706
 707       fNClusters++;
 708       counter++;
 709     }
 710 }
 711
 712 // STILL TO FIX  ----------------------------------------------------------------------------
 713
 714 #ifdef do_mc
 715 void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
 716 {
 717   //get mc id
 718   AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
 719
 720   trackID[0]=trackID[1]=trackID[2]=-2;
 721   for(Int_t i=fFirstRow; i<=fLastRow; i++){
 722     if(rowPt->fRow < (UInt_t)fCurrentRow){
 723       AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
 724       continue;
 725     }
 726     AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
 727     for(UInt_t j=0; j<rowPt->fNDigit; j++){
 728       Int_t cpad = digPt[j].fPad;
 729       Int_t ctime = digPt[j].fTime;
 730       if(cpad != pad) continue;
 731       if(ctime != time) continue;
 732
 733       trackID[0] = digPt[j].fTrackID[0];
 734       trackID[1] = digPt[j].fTrackID[1];
 735       trackID[2] = digPt[j].fTrackID[2];
 736
 737       break;
 738     }
 739     break;
 740   }
 741 }
 742 #endif
 743
 744 //----------------------------------Methods for the new unsorted way of reading the data --------------------------------
 745
 746 void AliHLTTPCClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size)
 747 {
 748   //set input pointer
 749   fPtr = (UChar_t*)ptr;
 750   fSize = size;
 751
 752   if(!fVectorInitialized){
 753     InitializePadArray();
 754   }
 755
 756   fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
 757
 758   while(fDigitReader->NextChannel()){
 759     UInt_t row=fDigitReader->GetRow();
 760     UInt_t pad=fDigitReader->GetPad();
 761     if(row==1000||pad==1000){
 762       HLTInfo("Corrupt data, hw address exceeded the maximum amount");
 763       continue;
 764     }
 765
 766     fRowPadVector[row][pad]->ClearCandidates();
 767     while(fDigitReader->NextBunch()){
 768       if(fDigitReader->GetBunchSize()>1){//to remove single timebin values, this will have to change at some point
 769         const UInt_t *bunchData= fDigitReader->GetSignals();
 770         UInt_t time = fDigitReader->GetTime();
 771         AliHLTTPCClusters candidate;
 772         for(Int_t i=0;i<fDigitReader->GetBunchSize();i++){
 773           candidate.fTotalCharge+=bunchData[i];
 774           candidate.fTime += time*bunchData[i];
 775           candidate.fTime2 += time*time*bunchData[i];
 776           time++;
 777         }
 778         if(candidate.fTotalCharge>0){
 779           candidate.fMean=candidate.fTime/candidate.fTotalCharge;
 780           candidate.fPad=candidate.fTotalCharge*pad;
 781           candidate.fPad2=candidate.fPad*pad;
 782           candidate.fLastMergedPad=pad;
 783           candidate.fRowNumber=row;
 784         }
 785         fRowPadVector[row][pad]->AddClusterCandidate(candidate);
 786       }
 787     }
 788   }
 789 }
 790
 791 Bool_t AliHLTTPCClusterFinder::ComparePads(AliHLTTPCPad *nextPad,AliHLTTPCClusters* cluster,Int_t nextPadToRead){
 792   //Checking if we have a match on the next pad
 793   for(UInt_t candidateNumber=0;candidateNumber<nextPad->fClusterCandidates.size();candidateNumber++){
 794     AliHLTTPCClusters *candidate =&nextPad->fClusterCandidates[candidateNumber];
 795     if(cluster->fMean-candidate->fMean==1 || candidate->fMean-cluster->fMean==1 || cluster->fMean-candidate->fMean==0){
 796       cluster->fMean=candidate->fMean;
 797       cluster->fTotalCharge+=candidate->fTotalCharge;
 798       cluster->fTime += candidate->fTime;
 799       cluster->fTime2 += candidate->fTime2;
 800       cluster->fPad+=candidate->fPad;
 801       cluster->fPad2=candidate->fPad2;
 802       cluster->fLastMergedPad=candidate->fPad;
 803
 804       //setting the matched pad to used
 805       nextPad->fUsedClusterCandidates[candidateNumber]=1;
 806       nextPadToRead++;
 807       if(nextPadToRead<(Int_t)fNumberOfPadsInRow[fRowOfFirstCandidate]){
 808         nextPad=fRowPadVector[fRowOfFirstCandidate][nextPadToRead];
 809         ComparePads(nextPad,cluster,nextPadToRead);
 810       }
 811       else{
 812         return kFALSE;
 813       }
 814     }
 815     else{
 816       return kFALSE;
 817     }
 818   }
 819   return kFALSE;
 820 }
 821
 822 void AliHLTTPCClusterFinder::FindClusters()
 823 {
 824   // see header file for function documentation
 825
 826   AliHLTTPCClusters* tmpCandidate=NULL;
 827   for(UInt_t row=0;row<fNumberOfRows;row++){
 828     fRowOfFirstCandidate=row;
 829     for(UInt_t pad=0;pad<fNumberOfPadsInRow[row]-1;pad++){
 830       AliHLTTPCPad *tmpPad=fRowPadVector[row][pad];
 831       for(size_t candidate=0;candidate<tmpPad->fClusterCandidates.size();candidate++){
 832         if(tmpPad->fUsedClusterCandidates[candidate]){
 833           continue;
 834         }
 835         tmpCandidate=&tmpPad->fClusterCandidates[candidate];
 836         UInt_t tmpTotalCharge=tmpCandidate->fTotalCharge;
 837         ComparePads(fRowPadVector[row][pad+1],tmpCandidate,pad+1);
 838         if(tmpCandidate->fTotalCharge>tmpTotalCharge){
 839           //we have a cluster
 840           fClusters.push_back(*tmpCandidate);
 841         }
 842       }
 843     }
 844   }
 845
 846   HLTInfo("Found %d clusters.",fClusters.size());
 847
 848   //TODO:  Change so it stores AliHLTTPCSpacePointData directly, instead of this copying
 849
 850   AliClusterData * clusterlist = new AliClusterData[fClusters.size()]; //Clusterlist
 851   for(unsigned int i=0;i<fClusters.size();i++){
 852     clusterlist[i].fTotalCharge = fClusters[i].fTotalCharge;
 853     clusterlist[i].fPad = fClusters[i].fPad;
 854     clusterlist[i].fPad2 = fClusters[i].fPad2;
 855     clusterlist[i].fTime = fClusters[i].fTime;
 856     clusterlist[i].fTime2 = fClusters[i].fTime2;
 857     clusterlist[i].fMean = fClusters[i].fMean;
 858     clusterlist[i].fFlags = fClusters[i].fFlags;
 859     clusterlist[i].fChargeFalling = fClusters[i].fChargeFalling;
 860     clusterlist[i].fLastCharge = fClusters[i].fLastCharge;
 861     clusterlist[i].fLastMergedPad = fClusters[i].fLastMergedPad;
 862     clusterlist[i].fRow = fClusters[i].fRowNumber;
 863   }
 864   //  PrintClusters();
 865   WriteClusters(fClusters.size(),clusterlist);
 866   delete [] clusterlist;
 867   fClusters.clear();
 868 }
 869
 870 void AliHLTTPCClusterFinder::PrintClusters()
 871 {
 872   // see header file for class documentation
 873   for(size_t i=0;i<fClusters.size();i++){
 874     HLTInfo("Cluster number: %d",i);
 875     HLTInfo("Row: %d \t Pad: %d",fClusters[i].fRowNumber,fClusters[i].fFirstPad);
 876     HLTInfo("Total Charge:   %d",fClusters[i].fTotalCharge);
 877     HLTInfo("fPad:           %d",fClusters[i].fPad);
 878     HLTInfo("PadError:       %d",fClusters[i].fPad2);
 879     HLTInfo("TimeMean:       %d",fClusters[i].fTime);
 880     HLTInfo("TimeError:      %d",fClusters[i].fTime2);
 881     HLTInfo("EndOfCluster:");
 882   }
 883 }
 884
 885 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list)//This is used when using the AliHLTTPCClusters class for cluster data
 886 {
 887   //write cluster to output pointer
 888   Int_t thisrow,thissector;
 889   UInt_t counter = fNClusters;
 890
 891   for(int j=0; j<nclusters; j++)
 892     {
 893       if(!list[j].fFlags) continue; //discard single pad clusters
 894       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
 895
 896       Float_t xyz[3];
 897       Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
 898       Float_t fpad2=fXYErr*fXYErr; //fixed given error
 899       Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
 900       Float_t ftime2=fZErr*fZErr;  //fixed given error
 901
 902
 903       if(fCalcerr) { //calc the errors, otherwice take the fixed error
 904         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
 905         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
 906         Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
 907         sy2/=q2;
 908         if(sy2 < 0) {
 909             LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
 910               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
 911             continue;
 912         } else {
 913           if(!fRawSP){
 914             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
 915             if(sy2 != 0){
 916               fpad2*=0.108; //constants are from offline studies
 917               if(patch<2)
 918                 fpad2*=2.07;
 919             }
 920           } else fpad2=sy2; //take the width not the error
 921         }
 922         Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
 923         sz2/=q2;
 924         if(sz2 < 0){
 925           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
 926             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
 927           continue;
 928         } else {
 929           if(!fRawSP){
 930             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
 931             if(sz2 != 0) {
 932               ftime2 *= 0.169; //constants are from offline studies
 933               if(patch<2)
 934                 ftime2 *= 1.77;
 935             }
 936           } else ftime2=sz2; //take the width, not the error
 937         }
 938       }
 939       if(fStdout==kTRUE)
 940         HLTInfo("WriteCluster: padrow %d pad %d +- %d time +- %d charge %d",fCurrentRow, fpad, fpad2, ftime, ftime2, list[j].fTotalCharge);
 941
 942       if(!fRawSP){
 943         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
 944         AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
 945
 946         if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
 947           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
 948         if(fNClusters >= fMaxNClusters)
 949           {
 950             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
 951               <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
 952             return;
 953           }
 954
 955         fSpacePointData[counter].fX = xyz[0];
 956         fSpacePointData[counter].fY = xyz[1];
 957         fSpacePointData[counter].fZ = xyz[2];
 958
 959       } else {
 960         fSpacePointData[counter].fX = fCurrentRow;
 961         fSpacePointData[counter].fY = fpad;
 962         fSpacePointData[counter].fZ = ftime;
 963       }
 964
 965       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
 966       fSpacePointData[counter].fPadRow = fCurrentRow;
 967       fSpacePointData[counter].fSigmaY2 = fpad2;
 968       fSpacePointData[counter].fSigmaZ2  = ftime2;
 969
 970       fSpacePointData[counter].fUsed = kFALSE;         // only used / set in AliHLTTPCDisplay
 971       fSpacePointData[counter].fTrackN = -1;           // only used / set in AliHLTTPCDisplay
 972
 973       Int_t patch=fCurrentPatch;
 974       if(patch==-1) patch=0; //never store negative patch number
 975       fSpacePointData[counter].fID = counter
 976         +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
 977
 978 #ifdef do_mc
 979       Int_t trackID[3];
 980       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
 981
 982       fSpacePointData[counter].fTrackID[0] = trackID[0];
 983       fSpacePointData[counter].fTrackID[1] = trackID[1];
 984       fSpacePointData[counter].fTrackID[2] = trackID[2];
 985
 986 #endif
 987
 988       fNClusters++;
 989       counter++;
 990     }
 991 }