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