Add outliers map to the statistical functions
authorcvetan <cvetan@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 3 May 2007 16:30:38 +0000 (16:30 +0000)
committercvetan <cvetan@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 3 May 2007 16:30:38 +0000 (16:30 +0000)
(Lars, Stefan)

TPC/AliTPCCalPad.cxx
TPC/AliTPCCalPad.h
TPC/AliTPCCalROC.cxx
TPC/AliTPCCalROC.h

index f62a2cb..27f81d2 100644 (file)
@@ -29,6 +29,8 @@
 #include <TGraph2D.h>
 #include <TH2F.h>
 #include "TTreeStream.h"
+#include "TFile.h"
+#include "TKey.h"
 
 ClassImp(AliTPCCalPad)
 
@@ -259,7 +261,7 @@ Double_t AliTPCCalPad::GetMeanRMS(Double_t &rms)
 
 
 //_____________________________________________________________________________
-Double_t AliTPCCalPad::GetMean()
+Double_t AliTPCCalPad::GetMean(AliTPCCalPad* outlierPad)
 {
     //
     // return mean of the mean of all ROCs
@@ -267,17 +269,19 @@ Double_t AliTPCCalPad::GetMean()
     Double_t arr[kNsec];
     Int_t n=0;
     for (Int_t isec = 0; isec < kNsec; isec++) {
-        AliTPCCalROC *calRoc = fROC[isec];
-       if ( calRoc ){
-           arr[n] = calRoc->GetMean();
-            n++;
-       }
+       AliTPCCalROC *calRoc = fROC[isec];
+       if ( calRoc ){
+          AliTPCCalROC* outlierROC = 0;
+          if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
+              arr[n] = calRoc->GetMean(outlierROC);
+          n++;
+       }
     }
     return TMath::Mean(n,arr);
 }
 
 //_____________________________________________________________________________
-Double_t AliTPCCalPad::GetRMS()
+Double_t AliTPCCalPad::GetRMS(AliTPCCalPad* outlierPad)
 {
     //
     // return mean of the RMS of all ROCs
@@ -285,17 +289,19 @@ Double_t AliTPCCalPad::GetRMS()
     Double_t arr[kNsec];
     Int_t n=0;
     for (Int_t isec = 0; isec < kNsec; isec++) {
-        AliTPCCalROC *calRoc = fROC[isec];
-       if ( calRoc ){
-           arr[n] = calRoc->GetRMS();
-            n++;
-       }
+       AliTPCCalROC *calRoc = fROC[isec];
+       if ( calRoc ){
+          AliTPCCalROC* outlierROC = 0;
+          if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
+          arr[n] = calRoc->GetRMS(outlierROC);
+          n++;
+       }
     }
     return TMath::Mean(n,arr);
 }
 
 //_____________________________________________________________________________
-Double_t AliTPCCalPad::GetMedian()
+Double_t AliTPCCalPad::GetMedian(AliTPCCalPad* outlierPad)
 {
     //
     // return mean of the median of all ROCs
@@ -303,17 +309,19 @@ Double_t AliTPCCalPad::GetMedian()
     Double_t arr[kNsec];
     Int_t n=0;
     for (Int_t isec = 0; isec < kNsec; isec++) {
-        AliTPCCalROC *calRoc = fROC[isec];
-       if ( calRoc ){
-           arr[n] = calRoc->GetMedian();
-            n++;
-       }
+       AliTPCCalROC *calRoc = fROC[isec];
+       if ( calRoc ){
+          AliTPCCalROC* outlierROC = 0;
+          if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
+          arr[n] = calRoc->GetMedian(outlierROC);
+          n++;
+       }
     }
     return TMath::Mean(n,arr);
 }
 
 //_____________________________________________________________________________
-Double_t AliTPCCalPad::GetLTM(Double_t *sigma, Double_t fraction)
+Double_t AliTPCCalPad::GetLTM(Double_t *sigma, Double_t fraction, AliTPCCalPad* outlierPad)
 {
     //
     // return mean of the LTM and sigma of all ROCs
@@ -327,7 +335,9 @@ Double_t AliTPCCalPad::GetLTM(Double_t *sigma, Double_t fraction)
         AliTPCCalROC *calRoc = fROC[isec];
        if ( calRoc ){
            if ( sigma ) sTemp=arrs+n;
-           arrm[n] = calRoc->GetLTM(sTemp,fraction);
+       AliTPCCalROC* outlierROC = 0;
+       if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
+           arrm[n] = calRoc->GetLTM(sTemp,fraction, outlierROC);
             n++;
        }
     }
