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