HLT/TPCLib/AliHLTTPCClusterFinder.cxx

   1 // @(#) $Id$
   2 // Original: AliHLTClustFinderNew.cxx,v 1.29 2005/06/14 10:55:21 cvetan Exp
   3
   4 /**************************************************************************
   5  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   6  *                                                                        *
   7  * Authors: Anders Vestbo <mailto:vestbo@fi.uib.no>,                      *
   8  *          Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>    *
   9  *          Jochen Thaeder <mailto:thaeder@kip.uni-heidelberg.de>         *
  10  *          for The ALICE Off-line 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 <mailto:vestbo@fi.uib.no>,
  23             Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
  24             Jochen Thaeder <mailto:thaeder@kip.uni-heidelberg.de>
  25     @date
  26     @brief  Cluster Finder for the TPC
  27 */
  28
  29 #include "AliHLTTPCDigitReader.h"
  30 #include "AliHLTTPCRootTypes.h"
  31 #include "AliHLTTPCLogging.h"
  32 #include "AliHLTTPCClusterFinder.h"
  33 #include "AliHLTTPCDigitData.h"
  34 #include "AliHLTTPCTransform.h"
  35 #include "AliHLTTPCSpacePointData.h"
  36 #include "AliHLTTPCMemHandler.h"
  37 #include "AliHLTTPCPad.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 // The cluster finder is initialized with the Init function,
  49 // providing the slice and patch information to work on.
  50 //
  51 // The input is a provided by the AliHLTTPCDigitReader class,
  52 // using the init() funktion, and the next() funktion in order
  53 // to get the next bin. Either packed or unpacked data can be
  54 // processed, dependent if one uses AliHLTTPCDigitReaderPacked
  55 // class or AliHLTTPCDigitReaderUnpacked class in the
  56 // Clusterfinder Component.
  57 // The resulting space points will be in the
  58 // array given by the SetOutputArray function.
  59 //
  60 // There are several setters which control the behaviour:
  61 //
  62 // - SetXYError(Float_t):   set fixed error in XY direction
  63 // - SetZError(Float_t):    set fixed error in Z  direction
  64 //                            (used if errors are not calculated)
  65 // - SetDeconv(Bool_t):     switch on/off deconvolution
  66 // - SetThreshold(UInt_t):  set charge threshold for cluster
  67 // - SetMatchWidth(UInt_t): set the match distance in
  68 //                            time for sequences to be merged
  69 // - SetSTDOutput(Bool_t):  switch on/off output about found clusters
  70 // - SetCalcErr(Bool_t):    switch on/off calculation of
  71 //                          space point errors (or widths in raw system)
  72 // - SetRawSP(Bool_t):      switch on/off convertion to raw system
  73 //
  74 //
  75 // Example Usage:
  76 //
  77 // AliHLTTPCFileHandler *file = new AliHLTTPCFileHandler();
  78 // file->SetAliInput(digitfile); //give some input file
  79 // for(int slice=0; slice<=35; slice++){
  80 //   for(int patch=0; pat<6; pat++){
  81 //     file->Init(slice,patch);
  82 //     UInt_t ndigits=0;
  83 //     UInt_t maxclusters=100000;
  84 //     UInt_t pointsize = maxclusters*sizeof(AliHLTTPCSpacePointData);
  85 //     AliHLTTPCSpacePointData *points = (AliHLTTPCSpacePointData*)memory->Allocate(pointsize);
  86 //     AliHLTTPCDigitRowData *digits = (AliHLTTPCDigitRowData*)file->AliAltroDigits2Memory(ndigits,event);
  87 //     AliHLTTPCClusterFinder *cf = new AliHLTTPCClusterFinder();
  88 //     cf->SetMatchWidth(2);
  89 //     cf->InitSlice( slice, patch, row[0], row[1], maxPoints );
  90 //     cf->SetSTDOutput(kTRUE);    //Some output to standard IO
  91 //     cf->SetRawSP(kFALSE);       //Convert space points to local system
  92 //     cf->SetThreshold(5);        //Threshold of cluster charge
  93 //     cf->SetDeconv(kTRUE);       //Deconv in pad and time direction
  94 //     cf->SetCalcErr(kTRUE);      //Calculate the errors of the spacepoints
  95 //     cf->SetOutputArray(points); //Move the spacepoints to the array
  96 //     cf->Read(iter->fPtr, iter->fSize ); //give the data to the cf
  97 //     cf->ProcessDigits();        //process the rows given by init
  98 //     Int_t npoints = cf->GetNumberOfClusters();
  99 //     AliHLTTPCMemHandler *out= new AliHLTTPCMemHandler();
 100 //     out->SetBinaryOutput(fname);
 101 //     out->Memory2Binary(npoints,points); //store the spacepoints
 102 //     out->CloseBinaryOutput();
 103 //     delete out;
 104 //     file->free();
 105 //     delete cf;
 106 //   }
 107 // }
 108 */
 109
 110 ClassImp(AliHLTTPCClusterFinder)
 111
 112 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
 113   :
 114   fSpacePointData(NULL),
 115   fDigitReader(NULL),
 116   fPtr(NULL),
 117   fSize(0),
 118   fDeconvTime(kTRUE),
 119   fDeconvPad(kTRUE),
 120   fStdout(kFALSE),
 121   fCalcerr(kTRUE),
 122   fRawSP(kFALSE),
 123   fFirstRow(0),
 124   fLastRow(0),
 125   fCurrentRow(0),
 126   fCurrentSlice(0),
 127   fCurrentPatch(0),
 128   fMatch(1),
 129   fThreshold(10),
 130   fSignalThreshold(-1),
 131   fNClusters(0),
 132   fMaxNClusters(0),
 133   fXYErr(0.2),
 134   fZErr(0.3),
 135   fOccupancyLimit(1.0)
 136 {
 137   //constructor
 138 }
 139
 140 AliHLTTPCClusterFinder::AliHLTTPCClusterFinder(const AliHLTTPCClusterFinder& src)
 141   :
 142   fSpacePointData(NULL),
 143   fDigitReader(src.fDigitReader),
 144   fPtr(NULL),
 145   fSize(src.fSize),
 146   fDeconvTime(src.fDeconvTime),
 147   fDeconvPad(src.fDeconvPad),
 148   fStdout(src.fStdout),
 149   fCalcerr(src.fCalcerr),
 150   fRawSP(src.fRawSP),
 151   fFirstRow(src.fFirstRow),
 152   fLastRow(src.fLastRow),
 153   fCurrentRow(src.fCurrentRow),
 154   fCurrentSlice(src.fCurrentSlice),
 155   fCurrentPatch(src.fCurrentPatch),
 156   fMatch(src.fMatch),
 157   fThreshold(src.fThreshold),
 158   fSignalThreshold(src.fSignalThreshold),
 159   fNClusters(src.fNClusters),
 160   fMaxNClusters(src.fMaxNClusters),
 161   fXYErr(src.fXYErr),
 162   fZErr(src.fZErr),
 163   fOccupancyLimit(src.fOccupancyLimit)
 164 {
 165 }
 166
 167 AliHLTTPCClusterFinder& AliHLTTPCClusterFinder::operator=(const AliHLTTPCClusterFinder& src)
 168 {
 169   fMatch=src.fMatch;
 170   fThreshold=src.fThreshold;
 171   fSignalThreshold=src.fSignalThreshold;
 172   fOccupancyLimit=src.fOccupancyLimit;
 173   fXYErr=src.fXYErr;
 174   fZErr=src.fZErr;
 175   fDeconvPad=src.fDeconvPad;
 176   fDeconvTime=src.fDeconvTime;
 177   fStdout=src.fStdout;
 178   fCalcerr=src.fCalcerr;
 179   fRawSP=src.fRawSP;
 180   fFirstRow=src.fFirstRow;
 181   fLastRow=src.fLastRow;
 182   fDigitReader=src.fDigitReader;
 183   return (*this);
 184 }
 185
 186 AliHLTTPCClusterFinder::~AliHLTTPCClusterFinder()
 187 {
 188   //destructor
 189 }
 190
 191 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
 192 {
 193   //init slice
 194   fNClusters = 0;
 195   fMaxNClusters = nmaxpoints;
 196   fCurrentSlice = slice;
 197   fCurrentPatch = patch;
 198   fFirstRow = firstrow;
 199   fLastRow = lastrow;
 200 }
 201
 202 void AliHLTTPCClusterFinder::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
 203 {
 204   //init slice
 205   fNClusters = 0;
 206   fMaxNClusters = nmaxpoints;
 207   fCurrentSlice = slice;
 208   fCurrentPatch = patch;
 209   fFirstRow=AliHLTTPCTransform::GetFirstRow(patch);
 210   fLastRow=AliHLTTPCTransform::GetLastRow(patch);
 211 }
 212
 213 void AliHLTTPCClusterFinder::SetOutputArray(AliHLTTPCSpacePointData *pt)
 214 {
 215   //set pointer to output
 216   fSpacePointData = pt;
 217 }
 218
 219 void AliHLTTPCClusterFinder::Read(void* ptr,unsigned long size){
 220   //set input pointer
 221   fPtr = (UChar_t*)ptr;
 222   fSize = size;
 223 }
 224
 225 void AliHLTTPCClusterFinder::ProcessDigits()
 226 {
 227   bool readValue = true;
 228   Int_t newRow = 0;
 229   Int_t rowOffset = 0;
 230   UShort_t time=0,newTime=0;
 231   UInt_t pad=0,newPad=0;
 232   AliHLTTPCSignal_t charge=0;
 233
 234   fNClusters = 0;
 235
 236   // initialize block for reading packed data
 237   fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
 238   readValue = fDigitReader->Next();
 239
 240   // Matthias 08.11.2006 the following return would cause termination without writing the
 241   // ClusterData and thus would block the component. I just want to have the commented line
 242   // here for information
 243   //if (!readValue)return;
 244
 245   pad = fDigitReader->GetPad();
 246   time = fDigitReader->GetTime();
 247   fCurrentRow = fDigitReader->GetRow();
 248
 249   if ( fCurrentPatch >= 2 ) // Outer sector, patches 2, 3, 4, 5
 250     rowOffset = AliHLTTPCTransform::GetFirstRow( 2 );
 251
 252   fCurrentRow += rowOffset;
 253
 254   UInt_t lastpad = 123456789;
 255   const UInt_t kPadArraySize=5000;
 256   const UInt_t kClusterListSize=10000;
 257   AliClusterData *pad1[kPadArraySize]; //2 lists for internal memory=2pads
 258   AliClusterData *pad2[kPadArraySize]; //2 lists for internal memory=2pads
 259   AliClusterData clusterlist[kClusterListSize]; //Clusterlist
 260
 261   AliClusterData **currentPt;  //List of pointers to the current pad
 262   AliClusterData **previousPt; //List of pointers to the previous pad
 263   currentPt = pad2;
 264   previousPt = pad1;
 265   UInt_t nprevious=0,ncurrent=0,ntotal=0;
 266
 267   /* quick implementation of baseline calculation and zero suppression
 268      open a pad object for each pad and delete it after processing.
 269      later a list of pad objects with base line history can be used
 270      The whole thing only works if we really get unprocessed raw data, if
 271      the data is already zero suppressed, there might be gaps in the time
 272      bins.
 273    */
 274   Int_t gatingGridOffset=50;
 275   AliHLTTPCPad baseline(gatingGridOffset, AliHLTTPCTransform::GetNTimeBins());
 276   // just to make later conversion to a list of objects easier
 277   AliHLTTPCPad* pCurrentPad=NULL;
 278   if (fSignalThreshold>=0) {
 279     pCurrentPad=&baseline;
 280     baseline.SetThreshold(fSignalThreshold);
 281   }
 282
 283   while ( readValue ){   // Reads through all digits in block
 284
 285     if(pad != lastpad){
 286       //This is a new pad
 287
 288       //Switch the lists:
 289       if(currentPt == pad2){
 290         currentPt = pad1;
 291         previousPt = pad2;
 292       }
 293       else {
 294         currentPt = pad2;
 295         previousPt = pad1;
 296       }
 297       nprevious = ncurrent;
 298       ncurrent = 0;
 299       if(pad != lastpad+1){
 300         //this happens if there is a pad with no signal.
 301         nprevious = ncurrent = 0;
 302       }
 303       lastpad = pad;
 304     }
 305
 306     Bool_t newcluster = kTRUE;
 307     UInt_t seqcharge=0,seqaverage=0,seqerror=0;
 308     AliHLTTPCSignal_t lastcharge=0;
 309     UInt_t bLastWasFalling=0;
 310     Int_t newbin=-1;
 311
 312
 313     if(fDeconvTime){
 314       redo: //This is a goto.
 315
 316       if(newbin > -1){
 317         //bin = newbin;
 318         newbin = -1;
 319       }
 320
 321       lastcharge=0;
 322       bLastWasFalling = 0;
 323     }
 324
 325     while(1){ //Loop over time bins of current pad
 326       // read all the values for one pad at once to calculate the base line
 327       if (pCurrentPad) {
 328         if (!pCurrentPad->IsStarted()) {
 329           //HLTDebug("reading data for pad %d, padrow %d", fDigitReader->GetPad(), fDigitReader->GetRow()+rowOffset);
 330           pCurrentPad->SetID(fDigitReader->GetRow()+rowOffset,fDigitReader->GetPad());
 331           if ((pCurrentPad->StartEvent())>=0) {
 332             do {
 333               if ((fDigitReader->GetRow()+rowOffset)!=pCurrentPad->GetRowNumber()) break;
 334               if (fDigitReader->GetPad()!=pCurrentPad->GetPadNumber()) break;
 335               pCurrentPad->SetRawData(fDigitReader->GetTime(), fDigitReader->GetSignal());
 336               //HLTDebug("set raw data to pad: bin %d charge %d", fDigitReader->GetTime(), fDigitReader->GetSignal());
 337             } while ((readValue = fDigitReader->Next())!=0);
 338           }
 339           pCurrentPad->CalculateBaseLine(AliHLTTPCTransform::GetNTimeBins()/2);
 340           if (pCurrentPad->Next(kTRUE/*do zero suppression*/)==0) {
 341             //HLTDebug("no data available after zero suppression");
 342             pCurrentPad->StopEvent();
 343             pCurrentPad->ResetHistory();
 344             break;
 345           }
 346           time=pCurrentPad->GetCurrentPosition();
 347           if (time>pCurrentPad->GetSize()) {
 348             HLTError("invalid time bin for pad");
 349             break;
 350           }
 351         }
 352       }
 353
 354       if (pCurrentPad) {
 355         Float_t occupancy=pCurrentPad->GetOccupancy();
 356         //HLTDebug("pad %d occupancy level: %f", pCurrentPad->GetPadNumber(), occupancy);
 357         if ( occupancy < fOccupancyLimit ) {
 358           charge = pCurrentPad->GetCorrectedData();
 359         } else {
 360           charge = 0;
 361           //HLTDebug("ignoring pad %d with occupancy level %f", pCurrentPad->GetPadNumber(), occupancy);
 362         }
 363       } else {
 364         charge = fDigitReader->GetSignal();
 365       }
 366       //HLTDebug("get next charge value: position %d charge %d", time, charge);
 367
 368
 369       // CHARGE DEBUG
 370       if (fDigitReader->GetRow() == 90){
 371 /////     LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::Row","row90")  << "PAD=" <<  fDigitReader->GetPad() << "  TIME=" <<  fDigitReader->GetTime()
 372           //                                       << "  SIGNAL=" <<  fDigitReader->GetSignal() << ENDLOG;
 373
 374       }
 375
 376       if(time >= AliHLTTPCTransform::GetNTimeBins()){
 377         HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
 378         break;
 379       }
 380
 381
 382       //Get the current ADC-value
 383       if(fDeconvTime){
 384
 385         //Check if the last pixel in the sequence is smaller than this
 386         if(charge > lastcharge){
 387           if(bLastWasFalling){
 388             newbin = 1;
 389             break;
 390           }
 391         }
 392         else bLastWasFalling = 1; //last pixel was larger than this
 393         lastcharge = charge;
 394       }
 395
 396       //Sum the total charge of this sequence
 397       seqcharge += charge;
 398       seqaverage += time*charge;
 399       seqerror += time*time*charge;
 400
 401       if (pCurrentPad) {
 402
 403         if((pCurrentPad->Next(kTRUE/*do zero suppression*/))==0) {
 404           pCurrentPad->StopEvent();
 405           pCurrentPad->ResetHistory();
 406           if(readValue) {
 407             newPad = fDigitReader->GetPad();
 408             newTime = fDigitReader->GetTime();
 409             newRow = fDigitReader->GetRow() + rowOffset;
 410           }
 411           break;
 412         }
 413
 414         newPad=pCurrentPad->GetPadNumber();
 415         newTime=pCurrentPad->GetCurrentPosition();
 416         newRow=pCurrentPad->GetRowNumber();
 417       } else {
 418       readValue = fDigitReader->Next();
 419       //Check where to stop:
 420       if(!readValue) break; //No more value
 421
 422       newPad = fDigitReader->GetPad();
 423       newTime = fDigitReader->GetTime();
 424       newRow = fDigitReader->GetRow() + rowOffset;
 425       }
 426
 427       if(newPad != pad)break; //new pad
 428       if(newTime != time+1) break; //end of sequence
 429
 430       // pad = newpad;    is equal
 431       time = newTime;
 432
 433     }//end loop over sequence
 434
 435     //HLTDebug("ended time bin sequence loop: seqcharge=%d readValue=%d", seqcharge, readValue);
 436     //HLTDebug("pad=%d newpad=%d current row=%d newrow=%d", pad, newPad, fCurrentRow, newRow);
 437     if (seqcharge<=0) {
 438       // with active zero suppression zero values are possible
 439       continue;
 440     }
 441
 442     //Calculate mean of sequence:
 443     Int_t seqmean=0;
 444     if(seqcharge)
 445       seqmean = seqaverage/seqcharge;
 446     else{
 447       LOG(AliHLTTPCLog::kFatal,"AliHLTTPCClusterFinder::ProcessRow","Data")
 448         <<"Error in data given to the cluster finder"<<ENDLOG;
 449       seqmean = 1;
 450       seqcharge = 1;
 451     }
 452
 453     //Calculate mean in pad direction:
 454     Int_t padmean = seqcharge*pad;
 455     Int_t paderror = pad*padmean;
 456
 457     //Compare with results on previous pad:
 458     for(UInt_t p=0; p<nprevious && p<kPadArraySize && ncurrent<kPadArraySize; p++){
 459
 460       //dont merge sequences on the same pad twice
 461       if(previousPt[p]->fLastMergedPad==pad) continue;
 462
 463       Int_t difference = seqmean - previousPt[p]->fMean;
 464       if(difference < -fMatch) break;
 465
 466       if(difference <= fMatch){ //There is a match here!!
 467         AliClusterData *local = previousPt[p];
 468
 469         if(fDeconvPad){
 470           if(seqcharge > local->fLastCharge){
 471             if(local->fChargeFalling){ //The previous pad was falling
 472               break; //create a new cluster
 473             }
 474           }
 475           else local->fChargeFalling = 1;
 476           local->fLastCharge = seqcharge;
 477         }
 478
 479         //Don't create a new cluster, because we found a match
 480         newcluster = kFALSE;
 481
 482         //Update cluster on current pad with the matching one:
 483         local->fTotalCharge += seqcharge;
 484         local->fPad += padmean;
 485         local->fPad2 += paderror;
 486         local->fTime += seqaverage;
 487         local->fTime2 += seqerror;
 488         local->fMean = seqmean;
 489         local->fFlags++; //means we have more than one pad
 490         local->fLastMergedPad = pad;
 491
 492         currentPt[ncurrent] = local;
 493         ncurrent++;
 494
 495         break;
 496       } //Checking for match at previous pad
 497     } //Loop over results on previous pad.
 498
 499     if(newcluster && ncurrent<kPadArraySize){
 500       //Start a new cluster. Add it to the clusterlist, and update
 501       //the list of pointers to clusters in current pad.
 502       //current pad will be previous pad on next pad.
 503
 504       //Add to the clusterlist:
 505       AliClusterData *tmp = &clusterlist[ntotal];
 506       tmp->fTotalCharge = seqcharge;
 507       tmp->fPad = padmean;
 508       tmp->fPad2 = paderror;
 509       tmp->fTime = seqaverage;
 510       tmp->fTime2 = seqerror;
 511       tmp->fMean = seqmean;
 512       tmp->fFlags = 0;  //flags for single pad clusters
 513       tmp->fLastMergedPad = pad;
 514
 515       if(fDeconvPad){
 516         tmp->fChargeFalling = 0;
 517         tmp->fLastCharge = seqcharge;
 518       }
 519
 520       //Update list of pointers to previous pad:
 521       currentPt[ncurrent] = &clusterlist[ntotal];
 522       ntotal++;
 523       ncurrent++;
 524     }
 525
 526     if(fDeconvTime)
 527       if(newbin >= 0) goto redo;
 528
 529     // to prevent endless loop
 530     if(time >= AliHLTTPCTransform::GetNTimeBins()){
 531       HLTWarning("Timebin (%d) out of range (%d)", time, AliHLTTPCTransform::GetNTimeBins());
 532       break;
 533     }
 534
 535
 536     if(!readValue) break; //No more value
 537
 538     if (ntotal>=kClusterListSize || ncurrent>=kPadArraySize) {
 539       HLTWarning("pad array size exceeded ntotal=%d ncurrent=%d, skip rest of the data", ntotal, ncurrent);
 540       break;
 541     }
 542
 543     if(fCurrentRow != newRow){
 544       WriteClusters(ntotal,clusterlist);
 545
 546       lastpad = 123456789;
 547
 548       currentPt = pad2;
 549       previousPt = pad1;
 550       nprevious=0;
 551       ncurrent=0;
 552       ntotal=0;
 553
 554       fCurrentRow = newRow;
 555     }
 556
 557     pad = newPad;
 558     time = newTime;
 559
 560   } // END while(readValue)
 561
 562   WriteClusters(ntotal,clusterlist);
 563
 564   HLTInfo("ClusterFinder found %d clusters in slice %d patch %d", fNClusters, fCurrentSlice, fCurrentPatch);
 565
 566 } // ENDEND
 567
 568 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliClusterData *list)
 569 {
 570   //write cluster to output pointer
 571   Int_t thisrow,thissector;
 572   UInt_t counter = fNClusters;
 573
 574   for(int j=0; j<nclusters; j++)
 575     {
 576       if(!list[j].fFlags) continue; //discard single pad clusters
 577       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
 578
 579       Float_t xyz[3];
 580       Float_t fpad =(Float_t)list[j].fPad / list[j].fTotalCharge;
 581       Float_t fpad2=fXYErr*fXYErr; //fixed given error
 582       Float_t ftime =(Float_t)list[j].fTime / list[j].fTotalCharge;
 583       Float_t ftime2=fZErr*fZErr;  //fixed given error
 584
 585
 586
 587
 588
 589
 590       if(fCalcerr) { //calc the errors, otherwice take the fixed error
 591         Int_t patch = AliHLTTPCTransform::GetPatch(fCurrentRow);
 592         UInt_t q2=list[j].fTotalCharge*list[j].fTotalCharge;
 593         Float_t sy2=list[j].fPad2 * list[j].fTotalCharge - list[j].fPad * list[j].fPad;
 594         sy2/=q2;
 595         if(sy2 < 0) {
 596             LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
 597               <<"SigmaY2 negative "<<sy2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
 598             continue;
 599         } else {
 600           if(!fRawSP){
 601             fpad2 = (sy2 + 1./12)*AliHLTTPCTransform::GetPadPitchWidth(patch)*AliHLTTPCTransform::GetPadPitchWidth(patch);
 602             if(sy2 != 0){
 603               fpad2*=0.108; //constants are from offline studies
 604               if(patch<2)
 605                 fpad2*=2.07;
 606             }
 607           } else fpad2=sy2; //take the width not the error
 608         }
 609         Float_t sz2=list[j].fTime2*list[j].fTotalCharge - list[j].fTime*list[j].fTime;
 610         sz2/=q2;
 611         if(sz2 < 0){
 612           LOG(AliHLTTPCLog::kError,"AliHLTTPCClusterFinder::WriteClusters","Cluster width")
 613             <<"SigmaZ2 negative "<<sz2<<" on row "<<fCurrentRow<<" "<<fpad<<" "<<ftime<<ENDLOG;
 614           continue;
 615         } else {
 616           if(!fRawSP){
 617             ftime2 = (sz2 + 1./12)*AliHLTTPCTransform::GetZWidth()*AliHLTTPCTransform::GetZWidth();
 618             if(sz2 != 0) {
 619               ftime2 *= 0.169; //constants are from offline studies
 620               if(patch<2)
 621                 ftime2 *= 1.77;
 622             }
 623           } else ftime2=sz2; //take the width, not the error
 624         }
 625       }
 626       if(fStdout==kTRUE)
 627         cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
 628
 629       if(!fRawSP){
 630         AliHLTTPCTransform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
 631         AliHLTTPCTransform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
 632
 633         if(xyz[0]==0) LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder","Cluster Finder")
 634           <<AliHLTTPCLog::kDec<<"Zero cluster"<<ENDLOG;
 635         if(fNClusters >= fMaxNClusters)
 636           {
 637             LOG(AliHLTTPCLog::kError,"AliHLTTPCClustFinder::WriteClusters","Cluster Finder")
 638               <<AliHLTTPCLog::kDec<<"Too many clusters "<<fNClusters<<ENDLOG;
 639             return;
 640           }
 641
 642         fSpacePointData[counter].fX = xyz[0];
 643         fSpacePointData[counter].fY = xyz[1];
 644         fSpacePointData[counter].fZ = xyz[2];
 645
 646       } else {
 647         fSpacePointData[counter].fX = fCurrentRow;
 648         fSpacePointData[counter].fY = fpad;
 649         fSpacePointData[counter].fZ = ftime;
 650       }
 651
 652       fSpacePointData[counter].fCharge = list[j].fTotalCharge;
 653       fSpacePointData[counter].fPadRow = fCurrentRow;
 654       fSpacePointData[counter].fSigmaY2 = fpad2;
 655       fSpacePointData[counter].fSigmaZ2  = ftime2;
 656
 657       fSpacePointData[counter].fUsed = kFALSE;         // only used / set in AliHLTTPCDisplay
 658       fSpacePointData[counter].fTrackN = -1;           // only used / set in AliHLTTPCDisplay
 659
 660       Int_t patch=fCurrentPatch;
 661       if(patch==-1) patch=0; //never store negative patch number
 662       fSpacePointData[counter].fID = counter
 663         +((fCurrentSlice&0x7f)<<25)+((patch&0x7)<<22);//Uli
 664
 665 #ifdef do_mc
 666       Int_t trackID[3];
 667       GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
 668
 669       fSpacePointData[counter].fTrackID[0] = trackID[0];
 670       fSpacePointData[counter].fTrackID[1] = trackID[1];
 671       fSpacePointData[counter].fTrackID[2] = trackID[2];
 672
 673       //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
 674 #endif
 675
 676       fNClusters++;
 677       counter++;
 678     }
 679 }
 680
 681 // STILL TO FIX  ----------------------------------------------------------------------------
 682
 683 #ifdef do_mc
 684 void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
 685 {
 686   //get mc id
 687   AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fDigitRowData;
 688
 689   trackID[0]=trackID[1]=trackID[2]=-2;
 690   //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
 691   for(Int_t i=fFirstRow; i<=fLastRow; i++){
 692     if(rowPt->fRow < (UInt_t)fCurrentRow){
 693       AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
 694       continue;
 695     }
 696     AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
 697     for(UInt_t j=0; j<rowPt->fNDigit; j++){
 698       Int_t cpad = digPt[j].fPad;
 699       Int_t ctime = digPt[j].fTime;
 700       if(cpad != pad) continue;
 701       if(ctime != time) continue;
 702
 703       trackID[0] = digPt[j].fTrackID[0];
 704       trackID[1] = digPt[j].fTrackID[1];
 705       trackID[2] = digPt[j].fTrackID[2];
 706
 707       //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
 708       break;
 709     }
 710     break;
 711   }
 712 }
 713 #endif