]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCcalibTracks.cxx
Print override: changed fot TObjArray in root v5-22-00. Now adapted.
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibTracks.cxx
index 0f61b56816d691a24369b42b39ed1c8866f30ffe..b8396289211cd60852e61110a4122b9f44ac4692 100644 (file)
                Offline/HLT             Offline/HLT                    OCDB entries (AliTPCClusterParam) 
 */            
 
+/*
+
+How to retrive it from file (created using calibration task):
+
+gSystem->Load("libANALYSIS");
+gSystem->Load("libTPCcalib");
+TFile fcalib("CalibObjects.root");
+TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
+AliTPCcalibTracks * calibTracks = ( AliTPCcalibTracks *)array->FindObject("calibTracks");
+
+
+//USAGE of debug stream example
+ gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
+  gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
+  AliXRDPROOFtoolkit tool;
+  TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
+  chainres->Lookup();
+*/
+
+
                                                                //
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -68,6 +88,7 @@ using namespace std;
 #include <TCollection.h>
 #include <iostream>
 #include <TLinearFitter.h>
+#include <TString.h>
 
 //
 // AliROOT includes 
@@ -93,11 +114,9 @@ using namespace std;
 #include "TText.h"
 #include "TPaveText.h"
 #include "TSystem.h"
+#include "TStatToolkit.h"
+#include "TCut.h"
 
-// Thread-stuff
-//#include "TThread.h"
-//#include "TMutex.h"
-//#include "TLockFile.h"
 
 
 ClassImp(AliTPCcalibTracks)
@@ -128,13 +147,7 @@ AliTPCcalibTracks::AliTPCcalibTracks():
   fHclusterPerPadrowRaw(0),
   fClusterCutHisto(0),
   fCalPadClusterPerPad(0),
-  fCalPadClusterPerPadRaw(0),
-  fFitterLinY1(0),   //!
-  fFitterLinZ1(0),   //! 
-  fFitterLinY2(0),   //! 
-  fFitterLinZ2(0),   //!
-  fFitterParY(0),    //! 
-  fFitterParZ(0)    //!  
+  fCalPadClusterPerPadRaw(0)
 { 
    // 
    // AliTPCcalibTracks default constructor
@@ -170,13 +183,7 @@ AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks):
   fHclusterPerPadrowRaw(0),
   fClusterCutHisto(0),
   fCalPadClusterPerPad(0),
-  fCalPadClusterPerPadRaw(0),
-  fFitterLinY1(0),   //!
-  fFitterLinZ1(0),   //! 
-  fFitterLinY2(0),   //! 
-  fFitterLinZ2(0),   //!
-  fFitterParY(0),    //! 
-  fFitterParZ(0)    //!  
+  fCalPadClusterPerPadRaw(0)
 {
    // 
    // AliTPCcalibTracks copy constructor
@@ -282,13 +289,7 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al
   fHclusterPerPadrowRaw(0),
   fClusterCutHisto(0),
   fCalPadClusterPerPad(0),
-  fCalPadClusterPerPadRaw(0),
-  fFitterLinY1(0),   //!
-  fFitterLinZ1(0),   //! 
-  fFitterLinY2(0),   //! 
-  fFitterLinZ2(0),   //!
-  fFitterParY(0),    //! 
-  fFitterParZ(0)    //!  
+  fCalPadClusterPerPadRaw(0)
  {
    // 
    // AliTPCcalibTracks constructor
@@ -323,7 +324,7 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al
    TProfile * prof1=0;
    TH1F     * his1 =0;
    fHclus = new TH1I("hclus","Number of clusters per track",160, 0, 160);     // valgrind 3
-   fRejectedTracksHisto    = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 101, 10);
+   fRejectedTracksHisto    = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 100, -1, 10);
    fHclusterPerPadrow      = new TH1I("fHclusterPerPadrow", " clusters per padRow, used for the resolution tree", 160, 0, 160);
    fHclusterPerPadrowRaw   = new TH1I("fHclusterPerPadrowRaw", " clusters per padRow, before cutting clusters", 160, 0, 160);
    fCalPadClusterPerPad    = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad");
@@ -420,18 +421,18 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al
       for (Int_t ipad = 0; ipad < 3; ipad++){
          Int_t   bin   = GetBin(iq, ipad);
          Float_t qmean = GetQ(bin);
-         char name[200];
-         sprintf(name,"ResolY Pad%d Qmiddle%f",ipad, qmean);
-         his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1);
+         char hname[200];
+         sprintf(hname,"ResolY Pad%d Qmiddle%f",ipad, qmean);
+         his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
          fArrayQDY->AddAt(his3D, bin);
-         sprintf(name,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
-         his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1);
+         sprintf(hname,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
+         his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
          fArrayQDZ->AddAt(his3D, bin);
-         sprintf(name,"RMSY Pad%d Qmiddle%f",ipad, qmean);
-         his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1);
+         sprintf(hname,"RMSY Pad%d Qmiddle%f",ipad, qmean);
+         his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
          fArrayQRMSY->AddAt(his3D, bin);
-         sprintf(name,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
-         his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1);
+         sprintf(hname,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
+         his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
          fArrayQRMSZ->AddAt(his3D, bin);
       }
    }
@@ -450,20 +451,6 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al
       }
    }
    
