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