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 to analyse tracks for calibration                               //
  20 //     to be used as a component in AliTPCSelectorTracks                     //
  21 //     In the constructor you have to specify name and title                 //
  22 //     to get the Object out of a file.                                      //
  23 //     The parameter 'clusterParam', a AliTPCClusterParam object             //
  24 //      (needed for TPC cluster error and shape parameterization)            //
  25 //     Normally you get this object out of the file 'TPCClusterParam.root'   //
  26 //     In the parameter 'cuts' the cuts are specified, that decide           //
  27 //     weather a track will be accepted for calibration or not.              //
  28 //                                                                           //
  29 //
  30 //     The data flow:
  31 //
  32 /*
  33    Raw Data -> Local Reconstruction -> Tracking ->     Calibration -> RefData (component itself)
  34                Offline/HLT             Offline/HLT                    OCDB entries (AliTPCClusterParam)
  35 */
  36
  37 /*
  38
  39 How to retrive it from file (created using calibration task):
  40
  41 gSystem->Load("libANALYSIS");
  42 gSystem->Load("libTPCcalib");
  43 TFile fcalib("CalibObjects.root");
  44 TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
  45 AliTPCcalibTracks * calibTracks = ( AliTPCcalibTracks *)array->FindObject("calibTracks");
  46
  47
  48 //USAGE of debug stream example
  49  gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  50   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  51   AliXRDPROOFtoolkit tool;
  52   TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
  53   chainres->Lookup();
  54 */
  55
  56
  57                                                                //
  58 ///////////////////////////////////////////////////////////////////////////////
  59
  60 //
  61 // ROOT includes
  62 //
  63 #include <iostream>
  64 #include <fstream>
  65 using namespace std;
  66
  67 #include <TROOT.h>
  68 #include <TChain.h>
  69 #include <TFile.h>
  70 #include <TH3F.h>
  71 #include <TProfile.h>
  72
  73 //
  74 //#include <TPDGCode.h>
  75 #include <TStyle.h>
  76 #include "TLinearFitter.h"
  77 //#include "TMatrixD.h"
  78 #include "TTreeStream.h"
  79 #include "TF1.h"
  80 #include <TCanvas.h>
  81 #include <TGraph2DErrors.h>
  82 #include "TPostScript.h"
  83 #include "TCint.h"
  84
  85 #include <TH2D.h>
  86 #include <TF2.h>
  87 #include <TSystem.h>
  88 #include <TCollection.h>
  89 #include <iostream>
  90 #include <TLinearFitter.h>
  91 #include <TString.h>
  92
  93 //
  94 // AliROOT includes
  95 //
  96 #include "AliMagF.h"
  97 #include "AliTracker.h"
  98 #include "AliESD.h"
  99 #include "AliESDtrack.h"
 100 #include "AliESDfriend.h"
 101 #include "AliESDfriendTrack.h"
 102 #include "AliTPCseed.h"
 103 #include "AliTPCclusterMI.h"
 104 #include "AliTPCROC.h"
 105
 106 #include "AliTPCParamSR.h"
 107 #include "AliTrackPointArray.h"
 108 #include "AliTPCcalibTracks.h"
 109 #include "AliTPCClusterParam.h"
 110 #include "AliTPCcalibTracksCuts.h"
 111 #include "AliTPCCalPadRegion.h"
 112 #include "AliTPCCalPad.h"
 113 #include "AliTPCCalROC.h"
 114 #include "TText.h"
 115 #include "TPaveText.h"
 116 #include "TSystem.h"
 117 #include "TStatToolkit.h"
 118 #include "TCut.h"
 119 #include "THnSparse.h"
 120
 121
 122
 123 ClassImp(AliTPCcalibTracks)
 124
 125
 126 AliTPCcalibTracks::AliTPCcalibTracks():
 127   AliTPCcalibBase(),
 128   fClusterParam(0),
 129   fROC(0),
 130   fHisDeltaY(0),    // THnSparse - delta Y
 131   fHisDeltaZ(0),    // THnSparse - delta Z
 132   fHisRMSY(0),      // THnSparse - rms Y
 133   fHisRMSZ(0),      // THnSparse - rms Z
 134   fHisQmax(0),      // THnSparse - qmax
 135   fHisQtot(0),      // THnSparse - qtot
 136   fArrayAmpRow(0),
 137   fArrayAmp(0),
 138   fArrayQDY(0),
 139   fArrayQDZ(0),
 140   fArrayQRMSY(0),
 141   fArrayQRMSZ(0),
 142   fArrayChargeVsDriftlength(0),
 143   fcalPadRegionChargeVsDriftlength(0),
 144   fDeltaY(0),
 145   fDeltaZ(0),
 146   fResolY(0),
 147   fResolZ(0),
 148   fRMSY(0),
 149   fRMSZ(0),
 150   fCuts(0),
 151   fHclus(0),
 152   fRejectedTracksHisto(0),
 153   fHclusterPerPadrow(0),
 154   fHclusterPerPadrowRaw(0),
 155   fClusterCutHisto(0),
 156   fCalPadClusterPerPad(0),
 157   fCalPadClusterPerPadRaw(0)
 158 {
 159    //
 160    // AliTPCcalibTracks default constructor
 161    //
 162   if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' default constructor called" << endl;
 163 }
 164
 165
 166
 167 AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks):
 168   AliTPCcalibBase(calibTracks),
 169   fClusterParam(0),
 170   fROC(0),
 171   fHisDeltaY(0),    // THnSparse - delta Y
 172   fHisDeltaZ(0),    // THnSparse - delta Z
 173   fHisRMSY(0),      // THnSparse - rms Y
 174   fHisRMSZ(0),      // THnSparse - rms Z
 175   fHisQmax(0),      // THnSparse - qmax
 176   fHisQtot(0),      // THnSparse - qtot
 177   fArrayAmpRow(0),
 178   fArrayAmp(0),
 179   fArrayQDY(0),
 180   fArrayQDZ(0),
 181   fArrayQRMSY(0),
 182   fArrayQRMSZ(0),
 183   fArrayChargeVsDriftlength(0),
 184   fcalPadRegionChargeVsDriftlength(0),
 185   fDeltaY(0),
 186   fDeltaZ(0),
 187   fResolY(0),
 188   fResolZ(0),
 189   fRMSY(0),
 190   fRMSZ(0),
 191   fCuts(0),
 192   fHclus(0),
 193   fRejectedTracksHisto(0),
 194   fHclusterPerPadrow(0),
 195   fHclusterPerPadrowRaw(0),
 196   fClusterCutHisto(0),
 197   fCalPadClusterPerPad(0),
 198   fCalPadClusterPerPadRaw(0)
 199 {
 200    //
 201    // AliTPCcalibTracks copy constructor
 202    //
 203   if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' copy constructor ***** " << endl;
 204
 205    Bool_t dirStatus = TH1::AddDirectoryStatus();
 206    TH1::AddDirectory(kFALSE);
 207
 208    Int_t length = -1;
 209    // backward compatibility: if the data member doesn't yet exist, it will not be merged
 210    (calibTracks.fArrayAmpRow) ? length = calibTracks.fArrayAmpRow->GetEntriesFast() : length = -1;
 211    fArrayAmpRow = new TObjArray(length);
 212    fArrayAmp = new TObjArray(length);
 213    for (Int_t i = 0; i < length; i++) {
 214       fArrayAmpRow->AddAt( (TProfile*)calibTracks.fArrayAmpRow->At(i)->Clone(), i);
 215       fArrayAmp->AddAt( ((TProfile*)calibTracks.fArrayAmp->At(i)->Clone()), i);
 216    }
 217
 218    (calibTracks.fArrayQDY) ? length = calibTracks.fArrayQDY->GetEntriesFast() : length = -1;
 219    fArrayQDY= new TObjArray(length);
 220    fArrayQDZ= new TObjArray(length);
 221    fArrayQRMSY= new TObjArray(length);
 222    fArrayQRMSZ= new TObjArray(length);
 223    for (Int_t i = 0; i < length; i++) {
 224       fArrayQDY->AddAt( ((TH1F*)calibTracks.fArrayQDY->At(i)->Clone()), i);
 225       fArrayQDZ->AddAt( ((TH1F*)calibTracks.fArrayQDZ->At(i)->Clone()), i);
 226       fArrayQRMSY->AddAt( ((TH1F*)calibTracks.fArrayQRMSY->At(i)->Clone()), i);
 227       fArrayQRMSZ->AddAt( ((TH1F*)calibTracks.fArrayQRMSZ->At(i)->Clone()), i);
 228    }
 229
 230    (calibTracks.fResolY) ? length = calibTracks.fResolY->GetEntriesFast() : length = -1;
 231    fResolY = new TObjArray(length);
 232    fResolZ = new TObjArray(length);
 233    fRMSY = new TObjArray(length);
 234    fRMSZ = new TObjArray(length);
 235    for (Int_t i = 0; i < length; i++) {
 236       fResolY->AddAt( ((TH1F*)calibTracks.fResolY->At(i)->Clone()), i);
 237       fResolZ->AddAt( ((TH1F*)calibTracks.fResolZ->At(i)->Clone()), i);
 238       fRMSY->AddAt( ((TH1F*)calibTracks.fRMSY->At(i)->Clone()), i);
 239       fRMSZ->AddAt( ((TH1F*)calibTracks.fRMSZ->At(i)->Clone()), i);
 240    }
 241
 242    (calibTracks.fArrayChargeVsDriftlength) ? length = calibTracks.fArrayChargeVsDriftlength->GetEntriesFast() : length = -1;
 243    (calibTracks.fArrayChargeVsDriftlength) ? fArrayChargeVsDriftlength = new TObjArray(length) : fArrayChargeVsDriftlength = 0;
 244    for (Int_t i = 0; i < length; i++) {
 245       fArrayChargeVsDriftlength->AddAt( ((TProfile*)calibTracks.fArrayChargeVsDriftlength->At(i)->Clone()), i);
 246    }
 247
 248    fDeltaY =  (TH1F*)calibTracks.fDeltaY->Clone();
 249    fDeltaZ =  (TH1F*)calibTracks.fDeltaZ->Clone();
 250    fHclus = (TH1I*)calibTracks.fHclus->Clone();
 251    fClusterCutHisto = (TH2I*)calibTracks.fClusterCutHisto->Clone();
 252    fRejectedTracksHisto    = (TH1I*)calibTracks.fRejectedTracksHisto->Clone();
 253    fHclusterPerPadrow      = (TH1I*)calibTracks.fHclusterPerPadrow->Clone();
 254    fHclusterPerPadrowRaw   = (TH1I*)calibTracks.fHclusterPerPadrowRaw->Clone();
 255    fcalPadRegionChargeVsDriftlength = (AliTPCCalPadRegion*)calibTracks.fcalPadRegionChargeVsDriftlength->Clone();
 256    fCalPadClusterPerPad    = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPad->Clone();
 257    fCalPadClusterPerPadRaw = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPadRaw->Clone();
 258
 259    fCuts = new AliTPCcalibTracksCuts(calibTracks.fCuts->GetMinClusters(), calibTracks.fCuts->GetMinRatio(),
 260       calibTracks.fCuts->GetMax1pt(), calibTracks.fCuts->GetEdgeYXCutNoise(), calibTracks.fCuts->GetEdgeThetaCutNoise());
 261    SetNameTitle(calibTracks.GetName(), calibTracks.GetTitle());
 262    TH1::AddDirectory(dirStatus); // set status back to original status
 263 //    cout << "+++++ end of copy constructor +++++" << endl;   // TO BE REMOVED
 264 }
 265
 266
 267 AliTPCcalibTracks & AliTPCcalibTracks::operator=(const AliTPCcalibTracks& calibTracks){
 268   //
 269   // assgnment operator
 270   //
 271   if (this != &calibTracks) {
 272     new (this) AliTPCcalibTracks(calibTracks);
 273   }
 274   return *this;
 275
 276 }
 277
 278
 279 AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam,  AliTPCcalibTracksCuts* cuts, Int_t logLevel) :
 280   AliTPCcalibBase(),
 281   fClusterParam(0),
 282   fROC(0),
 283   fHisDeltaY(0),    // THnSparse - delta Y
 284   fHisDeltaZ(0),    // THnSparse - delta Z
 285   fHisRMSY(0),      // THnSparse - rms Y
 286   fHisRMSZ(0),      // THnSparse - rms Z
 287   fHisQmax(0),      // THnSparse - qmax
 288   fHisQtot(0),      // THnSparse - qtot
 289   fArrayAmpRow(0),
 290   fArrayAmp(0),
 291   fArrayQDY(0),
 292   fArrayQDZ(0),
 293   fArrayQRMSY(0),
 294   fArrayQRMSZ(0),
 295   fArrayChargeVsDriftlength(0),
 296   fcalPadRegionChargeVsDriftlength(0),
 297   fDeltaY(0),
 298   fDeltaZ(0),
 299   fResolY(0),
 300   fResolZ(0),
 301   fRMSY(0),
 302   fRMSZ(0),
 303   fCuts(0),
 304   fHclus(0),
 305   fRejectedTracksHisto(0),
 306   fHclusterPerPadrow(0),
 307   fHclusterPerPadrowRaw(0),
 308   fClusterCutHisto(0),
 309   fCalPadClusterPerPad(0),
 310   fCalPadClusterPerPadRaw(0)
 311  {
 312    //
 313    // AliTPCcalibTracks constructor
 314    // specify 'name' and 'title' of your object
 315    // specify 'clusterParam', (needed for TPC cluster error and shape parameterization)
 316    // In the parameter 'cuts' the cuts are specified, that decide
 317    // weather a track will be accepted for calibration or not.
 318    //
 319    // fDebugLevel - debug output: -1: silence, 0: default, 1: things like constructor called, 5: write fDebugStreamer, 6: waste your screen
 320    //
 321    // All histograms are instatiated in this constructor.
 322    //
 323    this->SetName(name);
 324    this->SetTitle(title);
 325
 326    if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' main constructor ***** " << endl;
 327    G__SetCatchException(0);
 328
 329    fClusterParam = clusterParam;
 330    if (fClusterParam){
 331      fClusterParam->SetInstance(fClusterParam);
 332    }
 333    else {
 334      Error("AliTPCcalibTracks","No cluster parametrization found! A valid clusterParam object is needed in the constructor. (To be found in 'TPCClusterParam.root'.)");
 335    }
 336    fCuts = cuts;
 337    SetDebugLevel(logLevel);
 338    MakeHistos();
 339
 340    TH1::AddDirectory(kFALSE);
 341
 342    char chname[1000];
 343    TProfile * prof1=0;
 344    TH1F     * his1 =0;
 345    fHclus = new TH1I("hclus","Number of clusters per track",160, 0, 160);     // valgrind 3
 346    fRejectedTracksHisto    = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 100, -1, 10);
 347    fHclusterPerPadrow      = new TH1I("fHclusterPerPadrow", " clusters per padRow, used for the resolution tree", 160, 0, 160);
 348    fHclusterPerPadrowRaw   = new TH1I("fHclusterPerPadrowRaw", " clusters per padRow, before cutting clusters", 160, 0, 160);
 349    fCalPadClusterPerPad    = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad");
 350    fCalPadClusterPerPadRaw = new AliTPCCalPad("fCalPadClusterPerPadRaw", "clusters per pad, before cutting clusters");
 351    fClusterCutHisto = new TH2I("fClusterCutHisto", "Cutted cluster over padRow; Cut Criterium; PadRow", 5,1,5, 160,0,159);
 352
 353    // Amplitude  - sector - row histograms
 354    fArrayAmpRow = new TObjArray(72);
 355    fArrayAmp    = new TObjArray(72);
 356    fArrayChargeVsDriftlength = new TObjArray(72);
 357
 358    for (Int_t i = 0; i < 36; i++){
 359       sprintf(chname,"Amp_row_Sector%d",i);
 360       prof1 = new TProfile(chname,chname,63,0,64);          // valgrind 3   193,536 bytes in 354 blocks are still reachable
 361       prof1->SetXTitle("Pad row");
 362       prof1->SetYTitle("Mean Max amplitude");
 363       fArrayAmpRow->AddAt(prof1,i);
 364       sprintf(chname,"Amp_row_Sector%d",i+36);
 365       prof1 = new TProfile(chname,chname,96,0,97);       // valgrind 3   3,912 bytes in 6 blocks are possibly lost
 366       prof1->SetXTitle("Pad row");
 367       prof1->SetYTitle("Mean Max  amplitude");
 368       fArrayAmpRow->AddAt(prof1,i+36);
 369
 370       // amplitude
 371       sprintf(chname,"Amp_Sector%d",i);
 372       his1 = new TH1F(chname,chname,100,0,32);         // valgrind
 373       his1->SetXTitle("Max Amplitude (ADC)");
 374       fArrayAmp->AddAt(his1,i);
 375       sprintf(chname,"Amp_Sector%d",i+36);
 376       his1 = new TH1F(chname,chname,100,0,32);         // valgrind 3   13,408,208 bytes in 229 blocks are still reachable
 377       his1->SetXTitle("Max Amplitude (ADC)");
 378       fArrayAmp->AddAt(his1,i+36);
 379
 380       // driftlength
 381       sprintf(chname, "driftlengt vs. charge, ROC %i", i);
 382       prof1 = new TProfile(chname, chname, 25, 0, 250);
 383       prof1->SetYTitle("Charge");
 384       prof1->SetXTitle("Driftlength");
 385       fArrayChargeVsDriftlength->AddAt(prof1,i);
 386       sprintf(chname, "driftlengt vs. charge, ROC %i", i+36);
 387       prof1 = new TProfile(chname, chname, 25, 0, 250);
 388       prof1->SetYTitle("Charge");
 389       prof1->SetXTitle("Driftlength");
 390       fArrayChargeVsDriftlength->AddAt(prof1,i+36);
 391    }
 392
 393    TH1::AddDirectory(kFALSE);
 394
 395    fDeltaY = new TH1F("DeltaY","DeltaY",100,-1,1);
 396    fDeltaZ = new TH1F("DeltaZ","DeltaZ",100,-1,1);
 397
 398    fResolY = new TObjArray(3);
 399    fResolZ = new TObjArray(3);
 400    fRMSY   = new TObjArray(3);
 401    fRMSZ   = new TObjArray(3);
 402    TH3F * his3D;
 403    //
 404    his3D = new TH3F("Resol Y0","Resol Y0", 5,20,250, 4, 0,1., 50, -1,1);
 405    fResolY->AddAt(his3D,0);
 406    his3D = new TH3F("Resol Y1","Resol Y1", 5,20,250, 4, 0,1., 50, -1,1);
 407    fResolY->AddAt(his3D,1);
 408    his3D = new TH3F("Resol Y2","Resol Y2", 5,20,250, 4, 0,0.8, 50, -1,1);
 409    fResolY->AddAt(his3D,2);
 410    //
 411    his3D = new TH3F("Resol Z0","Resol Z0", 5,20,250, 4, 0,1, 50, -1,1);
 412    fResolZ->AddAt(his3D,0);
 413    his3D = new TH3F("Resol Z1","Resol Z1", 5,20,250, 4, 0,1, 50, -1,1);
 414    fResolZ->AddAt(his3D,1);
 415    his3D = new TH3F("Resol Z2","Resol Z2", 5,20,250, 4, 0,1, 50, -1,1);
 416    fResolZ->AddAt(his3D,2);
 417    //
 418    his3D = new TH3F("RMS Y0","RMS Y0", 5,20,250, 4, 0,1., 50, 0,0.8);
 419    fRMSY->AddAt(his3D,0);
 420    his3D = new TH3F("RMS Y1","RMS Y1", 5,20,250, 4, 0,1., 50, 0,0.8);
 421    fRMSY->AddAt(his3D,1);
 422    his3D = new TH3F("RMS Y2","RMS Y2", 5,20,250, 4, 0,0.8, 50, 0,0.8);
 423    fRMSY->AddAt(his3D,2);
 424    //
 425    his3D = new TH3F("RMS Z0","RMS Z0", 5,20,250, 4, 0,1, 50, 0,0.8);
 426    fRMSZ->AddAt(his3D,0);
 427    his3D = new TH3F("RMS Z1","RMS Z1", 5,20,250, 4, 0,1, 50, 0,0.8);
 428    fRMSZ->AddAt(his3D,1);
 429    his3D = new TH3F("RMS Z2","RMS Z2", 5,20,250, 4, 0,1, 50, 0,0.8);
 430    fRMSZ->AddAt(his3D,2);
 431    //
 432
 433    TH1::AddDirectory(kFALSE);
 434
 435    fArrayQDY = new TObjArray(300);
 436    fArrayQDZ = new TObjArray(300);
 437    fArrayQRMSY = new TObjArray(300);
 438    fArrayQRMSZ = new TObjArray(300);
 439    for (Int_t iq = 0; iq <= 10; iq++){
 440       for (Int_t ipad = 0; ipad < 3; ipad++){
 441          Int_t   bin   = GetBin(iq, ipad);
 442          Float_t qmean = GetQ(bin);
 443          char hname[200];
 444          sprintf(hname,"ResolY Pad%d Qmiddle%f",ipad, qmean);
 445          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
 446          fArrayQDY->AddAt(his3D, bin);
 447          sprintf(hname,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
 448          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
 449          fArrayQDZ->AddAt(his3D, bin);
 450          sprintf(hname,"RMSY Pad%d Qmiddle%f",ipad, qmean);
 451          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
 452          fArrayQRMSY->AddAt(his3D, bin);
 453          sprintf(hname,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
 454          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
 455          fArrayQRMSZ->AddAt(his3D, bin);
 456       }
 457    }
 458
 459    fcalPadRegionChargeVsDriftlength = new AliTPCCalPadRegion("fcalPadRegionChargeVsDriftlength", "TProfiles with charge vs driftlength for each pad region");
 460    TProfile *tempProf;
 461    for (UInt_t padSize = 0; padSize < 3; padSize++) {
 462       for (UInt_t isector = 0; isector < 36; isector++) {
 463          if (padSize == 0) sprintf(chname, "driftlengt vs. charge, sector %i, short pads", isector);
 464          if (padSize == 1) sprintf(chname, "driftlengt vs. charge, sector %i, medium  pads", isector);
 465          if (padSize == 2) sprintf(chname, "driftlengt vs. charge, sector %i, long  pads", isector);
 466          tempProf = new TProfile(chname, chname, 500, 0, 250);
 467          tempProf->SetYTitle("Charge");
 468          tempProf->SetXTitle("Driftlength");
 469          fcalPadRegionChargeVsDriftlength->SetObject(tempProf, isector, padSize);
 470       }
 471    }
 472
 473
 474    if (GetDebugLevel() > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl;
 475    cout << "end of main constructor" << endl; // TO BE REMOVED
 476 }
 477
 478
 479 AliTPCcalibTracks::~AliTPCcalibTracks() {
 480    //
 481    // AliTPCcalibTracks destructor
 482    //
 483
 484   if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' destuctor called." << endl;
 485    Int_t length = 0;
 486    if (fArrayAmpRow) length = fArrayAmpRow->GetEntriesFast();
 487    for (Int_t i = 0; i < length; i++){
 488       delete fArrayAmpRow->At(i);
 489       delete fArrayAmp->At(i);
 490    }
 491    delete fArrayAmpRow;
 492    delete fArrayAmp;
 493
 494    delete fDeltaY;
 495    delete fDeltaZ;
 496
 497    if (fResolY) length = fResolY->GetEntriesFast();
 498    for (Int_t i = 0; i < length; i++){
 499       delete fResolY->At(i);
 500       delete fResolZ->At(i);
 501       delete fRMSY->At(i);
 502       delete fRMSZ->At(i);
 503    }
 504    delete fResolY;
 505    delete fResolZ;
 506    delete fRMSY;
 507    delete fRMSZ;
 508
 509    if (fArrayQDY) length = fArrayQDY->GetEntriesFast();
 510    for (Int_t i = 0; i < length; i++){
 511       delete fArrayQDY->At(i);
 512       delete fArrayQDZ->At(i);
 513       delete fArrayQRMSY->At(i);
 514       delete fArrayQRMSZ->At(i);
 515    }
 516
 517    if (fArrayChargeVsDriftlength) length = fArrayChargeVsDriftlength->GetEntriesFast();
 518    for (Int_t i = 0; i < length; i++){
 519       delete fArrayChargeVsDriftlength->At(i);
 520    }
 521
 522
 523    delete fArrayQDY;
 524    delete fArrayQDZ;
 525    delete fArrayQRMSY;
 526    delete fArrayQRMSZ;
 527    delete fArrayChargeVsDriftlength;
 528
 529   delete fHclus;
 530   delete fRejectedTracksHisto;
 531   delete fClusterCutHisto;
 532   delete fHclusterPerPadrow;
 533   delete fHclusterPerPadrowRaw;
 534   if (fCalPadClusterPerPad)    delete fCalPadClusterPerPad;
 535   if (fCalPadClusterPerPadRaw) delete fCalPadClusterPerPadRaw;
 536   if(fcalPadRegionChargeVsDriftlength) {
 537      fcalPadRegionChargeVsDriftlength->Delete();
 538      delete fcalPadRegionChargeVsDriftlength;
 539   }
 540   delete fHisDeltaY;    // THnSparse - delta Y
 541   delete fHisDeltaZ;    // THnSparse - delta Z
 542   delete fHisRMSY;      // THnSparse - rms Y
 543   delete fHisRMSZ;      // THnSparse - rms Z
 544   delete fHisQmax;      // THnSparse - qmax
 545   delete fHisQtot;      // THnSparse - qtot
 546
 547 }
 548
 549
 550
 551 void AliTPCcalibTracks::Process(AliTPCseed *track){
 552    //
 553    // To be called in the selector
 554    // first AcceptTrack is evaluated, then calls all the following analyse functions:
 555    // FillResolutionHistoLocal(track)
 556
 557    //
 558   if (GetDebugLevel() > 5) Info("Process","Starting to process the track...");
 559    Int_t accpetStatus = AcceptTrack(track);
 560    if (accpetStatus == 0) {
 561       FillResolutionHistoLocal(track);
 562    }
 563    else fRejectedTracksHisto->Fill(accpetStatus);
 564 }
 565
 566
 567
 568 Int_t AliTPCcalibTracks::GetBin(Float_t q, Int_t pad){
 569   //
 570   // calculate bins for given q and pad type
 571   // used in TObjArray
 572   //
 573   Int_t res = TMath::Max( TMath::Nint((TMath::Sqrt(q) - 3.)), 0 );
 574   res *= 3;
 575   res += pad;
 576   return res;
 577 }
 578
 579
 580 Int_t AliTPCcalibTracks::GetBin(Int_t iq, Int_t pad){
 581   //
 582   // calculate bins for given iq and pad type
 583   // used in TObjArray
 584   //
 585   return iq * 3 + pad;;
 586 }
 587
 588
 589 Float_t AliTPCcalibTracks::GetQ(Int_t bin){
 590    //
 591    // returns to bin belonging charge
 592    // (bin / 3 + 3)^2
 593    //
 594    Int_t bin0 = bin / 3;
 595    bin0 += 3;
 596    return bin0 * bin0;
 597 }
 598
 599
 600 Float_t AliTPCcalibTracks::GetPad(Int_t bin){
 601    //
 602    // returns to bin belonging pad
 603    // bin % 3
 604    //
 605    return bin % 3;
 606 }
 607
 608
 609
 610 Int_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){
 611   //
 612   // Function, that decides wheather a given track is accepted for
 613   // the analysis or not.
 614   // The cuts are specified in the AliTPCcalibTracksCuts object 'fCuts'
 615   // Returns 0 if a track is accepted or an integer different from 0
 616   // to indicate the failed cut
 617   //
 618   const Int_t   kMinClusters  = fCuts->GetMinClusters();
 619   const Float_t kMinRatio     = fCuts->GetMinRatio();
 620   const Float_t kMax1pt       = fCuts->GetMax1pt();
 621   const Float_t kEdgeYXCutNoise    = fCuts->GetEdgeYXCutNoise();
 622   const Float_t kEdgeThetaCutNoise = fCuts->GetEdgeThetaCutNoise();
 623
 624   //
 625   // edge induced noise tracks - NEXT RELEASE will be removed during tracking
 626   if ( TMath::Abs(track->GetY() / track->GetX()) > kEdgeYXCutNoise )
 627     if ( TMath::Abs(track->GetTgl()) < kEdgeThetaCutNoise ) return 1;
 628   if (track->GetNumberOfClusters() < kMinClusters) return 2;
 629   Float_t ratio = track->GetNumberOfClusters() / (track->GetNFoundable() + 1.);
 630   if (ratio < kMinRatio) return 3;
 631 //   Float_t mpt = track->Get1Pt();       // Get1Pt() doesn't exist any more
 632   Float_t mpt = track->GetSigned1Pt();
 633   if (TMath::Abs(mpt) > kMax1pt) return 4;
 634   //if (TMath::Abs(track->GetZ())>240.) return kFALSE;
 635   //if (TMath::Abs(track->GetZ())<10.) return kFALSE;
 636   //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
 637
 638   if (GetDebugLevel() > 20) Info("AcceptTrack","Track has been accepted.");
 639   return 0;
 640 }
 641
 642
 643 void  AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
 644    //
 645    // fill resolution histograms - localy - tracklet in the neighborhood
 646    // write debug information to 'TPCSelectorDebug.root'
 647    //
 648    // _ the main function, called during track analysis _
 649    //
 650    // loop over all padrows along the track
 651    // fit tracklets (length: 13 clusters) calculate mean chi^2 for this track-fit in Y and Z direction
 652    //
 653    // loop again over all padrows along the track
 654    // fit tracklet (clusters in given padrow +- kDelta padrows)
 655    // with polynom of 2nd order and two polynoms of 1st order
 656    // take both polynoms of 1st order, calculate difference of their parameters
 657    // add covariance matrixes and calculate chi2 of this difference
 658    // if this chi2 is bigger than a given threshold, assume that the current cluster is
 659    // a kink an goto next padrow
 660    // if not:
 661    // fill fArrayAmpRow, array with amplitudes vs. row for given sector
 662    // fill fArrayAmp, array with amplitude histograms for give sector
 663    // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fDeltaY, fDeltaZ, fResolY, fResolZ, fArrayQDY, fArrayQDY
 664    //
 665    // write debug information to 'TPCSelectorDebug.root'
 666    // only for every kDeltaWriteDebugStream'th padrow to reduce data volume
 667    // and to avoid redundant data
 668    //
 669
 670   static TLinearFitter fFitterLinY1(2,"pol1");   //
 671   static TLinearFitter fFitterLinZ1(2,"pol1");   //
 672   static TLinearFitter fFitterLinY2(2,"pol1");   //
 673   static TLinearFitter fFitterLinZ2(2,"pol1");   //
 674   static TLinearFitter fFitterParY(3,"pol2");    //
 675   static TLinearFitter fFitterParZ(3,"pol2");    //
 676
 677   fFitterLinY1.StoreData(kFALSE);
 678   fFitterLinZ1.StoreData(kFALSE);
 679   fFitterLinY2.StoreData(kFALSE);
 680   fFitterLinZ2.StoreData(kFALSE);
 681   fFitterParY.StoreData(kFALSE);
 682   fFitterParZ.StoreData(kFALSE);
 683
 684
 685   if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****");
 686    const Int_t   kDelta    = 10;          // delta rows to fit
 687    const Float_t kMinRatio = 0.75;        // minimal ratio
 688    //   const Float_t kCutChi2  = 6.;          // cut chi2 - left right  - kink removal
 689    const Float_t kErrorFraction = 0.5;    // use only clusters with small interpolation error - for error param
 690    const Int_t   kFirstLargePad = 127;    // medium pads -> long pads
 691    const Float_t kLargePadSize  = 1.5;    // factor between medium and long pads' area
 692    const Int_t   kDeltaWriteDebugStream  = 5;  // only for every kDeltaWriteDebugStream'th padrow debug information is calulated and written to debugstream
 693    TVectorD paramY0(2);
 694    TVectorD paramZ0(2);
 695    TVectorD paramY1(2);
 696    TVectorD paramZ1(2);
 697    TVectorD paramY2(3);
 698    TVectorD paramZ2(3);
 699    TMatrixD matrixY0(2,2);
 700    TMatrixD matrixZ0(2,2);
 701    TMatrixD matrixY1(2,2);
 702    TMatrixD matrixZ1(2,2);
 703
 704    // estimate mean error
 705    Int_t nTrackletsAll = 0;       // number of tracklets for given track
 706    Float_t csigmaY     = 0;       // mean sigma for tracklet refit in Y direction
 707    Float_t csigmaZ     = 0;       // mean sigma for tracklet refit in Z direction
 708    Int_t nClusters     = 0;       // working variable, number of clusters per tracklet
 709    Int_t sectorG       = -1;      // working variable, sector of tracklet, has to stay constant for one tracklet
 710
 711    fHclus->Fill(track->GetNumberOfClusters());      // for statistics overview
 712    // ---------------------------------------------------------------------
 713    for (Int_t irow = 0; irow < 159; irow++){
 714       // loop over all rows along the track
 715       // fit tracklets (length: 13 rows) with pol2 in Y and Z direction
 716       // calculate mean chi^2 for this track-fit in Y and Z direction
 717       AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
 718       if (!cluster0) continue;  // no cluster found
 719       Int_t sector = cluster0->GetDetector();
 720       fHclusterPerPadrowRaw->Fill(irow);
 721
 722       Int_t ipad = TMath::Nint(cluster0->GetPad());
 723       Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
 724       fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
 725
 726       if (sector != sectorG){
 727          // track leaves sector before it crossed enough rows to fit / initialization
 728          nClusters = 0;
 729          fFitterParY.ClearPoints();
 730          fFitterParZ.ClearPoints();
 731          sectorG = sector;
 732       }
 733       else {
 734          nClusters++;
 735          Double_t x = cluster0->GetX();
 736          fFitterParY.AddPoint(&x, cluster0->GetY(), 1);
 737          fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1);
 738          //
 739          if ( nClusters >= kDelta + 3 ){
 740          // if more than 13 (kDelta+3) clusters were added to the fitters
 741          // fit the tracklet, increase trackletCounter
 742          fFitterParY.Eval();
 743          fFitterParZ.Eval();
 744          nTrackletsAll++;
 745          csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.);
 746          csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.);
 747          nClusters = -1;
 748          fFitterParY.ClearPoints();
 749          fFitterParZ.ClearPoints();
 750          }
 751       }
 752    }      // for (Int_t irow = 0; irow < 159; irow++)
 753    // mean chi^2 for all tracklet fits in Y and in Z direction:
 754    csigmaY = TMath::Sqrt(csigmaY / nTrackletsAll);
 755    csigmaZ = TMath::Sqrt(csigmaZ / nTrackletsAll);
 756    // ---------------------------------------------------------------------
 757
 758    for (Int_t irow = 0; irow < 159; irow++){
 759       // loop again over all rows along the track
 760       // do analysis
 761       //
 762       Int_t nclFound = 0;  // number of clusters in the neighborhood
 763       Int_t ncl0 = 0;      // number of clusters in rows < rowOfCenterCluster
 764       Int_t ncl1 = 0;      // number of clusters in rows > rowOfCenterCluster
 765       AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
 766       if (!cluster0) continue;
 767       Int_t sector = cluster0->GetDetector();
 768       Float_t xref = cluster0->GetX();
 769
 770       // Make Fit
 771       fFitterParY.ClearPoints();
 772       fFitterParZ.ClearPoints();
 773       fFitterLinY1.ClearPoints();
 774       fFitterLinZ1.ClearPoints();
 775       fFitterLinY2.ClearPoints();
 776       fFitterLinZ2.ClearPoints();
 777
 778       // fit tracklet (clusters in given padrow +- kDelta padrows)
 779       // with polynom of 2nd order and two polynoms of 1st order
 780       // take both polynoms of 1st order, calculate difference of their parameters
 781       // add covariance matrixes and calculate chi2 of this difference
 782       // if this chi2 is bigger than a given threshold, assume that the current cluster is
 783       // a kink an goto next padrow
 784
 785       for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
 786          // loop over irow +- kDelta rows (neighboured rows)
 787          //
 788          //
 789          if (idelta == 0) continue;                                // don't use center cluster
 790          if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
 791          AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
 792          if (!currentCluster) continue;
 793          if (currentCluster->GetType() < 0) continue;
 794          if (currentCluster->GetDetector() != sector) continue;
 795          Double_t x = currentCluster->GetX() - xref;  // x = differece: current cluster - cluster @ irow
 796          nclFound++;
 797          if (idelta < 0){
 798          ncl0++;
 799          fFitterLinY1.AddPoint(&x, currentCluster->GetY(), csigmaY);
 800          fFitterLinZ1.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
 801          }
 802          if (idelta > 0){
 803          ncl1++;
 804          fFitterLinY2.AddPoint(&x, currentCluster->GetY(), csigmaY);
 805          fFitterLinZ2.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
 806          }
 807          fFitterParY.AddPoint(&x, currentCluster->GetY(), csigmaY);
 808          fFitterParZ.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
 809       }  // loop over neighbourhood for fitter filling
 810
 811
 812
 813       if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
 814       if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
 815       if (nclFound < kDelta * kMinRatio) continue;    // if not enough clusters (7.5) found in neighbourhood goto next padrow
 816       fFitterParY.Eval();
 817       fFitterParZ.Eval();
 818       Double_t chi2 = (fFitterParY.GetChisquare() + fFitterParZ.GetChisquare()) / (2. * nclFound - 6.);
 819       TTreeSRedirector *cstream = GetDebugStreamer();
 820       if (cstream){
 821         (*cstream)<<"Cut9"<<
 822           "chi2="<<chi2<<
 823           "\n";
 824       }
 825       // REMOVE KINK
 826       // only when there are enough clusters (4) in each direction
 827       if (ncl0 > 4){
 828          fFitterLinY1.Eval();
 829          fFitterLinZ1.Eval();
 830       }
 831       if (ncl1 > 4){
 832          fFitterLinY2.Eval();
 833          fFitterLinZ2.Eval();
 834       }
 835
 836       if (ncl0 > 4 && ncl1 > 4){
 837          fFitterLinY1.GetCovarianceMatrix(matrixY0);
 838          fFitterLinY2.GetCovarianceMatrix(matrixY1);
 839          fFitterLinZ1.GetCovarianceMatrix(matrixZ0);
 840          fFitterLinZ2.GetCovarianceMatrix(matrixZ1);
 841          fFitterLinY2.GetParameters(paramY1);
 842          fFitterLinZ2.GetParameters(paramZ1);
 843          fFitterLinY1.GetParameters(paramY0);
 844          fFitterLinZ1.GetParameters(paramZ0);
 845          paramY0 -= paramY1;
 846          paramZ0 -= paramZ1;
 847          matrixY0 += matrixY1;
 848          matrixZ0 += matrixZ1;
 849          Double_t cchi2 = 0;
 850
 851          TMatrixD difY(2, 1, paramY0.GetMatrixArray());
 852          TMatrixD difYT(1, 2, paramY0.GetMatrixArray());
 853          matrixY0.Invert();
 854          TMatrixD mulY(matrixY0, TMatrixD::kMult, difY);
 855          TMatrixD chi2Y(difYT, TMatrixD::kMult, mulY);
 856          cchi2 += chi2Y(0, 0);
 857
 858          TMatrixD difZ(2, 1, paramZ0.GetMatrixArray());
 859          TMatrixD difZT(1, 2, paramZ0.GetMatrixArray());
 860          matrixZ0.Invert();
 861          TMatrixD mulZ(matrixZ0, TMatrixD::kMult, difZ);
 862          TMatrixD chi2Z(difZT, TMatrixD::kMult, mulZ);
 863          cchi2 += chi2Z(0, 0);
 864
 865          if (cstream){
 866            (*cstream)<<"Cut8"<<
 867              "chi2="<<cchi2<<
 868              "\n";
 869          }
 870       }
 871
 872       // current padrow has no kink
 873
 874       // get fit parameters from pol2 fit:
 875       Double_t paramY[4], paramZ[4];
 876       paramY[0] = fFitterParY.GetParameter(0);
 877       paramY[1] = fFitterParY.GetParameter(1);
 878       paramY[2] = fFitterParY.GetParameter(2);
 879       paramZ[0] = fFitterParZ.GetParameter(0);
 880       paramZ[1] = fFitterParZ.GetParameter(1);
 881       paramZ[2] = fFitterParZ.GetParameter(2);
 882
 883       Double_t tracky = paramY[0];
 884       Double_t trackz = paramZ[0];
 885       Float_t  deltay = tracky - cluster0->GetY();
 886       Float_t  deltaz = trackz - cluster0->GetZ();
 887       Float_t  angley = paramY[1] - paramY[0] / xref;
 888       Float_t  anglez = paramZ[1];
 889
 890       Float_t max = cluster0->GetMax();
 891       UInt_t isegment = cluster0->GetDetector() % 36;
 892       Int_t padSize = 0;                          // short pads
 893       if (cluster0->GetDetector() >= 36) {
 894          padSize = 1;                              // medium pads
 895          if (cluster0->GetRow() > 63) padSize = 2; // long pads
 896       }
 897
 898       // =========================================
 899       // wirte collected information to histograms
 900       // =========================================
 901
 902       TProfile *profAmpRow =  (TProfile*)fArrayAmpRow->At(sector);
 903       if ( irow >= kFirstLargePad) max /= kLargePadSize;
 904       Double_t smax = TMath::Sqrt(max);
 905       profAmpRow->Fill( (Double_t)cluster0->GetRow(), smax );
 906       TH1F *hisAmp =  (TH1F*)fArrayAmp->At(sector);
 907       hisAmp->Fill(smax);
 908
 909       // remove the following two lines one day:
 910       TProfile *profDriftLength = (TProfile*)fArrayChargeVsDriftlength->At(sector);
 911       profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax );
 912
 913       TProfile *profDriftLengthTmp = (TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isegment, padSize));
 914       profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax );
 915
 916       fHclusterPerPadrow->Fill(irow);   // fill histogram showing clusters per padrow
 917       Int_t ipad = TMath::Nint(cluster0->GetPad());
 918       Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
 919       fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
 920
 921
 922       TH3F * his3 = 0;
 923       his3 = (TH3F*)fRMSY->At(padSize);
 924       if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
 925       his3 = (TH3F*)fRMSZ->At(padSize);
 926       if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
 927
 928       his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
 929       if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
 930       his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
 931       if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
 932
 933
 934       // Fill resolution histograms
 935       Bool_t useForResol = kTRUE;
 936       if (fFitterParY.GetParError(0) > kErrorFraction * csigmaY) useForResol = kFALSE;
 937
 938       if (cstream){
 939         Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
 940         Float_t sy = cluster0->GetSigmaY2();
 941         Float_t sz = cluster0->GetSigmaZ2();
 942         (*cstream)<<"Resol0"<<
 943           "run="<<fRun<<              //  run number
 944           "event="<<fEvent<<          //  event number
 945           "time="<<fTime<<            //  time stamp of event
 946           "trigger="<<fTrigger<<      //  trigger
 947           "mag="<<fMagF<<             //  magnetic field
 948           "padSize="<<padSize<<
 949           "angley="<<angley<<
 950           "anglez="<<anglez<<
 951           "zdr="<<zdrift<<
 952           "dy="<<deltay<<
 953           "dz="<<deltaz<<
 954           "sy="<<sy<<
 955           "sz="<<sz<<
 956           "\n";
 957       }
 958
 959       if (useForResol){
 960          fDeltaY->Fill(deltay);
 961          fDeltaZ->Fill(deltaz);
 962          his3 = (TH3F*)fResolY->At(padSize);
 963          if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
 964          his3 = (TH3F*)fResolZ->At(padSize);
 965          if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
 966          his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
 967          if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
 968          his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
 969          if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );
 970       }
 971
 972       //=============================================================================================
 973
 974       if (useForResol && nclFound > 2 * kMinRatio * kDelta
 975           && irow % kDeltaWriteDebugStream == 0 && GetDebugLevel() > 4){
 976         if (GetDebugLevel() > 20) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow);
 977          FillResolutionHistoLocalDebugPart(track, cluster0, irow, angley, anglez, nclFound, kDelta);
 978       }  // if (useForResol && nclFound > 2 * kMinRatio * kDelta)
 979       //
 980       // Fill THN histograms
 981       //
 982       Double_t xvar[9];
 983       xvar[1]=padSize;
 984       xvar[2]=1.-TMath::Abs(cluster0->GetZ())/250.;
 985       xvar[3]=cluster0->GetMax();
 986       xvar[5]=angley;
 987       xvar[6]=anglez;
 988
 989       xvar[4]=cluster0->GetPad()-TMath::Nint(cluster0->GetPad());
 990       xvar[0]=deltay;
 991       fHisDeltaY->Fill(xvar);
 992       xvar[0]=TMath::Sqrt(cluster0->GetSigmaY2());
 993       fHisRMSY->Fill(xvar);
 994
 995       xvar[4]=cluster0->GetTimeBin()-TMath::Nint(cluster0->GetTimeBin());
 996       xvar[0]=deltaz;
 997       fHisDeltaZ->Fill(xvar);
 998       xvar[0]=TMath::Sqrt(cluster0->GetSigmaZ2());
 999       fHisRMSZ->Fill(xvar);
