-// @(#) $Id$
-
-/**************************************************************************
- * This file is property of and copyright by the ALICE HLT Project *
- * ALICE Experiment at CERN, All rights reserved. *
- * *
- * Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no> *
- * Kenneth Aamodt <Kenneth.aamodt@ift.uib.no> *
- * for The ALICE HLT Project. *
- * *
- * Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes is hereby granted *
- * without fee, provided that the above copyright notice appears in all *
- * copies and that both the copyright notice and this permission notice *
- * appear in the supporting documentation. The authors make no claims *
- * about the suitability of this software for any purpose. It is *
- * provided "as is" without express or implied warranty. *
- **************************************************************************/
-
-/** @file AliHLTTPCPad.cxx
- @author Matthias Richter, Kenneth Aamodt
- @date
- @brief Container Class for TPC Pads.
-*/
-
-// see header file for class documentation
-// or
-// refer to README to build package
-// or
-// visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
-
-#if __GNUC__>= 3
-using namespace std;
-#endif
+// $Id$
+
+//**************************************************************************
+//* This file is property of and copyright by the ALICE HLT Project *
+//* ALICE Experiment at CERN, All rights reserved. *
+//* *
+//* Primary Authors: Timm Steinbeck, Matthias Richter *
+//* Developers: Kenneth Aamodt <kenneth.aamodt@student.uib.no> *
+//* for The ALICE HLT Project. *
+//* *
+//* Permission to use, copy, modify and distribute this software and its *
+//* documentation strictly for non-commercial purposes is hereby granted *
+//* without fee, provided that the above copyright notice appears in all *
+//* copies and that both the copyright notice and this permission notice *
+//* appear in the supporting documentation. The authors make no claims *
+//* about the suitability of this software for any purpose. It is *
+//* provided "as is" without express or implied warranty. *
+//**************************************************************************
+
+/// @file AliHLTTPCPad.cxx
+/// @author Matthias Richter, Kenneth Aamodt
+/// @date
+/// @brief Container Class for TPC Pads.
+///
#include <cerrno>
#include "AliHLTTPCPad.h"
#include "AliHLTTPCTransform.h"
#include "AliHLTTPCClusters.h"
#include <sys/time.h>
+#include "TMath.h"
+#include "TFile.h"
//------------------------------
/** margin for the base line be re-avaluated */
:
fClusterCandidates(),
fUsedClusterCandidates(),
+ fSelectedPad(kFALSE),
+ fHWAddress(0),
fRowNo(-1),
fPadNo(-1),
fThreshold(0),
fDataSignals(NULL),
fSignalPositionArray(NULL),
fSizeOfSignalPositionArray(0),
- fNSigmaThreshold(0),
- fSignalThreshold(0)
+ fNGoodSignalsSent(0),
+ fCandidateDigitsVector()
{
// see header file for class documentation
// or
// HLTInfo("Entering default constructor");
fDataSignals= new AliHLTTPCSignal_t[AliHLTTPCTransform::GetNTimeBins()];
memset( fDataSignals, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
-
- fSignalPositionArray= new AliHLTTPCSignal_t[AliHLTTPCTransform::GetNTimeBins()];
+
+ fSignalPositionArray= new Int_t[AliHLTTPCTransform::GetNTimeBins()];
memset( fSignalPositionArray, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
fSizeOfSignalPositionArray=0;
+
+}
+
+AliHLTTPCPad::AliHLTTPCPad(Int_t /*dummy*/)
+ :
+ fClusterCandidates(),
+ fUsedClusterCandidates(),
+ fSelectedPad(kFALSE),
+ fHWAddress(0),
+ fRowNo(-1),
+ fPadNo(-1),
+ fThreshold(0),
+ fAverage(-1),
+ fNofEvents(0),
+ fSum(0),
+ fCount(0),
+ fTotal(0),
+ fBLMax(-1),
+ fBLMaxBin(-1),
+ fBLMin(-1),
+ fBLMinBin(-1),
+ fFirstBLBin(0),
+ fNofBins(0),
+ fReadPos(0),
+ fpRawData(NULL),
+ fDataSignals(NULL),
+ fSignalPositionArray(NULL),
+ fSizeOfSignalPositionArray(0),
+ fNGoodSignalsSent(0),
+ fCandidateDigitsVector()
+{
+ // see header file for class documentation
+ // or
+ // refer to README to build package
+ // or
+ // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
}
AliHLTTPCPad::AliHLTTPCPad(Int_t offset, Int_t nofBins)
:
fClusterCandidates(),
fUsedClusterCandidates(),
+ fSelectedPad(kFALSE),
+ fHWAddress(0),
fRowNo(-1),
fPadNo(-1),
fThreshold(0),
fDataSignals(NULL),
fSignalPositionArray(NULL),
fSizeOfSignalPositionArray(0),
- fNSigmaThreshold(0),
- fSignalThreshold(0)
+ fNGoodSignalsSent(0),
+ fCandidateDigitsVector()
{
// see header file for class documentation
}
StopEvent();
}
if (fDataSignals) {
- AliHLTTPCSignal_t* pData=fDataSignals;
+ delete [] fDataSignals;
fDataSignals=NULL;
- delete [] pData;
}
- if (fSignalPositionArray) {
- //AliHLTTPCSignal_t* pData=fSignalPositionArray;
+ if (fSignalPositionArray!=NULL) {
+ delete [] fSignalPositionArray;
fSignalPositionArray=NULL;
- // delete [] pData;
}
-
}
Int_t AliHLTTPCPad::SetID(Int_t rowno, Int_t padno)
// see header file for class documentation
fRowNo=rowno;
fPadNo=padno;
+
return 0;
}
Int_t AliHLTTPCPad::SetRawData(Int_t bin, AliHLTTPCSignal_t value)
{
// see header file for class documentation
+ // printf("Row: %d Pad: %d Time: %d Charge %d", fRowNo, fPadNo, bin, value);
Int_t iResult=0;
if (fpRawData) {
if (bin<fNofBins) {
// see header file for class documentation
for(Int_t bin=0;bin<AliHLTTPCTransform::GetNTimeBins();bin++){
if(GetDataSignal(bin)>0)
- cout<<fRowNo<<"\t"<<fPadNo<<"\t"<<bin<<"\t"<<GetDataSignal(bin)<<endl;;
+ //This cout should be here since using logging produces output that is much more difficult to read
+ cout<<fRowNo<<"\t"<<fPadNo<<"\t"<<bin<<"\t"<<GetDataSignal(bin)<<endl;
}
- cout<<"bins: "<<AliHLTTPCTransform::GetNTimeBins()<<endl;
+}
+
+void AliHLTTPCPad::ClearCandidates(){
+ fClusterCandidates.clear();
+ fUsedClusterCandidates.clear();
+ fCandidateDigitsVector.clear();
}
void AliHLTTPCPad::SetDataToDefault()
{
// see header file for class documentation
- if(fpRawData){
- memset( fDataSignals, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
- memset( fSignalPositionArray, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
+ // if(fDataSignals && fSignalPositionArray){
+ for(Int_t i =0;i<fSizeOfSignalPositionArray;i++){
+ fDataSignals[fSignalPositionArray[i]]=-1;
+ }
fSizeOfSignalPositionArray=0;
- }
+ fNGoodSignalsSent = 0;
+ // }
}
void AliHLTTPCPad::SetDataSignal(Int_t bin,Int_t signal)
fSizeOfSignalPositionArray++;
}
+Bool_t AliHLTTPCPad::GetNextGoodSignal(Int_t &time,Int_t &bunchSize){
+
+ if(fNGoodSignalsSent<fSizeOfSignalPositionArray&&fSizeOfSignalPositionArray>0){
+ time = fSignalPositionArray[fNGoodSignalsSent];
+ bunchSize=1;
+ fNGoodSignalsSent++;
+ while(fNGoodSignalsSent<fSizeOfSignalPositionArray){
+ if(fDataSignals[time+bunchSize+1]>0){
+ bunchSize++;
+ fNGoodSignalsSent++;
+ }
+ else{
+ break;
+ }
+ }
+ // fNGoodSignalsSent++;
+ return kTRUE;
+ }
+ return kFALSE;
+}
+
Int_t AliHLTTPCPad::GetDataSignal(Int_t bin) const
{
// see header file for class documentation
return fDataSignals[bin];
}
-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){
+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, bool speedup){
//see headerfile for documentation
- Bool_t useSigma= kFALSE;
- if(nSigma>0){
- useSigma=kTRUE;
+ //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);
+
+ Bool_t useRMS= kFALSE;
+ if(nRMS>0){
+ useRMS=kTRUE;
+ if(threshold>0){
+ HLTInfo("Both RMSThreshold and SignalThreshold defined, using RMSThreshold");
+ }
}
- if(threshold<1 && nSigma<=0){
+ if(threshold<1 && nRMS<=0){
//setting the data to -1 for this pad
- memset( fDataSignals, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
- fSizeOfSignalPositionArray=0;
+ HLTInfo("Neither of RMSThreshold and SignalThreshold set, zerosuppression aborted");
return;
}
- if(endTime>=AliHLTTPCTransform::GetNTimeBins()){
- endTime=AliHLTTPCTransform::GetNTimeBins()-1;
- }
Int_t fThresholdUsed=threshold;
-
+
+ Int_t maxVal=0;
Int_t nAdded=0;
Int_t sumNAdded=0;
fSizeOfSignalPositionArray=0;
- for(Int_t i=beginTime;i<endTime+1;i++){
- if(fDataSignals[i]>0){
- nAdded++;
- sumNAdded+=fDataSignals[i];
+ if(useRMS){
+ for(Int_t i=beginTime;i<endTime+1;i++){
+ if(fDataSignals[i]>0){
+ nAdded++;
+ sumNAdded+=fDataSignals[i]*fDataSignals[i];
+ if (maxVal<fDataSignals[i]) maxVal=fDataSignals[i];
+ }
}
}
-
- if(nAdded<reqMinPoint){
- return; //This will ensure that no data is read in FindClusterCandidates() (since fSizeOfSignalPositionArray=0)
- }
-
- Double_t averageValue=sumNAdded/nAdded;
-
- Double_t sigma=0;
- if(useSigma){
- //Calculate the sigma
- Double_t sumOfDifferenceSquared=0;
- for(Int_t i=endTime;i>=beginTime;i--){
+ else if(threshold>0){
+ for(Int_t i=beginTime;i<endTime+1;i++){
if(fDataSignals[i]>0){
- if(fDataSignals[i]-averageValue<50){
- sumOfDifferenceSquared+=(fDataSignals[i]-averageValue)*(fDataSignals[i]-averageValue);
- }
- else{
- nAdded--;
- }
+ nAdded++;
+ sumNAdded+=fDataSignals[i];
+ if (maxVal<fDataSignals[i]) maxVal=fDataSignals[i];
}
}
- sigma=sumOfDifferenceSquared/nAdded;
- fThresholdUsed=(Int_t)(nSigma*sigma);
}
-
- //For now just set the adc value outside [beginTime,endTime] to -1
- for(Int_t i=0;i<beginTime;i++){
- fDataSignals[i]=-1;
+ else{
+ HLTFatal("This should never happen because this is tested earlier in the code.(nRMSThreshold<1&&signal-threshold<1)");
+ }
+ if(nAdded<reqMinPoint){
+ HLTInfo("Number of signals is less than required, zero suppression aborted");
+ return;
}
- for(Int_t i=endTime+1;i<AliHLTTPCTransform::GetNTimeBins();i++){
- fDataSignals[i]=-1;
+
+ if(nAdded==0){
+ HLTInfo("No signals added for this pad, zerosuppression aborted: pad %d row %d",fPadNo,fRowNo);
+ return;
}
+
+ Double_t averageValue=(Double_t)sumNAdded/nAdded;//true average for threshold approach, average of signals squared for rms approach
- // Do zero suppression on the adc values within [beginTime,endTime]
- for(Int_t i=endTime;i>=beginTime;i--){
- //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
- //(better to set it to -1 which is default for no signal)
- if(fDataSignals[i]>(Int_t)(averageValue+fThresholdUsed+1) && fDataSignals[i-1]>(Int_t)(averageValue+fThresholdUsed+1)){
- //here the signals below threshold to the right of the candidate is added
- Bool_t contRight=kTRUE;
- Int_t endRight=i;
- Int_t nToAddRight=0;
- while(contRight){
- if(endRight+1<endTime){
- // cout<<fDataSignals[endRight+1]<<" "<<fDataSignals[endRight+2]<<endl;;
- if(fDataSignals[endRight+1]>=fDataSignals[endRight+2] && fDataSignals[endRight+1]>averageValue){
- nToAddRight++;
- }
- else{
- if(fDataSignals[endRight+1]> averageValue){
- nToAddRight++;
- }
- contRight=kFALSE;
- }
+ if(useRMS){
+ //Calculate the RMS
+ if(averageValue>0){
+ fThresholdUsed =(Int_t)(TMath::Sqrt(averageValue)*nRMS);
+ }
+ else{
+ HLTFatal("average value in ZeroSuppression less than 0, investigation needed. This should never happen");
+ }
+ }
+ else{
+ fThresholdUsed = (Int_t)(averageValue + threshold);
+ }
+ if (maxVal<fThresholdUsed) return;
+
+ // Do zero suppression on the adc values within [beginTime,endTime](add the good values)
+ for(Int_t i=beginTime;i<endTime;i++){
+ if(fDataSignals[i]>fThresholdUsed){
+ Int_t firstSignalTime=i;
+ for(Int_t left=1;left<timebinsLeft;left++){//looking 5 to the left of the signal to add tail
+ if(fDataSignals[i-left]-averageValue+valueUnderAverage>0 && i-left>=beginTime){
+ firstSignalTime--;
}
- else if(endRight>endTime+1){
- contRight=kFALSE;
+ else{
+ break;
}
- endRight++;
}
- for(int j=i+nToAddRight;j>i;j--){
- fDataSignals[j]=(Int_t)(fDataSignals[j]-averageValue);
- fSignalPositionArray[fSizeOfSignalPositionArray]=j;
- fSizeOfSignalPositionArray++;
+ Int_t lastSignalTime=i;
+ while(fDataSignals[lastSignalTime+1]>fThresholdUsed && lastSignalTime+1<endTime){
+ lastSignalTime++;
}
-
-
- //before while the two consecutive timebin values are added
- fDataSignals[i]=(Int_t)(fDataSignals[i]-averageValue);
- fSignalPositionArray[fSizeOfSignalPositionArray]=i;
- fSizeOfSignalPositionArray++;
- fDataSignals[i-1]=(Int_t)(fDataSignals[i-1]-averageValue);
- fSignalPositionArray[fSizeOfSignalPositionArray]=i-1;
- fSizeOfSignalPositionArray++;
- i--;
- // cout<<""<<endl;
- //Here the consecutive pads after the two first are added
- if(i-1>0){
- while(fDataSignals[i-1]>(Int_t)(averageValue+fThresholdUsed+1)){
- fDataSignals[i-1]=(Int_t)(fDataSignals[i-1]-averageValue);
- fSignalPositionArray[fSizeOfSignalPositionArray]=i-1;
- fSizeOfSignalPositionArray++;
- i--;
- }
- }
- //adding the signal below threshold belonging to the total signal
- Bool_t contLeft=kTRUE;
- while(contLeft){
- if(i-2>0){
- if(fDataSignals[i-1]>=fDataSignals[i-2] && fDataSignals[i-1]>averageValue){
- fDataSignals[i-1]=(Int_t)(fDataSignals[i-1]-averageValue);
- fSignalPositionArray[fSizeOfSignalPositionArray]=i-1;
- fSizeOfSignalPositionArray++;
- i--;
- }
- else{
- if(fDataSignals[i-1]> averageValue){
- fDataSignals[i-1]=(Int_t)(fDataSignals[i-1]-averageValue);
- fSignalPositionArray[fSizeOfSignalPositionArray]=i-1;
- fSizeOfSignalPositionArray++;
- i--;
- }
- contLeft=kFALSE;
- }
+ for(Int_t right=1;right<timebinsRight;right++){//looking 5 to the left of the signal to add tail
+ if(fDataSignals[i+right]-averageValue+valueUnderAverage>0&&i+right<endTime){
+ lastSignalTime++;
}
else{
- contLeft=kFALSE;
- }
-
+ break;
+ }
+ }
+
+ for(Int_t t=firstSignalTime;t<lastSignalTime;t++){
+ fDataSignals[t]=(AliHLTTPCSignal_t)(fDataSignals[t]-averageValue + valueUnderAverage);
+ fSignalPositionArray[fSizeOfSignalPositionArray]=t;
+ fSizeOfSignalPositionArray++;
+ // Matthias Oct 10 2008: trying hard to make the code faster for the
+ // AltroChannelSelection. For that we only need to know there is data
+ if (speedup) return;
}
+ i+=lastSignalTime;
}
}
- Int_t nReadFromPositionArray=0;
- for(Int_t i=endTime;i>=beginTime;i--){
- if(i==fSignalPositionArray[nReadFromPositionArray]){
- nReadFromPositionArray++;
+ //reset the rest of the data
+ Int_t counterSize=fSizeOfSignalPositionArray;
+
+ for(Int_t d=endTime;d>=beginTime;d--){
+ if(d==fSignalPositionArray[counterSize-1]&&counterSize-1>=0){
+ counterSize--;
}
else{
- fDataSignals[i]=-1;
- }
+ fDataSignals[d]=-1;
+ }
+ }
+ if(fDataSignals[beginTime+1]<1){
+ fDataSignals[beginTime]=0;
}
}
-void AliHLTTPCPad::FindClusterCandidates()
-{
- // see header file for class documentation
+void AliHLTTPCPad::AddClusterCandidate(const AliHLTTPCClusters& candidate){
+ fClusterCandidates.push_back(candidate);
+ fUsedClusterCandidates.push_back(0);
+}
- if(fNSigmaThreshold>0){
- ZeroSuppress(fNSigmaThreshold);
- }
- else if(fSignalThreshold>0){
- ZeroSuppress((Double_t)0,(Int_t)fSignalThreshold);
- }
- UInt_t seqcharge=0;
- UInt_t seqaverage=0;
- UInt_t seqerror=0;
- vector<Int_t> tmpPos;
- vector<Int_t> tmpSig;
- UInt_t isFalling=0;
-
- for(Int_t pos=fSizeOfSignalPositionArray-2;pos>=0;pos--){
- if(fSignalPositionArray[pos]==fSignalPositionArray[pos+1]+1){
- seqcharge+=fDataSignals[fSignalPositionArray[pos+1]];
- seqaverage += fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
- seqerror += fSignalPositionArray[pos+1]*fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
-
- tmpPos.push_back(fSignalPositionArray[pos+1]);
- tmpSig.push_back(fDataSignals[fSignalPositionArray[pos+1]]);
-
- if(fDataSignals[fSignalPositionArray[pos+1]]>fDataSignals[fSignalPositionArray[pos]]){
- isFalling=1;
- }
- if(fDataSignals[fSignalPositionArray[pos+1]]<fDataSignals[fSignalPositionArray[pos]]&&isFalling){
- Int_t seqmean=0;
- seqmean = seqaverage/seqcharge;
-
- //Calculate mean in pad direction:
- Int_t padmean = seqcharge*fPadNo;
- Int_t paderror = fPadNo*padmean;
- AliHLTTPCClusters candidate;
- candidate.fTotalCharge = seqcharge;
- candidate.fPad = padmean;
- candidate.fPad2 = paderror;
- candidate.fTime = seqaverage;
- candidate.fTime2 = seqerror;
- candidate.fMean = seqmean;
- candidate.fLastMergedPad = fPadNo;
- fClusterCandidates.push_back(candidate);
- fUsedClusterCandidates.push_back(0);
- isFalling=0;
- seqcharge=0;
- seqaverage=0;
- seqerror=0;
-
- tmpPos.clear();
- tmpSig.clear();
-
- continue;
- }
-
- if(pos<1){
- seqcharge+=fDataSignals[fSignalPositionArray[0]];
- seqaverage += fSignalPositionArray[0]*fDataSignals[fSignalPositionArray[0]];
- seqerror += fSignalPositionArray[0]*fSignalPositionArray[0]*fDataSignals[fSignalPositionArray[0]];
- tmpPos.push_back(fSignalPositionArray[0]);
- tmpSig.push_back(fDataSignals[fSignalPositionArray[0]]);
-
- //Calculate mean of sequence:
- Int_t seqmean=0;
- seqmean = seqaverage/seqcharge;
-
- //Calculate mean in pad direction:
- Int_t padmean = seqcharge*fPadNo;
- Int_t paderror = fPadNo*padmean;
- AliHLTTPCClusters candidate;
- candidate.fTotalCharge = seqcharge;
- candidate.fPad = padmean;
- candidate.fPad2 = paderror;
- candidate.fTime = seqaverage;
- candidate.fTime2 = seqerror;
- candidate.fMean = seqmean;
- candidate.fLastMergedPad = fPadNo;
- fClusterCandidates.push_back(candidate);
- fUsedClusterCandidates.push_back(0);
- isFalling=0;
- seqcharge=0;
- seqaverage=0;
- seqerror=0;
-
- tmpPos.clear();
- tmpSig.clear();
- }
- }
- else if(seqcharge>0){
- seqcharge+=fDataSignals[fSignalPositionArray[pos+1]];
- seqaverage += fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
- seqerror += fSignalPositionArray[pos+1]*fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
- tmpPos.push_back(fSignalPositionArray[pos+1]);
- tmpSig.push_back(fDataSignals[fSignalPositionArray[pos+1]]);
-
- //Calculate mean of sequence:
- Int_t seqmean=0;
- seqmean = seqaverage/seqcharge;
-
- //Calculate mean in pad direction:
- Int_t padmean = seqcharge*fPadNo;
- Int_t paderror = fPadNo*padmean;
- AliHLTTPCClusters candidate;
- candidate.fTotalCharge = seqcharge;
- candidate.fPad = padmean;
- candidate.fPad2 = paderror;
- candidate.fTime = seqaverage;
- candidate.fTime2 = seqerror;
- candidate.fMean = seqmean;
- candidate.fLastMergedPad = fPadNo;
- fClusterCandidates.push_back(candidate);
- fUsedClusterCandidates.push_back(0);
- isFalling=0;
- seqcharge=0;
- seqaverage=0;
- seqerror=0;
-
- tmpPos.clear();
- tmpSig.clear();
- }
- }
+void AliHLTTPCPad::AddCandidateDigits(const vector<AliHLTTPCDigitData>& candidateDigits){
+ fCandidateDigitsVector.push_back(candidateDigits);
}
+vector<AliHLTTPCDigitData> *AliHLTTPCPad::GetCandidateDigits(Int_t candidateIndex){
+ return (size_t(candidateIndex) < fCandidateDigitsVector.size()) ? &fCandidateDigitsVector.at(candidateIndex) :0;
+}