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 #if __GNUC__>= 3
  27 using namespace std;
  28 #endif
  29
  30 #include <cerrno>
  31 #include "AliHLTTPCPad.h"
  32 #include "AliHLTStdIncludes.h"
  33
  34
  35 //added by kenneth
  36 #include "AliHLTTPCTransform.h"
  37 #include "AliHLTTPCClusters.h"
  38 #include <sys/time.h>
  39 //------------------------------
  40
  41 /** margin for the base line be re-avaluated */
  42 #define ALIHLTPAD_BASELINE_MARGIN (2*fAverage)
  43
  44 /** ROOT macro for the implementation of ROOT specific class methods */
  45 ClassImp(AliHLTTPCPad)
  46
  47 AliHLTTPCPad::AliHLTTPCPad()
  48   :
  49   fRowNo(-1),
  50   fPadNo(-1),
  51   fThreshold(0),
  52   fAverage(-1),
  53   fNofEvents(0),
  54   fSum(0),
  55   fCount(0),
  56   fTotal(0),
  57   fBLMax(-1),
  58   fBLMaxBin(-1),
  59   fBLMin(-1),
  60   fBLMinBin(-1),
  61   fFirstBLBin(0),
  62   fNofBins(0),
  63   fReadPos(0),
  64   fpRawData(NULL),
  65   fClusterCandidates(),
  66   fUsedClusterCandidates(0),
  67   fDataSignals(NULL),
  68   fSignalPositionArray(NULL),
  69   fSizeOfSignalPositionArray(0)
  70 {
  71   // see header file for class documentation
  72   // or
  73   // refer to README to build package
  74   // or
  75   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  76   //  HLTInfo("Entering default constructor");
  77   fDataSignals= new AliHLTTPCSignal_t[AliHLTTPCTransform::GetNTimeBins()];
  78   memset( fDataSignals, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
  79
  80   fSignalPositionArray= new AliHLTTPCSignal_t[AliHLTTPCTransform::GetNTimeBins()];
  81   memset( fSignalPositionArray, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
  82   fSizeOfSignalPositionArray=0;
  83 }
  84
  85 AliHLTTPCPad::AliHLTTPCPad(Int_t offset, Int_t nofBins)
  86   :
  87   fRowNo(-1),
  88   fPadNo(-1),
  89   fThreshold(0),
  90   fAverage(-1),
  91   fNofEvents(0),
  92   fSum(0),
  93   fCount(0),
  94   fTotal(0),
  95   fBLMax(-1),
  96   fBLMaxBin(-1),
  97   fBLMin(-1),
  98   fBLMinBin(-1),
  99   fFirstBLBin(offset),
 100   fNofBins(nofBins),
 101   fReadPos(0),
 102   fpRawData(NULL),
 103   fClusterCandidates(),
 104   fUsedClusterCandidates(0),
 105   fDataSignals(NULL),
 106   fSignalPositionArray(NULL),
 107   fSizeOfSignalPositionArray(0)
 108 {
 109   // see header file for class documentation
 110 }
 111
 112 AliHLTTPCPad::AliHLTTPCPad(const AliHLTTPCPad& srcPad)
 113   :
 114   fRowNo(srcPad.fRowNo),
 115   fPadNo(srcPad.fPadNo),
 116   fThreshold(0),
 117   fAverage(-1),
 118   fNofEvents(0),
 119   fSum(0),
 120   fCount(0),
 121   fTotal(0),
 122   fBLMax(-1),
 123   fBLMaxBin(-1),
 124   fBLMin(-1),
 125   fBLMinBin(-1),
 126   fFirstBLBin(0),
 127   fNofBins(0),
 128   fReadPos(0),
 129   fpRawData(NULL),
 130   fClusterCandidates(),
 131   fUsedClusterCandidates(0),
 132   fDataSignals(NULL),
 133   fSignalPositionArray(NULL),
 134   fSizeOfSignalPositionArray(0)
 135 {
 136   // see header file for class documentation
 137   HLTFatal("copy constructor not implemented");
 138 }
 139
 140 AliHLTTPCPad& AliHLTTPCPad::operator=(const AliHLTTPCPad&)
 141 {
 142   // see header file for class documentation
 143   HLTFatal("assignment operator not implemented");
 144   return (*this);
 145 }
 146
 147 AliHLTTPCPad::~AliHLTTPCPad()
 148 {
 149   // see header file for class documentation
 150   if (fpRawData) {
 151     HLTWarning("event data acquisition not stopped");
 152     StopEvent();
 153   }
 154   if (fDataSignals) {
 155     AliHLTTPCSignal_t* pData=fDataSignals;
 156     fDataSignals=NULL;
 157    delete [] pData;
 158   }
 159   if (fSignalPositionArray) {
 160     AliHLTTPCSignal_t* pData=fSignalPositionArray;
 161     fSignalPositionArray=NULL;
 162     //   delete [] pData;
 163   }
 164
 165 }
 166
 167 Int_t AliHLTTPCPad::SetID(Int_t rowno, Int_t padno)
 168 {
 169   // see header file for class documentation
 170   fRowNo=rowno;
 171   fPadNo=padno;
 172   return 0;
 173 }
 174
 175 Int_t AliHLTTPCPad::StartEvent()
 176 {
 177   // see header file for class documentation
 178   Int_t iResult=0;
 179   if (fpRawData==NULL) {
 180     fBLMax=-1;
 181     fBLMaxBin=-1;
 182     fBLMin=-1;
 183     fBLMinBin=-1;
 184     fSum=0;
 185     fCount=0;
 186     fTotal=0;
 187     if (fNofBins>0) {
 188       fpRawData=new AliHLTTPCSignal_t[fNofBins];
 189       if (fpRawData) {
 190         for (int i=0; i<fNofBins; i++) fpRawData[i]=-1;
 191       } else {
 192         HLTError("memory allocation failed");
 193         iResult=-ENOMEM;
 194       }
 195     }
 196   } else {
 197     HLTWarning("event data acquisition already started");
 198     iResult=-EALREADY;
 199   }
 200   return iResult;
 201 }
 202
 203 Int_t AliHLTTPCPad::CalculateBaseLine(Int_t reqMinCount)
 204 {
 205   // see header file for class documentation
 206   Int_t iResult=0;
 207   AliHLTTPCSignal_t avBackup=fAverage;
 208   //HLTDebug("reqMinCount=%d fCount=%d fTotal=%d fSum=%d fBLMax=%d fBLMin=%d", reqMinCount, fCount, fTotal, fSum, fBLMax, fBLMin);
 209   if (fCount>=reqMinCount && fCount>=fTotal/2) {
 210     fAverage=fCount>0?fSum/fCount:0;
 211     if (fAverage>0) {
 212       //HLTDebug("average for current event %d (%d - %d)", fAverage, fBLMax, fBLMin);
 213       fCount=0;fSum=-1;
 214       if (fBLMax>ALIHLTPAD_BASELINE_MARGIN) {
 215         // calculate again
 216         //HLTDebug("maximum value %d exceeds margin for base line (%d) "
 217         //       "-> re-evaluate base line", fBLMax, ALIHLTPAD_BASELINE_MARGIN);
 218         if (fpRawData) {
 219           for (Int_t i=fFirstBLBin; i<fNofBins; i++)
 220             if (fpRawData[i]>=0) AddBaseLineValue(i, fpRawData[i]);
 221           if (fCount>0 && fCount>=reqMinCount && fCount>=fTotal/2) {
 222             fAverage=fSum/fCount;
 223             //HLTDebug("new average %d", fAverage);
 224           } else {
 225 //          HLTDebug("baseline re-eveluation skipped because of to few "
 226 //                     "contributing bins: total=%d, contributing=%d, req=%d"
 227 //                     "\ndata might be already zero suppressed"
 228 //                     , fTotal, fCount, reqMinCount);
 229             iResult=-ENODATA;
 230           }
 231           fCount=0;fSum=-1;
 232         } else {
 233           HLTError("missing raw data for base line calculation");
 234           iResult=-ENOBUFS;
 235         }
 236       }
 237       if (iResult>=0) {
 238         // calculate average for all events
 239         fAverage=((avBackup*fNofEvents)+fAverage)/(fNofEvents+1);
 240         //HLTDebug("base line average for %d event(s): %d", fNofEvents+1, fAverage);
 241       } else {
 242         fAverage=avBackup;
 243       }
 244     } else {
 245       fAverage=avBackup;
 246     }
 247   } else {
 248 //     HLTDebug("baseline calculation skipped because of to few contributing "
 249 //             "bins: total=%d, contributing=%d, required=%d \ndata might be "
 250 //             "already zero suppressed", fTotal, fCount, reqMinCount);
 251   }
 252
 253   return iResult;
 254 }
 255
 256 Int_t AliHLTTPCPad::StopEvent()
 257 {
 258   // see header file for class documentation
 259   Int_t iResult=0;
 260   if (fpRawData) {
 261     AliHLTTPCSignal_t* pData=fpRawData;
 262     fpRawData=NULL;
 263     delete [] pData;
 264     fTotal=0;
 265     fNofEvents++;
 266     Rewind();
 267   } else if (fNofBins>0) {
 268     HLTError("event data acquisition not started");
 269     iResult=-EBADF;
 270   }
 271   return iResult;
 272 }
 273
 274 Int_t AliHLTTPCPad::ResetHistory()
 275 {
 276   // see header file for class documentation
 277   Int_t iResult=0;
 278   fAverage=-1;
 279   fNofEvents=0;
 280   return iResult;
 281 }
 282
 283 Int_t AliHLTTPCPad::SetThreshold(AliHLTTPCSignal_t thresh)
 284 {
 285   // see header file for class documentation
 286   Int_t iResult=0;
 287   fThreshold=thresh;
 288   return iResult;
 289 }
 290
 291 Int_t AliHLTTPCPad::AddBaseLineValue(Int_t bin, AliHLTTPCSignal_t value)
 292 {
 293   // see header file for class documentation
 294   Int_t iResult=0;
 295   if (bin>=fFirstBLBin) {
 296     if (fAverage<0 || value<ALIHLTPAD_BASELINE_MARGIN) {
 297       // add to the current sum and count
 298       fSum+=value;
 299       fCount++;
 300       if (fBLMax<value) {
 301         // keep the maximum value for later quality control of the base
 302         // line calculation
 303         fBLMax=value;
 304         fBLMaxBin=bin;
 305       }
 306       if (fBLMin<0 || fBLMin>value) {
 307         // keep the minimum value for later quality control of the base
 308         // line calculation
 309         fBLMin=value;
 310         fBLMinBin=bin;
 311       }
 312     } else {
 313 //       HLTDebug("ignoring value %d (bin %d) for base line calculation "
 314 //             "(current average is %d)",
 315 //             value, bin, fAverage);
 316     }
 317   }
 318   return iResult;
 319 }
 320
 321 Int_t AliHLTTPCPad::SetRawData(Int_t bin, AliHLTTPCSignal_t value)
 322 {
 323   // see header file for class documentation
 324   Int_t iResult=0;
 325   if (fpRawData) {
 326     if (bin<fNofBins) {
 327       if (value>=0) {
 328         if (fpRawData[bin]<0) {
 329           AddBaseLineValue(bin, value);
 330           fTotal++;
 331         } else {
 332           // ignore value for average calculation
 333           HLTWarning("overriding content of bin %d (%d)", bin, fpRawData[bin]);
 334         }
 335         fpRawData[bin]=value;
 336       } else {
 337         HLTWarning("ignoring neg. raw data");
 338       }
 339     } else {
 340       HLTWarning("bin %d out of range (%d)", bin, fNofBins);
 341       iResult=-ERANGE;
 342     }
 343   } else if (fNofBins>0) {
 344     HLTError("event cycle not started");
 345     iResult=-EBADF;
 346   }
 347   return iResult;
 348 }
 349
 350 Int_t AliHLTTPCPad::Next(Int_t bZeroSuppression)
 351 {
 352   // see header file for class documentation
 353   if (fpRawData==NULL) return 0;
 354   Int_t iResult=fReadPos<fNofBins;
 355   if (iResult>0 && (iResult=(++fReadPos<fNofBins))>0) {
 356     if (bZeroSuppression) {
 357       while ((iResult=(fReadPos<fNofBins))>0 &&
 358              GetCorrectedData(fReadPos)<=0)
 359         fReadPos++;
 360     }
 361   }
 362   return iResult;
 363 }
 364
 365 Int_t AliHLTTPCPad::Rewind(Int_t bZeroSuppression)
 366 {
 367   // see header file for class documentation
 368   fReadPos=(bZeroSuppression>0?0:fFirstBLBin)-1;
 369   return Next(bZeroSuppression);
 370 }
 371
 372 AliHLTTPCSignal_t AliHLTTPCPad::GetRawData(Int_t bin) const
 373 {
 374   // see header file for class documentation
 375   AliHLTTPCSignal_t data=0;
 376   if (fpRawData) {
 377     if (bin<fNofBins) {
 378       data=fpRawData[bin];
 379     } else {
 380       HLTWarning("requested bin %d out of range (%d)", bin, fNofBins);
 381     }
 382   } else if (fNofBins>0) {
 383     HLTWarning("data only available within event cycle");
 384   }
 385   return data;
 386 }
 387
 388 AliHLTTPCSignal_t AliHLTTPCPad::GetCorrectedData(Int_t bin) const
 389 {
 390   // see header file for class documentation
 391   AliHLTTPCSignal_t data=GetRawData(bin)-GetBaseLine(bin);
 392   AliHLTTPCSignal_t prev=0;
 393   if (bin>1) prev=GetRawData(bin-1)-GetBaseLine(bin-1);
 394   AliHLTTPCSignal_t succ=0;
 395   if (bin+1<GetSize()) succ=GetRawData(bin+1)-GetBaseLine(bin+1);
 396   if (fThreshold>0) {
 397     data-=fThreshold;
 398     prev-=fThreshold;
 399     succ-=fThreshold;
 400   }
 401
 402   // case 1:
 403   // the signal is below the base-line and threshold
 404   if (data<0) data=0;
 405
 406   //case 2:
 407   // the neighboring bins are both below base-line/threshold
 408   // a real signal is always more than one bin wide because of the shaper
 409   if (prev<=0 && succ<=0) data=0;
 410
 411   // case 3:
 412   // the bin is inside the range of ignored bins
 413   if (bin<fFirstBLBin) data=0;
 414   //HLTDebug("fReadPos=%d data=%d threshold=%d raw data=%d base line=%d", fReadPos, data, fThreshold, GetRawData(bin), GetBaseLine(bin));
 415   return data;
 416 }
 417
 418 AliHLTTPCSignal_t AliHLTTPCPad::GetBaseLine(Int_t bin) const
 419 {
 420   // see header file for class documentation
 421   AliHLTTPCSignal_t val=0;
 422   if (fAverage>0) {
 423     // we take the minumum value as the base line if it doesn't differ from
 424     // the average to much
 425     val=fAverage;
 426 #ifdef KEEP_NOISE
 427     const AliHLTTPCSignal_t kMaxDifference=15;
 428     if ((fAverage-fBLMin)<=kMaxDifference) val=fBLMin;
 429     else val>kMaxDifference?val-=kMaxDifference:0;
 430 #endif
 431   }
 432   if (val<0) {
 433     // here we should never get
 434     val=0;
 435     HLTFatal("wrong base line value");
 436   }
 437   return val;
 438 }
 439
 440 AliHLTTPCSignal_t AliHLTTPCPad::GetAverage() const
 441 {
 442   // see header file for class documentation
 443   return fAverage>0?fAverage:0;
 444 }
 445
 446 Float_t AliHLTTPCPad::GetOccupancy() const
 447 {
 448   // see header file for class documentation
 449   Float_t occupancy=0;
 450   if (fpRawData && fNofBins>0) {
 451     for (Int_t i=fFirstBLBin; i<fNofBins; i++) {
 452       if (GetCorrectedData(i)>0) occupancy+=1;
 453     }
 454     if (fNofBins-fFirstBLBin>0)
 455       occupancy/=fNofBins-fFirstBLBin;
 456   }
 457   return occupancy;
 458 }
 459
 460 Float_t AliHLTTPCPad::GetAveragedOccupancy() const
 461 {
 462   // see header file for class documentation
 463
 464   // history is not yet implemented
 465   return GetOccupancy();
 466 }
 467 void AliHLTTPCPad::PrintRawData(){
 468   for(Int_t bin=0;bin<AliHLTTPCTransform::GetNTimeBins();bin++){
 469     if(GetDataSignal(bin)>0)
 470       cout<<fRowNo<<"\t"<<fPadNo<<"\t"<<bin<<"\t"<<GetDataSignal(bin)<<endl;;
 471   }
 472   cout<<"bins: "<<AliHLTTPCTransform::GetNTimeBins()<<endl;
 473 }
 474
 475 void AliHLTTPCPad::SetDataToDefault(){
 476   if(fpRawData){
 477     memset( fDataSignals, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
 478     memset( fSignalPositionArray, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
 479     fSizeOfSignalPositionArray=0;
 480   }
 481 }
 482 void AliHLTTPCPad::SetDataSignal(Int_t bin,Int_t signal){
 483   fDataSignals[bin]=signal;
 484   fSignalPositionArray[fSizeOfSignalPositionArray]=bin;
 485   fSizeOfSignalPositionArray++;
 486 }
 487 Int_t AliHLTTPCPad::GetDataSignal(Int_t bin){
 488   return fDataSignals[bin];
 489 }
 490 void AliHLTTPCPad::FindClusterCandidates(){
 491   UInt_t seqcharge=0;
 492   UInt_t seqaverage=0;
 493   UInt_t seqerror=0;
 494   vector<Int_t> tmpPos;
 495   vector<Int_t> tmpSig;
 496   UInt_t isFalling=0;
 497
 498   for(Int_t pos=fSizeOfSignalPositionArray-2;pos>=0;pos--){
 499
 500     if(fSignalPositionArray[pos]==fSignalPositionArray[pos+1]+1){
 501       seqcharge+=fDataSignals[fSignalPositionArray[pos+1]];
 502       seqaverage += fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
 503       seqerror += fSignalPositionArray[pos+1]*fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
 504
 505       tmpPos.push_back(fSignalPositionArray[pos+1]);
 506       tmpSig.push_back(fDataSignals[fSignalPositionArray[pos+1]]);
 507
 508       if(fDataSignals[fSignalPositionArray[pos+1]]>fDataSignals[fSignalPositionArray[pos]]){
 509         isFalling=1;
 510       }
 511       if(fDataSignals[fSignalPositionArray[pos+1]]<fDataSignals[fSignalPositionArray[pos]]&&isFalling){
 512         Int_t seqmean=0;
 513         seqmean = seqaverage/seqcharge;
 514
 515         //Calculate mean in pad direction:
 516         Int_t padmean = seqcharge*fPadNo;
 517         Int_t paderror = fPadNo*padmean;
 518         AliHLTTPCClusters candidate;
 519         candidate.fTotalCharge   = seqcharge;
 520         candidate.fPad       = padmean;
 521         candidate.fPad2      = paderror;
 522         candidate.fTime      = seqaverage;
 523         candidate.fTime2     = seqerror;
 524         candidate.fMean          = seqmean;
 525         candidate.fLastMergedPad = fPadNo;
 526         fClusterCandidates.push_back(candidate);
 527         fUsedClusterCandidates.push_back(0);
 528         isFalling=0;
 529         seqcharge=0;
 530         seqaverage=0;
 531         seqerror=0;
 532
 533         tmpPos.clear();
 534         tmpSig.clear();
 535
 536         continue;
 537       }
 538
 539       if(pos<1){
 540         seqcharge+=fDataSignals[fSignalPositionArray[0]];
 541         seqaverage += fSignalPositionArray[0]*fDataSignals[fSignalPositionArray[0]];
 542         seqerror += fSignalPositionArray[0]*fSignalPositionArray[0]*fDataSignals[fSignalPositionArray[0]];
 543         tmpPos.push_back(fSignalPositionArray[0]);
 544         tmpSig.push_back(fDataSignals[fSignalPositionArray[0]]);
 545
 546         //Calculate mean of sequence:
 547         Int_t seqmean=0;
 548         seqmean = seqaverage/seqcharge;
 549
 550         //Calculate mean in pad direction:
 551         Int_t padmean = seqcharge*fPadNo;
 552         Int_t paderror = fPadNo*padmean;
 553         AliHLTTPCClusters candidate;
 554         candidate.fTotalCharge   = seqcharge;
 555         candidate.fPad       = padmean;
 556         candidate.fPad2      = paderror;
 557         candidate.fTime      = seqaverage;
 558         candidate.fTime2     = seqerror;
 559         candidate.fMean          = seqmean;
 560         candidate.fLastMergedPad = fPadNo;
 561         fClusterCandidates.push_back(candidate);
 562         fUsedClusterCandidates.push_back(0);
 563         isFalling=0;
 564         seqcharge=0;
 565         seqaverage=0;
 566         seqerror=0;
 567
 568         tmpPos.clear();
 569         tmpSig.clear();
 570       }
 571     }
 572     else if(seqcharge>0){
 573       seqcharge+=fDataSignals[fSignalPositionArray[pos+1]];
 574       seqaverage += fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
 575         seqerror += fSignalPositionArray[pos+1]*fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
 576         tmpPos.push_back(fSignalPositionArray[pos+1]);
 577         tmpSig.push_back(fDataSignals[fSignalPositionArray[pos+1]]);
 578
 579         //Calculate mean of sequence:
 580         Int_t seqmean=0;
 581         seqmean = seqaverage/seqcharge;
 582
 583         //Calculate mean in pad direction:
 584         Int_t padmean = seqcharge*fPadNo;
 585         Int_t paderror = fPadNo*padmean;
 586         AliHLTTPCClusters candidate;
 587         candidate.fTotalCharge   = seqcharge;
 588         candidate.fPad       = padmean;
 589         candidate.fPad2      = paderror;
 590         candidate.fTime      = seqaverage;
 591         candidate.fTime2     = seqerror;
 592         candidate.fMean          = seqmean;
 593         candidate.fLastMergedPad = fPadNo;
 594         fClusterCandidates.push_back(candidate);
 595         fUsedClusterCandidates.push_back(0);
 596         isFalling=0;
 597         seqcharge=0;
 598         seqaverage=0;
 599         seqerror=0;
 600
 601         tmpPos.clear();
 602         tmpSig.clear();
 603     }
 604   }
 605 }