-   fFitterLinY1 = new TLinearFitter (2,"pol1");
-   fFitterLinZ1 = new TLinearFitter (2,"pol1");
-   fFitterLinY2 = new TLinearFitter (2,"pol1");
-   fFitterLinZ2 = new TLinearFitter (2,"pol1");  
-   fFitterParY  = new TLinearFitter (3,"pol2");
-   fFitterParZ  = new TLinearFitter (3,"pol2");
-   //
-   fFitterLinY1->StoreData(kFALSE);
-   fFitterLinZ1->StoreData(kFALSE);
-   fFitterLinY2->StoreData(kFALSE);
-   fFitterLinZ2->StoreData(kFALSE);
-   fFitterParY->StoreData(kFALSE);
-   fFitterParZ->StoreData(kFALSE);
-
 
    if (GetDebugLevel() > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl; 
    cout << "end of main constructor" << endl; // TO BE REMOVED
@@ -513,13 +500,7 @@ AliTPCcalibTracks::~AliTPCcalibTracks() {
       delete fArrayChargeVsDriftlength->At(i);
    }
    
-   delete fFitterLinY1;
-   delete fFitterLinZ1;
-   delete fFitterLinY2;
-   delete fFitterLinZ2;
-   delete fFitterParY;
-   delete fFitterParZ;
-   
+    
    delete fArrayQDY;
    delete fArrayQDZ;
    delete fArrayQRMSY;
@@ -540,16 +521,6 @@ AliTPCcalibTracks::~AliTPCcalibTracks() {
 }
    
   
-void AliTPCcalibTracks::AddInfo(TChain * chain, char* fileName){
-   // 
-   // Add the neccessary information for processing to the chain 
-   // (cluster parametrization)
-   // 
-   TFile clusterParamFile(fileName);
-   AliTPCClusterParam *clusterParam  =  (AliTPCClusterParam *) clusterParamFile.Get("Param");
-   chain->GetUserInfo()->AddLast((TObject*)clusterParam);
-   cout << "Clusterparametrization added to the chain." << endl;
-}
 
 void AliTPCcalibTracks::Process(AliTPCseed *track){
    // 
@@ -639,7 +610,7 @@ Int_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){
   //if (TMath::Abs(track->GetZ())<10.) return kFALSE;
   //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
   
-  if (GetDebugLevel() > 5) Info("AcceptTrack","Track has been accepted.");  
+  if (GetDebugLevel() > 20) Info("AcceptTrack","Track has been accepted.");  
   return 0;
 }
 
@@ -671,10 +642,25 @@ void  AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
    // and to avoid redundant data
    // 
 
+  static TLinearFitter fFitterLinY1(2,"pol1");   //
+  static TLinearFitter fFitterLinZ1(2,"pol1");   // 
+  static TLinearFitter fFitterLinY2(2,"pol1");   // 
+  static TLinearFitter fFitterLinZ2(2,"pol1");   //
+  static TLinearFitter fFitterParY(3,"pol2");    // 
+  static TLinearFitter fFitterParZ(3,"pol2");    //
+
+  fFitterLinY1.StoreData(kFALSE);
+  fFitterLinZ1.StoreData(kFALSE);
+  fFitterLinY2.StoreData(kFALSE);
+  fFitterLinZ2.StoreData(kFALSE);
+  fFitterParY.StoreData(kFALSE);
+  fFitterParZ.StoreData(kFALSE);
+
+
   if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****");
    const Int_t   kDelta    = 10;          // delta rows to fit
    const Float_t kMinRatio = 0.75;        // minimal ratio
-   const Float_t kCutChi2  = 6.;          // cut chi2 - left right  - kink removal
+   //   const Float_t kCutChi2  = 6.;          // cut chi2 - left right  - kink removal
    const Float_t kErrorFraction = 0.5;    // use only clusters with small interpolation error - for error param
    const Int_t   kFirstLargePad = 127;    // medium pads -> long pads
    const Float_t kLargePadSize  = 1.5;    // factor between medium and long pads' area
@@ -715,27 +701,27 @@ void  AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
       if (sector != sectorG){
          // track leaves sector before it crossed enough rows to fit / initialization
          nClusters = 0;
-         fFitterParY->ClearPoints();
-         fFitterParZ->ClearPoints();
+         fFitterParY.ClearPoints();
+         fFitterParZ.ClearPoints();
          sectorG = sector;
       }
       else {
          nClusters++;
          Double_t x = cluster0->GetX();
-         fFitterParY->AddPoint(&x, cluster0->GetY(), 1);
-         fFitterParZ->AddPoint(&x, cluster0->GetZ(), 1);
+         fFitterParY.AddPoint(&x, cluster0->GetY(), 1);
+         fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1);
          //
          if ( nClusters >= kDelta + 3 ){  
          // if more than 13 (kDelta+3) clusters were added to the fitters
          // fit the tracklet, increase trackletCounter
-         fFitterParY->Eval();
-         fFitterParZ->Eval();
+         fFitterParY.Eval();
+         fFitterParZ.Eval();
          nTrackletsAll++;
-         csigmaY += fFitterParY->GetChisquare() / (nClusters - 3.);
-         csigmaZ += fFitterParZ->GetChisquare() / (nClusters - 3.);
+         csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.);
+         csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.);
          nClusters = -1;
-         fFitterParY->ClearPoints();
-         fFitterParZ->ClearPoints();
+         fFitterParY.ClearPoints();
+         fFitterParZ.ClearPoints();
          }
       }
    }      // for (Int_t irow = 0; irow < 159; irow++)
@@ -757,12 +743,12 @@ void  AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
       Float_t xref = cluster0->GetX();
          
       // Make Fit
-      fFitterParY->ClearPoints();
-      fFitterParZ->ClearPoints();
-      fFitterLinY1->ClearPoints();
-      fFitterLinZ1->ClearPoints();
-      fFitterLinY2->ClearPoints();
-      fFitterLinZ2->ClearPoints();
+      fFitterParY.ClearPoints();
+      fFitterParZ.ClearPoints();
+      fFitterLinY1.ClearPoints();
+      fFitterLinZ1.ClearPoints();
+      fFitterLinY2.ClearPoints();
+      fFitterLinZ2.ClearPoints();
       
       // fit tracklet (clusters in given padrow +- kDelta padrows) 
       // with polynom of 2nd order and two polynoms of 1st order
@@ -785,89 +771,102 @@ void  AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
          nclFound++;
          if (idelta < 0){
          ncl0++;
-         fFitterLinY1->AddPoint(&x, currentCluster->GetY(), csigmaY);
-         fFitterLinZ1->AddPoint(&x, currentCluster->GetZ(), csigmaZ);
+         fFitterLinY1.AddPoint(&x, currentCluster->GetY(), csigmaY);
+         fFitterLinZ1.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
          }
          if (idelta > 0){
          ncl1++;
-         fFitterLinY2->AddPoint(&x, currentCluster->GetY(), csigmaY);
-         fFitterLinZ2->AddPoint(&x, currentCluster->GetZ(), csigmaZ);
+         fFitterLinY2.AddPoint(&x, currentCluster->GetY(), csigmaY);
+         fFitterLinZ2.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
          }
-         fFitterParY->AddPoint(&x, currentCluster->GetY(), csigmaY);  
-         fFitterParZ->AddPoint(&x, currentCluster->GetZ(), csigmaZ);  
+         fFitterParY.AddPoint(&x, currentCluster->GetY(), csigmaY);  
+         fFitterParZ.AddPoint(&x, currentCluster->GetZ(), csigmaZ);  
       }  // loop over neighbourhood for fitter filling 
+
+
       
       if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
       if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
       if (nclFound < kDelta * kMinRatio) continue;    // if not enough clusters (7.5) found in neighbourhood goto next padrow
