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