TPC/AliTPCcalibTracks.cxx

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