ANALYSIS/AliESDtrackCuts.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 /* $Id: AliESDtrackCuts.cxx 24534 2008-03-16 22:22:11Z fca $ */
  17
  18 #include "AliESDtrackCuts.h"
  19
  20 #include <AliESDtrack.h>
  21 #include <AliESDVertex.h>
  22 #include <AliESDEvent.h>
  23 #include <AliLog.h>
  24
  25 #include <TTree.h>
  26 #include <TCanvas.h>
  27 #include <TDirectory.h>
  28
  29 //____________________________________________________________________
  30 ClassImp(AliESDtrackCuts)
  31
  32 // Cut names
  33 const Char_t* AliESDtrackCuts::fgkCutNames[kNCuts] = {
  34  "require TPC refit",
  35  "require ITS refit",
  36  "n clusters TPC",
  37  "n clusters ITS",
  38  "#Chi^{2}/clusters TPC",
  39  "#Chi^{2}/clusters ITS",
  40  "cov 11",
  41  "cov 22",
  42  "cov 33",
  43  "cov 44",
  44  "cov 55",
  45  "trk-to-vtx",
  46  "trk-to-vtx failed",
  47  "kink daughters",
  48  "p",
  49  "p_{T}",
  50  "p_{x}",
  51  "p_{y}",
  52  "p_{z}",
  53  "y",
  54  "eta",
  55  "trk-to-vtx dca absolute",
  56  "trk-to-vtx dca xy absolute",
  57  "trk-to-vtx dca z absolute",
  58  "SPD cluster requirement",
  59  "SDD cluster requirement",
  60  "SSD cluster requirement"
  61 };
  62
  63 //____________________________________________________________________
  64 AliESDtrackCuts::AliESDtrackCuts(const Char_t* name, const Char_t* title) : AliAnalysisCuts(name,title),
  65   fCutMinNClusterTPC(0),
  66   fCutMinNClusterITS(0),
  67   fCutMaxChi2PerClusterTPC(0),
  68   fCutMaxChi2PerClusterITS(0),
  69   fCutMaxC11(0),
  70   fCutMaxC22(0),
  71   fCutMaxC33(0),
  72   fCutMaxC44(0),
  73   fCutMaxC55(0),
  74   fCutAcceptKinkDaughters(0),
  75   fCutRequireTPCRefit(0),
  76   fCutRequireITSRefit(0),
  77   fCutNsigmaToVertex(0),
  78   fCutSigmaToVertexRequired(0),
  79   fCutDCAToVertex(0),
  80   fCutDCAToVertexXY(0),
  81   fCutDCAToVertexZ(0),
  82   fPMin(0),
  83   fPMax(0),
  84   fPtMin(0),
  85   fPtMax(0),
  86   fPxMin(0),
  87   fPxMax(0),
  88   fPyMin(0),
  89   fPyMax(0),
  90   fPzMin(0),
  91   fPzMax(0),
  92   fEtaMin(0),
  93   fEtaMax(0),
  94   fRapMin(0),
  95   fRapMax(0),
  96   fHistogramsOn(0),
  97   ffDTheoretical(0),
  98   fhCutStatistics(0),
  99   fhCutCorrelation(0)
 100 {
 101   //
 102   // constructor
 103   //
 104
 105   Init();
 106
 107   //##############################################################################
 108   // setting default cuts
 109   SetMinNClustersTPC();
 110   SetMinNClustersITS();
 111   SetMaxChi2PerClusterTPC();
 112   SetMaxChi2PerClusterITS();
 113   SetMaxCovDiagonalElements();
 114   SetRequireTPCRefit();
 115   SetRequireITSRefit();
 116   SetAcceptKingDaughters();
 117   SetMaxNsigmaToVertex();
 118   SetMaxDCAToVertex();
 119   SetMaxDCAToVertexXY();
 120   SetMaxDCAToVertexZ();
 121   SetPRange();
 122   SetPtRange();
 123   SetPxRange();
 124   SetPyRange();
 125   SetPzRange();
 126   SetEtaRange();
 127   SetRapRange();
 128   SetClusterRequirementITS(kSPD);
 129   SetClusterRequirementITS(kSDD);
 130   SetClusterRequirementITS(kSSD);
 131
 132   SetHistogramsOn();
 133 }
 134
 135 //_____________________________________________________________________________
 136 AliESDtrackCuts::AliESDtrackCuts(const AliESDtrackCuts &c) : AliAnalysisCuts(c),
 137   fCutMinNClusterTPC(0),
 138   fCutMinNClusterITS(0),
 139   fCutMaxChi2PerClusterTPC(0),
 140   fCutMaxChi2PerClusterITS(0),
 141   fCutMaxC11(0),
 142   fCutMaxC22(0),
 143   fCutMaxC33(0),
 144   fCutMaxC44(0),
 145   fCutMaxC55(0),
 146   fCutAcceptKinkDaughters(0),
 147   fCutRequireTPCRefit(0),
 148   fCutRequireITSRefit(0),
 149   fCutNsigmaToVertex(0),
 150   fCutSigmaToVertexRequired(0),
 151   fCutDCAToVertex(0),
 152   fCutDCAToVertexXY(0),
 153   fCutDCAToVertexZ(0),
 154   fPMin(0),
 155   fPMax(0),
 156   fPtMin(0),
 157   fPtMax(0),
 158   fPxMin(0),
 159   fPxMax(0),
 160   fPyMin(0),
 161   fPyMax(0),
 162   fPzMin(0),
 163   fPzMax(0),
 164   fEtaMin(0),
 165   fEtaMax(0),
 166   fRapMin(0),
 167   fRapMax(0),
 168   fHistogramsOn(0),
 169   ffDTheoretical(0),
 170   fhCutStatistics(0),
 171   fhCutCorrelation(0)
 172 {
 173   //
 174   // copy constructor
 175   //
 176
 177   ((AliESDtrackCuts &) c).Copy(*this);
 178 }
 179
 180 AliESDtrackCuts::~AliESDtrackCuts()
 181 {
 182   //
 183   // destructor
 184   //
 185
 186   for (Int_t i=0; i<2; i++) {
 187
 188     if (fhNClustersITS[i])
 189       delete fhNClustersITS[i];
 190     if (fhNClustersTPC[i])
 191       delete fhNClustersTPC[i];
 192     if (fhChi2PerClusterITS[i])
 193       delete fhChi2PerClusterITS[i];
 194     if (fhChi2PerClusterTPC[i])
 195       delete fhChi2PerClusterTPC[i];
 196     if (fhC11[i])
 197       delete fhC11[i];
 198     if (fhC22[i])
 199       delete fhC22[i];
 200     if (fhC33[i])
 201       delete fhC33[i];
 202     if (fhC44[i])
 203       delete fhC44[i];
 204     if (fhC55[i])
 205     delete fhC55[i];
 206
 207     if (fhDXY[i])
 208       delete fhDXY[i];
 209     if (fhDZ[i])
 210       delete fhDZ[i];
 211     if (fhDXYDZ[i])
 212       delete fhDXYDZ[i];
 213     if (fhDXYvsDZ[i])
 214       delete fhDXYvsDZ[i];
 215
 216     if (fhDXYNormalized[i])
 217       delete fhDXYNormalized[i];
 218     if (fhDZNormalized[i])
 219       delete fhDZNormalized[i];
 220     if (fhDXYvsDZNormalized[i])
 221       delete fhDXYvsDZNormalized[i];
 222     if (fhNSigmaToVertex[i])
 223       delete fhNSigmaToVertex[i];
 224     if (fhPt[i])
 225       delete fhPt[i];
 226     if (fhEta[i])
 227       delete fhEta[i];
 228   }
 229
 230   if (ffDTheoretical)
 231     delete ffDTheoretical;
 232
 233   if (fhCutStatistics)
 234     delete fhCutStatistics;
 235   if (fhCutCorrelation)
 236     delete fhCutCorrelation;
 237 }
 238
 239 void AliESDtrackCuts::Init()
 240 {
 241   //
 242   // sets everything to zero
 243   //
 244
 245   fCutMinNClusterTPC = 0;
 246   fCutMinNClusterITS = 0;
 247
 248   fCutMaxChi2PerClusterTPC = 0;
 249   fCutMaxChi2PerClusterITS = 0;
 250
 251   for (Int_t i = 0; i < 3; i++)
 252         fCutClusterRequirementITS[i] = kOff;
 253
 254   fCutMaxC11 = 0;
 255   fCutMaxC22 = 0;
 256   fCutMaxC33 = 0;
 257   fCutMaxC44 = 0;
 258   fCutMaxC55 = 0;
 259
 260   fCutAcceptKinkDaughters = 0;
 261   fCutRequireTPCRefit = 0;
 262   fCutRequireITSRefit = 0;
 263
 264   fCutNsigmaToVertex = 0;
 265   fCutSigmaToVertexRequired = 0;
 266   fCutDCAToVertex = 0;
 267   fCutDCAToVertexXY = 0;
 268   fCutDCAToVertexZ = 0;
 269
 270   fPMin = 0;
 271   fPMax = 0;
 272   fPtMin = 0;
 273   fPtMax = 0;
 274   fPxMin = 0;
 275   fPxMax = 0;
 276   fPyMin = 0;
 277   fPyMax = 0;
 278   fPzMin = 0;
 279   fPzMax = 0;
 280   fEtaMin = 0;
 281   fEtaMax = 0;
 282   fRapMin = 0;
 283   fRapMax = 0;
 284
 285   fHistogramsOn = kFALSE;
 286
 287   for (Int_t i=0; i<2; ++i)
 288   {
 289     fhNClustersITS[i] = 0;
 290     fhNClustersTPC[i] = 0;
 291
 292     fhChi2PerClusterITS[i] = 0;
 293     fhChi2PerClusterTPC[i] = 0;
 294
 295     fhC11[i] = 0;
 296     fhC22[i] = 0;
 297     fhC33[i] = 0;
 298     fhC44[i] = 0;
 299     fhC55[i] = 0;
 300
 301     fhDXY[i] = 0;
 302     fhDZ[i] = 0;
 303     fhDXYDZ[i] = 0;
 304     fhDXYvsDZ[i] = 0;
 305
 306     fhDXYNormalized[i] = 0;
 307     fhDZNormalized[i] = 0;
 308     fhDXYvsDZNormalized[i] = 0;
 309     fhNSigmaToVertex[i] = 0;
 310
 311     fhPt[i] = 0;
 312     fhEta[i] = 0;
 313   }
 314   ffDTheoretical = 0;
 315
 316   fhCutStatistics = 0;
 317   fhCutCorrelation = 0;
 318 }
 319
 320 //_____________________________________________________________________________
 321 AliESDtrackCuts &AliESDtrackCuts::operator=(const AliESDtrackCuts &c)
 322 {
 323   //
 324   // Assignment operator
 325   //
 326
 327   if (this != &c) ((AliESDtrackCuts &) c).Copy(*this);
 328   return *this;
 329 }
 330
 331 //_____________________________________________________________________________
 332 void AliESDtrackCuts::Copy(TObject &c) const
 333 {
 334   //
 335   // Copy function
 336   //
 337
 338   AliESDtrackCuts& target = (AliESDtrackCuts &) c;
 339
 340   target.Init();
 341
 342   target.fCutMinNClusterTPC = fCutMinNClusterTPC;
 343   target.fCutMinNClusterITS = fCutMinNClusterITS;
 344
 345   target.fCutMaxChi2PerClusterTPC = fCutMaxChi2PerClusterTPC;
 346   target.fCutMaxChi2PerClusterITS = fCutMaxChi2PerClusterITS;
 347
 348   for (Int_t i = 0; i < 3; i++)
 349     target.fCutClusterRequirementITS[i] = fCutClusterRequirementITS[i];
 350
 351   target.fCutMaxC11 = fCutMaxC11;
 352   target.fCutMaxC22 = fCutMaxC22;
 353   target.fCutMaxC33 = fCutMaxC33;
 354   target.fCutMaxC44 = fCutMaxC44;
 355   target.fCutMaxC55 = fCutMaxC55;
 356
 357   target.fCutAcceptKinkDaughters = fCutAcceptKinkDaughters;
 358   target.fCutRequireTPCRefit = fCutRequireTPCRefit;
 359   target.fCutRequireITSRefit = fCutRequireITSRefit;
 360
 361   target.fCutNsigmaToVertex = fCutNsigmaToVertex;
 362   target.fCutSigmaToVertexRequired = fCutSigmaToVertexRequired;
 363   target.fCutDCAToVertex = fCutDCAToVertex;
 364   target.fCutDCAToVertexXY = fCutDCAToVertexXY;
 365   target.fCutDCAToVertexZ = fCutDCAToVertexZ;
 366
 367   target.fPMin = fPMin;
 368   target.fPMax = fPMax;
 369   target.fPtMin = fPtMin;
 370   target.fPtMax = fPtMax;
 371   target.fPxMin = fPxMin;
 372   target.fPxMax = fPxMax;
 373   target.fPyMin = fPyMin;
 374   target.fPyMax = fPyMax;
 375   target.fPzMin = fPzMin;
 376   target.fPzMax = fPzMax;
 377   target.fEtaMin = fEtaMin;
 378   target.fEtaMax = fEtaMax;
 379   target.fRapMin = fRapMin;
 380   target.fRapMax = fRapMax;
 381
 382   target.fHistogramsOn = fHistogramsOn;
 383
 384   for (Int_t i=0; i<2; ++i)
 385   {
 386     if (fhNClustersITS[i]) target.fhNClustersITS[i] = (TH1F*) fhNClustersITS[i]->Clone();
 387     if (fhNClustersTPC[i]) target.fhNClustersTPC[i] = (TH1F*) fhNClustersTPC[i]->Clone();
 388
 389     if (fhChi2PerClusterITS[i]) target.fhChi2PerClusterITS[i] = (TH1F*) fhChi2PerClusterITS[i]->Clone();
 390     if (fhChi2PerClusterTPC[i]) target.fhChi2PerClusterTPC[i] = (TH1F*) fhChi2PerClusterTPC[i]->Clone();
 391
 392     if (fhC11[i]) target.fhC11[i] = (TH1F*) fhC11[i]->Clone();
 393     if (fhC22[i]) target.fhC22[i] = (TH1F*) fhC22[i]->Clone();
 394     if (fhC33[i]) target.fhC33[i] = (TH1F*) fhC33[i]->Clone();
 395     if (fhC44[i]) target.fhC44[i] = (TH1F*) fhC44[i]->Clone();
 396     if (fhC55[i]) target.fhC55[i] = (TH1F*) fhC55[i]->Clone();
 397
 398     if (fhDXY[i]) target.fhDXY[i] = (TH1F*) fhDXY[i]->Clone();
 399     if (fhDZ[i]) target.fhDZ[i] = (TH1F*) fhDZ[i]->Clone();
 400     if (fhDXYDZ[i]) target.fhDXYDZ[i] = (TH1F*) fhDXYDZ[i]->Clone();
 401     if (fhDXYvsDZ[i]) target.fhDXYvsDZ[i] = (TH2F*) fhDXYvsDZ[i]->Clone();
 402
 403     if (fhDXYNormalized[i]) target.fhDXYNormalized[i] = (TH1F*) fhDXYNormalized[i]->Clone();
 404     if (fhDZNormalized[i]) target.fhDZNormalized[i] = (TH1F*) fhDZNormalized[i]->Clone();
 405     if (fhDXYvsDZNormalized[i]) target.fhDXYvsDZNormalized[i] = (TH2F*) fhDXYvsDZNormalized[i]->Clone();
 406     if (fhNSigmaToVertex[i]) target.fhNSigmaToVertex[i] = (TH1F*) fhNSigmaToVertex[i]->Clone();
 407
 408     if (fhPt[i]) target.fhPt[i] = (TH1F*) fhPt[i]->Clone();
 409     if (fhEta[i]) target.fhEta[i] = (TH1F*) fhEta[i]->Clone();
 410   }
 411   if (ffDTheoretical) target.ffDTheoretical = (TF1*) ffDTheoretical->Clone();
 412
 413   if (fhCutStatistics) target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
 414   if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();
 415
 416   TNamed::Copy(c);
 417 }
 418
 419 //_____________________________________________________________________________
 420 Long64_t AliESDtrackCuts::Merge(TCollection* list) {
 421   // Merge a list of AliESDtrackCuts objects with this (needed for PROOF)
 422   // Returns the number of merged objects (including this)
 423   if (!list)
 424     return 0;
 425   if (list->IsEmpty())
 426     return 1;
 427   if (!fHistogramsOn)
 428     return 0;
 429   TIterator* iter = list->MakeIterator();
 430   TObject* obj;
 431
 432   // collection of measured and generated histograms
 433   Int_t count = 0;
 434   while ((obj = iter->Next())) {
 435
 436     AliESDtrackCuts* entry = dynamic_cast<AliESDtrackCuts*>(obj);
 437     if (entry == 0)
 438       continue;
 439
 440     if (!entry->fHistogramsOn)
 441       continue;
 442
 443     for (Int_t i=0; i<2; i++) {
 444
 445       fhNClustersITS[i]      ->Add(entry->fhNClustersITS[i]     );
 446       fhNClustersTPC[i]      ->Add(entry->fhNClustersTPC[i]     );
 447
 448       fhChi2PerClusterITS[i] ->Add(entry->fhChi2PerClusterITS[i]);
 449       fhChi2PerClusterTPC[i] ->Add(entry->fhChi2PerClusterTPC[i]);
 450
 451       fhC11[i]               ->Add(entry->fhC11[i]              );
 452       fhC22[i]               ->Add(entry->fhC22[i]              );
 453       fhC33[i]               ->Add(entry->fhC33[i]              );
 454       fhC44[i]               ->Add(entry->fhC44[i]              );
 455       fhC55[i]               ->Add(entry->fhC55[i]              );
 456
 457       fhDXY[i]               ->Add(entry->fhDXY[i]              );
 458       fhDZ[i]                ->Add(entry->fhDZ[i]               );
 459       fhDXYDZ[i]             ->Add(entry->fhDXYDZ[i]          );
 460       fhDXYvsDZ[i]           ->Add(entry->fhDXYvsDZ[i]          );
 461
 462       fhDXYNormalized[i]     ->Add(entry->fhDXYNormalized[i]    );
 463       fhDZNormalized[i]      ->Add(entry->fhDZNormalized[i]     );
 464       fhDXYvsDZNormalized[i] ->Add(entry->fhDXYvsDZNormalized[i]);
 465       fhNSigmaToVertex[i]    ->Add(entry->fhNSigmaToVertex[i]);
 466
 467       fhPt[i]                ->Add(entry->fhPt[i]);
 468       fhEta[i]               ->Add(entry->fhEta[i]);
 469     }
 470
 471     fhCutStatistics  ->Add(entry->fhCutStatistics);
 472     fhCutCorrelation ->Add(entry->fhCutCorrelation);
 473
 474     count++;
 475   }
 476   return count+1;
 477 }
 478
 479
 480 //____________________________________________________________________
 481 Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
 482 {
 483   // Calculates the number of sigma to the vertex.
 484
 485   Float_t b[2];
 486   Float_t bRes[2];
 487   Float_t bCov[3];
 488   esdTrack->GetImpactParameters(b,bCov);
 489
 490   if (bCov[0]<=0 || bCov[2]<=0) {
 491     AliDebugClass(1, "Estimated b resolution lower or equal zero!");
 492     bCov[0]=0; bCov[2]=0;
 493   }
 494   bRes[0] = TMath::Sqrt(bCov[0]);
 495   bRes[1] = TMath::Sqrt(bCov[2]);
 496
 497   // -----------------------------------
 498   // How to get to a n-sigma cut?
 499   //
 500   // The accumulated statistics from 0 to d is
 501   //
 502   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 503   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 504   //
 505   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 506   // Can this be expressed in a different way?
 507
 508   if (bRes[0] == 0 || bRes[1] ==0)
 509     return -1;
 510
 511   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 512
 513   // work around precision problem
 514   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 515   // 1e-15 corresponds to nsigma ~ 7.7
 516   if (TMath::Exp(-d * d / 2) < 1e-15)
 517     return 1000;
 518
 519   Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 520   return nSigma;
 521 }
 522
 523 void AliESDtrackCuts::EnableNeededBranches(TTree* tree)
 524 {
 525   // enables the branches needed by AcceptTrack, for a list see comment of AcceptTrack
 526
 527   tree->SetBranchStatus("fTracks.fFlags", 1);
 528   tree->SetBranchStatus("fTracks.fITSncls", 1);
 529   tree->SetBranchStatus("fTracks.fTPCncls", 1);
 530   tree->SetBranchStatus("fTracks.fITSchi2", 1);
 531   tree->SetBranchStatus("fTracks.fTPCchi2", 1);
 532   tree->SetBranchStatus("fTracks.fC*", 1);
 533   tree->SetBranchStatus("fTracks.fD", 1);
 534   tree->SetBranchStatus("fTracks.fZ", 1);
 535   tree->SetBranchStatus("fTracks.fCdd", 1);
 536   tree->SetBranchStatus("fTracks.fCdz", 1);
 537   tree->SetBranchStatus("fTracks.fCzz", 1);
 538   tree->SetBranchStatus("fTracks.fP*", 1);
 539   tree->SetBranchStatus("fTracks.fR*", 1);
 540   tree->SetBranchStatus("fTracks.fKinkIndexes*", 1);
 541 }
 542
 543 //____________________________________________________________________
 544 Bool_t
 545 AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
 546   //
 547   // figure out if the tracks survives all the track cuts defined
 548   //
 549   // the different quality parameter and kinematic values are first
 550   // retrieved from the track. then it is found out what cuts the
 551   // track did not survive and finally the cuts are imposed.
 552
 553   // this function needs the following branches:
 554   // fTracks.fFlags
 555   // fTracks.fITSncls
 556   // fTracks.fTPCncls
 557   // fTracks.fITSchi2
 558   // fTracks.fTPCchi2
 559   // fTracks.fC   //GetExternalCovariance
 560   // fTracks.fD   //GetImpactParameters
 561   // fTracks.fZ   //GetImpactParameters
 562   // fTracks.fCdd //GetImpactParameters
 563   // fTracks.fCdz //GetImpactParameters
 564   // fTracks.fCzz //GetImpactParameters
 565   // fTracks.fP   //GetPxPyPz
 566   // fTracks.fR   //GetMass
 567   // fTracks.fP   //GetMass
 568   // fTracks.fKinkIndexes
 569
 570
 571   UInt_t status = esdTrack->GetStatus();
 572
 573   // getting quality parameters from the ESD track
 574   Int_t nClustersITS = esdTrack->GetITSclusters(0);
 575   Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
 576
 577   Float_t chi2PerClusterITS = -1;
 578   Float_t chi2PerClusterTPC = -1;
 579   if (nClustersITS!=0)
 580     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
 581   if (nClustersTPC!=0)
 582     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
 583   Double_t extCov[15];
 584   esdTrack->GetExternalCovariance(extCov);
 585
 586   // getting the track to vertex parameters
 587   Float_t nSigmaToVertex = GetSigmaToVertex(esdTrack);
 588
 589   Float_t b[2];
 590   Float_t bCov[3];
 591   esdTrack->GetImpactParameters(b,bCov);
 592   if (bCov[0]<=0 || bCov[2]<=0) {
 593     AliDebug(1, "Estimated b resolution lower or equal zero!");
 594     bCov[0]=0; bCov[2]=0;
 595   }
 596
 597   Float_t dcaToVertexXY = b[0];
 598   Float_t dcaToVertexZ = b[1];
 599
 600   Float_t dcaToVertex   = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
 601
 602   // getting the kinematic variables of the track
 603   // (assuming the mass is known)
 604   Double_t p[3];
 605   esdTrack->GetPxPyPz(p);
 606
 607   Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
 608   Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
 609   Float_t energy   = TMath::Sqrt(TMath::Power(esdTrack->GetMass(),2) + TMath::Power(momentum,2));
 610
 611
 612   //y-eta related calculations
 613   Float_t eta = -100.;
 614   Float_t y   = -100.;
 615   if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
 616     eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
 617   if((energy != TMath::Abs(p[2]))&&(momentum != 0))
 618     y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
 619
 620
 621   //########################################################################
 622   // cut the track?
 623
 624   Bool_t cuts[kNCuts];
 625   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
 626
 627   // track quality cuts
 628   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
 629     cuts[0]=kTRUE;
 630   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
 631     cuts[1]=kTRUE;
 632   if (nClustersTPC<fCutMinNClusterTPC)
 633     cuts[2]=kTRUE;
 634   if (nClustersITS<fCutMinNClusterITS)
 635     cuts[3]=kTRUE;
 636   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
 637     cuts[4]=kTRUE;
 638   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
 639     cuts[5]=kTRUE;
 640   if (extCov[0]  > fCutMaxC11)
 641     cuts[6]=kTRUE;
 642   if (extCov[2]  > fCutMaxC22)
 643     cuts[7]=kTRUE;
 644   if (extCov[5]  > fCutMaxC33)
 645     cuts[8]=kTRUE;
 646   if (extCov[9]  > fCutMaxC44)
 647     cuts[9]=kTRUE;
 648   if (extCov[14]  > fCutMaxC55)
 649     cuts[10]=kTRUE;
 650   if (nSigmaToVertex > fCutNsigmaToVertex && fCutSigmaToVertexRequired)
 651     cuts[11] = kTRUE;
 652   // if n sigma could not be calculated
 653   if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)
 654     cuts[12]=kTRUE;
 655   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
 656     cuts[13]=kTRUE;
 657   // track kinematics cut
 658   if((momentum < fPMin) || (momentum > fPMax))
 659     cuts[14]=kTRUE;
 660   if((pt < fPtMin) || (pt > fPtMax))
 661     cuts[15] = kTRUE;
 662   if((p[0] < fPxMin) || (p[0] > fPxMax))
 663     cuts[16] = kTRUE;
 664   if((p[1] < fPyMin) || (p[1] > fPyMax))
 665     cuts[17] = kTRUE;
 666   if((p[2] < fPzMin) || (p[2] > fPzMax))
 667     cuts[18] = kTRUE;
 668   if((eta < fEtaMin) || (eta > fEtaMax))
 669     cuts[19] = kTRUE;
 670   if((y < fRapMin) || (y > fRapMax))
 671     cuts[20] = kTRUE;
 672   if (dcaToVertex > fCutDCAToVertex)
 673     cuts[21] = kTRUE;
 674   if (TMath::Abs(dcaToVertexXY) > fCutDCAToVertexXY)
 675     cuts[22] = kTRUE;
 676   if (TMath::Abs(dcaToVertexZ) > fCutDCAToVertexZ)
 677     cuts[23] = kTRUE;
 678
 679   for (Int_t i = 0; i < 3; i++)
 680     cuts[24+i] = !CheckITSClusterRequirement(fCutClusterRequirementITS[i], esdTrack->HasPointOnITSLayer(i*2), esdTrack->HasPointOnITSLayer(i*2+1));
 681
 682   Bool_t cut=kFALSE;
 683   for (Int_t i=0; i<kNCuts; i++)
 684     if (cuts[i]) cut = kTRUE;
 685
 686
 687
 688   //########################################################################
 689   // filling histograms
 690   if (fHistogramsOn) {
 691     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
 692     if (cut)
 693       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
 694
 695     for (Int_t i=0; i<kNCuts; i++) {
 696       if (cuts[i])
 697         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
 698
 699       for (Int_t j=i; j<kNCuts; j++) {
 700         if (cuts[i] && cuts[j]) {
 701           Float_t xC = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
 702           Float_t yC = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
 703           fhCutCorrelation->Fill(xC, yC);
 704         }
 705       }
 706     }
 707   }
 708
 709   // now we loop over the filling of the histograms twice: once "before" the cut, once "after"
 710   // the code is not in a function due to too many local variables that would need to be passed
 711
 712   for (Int_t id = 0; id < 2; id++)
 713   {
 714     // id = 0 --> before cut
 715     // id = 1 --> after cut
 716
 717     if (fHistogramsOn)
 718     {
 719       fhNClustersITS[id]->Fill(nClustersITS);
 720       fhNClustersTPC[id]->Fill(nClustersTPC);
 721       fhChi2PerClusterITS[id]->Fill(chi2PerClusterITS);
 722       fhChi2PerClusterTPC[id]->Fill(chi2PerClusterTPC);
 723
 724       fhC11[id]->Fill(extCov[0]);
 725       fhC22[id]->Fill(extCov[2]);
 726       fhC33[id]->Fill(extCov[5]);
 727       fhC44[id]->Fill(extCov[9]);
 728       fhC55[id]->Fill(extCov[14]);
 729
 730       fhPt[id]->Fill(pt);
 731       fhEta[id]->Fill(eta);
 732
 733       Float_t bRes[2];
 734       bRes[0] = TMath::Sqrt(bCov[0]);
 735       bRes[1] = TMath::Sqrt(bCov[2]);
 736
 737       fhDZ[id]->Fill(b[1]);
 738       fhDXY[id]->Fill(b[0]);
 739       fhDXYDZ[id]->Fill(dcaToVertex);
 740       fhDXYvsDZ[id]->Fill(b[1],b[0]);
 741
 742       if (bRes[0]!=0 && bRes[1]!=0) {
 743         fhDZNormalized[id]->Fill(b[1]/bRes[1]);
 744         fhDXYNormalized[id]->Fill(b[0]/bRes[0]);
 745         fhDXYvsDZNormalized[id]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
 746         fhNSigmaToVertex[id]->Fill(nSigmaToVertex);
 747       }
 748     }
 749
 750     // cut the track
 751     if (cut)
 752       return kFALSE;
 753   }
 754
 755   return kTRUE;
 756 }
 757
 758 //____________________________________________________________________
 759 Bool_t AliESDtrackCuts::CheckITSClusterRequirement(ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2)
 760 {
 761   // checks if the cluster requirement is fullfilled (in this case: return kTRUE)
 762
 763   switch (req)
 764   {
 765         case kOff:        return kTRUE;
 766         case kNone:       return !clusterL1 && !clusterL2;
 767         case kAny:        return clusterL1 || clusterL2;
 768         case kFirst:      return clusterL1;
 769         case kOnlyFirst:  return clusterL1 && !clusterL2;
 770         case kSecond:     return clusterL2;
 771         case kOnlySecond: return clusterL2 && !clusterL1;
 772         case kBoth:       return clusterL1 && clusterL2;
 773   }
 774
 775   return kFALSE;
 776 }
 777
 778 //____________________________________________________________________
 779 AliESDtrack* AliESDtrackCuts::GetTPCOnlyTrack(AliESDEvent* esd, Int_t iTrack)
 780 {
 781   // creates a TPC only track from the given esd track
 782   // the track has to be deleted by the user
 783   //
 784   // NB. most of the functionality to get a TPC only track from an ESD track is in AliESDtrack, where it should be
 785   // there are only missing propagations here that are needed for old data
 786   // this function will therefore become obsolete
 787   //
 788   // adapted from code provided by CKB
 789
 790   if (!esd->GetPrimaryVertexTPC())
 791     return 0; // No TPC vertex no TPC tracks
 792
 793   if(!esd->GetPrimaryVertexTPC()->GetStatus())
 794     return 0; // TPC Vertex is created by default in AliESDEvent, do not use in this case
 795
 796   AliESDtrack* track = esd->GetTrack(iTrack);
 797   if (!track)
 798     return 0;
 799
 800   AliESDtrack *tpcTrack = new AliESDtrack();
 801
 802   // This should have been done during the reconstruction
 803   // fixed by Juri in r26675
 804   // but recalculate for older data CKB
 805   Float_t p[2],cov[3];
 806   track->GetImpactParametersTPC(p,cov);
 807   if(p[0]==0&&p[1]==0)
 808     track->RelateToVertexTPC(esd->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig);
 809   // BKC
 810
 811   // only true if we have a tpc track
 812   if (!track->FillTPCOnlyTrack(*tpcTrack))
 813   {
 814     delete tpcTrack;
 815     return 0;
 816   }
 817
 818   // propagate to Vertex
 819   // not needed for normal reconstructed ESDs...
 820   // Double_t pTPC[2],covTPC[3];
 821   // tpcTrack->PropagateToDCA(esd->GetPrimaryVertexTPC(), esd->GetMagneticField(), 10000,  pTPC, covTPC);
 822
 823   return tpcTrack;
 824 }
 825
 826 //____________________________________________________________________
 827 TObjArray* AliESDtrackCuts::GetAcceptedTracks(AliESDEvent* esd,Bool_t bTPC)
 828 {
 829   //
 830   // returns an array of all tracks that pass the cuts
 831   // or an array of TPC only tracks (propagated to the TPC vertex during reco)
 832   // tracks that pass the cut
 833
 834   TObjArray* acceptedTracks = new TObjArray();
 835
 836   // loop over esd tracks
 837   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 838     if(bTPC){
 839       if(!esd->GetPrimaryVertexTPC())return acceptedTracks; // No TPC vertex no TPC tracks
 840       if(!esd->GetPrimaryVertexTPC()->GetStatus())return acceptedTracks; // No proper TPC vertex, only the default
 841
 842       AliESDtrack *tpcTrack = GetTPCOnlyTrack(esd, iTrack);
 843       if (!tpcTrack)
 844         continue;
 845
 846       if (AcceptTrack(tpcTrack)) {
 847         acceptedTracks->Add(tpcTrack);
 848       }
 849       else
 850         delete tpcTrack;
 851     }
 852     else
 853     {
 854       AliESDtrack* track = esd->GetTrack(iTrack);
 855       if(AcceptTrack(track))
 856         acceptedTracks->Add(track);
 857     }
 858   }
 859   if(bTPC)acceptedTracks->SetOwner(kTRUE);
 860   return acceptedTracks;
 861 }
 862
 863 //____________________________________________________________________
 864 Int_t AliESDtrackCuts::CountAcceptedTracks(AliESDEvent* esd)
 865 {
 866   //
 867   // returns an the number of tracks that pass the cuts
 868   //
 869
 870   Int_t count = 0;
 871
 872   // loop over esd tracks
 873   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 874     AliESDtrack* track = esd->GetTrack(iTrack);
 875     if (AcceptTrack(track))
 876       count++;
 877   }
 878
 879   return count;
 880 }
 881
 882 //____________________________________________________________________
 883  void AliESDtrackCuts::DefineHistograms(Int_t color) {
 884    //
 885    // diagnostics histograms are defined
 886    //
 887
 888    fHistogramsOn=kTRUE;
 889
 890    Bool_t oldStatus = TH1::AddDirectoryStatus();
 891    TH1::AddDirectory(kFALSE);
 892
 893    //###################################################################################
 894    // defining histograms
 895
 896    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
 897
 898    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
 899    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
 900
 901    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
 902
 903    for (Int_t i=0; i<kNCuts; i++) {
 904      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
 905      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
 906      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
 907    }
 908
 909   fhCutStatistics  ->SetLineColor(color);
 910   fhCutCorrelation ->SetLineColor(color);
 911   fhCutStatistics  ->SetLineWidth(2);
 912   fhCutCorrelation ->SetLineWidth(2);
 913
 914   for (Int_t i=0; i<2; i++) {
 915     fhNClustersITS[i]        = new TH1F("nClustersITS"    ,"",8,-0.5,7.5);
 916     fhNClustersTPC[i]        = new TH1F("nClustersTPC"     ,"",165,-0.5,164.5);
 917     fhChi2PerClusterITS[i]   = new TH1F("chi2PerClusterITS","",500,0,10);
 918     fhChi2PerClusterTPC[i]   = new TH1F("chi2PerClusterTPC","",500,0,10);
 919
 920     fhC11[i]                 = new TH1F("covMatrixDiagonal11","",2000,0,20);
 921     fhC22[i]                 = new TH1F("covMatrixDiagonal22","",2000,0,20);
 922     fhC33[i]                 = new TH1F("covMatrixDiagonal33","",1000,0,0.1);
 923     fhC44[i]                 = new TH1F("covMatrixDiagonal44","",1000,0,0.1);
 924     fhC55[i]                 = new TH1F("covMatrixDiagonal55","",1000,0,5);
 925
 926     fhDXY[i]                 = new TH1F("dXY"    ,"",500,-10,10);
 927     fhDZ[i]                  = new TH1F("dZ"     ,"",500,-10,10);
 928     fhDXYDZ[i]               = new TH1F("dXYDZ"  ,"",500,0,10);
 929     fhDXYvsDZ[i]             = new TH2F("dXYvsDZ","",200,-10,10,200,-10,10);
 930
 931     fhDXYNormalized[i]       = new TH1F("dXYNormalized"    ,"",500,-10,10);
 932     fhDZNormalized[i]        = new TH1F("dZNormalized"     ,"",500,-10,10);
 933     fhDXYvsDZNormalized[i]   = new TH2F("dXYvsDZNormalized","",200,-10,10,200,-10,10);
 934
 935     fhNSigmaToVertex[i]      = new TH1F("nSigmaToVertex","",500,0,10);
 936
 937     fhPt[i]                  = new TH1F("pt"     ,"p_{T} distribution;p_{T} (GeV/c)",500,0.0,100.0);
 938     fhEta[i]                 = new TH1F("eta"     ,"#eta distribution;#eta",40,-2.0,2.0);
 939
 940     fhNClustersITS[i]->SetTitle("n ITS clusters");
 941     fhNClustersTPC[i]->SetTitle("n TPC clusters");
 942     fhChi2PerClusterITS[i]->SetTitle("#Chi^{2} per ITS cluster");
 943     fhChi2PerClusterTPC[i]->SetTitle("#Chi^{2} per TPC cluster");
 944
 945     fhC11[i]->SetTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
 946     fhC22[i]->SetTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
 947     fhC33[i]->SetTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
 948     fhC44[i]->SetTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
 949     fhC55[i]->SetTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
 950
 951     fhDXY[i]->SetTitle("transverse impact parameter");
 952     fhDZ[i]->SetTitle("longitudinal impact parameter");
 953     fhDXYDZ[i]->SetTitle("absolute impact parameter;sqrt(dXY**2 + dZ**2) in cm");
 954     fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter");
 955     fhDXYvsDZ[i]->SetYTitle("transverse impact parameter");
 956
 957     fhDXYNormalized[i]->SetTitle("normalized trans impact par");
 958     fhDZNormalized[i]->SetTitle("normalized long impact par");
 959     fhDXYvsDZNormalized[i]->SetTitle("normalized long impact par");
 960     fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par");
 961     fhNSigmaToVertex[i]->SetTitle("n #sigma to vertex");
 962
 963     fhNClustersITS[i]->SetLineColor(color);   fhNClustersITS[i]->SetLineWidth(2);
 964     fhNClustersTPC[i]->SetLineColor(color);   fhNClustersTPC[i]->SetLineWidth(2);
 965     fhChi2PerClusterITS[i]->SetLineColor(color);   fhChi2PerClusterITS[i]->SetLineWidth(2);
 966     fhChi2PerClusterTPC[i]->SetLineColor(color);   fhChi2PerClusterTPC[i]->SetLineWidth(2);
 967
 968     fhC11[i]->SetLineColor(color);   fhC11[i]->SetLineWidth(2);
 969     fhC22[i]->SetLineColor(color);   fhC22[i]->SetLineWidth(2);
 970     fhC33[i]->SetLineColor(color);   fhC33[i]->SetLineWidth(2);
 971     fhC44[i]->SetLineColor(color);   fhC44[i]->SetLineWidth(2);
 972     fhC55[i]->SetLineColor(color);   fhC55[i]->SetLineWidth(2);
 973
 974     fhDXY[i]->SetLineColor(color);   fhDXY[i]->SetLineWidth(2);
 975     fhDZ[i]->SetLineColor(color);    fhDZ[i]->SetLineWidth(2);
 976     fhDXYDZ[i]->SetLineColor(color); fhDXYDZ[i]->SetLineWidth(2);
 977
 978     fhDXYNormalized[i]->SetLineColor(color);   fhDXYNormalized[i]->SetLineWidth(2);
 979     fhDZNormalized[i]->SetLineColor(color);    fhDZNormalized[i]->SetLineWidth(2);
 980     fhNSigmaToVertex[i]->SetLineColor(color);  fhNSigmaToVertex[i]->SetLineWidth(2);
 981   }
 982
 983   // The number of sigmas to the vertex is per definition gaussian
 984   ffDTheoretical = new TF1("nSigmaToVertexTheoretical","([0]/2.506628274)*exp(-(x**2)/2)",0,50);
 985   ffDTheoretical->SetParameter(0,1);
 986
 987   TH1::AddDirectory(oldStatus);
 988 }
 989
 990 //____________________________________________________________________
 991 Bool_t AliESDtrackCuts::LoadHistograms(const Char_t* dir)
 992 {
 993   //
 994   // loads the histograms from a file
 995   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
 996   //
 997
 998   if (!dir)
 999     dir = GetName();
1000
1001   if (!gDirectory->cd(dir))
1002     return kFALSE;
1003
1004   ffDTheoretical = dynamic_cast<TF1*> (gDirectory->Get("nSigmaToVertexTheory"));
1005
1006   fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
1007   fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
1008
1009   for (Int_t i=0; i<2; i++) {
1010     if (i==0)
1011     {
1012       gDirectory->cd("before_cuts");
1013     }
1014     else
1015       gDirectory->cd("after_cuts");
1016
1017     fhNClustersITS[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersITS"     ));
1018     fhNClustersTPC[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersTPC"     ));
1019     fhChi2PerClusterITS[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterITS"));
1020     fhChi2PerClusterTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterTPC"));
1021
1022     fhC11[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal11"));
1023     fhC22[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal22"));
1024     fhC33[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal33"));
1025     fhC44[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal44"));
1026     fhC55[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal55"));
1027
1028     fhDXY[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXY"    ));
1029     fhDZ[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZ"     ));
1030     fhDXYDZ[i] =   dynamic_cast<TH1F*> (gDirectory->Get("dXYDZ"));
1031     fhDXYvsDZ[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZ"));
1032
1033     fhDXYNormalized[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXYNormalized"    ));
1034     fhDZNormalized[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZNormalized"     ));
1035     fhDXYvsDZNormalized[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZNormalized"));
1036     fhNSigmaToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSigmaToVertex"));
1037
1038     fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get("pt"));
1039     fhEta[i] = dynamic_cast<TH1F*> (gDirectory->Get("eta"));
1040
1041     gDirectory->cd("../");
1042   }
1043
1044   gDirectory->cd("..");
1045
1046   return kTRUE;
1047 }
1048
1049 //____________________________________________________________________
1050 void AliESDtrackCuts::SaveHistograms(const Char_t* dir) {
1051   //
1052   // saves the histograms in a directory (dir)
1053   //
1054
1055   if (!fHistogramsOn) {
1056     AliDebug(0, "Histograms not on - cannot save histograms!!!");
1057     return;
1058   }
1059
1060   if (!dir)
1061     dir = GetName();
1062
1063   gDirectory->mkdir(dir);
1064   gDirectory->cd(dir);
1065
1066   gDirectory->mkdir("before_cuts");
1067   gDirectory->mkdir("after_cuts");
1068
1069   // a factor of 2 is needed since n sigma is positive
1070   ffDTheoretical->SetParameter(0,2*fhNSigmaToVertex[0]->Integral("width"));
1071   ffDTheoretical->Write("nSigmaToVertexTheory");
1072
1073   fhCutStatistics->Write();
1074   fhCutCorrelation->Write();
1075
1076   for (Int_t i=0; i<2; i++) {
1077     if (i==0)
1078       gDirectory->cd("before_cuts");
1079     else
1080       gDirectory->cd("after_cuts");
1081
1082     fhNClustersITS[i]        ->Write();
1083     fhNClustersTPC[i]        ->Write();
1084     fhChi2PerClusterITS[i]   ->Write();
1085     fhChi2PerClusterTPC[i]   ->Write();
1086
1087     fhC11[i]                 ->Write();
1088     fhC22[i]                 ->Write();
1089     fhC33[i]                 ->Write();
1090     fhC44[i]                 ->Write();
1091     fhC55[i]                 ->Write();
1092
1093     fhDXY[i]                 ->Write();
1094     fhDZ[i]                  ->Write();
1095     fhDXYDZ[i]               ->Write();
1096     fhDXYvsDZ[i]             ->Write();
1097
1098     fhDXYNormalized[i]       ->Write();
1099     fhDZNormalized[i]        ->Write();
1100     fhDXYvsDZNormalized[i]   ->Write();
1101     fhNSigmaToVertex[i]      ->Write();
1102
1103     fhPt[i]                  ->Write();
1104     fhEta[i]                 ->Write();
1105
1106     gDirectory->cd("../");
1107   }
1108
1109   gDirectory->cd("../");
1110 }
1111
1112 //____________________________________________________________________
1113 void AliESDtrackCuts::DrawHistograms()
1114 {
1115   // draws some histograms
1116
1117   TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "Track Quality Results1", 800, 800);
1118   canvas1->Divide(2, 2);
1119
1120   canvas1->cd(1);
1121   fhNClustersTPC[0]->SetStats(kFALSE);
1122   fhNClustersTPC[0]->Draw();
1123
1124   canvas1->cd(2);
1125   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1126   fhChi2PerClusterTPC[0]->Draw();
1127
1128   canvas1->cd(3);
1129   fhNSigmaToVertex[0]->SetStats(kFALSE);
1130   fhNSigmaToVertex[0]->GetXaxis()->SetRangeUser(0, 10);
1131   fhNSigmaToVertex[0]->Draw();
1132
1133   canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
1134
1135   TCanvas* canvas2 = new TCanvas(Form("%s_2", GetName()), "Track Quality Results2", 1200, 800);
1136   canvas2->Divide(3, 2);
1137
1138   canvas2->cd(1);
1139   fhC11[0]->SetStats(kFALSE);
1140   gPad->SetLogy();
1141   fhC11[0]->Draw();
1142
1143   canvas2->cd(2);
1144   fhC22[0]->SetStats(kFALSE);
1145   gPad->SetLogy();
1146   fhC22[0]->Draw();
1147
1148   canvas2->cd(3);
1149   fhC33[0]->SetStats(kFALSE);
1150   gPad->SetLogy();
1151   fhC33[0]->Draw();
1152
1153   canvas2->cd(4);
1154   fhC44[0]->SetStats(kFALSE);
1155   gPad->SetLogy();
1156   fhC44[0]->Draw();
1157
1158   canvas2->cd(5);
1159   fhC55[0]->SetStats(kFALSE);
1160   gPad->SetLogy();
1161   fhC55[0]->Draw();
1162
1163   canvas2->SaveAs(Form("%s_%s.gif", GetName(), canvas2->GetName()));
1164
1165   TCanvas* canvas3 = new TCanvas(Form("%s_3", GetName()), "Track Quality Results3", 1200, 800);
1166   canvas3->Divide(3, 2);
1167
1168   canvas3->cd(1);
1169   fhDXY[0]->SetStats(kFALSE);
1170   gPad->SetLogy();
1171   fhDXY[0]->Draw();
1172
1173   canvas3->cd(2);
1174   fhDZ[0]->SetStats(kFALSE);
1175   gPad->SetLogy();
1176   fhDZ[0]->Draw();
1177
1178   canvas3->cd(3);
1179   fhDXYvsDZ[0]->SetStats(kFALSE);
1180   gPad->SetLogz();
1181   gPad->SetRightMargin(0.15);
1182   fhDXYvsDZ[0]->Draw("COLZ");
1183
1184   canvas3->cd(4);
1185   fhDXYNormalized[0]->SetStats(kFALSE);
1186   gPad->SetLogy();
1187   fhDXYNormalized[0]->Draw();
1188
1189   canvas3->cd(5);
1190   fhDZNormalized[0]->SetStats(kFALSE);
1191   gPad->SetLogy();
1192   fhDZNormalized[0]->Draw();
1193
1194   canvas3->cd(6);
1195   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1196   gPad->SetLogz();
1197   gPad->SetRightMargin(0.15);
1198   fhDXYvsDZNormalized[0]->Draw("COLZ");
1199
1200   canvas3->SaveAs(Form("%s_%s.gif", GetName(), canvas3->GetName()));
1201
1202   TCanvas* canvas4 = new TCanvas(Form("%s_4", GetName()), "Track Quality Results4", 800, 500);
1203   canvas4->Divide(2, 1);
1204
1205   canvas4->cd(1);
1206   fhCutStatistics->SetStats(kFALSE);
1207   fhCutStatistics->LabelsOption("v");
1208   gPad->SetBottomMargin(0.3);
1209   fhCutStatistics->Draw();
1210
1211   canvas4->cd(2);
1212   fhCutCorrelation->SetStats(kFALSE);
1213   fhCutCorrelation->LabelsOption("v");
1214   gPad->SetBottomMargin(0.3);
1215   gPad->SetLeftMargin(0.3);
1216   fhCutCorrelation->Draw("COLZ");
1217
1218   canvas4->SaveAs(Form("%s_%s.gif", GetName(), canvas4->GetName()));
1219
1220   /*canvas->cd(1);
1221   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1222   fhDXYvsDZNormalized[0]->DrawCopy("COLZ");
1223
1224   canvas->cd(2);
1225   fhNClustersTPC[0]->SetStats(kFALSE);
1226   fhNClustersTPC[0]->DrawCopy();
1227
1228   canvas->cd(3);
1229   fhChi2PerClusterITS[0]->SetStats(kFALSE);
1230   fhChi2PerClusterITS[0]->DrawCopy();
1231   fhChi2PerClusterITS[1]->SetLineColor(2);
1232   fhChi2PerClusterITS[1]->DrawCopy("SAME");
1233
1234   canvas->cd(4);
1235   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1236   fhChi2PerClusterTPC[0]->DrawCopy();
1237   fhChi2PerClusterTPC[1]->SetLineColor(2);
1238   fhChi2PerClusterTPC[1]->DrawCopy("SAME");*/
1239 }
1240
1241 Float_t AliESDtrackCuts::GetMinNsigmaToVertex() const
1242 {
1243   // deprecated, please use GetMaxNsigmaToVertex
1244
1245   Printf("WARNING: AliESDtrackCuts::GetMinNsigmaToVertex is DEPRECATED and will be removed in the next release. Please use GetMaxNsigmaToVertex instead. Renaming was done to improve code readability.");
1246
1247   return GetMaxNsigmaToVertex();
1248 }
1249
1250 void AliESDtrackCuts::SetMinNsigmaToVertex(Float_t sigma)
1251 {
1252   // deprecated, will be removed in next release
1253
1254   Printf("WARNING: AliESDtrackCuts::SetMinNsigmaToVertex is DEPRECATED and will be removed in the next release. Please use SetMaxNsigmaToVertex instead. Renaming was done to improve code readability.");
1255
1256   SetMaxNsigmaToVertex(sigma);
1257 }
1258
1259 void AliESDtrackCuts::SetDCAToVertex(Float_t dist)
1260 {
1261   // deprecated, will be removed in next release
1262
1263   Printf("WARNING: AliESDtrackCuts::SetDCAToVertex is DEPRECATED and will be removed in the next release. Please use SetMaxDCAToVertex instead. Renaming was done to improve code readability.");
1264
1265   SetMaxDCAToVertex(dist);
1266 }
1267
1268 void AliESDtrackCuts::SetDCAToVertexXY(Float_t dist)
1269 {
1270   // deprecated, will be removed in next release
1271
1272   Printf("WARNING: AliESDtrackCuts::SetDCAToVertexXY is DEPRECATED and will be removed in the next release. Please use SetMaxDCAToVertexXY instead. Renaming was done to improve code readability.");
1273
1274   SetMaxDCAToVertexXY(dist);
1275 }
1276
1277 void AliESDtrackCuts::SetDCAToVertexZ(Float_t dist)
1278 {
1279   // deprecated, will be removed in next release
1280
1281   Printf("WARNING: AliESDtrackCuts::SetDCAToVertexZ is DEPRECATED and will be removed in the next release. Please use SetMaxDCAToVertexZ instead. Renaming was done to improve code readability.");
1282
1283   SetMaxDCAToVertexZ(dist);
1284 }