HLT/TPCLib/AliHLTTPCClusterFinder.cxx

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