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,100,0,32);         // 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,100,0,32);         // 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, 25, 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, 25, 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       Double_t smax = TMath::Sqrt(max);
 893       profAmpRow->Fill( (Double_t)cluster0->GetRow(), smax );
 894       TH1F *hisAmp =  (TH1F*)fArrayAmp->At(sector);
 895       hisAmp->Fill(smax);
 896
 897       // remove the following two lines one day:
 898       TProfile *profDriftLength = (TProfile*)fArrayChargeVsDriftlength->At(sector);
 899       profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax );
 900
 901       TProfile *profDriftLengthTmp = (TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isegment, padSize));
 902       profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax );
 903
 904       fHclusterPerPadrow->Fill(irow);   // fill histogram showing clusters per padrow
 905       Int_t ipad = TMath::Nint(cluster0->GetPad());
 906       Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
 907       fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
 908
 909
 910       TH3F * his3 = 0;
 911       his3 = (TH3F*)fRMSY->At(padSize);
 912       if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
 913       his3 = (TH3F*)fRMSZ->At(padSize);
 914       if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
 915
 916       his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
 917       if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
 918       his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
 919       if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
 920
 921
 922       // Fill resolution histograms
 923       Bool_t useForResol = kTRUE;
 924       if (fFitterParY.GetParError(0) > kErrorFraction * csigmaY) useForResol = kFALSE;
 925
 926       if (cstream){
 927         Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
 928         Float_t sy = cluster0->GetSigmaY2();
 929         Float_t sz = cluster0->GetSigmaZ2();
 930         (*cstream)<<"Resol0"<<
 931           "run="<<fRun<<              //  run number
 932           "event="<<fEvent<<          //  event number
 933           "time="<<fTime<<            //  time stamp of event
 934           "trigger="<<fTrigger<<      //  trigger
 935           "mag="<<fMagF<<             //  magnetic field
 936           "padSize="<<padSize<<
 937           "angley="<<angley<<
 938           "anglez="<<anglez<<
 939           "zdr="<<zdrift<<
 940           "dy="<<deltay<<
 941           "dz="<<deltaz<<
 942           "sy="<<sy<<
 943           "sz="<<sz<<
 944           "\n";
 945       }
 946
 947       if (useForResol){
 948          fDeltaY->Fill(deltay);
 949          fDeltaZ->Fill(deltaz);
 950          his3 = (TH3F*)fResolY->At(padSize);
 951          if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
 952          his3 = (TH3F*)fResolZ->At(padSize);
 953          if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
 954          his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
 955          if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
 956          his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
 957          if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );
 958       }
 959
 960       //=============================================================================================
 961
 962       if (useForResol && nclFound > 2 * kMinRatio * kDelta
 963           && irow % kDeltaWriteDebugStream == 0 && GetDebugLevel() > 4){
 964         if (GetDebugLevel() > 20) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow);
 965          FillResolutionHistoLocalDebugPart(track, cluster0, irow, angley, anglez, nclFound, kDelta);
 966       }  // if (useForResol && nclFound > 2 * kMinRatio * kDelta)
 967
 968    }    // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
 969 }  // FillResolutionHistoLocal(...)
 970
 971
 972
 973 void AliTPCcalibTracks::FillResolutionHistoLocalDebugPart(AliTPCseed *track, AliTPCclusterMI *cluster0, Int_t irow, Float_t  angley, Float_t  anglez, Int_t nclFound, Int_t kDelta) {
 974    //
 975    //  - debug part of FillResolutionHistoLocal -
 976    // called only for every kDeltaWriteDebugStream'th padrow, to avoid to much redundant data
 977    // called only for GetStreamLevel() > 4
 978    // fill resolution trees
 979    //
 980
 981    Int_t sector = cluster0->GetDetector();
 982    Float_t xref = cluster0->GetX();
 983    Int_t padSize = 0;                          // short pads
 984    if (cluster0->GetDetector() >= 36) {
 985       padSize = 1;                              // medium pads
 986       if (cluster0->GetRow() > 63) padSize = 2; // long pads
 987    }
 988
 989    static TLinearFitter fitY0(3, "pol2");
 990    static TLinearFitter fitZ0(3, "pol2");
 991    static TLinearFitter fitY2(5, "hyp4");
 992    static TLinearFitter fitZ2(5, "hyp4");
 993    static TLinearFitter fitY2Q(5, "hyp4");
 994    static TLinearFitter fitZ2Q(5, "hyp4");
 995    static TLinearFitter fitY2S(5, "hyp4");
 996    static TLinearFitter fitZ2S(5, "hyp4");
 997    fitY0.ClearPoints();
 998    fitZ0.ClearPoints();
 999    fitY2.ClearPoints();
1000    fitZ2.ClearPoints();
1001    fitY2Q.ClearPoints();
1002    fitZ2Q.ClearPoints();
1003    fitY2S.ClearPoints();
1004    fitZ2S.ClearPoints();
1005
1006    for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
1007      // loop over irow +- kDelta rows (neighboured rows)
1008      //
1009      //
1010      if (idelta == 0) continue;
1011      if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
1012      AliTPCclusterMI * cluster = track->GetClusterPointer(irow + idelta);
1013      if (!cluster) continue;
1014      if (cluster->GetType() < 0) continue;
1015      if (cluster->GetDetector() != sector) continue;
1016      Double_t x = cluster->GetX() - xref;
1017      Double_t sigmaY0 = fClusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley) );
1018      Double_t sigmaZ0 = fClusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez) );
1019      //
1020      Double_t sigmaYQ = fClusterParam->GetErrorQPar( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
1021      Double_t sigmaZQ = fClusterParam->GetErrorQPar( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
1022      Double_t sigmaYS = fClusterParam->GetErrorQParScaled( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())),
1023                                                            TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
1024      Double_t sigmaZS = fClusterParam->GetErrorQParScaled( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())),
1025                                                            TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
1026      Float_t rmsYFactor = fClusterParam->GetShapeFactor( 0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
1027                                                          TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
1028                                                          TMath::Sqrt(cluster0->GetSigmaY2()), 0 );
1029      Float_t rmsZFactor = fClusterParam->GetShapeFactor(0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
1030                                                         TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
1031                                                         TMath::Sqrt(cluster0->GetSigmaZ2()),0 );
1032      sigmaYS  = TMath::Sqrt(sigmaYS * sigmaYS + rmsYFactor * rmsYFactor / 12.);
1033      sigmaZS  = TMath::Sqrt(sigmaZS * sigmaZS + rmsZFactor * rmsZFactor / 12. + rmsYFactor * rmsYFactor / 24.);
1034      //
1035      if (kDelta != 0){
1036        fitY0.AddPoint(&x, cluster->GetY(), sigmaY0);
1037        fitZ0.AddPoint(&x, cluster->GetZ(), sigmaZ0);
1038      }
1039      Double_t xxx[4];
1040      xxx[0] = ( (idelta+irow) % 2 == 0 ) ? 1 : 0;
1041      xxx[1] = x;
1042      xxx[2] = ( (idelta+irow) % 2 == 0 ) ? x : 0;
1043      xxx[3] = x * x;
1044      fitY2.AddPoint(xxx, cluster->GetY(), sigmaY0);
1045      fitY2Q.AddPoint(xxx, cluster->GetY(), sigmaYQ);
1046      fitY2S.AddPoint(xxx, cluster->GetY(), sigmaYS);
1047      fitZ2.AddPoint(xxx, cluster->GetZ(), sigmaZ0);
1048      fitZ2Q.AddPoint(xxx, cluster->GetZ(), sigmaZQ);
1049      fitZ2S.AddPoint(xxx, cluster->GetZ(), sigmaZS);
1050      //
1051    }  // neigbouhood-loop
1052       //
1053    fitY0.Eval();
1054    fitZ0.Eval();
1055    fitY2.Eval();
1056    fitZ2.Eval();
1057    fitY2Q.Eval();
1058    fitZ2Q.Eval();
1059    fitY2S.Eval();
1060    fitZ2S.Eval();
1061    Float_t chi2Y0 = fitY0.GetChisquare() / (nclFound-3.);
1062    Float_t chi2Z0 = fitZ0.GetChisquare() / (nclFound-3.);
1063    Float_t chi2Y2 = fitY2.GetChisquare() / (nclFound-5.);
1064    Float_t chi2Z2 = fitZ2.GetChisquare() / (nclFound-5.);
1065    Float_t chi2Y2Q = fitY2Q.GetChisquare() / (nclFound-5.);
1066    Float_t chi2Z2Q = fitZ2Q.GetChisquare() / (nclFound-5.);
1067    Float_t chi2Y2S = fitY2S.GetChisquare() / (nclFound-5.);
1068    Float_t chi2Z2S = fitZ2S.GetChisquare() / (nclFound-5.);
1069    //
1070    static  TVectorD    parY0(3);
1071    static  TMatrixD    matY0(3, 3);
1072    static  TVectorD    parZ0(3);
1073    static  TMatrixD    matZ0(3, 3);
1074    fitY0.GetParameters(parY0);
1075    fitY0.GetCovarianceMatrix(matY0);
1076    fitZ0.GetParameters(parZ0);
1077    fitZ0.GetCovarianceMatrix(matZ0);
1078    //
1079    static  TVectorD    parY2(5);
1080    static  TMatrixD    matY2(5,5);
1081    static  TVectorD    parZ2(5);
1082    static  TMatrixD    matZ2(5,5);
1083    fitY2.GetParameters(parY2);
1084    fitY2.GetCovarianceMatrix(matY2);
1085    fitZ2.GetParameters(parZ2);
1086    fitZ2.GetCovarianceMatrix(matZ2);
1087    //
1088    static  TVectorD    parY2Q(5);
1089    static  TMatrixD    matY2Q(5,5);
1090    static  TVectorD    parZ2Q(5);
1091    static  TMatrixD    matZ2Q(5,5);
1092    fitY2Q.GetParameters(parY2Q);
1093    fitY2Q.GetCovarianceMatrix(matY2Q);
1094    fitZ2Q.GetParameters(parZ2Q);
1095    fitZ2Q.GetCovarianceMatrix(matZ2Q);
1096    static  TVectorD    parY2S(5);
1097    static  TMatrixD    matY2S(5,5);
1098    static  TVectorD    parZ2S(5);
1099    static  TMatrixD    matZ2S(5,5);
1100    fitY2S.GetParameters(parY2S);
1101    fitY2S.GetCovarianceMatrix(matY2S);
1102    fitZ2S.GetParameters(parZ2S);
1103    fitZ2S.GetCovarianceMatrix(matZ2S);
1104    Float_t sigmaY0   = TMath::Sqrt(matY0(0,0));
1105    Float_t sigmaZ0   = TMath::Sqrt(matZ0(0,0));
1106    Float_t sigmaDY0  = TMath::Sqrt(matY0(1,1));
1107    Float_t sigmaDZ0  = TMath::Sqrt(matZ0(1,1));
1108    Float_t sigmaY2   = TMath::Sqrt(matY2(1,1));
1109    Float_t sigmaZ2   = TMath::Sqrt(matZ2(1,1));
1110    Float_t sigmaDY2  = TMath::Sqrt(matY2(3,3));
1111    Float_t sigmaDZ2  = TMath::Sqrt(matZ2(3,3));
1112    Float_t sigmaY2Q  = TMath::Sqrt(matY2Q(1,1));
1113    Float_t sigmaZ2Q  = TMath::Sqrt(matZ2Q(1,1));
1114    Float_t sigmaDY2Q = TMath::Sqrt(matY2Q(3,3));
1115    Float_t sigmaDZ2Q = TMath::Sqrt(matZ2Q(3,3));
1116    Float_t sigmaY2S  = TMath::Sqrt(matY2S(1,1));
1117    Float_t sigmaZ2S  = TMath::Sqrt(matZ2S(1,1));
1118    Float_t sigmaDY2S = TMath::Sqrt(matY2S(3,3));
1119    Float_t sigmaDZ2S = TMath::Sqrt(matZ2S(3,3));
1120
1121    // Error parameters
1122    Float_t csigmaY0 = fClusterParam->GetError0Par(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(angley));
1123    Float_t csigmaZ0 = fClusterParam->GetError0Par(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(anglez));
1124    //
1125    Float_t csigmaYQ = fClusterParam->GetErrorQPar(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1126                                                   TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1127    Float_t csigmaZQ = fClusterParam->GetErrorQPar(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1128                                                   TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
1129    Float_t csigmaYS = fClusterParam->GetErrorQParScaled(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1130                                                         TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1131    Float_t csigmaZS = fClusterParam->GetErrorQParScaled(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1132                                                         TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
1133
1134
1135    // RMS parameters
1136    Float_t meanRMSY = 0;
1137    Float_t meanRMSZ = 0;
1138    Int_t   nclRMS = 0;
1139    for (Int_t idelta = -2; idelta <= 2; idelta++){
1140      // loop over neighbourhood
1141      if (idelta+irow < 0 || idelta+irow > 159) continue;
1142      //         if (idelta+irow>159) continue;
1143      AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta);
1144      if (!cluster) continue;
1145      meanRMSY += TMath::Sqrt(cluster->GetSigmaY2());
1146      meanRMSZ += TMath::Sqrt(cluster->GetSigmaZ2());
1147      nclRMS++;
1148    }
1149    meanRMSY /= nclRMS;
1150    meanRMSZ /= nclRMS;
1151
1152    Float_t rmsY      = TMath::Sqrt(cluster0->GetSigmaY2());
1153    Float_t rmsZ      = TMath::Sqrt(cluster0->GetSigmaZ2());
1154    Float_t rmsYT     = fClusterParam->GetRMSQ(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1155                                               TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1156    Float_t rmsZT     = fClusterParam->GetRMSQ(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1157                                               TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1158    Float_t rmsYT0    = fClusterParam->GetRMS0(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1159                                               TMath::Abs(angley));
1160    Float_t rmsZT0    = fClusterParam->GetRMS0(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1161                                                  TMath::Abs(anglez));
1162    Float_t rmsYSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1163                                                   TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1164    Float_t rmsZSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1165                                                   TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1166    Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1167                                                       TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
1168                                                       rmsY,meanRMSY);
1169    Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1170                                                       TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
1171                                                       rmsZ,meanRMSZ);
1172    //
1173    // cluster debug
1174    TTreeSRedirector *cstream = GetDebugStreamer();
1175    if (cstream){
1176      (*cstream)<<"ResolCl"<<    // valgrind 3   40,000 bytes in 1 blocks are possibly
1177        "run="<<fRun<<              //  run number
1178        "event="<<fEvent<<          //  event number
1179        "time="<<fTime<<            //  time stamp of event
1180        "trigger="<<fTrigger<<      //  trigger
1181        "mag="<<fMagF<<             //  magnetic field
1182        "Sector="<<sector<<
1183        "Cl.="<<cluster0<<
1184        "CSigmaY0="<<csigmaY0<<   // cluster errorY
1185        "CSigmaYQ="<<csigmaYQ<<   // cluster errorY - q scaled
1186        "CSigmaYS="<<csigmaYS<<   // cluster errorY - q scaled
1187        "CSigmaZ0="<<csigmaZ0<<   //
1188        "CSigmaZQ="<<csigmaZQ<<
1189        "CSigmaZS="<<csigmaZS<<
1190        "shapeYF="<<rmsYFactor<<
1191        "shapeZF="<<rmsZFactor<<
1192        "rmsY="<<rmsY<<
1193        "rmsZ="<<rmsZ<<
1194        "rmsYM="<<meanRMSY<<
1195        "rmsZM="<<meanRMSZ<<
1196        "rmsYT="<<rmsYT<<
1197        "rmsZT="<<rmsZT<<
1198        "rmsYT0="<<rmsYT0<<
1199        "rmsZT0="<<rmsZT0<<
1200        "rmsYS="<<rmsYSigma<<
1201        "rmsZS="<<rmsZSigma<<
1202        "padSize="<<padSize<<
1203        "Ncl="<<nclFound<<
1204        "PY0.="<<&parY0<<
1205        "PZ0.="<<&parZ0<<
1206        "SigmaY0="<<sigmaY0<<
1207        "SigmaZ0="<<sigmaZ0<<
1208        "angley="<<angley<<
1209        "anglez="<<anglez<<
1210        "\n";
1211
1212      //       tracklet dubug
1213      (*cstream)<<"ResolTr"<<
1214        "run="<<fRun<<              //  run number
1215        "event="<<fEvent<<          //  event number
1216        "time="<<fTime<<            //  time stamp of event
1217        "trigger="<<fTrigger<<      //  trigger
1218        "mag="<<fMagF<<             //  magnetic field
1219        "padSize="<<padSize<<
1220        "IPad="<<padSize<<
1221        "Sector="<<sector<<
1222        "Ncl="<<nclFound<<
1223        "chi2Y0="<<chi2Y0<<
1224        "chi2Z0="<<chi2Z0<<
1225        "chi2Y2="<<chi2Y2<<
1226        "chi2Z2="<<chi2Z2<<
1227        "chi2Y2Q="<<chi2Y2Q<<
1228        "chi2Z2Q="<<chi2Z2Q<<
1229        "chi2Y2S="<<chi2Y2S<<
1230        "chi2Z2S="<<chi2Z2S<<
1231        "PY0.="<<&parY0<<
1232        "PZ0.="<<&parZ0<<
1233        "PY2.="<<&parY2<<
1234        "PZ2.="<<&parZ2<<
1235        "PY2Q.="<<&parY2Q<<
1236        "PZ2Q.="<<&parZ2Q<<
1237        "PY2S.="<<&parY2S<<
1238        "PZ2S.="<<&parZ2S<<
1239        "SigmaY0="<<sigmaY0<<
1240        "SigmaZ0="<<sigmaZ0<<
1241        "SigmaDY0="<<sigmaDY0<<
1242        "SigmaDZ0="<<sigmaDZ0<<
1243        "SigmaY2="<<sigmaY2<<
1244        "SigmaZ2="<<sigmaZ2<<
1245        "SigmaDY2="<<sigmaDY2<<
1246        "SigmaDZ2="<<sigmaDZ2<<
1247        "SigmaY2Q="<<sigmaY2Q<<
1248        "SigmaZ2Q="<<sigmaZ2Q<<
1249        "SigmaDY2Q="<<sigmaDY2Q<<
1250        "SigmaDZ2Q="<<sigmaDZ2Q<<
1251        "SigmaY2S="<<sigmaY2S<<
1252        "SigmaZ2S="<<sigmaZ2S<<
1253        "SigmaDY2S="<<sigmaDY2S<<
1254        "SigmaDZ2S="<<sigmaDZ2S<<
1255        "angley="<<angley<<
1256        "anglez="<<anglez<<
1257        "\n";
1258    }
1259 }
1260
1261
1262
1263
1264
1265 TH2D * AliTPCcalibTracks::MakeDiff(TH2D * hfit, TF2 * func){
1266    //
1267    // creates a new histogram which contains the difference between
1268    // the histogram hfit and the function func
1269    //
1270   TH2D * result = (TH2D*)hfit->Clone();      // valgrind 3   40,139 bytes in 11 blocks are still reachable
1271   result->SetTitle(Form("%s fit residuals",result->GetTitle()));
1272   result->SetName(Form("%s fit residuals",result->GetName()));
1273   TAxis *xaxis = hfit->GetXaxis();
1274   TAxis *yaxis = hfit->GetYaxis();
1275   Double_t x[2];
1276   for (Int_t biny = 0; biny <= yaxis->GetNbins(); biny++) {
1277     x[1]  = yaxis->GetBinCenter(biny);
1278     for (Int_t binx = 0; binx <= xaxis->GetNbins(); binx++) {
1279       x[0]  = xaxis->GetBinCenter(binx);
1280       Int_t bin = hfit->GetBin(binx, biny);
1281       Double_t val = hfit->GetBinContent(bin);
1282 //      result->SetBinContent( bin, (val - func->Eval(x[0], x[1])) / func->Eval(x[0], x[1]) );
1283       result->SetBinContent( bin, (val / func->Eval(x[0], x[1])) - 1 );
1284     }
1285   }
1286   return result;
1287 }
1288
1289
1290 void  AliTPCcalibTracks::SetStyle() const {
1291    //
1292    // set style, can be called by all draw functions
1293    //
1294    gROOT->SetStyle("Plain");
1295    gStyle->SetFillColor(10);
1296    gStyle->SetPadColor(10);
1297    gStyle->SetCanvasColor(10);
1298    gStyle->SetStatColor(10);
1299    gStyle->SetPalette(1,0);
1300    gStyle->SetNumberContours(60);
1301 }
1302
1303
1304 void AliTPCcalibTracks::Draw(Option_t* opt){
1305    //
1306    // draw-function of AliTPCcalibTracks
1307    // will draws some exemplaric pictures
1308    //
1309
1310   if (GetDebugLevel() > 6) Info("Draw", "Drawing an exemplaric picture.");
1311    SetStyle();
1312    Double_t min = 0;
1313    Double_t max = 0;
1314    TCanvas *c1 = new TCanvas();
1315    c1->Divide(0, 3);
1316    TVirtualPad *upperThird = c1->GetPad(1);
1317    TVirtualPad *middleThird = c1->GetPad(2);
1318    TVirtualPad *lowerThird = c1->GetPad(3);
1319    upperThird->Divide(2,0);
1320    TVirtualPad *upleft  = upperThird->GetPad(1);
1321    TVirtualPad *upright = upperThird->GetPad(2);
1322    middleThird->Divide(2,0);
1323    TVirtualPad *middleLeft  = middleThird->GetPad(1);
1324    TVirtualPad *middleRight = middleThird->GetPad(2);
1325    lowerThird->Divide(2,0);
1326    TVirtualPad *downLeft  = lowerThird->GetPad(1);
1327    TVirtualPad *downRight = lowerThird->GetPad(2);
1328
1329
1330    upleft->cd(0);
1331    min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
1332    max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
1333    fDeltaY->SetAxisRange(min, max);
1334    fDeltaY->Fit("gaus","q","",min, max);        // valgrind 3  7 block possibly lost   2,400 bytes in 1 blocks are still reachable
1335    c1->Update();
1336
1337    upright->cd(0);
1338    max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
1339    min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
1340    fDeltaZ->SetAxisRange(min, max);
1341    fDeltaZ->Fit("gaus","q","",min, max);
1342    c1->Update();
1343
1344    middleLeft->cd();
1345    fHclus->Draw(opt);
1346
1347    middleRight->cd();
1348    fRejectedTracksHisto->Draw(opt);
1349    TPaveText *pt = new TPaveText(0.6,0.6, 0.8,0.8, "NDC");
1350    TText *t1 = pt->AddText("1: kEdgeThetaCutNoise");
1351    TText *t2 = pt->AddText("2: kMinClusters");
1352    TText *t3 = pt->AddText("3: kMinRatio");
1353    TText *t4 = pt->AddText("4: kMax1pt");
1354    t1 = t1; t2 = t2; t3 = t3; t4 = t4;    // avoid compiler warnings
1355    pt->SetToolTipText("Legend for failed cuts");
1356    pt->Draw();
1357
1358    downLeft->cd();
1359    fHclusterPerPadrowRaw->Draw(opt);
1360
1361    downRight->cd();
1362    fHclusterPerPadrow->Draw(opt);
1363 }
1364
1365
1366 void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){
1367    //
1368    // all functions are called, that produce pictures
1369    // the histograms are written to the directory 'pathName'
1370    // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file
1371    // 'stat' is also the number of minEntries for MakeResPlotsQTree
1372    //
1373
1374   if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
1375    MakeAmpPlots(stat, pathName);
1376    MakeDeltaPlots(pathName);
1377    FitResolutionNew(pathName);
1378    FitRMSNew(pathName);
1379    MakeChargeVsDriftLengthPlots(pathName);
1380 //    MakeResPlotsQ(1, 1);
1381    MakeResPlotsQTree(stat, pathName);
1382 }
1383
1384
1385 void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, const char* pathName){
1386    //
1387    // creates several plots:
1388    // fArrayAmp.ps, fArrayAmpRow.ps and DeltaYZ.ps
1389    // fArrayAmp.ps: one histogram per sector, the histogram shows the charge per cluster
1390    // fArrayAmpRow.ps: one histogram per sector, mean max. amplitude vs. pad row with landau fit
1391    // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
1392    // Empty histograms (sectors without data) are not written to file
1393    // the ps-files are written to the directory 'pathName', that is created if it does not exist
1394    // 'stat': only histograms with more than 'stat' entries are written to file.
1395    //
1396
1397    SetStyle();
1398    gSystem->MakeDirectory(pathName);
1399    gSystem->ChangeDirectory(pathName);
1400
1401    TCanvas* c1 = new TCanvas();     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
1402    TPostScript *ps;
1403    // histograms with accumulated amplitude for all IROCs and OROCs
1404    TH1F *allAmpHisIROC = ((TH1F*)(fArrayAmp->At(0))->Clone());
1405    allAmpHisIROC->SetName("Amp all IROCs");
1406    allAmpHisIROC->SetTitle("Amp all IROCs");
1407    TH1F *allAmpHisOROC = ((TH1F*)(fArrayAmp->At(36))->Clone());
1408    allAmpHisOROC->SetName("Amp all OROCs");
1409    allAmpHisOROC->SetTitle("Amp all OROCs");
1410
1411    ps = new TPostScript("fArrayAmp.ps", 112);
1412    if (GetDebugLevel() > -1) cout << "creating fArrayAmp.ps..." << endl;
1413    for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++){
1414       if ( ((TH1F*)fArrayAmp->At(i))->GetEntries() < stat  ) continue;
1415       ps->NewPage();
1416       ((TH1F*)fArrayAmp->At(i))->Draw();
1417       c1->Update();              // valgrind 3
1418       if (i > 0 && i < 36) {
1419          allAmpHisIROC->Add(((TH1F*)fArrayAmp->At(i)));
1420          allAmpHisOROC->Add(((TH1F*)fArrayAmp->At(i+36)));
1421       }
1422    }
1423    ps->NewPage();
1424    allAmpHisIROC->Draw();
1425    c1->Update();              // valgrind
1426    ps->NewPage();
1427    allAmpHisOROC->Draw();
1428    c1->Update();
1429    ps->Close();
1430    delete ps;
1431
1432    TH1F *his = 0;
1433    Double_t min = 0;
1434    Double_t max = 0;
1435    ps = new TPostScript("fArrayAmpRow.ps", 112);
1436    if (GetDebugLevel() > -1) cout << "creating fArrayAmpRow.ps..." << endl;
1437    for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++){
1438       his = (TH1F*)fArrayAmpRow->At(i);
1439       if (his->GetEntries() < stat) continue;
1440       ps->NewPage();
1441       min = TMath::Max( his->GetBinCenter(his->GetMaximumBin() )-100., 0.);
1442       max = his->GetBinCenter(5*his->GetMaximumBin()) + 100;
1443       his->SetAxisRange(min, max);
1444       his->Fit("pol3", "q", "", min, max);
1445       // his->Draw("error");    // don't use this line when you don't want to have empty pages in the ps-file
1446       c1->Update();
1447    }
1448    ps->Close();
1449    delete ps;
1450    delete c1;
1451    gSystem->ChangeDirectory("..");
1452 }
1453
1454
1455 void AliTPCcalibTracks::MakeDeltaPlots(const char* pathName){
1456    //
1457    // creates several plots:
1458    // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
1459    // the ps-files are written to the directory 'pathName', that is created if it does not exist
1460    //
1461
1462    SetStyle();
1463    gSystem->MakeDirectory(pathName);
1464    gSystem->ChangeDirectory(pathName);
1465
1466    TCanvas* c1 = new TCanvas();     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
1467    TPostScript *ps;
1468    Double_t min = 0;
1469    Double_t max = 0;
1470
1471    ps = new TPostScript("DeltaYZ.ps", 112);
1472    if (GetDebugLevel() > -1) cout << "creating DeltaYZ.ps..." << endl;
1473    min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
1474    max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
1475    fDeltaY->SetAxisRange(min, max);
1476    ps->NewPage();
1477    fDeltaY->Fit("gaus","q","",min, max);        // valgrind 3  7 block possibly lost   2,400 bytes in 1 blocks are still reachable
1478    c1->Update();
1479    ps->NewPage();
1480    max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
1481    min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
1482    fDeltaZ->SetAxisRange(min, max);
1483    fDeltaZ->Fit("gaus","q","",min, max);
1484    c1->Update();
1485    ps->Close();
1486    delete ps;
1487    delete c1;
1488    gSystem->ChangeDirectory("..");
1489 }
1490
1491
1492 void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(const char* pathName){
1493    //
1494    // creates charge vs. driftlength plots, one TProfile for each ROC
1495    // is not correct like this, should be one TProfile for each sector and padsize
1496    //
1497
1498    SetStyle();
1499    gSystem->MakeDirectory(pathName);
1500    gSystem->ChangeDirectory(pathName);
1501
1502    TCanvas* c1 = new TCanvas();     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
1503    TPostScript *ps;
1504    ps = new TPostScript("chargeVsDriftlengthOld.ps", 112);
1505    if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlength.ps..." << endl;
1506    TProfile *chargeVsDriftlengthAllIROCs = ((TProfile*)fArrayChargeVsDriftlength->At(0)->Clone());
1507    TProfile *chargeVsDriftlengthAllOROCs = ((TProfile*)fArrayChargeVsDriftlength->At(36)->Clone());
1508    chargeVsDriftlengthAllIROCs->SetName("allAmpHisIROC");
1509    chargeVsDriftlengthAllIROCs->SetTitle("charge vs. driftlength, all IROCs");
1510    chargeVsDriftlengthAllOROCs->SetName("allAmpHisOROC");
1511    chargeVsDriftlengthAllOROCs->SetTitle("charge vs. driftlength, all OROCs");
1512
1513    for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) {
1514       ((TProfile*)fArrayChargeVsDriftlength->At(i))->Draw();
1515       c1->Update();
1516       if (i > 0 && i < 36) {
1517          chargeVsDriftlengthAllIROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i)));
1518          chargeVsDriftlengthAllOROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i+36)));
1519       }
1520       ps->NewPage();
1521    }
1522    chargeVsDriftlengthAllIROCs->Draw();
1523    c1->Update();              // valgrind
1524    ps->NewPage();
1525    chargeVsDriftlengthAllOROCs->Draw();
1526    c1->Update();
1527    ps->Close();
1528    delete ps;
1529    delete c1;
1530    gSystem->ChangeDirectory("..");
1531 }
1532
1533
1534 void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(const char* pathName){
1535    //
1536    // creates charge vs. driftlength plots, one TProfile for each ROC
1537    // under development....
1538    //
1539
1540    SetStyle();
1541    gSystem->MakeDirectory(pathName);
1542    gSystem->ChangeDirectory(pathName);
1543
1544    TCanvas* c1 = new TCanvas("c1", "c1", 700,(Int_t)(TMath::Sqrt(2)*700));     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
1545 //    TCanvas c1("c1", "c1", 500,(sqrt(2)*500))
1546    c1->Divide(0,3);
1547    TPostScript *ps;
1548    ps = new TPostScript("chargeVsDriftlength.ps", 111);
1549    if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlengthNew.ps..." << endl;
1550
1551    TProfile *chargeVsDriftlengthAllShortPads  = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0)->Clone());
1552    TProfile *chargeVsDriftlengthAllMediumPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1)->Clone());
1553    TProfile *chargeVsDriftlengthAllLongPads   = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2)->Clone());
1554    chargeVsDriftlengthAllShortPads->SetName("allAmpHisShortPads");
1555    chargeVsDriftlengthAllShortPads->SetTitle("charge vs. driftlength, all sectors, short pads");
1556    chargeVsDriftlengthAllMediumPads->SetName("allAmpHisMediumPads");
1557    chargeVsDriftlengthAllMediumPads->SetTitle("charge vs. driftlength, all sectors, medium pads");
1558    chargeVsDriftlengthAllLongPads->SetName("allAmpHisLongPads");
1559    chargeVsDriftlengthAllLongPads->SetTitle("charge vs. driftlength, all sectors, long pads");
1560
1561    for (Int_t i = 0; i < 36; i++) {
1562       c1->cd(1)->SetGridx();
1563       c1->cd(1)->SetGridy();
1564       ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,0))->Draw();
1565       c1->cd(2)->SetGridx();
1566       c1->cd(2)->SetGridy();
1567       ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,1))->Draw();
1568       c1->cd(3)->SetGridx();
1569       c1->cd(3)->SetGridy();
1570       ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,2))->Draw();
1571       c1->Update();
1572       chargeVsDriftlengthAllShortPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0));
1573       chargeVsDriftlengthAllMediumPads->Add((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1));
1574       chargeVsDriftlengthAllLongPads->Add(  (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2));
1575       ps->NewPage();
1576    }
1577    c1->cd(1)->SetGridx();
1578    c1->cd(1)->SetGridy();
1579    chargeVsDriftlengthAllShortPads->Draw();
1580    c1->cd(2)->SetGridx();
1581    c1->cd(2)->SetGridy();
1582    chargeVsDriftlengthAllMediumPads->Draw();
1583    c1->cd(3)->SetGridx();
1584    c1->cd(3)->SetGridy();
1585    chargeVsDriftlengthAllLongPads->Draw();
1586    c1->Update();              // valgrind
1587 //    ps->NewPage();
1588    ps->Close();
1589    delete ps;
1590    delete c1;
1591    gSystem->ChangeDirectory("..");
1592 }
1593
1594
1595
1596 void AliTPCcalibTracks::FitResolutionNew(const char* pathName){
1597    //
1598    // calculates different resulution fits in Y and Z direction
1599    // the histograms are written to 'ResolutionYZ.ps'
1600    // writes calculated resolution to 'resol.txt'
1601    // all files are stored in the directory pathName
1602    //
1603
1604    SetStyle();
1605    gSystem->MakeDirectory(pathName);
1606    gSystem->ChangeDirectory(pathName);
1607
1608    TCanvas c;
1609    c.Divide(2,1);
1610    if (GetDebugLevel() > -1) cout << "creating ResolutionYZ.ps..." << endl;
1611    TPostScript *ps = new TPostScript("ResolutionYZ.ps", 112);
1612    TF2 *fres = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
1613    fres->SetParameter(0,0.02);
1614    fres->SetParameter(1,0.0054);
1615    fres->SetParameter(2,0.13);
1616
1617    TH1::AddDirectory(kTRUE);  // TH3F::FitSlicesZ() writes histograms into the current directory
1618
1619    // create histogramw for Y-resolution
1620    TH3F * hisResY0 = (TH3F*)fResolY->At(0);
1621    hisResY0->FitSlicesZ();
1622    TH2D * hisResY02 = (TH2D*)gDirectory->Get("Resol Y0_2");
1623    TH3F * hisResY1 = (TH3F*)fResolY->At(1);
1624    hisResY1->FitSlicesZ();
1625    TH2D * hisResY12 = (TH2D*)gDirectory->Get("Resol Y1_2");
1626    TH3F * hisResY2 = (TH3F*)fResolY->At(2);
1627    hisResY2->FitSlicesZ();
1628    TH2D * hisResY22 = (TH2D*)gDirectory->Get("Resol Y2_2");
1629     //
1630    ps->NewPage();
1631    c.cd(1);
1632    hisResY02->Fit(fres, "q");      // valgrind    132,072 bytes in 6 blocks are indirectly lost
1633    hisResY02->Draw("surf1");
1634    c.cd(2);
1635    MakeDiff(hisResY02,fres)->Draw("surf1");
1636    c.Update();
1637    //   c.SaveAs("ResolutionYPad0.eps");
1638    ps->NewPage();
1639    c.cd(1);
1640    hisResY12->Fit(fres, "q");
1641    hisResY12->Draw("surf1");
1642    c.cd(2);
1643    MakeDiff(hisResY12,fres)->Draw("surf1");
1644    c.Update();
1645    //   c.SaveAs("ResolutionYPad1.eps");
1646    ps->NewPage();
1647    c.cd(1);
1648    hisResY22->Fit(fres, "q");
1649    hisResY22->Draw("surf1");
1650    c.cd(2);
1651    MakeDiff(hisResY22,fres)->Draw("surf1");
1652    c.Update();
1653 //    c.SaveAs("ResolutionYPad2.eps");
1654
1655    // create histogramw for Z-resolution
1656    TH3F * hisResZ0 = (TH3F*)fResolZ->At(0);
1657    hisResZ0->FitSlicesZ();
1658    TH2D * hisResZ02 = (TH2D*)gDirectory->Get("Resol Z0_2");
1659    TH3F * hisResZ1 = (TH3F*)fResolZ->At(1);
1660    hisResZ1->FitSlicesZ();
1661    TH2D * hisResZ12 = (TH2D*)gDirectory->Get("Resol Z1_2");
1662    TH3F * hisResZ2 = (TH3F*)fResolZ->At(2);
1663    hisResZ2->FitSlicesZ();
1664    TH2D * hisResZ22 = (TH2D*)gDirectory->Get("Resol Z2_2");
1665
1666    ps->NewPage();
1667    c.cd(1);
1668    hisResZ02->Fit(fres, "q");
1669    hisResZ02->Draw("surf1");
1670    c.cd(2);
1671    MakeDiff(hisResZ02,fres)->Draw("surf1");
1672    c.Update();
1673 //    c.SaveAs("ResolutionZPad0.eps");
1674    ps->NewPage();
1675    c.cd(1);
1676    hisResZ12->Fit(fres, "q");
1677    hisResZ12->Draw("surf1");
1678    c.cd(2);
1679    MakeDiff(hisResZ12,fres)->Draw("surf1");
1680    c.Update();
1681 //    c.SaveAs("ResolutionZPad1.eps");
1682    ps->NewPage();
1683    c.cd(1);
1684    hisResZ22->Fit(fres, "q");
1685    hisResZ22->Draw("surf1");
1686    c.cd(2);
1687    MakeDiff(hisResZ22,fres)->Draw("surf1");
1688    c.Update();
1689 //    c.SaveAs("ResolutionZPad2.eps");
1690    ps->Close();
1691    delete ps;
1692
1693    // write calculated resoltuions to 'resol.txt'
1694    ofstream fresol("resol.txt");
1695    fresol<<"Pad 0.75 cm"<<"\n";
1696    hisResY02->Fit(fres, "q");                     // valgrind
1697    fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1698    hisResZ02->Fit(fres, "q");
1699    fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1700    //
1701    fresol<<"Pad 1.00 cm"<<1<<"\n";
1702    hisResY12->Fit(fres, "q");                     // valgrind
1703    fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1704    hisResZ12->Fit(fres, "q");
1705    fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1706    //
1707    fresol<<"Pad 1.50 cm"<<0<<"\n";
1708    hisResY22->Fit(fres, "q");
1709    fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1710    hisResZ22->Fit(fres, "q");
1711    fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1712
1713    TH1::AddDirectory(kFALSE);
1714    gSystem->ChangeDirectory("..");
1715    delete fres;
1716 }
1717
1718
1719 void AliTPCcalibTracks::FitRMSNew(const char* pathName){
1720    //
1721    // calculates different resulution-rms fits in Y and Z direction
1722    // the histograms are written to 'RMS_YZ.ps'
1723    // writes calculated resolution rms to 'rms.txt'
1724    // all files are stored in the directory pathName
1725    //
1726
1727    SetStyle();
1728    gSystem->MakeDirectory(pathName);
1729    gSystem->ChangeDirectory(pathName);
1730
1731    TCanvas c;        // valgrind 3   42,120 bytes in 405 blocks are still reachable   23,816 bytes in 229 blocks are still reachable
1732    c.Divide(2,1);
1733    if (GetDebugLevel() > -1) cout << "creating RMS_YZ.ps..." << endl;
1734    TPostScript *ps = new TPostScript("RMS_YZ.ps", 112);
1735    TF2 *frms = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
1736    frms->SetParameter(0,0.02);
1737    frms->SetParameter(1,0.0054);
1738    frms->SetParameter(2,0.13);
1739
1740    TH1::AddDirectory(kTRUE);  // TH3F::FitSlicesZ() writes histograms into the current directory
1741
1742    // create histogramw for Y-RMS
1743    TH3F * hisResY0 = (TH3F*)fRMSY->At(0);
1744    hisResY0->FitSlicesZ();
1745    TH2D * hisResY02 = (TH2D*)gDirectory->Get("RMS Y0_1");
1746    TH3F * hisResY1 = (TH3F*)fRMSY->At(1);
1747    hisResY1->FitSlicesZ();
1748    TH2D * hisResY12 = (TH2D*)gDirectory->Get("RMS Y1_1");
1749    TH3F * hisResY2 = (TH3F*)fRMSY->At(2);
1750    hisResY2->FitSlicesZ();
1751    TH2D * hisResY22 = (TH2D*)gDirectory->Get("RMS Y2_1");
1752    //
1753    ps->NewPage();
1754    c.cd(1);
1755    hisResY02->Fit(frms, "qn0");
1756    hisResY02->Draw("surf1");
1757    c.cd(2);
1758    MakeDiff(hisResY02,frms)->Draw("surf1");
1759    c.Update();
1760 //    c.SaveAs("RMSYPad0.eps");
1761    ps->NewPage();
1762    c.cd(1);
1763    hisResY12->Fit(frms, "qn0");               // valgrind   several blocks possibly lost
1764    hisResY12->Draw("surf1");
1765    c.cd(2);
1766    MakeDiff(hisResY12,frms)->Draw("surf1");
1767    c.Update();
1768 //    c.SaveAs("RMSYPad1.eps");
1769    ps->NewPage();
1770    c.cd(1);
1771    hisResY22->Fit(frms, "qn0");
1772    hisResY22->Draw("surf1");
1773    c.cd(2);
1774    MakeDiff(hisResY22,frms)->Draw("surf1");
1775    c.Update();
1776 //    c.SaveAs("RMSYPad2.eps");
1777
1778    // create histogramw for Z-RMS
1779    TH3F * hisResZ0 = (TH3F*)fRMSZ->At(0);
1780    hisResZ0->FitSlicesZ();
1781    TH2D * hisResZ02 = (TH2D*)gDirectory->Get("RMS Z0_1");
1782    TH3F * hisResZ1 = (TH3F*)fRMSZ->At(1);
1783    hisResZ1->FitSlicesZ();
1784    TH2D * hisResZ12 = (TH2D*)gDirectory->Get("RMS Z1_1");
1785    TH3F * hisResZ2 = (TH3F*)fRMSZ->At(2);
1786    hisResZ2->FitSlicesZ();
1787    TH2D * hisResZ22 = (TH2D*)gDirectory->Get("RMS Z2_1");
1788    //
1789    ps->NewPage();
1790    c.cd(1);
1791    hisResZ02->Fit(frms, "qn0");         // valgrind
1792    hisResZ02->Draw("surf1");
1793    c.cd(2);
1794    MakeDiff(hisResZ02,frms)->Draw("surf1");
1795    c.Update();
1796 //    c.SaveAs("RMSZPad0.eps");
1797    ps->NewPage();
1798    c.cd(1);
1799    hisResZ12->Fit(frms, "qn0");
1800    hisResZ12->Draw("surf1");
1801    c.cd(2);
1802    MakeDiff(hisResZ12,frms)->Draw("surf1");
1803    c.Update();
1804 //    c.SaveAs("RMSZPad1.eps");
1805    ps->NewPage();
1806    c.cd(1);
1807    hisResZ22->Fit(frms, "qn0");         // valgrind  1 block possibly lost
1808    hisResZ22->Draw("surf1");
1809    c.cd(2);
1810    MakeDiff(hisResZ22,frms)->Draw("surf1");
1811    c.Update();
1812 //    c.SaveAs("RMSZPad2.eps");
1813
1814    // write calculated resoltuion rms to 'rms.txt'
1815    ofstream filerms("rms.txt");
1816    filerms<<"Pad 0.75 cm"<<"\n";
1817    hisResY02->Fit(frms, "qn0");         // valgrind   23 blocks indirectly lost
1818    filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1819    hisResZ02->Fit(frms, "qn0");         // valgrind   23 blocks indirectly lost
1820    filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1821    //
1822    filerms<<"Pad 1.00 cm"<<1<<"\n";
1823    hisResY12->Fit(frms, "qn0");         // valgrind      3,256 bytes in 22 blocks are indirectly lost
1824    filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1825    hisResZ12->Fit(frms, "qn0");         // valgrind    66,036 bytes in 3 blocks are still reachable
1826    filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1827    //
1828    filerms<<"Pad 1.50 cm"<<0<<"\n";
1829    hisResY22->Fit(frms, "qn0");      // valgrind   40,139 bytes in 11 blocks are still reachable   330,180 bytes in 15 blocks are possibly lost
1830    filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1831    hisResZ22->Fit(frms, "qn0");
1832    filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1833
1834    TH1::AddDirectory(kFALSE);
1835    gSystem->ChangeDirectory("..");
1836    ps->Close();
1837    delete ps;
1838    delete frms;
1839 }
1840
1841
1842 void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){
1843    //
1844    // Make tree with resolution parameters
1845    // the result is written to 'resol.root' in directory 'pathname'
1846    // file information are available in fileInfo
1847    // available variables in the tree 'Resol':
1848    //  Entries: number of entries for this resolution point
1849    //  nbins:   number of bins that were accumulated
1850    //  Dim:     direction, Dim==0: y-direction, Dim==1: z-direction
1851    //  Pad:     padSize; short, medium and long
1852    //  Length:  pad length, 0.75, 1, 1.5
1853    //  QMean:   mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
1854    //  Zc:      center of middle bin in drift direction
1855    //  Zm:      mean dirftlength for accumulated Delta-Histograms
1856    //  Zs:      width of driftlength bin
1857    //  AngleC:  center of middle bin in Angle-Direction
1858    //  AngleM:  mean angle for accumulated Delta-Histograms
1859    //  AngleS:  width of Angle-bin
1860    //  Resol:   sigma for gaus fit through Delta-Histograms
1861    //  Sigma:   error of sigma for gaus fit through Delta Histograms
1862    //  MeanR:   mean of the Delta-Histogram
1863    //  SigmaR:  rms of the Delta-Histogram
1864    //  RMSm:    mean of the gaus fit through RMS-Histogram
1865    //  RMS:     sigma of the gaus fit through RMS-Histogram
1866    //  RMSe0:   error of mean of gaus fit in RMS-Histogram
1867    //  RMSe1:   error of sigma of gaus fit in RMS-Histogram
1868    //
1869
1870   if (GetDebugLevel() > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl;
1871   if (GetDebugLevel() > -1) cout << "    relax, the calculation will take a while..." << endl;
1872
1873    gSystem->MakeDirectory(pathName);
1874    gSystem->ChangeDirectory(pathName);
1875    TString kFileName = "resol.root";
1876    TTreeSRedirector fTreeResol(kFileName.Data());
1877
1878    TH3F *resArray[2][3][11];
1879    TH3F *rmsArray[2][3][11];
1880
1881    // load histograms from fArrayQDY and fArrayQDZ
1882    // into resArray and rmsArray
1883    // that is all we need here
1884    for (Int_t idim = 0; idim < 2; idim++){
1885       for (Int_t ipad = 0; ipad < 3; ipad++){
1886          for (Int_t iq = 0; iq <= 10; iq++){
1887             rmsArray[idim][ipad][iq]=0;
1888             resArray[idim][ipad][iq]=0;
1889             Int_t bin = GetBin(iq,ipad);
1890             TH3F *hresl = 0;
1891             if (idim == 0) hresl = (TH3F*)fArrayQDY->At(bin);
1892             if (idim == 1) hresl = (TH3F*)fArrayQDZ->At(bin);
1893             if (!hresl) continue;
1894             resArray[idim][ipad][iq] = (TH3F*) hresl->Clone();
1895             resArray[idim][ipad][iq]->SetDirectory(0);
1896             TH3F * hreslRMS = 0;
1897             if (idim == 0) hreslRMS = (TH3F*)fArrayQRMSY->At(bin);
1898             if (idim == 1) hreslRMS = (TH3F*)fArrayQRMSZ->At(bin);
1899             if (!hreslRMS) continue;
1900             rmsArray[idim][ipad][iq] = (TH3F*) hreslRMS->Clone();
1901             rmsArray[idim][ipad][iq]->SetDirectory(0);
1902          }
1903       }
1904    }
1905    if (GetDebugLevel() > -1) cout << "Histograms loaded, starting to proces..." << endl;
1906
1907    //--------------------------------------------------------------------------------------------
1908
1909    char name[200];
1910    Double_t qMean;
1911    Double_t zMean, angleMean, zCenter, angleCenter;
1912    Double_t zSigma, angleSigma;
1913    TH1D *projectionRes = new TH1D("projectionRes", "projectionRes", 50, -1, 1);
1914    TH1D *projectionRms = new TH1D("projectionRms", "projectionRms", 50, -1, 1);
1915    TF1 *fitFunction = new TF1("fitFunction", "gaus");
1916    Float_t entriesQ = 0;
1917    Int_t loopCounter = 1;
1918
1919    for (Int_t idim = 0; idim < 2; idim++){
1920       // Loop y-z corrdinate
1921       for (Int_t ipad = 0; ipad < 3; ipad++){
1922          // loop pad type
1923          for (Int_t iq = -1; iq < 10; iq++){
1924             // LOOP Q
1925            if (GetDebugLevel() > -1)
1926                cout << "Loop-counter, this is loop " << loopCounter << " of 66, ("
1927                   << (Int_t)((loopCounter)/66.*100) << "% done), "
1928                   << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << "  \r" << std::flush;
1929             loopCounter++;
1930             TH3F *hres = 0;
1931             TH3F *hrms = 0;
1932             qMean    = 0;
1933             entriesQ = 0;
1934
1935             // calculate qMean
1936             if (iq == -1){
1937                // integrated spectra
1938                for (Int_t iql = 0; iql < 10; iql++){
1939                   Int_t bin = GetBin(iql,ipad);
1940                   TH3F *hresl = resArray[idim][ipad][iql];
1941                   TH3F *hrmsl = rmsArray[idim][ipad][iql];
1942                   if (!hresl) continue;
1943                   if (!hrmsl) continue;
1944                   entriesQ += hresl->GetEntries();
1945                   qMean += hresl->GetEntries() * GetQ(bin);
1946                   if (!hres) {
1947                      hres = (TH3F*)hresl->Clone();
1948                      hrms = (TH3F*)hrmsl->Clone();
1949                   }
1950                   else{
1951                      hres->Add(hresl);
1952                      hrms->Add(hrmsl);
1953                   }
1954                }
1955                qMean /= entriesQ;
1956                qMean *= -1.;  // integral mean charge
1957             }
1958             else {
1959                // loop over neighboured Q-bins
1960                // accumulate entries from neighboured Q-bins
1961                for (Int_t iql = iq - 1; iql <= iq + 1; iql++){
1962                   if (iql < 0) continue;
1963                   Int_t bin = GetBin(iql,ipad);
1964                   TH3F * hresl = resArray[idim][ipad][iql];
1965                   TH3F * hrmsl = rmsArray[idim][ipad][iql];
1966                   if (!hresl) continue;
1967                   if (!hrmsl) continue;
1968                   entriesQ += hresl->GetEntries();
1969                   qMean += hresl->GetEntries() * GetQ(bin);
1970                   if (!hres) {
1971                      hres = (TH3F*) hresl->Clone();
1972                      hrms = (TH3F*) hrmsl->Clone();
1973                   }
1974                   else{
1975                      hres->Add(hresl);
1976                      hrms->Add(hrmsl);
1977                   }
1978                }
1979                qMean/=entriesQ;
1980             }
1981
1982             TAxis *xAxisDriftLength = hres->GetXaxis();   // driftlength / z - axis
1983             TAxis *yAxisAngle       = hres->GetYaxis();   // angle axis
1984             TAxis *zAxisDelta       = hres->GetZaxis();   // delta axis
1985             TAxis *zAxisRms         = hrms->GetZaxis();   // rms axis
1986
1987             // loop over all angle bins
1988             for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) {
1989                angleCenter = yAxisAngle->GetBinCenter(ibinyAngle);
1990                // loop over all driftlength bins
1991                for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) {
1992                   zCenter    = xAxisDriftLength->GetBinCenter(ibinxDL);
1993                   zSigma     = xAxisDriftLength->GetBinWidth(ibinxDL);
1994                   angleSigma = yAxisAngle->GetBinWidth(ibinyAngle);
1995                   zMean      = zCenter;      // changens, when more statistic is accumulated
1996                   angleMean  = angleCenter;  // changens, when more statistic is accumulated
1997
1998                   // create 2 1D-Histograms, projectionRes and projectionRms
1999                   // these histograms are delta histograms for given direction, padSize, chargeBin,
2000                   // angleBin and driftLengthBin
2001                   // later on they will be fitted with a gausian, its sigma is the resoltuion...
2002                   sprintf(name,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter);
2003                   // TH1D * projectionRes = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
2004                   projectionRes->SetNameTitle(name, name);
2005                   sprintf(name,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter);
2006                   // TH1D * projectionRms =  new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
2007                   projectionRms->SetNameTitle(name, name);
2008
2009                   projectionRes->Reset();
2010                   projectionRes->SetBins(zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
2011                   projectionRms->Reset();
2012                   projectionRms->SetBins(zAxisRms->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
2013                   projectionRes->SetDirectory(0);
2014                   projectionRms->SetDirectory(0);
2015
2016                   Double_t entries = 0;
2017                   Int_t    nbins   = 0;   // counts, how many bins were accumulated
2018                   zMean     = 0;
2019                   angleMean = 0;
2020                   entries   = 0;
2021
2022                   // fill projectionRes and projectionRms for given dim, ipad and iq,
2023                   // as well as for given angleBin and driftlengthBin
2024                   // if this gives not enough statistic, include neighbourhood
2025                   // (angle and driftlength) successifely
2026                   for (Int_t dbin = 0; dbin <= 8; dbin++){              // delta-bins around centered angleBin and driftlengthBin
2027                      for (Int_t dbiny2 = -1; dbiny2 <= 1; dbiny2++) {   // delta-bins in angle direction
2028                         for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ // delta-bins in driftlength direction
2029                            if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue;   // add each bin only one time !
2030                            Int_t binx2 = ibinxDL + dbinx2;                       // position variable in x (driftlength) direction
2031                            Int_t biny2 = ibinyAngle + dbiny2;                    // position variable in y (angle)  direction
2032                            if (binx2 < 1 || biny2 < 1) continue;                 // don't go out of the histogram!
2033                            if (binx2 >= xAxisDriftLength->GetNbins()) continue;  // don't go out of the histogram!
2034                            if (biny2 >= yAxisAngle->GetNbins()) continue;        // don't go out of the histogram!
2035                            nbins++;                                              // count the number of accumulated bins
2036                            // Fill resolution histo
2037                            for (Int_t ibin3 = 1; ibin3 < zAxisDelta->GetNbins(); ibin3++) {
2038                               // Int_t content = (Int_t)hres->GetBinContent(binx2, biny2, ibin3);     // unused variable
2039                               projectionRes->Fill(zAxisDelta->GetBinCenter(ibin3), hres->GetBinContent(binx2, biny2, ibin3));
2040                               entries   += hres->GetBinContent(binx2, biny2, ibin3);
2041                               zMean     += hres->GetBinContent(binx2, biny2, ibin3) * xAxisDriftLength->GetBinCenter(binx2);
2042                               angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yAxisAngle->GetBinCenter(biny2);
2043                            }  // ibin3 loop
2044                            // fill RMS histo
2045                            for (Int_t ibin3 = 1; ibin3 < zAxisRms->GetNbins(); ibin3++) {
2046                               projectionRms->Fill(zAxisRms->GetBinCenter(ibin3), hrms->GetBinContent(binx2, biny2, ibin3));
2047                            }
2048                         }  //dbinx2 loop
2049                         if (entries > minEntries) break; // enough statistic accumulated
2050                      }  // dbiny2 loop
2051                      if (entries > minEntries) break;    // enough statistic accumulated
2052                   }  // dbin loop
2053                   if ( entries< minEntries) continue;  // when it was absolutly impossible to get enough statistic, don't write this point into the resolution tree
2054                   zMean /= entries;
2055                   angleMean /= entries;
2056
2057                   if (entries > minEntries) {
2058                      //  when enough statistic is accumulated
2059                      //  fit Delta histograms with a gausian
2060                      //  of the gausian is the resolution (resol), its fit error is sigma
2061                      //  store also mean and RMS of the histogram
2062                      Float_t xmin     = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
2063                      Float_t xmax     = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
2064
2065 //                      projectionRes->Fit("gaus", "q0", "", xmin, xmax);
2066 //                      Float_t resol    = projectionRes->GetFunction("gaus")->GetParameter(2);
2067 //                      Float_t sigma    = projectionRes->GetFunction("gaus")->GetParError(2);
2068                      fitFunction->SetMaximum(xmax);
2069                      fitFunction->SetMinimum(xmin);
2070                      projectionRes->Fit("fitFunction", "qN0", "", xmin, xmax);
2071                      Float_t resol    = fitFunction->GetParameter(2);
2072                      Float_t sigma    = fitFunction->GetParError(2);
2073
2074                      Float_t meanR    = projectionRes->GetMean();
2075                      Float_t sigmaR   = projectionRes->GetRMS();
2076                      // fit also RMS histograms with a gausian
2077                      // store mean and sigma of the gausian in rmsMean and rmsSigma
2078                      // store also the fit errors in errorRMS and errorSigma
2079                      xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
2080                      xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
2081
2082 //                      projectionRms->Fit("gaus","q0","",xmin,xmax);
2083 //                      Float_t rmsMean    = projectionRms->GetFunction("gaus")->GetParameter(1);
2084 //                      Float_t rmsSigma   = projectionRms->GetFunction("gaus")->GetParameter(2);
2085 //                      Float_t errorRMS   = projectionRms->GetFunction("gaus")->GetParError(1);
2086 //                      Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2);
2087                      projectionRms->Fit("fitFunction", "qN0", "", xmin, xmax);
2088                      Float_t rmsMean    = fitFunction->GetParameter(1);
2089                      Float_t rmsSigma   = fitFunction->GetParameter(2);
2090                      Float_t errorRMS   = fitFunction->GetParError(1);
2091                      Float_t errorSigma = fitFunction->GetParError(2);
2092
2093                      Float_t length = 0.75;
2094                      if (ipad == 1) length = 1;
2095                      if (ipad == 2) length = 1.5;
2096
2097                      fTreeResol<<"Resol"<<
2098                         "Entries="<<entries<<      // number of entries for this resolution point
2099                         "nbins="<<nbins<<          // number of bins that were accumulated
2100                         "Dim="<<idim<<             // direction, Dim==0: y-direction, Dim==1: z-direction
2101                         "Pad="<<ipad<<             // padSize; short, medium and long
2102                         "Length="<<length<<        // pad length, 0.75, 1, 1.5
2103                         "QMean="<<qMean<<          // mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
2104                         "Zc="<<zCenter<<           // center of middle bin in drift direction
2105                         "Zm="<<zMean<<             // mean dirftlength for accumulated Delta-Histograms
2106                         "Zs="<<zSigma<<            // width of driftlength bin
2107                         "AngleC="<<angleCenter<<   // center of middle bin in Angle-Direction
2108                         "AngleM="<<angleMean<<     // mean angle for accumulated Delta-Histograms
2109                         "AngleS="<<angleSigma<<    // width of Angle-bin
2110                         "Resol="<<resol<<          // sigma for gaus fit through Delta-Histograms
2111                         "Sigma="<<sigma<<          // error of sigma for gaus fit through Delta Histograms
2112                         "MeanR="<<meanR<<          // mean of the Delta-Histogram
2113                         "SigmaR="<<sigmaR<<        // rms of the Delta-Histogram
2114                         "RMSm="<<rmsMean<<         // mean of the gaus fit through RMS-Histogram
2115                         "RMSs="<<rmsSigma<<        // sigma of the gaus fit through RMS-Histogram
2116                         "RMSe0="<<errorRMS<<       // error of mean of gaus fit in RMS-Histogram
2117                         "RMSe1="<<errorSigma<<     // error of sigma of gaus fit in RMS-Histogram
2118                         "\n";
2119                      if (GetDebugLevel() > 5) {
2120                         projectionRes->SetDirectory(fTreeResol.GetFile());
2121                         projectionRes->Write(projectionRes->GetName());
2122                         projectionRes->SetDirectory(0);
2123                         projectionRms->SetDirectory(fTreeResol.GetFile());
2124                         projectionRms->Write(projectionRms->GetName());
2125                         projectionRes->SetDirectory(0);
2126                      }
2127                   }  // if (projectionRes->GetSum() > minEntries)
2128                }  // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++)
2129             }  // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++)
2130
2131          }  // iq-loop
2132       }  // ipad-loop
2133    }  // idim-loop
2134    delete projectionRes;
2135    delete projectionRms;
2136
2137 //    TFile resolFile(fTreeResol.GetFile());
2138    TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries));
2139    fileInfo.Write("fileInfo");
2140 //    resolFile.Close();
2141 //    fTreeResol.GetFile()->Close();
2142    if (GetDebugLevel() > -1) cout << endl;
2143    if (GetDebugLevel() > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl;
2144    gSystem->ChangeDirectory("..");
2145 }
2146
2147
2148
2149
2150
2151 Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
2152    //
2153    // function to merge several AliTPCcalibTracks objects after PROOF calculation
2154    // The object's histograms are merged via their merge functions
2155    // Be carefull: histograms are linked to a file, switch this off by TH1::AddDirectory(kFALSE) !!!
2156    //
2157
2158   if (GetDebugLevel() > 0) cout << " *****  this is AliTPCcalibTracks::Merge(TCollection *collectionList)  *****"<< endl;
2159    if (!collectionList) return 0;
2160    if (collectionList->IsEmpty()) return -1;
2161
2162    if (GetDebugLevel() > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl;     //    REMOVE THIS LINE!!!!!!!!!!!!!!!!!1
2163    if (GetDebugLevel() > 5) cout << " the list in the merge-function looks as follows: " << endl;
2164    collectionList->Print();
2165
2166    // create a list for each data member
2167    TList* deltaYList = new TList;
2168    TList* deltaZList = new TList;
2169    TList* arrayAmpRowList = new TList;
2170    TList* rejectedTracksList = new TList;
2171    TList* hclusList = new TList;
2172    TList* clusterCutHistoList = new TList;
2173    TList* arrayAmpList = new TList;
2174    TList* arrayQDYList = new TList;
2175    TList* arrayQDZList = new TList;
2176    TList* arrayQRMSYList = new TList;
2177    TList* arrayQRMSZList = new TList;
2178    TList* arrayChargeVsDriftlengthList = new TList;
2179    TList* calPadRegionChargeVsDriftlengthList = new TList;
2180    TList* hclusterPerPadrowList = new TList;
2181    TList* hclusterPerPadrowRawList = new TList;
2182    TList* resolYList = new TList;
2183    TList* resolZList = new TList;
2184    TList* rMSYList = new TList;
2185    TList* rMSZList = new TList;
2186
2187 //    TList* nRowsList = new TList;
2188 //    TList* nSectList = new TList;
2189 //    TList* fileNoList = new TList;
2190
2191    TIterator *listIterator = collectionList->MakeIterator();
2192    AliTPCcalibTracks *calibTracks = 0;
2193    if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;
2194    Int_t counter = 0;
2195    while ( (calibTracks = dynamic_cast<AliTPCcalibTracks*> (listIterator->Next())) ){
2196       // loop over all entries in the collectionList and get dataMembers into lists
2197       if (!calibTracks) continue;
2198
2199       deltaYList->Add( calibTracks->GetfDeltaY() );
2200       deltaZList->Add( calibTracks->GetfDeltaZ() );
2201       arrayAmpRowList->Add(calibTracks->GetfArrayAmpRow());
2202       arrayAmpList->Add(calibTracks->GetfArrayAmp());
2203       arrayQDYList->Add(calibTracks->GetfArrayQDY());
2204       arrayQDZList->Add(calibTracks->GetfArrayQDZ());
2205       arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
2206       arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ());
2207       resolYList->Add(calibTracks->GetfResolY());
2208       resolZList->Add(calibTracks->GetfResolZ());
2209       rMSYList->Add(calibTracks->GetfRMSY());
2210       rMSZList->Add(calibTracks->GetfRMSZ());
2211       arrayChargeVsDriftlengthList->Add(calibTracks->GetfArrayChargeVsDriftlength());
2212       calPadRegionChargeVsDriftlengthList->Add(calibTracks->GetCalPadRegionchargeVsDriftlength());
2213       hclusList->Add(calibTracks->GetfHclus());
2214       rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
2215       clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
2216       hclusterPerPadrowList->Add(calibTracks->GetfHclusterPerPadrow());
2217       hclusterPerPadrowRawList->Add(calibTracks->GetfHclusterPerPadrowRaw());
2218       //
2219       if (fCalPadClusterPerPad && calibTracks->GetfCalPadClusterPerPad())
2220         fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
2221       //      fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
2222       counter++;
2223       if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
2224    }
2225
2226
2227    // merge data members
2228    if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl;
2229    fDeltaY->Merge(deltaYList);
2230    fDeltaZ->Merge(deltaZList);
2231    fHclus->Merge(hclusList);
2232    fClusterCutHisto->Merge(clusterCutHistoList);
2233    fRejectedTracksHisto->Merge(rejectedTracksList);
2234    fHclusterPerPadrow->Merge(hclusterPerPadrowList);
2235    fHclusterPerPadrowRaw->Merge(hclusterPerPadrowRawList);
2236
2237    TObjArray* objarray = 0;
2238    TH1* hist = 0;
2239    TList* histList = 0;
2240    TIterator *objListIterator = 0;
2241
2242    if (GetDebugLevel() > 0) cout << "merging fArrayAmpRows..." << endl;
2243    // merge fArrayAmpRows
2244    for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++ ) {  // loop over data member, i<72
2245       objListIterator = arrayAmpRowList->MakeIterator();
2246       histList = new TList;
2247       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2248          // loop over arrayAmpRowList, get TObjArray, get object at position i, cast it into TProfile
2249          if (!objarray) continue;
2250          hist = (TProfile*)(objarray->At(i));
2251          histList->Add(hist);
2252       }
2253       ((TProfile*)(fArrayAmpRow->At(i)))->Merge(histList);
2254       delete histList;
2255       delete objListIterator;
2256    }
2257
2258    if (GetDebugLevel() > 0) cout << "merging fArrayAmps..." << endl;
2259    // merge fArrayAmps
2260    for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++ ) {  // loop over data member, i<72
2261       TIterator *cobjListIterator = arrayAmpList->MakeIterator();
2262       histList = new TList;
2263       while (( objarray =  (TObjArray*)cobjListIterator->Next() )) {
2264          // loop over arrayAmpList, get TObjArray, get object at position i, cast it into TH1F
2265          if (!objarray) continue;
2266          hist = (TH1F*)(objarray->At(i));
2267          histList->Add(hist);
2268       }
2269       ((TH1F*)(fArrayAmp->At(i)))->Merge(histList);
2270       delete histList;
2271       delete cobjListIterator;
2272    }
2273
2274    if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
2275    // merge fArrayQDY
2276    for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
2277       objListIterator = arrayQDYList->MakeIterator();
2278       histList = new TList;
2279       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2280          // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F
2281          if (!objarray) continue;
2282          hist = (TH3F*)(objarray->At(i));
2283          histList->Add(hist);
2284       }
2285       ((TH3F*)(fArrayQDY->At(i)))->Merge(histList);
2286       delete histList;
2287       delete objListIterator;
2288    }
2289
2290    if (GetDebugLevel() > 0) cout << "merging fArrayQDZ..." << endl;
2291    // merge fArrayQDZ
2292    for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300
2293       objListIterator = arrayQDZList->MakeIterator();
2294       histList = new TList;
2295       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2296          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2297          if (!objarray) continue;
2298          hist = (TH3F*)(objarray->At(i));
2299          histList->Add(hist);
2300       }
2301       ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList);
2302       delete histList;
2303       delete objListIterator;
2304    }
2305
2306    if (GetDebugLevel() > 0) cout << "merging fArrayQRMSY..." << endl;
2307    // merge fArrayQRMSY
2308    for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300
2309       objListIterator = arrayQRMSYList->MakeIterator();
2310       histList = new TList;
2311       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2312          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2313          if (!objarray) continue;
2314          hist = (TH3F*)(objarray->At(i));
2315          histList->Add(hist);
2316       }
2317       ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList);
2318       delete histList;
2319       delete objListIterator;
2320    }
2321
2322    if (GetDebugLevel() > 0) cout << "merging fArrayQRMSZ..." << endl;
2323    // merge fArrayQRMSZ
2324    for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300
2325       objListIterator = arrayQRMSZList->MakeIterator();
2326       histList = new TList;
2327       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2328          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2329          if (!objarray) continue;
2330          hist = (TH3F*)(objarray->At(i));
2331          histList->Add(hist);
2332       }
2333       ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList);
2334       delete histList;
2335       delete objListIterator;
2336    }
2337
2338    if (GetDebugLevel() > 0) cout << "merging fArrayChargeVsDriftlength..." << endl;
2339    // merge fArrayChargeVsDriftlength
2340    for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) { // loop over data member, i < 300
2341       objListIterator = arrayChargeVsDriftlengthList->MakeIterator();
2342       histList = new TList;
2343       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2344          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TProfile
2345          if (!objarray) continue;
2346          hist = (TProfile*)(objarray->At(i));
2347          histList->Add(hist);
2348       }
2349       ((TProfile*)(fArrayChargeVsDriftlength->At(i)))->Merge(histList);
2350       delete histList;
2351       delete objListIterator;
2352    }
2353
2354    if (GetDebugLevel() > 0) cout << "merging fcalPadRegionChargeVsDriftlength..." << endl;
2355    // merge fcalPadRegionChargeVsDriftlength
2356    AliTPCCalPadRegion *cpr = 0x0;
2357
2358    /*
2359    TIterator *regionIterator = fcalPadRegionChargeVsDriftlength->MakeIterator();
2360    while (hist = (TProfile*)regionIterator->Next()) {
2361       // loop over all calPadRegion's in destination calibTracks object
2362          objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
2363          while (( cpr =  (AliTPCCalPadRegion*)objListIterator->Next() )) {
2364
2365
2366       hist->Merge(...);
2367    }
2368    */
2369
2370    for (UInt_t isec = 0; isec < 36; isec++) {
2371       for (UInt_t padSize = 0; padSize < 3; padSize++){
2372          objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
2373          histList = new TList;
2374          while (( cpr =  (AliTPCCalPadRegion*)objListIterator->Next() )) {
2375             // loop over calPadRegionChargeVsDriftlengthList, get AliTPCCalPadRegion, get object
2376             if (!cpr) continue;
2377             hist = (TProfile*)cpr->GetObject(isec, padSize);
2378             histList->Add(hist);
2379          }
2380          ((TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isec, padSize)))->Merge(histList);
2381          delete histList;
2382          delete objListIterator;
2383       }
2384    }
2385
2386
2387
2388
2389    if (GetDebugLevel() > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl;
2390    // merge fResolY
2391    for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3
2392       objListIterator = resolYList->MakeIterator();
2393       histList = new TList;
2394       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2395          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2396          if (!objarray) continue;
2397          hist = (TH3F*)(objarray->At(i));
2398          histList->Add(hist);
2399       }
2400       ((TH3F*)(fResolY->At(i)))->Merge(histList);
2401       delete histList;
2402       delete objListIterator;
2403    }
2404
2405    // merge fResolZ
2406    for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3
2407       objListIterator = resolZList->MakeIterator();
2408       histList = new TList;
2409       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2410          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2411          if (!objarray) continue;
2412          hist = (TH3F*)(objarray->At(i));
2413          histList->Add(hist);
2414       }
2415       ((TH3F*)(fResolZ->At(i)))->Merge(histList);
2416       delete histList;
2417       delete objListIterator;
2418    }
2419
2420    // merge fRMSY
2421    for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3
2422       objListIterator = rMSYList->MakeIterator();
2423       histList = new TList;
2424       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2425          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2426          if (!objarray) continue;
2427          hist = (TH3F*)(objarray->At(i));
2428          histList->Add(hist);
2429       }
2430       ((TH3F*)(fRMSY->At(i)))->Merge(histList);
2431       delete histList;
2432       delete objListIterator;
2433    }
2434
2435    // merge fRMSZ
2436    for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3
2437       objListIterator = rMSZList->MakeIterator();
2438       histList = new TList;
2439       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
2440          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2441          if (!objarray) continue;
2442          hist = (TH3F*)(objarray->At(i));
2443          histList->Add(hist);
2444       }
2445       ((TH3F*)(fRMSZ->At(i)))->Merge(histList);
2446       delete histList;
2447       delete objListIterator;
2448    }
2449
2450    delete deltaYList;
2451    delete deltaZList;
2452    delete arrayAmpRowList;
2453    delete arrayAmpList;
2454    delete arrayQDYList;
2455    delete arrayQDZList;
2456    delete arrayQRMSYList;
2457    delete arrayQRMSZList;
2458    delete resolYList;
2459    delete resolZList;
2460    delete rMSYList;
2461    delete rMSZList;
2462    delete listIterator;
2463
2464    if (GetDebugLevel() > 0) cout << "merging done!" << endl;
2465
2466    return 1;
2467 }
2468
2469
2470 AliTPCcalibTracks* AliTPCcalibTracks::TestMerge(AliTPCcalibTracks *ct, AliTPCClusterParam *clusterParam, Int_t nCalTracks){
2471    //
2472    // function to test AliTPCcalibTrack::Merge:
2473    // in the file 'f' is a AliTPCcalibTrack object with name "calibTracks"
2474    // this object is appended 'nCalTracks' times to a TList
2475    // A new AliTPCcalibTrack object is created which merges the list
2476    // this object is returned
2477    //
2478    /*
2479    // .L AliTPCcalibTracks.cxx+g
2480    .L libTPCcalib.so
2481    TFile f("Output.root");
2482    AliTPCcalibTracks* calTracks = (AliTPCcalibTracks*)f.Get("calibTracks");
2483    //f.Close();
2484    TFile clusterParamFile("/u/lbozyk/calibration/workdir/calibTracks/TPCClusterParam.root");
2485    AliTPCClusterParam *clusterParam  =  (AliTPCClusterParam *) clusterParamFile.Get("Param");
2486    clusterParamFile.Close();
2487
2488    AliTPCcalibTracks::TestMerge(calTracks, clusterParam);
2489    */
2490    TList *list = new TList();
2491    if (ct == 0 || clusterParam == 0) return 0;
2492    cout << "making list with " << nCalTracks << " AliTPCcalibTrack objects" << endl;
2493    for (Int_t i = 0; i < nCalTracks; i++) {
2494       if (i%10==0) cout << "Adding element " << i << " of " << nCalTracks << endl;
2495       list->Add(new AliTPCcalibTracks(*ct));
2496    }
2497
2498    // only for check at the end
2499    AliTPCcalibTracks* cal1 = new AliTPCcalibTracks(*ct);
2500    Double_t cal1Entries = ((TH1F*)cal1->GetfArrayAmpRow()->At(5))->GetEntries();
2501 //    Double_t cal1Entries = 5; //((TH1F*)ct->GetfArrayAmpRow()->At(5))->GetEntries();
2502
2503    cout  << "The list contains " << list->GetEntries() << " entries. " << endl;
2504
2505
2506    AliTPCcalibTracksCuts *cuts = new AliTPCcalibTracksCuts(20, 0.4, 0.5, 0.13, 0.018);
2507    AliTPCcalibTracks* cal = new AliTPCcalibTracks("calTracksMerged", "calTracksMerged", clusterParam, cuts, 5);
2508    cal->Merge(list);
2509
2510    cout << "cal->GetfArrayAmpRow()->At(5)->Print():" << endl;
2511    cal->GetfArrayAmpRow()->At(5)->Print();
2512    Double_t calEntries = ((TH1F*)cal->GetfArrayAmpRow()->At(5))->GetEntries();
2513
2514    cout << "cal1->GetfArrayAmpRow()->At(5))->GetEntries() = " << cal1Entries << endl;
2515    cout << " cal->GetfArrayAmpRow()->At(5))->GetEntries() = " <<  calEntries << endl;
2516    printf("That means there were %f / %f = %f AliTPCcalibTracks-Objects merged. \n",
2517       calEntries, cal1Entries, ((Double_t)calEntries/cal1Entries));
2518
2519    return cal;
2520
2521 }
2522
2523
2524 void  AliTPCcalibTracks::MakeQPosNormAll(TTree * chainres, AliTPCClusterParam * param, Int_t maxPoints){
2525   //
2526   // Make position corrections
2527   // for the moment Only using debug streamer
2528   // chainres  - debug tree
2529   // param     - parameters to be updated
2530   // maxPoints - maximal number of points using for fit
2531   // verbose   - print info flag
2532   //
2533   // Current implementation - only using debug streamers
2534   //
2535
2536   /*
2537     //Defaults
2538     Int_t maxPoints=100000;
2539   */
2540   /*
2541     Usage:
2542     //0. Load libraries
2543     gSystem->Load("libANALYSIS");
2544     gSystem->Load("libSTAT");
2545     gSystem->Load("libTPCcalib");
2546
2547
2548     //1. Load Parameters to be modified:
2549     //e.g:
2550
2551     AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
2552     AliCDBManager::Instance()->SetRun(0)
2553     AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
2554
2555     //2. Load chain from debug streamers
2556     //
2557     //e.g
2558     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
2559     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
2560     AliXRDPROOFtoolkit tool;
2561     TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
2562     chainres->Lookup();
2563     //3. Do fits and store results
2564     //
2565     AliTPCcalibTracks::MakeQPosNormAll(chainres,param,200000,0) ;
2566     TFile f("paramout.root","recreate");
2567     param->Write("clusterParam");
2568     f.Close();
2569     //4. Verification
2570     TFile f2("paramout.root");
2571     AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
2572     param2->SetInstance(param2);
2573     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
2574
2575    */
2576
2577
2578   TStatToolkit toolkit;
2579   Double_t chi2;
2580   TVectorD fitParamY0;
2581   TVectorD fitParamY1;
2582   TVectorD fitParamZ0;
2583   TVectorD fitParamZ1;
2584   TMatrixD covMatrix;
2585   Int_t npoints;
2586
2587   chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
2588   chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
2589   chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
2590   chainres->SetAlias("st","(sin(dt)-dt)");
2591   //
2592   chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
2593   //
2594   //
2595   //
2596   TCut cutA("1");
2597   TString fstringY="";
2598   //
2599   fstringY+="(dp)++";            //1
2600   fstringY+="(dp)*di++";         //2
2601   fstringY+="(sp)++";            //3
2602   fstringY+="(sp)*di++";         //4
2603   TString fstringZ="";
2604   fstringZ+="(dt)++";            //1
2605   fstringZ+="(dt)*di++";         //2
2606   fstringZ+="(st)++";            //3
2607   fstringZ+="(st)*di++";         //4
2608   //
2609   // Z corrections
2610   //
2611   TString *strZ0 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
2612   printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2613   param->fPosZcor[0] = (TVectorD*) fitParamZ0.Clone();
2614   //
2615   TString *strZ1 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
2616   //
2617   printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2618   param->fPosZcor[1] = (TVectorD*) fitParamZ1.Clone();
2619   param->fPosZcor[2] = (TVectorD*) fitParamZ1.Clone();
2620   //
2621   // Y corrections
2622   //
2623   TString *strY0 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
2624   printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2625   param->fPosYcor[0] = (TVectorD*) fitParamY0.Clone();
2626
2627
2628   TString *strY1 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
2629   //
2630   printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2631   param->fPosYcor[1] = (TVectorD*) fitParamY1.Clone();
2632   param->fPosYcor[2] = (TVectorD*) fitParamY1.Clone();
2633   //
2634   //
2635   //
2636   chainres->SetAlias("fitZ0",strZ0->Data());
2637   chainres->SetAlias("fitZ1",strZ1->Data());
2638   chainres->SetAlias("fitY0",strY0->Data());
2639   chainres->SetAlias("fitY1",strY1->Data());
2640   //  chainres->Draw("Cl.fZ-PZ0.fElements[0]","CSigmaY0<0.7&&CSigmaZ0<0.7"+cutA,"",10000);
2641 }
2642
2643
2644