ATO-78 - Robust filters - modifications for boundary values

AliTPCCalROC.cxx AliTPCCalROC.h test/UnitTest.C

AliTPCCalROC.cxx  - handling of edge effects for the median and LTM filters (http://en.wikipedia.org/wiki/Median_filter#Edge_preservation_properties) + mirroring values at the boundaries
AliTPCCalROC.h    - check the ranges
test/UnitTest.C   - check invariants of median and LTM filters - assert in case of failures

diff --git a/TPC/Base/AliTPCCalROC.cxx b/TPC/Base/AliTPCCalROC.cxx
index 933b4aed0b6975ef499ed47692c610b9692c1c34..20ab069542e255d1cc56b11168e2cbbaebf4d043 100644 (file)
--- a/TPC/Base/AliTPCCalROC.cxx
+++ b/TPC/Base/AliTPCCalROC.cxx
@@ -143,29 +143,47 @@ void AliTPCCalROC::Streamer(TBuffer &R__b)
    }
 }
 
+
+//
+
 Bool_t AliTPCCalROC::MedianFilter(Int_t deltaRow, Int_t deltaPad){
   //
-  //   Modify content of the class
-  //   write median
+  //   Modify content of the object - raplace value by median in neighorhood
+  //
   Float_t *newBuffer=new Float_t[fNChannels] ;
-  Float_t *cacheBuffer=new Float_t[fNChannels] ;
+  Double_t *cacheBuffer=new Double_t[fNChannels];
   //
-  for (UInt_t iRow=0; iRow<fNRows; iRow++){
+  for (Int_t iRow=0; iRow< Int_t(fNRows); iRow++){
     Int_t nPads=GetNPads(iRow); // number of rows in current row
     for (Int_t iPad=0; iPad<nPads; iPad++){
+      Double_t value=GetValue(iRow,iPad);
       Int_t counter=0;
+      //
       for (Int_t dRow=-deltaRow; dRow<=deltaRow; dRow++){
-       Int_t jRow=iRow+dRow;
-       Int_t jPads= GetNPads(jRow);
-       if (jPads==0) continue;
+       Int_t jRow=iRow+dRow;  //take the row - mirror if ouside of range
+       Float_t sign0=1.;
+       if (jRow<0) sign0=-1.;
+       if (UInt_t(jRow)>=fNRows) sign0=-1.;    
+       Int_t jPads= GetNPads(iRow+sign0*dRow);
        Int_t offset=(nPads-jPads)/2.;
+       //
        for (Int_t dPad=-deltaPad; dPad<=deltaPad; dPad++){
-         Int_t jPad=iPad+dPad+offset;
-         if (jPad<0 || jPad>=jPads) continue;
-         cacheBuffer[counter++]=GetValue(jRow,jPad);
+         Float_t sign=sign0;
+         Int_t jPad=iPad+dPad+offset;  //take the pad - mirror if ouside of range
+         Int_t kRow=jRow;
+         if (jPad<0) sign=-1;
+         if (jPad>=jPads) sign=-1;
+         if (sign<0){
+           kRow=iRow-dRow;
+           jPad=iPad-dPad+offset;
+         }       
+         if (IsInRange(UInt_t(kRow),jPad)){
+           cacheBuffer[counter]=sign*(GetValue(kRow,jPad)-value);
+           counter++;
+         }
        }
       }
-      newBuffer[fkIndexes[iRow]+iPad] = TMath::Median(counter,cacheBuffer);
+      newBuffer[fkIndexes[iRow]+iPad] = TMath::Median(counter, cacheBuffer)+value;
     }
   }
   memcpy(fData, newBuffer,GetNchannels()*sizeof(Float_t));
@@ -174,33 +192,51 @@ Bool_t AliTPCCalROC::MedianFilter(Int_t deltaRow, Int_t deltaPad){
   return kTRUE;
 }
 
+
+
 Bool_t AliTPCCalROC::LTMFilter(Int_t deltaRow, Int_t deltaPad, Float_t fraction, Int_t type){
   //
   //
   // //
   //   Modify content of the class
-  //   write median
+  //   write LTM mean or median
   if (fraction<0 || fraction>1) return kFALSE;
   Float_t *newBuffer=new Float_t[fNChannels] ;
   Double_t *cacheBuffer=new Double_t[fNChannels];
   //
-  for (UInt_t iRow=0; iRow<fNRows; iRow++){
+  for (Int_t iRow=0; iRow< Int_t(fNRows); iRow++){
     Int_t nPads=GetNPads(iRow); // number of rows in current row
     for (Int_t iPad=0; iPad<nPads; iPad++){
+      Double_t value=GetValue(iRow,iPad);
       Int_t counter=0;
+      //
       for (Int_t dRow=-deltaRow; dRow<=deltaRow; dRow++){
-       Int_t jRow=iRow+dRow;
-       Int_t jPads= GetNPads(jRow);
-       if (jPads==0) continue;
+       Int_t jRow=iRow+dRow;  //take the row - mirror if ouside of range
+       Float_t sign0=1.;
+       if (jRow<0) sign0=-1.;
+       if (UInt_t(jRow)>=fNRows) sign0=-1.;    
+       Int_t jPads= GetNPads(iRow+sign0*dRow);
        Int_t offset=(nPads-jPads)/2.;
+       //
        for (Int_t dPad=-deltaPad; dPad<=deltaPad; dPad++){
-         Int_t jPad=iPad+dPad+offset;
-         if (jPad<0 || jPad>=jPads) continue;
-         cacheBuffer[counter++]=GetValue(jRow,jPad);
+         Float_t sign=sign0;
+         Int_t jPad=iPad+dPad+offset;  //take the pad - mirror if ouside of range
+         Int_t kRow=jRow;
+         if (jPad<0) sign=-1;
+         if (jPad>=jPads) sign=-1;
+         if (sign<0){
+           kRow=iRow-dRow;
+           jPad=iPad-dPad+offset;
+         } 
+         if (IsInRange(UInt_t(kRow),jPad)){
+           cacheBuffer[counter]=sign*(GetValue(kRow,jPad)-value);
+           counter++;
+         }
        }
       }
       Double_t mean,rms;
-      AliMathBase::EvaluateUni(counter,cacheBuffer,mean,rms,fraction*Double_t(counter));
+      AliMathBase::EvaluateUni(counter,cacheBuffer,mean,rms,TMath::Nint(fraction*Double_t(counter)));
+      mean+=value;
       newBuffer[fkIndexes[iRow]+iPad] = (type==0)? mean:rms;
     }
   }

diff --git a/TPC/Base/AliTPCCalROC.h b/TPC/Base/AliTPCCalROC.h
index 374ba8523a2c4dc949d40e67a85768707310e985..1baee19b846719a0c0d5f9a5ef4a2155fa51e12d 100644 (file)
--- a/TPC/Base/AliTPCCalROC.h
+++ b/TPC/Base/AliTPCCalROC.h
@@ -31,6 +31,7 @@ class AliTPCCalROC : public TNamed {
   UInt_t        GetSector() const { return fSector;}
   UInt_t        GetNrows() const               { return fNRows;};
   UInt_t        GetNchannels()       const     { return fNChannels;};
+  Bool_t        IsInRange(UInt_t row, UInt_t pad) { return (row<fNRows)? (fkIndexes[row]+pad)<fNChannels:kFALSE;}
   UInt_t        GetNPads(UInt_t row)  const     { return (row<fNRows)? AliTPCROC::Instance()->GetNPads(fSector,row):0;};
   Float_t      GetValue(UInt_t row, UInt_t pad) const { return ( (row<fNRows) && (fkIndexes[row]+pad)<fNChannels)? fData[fkIndexes[row]+pad]: 0; };
   Float_t      GetValue(UInt_t channel) const { return  fData[channel]; };

diff --git a/TPC/Base/test/UnitTest.C b/TPC/Base/test/UnitTest.C
index acdeae6931d5da61100beb9c396018fc0fc42cd4..780f37247bf54d3c88d8ae67708fdcc04efc2717 100644 (file)
--- a/TPC/Base/test/UnitTest.C
+++ b/TPC/Base/test/UnitTest.C
@@ -1,10 +1,11 @@
 /*
-  Unit test for fucntions used in the calibrations
+  Unit test for fucntions used in the Base directory:
   gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/ -I$ALICE_ROOT/include -I$ALICE_ROOT/STEER\
    -I$ALICE_ROOT/TPC -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TRD -I$ALICE_ROOT/TOF -I$ALICE_ROOT/RAW -I$ALICE_ROOT/PWG1 -I$ALICE_ROOT/STAT -I$ALICE_ROOT/TPC/Base -I$ALICE_ROOT/TPC/Calib");
    
   .L $ALICE_ROOT/TPC/Base/test/UnitTest.C+ 
- 
+  UnitTestAliTPCCalPadTree();
+
 */
 
 #include "TF1.h"
@@ -21,7 +22,7 @@
 #include "AliTPCcalibDButil.h"
 
 //
-// PARAMETERS:
+// PARAMETERS to set from outside:
 //
 TString baseDir="/hera/alice/wiechula/calib/guiTrees";
 //
@@ -30,33 +31,64 @@ TString baseDir="/hera/alice/wiechula/calib/guiTrees";
 
 void  UnitTestAliTPCCalPadTree(){
   //
-  //
-  //
+  //  Make a UnitTest of the AliTPCCalPad 
+  //   a.) TTree functionaility
+  //   b.) MedianFilterFunctionality
+  //   c.) LTMFilterFunctionality
+  // 
   TObjArray *fArray = new TObjArray(100);
   TTree * treePad=AliTPCcalibDButil::ConnectGainTrees(baseDir);
-  for (Int_t i=0; i<3; i++){
-    AliTPCCalPad * padLx = AliTPCCalPad::MakePadFromTree(treePad,"lx.fElements","QMax");
-    AliTPCCalPad * padLy = AliTPCCalPad::MakePadFromTree(treePad,"ly.fElements","QMax");
+  for (Int_t i=0; i<5; i+=2){
+    AliTPCCalPad * padLx = AliTPCCalPad::MakePadFromTree(treePad,"lx.fElements","Lx");
+    AliTPCCalPad * padLy = AliTPCCalPad::MakePadFromTree(treePad,"ly.fElements","Ly");
+    AliTPCCalPad * padLLx = AliTPCCalPad::MakePadFromTree(treePad,"lx.fElements","LLx");
+    AliTPCCalPad * padLLy = AliTPCCalPad::MakePadFromTree(treePad,"ly.fElements","LLy");
     AliTPCCalPad * padMax = AliTPCCalPad::MakePadFromTree(treePad,"QA.2010.LHC10d.MaxCharge.fElements","QMax");
     AliTPCCalPad * padMean = AliTPCCalPad::MakePadFromTree(treePad,"QA.2010.LHC10d.MeanCharge.fElements","QTot");
+    AliTPCCalPad * padMaxL = AliTPCCalPad::MakePadFromTree(treePad,"QA.2010.LHC10d.MaxCharge.fElements","QMax");
+    AliTPCCalPad * padMeanL = AliTPCCalPad::MakePadFromTree(treePad,"QA.2010.LHC10d.MeanCharge.fElements","QTot");
     if (i>0) {
       padLx->MedianFilter(i,2*i);
       padLy->MedianFilter(i,2*i);
+      padLLx->LTMFilter(i,2*i,1.00, 0);
+      padLLy->LTMFilter(i,2*i,1.00, 0);
       padMax->MedianFilter(i,2*i);
       padMean->MedianFilter(i,2*i);
+      padMaxL->LTMFilter(i,2*i,0.8,0);
+      padMeanL->LTMFilter(i,2*i,0.8,0);
     }
-    padLx->SetName(TString::Format("QLx%d",i).Data());
-    padLy->SetName(TString::Format("QLy%d",i).Data());
+    padLx->SetName(TString::Format("Lx%d",i).Data());
+    padLy->SetName(TString::Format("Ly%d",i).Data());
+    padLLx->SetName(TString::Format("LLx%d",i).Data());
+    padLLy->SetName(TString::Format("LLy%d",i).Data());
     padMax->SetName(TString::Format("QMax%d",i).Data());
     padMean->SetName(TString::Format("QTot%d",i).Data());
+    padMaxL->SetName(TString::Format("QMaxL%d",i).Data());
+    padMeanL->SetName(TString::Format("QTotL%d",i).Data());
     fArray->AddLast(padLx);
     fArray->AddLast(padLy);
+    fArray->AddLast(padLLx);
+    fArray->AddLast(padLLy);
     fArray->AddLast(padMax);
     fArray->AddLast(padMean);
+    fArray->AddLast(padMaxL);
+    fArray->AddLast(padMeanL);
   }
   AliTPCCalibViewer::MakeTree("QAtest.root", fArray,0);
   //
+  // 2.) Check invariants
+  //
   TFile*fout= TFile::Open("QAtest.root");
   TTree * tree  = (TTree*)fout->Get("calPads");
-  
+  Int_t isOutM0 = tree->Draw("(Ly2.fElements-Ly0.fElements)>>his0(100,-10,10)","abs((Ly2.fElements-Ly0.fElements))>2","goff");
+  Int_t isOutM1=tree->Draw("(Lx2.fElements-Lx0.fElements)/0.75>>his1(100,-10,10)","abs((Lx2.fElements-Lx0.fElements))>0","goff");
+  printf("IsOut=%d\t%d\n",isOutM0,isOutM1);
+  if ((isOutM0+isOutM1)==0) ::Info("UnitTestAliTPCCalPadTree","MedianTest OK");
+  if (isOutM0||isOutM1) ::Fatal("UnitTestAliTPCCalPadTree","MedianTest FAILED");
+  //
+  Int_t isOutL0 = tree->Draw("(LLy2.fElements-Ly0.fElements)>>his0(100,-10,10)","abs((LLy2.fElements-LLy0.fElements))>0","goff");
+  Int_t isOutL1=tree->Draw("(LLx2.fElements-Lx0.fElements)/0.75>>his1(100,-10,10)","abs((LLx2.fElements-LLx0.fElements))>0","goff");
+  printf("IsOut=%d\t%d\n",isOutL0,isOutL1);
+  if ((isOutL0+isOutL1)==0) ::Info("UnitTestAliTPCCalPadTree","LTMTest OK");
+  if (isOutL0||isOutL1) ::Fatal("UnitTestAliTPCCalPadTree","LTMTest FAILED");
 }