1000
1001    }    // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
1002 }  // FillResolutionHistoLocal(...)
1003
1004
1005
1006 void AliTPCcalibTracks::FillResolutionHistoLocalDebugPart(AliTPCseed *track, AliTPCclusterMI *cluster0, Int_t irow, Float_t  angley, Float_t  anglez, Int_t nclFound, Int_t kDelta) {
1007    //
1008    //  - debug part of FillResolutionHistoLocal -
1009    // called only for every kDeltaWriteDebugStream'th padrow, to avoid to much redundant data
1010    // called only for GetStreamLevel() > 4
1011    // fill resolution trees
1012    //
1013
1014    Int_t sector = cluster0->GetDetector();
1015    Float_t xref = cluster0->GetX();
1016    Int_t padSize = 0;                          // short pads
1017    if (cluster0->GetDetector() >= 36) {
1018       padSize = 1;                              // medium pads
1019       if (cluster0->GetRow() > 63) padSize = 2; // long pads
1020    }
1021
1022    static TLinearFitter fitY0(3, "pol2");
1023    static TLinearFitter fitZ0(3, "pol2");
1024    static TLinearFitter fitY2(5, "hyp4");
1025    static TLinearFitter fitZ2(5, "hyp4");
1026    static TLinearFitter fitY2Q(5, "hyp4");
1027    static TLinearFitter fitZ2Q(5, "hyp4");
1028    static TLinearFitter fitY2S(5, "hyp4");
1029    static TLinearFitter fitZ2S(5, "hyp4");
1030    fitY0.ClearPoints();
1031    fitZ0.ClearPoints();
1032    fitY2.ClearPoints();
1033    fitZ2.ClearPoints();
1034    fitY2Q.ClearPoints();
1035    fitZ2Q.ClearPoints();
1036    fitY2S.ClearPoints();
1037    fitZ2S.ClearPoints();
1038
1039    for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
1040      // loop over irow +- kDelta rows (neighboured rows)
1041      //
1042      //
1043      if (idelta == 0) continue;
1044      if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
1045      AliTPCclusterMI * cluster = track->GetClusterPointer(irow + idelta);
1046      if (!cluster) continue;
1047      if (cluster->GetType() < 0) continue;
1048      if (cluster->GetDetector() != sector) continue;
1049      Double_t x = cluster->GetX() - xref;
1050      Double_t sigmaY0 = fClusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley) );
1051      Double_t sigmaZ0 = fClusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez) );
1052      //
1053      Double_t sigmaYQ = fClusterParam->GetErrorQPar( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
1054      Double_t sigmaZQ = fClusterParam->GetErrorQPar( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
1055      Double_t sigmaYS = fClusterParam->GetErrorQParScaled( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())),
1056                                                            TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
1057      Double_t sigmaZS = fClusterParam->GetErrorQParScaled( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())),
1058                                                            TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
1059      Float_t rmsYFactor = fClusterParam->GetShapeFactor( 0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
1060                                                          TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
1061                                                          TMath::Sqrt(cluster0->GetSigmaY2()), 0 );
1062      Float_t rmsZFactor = fClusterParam->GetShapeFactor(0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
1063                                                         TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
1064                                                         TMath::Sqrt(cluster0->GetSigmaZ2()),0 );
1065      sigmaYS  = TMath::Sqrt(sigmaYS * sigmaYS + rmsYFactor * rmsYFactor / 12.);
1066      sigmaZS  = TMath::Sqrt(sigmaZS * sigmaZS + rmsZFactor * rmsZFactor / 12. + rmsYFactor * rmsYFactor / 24.);
1067      //
1068      if (kDelta != 0){
1069        fitY0.AddPoint(&x, cluster->GetY(), sigmaY0);
1070        fitZ0.AddPoint(&x, cluster->GetZ(), sigmaZ0);
1071      }
1072      Double_t xxx[4];
1073      xxx[0] = ( (idelta+irow) % 2 == 0 ) ? 1 : 0;
1074      xxx[1] = x;
1075      xxx[2] = ( (idelta+irow) % 2 == 0 ) ? x : 0;
1076      xxx[3] = x * x;
1077      fitY2.AddPoint(xxx, cluster->GetY(), sigmaY0);
1078      fitY2Q.AddPoint(xxx, cluster->GetY(), sigmaYQ);
1079      fitY2S.AddPoint(xxx, cluster->GetY(), sigmaYS);
1080      fitZ2.AddPoint(xxx, cluster->GetZ(), sigmaZ0);
1081      fitZ2Q.AddPoint(xxx, cluster->GetZ(), sigmaZQ);
1082      fitZ2S.AddPoint(xxx, cluster->GetZ(), sigmaZS);
1083      //
1084    }  // neigbouhood-loop
1085       //
1086    fitY0.Eval();
1087    fitZ0.Eval();
1088    fitY2.Eval();
1089    fitZ2.Eval();
1090    fitY2Q.Eval();
1091    fitZ2Q.Eval();
1092    fitY2S.Eval();
1093    fitZ2S.Eval();
1094    Float_t chi2Y0 = fitY0.GetChisquare() / (nclFound-3.);
1095    Float_t chi2Z0 = fitZ0.GetChisquare() / (nclFound-3.);
1096    Float_t chi2Y2 = fitY2.GetChisquare() / (nclFound-5.);
1097    Float_t chi2Z2 = fitZ2.GetChisquare() / (nclFound-5.);
1098    Float_t chi2Y2Q = fitY2Q.GetChisquare() / (nclFound-5.);
1099    Float_t chi2Z2Q = fitZ2Q.GetChisquare() / (nclFound-5.);
1100    Float_t chi2Y2S = fitY2S.GetChisquare() / (nclFound-5.);
1101    Float_t chi2Z2S = fitZ2S.GetChisquare() / (nclFound-5.);
1102    //
1103    static  TVectorD    parY0(3);
1104    static  TMatrixD    matY0(3, 3);
1105    static  TVectorD    parZ0(3);
1106    static  TMatrixD    matZ0(3, 3);
1107    fitY0.GetParameters(parY0);
1108    fitY0.GetCovarianceMatrix(matY0);
1109    fitZ0.GetParameters(parZ0);
1110    fitZ0.GetCovarianceMatrix(matZ0);
1111    //
1112    static  TVectorD    parY2(5);
1113    static  TMatrixD    matY2(5,5);
1114    static  TVectorD    parZ2(5);
1115    static  TMatrixD    matZ2(5,5);
1116    fitY2.GetParameters(parY2);
1117    fitY2.GetCovarianceMatrix(matY2);
1118    fitZ2.GetParameters(parZ2);
1119    fitZ2.GetCovarianceMatrix(matZ2);
1120    //
1121    static  TVectorD    parY2Q(5);
1122    static  TMatrixD    matY2Q(5,5);
1123    static  TVectorD    parZ2Q(5);
1124    static  TMatrixD    matZ2Q(5,5);
1125    fitY2Q.GetParameters(parY2Q);
1126    fitY2Q.GetCovarianceMatrix(matY2Q);
1127    fitZ2Q.GetParameters(parZ2Q);
1128    fitZ2Q.GetCovarianceMatrix(matZ2Q);
1129    static  TVectorD    parY2S(5);
1130    static  TMatrixD    matY2S(5,5);
1131    static  TVectorD    parZ2S(5);
1132    static  TMatrixD    matZ2S(5,5);
1133    fitY2S.GetParameters(parY2S);
1134    fitY2S.GetCovarianceMatrix(matY2S);
1135    fitZ2S.GetParameters(parZ2S);
1136    fitZ2S.GetCovarianceMatrix(matZ2S);
1137    Float_t sigmaY0   = TMath::Sqrt(matY0(0,0));
1138    Float_t sigmaZ0   = TMath::Sqrt(matZ0(0,0));
1139    Float_t sigmaDY0  = TMath::Sqrt(matY0(1,1));
1140    Float_t sigmaDZ0  = TMath::Sqrt(matZ0(1,1));
1141    Float_t sigmaY2   = TMath::Sqrt(matY2(1,1));
1142    Float_t sigmaZ2   = TMath::Sqrt(matZ2(1,1));
1143    Float_t sigmaDY2  = TMath::Sqrt(matY2(3,3));
1144    Float_t sigmaDZ2  = TMath::Sqrt(matZ2(3,3));
1145    Float_t sigmaY2Q  = TMath::Sqrt(matY2Q(1,1));
1146    Float_t sigmaZ2Q  = TMath::Sqrt(matZ2Q(1,1));
1147    Float_t sigmaDY2Q = TMath::Sqrt(matY2Q(3,3));
1148    Float_t sigmaDZ2Q = TMath::Sqrt(matZ2Q(3,3));
1149    Float_t sigmaY2S  = TMath::Sqrt(matY2S(1,1));
1150    Float_t sigmaZ2S  = TMath::Sqrt(matZ2S(1,1));
1151    Float_t sigmaDY2S = TMath::Sqrt(matY2S(3,3));
1152    Float_t sigmaDZ2S = TMath::Sqrt(matZ2S(3,3));
1153
1154    // Error parameters
1155    Float_t csigmaY0 = fClusterParam->GetError0Par(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(angley));
1156    Float_t csigmaZ0 = fClusterParam->GetError0Par(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(anglez));
1157    //
1158    Float_t csigmaYQ = fClusterParam->GetErrorQPar(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1159                                                   TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1160    Float_t csigmaZQ = fClusterParam->GetErrorQPar(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1161                                                   TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
1162    Float_t csigmaYS = fClusterParam->GetErrorQParScaled(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1163                                                         TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1164    Float_t csigmaZS = fClusterParam->GetErrorQParScaled(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1165                                                         TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
1166
1167
1168    // RMS parameters
1169    Float_t meanRMSY = 0;
1170    Float_t meanRMSZ = 0;
1171    Int_t   nclRMS = 0;
1172    for (Int_t idelta = -2; idelta <= 2; idelta++){
1173      // loop over neighbourhood
1174      if (idelta+irow < 0 || idelta+irow > 159) continue;
1175      //         if (idelta+irow>159) continue;
1176      AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta);
1177      if (!cluster) continue;
1178      meanRMSY += TMath::Sqrt(cluster->GetSigmaY2());
1179      meanRMSZ += TMath::Sqrt(cluster->GetSigmaZ2());
1180      nclRMS++;
1181    }
1182    meanRMSY /= nclRMS;
1183    meanRMSZ /= nclRMS;
1184
1185    Float_t rmsY      = TMath::Sqrt(cluster0->GetSigmaY2());
1186    Float_t rmsZ      = TMath::Sqrt(cluster0->GetSigmaZ2());
1187    Float_t rmsYT     = fClusterParam->GetRMSQ(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1188                                               TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1189    Float_t rmsZT     = fClusterParam->GetRMSQ(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1190                                               TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1191    Float_t rmsYT0    = fClusterParam->GetRMS0(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1192                                               TMath::Abs(angley));
1193    Float_t rmsZT0    = fClusterParam->GetRMS0(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1194                                                  TMath::Abs(anglez));
1195    Float_t rmsYSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1196                                                   TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1197    Float_t rmsZSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1198                                                   TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1199    Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1200                                                       TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
1201                                                       rmsY,meanRMSY);
1202    Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1203                                                       TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
1204                                                       rmsZ,meanRMSZ);
1205    //
1206    // cluster debug
1207    TTreeSRedirector *cstream = GetDebugStreamer();
1208    if (cstream){
1209      (*cstream)<<"ResolCl"<<    // valgrind 3   40,000 bytes in 1 blocks are possibly
1210        "run="<<fRun<<              //  run number
1211        "event="<<fEvent<<          //  event number
1212        "time="<<fTime<<            //  time stamp of event
1213        "trigger="<<fTrigger<<      //  trigger
1214        "mag="<<fMagF<<             //  magnetic field
1215        "Sector="<<sector<<
1216        "Cl.="<<cluster0<<
1217        "CSigmaY0="<<csigmaY0<<   // cluster errorY
1218        "CSigmaYQ="<<csigmaYQ<<   // cluster errorY - q scaled
1219        "CSigmaYS="<<csigmaYS<<   // cluster errorY - q scaled
1220        "CSigmaZ0="<<csigmaZ0<<   //
1221        "CSigmaZQ="<<csigmaZQ<<
1222        "CSigmaZS="<<csigmaZS<<
1223        "shapeYF="<<rmsYFactor<<
1224        "shapeZF="<<rmsZFactor<<
1225        "rmsY="<<rmsY<<
1226        "rmsZ="<<rmsZ<<
1227        "rmsYM="<<meanRMSY<<
1228        "rmsZM="<<meanRMSZ<<
1229        "rmsYT="<<rmsYT<<
1230        "rmsZT="<<rmsZT<<
1231        "rmsYT0="<<rmsYT0<<
1232        "rmsZT0="<<rmsZT0<<
1233        "rmsYS="<<rmsYSigma<<
1234        "rmsZS="<<rmsZSigma<<
1235        "padSize="<<padSize<<
1236        "Ncl="<<nclFound<<
1237        "PY0.="<<&parY0<<
1238        "PZ0.="<<&parZ0<<
1239        "SigmaY0="<<sigmaY0<<
1240        "SigmaZ0="<<sigmaZ0<<
1241        "angley="<<angley<<
1242        "anglez="<<anglez<<
1243        "\n";
1244
1245      //       tracklet dubug
1246      (*cstream)<<"ResolTr"<<
1247        "run="<<fRun<<              //  run number
1248        "event="<<fEvent<<          //  event number
1249        "time="<<fTime<<            //  time stamp of event
1250        "trigger="<<fTrigger<<      //  trigger
1251        "mag="<<fMagF<<             //  magnetic field
1252        "padSize="<<padSize<<
1253        "IPad="<<padSize<<
1254        "Sector="<<sector<<
1255        "Ncl="<<nclFound<<
1256        "chi2Y0="<<chi2Y0<<
1257        "chi2Z0="<<chi2Z0<<
1258        "chi2Y2="<<chi2Y2<<
1259        "chi2Z2="<<chi2Z2<<
1260        "chi2Y2Q="<<chi2Y2Q<<
1261        "chi2Z2Q="<<chi2Z2Q<<
1262        "chi2Y2S="<<chi2Y2S<<
1263        "chi2Z2S="<<chi2Z2S<<
1264        "PY0.="<<&parY0<<
1265        "PZ0.="<<&parZ0<<
1266        "PY2.="<<&parY2<<
1267        "PZ2.="<<&parZ2<<
1268        "PY2Q.="<<&parY2Q<<
1269        "PZ2Q.="<<&parZ2Q<<
1270        "PY2S.="<<&parY2S<<
1271        "PZ2S.="<<&parZ2S<<
1272        "SigmaY0="<<sigmaY0<<
1273        "SigmaZ0="<<sigmaZ0<<
1274        "SigmaDY0="<<sigmaDY0<<
1275        "SigmaDZ0="<<sigmaDZ0<<
1276        "SigmaY2="<<sigmaY2<<
1277        "SigmaZ2="<<sigmaZ2<<
1278        "SigmaDY2="<<sigmaDY2<<
1279        "SigmaDZ2="<<sigmaDZ2<<
1280        "SigmaY2Q="<<sigmaY2Q<<
1281        "SigmaZ2Q="<<sigmaZ2Q<<
1282        "SigmaDY2Q="<<sigmaDY2Q<<
1283        "SigmaDZ2Q="<<sigmaDZ2Q<<
1284        "SigmaY2S="<<sigmaY2S<<
1285        "SigmaZ2S="<<sigmaZ2S<<
1286        "SigmaDY2S="<<sigmaDY2S<<
1287        "SigmaDZ2S="<<sigmaDZ2S<<
1288        "angley="<<angley<<
1289        "anglez="<<anglez<<
1290        "\n";
1291    }
1292 }
1293
1294
1295
1296
1297
1298 TH2D * AliTPCcalibTracks::MakeDiff(TH2D * hfit, TF2 * func){
1299    //
1300    // creates a new histogram which contains the difference between
1301    // the histogram hfit and the function func
1302    //
1303   TH2D * result = (TH2D*)hfit->Clone();      // valgrind 3   40,139 bytes in 11 blocks are still reachable
1304   result->SetTitle(Form("%s fit residuals",result->GetTitle()));
1305   result->SetName(Form("%s fit residuals",result->GetName()));
1306   TAxis *xaxis = hfit->GetXaxis();
1307   TAxis *yaxis = hfit->GetYaxis();
1308   Double_t x[2];
1309   for (Int_t biny = 0; biny <= yaxis->GetNbins(); biny++) {
1310     x[1]  = yaxis->GetBinCenter(biny);
1311     for (Int_t binx = 0; binx <= xaxis->GetNbins(); binx++) {
1312       x[0]  = xaxis->GetBinCenter(binx);
1313       Int_t bin = hfit->GetBin(binx, biny);
1314       Double_t val = hfit->GetBinContent(bin);
1315 //      result->SetBinContent( bin, (val - func->Eval(x[0], x[1])) / func->Eval(x[0], x[1]) );
1316       result->SetBinContent( bin, (val / func->Eval(x[0], x[1])) - 1 );
1317     }
1318   }
1319   return result;
1320 }
1321
1322
1323 void  AliTPCcalibTracks::SetStyle() const {
1324    //
1325    // set style, can be called by all draw functions
1326    //
1327    gROOT->SetStyle("Plain");
1328    gStyle->SetFillColor(10);
1329    gStyle->SetPadColor(10);
1330    gStyle->SetCanvasColor(10);
1331    gStyle->SetStatColor(10);
1332    gStyle->SetPalette(1,0);
1333    gStyle->SetNumberContours(60);
1334 }
1335
1336
1337 void AliTPCcalibTracks::Draw(Option_t* opt){
1338    //
1339    // draw-function of AliTPCcalibTracks
1340    // will draws some exemplaric pictures
1341    //
1342
1343   if (GetDebugLevel() > 6) Info("Draw", "Drawing an exemplaric picture.");
1344    SetStyle();
1345    Double_t min = 0;
1346    Double_t max = 0;
1347    TCanvas *c1 = new TCanvas();
1348    c1->Divide(0, 3);
1349    TVirtualPad *upperThird = c1->GetPad(1);
1350    TVirtualPad *middleThird = c1->GetPad(2);
1351    TVirtualPad *lowerThird = c1->GetPad(3);
1352    upperThird->Divide(2,0);
1353    TVirtualPad *upleft  = upperThird->GetPad(1);
1354    TVirtualPad *upright = upperThird->GetPad(2);
1355    middleThird->Divide(2,0);
1356    TVirtualPad *middleLeft  = middleThird->GetPad(1);
1357    TVirtualPad *middleRight = middleThird->GetPad(2);
1358    lowerThird->Divide(2,0);
1359    TVirtualPad *downLeft  = lowerThird->GetPad(1);
1360    TVirtualPad *downRight = lowerThird->GetPad(2);
1361
1362
1363    upleft->cd(0);
1364    min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
1365    max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
1366    fDeltaY->SetAxisRange(min, max);
1367    fDeltaY->Fit("gaus","q","",min, max);        // valgrind 3  7 block possibly lost   2,400 bytes in 1 blocks are still reachable
1368    c1->Update();
1369
1370    upright->cd(0);
1371    max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
1372    min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
1373    fDeltaZ->SetAxisRange(min, max);
1374    fDeltaZ->Fit("gaus","q","",min, max);
1375    c1->Update();
1376
1377    middleLeft->cd();
1378    fHclus->Draw(opt);
1379
1380    middleRight->cd();
1381    fRejectedTracksHisto->Draw(opt);
1382    TPaveText *pt = new TPaveText(0.6,0.6, 0.8,0.8, "NDC");
1383    TText *t1 = pt->AddText("1: kEdgeThetaCutNoise");
1384    TText *t2 = pt->AddText("2: kMinClusters");
1385    TText *t3 = pt->AddText("3: kMinRatio");
1386    TText *t4 = pt->AddText("4: kMax1pt");
1387    t1 = t1; t2 = t2; t3 = t3; t4 = t4;    // avoid compiler warnings
1388    pt->SetToolTipText("Legend for failed cuts");
1389    pt->Draw();
1390
1391    downLeft->cd();
1392    fHclusterPerPadrowRaw->Draw(opt);
1393
1394    downRight->cd();
1395    fHclusterPerPadrow->Draw(opt);
1396 }
1397
1398
1399 void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){
1400    //
1401    // all functions are called, that produce pictures
1402    // the histograms are written to the directory 'pathName'
1403    // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file
1404    // 'stat' is also the number of minEntries for MakeResPlotsQTree
1405    //
1406
1407   if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
1408    MakeAmpPlots(stat, pathName);
1409    MakeDeltaPlots(pathName);
1410    FitResolutionNew(pathName);
1411    FitRMSNew(pathName);
1412    MakeChargeVsDriftLengthPlots(pathName);
1413 //    MakeResPlotsQ(1, 1);
1414    MakeResPlotsQTree(stat, pathName);
1415 }
1416
1417
1418 void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, const char* pathName){
1419    //
1420    // creates several plots:
1421    // fArrayAmp.ps, fArrayAmpRow.ps and DeltaYZ.ps
1422    // fArrayAmp.ps: one histogram per sector, the histogram shows the charge per cluster
1423    // fArrayAmpRow.ps: one histogram per sector, mean max. amplitude vs. pad row with landau fit
1424    // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
1425    // Empty histograms (sectors without data) are not written to file
1426    // the ps-files are written to the directory 'pathName', that is created if it does not exist
1427    // 'stat': only histograms with more than 'stat' entries are written to file.
1428    //
1429
1430    SetStyle();
1431    gSystem->MakeDirectory(pathName);
1432    gSystem->ChangeDirectory(pathName);
1433
1434    TCanvas* c1 = new TCanvas();     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
1435    TPostScript *ps;
1436    // histograms with accumulated amplitude for all IROCs and OROCs
1437    TH1F *allAmpHisIROC = ((TH1F*)(fArrayAmp->At(0))->Clone());
1438    allAmpHisIROC->SetName("Amp all IROCs");
1439    allAmpHisIROC->SetTitle("Amp all IROCs");
1440    TH1F *allAmpHisOROC = ((TH1F*)(fArrayAmp->At(36))->Clone());
1441    allAmpHisOROC->SetName("Amp all OROCs");
1442    allAmpHisOROC->SetTitle("Amp all OROCs");
1443
1444    ps = new TPostScript("fArrayAmp.ps", 112);
1445    if (GetDebugLevel() > -1) cout << "creating fArrayAmp.ps..." << endl;
1446    for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++){
1447       if ( ((TH1F*)fArrayAmp->At(i))->GetEntries() < stat  ) continue;
1448       ps->NewPage();
1449       ((TH1F*)fArrayAmp->At(i))->Draw();
1450       c1->Update();              // valgrind 3
1451       if (i > 0 && i < 36) {
1452          allAmpHisIROC->Add(((TH1F*)fArrayAmp->At(i)));
1453          allAmpHisOROC->Add(((TH1F*)fArrayAmp->At(i+36)));
1454       }
1455    }
1456    ps->NewPage();
1457    allAmpHisIROC->Draw();
1458    c1->Update();              // valgrind
1459    ps->NewPage();
1460    allAmpHisOROC->Draw();
1461    c1->Update();
1462    ps->Close();
1463    delete ps;
1464
1465    TH1F *his = 0;
1466    Double_t min = 0;
1467    Double_t max = 0;
1468    ps = new TPostScript("fArrayAmpRow.ps", 112);
1469    if (GetDebugLevel() > -1) cout << "creating fArrayAmpRow.ps..." << endl;
1470    for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++){
1471       his = (TH1F*)fArrayAmpRow->At(i);
1472       if (his->GetEntries() < stat) continue;
1473       ps->NewPage();
1474       min = TMath::Max( his->GetBinCenter(his->GetMaximumBin() )-100., 0.);
1475       max = his->GetBinCenter(5*his->GetMaximumBin()) + 100;
1476       his->SetAxisRange(min, max);
1477       his->Fit("pol3", "q", "", min, max);
1478       // his->Draw("error");    // don't use this line when you don't want to have empty pages in the ps-file
1479       c1->Update();
1480    }
1481    ps->Close();
1482    delete ps;
1483    delete c1;
1484    gSystem->ChangeDirectory("..");
1485 }
1486
1487
1488 void AliTPCcalibTracks::MakeDeltaPlots(const char* pathName){
1489    //
1490    // creates several plots:
1491    // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
1492    // the ps-files are written to the directory 'pathName', that is created if it does not exist
1493    //
1494
1495    SetStyle();
1496    gSystem->MakeDirectory(pathName);
1497    gSystem->ChangeDirectory(pathName);
1498
1499    TCanvas* c1 = new TCanvas();     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
1500    TPostScript *ps;
1501    Double_t min = 0;
1502    Double_t max = 0;
1503
1504    ps = new TPostScript("DeltaYZ.ps", 112);
1505    if (GetDebugLevel() > -1) cout << "creating DeltaYZ.ps..." << endl;
1506    min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
1507    max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
1508    fDeltaY->SetAxisRange(min, max);
1509    ps->NewPage();
1510    fDeltaY->Fit("gaus","q","",min, max);        // valgrind 3  7 block possibly lost   2,400 bytes in 1 blocks are still reachable
1511    c1->Update();
1512    ps->NewPage();
1513    max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
1514    min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
1515    fDeltaZ->SetAxisRange(min, max);
1516    fDeltaZ->Fit("gaus","q","",min, max);
1517    c1->Update();
1518    ps->Close();
1519    delete ps;
1520    delete c1;
1521    gSystem->ChangeDirectory("..");
1522 }
1523
1524
1525 void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(const char* pathName){
1526    //
1527    // creates charge vs. driftlength plots, one TProfile for each ROC
1528    // is not correct like this, should be one TProfile for each sector and padsize
1529    //
1530
1531    SetStyle();
1532    gSystem->MakeDirectory(pathName);
1533    gSystem->ChangeDirectory(pathName);
1534
1535    TCanvas* c1 = new TCanvas();     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
1536    TPostScript *ps;
1537    ps = new TPostScript("chargeVsDriftlengthOld.ps", 112);
1538    if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlength.ps..." << endl;
1539    TProfile *chargeVsDriftlengthAllIROCs = ((TProfile*)fArrayChargeVsDriftlength->At(0)->Clone());
1540    TProfile *chargeVsDriftlengthAllOROCs = ((TProfile*)fArrayChargeVsDriftlength->At(36)->Clone());
1541    chargeVsDriftlengthAllIROCs->SetName("allAmpHisIROC");
1542    chargeVsDriftlengthAllIROCs->SetTitle("charge vs. driftlength, all IROCs");
1543    chargeVsDriftlengthAllOROCs->SetName("allAmpHisOROC");
1544    chargeVsDriftlengthAllOROCs->SetTitle("charge vs. driftlength, all OROCs");
1545
1546    for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) {
1547       ((TProfile*)fArrayChargeVsDriftlength->At(i))->Draw();
1548       c1->Update();
1549       if (i > 0 && i < 36) {
1550          chargeVsDriftlengthAllIROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i)));
1551          chargeVsDriftlengthAllOROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i+36)));
1552       }
1553       ps->NewPage();
1554    }
1555    chargeVsDriftlengthAllIROCs->Draw();
1556    c1->Update();              // valgrind
1557    ps->NewPage();
1558    chargeVsDriftlengthAllOROCs->Draw();
1559    c1->Update();
1560    ps->Close();
1561    delete ps;
1562    delete c1;
1563    gSystem->ChangeDirectory("..");
1564 }
1565
1566
1567 void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(const char* pathName){
1568    //
1569    // creates charge vs. driftlength plots, one TProfile for each ROC
1570    // under development....
1571    //
1572
1573    SetStyle();
1574    gSystem->MakeDirectory(pathName);
1575    gSystem->ChangeDirectory(pathName);
1576
1577    TCanvas* c1 = new TCanvas("c1", "c1", 700,(Int_t)(TMath::Sqrt(2)*700));     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
1578 //    TCanvas c1("c1", "c1", 500,(sqrt(2)*500))
1579    c1->Divide(0,3);
1580    TPostScript *ps;
1581    ps = new TPostScript("chargeVsDriftlength.ps", 111);
1582    if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlengthNew.ps..." << endl;
1583
1584    TProfile *chargeVsDriftlengthAllShortPads  = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0)->Clone());
1585    TProfile *chargeVsDriftlengthAllMediumPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1)->Clone());
1586    TProfile *chargeVsDriftlengthAllLongPads   = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2)->Clone());
1587    chargeVsDriftlengthAllShortPads->SetName("allAmpHisShortPads");
1588    chargeVsDriftlengthAllShortPads->SetTitle("charge vs. driftlength, all sectors, short pads");
1589    chargeVsDriftlengthAllMediumPads->SetName("allAmpHisMediumPads");
1590    chargeVsDriftlengthAllMediumPads->SetTitle("charge vs. driftlength, all sectors, medium pads");
1591    chargeVsDriftlengthAllLongPads->SetName("allAmpHisLongPads");
1592    chargeVsDriftlengthAllLongPads->SetTitle("charge vs. driftlength, all sectors, long pads");
1593
1594    for (Int_t i = 0; i < 36; i++) {
1595       c1->cd(1)->SetGridx();
1596       c1->cd(1)->SetGridy();
1597       ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,0))->Draw();
1598       c1->cd(2)->SetGridx();
1599       c1->cd(2)->SetGridy();
1600       ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,1))->Draw();
1601       c1->cd(3)->SetGridx();
1602       c1->cd(3)->SetGridy();
1603       ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,2))->Draw();
1604       c1->Update();
1605       chargeVsDriftlengthAllShortPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0));
1606       chargeVsDriftlengthAllMediumPads->Add((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1));
1607       chargeVsDriftlengthAllLongPads->Add(  (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2));
1608       ps->NewPage();
1609    }
1610    c1->cd(1)->SetGridx();
1611    c1->cd(1)->SetGridy();
1612    chargeVsDriftlengthAllShortPads->Draw();
1613    c1->cd(2)->SetGridx();
1614    c1->cd(2)->SetGridy();
1615    chargeVsDriftlengthAllMediumPads->Draw();
1616    c1->cd(3)->SetGridx();
1617    c1->cd(3)->SetGridy();
1618    chargeVsDriftlengthAllLongPads->Draw();
1619    c1->Update();              // valgrind
1620 //    ps->NewPage();
1621    ps->Close();
1622    delete ps;
1623    delete c1;
1624    gSystem->ChangeDirectory("..");
1625 }
1626
1627
1628
1629 void AliTPCcalibTracks::FitResolutionNew(const char* pathName){
1630    //
1631    // calculates different resulution fits in Y and Z direction
1632    // the histograms are written to 'ResolutionYZ.ps'
1633    // writes calculated resolution to 'resol.txt'
1634    // all files are stored in the directory pathName
1635    //
1636
1637    SetStyle();
1638    gSystem->MakeDirectory(pathName);
1639    gSystem->ChangeDirectory(pathName);
1640
1641    TCanvas c;
1642    c.Divide(2,1);
1643    if (GetDebugLevel() > -1) cout << "creating ResolutionYZ.ps..." << endl;
1644    TPostScript *ps = new TPostScript("ResolutionYZ.ps", 112);
1645    TF2 *fres = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
1646    fres->SetParameter(0,0.02);
1647    fres->SetParameter(1,0.0054);
1648    fres->SetParameter(2,0.13);
1649
1650    TH1::AddDirectory(kTRUE);  // TH3F::FitSlicesZ() writes histograms into the current directory
1651
1652    // create histogramw for Y-resolution
1653    TH3F * hisResY0 = (TH3F*)fResolY->At(0);
1654    hisResY0->FitSlicesZ();
1655    TH2D * hisResY02 = (TH2D*)gDirectory->Get("Resol Y0_2");
1656    TH3F * hisResY1 = (TH3F*)fResolY->At(1);
1657    hisResY1->FitSlicesZ();
1658    TH2D * hisResY12 = (TH2D*)gDirectory->Get("Resol Y1_2");
1659    TH3F * hisResY2 = (TH3F*)fResolY->At(2);
1660    hisResY2->FitSlicesZ();
1661    TH2D * hisResY22 = (TH2D*)gDirectory->Get("Resol Y2_2");
1662     //
1663    ps->NewPage();
1664    c.cd(1);
1665    hisResY02->Fit(fres, "q");      // valgrind    132,072 bytes in 6 blocks are indirectly lost
1666    hisResY02->Draw("surf1");
1667    c.cd(2);
1668    MakeDiff(hisResY02,fres)->Draw("surf1");
1669    c.Update();
1670    //   c.SaveAs("ResolutionYPad0.eps");
1671    ps->NewPage();
1672    c.cd(1);
1673    hisResY12->Fit(fres, "q");
1674    hisResY12->Draw("surf1");
1675    c.cd(2);
1676    MakeDiff(hisResY12,fres)->Draw("surf1");
1677    c.Update();
1678    //   c.SaveAs("ResolutionYPad1.eps");
1679    ps->NewPage();
1680    c.cd(1);
1681    hisResY22->Fit(fres, "q");
1682    hisResY22->Draw("surf1");
1683    c.cd(2);
1684    MakeDiff(hisResY22,fres)->Draw("surf1");
1685    c.Update();
1686 //    c.SaveAs("ResolutionYPad2.eps");
1687
1688    // create histogramw for Z-resolution
1689    TH3F * hisResZ0 = (TH3F*)fResolZ->At(0);
1690    hisResZ0->FitSlicesZ();
1691    TH2D * hisResZ02 = (TH2D*)gDirectory->Get("Resol Z0_2");
1692    TH3F * hisResZ1 = (TH3F*)fResolZ->At(1);
1693    hisResZ1->FitSlicesZ();
1694    TH2D * hisResZ12 = (TH2D*)gDirectory->Get("Resol Z1_2");
1695    TH3F * hisResZ2 = (TH3F*)fResolZ->At(2);
1696    hisResZ2->FitSlicesZ();
1697    TH2D * hisResZ22 = (TH2D*)gDirectory->Get("Resol Z2_2");
1698
1699    ps->NewPage();
1700    c.cd(1);
1701    hisResZ02->Fit(fres, "q");
1702    hisResZ02->Draw("surf1");
1703    c.cd(2);
1704    MakeDiff(hisResZ02,fres)->Draw("surf1");
1705    c.Update();
1706 //    c.SaveAs("ResolutionZPad0.eps");
1707    ps->NewPage();
1708    c.cd(1);
1709    hisResZ12->Fit(fres, "q");
1710    hisResZ12->Draw("surf1");
1711    c.cd(2);
1712    MakeDiff(hisResZ12,fres)->Draw("surf1");
1713    c.Update();
1714 //    c.SaveAs("ResolutionZPad1.eps");
1715    ps->NewPage();
1716    c.cd(1);
1717    hisResZ22->Fit(fres, "q");
1718    hisResZ22->Draw("surf1");
1719    c.cd(2);
1720    MakeDiff(hisResZ22,fres)->Draw("surf1");
1721    c.Update();
1722 //    c.SaveAs("ResolutionZPad2.eps");
1723    ps->Close();
1724    delete ps;
1725
1726    // write calculated resoltuions to 'resol.txt'
1727    ofstream fresol("resol.txt");
1728    fresol<<"Pad 0.75 cm"<<"\n";
1729    hisResY02->Fit(fres, "q");                     // valgrind
1730    fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1731    hisResZ02->Fit(fres, "q");
1732    fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1733    //
1734    fresol<<"Pad 1.00 cm"<<1<<"\n";
1735    hisResY12->Fit(fres, "q");                     // valgrind
1736    fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1737    hisResZ12->Fit(fres, "q");
1738    fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1739    //
1740    fresol<<"Pad 1.50 cm"<<0<<"\n";
1741    hisResY22->Fit(fres, "q");
1742    fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1743    hisResZ22->Fit(fres, "q");
1744    fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1745
1746    TH1::AddDirectory(kFALSE);
1747    gSystem->ChangeDirectory("..");
1748    delete fres;
1749 }
1750
1751
1752 void AliTPCcalibTracks::FitRMSNew(const char* pathName){
1753    //
1754    // calculates different resulution-rms fits in Y and Z direction
1755    // the histograms are written to 'RMS_YZ.ps'
1756    // writes calculated resolution rms to 'rms.txt'
1757    // all files are stored in the directory pathName
1758    //
1759
1760    SetStyle();
1761    gSystem->MakeDirectory(pathName);
1762    gSystem->ChangeDirectory(pathName);
1763
1764    TCanvas c;        // valgrind 3   42,120 bytes in 405 blocks are still reachable   23,816 bytes in 229 blocks are still reachable
1765    c.Divide(2,1);
1766    if (GetDebugLevel() > -1) cout << "creating RMS_YZ.ps..." << endl;
1767    TPostScript *ps = new TPostScript("RMS_YZ.ps", 112);
1768    TF2 *frms = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
1769    frms->SetParameter(0,0.02);
1770    frms->SetParameter(1,0.0054);
1771    frms->SetParameter(2,0.13);
1772
1773    TH1::AddDirectory(kTRUE);  // TH3F::FitSlicesZ() writes histograms into the current directory
1774
1775    // create histogramw for Y-RMS
1776    TH3F * hisResY0 = (TH3F*)fRMSY->At(0);
1777    hisResY0->FitSlicesZ();
1778    TH2D * hisResY02 = (TH2D*)gDirectory->Get("RMS Y0_1");
1779    TH3F * hisResY1 = (TH3F*)fRMSY->At(1);
1780    hisResY1->FitSlicesZ();
1781    TH2D * hisResY12 = (TH2D*)gDirectory->Get("RMS Y1_1");
1782    TH3F * hisResY2 = (TH3F*)fRMSY->At(2);
1783    hisResY2->FitSlicesZ();
1784    TH2D * hisResY22 = (TH2D*)gDirectory->Get("RMS Y2_1");
1785    //
1786    ps->NewPage();
1787    c.cd(1);
1788    hisResY02->Fit(frms, "qn0");
1789    hisResY02->Draw("surf1");
1790    c.cd(2);
1791    MakeDiff(hisResY02,frms)->Draw("surf1");
1792    c.Update();
1793 //    c.SaveAs("RMSYPad0.eps");
1794    ps->NewPage();
1795    c.cd(1);
1796    hisResY12->Fit(frms, "qn0");               // valgrind   several blocks possibly lost
1797    hisResY12->Draw("surf1");
1798    c.cd(2);
1799    MakeDiff(hisResY12,frms)->Draw("surf1");
1800    c.Update();
1801 //    c.SaveAs("RMSYPad1.eps");
1802    ps->NewPage();
1803    c.cd(1);
1804    hisResY22->Fit(frms, "qn0");
1805    hisResY22->Draw("surf1");
1806    c.cd(2);
1807    MakeDiff(hisResY22,frms)->Draw("surf1");
1808    c.Update();
1809 //    c.SaveAs("RMSYPad2.eps");
1810
1811    // create histogramw for Z-RMS
1812    TH3F * hisResZ0 = (TH3F*)fRMSZ->At(0);
1813    hisResZ0->FitSlicesZ();
1814    TH2D * hisResZ02 = (TH2D*)gDirectory->Get("RMS Z0_1");
1815    TH3F * hisResZ1 = (TH3F*)fRMSZ->At(1);
1816    hisResZ1->FitSlicesZ();
1817    TH2D * hisResZ12 = (TH2D*)gDirectory->Get("RMS Z1_1");
1818    TH3F * hisResZ2 = (TH3F*)fRMSZ->At(2);
1819    hisResZ2->FitSlicesZ();
1820    TH2D * hisResZ22 = (TH2D*)gDirectory->Get("RMS Z2_1");
1821    //
1822    ps->NewPage();
1823    c.cd(1);
1824    hisResZ02->Fit(frms, "qn0");         // valgrind
1825    hisResZ02->Draw("surf1");
1826    c.cd(2);
1827    MakeDiff(hisResZ02,frms)->Draw("surf1");
1828    c.Update();
1829 //    c.SaveAs("RMSZPad0.eps");
1830    ps->NewPage();
1831    c.cd(1);
1832    hisResZ12->Fit(frms, "qn0");
1833    hisResZ12->Draw("surf1");
1834    c.cd(2);
1835    MakeDiff(hisResZ12,frms)->Draw("surf1");
1836    c.Update();
1837 //    c.SaveAs("RMSZPad1.eps");
1838    ps->NewPage();
1839    c.cd(1);
1840    hisResZ22->Fit(frms, "qn0");         // valgrind  1 block possibly lost
1841    hisResZ22->Draw("surf1");
1842    c.cd(2);
1843    MakeDiff(hisResZ22,frms)->Draw("surf1");
1844    c.Update();
1845 //    c.SaveAs("RMSZPad2.eps");
1846
1847    // write calculated resoltuion rms to 'rms.txt'
1848    ofstream filerms("rms.txt");
1849    filerms<<"Pad 0.75 cm"<<"\n";
1850    hisResY02->Fit(frms, "qn0");         // valgrind   23 blocks indirectly lost
1851    filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1852    hisResZ02->Fit(frms, "qn0");         // valgrind   23 blocks indirectly lost
1853    filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1854    //
1855    filerms<<"Pad 1.00 cm"<<1<<"\n";
1856    hisResY12->Fit(frms, "qn0");         // valgrind      3,256 bytes in 22 blocks are indirectly lost
1857    filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1858    hisResZ12->Fit(frms, "qn0");         // valgrind    66,036 bytes in 3 blocks are still reachable
1859    filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1860    //
1861    filerms<<"Pad 1.50 cm"<<0<<"\n";
1862    hisResY22->Fit(frms, "qn0");      // valgrind   40,139 bytes in 11 blocks are still reachable   330,180 bytes in 15 blocks are possibly lost
1863    filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1864    hisResZ22->Fit(frms, "qn0");
1865    filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1866
1867    TH1::AddDirectory(kFALSE);
1868    gSystem->ChangeDirectory("..");
1869    ps->Close();
1870    delete ps;
1871    delete frms;
1872 }
1873
1874
1875 void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){
1876    //
1877    // Make tree with resolution parameters
1878    // the result is written to 'resol.root' in directory 'pathname'
1879    // file information are available in fileInfo
1880    // available variables in the tree 'Resol':
1881    //  Entries: number of entries for this resolution point
1882    //  nbins:   number of bins that were accumulated
1883    //  Dim:     direction, Dim==0: y-direction, Dim==1: z-direction
1884    //  Pad:     padSize; short, medium and long
1885    //  Length:  pad length, 0.75, 1, 1.5
1886    //  QMean:   mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
1887    //  Zc:      center of middle bin in drift direction
1888    //  Zm:      mean dirftlength for accumulated Delta-Histograms
1889    //  Zs:      width of driftlength bin
1890    //  AngleC:  center of middle bin in Angle-Direction
1891    //  AngleM:  mean angle for accumulated Delta-Histograms
1892    //  AngleS:  width of Angle-bin
1893    //  Resol:   sigma for gaus fit through Delta-Histograms
1894    //  Sigma:   error of sigma for gaus fit through Delta Histograms
1895    //  MeanR:   mean of the Delta-Histogram
1896    //  SigmaR:  rms of the Delta-Histogram
1897    //  RMSm:    mean of the gaus fit through RMS-Histogram
1898    //  RMS:     sigma of the gaus fit through RMS-Histogram
1899    //  RMSe0:   error of mean of gaus fit in RMS-Histogram
1900    //  RMSe1:   error of sigma of gaus fit in RMS-Histogram
1901    //
1902
1903   if (GetDebugLevel() > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl;
1904   if (GetDebugLevel() > -1) cout << "    relax, the calculation will take a while..." << endl;
1905
1906    gSystem->MakeDirectory(pathName);
1907    gSystem->ChangeDirectory(pathName);
1908    TString kFileName = "resol.root";
1909    TTreeSRedirector fTreeResol(kFileName.Data());
1910
1911    TH3F *resArray[2][3][11];
1912    TH3F *rmsArray[2][3][11];
1913
1914    // load histograms from fArrayQDY and fArrayQDZ
1915    // into resArray and rmsArray
1916    // that is all we need here
1917    for (Int_t idim = 0; idim < 2; idim++){
1918       for (Int_t ipad = 0; ipad < 3; ipad++){
1919          for (Int_t iq = 0; iq <= 10; iq++){
1920             rmsArray[idim][ipad][iq]=0;
1921             resArray[idim][ipad][iq]=0;
1922             Int_t bin = GetBin(iq,ipad);
1923             TH3F *hresl = 0;
1924             if (idim == 0) hresl = (TH3F*)fArrayQDY->At(bin);
1925             if (idim == 1) hresl = (TH3F*)fArrayQDZ->At(bin);
1926             if (!hresl) continue;
1927             resArray[idim][ipad][iq] = (TH3F*) hresl->Clone();
1928             resArray[idim][ipad][iq]->SetDirectory(0);
1929             TH3F * hreslRMS = 0;
1930             if (idim == 0) hreslRMS = (TH3F*)fArrayQRMSY->At(bin);
1931             if (idim == 1) hreslRMS = (TH3F*)fArrayQRMSZ->At(bin);
1932             if (!hreslRMS) continue;
1933             rmsArray[idim][ipad][iq] = (TH3F*) hreslRMS->Clone();
1934             rmsArray[idim][ipad][iq]->SetDirectory(0);
1935          }
1936       }
1937    }
1938    if (GetDebugLevel() > -1) cout << "Histograms loaded, starting to proces..." << endl;
1939
1940    //--------------------------------------------------------------------------------------------
1941
1942    char name[200];
1943    Double_t qMean;
1944    Double_t zMean, angleMean, zCenter, angleCenter;
1945    Double_t zSigma, angleSigma;
1946    TH1D *projectionRes = new TH1D("projectionRes", "projectionRes", 50, -1, 1);
1947    TH1D *projectionRms = new TH1D("projectionRms", "projectionRms", 50, -1, 1);
1948    TF1 *fitFunction = new TF1("fitFunction", "gaus");
1949    Float_t entriesQ = 0;
1950    Int_t loopCounter = 1;
1951
1952    for (Int_t idim = 0; idim < 2; idim++){
1953       // Loop y-z corrdinate
1954       for (Int_t ipad = 0; ipad < 3; ipad++){
1955          // loop pad type
1956          for (Int_t iq = -1; iq < 10; iq++){
1957             // LOOP Q
1958            if (GetDebugLevel() > -1)
1959                cout << "Loop-counter, this is loop " << loopCounter << " of 66, ("
1960                   << (Int_t)((loopCounter)/66.*100) << "% done), "
1961                   << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << "  \r" << std::flush;
1962             loopCounter++;
1963             TH3F *hres = 0;
1964             TH3F *hrms = 0;
1965             qMean    = 0;
1966             entriesQ = 0;
1967
1968             // calculate qMean
1969             if (iq == -1){
1970                // integrated spectra
1971                for (Int_t iql = 0; iql < 10; iql++){
1972                   Int_t bin = GetBin(iql,ipad);
1973                   TH3F *hresl = resArray[idim][ipad][iql];
1974                   TH3F *hrmsl = rmsArray[idim][ipad][iql];
1975                   if (!hresl) continue;
1976                   if (!hrmsl) continue;
1977                   entriesQ += hresl->GetEntries();
1978                   qMean += hresl->GetEntries() * GetQ(bin);
1979                   if (!hres) {
1980                      hres = (TH3F*)hresl->Clone();
1981                      hrms = (TH3F*)hrmsl->Clone();
1982                   }
1983                   else{
1984                      hres->Add(hresl);
1985                      hrms->Add(hrmsl);
1986                   }
1987                }
1988                qMean /= entriesQ;
1989                qMean *= -1.;  // integral mean charge
1990             }
1991             else {
1992                // loop over neighboured Q-bins
1993                // accumulate entries from neighboured Q-bins
1994                for (Int_t iql = iq - 1; iql <= iq + 1; iql++){
1995                   if (iql < 0) continue;
1996                   Int_t bin = GetBin(iql,ipad);
1997                   TH3F * hresl = resArray[idim][ipad][iql];
1998                   TH3F * hrmsl = rmsArray[idim][ipad][iql];
1999                   if (!hresl) continue;
2000                   if (!hrmsl) continue;
2001                   entriesQ += hresl->GetEntries();
2002                   qMean += hresl->GetEntries() * GetQ(bin);
2003                   if (!hres) {
2004                      hres = (TH3F*) hresl->Clone();
2005                      hrms = (TH3F*) hrmsl->Clone();
2006                   }
2007                   else{
2008                      hres->Add(hresl);
2009                      hrms->Add(hrmsl);
2010                   }
2011                }
2012                qMean/=entriesQ;
2013             }
2014
2015             TAxis *xAxisDriftLength = hres->GetXaxis();   // driftlength / z - axis
2016             TAxis *yAxisAngle       = hres->GetYaxis();   // angle axis
2017             TAxis *zAxisDelta       = hres->GetZaxis();   // delta axis
2018             TAxis *zAxisRms         = hrms->GetZaxis();   // rms axis
2019
2020             // loop over all angle bins
2021             for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) {
2022                angleCenter = yAxisAngle->GetBinCenter(ibinyAngle);
2023                // loop over all driftlength bins
2024                for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) {
2025                   zCenter    = xAxisDriftLength->GetBinCenter(ibinxDL);
2026                   zSigma     = xAxisDriftLength->GetBinWidth(ibinxDL);
2027                   angleSigma = yAxisAngle->GetBinWidth(ibinyAngle);
2028                   zMean      = zCenter;      // changens, when more statistic is accumulated
2029                   angleMean  = angleCenter;  // changens, when more statistic is accumulated
2030
2031                   // create 2 1D-Histograms, projectionRes and projectionRms
2032                   // these histograms are delta histograms for given direction, padSize, chargeBin,
2033                   // angleBin and driftLengthBin
2034                   // later on they will be fitted with a gausian, its sigma is the resoltuion...
2035                   sprintf(name,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter);
2036                   // TH1D * projectionRes = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
2037                   projectionRes->SetNameTitle(name, name);
2038                   sprintf(name,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter);
2039                   // TH1D * projectionRms =  new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
2040                   projectionRms->SetNameTitle(name, name);
2041
2042                   projectionRes->Reset();
2043                   projectionRes->SetBins(zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
2044                   projectionRms->Reset();
2045                   projectionRms->SetBins(zAxisRms->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
2046                   projectionRes->SetDirectory(0);
2047                   projectionRms->SetDirectory(0);
2048
2049                   Double_t entries = 0;
2050                   Int_t    nbins   = 0;   // counts, how many bins were accumulated
2051                   zMean     = 0;
2052                   angleMean = 0;
2053                   entries   = 0;
2054
2055                   // fill projectionRes and projectionRms for given dim, ipad and iq,
2056                   // as well as for given angleBin and driftlengthBin
2057                   // if this gives not enough statistic, include neighbourhood
2058                   // (angle and driftlength) successifely
2059                   for (Int_t dbin = 0; dbin <= 8; dbin++){              // delta-bins around centered angleBin and driftlengthBin
2060                      for (Int_t dbiny2 = -1; dbiny2 <= 1; dbiny2++) {   // delta-bins in angle direction
2061                         for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ // delta-bins in driftlength direction
2062                            if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue;   // add each bin only one time !
2063                            Int_t binx2 = ibinxDL + dbinx2;                       // position variable in x (driftlength) direction
2064                            Int_t biny2 = ibinyAngle + dbiny2;                    // position variable in y (angle)  direction
2065                            if (binx2 < 1 || biny2 < 1) continue;                 // don't go out of the histogram!
2066                            if (binx2 >= xAxisDriftLength->GetNbins()) continue;  // don't go out of the histogram!
2067                            if (biny2 >= yAxisAngle->GetNbins()) continue;        // don't go out of the histogram!
2068                            nbins++;                                              // count the number of accumulated bins
2069                            // Fill resolution histo
2070                            for (Int_t ibin3 = 1; ibin3 < zAxisDelta->GetNbins(); ibin3++) {
2071                               // Int_t content = (Int_t)hres->GetBinContent(binx2, biny2, ibin3);     // unused variable
2072                               projectionRes->Fill(zAxisDelta->GetBinCenter(ibin3), hres->GetBinContent(binx2, biny2, ibin3));
2073                               entries   += hres->GetBinContent(binx2, biny2, ibin3);
2074                               zMean     += hres->GetBinContent(binx2, biny2, ibin3) * xAxisDriftLength->GetBinCenter(binx2);
2075                               angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yAxisAngle->GetBinCenter(biny2);
2076                            }  // ibin3 loop
2077                            // fill RMS histo
2078                            for (Int_t ibin3 = 1; ibin3 < zAxisRms->GetNbins(); ibin3++) {
2079                               projectionRms->Fill(zAxisRms->GetBinCenter(ibin3), hrms->GetBinContent(binx2, biny2, ibin3));
2080                            }
2081                         }  //dbinx2 loop
2082                         if (entries > minEntries) break; // enough statistic accumulated
2083                      }  // dbiny2 loop
2084                      if (entries > minEntries) break;    // enough statistic accumulated
2085                   }  // dbin loop
2086                   if ( entries< minEntries) continue;  // when it was absolutly impossible to get enough statistic, don't write this point into the resolution tree
2087                   zMean /= entries;
2088                   angleMean /= entries;
2089
2090                   if (entries > minEntries) {
2091                      //  when enough statistic is accumulated
2092                      //  fit Delta histograms with a gausian
2093                      //  of the gausian is the resolution (resol), its fit error is sigma
2094                      //  store also mean and RMS of the histogram
2095                      Float_t xmin     = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
2096                      Float_t xmax     = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
2097
2098 //                      projectionRes->Fit("gaus", "q0", "", xmin, xmax);
2099 //                      Float_t resol    = projectionRes->GetFunction("gaus")->GetParameter(2);
2100 //                      Float_t sigma    = projectionRes->GetFunction("gaus")->GetParError(2);
2101                      fitFunction->SetMaximum(xmax);
2102                      fitFunction->SetMinimum(xmin);
2103                      projectionRes->Fit("fitFunction", "qN0", "", xmin, xmax);
2104                      Float_t resol    = fitFunction->GetParameter(2);
2105                      Float_t sigma    = fitFunction->GetParError(2);
2106
2107                      Float_t meanR    = projectionRes->GetMean();
2108                      Float_t sigmaR   = projectionRes->GetRMS();
2109                      // fit also RMS histograms with a gausian
2110                      // store mean and sigma of the gausian in rmsMean and rmsSigma
2111                      // store also the fit errors in errorRMS and errorSigma
2112                      xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
2113                      xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
2114
2115 //                      projectionRms->Fit("gaus","q0","",xmin,xmax);
2116 //                      Float_t rmsMean    = projectionRms->GetFunction("gaus")->GetParameter(1);
2117 //                      Float_t rmsSigma   = projectionRms->GetFunction("gaus")->GetParameter(2);
2118 //                      Float_t errorRMS   = projectionRms->GetFunction("gaus")->GetParError(1);
2119 //                      Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2);
2120                      projectionRms->Fit("fitFunction", "qN0", "", xmin, xmax);
2121                      Float_t rmsMean    = fitFunction->GetParameter(1);
2122                      Float_t rmsSigma   = fitFunction->GetParameter(2);
2123                      Float_t errorRMS   = fitFunction->GetParError(1);
2124                      Float_t errorSigma = fitFunction->GetParError(2);
2125
2126                      Float_t length = 0.75;
2127                      if (ipad == 1) length = 1;
2128                      if (ipad == 2) length = 1.5;
2129
2130                      fTreeResol<<"Resol"<<
2131                         "Entries="<<entries<<      // number of entries for this resolution point
2132                         "nbins="<<nbins<<          // number of bins that were accumulated
2133                         "Dim="<<idim<<             // direction, Dim==0: y-direction, Dim==1: z-direction
2134                         "Pad="<<ipad<<             // padSize; short, medium and long
2135                         "Length="<<length<<        // pad length, 0.75, 1, 1.5
2136                         "QMean="<<qMean<<          // mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
2137                         "Zc="<<zCenter<<           // center of middle bin in drift direction
2138                         "Zm="<<zMean<<             // mean dirftlength for accumulated Delta-Histograms
2139                         "Zs="<<zSigma<<            // width of driftlength bin
2140                         "AngleC="<<angleCenter<<   // center of middle bin in Angle-Direction
2141                         "AngleM="<<angleMean<<     // mean angle for accumulated Delta-Histograms
2142                         "AngleS="<<angleSigma<<    // width of Angle-bin
2143                         "Resol="<<resol<<          // sigma for gaus fit through Delta-Histograms
2144                         "Sigma="<<sigma<<          // error of sigma for gaus fit through Delta Histograms
2145                         "MeanR="<<meanR<<          // mean of the Delta-Histogram
2146                         "SigmaR="<<sigmaR<<        // rms of the Delta-Histogram
2147                         "RMSm="<<rmsMean<<         // mean of the gaus fit through RMS-Histogram
2148                         "RMSs="<<rmsSigma<<        // sigma of the gaus fit through RMS-Histogram
2149                         "RMSe0="<<errorRMS<<       // error of mean of gaus fit in RMS-Histogram
2150                         "RMSe1="<<errorSigma<<     // error of sigma of gaus fit in RMS-Histogram
2151                         "\n";
2152                      if (GetDebugLevel() > 5) {
2153                         projectionRes->SetDirectory(fTreeResol.GetFile());
2154                         projectionRes->Write(projectionRes->GetName());
2155                         projectionRes->SetDirectory(0);
2156                         projectionRms->SetDirectory(fTreeResol.GetFile());
2157                         projectionRms->Write(projectionRms->GetName());
2158                         projectionRes->SetDirectory(0);
2159                      }
2160                   }  // if (projectionRes->GetSum() > minEntries)
2161                }  // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++)
2162             }  // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++)
2163
2164          }  // iq-loop
2165       }  // ipad-loop
2166    }  // idim-loop
2167    delete projectionRes;
2168    delete projectionRms;
2169
2170 //    TFile resolFile(fTreeResol.GetFile());
2171    TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries));
2172    fileInfo.Write("fileInfo");
2173 //    resolFile.Close();
2174 //    fTreeResol.GetFile()->Close();
2175    if (GetDebugLevel() > -1) cout << endl;
2176    if (GetDebugLevel() > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl;
2177    gSystem->ChangeDirectory("..");
2178 }
2179
2180
2181
2182
2183
2184 Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
2185    //
2186    // function to merge several AliTPCcalibTracks objects after PROOF calculation
2187    // The object's histograms are merged via their merge functions
2188    // Be carefull: histograms are linked to a file, switch this off by TH1::AddDirectory(kFALSE) !!!
2189    //
2190
2191   if (GetDebugLevel() > 0) cout << " *****  this is AliTPCcalibTracks::Merge(TCollection *collectionList)  *****"<< endl;
2192    if (!collectionList) return 0;
2193    if (collectionList->IsEmpty()) return -1;
2194
2195    if (GetDebugLevel() > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl;     //    REMOVE THIS LINE!!!!!!!!!!!!!!!!!1
2196    if (GetDebugLevel() > 5) cout << " the list in the merge-function looks as follows: " << endl;
2197    collectionList->Print();
2198
2199    // create a list for each data member
2200    TList* deltaYList = new TList;
2201    TList* deltaZList = new TList;
2202    TList* arrayAmpRowList = new TList;
2203    TList* rejectedTracksList = new TList;
2204    TList* hclusList = new TList;
2205    TList* clusterCutHistoList = new TList;
2206    TList* arrayAmpList = new TList;
2207    TList* arrayQDYList = new TList;
2208    TList* arrayQDZList = new TList;
2209    TList* arrayQRMSYList = new TList;
2210    TList* arrayQRMSZList = new TList;
2211    TList* arrayChargeVsDriftlengthList = new TList;
2212    TList* calPadRegionChargeVsDriftlengthList = new TList;
2213    TList* hclusterPerPadrowList = new TList;
2214    TList* hclusterPerPadrowRawList = new TList;
2215    TList* resolYList = new TList;
2216    TList* resolZList = new TList;
2217    TList* rMSYList = new TList;
2218    TList* rMSZList = new TList;
2219
2220 //    TList* nRowsList = new TList;
2221 //    TList* nSectList = new TList;
2222 //    TList* fileNoList = new TList;
2223
2224    TIterator *listIterator = collectionList->MakeIterator();
2225    AliTPCcalibTracks *calibTracks = 0;
2226    if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;
2227    Int_t counter = 0;
2228    while ( (calibTracks = dynamic_cast<AliTPCcalibTracks*> (listIterator->Next())) ){
2229       // loop over all entries in the collectionList and get dataMembers into lists
2230       if (!calibTracks) continue;
2231
2232       deltaYList->Add( calibTracks->GetfDeltaY() );
2233       deltaZList->Add( calibTracks->GetfDeltaZ() );
2234       arrayAmpRowList->Add(calibTracks->GetfArrayAmpRow());
2235       arrayAmpList->Add(calibTracks->GetfArrayAmp());
2236       arrayQDYList->Add(calibTracks->GetfArrayQDY());
2237       arrayQDZList->Add(calibTracks->GetfArrayQDZ());
2238       arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
2239       arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ());
2240       resolYList->Add(calibTracks->GetfResolY());
2241       resolZList->Add(calibTracks->GetfResolZ());
2242       rMSYList->Add(calibTracks->GetfRMSY());
2243       rMSZList->Add(calibTracks->GetfRMSZ());
2244       arrayChargeVsDriftlengthList->Add(calibTracks->GetfArrayChargeVsDriftlength());
2245       calPadRegionChargeVsDriftlengthList->Add(calibTracks->GetCalPadRegionchargeVsDriftlength());
2246       hclusList->Add(calibTracks->GetfHclus());
2247       rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
2248       clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
2249       hclusterPerPadrowList->Add(calibTracks->GetfHclusterPerPadrow());
2250       hclusterPerPadrowRawList->Add(calibTracks->GetfHclusterPerPadrowRaw());
2251       //
2252       if (fCalPadClusterPerPad && calibTracks->GetfCalPadClusterPerPad())
2253         fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
2254       //      fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
2255       counter++;
2256       if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
2257       AddHistos(calibTracks);
2258    }
2259
2260
2261    // merge data members
2262    if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl;
2263    fDeltaY->Merge(deltaYList);
2264    fDeltaZ->Merge(deltaZList);
2265    fHclus->Merge(hclusList);
2266    fClusterCutHisto->Merge(clusterCutHistoList);
2267    fRejectedTracksHisto->Merge(rejectedTracksList);
2268    fHclusterPerPadrow->Merge(hclusterPerPadrowList);
2269    fHclusterPerPadrowRaw->Merge(hclusterPerPadrowRawList);
2270
2271    TObjArray* objarray = 0;
2272    TH1* hist = 0;
2273    TList* histList = 0;
2274    TIterator *objListIterator = 0;
2275
2276    if (GetDebugLevel() > 0) cout << "merging fArrayAmpRows..." << endl;
2277    // merge fArrayAmpRows
2278    for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++ ) {  // loop over data member, i<72
2279       objListIterator = arrayAmpRowList->MakeIterator();
2280       histList = new TList;
2281       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2282          // loop over arrayAmpRowList, get TObjArray, get object at position i, cast it into TProfile
2283          if (!objarray) continue;
2284          hist = (TProfile*)(objarray->At(i));
2285          histList->Add(hist);
2286       }
2287       ((TProfile*)(fArrayAmpRow->At(i)))->Merge(histList);
2288       delete histList;
2289       delete objListIterator;
2290    }
2291
2292    if (GetDebugLevel() > 0) cout << "merging fArrayAmps..." << endl;
2293    // merge fArrayAmps
2294    for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++ ) {  // loop over data member, i<72
2295       TIterator *cobjListIterator = arrayAmpList->MakeIterator();
2296       histList = new TList;
2297       while (( objarray =  (TObjArray*)cobjListIterator->Next() )) {
2298          // loop over arrayAmpList, get TObjArray, get object at position i, cast it into TH1F
2299          if (!objarray) continue;
2300          hist = (TH1F*)(objarray->At(i));
2301          histList->Add(hist);
2302       }
2303       ((TH1F*)(fArrayAmp->At(i)))->Merge(histList);
2304       delete histList;
2305       delete cobjListIterator;
2306    }
2307
2308    if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
2309    // merge fArrayQDY
2310    for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
2311       objListIterator = arrayQDYList->MakeIterator();
2312       histList = new TList;
2313       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2314          // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F
2315          if (!objarray) continue;
2316          hist = (TH3F*)(objarray->At(i));
2317          histList->Add(hist);
2318       }
2319       ((TH3F*)(fArrayQDY->At(i)))->Merge(histList);
2320       delete histList;
2321       delete objListIterator;
2322    }
2323
2324    if (GetDebugLevel() > 0) cout << "merging fArrayQDZ..." << endl;
2325    // merge fArrayQDZ
2326    for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300
2327       objListIterator = arrayQDZList->MakeIterator();
2328       histList = new TList;
2329       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2330          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2331          if (!objarray) continue;
2332          hist = (TH3F*)(objarray->At(i));
2333          histList->Add(hist);
2334       }
2335       ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList);
2336       delete histList;
2337       delete objListIterator;
2338    }
2339
2340    if (GetDebugLevel() > 0) cout << "merging fArrayQRMSY..." << endl;
2341    // merge fArrayQRMSY
2342    for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300
2343       objListIterator = arrayQRMSYList->MakeIterator();
2344       histList = new TList;
2345       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2346          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2347          if (!objarray) continue;
2348          hist = (TH3F*)(objarray->At(i));
2349          histList->Add(hist);
2350       }
2351       ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList);
2352       delete histList;
2353       delete objListIterator;
2354    }
2355
2356    if (GetDebugLevel() > 0) cout << "merging fArrayQRMSZ..." << endl;
2357    // merge fArrayQRMSZ
2358    for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300
2359       objListIterator = arrayQRMSZList->MakeIterator();
2360       histList = new TList;
2361       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2362          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2363          if (!objarray) continue;
2364          hist = (TH3F*)(objarray->At(i));
2365          histList->Add(hist);
2366       }
2367       ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList);
2368       delete histList;
2369       delete objListIterator;
2370    }
2371
2372    if (GetDebugLevel() > 0) cout << "merging fArrayChargeVsDriftlength..." << endl;
2373    // merge fArrayChargeVsDriftlength
2374    for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) { // loop over data member, i < 300
2375       objListIterator = arrayChargeVsDriftlengthList->MakeIterator();
2376       histList = new TList;
2377       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2378          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TProfile
2379          if (!objarray) continue;
2380          hist = (TProfile*)(objarray->At(i));
2381          histList->Add(hist);
2382       }
2383       ((TProfile*)(fArrayChargeVsDriftlength->At(i)))->Merge(histList);
2384       delete histList;
2385       delete objListIterator;
2386    }
2387
2388    if (GetDebugLevel() > 0) cout << "merging fcalPadRegionChargeVsDriftlength..." << endl;
2389    // merge fcalPadRegionChargeVsDriftlength
2390    AliTPCCalPadRegion *cpr = 0x0;
2391
2392    /*
2393    TIterator *regionIterator = fcalPadRegionChargeVsDriftlength->MakeIterator();
2394    while (hist = (TProfile*)regionIterator->Next()) {
2395       // loop over all calPadRegion's in destination calibTracks object
2396          objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
2397          while (( cpr =  (AliTPCCalPadRegion*)objListIterator->Next() )) {
2398
2399
2400       hist->Merge(...);
2401    }
2402    */
2403
2404    for (UInt_t isec = 0; isec < 36; isec++) {
2405       for (UInt_t padSize = 0; padSize < 3; padSize++){
2406          objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
2407          histList = new TList;
2408          while (( cpr =  (AliTPCCalPadRegion*)objListIterator->Next() )) {
2409             // loop over calPadRegionChargeVsDriftlengthList, get AliTPCCalPadRegion, get object
2410             if (!cpr) continue;
2411             hist = (TProfile*)cpr->GetObject(isec, padSize);
2412             histList->Add(hist);
2413          }
2414          ((TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isec, padSize)))->Merge(histList);
2415          delete histList;
2416          delete objListIterator;
2417       }
2418    }
2419
2420
2421
2422
2423    if (GetDebugLevel() > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl;
2424    // merge fResolY
2425    for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3
2426       objListIterator = resolYList->MakeIterator();
2427       histList = new TList;
2428       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2429          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2430          if (!objarray) continue;
2431          hist = (TH3F*)(objarray->At(i));
2432          histList->Add(hist);
2433       }
2434       ((TH3F*)(fResolY->At(i)))->Merge(histList);
2435       delete histList;
2436       delete objListIterator;
2437    }
2438
2439    // merge fResolZ
2440    for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3
2441       objListIterator = resolZList->MakeIterator();
2442       histList = new TList;
2443       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2444          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2445          if (!objarray) continue;
2446          hist = (TH3F*)(objarray->At(i));
2447          histList->Add(hist);
2448       }
2449       ((TH3F*)(fResolZ->At(i)))->Merge(histList);
2450       delete histList;
2451       delete objListIterator;
2452    }
2453
2454    // merge fRMSY
2455    for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3
2456       objListIterator = rMSYList->MakeIterator();
2457       histList = new TList;
2458       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2459          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2460          if (!objarray) continue;
2461          hist = (TH3F*)(objarray->At(i));
2462          histList->Add(hist);
2463       }
2464       ((TH3F*)(fRMSY->At(i)))->Merge(histList);
2465       delete histList;
2466       delete objListIterator;
2467    }
2468
2469    // merge fRMSZ
2470    for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3
2471       objListIterator = rMSZList->MakeIterator();
2472       histList = new TList;
2473       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2474          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2475          if (!objarray) continue;
2476          hist = (TH3F*)(objarray->At(i));
2477          histList->Add(hist);
2478       }
2479       ((TH3F*)(fRMSZ->At(i)))->Merge(histList);
2480       delete histList;
2481       delete objListIterator;
2482    }
2483
2484    delete deltaYList;
2485    delete deltaZList;
2486    delete arrayAmpRowList;
2487    delete arrayAmpList;
2488    delete arrayQDYList;
2489    delete arrayQDZList;
2490    delete arrayQRMSYList;
2491    delete arrayQRMSZList;
2492    delete resolYList;
2493    delete resolZList;
2494    delete rMSYList;
2495    delete rMSZList;
2496    delete listIterator;
2497
2498    if (GetDebugLevel() > 0) cout << "merging done!" << endl;
2499
2500    return 1;
2501 }
2502
2503
2504 AliTPCcalibTracks* AliTPCcalibTracks::TestMerge(AliTPCcalibTracks *ct, AliTPCClusterParam *clusterParam, Int_t nCalTracks){
2505    //
2506    // function to test AliTPCcalibTrack::Merge:
2507    // in the file 'f' is a AliTPCcalibTrack object with name "calibTracks"
2508    // this object is appended 'nCalTracks' times to a TList
2509    // A new AliTPCcalibTrack object is created which merges the list
2510    // this object is returned
2511    //
2512    /*
2513    // .L AliTPCcalibTracks.cxx+g
2514    .L libTPCcalib.so
2515    TFile f("Output.root");
2516    AliTPCcalibTracks* calTracks = (AliTPCcalibTracks*)f.Get("calibTracks");
2517    //f.Close();
2518    TFile clusterParamFile("/u/lbozyk/calibration/workdir/calibTracks/TPCClusterParam.root");
2519    AliTPCClusterParam *clusterParam  =  (AliTPCClusterParam *) clusterParamFile.Get("Param");
2520    clusterParamFile.Close();
2521
2522    AliTPCcalibTracks::TestMerge(calTracks, clusterParam);
2523    */
2524    TList *list = new TList();
2525    if (ct == 0 || clusterParam == 0) return 0;
2526    cout << "making list with " << nCalTracks << " AliTPCcalibTrack objects" << endl;
2527    for (Int_t i = 0; i < nCalTracks; i++) {
2528       if (i%10==0) cout << "Adding element " << i << " of " << nCalTracks << endl;
2529       list->Add(new AliTPCcalibTracks(*ct));
2530    }
2531
2532    // only for check at the end
2533    AliTPCcalibTracks* cal1 = new AliTPCcalibTracks(*ct);
2534    Double_t cal1Entries = ((TH1F*)cal1->GetfArrayAmpRow()->At(5))->GetEntries();
2535 //    Double_t cal1Entries = 5; //((TH1F*)ct->GetfArrayAmpRow()->At(5))->GetEntries();
2536
2537    cout  << "The list contains " << list->GetEntries() << " entries. " << endl;
2538
2539
2540    AliTPCcalibTracksCuts *cuts = new AliTPCcalibTracksCuts(20, 0.4, 0.5, 0.13, 0.018);
2541    AliTPCcalibTracks* cal = new AliTPCcalibTracks("calTracksMerged", "calTracksMerged", clusterParam, cuts, 5);
2542    cal->Merge(list);
2543
2544    cout << "cal->GetfArrayAmpRow()->At(5)->Print():" << endl;
2545    cal->GetfArrayAmpRow()->At(5)->Print();
2546    Double_t calEntries = ((TH1F*)cal->GetfArrayAmpRow()->At(5))->GetEntries();
2547
2548    cout << "cal1->GetfArrayAmpRow()->At(5))->GetEntries() = " << cal1Entries << endl;
2549    cout << " cal->GetfArrayAmpRow()->At(5))->GetEntries() = " <<  calEntries << endl;
2550    printf("That means there were %f / %f = %f AliTPCcalibTracks-Objects merged. \n",
2551       calEntries, cal1Entries, ((Double_t)calEntries/cal1Entries));
2552
2553    return cal;
2554
2555 }
2556
2557
2558 void  AliTPCcalibTracks::MakeQPosNormAll(TTree * chainres, AliTPCClusterParam * param, Int_t maxPoints){
2559   //
2560   // Make position corrections
2561   // for the moment Only using debug streamer
2562   // chainres  - debug tree
2563   // param     - parameters to be updated
2564   // maxPoints - maximal number of points using for fit
2565   // verbose   - print info flag
2566   //
2567   // Current implementation - only using debug streamers
2568   //
2569
2570   /*
2571     //Defaults
2572     Int_t maxPoints=100000;
2573   */
2574   /*
2575     Usage:
2576     //0. Load libraries
2577     gSystem->Load("libANALYSIS");
2578     gSystem->Load("libSTAT");
2579     gSystem->Load("libTPCcalib");
2580
2581
2582     //1. Load Parameters to be modified:
2583     //e.g:
2584
2585     AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
2586     AliCDBManager::Instance()->SetRun(0)
2587     AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
2588
2589     //2. Load chain from debug streamers
2590     //
2591     //e.g
2592     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
2593     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
2594     AliXRDPROOFtoolkit tool;
2595     TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
2596     chainres->Lookup();
2597     //3. Do fits and store results
2598     //
2599     AliTPCcalibTracks::MakeQPosNormAll(chainres,param,200000,0) ;
2600     TFile f("paramout.root","recreate");
2601     param->Write("clusterParam");
2602     f.Close();
2603     //4. Verification
2604     TFile f2("paramout.root");
2605     AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
2606     param2->SetInstance(param2);
2607     chainres->Draw("fitZ0:AliTPCClusterParam::SPosCorrection(1,0,Cl.fPad,Cl.fTimeBin,Cl.fZ,Cl.fSigmaY2,Cl.fSigmaZ2,Cl.fMax)","Cl.fDetector<36","",10000); // should be line
2608
2609    */
2610
2611
2612   TStatToolkit toolkit;
2613   Double_t chi2;
2614   TVectorD fitParamY0;
2615   TVectorD fitParamY1;
2616   TVectorD fitParamZ0;
2617   TVectorD fitParamZ1;
2618   TMatrixD covMatrix;
2619   Int_t npoints;
2620
2621   chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
2622   chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
2623   chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
2624   chainres->SetAlias("st","(sin(dt)-dt)");
2625   //
2626   chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
2627   //
2628   //
2629   //
2630   TCut cutA("1");
2631   TString fstringY="";
2632   //
2633   fstringY+="(dp)++";            //1
2634   fstringY+="(dp)*di++";         //2
2635   fstringY+="(sp)++";            //3
2636   fstringY+="(sp)*di++";         //4
2637   TString fstringZ="";
2638   fstringZ+="(dt)++";            //1
2639   fstringZ+="(dt)*di++";         //2
2640   fstringZ+="(st)++";            //3
2641   fstringZ+="(st)*di++";         //4
2642   //
2643   // Z corrections
2644   //
2645   TString *strZ0 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
2646   printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2647   param->PosZcor(0) = (TVectorD*) fitParamZ0.Clone();
2648   //
2649   TString *strZ1 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
2650   //
2651   printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2652   param->PosZcor(1) = (TVectorD*) fitParamZ1.Clone();
2653   param->PosZcor(2) = (TVectorD*) fitParamZ1.Clone();
2654   //
2655   // Y corrections
2656   //
2657   TString *strY0 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
2658   printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2659   param->PosYcor(0) = (TVectorD*) fitParamY0.Clone();
2660
2661
2662   TString *strY1 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
2663   //
2664   printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2665   param->PosYcor(1) = (TVectorD*) fitParamY1.Clone();
2666   param->PosYcor(2) = (TVectorD*) fitParamY1.Clone();
2667   //
2668   //
2669   //
2670   chainres->SetAlias("fitZ0",strZ0->Data());
2671   chainres->SetAlias("fitZ1",strZ1->Data());
2672   chainres->SetAlias("fitY0",strY0->Data());
2673   chainres->SetAlias("fitY1",strY1->Data());
2674   //  chainres->Draw("Cl.fZ-PZ0.fElements[0]","CSigmaY0<0.7&&CSigmaZ0<0.7"+cutA,"",10000);
2675 }
2676
2677
2678
2679 void AliTPCcalibTracks::MakeHistos(){
2680   //
2681   ////make THnSparse
2682   //
2683   //THnSparse  *fHisDeltaY;    // THnSparse - delta Y
2684   //THnSparse  *fHisDeltaZ;    // THnSparse - delta Z
2685   //THnSparse  *fHisRMSY;      // THnSparse - rms Y
2686   //THnSparse  *fHisRMSZ;      // THnSparse - rms Z
2687   //THnSparse  *fHisQmax;      // THnSparse - qmax
2688   //THnSparse  *fHisQtot;      // THnSparse - qtot
2689   // cluster  performance bins
2690   // 0 - variable of interest
2691   // 1 - pad type   - 0- short 1-medium 2-long pads
2692   // 2 - drift length - drift length -0-1
2693   // 3 - Qmax         - Qmax  - 2- 400
2694   // 4 - cog          - COG position - 0-1
2695   // 5 - tan(phi)     - local y angle
2696   // 6 - tan(theta)   - local z angle
2697   // 7 - sector       - sector number
2698   Double_t xminTrack[8], xmaxTrack[8];
2699   Int_t binsTrack[8];
2700   TString axisName[8];
2701
2702   //
2703   binsTrack[0]=100;
2704   axisName[0]  ="var";
2705
2706   binsTrack[1] =3;
2707   xminTrack[1] =0; xmaxTrack[1]=3;
2708   axisName[1]  ="pad type";
2709   //
2710   binsTrack[2] =10;
2711   xminTrack[2] =0; xmaxTrack[2]=1;
2712   axisName[2]  ="drift length";
2713   //
2714   binsTrack[3] =10;
2715   xminTrack[3] =1; xmaxTrack[3]=400;
2716   axisName[3]  ="Qmax";
2717   //
2718   binsTrack[4] =10;
2719   xminTrack[4] =0; xmaxTrack[4]=1;
2720   axisName[4]  ="cog";
2721   //
2722   binsTrack[5] =10;
2723   xminTrack[5] =0; xmaxTrack[5]=2;
2724   axisName[5]  ="tan(phi)";
2725   //
2726   binsTrack[6] =10;
2727   xminTrack[6] =0; xmaxTrack[6]=2;
2728   axisName[6]  ="tan(theta)";
2729   //
2730   xminTrack[0] =-0.5; xmaxTrack[0]=0.5;
2731   fHisDeltaY=new THnSparseS("#Delta_{y} (cm)","#Delta_{y} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
2732   xminTrack[0] =-0.5; xmaxTrack[0]=0.5;
2733   fHisDeltaZ=new THnSparseS("#Delta_{z} (cm)","#Delta_{z} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
2734   xminTrack[0] =0.; xmaxTrack[0]=0.5;
2735   fHisRMSY=new THnSparseS("#RMS_{y} (cm)","#RMS_{y} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
2736   xminTrack[0] =0.; xmaxTrack[0]=0.5;
2737   fHisRMSZ=new THnSparseS("#RMS_{z} (cm)","#RMS_{z} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
2738   xminTrack[0] =0.; xmaxTrack[0]=100;
2739   fHisQmax=new THnSparseS("Qmax (ADC)","Qmax (ADC)", 7, binsTrack,xminTrack, xmaxTrack);
2740
2741   xminTrack[0] =0.; xmaxTrack[0]=250;
2742   fHisQtot=new THnSparseS("Qtot (ADC)","Qtot (ADC)", 7, binsTrack,xminTrack, xmaxTrack);
2743   BinLogX(fHisDeltaY,3);
2744   BinLogX(fHisDeltaZ,3);
2745   BinLogX(fHisRMSY,3);
2746   BinLogX(fHisRMSZ,3);
2747   BinLogX(fHisQmax,3);
2748   BinLogX(fHisQtot,3);
2749
2750 }
2751
2752 void    AliTPCcalibTracks::AddHistos(AliTPCcalibTracks* calib){
2753   //
2754   // Add histograms
2755   //
2756   if (calib->fHisDeltaY) fHisDeltaY->Add(calib->fHisDeltaY);
2757   if (calib->fHisDeltaZ) fHisDeltaZ->Add(calib->fHisDeltaZ);
2758   if (calib->fHisRMSY)   fHisRMSY->Add(calib->fHisRMSY);
2759   if (calib->fHisRMSZ)   fHisRMSZ->Add(calib->fHisRMSZ);
2760 }