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