-      fFitterParY->Eval();
-      fFitterParZ->Eval();
-      Double_t chi2 = (fFitterParY->GetChisquare() + fFitterParZ->GetChisquare()) / (2. * nclFound - 6.);
-      if (chi2 > kCutChi2) fRejectedTracksHisto->Fill(9);
-      if (chi2 > kCutChi2) fClusterCutHisto->Fill(2, irow);
-      if (chi2 > kCutChi2) continue;   // if chi^2 is too big goto next padrow
-      
+      fFitterParY.Eval();
+      fFitterParZ.Eval();
+      Double_t chi2 = (fFitterParY.GetChisquare() + fFitterParZ.GetChisquare()) / (2. * nclFound - 6.);
+      //if (chi2 > kCutChi2) fRejectedTracksHisto->Fill(9);
+      //if (chi2 > kCutChi2) fClusterCutHisto->Fill(2, irow);
+      //if (chi2 > kCutChi2) continue;   // if chi^2 is too big goto next padrow
+      TTreeSRedirector *cstream = GetDebugStreamer();
+      if (cstream){
+       (*cstream)<<"Cut9"<<
+         "chi2="<<chi2<<
+         "\n";
+      }
       // REMOVE KINK
       // only when there are enough clusters (4) in each direction
       if (ncl0 > 4){
-         fFitterLinY1->Eval();
-         fFitterLinZ1->Eval();
+         fFitterLinY1.Eval();
+         fFitterLinZ1.Eval();
       }
       if (ncl1 > 4){
-         fFitterLinY2->Eval();
-         fFitterLinZ2->Eval();
+         fFitterLinY2.Eval();
+         fFitterLinZ2.Eval();
       }
       
       if (ncl0 > 4 && ncl1 > 4){
-         fFitterLinY1->GetCovarianceMatrix(matrixY0);
-         fFitterLinY2->GetCovarianceMatrix(matrixY1);
-         fFitterLinZ1->GetCovarianceMatrix(matrixZ0);
-         fFitterLinZ2->GetCovarianceMatrix(matrixZ1);
-         fFitterLinY2->GetParameters(paramY1);
-         fFitterLinZ2->GetParameters(paramZ1);
-         fFitterLinY1->GetParameters(paramY0);
-         fFitterLinZ1->GetParameters(paramZ0);
+         fFitterLinY1.GetCovarianceMatrix(matrixY0);
+         fFitterLinY2.GetCovarianceMatrix(matrixY1);
+         fFitterLinZ1.GetCovarianceMatrix(matrixZ0);
+         fFitterLinZ2.GetCovarianceMatrix(matrixZ1);
+         fFitterLinY2.GetParameters(paramY1);
+         fFitterLinZ2.GetParameters(paramZ1);
+         fFitterLinY1.GetParameters(paramY0);
+         fFitterLinZ1.GetParameters(paramZ0);
          paramY0 -= paramY1;
          paramZ0 -= paramZ1;
          matrixY0 += matrixY1;
          matrixZ0 += matrixZ1;
-         Double_t chi2 = 0;
+         Double_t cchi2 = 0;
          
          TMatrixD difY(2, 1, paramY0.GetMatrixArray());
          TMatrixD difYT(1, 2, paramY0.GetMatrixArray());
          matrixY0.Invert();
          TMatrixD mulY(matrixY0, TMatrixD::kMult, difY);
          TMatrixD chi2Y(difYT, TMatrixD::kMult, mulY);
-         chi2 += chi2Y(0, 0);
+         cchi2 += chi2Y(0, 0);
          
          TMatrixD difZ(2, 1, paramZ0.GetMatrixArray());
          TMatrixD difZT(1, 2, paramZ0.GetMatrixArray());
          matrixZ0.Invert();
          TMatrixD mulZ(matrixZ0, TMatrixD::kMult, difZ);
          TMatrixD chi2Z(difZT, TMatrixD::kMult, mulZ);
-         chi2 += chi2Z(0, 0);      
+         cchi2 += chi2Z(0, 0);      
          
-         // REMOVE KINK
-         if (chi2 * 0.25 > kCutChi2) fRejectedTracksHisto->Fill(8);
-         if (chi2 * 0.25 > kCutChi2) fClusterCutHisto->Fill(3, irow);
-         if (chi2 * 0.25 > kCutChi2) continue;   // if chi2 is too big goto next padrow
+         // REMOVE KINK - TO be fixed - proper chi2 calculation for curved track to be implemented
+         //if (chi2 * 0.25 > kCutChi2) fRejectedTracksHisto->Fill(8);
+         //if (chi2 * 0.25 > kCutChi2) fClusterCutHisto->Fill(3, irow);
+         //if (chi2 * 0.25 > kCutChi2) continue;   // if chi2 is too big goto next padrow
          // fit tracklet with polynom of 2nd order and two polynoms of 1st order
          // take both polynoms of 1st order, calculate difference of their parameters
          // add covariance matrixes and calculate chi2 of this difference
          // if this chi2 is bigger than a given threshold, assume that the current cluster is
          // a kink an goto next padrow
+
+        if (cstream){
+          (*cstream)<<"Cut8"<<
+            "chi2="<<cchi2<<
+            "\n";
+        }       
       }
       
       // current padrow has no kink
       
       // get fit parameters from pol2 fit: 
       Double_t paramY[4], paramZ[4];
-      paramY[0] = fFitterParY->GetParameter(0);
-      paramY[1] = fFitterParY->GetParameter(1);
-      paramY[2] = fFitterParY->GetParameter(2);
-      paramZ[0] = fFitterParZ->GetParameter(0);
-      paramZ[1] = fFitterParZ->GetParameter(1);
-      paramZ[2] = fFitterParZ->GetParameter(2);    
+      paramY[0] = fFitterParY.GetParameter(0);
+      paramY[1] = fFitterParY.GetParameter(1);
+      paramY[2] = fFitterParY.GetParameter(2);
+      paramZ[0] = fFitterParZ.GetParameter(0);
+      paramZ[1] = fFitterParZ.GetParameter(1);
+      paramZ[2] = fFitterParZ.GetParameter(2);    
       
       Double_t tracky = paramY[0];
       Double_t trackz = paramZ[0];
@@ -921,8 +920,29 @@ void  AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
       
       // Fill resolution histograms
       Bool_t useForResol = kTRUE;
-      if (fFitterParY->GetParError(0) > kErrorFraction * csigmaY) useForResol = kFALSE;
+      if (fFitterParY.GetParError(0) > kErrorFraction * csigmaY) useForResol = kFALSE;
    
+      if (cstream){
+       Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
+       Float_t sy = cluster0->GetSigmaY2();
+       Float_t sz = cluster0->GetSigmaZ2();
+       (*cstream)<<"Resol0"<<
+         "run="<<fRun<<              //  run number
+         "event="<<fEvent<<          //  event number
+         "time="<<fTime<<            //  time stamp of event
+         "trigger="<<fTrigger<<      //  trigger
+         "mag="<<fMagF<<             //  magnetic field              
+         "padSize="<<padSize<<
+         "angley="<<angley<<
+         "anglez="<<anglez<<
+         "zdr="<<zdrift<<
+         "dy="<<deltay<<
+         "dz="<<deltaz<<
+         "sy="<<sy<<
+         "sz="<<sz<<
+         "\n";
+      }
+
       if (useForResol){
          fDeltaY->Fill(deltay);
          fDeltaZ->Fill(deltaz);
@@ -940,7 +960,7 @@ void  AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
       
       if (useForResol && nclFound > 2 * kMinRatio * kDelta 
          && irow % kDeltaWriteDebugStream == 0 && GetDebugLevel() > 4){
-       if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow);
+       if (GetDebugLevel() > 20) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow);
          FillResolutionHistoLocalDebugPart(track, cluster0, irow, angley, anglez, nclFound, kDelta);
       }  // if (useForResol && nclFound > 2 * kMinRatio * kDelta)
    
@@ -965,267 +985,276 @@ void AliTPCcalibTracks::FillResolutionHistoLocalDebugPart(AliTPCseed *track, Ali
       if (cluster0->GetRow() > 63) padSize = 2; // long pads
    }
       
