zero suppression and selction of active pads (Kenneth)
authorrichterm <richterm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 24 Nov 2007 09:35:11 +0000 (09:35 +0000)
committerrichterm <richterm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 24 Nov 2007 09:35:11 +0000 (09:35 +0000)
HLT/TPCLib/AliHLTTPCClusterFinder.cxx
HLT/TPCLib/AliHLTTPCClusterFinder.h
HLT/TPCLib/AliHLTTPCClusterFinderComponent.cxx
HLT/TPCLib/AliHLTTPCClusterFinderComponent.h
HLT/TPCLib/AliHLTTPCDefinitions.cxx
HLT/TPCLib/AliHLTTPCDefinitions.h
HLT/TPCLib/AliHLTTPCDigitReaderPacked.h
HLT/TPCLib/AliHLTTPCPad.cxx
HLT/TPCLib/AliHLTTPCPad.h
HLT/TPCLib/AliHLTTPCPadArray.cxx
HLT/TPCLib/AliHLTTPCPadArray.h

index 6f030d9..92d811f 100644 (file)
@@ -141,7 +141,8 @@ AliHLTTPCClusterFinder::AliHLTTPCClusterFinder()
   fXYErr(0.2),
   fZErr(0.3),
   fOccupancyLimit(1.0),
-  fPadArray(NULL)
+  fPadArray(NULL),
+  fNSigmaThreshold(0)
 {
   //constructor
 }
@@ -684,6 +685,7 @@ void AliHLTTPCClusterFinder::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
 //----------------------------------Methods for the new unsorted way of reading the data --------------------------------
 
 void AliHLTTPCClusterFinder::SetPadArray(AliHLTTPCPadArray * padArray){
+  // see header file for function documentation
   fPadArray=padArray;
 }
 
