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