-      static TLinearFitter fitY0(3, "pol2");
-      static TLinearFitter fitZ0(3, "pol2");
-      static TLinearFitter fitY2(5, "hyp4");
-      static TLinearFitter fitZ2(5, "hyp4");
-      static TLinearFitter fitY2Q(5, "hyp4");
-      static TLinearFitter fitZ2Q(5, "hyp4");
-      static TLinearFitter fitY2S(5, "hyp4");
-      static TLinearFitter fitZ2S(5, "hyp4");
-      fitY0.ClearPoints();
-      fitZ0.ClearPoints();
-      fitY2.ClearPoints();
-      fitZ2.ClearPoints();
-      fitY2Q.ClearPoints();
-      fitZ2Q.ClearPoints();
-      fitY2S.ClearPoints();
-      fitZ2S.ClearPoints();
-      
-      for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
-         // loop over irow +- kDelta rows (neighboured rows)
-         // 
-         // 
-       if (idelta == 0) continue;
-       if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
-       AliTPCclusterMI * cluster = track->GetClusterPointer(irow + idelta);
-       if (!cluster) continue;
-       if (cluster->GetType() < 0) continue;
-       if (cluster->GetDetector() != sector) continue;
-       Double_t x = cluster->GetX() - xref;
-       Double_t sigmaY0 = fClusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley) );
-       Double_t sigmaZ0 = fClusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez) );
-       //
-       Double_t sigmaYQ = fClusterParam->GetErrorQPar( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
-       Double_t sigmaZQ = fClusterParam->GetErrorQPar( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
-       Double_t sigmaYS = fClusterParam->GetErrorQParScaled( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), 
+   static TLinearFitter fitY0(3, "pol2");
+   static TLinearFitter fitZ0(3, "pol2");
+   static TLinearFitter fitY2(5, "hyp4");
+   static TLinearFitter fitZ2(5, "hyp4");
+   static TLinearFitter fitY2Q(5, "hyp4");
+   static TLinearFitter fitZ2Q(5, "hyp4");
+   static TLinearFitter fitY2S(5, "hyp4");
+   static TLinearFitter fitZ2S(5, "hyp4");
+   fitY0.ClearPoints();
+   fitZ0.ClearPoints();
+   fitY2.ClearPoints();
+   fitZ2.ClearPoints();
+   fitY2Q.ClearPoints();
+   fitZ2Q.ClearPoints();
+   fitY2S.ClearPoints();
+   fitZ2S.ClearPoints();
+   
+   for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
+     // loop over irow +- kDelta rows (neighboured rows)
+     // 
+     // 
+     if (idelta == 0) continue;
+     if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
+     AliTPCclusterMI * cluster = track->GetClusterPointer(irow + idelta);
+     if (!cluster) continue;
+     if (cluster->GetType() < 0) continue;
+     if (cluster->GetDetector() != sector) continue;
+     Double_t x = cluster->GetX() - xref;
+     Double_t sigmaY0 = fClusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley) );
+     Double_t sigmaZ0 = fClusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez) );
+     //
+     Double_t sigmaYQ = fClusterParam->GetErrorQPar( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
+     Double_t sigmaZQ = fClusterParam->GetErrorQPar( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
+     Double_t sigmaYS = fClusterParam->GetErrorQParScaled( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), 
                                                            TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
-       Double_t sigmaZS = fClusterParam->GetErrorQParScaled( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), 
+     Double_t sigmaZS = fClusterParam->GetErrorQParScaled( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), 
                                                            TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
-       Float_t rmsYFactor = fClusterParam->GetShapeFactor( 0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
-                                                          TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
-                                                          TMath::Sqrt(cluster0->GetSigmaY2()), 0 );
-       Float_t rmsZFactor = fClusterParam->GetShapeFactor(0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
-                                                          TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
-                                                          TMath::Sqrt(cluster0->GetSigmaZ2()),0 );
-       sigmaYS  = TMath::Sqrt(sigmaYS * sigmaYS + rmsYFactor * rmsYFactor / 12.);
-       sigmaZS  = TMath::Sqrt(sigmaZS * sigmaZS + rmsZFactor * rmsZFactor / 12. + rmsYFactor * rmsYFactor / 24.);
-       //
-       if (kDelta != 0){
-         fitY0.AddPoint(&x, cluster->GetY(), sigmaY0);
-         fitZ0.AddPoint(&x, cluster->GetZ(), sigmaZ0);
-       }
-       Double_t xxx[4];
-       xxx[0] = ( (idelta+irow) % 2 == 0 ) ? 1 : 0;
-       xxx[1] = x;
-       xxx[2] = ( (idelta+irow) % 2 == 0 ) ? x : 0;
-       xxx[3] = x * x; 
-       fitY2.AddPoint(xxx, cluster->GetY(), sigmaY0);
-       fitY2Q.AddPoint(xxx, cluster->GetY(), sigmaYQ);
-       fitY2S.AddPoint(xxx, cluster->GetY(), sigmaYS);
-       fitZ2.AddPoint(xxx, cluster->GetZ(), sigmaZ0);
-       fitZ2Q.AddPoint(xxx, cluster->GetZ(), sigmaZQ);
-       fitZ2S.AddPoint(xxx, cluster->GetZ(), sigmaZS);
-       //
-      }  // neigbouhood-loop
-      //
-      fitY0.Eval();
-      fitZ0.Eval();
-      fitY2.Eval();
-      fitZ2.Eval();
-      fitY2Q.Eval();
-      fitZ2Q.Eval();
-      fitY2S.Eval();
-      fitZ2S.Eval();
-      Float_t chi2Y0 = fitY0.GetChisquare() / (nclFound-3.);
-      Float_t chi2Z0 = fitZ0.GetChisquare() / (nclFound-3.);
-      Float_t chi2Y2 = fitY2.GetChisquare() / (nclFound-5.);
-      Float_t chi2Z2 = fitZ2.GetChisquare() / (nclFound-5.);
-      Float_t chi2Y2Q = fitY2Q.GetChisquare() / (nclFound-5.);
-      Float_t chi2Z2Q = fitZ2Q.GetChisquare() / (nclFound-5.);
-      Float_t chi2Y2S = fitY2S.GetChisquare() / (nclFound-5.);
-      Float_t chi2Z2S = fitZ2S.GetChisquare() / (nclFound-5.);
+     Float_t rmsYFactor = fClusterParam->GetShapeFactor( 0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
+                                                        TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
+                                                        TMath::Sqrt(cluster0->GetSigmaY2()), 0 );
+     Float_t rmsZFactor = fClusterParam->GetShapeFactor(0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
+                                                       TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
+                                                       TMath::Sqrt(cluster0->GetSigmaZ2()),0 );
+     sigmaYS  = TMath::Sqrt(sigmaYS * sigmaYS + rmsYFactor * rmsYFactor / 12.);
+     sigmaZS  = TMath::Sqrt(sigmaZS * sigmaZS + rmsZFactor * rmsZFactor / 12. + rmsYFactor * rmsYFactor / 24.);
+     //
+     if (kDelta != 0){
+       fitY0.AddPoint(&x, cluster->GetY(), sigmaY0);
+       fitZ0.AddPoint(&x, cluster->GetZ(), sigmaZ0);
+     }
+     Double_t xxx[4];
+     xxx[0] = ( (idelta+irow) % 2 == 0 ) ? 1 : 0;
+     xxx[1] = x;
+     xxx[2] = ( (idelta+irow) % 2 == 0 ) ? x : 0;
+     xxx[3] = x * x;   
+     fitY2.AddPoint(xxx, cluster->GetY(), sigmaY0);
+     fitY2Q.AddPoint(xxx, cluster->GetY(), sigmaYQ);
+     fitY2S.AddPoint(xxx, cluster->GetY(), sigmaYS);
+     fitZ2.AddPoint(xxx, cluster->GetZ(), sigmaZ0);
+     fitZ2Q.AddPoint(xxx, cluster->GetZ(), sigmaZQ);
+     fitZ2S.AddPoint(xxx, cluster->GetZ(), sigmaZS);
+     //
+   }  // neigbouhood-loop
       //