@@ -695,9 +697,17 @@ void AliHLTTPCClusterFinder::ReadDataUnsorted(void* ptr,unsigned long size){
   fDigitReader->InitBlock(fPtr,fSize,fFirstRow,fLastRow,fCurrentPatch,fCurrentSlice);
 
   fPadArray->SetDigitReader(fDigitReader);
+
+  if(fSignalThreshold>0){
+    fPadArray->SetSignalThreshold(fSignalThreshold);
+  }
+  if(fNSigmaThreshold>0){
+    fPadArray->SetNSigmaThreshold(fNSigmaThreshold);
+  }
   fPadArray->ReadData();
 }
 void AliHLTTPCClusterFinder::FindClusters(){
+  // see header file for function documentation
   fPadArray->FindClusterCandidates();
   fPadArray->FindClusters(fMatch);
 
@@ -719,6 +729,10 @@ void AliHLTTPCClusterFinder::FindClusters(){
   delete [] clusterlist;
   fPadArray->DataToDefault();
 }
+Int_t AliHLTTPCClusterFinder::GetActivePads(AliHLTTPCPadArray::AliHLTTPCActivePads* activePads,Int_t maxActivePads){
+  // see header file for function documentation
+  return fPadArray->GetActivePads((AliHLTTPCPadArray::AliHLTTPCActivePads*)activePads , maxActivePads);
+}
 void AliHLTTPCClusterFinder::WriteClusters(Int_t nclusters,AliHLTTPCClusters *list)//This is used when using the AliHLTTPCClusters class for cluster data
 {
   //write cluster to output pointer
index 202c622..ad1b3e3 100644 (file)
@@ -20,7 +20,7 @@
 // or
 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
 
-#define UNSORTED 0
+#define UNSORTED 1
 
 #include "AliHLTLogging.h"
 #include "AliHLTTPCPadArray.h"
@@ -66,6 +66,7 @@ class AliHLTTPCClusterFinder : public AliHLTLogging {
   void SetThreshold(UInt_t i) {fThreshold=i;}
   void SetOccupancyLimit(Float_t f) {fOccupancyLimit=f;}
   void SetSignalThreshold(Int_t i) {fSignalThreshold=i;}
+  void SetNSigmaThreshold(Double_t d) {fNSigmaThreshold=d;}
   void SetMatchWidth(UInt_t i) {fMatch=i;}
   void SetSTDOutput(Bool_t f=kFALSE) {fStdout=f;}  
   void SetCalcErr(Bool_t f=kTRUE) {fCalcerr=f;}
@@ -77,6 +78,7 @@ class AliHLTTPCClusterFinder : public AliHLTLogging {
   void SetPadArray(AliHLTTPCPadArray *padArray);
   void ReadDataUnsorted(void* ptr,unsigned long size);
   void FindClusters();
+  Int_t GetActivePads(AliHLTTPCPadArray::AliHLTTPCActivePads* activePads,Int_t maxActivePads);
   void WriteClusters(Int_t nclusters,AliHLTTPCClusters *list);
   
  private: 
@@ -105,7 +107,9 @@ class AliHLTTPCClusterFinder : public AliHLTLogging {
   Int_t fMatch;          //size of match
   UInt_t fThreshold;     //threshold for clusters
   /** threshold for zero suppression (applied per bin) */
-  Int_t fSignalThreshold;
+  Int_t fSignalThreshold;// see above
+  /** threshold for zero suppression 2007 December run */
+  Double_t fNSigmaThreshold; // see above
   Int_t fNClusters;      //number of found clusters
   Int_t fMaxNClusters;   //max. number of clusters
   Float_t fXYErr;        //fixed error in XY
index e6d9195..2dd63af 100644 (file)
@@ -66,7 +66,8 @@ AliHLTTPCClusterFinderComponent::AliHLTTPCClusterFinderComponent(bool packed)
   fPackedSwitch(packed),
   fUnsorted(0),
   fPatch(0),
-  fPadArray(NULL)
+  fPadArray(NULL),
+  fGetActivePads(0)
 {
   // see header file for class documentation
   // or
@@ -130,6 +131,7 @@ int AliHLTTPCClusterFinderComponent::DoInit( int argc, const char** argv )
 
   Int_t rawreadermode =  -1;
   Int_t sigthresh = -1;
+  Double_t sigmathresh= -1;
   Float_t occulimit = 1.0;
   Int_t oldRCUFormat=0;
   // Data Format version numbers:
@@ -241,6 +243,28 @@ int AliHLTTPCClusterFinderComponent::DoInit( int argc, const char** argv )
       continue;
     }
 
+    // -- checking for active pads, used in 2007 December run
+    if ( !strcmp( argv[i], "activepads" ) ) {
+      fGetActivePads = strtoul( argv[i+1], &cpErr ,0);
+      if ( *cpErr ){
+       HLTError("Cannot convert activepads specifier '%s'. Should  be 0(off) or 1(on), must be integer", argv[i+1]);
+       return EINVAL;
+      }
+      i+=2;
+      continue;
+    }
+
+    // -- checking for nsigma-threshold, used in 2007 December run in ZeroSuppression
+    if ( !strcmp( argv[i], "nsigma-threshold" ) ) {
+       sigmathresh = strtoul( argv[i+1], &cpErr ,0);
+      if ( *cpErr ){
+       HLTError("Cannot convert nsigma-threshold specifier '%s'. Must be integer", argv[i+1]);
+       return EINVAL;
+      }
+      i+=2;
+      continue;
+    }
+
     Logging(kHLTLogError, "HLT::TPCClusterFinder::DoInit", "Unknown Option", "Unknown option '%s'", argv[i] );
     return EINVAL;
 
@@ -304,7 +328,7 @@ int AliHLTTPCClusterFinderComponent::DoInit( int argc, const char** argv )
   if ( (fXYClusterError>0) && (fZClusterError>0) )
     fClusterFinder->SetCalcErr( false );
   fClusterFinder->SetSignalThreshold(sigthresh);
-    
+  fClusterFinder->SetNSigmaThreshold(sigmathresh);
   if(fUnsorted&&fPatch>-1&&fPatch<6){
     fPadArray = new AliHLTTPCPadArray(fPatch);
     fPadArray->InitializeVector();
@@ -454,6 +478,28 @@ int AliHLTTPCClusterFinderComponent::DoEvent( const AliHLTComponentEventData& ev
       outBPtr += mysize;
       outPtr = (AliHLTTPCClusterData*)outBPtr;
        
+      if(fGetActivePads){
+       AliHLTTPCPadArray::AliHLTTPCActivePads* outPtrActive;
+       UInt_t activePadsSize, activePadsN = 0;
+       outPtrActive = (AliHLTTPCPadArray::AliHLTTPCActivePads*)outBPtr;
+       offset=tSize;
+       Int_t maxActivePads = (size-tSize)/sizeof(AliHLTTPCPadArray::AliHLTTPCActivePads);
+       activePadsSize= fClusterFinder->GetActivePads((AliHLTTPCPadArray::AliHLTTPCActivePads*)outPtrActive,maxActivePads)*sizeof(AliHLTTPCPadArray::AliHLTTPCActivePads);
+       
+       AliHLTComponentBlockData bdActive;
+       FillBlockData( bdActive );
+       bdActive.fOffset = offset;
+       bdActive.fSize = activePadsSize;
+       bdActive.fSpecification = iter->fSpecification;
+       bdActive.fDataType = AliHLTTPCDefinitions::fgkActivePadsDataType;
+       outputBlocks.push_back( bdActive );
+       
+       tSize+=activePadsSize;
+       outBPtr += activePadsSize;
+       outPtrActive = (AliHLTTPCPadArray::AliHLTTPCActivePads*)outBPtr;
+      }
+
       if ( tSize > size )
        {
          Logging( kHLTLogFatal, "HLT::TPCClusterFinder::DoEvent", "Too much data", 
index 192a2a1..3b024b0 100644 (file)
@@ -119,12 +119,18 @@ class AliHLTTPCClusterFinderComponent : public AliHLTProcessor
       Int_t fPatch;                                                                  //!transient
 
       /*
+       * Switch to specify if one ship out a list of active pads.
+       * Used for the 2007 December run. 
+       */
+      Int_t fGetActivePads;                                                          //!transient
+
+      /*
        * Pointer to a PadArray object containing a double array of all the pads in
        * the current patch.
        */
       AliHLTTPCPadArray * fPadArray;                                                 //!transient
 
-      ClassDef(AliHLTTPCClusterFinderComponent, 1)
+      ClassDef(AliHLTTPCClusterFinderComponent, 2)
 
     };
 #endif
index ce6f590..15ec2b9 100644 (file)
@@ -51,6 +51,8 @@ const AliHLTComponentDataType AliHLTTPCDefinitions::fgkCalibPedestalDataType =
   (AliHLTComponentDataType){sizeof(AliHLTComponentDataType),                          {'C','A','L','_','P','E','D',' '},  kAliHLTDataOriginAny} | kAliHLTDataOriginTPC;
 const AliHLTComponentDataType AliHLTTPCDefinitions::fgkCalibPulserDataType =                                                                                 
   (AliHLTComponentDataType){sizeof(AliHLTComponentDataType),                          {'C','A','L','_','P','U','L','S'},  kAliHLTDataOriginAny} | kAliHLTDataOriginTPC;
+const AliHLTComponentDataType AliHLTTPCDefinitions::fgkActivePadsDataType =                                                                                  
+  (AliHLTComponentDataType){sizeof(AliHLTComponentDataType),                          {'A','C','T','I','V','E','P','A'},  kAliHLTDataOriginAny} | kAliHLTDataOriginTPC;
 
 
 AliHLTTPCDefinitions::AliHLTTPCDefinitions()
index 02c2110..7366c51 100644 (file)
@@ -80,8 +80,10 @@ public:
   static const AliHLTComponentDataType fgkCalibPedestalDataType;   // see above
   /** signal calibration data */
   static const AliHLTComponentDataType fgkCalibPulserDataType;     // see above
+  /** active pads data type, Used for cosmics test december 2007 */
+  static const AliHLTComponentDataType fgkActivePadsDataType;     // see above
 
-  ClassDef(AliHLTTPCDefinitions, 0);
+  ClassDef(AliHLTTPCDefinitions, 1);
 
 };
 
index 509c641..8303df9 100644 (file)
@@ -13,7 +13,7 @@
     @brief  A digit reader implementation for simulated, packed TPC 'raw' data.
 */
 
-#define ENABLE_PAD_SORTING 1
+#define ENABLE_PAD_SORTING 0
 
 #include "AliHLTTPCDigitReader.h"
 
index 24e8c3f..ecdd538 100644 (file)
@@ -72,7 +72,9 @@ AliHLTTPCPad::AliHLTTPCPad()
   fpRawData(NULL),
   fDataSignals(NULL),
   fSignalPositionArray(NULL),
-  fSizeOfSignalPositionArray(0)
+  fSizeOfSignalPositionArray(0),
+  fNSigmaThreshold(0),
+  fSignalThreshold(0)
 {
   // see header file for class documentation
   // or
@@ -110,7 +112,9 @@ AliHLTTPCPad::AliHLTTPCPad(Int_t offset, Int_t nofBins)
   fpRawData(NULL),
   fDataSignals(NULL),
   fSignalPositionArray(NULL),
-  fSizeOfSignalPositionArray(0)
+  fSizeOfSignalPositionArray(0),
+  fNSigmaThreshold(0),
+  fSignalThreshold(0)
 {
   // see header file for class documentation
 }
@@ -125,7 +129,7 @@ AliHLTTPCPad::~AliHLTTPCPad()
   if (fDataSignals) {
     AliHLTTPCSignal_t* pData=fDataSignals;
     fDataSignals=NULL;
-   delete [] pData;
+    delete [] pData;
   }
   if (fSignalPositionArray) {
     //AliHLTTPCSignal_t* pData=fSignalPositionArray;
@@ -193,10 +197,10 @@ Int_t AliHLTTPCPad::CalculateBaseLine(Int_t reqMinCount)
            fAverage=fSum/fCount;
            //HLTDebug("new average %d", fAverage);
          } else {
-//         HLTDebug("baseline re-eveluation skipped because of to few "
-//                    "contributing bins: total=%d, contributing=%d, req=%d"
-//                    "\ndata might be already zero suppressed"
-//                    , fTotal, fCount, reqMinCount);
+           //      HLTDebug("baseline re-eveluation skipped because of to few "
+           //                 "contributing bins: total=%d, contributing=%d, req=%d"
+           //                 "\ndata might be already zero suppressed"
+           //                 , fTotal, fCount, reqMinCount);
            iResult=-ENODATA;
          }
          fCount=0;fSum=-1;
@@ -216,9 +220,9 @@ Int_t AliHLTTPCPad::CalculateBaseLine(Int_t reqMinCount)
       fAverage=avBackup;
     }
   } else {
-//     HLTDebug("baseline calculation skipped because of to few contributing "
-//            "bins: total=%d, contributing=%d, required=%d \ndata might be "
-//            "already zero suppressed", fTotal, fCount, reqMinCount);
+    //     HLTDebug("baseline calculation skipped because of to few contributing "
+    //                "bins: total=%d, contributing=%d, required=%d \ndata might be "
+    //                "already zero suppressed", fTotal, fCount, reqMinCount);
   }
 
   return iResult;
@@ -281,9 +285,9 @@ Int_t AliHLTTPCPad::AddBaseLineValue(Int_t bin, AliHLTTPCSignal_t value)
        fBLMinBin=bin;
       }
     } else {
-//       HLTDebug("ignoring value %d (bin %d) for base line calculation "
-//            "(current average is %d)",
-//            value, bin, fAverage);
+      //       HLTDebug("ignoring value %d (bin %d) for base line calculation "
+      //              "(current average is %d)",
+      //              value, bin, fAverage);
     }
   }
   return iResult;
@@ -469,9 +473,167 @@ Int_t AliHLTTPCPad::GetDataSignal(Int_t bin) const
   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){
+  //see headerfile for documentation
+  Bool_t useSigma= kFALSE;
+  if(nSigma>0){
+    useSigma=kTRUE;
+  }
+  if(threshold<1 && nSigma<=0){
+    //setting the data to -1 for this pad
+    memset( fDataSignals, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
+    fSizeOfSignalPositionArray=0;
+    return;
+  }
+  if(endTime>=AliHLTTPCTransform::GetNTimeBins()){
+    endTime=AliHLTTPCTransform::GetNTimeBins()-1;
+  }
+  Int_t fThresholdUsed=threshold;
+  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(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--){
+      if(fDataSignals[i]>0){
+       if(fDataSignals[i]-averageValue<50){
+         sumOfDifferenceSquared+=(fDataSignals[i]-averageValue)*(fDataSignals[i]-averageValue);
+       }
+       else{
+         nAdded--;
+       }
+      }
+    }
+    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;
+  }
+  for(Int_t i=endTime+1;i<AliHLTTPCTransform::GetNTimeBins();i++){
+    fDataSignals[i]=-1;
+  }
+  // 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;
+         }
+       }
+       else if(endRight>endTime+1){
+         contRight=kFALSE;
+       }
+       endRight++;
+      }
+      for(int j=i+nToAddRight;j>i;j--){
+       fDataSignals[j]=(Int_t)(fDataSignals[j]-averageValue);
+       fSignalPositionArray[fSizeOfSignalPositionArray]=j;
+       fSizeOfSignalPositionArray++;
+      }
+      //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;
+         }
+       }
+       else{
+         contLeft=kFALSE;
+       }
+      }
+    }
+  }
+  Int_t nReadFromPositionArray=0;
+  for(Int_t i=endTime;i>=beginTime;i--){
+    if(i==fSignalPositionArray[nReadFromPositionArray]){
+      nReadFromPositionArray++;
+    }
+    else{
+      fDataSignals[i]=-1;
+    } 
+  }
+}
+
 void AliHLTTPCPad::FindClusterCandidates()
 {
   // see header file for class documentation
+
+  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;
@@ -480,7 +642,6 @@ void AliHLTTPCPad::FindClusterCandidates()
   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]];
