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