@@ -418,12 +428,51 @@ TH2F *AliTPCCalPad::MakeHisto2D(Int_t side){
 
 
 
-void AliTPCCalPad::MakeTree(const char * fileName, TObjArray * array) {
+void AliTPCCalPad::MakeTree(const char * fileName, TObjArray * array, const char * mapFileName, AliTPCCalPad* outlierPad, Float_t ltmFraction) {
   //
   // Write tree with all available information
   //
-   TTreeSRedirector cstream(fileName);
    AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
+
+   TObjArray* mapIROCs = 0;
+   TObjArray* mapOROCs = 0;
+   TVectorF *mapIROCArray = 0;
+   TVectorF *mapOROCArray = 0;
+   Int_t mapEntries = 0;
+   TString* mapNames = 0;
+   
+   if (mapFileName) {
+      TFile mapFile(mapFileName, "read");
+      
+      TList* listOfROCs = mapFile.GetListOfKeys();
+      mapEntries = listOfROCs->GetEntries()/2;
+      mapIROCs = new TObjArray(mapEntries*2);
+      mapOROCs = new TObjArray(mapEntries*2);
+      mapIROCArray = new TVectorF[mapEntries];
+      mapOROCArray = new TVectorF[mapEntries];
+      
+      mapNames = new TString[mapEntries];
+      for (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
+         TString ROCname(((TKey*)(listOfROCs->At(ivalue*2)))->GetName());
+         ROCname.Remove(ROCname.Length()-4, 4);
+         mapIROCs->AddAt((AliTPCCalROC*)mapFile.Get((ROCname + "IROC").Data()), ivalue);
+         mapOROCs->AddAt((AliTPCCalROC*)mapFile.Get((ROCname + "OROC").Data()), ivalue);
+         mapNames[ivalue].Append(ROCname);
+      }
+      
+      for (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
+         mapIROCArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(0));
+         mapOROCArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(36));
+      
+         for (UInt_t ichannel = 0; ichannel < tpcROCinstance->GetNChannels(0); ichannel++)
+            (mapIROCArray[ivalue])[ichannel] = ((AliTPCCalROC*)(mapIROCs->At(ivalue)))->GetValue(ichannel);
+         for (UInt_t ichannel = 0; ichannel < tpcROCinstance->GetNChannels(36); ichannel++)
+            (mapOROCArray[ivalue])[ichannel] = ((AliTPCCalROC*)(mapOROCs->At(ivalue)))->GetValue(ichannel);
+      }
+
+   } //  if (mapFileName)
+  
+   TTreeSRedirector cstream(fileName);
    Int_t arrayEntries = array->GetEntries();
    
    TString* names = new TString[arrayEntries];
@@ -438,6 +487,12 @@ void AliTPCCalPad::MakeTree(const char * fileName, TObjArray * array) {
       TVectorF mean(arrayEntries);
       TVectorF rms(arrayEntries);
       TVectorF ltm(arrayEntries);
+      TVectorF ltmrms(arrayEntries);
+      TVectorF medianWithOut(arrayEntries);
+      TVectorF meanWithOut(arrayEntries);
+      TVectorF rmsWithOut(arrayEntries);
+      TVectorF ltmWithOut(arrayEntries);
+      TVectorF ltmrmsWithOut(arrayEntries);
       
       TVectorF *vectorArray = new TVectorF[arrayEntries];
       for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++)
@@ -446,17 +501,35 @@ void AliTPCCalPad::MakeTree(const char * fileName, TObjArray * array) {
       for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
          AliTPCCalPad* calPad = (AliTPCCalPad*) array->At(ivalue);
          AliTPCCalROC* calROC = calPad->GetCalROC(isector);
+         AliTPCCalROC* outlierROC = 0;
+         if (outlierPad) outlierROC = outlierPad->GetCalROC(isector);
          if (calROC) {
             median[ivalue] = calROC->GetMedian();
             mean[ivalue] = calROC->GetMean();
             rms[ivalue] = calROC->GetRMS();
-            ltm[ivalue] = calROC->GetLTM();
+            Double_t ltmrmsValue = 0;
+            ltm[ivalue] = calROC->GetLTM(&ltmrmsValue, ltmFraction);
+            ltmrms[ivalue] = ltmrmsValue;
+            if (outlierROC) {
+               medianWithOut[ivalue] = calROC->GetMedian(outlierROC);
+               meanWithOut[ivalue] = calROC->GetMean(outlierROC);
+               rmsWithOut[ivalue] = calROC->GetRMS(outlierROC);
+               ltmrmsValue = 0;
+               ltmWithOut[ivalue] = calROC->GetLTM(&ltmrmsValue, ltmFraction, outlierROC);
+               ltmrmsWithOut[ivalue] = ltmrmsValue;
+            }
          }
          else {
             median[ivalue] = 0.;
             mean[ivalue] = 0.;
             rms[ivalue] = 0.;
             ltm[ivalue] = 0.;
+            ltmrms[ivalue] = 0.;
+            medianWithOut[ivalue] = 0.;
+            meanWithOut[ivalue] = 0.;
+            rmsWithOut[ivalue] = 0.;
+            ltmWithOut[ivalue] = 0.;
+            ltmrmsWithOut[ivalue] = 0.;
          }
       }
       
@@ -502,7 +575,16 @@ void AliTPCCalPad::MakeTree(const char * fileName, TObjArray * array) {
             (Char_t*)((names[ivalue] + "_Median=").Data()) << median[ivalue] <<
             (Char_t*)((names[ivalue] + "_Mean=").Data()) << mean[ivalue] <<
             (Char_t*)((names[ivalue] + "_RMS=").Data()) << rms[ivalue] <<
-            (Char_t*)((names[ivalue] + "_LTM=").Data()) << ltm[ivalue];
+            (Char_t*)((names[ivalue] + "_LTM=").Data()) << ltm[ivalue] <<
+            (Char_t*)((names[ivalue] + "_RMS_LTM=").Data()) << ltmrms[ivalue];
+         if (outlierPad) {
+            cstream << "calPads" <<
+               (Char_t*)((names[ivalue] + "_Median_OutlierCutted=").Data()) << medianWithOut[ivalue] <<
+               (Char_t*)((names[ivalue] + "_Mean_OutlierCutted=").Data()) << meanWithOut[ivalue] <<
+               (Char_t*)((names[ivalue] + "_RMS_OutlierCutted=").Data()) << rmsWithOut[ivalue] <<
+               (Char_t*)((names[ivalue] + "_LTM_OutlierCutted=").Data()) << ltmWithOut[ivalue] <<
+               (Char_t*)((names[ivalue] + "_RMS_LTM_OutlierCutted=").Data()) << ltmrmsWithOut[ivalue];
+         }
       }
 
       for  (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
@@ -510,6 +592,17 @@ void AliTPCCalPad::MakeTree(const char * fileName, TObjArray * array) {
             (Char_t*)((names[ivalue] + ".=").Data()) << &vectorArray[ivalue];
       }
 
+      if (mapFileName) {
+         for  (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
+            if (isector < 36)
+               cstream << "calPads" <<
+                  (Char_t*)((mapNames[ivalue] + ".=").Data()) << &mapIROCArray[ivalue];
+            else
+               cstream << "calPads" <<
+                  (Char_t*)((mapNames[ivalue] + ".=").Data()) << &mapOROCArray[ivalue];
+         }
+      }
+
       cstream << "calPads" <<
          "row.=" << &posArray[0] <<
          "pad.=" << &posArray[1] <<
@@ -525,6 +618,13 @@ void AliTPCCalPad::MakeTree(const char * fileName, TObjArray * array) {
       delete[] vectorArray;
    }
    delete[] names;
+   if (mapFileName) {
+      delete mapIROCs;
+      delete mapOROCs;
+      delete[] mapIROCArray;
+      delete[] mapOROCArray;
+      delete[] mapNames;
+   }
 }
 
 
index 86fc7fd..095d143 100644 (file)
@@ -40,14 +40,14 @@ class AliTPCCalPad : public TNamed {
   void Divide(const AliTPCCalPad * pad);
   //
   Double_t GetMeanRMS(Double_t &rms);
-  Double_t GetMean();
-  Double_t GetRMS() ;
-  Double_t GetMedian() ;
-  Double_t GetLTM(Double_t *sigma=0, Double_t fraction=0.9);
+  Double_t GetMean(AliTPCCalPad* outlierPad = 0);
+  Double_t GetRMS(AliTPCCalPad* outlierPad = 0) ;
+  Double_t GetMedian(AliTPCCalPad* outlierPad = 0) ;
+  Double_t GetLTM(Double_t *sigma=0, Double_t fraction=0.9, AliTPCCalPad* outlierPad = 0);
   TGraph  *MakeGraph(Int_t type=0, Float_t ratio=0.7);
   TH2F    *MakeHisto2D(Int_t side=0);
   TH1F    *MakeHisto1D(Float_t min=4, Float_t max=-4, Int_t type=0);  
-  static void MakeTree(const char * fileName, TObjArray * array);
+  static void MakeTree(const char * fileName, TObjArray * array, const char * mapFileName = 0, AliTPCCalPad* outlierPad = 0, Float_t ltmFraction = 0.9);
  protected:
   AliTPCCalROC *fROC[kNsec];                    //  Array of ROC objects which contain the values per pad
   ClassDef(AliTPCCalPad,1)                      //  TPC calibration class for parameters which are saved per pad
index 831c3e1..f8e23dc 100644 (file)
@@ -36,6 +36,8 @@
 #include "AliTPCCalROC.h"
 #include "AliMathBase.h"
 
+#include "TRandom3.h"      // only needed by the AliTPCCalROCTest() method
+
 ClassImp(AliTPCCalROC)
 
 
@@ -181,18 +183,66 @@ void AliTPCCalROC::Divide(const AliTPCCalROC*  roc) {
   }
 }
 
+Double_t AliTPCCalROC::GetMean(AliTPCCalROC* outlierROC) {
+   if (!outlierROC) return TMath::Mean(fNChannels, fData);
+   Double_t *ddata = new Double_t[fNChannels];
+   Int_t NPoints = 0;
+   for (UInt_t i=0;i<fNChannels;i++) {
+      if (!(outlierROC->GetValue(i))) {
+         ddata[NPoints]= fData[NPoints];
+         NPoints++;
+      }
+   }
+   Double_t mean = TMath::Mean(NPoints, ddata);
+   delete [] ddata;
+   return mean;
+}
 
+Double_t AliTPCCalROC::GetMedian(AliTPCCalROC* outlierROC) {
+   if (!outlierROC) return TMath::Median(fNChannels, fData);
+   Double_t *ddata = new Double_t[fNChannels];
+   Int_t NPoints = 0;
+   for (UInt_t i=0;i<fNChannels;i++) {
+      if (!(outlierROC->GetValue(i))) {
+         ddata[NPoints]= fData[NPoints];
+         NPoints++;
+      }
+   }
+   Double_t mean = TMath::Median(NPoints, ddata);
+   delete [] ddata;
+   return mean;
+}
 
+Double_t AliTPCCalROC::GetRMS(AliTPCCalROC* outlierROC) {
+   if (!outlierROC) return TMath::RMS(fNChannels, fData);
+   Double_t *ddata = new Double_t[fNChannels];
+   Int_t NPoints = 0;
+   for (UInt_t i=0;i<fNChannels;i++) {
+      if (!(outlierROC->GetValue(i))) {
+         ddata[NPoints]= fData[NPoints];
+         NPoints++;
+      }
+   }
+   Double_t mean = TMath::RMS(NPoints, ddata);
+   delete [] ddata;
+   return mean;
+}
 
-Double_t AliTPCCalROC::GetLTM(Double_t *sigma, Double_t fraction){
+Double_t AliTPCCalROC::GetLTM(Double_t *sigma, Double_t fraction, AliTPCCalROC* outlierROC){
   //
   //  Calculate LTM mean and sigma
   //
   Double_t *ddata = new Double_t[fNChannels];
   Double_t mean=0, lsigma=0;
-  Int_t hh = TMath::Min(TMath::Nint(fraction *fNChannels), Int_t(fNChannels));
-  for (UInt_t i=0;i<fNChannels;i++) ddata[i]= fData[i];
-  AliMathBase::EvaluateUni(UInt_t(fNChannels),ddata, mean, lsigma, hh);
+  UInt_t NPoints = 0;
+  for (UInt_t i=0;i<fNChannels;i++) {
+     if (!outlierROC || !(outlierROC->GetValue(i))) {
+        ddata[NPoints]= fData[NPoints];
+        NPoints++;
+     }
+  }
+  Int_t hh = TMath::Min(TMath::Nint(fraction *NPoints), Int_t(NPoints));
+  AliMathBase::EvaluateUni(NPoints,ddata, mean, lsigma, hh);
   if (sigma) *sigma=lsigma;
   delete [] ddata;
   return mean;
@@ -345,10 +395,14 @@ void AliTPCCalROC::Draw(Option_t* opt){
 
 
 
-void AliTPCCalROC::Test(){
+void AliTPCCalROC::Test() {
   //
-  // example function to show functionality and tes AliTPCCalROC
+  // example function to show functionality and test AliTPCCalROC
   //
+
+  Float_t kEpsilon=0.00001;
+  
+  // create CalROC with defined values
   AliTPCCalROC  roc0(0);  
   for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
     for (UInt_t ipad = 0; ipad <roc0.GetNPads(irow); ipad++){
@@ -356,31 +410,124 @@ void AliTPCCalROC::Test(){
       roc0.SetValue(irow,ipad,value);
     }
   }
-  //
+  
+  // copy CalROC, readout values and compare to original
   AliTPCCalROC roc1(roc0);
   for (UInt_t irow = 0; irow <roc1.GetNrows(); irow++){
     for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
       Float_t value  = irow+ipad/1000.;
       if (roc1.GetValue(irow,ipad)!=value){
-       printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
+        printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
       }
     }
-  }  
+  }
+
+  // write original CalROC to file
   TFile f("calcTest.root","recreate");
   roc0.Write("Roc0");
   AliTPCCalROC * roc2 = (AliTPCCalROC*)f.Get("Roc0");
   f.Close();
-  //
+  
+  // read CalROC from file and compare to original CalROC
   for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
     if (roc0.GetNPads(irow)!=roc2->GetNPads(irow))
       printf("NPads - Read/Write error\trow=%d\n",irow);
     for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
       Float_t value  = irow+ipad/1000.;
       if (roc2->GetValue(irow,ipad)!=value){
-       printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
+        printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
       }
     }
-  }   
+  }
+
+  // 
+  // Algebra Tests
+  //
+  
+  // Add constant
+  AliTPCCalROC roc3(roc0);
+  roc3.Add(1.5);
+  for (UInt_t irow = 0; irow <roc3.GetNrows(); irow++){
+    for (UInt_t ipad = 0; ipad <roc3.GetNPads(irow); ipad++){
+      Float_t value  = irow+ipad/1000. + 1.5;
+      if (TMath::Abs(roc3.GetValue(irow,ipad)-value) > kEpsilon){
+        printf("Add constant - error\trow=%d\tpad=%d\n",irow,ipad);
+      }
+    }
+  }
+
+  // Add another CalROC
+  AliTPCCalROC roc4(roc0);
+  roc4.Add(&roc0, -1.5);
+  for (UInt_t irow = 0; irow <roc4.GetNrows(); irow++){
+    for (UInt_t ipad = 0; ipad <roc4.GetNPads(irow); ipad++){
+      Float_t value  = irow+ipad/1000. - 1.5 * (irow+ipad/1000.);
+      if (TMath::Abs(roc4.GetValue(irow,ipad)-value) > kEpsilon){
+        printf("Add CalROC - error\trow=%d\tpad=%d\n",irow,ipad);
+      }
+    }
+  }
+  
+  // Multiply with constant
+  AliTPCCalROC roc5(roc0);
+  roc5.Multiply(-1.4);
+  for (UInt_t irow = 0; irow <roc5.GetNrows(); irow++){
+    for (UInt_t ipad = 0; ipad <roc5.GetNPads(irow); ipad++){
+      Float_t value  = (irow+ipad/1000.) * (-1.4);
+      if (TMath::Abs(roc5.GetValue(irow,ipad)-value) > kEpsilon){
+        printf("Multiply with constant - error\trow=%d\tpad=%d\n",irow,ipad);
+      }
+    }
+  }
+
+  // Multiply another CalROC
+  AliTPCCalROC roc6(roc0);
+  roc6.Multiply(&roc0);
+  for (UInt_t irow = 0; irow <roc6.GetNrows(); irow++){
+    for (UInt_t ipad = 0; ipad <roc6.GetNPads(irow); ipad++){
+      Float_t value  = (irow+ipad/1000.) * (irow+ipad/1000.);
+      if (TMath::Abs(roc6.GetValue(irow,ipad)-value) > kEpsilon*100){
+        printf("Multiply with CalROC - error\trow=%d\tpad=%d\n",irow,ipad);
+      }
+    }
+  }
+
+
+  // Divide by CalROC
+  AliTPCCalROC roc7(roc0);
+  roc7.Divide(&roc0);
+  for (UInt_t irow = 0; irow <roc7.GetNrows(); irow++){
+    for (UInt_t ipad = 0; ipad <roc7.GetNPads(irow); ipad++){
+      Float_t value  = 1.;
+      if (irow+ipad == 0) value = roc0.GetValue(irow,ipad);
+      if (TMath::Abs(roc7.GetValue(irow,ipad)-value) > kEpsilon){
+        printf("Multiply with CalROC - error\trow=%d\tpad=%d\n",irow,ipad);
+      }
+    }
+  }
+
+  //
+  // Statistics Test
+  //
+  
+  // create CalROC with defined values
+  TRandom3 rnd(0);
+  AliTPCCalROC sroc0(0);
+  for (UInt_t ichannel = 0; ichannel < sroc0.GetNchannels(); ichannel++){
+    Float_t value  = rnd.Gaus(10., 2.);
+    sroc0.SetValue(ichannel,value);
+  }
+
+  printf("Mean   (should be close to 10): %f\n", sroc0.GetMean());
+  printf("RMS    (should be close to  2): %f\n", sroc0.GetRMS());
+  printf("Median (should be close to 10): %f\n", sroc0.GetMedian());
+  printf("LTM    (should be close to 10): %f\n", sroc0.GetLTM());
+
+  //AliTPCCalROC* sroc1 = sroc0.LocalFit(4, 8);
+  
+  //delete sroc1;
+
+//        std::cout << TMath::Abs(roc5.GetValue(irow,ipad)-value) << std::endl;
 }
 
 
