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