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