TPC/AliTPCcalibTracks.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 //     Class for calibration of the cluster properties:                       //
  20 //     Cluster position resolution (RMS)  and short term distortions (within pad or within time bin)
  21 //
  22 //  Algorithm:
  23 //    1. Creation of the residual/properties  histograms
  24 //    2. Filling of the histograms.
  25 //       2.a Traklet fitting
  26 //       2.b Resudual filling
  27 //    3. Postprocessing: Creation of the tree with the mean values and RMS at different bins
  28 //    4.               : Paramaterization fitting.
  29 //                       In old AliRoot the RMS paramterization was used to characterize cluster resolution
  30 //                       Mean value /bias was neglected
  31 //    5. To be implemented
  32 //          a.) lookup table for the distortion parmaterization - THn
  33 //          b.) Usage in the reconstruction
  34 //
  35 //
  36 //   1. Creation of the histograms:   MakeHistos()
  37 //
  38 //         6 n dimensional histograms  are filled during the calibration:
  39 //         cluster  performance bins
  40 //         0 - variable of interest
  41 //         1 - pad type   - 0- IROC Short 1- OCROC medium 2 - OROC long pads
  42 //         2 - drift length - drift length -0-1
  43 //         3 - Qmax         - Qmax  - 2- 400
  44 //         4 - cog          - COG position - 0-1
  45 //         5 - tan(phi)     - local  angle
  46 //        Histograms:
  47 //        THnSparse  *fHisDeltaY;    // THnSparse - delta Y between the cluster and tracklet interpolation (+-X(5?) rows)
  48 //        THnSparse  *fHisDeltaZ;    // THnSparse - delta Z
  49 //        THnSparse  *fHisRMSY;      // THnSparse - rms Y
  50 //        THnSparse  *fHisRMSZ;      // THnSparse - rms Z
  51 //        THnSparse  *fHisQmax;      // THnSparse - qmax
  52 //        THnSparse  *fHisQtot;      // THnSparse - qtot
  53 //
  54 //
  55 //
  56
  57
  58                      //
  59 //     The parameter 'clusterParam', a AliTPCClusterParam object             //
  60 //      (needed for TPC cluster error and shape parameterization)            //
  61 //
  62 //     Normally you get this object out of the file 'TPCClusterParam.root'   //
  63 //     In the parameter 'cuts' the cuts are specified, that decide           //
  64 //     weather a track will be accepted for calibration or not.              //
  65 //                                                                           //
  66 //
  67 //     The data flow:
  68 //
  69 /*
  70    Raw Data -> Local Reconstruction -> Tracking ->     Calibration -> RefData (component itself)
  71                Offline/HLT             Offline/HLT                    OCDB entries (AliTPCClusterParam)
  72 */
  73
  74 /*
  75   Example usage - creation of the residual trees:
  76   TFile f("CalibObjects.root");
  77   AliTPCcalibTracks *calibTracks = ( AliTPCcalibTracks *)TPCCalib->FindObject("calibTracks");
  78   TH2 * his2 = calibTracks->fHisDeltaY->Projection(0,4);
  79   his2->SetName("his2");
  80   his2->FitSlicesY();
  81
  82
  83   TTreeSRedirector *pcstream = new TTreeSRedirector("clusterLookup.root");
  84   AliTPCcalibTracks::MakeSummaryTree(calibTracks->fHisDeltaY, pcstream, 0);
  85   AliTPCcalibTracks::MakeSummaryTree(calibTracks->fHisDeltaZ, pcstream, 1);
  86   delete pcstream;
  87   TFile fl("clusterLookup.root");
  88   TCut cutNI="cogA==0&&angleA==0&&qmaxA==0&&zA==0&&ipadA==0"
  89
  90 */
  91
  92                                                                //
  93 ///////////////////////////////////////////////////////////////////////////////
  94
  95 //
  96 // ROOT includes
  97 //
  98 #include <iostream>
  99 #include <fstream>
 100 using namespace std;
 101
 102 #include <TROOT.h>
 103 #include <TChain.h>
 104 #include <TFile.h>
 105 #include <TH3F.h>
 106 #include <TProfile.h>
 107
 108 //
 109 //#include <TPDGCode.h>
 110 #include <TStyle.h>
 111 #include "TLinearFitter.h"
 112 //#include "TMatrixD.h"
 113 #include "TTreeStream.h"
 114 #include "TF1.h"
 115 #include <TCanvas.h>
 116 #include <TGraph2DErrors.h>
 117 #include "TPostScript.h"
 118 #include "TCint.h"
 119
 120 #include <TH2D.h>
 121 #include <TF2.h>
 122 #include <TSystem.h>
 123 #include <TCollection.h>
 124 #include <iostream>
 125 #include <TLinearFitter.h>
 126 #include <TString.h>
 127
 128 //
 129 // AliROOT includes
 130 //
 131 #include "AliMagF.h"
 132 #include "AliTracker.h"
 133 #include "AliESD.h"
 134 #include "AliESDtrack.h"
 135 #include "AliESDfriend.h"
 136 #include "AliESDfriendTrack.h"
 137 #include "AliTPCseed.h"
 138 #include "AliTPCclusterMI.h"
 139 #include "AliTPCROC.h"
 140
 141 #include "AliTPCParamSR.h"
 142 #include "AliTrackPointArray.h"
 143 #include "AliTPCcalibTracks.h"
 144 #include "AliTPCClusterParam.h"
 145 #include "AliTPCcalibTracksCuts.h"
 146 #include "AliTPCCalPad.h"
 147 #include "AliTPCCalROC.h"
 148 #include "TText.h"
 149 #include "TPaveText.h"
 150 #include "TSystem.h"
 151 #include "TStatToolkit.h"
 152 #include "TCut.h"
 153 #include "THnSparse.h"
 154 #include "AliRieman.h"
 155
 156
 157
 158 ClassImp(AliTPCcalibTracks)
 159
 160
 161 AliTPCcalibTracks::AliTPCcalibTracks():
 162   AliTPCcalibBase(),
 163   fClusterParam(0),
 164   fROC(0),
 165   fHisDeltaY(0),    // THnSparse - delta Y
 166   fHisDeltaZ(0),    // THnSparse - delta Z
 167   fHisRMSY(0),      // THnSparse - rms Y
 168   fHisRMSZ(0),      // THnSparse - rms Z
 169   fHisQmax(0),      // THnSparse - qmax
 170   fHisQtot(0),      // THnSparse - qtot
 171   fArrayQDY(0),
 172   fArrayQDZ(0),
 173   fArrayQRMSY(0),
 174   fArrayQRMSZ(0),
 175   fResolY(0),
 176   fResolZ(0),
 177   fRMSY(0),
 178   fRMSZ(0),
 179   fCuts(0),
 180   fRejectedTracksHisto(0),
 181   fClusterCutHisto(0),
 182   fCalPadClusterPerPad(0),
 183   fCalPadClusterPerPadRaw(0)
 184 {
 185    //
 186    // AliTPCcalibTracks default constructor
 187    //
 188   if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' default constructor called" << endl;
 189 }
 190
 191
 192
 193 AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks):
 194   AliTPCcalibBase(calibTracks),
 195   fClusterParam(0),
 196   fROC(0),
 197   fHisDeltaY(0),    // THnSparse - delta Y
 198   fHisDeltaZ(0),    // THnSparse - delta Z
 199   fHisRMSY(0),      // THnSparse - rms Y
 200   fHisRMSZ(0),      // THnSparse - rms Z
 201   fHisQmax(0),      // THnSparse - qmax
 202   fHisQtot(0),      // THnSparse - qtot
 203   fArrayQDY(0),
 204   fArrayQDZ(0),
 205   fArrayQRMSY(0),
 206   fArrayQRMSZ(0),
 207   fResolY(0),
 208   fResolZ(0),
 209   fRMSY(0),
 210   fRMSZ(0),
 211   fCuts(0),
 212   fRejectedTracksHisto(0),
 213   fClusterCutHisto(0),
 214   fCalPadClusterPerPad(0),
 215   fCalPadClusterPerPadRaw(0)
 216 {
 217    //
 218    // AliTPCcalibTracks copy constructor
 219    //
 220   if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' copy constructor ***** " << endl;
 221
 222    Bool_t dirStatus = TH1::AddDirectoryStatus();
 223    TH1::AddDirectory(kFALSE);
 224
 225    Int_t length = -1;
 226
 227    (calibTracks.fArrayQDY) ? length = calibTracks.fArrayQDY->GetEntriesFast() : length = -1;
 228    fArrayQDY= new TObjArray(length);
 229    fArrayQDZ= new TObjArray(length);
 230    fArrayQRMSY= new TObjArray(length);
 231    fArrayQRMSZ= new TObjArray(length);
 232    for (Int_t i = 0; i < length; i++) {
 233       fArrayQDY->AddAt( ((TH1F*)calibTracks.fArrayQDY->At(i)->Clone()), i);
 234       fArrayQDZ->AddAt( ((TH1F*)calibTracks.fArrayQDZ->At(i)->Clone()), i);
 235       fArrayQRMSY->AddAt( ((TH1F*)calibTracks.fArrayQRMSY->At(i)->Clone()), i);
 236       fArrayQRMSZ->AddAt( ((TH1F*)calibTracks.fArrayQRMSZ->At(i)->Clone()), i);
 237    }
 238
 239    (calibTracks.fResolY) ? length = calibTracks.fResolY->GetEntriesFast() : length = -1;
 240    fResolY = new TObjArray(length);
 241    fResolZ = new TObjArray(length);
 242    fRMSY = new TObjArray(length);
 243    fRMSZ = new TObjArray(length);
 244    for (Int_t i = 0; i < length; i++) {
 245       fResolY->AddAt( ((TH1F*)calibTracks.fResolY->At(i)->Clone()), i);
 246       fResolZ->AddAt( ((TH1F*)calibTracks.fResolZ->At(i)->Clone()), i);
 247       fRMSY->AddAt( ((TH1F*)calibTracks.fRMSY->At(i)->Clone()), i);
 248       fRMSZ->AddAt( ((TH1F*)calibTracks.fRMSZ->At(i)->Clone()), i);
 249    }
 250
 251
 252    fClusterCutHisto = (TH2I*)calibTracks.fClusterCutHisto->Clone();
 253    fRejectedTracksHisto    = (TH1I*)calibTracks.fRejectedTracksHisto->Clone();
 254    fCalPadClusterPerPad    = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPad->Clone();
 255    fCalPadClusterPerPadRaw = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPadRaw->Clone();
 256
 257    fCuts = new AliTPCcalibTracksCuts(calibTracks.fCuts->GetMinClusters(), calibTracks.fCuts->GetMinRatio(),
 258       calibTracks.fCuts->GetMax1pt(), calibTracks.fCuts->GetEdgeYXCutNoise(), calibTracks.fCuts->GetEdgeThetaCutNoise());
 259    SetNameTitle(calibTracks.GetName(), calibTracks.GetTitle());
 260    TH1::AddDirectory(dirStatus); // set status back to original status
 261 //    cout << "+++++ end of copy constructor +++++" << endl;   // TO BE REMOVED
 262 }
 263
 264
 265 AliTPCcalibTracks & AliTPCcalibTracks::operator=(const AliTPCcalibTracks& calibTracks){
 266   //
 267   // assgnment operator
 268   //
 269   if (this != &calibTracks) {
 270     new (this) AliTPCcalibTracks(calibTracks);
 271   }
 272   return *this;
 273
 274 }
 275
 276
 277 AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam,  AliTPCcalibTracksCuts* cuts, Int_t logLevel) :
 278   AliTPCcalibBase(),
 279   fClusterParam(0),
 280   fROC(0),
 281   fHisDeltaY(0),    // THnSparse - delta Y
 282   fHisDeltaZ(0),    // THnSparse - delta Z
 283   fHisRMSY(0),      // THnSparse - rms Y
 284   fHisRMSZ(0),      // THnSparse - rms Z
 285   fHisQmax(0),      // THnSparse - qmax
 286   fHisQtot(0),      // THnSparse - qtot
 287   fArrayQDY(0),
 288   fArrayQDZ(0),
 289   fArrayQRMSY(0),
 290   fArrayQRMSZ(0),
 291   fResolY(0),
 292   fResolZ(0),
 293   fRMSY(0),
 294   fRMSZ(0),
 295   fCuts(0),
 296   fRejectedTracksHisto(0),
 297   fClusterCutHisto(0),
 298   fCalPadClusterPerPad(0),
 299   fCalPadClusterPerPadRaw(0)
 300  {
 301    //
 302    // AliTPCcalibTracks constructor
 303    // specify 'name' and 'title' of your object
 304    // specify 'clusterParam', (needed for TPC cluster error and shape parameterization)
 305    // In the parameter 'cuts' the cuts are specified, that decide
 306    // weather a track will be accepted for calibration or not.
 307    //
 308    // fDebugLevel - debug output: -1: silence, 0: default, 1: things like constructor called, 5: write fDebugStreamer, 6: waste your screen
 309    //
 310    // All histograms are instatiated in this constructor.
 311    //
 312    this->SetName(name);
 313    this->SetTitle(title);
 314
 315    if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' main constructor ***** " << endl;
 316    G__SetCatchException(0);
 317
 318    fClusterParam = clusterParam;
 319    if (fClusterParam){
 320      fClusterParam->SetInstance(fClusterParam);
 321    }
 322    else {
 323      Error("AliTPCcalibTracks","No cluster parametrization found! A valid clusterParam object is needed in the constructor. (To be found in 'TPCClusterParam.root'.)");
 324    }
 325    fCuts = cuts;
 326    SetDebugLevel(logLevel);
 327    MakeHistos();
 328
 329    TH1::AddDirectory(kFALSE);
 330
 331    fRejectedTracksHisto    = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 100, -1, 10);
 332    fCalPadClusterPerPad    = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad");
 333    fCalPadClusterPerPadRaw = new AliTPCCalPad("fCalPadClusterPerPadRaw", "clusters per pad, before cutting clusters");
 334    fClusterCutHisto = new TH2I("fClusterCutHisto", "Cutted cluster over padRow; Cut Criterium; PadRow", 5,1,5, 160,0,159);
 335
 336
 337    TH1::AddDirectory(kFALSE);
 338
 339
 340    fResolY = new TObjArray(3);
 341    fResolZ = new TObjArray(3);
 342    fRMSY   = new TObjArray(3);
 343    fRMSZ   = new TObjArray(3);
 344    TH3F * his3D;
 345    //
 346    his3D = new TH3F("Resol Y0","Resol Y0", 5,20,250, 4, 0,1., 50, -1,1);
 347    fResolY->AddAt(his3D,0);
 348    his3D = new TH3F("Resol Y1","Resol Y1", 5,20,250, 4, 0,1., 50, -1,1);
 349    fResolY->AddAt(his3D,1);
 350    his3D = new TH3F("Resol Y2","Resol Y2", 5,20,250, 4, 0,0.8, 50, -1,1);
 351    fResolY->AddAt(his3D,2);
 352    //
 353    his3D = new TH3F("Resol Z0","Resol Z0", 5,20,250, 4, 0,1, 50, -1,1);
 354    fResolZ->AddAt(his3D,0);
 355    his3D = new TH3F("Resol Z1","Resol Z1", 5,20,250, 4, 0,1, 50, -1,1);
 356    fResolZ->AddAt(his3D,1);
 357    his3D = new TH3F("Resol Z2","Resol Z2", 5,20,250, 4, 0,1, 50, -1,1);
 358    fResolZ->AddAt(his3D,2);
 359    //
 360    his3D = new TH3F("RMS Y0","RMS Y0", 5,20,250, 4, 0,1., 50, 0,0.8);
 361    fRMSY->AddAt(his3D,0);
 362    his3D = new TH3F("RMS Y1","RMS Y1", 5,20,250, 4, 0,1., 50, 0,0.8);
 363    fRMSY->AddAt(his3D,1);
 364    his3D = new TH3F("RMS Y2","RMS Y2", 5,20,250, 4, 0,0.8, 50, 0,0.8);
 365    fRMSY->AddAt(his3D,2);
 366    //
 367    his3D = new TH3F("RMS Z0","RMS Z0", 5,20,250, 4, 0,1, 50, 0,0.8);
 368    fRMSZ->AddAt(his3D,0);
 369    his3D = new TH3F("RMS Z1","RMS Z1", 5,20,250, 4, 0,1, 50, 0,0.8);
 370    fRMSZ->AddAt(his3D,1);
 371    his3D = new TH3F("RMS Z2","RMS Z2", 5,20,250, 4, 0,1, 50, 0,0.8);
 372    fRMSZ->AddAt(his3D,2);
 373    //
 374
 375    TH1::AddDirectory(kFALSE);
 376
 377    fArrayQDY = new TObjArray(300);
 378    fArrayQDZ = new TObjArray(300);
 379    fArrayQRMSY = new TObjArray(300);
 380    fArrayQRMSZ = new TObjArray(300);
 381    for (Int_t iq = 0; iq <= 10; iq++){
 382       for (Int_t ipad = 0; ipad < 3; ipad++){
 383          Int_t   bin   = GetBin(iq, ipad);
 384          Float_t qmean = GetQ(bin);
 385          char hname[200];
 386          snprintf(hname,100,"ResolY Pad%d Qmiddle%f",ipad, qmean);
 387          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
 388          fArrayQDY->AddAt(his3D, bin);
 389          snprintf(hname,100,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
 390          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
 391          fArrayQDZ->AddAt(his3D, bin);
 392          snprintf(hname,100,"RMSY Pad%d Qmiddle%f",ipad, qmean);
 393          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
 394          fArrayQRMSY->AddAt(his3D, bin);
 395          snprintf(hname,100,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
 396          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
 397          fArrayQRMSZ->AddAt(his3D, bin);
 398       }
 399    }
 400
 401
 402
 403    if (GetDebugLevel() > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl;
 404    cout << "end of main constructor" << endl; // TO BE REMOVED
 405 }
 406
 407
 408 AliTPCcalibTracks::~AliTPCcalibTracks() {
 409    //
 410    // AliTPCcalibTracks destructor
 411    //
 412
 413   if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' destuctor called." << endl;
 414    Int_t length = 0;
 415
 416
 417    if (fResolY) {
 418      length = fResolY->GetEntriesFast();
 419      for (Int_t i = 0; i < length; i++){
 420        delete fResolY->At(i);
 421        delete fResolZ->At(i);
 422        delete fRMSY->At(i);
 423        delete fRMSZ->At(i);
 424      }
 425      delete fResolY;
 426      delete fResolZ;
 427      delete fRMSY;
 428      delete fRMSZ;
 429    }
 430
 431    if (fArrayQDY) {
 432      length = fArrayQDY->GetEntriesFast();
 433      for (Int_t i = 0; i < length; i++){
 434        delete fArrayQDY->At(i);
 435        delete fArrayQDZ->At(i);
 436        delete fArrayQRMSY->At(i);
 437        delete fArrayQRMSZ->At(i);
 438      }
 439    }
 440
 441
 442
 443    delete fArrayQDY;
 444    delete fArrayQDZ;
 445    delete fArrayQRMSY;
 446    delete fArrayQRMSZ;
 447
 448   delete fRejectedTracksHisto;
 449   delete fClusterCutHisto;
 450   if (fCalPadClusterPerPad)    delete fCalPadClusterPerPad;
 451   if (fCalPadClusterPerPadRaw) delete fCalPadClusterPerPadRaw;
 452   delete fHisDeltaY;    // THnSparse - delta Y
 453   delete fHisDeltaZ;    // THnSparse - delta Z
 454   delete fHisRMSY;      // THnSparse - rms Y
 455   delete fHisRMSZ;      // THnSparse - rms Z
 456   delete fHisQmax;      // THnSparse - qmax
 457   delete fHisQtot;      // THnSparse - qtot
 458
 459 }
 460
 461
 462
 463 void AliTPCcalibTracks::Process(AliTPCseed *track){
 464    //
 465    // To be called in the selector
 466    // first AcceptTrack is evaluated, then calls all the following analyse functions:
 467    // FillResolutionHistoLocal(track)
 468
 469    //
 470   if (GetDebugLevel() > 5) Info("Process","Starting to process the track...");
 471    Int_t accpetStatus = AcceptTrack(track);
 472    if (accpetStatus == 0) {
 473       FillResolutionHistoLocal(track);
 474    }
 475    else fRejectedTracksHisto->Fill(accpetStatus);
 476 }
 477
 478
 479
 480 Int_t AliTPCcalibTracks::GetBin(Float_t q, Int_t pad){
 481   //
 482   // calculate bins for given q and pad type
 483   // used in TObjArray
 484   //
 485   Int_t res = TMath::Max( TMath::Nint((TMath::Sqrt(q) - 3.)), 0 );
 486   res *= 3;
 487   res += pad;
 488   return res;
 489 }
 490
 491
 492 Int_t AliTPCcalibTracks::GetBin(Int_t iq, Int_t pad){
 493   //
 494   // calculate bins for given iq and pad type
 495   // used in TObjArray
 496   //
 497   return iq * 3 + pad;;
 498 }
 499
 500
 501 Float_t AliTPCcalibTracks::GetQ(Int_t bin){
 502    //
 503    // returns to bin belonging charge
 504    // (bin / 3 + 3)^2
 505    //
 506    Int_t bin0 = bin / 3;
 507    bin0 += 3;
 508    return bin0 * bin0;
 509 }
 510
 511
 512 Float_t AliTPCcalibTracks::GetPad(Int_t bin){
 513    //
 514    // returns to bin belonging pad
 515    // bin % 3
 516    //
 517    return bin % 3;
 518 }
 519
 520
 521
 522 Int_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){
 523   //
 524   // Function, that decides wheather a given track is accepted for
 525   // the analysis or not.
 526   // The cuts are specified in the AliTPCcalibTracksCuts object 'fCuts'
 527   // Returns 0 if a track is accepted or an integer different from 0
 528   // to indicate the failed cut
 529   //
 530   const Int_t   kMinClusters  = fCuts->GetMinClusters();
 531   const Float_t kMinRatio     = fCuts->GetMinRatio();
 532   const Float_t kMax1pt       = fCuts->GetMax1pt();
 533   const Float_t kEdgeYXCutNoise    = fCuts->GetEdgeYXCutNoise();
 534   const Float_t kEdgeThetaCutNoise = fCuts->GetEdgeThetaCutNoise();
 535
 536   //
 537   // edge induced noise tracks - NEXT RELEASE will be removed during tracking
 538   if ( TMath::Abs(track->GetY() / track->GetX()) > kEdgeYXCutNoise )
 539     if ( TMath::Abs(track->GetTgl()) < kEdgeThetaCutNoise ) return 1;
 540   if (track->GetNumberOfClusters() < kMinClusters) return 2;
 541   Float_t ratio = track->GetNumberOfClusters() / (track->GetNFoundable() + 1.);
 542   if (ratio < kMinRatio) return 3;
 543 //   Float_t mpt = track->Get1Pt();       // Get1Pt() doesn't exist any more
 544   Float_t mpt = track->GetSigned1Pt();
 545   if (TMath::Abs(mpt) > kMax1pt) return 4;
 546   //if (TMath::Abs(track->GetZ())>240.) return kFALSE;
 547   //if (TMath::Abs(track->GetZ())<10.) return kFALSE;
 548   //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
 549
 550   if (GetDebugLevel() > 20) Info("AcceptTrack","Track has been accepted.");
 551   return 0;
 552 }
 553
 554
 555 void  AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
 556    //
 557    // fill resolution histograms - localy - tracklet in the neighborhood
 558    // write debug information to 'TPCSelectorDebug.root'
 559    //
 560    // _ the main function, called during track analysis _
 561    //
 562    // loop over all padrows along the track
 563    // fit tracklets (length: 13 clusters) calculate mean chi^2 for this track-fit in Y and Z direction
 564    //
 565    // loop again over all padrows along the track
 566    // fit tracklet (clusters in given padrow +- kDelta padrows)
 567    // with polynom of 2nd order and two polynoms of 1st order
 568    // take both polynoms of 1st order, calculate difference of their parameters
 569    // add covariance matrixes and calculate chi2 of this difference
 570    // if this chi2 is bigger than a given threshold, assume that the current cluster is
 571    // a kink an goto next padrow
 572    // if not:
 573    // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fResolY, fResolZ, fArrayQDY, fArrayQDY
 574    //
 575    // write debug information to 'TPCSelectorDebug.root'
 576    // only for every kDeltaWriteDebugStream'th padrow to reduce data volume
 577    // and to avoid redundant data
 578    //
 579
 580   /*  {//SG
 581     static long n=0;
 582     if( n==0 ){
 583       if( fClusterParam && fClusterParam->GetWaveCorrectionMap() ){
 584         cout<<"Map found "<<endl;
 585       }else{
 586         cout<<"Can't find the map "<<endl;
 587       }
 588     }
 589     if( n%100 ==0 )cout<<n<<" Tracks processed"<<endl;
 590     n++;
 591   }*/
 592   static TLinearFitter fFitterParY(3,"pol2");    //
 593   static TLinearFitter fFitterParZ(3,"pol2");    //
 594   static TLinearFitter fFitterParYWeight(3,"pol2");    //
 595   static TLinearFitter fFitterParZWeight(3,"pol2");    //
 596   fFitterParY.StoreData(kTRUE);
 597   fFitterParZ.StoreData(kTRUE);
 598   fFitterParYWeight.StoreData(kTRUE);
 599   fFitterParZWeight.StoreData(kTRUE);
 600   if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****");
 601   const Int_t   kDelta     = 10;          // delta rows to fit
 602   const Float_t kMinRatio  = 0.75;        // minimal ratio
 603   const Float_t kChi2Cut   =  10.;          // cut chi2 - left right
 604   const Float_t kSigmaCut  = 3.;        //sigma cut
 605   const Float_t kdEdxCut=300;
 606   const Float_t kPtCut=0.300;
 607
 608   if (track->GetTPCsignal()>kdEdxCut) return;
 609   if (TMath::Abs(AliTracker::GetBz())>0.1 &&TMath::Abs(track->Pt())<kPtCut) return;
 610
 611   // estimate mean error
 612   Int_t nTrackletsAll = 0;       // number of tracklets for given track
 613   Float_t csigmaY     = 0;       // mean sigma for tracklet refit in Y direction
 614   Float_t csigmaZ     = 0;       // mean sigma for tracklet refit in Z direction
 615   Int_t nClusters     = 0;       // working variable, number of clusters per tracklet
 616   Int_t sectorG       = -1;      // working variable, sector of tracklet, has to stay constant for one tracklet
 617   Double_t refX=0;
 618   // ---------------------------------------------------------------------
 619   for (Int_t irow = 0; irow < 159; irow++){
 620     // loop over all rows along the track
 621     // fit tracklets (length: 13 rows) with pol2 in Y and Z direction
 622     // calculate mean chi^2 for this track-fit in Y and Z direction
 623     AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
 624     if (!cluster0) continue;  // no cluster found
 625     Int_t sector = cluster0->GetDetector();
 626
 627     Int_t ipad = TMath::Nint(cluster0->GetPad());
 628     Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
 629     fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
 630
 631     if (sector != sectorG){
 632       // track leaves sector before it crossed enough rows to fit / initialization
 633       nClusters = 0;
 634       fFitterParY.ClearPoints();
 635       fFitterParZ.ClearPoints();
 636       sectorG = sector;
 637     }
 638     else {
 639       nClusters++;
 640       if (refX<1) refX=cluster0->GetX()+kDelta*0.5;
 641       Double_t x = cluster0->GetX()-refX;
 642       fFitterParY.AddPoint(&x, cluster0->GetY(), 1);
 643       fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1);
 644       //
 645       if ( nClusters >= kDelta + 3 ){
 646         // if more than 13 (kDelta+3) clusters were added to the fitters
 647         // fit the tracklet, increase trackletCounter
 648         fFitterParY.Eval();
 649         fFitterParZ.Eval();
 650         nTrackletsAll++;
 651         csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.);
 652         csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.);
 653         nClusters = -1;
 654         fFitterParY.ClearPoints();
 655         fFitterParZ.ClearPoints();
 656         refX=0;
 657       }
 658     }
 659   }      // for (Int_t irow = 0; irow < 159; irow++)
 660   // mean chi^2 for all tracklet fits in Y and in Z direction:
 661   csigmaY = TMath::Sqrt(TMath::Abs(csigmaY) / (nTrackletsAll+0.1));
 662   csigmaZ = TMath::Sqrt(TMath::Abs(csigmaZ) / (nTrackletsAll+0.1));
 663   if (csigmaY<=TMath::KUncertainty() || csigmaZ<= TMath::KUncertainty()) return;
 664   // ---------------------------------------------------------------------
 665   //
 666   //
 667
 668   for (Int_t irow = kDelta; irow < 159-kDelta; irow++){
 669     // loop again over all rows along the track
 670     // do analysis
 671     //
 672     Int_t nclFound = 0;  // number of clusters in the neighborhood
 673     Int_t ncl0 = 0;      // number of clusters in rows < rowOfCenterCluster
 674     Int_t ncl1 = 0;      // number of clusters in rows > rowOfCenterCluster
 675     AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
 676     if (!cluster0) continue;
 677     Int_t sector = cluster0->GetDetector();
 678     Float_t xref = cluster0->GetX();
 679
 680     // Make Fit
 681     fFitterParY.ClearPoints();
 682     fFitterParZ.ClearPoints();
 683     fFitterParYWeight.ClearPoints();
 684     fFitterParZWeight.ClearPoints();
 685     // fit tracklet (clusters in given padrow +- kDelta padrows)
 686     // with polynom of 2nd order and two polynoms of 1st order
 687     // take both polynoms of 1st order, calculate difference of their parameters
 688     // add covariance matrixes and calculate chi2 of this difference
 689     // if this chi2 is bigger than a given threshold, assume that the current cluster is
 690     // a kink an goto next padrow
 691
 692     // first fit the track angle for wave correction
 693
 694     AliRieman riemanFitAngle( 2*kDelta ); //SG
 695
 696     if( fClusterParam && fClusterParam->GetWaveCorrectionMap() ){
 697       for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
 698         // loop over irow +- kDelta rows (neighboured rows)
 699         if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
 700         AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
 701         if (!currentCluster) continue;
 702         if (currentCluster->GetType() < 0) continue;
 703         if (currentCluster->GetDetector() != sector) continue;
 704         riemanFitAngle.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY,csigmaZ);
 705       }
 706       if( riemanFitAngle.GetN()>3 ) riemanFitAngle.Update();
 707     }
 708
 709     // do fit
 710
 711     AliRieman riemanFit(2*kDelta);
 712     AliRieman riemanFitW(2*kDelta);
 713     for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
 714       // loop over irow +- kDelta rows (neighboured rows)
 715       //
 716       //
 717       if (idelta == 0) continue;                                // don't use center cluster
 718       if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
 719       AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
 720       if (!currentCluster) continue;
 721       if (currentCluster->GetType() < 0) continue;
 722       if (currentCluster->GetDetector() != sector) continue;
 723       nclFound++;
 724       if (idelta < 0){
 725         ncl0++;
 726       }
 727       if (idelta > 0){
 728         ncl1++;
 729       }
 730       //SG!!!
 731       double dY = 0;
 732       if( fClusterParam ){
 733         Int_t padSize = 0;                          // short pads
 734         if (currentCluster->GetDetector() >= 36) {
 735           padSize = 1;                              // medium pads
 736           if (currentCluster->GetRow() > 63) padSize = 2; // long pads
 737         }
 738         dY = fClusterParam->GetWaveCorrection( padSize,
 739                                                currentCluster->GetZ(),
 740                                                currentCluster->GetMax(),
 741                                                currentCluster->GetPad(),
 742                                                riemanFitAngle.GetDYat(currentCluster->GetX())
 743                                                );
 744       }
 745       riemanFit.AddPoint(currentCluster->GetX(), currentCluster->GetY()-dY,currentCluster->GetZ(), csigmaY,csigmaZ);
 746       riemanFitW.AddPoint(currentCluster->GetX(), currentCluster->GetY()-dY,currentCluster->GetZ(), TMath::Abs(csigmaY*TMath::Sqrt(1+TMath::Abs(idelta))),TMath::Abs(csigmaZ*TMath::Sqrt(1+TMath::Abs(idelta))));
 747     }  // loop over neighbourhood for fitter filling
 748     if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
 749     if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
 750     if (nclFound < kDelta * kMinRatio) continue;    // if not enough clusters (7.5) found in neighbourhood goto next padrow
 751     riemanFit.Update();
 752     riemanFitW.Update();
 753     Double_t chi2R=TMath::Sqrt(riemanFit.CalcChi2()/(2*nclFound-5));
 754     Double_t chi2RW=TMath::Sqrt(riemanFitW.CalcChi2()/(2*nclFound-5));
 755     Double_t paramYR[3], paramZR[3];
 756     Double_t paramYRW[3], paramZRW[3];
 757     //
 758     paramYR[0]    = riemanFit.GetYat(xref);
 759     paramZR[0]    = riemanFit.GetZat(xref);
 760     paramYRW[0]   = riemanFitW.GetYat(xref);
 761     paramZRW[0]   = riemanFitW.GetZat(xref);
 762     //
 763     paramYR[1]    = riemanFit.GetDYat(xref);
 764     paramZR[1]    = riemanFit.GetDZat(xref);
 765     paramYRW[1]   = riemanFitW.GetDYat(xref);
 766     paramZRW[1]   = riemanFitW.GetDZat(xref);
 767     //
 768     Int_t reject=0;
 769     if (chi2R>kChi2Cut) reject+=1;
 770     if (chi2RW>kChi2Cut) reject+=2;
 771     if (TMath::Abs(paramYR[0]-paramYRW[0])>kSigmaCut*csigmaY) reject+=4;
 772     if (TMath::Abs(paramZR[0]-paramZRW[0])>kSigmaCut*csigmaZ) reject+=8;
 773     if (TMath::Abs(paramYR[1]-paramYRW[1])>csigmaY) reject+=16;
 774     if (TMath::Abs(paramZR[1]-paramZRW[1])>csigmaZ) reject+=32;
 775     //
 776     TTreeSRedirector *cstream = GetDebugStreamer();
 777     // get fit parameters from pol2 fit:
 778     Double_t tracky = paramYR[0];
 779     Double_t trackz = paramZR[0];
 780     Float_t  deltay = cluster0->GetY()-tracky;
 781     Float_t  deltaz = cluster0->GetZ()-trackz;
 782     Float_t  angley = paramYR[1];
 783     Float_t  anglez = paramZR[1];
 784
 785     Int_t padSize = 0;                          // short pads
 786     if (cluster0->GetDetector() >= 36) {
 787       padSize = 1;                              // medium pads
 788       if (cluster0->GetRow() > 63) padSize = 2; // long pads
 789     }
 790     Int_t ipad = TMath::Nint(cluster0->GetPad());
 791     if (cstream){
 792       Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
 793       (*cstream)<<"Resol0"<<
 794         "run="<<fRun<<              //  run number
 795         "event="<<fEvent<<          //  event number
 796         "time="<<fTime<<            //  time stamp of event
 797         "trigger="<<fTrigger<<      //  trigger
 798         "mag="<<fMagF<<             //  magnetic field
 799         "padSize="<<padSize<<
 800         //
 801         "reject="<<reject<<
 802         "cl.="<<cluster0<<          // cluster info
 803         "xref="<<xref<<             // cluster reference X
 804         //rieman fit
 805         "yr="<<paramYR[0]<<         // track position y - rieman fit
 806         "zr="<<paramZR[0]<<         // track position z - rieman fit
 807         "yrW="<<paramYRW[0]<<         // track position y - rieman fit - weight
 808         "zrW="<<paramZRW[0]<<         // track position z - rieman fit - weight
 809         "dyr="<<paramYR[1]<<         // track position y - rieman fit
 810         "dzr="<<paramZR[1]<<         // track position z - rieman fit
 811         "dyrW="<<paramYRW[1]<<         // track position y - rieman fit - weight
 812         "dzrW="<<paramZRW[1]<<         // track position z - rieman fit - weight
 813         //
 814         "angley="<<angley<<         // angle par fit
 815         "anglez="<<anglez<<         // angle par fit
 816         "zdr="<<zdrift<<            //
 817         "dy="<<deltay<<
 818         "dz="<<deltaz<<
 819         "sy="<<csigmaY<<
 820         "sz="<<csigmaZ<<
 821         "chi2R="<<chi2R<<
 822         "chi2RW="<<chi2RW<<
 823         "\n";
 824     }
 825     if (reject) continue;
 826
 827
 828     // =========================================
 829     // wirte collected information to histograms
 830     // =========================================
 831
 832     Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
 833     fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
 834
 835
 836     TH3F * his3 = 0;
 837     his3 = (TH3F*)fRMSY->At(padSize);
 838     if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaY2()) ));
 839     his3 = (TH3F*)fRMSZ->At(padSize);
 840     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaZ2()) ));
 841
 842     his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
 843     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaY2()) ));
 844     his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
 845     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaZ2()) ));
 846
 847
 848     his3 = (TH3F*)fResolY->At(padSize);
 849     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
 850     his3 = (TH3F*)fResolZ->At(padSize);
 851     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
 852     his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
 853     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
 854     his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
 855     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );
 856     //=============================================================================================
 857     //
 858     // Fill THN histograms
 859     //
 860     Double_t xvar[9];
 861     xvar[1]=padSize;   // pad type
 862     xvar[2]=cluster0->GetZ();  //
 863     xvar[3]=cluster0->GetMax();
 864
 865     xvar[0]=deltay;
 866     double cog = cluster0->GetPad()-Int_t(cluster0->GetPad());// distance to the center of the pad
 867     xvar[4] = cog;
 868     xvar[5]=angley;
 869
 870     if( TMath::Abs(cog-0.5)<1.e-8 ) xvar[4]=-1; // fill one-pad clusters in the underflow bin
 871     fHisDeltaY->Fill(xvar);
 872
 873     xvar[4]=cog;
 874     xvar[0]=TMath::Sqrt(cluster0->GetSigmaY2());
 875     fHisRMSY->Fill(xvar);
 876
 877     xvar[0]=deltaz;
 878     xvar[4]=cluster0->GetTimeBin()-Int_t(cluster0->GetTimeBin()); // distance to the center of the time bin
 879     xvar[5]=anglez;
 880     fHisDeltaZ->Fill(xvar);
 881     xvar[0]=TMath::Sqrt(cluster0->GetSigmaZ2());
 882     fHisRMSZ->Fill(xvar);
 883
 884   }    // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
 885 }  // FillResolutionHistoLocal(...)
 886
 887
 888
 889
 890
 891
 892
 893
 894 void  AliTPCcalibTracks::SetStyle() const {
 895    //
 896    // set style, can be called by all draw functions
 897    //
 898    gROOT->SetStyle("Plain");
 899    gStyle->SetFillColor(10);
 900    gStyle->SetPadColor(10);
 901    gStyle->SetCanvasColor(10);
 902    gStyle->SetStatColor(10);
 903    gStyle->SetPalette(1,0);
 904    gStyle->SetNumberContours(60);
 905 }
 906
 907
 908
 909 void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){
 910    //
 911    // all functions are called, that produce pictures
 912    // the histograms are written to the directory 'pathName'
 913    // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file
 914    // 'stat' is also the number of minEntries for MakeResPlotsQTree
 915    //
 916
 917   if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
 918   MakeResPlotsQTree(stat, pathName);
 919 }
 920
 921
 922
 923
 924 void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){
 925    //
 926    // Make tree with resolution parameters
 927    // the result is written to 'resol.root' in directory 'pathname'
 928    // file information are available in fileInfo
 929    // available variables in the tree 'Resol':
 930    //  Entries: number of entries for this resolution point
 931    //  nbins:   number of bins that were accumulated
 932    //  Dim:     direction, Dim==0: y-direction, Dim==1: z-direction
 933    //  Pad:     padSize; short, medium and long
 934    //  Length:  pad length, 0.75, 1, 1.5
 935    //  QMean:   mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
 936    //  Zc:      center of middle bin in drift direction
 937    //  Zm:      mean dirftlength for accumulated Delta-Histograms
 938    //  Zs:      width of driftlength bin
 939    //  AngleC:  center of middle bin in Angle-Direction
 940    //  AngleM:  mean angle for accumulated Delta-Histograms
 941    //  AngleS:  width of Angle-bin
 942    //  Resol:   sigma for gaus fit through Delta-Histograms
 943    //  Sigma:   error of sigma for gaus fit through Delta Histograms
 944    //  MeanR:   mean of the Delta-Histogram
 945    //  SigmaR:  rms of the Delta-Histogram
 946    //  RMSm:    mean of the gaus fit through RMS-Histogram
 947    //  RMS:     sigma of the gaus fit through RMS-Histogram
 948    //  RMSe0:   error of mean of gaus fit in RMS-Histogram
 949    //  RMSe1:   error of sigma of gaus fit in RMS-Histogram
 950    //
 951
 952   if (GetDebugLevel() > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl;
 953   if (GetDebugLevel() > -1) cout << "    relax, the calculation will take a while..." << endl;
 954
 955    gSystem->MakeDirectory(pathName);
 956    gSystem->ChangeDirectory(pathName);
 957    TString kFileName = "resol.root";
 958    TTreeSRedirector fTreeResol(kFileName.Data());
 959
 960    TH3F *resArray[2][3][11];
 961    TH3F *rmsArray[2][3][11];
 962
 963    // load histograms from fArrayQDY and fArrayQDZ
 964    // into resArray and rmsArray
 965    // that is all we need here
 966    for (Int_t idim = 0; idim < 2; idim++){
 967       for (Int_t ipad = 0; ipad < 3; ipad++){
 968          for (Int_t iq = 0; iq <= 10; iq++){
 969             rmsArray[idim][ipad][iq]=0;
 970             resArray[idim][ipad][iq]=0;
 971             Int_t bin = GetBin(iq,ipad);
 972             TH3F *hresl = 0;
 973             if (idim == 0) hresl = (TH3F*)fArrayQDY->At(bin);
 974             if (idim == 1) hresl = (TH3F*)fArrayQDZ->At(bin);
 975             if (!hresl) continue;
 976             resArray[idim][ipad][iq] = (TH3F*) hresl->Clone();
 977             resArray[idim][ipad][iq]->SetDirectory(0);
 978             TH3F * hreslRMS = 0;
 979             if (idim == 0) hreslRMS = (TH3F*)fArrayQRMSY->At(bin);
 980             if (idim == 1) hreslRMS = (TH3F*)fArrayQRMSZ->At(bin);
 981             if (!hreslRMS) continue;
 982             rmsArray[idim][ipad][iq] = (TH3F*) hreslRMS->Clone();
 983             rmsArray[idim][ipad][iq]->SetDirectory(0);
 984          }
 985       }
 986    }
 987    if (GetDebugLevel() > -1) cout << "Histograms loaded, starting to proces..." << endl;
 988
 989    //--------------------------------------------------------------------------------------------
 990
 991    char name[200];
 992    Double_t qMean;
 993    Double_t zMean, angleMean, zCenter, angleCenter;
 994    Double_t zSigma, angleSigma;
 995    TH1D *projectionRes = new TH1D("projectionRes", "projectionRes", 50, -1, 1);
 996    TH1D *projectionRms = new TH1D("projectionRms", "projectionRms", 50, -1, 1);
 997    TF1 *fitFunction = new TF1("fitFunction", "gaus");
 998    Float_t entriesQ = 0;
 999    Int_t loopCounter = 1;
1000
1001    for (Int_t idim = 0; idim < 2; idim++){
1002       // Loop y-z corrdinate
1003       for (Int_t ipad = 0; ipad < 3; ipad++){
1004          // loop pad type
1005          for (Int_t iq = -1; iq < 10; iq++){
1006             // LOOP Q
1007            if (GetDebugLevel() > -1)
1008                cout << "Loop-counter, this is loop " << loopCounter << " of 66, ("
1009                   << (Int_t)((loopCounter)/66.*100) << "% done), "
1010                   << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << "  \r" << std::flush;
1011             loopCounter++;
1012             TH3F *hres = 0;
1013             TH3F *hrms = 0;
1014             qMean    = 0;
1015             entriesQ = 0;
1016
1017             // calculate qMean
1018             if (iq == -1){
1019                // integrated spectra
1020                for (Int_t iql = 0; iql < 10; iql++){
1021                   Int_t bin = GetBin(iql,ipad);
1022                   TH3F *hresl = resArray[idim][ipad][iql];
1023                   TH3F *hrmsl = rmsArray[idim][ipad][iql];
1024                   if (!hresl) continue;
1025                   if (!hrmsl) continue;
1026                   entriesQ += hresl->GetEntries();
1027                   qMean += hresl->GetEntries() * GetQ(bin);
1028                   if (!hres) {
1029                      hres = (TH3F*)hresl->Clone();
1030                      hrms = (TH3F*)hrmsl->Clone();
1031                   }
1032                   else{
1033                      hres->Add(hresl);
1034                      hrms->Add(hrmsl);
1035                   }
1036                }
1037                qMean /= entriesQ;
1038                qMean *= -1.;  // integral mean charge
1039             }
1040             else {
1041                // loop over neighboured Q-bins
1042                // accumulate entries from neighboured Q-bins
1043                for (Int_t iql = iq - 1; iql <= iq + 1; iql++){
1044                   if (iql < 0) continue;
1045                   Int_t bin = GetBin(iql,ipad);
1046                   TH3F * hresl = resArray[idim][ipad][iql];
1047                   TH3F * hrmsl = rmsArray[idim][ipad][iql];
1048                   if (!hresl) continue;
1049                   if (!hrmsl) continue;
1050                   entriesQ += hresl->GetEntries();
1051                   qMean += hresl->GetEntries() * GetQ(bin);
1052                   if (!hres) {
1053                      hres = (TH3F*) hresl->Clone();
1054                      hrms = (TH3F*) hrmsl->Clone();
1055                   }
1056                   else{
1057                      hres->Add(hresl);
1058                      hrms->Add(hrmsl);
1059                   }
1060                }
1061                qMean/=entriesQ;
1062             }
1063             if (!hres) continue;
1064             if (!hrms) continue;
1065
1066             TAxis *xAxisDriftLength = hres->GetXaxis();   // driftlength / z - axis
1067             TAxis *yAxisAngle       = hres->GetYaxis();   // angle axis
1068             TAxis *zAxisDelta       = hres->GetZaxis();   // delta axis
1069             TAxis *zAxisRms         = hrms->GetZaxis();   // rms axis
1070
1071             // loop over all angle bins
1072             for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) {
1073                angleCenter = yAxisAngle->GetBinCenter(ibinyAngle);
1074                // loop over all driftlength bins
1075                for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) {
1076                   zCenter    = xAxisDriftLength->GetBinCenter(ibinxDL);
1077                   zSigma     = xAxisDriftLength->GetBinWidth(ibinxDL);
1078                   angleSigma = yAxisAngle->GetBinWidth(ibinyAngle);
1079                   zMean      = zCenter;      // changens, when more statistic is accumulated
1080                   angleMean  = angleCenter;  // changens, when more statistic is accumulated
1081
1082                   // create 2 1D-Histograms, projectionRes and projectionRms
1083                   // these histograms are delta histograms for given direction, padSize, chargeBin,
1084                   // angleBin and driftLengthBin
1085                   // later on they will be fitted with a gausian, its sigma is the resoltuion...
1086                   snprintf(name,200,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter);
1087                   // TH1D * projectionRes = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
1088                   projectionRes->SetNameTitle(name, name);
1089                   snprintf(name,200,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter);
1090                   // TH1D * projectionRms =  new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
1091                   projectionRms->SetNameTitle(name, name);
1092
1093                   projectionRes->Reset();
1094                   projectionRes->SetBins(zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
1095                   projectionRms->Reset();
1096                   projectionRms->SetBins(zAxisRms->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
1097                   projectionRes->SetDirectory(0);
1098                   projectionRms->SetDirectory(0);
1099
1100                   Double_t entries = 0;
1101                   Int_t    nbins   = 0;   // counts, how many bins were accumulated
1102                   zMean     = 0;
1103                   angleMean = 0;
1104                   entries   = 0;
1105
1106                   // fill projectionRes and projectionRms for given dim, ipad and iq,
1107                   // as well as for given angleBin and driftlengthBin
1108                   // if this gives not enough statistic, include neighbourhood
1109                   // (angle and driftlength) successifely
1110                   for (Int_t dbin = 0; dbin <= 8; dbin++){              // delta-bins around centered angleBin and driftlengthBin
1111                      for (Int_t dbiny2 = -1; dbiny2 <= 1; dbiny2++) {   // delta-bins in angle direction
1112                         for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ // delta-bins in driftlength direction
1113                            if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue;   // add each bin only one time !
1114                            Int_t binx2 = ibinxDL + dbinx2;                       // position variable in x (driftlength) direction
1115                            Int_t biny2 = ibinyAngle + dbiny2;                    // position variable in y (angle)  direction
1116                            if (binx2 < 1 || biny2 < 1) continue;                 // don't go out of the histogram!
1117                            if (binx2 >= xAxisDriftLength->GetNbins()) continue;  // don't go out of the histogram!
1118                            if (biny2 >= yAxisAngle->GetNbins()) continue;        // don't go out of the histogram!
1119                            nbins++;                                              // count the number of accumulated bins
1120                            // Fill resolution histo
1121                            for (Int_t ibin3 = 1; ibin3 < zAxisDelta->GetNbins(); ibin3++) {
1122                               // Int_t content = (Int_t)hres->GetBinContent(binx2, biny2, ibin3);     // unused variable
1123                               projectionRes->Fill(zAxisDelta->GetBinCenter(ibin3), hres->GetBinContent(binx2, biny2, ibin3));
1124                               entries   += hres->GetBinContent(binx2, biny2, ibin3);
1125                               zMean     += hres->GetBinContent(binx2, biny2, ibin3) * xAxisDriftLength->GetBinCenter(binx2);
1126                               angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yAxisAngle->GetBinCenter(biny2);
1127                            }  // ibin3 loop
1128                            // fill RMS histo
1129                            for (Int_t ibin3 = 1; ibin3 < zAxisRms->GetNbins(); ibin3++) {
1130                               projectionRms->Fill(zAxisRms->GetBinCenter(ibin3), hrms->GetBinContent(binx2, biny2, ibin3));
1131                            }
1132                         }  //dbinx2 loop
1133                         if (entries > minEntries) break; // enough statistic accumulated
1134                      }  // dbiny2 loop
1135                      if (entries > minEntries) break;    // enough statistic accumulated
1136                   }  // dbin loop
1137                   if ( entries< minEntries) continue;  // when it was absolutly impossible to get enough statistic, don't write this point into the resolution tree
1138                   zMean /= entries;
1139                   angleMean /= entries;
1140
1141                   if (entries > minEntries) {
1142                      //  when enough statistic is accumulated
1143                      //  fit Delta histograms with a gausian
1144                      //  of the gausian is the resolution (resol), its fit error is sigma
1145                      //  store also mean and RMS of the histogram
1146                      Float_t xmin     = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
1147                      Float_t xmax     = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
1148
1149 //                      projectionRes->Fit("gaus", "q0", "", xmin, xmax);
1150 //                      Float_t resol    = projectionRes->GetFunction("gaus")->GetParameter(2);
1151 //                      Float_t sigma    = projectionRes->GetFunction("gaus")->GetParError(2);
1152                      fitFunction->SetMaximum(xmax);
1153                      fitFunction->SetMinimum(xmin);
1154                      projectionRes->Fit("fitFunction", "qN0", "", xmin, xmax);
1155                      Float_t resol    = fitFunction->GetParameter(2);
1156                      Float_t sigma    = fitFunction->GetParError(2);
1157
1158                      Float_t meanR    = projectionRes->GetMean();
1159                      Float_t sigmaR   = projectionRes->GetRMS();
1160                      // fit also RMS histograms with a gausian
1161                      // store mean and sigma of the gausian in rmsMean and rmsSigma
1162                      // store also the fit errors in errorRMS and errorSigma
1163                      xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
1164                      xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
1165
1166 //                      projectionRms->Fit("gaus","q0","",xmin,xmax);
1167 //                      Float_t rmsMean    = projectionRms->GetFunction("gaus")->GetParameter(1);
1168 //                      Float_t rmsSigma   = projectionRms->GetFunction("gaus")->GetParameter(2);
1169 //                      Float_t errorRMS   = projectionRms->GetFunction("gaus")->GetParError(1);
1170 //                      Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2);
1171                      projectionRms->Fit("fitFunction", "qN0", "", xmin, xmax);
1172                      Float_t rmsMean    = fitFunction->GetParameter(1);
1173                      Float_t rmsSigma   = fitFunction->GetParameter(2);
1174                      Float_t errorRMS   = fitFunction->GetParError(1);
1175                      Float_t errorSigma = fitFunction->GetParError(2);
1176
1177                      Float_t length = 0.75;
1178                      if (ipad == 1) length = 1;
1179                      if (ipad == 2) length = 1.5;
1180
1181                      fTreeResol<<"Resol"<<
1182                         "Entries="<<entries<<      // number of entries for this resolution point
1183                         "nbins="<<nbins<<          // number of bins that were accumulated
1184                         "Dim="<<idim<<             // direction, Dim==0: y-direction, Dim==1: z-direction
1185                         "Pad="<<ipad<<             // padSize; short, medium and long
1186                         "Length="<<length<<        // pad length, 0.75, 1, 1.5
1187                         "QMean="<<qMean<<          // mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
1188                         "Zc="<<zCenter<<           // center of middle bin in drift direction
1189                         "Zm="<<zMean<<             // mean dirftlength for accumulated Delta-Histograms
1190                         "Zs="<<zSigma<<            // width of driftlength bin
1191                         "AngleC="<<angleCenter<<   // center of middle bin in Angle-Direction
1192                         "AngleM="<<angleMean<<     // mean angle for accumulated Delta-Histograms
1193                         "AngleS="<<angleSigma<<    // width of Angle-bin
1194                         "Resol="<<resol<<          // sigma for gaus fit through Delta-Histograms
1195                         "Sigma="<<sigma<<          // error of sigma for gaus fit through Delta Histograms
1196                         "MeanR="<<meanR<<          // mean of the Delta-Histogram
1197                         "SigmaR="<<sigmaR<<        // rms of the Delta-Histogram
1198                         "RMSm="<<rmsMean<<         // mean of the gaus fit through RMS-Histogram
1199                         "RMSs="<<rmsSigma<<        // sigma of the gaus fit through RMS-Histogram
1200                         "RMSe0="<<errorRMS<<       // error of mean of gaus fit in RMS-Histogram
1201                         "RMSe1="<<errorSigma<<     // error of sigma of gaus fit in RMS-Histogram
1202                         "\n";
1203                      if (GetDebugLevel() > 5) {
1204                         projectionRes->SetDirectory(fTreeResol.GetFile());
1205                         projectionRes->Write(projectionRes->GetName());
1206                         projectionRes->SetDirectory(0);
1207                         projectionRms->SetDirectory(fTreeResol.GetFile());
1208                         projectionRms->Write(projectionRms->GetName());
1209                         projectionRes->SetDirectory(0);
1210                      }
1211                   }  // if (projectionRes->GetSum() > minEntries)
1212                }  // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++)
1213             }  // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++)
1214
1215          }  // iq-loop
1216       }  // ipad-loop
1217    }  // idim-loop
1218    delete projectionRes;
1219    delete projectionRms;
1220
1221 //    TFile resolFile(fTreeResol.GetFile());
1222    TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries));
1223    fileInfo.Write("fileInfo");
1224 //    resolFile.Close();
1225 //    fTreeResol.GetFile()->Close();
1226    if (GetDebugLevel() > -1) cout << endl;
1227    if (GetDebugLevel() > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl;
1228    gSystem->ChangeDirectory("..");
1229 }
1230
1231
1232
1233
1234
1235 Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
1236    //
1237    // function to merge several AliTPCcalibTracks objects after PROOF calculation
1238    // The object's histograms are merged via their merge functions
1239    // Be carefull: histograms are linked to a file, switch this off by TH1::AddDirectory(kFALSE) !!!
1240    //
1241
1242   if (GetDebugLevel() > 0) cout << " *****  this is AliTPCcalibTracks::Merge(TCollection *collectionList)  *****"<< endl;
1243    if (!collectionList) return 0;
1244    if (collectionList->IsEmpty()) return -1;
1245
1246    if (GetDebugLevel() > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl;     //    REMOVE THIS LINE!!!!!!!!!!!!!!!!!1
1247    if (GetDebugLevel() > 5) cout << " the list in the merge-function looks as follows: " << endl;
1248    collectionList->Print();
1249
1250    // create a list for each data member
1251    TList* deltaYList = new TList;
1252    TList* deltaZList = new TList;
1253    TList* arrayAmpRowList = new TList;
1254    TList* rejectedTracksList = new TList;
1255    TList* clusterCutHistoList = new TList;
1256    TList* arrayAmpList = new TList;
1257    TList* arrayQDYList = new TList;
1258    TList* arrayQDZList = new TList;
1259    TList* arrayQRMSYList = new TList;
1260    TList* arrayQRMSZList = new TList;
1261    TList* resolYList = new TList;
1262    TList* resolZList = new TList;
1263    TList* rMSYList = new TList;
1264    TList* rMSZList = new TList;
1265
1266 //    TList* nRowsList = new TList;
1267 //    TList* nSectList = new TList;
1268 //    TList* fileNoList = new TList;
1269
1270    TIterator *listIterator = collectionList->MakeIterator();
1271    AliTPCcalibTracks *calibTracks = 0;
1272    if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;
1273    Int_t counter = 0;
1274    while ( (calibTracks = dynamic_cast<AliTPCcalibTracks*> (listIterator->Next())) ){
1275       // loop over all entries in the collectionList and get dataMembers into lists
1276
1277       arrayQDYList->Add(calibTracks->GetfArrayQDY());
1278       arrayQDZList->Add(calibTracks->GetfArrayQDZ());
1279       arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
1280       arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ());
1281       resolYList->Add(calibTracks->GetfResolY());
1282       resolZList->Add(calibTracks->GetfResolZ());
1283       rMSYList->Add(calibTracks->GetfRMSY());
1284       rMSZList->Add(calibTracks->GetfRMSZ());
1285       rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
1286       clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
1287       //
1288       if (fCalPadClusterPerPad && calibTracks->GetfCalPadClusterPerPad())
1289         fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
1290       //      fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
1291       counter++;
1292       if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
1293       AddHistos(calibTracks);
1294    }
1295
1296
1297    // merge data members
1298    if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl;
1299    fClusterCutHisto->Merge(clusterCutHistoList);
1300    fRejectedTracksHisto->Merge(rejectedTracksList);
1301
1302    TObjArray* objarray = 0;
1303    TH1* hist = 0;
1304    TList* histList = 0;
1305    TIterator *objListIterator = 0;
1306
1307
1308    if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
1309    // merge fArrayQDY
1310    for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
1311       objListIterator = arrayQDYList->MakeIterator();
1312       histList = new TList;
1313       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
1314          // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F
1315          hist = (TH3F*)(objarray->At(i));
1316          histList->Add(hist);
1317       }
1318       ((TH3F*)(fArrayQDY->At(i)))->Merge(histList);
1319       delete histList;
1320       delete objListIterator;
1321    }
1322
1323    if (GetDebugLevel() > 0) cout << "merging fArrayQDZ..." << endl;
1324    // merge fArrayQDZ
1325    for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300
1326       objListIterator = arrayQDZList->MakeIterator();
1327       histList = new TList;
1328       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
1329          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1330          hist = (TH3F*)(objarray->At(i));
1331          histList->Add(hist);
1332       }
1333       ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList);
1334       delete histList;
1335       delete objListIterator;
1336    }
1337
1338    if (GetDebugLevel() > 0) cout << "merging fArrayQRMSY..." << endl;
1339    // merge fArrayQRMSY
1340    for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300
1341       objListIterator = arrayQRMSYList->MakeIterator();
1342       histList = new TList;
1343       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
1344          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1345          //   if (!objarray) continue; // removed for coverity -> JMT
1346          hist = (TH3F*)(objarray->At(i));
1347          histList->Add(hist);
1348       }
1349       ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList);
1350       delete histList;
1351       delete objListIterator;
1352    }
1353
1354    if (GetDebugLevel() > 0) cout << "merging fArrayQRMSZ..." << endl;
1355    // merge fArrayQRMSZ
1356    for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300
1357       objListIterator = arrayQRMSZList->MakeIterator();
1358       histList = new TList;
1359       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
1360          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1361          hist = (TH3F*)(objarray->At(i));
1362          histList->Add(hist);
1363       }
1364       ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList);
1365       delete histList;
1366       delete objListIterator;
1367    }
1368
1369
1370
1371
1372
1373
1374    if (GetDebugLevel() > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl;
1375    // merge fResolY
1376    for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3
1377       objListIterator = resolYList->MakeIterator();
1378       histList = new TList;
1379       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
1380          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1381          hist = (TH3F*)(objarray->At(i));
1382          histList->Add(hist);
1383       }
1384       ((TH3F*)(fResolY->At(i)))->Merge(histList);
1385       delete histList;
1386       delete objListIterator;
1387    }
1388
1389    // merge fResolZ
1390    for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3
1391       objListIterator = resolZList->MakeIterator();
1392       histList = new TList;
1393       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
1394          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1395          hist = (TH3F*)(objarray->At(i));
1396          histList->Add(hist);
1397       }
1398       ((TH3F*)(fResolZ->At(i)))->Merge(histList);
1399       delete histList;
1400       delete objListIterator;
1401    }
1402
1403    // merge fRMSY
1404    for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3
1405       objListIterator = rMSYList->MakeIterator();
1406       histList = new TList;
1407       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
1408          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1409          hist = (TH3F*)(objarray->At(i));
1410          histList->Add(hist);
1411       }
1412       ((TH3F*)(fRMSY->At(i)))->Merge(histList);
1413       delete histList;
1414       delete objListIterator;
1415    }
1416
1417    // merge fRMSZ
1418    for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3
1419       objListIterator = rMSZList->MakeIterator();
1420       histList = new TList;
1421       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
1422          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1423          hist = (TH3F*)(objarray->At(i));
1424          histList->Add(hist);
1425       }
1426       ((TH3F*)(fRMSZ->At(i)))->Merge(histList);
1427       delete histList;
1428       delete objListIterator;
1429    }
1430
1431    delete deltaYList;
1432    delete deltaZList;
1433    delete arrayAmpRowList;
1434    delete arrayAmpList;
1435    delete arrayQDYList;
1436    delete arrayQDZList;
1437    delete arrayQRMSYList;
1438    delete arrayQRMSZList;
1439    delete resolYList;
1440    delete resolZList;
1441    delete rMSYList;
1442    delete rMSZList;
1443    delete listIterator;
1444
1445    if (GetDebugLevel() > 0) cout << "merging done!" << endl;
1446
1447    return 1;
1448 }
1449
1450
1451
1452
1453 void AliTPCcalibTracks::MakeHistos(){
1454   //
1455   ////make THnSparse
1456   //
1457   //THnSparse  *fHisDeltaY;    // THnSparse - delta Y
1458   //THnSparse  *fHisDeltaZ;    // THnSparse - delta Z
1459   //THnSparse  *fHisRMSY;      // THnSparse - rms Y
1460   //THnSparse  *fHisRMSZ;      // THnSparse - rms Z
1461   //THnSparse  *fHisQmax;      // THnSparse - qmax
1462   //THnSparse  *fHisQtot;      // THnSparse - qtot
1463   // cluster  performance bins
1464   // 0 - variable of interest
1465   // 1 - pad type   - 0- short 1-medium 2-long pads
1466   // 2 - drift length - drift length -0-1
1467   // 3 - Qmax         - Qmax  - 2- 400
1468   // 4 - cog          - COG position - 0-1
1469   // 5 - tan(phi)     - local y angle
1470   // 6 - tan(theta)   - local z angle
1471   // 7 - sector       - sector number
1472   Double_t xminTrack[8], xmaxTrack[8];
1473   Int_t binsTrack[8];
1474   TString axisName[8];
1475
1476   //
1477   binsTrack[0]=200;
1478   axisName[0]  ="var";
1479
1480   binsTrack[1] =3;
1481   xminTrack[1] =0; xmaxTrack[1]=3;
1482   axisName[1]  ="pad type";
1483   //
1484   binsTrack[2] =20;
1485   xminTrack[2] =-250; xmaxTrack[2]=250;
1486   axisName[2]  ="z";
1487   //
1488   binsTrack[3] =10;
1489   xminTrack[3] =1; xmaxTrack[3]=400;
1490   axisName[3]  ="Qmax";
1491   //
1492   binsTrack[4] =20;
1493   xminTrack[4] =0; xmaxTrack[4]=1;
1494   axisName[4]  ="cog";
1495   //
1496   binsTrack[5] =15;
1497   xminTrack[5] =-1.5; xmaxTrack[5]=1.5;
1498   axisName[5]  ="tan(angle)";
1499   //
1500   //
1501   xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
1502   fHisDeltaY=new THnSparseF("#Delta_{y} (cm)","#Delta_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1503   xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
1504   fHisDeltaZ=new THnSparseF("#Delta_{z} (cm)","#Delta_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1505   xminTrack[0] =0.; xmaxTrack[0]=0.5;
1506   fHisRMSY=new THnSparseF("#RMS_{y} (cm)","#RMS_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1507   xminTrack[0] =0.; xmaxTrack[0]=0.5;
1508   fHisRMSZ=new THnSparseF("#RMS_{z} (cm)","#RMS_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1509   xminTrack[0] =0.; xmaxTrack[0]=100;
1510   fHisQmax=new THnSparseF("Qmax (ADC)","Qmax (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
1511
1512   xminTrack[0] =0.; xmaxTrack[0]=250;
1513   fHisQtot=new THnSparseF("Qtot (ADC)","Qtot (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
1514
1515
1516   for (Int_t ivar=0;ivar<6;ivar++){
1517     fHisDeltaY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1518     fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1519     fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1520     fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1521     fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
1522     fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
1523     fHisDeltaY->GetAxis(ivar)->SetTitle(axisName[ivar].Data());
1524     fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1525     fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1526     fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1527     fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
1528     fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
1529   }
1530
1531
1532   BinLogX(fHisDeltaY,3);
1533   BinLogX(fHisDeltaZ,3);
1534   BinLogX(fHisRMSY,3);
1535   BinLogX(fHisRMSZ,3);
1536   BinLogX(fHisQmax,3);
1537   BinLogX(fHisQtot,3);
1538
1539 }
1540
1541 void    AliTPCcalibTracks::AddHistos(AliTPCcalibTracks* calib){
1542   //
1543   // Add histograms
1544   //
1545   if (calib->fHisDeltaY) fHisDeltaY->Add(calib->fHisDeltaY);
1546   if (calib->fHisDeltaZ) fHisDeltaZ->Add(calib->fHisDeltaZ);
1547   if (calib->fHisRMSY)   fHisRMSY->Add(calib->fHisRMSY);
1548   if (calib->fHisRMSZ)   fHisRMSZ->Add(calib->fHisRMSZ);
1549 }
1550
1551
1552
1553 void AliTPCcalibTracks::MakeSummaryTree(THnSparse *hisInput, TTreeSRedirector *pcstream, Int_t ptype){
1554   //
1555   // Dump summary info
1556   //
1557   //  0.OBJ: TAxis     var
1558   //  1.OBJ: TAxis     pad type
1559   //  2.OBJ: TAxis     z
1560   //  3.OBJ: TAxis     Qmax
1561   //  4.OBJ: TAxis     cog
1562   //  5.OBJ: TAxis     tan(angle)
1563   //
1564   if (ptype>3) return;
1565   Int_t idim[6]={0,1,2,3,4,5};
1566   TString hname[4]={"dy","dz","rmsy","rmsz"};
1567   //
1568   Int_t nbins5=hisInput->GetAxis(5)->GetNbins();
1569   Int_t first5=hisInput->GetAxis(5)->GetFirst();
1570   Int_t last5 =hisInput->GetAxis(5)->GetLast();
1571   //
1572   for (Int_t ibin5=first5-1; ibin5<=last5; ibin5+=1){   // axis 5 - local angle
1573     Bool_t bin5All=kTRUE;
1574     if (ibin5>=first5){
1575       hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5-1,first5),TMath::Min(ibin5+1,nbins5));
1576       if (ibin5==first5) hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5,first5),TMath::Min(ibin5,nbins5));
1577       bin5All=kFALSE;
1578     }
1579     Double_t      x5= hisInput->GetAxis(5)->GetBinCenter(ibin5);
1580     THnSparse * his5= hisInput->Projection(5,idim);         //projected histogram according selection 5
1581     Int_t nbins4=his5->GetAxis(4)->GetNbins();
1582     Int_t first4=his5->GetAxis(4)->GetFirst();
1583     Int_t last4 =his5->GetAxis(4)->GetLast();
1584
1585     for (Int_t ibin4=first4-1; ibin4<=last4; ibin4+=1){   // axis 4 - cog
1586       Bool_t bin4All=kTRUE;
1587       if (ibin4>=first4){
1588         his5->GetAxis(4)->SetRange(TMath::Max(ibin4+1,first4),TMath::Min(ibin4-1,nbins4));
1589         if (ibin4==first4||ibin4==last4) his5->GetAxis(4)->SetRange(TMath::Max(ibin4,first4),TMath::Min(ibin4,nbins4));
1590         bin4All=kFALSE;
1591       }
1592       Double_t      x4= his5->GetAxis(4)->GetBinCenter(ibin4);
1593       THnSparse * his4= his5->Projection(4,idim);         //projected histogram according selection 4
1594       //
1595       Int_t nbins3=his4->GetAxis(3)->GetNbins();
1596       Int_t first3=his4->GetAxis(3)->GetFirst();
1597       Int_t last3 =his4->GetAxis(3)->GetLast();
1598       //
1599       for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){   // axis 3 - Qmax
1600         Bool_t bin3All=kTRUE;
1601         if (ibin3>=first3){
1602           his4->GetAxis(3)->SetRange(TMath::Max(ibin3,first3),TMath::Min(ibin3,nbins3));
1603           bin3All=kFALSE;
1604         }
1605         Double_t      x3= his4->GetAxis(3)->GetBinCenter(ibin3);
1606         THnSparse * his3= his4->Projection(3,idim);         //projected histogram according selection 3
1607         //
1608         Int_t nbins2    = his3->GetAxis(2)->GetNbins();
1609         Int_t first2    = his3->GetAxis(2)->GetFirst();
1610         Int_t last2     = his3->GetAxis(2)->GetLast();
1611         //
1612         for (Int_t ibin2=first2-1; ibin2<=last2; ibin2+=1){   // axis 2 - z
1613           Bool_t bin2All=kTRUE;
1614           Double_t      x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1615           if (ibin2>=first2){
1616             his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,first2),TMath::Min(ibin2+1,nbins2));
1617             if (ibin2==first2||ibin2==last2||TMath::Abs(x2)<20) his3->GetAxis(2)->SetRange(TMath::Max(ibin2,first2),TMath::Min(ibin2,nbins2));
1618             bin2All=kFALSE;
1619           }
1620           THnSparse * his2= his3->Projection(2,idim);         //projected histogram according selection 2
1621           //
1622           Int_t nbins1     = his2->GetAxis(1)->GetNbins();
1623           Int_t first1     = his2->GetAxis(1)->GetFirst();
1624           Int_t last1      = his2->GetAxis(1)->GetLast();
1625           for (Int_t ibin1=first1-1; ibin1<=last1; ibin1++){   //axis 1 - pad type
1626             Bool_t bin1All=kTRUE;
1627             if (ibin1>=first1){
1628               his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1629               bin1All=kFALSE;
1630             }
1631             Double_t       x1= TMath::Nint(his2->GetAxis(1)->GetBinCenter(ibin1)-0.5);
1632             TH1 * hisDelta = his2->Projection(0);
1633             Double_t entries = hisDelta->GetEntries();
1634             Double_t mean=0, rms=0;
1635             if (entries>10){
1636               mean    = hisDelta->GetMean();
1637               rms = hisDelta->GetRMS();
1638               hisDelta->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
1639               mean    = hisDelta->GetMean();
1640               rms = hisDelta->GetRMS();
1641             }
1642             //
1643             (*pcstream)<<hname[ptype].Data()<<
1644               // flag - intgrated values for given bin
1645               "angleA="<<bin5All<<
1646               "cogA="<<bin4All<<
1647               "qmaxA="<<bin3All<<
1648               "zA="<<bin2All<<
1649               "ipadA="<<bin1All<<
1650               // center of bin value
1651               "angle="<<x5<<       // local angle
1652               "cog="<<x4<<         // distance cluster to center
1653               "qmax="<<x3<<        // qmax
1654               "z="<<x2<<           // z of the cluster
1655               "ipad="<<x1<<        // type of the region
1656               // mean values
1657               //
1658               "entries="<<entries<<
1659               "mean="<<mean<<
1660               "rms="<<rms<<
1661               "ptype="<<ptype<<   //
1662               "\n";
1663             delete hisDelta;
1664             printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n",x5,x4,x3,x2,x1, entries,mean);
1665           }//loop z
1666           delete his2;
1667         } // loop Q max
1668         delete his3;
1669       } // loop COG
1670       delete his4;
1671     }//loop local angle
1672     delete his5;
1673   }
1674 }
1675
1676
1677 int AliTPCcalibTracks::GetTHnStat(const  THnBase *H, THnBase *&Mean, THnBase *&Sigma, THnBase *&Entr )
1678 {
1679   // H is THnSparse:
1680   //
1681   // OBJ: TAxis     var     var
1682   // OBJ: TAxis     axis 1
1683   // OBJ: TAxis     axis 2
1684   // ...
1685
1686   // Output:
1687
1688   // Mean, Sigma and Entr is THnF of mean, RMS and entries of var:
1689   //
1690   // OBJ: TAxis     axis 1
1691   // OBJ: TAxis     axis 2
1692   // ..
1693
1694   Mean = 0;
1695   Sigma = 0;
1696   Entr = 0;
1697   if( !H ) return -1;
1698
1699   Bool_t ok = kTRUE;
1700
1701   int nDim = H->GetNdimensions();
1702   Long64_t nFilledBins = H->GetNbins();
1703   Long64_t nStatBins = 0;
1704
1705   Int_t *nBins = new Int_t [nDim];
1706   Double_t *xMin = new Double_t [nDim];
1707   Double_t *xMax = new Double_t [nDim];
1708   Int_t *idx = new Int_t [nDim];
1709
1710   TString nameMean = H->GetName();
1711   TString nameSigma = H->GetName();
1712   TString nameEntr = H->GetName();
1713   nameMean+="_Mean";
1714   nameSigma+="_Sigma";
1715   nameEntr+="_Entr";
1716
1717   ok = ok && ( nBins && xMin && xMax && idx );
1718
1719   if( ok ){
1720     for( int i=0; i<nDim; i++){
1721       xMin[i] = H->GetAxis(i)->GetXmin();
1722       xMax[i] = H->GetAxis(i)->GetXmax();
1723       nBins[i] = H->GetAxis(i)->GetNbins();
1724     }
1725
1726     Mean  = new THnSparseF( nameMean.Data(),  nameMean.Data(),  nDim-1, nBins+1, xMin+1, xMax+1 );
1727     Sigma = new THnSparseF( nameSigma.Data(), nameSigma.Data(), nDim-1, nBins+1, xMin+1, xMax+1 );
1728     Entr  = new THnSparseF( nameEntr.Data(),  nameEntr.Data(),  nDim-1, nBins+1, xMin+1, xMax+1 );
1729   }
1730
1731   ok = ok && ( Mean && Sigma && Entr );
1732
1733   if( ok ){
1734     for( int i=0; i<nDim-1; i++){
1735       const TAxis *ax = H->GetAxis(i+1);
1736       Mean->GetAxis(i)->SetName(ax->GetName());
1737       Mean->GetAxis(i)->SetTitle(ax->GetTitle());
1738       Sigma->GetAxis(i)->SetName(ax->GetName());
1739       Sigma->GetAxis(i)->SetTitle(ax->GetTitle());
1740       Entr->GetAxis(i)->SetName(ax->GetName());
1741       Entr->GetAxis(i)->SetTitle(ax->GetTitle());
1742       if( ax->GetXbins()->fN>0 ){
1743         Mean->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
1744         Sigma->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
1745         Entr->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
1746       }
1747     }
1748
1749     // Allocate bins
1750
1751     for( Long64_t binS=0; binS<nFilledBins; binS++){
1752       H->GetBinContent(binS,idx);
1753       Mean->GetBin(idx+1,kTRUE);
1754       Sigma->GetBin(idx+1,kTRUE);
1755       Entr->GetBin(idx+1,kTRUE);
1756     }
1757
1758     nStatBins = Mean->GetNbins();
1759
1760   }
1761
1762
1763   Long64_t *hMap = new Long64_t[nFilledBins];
1764   Double_t *hVal  = new Double_t[nFilledBins];
1765   Double_t *hEntr = new Double_t[nFilledBins];
1766   Double_t *meanD = new Double_t[nStatBins];
1767   Double_t *sigmaD = new Double_t[nStatBins];
1768
1769   ok = ok && ( hMap && hVal && hEntr && meanD && sigmaD );
1770
1771   // Map bins to Stat bins
1772
1773   if( ok ){
1774     for( Long64_t binS=0; binS<nFilledBins; binS++){
1775       double val = H->GetBinContent(binS,idx);
1776       Long64_t bin = Mean->GetBin(idx+1,kFALSE);
1777       ok = ok && ( bin>=0 && bin<nStatBins && bin==Sigma->GetBin(idx+1,kFALSE) && bin== Entr->GetBin(idx+1,kFALSE) );
1778       if( !ok ) break;
1779       if( val < 0 ){
1780         cout << "AliTPCcalibTracks: GetSparseStat : Unexpected content of the input THn"<<endl;
1781         bin = -1;
1782       }
1783       if( idx[0]<1 || idx[0]>nBins[0] ) bin = -1;
1784       hMap[binS] = bin;
1785       hVal[binS] = idx[0];
1786       hEntr[binS] = val;
1787     }
1788   }
1789
1790   // do 2 iteration with cut at 4 Sigma
1791
1792   for( int iter=0; ok && iter<2; iter++ ){
1793
1794     // clean up entries
1795
1796     for( Long64_t bin=0; bin < nStatBins; bin++){
1797       Entr->SetBinContent(bin, 0);
1798       meanD[bin]=0;
1799       sigmaD[bin]=0;
1800     }
1801
1802     // get Entries and Mean
1803
1804     for( Long64_t binS=0; binS<nFilledBins; binS++){
1805       Long64_t bin = hMap[binS];
1806       Double_t val  = hVal[binS];
1807       Double_t entr = hEntr[binS];
1808       if( bin < 0 ) continue;
1809       if( iter!=0 ){ // cut
1810         double s2 = Sigma->GetBinContent(bin);
1811         double d = val - Mean->GetBinContent(bin);
1812         if( d*d>s2*16 ) continue;
1813       }
1814       meanD[bin]+= entr*val;
1815       Entr->AddBinContent(bin,entr);
1816     }
1817
1818     for( Long64_t bin=0; bin<nStatBins; bin++){
1819       double val = meanD[bin];
1820       double sum = Entr->GetBinContent(bin);
1821       if( sum>=1 ){
1822         Mean->SetBinContent(bin, val/sum );
1823       }
1824       else Mean->SetBinContent(bin, 0);
1825       Entr->SetBinContent(bin, 0);
1826     }
1827
1828     // get RMS
1829
1830     for( Long64_t binS=0; binS<nFilledBins; binS++){
1831       Long64_t bin = hMap[binS];
1832       Double_t val  = hVal[binS];
1833       Double_t entr = hEntr[binS];
1834       if( bin < 0 ) continue;
1835       double d2 = val - Mean->GetBinContent(bin);
1836       d2 *= d2;
1837       if( iter!=0 ){ // cut
1838         double s2 = Sigma->GetBinContent(bin);
1839         if( d2>s2*16 ) continue;
1840       }
1841       sigmaD[bin] += entr*d2;
1842       Entr->AddBinContent(bin,entr);
1843     }
1844
1845     for( Long64_t bin=0; bin<nStatBins; bin++ ){
1846       double val = sigmaD[bin];
1847       double sum = Entr->GetBinContent(bin);
1848       if( sum>=1 && val>=0 ){
1849         Sigma->SetBinContent(bin, val/sum );
1850       }
1851       else Sigma->SetBinContent(bin, 0);
1852     }
1853   }
1854
1855   // scale the Mean and the Sigma
1856
1857   if( ok ){
1858     double v0 = H->GetAxis(0)->GetBinCenter(1);
1859     double vscale = H->GetAxis(0)->GetBinWidth(1);
1860
1861     for( Long64_t bin=0; bin<nStatBins; bin++){
1862       double m = Mean->GetBinContent(bin);
1863       double s2 = Sigma->GetBinContent(bin);
1864       double sum = Entr->GetBinContent(bin);
1865       if( sum>0 && s2>=0 ){
1866         Mean->SetBinContent(bin, v0 + (m-1)*vscale );
1867         Sigma->SetBinContent(bin, sqrt(s2)*vscale );
1868       }
1869       else{
1870         Mean->SetBinContent(bin, 0);
1871         Sigma->SetBinContent(bin, 0);
1872         Entr->SetBinContent(bin, 0);
1873       }
1874     }
1875   }
1876
1877   delete[] nBins;
1878   delete[] xMin;
1879   delete[] xMax;
1880   delete[] idx;
1881   delete[] hMap;
1882   delete[] hVal;
1883   delete[] hEntr;
1884   delete[] meanD;
1885   delete[] sigmaD;
1886
1887   if( !ok ){
1888     cout << "AliTPCcalibTracks: GetSparseStat : No memory or internal Error "<<endl;
1889     delete Mean;
1890     delete Sigma;
1891     delete Entr;
1892     Mean =0;
1893     Sigma = 0;
1894     Entr = 0;
1895     return -1;
1896   }
1897
1898   return 0;
1899 }
1900
1901
1902
1903 int AliTPCcalibTracks::CreateWaveCorrection(const  THnBase *DeltaY, THnBase *&MeanY, THnBase *&SigmaY, THnBase *&EntrY,
1904                                             Bool_t MirrorZ, Bool_t MirrorPad, Bool_t MirrorAngle, Int_t MinStat  )
1905 {
1906   // DeltaY is THnSparse:
1907   //
1908   // OBJ: TAxis 0    var     var
1909   // OBJ: TAxis 1    pad type        pad type
1910   // OBJ: TAxis 2    z       z
1911   // OBJ: TAxis 3    Qmax    Qmax
1912   // OBJ: TAxis 4    cog     cog
1913   // OBJ: TAxis 5    tan(angle)      tan(angle)
1914   // ...
1915
1916   MeanY = 0;
1917   SigmaY = 0;
1918   EntrY = 0;
1919
1920   int nDim = DeltaY->GetNdimensions();
1921   if( nDim<6 ){
1922     cout << "AliTPCcalibTracks: CreateWaveCorrection: Unknown input"<<endl;
1923     return -1;
1924   }
1925
1926   Int_t *nBins = new Int_t[nDim];
1927   Int_t *nBinsNew = new Int_t[nDim];
1928   Double_t *xMin = new Double_t[nDim];
1929   Double_t *xMax = new Double_t[nDim];
1930   Int_t *idx = new Int_t[nDim];
1931   THnSparseD *mergedDeltaY = 0;
1932
1933   int ret = 0;
1934
1935   if( !nBins || !nBinsNew || !xMin || !xMax || !idx ){
1936     ret = -1;
1937     cout << "AliTPCcalibTracks: CreateWaveCorrection: Out of memory"<<endl;
1938   }
1939
1940   if( ret==0 ){
1941
1942     for( int i=0; i<nDim; i++){
1943       xMin[i] = DeltaY->GetAxis(i)->GetXmin();
1944       xMax[i] = DeltaY->GetAxis(i)->GetXmax();
1945       nBins[i] = DeltaY->GetAxis(i)->GetNbins();
1946       nBinsNew[i] = nBins[i];
1947     }
1948
1949     // Merge cog axis
1950     if( MirrorPad ){
1951       Int_t centralBin = DeltaY->GetAxis(4)->FindFixBin(0.5);
1952       xMin[4] = DeltaY->GetAxis(4)->GetBinLowEdge(centralBin);
1953       nBinsNew[4] = nBinsNew[4]-centralBin+1;
1954     }
1955
1956     // Merge Z axis
1957     if( MirrorZ ){
1958       Int_t centralBin = DeltaY->GetAxis(2)->FindFixBin(0.0);
1959       xMin[2] = DeltaY->GetAxis(2)->GetBinLowEdge(centralBin)-0.0;
1960       nBinsNew[2] = nBinsNew[2]-centralBin+1;
1961     }
1962
1963     // Merge Angle axis
1964     if( MirrorAngle ){
1965       Int_t centralBin = DeltaY->GetAxis(5)->FindFixBin(0.0);
1966       xMin[5] = DeltaY->GetAxis(5)->GetBinLowEdge(centralBin)-0.0;
1967       nBinsNew[5] = nBinsNew[5]-centralBin+1;
1968     }
1969
1970     // Merge Sparse
1971
1972     mergedDeltaY = new THnSparseD("mergedDeltaY", "mergedDeltaY",nDim, nBinsNew, xMin, xMax );
1973
1974     if( !mergedDeltaY ){
1975       cout << "AliTPCcalibTracks: CreateWaveCorrection: Can not copy a Sparse"<<endl;
1976       ret = -1;
1977     }
1978   }
1979
1980   if( ret == 0 ){
1981
1982     for( int i=0; i<nDim; i++){
1983       const TAxis *ax = DeltaY->GetAxis(i);
1984       mergedDeltaY->GetAxis(i)->SetName(ax->GetName());
1985       mergedDeltaY->GetAxis(i)->SetTitle(ax->GetTitle());
1986       if( ax->GetXbins()->fN>0 ){
1987         mergedDeltaY->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
1988       }
1989     }
1990
1991     for( Long64_t binS=0; binS<DeltaY->GetNbins(); binS++){
1992       Double_t stat = DeltaY->GetBinContent(binS,idx);
1993       if( stat < 1 ) continue;
1994       bool swap0=0;
1995
1996       if( MirrorPad && idx[4]>0){ // underflow reserved for contains one-pad clusters, don't mirror
1997         Double_t v = DeltaY->GetAxis(4)->GetBinCenter(idx[4]);
1998         if( v < 0.5 ) swap0 = !swap0;
1999         idx[4] = mergedDeltaY->GetAxis(4)->FindFixBin( 0.5 + TMath::Abs(0.5 - v) );
2000       }
2001
2002       if( MirrorZ ){
2003         Double_t v = DeltaY->GetAxis(2)->GetBinCenter(idx[2]);
2004         if( v < 0.0 ) swap0 = !swap0;
2005         if( idx[2]<=0 ) idx[2] = nBinsNew[2]+1;
2006         else idx[2] = mergedDeltaY->GetAxis(2)->FindFixBin( TMath::Abs(v) );
2007       }
2008
2009       if( MirrorAngle ){
2010         Double_t v = DeltaY->GetAxis(5)->GetBinCenter(idx[5]);
2011         if( idx[5]<=0 ) idx[5] = nBinsNew[5]+1;
2012         else idx[5] = mergedDeltaY->GetAxis(5)->FindFixBin( TMath::Abs(v) );
2013       }
2014
2015       if( swap0 ){
2016         if( idx[0]<=0 ) idx[0] = nBinsNew[0]+1;
2017         else if( idx[0] >= nBins[0]+1 ) idx[0] = 0;
2018         else {
2019           Double_t v = DeltaY->GetAxis(0)->GetBinCenter(idx[0]);
2020           idx[0] = mergedDeltaY->GetAxis(0)->FindFixBin(-v);
2021         }
2022       }
2023
2024       Long64_t bin = mergedDeltaY->GetBin(idx,kTRUE);
2025       if( bin<0 ){
2026         cout << "AliTPCcalibTracks: CreateWaveCorrection : wrong bining"<<endl;
2027         continue;
2028       }
2029       mergedDeltaY->AddBinContent(bin,stat);
2030     }
2031
2032     ret = GetTHnStat(mergedDeltaY, MeanY, SigmaY, EntrY );
2033   }
2034
2035   if( ret==0 ){
2036
2037     MeanY->SetName("TPCWaveCorrectionMap");
2038     MeanY->SetTitle("TPCWaveCorrectionMap");
2039     SigmaY->SetName("TPCResolutionMap");
2040     SigmaY->SetTitle("TPCResolutionMap");
2041     EntrY->SetName("TPCWaveCorrectionEntr");
2042     EntrY->SetTitle("TPCWaveCorrectionEntr");
2043
2044     for( Long64_t bin=0; bin<MeanY->GetNbins(); bin++){
2045       Double_t stat = EntrY->GetBinContent(bin,idx);
2046
2047       // Normalize : Set no correction for one pad clusters
2048
2049       if( idx[3]<=0 ) MeanY->SetBinContent(bin,0);
2050
2051       // Suppress bins with low statistic
2052
2053       if( stat<MinStat ){
2054         EntrY->SetBinContent(bin,0);
2055         MeanY->SetBinContent(bin,0);
2056         SigmaY->SetBinContent(bin,-1);
2057       }
2058     }
2059
2060   }
2061
2062   delete[] nBins;
2063   delete[] nBinsNew;
2064   delete[] xMin;
2065   delete[] xMax;
2066   delete[] idx;
2067   delete mergedDeltaY;
2068
2069   if( ret!=0 ){
2070     delete MeanY;
2071     delete SigmaY;
2072     delete EntrY;
2073     MeanY = 0;
2074     SigmaY = 0;
2075     EntrY = 0;
2076   }
2077
2078   return ret;
2079 }
2080
2081 int AliTPCcalibTracks::UpdateClusterParam( AliTPCClusterParam *cParam, Bool_t MirrorZ, Bool_t MirrorPad, Bool_t MirrorAngle, Int_t MinStat )
2082 {
2083   if( !cParam || !fHisDeltaY) return -1;
2084
2085   THnBase *meanY  = 0;
2086   THnBase *sigmaY = 0;
2087   THnBase *entrY  = 0;
2088   int ret = CreateWaveCorrection(fHisDeltaY, meanY, sigmaY, entrY,  MirrorZ, MirrorAngle, MirrorPad, MinStat  );
2089   if( ret<0 ) return ret;
2090   cParam->SetWaveCorrectionMap(meanY);
2091   cParam->SetResolutionYMap(sigmaY);
2092   delete meanY;
2093   delete sigmaY;
2094   delete entrY;
2095   return 0;
2096 }