@@ -556,34 +717,35 @@ void AliHLTTPCPad::FindClusterCandidates()
     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]]);
+      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 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();
+      //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();
     }
   }
 }
+
index 3e4baff..3ba1fa6 100644 (file)
@@ -264,6 +264,25 @@ public:
   Int_t GetDataSignal(Int_t bin) const;
 
   /**
+   * Zerosuppression where one can choose wether one want to cut on sigma(default 3) or a given adc threshold above baseline
+   * It works like this: Finds the signals above baseline+threshold, it then looks to the right of the signals adding the 
+   * the signal below threshold but over the average. It stops when the adc value rises (since you should expect a fall) 
+   * or if the signal in that bin is below average. It now does the same thing to the left. 
+   * This is a very timeconsuming approach(but good), and will likely only be used only for cosmics and laser.
+   * If you use sigma approach you can give as input n sigma or stick with default=3. For channels with very large signals 
+   * the (value-average)² is not calculated when value-average>50. This is due to sigma shooting sky high for only a few values.
+   * The method is checked with 2006 cosmics data, and it looks good.
+   * If you want to use the threshold approach you HAVE to set nSigma=-1 and threshold>0. For example: If you want all signals 
+   * 30 adc counts above threshold you should call the function like this: ZeroSuppress(-1,30)
+   * @param nSigma          Specify nSigma above threshold default=3
+   * @param threshold       Specify what adc threshold above average default=20 (remember to give nSigma=-1 if you want to use this approach)
+   * @param reqMinPoint     Required minimum number of points to do zerosuppression default AliHLTTPCTransform::GetNTimeBins/2 (1024/2).
+   * @param beginTime       Lowest timebin value. Gating grid causes some problems in the first timebins. default= 50
+   * @param endTime         Highest timebin value. Default AliHLTTPCTransform::GetNTimeBins-1
+   */
+  void ZeroSuppress(Double_t nSigma,Int_t threshold,Int_t reqMinPoint,Int_t beginTime,Int_t endTime);
+
+  /**
    * Finds the cluster candidate. If atleast two signals in the data array are neighbours
    * they are stored in a cluster candidate vector.
    */
