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