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