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