HLT/TPCLib/AliHLTTPCClusterFinder.cxx

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