-      static  TVectorD    parY0(3);
-      static  TMatrixD    matY0(3, 3);
-      static  TVectorD    parZ0(3);
-      static  TMatrixD    matZ0(3, 3);
-      fitY0.GetParameters(parY0);
-      fitY0.GetCovarianceMatrix(matY0);
-      fitZ0.GetParameters(parZ0);
-      fitZ0.GetCovarianceMatrix(matZ0);
-      //
-      static  TVectorD    parY2(5);
-      static  TMatrixD    matY2(5,5);
-      static  TVectorD    parZ2(5);
-      static  TMatrixD    matZ2(5,5);
-      fitY2.GetParameters(parY2);
-      fitY2.GetCovarianceMatrix(matY2);
-      fitZ2.GetParameters(parZ2);
-      fitZ2.GetCovarianceMatrix(matZ2);
-      //
-      static  TVectorD    parY2Q(5);
-      static  TMatrixD    matY2Q(5,5);
-      static  TVectorD    parZ2Q(5);
-      static  TMatrixD    matZ2Q(5,5);
-      fitY2Q.GetParameters(parY2Q);
-      fitY2Q.GetCovarianceMatrix(matY2Q);
-      fitZ2Q.GetParameters(parZ2Q);
-      fitZ2Q.GetCovarianceMatrix(matZ2Q);
-      static  TVectorD    parY2S(5);
-      static  TMatrixD    matY2S(5,5);
-      static  TVectorD    parZ2S(5);
-      static  TMatrixD    matZ2S(5,5);
-      fitY2S.GetParameters(parY2S);
-      fitY2S.GetCovarianceMatrix(matY2S);
-      fitZ2S.GetParameters(parZ2S);
-      fitZ2S.GetCovarianceMatrix(matZ2S);
-      Float_t sigmaY0   = TMath::Sqrt(matY0(0,0));
-      Float_t sigmaZ0   = TMath::Sqrt(matZ0(0,0));
-      Float_t sigmaDY0  = TMath::Sqrt(matY0(1,1));
-      Float_t sigmaDZ0  = TMath::Sqrt(matZ0(1,1));
-      Float_t sigmaY2   = TMath::Sqrt(matY2(1,1));
-      Float_t sigmaZ2   = TMath::Sqrt(matZ2(1,1));
-      Float_t sigmaDY2  = TMath::Sqrt(matY2(3,3));
-      Float_t sigmaDZ2  = TMath::Sqrt(matZ2(3,3));
-      Float_t sigmaY2Q  = TMath::Sqrt(matY2Q(1,1));
-      Float_t sigmaZ2Q  = TMath::Sqrt(matZ2Q(1,1));
-      Float_t sigmaDY2Q = TMath::Sqrt(matY2Q(3,3));
-      Float_t sigmaDZ2Q = TMath::Sqrt(matZ2Q(3,3));
-      Float_t sigmaY2S  = TMath::Sqrt(matY2S(1,1));
-      Float_t sigmaZ2S  = TMath::Sqrt(matZ2S(1,1));
-      Float_t sigmaDY2S = TMath::Sqrt(matY2S(3,3));
-      Float_t sigmaDZ2S = TMath::Sqrt(matZ2S(3,3));
-      
-      // Error parameters
-      Float_t csigmaY0 = fClusterParam->GetError0Par(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(angley));
-      Float_t csigmaZ0 = fClusterParam->GetError0Par(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(anglez));
-      //
-      Float_t csigmaYQ = fClusterParam->GetErrorQPar(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                    TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
-      Float_t csigmaZQ = fClusterParam->GetErrorQPar(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                      TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
-      Float_t csigmaYS = fClusterParam->GetErrorQParScaled(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                    TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
-      Float_t csigmaZS = fClusterParam->GetErrorQParScaled(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                      TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
-
-
-      // RMS parameters
-      Float_t meanRMSY = 0;
-      Float_t meanRMSZ = 0;
-      Int_t   nclRMS = 0;
-      for (Int_t idelta = -2; idelta <= 2; idelta++){
-        // loop over neighbourhood
-       if (idelta+irow < 0 || idelta+irow > 159) continue;
-//     if (idelta+irow>159) continue;
-       AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta);
-       if (!cluster) continue;
-       meanRMSY += TMath::Sqrt(cluster->GetSigmaY2());
-       meanRMSZ += TMath::Sqrt(cluster->GetSigmaZ2());
-       nclRMS++;
-      }
-      meanRMSY /= nclRMS; 
-      meanRMSZ /= nclRMS; 
-
-      Float_t rmsY      = TMath::Sqrt(cluster0->GetSigmaY2());  
-      Float_t rmsZ      = TMath::Sqrt(cluster0->GetSigmaZ2());
-      Float_t rmsYT     = fClusterParam->GetRMSQ(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                               TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
-      Float_t rmsZT     = fClusterParam->GetRMSQ(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                               TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
-      Float_t rmsYT0    = fClusterParam->GetRMS0(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                TMath::Abs(angley));
-      Float_t rmsZT0    = fClusterParam->GetRMS0(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
+   fitY0.Eval();
+   fitZ0.Eval();
+   fitY2.Eval();
+   fitZ2.Eval();
+   fitY2Q.Eval();
+   fitZ2Q.Eval();
+   fitY2S.Eval();
+   fitZ2S.Eval();
+   Float_t chi2Y0 = fitY0.GetChisquare() / (nclFound-3.);
+   Float_t chi2Z0 = fitZ0.GetChisquare() / (nclFound-3.);
+   Float_t chi2Y2 = fitY2.GetChisquare() / (nclFound-5.);
+   Float_t chi2Z2 = fitZ2.GetChisquare() / (nclFound-5.);
+   Float_t chi2Y2Q = fitY2Q.GetChisquare() / (nclFound-5.);
+   Float_t chi2Z2Q = fitZ2Q.GetChisquare() / (nclFound-5.);
+   Float_t chi2Y2S = fitY2S.GetChisquare() / (nclFound-5.);
+   Float_t chi2Z2S = fitZ2S.GetChisquare() / (nclFound-5.);
+   //
+   static  TVectorD    parY0(3);
+   static  TMatrixD    matY0(3, 3);
+   static  TVectorD    parZ0(3);
+   static  TMatrixD    matZ0(3, 3);
+   fitY0.GetParameters(parY0);
+   fitY0.GetCovarianceMatrix(matY0);
+   fitZ0.GetParameters(parZ0);
+   fitZ0.GetCovarianceMatrix(matZ0);
+   //
+   static  TVectorD    parY2(5);
+   static  TMatrixD    matY2(5,5);
+   static  TVectorD    parZ2(5);
+   static  TMatrixD    matZ2(5,5);
+   fitY2.GetParameters(parY2);
+   fitY2.GetCovarianceMatrix(matY2);
+   fitZ2.GetParameters(parZ2);
+   fitZ2.GetCovarianceMatrix(matZ2);
+   //
+   static  TVectorD    parY2Q(5);
+   static  TMatrixD    matY2Q(5,5);
+   static  TVectorD    parZ2Q(5);
+   static  TMatrixD    matZ2Q(5,5);
+   fitY2Q.GetParameters(parY2Q);
+   fitY2Q.GetCovarianceMatrix(matY2Q);
+   fitZ2Q.GetParameters(parZ2Q);
+   fitZ2Q.GetCovarianceMatrix(matZ2Q);
+   static  TVectorD    parY2S(5);
+   static  TMatrixD    matY2S(5,5);
+   static  TVectorD    parZ2S(5);
+   static  TMatrixD    matZ2S(5,5);
+   fitY2S.GetParameters(parY2S);
+   fitY2S.GetCovarianceMatrix(matY2S);
+   fitZ2S.GetParameters(parZ2S);
+   fitZ2S.GetCovarianceMatrix(matZ2S);
+   Float_t sigmaY0   = TMath::Sqrt(matY0(0,0));
+   Float_t sigmaZ0   = TMath::Sqrt(matZ0(0,0));
+   Float_t sigmaDY0  = TMath::Sqrt(matY0(1,1));
+   Float_t sigmaDZ0  = TMath::Sqrt(matZ0(1,1));
+   Float_t sigmaY2   = TMath::Sqrt(matY2(1,1));
+   Float_t sigmaZ2   = TMath::Sqrt(matZ2(1,1));
+   Float_t sigmaDY2  = TMath::Sqrt(matY2(3,3));
+   Float_t sigmaDZ2  = TMath::Sqrt(matZ2(3,3));
+   Float_t sigmaY2Q  = TMath::Sqrt(matY2Q(1,1));
+   Float_t sigmaZ2Q  = TMath::Sqrt(matZ2Q(1,1));
+   Float_t sigmaDY2Q = TMath::Sqrt(matY2Q(3,3));
+   Float_t sigmaDZ2Q = TMath::Sqrt(matZ2Q(3,3));
+   Float_t sigmaY2S  = TMath::Sqrt(matY2S(1,1));
+   Float_t sigmaZ2S  = TMath::Sqrt(matZ2S(1,1));
+   Float_t sigmaDY2S = TMath::Sqrt(matY2S(3,3));
+   Float_t sigmaDZ2S = TMath::Sqrt(matZ2S(3,3));
+   
+   // Error parameters
+   Float_t csigmaY0 = fClusterParam->GetError0Par(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(angley));
+   Float_t csigmaZ0 = fClusterParam->GetError0Par(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(anglez));
+   //
+   Float_t csigmaYQ = fClusterParam->GetErrorQPar(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
+                                                 TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
+   Float_t csigmaZQ = fClusterParam->GetErrorQPar(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
+                                                 TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
+   Float_t csigmaYS = fClusterParam->GetErrorQParScaled(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
+                                                       TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
+   Float_t csigmaZS = fClusterParam->GetErrorQParScaled(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
+                                                       TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
+   
+   
+   // RMS parameters
+   Float_t meanRMSY = 0;
+   Float_t meanRMSZ = 0;
+   Int_t   nclRMS = 0;
+   for (Int_t idelta = -2; idelta <= 2; idelta++){
+     // loop over neighbourhood
+     if (idelta+irow < 0 || idelta+irow > 159) continue;
+     //        if (idelta+irow>159) continue;
+     AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta);
+     if (!cluster) continue;
+     meanRMSY += TMath::Sqrt(cluster->GetSigmaY2());
+     meanRMSZ += TMath::Sqrt(cluster->GetSigmaZ2());
+     nclRMS++;
+   }
+   meanRMSY /= nclRMS; 
+   meanRMSZ /= nclRMS; 
+   
+   Float_t rmsY      = TMath::Sqrt(cluster0->GetSigmaY2());  
+   Float_t rmsZ      = TMath::Sqrt(cluster0->GetSigmaZ2());
+   Float_t rmsYT     = fClusterParam->GetRMSQ(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
+                                             TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
+   Float_t rmsZT     = fClusterParam->GetRMSQ(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
+                                             TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
+   Float_t rmsYT0    = fClusterParam->GetRMS0(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
+                                             TMath::Abs(angley));
+   Float_t rmsZT0    = fClusterParam->GetRMS0(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
                                                 TMath::Abs(anglez));
-      Float_t rmsYSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                    TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
-      Float_t rmsZSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                    TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
-      Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                        TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
-                                                        rmsY,meanRMSY);
-      Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                        TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
-                                                        rmsZ,meanRMSZ);
-      //
-      // cluster debug
-      TTreeSRedirector *cstream = GetDebugStreamer();
-      if (cstream){
-       (*cstream)<<"ResolCl"<< // valgrind 3   40,000 bytes in 1 blocks are possibly 
-         "Sector="<<sector<<
-         "Cl.="<<cluster0<<
-         "CSigmaY0="<<csigmaY0<<   // cluster errorY
-         "CSigmaYQ="<<csigmaYQ<<   // cluster errorY - q scaled
-         "CSigmaYS="<<csigmaYS<<   // cluster errorY - q scaled
-         "CSigmaZ0="<<csigmaZ0<<   // 
-         "CSigmaZQ="<<csigmaZQ<<
-         "CSigmaZS="<<csigmaZS<<
-         "shapeYF="<<rmsYFactor<<
-         "shapeZF="<<rmsZFactor<<
-         "rmsY="<<rmsY<<
-         "rmsZ="<<rmsZ<<
-         "rmsYM="<<meanRMSY<<
-         "rmsZM="<<meanRMSZ<<
-         "rmsYT="<<rmsYT<<
-         "rmsZT="<<rmsZT<<
-         "rmsYT0="<<rmsYT0<<
-         "rmsZT0="<<rmsZT0<<
-         "rmsYS="<<rmsYSigma<<  
-         "rmsZS="<<rmsZSigma<<
-         "padSize="<<padSize<<
-         "Ncl="<<nclFound<<    
-         "PY0.="<<&parY0<<
-         "PZ0.="<<&parZ0<<
-         "SigmaY0="<<sigmaY0<< 
-         "SigmaZ0="<<sigmaZ0<< 
-         "angley="<<angley<<
-         "anglez="<<anglez<<
-         "\n";
-       
-       //       tracklet dubug
-       (*cstream)<<"ResolTr"<< 
-         "padSize="<<padSize<<
-         "IPad="<<padSize<<
-         "Sector="<<sector<<
-         "Ncl="<<nclFound<<    
-         "chi2Y0="<<chi2Y0<<
-         "chi2Z0="<<chi2Z0<<
-         "chi2Y2="<<chi2Y2<<
-         "chi2Z2="<<chi2Z2<<
-         "chi2Y2Q="<<chi2Y2Q<<
-         "chi2Z2Q="<<chi2Z2Q<<
-         "chi2Y2S="<<chi2Y2S<<
-         "chi2Z2S="<<chi2Z2S<<
-         "PY0.="<<&parY0<<
-         "PZ0.="<<&parZ0<<
-         "PY2.="<<&parY2<<
-         "PZ2.="<<&parZ2<<
-         "PY2Q.="<<&parY2Q<<
-         "PZ2Q.="<<&parZ2Q<<
-         "PY2S.="<<&parY2S<<
-         "PZ2S.="<<&parZ2S<<
-         "SigmaY0="<<sigmaY0<< 
-         "SigmaZ0="<<sigmaZ0<< 
-         "SigmaDY0="<<sigmaDY0<< 
-         "SigmaDZ0="<<sigmaDZ0<< 
-         "SigmaY2="<<sigmaY2<< 
-         "SigmaZ2="<<sigmaZ2<< 
-         "SigmaDY2="<<sigmaDY2<< 
-         "SigmaDZ2="<<sigmaDZ2<< 
-         "SigmaY2Q="<<sigmaY2Q<< 
-         "SigmaZ2Q="<<sigmaZ2Q<< 
-         "SigmaDY2Q="<<sigmaDY2Q<< 
-         "SigmaDZ2Q="<<sigmaDZ2Q<< 
-         "SigmaY2S="<<sigmaY2S<< 
-         "SigmaZ2S="<<sigmaZ2S<< 
-         "SigmaDY2S="<<sigmaDY2S<< 
-         "SigmaDZ2S="<<sigmaDZ2S<< 
-         "angley="<<angley<<
-         "anglez="<<anglez<<
-         "\n";
-      }
-
+   Float_t rmsYSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
+                                                 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
+   Float_t rmsZSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
+                                                 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
+   Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
+                                                     TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
+                                                     rmsY,meanRMSY);
+   Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
+                                                     TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
+                                                     rmsZ,meanRMSZ);
+   //
+   // cluster debug
+   TTreeSRedirector *cstream = GetDebugStreamer();
+   if (cstream){
+     (*cstream)<<"ResolCl"<<   // valgrind 3   40,000 bytes in 1 blocks are possibly 
+       "run="<<fRun<<              //  run number
+       "event="<<fEvent<<          //  event number
+       "time="<<fTime<<            //  time stamp of event
+       "trigger="<<fTrigger<<      //  trigger
+       "mag="<<fMagF<<             //  magnetic field        
+       "Sector="<<sector<<
+       "Cl.="<<cluster0<<
+       "CSigmaY0="<<csigmaY0<<   // cluster errorY
+       "CSigmaYQ="<<csigmaYQ<<   // cluster errorY - q scaled
+       "CSigmaYS="<<csigmaYS<<   // cluster errorY - q scaled
+       "CSigmaZ0="<<csigmaZ0<<   // 
+       "CSigmaZQ="<<csigmaZQ<<
+       "CSigmaZS="<<csigmaZS<<
+       "shapeYF="<<rmsYFactor<<
+       "shapeZF="<<rmsZFactor<<
+       "rmsY="<<rmsY<<
+       "rmsZ="<<rmsZ<<
+       "rmsYM="<<meanRMSY<<
+       "rmsZM="<<meanRMSZ<<
+       "rmsYT="<<rmsYT<<
+       "rmsZT="<<rmsZT<<
+       "rmsYT0="<<rmsYT0<<
+       "rmsZT0="<<rmsZT0<<
+       "rmsYS="<<rmsYSigma<<  
+       "rmsZS="<<rmsZSigma<<
+       "padSize="<<padSize<<
+       "Ncl="<<nclFound<<      
+       "PY0.="<<&parY0<<
+       "PZ0.="<<&parZ0<<
+       "SigmaY0="<<sigmaY0<< 
+       "SigmaZ0="<<sigmaZ0<< 
+       "angley="<<angley<<
+       "anglez="<<anglez<<
+       "\n";
+     
+     //       tracklet dubug
+     (*cstream)<<"ResolTr"<<   
+       "run="<<fRun<<              //  run number
+       "event="<<fEvent<<          //  event number
+       "time="<<fTime<<            //  time stamp of event
+       "trigger="<<fTrigger<<      //  trigger
+       "mag="<<fMagF<<             //  magnetic field        
+       "padSize="<<padSize<<
+       "IPad="<<padSize<<
+       "Sector="<<sector<<
+       "Ncl="<<nclFound<<      
+       "chi2Y0="<<chi2Y0<<
+       "chi2Z0="<<chi2Z0<<
+       "chi2Y2="<<chi2Y2<<
+       "chi2Z2="<<chi2Z2<<
+       "chi2Y2Q="<<chi2Y2Q<<
+       "chi2Z2Q="<<chi2Z2Q<<
+       "chi2Y2S="<<chi2Y2S<<
+       "chi2Z2S="<<chi2Z2S<<
+       "PY0.="<<&parY0<<
+       "PZ0.="<<&parZ0<<
+       "PY2.="<<&parY2<<
+       "PZ2.="<<&parZ2<<
+       "PY2Q.="<<&parY2Q<<
+       "PZ2Q.="<<&parZ2Q<<
+       "PY2S.="<<&parY2S<<
+       "PZ2S.="<<&parZ2S<<
+       "SigmaY0="<<sigmaY0<< 
+       "SigmaZ0="<<sigmaZ0<< 
+       "SigmaDY0="<<sigmaDY0<< 
+       "SigmaDZ0="<<sigmaDZ0<< 
+       "SigmaY2="<<sigmaY2<< 
+       "SigmaZ2="<<sigmaZ2<< 
+       "SigmaDY2="<<sigmaDY2<< 
+       "SigmaDZ2="<<sigmaDZ2<< 
+       "SigmaY2Q="<<sigmaY2Q<< 
+       "SigmaZ2Q="<<sigmaZ2Q<< 
+       "SigmaDY2Q="<<sigmaDY2Q<< 
+       "SigmaDZ2Q="<<sigmaDZ2Q<< 
+       "SigmaY2S="<<sigmaY2S<< 
+       "SigmaZ2S="<<sigmaZ2S<< 
+       "SigmaDY2S="<<sigmaDY2S<< 
+       "SigmaDZ2S="<<sigmaDZ2S<< 
+       "angley="<<angley<<
+       "anglez="<<anglez<<
+       "\n";
+   }  
 }
 
 
@@ -1333,7 +1362,7 @@ void AliTPCcalibTracks::Draw(Option_t* opt){
 }
 
 
-void AliTPCcalibTracks::MakeReport(Int_t stat, char* pathName){ 
+void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){ 
    // 
    // all functions are called, that produce pictures
    // the histograms are written to the directory 'pathName'
@@ -1352,7 +1381,7 @@ void AliTPCcalibTracks::MakeReport(Int_t stat, char* pathName){
 }
    
 
-void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, char* pathName){ 
+void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, const char* pathName){ 
    // 
    // creates several plots:
    // fArrayAmp.ps, fArrayAmpRow.ps and DeltaYZ.ps
@@ -1422,7 +1451,7 @@ void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, char* pathName){
 }
 
 
-void AliTPCcalibTracks::MakeDeltaPlots(char* pathName){
+void AliTPCcalibTracks::MakeDeltaPlots(const char* pathName){
    // 
    // creates several plots:
    // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
@@ -1459,7 +1488,7 @@ void AliTPCcalibTracks::MakeDeltaPlots(char* pathName){
 }
 
 
-void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(char* pathName){
+void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(const char* pathName){
    // 
    // creates charge vs. driftlength plots, one TProfile for each ROC
    // is not correct like this, should be one TProfile for each sector and padsize
@@ -1501,7 +1530,7 @@ void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(char* pathName){
 }   
 
 
-void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(char* pathName){
+void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(const char* pathName){
    // 
    // creates charge vs. driftlength plots, one TProfile for each ROC
    // under development....
@@ -1563,7 +1592,7 @@ void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(char* pathName){
 
 
 
-void AliTPCcalibTracks::FitResolutionNew(char* pathName){
+void AliTPCcalibTracks::FitResolutionNew(const char* pathName){
    // 
    // calculates different resulution fits in Y and Z direction
    // the histograms are written to 'ResolutionYZ.ps'
@@ -1686,7 +1715,7 @@ void AliTPCcalibTracks::FitResolutionNew(char* pathName){
 }
 
 
-void AliTPCcalibTracks::FitRMSNew(char* pathName){
+void AliTPCcalibTracks::FitRMSNew(const char* pathName){
    // 
    // calculates different resulution-rms fits in Y and Z direction
    // the histograms are written to 'RMS_YZ.ps'
@@ -1809,7 +1838,7 @@ void AliTPCcalibTracks::FitRMSNew(char* pathName){
 }
 
 
-void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, char* pathName){
+void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){
    //
    // Make tree with resolution parameters
    // the result is written to 'resol.root' in directory 'pathname'
@@ -2162,7 +2191,7 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
    AliTPCcalibTracks *calibTracks = 0;
    if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;    
    Int_t counter = 0;
-   while ( (calibTracks = (AliTPCcalibTracks*)listIterator->Next()) ){
+   while ( (calibTracks = dynamic_cast<AliTPCcalibTracks*> (listIterator->Next())) ){
       // loop over all entries in the collectionList and get dataMembers into lists
       if (!calibTracks) continue;
       
@@ -2185,8 +2214,10 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
       clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
       hclusterPerPadrowList->Add(calibTracks->GetfHclusterPerPadrow());
       hclusterPerPadrowRawList->Add(calibTracks->GetfHclusterPerPadrowRaw());
-      fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
-      fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
+      //
+      if (fCalPadClusterPerPad && calibTracks->GetfCalPadClusterPerPad())
+       fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());      
+      //      fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
       counter++;
       if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
    }
@@ -2226,9 +2257,9 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
    if (GetDebugLevel() > 0) cout << "merging fArrayAmps..." << endl;
    // merge fArrayAmps
    for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++ ) {  // loop over data member, i<72
-      TIterator *objListIterator = arrayAmpList->MakeIterator();
+      TIterator *cobjListIterator = arrayAmpList->MakeIterator();
       histList = new TList;
-      while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
+      while (( objarray =  (TObjArray*)cobjListIterator->Next() )) { 
          // loop over arrayAmpList, get TObjArray, get object at position i, cast it into TH1F
          if (!objarray) continue;
          hist = (TH1F*)(objarray->At(i));
@@ -2236,7 +2267,7 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
       }
       ((TH1F*)(fArrayAmp->At(i)))->Merge(histList);
       delete histList;
-      delete objListIterator;
+      delete cobjListIterator;
    }
    
    if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
@@ -2489,5 +2520,124 @@ AliTPCcalibTracks* AliTPCcalibTracks::TestMerge(AliTPCcalibTracks *ct, AliTPCClu
 }
 
 
+void  AliTPCcalibTracks::MakeQPosNormAll(TTree * chainres, AliTPCClusterParam * param, Int_t maxPoints){
+  //
+  // Make position corrections
+  // for the moment Only using debug streamer 
+  // chainres  - debug tree
+  // param     - parameters to be updated
+  // maxPoints - maximal number of points using for fit
+  // verbose   - print info flag
+  //
+  // Current implementation - only using debug streamers
+  // 
+  
+  /*    
+    //Defaults
+    Int_t maxPoints=100000;
+  */
+  /*
+    Usage: 
+    //0. Load libraries
+    gSystem->Load("libANALYSIS");
+    gSystem->Load("libSTAT");
+    gSystem->Load("libTPCcalib");
+    
+
+    //1. Load Parameters to be modified:
+    //e.g:
+    
+    AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
+    AliCDBManager::Instance()->SetRun(0) 
+    AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
+
+    //2. Load chain from debug streamers
+    //
+    //e.g
+    gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
+    gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
+    AliXRDPROOFtoolkit tool;  
+    TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
+    chainres->Lookup();
+    //3. Do fits and store results
+    // 
+    AliTPCcalibTracks::MakeQPosNormAll(chainres,param,200000,0) ;
+    TFile f("paramout.root","recreate");
+    param->Write("clusterParam");
+    f.Close();
+    //4. Verification
+    TFile f2("paramout.root");
+    AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
+    param2->SetInstance(param2);
+    chainres->Draw("fitZ0:AliTPCClusterParam::SPosCorrection(1,0,Cl.fPad,Cl.fTimeBin,Cl.fZ,Cl.fSigmaY2,Cl.fSigmaZ2,Cl.fMax)","Cl.fDetector<36","",10000); // should be line 
+    
+   */
+
+
+  TStatToolkit toolkit;
+  Double_t chi2;
+  TVectorD fitParamY0;
+  TVectorD fitParamY1;
+  TVectorD fitParamZ0;
+  TVectorD fitParamZ1;
+  TMatrixD covMatrix;
+  Int_t npoints;
+  
+  chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
+  chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
+  chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
+  chainres->SetAlias("st","(sin(dt)-dt)");
+  //
+  chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
+  //
+  //
+  //
+  TCut cutA("1");
+  TString fstringY="";  
+  //
+  fstringY+="(dp)++";            //1
+  fstringY+="(dp)*di++";         //2
+  fstringY+="(sp)++";            //3
+  fstringY+="(sp)*di++";         //4
+  TString fstringZ="";  
+  fstringZ+="(dt)++";            //1
+  fstringZ+="(dt)*di++";         //2
+  fstringZ+="(st)++";            //3
+  fstringZ+="(st)*di++";         //4
+  //
+  // Z corrections
+  //
+  TString *strZ0 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
+  printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
+  param->fPosZcor[0] = (TVectorD*) fitParamZ0.Clone();
+  //
+  TString *strZ1 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
+  //
+  printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
+  param->fPosZcor[1] = (TVectorD*) fitParamZ1.Clone();
+  param->fPosZcor[2] = (TVectorD*) fitParamZ1.Clone();
+  //
+  // Y corrections
+  //   
+  TString *strY0 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
+  printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
+  param->fPosYcor[0] = (TVectorD*) fitParamY0.Clone();
+  
+
+  TString *strY1 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
+  //
+  printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
+  param->fPosYcor[1] = (TVectorD*) fitParamY1.Clone();
+  param->fPosYcor[2] = (TVectorD*) fitParamY1.Clone();
+  //
+  //
+  //
+  chainres->SetAlias("fitZ0",strZ0->Data());
+  chainres->SetAlias("fitZ1",strZ1->Data());
+  chainres->SetAlias("fitY0",strY0->Data());
+  chainres->SetAlias("fitY1",strY1->Data());
+  //  chainres->Draw("Cl.fZ-PZ0.fElements[0]","CSigmaY0<0.7&&CSigmaZ0<0.7"+cutA,"",10000);   
+}
+