@@ -272,7 +291,18 @@ public:
    * Prints the raw data og this pad.
    */
   void PrintRawData();
+
+   /**
+   * Set the Signal Threshold
+   */
+  void SetSignalThreshold(Int_t i){fSignalThreshold=i;}
+
+  /**
+   * Set the nSigma threshold
+   */
+  void SetNSigmaThreshold(Double_t i){fNSigmaThreshold=i;}
+
+
   /**
    * Vector of cluster candidates
    */
@@ -343,7 +373,10 @@ public:
    */
   Int_t *fSignalPositionArray;                                     //! transient
   Int_t fSizeOfSignalPositionArray;                                //! transient
+  
+  Double_t fNSigmaThreshold;                                       //! transient
+  Double_t fSignalThreshold;                                       //! transient
 
-  ClassDef(AliHLTTPCPad, 1)
+  ClassDef(AliHLTTPCPad, 3)
 };
 #endif // ALIHLTTPCPAD_H
index 8152b3e..8c2fdec 100644 (file)
@@ -55,7 +55,9 @@ AliHLTTPCPadArray::AliHLTTPCPadArray()
   fThreshold(10),
   fNumberOfPadsInRow(NULL),
   fNumberOfRows(0),
-  fDigitReader(NULL)
+  fDigitReader(NULL),
+  fSignalThreshold(0),
+  fNSigmaThreshold(0)
 {
   // see header file for class documentation
   // or
@@ -74,7 +76,9 @@ AliHLTTPCPadArray::AliHLTTPCPadArray(Int_t patch)
   fThreshold(10),
   fNumberOfPadsInRow(NULL),
   fNumberOfRows(0),
-  fDigitReader(NULL)
+  fDigitReader(NULL),
+  fSignalThreshold(0),
+  fNSigmaThreshold(0)
 {
   // see header file for class documentation
 }
@@ -451,6 +455,7 @@ void AliHLTTPCPadArray::PrintClusters()
 
 void AliHLTTPCPadArray::DataToDefault()
 {
+  //see header file for documentation
   for(Int_t i=0;i<fNumberOfRows;i++){
     for(Int_t j=0;j<fNumberOfPadsInRow[i];j++){
        fRowPadVector[i][j]->SetDataToDefault();
@@ -460,9 +465,50 @@ void AliHLTTPCPadArray::DataToDefault()
 
 void AliHLTTPCPadArray::FindClusterCandidates()
 {
-  for(Int_t row=0;row<fNumberOfRows;row++){
-    for(Int_t pad=0;pad<fNumberOfPadsInRow[row];pad++){
+  //see header file for documentation
+  if(fNSigmaThreshold>0){
+    for(Int_t row=0;row<fNumberOfRows;row++){
+      for(Int_t pad=0;pad<fNumberOfPadsInRow[row];pad++){
+       fRowPadVector[row][pad]->SetNSigmaThreshold(fNSigmaThreshold);
+       fRowPadVector[row][pad]->FindClusterCandidates();
+      }
+    }
+  }
+  else if(fSignalThreshold>0){
+    for(Int_t row=0;row<fNumberOfRows;row++){
+      for(Int_t pad=0;pad<fNumberOfPadsInRow[row];pad++){
+       fRowPadVector[row][pad]->SetSignalThreshold(fSignalThreshold);
        fRowPadVector[row][pad]->FindClusterCandidates();
+      }
     }
   }
+  else{
+    for(Int_t row=0;row<fNumberOfRows;row++){
+      for(Int_t pad=0;pad<fNumberOfPadsInRow[row];pad++){
+       fRowPadVector[row][pad]->FindClusterCandidates();
+      }
+    }
+  }
+}
+Int_t AliHLTTPCPadArray::GetActivePads(AliHLTTPCActivePads * activePads,Int_t maxActivePads){
+  //see header file for documentation
+
+  Int_t counter=0;
+  for(Int_t row=0;row<fNumberOfRows;row++){
+    for(Int_t pad=0;pad<fNumberOfPadsInRow[row]-1;pad++){
+      if(fRowPadVector[row][pad]->fClusterCandidates.size()>0){
+       AliHLTTPCActivePads tmpAP;
+       tmpAP.fRow=row;
+       tmpAP.fPad=pad;
+       activePads[counter]= tmpAP;
+       counter++;
+       if(counter>=maxActivePads){
+         return counter;
+       }
+      }
+       
+    }
+  }  
+  return counter;
 }
index 82aaed9..d4bacca 100644 (file)
@@ -37,6 +37,14 @@ class AliHLTTPCPadArray : public AliHLTLogging {
 
 public:
 
+  struct AliHLTTPCActivePads
+  {
+    UInt_t fPad;           //pad value
+    UInt_t fRow;           //row value
+  };
+  typedef struct AliHLTTPCActivePads AliHLTTPCActivePads; //!
+
+
   /** standard constructor */
   AliHLTTPCPadArray();
 
@@ -104,6 +112,24 @@ public:
    */
   void PrintClusters();
 
+  /**
+   * Set the Signal Threshold
+   */
+  void SetSignalThreshold(Int_t i){fSignalThreshold=i;}
+
+  /**
+   * Set the nSigma threshold
+   */
+  void SetNSigmaThreshold(Double_t i){fNSigmaThreshold=i;}
+
+
+  /**
+   * Loop over all pads adding pads with signal to active pads
+   * Returns number of active pads
+   *
+   */
+  Int_t GetActivePads(AliHLTTPCActivePads* activePads, Int_t maxActivePads);
+
   typedef vector<AliHLTTPCPad*> AliHLTTPCPadVector;
 
   vector<AliHLTTPCPadVector> fRowPadVector;                        //! transient
@@ -123,9 +149,16 @@ private:
 
   Int_t fLastRow;                                                  //! transient
 
-  //TODO: I suggest making the following UInt_t if it is never supposed to be negative.
+  //TODO: I suggest making the following UInt_t if it is never supposed to be negative. Will do!
+  /* total charge of Cluster threshold*/
   Int_t fThreshold;                                                //! transient
 
+  //TODO: I suggest making the following UInt_t if it is never supposed to be negative.
+  Int_t fSignalThreshold;                                                //! transient
+
+  //TODO: I suggest making the following UInt_t if it is never supposed to be negative.
+  Double_t fNSigmaThreshold;                                                //! transient
+
   Int_t* fNumberOfPadsInRow;                                       //! transient
 
   Int_t fNumberOfRows;                                             //! transient