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 #include "TMath.h"
  46 #include "TFile.h"
  47 //------------------------------
  48
  49 /** margin for the base line be re-avaluated */
  50 #define ALIHLTPAD_BASELINE_MARGIN (2*fAverage)
  51
  52 /** ROOT macro for the implementation of ROOT specific class methods */
  53 ClassImp(AliHLTTPCPad)
  54
  55 AliHLTTPCPad::AliHLTTPCPad()
  56   :
  57   fClusterCandidates(),
  58   fUsedClusterCandidates(),
  59   fRowNo(-1),
  60   fPadNo(-1),
  61   fThreshold(0),
  62   fAverage(-1),
  63   fNofEvents(0),
  64   fSum(0),
  65   fCount(0),
  66   fTotal(0),
  67   fBLMax(-1),
  68   fBLMaxBin(-1),
  69   fBLMin(-1),
  70   fBLMinBin(-1),
  71   fFirstBLBin(0),
  72   fNofBins(0),
  73   fReadPos(0),
  74   fpRawData(NULL),
  75   fDataSignals(NULL),
  76   fSignalPositionArray(NULL),
  77   fSizeOfSignalPositionArray(0),
  78   fNSigmaThreshold(0),
  79   fSignalThreshold(0),
  80   fModeSwitch(0),
  81   fNGoodSignalsSent(0),
  82   fDebugHistoBeforeZS(NULL),
  83   fDebugHistoAfterZS(NULL)
  84 {
  85   // see header file for class documentation
  86   // or
  87   // refer to README to build package
  88   // or
  89   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  90   //  HLTInfo("Entering default constructor");
  91   fDataSignals= new AliHLTTPCSignal_t[AliHLTTPCTransform::GetNTimeBins()];
  92   memset( fDataSignals, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
  93
  94   fSignalPositionArray= new AliHLTTPCSignal_t[AliHLTTPCTransform::GetNTimeBins()];
  95   memset( fSignalPositionArray, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
  96   fSizeOfSignalPositionArray=0;
  97
  98 }
  99
 100 AliHLTTPCPad::AliHLTTPCPad(Int_t mode)
 101   :
 102   fClusterCandidates(0),
 103   fUsedClusterCandidates(0),
 104   fRowNo(-1),
 105   fPadNo(-1),
 106   fThreshold(0),
 107   fAverage(-1),
 108   fNofEvents(0),
 109   fSum(0),
 110   fCount(0),
 111   fTotal(0),
 112   fBLMax(-1),
 113   fBLMaxBin(-1),
 114   fBLMin(-1),
 115   fBLMinBin(-1),
 116   fFirstBLBin(0),
 117   fNofBins(0),
 118   fReadPos(0),
 119   fpRawData(NULL),
 120   fDataSignals(NULL),
 121   fSignalPositionArray(NULL),
 122   fSizeOfSignalPositionArray(0),
 123   fNSigmaThreshold(0),
 124   fSignalThreshold(0),
 125   fModeSwitch(mode),
 126   fNGoodSignalsSent(0),
 127   fDebugHistoBeforeZS(NULL),
 128   fDebugHistoAfterZS(NULL)
 129 {
 130   // see header file for class documentation
 131   // or
 132   // refer to README to build package
 133   // or
 134   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
 135 }
 136
 137 AliHLTTPCPad::AliHLTTPCPad(Int_t offset, Int_t nofBins)
 138   :
 139   fClusterCandidates(),
 140   fUsedClusterCandidates(),
 141   fRowNo(-1),
 142   fPadNo(-1),
 143   fThreshold(0),
 144   fAverage(-1),
 145   fNofEvents(0),
 146   fSum(0),
 147   fCount(0),
 148   fTotal(0),
 149   fBLMax(-1),
 150   fBLMaxBin(-1),
 151   fBLMin(-1),
 152   fBLMinBin(-1),
 153   fFirstBLBin(offset),
 154   fNofBins(nofBins),
 155   fReadPos(0),
 156   fpRawData(NULL),
 157   fDataSignals(NULL),
 158   fSignalPositionArray(NULL),
 159   fSizeOfSignalPositionArray(0),
 160   fNSigmaThreshold(0),
 161   fSignalThreshold(0),
 162   fModeSwitch(0),
 163   fNGoodSignalsSent(0),
 164   fDebugHistoBeforeZS(NULL),
 165   fDebugHistoAfterZS(NULL)
 166 {
 167   // see header file for class documentation
 168 }
 169
 170 AliHLTTPCPad::~AliHLTTPCPad()
 171 {
 172   // see header file for class documentation
 173   if (fpRawData) {
 174     HLTWarning("event data acquisition not stopped");
 175     StopEvent();
 176   }
 177   if (fDataSignals) {
 178     delete [] fDataSignals;
 179     fDataSignals=NULL;
 180   }
 181   if (fSignalPositionArray) {
 182     delete [] fSignalPositionArray;
 183     fSignalPositionArray=NULL;
 184   }
 185   if(fDebugHistoBeforeZS){
 186     delete fDebugHistoBeforeZS;
 187     fDebugHistoBeforeZS=NULL;
 188   }
 189   if(fDebugHistoAfterZS){
 190     delete fDebugHistoAfterZS;
 191     fDebugHistoAfterZS=NULL;
 192   }
 193 }
 194
 195 Int_t AliHLTTPCPad::SetID(Int_t rowno, Int_t padno)
 196 {
 197   // see header file for class documentation
 198   fRowNo=rowno;
 199   fPadNo=padno;
 200
 201 #if DebugHisto
 202   char *nameBefore;
 203   sprintf(nameBefore,"beforeRow%dPad%d",fRowNo,fPadNo);
 204   char *nameAfter;
 205   sprintf(nameAfter,"afterRow%dPad%d",fRowNo,fPadNo);
 206   fDebugHistoBeforeZS = new TH1F(nameBefore,nameBefore,1024,0,1024);
 207   fDebugHistoAfterZS = new TH1F(nameAfter,nameAfter,1024,0,1024);
 208 #endif
 209
 210   return 0;
 211 }
 212
 213 Int_t AliHLTTPCPad::StartEvent()
 214 {
 215   // see header file for class documentation
 216   Int_t iResult=0;
 217   if (fpRawData==NULL) {
 218     fBLMax=-1;
 219     fBLMaxBin=-1;
 220     fBLMin=-1;
 221     fBLMinBin=-1;
 222     fSum=0;
 223     fCount=0;
 224     fTotal=0;
 225     if (fNofBins>0) {
 226       fpRawData=new AliHLTTPCSignal_t[fNofBins];
 227       if (fpRawData) {
 228         for (int i=0; i<fNofBins; i++) fpRawData[i]=-1;
 229       } else {
 230         HLTError("memory allocation failed");
 231         iResult=-ENOMEM;
 232       }
 233     }
 234   } else {
 235     HLTWarning("event data acquisition already started");
 236     iResult=-EALREADY;
 237   }
 238   return iResult;
 239 }
 240
 241 Int_t AliHLTTPCPad::CalculateBaseLine(Int_t reqMinCount)
 242 {
 243   // see header file for class documentation
 244   Int_t iResult=0;
 245   AliHLTTPCSignal_t avBackup=fAverage;
 246   //HLTDebug("reqMinCount=%d fCount=%d fTotal=%d fSum=%d fBLMax=%d fBLMin=%d", reqMinCount, fCount, fTotal, fSum, fBLMax, fBLMin);
 247   if (fCount>=reqMinCount && fCount>=fTotal/2) {
 248     fAverage=fCount>0?fSum/fCount:0;
 249     if (fAverage>0) {
 250       //HLTDebug("average for current event %d (%d - %d)", fAverage, fBLMax, fBLMin);
 251       fCount=0;fSum=-1;
 252       if (fBLMax>ALIHLTPAD_BASELINE_MARGIN) {
 253         // calculate again
 254         //HLTDebug("maximum value %d exceeds margin for base line (%d) "
 255         //       "-> re-evaluate base line", fBLMax, ALIHLTPAD_BASELINE_MARGIN);
 256         if (fpRawData) {
 257           for (Int_t i=fFirstBLBin; i<fNofBins; i++)
 258             if (fpRawData[i]>=0) AddBaseLineValue(i, fpRawData[i]);
 259           if (fCount>0 && fCount>=reqMinCount && fCount>=fTotal/2) {
 260             fAverage=fSum/fCount;
 261             //HLTDebug("new average %d", fAverage);
 262           } else {
 263             //      HLTDebug("baseline re-eveluation skipped because of to few "
 264             //                 "contributing bins: total=%d, contributing=%d, req=%d"
 265             //                 "\ndata might be already zero suppressed"
 266             //                 , fTotal, fCount, reqMinCount);
 267             iResult=-ENODATA;
 268           }
 269           fCount=0;fSum=-1;
 270         } else {
 271           HLTError("missing raw data for base line calculation");
 272           iResult=-ENOBUFS;
 273         }
 274       }
 275       if (iResult>=0) {
 276         // calculate average for all events
 277         fAverage=((avBackup*fNofEvents)+fAverage)/(fNofEvents+1);
 278         //HLTDebug("base line average for %d event(s): %d", fNofEvents+1, fAverage);
 279       } else {
 280         fAverage=avBackup;
 281       }
 282     } else {
 283       fAverage=avBackup;
 284     }
 285   } else {
 286     //     HLTDebug("baseline calculation skipped because of to few contributing "
 287     //         "bins: total=%d, contributing=%d, required=%d \ndata might be "
 288     //         "already zero suppressed", fTotal, fCount, reqMinCount);
 289   }
 290
 291   return iResult;
 292 }
 293
 294 Int_t AliHLTTPCPad::StopEvent()
 295 {
 296   // see header file for class documentation
 297   Int_t iResult=0;
 298   if (fpRawData) {
 299     AliHLTTPCSignal_t* pData=fpRawData;
 300     fpRawData=NULL;
 301     delete [] pData;
 302     fTotal=0;
 303     fNofEvents++;
 304     Rewind();
 305   } else if (fNofBins>0) {
 306     HLTError("event data acquisition not started");
 307     iResult=-EBADF;
 308   }
 309   return iResult;
 310 }
 311
 312 Int_t AliHLTTPCPad::ResetHistory()
 313 {
 314   // see header file for class documentation
 315   Int_t iResult=0;
 316   fAverage=-1;
 317   fNofEvents=0;
 318   return iResult;
 319 }
 320
 321 Int_t AliHLTTPCPad::SetThreshold(AliHLTTPCSignal_t thresh)
 322 {
 323   // see header file for class documentation
 324   Int_t iResult=0;
 325   fThreshold=thresh;
 326   return iResult;
 327 }
 328
 329 Int_t AliHLTTPCPad::AddBaseLineValue(Int_t bin, AliHLTTPCSignal_t value)
 330 {
 331   // see header file for class documentation
 332   Int_t iResult=0;
 333   if (bin>=fFirstBLBin) {
 334     if (fAverage<0 || value<ALIHLTPAD_BASELINE_MARGIN) {
 335       // add to the current sum and count
 336       fSum+=value;
 337       fCount++;
 338       if (fBLMax<value) {
 339         // keep the maximum value for later quality control of the base
 340         // line calculation
 341         fBLMax=value;
 342         fBLMaxBin=bin;
 343       }
 344       if (fBLMin<0 || fBLMin>value) {
 345         // keep the minimum value for later quality control of the base
 346         // line calculation
 347         fBLMin=value;
 348         fBLMinBin=bin;
 349       }
 350     } else {
 351       //       HLTDebug("ignoring value %d (bin %d) for base line calculation "
 352       //               "(current average is %d)",
 353       //               value, bin, fAverage);
 354     }
 355   }
 356   return iResult;
 357 }
 358
 359 Int_t AliHLTTPCPad::SetRawData(Int_t bin, AliHLTTPCSignal_t value)
 360 {
 361   // see header file for class documentation
 362   //  printf("Row: %d    Pad: %d  Time: %d Charge %d", fRowNo, fPadNo, bin, value);
 363   Int_t iResult=0;
 364   if (fpRawData) {
 365     if (bin<fNofBins) {
 366       if (value>=0) {
 367         if (fpRawData[bin]<0) {
 368           AddBaseLineValue(bin, value);
 369           fTotal++;
 370         } else {
 371           // ignore value for average calculation
 372           HLTWarning("overriding content of bin %d (%d)", bin, fpRawData[bin]);
 373         }
 374         fpRawData[bin]=value;
 375       } else {
 376         HLTWarning("ignoring neg. raw data");
 377       }
 378     } else {
 379       HLTWarning("bin %d out of range (%d)", bin, fNofBins);
 380       iResult=-ERANGE;
 381     }
 382   } else if (fNofBins>0) {
 383     HLTError("event cycle not started");
 384     iResult=-EBADF;
 385   }
 386   return iResult;
 387 }
 388
 389 Int_t AliHLTTPCPad::Next(Int_t bZeroSuppression)
 390 {
 391   // see header file for class documentation
 392   if (fpRawData==NULL) return 0;
 393   Int_t iResult=fReadPos<fNofBins;
 394   if (iResult>0 && (iResult=(++fReadPos<fNofBins))>0) {
 395     if (bZeroSuppression) {
 396       while ((iResult=(fReadPos<fNofBins))>0 &&
 397              GetCorrectedData(fReadPos)<=0)
 398         fReadPos++;
 399     }
 400   }
 401   return iResult;
 402 }
 403
 404 Int_t AliHLTTPCPad::Rewind(Int_t bZeroSuppression)
 405 {
 406   // see header file for class documentation
 407   fReadPos=(bZeroSuppression>0?0:fFirstBLBin)-1;
 408   return Next(bZeroSuppression);
 409 }
 410
 411 AliHLTTPCSignal_t AliHLTTPCPad::GetRawData(Int_t bin) const
 412 {
 413   // see header file for class documentation
 414   AliHLTTPCSignal_t data=0;
 415   if (fpRawData) {
 416     if (bin<fNofBins) {
 417       data=fpRawData[bin];
 418     } else {
 419       HLTWarning("requested bin %d out of range (%d)", bin, fNofBins);
 420     }
 421   } else if (fNofBins>0) {
 422     HLTWarning("data only available within event cycle");
 423   }
 424   return data;
 425 }
 426
 427 AliHLTTPCSignal_t AliHLTTPCPad::GetCorrectedData(Int_t bin) const
 428 {
 429   // see header file for class documentation
 430   AliHLTTPCSignal_t data=GetRawData(bin)-GetBaseLine(bin);
 431   AliHLTTPCSignal_t prev=0;
 432   if (bin>1) prev=GetRawData(bin-1)-GetBaseLine(bin-1);
 433   AliHLTTPCSignal_t succ=0;
 434   if (bin+1<GetSize()) succ=GetRawData(bin+1)-GetBaseLine(bin+1);
 435   if (fThreshold>0) {
 436     data-=fThreshold;
 437     prev-=fThreshold;
 438     succ-=fThreshold;
 439   }
 440
 441   // case 1:
 442   // the signal is below the base-line and threshold
 443   if (data<0) data=0;
 444
 445   //case 2:
 446   // the neighboring bins are both below base-line/threshold
 447   // a real signal is always more than one bin wide because of the shaper
 448   if (prev<=0 && succ<=0) data=0;
 449
 450   // case 3:
 451   // the bin is inside the range of ignored bins
 452   if (bin<fFirstBLBin) data=0;
 453   //HLTDebug("fReadPos=%d data=%d threshold=%d raw data=%d base line=%d", fReadPos, data, fThreshold, GetRawData(bin), GetBaseLine(bin));
 454   return data;
 455 }
 456
 457 AliHLTTPCSignal_t AliHLTTPCPad::GetBaseLine(Int_t /*bin*/) const //TODO: Why is bin being ignored?
 458 {
 459   // see header file for class documentation
 460   AliHLTTPCSignal_t val=0;
 461   if (fAverage>0) {
 462     // we take the minumum value as the base line if it doesn't differ from
 463     // the average to much
 464     val=fAverage;
 465 #ifdef KEEP_NOISE
 466     const AliHLTTPCSignal_t kMaxDifference=15;
 467     if ((fAverage-fBLMin)<=kMaxDifference) val=fBLMin;
 468     else val>kMaxDifference?val-=kMaxDifference:0;
 469 #endif
 470   }
 471   if (val<0) {
 472     // here we should never get
 473     val=0;
 474     HLTFatal("wrong base line value");
 475   }
 476   return val;
 477 }
 478
 479 AliHLTTPCSignal_t AliHLTTPCPad::GetAverage() const
 480 {
 481   // see header file for class documentation
 482   return fAverage>0?fAverage:0;
 483 }
 484
 485 Float_t AliHLTTPCPad::GetOccupancy() const
 486 {
 487   // see header file for class documentation
 488   Float_t occupancy=0;
 489   if (fpRawData && fNofBins>0) {
 490     for (Int_t i=fFirstBLBin; i<fNofBins; i++) {
 491       if (GetCorrectedData(i)>0) occupancy+=1;
 492     }
 493     if (fNofBins-fFirstBLBin>0)
 494       occupancy/=fNofBins-fFirstBLBin;
 495   }
 496   return occupancy;
 497 }
 498
 499 Float_t AliHLTTPCPad::GetAveragedOccupancy() const
 500 {
 501   // see header file for class documentation
 502
 503   // history is not yet implemented
 504   return GetOccupancy();
 505 }
 506 void AliHLTTPCPad::PrintRawData()
 507 {
 508   // see header file for class documentation
 509   for(Int_t bin=0;bin<AliHLTTPCTransform::GetNTimeBins();bin++){
 510     if(GetDataSignal(bin)>0)
 511       cout<<fRowNo<<"\t"<<fPadNo<<"\t"<<bin<<"\t"<<GetDataSignal(bin)<<endl;;
 512   }
 513   //  cout<<"bins: "<<AliHLTTPCTransform::GetNTimeBins()<<endl;
 514 }
 515
 516 void AliHLTTPCPad::ClearCandidates(){
 517   fClusterCandidates.clear();
 518   fUsedClusterCandidates.clear();
 519 }
 520
 521 void AliHLTTPCPad::SetDataToDefault()
 522 {
 523   // see header file for class documentation
 524   if(fDataSignals && fSignalPositionArray){
 525     for(Int_t i =0;i<fSizeOfSignalPositionArray;i++){
 526       fDataSignals[fSignalPositionArray[i]]=-1;
 527     }
 528     fSizeOfSignalPositionArray=0;
 529   }
 530 }
 531
 532 void AliHLTTPCPad::SetDataSignal(Int_t bin,Int_t signal)
 533 {
 534   // see header file for class documentation
 535   fDataSignals[bin]=signal;
 536   fSignalPositionArray[fSizeOfSignalPositionArray]=bin;
 537   fSizeOfSignalPositionArray++;
 538 #if DebugHisto
 539   fDebugHistoBeforeZS->Fill(bin,signal);
 540 #endif
 541 }
 542
 543 Bool_t AliHLTTPCPad::GetNextGoodSignal(Int_t &time, Int_t &signal ){
 544   /*  for(Int_t i=70;i<900;i++){
 545     if(fDataSignals[i]>0){
 546       printf("Signals which are good: Bin: %d Signal: %d\n",i,fDataSignals[i]);
 547     }
 548     }*/
 549   if(fNGoodSignalsSent<fSizeOfSignalPositionArray&&fSizeOfSignalPositionArray>0){
 550     time = fSignalPositionArray[fNGoodSignalsSent];
 551     signal = GetDataSignal(time);
 552     //    printf("GoodSignal: Row: %d Pad: %d time %d  signal %d  signalsSent: %d\n",fRowNo,fPadNo,fSignalPositionArray[fNGoodSignalsSent],GetDataSignal(time), fNGoodSignalsSent);
 553     fNGoodSignalsSent++;
 554     return kTRUE;
 555   }
 556   return kFALSE;
 557 }
 558
 559 Int_t AliHLTTPCPad::GetDataSignal(Int_t bin) const
 560 {
 561   // see header file for class documentation
 562   return fDataSignals[bin];
 563 }
 564
 565 void AliHLTTPCPad::ZeroSuppress(Double_t nRMS, Int_t threshold, Int_t reqMinPoint, Int_t beginTime, Int_t endTime, Int_t timebinsLeft, Int_t timebinsRight, Int_t valueUnderAverage){
 566   //see headerfile for documentation
 567
 568   //HLTDebug("In Pad: nRMS=%d, threshold=%d, reqMinPoint=%d, beginTime=%d, endTime=%d, timebinsLeft=%d timebinsRight=%d valueUnderAverage=%d \n",nRMS,threshold,reqMinPoint,beginTime,endTime,timebinsLeft,timebinsRight,valueUnderAverage);
 569
 570   Bool_t useRMS= kFALSE;
 571   if(nRMS>0){
 572     useRMS=kTRUE;
 573     if(threshold>0){
 574       HLTInfo("Both RMSThreshold and SignalThreshold defined, using RMSThreshold");
 575     }
 576   }
 577   if(threshold<1 && nRMS<=0){
 578     //setting the data to -1 for this pad
 579     HLTInfo("Neither of RMSThreshold and SignalThreshold set, zerosuppression aborted");
 580     return;
 581   }
 582
 583   Int_t fThresholdUsed=threshold;
 584
 585   Int_t nAdded=0;
 586   Int_t sumNAdded=0;
 587   fSizeOfSignalPositionArray=0;
 588   if(useRMS){
 589     for(Int_t i=beginTime;i<endTime+1;i++){
 590       if(fDataSignals[i]>0){
 591         nAdded++;
 592         sumNAdded+=fDataSignals[i]*fDataSignals[i];
 593       }
 594     }
 595   }
 596   else if(threshold>0){
 597     for(Int_t i=beginTime;i<endTime+1;i++){
 598       if(fDataSignals[i]>0){
 599         nAdded++;
 600         sumNAdded+=fDataSignals[i];
 601       }
 602     }
 603   }
 604   else{
 605     HLTFatal("This should never happen because this is tested earlier in the code.(nRMSThreshold<1&&signal-threshold<1)");
 606   }
 607   if(nAdded<reqMinPoint){
 608     HLTInfo("Number of signals is less than required, zero suppression aborted");
 609     return;
 610   }
 611
 612   if(nAdded==0){
 613     HLTInfo("No signals added for this pad, zerosuppression aborted: pad %d row %d",fPadNo,fRowNo);
 614     return;
 615   }
 616   // HLTInfo("sumNAdded=%d    nAdded=%d pad %d ",sumNAdded,nAdded,fPadNo);
 617   Double_t averageValue=(Double_t)sumNAdded/nAdded;//true average for threshold approach, average of signals squared for rms approach
 618
 619   //  Double_t rms=0;
 620   if(useRMS){
 621     //Calculate the RMS
 622     if(averageValue>0){
 623       fThresholdUsed =(Int_t)(TMath::Sqrt(averageValue)*nRMS);
 624     }
 625     else{
 626       HLTFatal("average value in ZeroSuppression less than 0, investigation needed. This should never happen");
 627     }
 628   }
 629   else{
 630     fThresholdUsed = (Int_t)(averageValue + threshold);
 631   }
 632
 633   averageValue = 55.3;
 634   // Do zero suppression on the adc values within [beginTime,endTime]
 635   for(Int_t i=beginTime;i<endTime;i++){
 636     if(fDataSignals[i]>fThresholdUsed){
 637       //  HLTInfo("Signal Larger in pad %d time %d signal %d  ,   threshold: %d  averageValue %e",fPadNo,i,fDataSignals[i],fThresholdUsed, averageValue);
 638       Int_t firstSignalTime=i;
 639       for(Int_t left=1;left<timebinsLeft;left++){//looking 5 to the left of the signal to add tail
 640         if(fDataSignals[i-left]-averageValue+valueUnderAverage>0&&i-left>=beginTime){
 641           firstSignalTime--;
 642         }
 643         else{
 644           break;
 645         }
 646       }
 647       Int_t lastSignalTime=i;
 648       for(Int_t right=1;right<timebinsRight;right++){//looking 5 to the left of the signal to add tail
 649         if(fDataSignals[i+right]-averageValue+valueUnderAverage>0&&i+right<endTime){
 650           lastSignalTime++;
 651         }
 652         else{
 653           break;
 654         }
 655       }
 656       for(Int_t t=firstSignalTime;t<lastSignalTime;t++){
 657         //      cout<<"Row: "<<fRowNo<<" Pad: "<<fPadNo<<"   Adding to tmebin: "<<t<<" signal: "<<(AliHLTTPCSignal_t)(fDataSignals[t]-averageValue + valueUnderAverage)<<endl;
 658         fDataSignals[t]=(AliHLTTPCSignal_t)(fDataSignals[t]-averageValue + valueUnderAverage);
 659         //      cout<<"Adding to signalPosition array bin number: "<<fSizeOfSignalPositionArray<<"    timebin number: "<<t<<endl;
 660         fSignalPositionArray[fSizeOfSignalPositionArray]=t;
 661         fSizeOfSignalPositionArray++;
 662         //      cout<<"Number of signals added so far: "<<fSizeOfSignalPositionArray<<"     firstSignalTimeBin: "<<firstSignalTime<<"     lastSignalTimeBin: "<<lastSignalTime<<endl;
 663         /*      if(fRowNo==29&&fPadNo==58){
 664           cout<<"Signal added: Row: "<<fRowNo<<" Pad: "<<fPadNo<<"  Time: "<<t<<" signal: "<<fDataSignals[t]<<"  #signals: "<<fSizeOfSignalPositionArray<<endl;
 665         }
 666         */
 667 #if DebugHisto
 668       fDebugHistoAfterZS->Fill(t,fDataSignals[t]);
 669 #endif
 670       }
 671       i+=lastSignalTime;
 672     }
 673     else{
 674       fDataSignals[i]=-1;
 675     }
 676   }
 677 }
 678
 679 void AliHLTTPCPad::AddClusterCandidate(AliHLTTPCClusters candidate){
 680   fClusterCandidates.push_back(candidate);
 681   fUsedClusterCandidates.push_back(0);
 682 }
 683
 684 void AliHLTTPCPad::SaveHistograms(){
 685 #if DebugHisto
 686   if(fSizeOfSignalPositionArray==0){
 687     return;
 688   }
 689   char* filename;
 690   sprintf(filename,"/afsuser/kenneth/SimpleComponentWrapper/histos/HistogramsRow%dPad%d.root",fRowNo,fPadNo);
 691   TFile file(filename,"RECREATE");
 692   fDebugHistoBeforeZS->Write();
 693   fDebugHistoAfterZS->Write();
 694   file.Close();
 695 #endif
 696 }
 697
 698 void AliHLTTPCPad::FindClusterCandidates()
 699 {
 700   // see header file for class documentation
 701   /*
 702   if(fSizeOfSignalPositionArray<2){
 703     return;
 704   }
 705
 706   if(fNSigmaThreshold>0){
 707     ZeroSuppress(fNSigmaThreshold);
 708   }
 709   else if(fSignalThreshold>0){
 710     ZeroSuppress((Double_t)0,(Int_t)fSignalThreshold);
 711   }
 712   UInt_t seqcharge=0;
 713   UInt_t seqaverage=0;
 714   UInt_t seqerror=0;
 715   vector<Int_t> tmpPos;
 716   vector<Int_t> tmpSig;
 717   UInt_t isFalling=0;
 718
 719   for(Int_t pos=fSizeOfSignalPositionArray-2;pos>=0;pos--){
 720     if(fSignalPositionArray[pos]==fSignalPositionArray[pos+1]+1){
 721       seqcharge+=fDataSignals[fSignalPositionArray[pos+1]];
 722       seqaverage += fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
 723       seqerror += fSignalPositionArray[pos+1]*fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
 724
 725       tmpPos.push_back(fSignalPositionArray[pos+1]);
 726       tmpSig.push_back(fDataSignals[fSignalPositionArray[pos+1]]);
 727
 728       if(fDataSignals[fSignalPositionArray[pos+1]]>fDataSignals[fSignalPositionArray[pos]]){
 729         isFalling=1;
 730       }
 731       if(fDataSignals[fSignalPositionArray[pos+1]]<fDataSignals[fSignalPositionArray[pos]]&&isFalling){
 732         Int_t seqmean=0;
 733         seqmean = seqaverage/seqcharge;
 734
 735         //Calculate mean in pad direction:
 736         Int_t padmean = seqcharge*fPadNo;
 737         Int_t paderror = fPadNo*padmean;
 738         AliHLTTPCClusters candidate;
 739         candidate.fTotalCharge   = seqcharge;
 740         candidate.fPad       = padmean;
 741         candidate.fPad2      = paderror;
 742         candidate.fTime      = seqaverage;
 743         candidate.fTime2     = seqerror;
 744         candidate.fMean          = seqmean;
 745         candidate.fLastMergedPad = fPadNo;
 746         fClusterCandidates.push_back(candidate);
 747         fUsedClusterCandidates.push_back(0);
 748         isFalling=0;
 749         seqcharge=0;
 750         seqaverage=0;
 751         seqerror=0;
 752
 753         tmpPos.clear();
 754         tmpSig.clear();
 755
 756         continue;
 757       }
 758
 759       if(pos<1){
 760         seqcharge+=fDataSignals[fSignalPositionArray[0]];
 761         seqaverage += fSignalPositionArray[0]*fDataSignals[fSignalPositionArray[0]];
 762         seqerror += fSignalPositionArray[0]*fSignalPositionArray[0]*fDataSignals[fSignalPositionArray[0]];
 763         tmpPos.push_back(fSignalPositionArray[0]);
 764         tmpSig.push_back(fDataSignals[fSignalPositionArray[0]]);
 765
 766         //Calculate mean of sequence:
 767         Int_t seqmean=0;
 768         seqmean = seqaverage/seqcharge;
 769
 770         //Calculate mean in pad direction:
 771         Int_t padmean = seqcharge*fPadNo;
 772         Int_t paderror = fPadNo*padmean;
 773         AliHLTTPCClusters candidate;
 774         candidate.fTotalCharge   = seqcharge;
 775         candidate.fPad       = padmean;
 776         candidate.fPad2      = paderror;
 777         candidate.fTime      = seqaverage;
 778         candidate.fTime2     = seqerror;
 779         candidate.fMean          = seqmean;
 780         candidate.fLastMergedPad = fPadNo;
 781         fClusterCandidates.push_back(candidate);
 782         fUsedClusterCandidates.push_back(0);
 783         isFalling=0;
 784         seqcharge=0;
 785         seqaverage=0;
 786         seqerror=0;
 787
 788         tmpPos.clear();
 789         tmpSig.clear();
 790       }
 791     }
 792     else if(seqcharge>0){
 793       seqcharge+=fDataSignals[fSignalPositionArray[pos+1]];
 794       seqaverage += fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
 795       seqerror += fSignalPositionArray[pos+1]*fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
 796       tmpPos.push_back(fSignalPositionArray[pos+1]);
 797       tmpSig.push_back(fDataSignals[fSignalPositionArray[pos+1]]);
 798
 799       //Calculate mean of sequence:
 800       Int_t seqmean=0;
 801       seqmean = seqaverage/seqcharge;
 802
 803       //Calculate mean in pad direction:
 804       Int_t padmean = seqcharge*fPadNo;
 805       Int_t paderror = fPadNo*padmean;
 806       AliHLTTPCClusters candidate;
 807       candidate.fTotalCharge   = seqcharge;
 808       candidate.fPad       = padmean;
 809       candidate.fPad2      = paderror;
 810       candidate.fTime      = seqaverage;
 811       candidate.fTime2     = seqerror;
 812       candidate.fMean          = seqmean;
 813       candidate.fLastMergedPad = fPadNo;
 814       fClusterCandidates.push_back(candidate);
 815       fUsedClusterCandidates.push_back(0);
 816       isFalling=0;
 817       seqcharge=0;
 818       seqaverage=0;
 819       seqerror=0;
 820
 821       tmpPos.clear();
 822       tmpSig.clear();
 823     }
 824   }
 825   */
 826 }
 827