HLT/TPCLib/AliHLTTPCPad.cxx

   1 // @(#) $Id$
   2
   3 /**************************************************************************
   4  * This file is property of and copyright by the ALICE HLT Project        *
   5  * ALICE Experiment at CERN, All rights reserved.                         *
   6  *                                                                        *
   7  * Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
   8  *                  Kenneth Aamodt   <Kenneth.aamodt@ift.uib.no>          *
   9  *                  for The ALICE HLT Project.                            *
  10  *                                                                        *
  11  * Permission to use, copy, modify and distribute this software and its   *
  12  * documentation strictly for non-commercial purposes is hereby granted   *
  13  * without fee, provided that the above copyright notice appears in all   *
  14  * copies and that both the copyright notice and this permission notice   *
  15  * appear in the supporting documentation. The authors make no claims     *
  16  * about the suitability of this software for any purpose. It is          *
  17  * provided "as is" without express or implied warranty.                  *
  18  **************************************************************************/
  19
  20 /** @file   AliHLTTPCPad.cxx
  21     @author Matthias Richter, Kenneth Aamodt
  22     @date
  23     @brief  Container Class for TPC Pads.
  24 */
  25
  26 // see header file for class documentation
  27 // or
  28 // refer to README to build package
  29 // or
  30 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  31
  32 #if __GNUC__>= 3
  33 using namespace std;
  34 #endif
  35
  36 #include <cerrno>
  37 #include "AliHLTTPCPad.h"
  38 #include "AliHLTStdIncludes.h"
  39
  40
  41 //added by kenneth
  42 #include "AliHLTTPCTransform.h"
  43 #include "AliHLTTPCClusters.h"
  44 #include <sys/time.h>
  45 //------------------------------
  46
  47 /** margin for the base line be re-avaluated */
  48 #define ALIHLTPAD_BASELINE_MARGIN (2*fAverage)
  49
  50 /** ROOT macro for the implementation of ROOT specific class methods */
  51 ClassImp(AliHLTTPCPad)
  52
  53 AliHLTTPCPad::AliHLTTPCPad()
  54   :
  55   fClusterCandidates(),
  56   fUsedClusterCandidates(),
  57   fRowNo(-1),
  58   fPadNo(-1),
  59   fThreshold(0),
  60   fAverage(-1),
  61   fNofEvents(0),
  62   fSum(0),
  63   fCount(0),
  64   fTotal(0),
  65   fBLMax(-1),
  66   fBLMaxBin(-1),
  67   fBLMin(-1),
  68   fBLMinBin(-1),
  69   fFirstBLBin(0),
  70   fNofBins(0),
  71   fReadPos(0),
  72   fpRawData(NULL),
  73   fDataSignals(NULL),
  74   fSignalPositionArray(NULL),
  75   fSizeOfSignalPositionArray(0),
  76   fNSigmaThreshold(0),
  77   fSignalThreshold(0)
  78 {
  79   // see header file for class documentation
  80   // or
  81   // refer to README to build package
  82   // or
  83   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  84   //  HLTInfo("Entering default constructor");
  85   fDataSignals= new AliHLTTPCSignal_t[AliHLTTPCTransform::GetNTimeBins()];
  86   memset( fDataSignals, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
  87
  88   fSignalPositionArray= new AliHLTTPCSignal_t[AliHLTTPCTransform::GetNTimeBins()];
  89   memset( fSignalPositionArray, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
  90   fSizeOfSignalPositionArray=0;
  91 }
  92
  93 AliHLTTPCPad::AliHLTTPCPad(Int_t offset, Int_t nofBins)
  94   :
  95   fClusterCandidates(),
  96   fUsedClusterCandidates(),
  97   fRowNo(-1),
  98   fPadNo(-1),
  99   fThreshold(0),
 100   fAverage(-1),
 101   fNofEvents(0),
 102   fSum(0),
 103   fCount(0),
 104   fTotal(0),
 105   fBLMax(-1),
 106   fBLMaxBin(-1),
 107   fBLMin(-1),
 108   fBLMinBin(-1),
 109   fFirstBLBin(offset),
 110   fNofBins(nofBins),
 111   fReadPos(0),
 112   fpRawData(NULL),
 113   fDataSignals(NULL),
 114   fSignalPositionArray(NULL),
 115   fSizeOfSignalPositionArray(0),
 116   fNSigmaThreshold(0),
 117   fSignalThreshold(0)
 118 {
 119   // see header file for class documentation
 120 }
 121
 122 AliHLTTPCPad::~AliHLTTPCPad()
 123 {
 124   // see header file for class documentation
 125   if (fpRawData) {
 126     HLTWarning("event data acquisition not stopped");
 127     StopEvent();
 128   }
 129   if (fDataSignals) {
 130     AliHLTTPCSignal_t* pData=fDataSignals;
 131     fDataSignals=NULL;
 132     delete [] pData;
 133   }
 134   if (fSignalPositionArray) {
 135     //AliHLTTPCSignal_t* pData=fSignalPositionArray;
 136     fSignalPositionArray=NULL;
 137     //   delete [] pData;
 138   }
 139
 140 }
 141
 142 Int_t AliHLTTPCPad::SetID(Int_t rowno, Int_t padno)
 143 {
 144   // see header file for class documentation
 145   fRowNo=rowno;
 146   fPadNo=padno;
 147   return 0;
 148 }
 149
 150 Int_t AliHLTTPCPad::StartEvent()
 151 {
 152   // see header file for class documentation
 153   Int_t iResult=0;
 154   if (fpRawData==NULL) {
 155     fBLMax=-1;
 156     fBLMaxBin=-1;
 157     fBLMin=-1;
 158     fBLMinBin=-1;
 159     fSum=0;
 160     fCount=0;
 161     fTotal=0;
 162     if (fNofBins>0) {
 163       fpRawData=new AliHLTTPCSignal_t[fNofBins];
 164       if (fpRawData) {
 165         for (int i=0; i<fNofBins; i++) fpRawData[i]=-1;
 166       } else {
 167         HLTError("memory allocation failed");
 168         iResult=-ENOMEM;
 169       }
 170     }
 171   } else {
 172     HLTWarning("event data acquisition already started");
 173     iResult=-EALREADY;
 174   }
 175   return iResult;
 176 }
 177
 178 Int_t AliHLTTPCPad::CalculateBaseLine(Int_t reqMinCount)
 179 {
 180   // see header file for class documentation
 181   Int_t iResult=0;
 182   AliHLTTPCSignal_t avBackup=fAverage;
 183   //HLTDebug("reqMinCount=%d fCount=%d fTotal=%d fSum=%d fBLMax=%d fBLMin=%d", reqMinCount, fCount, fTotal, fSum, fBLMax, fBLMin);
 184   if (fCount>=reqMinCount && fCount>=fTotal/2) {
 185     fAverage=fCount>0?fSum/fCount:0;
 186     if (fAverage>0) {
 187       //HLTDebug("average for current event %d (%d - %d)", fAverage, fBLMax, fBLMin);
 188       fCount=0;fSum=-1;
 189       if (fBLMax>ALIHLTPAD_BASELINE_MARGIN) {
 190         // calculate again
 191         //HLTDebug("maximum value %d exceeds margin for base line (%d) "
 192         //       "-> re-evaluate base line", fBLMax, ALIHLTPAD_BASELINE_MARGIN);
 193         if (fpRawData) {
 194           for (Int_t i=fFirstBLBin; i<fNofBins; i++)
 195             if (fpRawData[i]>=0) AddBaseLineValue(i, fpRawData[i]);
 196           if (fCount>0 && fCount>=reqMinCount && fCount>=fTotal/2) {
 197             fAverage=fSum/fCount;
 198             //HLTDebug("new average %d", fAverage);
 199           } else {
 200             //      HLTDebug("baseline re-eveluation skipped because of to few "
 201             //                 "contributing bins: total=%d, contributing=%d, req=%d"
 202             //                 "\ndata might be already zero suppressed"
 203             //                 , fTotal, fCount, reqMinCount);
 204             iResult=-ENODATA;
 205           }
 206           fCount=0;fSum=-1;
 207         } else {
 208           HLTError("missing raw data for base line calculation");
 209           iResult=-ENOBUFS;
 210         }
 211       }
 212       if (iResult>=0) {
 213         // calculate average for all events
 214         fAverage=((avBackup*fNofEvents)+fAverage)/(fNofEvents+1);
 215         //HLTDebug("base line average for %d event(s): %d", fNofEvents+1, fAverage);
 216       } else {
 217         fAverage=avBackup;
 218       }
 219     } else {
 220       fAverage=avBackup;
 221     }
 222   } else {
 223     //     HLTDebug("baseline calculation skipped because of to few contributing "
 224     //         "bins: total=%d, contributing=%d, required=%d \ndata might be "
 225     //         "already zero suppressed", fTotal, fCount, reqMinCount);
 226   }
 227
 228   return iResult;
 229 }
 230
 231 Int_t AliHLTTPCPad::StopEvent()
 232 {
 233   // see header file for class documentation
 234   Int_t iResult=0;
 235   if (fpRawData) {
 236     AliHLTTPCSignal_t* pData=fpRawData;
 237     fpRawData=NULL;
 238     delete [] pData;
 239     fTotal=0;
 240     fNofEvents++;
 241     Rewind();
 242   } else if (fNofBins>0) {
 243     HLTError("event data acquisition not started");
 244     iResult=-EBADF;
 245   }
 246   return iResult;
 247 }
 248
 249 Int_t AliHLTTPCPad::ResetHistory()
 250 {
 251   // see header file for class documentation
 252   Int_t iResult=0;
 253   fAverage=-1;
 254   fNofEvents=0;
 255   return iResult;
 256 }
 257
 258 Int_t AliHLTTPCPad::SetThreshold(AliHLTTPCSignal_t thresh)
 259 {
 260   // see header file for class documentation
 261   Int_t iResult=0;
 262   fThreshold=thresh;
 263   return iResult;
 264 }
 265
 266 Int_t AliHLTTPCPad::AddBaseLineValue(Int_t bin, AliHLTTPCSignal_t value)
 267 {
 268   // see header file for class documentation
 269   Int_t iResult=0;
 270   if (bin>=fFirstBLBin) {
 271     if (fAverage<0 || value<ALIHLTPAD_BASELINE_MARGIN) {
 272       // add to the current sum and count
 273       fSum+=value;
 274       fCount++;
 275       if (fBLMax<value) {
 276         // keep the maximum value for later quality control of the base
 277         // line calculation
 278         fBLMax=value;
 279         fBLMaxBin=bin;
 280       }
 281       if (fBLMin<0 || fBLMin>value) {
 282         // keep the minimum value for later quality control of the base
 283         // line calculation
 284         fBLMin=value;
 285         fBLMinBin=bin;
 286       }
 287     } else {
 288       //       HLTDebug("ignoring value %d (bin %d) for base line calculation "
 289       //               "(current average is %d)",
 290       //               value, bin, fAverage);
 291     }
 292   }
 293   return iResult;
 294 }
 295
 296 Int_t AliHLTTPCPad::SetRawData(Int_t bin, AliHLTTPCSignal_t value)
 297 {
 298   // see header file for class documentation
 299   Int_t iResult=0;
 300   if (fpRawData) {
 301     if (bin<fNofBins) {
 302       if (value>=0) {
 303         if (fpRawData[bin]<0) {
 304           AddBaseLineValue(bin, value);
 305           fTotal++;
 306         } else {
 307           // ignore value for average calculation
 308           HLTWarning("overriding content of bin %d (%d)", bin, fpRawData[bin]);
 309         }
 310         fpRawData[bin]=value;
 311       } else {
 312         HLTWarning("ignoring neg. raw data");
 313       }
 314     } else {
 315       HLTWarning("bin %d out of range (%d)", bin, fNofBins);
 316       iResult=-ERANGE;
 317     }
 318   } else if (fNofBins>0) {
 319     HLTError("event cycle not started");
 320     iResult=-EBADF;
 321   }
 322   return iResult;
 323 }
 324
 325 Int_t AliHLTTPCPad::Next(Int_t bZeroSuppression)
 326 {
 327   // see header file for class documentation
 328   if (fpRawData==NULL) return 0;
 329   Int_t iResult=fReadPos<fNofBins;
 330   if (iResult>0 && (iResult=(++fReadPos<fNofBins))>0) {
 331     if (bZeroSuppression) {
 332       while ((iResult=(fReadPos<fNofBins))>0 &&
 333              GetCorrectedData(fReadPos)<=0)
 334         fReadPos++;
 335     }
 336   }
 337   return iResult;
 338 }
 339
 340 Int_t AliHLTTPCPad::Rewind(Int_t bZeroSuppression)
 341 {
 342   // see header file for class documentation
 343   fReadPos=(bZeroSuppression>0?0:fFirstBLBin)-1;
 344   return Next(bZeroSuppression);
 345 }
 346
 347 AliHLTTPCSignal_t AliHLTTPCPad::GetRawData(Int_t bin) const
 348 {
 349   // see header file for class documentation
 350   AliHLTTPCSignal_t data=0;
 351   if (fpRawData) {
 352     if (bin<fNofBins) {
 353       data=fpRawData[bin];
 354     } else {
 355       HLTWarning("requested bin %d out of range (%d)", bin, fNofBins);
 356     }
 357   } else if (fNofBins>0) {
 358     HLTWarning("data only available within event cycle");
 359   }
 360   return data;
 361 }
 362
 363 AliHLTTPCSignal_t AliHLTTPCPad::GetCorrectedData(Int_t bin) const
 364 {
 365   // see header file for class documentation
 366   AliHLTTPCSignal_t data=GetRawData(bin)-GetBaseLine(bin);
 367   AliHLTTPCSignal_t prev=0;
 368   if (bin>1) prev=GetRawData(bin-1)-GetBaseLine(bin-1);
 369   AliHLTTPCSignal_t succ=0;
 370   if (bin+1<GetSize()) succ=GetRawData(bin+1)-GetBaseLine(bin+1);
 371   if (fThreshold>0) {
 372     data-=fThreshold;
 373     prev-=fThreshold;
 374     succ-=fThreshold;
 375   }
 376
 377   // case 1:
 378   // the signal is below the base-line and threshold
 379   if (data<0) data=0;
 380
 381   //case 2:
 382   // the neighboring bins are both below base-line/threshold
 383   // a real signal is always more than one bin wide because of the shaper
 384   if (prev<=0 && succ<=0) data=0;
 385
 386   // case 3:
 387   // the bin is inside the range of ignored bins
 388   if (bin<fFirstBLBin) data=0;
 389   //HLTDebug("fReadPos=%d data=%d threshold=%d raw data=%d base line=%d", fReadPos, data, fThreshold, GetRawData(bin), GetBaseLine(bin));
 390   return data;
 391 }
 392
 393 AliHLTTPCSignal_t AliHLTTPCPad::GetBaseLine(Int_t /*bin*/) const //TODO: Why is bin being ignored?
 394 {
 395   // see header file for class documentation
 396   AliHLTTPCSignal_t val=0;
 397   if (fAverage>0) {
 398     // we take the minumum value as the base line if it doesn't differ from
 399     // the average to much
 400     val=fAverage;
 401 #ifdef KEEP_NOISE
 402     const AliHLTTPCSignal_t kMaxDifference=15;
 403     if ((fAverage-fBLMin)<=kMaxDifference) val=fBLMin;
 404     else val>kMaxDifference?val-=kMaxDifference:0;
 405 #endif
 406   }
 407   if (val<0) {
 408     // here we should never get
 409     val=0;
 410     HLTFatal("wrong base line value");
 411   }
 412   return val;
 413 }
 414
 415 AliHLTTPCSignal_t AliHLTTPCPad::GetAverage() const
 416 {
 417   // see header file for class documentation
 418   return fAverage>0?fAverage:0;
 419 }
 420
 421 Float_t AliHLTTPCPad::GetOccupancy() const
 422 {
 423   // see header file for class documentation
 424   Float_t occupancy=0;
 425   if (fpRawData && fNofBins>0) {
 426     for (Int_t i=fFirstBLBin; i<fNofBins; i++) {
 427       if (GetCorrectedData(i)>0) occupancy+=1;
 428     }
 429     if (fNofBins-fFirstBLBin>0)
 430       occupancy/=fNofBins-fFirstBLBin;
 431   }
 432   return occupancy;
 433 }
 434
 435 Float_t AliHLTTPCPad::GetAveragedOccupancy() const
 436 {
 437   // see header file for class documentation
 438
 439   // history is not yet implemented
 440   return GetOccupancy();
 441 }
 442 void AliHLTTPCPad::PrintRawData()
 443 {
 444   // see header file for class documentation
 445   for(Int_t bin=0;bin<AliHLTTPCTransform::GetNTimeBins();bin++){
 446     if(GetDataSignal(bin)>0)
 447       cout<<fRowNo<<"\t"<<fPadNo<<"\t"<<bin<<"\t"<<GetDataSignal(bin)<<endl;;
 448   }
 449   cout<<"bins: "<<AliHLTTPCTransform::GetNTimeBins()<<endl;
 450 }
 451
 452 void AliHLTTPCPad::SetDataToDefault()
 453 {
 454   // see header file for class documentation
 455   if(fpRawData){
 456     memset( fDataSignals, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
 457     memset( fSignalPositionArray, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
 458     fSizeOfSignalPositionArray=0;
 459   }
 460 }
 461
 462 void AliHLTTPCPad::SetDataSignal(Int_t bin,Int_t signal)
 463 {
 464   // see header file for class documentation
 465   fDataSignals[bin]=signal;
 466   fSignalPositionArray[fSizeOfSignalPositionArray]=bin;
 467   fSizeOfSignalPositionArray++;
 468 }
 469
 470 Int_t AliHLTTPCPad::GetDataSignal(Int_t bin) const
 471 {
 472   // see header file for class documentation
 473   return fDataSignals[bin];
 474 }
 475
 476 void AliHLTTPCPad::ZeroSuppress(Double_t nSigma = 3,Int_t threshold = 20 ,Int_t reqMinPoint = AliHLTTPCTransform::GetNTimeBins()/2, Int_t beginTime = 50,Int_t endTime = AliHLTTPCTransform::GetNTimeBins()-1){
 477   //see headerfile for documentation
 478
 479   Bool_t useSigma= kFALSE;
 480   if(nSigma>0){
 481     useSigma=kTRUE;
 482   }
 483   if(threshold<1 && nSigma<=0){
 484     //setting the data to -1 for this pad
 485     memset( fDataSignals, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
 486     fSizeOfSignalPositionArray=0;
 487     return;
 488   }
 489   if(endTime>=AliHLTTPCTransform::GetNTimeBins()){
 490     endTime=AliHLTTPCTransform::GetNTimeBins()-1;
 491   }
 492
 493   Int_t fThresholdUsed=threshold;
 494
 495   Int_t nAdded=0;
 496   Int_t sumNAdded=0;
 497   fSizeOfSignalPositionArray=0;
 498   for(Int_t i=beginTime;i<endTime+1;i++){
 499     if(fDataSignals[i]>0){
 500       nAdded++;
 501       sumNAdded+=fDataSignals[i];
 502     }
 503   }
 504
 505   if(nAdded<reqMinPoint){
 506     return;      //This will ensure that no data is read in FindClusterCandidates() (since fSizeOfSignalPositionArray=0)
 507   }
 508
 509   Double_t averageValue=sumNAdded/nAdded;
 510
 511   Double_t sigma=0;
 512   if(useSigma){
 513     //Calculate the sigma
 514     Double_t sumOfDifferenceSquared=0;
 515     for(Int_t i=endTime;i>=beginTime;i--){
 516       if(fDataSignals[i]>0){
 517         if(fDataSignals[i]-averageValue<50){
 518           sumOfDifferenceSquared+=(fDataSignals[i]-averageValue)*(fDataSignals[i]-averageValue);
 519         }
 520         else{
 521           nAdded--;
 522         }
 523       }
 524     }
 525     sigma=sumOfDifferenceSquared/nAdded;
 526     fThresholdUsed=(Int_t)(nSigma*sigma);
 527   }
 528
 529   //For now just set the adc value outside [beginTime,endTime] to -1
 530   for(Int_t i=0;i<beginTime;i++){
 531     fDataSignals[i]=-1;
 532   }
 533   for(Int_t i=endTime+1;i<AliHLTTPCTransform::GetNTimeBins();i++){
 534     fDataSignals[i]=-1;
 535   }
 536
 537   // Do zero suppression on the adc values within [beginTime,endTime]
 538   for(Int_t i=endTime;i>=beginTime;i--){
 539     //the +1 in the if below is there to avoid getting a signal which is 0, adding to the numbers you have to loop over in the end
 540     //(better to set it to -1 which is default for no signal)
 541     if(fDataSignals[i]>(Int_t)(averageValue+fThresholdUsed+1) && fDataSignals[i-1]>(Int_t)(averageValue+fThresholdUsed+1)){
 542       //here the signals below threshold to the right of the candidate is added
 543       Bool_t contRight=kTRUE;
 544       Int_t endRight=i;
 545       Int_t nToAddRight=0;
 546       while(contRight){
 547         if(endRight+1<endTime){
 548           //      cout<<fDataSignals[endRight+1]<<"    "<<fDataSignals[endRight+2]<<endl;;
 549           if(fDataSignals[endRight+1]>=fDataSignals[endRight+2] && fDataSignals[endRight+1]>averageValue){
 550             nToAddRight++;
 551           }
 552           else{
 553             if(fDataSignals[endRight+1]> averageValue){
 554               nToAddRight++;
 555             }
 556             contRight=kFALSE;
 557           }
 558         }
 559         else if(endRight>endTime+1){
 560           contRight=kFALSE;
 561         }
 562         endRight++;
 563       }
 564       for(int j=i+nToAddRight;j>i;j--){
 565         fDataSignals[j]=(Int_t)(fDataSignals[j]-averageValue);
 566         fSignalPositionArray[fSizeOfSignalPositionArray]=j;
 567         fSizeOfSignalPositionArray++;
 568       }
 569
 570
 571       //before while the two consecutive timebin values are added
 572       fDataSignals[i]=(Int_t)(fDataSignals[i]-averageValue);
 573       fSignalPositionArray[fSizeOfSignalPositionArray]=i;
 574       fSizeOfSignalPositionArray++;
 575       fDataSignals[i-1]=(Int_t)(fDataSignals[i-1]-averageValue);
 576       fSignalPositionArray[fSizeOfSignalPositionArray]=i-1;
 577       fSizeOfSignalPositionArray++;
 578       i--;
 579       //      cout<<""<<endl;
 580       //Here the consecutive pads after the two first are added
 581       if(i-1>0){
 582         while(fDataSignals[i-1]>(Int_t)(averageValue+fThresholdUsed+1)){
 583           fDataSignals[i-1]=(Int_t)(fDataSignals[i-1]-averageValue);
 584           fSignalPositionArray[fSizeOfSignalPositionArray]=i-1;
 585           fSizeOfSignalPositionArray++;
 586           i--;
 587         }
 588       }
 589       //adding the signal below threshold belonging to the total signal
 590       Bool_t contLeft=kTRUE;
 591       while(contLeft){
 592         if(i-2>0){
 593           if(fDataSignals[i-1]>=fDataSignals[i-2] && fDataSignals[i-1]>averageValue){
 594             fDataSignals[i-1]=(Int_t)(fDataSignals[i-1]-averageValue);
 595             fSignalPositionArray[fSizeOfSignalPositionArray]=i-1;
 596             fSizeOfSignalPositionArray++;
 597             i--;
 598           }
 599           else{
 600             if(fDataSignals[i-1]> averageValue){
 601               fDataSignals[i-1]=(Int_t)(fDataSignals[i-1]-averageValue);
 602               fSignalPositionArray[fSizeOfSignalPositionArray]=i-1;
 603               fSizeOfSignalPositionArray++;
 604               i--;
 605             }
 606             contLeft=kFALSE;
 607           }
 608         }
 609         else{
 610           contLeft=kFALSE;
 611         }
 612
 613       }
 614     }
 615   }
 616   Int_t nReadFromPositionArray=0;
 617   for(Int_t i=endTime;i>=beginTime;i--){
 618     if(i==fSignalPositionArray[nReadFromPositionArray]){
 619       nReadFromPositionArray++;
 620     }
 621     else{
 622       fDataSignals[i]=-1;
 623     }
 624   }
 625 }
 626
 627 void AliHLTTPCPad::FindClusterCandidates()
 628 {
 629   // see header file for class documentation
 630
 631   if(fSizeOfSignalPositionArray<2){
 632     return;
 633   }
 634
 635   if(fNSigmaThreshold>0){
 636     ZeroSuppress(fNSigmaThreshold);
 637   }
 638   else if(fSignalThreshold>0){
 639     ZeroSuppress((Double_t)0,(Int_t)fSignalThreshold);
 640   }
 641   UInt_t seqcharge=0;
 642   UInt_t seqaverage=0;
 643   UInt_t seqerror=0;
 644   vector<Int_t> tmpPos;
 645   vector<Int_t> tmpSig;
 646   UInt_t isFalling=0;
 647
 648   for(Int_t pos=fSizeOfSignalPositionArray-2;pos>=0;pos--){
 649     if(fSignalPositionArray[pos]==fSignalPositionArray[pos+1]+1){
 650       seqcharge+=fDataSignals[fSignalPositionArray[pos+1]];
 651       seqaverage += fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
 652       seqerror += fSignalPositionArray[pos+1]*fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
 653
 654       tmpPos.push_back(fSignalPositionArray[pos+1]);
 655       tmpSig.push_back(fDataSignals[fSignalPositionArray[pos+1]]);
 656
 657       if(fDataSignals[fSignalPositionArray[pos+1]]>fDataSignals[fSignalPositionArray[pos]]){
 658         isFalling=1;
 659       }
 660       if(fDataSignals[fSignalPositionArray[pos+1]]<fDataSignals[fSignalPositionArray[pos]]&&isFalling){
 661         Int_t seqmean=0;
 662         seqmean = seqaverage/seqcharge;
 663
 664         //Calculate mean in pad direction:
 665         Int_t padmean = seqcharge*fPadNo;
 666         Int_t paderror = fPadNo*padmean;
 667         AliHLTTPCClusters candidate;
 668         candidate.fTotalCharge   = seqcharge;
 669         candidate.fPad       = padmean;
 670         candidate.fPad2      = paderror;
 671         candidate.fTime      = seqaverage;
 672         candidate.fTime2     = seqerror;
 673         candidate.fMean          = seqmean;
 674         candidate.fLastMergedPad = fPadNo;
 675         fClusterCandidates.push_back(candidate);
 676         fUsedClusterCandidates.push_back(0);
 677         isFalling=0;
 678         seqcharge=0;
 679         seqaverage=0;
 680         seqerror=0;
 681
 682         tmpPos.clear();
 683         tmpSig.clear();
 684
 685         continue;
 686       }
 687
 688       if(pos<1){
 689         seqcharge+=fDataSignals[fSignalPositionArray[0]];
 690         seqaverage += fSignalPositionArray[0]*fDataSignals[fSignalPositionArray[0]];
 691         seqerror += fSignalPositionArray[0]*fSignalPositionArray[0]*fDataSignals[fSignalPositionArray[0]];
 692         tmpPos.push_back(fSignalPositionArray[0]);
 693         tmpSig.push_back(fDataSignals[fSignalPositionArray[0]]);
 694
 695         //Calculate mean of sequence:
 696         Int_t seqmean=0;
 697         seqmean = seqaverage/seqcharge;
 698
 699         //Calculate mean in pad direction:
 700         Int_t padmean = seqcharge*fPadNo;
 701         Int_t paderror = fPadNo*padmean;
 702         AliHLTTPCClusters candidate;
 703         candidate.fTotalCharge   = seqcharge;
 704         candidate.fPad       = padmean;
 705         candidate.fPad2      = paderror;
 706         candidate.fTime      = seqaverage;
 707         candidate.fTime2     = seqerror;
 708         candidate.fMean          = seqmean;
 709         candidate.fLastMergedPad = fPadNo;
 710         fClusterCandidates.push_back(candidate);
 711         fUsedClusterCandidates.push_back(0);
 712         isFalling=0;
 713         seqcharge=0;
 714         seqaverage=0;
 715         seqerror=0;
 716
 717         tmpPos.clear();
 718         tmpSig.clear();
 719       }
 720     }
 721     else if(seqcharge>0){
 722       seqcharge+=fDataSignals[fSignalPositionArray[pos+1]];
 723       seqaverage += fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
 724       seqerror += fSignalPositionArray[pos+1]*fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
 725       tmpPos.push_back(fSignalPositionArray[pos+1]);
 726       tmpSig.push_back(fDataSignals[fSignalPositionArray[pos+1]]);
 727
 728       //Calculate mean of sequence:
 729       Int_t seqmean=0;
 730       seqmean = seqaverage/seqcharge;
 731
 732       //Calculate mean in pad direction:
 733       Int_t padmean = seqcharge*fPadNo;
 734       Int_t paderror = fPadNo*padmean;
 735       AliHLTTPCClusters candidate;
 736       candidate.fTotalCharge   = seqcharge;
 737       candidate.fPad       = padmean;
 738       candidate.fPad2      = paderror;
 739       candidate.fTime      = seqaverage;
 740       candidate.fTime2     = seqerror;
 741       candidate.fMean          = seqmean;
 742       candidate.fLastMergedPad = fPadNo;
 743       fClusterCandidates.push_back(candidate);
 744       fUsedClusterCandidates.push_back(0);
 745       isFalling=0;
 746       seqcharge=0;
 747       seqaverage=0;
 748       seqerror=0;
 749
 750       tmpPos.clear();
 751       tmpSig.clear();
 752     }
 753   }
 754 }
 755