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