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 "AliHLTTPCPadArray.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   fNClusters(0),
 140   fMaxNClusters(0),
 141   fXYErr(0.2),
 142   fZErr(0.3),
 143   fOccupancyLimit(1.0),
 144   fPadArray(NULL)
 145 {
 146   //constructor
 147 }
 148
 149 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder(const AliHLTTPCClusterFinder& src)
 150   :
 151   fSpacePointData(NULL),
 152   fDigitReader(src.fDigitReader),
 153   fPtr(NULL),
 154   fSize(src.fSize),
 155   fDeconvTime(src.fDeconvTime),
 156   fDeconvPad(src.fDeconvPad),
 157   fStdout(src.fStdout),
 158   fCalcerr(src.fCalcerr),
 159   fRawSP(src.fRawSP),
 160   fFirstRow(src.fFirstRow),
 161   fLastRow(src.fLastRow),
 162   fCurrentRow(src.fCurrentRow),
 163   fCurrentSlice(src.fCurrentSlice),
 164   fCurrentPatch(src.fCurrentPatch),
 165   fMatch(src.fMatch),
 166   fThreshold(src.fThreshold),
 167   fSignalThreshold(src.fSignalThreshold),
 168   fNClusters(src.fNClusters),
 169   fMaxNClusters(src.fMaxNClusters),
 170   fXYErr(src.fXYErr),
 171   fZErr(src.fZErr),
 172   fOccupancyLimit(src.fOccupancyLimit),
 173   fPadArray(NULL)
 174 {
 175 }
 176
 177 AliHLTTPCClusterFinder& AliHLTTPCClusterFinder::operator=(const AliHLTTPCClusterFinder& src)
 178 {
 179   fMatch=src.fMatch;
 180   fThreshold=src.fThreshold;
 181   fSignalThreshold=src.fSignalThreshold;
 182   fOccupancyLimit=src.fOccupancyLimit;
 183   fXYErr=src.fXYErr;
 184   fZErr=src.fZErr;
 185   fDeconvPad=src.fDeconvPad;
 186   fDeconvTime=src.fDeconvTime;
 187   fStdout=src.fStdout;
 188   fCalcerr=src.fCalcerr;
 189   fRawSP=src.fRawSP;
 190   fFirstRow=src.fFirstRow;
 191   fLastRow=src.fLastRow;
 192   fDigitReader=src.fDigitReader;
 193   return (*this);
 194 }
 195
 196 AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder()
 197 {
 198   //destructor
 199 }
 200
 201 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
 202 {
 203   //init slice
 204   fNClusters = 0;
 205   fMaxNClusters = nmaxpoints;
 206   fCurrentSlice = slice;
 207   fCurrentPatch = patch;
 208   fFirstRow = firstrow;
 209   fLastRow = lastrow;
 210 }
 211
 212 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
 213 {
 214   //init slice
 215   fNClusters = 0;
 216   fMaxNClusters = nmaxpoints;
 217   fCurrentSlice = slice;
 218   fCurrentPatch = patch;
 219   fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
 220   fLastRow=AliHLTTPCTransform::GetLastRow(patch);
 221 }
 222
 223 void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt)
 224 {
 225   //set pointer to output
 226   fSpacePointData = pt;
 227 }
 228
 229 void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
 230   //set input pointer
 231   fPtr = (UChar_t*)ptr;
 232   fSize = size;
 233 }
 234
 235 void AliHLTTPCClusterFinder::ProcessDigits()
 236 {
 237   int iResult=0;
 238   bool readValue = true;
 239   Int_t newRow = 0;
 240   Int_t rowOffset = 0;
 241   UShort_t time=0,newTime=0;
 242   UInt_t pad=0,newPad=0;
 243   AliHLTTPCSignal_t charge=0;
 244
 245   fNClusters = 0;
 246
 247   // initialize block for reading packed data
 248   fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
 249   readValue = fDigitReader->Next();
 250
 251   // Matthias 08.11.2006 the following return would cause termination without writing the
 252   // ClusterData and thus would block the component. I just want to have the commented line
 253   // here for information
 254   //if (!readValue)return;
 255
 256   pad = fDigitReader->GetPad();
 257   time = fDigitReader->GetTime();
 258   fCurrentRow = fDigitReader->GetRow();
 259
 260   if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
 261     rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
 262
 263   fCurrentRow += rowOffset;
 264
 265   UInt_t lastpad = 123456789;
 266   const UInt_t kPadArraySize=5000;
 267   const UInt_t kClusterListSize=10000;
 268   AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
 269   AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
 270   AliClusterData clusterlist[kClusterListSize]; //Clusterlist
 271
 272   AliClusterData **currentPt;  //List of pointers to the current pad
 273   AliClusterData **previousPt; //List of pointers to the previous pad
 274   currentPt = pad2;
 275   previousPt = pad1;
 276   UInt_t nprevious=0,ncurrent=0,ntotal=0;
 277
 278   /* quick implementation of baseline calculation and zero suppression
 279      open a pad object for each pad and delete it after processing.
 280      later a list of pad objects with base line history can be used
 281      The whole thing only works if we really get unprocessed raw data, if
 282      the data is already zero suppressed, there might be gaps in the time
 283      bins.
 284    */
 285   Int_t gatingGridOffset=50;
 286   AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
 287   // just to make later conversion to a list of objects easier
 288   AliHLTTPCPad* pCurrentPad=NULL;
 289   if (fSignalThreshold>=0) {
 290     pCurrentPad=&baseline;
 291     baseline.SetThreshold(fSignalThreshold);
 292   }
 293
 294   while ( readValue!=0 && iResult>=0){   // Reads through all digits in block
 295     iResult=0;
 296
 297     if(pad != lastpad){
 298       //This is a new pad
 299
 300       //Switch the lists:
 301       if(currentPt == pad2){
 302         currentPt = pad1;
 303         previousPt = pad2;
 304       }
 305       else {
 306         currentPt = pad2;
 307         previousPt = pad1;
 308       }
 309       nprevious = ncurrent;
 310       ncurrent = 0;
 311       if(pad != lastpad+1){
 312         //this happens if there is a pad with no signal.
 313         nprevious = ncurrent = 0;
 314       }
 315       lastpad = pad;
 316     }
 317
 318     Bool_t newcluster = kTRUE;
 319     UInt_t seqcharge=0,seqaverage=0,seqerror=0;
 320     AliHLTTPCSignal_t lastcharge=0;
 321     UInt_t bLastWasFalling=0;
 322     Int_t newbin=-1;
 323
 324
 325     if(fDeconvTime){
 326       redo: //This is a goto.
 327
 328       if(newbin > -1){
 329         //bin = newbin;
 330         newbin = -1;
 331       }
 332
 333       lastcharge=0;
 334       bLastWasFalling = 0;
 335     }
 336
 337     while(iResult>=0){ //Loop over time bins of current pad
 338       // read all the values for one pad at once to calculate the base line
 339       if (pCurrentPad) {
 340         if (!pCurrentPad->IsStarted()) {
 341           //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
 342           pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
 343           if ((pCurrentPad->StartEvent())>=0) {
 344             do {
 345               if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
 346               if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
 347               pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
 348               //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
 349             } while ((readValue = fDigitReader->Next())!=0);
 350           }
 351           pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
 352           if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
 353             //HLTDebug("no data available after zero suppression");
 354             pCurrentPad->StopEvent();
 355             pCurrentPad->ResetHistory();
 356             break;
 357           }
 358           time=pCurrentPad->GetCurrentPosition();
 359           if (time>pCurrentPad->GetSize()) {
 360             HLTError("invalid time bin for pad");
 361             break;
 362           }
 363         }
 364       }
 365
 366       if (pCurrentPad) {
 367         Float_t occupancy=pCurrentPad->GetOccupancy();
 368         //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
 369         if ( occupancy < fOccupancyLimit ) {
 370           charge = pCurrentPad->GetCorrectedData();
 371         } else {
 372           charge = 0;
 373           //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
 374         }
 375       } else {
 376         charge = fDigitReader->GetSignal();
 377       }
 378       //HLTDebug("get next charge value: position %d charge %d", time, charge);
 379
 380
 381       // CHARGE DEBUG
 382       if (fDigitReader->GetRow() == 90){
 383 /////     LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90")  << "PAD=" <<  fDigitReader->GetPad() << "  TIME=" <<  fDigitReader->GetTime()
 384           //                                       << "  SIGNAL=" <<  fDigitReader->GetSignal() << ENDLOG;
 385
 386       }
 387
 388       if(time >= AliHLTTPCTransform::GetNTimeBins()){
 389         HLTWarning("Pad %d: Timebin (%d) out of range (%d)", pad, time, AliHLTTPCTransform::GetNTimeBins());
 390         iResult=-ERANGE;
 391       }
 392
 393
 394       //Get the current ADC-value
 395       if(fDeconvTime){
 396
 397         //Check if the last pixel in the sequence is smaller than this
 398         if(charge > lastcharge){
 399           if(bLastWasFalling){
 400             newbin = 1;
 401             break;
 402           }
 403         }
 404         else bLastWasFalling = 1; //last pixel was larger than this
 405         lastcharge = charge;
 406       }
 407
 408       //Sum the total charge of this sequence
 409       seqcharge += charge;
 410       seqaverage += time*charge;
 411       seqerror += time*time*charge;
 412
 413       if (pCurrentPad) {
 414
 415         if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
 416           pCurrentPad->StopEvent();
 417           pCurrentPad->ResetHistory();
 418           if(readValue) {
 419             newPad = fDigitReader->GetPad();
 420             newTime = fDigitReader->GetTime();
 421             newRow = fDigitReader->GetRow() + rowOffset;
 422           }
 423           break;
 424         }
 425
 426         newPad=pCurrentPad->GetPadNumber();
 427         newTime=pCurrentPad->GetCurrentPosition();
 428         newRow=pCurrentPad->GetRowNumber();
 429       } else {
 430       readValue = fDigitReader->Next();
 431       //Check where to stop:
 432       if(!readValue) break; //No more value
 433
 434       newPad = fDigitReader->GetPad();
 435       newTime = fDigitReader->GetTime();
 436       newRow = fDigitReader->GetRow() + rowOffset;
 437       }
 438
 439       if(newPad != pad)break; //new pad
 440       if(newTime != time+1) break; //end of sequence
 441       if(iResult<0) break;
 442
 443       // pad = newpad;    is equal
 444       time = newTime;
 445
 446     }//end loop over sequence
 447
 448     //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
 449     //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
 450     if (seqcharge<=0) {
 451       // with active zero suppression zero values are possible
 452       continue;
 453     }
 454
 455     //Calculate mean of sequence:
 456     Int_t seqmean=0;
 457     if(seqcharge)
 458       seqmean = seqaverage/seqcharge;
 459     else{
 460       LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
 461         <<"Error in data given to the cluster finder"<<ENDLOG;
 462       seqmean = 1;
 463       seqcharge = 1;
 464     }
 465
 466     //Calculate mean in pad direction:
 467     Int_t padmean = seqcharge*pad;
 468     Int_t paderror = pad*padmean;
 469
 470     //Compare with results on previous pad:
 471     for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
 472
 473       //dont merge sequences on the same pad twice
 474       if(previousPt[p]->fLastMergedPad==pad) continue;
 475
 476       Int_t difference = seqmean - previousPt[p]->fMean;
 477       if(difference < -fMatch) break;
 478
 479       if(difference <= fMatch){ //There is a match here!!
 480         AliClusterData *local = previousPt[p];
 481
 482         if(fDeconvPad){
 483           if(seqcharge > local->fLastCharge){
 484             if(local->fChargeFalling){ //The previous pad was falling
 485               break; //create a new cluster
 486             }
 487           }
 488           else local->fChargeFalling = 1;
 489           local->fLastCharge = seqcharge;
 490         }
 491
 492         //Don't create a new cluster, because we found a match
 493         newcluster = kFALSE;
 494
 495         //Update cluster on current pad with the matching one:
 496         local->fTotalCharge += seqcharge;
 497         local->fPad += padmean;
 498         local->fPad2 += paderror;
 499         local->fTime += seqaverage;
 500         local->fTime2 += seqerror;
 501         local->fMean = seqmean;
 502         local->fFlags++; //means we have more than one pad
 503         local->fLastMergedPad = pad;
 504
 505         currentPt[ncurrent] = local;
 506         ncurrent++;
 507
 508         break;
 509       } //Checking for match at previous pad
 510     } //Loop over results on previous pad.
 511
 512     if(newcluster && ncurrent<kPadArraySize){
 513       //Start a new cluster. Add it to the clusterlist, and update
 514       //the list of pointers to clusters in current pad.
 515       //current pad will be previous pad on next pad.
 516
 517       //Add to the clusterlist:
 518       AliClusterData *tmp = &clusterlist[ntotal];
 519       tmp->fTotalCharge = seqcharge;
 520       tmp->fPad = padmean;
 521       tmp->fPad2 = paderror;
 522       tmp->fTime = seqaverage;
 523       tmp->fTime2 = seqerror;
 524       tmp->fMean = seqmean;
 525       tmp->fFlags = 0;  //flags for single pad clusters
 526       tmp->fLastMergedPad = pad;
 527
 528       if(fDeconvPad){
 529         tmp->fChargeFalling = 0;
 530         tmp->fLastCharge = seqcharge;
 531       }
 532
 533       //Update list of pointers to previous pad:
 534       currentPt[ncurrent] = &clusterlist[ntotal];
 535       ntotal++;
 536       ncurrent++;
 537     }
 538
 539     if(fDeconvTime)
 540       if(newbin >= 0) goto redo;
 541
 542     // to prevent endless loop
 543     if(time >= AliHLTTPCTransform::GetNTimeBins()){
 544       HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
 545       break;
 546     }
 547
 548
 549     if(!readValue) break; //No more value
 550
 551     if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
 552       HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
 553       break;
 554     }
 555
 556     if(fCurrentRow != newRow){
 557       WriteClusters(ntotal,clusterlist);
 558
 559       lastpad = 123456789;
 560
 561       currentPt = pad2;
 562       previousPt = pad1;
 563       nprevious=0;
 564       ncurrent=0;
 565       ntotal=0;
 566
 567       fCurrentRow = newRow;
 568     }
 569
 570     pad = newPad;
 571     time = newTime;
 572
 573   } // END while(readValue)
 574
 575   WriteClusters(ntotal,clusterlist);
 576
 577   HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
 578
 579 } // ENDEND
 580
 581 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list)
 582 {
 583   //write cluster to output pointer
 584   Int_t thisrow,thissector;
 585   UInt_t counter = fNClusters;
 586
 587   for(int j=0; j<nclusters; j++)
 588     {
 589
 590
 591
 592       if(!list[j].fFlags) continue; //discard single pad clusters
 593       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
 594
 595       Float_t xyz[3];
 596       Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
 597       Float_t fpad2=fXYErr*fXYErr; //fixed given error
 598       Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
 599       Float_t ftime2=fZErr*fZErr;  //fixed given error
 600
 601
 602 #if UNSORTED
 603       fCurrentRow=list[j].fRow;
 604 #endif
 605
 606       if(fCalcerr) { //calc the errors, otherwice take the fixed error
 607         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
 608         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
 609         Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
 610         sy2/=q2;
 611         if(sy2 < 0) {
 612             LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
 613               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
 614             continue;
 615         } else {
 616           if(!fRawSP){
 617             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
 618             if(sy2 != 0){
 619               fpad2*=0.108; //constants are from offline studies
 620               if(patch<2)
 621                 fpad2*=2.07;
 622             }
 623           } else fpad2=sy2; //take the width not the error
 624         }
 625         Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
 626         sz2/=q2;
 627         if(sz2 < 0){
 628           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
 629             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
 630           continue;
 631         } else {
 632           if(!fRawSP){
 633             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
 634             if(sz2 != 0) {
 635               ftime2 *= 0.169; //constants are from offline studies
 636               if(patch<2)
 637                 ftime2 *= 1.77;
 638             }
 639           } else ftime2=sz2; //take the width, not the error
 640         }
 641       }
 642       if(fStdout==kTRUE)
 643         cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
 644
 645       if(!fRawSP){
 646         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
 647         AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
 648
 649         if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
 650           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
 651         if(fNClusters >= fMaxNClusters)
 652           {
 653             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
 654               <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
 655             return;
 656           }
 657
 658         fSpacePointData[counter].fX = xyz[0];
 659         fSpacePointData[counter].fY = xyz[1];
 660         fSpacePointData[counter].fZ = xyz[2];
 661
 662       } else {
 663         fSpacePointData[counter].fX = fCurrentRow;
 664         fSpacePointData[counter].fY = fpad;
 665         fSpacePointData[counter].fZ = ftime;
 666       }
 667
 668       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
 669       fSpacePointData[counter].fPadRow = fCurrentRow;
 670       fSpacePointData[counter].fSigmaY2 = fpad2;
 671       fSpacePointData[counter].fSigmaZ2  = ftime2;
 672
 673       fSpacePointData[counter].fUsed = kFALSE;         // only used / set in AliHLTTPCDisplay
 674       fSpacePointData[counter].fTrackN = -1;           // only used / set in AliHLTTPCDisplay
 675
 676       Int_t patch=fCurrentPatch;
 677       if(patch==-1) patch=0; //never store negative patch number
 678       fSpacePointData[counter].fID = counter
 679         +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
 680
 681 #ifdef do_mc
 682       Int_t trackID[3];
 683       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
 684
 685       fSpacePointData[counter].fTrackID[0] = trackID[0];
 686       fSpacePointData[counter].fTrackID[1] = trackID[1];
 687       fSpacePointData[counter].fTrackID[2] = trackID[2];
 688
 689       //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
 690 #endif
 691
 692       fNClusters++;
 693       counter++;
 694     }
 695 }
 696
 697 // STILL TO FIX  ----------------------------------------------------------------------------
 698
 699 #ifdef do_mc
 700 void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
 701 {
 702   //get mc id
 703   AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
 704
 705   trackID[0]=trackID[1]=trackID[2]=-2;
 706   //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
 707   for(Int_t i=fFirstRow; i<=fLastRow; i++){
 708     if(rowPt->fRow < (UInt_t)fCurrentRow){
 709       AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
 710       continue;
 711     }
 712     AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
 713     for(UInt_t j=0; j<rowPt->fNDigit; j++){
 714       Int_t cpad = digPt[j].fPad;
 715       Int_t ctime = digPt[j].fTime;
 716       if(cpad != pad) continue;
 717       if(ctime != time) continue;
 718
 719       trackID[0] = digPt[j].fTrackID[0];
 720       trackID[1] = digPt[j].fTrackID[1];
 721       trackID[2] = digPt[j].fTrackID[2];
 722
 723       //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
 724       break;
 725     }
 726     break;
 727   }
 728 }
 729 #endif
 730
 731 //----------------------------------Methods for the new unsorted way of reading the data --------------------------------
 732
 733 void AliHLTTPCClusterFinder::SetPadArray(AliHLTTPCPadArray * padArray){
 734   fPadArray=padArray;
 735 }
 736
 737 void AliHLTTPCClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size){
 738   //set input pointer
 739   fPtr = (UChar_t*)ptr;
 740   fSize = size;
 741
 742   fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
 743
 744   fPadArray->SetDigitReader(fDigitReader);
 745   fPadArray->ReadData();
 746 }
 747 void AliHLTTPCClusterFinder::FindClusters(){
 748   fPadArray->FindClusterCandidates();
 749   fPadArray->FindClusters(fMatch);
 750
 751   AliClusterData * clusterlist = new AliClusterData[fPadArray->fClusters.size()]; //Clusterlist
 752   for(int i=0;i<fPadArray->fClusters.size();i++){
 753     clusterlist[i].fTotalCharge = fPadArray->fClusters[i].fTotalCharge;
 754     clusterlist[i].fPad = fPadArray->fClusters[i].fPad;
 755     clusterlist[i].fPad2 = fPadArray->fClusters[i].fPad2;
 756     clusterlist[i].fTime = fPadArray->fClusters[i].fTime;
 757     clusterlist[i].fTime2 = fPadArray->fClusters[i].fTime2;
 758     clusterlist[i].fMean = fPadArray->fClusters[i].fMean;
 759     clusterlist[i].fFlags = fPadArray->fClusters[i].fFlags;
 760     clusterlist[i].fChargeFalling = fPadArray->fClusters[i].fChargeFalling;
 761     clusterlist[i].fLastCharge = fPadArray->fClusters[i].fLastCharge;
 762     clusterlist[i].fLastMergedPad = fPadArray->fClusters[i].fLastMergedPad;
 763     clusterlist[i].fRow = fPadArray->fClusters[i].fRowNumber;
 764   }
 765   WriteClusters(fPadArray->fClusters.size(),clusterlist);
 766   delete [] clusterlist;
 767   fPadArray->DataToDefault();
 768 }
 769 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list)//This is used when using the AliHLTTPCClusters class for cluster data
 770 {
 771   //write cluster to output pointer
 772   Int_t thisrow,thissector;
 773   UInt_t counter = fNClusters;
 774
 775   for(int j=0; j<nclusters; j++)
 776     {
 777       if(!list[j].fFlags) continue; //discard single pad clusters
 778       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
 779
 780       Float_t xyz[3];
 781       Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
 782       Float_t fpad2=fXYErr*fXYErr; //fixed given error
 783       Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
 784       Float_t ftime2=fZErr*fZErr;  //fixed given error
 785
 786
 787       if(fCalcerr) { //calc the errors, otherwice take the fixed error
 788         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
 789         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
 790         Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
 791         sy2/=q2;
 792         if(sy2 < 0) {
 793             LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
 794               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
 795             continue;
 796         } else {
 797           if(!fRawSP){
 798             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
 799             if(sy2 != 0){
 800               fpad2*=0.108; //constants are from offline studies
 801               if(patch<2)
 802                 fpad2*=2.07;
 803             }
 804           } else fpad2=sy2; //take the width not the error
 805         }
 806         Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
 807         sz2/=q2;
 808         if(sz2 < 0){
 809           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
 810             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
 811           continue;
 812         } else {
 813           if(!fRawSP){
 814             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
 815             if(sz2 != 0) {
 816               ftime2 *= 0.169; //constants are from offline studies
 817               if(patch<2)
 818                 ftime2 *= 1.77;
 819             }
 820           } else ftime2=sz2; //take the width, not the error
 821         }
 822       }
 823       if(fStdout==kTRUE)
 824         cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
 825
 826       if(!fRawSP){
 827         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
 828         AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
 829
 830         if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
 831           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
 832         if(fNClusters >= fMaxNClusters)
 833           {
 834             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
 835               <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
 836             return;
 837           }
 838
 839         fSpacePointData[counter].fX = xyz[0];
 840         fSpacePointData[counter].fY = xyz[1];
 841         fSpacePointData[counter].fZ = xyz[2];
 842
 843       } else {
 844         fSpacePointData[counter].fX = fCurrentRow;
 845         fSpacePointData[counter].fY = fpad;
 846         fSpacePointData[counter].fZ = ftime;
 847       }
 848
 849       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
 850       fSpacePointData[counter].fPadRow = fCurrentRow;
 851       fSpacePointData[counter].fSigmaY2 = fpad2;
 852       fSpacePointData[counter].fSigmaZ2  = ftime2;
 853
 854       fSpacePointData[counter].fUsed = kFALSE;         // only used / set in AliHLTTPCDisplay
 855       fSpacePointData[counter].fTrackN = -1;           // only used / set in AliHLTTPCDisplay
 856
 857       Int_t patch=fCurrentPatch;
 858       if(patch==-1) patch=0; //never store negative patch number
 859       fSpacePointData[counter].fID = counter
 860         +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
 861
 862 #ifdef do_mc
 863       Int_t trackID[3];
 864       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
 865
 866       fSpacePointData[counter].fTrackID[0] = trackID[0];
 867       fSpacePointData[counter].fTrackID[1] = trackID[1];
 868       fSpacePointData[counter].fTrackID[2] = trackID[2];
 869
 870       //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
 871 #endif
 872
 873       fNClusters++;
 874       counter++;
 875     }
 876 }