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