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