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