@@ -427,13 +574,16 @@ Double_t AliTPCCalROC::GetNeighbourhoodValue(TLinearFitter* fitterQ, Int_t row,
   
   TArrayI *neighbourhoodRows = 0;
   TArrayI *neighbourhoodPads = 0;
+  
+  //std::cerr << "Trying to get neighbourhood for row " << row << ", pad " << pad << std::endl;
   GetNeighbourhood(neighbourhoodRows, neighbourhoodPads, row, pad, rRadius, pRadius);
+  //std::cerr << "Got neighbourhood for row " << row << ", pad " << pad << std::endl;
   
   Int_t r, p;
   for (Int_t i=0; i < (2*rRadius+1)*(2*pRadius+1); i++) {
     r = neighbourhoodRows->At(i);
     p = neighbourhoodPads->At(i);
-    if (r == -1 || p == -1) continue;
+    if (r == -1 || p == -1) continue;   // window is bigger than ROC
     tpcROCinstance->GetPositionLocal(fSector, r, p, localXY);   // calculate position localXY by pad and row number
     dlx = lPad[0] - localXY[0];
     dly = lPad[1] - localXY[1];
@@ -443,11 +593,15 @@ Double_t AliTPCCalROC::GetNeighbourhoodValue(TLinearFitter* fitterQ, Int_t row,
     xx[3] = dlx*dlx;
     xx[4] = dly*dly;
     xx[5] = dlx*dly;
-    if (ROCoutliers && ROCoutliers->GetValue(r,p) != 1) {
+    if (!ROCoutliers || ROCoutliers->GetValue(r,p) != 1) {
       fitterQ->AddPoint(xx, GetValue(r, p), 1);
       npoints++;
     }
   }
+  
+  delete neighbourhoodRows;
+  delete neighbourhoodPads;
+  
   if (npoints < 0.5 * ((2*rRadius+1)*(2*pRadius+1)) ) {
     // std::cerr << "Too few data points for fitting @ row " << row << ", pad " << pad << " in sector " << fSector << std::endl;
     return 0.;  // for diagnostic
@@ -464,9 +618,6 @@ Double_t AliTPCCalROC::GetNeighbourhoodValue(TLinearFitter* fitterQ, Int_t row,
   }
   Double_t value = fitParam[0];
   
-  delete neighbourhoodRows;
-  delete neighbourhoodPads;
-  
   //if (value < 0) std::cerr << "negative fit-value " << value << " in sector "<< this->fSector << " @ row: " << row << " and pad: " << pad << ", with fitter Chi2 = " << chi2Q <<  std::endl;
   
   return value;
@@ -481,8 +632,8 @@ void AliTPCCalROC::GetNeighbourhood(TArrayI* &rowArray, TArrayI* &padArray, Int_
   //
   rowArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
   padArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
-  Int_t* rowArrayTemp = rowArray->GetArray();
-  Int_t* padArrayTemp = padArray->GetArray();
+  //Int_t* rowArrayTemp = rowArray->GetArray();
+  //Int_t* padArrayTemp = padArray->GetArray();
   
   Int_t rmin = row - rRadius;
   UInt_t rmax = row + rRadius;
@@ -514,17 +665,17 @@ void AliTPCCalROC::GetNeighbourhood(TArrayI* &rowArray, TArrayI* &padArray, Int_
       if (pmin  < 0 ) pmin = 0; // if the window is bigger than the ROC
     }
     for (Int_t p = pmin; p <= pmax; p++) {
-      rowArrayTemp[i] = r;
-      padArrayTemp[i] = p;
+      (*rowArray)[i] = r;
+      (*padArray)[i] = p;
       i++;
     }
   }
   for (Int_t j = i; j < rowArray->GetSize(); j++){  // unused padArray-entries, in the case that the window is bigger than the ROC
     //std::cout << "trying to write -1" << std::endl;
-    rowArrayTemp[j] = -1;
-    padArrayTemp[j] = -1;
+    (*rowArray)[j] = -1;
+    (*padArray)[j] = -1;
     //std::cout << "writing -1" << std::endl;
-  } 
+  }
 }
 
 
@@ -571,7 +722,7 @@ void AliTPCCalROC::GlobalFit(const AliTPCCalROC* ROCoutliers, Bool_t robust, TVe
        xx[3] = dlx*dlx;
        xx[4] = dly*dly;
        xx[5] = dlx*dly;
-       if (ROCoutliers && ROCoutliers->GetValue(irow, ipad) != 1) {
+       if (!ROCoutliers || ROCoutliers->GetValue(irow, ipad) != 1) {
            npoints++;
           fitterG->AddPoint(xx, GetValue(irow, ipad), 1);  
         }
@@ -590,7 +741,7 @@ void AliTPCCalROC::GlobalFit(const AliTPCCalROC* ROCoutliers, Bool_t robust, TVe
        xx[0] = 1;
        xx[1] = dlx;
        xx[2] = dly;
-       if (ROCoutliers && ROCoutliers->GetValue(irow, ipad) != 1) {
+       if (!ROCoutliers || ROCoutliers->GetValue(irow, ipad) != 1) {
            npoints++;
           fitterG->AddPoint(xx, GetValue(irow, ipad), 1);  
         }
index 2939c76..c633868 100644 (file)
@@ -46,10 +46,10 @@ class AliTPCCalROC : public TObject {
   void Divide(const AliTPCCalROC * roc);   
   // statistic
   //
-  Double_t GetMean(){return TMath::Mean(fNChannels, fData);}
-  Double_t GetRMS() {return TMath::RMS(fNChannels, fData);}
-  Double_t GetMedian() {return TMath::Median(fNChannels, fData);}
-  Double_t GetLTM(Double_t *sigma=0, Double_t fraction=0.9);
+  Double_t GetMean(AliTPCCalROC* outlierROC = 0);
+  Double_t GetRMS(AliTPCCalROC* outlierROC = 0);
+  Double_t GetMedian(AliTPCCalROC* outlierROC = 0) ;
+  Double_t GetLTM(Double_t *sigma=0, Double_t fraction=0.9, AliTPCCalROC* outlierROC = 0);
   TH1F * MakeHisto1D(Float_t min=4, Float_t max=-4, Int_t type=0);     
   TH2F * MakeHisto2D(Float_t min=4, Float_t max=-4, Int_t type=0);   
   TH2F * MakeHistoOutliers(Float_t delta=4, Float_t fraction=0.7, Int_t mode=0);