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