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