X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCCalROC.cxx;h=5928d9df70f13672b50e2439267b4f7924da3270;hb=36660aadb0ce6be85d5a29aa0a06e8b15ba0db5b;hp=c8810c6a684a5927df3c4e6a8826aff997621fb7;hpb=43b569c61a8b50571e412584e82c9f33742bbedf;p=u%2Fmrichter%2FAliRoot.git

diff --git a/TPC/AliTPCCalROC.cxx b/TPC/AliTPCCalROC.cxx
index c8810c6a684..5928d9df70f 100644
--- a/TPC/AliTPCCalROC.cxx
+++ b/TPC/AliTPCCalROC.cxx
@@ -16,29 +16,39 @@
 
 ///////////////////////////////////////////////////////////////////////////////
 //                                                                           //
-//  Calibration base class for a single ROC                                  //
-//  Contains one float value per pad                                         //
+//     Calibration base class for a single ROC                               //
+//     Contains one float value per pad                                      //
 //     mapping of the pads taken form AliTPCROC                              //
 //                                                                           //
 ///////////////////////////////////////////////////////////////////////////////
 
-#include "AliTPCCalROC.h"
+//
+// ROOT includes 
+//
 #include "TMath.h"
 #include "TClass.h"
 #include "TFile.h"
 #include "TH1F.h"
 #include "TH2F.h"
+#include "TArrayI.h"
+
+//
+//
+#include "AliTPCCalROC.h"
 #include "AliMathBase.h"
+
+#include "TRandom3.h"      // only needed by the AliTPCCalROCTest() method
+
 ClassImp(AliTPCCalROC)
 
 
 //_____________________________________________________________________________
 AliTPCCalROC::AliTPCCalROC()
-             :TObject(),
+             :TNamed(),
 	      fSector(0),
 	      fNChannels(0),
 	      fNRows(0),
-	      fIndexes(0),
+	      fkIndexes(0),
 	      fData(0)
 {
   //
@@ -49,11 +59,11 @@ AliTPCCalROC::AliTPCCalROC()
 
 //_____________________________________________________________________________
 AliTPCCalROC::AliTPCCalROC(UInt_t  sector)
-             :TObject(),
+             :TNamed(),
 	      fSector(0),
 	      fNChannels(0),
 	      fNRows(0),
-	      fIndexes(0),
+	      fkIndexes(0),
 	      fData(0)
 {
   //
@@ -62,18 +72,18 @@ AliTPCCalROC::AliTPCCalROC(UInt_t  sector)
   fSector = sector;
   fNChannels    =  AliTPCROC::Instance()->GetNChannels(fSector);
   fNRows        =  AliTPCROC::Instance()->GetNRows(fSector);
-  fIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
+  fkIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
   fData = new Float_t[fNChannels];
   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata] = 0.;
 }
 
 //_____________________________________________________________________________
 AliTPCCalROC::AliTPCCalROC(const AliTPCCalROC &c)
-             :TObject(c),
+             :TNamed(c),
 	      fSector(0),
 	      fNChannels(0),
 	      fNRows(0),
-	      fIndexes(0),
+	      fkIndexes(0),
 	      fData(0)
 {
   //
@@ -82,7 +92,7 @@ AliTPCCalROC::AliTPCCalROC(const AliTPCCalROC &c)
   fSector = c.fSector;
   fNChannels    =  AliTPCROC::Instance()->GetNChannels(fSector);
   fNRows        =  AliTPCROC::Instance()->GetNRows(fSector);
-  fIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
+  fkIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
   //
   fData   = new Float_t[fNChannels];
   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata] = c.fData[idata];
@@ -114,38 +124,40 @@ AliTPCCalROC::~AliTPCCalROC()
 
 void AliTPCCalROC::Streamer(TBuffer &R__b)
 {
+   //
    // Stream an object of class AliTPCCalROC.
+   //
    if (R__b.IsReading()) {
       AliTPCCalROC::Class()->ReadBuffer(R__b, this);
-      fIndexes =  AliTPCROC::Instance()->GetRowIndexes(fSector);
+      fkIndexes =  AliTPCROC::Instance()->GetRowIndexes(fSector);
    } else {
       AliTPCCalROC::Class()->WriteBuffer(R__b,this);
    }
 }
 
-//  //
-//   // algebra
-//   void Add(Float_t c1);
-//   void Multiply(Float_t c1);
-//   void Add(const AliTPCCalROC * roc, Double_t c1 = 1);
-//   void Divide(const AliTPCCalROC * roc);   
+
+// algebra fuctions:
 
 void AliTPCCalROC::Add(Float_t c1){
   //
-  // add constant
+  // add c1 to each channel of the ROC  
   //
   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata]+=c1;
 }
