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