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