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