+
+
 void AliTPCCalROC::Multiply(Float_t c1){
   //
-  // add constant
+  // multiply each channel of the ROC with c1
   //
   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata]*=c1;
 }
 
+
 void AliTPCCalROC::Add(const AliTPCCalROC * roc, Double_t c1){
   //
-  // add values 
+  // multiply AliTPCCalROC roc by c1 and add each channel to the coresponing channel in the ROC
+  //  - pad by pad 
   //
   for (UInt_t  idata = 0; idata< fNChannels; idata++){
     fData[idata]+=roc->fData[idata]*c1;
@@ -155,7 +167,8 @@ void AliTPCCalROC::Add(const AliTPCCalROC * roc, Double_t c1){
 
 void AliTPCCalROC::Multiply(const AliTPCCalROC*  roc) {
   //
-  // multiply values - per by pad
+  // multiply each channel of the ROC with the coresponding channel of 'roc'
+  //     - pad by pad -
   //
   for (UInt_t  idata = 0; idata< fNChannels; idata++){
     fData[idata]*=roc->fData[idata];
@@ -165,7 +178,8 @@ void AliTPCCalROC::Multiply(const AliTPCCalROC*  roc) {
 
 void AliTPCCalROC::Divide(const AliTPCCalROC*  roc) {
   //
-  // divide values 
+  // divide each channel of the ROC by the coresponding channel of 'roc'
+  //     - pad by pad -
   //
   Float_t kEpsilon=0.00000000000000001;
   for (UInt_t  idata = 0; idata< fNChannels; idata++){
@@ -174,21 +188,96 @@ void AliTPCCalROC::Divide(const AliTPCCalROC*  roc) {
   }
 }
 
+Double_t AliTPCCalROC::GetMean(AliTPCCalROC *const outlierROC) const {
+   //
+   //  returns the mean value of the ROC
+   //  pads with value != 0 in outlierROC are not used for the calculation
+   //  return 0 if no data is accepted by the outlier cuts 
+   //
+   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[i];
+         nPoints++;
+      }
+   }
+   Double_t mean = 0;
+   if(nPoints>0)
+     mean = TMath::Mean(nPoints, ddata);
+   delete [] ddata;
+   return mean;
+}
 
+Double_t AliTPCCalROC::GetMedian(AliTPCCalROC *const outlierROC) const {
+   //
+   //  returns the median value of the ROC
+   //  pads with value != 0 in outlierROC are not used for the calculation
+   //  return 0 if no data is accepted by the outlier cuts 
+   //
+   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[i];
+         nPoints++;
+      }
+   }
+   Double_t median = 0;
+   if(nPoints>0)
+     median = TMath::Median(nPoints, ddata);
+   delete [] ddata;
+   return median;
+}
 
+Double_t AliTPCCalROC::GetRMS(AliTPCCalROC *const outlierROC) const {
+   //
+   //  returns the RMS value of the ROC
+   //  pads with value != 0 in outlierROC are not used for the calculation
+   //  return 0 if no data is accepted by the outlier cuts 
+   //
+   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[i];
+         nPoints++;
+      }
+   }
+   Double_t rms = 0;
+   if(nPoints>0)
+     rms = TMath::RMS(nPoints, ddata);
+   delete [] ddata;
+   return rms;
+}
 
-Double_t AliTPCCalROC::GetLTM(Double_t *sigma, Double_t fraction){
+Double_t AliTPCCalROC::GetLTM(Double_t *const sigma, Double_t fraction, AliTPCCalROC *const outlierROC){
   //
-  //  Calculate LTM mean and sigma
+  //  returns the LTM and sigma
+  //  pads with value != 0 in outlierROC are not used for the calculation
+   //  return 0 if no data is accepted by the outlier cuts 
   //
   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);
-  if (sigma) *sigma=lsigma;
+  UInt_t nPoints = 0;
+  for (UInt_t i=0;i<fNChannels;i++) {
+     if (!outlierROC || !(outlierROC->GetValue(i))) {
+        ddata[nPoints]= fData[i];
+        nPoints++;
+     }
+  }
+
+  Double_t ltm =0, lsigma=0;
+  if(nPoints>0) {
+    Int_t hh = TMath::Min(TMath::Nint(fraction *nPoints), Int_t(nPoints));
+    AliMathBase::EvaluateUni(nPoints,ddata, ltm, lsigma, hh);
+    if (sigma) *sigma=lsigma;
+  }
+  
   delete [] ddata;
-  return mean;
+  return ltm;
 }
 
 TH1F * AliTPCCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type){
@@ -197,6 +286,7 @@ TH1F * AliTPCCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type){
   // type -1 = user defined range
   //       0 = nsigma cut nsigma=min
   //       1 = delta cut around median delta=min
+  //
   if (type>=0){
     if (type==0){
       // nsigma range
@@ -229,7 +319,7 @@ TH1F * AliTPCCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type){
   TH1F * his = new TH1F(name,name,100, min,max);
   for (UInt_t irow=0; irow<fNRows; irow++){
     UInt_t npads = (Int_t)GetNPads(irow);
-    for (UInt_t ipad=0; ipad<=npads; ipad++){
+    for (UInt_t ipad=0; ipad<npads; ipad++){
       his->Fill(GetValue(irow,ipad));
     }
   }
@@ -244,6 +334,7 @@ TH2F * AliTPCCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type){
   // type -1 = user defined range
   //       0 = nsigma cut nsigma=min
   //       1 = delta cut around median delta=min
+  //
   if (type>=0){
     if (type==0){
       // nsigma range
@@ -278,7 +369,7 @@ TH2F * AliTPCCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type){
   TH2F * his = new TH2F(name,name,fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
   for (UInt_t irow=0; irow<fNRows; irow++){
     UInt_t npads = (Int_t)GetNPads(irow);
-    for (UInt_t ipad=0; ipad<=npads; ipad++){
+    for (UInt_t ipad=0; ipad<npads; ipad++){
       his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,GetValue(irow,ipad));
     }
   }
@@ -295,6 +386,7 @@ TH2F * AliTPCCalROC::MakeHistoOutliers(Float_t delta, Float_t fraction, Int_t ty
   // fraction - fraction of values used to define sigma
   // delta - in mode 0 - nsigma cut
   //            mode 1 - delta cut
+  //
   Double_t sigma;
   Float_t mean  = GetLTM(&sigma,fraction);  
   if (type==0) delta*=sigma; 
@@ -308,7 +400,7 @@ TH2F * AliTPCCalROC::MakeHistoOutliers(Float_t delta, Float_t fraction, Int_t ty
   TH2F * his = new TH2F(name,name,fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
   for (UInt_t irow=0; irow<fNRows; irow++){
     UInt_t npads = (Int_t)GetNPads(irow);
-    for (UInt_t ipad=0; ipad<=npads; ipad++){
+    for (UInt_t ipad=0; ipad<npads; ipad++){
       if (TMath::Abs(GetValue(irow,ipad)-mean)>delta)
 	his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,1);
     }
@@ -338,10 +430,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++){
@@ -349,30 +445,394 @@ 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;
+}
+
+
+AliTPCCalROC * AliTPCCalROC::LocalFit(Int_t rowRadius, Int_t padRadius, AliTPCCalROC* ROCoutliers, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction) {
+  //
+  // MakeLocalFit - smoothing
+  // returns a AliTPCCalROC with smoothed data
+  // rowRadius and padRadius specifies a window around a given pad. 
+  // The data of this window are fitted with a parabolic function. 
+  // This function is evaluated at the pad's position.
+  // At the edges the window is shifted, so that the specified pad is not anymore in the center of the window. 
+  // rowRadius  -  radius - rows to be used for smoothing
+  // padradius  -  radius - pads to be used for smoothing
+  // ROCoutlier -  map of outliers - pads not to be used for local smoothing
+  // robust     -  robust method of fitting  - (much slower)
+  // chi2Threshold: Threshold for chi2 when EvalRobust is called
+  // robustFraction: Fraction of data that will be used in EvalRobust
+  //
+  AliTPCCalROC * xROCfitted = new AliTPCCalROC(fSector);
+  TLinearFitter fitterQ(6,"hyp5");
+   // TLinearFitter fitterQ(6,"x0++x1++x2++x3++x4++x5");
+  fitterQ.StoreData(kTRUE);
+  for (UInt_t row=0; row < GetNrows(); row++) {
+    //std::cout << "Entering row " << row << " of " << GetNrows() << " @ sector "<< fSector << " for local fitting... "<< std::endl;
+    for (UInt_t pad=0; pad < GetNPads(row); pad++)
+      xROCfitted->SetValue(row, pad, GetNeighbourhoodValue(&fitterQ, row, pad, rowRadius, padRadius, ROCoutliers, robust, chi2Threshold, robustFraction));
+  }
+  return xROCfitted;
+}
+
+
+Double_t AliTPCCalROC::GetNeighbourhoodValue(TLinearFitter* fitterQ, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius, AliTPCCalROC *const ROCoutliers, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction) {
+  //
+  // AliTPCCalROC::GetNeighbourhoodValue - smoothing - PRIVATE
+  // in this function the fit for LocalFit is done
+  //
+
+  fitterQ->ClearPoints();
+  TVectorD fitParam(6);
+  Int_t    npoints=0;
+  Double_t xx[6];
+  Float_t dlx, dly;
+  Float_t lPad[3] = {0};
+  Float_t localXY[3] = {0};
+  
+  AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
+  tpcROCinstance->GetPositionLocal(fSector, row, pad, lPad);  // calculate position lPad by pad and row number
+  
+  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;   // 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];
+    //xx[0] = 1;
+    xx[1] = dlx;
+    xx[2] = dly;
+    xx[3] = dlx*dlx;
+    xx[4] = dly*dly;
+    xx[5] = dlx*dly;
+    if (!ROCoutliers || ROCoutliers->GetValue(r,p) != 1) {
+      fitterQ->AddPoint(&xx[1], 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
+  }
+  fitterQ->Eval();
+  fitterQ->GetParameters(fitParam);
+  Float_t chi2Q = 0;
+  if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
+  //if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
+  if (robust && chi2Q > chi2Threshold) {
+    //std::cout << "robust fitter called... " << std::endl;
+    fitterQ->EvalRobust(robustFraction);
+    fitterQ->GetParameters(fitParam);
+  }
+  Double_t value = fitParam[0];
+  
+  //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;
+}
+
+
+
+
+void AliTPCCalROC::GetNeighbourhood(TArrayI* &rowArray, TArrayI* &padArray, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius) {
+  //
+  // AliTPCCalROC::GetNeighbourhood - PRIVATE
+  // in this function the window for LocalFit is determined
+  //
+  rowArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
+  padArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
+  
+  Int_t rmin = row - rRadius;
+  UInt_t rmax = row + rRadius;
+  
+  // if window goes out of ROC
+  if (rmin < 0) {
+    rmax = rmax - rmin;
+    rmin = 0;
+  }
+  if (rmax >= GetNrows()) {
+    rmin = rmin - (rmax - GetNrows()+1);
+    rmax = GetNrows() - 1;
+      if (rmin  < 0 ) rmin = 0; // if the window is bigger than the ROC
+  }
+  
+  Int_t pmin, pmax;
+  Int_t i = 0;
+  
+  for (UInt_t r = rmin; r <= rmax; r++) {
+    pmin = pad - pRadius;
+    pmax = pad + pRadius;
+    if (pmin < 0) {
+      pmax = pmax - pmin;
+      pmin = 0;
+    }
+    if (pmax >= (Int_t)GetNPads(r)) {
+      pmin = pmin - (pmax - GetNPads(r)+1);
+      pmax = GetNPads(r) - 1;
+      if (pmin  < 0 ) pmin = 0; // if the window is bigger than the ROC
+    }
+    for (Int_t p = pmin; p <= pmax; 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;
+    (*rowArray)[j] = -1;
+    (*padArray)[j] = -1;
+    //std::cout << "writing -1" << std::endl;
+  }
+}
+
+
+
+void AliTPCCalROC::GlobalFit(const AliTPCCalROC* ROCoutliers, Bool_t robust, TVectorD &fitParam, TMatrixD &covMatrix, Float_t & chi2, Int_t fitType, Double_t chi2Threshold, Double_t robustFraction, Double_t err){
+  //
+  // Makes a  GlobalFit for the given secotr and return fit-parameters, covariance and chi2
+  // The origin of the fit function is the center of the ROC!
+  // fitType == 0: fit plane function
+  // fitType == 1: fit parabolic function
+  // ROCoutliers - pads with value !=0 are not used in fitting procedure
+  // chi2Threshold: Threshold for chi2 when EvalRobust is called
+  // robustFraction: Fraction of data that will be used in EvalRobust
+  // err: error of the data points
+  //
+  TLinearFitter* fitterG = 0;
+  Double_t xx[6];
+  
+  if (fitType  == 1) {
+    fitterG = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
+    fitParam.ResizeTo(6);
+    covMatrix.ResizeTo(6,6);
+  } else {
+    fitterG = new TLinearFitter(3,"x0++x1++x2");
+    fitParam.ResizeTo(3);
+    covMatrix.ResizeTo(3,3);
+  }
+  fitterG->StoreData(kTRUE);   
+  fitterG->ClearPoints();
+  Int_t    npoints=0;
+  
+  Float_t dlx, dly;
+  Float_t centerPad[3] = {0};
+  Float_t localXY[3] = {0};
+  
+  AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
+  tpcROCinstance->GetPositionLocal(fSector, GetNrows()/2, GetNPads(GetNrows()/2)/2, centerPad);  // calculate center of ROC 
+  
+  // loop over all channels and read data into fitterG
+  for (UInt_t irow = 0; irow < GetNrows(); irow++) {
+    for (UInt_t ipad = 0; ipad < GetNPads(irow); ipad++) {
+      // fill fitterG
+      if (ROCoutliers && ROCoutliers->GetValue(irow, ipad) != 0) continue;
+      tpcROCinstance->GetPositionLocal(fSector, irow, ipad, localXY);   // calculate position localXY by pad and row number
+      dlx = localXY[0] - centerPad[0];
+      dly = localXY[1] - centerPad[1];
+      xx[0] = 1;
+      xx[1] = dlx;
+      xx[2] = dly;
+      xx[3] = dlx*dlx;
+      xx[4] = dly*dly;
+      xx[5] = dlx*dly;
+      npoints++;
+      fitterG->AddPoint(xx, GetValue(irow, ipad), err);
+    }
+  }
+  if(npoints>10) { // make sure there is something to fit
+    fitterG->Eval();
+    fitterG->GetParameters(fitParam);
+    fitterG->GetCovarianceMatrix(covMatrix);
+    if (fitType == 1)
+      chi2 = fitterG->GetChisquare()/(npoints-6.);
+    else chi2 = fitterG->GetChisquare()/(npoints-3.);
+    if (robust && chi2 > chi2Threshold) {
+      //    std::cout << "robust fitter called... " << std::endl;
+      fitterG->EvalRobust(robustFraction);
+      fitterG->GetParameters(fitParam);
+    }
+  } else {
+    // set parameteres to 0
+    Int_t nParameters = 3;
+    if (fitType  == 1)
+      nParameters = 6;
+
+    for(Int_t i = 0; i < nParameters; i++)
+      fitParam[i] = 0;
+  }
+  
+  delete fitterG;
+}
+
+
+AliTPCCalROC* AliTPCCalROC::CreateGlobalFitCalROC(TVectorD &fitParam, Int_t sector){
+  //
+  // Create ROC with global fit parameters
+  // The origin of the fit function is the center of the ROC!
+  // loop over all channels, write fit values into new ROC and return it
+  //
+  Float_t dlx, dly;
+  Float_t centerPad[3] = {0};
+  Float_t localXY[3] = {0};
+  AliTPCCalROC * xROCfitted = new AliTPCCalROC(sector);
+  AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
+  tpcROCinstance->GetPositionLocal(sector, xROCfitted->GetNrows()/2, xROCfitted->GetNPads(xROCfitted->GetNrows()/2)/2, centerPad);  // calculate center of ROC 
+  Int_t fitType = 1;
+  if (fitParam.GetNoElements() == 6) fitType = 1;
+  else fitType = 0;
+  Double_t value = 0;
+  if (fitType == 1) { // parabolic fit
+    for (UInt_t irow = 0; irow < xROCfitted->GetNrows(); irow++) {
+      for (UInt_t ipad = 0; ipad < xROCfitted->GetNPads(irow); ipad++) {
+	tpcROCinstance->GetPositionLocal(sector, irow, ipad, localXY);   // calculate position localXY by pad and row number
+	dlx = localXY[0] - centerPad[0];
+	dly = localXY[1] - centerPad[1];
+	value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly + fitParam[3]*dlx*dlx + fitParam[4]*dly*dly + fitParam[5]*dlx*dly;
+	xROCfitted->SetValue(irow, ipad, value);
+      }
+    }   
+  }
+  else {  // linear fit
+    for (UInt_t irow = 0; irow < xROCfitted->GetNrows(); irow++) {
+      for (UInt_t ipad = 0; ipad < xROCfitted->GetNPads(irow); ipad++) {
+	tpcROCinstance->GetPositionLocal(sector, irow, ipad, localXY);   // calculate position localXY by pad and row number
+	dlx = localXY[0] - centerPad[0];
+	dly = localXY[1] - centerPad[1];
+	value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly;
+	xROCfitted->SetValue(irow, ipad, value);
+      }
+    }   
+  }
+  return xROCfitted;
 }