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 TPC standalone",
  38  "require ITS refit",
  39  "n clusters TPC",
  40  "n clusters ITS",
  41  "#Chi^{2}/cluster TPC",
  42  "#Chi^{2}/cluster ITS",
  43  "cov 11",
  44  "cov 22",
  45  "cov 33",
  46  "cov 44",
  47  "cov 55",
  48  "trk-to-vtx",
  49  "trk-to-vtx failed",
  50  "kink daughters",
  51  "p",
  52  "p_{T}",
  53  "p_{x}",
  54  "p_{y}",
  55  "p_{z}",
  56  "eta",
  57  "y",
  58  "trk-to-vtx max dca 2D absolute",
  59  "trk-to-vtx max dca xy absolute",
  60  "trk-to-vtx max dca z absolute",
  61  "trk-to-vtx min dca 2D absolute",
  62  "trk-to-vtx min dca xy absolute",
  63  "trk-to-vtx min dca z absolute",
  64  "SPD cluster requirement",
  65  "SDD cluster requirement",
  66  "SSD cluster requirement",
  67  "require ITS stand-alone",
  68  "rel 1/pt uncertainty",
  69  "TPC n shared clusters",
  70  "TPC rel shared clusters",
  71  "require ITS Pid",
  72  "n crossed rows TPC",
  73  "n crossed rows / n findable clusters",
  74  "missing ITS points"
  75 };
  76
  77 //____________________________________________________________________
  78 AliESDtrackCuts::AliESDtrackCuts(const Char_t* name, const Char_t* title) : AliAnalysisCuts(name,title),
  79   fCutMinNClusterTPC(0),
  80   fCutMinNClusterITS(0),
  81   fCutMinNCrossedRowsTPC(0),
  82   fCutMinRatioCrossedRowsOverFindableClustersTPC(0),
  83   fCutMaxChi2PerClusterTPC(0),
  84   fCutMaxChi2PerClusterITS(0),
  85   fCutMaxMissingITSPoints(0),
  86   fCutMaxC11(0),
  87   fCutMaxC22(0),
  88   fCutMaxC33(0),
  89   fCutMaxC44(0),
  90   fCutMaxC55(0),
  91   fCutMaxRel1PtUncertainty(0),
  92   fCutAcceptKinkDaughters(0),
  93   fCutAcceptSharedTPCClusters(0),
  94   fCutMaxFractionSharedTPCClusters(0),
  95   fCutRequireTPCRefit(0),
  96   fCutRequireTPCStandAlone(0),
  97   fCutRequireITSRefit(0),
  98   fCutRequireITSPid(0),
  99   fCutRequireITSStandAlone(0),
 100   fCutRequireITSpureSA(0),
 101   fCutNsigmaToVertex(0),
 102   fCutSigmaToVertexRequired(0),
 103   fCutMaxDCAToVertexXY(0),
 104   fCutMaxDCAToVertexZ(0),
 105   fCutMinDCAToVertexXY(0),
 106   fCutMinDCAToVertexZ(0),
 107   fCutMaxDCAToVertexXYPtDep(""),
 108   fCutMaxDCAToVertexZPtDep(""),
 109   fCutMinDCAToVertexXYPtDep(""),
 110   fCutMinDCAToVertexZPtDep(""),
 111   f1CutMaxDCAToVertexXYPtDep(0x0),
 112   f1CutMaxDCAToVertexZPtDep(0x0),
 113   f1CutMinDCAToVertexXYPtDep(0x0),
 114   f1CutMinDCAToVertexZPtDep(0x0),
 115   fCutDCAToVertex2D(0),
 116   fPMin(0),
 117   fPMax(0),
 118   fPtMin(0),
 119   fPtMax(0),
 120   fPxMin(0),
 121   fPxMax(0),
 122   fPyMin(0),
 123   fPyMax(0),
 124   fPzMin(0),
 125   fPzMax(0),
 126   fEtaMin(0),
 127   fEtaMax(0),
 128   fRapMin(0),
 129   fRapMax(0),
 130   fHistogramsOn(0),
 131   ffDTheoretical(0),
 132   fhCutStatistics(0),
 133   fhCutCorrelation(0)
 134 {
 135   //
 136   // constructor
 137   //
 138
 139   Init();
 140
 141   //##############################################################################
 142   // setting default cuts
 143   SetMinNClustersTPC();
 144   SetMinNClustersITS();
 145   SetMinNCrossedRowsTPC();
 146   SetMinRatioCrossedRowsOverFindableClustersTPC();
 147   SetMaxChi2PerClusterTPC();
 148   SetMaxChi2PerClusterITS();
 149   SetMaxNOfMissingITSPoints();
 150   SetMaxCovDiagonalElements();
 151   SetMaxRel1PtUncertainty();
 152   SetRequireTPCRefit();
 153   SetRequireTPCStandAlone();
 154   SetRequireITSRefit();
 155   SetRequireITSPid(kFALSE);
 156   SetRequireITSStandAlone(kFALSE);
 157   SetRequireITSPureStandAlone(kFALSE);
 158   SetAcceptKinkDaughters();
 159   SetAcceptSharedTPCClusters();
 160   SetMaxFractionSharedTPCClusters();
 161   SetMaxNsigmaToVertex();
 162   SetMaxDCAToVertexXY();
 163   SetMaxDCAToVertexZ();
 164   SetDCAToVertex2D();
 165   SetMinDCAToVertexXY();
 166   SetMinDCAToVertexZ();
 167   SetPRange();
 168   SetPtRange();
 169   SetPxRange();
 170   SetPyRange();
 171   SetPzRange();
 172   SetEtaRange();
 173   SetRapRange();
 174   SetClusterRequirementITS(kSPD);
 175   SetClusterRequirementITS(kSDD);
 176   SetClusterRequirementITS(kSSD);
 177
 178   SetHistogramsOn();
 179 }
 180
 181 //_____________________________________________________________________________
 182 AliESDtrackCuts::AliESDtrackCuts(const AliESDtrackCuts &c) : AliAnalysisCuts(c),
 183   fCutMinNClusterTPC(0),
 184   fCutMinNClusterITS(0),
 185   fCutMinNCrossedRowsTPC(0),
 186   fCutMinRatioCrossedRowsOverFindableClustersTPC(0),
 187   fCutMaxChi2PerClusterTPC(0),
 188   fCutMaxChi2PerClusterITS(0),
 189   fCutMaxMissingITSPoints(0),
 190   fCutMaxC11(0),
 191   fCutMaxC22(0),
 192   fCutMaxC33(0),
 193   fCutMaxC44(0),
 194   fCutMaxC55(0),
 195   fCutMaxRel1PtUncertainty(0),
 196   fCutAcceptKinkDaughters(0),
 197   fCutAcceptSharedTPCClusters(0),
 198   fCutMaxFractionSharedTPCClusters(0),
 199   fCutRequireTPCRefit(0),
 200   fCutRequireTPCStandAlone(0),
 201   fCutRequireITSRefit(0),
 202   fCutRequireITSPid(0),
 203   fCutRequireITSStandAlone(0),
 204   fCutRequireITSpureSA(0),
 205   fCutNsigmaToVertex(0),
 206   fCutSigmaToVertexRequired(0),
 207   fCutMaxDCAToVertexXY(0),
 208   fCutMaxDCAToVertexZ(0),
 209   fCutMinDCAToVertexXY(0),
 210   fCutMinDCAToVertexZ(0),
 211   fCutMaxDCAToVertexXYPtDep(""),
 212   fCutMaxDCAToVertexZPtDep(""),
 213   fCutMinDCAToVertexXYPtDep(""),
 214   fCutMinDCAToVertexZPtDep(""),
 215   f1CutMaxDCAToVertexXYPtDep(0x0),
 216   f1CutMaxDCAToVertexZPtDep(0x0),
 217   f1CutMinDCAToVertexXYPtDep(0x0),
 218   f1CutMinDCAToVertexZPtDep(0x0),
 219   fCutDCAToVertex2D(0),
 220   fPMin(0),
 221   fPMax(0),
 222   fPtMin(0),
 223   fPtMax(0),
 224   fPxMin(0),
 225   fPxMax(0),
 226   fPyMin(0),
 227   fPyMax(0),
 228   fPzMin(0),
 229   fPzMax(0),
 230   fEtaMin(0),
 231   fEtaMax(0),
 232   fRapMin(0),
 233   fRapMax(0),
 234   fHistogramsOn(0),
 235   ffDTheoretical(0),
 236   fhCutStatistics(0),
 237   fhCutCorrelation(0)
 238 {
 239   //
 240   // copy constructor
 241   //
 242
 243   ((AliESDtrackCuts &) c).Copy(*this);
 244 }
 245
 246 AliESDtrackCuts::~AliESDtrackCuts()
 247 {
 248   //
 249   // destructor
 250   //
 251
 252   for (Int_t i=0; i<2; i++) {
 253
 254     if (fhNClustersITS[i])
 255       delete fhNClustersITS[i];
 256     if (fhNClustersTPC[i])
 257       delete fhNClustersTPC[i];
 258     if (fhNSharedClustersTPC[i])
 259       delete fhNSharedClustersTPC[i];
 260     if (fhNCrossedRowsTPC[i])
 261       delete fhNCrossedRowsTPC[i];
 262     if (fhRatioCrossedRowsOverFindableClustersTPC[i])
 263       delete fhRatioCrossedRowsOverFindableClustersTPC[i];
 264     if (fhChi2PerClusterITS[i])
 265       delete fhChi2PerClusterITS[i];
 266     if (fhChi2PerClusterTPC[i])
 267       delete fhChi2PerClusterTPC[i];
 268     if (fhC11[i])
 269       delete fhC11[i];
 270     if (fhC22[i])
 271       delete fhC22[i];
 272     if (fhC33[i])
 273       delete fhC33[i];
 274     if (fhC44[i])
 275       delete fhC44[i];
 276     if (fhC55[i])
 277       delete fhC55[i];
 278
 279     if (fhRel1PtUncertainty[i])
 280       delete fhRel1PtUncertainty[i];
 281
 282     if (fhDXY[i])
 283       delete fhDXY[i];
 284     if (fhDZ[i])
 285       delete fhDZ[i];
 286     if (fhDXYDZ[i])
 287       delete fhDXYDZ[i];
 288     if (fhDXYvsDZ[i])
 289       delete fhDXYvsDZ[i];
 290
 291     if (fhDXYNormalized[i])
 292       delete fhDXYNormalized[i];
 293     if (fhDZNormalized[i])
 294       delete fhDZNormalized[i];
 295     if (fhDXYvsDZNormalized[i])
 296       delete fhDXYvsDZNormalized[i];
 297     if (fhNSigmaToVertex[i])
 298       delete fhNSigmaToVertex[i];
 299     if (fhPt[i])
 300       delete fhPt[i];
 301     if (fhEta[i])
 302       delete fhEta[i];
 303   }
 304
 305   if(f1CutMaxDCAToVertexXYPtDep)delete f1CutMaxDCAToVertexXYPtDep;
 306   f1CutMaxDCAToVertexXYPtDep = 0;
 307   if( f1CutMaxDCAToVertexZPtDep) delete  f1CutMaxDCAToVertexZPtDep;
 308   f1CutMaxDCAToVertexZPtDep = 0;
 309   if( f1CutMinDCAToVertexXYPtDep)delete  f1CutMinDCAToVertexXYPtDep;
 310   f1CutMinDCAToVertexXYPtDep = 0;
 311   if(f1CutMinDCAToVertexZPtDep)delete  f1CutMinDCAToVertexZPtDep;
 312   f1CutMinDCAToVertexZPtDep = 0;
 313
 314
 315   if (ffDTheoretical)
 316     delete ffDTheoretical;
 317
 318   if (fhCutStatistics)
 319     delete fhCutStatistics;
 320   if (fhCutCorrelation)
 321     delete fhCutCorrelation;
 322 }
 323
 324 void AliESDtrackCuts::Init()
 325 {
 326   //
 327   // sets everything to zero
 328   //
 329
 330   fCutMinNClusterTPC = 0;
 331   fCutMinNClusterITS = 0;
 332
 333   fCutMaxChi2PerClusterTPC = 0;
 334   fCutMaxChi2PerClusterITS = 0;
 335   fCutMaxMissingITSPoints  = 0;
 336
 337   for (Int_t i = 0; i < 3; i++)
 338         fCutClusterRequirementITS[i] = kOff;
 339
 340   fCutMaxC11 = 0;
 341   fCutMaxC22 = 0;
 342   fCutMaxC33 = 0;
 343   fCutMaxC44 = 0;
 344   fCutMaxC55 = 0;
 345
 346   fCutMaxRel1PtUncertainty = 0;
 347
 348   fCutAcceptKinkDaughters = 0;
 349   fCutAcceptSharedTPCClusters = 0;
 350   fCutMaxFractionSharedTPCClusters = 0;
 351   fCutRequireTPCRefit = 0;
 352   fCutRequireTPCStandAlone = 0;
 353   fCutRequireITSRefit = 0;
 354   fCutRequireITSPid = 0;
 355   fCutRequireITSStandAlone = 0;
 356   fCutRequireITSpureSA = 0;
 357
 358   fCutNsigmaToVertex = 0;
 359   fCutSigmaToVertexRequired = 0;
 360   fCutMaxDCAToVertexXY = 0;
 361   fCutMaxDCAToVertexZ = 0;
 362   fCutDCAToVertex2D = 0;
 363   fCutMinDCAToVertexXY = 0;
 364   fCutMinDCAToVertexZ = 0;
 365   fCutMaxDCAToVertexXYPtDep = "";
 366   fCutMaxDCAToVertexZPtDep = "";
 367   fCutMinDCAToVertexXYPtDep = "";
 368   fCutMinDCAToVertexZPtDep = "";
 369
 370   if(f1CutMaxDCAToVertexXYPtDep)delete f1CutMaxDCAToVertexXYPtDep;
 371   f1CutMaxDCAToVertexXYPtDep = 0;
 372   if( f1CutMaxDCAToVertexXYPtDep) delete  f1CutMaxDCAToVertexXYPtDep;
 373   f1CutMaxDCAToVertexXYPtDep = 0;
 374   if( f1CutMaxDCAToVertexZPtDep) delete  f1CutMaxDCAToVertexZPtDep;
 375   f1CutMaxDCAToVertexZPtDep = 0;
 376   if( f1CutMinDCAToVertexXYPtDep)delete  f1CutMinDCAToVertexXYPtDep;
 377   f1CutMinDCAToVertexXYPtDep = 0;
 378   if(f1CutMinDCAToVertexZPtDep)delete f1CutMinDCAToVertexZPtDep;
 379   f1CutMinDCAToVertexZPtDep = 0;
 380
 381
 382   fPMin = 0;
 383   fPMax = 0;
 384   fPtMin = 0;
 385   fPtMax = 0;
 386   fPxMin = 0;
 387   fPxMax = 0;
 388   fPyMin = 0;
 389   fPyMax = 0;
 390   fPzMin = 0;
 391   fPzMax = 0;
 392   fEtaMin = 0;
 393   fEtaMax = 0;
 394   fRapMin = 0;
 395   fRapMax = 0;
 396
 397   fHistogramsOn = kFALSE;
 398
 399   for (Int_t i=0; i<2; ++i)
 400   {
 401     fhNClustersITS[i] = 0;
 402     fhNClustersTPC[i] = 0;
 403     fhNSharedClustersTPC[i] = 0;
 404     fhNCrossedRowsTPC[i] = 0;
 405     fhRatioCrossedRowsOverFindableClustersTPC[i] = 0;
 406
 407     fhChi2PerClusterITS[i] = 0;
 408     fhChi2PerClusterTPC[i] = 0;
 409
 410     fhC11[i] = 0;
 411     fhC22[i] = 0;
 412     fhC33[i] = 0;
 413     fhC44[i] = 0;
 414     fhC55[i] = 0;
 415
 416     fhRel1PtUncertainty[i] = 0;
 417
 418     fhDXY[i] = 0;
 419     fhDZ[i] = 0;
 420     fhDXYDZ[i] = 0;
 421     fhDXYvsDZ[i] = 0;
 422
 423     fhDXYNormalized[i] = 0;
 424     fhDZNormalized[i] = 0;
 425     fhDXYvsDZNormalized[i] = 0;
 426     fhNSigmaToVertex[i] = 0;
 427
 428     fhPt[i] = 0;
 429     fhEta[i] = 0;
 430   }
 431   ffDTheoretical = 0;
 432
 433   fhCutStatistics = 0;
 434   fhCutCorrelation = 0;
 435 }
 436
 437 //_____________________________________________________________________________
 438 AliESDtrackCuts &AliESDtrackCuts::operator=(const AliESDtrackCuts &c)
 439 {
 440   //
 441   // Assignment operator
 442   //
 443
 444   if (this != &c) ((AliESDtrackCuts &) c).Copy(*this);
 445   return *this;
 446 }
 447
 448 //_____________________________________________________________________________
 449 void AliESDtrackCuts::Copy(TObject &c) const
 450 {
 451   //
 452   // Copy function
 453   //
 454
 455   AliESDtrackCuts& target = (AliESDtrackCuts &) c;
 456
 457   target.Init();
 458
 459   target.fCutMinNClusterTPC = fCutMinNClusterTPC;
 460   target.fCutMinNClusterITS = fCutMinNClusterITS;
 461   target.fCutMinNCrossedRowsTPC = fCutMinNCrossedRowsTPC;
 462   target.fCutMinRatioCrossedRowsOverFindableClustersTPC = fCutMinRatioCrossedRowsOverFindableClustersTPC;
 463
 464
 465   target.fCutMaxChi2PerClusterTPC = fCutMaxChi2PerClusterTPC;
 466   target.fCutMaxChi2PerClusterITS = fCutMaxChi2PerClusterITS;
 467   target.fCutMaxMissingITSPoints = fCutMaxMissingITSPoints;
 468
 469   for (Int_t i = 0; i < 3; i++)
 470     target.fCutClusterRequirementITS[i] = fCutClusterRequirementITS[i];
 471
 472   target.fCutMaxC11 = fCutMaxC11;
 473   target.fCutMaxC22 = fCutMaxC22;
 474   target.fCutMaxC33 = fCutMaxC33;
 475   target.fCutMaxC44 = fCutMaxC44;
 476   target.fCutMaxC55 = fCutMaxC55;
 477
 478   target.fCutMaxRel1PtUncertainty = fCutMaxRel1PtUncertainty;
 479
 480   target.fCutAcceptKinkDaughters = fCutAcceptKinkDaughters;
 481   target.fCutAcceptSharedTPCClusters = fCutAcceptSharedTPCClusters;
 482   target.fCutMaxFractionSharedTPCClusters = fCutMaxFractionSharedTPCClusters;
 483   target.fCutRequireTPCRefit = fCutRequireTPCRefit;
 484   target.fCutRequireTPCStandAlone = fCutRequireTPCStandAlone;
 485   target.fCutRequireITSRefit = fCutRequireITSRefit;
 486   target.fCutRequireITSPid = fCutRequireITSPid;
 487   target.fCutRequireITSStandAlone = fCutRequireITSStandAlone;
 488   target.fCutRequireITSpureSA = fCutRequireITSpureSA;
 489
 490   target.fCutNsigmaToVertex = fCutNsigmaToVertex;
 491   target.fCutSigmaToVertexRequired = fCutSigmaToVertexRequired;
 492   target.fCutMaxDCAToVertexXY = fCutMaxDCAToVertexXY;
 493   target.fCutMaxDCAToVertexZ = fCutMaxDCAToVertexZ;
 494   target.fCutDCAToVertex2D = fCutDCAToVertex2D;
 495   target.fCutMinDCAToVertexXY = fCutMinDCAToVertexXY;
 496   target.fCutMinDCAToVertexZ = fCutMinDCAToVertexZ;
 497
 498   target.fCutMaxDCAToVertexXYPtDep = fCutMaxDCAToVertexXYPtDep;
 499   target.SetMaxDCAToVertexXYPtDep(fCutMaxDCAToVertexXYPtDep.Data());
 500
 501   target.fCutMaxDCAToVertexZPtDep = fCutMaxDCAToVertexZPtDep;
 502   target.SetMaxDCAToVertexZPtDep(fCutMaxDCAToVertexZPtDep.Data());
 503
 504   target.fCutMinDCAToVertexXYPtDep = fCutMinDCAToVertexXYPtDep;
 505   target.SetMinDCAToVertexXYPtDep(fCutMinDCAToVertexXYPtDep.Data());
 506
 507   target.fCutMinDCAToVertexZPtDep = fCutMinDCAToVertexZPtDep;
 508   target.SetMinDCAToVertexZPtDep(fCutMinDCAToVertexZPtDep.Data());
 509
 510   target.fPMin = fPMin;
 511   target.fPMax = fPMax;
 512   target.fPtMin = fPtMin;
 513   target.fPtMax = fPtMax;
 514   target.fPxMin = fPxMin;
 515   target.fPxMax = fPxMax;
 516   target.fPyMin = fPyMin;
 517   target.fPyMax = fPyMax;
 518   target.fPzMin = fPzMin;
 519   target.fPzMax = fPzMax;
 520   target.fEtaMin = fEtaMin;
 521   target.fEtaMax = fEtaMax;
 522   target.fRapMin = fRapMin;
 523   target.fRapMax = fRapMax;
 524
 525   target.fHistogramsOn = fHistogramsOn;
 526
 527   for (Int_t i=0; i<2; ++i)
 528   {
 529     if (fhNClustersITS[i]) target.fhNClustersITS[i] = (TH1F*) fhNClustersITS[i]->Clone();
 530     if (fhNClustersTPC[i]) target.fhNClustersTPC[i] = (TH1F*) fhNClustersTPC[i]->Clone();
 531     if (fhNSharedClustersTPC[i]) target.fhNSharedClustersTPC[i] = (TH1F*) fhNSharedClustersTPC[i]->Clone();
 532     if (fhNCrossedRowsTPC[i]) target.fhNCrossedRowsTPC[i] = (TH1F*) fhNCrossedRowsTPC[i]->Clone();
 533     if (fhRatioCrossedRowsOverFindableClustersTPC[i]) target.fhRatioCrossedRowsOverFindableClustersTPC[i] = (TH1F*) fhRatioCrossedRowsOverFindableClustersTPC[i]->Clone();
 534
 535     if (fhChi2PerClusterITS[i]) target.fhChi2PerClusterITS[i] = (TH1F*) fhChi2PerClusterITS[i]->Clone();
 536     if (fhChi2PerClusterTPC[i]) target.fhChi2PerClusterTPC[i] = (TH1F*) fhChi2PerClusterTPC[i]->Clone();
 537
 538     if (fhC11[i]) target.fhC11[i] = (TH1F*) fhC11[i]->Clone();
 539     if (fhC22[i]) target.fhC22[i] = (TH1F*) fhC22[i]->Clone();
 540     if (fhC33[i]) target.fhC33[i] = (TH1F*) fhC33[i]->Clone();
 541     if (fhC44[i]) target.fhC44[i] = (TH1F*) fhC44[i]->Clone();
 542     if (fhC55[i]) target.fhC55[i] = (TH1F*) fhC55[i]->Clone();
 543
 544     if (fhRel1PtUncertainty[i]) target.fhRel1PtUncertainty[i] = (TH1F*) fhRel1PtUncertainty[i]->Clone();
 545
 546     if (fhDXY[i]) target.fhDXY[i] = (TH1F*) fhDXY[i]->Clone();
 547     if (fhDZ[i]) target.fhDZ[i] = (TH1F*) fhDZ[i]->Clone();
 548     if (fhDXYDZ[i]) target.fhDXYDZ[i] = (TH1F*) fhDXYDZ[i]->Clone();
 549     if (fhDXYvsDZ[i]) target.fhDXYvsDZ[i] = (TH2F*) fhDXYvsDZ[i]->Clone();
 550
 551     if (fhDXYNormalized[i]) target.fhDXYNormalized[i] = (TH1F*) fhDXYNormalized[i]->Clone();
 552     if (fhDZNormalized[i]) target.fhDZNormalized[i] = (TH1F*) fhDZNormalized[i]->Clone();
 553     if (fhDXYvsDZNormalized[i]) target.fhDXYvsDZNormalized[i] = (TH2F*) fhDXYvsDZNormalized[i]->Clone();
 554     if (fhNSigmaToVertex[i]) target.fhNSigmaToVertex[i] = (TH1F*) fhNSigmaToVertex[i]->Clone();
 555
 556     if (fhPt[i]) target.fhPt[i] = (TH1F*) fhPt[i]->Clone();
 557     if (fhEta[i]) target.fhEta[i] = (TH1F*) fhEta[i]->Clone();
 558   }
 559   if (ffDTheoretical) target.ffDTheoretical = (TF1*) ffDTheoretical->Clone();
 560
 561   if (fhCutStatistics) target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
 562   if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();
 563
 564   TNamed::Copy(c);
 565 }
 566
 567 //_____________________________________________________________________________
 568 Long64_t AliESDtrackCuts::Merge(TCollection* list) {
 569   // Merge a list of AliESDtrackCuts objects with this (needed for PROOF)
 570   // Returns the number of merged objects (including this)
 571   if (!list)
 572     return 0;
 573   if (list->IsEmpty())
 574     return 1;
 575   if (!fHistogramsOn)
 576     return 0;
 577   TIterator* iter = list->MakeIterator();
 578   TObject* obj;
 579
 580   // collection of measured and generated histograms
 581   Int_t count = 0;
 582   while ((obj = iter->Next())) {
 583
 584     AliESDtrackCuts* entry = dynamic_cast<AliESDtrackCuts*>(obj);
 585     if (entry == 0)
 586       continue;
 587
 588     if (!entry->fHistogramsOn)
 589       continue;
 590
 591     for (Int_t i=0; i<2; i++) {
 592
 593       fhNClustersITS[i]      ->Add(entry->fhNClustersITS[i]     );
 594       fhNClustersTPC[i]      ->Add(entry->fhNClustersTPC[i]     );
 595       if (fhNSharedClustersTPC[i])
 596         fhNSharedClustersTPC[i]      ->Add(entry->fhNSharedClustersTPC[i]     );
 597       if (fhNCrossedRowsTPC[i])
 598         fhNCrossedRowsTPC[i]   ->Add(entry->fhNCrossedRowsTPC[i]     );
 599       if (fhRatioCrossedRowsOverFindableClustersTPC[i])
 600         fhRatioCrossedRowsOverFindableClustersTPC[i]      ->Add(entry->fhRatioCrossedRowsOverFindableClustersTPC[i]     );
 601
 602       fhChi2PerClusterITS[i] ->Add(entry->fhChi2PerClusterITS[i]);
 603       fhChi2PerClusterTPC[i] ->Add(entry->fhChi2PerClusterTPC[i]);
 604
 605       fhC11[i]               ->Add(entry->fhC11[i]              );
 606       fhC22[i]               ->Add(entry->fhC22[i]              );
 607       fhC33[i]               ->Add(entry->fhC33[i]              );
 608       fhC44[i]               ->Add(entry->fhC44[i]              );
 609       fhC55[i]               ->Add(entry->fhC55[i]              );
 610
 611       fhRel1PtUncertainty[i] ->Add(entry->fhRel1PtUncertainty[i]);
 612
 613       fhDXY[i]               ->Add(entry->fhDXY[i]              );
 614       fhDZ[i]                ->Add(entry->fhDZ[i]               );
 615       fhDXYDZ[i]             ->Add(entry->fhDXYDZ[i]          );
 616       fhDXYvsDZ[i]           ->Add(entry->fhDXYvsDZ[i]          );
 617
 618       fhDXYNormalized[i]     ->Add(entry->fhDXYNormalized[i]    );
 619       fhDZNormalized[i]      ->Add(entry->fhDZNormalized[i]     );
 620       fhDXYvsDZNormalized[i] ->Add(entry->fhDXYvsDZNormalized[i]);
 621       fhNSigmaToVertex[i]    ->Add(entry->fhNSigmaToVertex[i]);
 622
 623       fhPt[i]                ->Add(entry->fhPt[i]);
 624       fhEta[i]               ->Add(entry->fhEta[i]);
 625     }
 626
 627     fhCutStatistics  ->Add(entry->fhCutStatistics);
 628     fhCutCorrelation ->Add(entry->fhCutCorrelation);
 629
 630     count++;
 631   }
 632   return count+1;
 633 }
 634
 635 //____________________________________________________________________
 636 AliESDtrackCuts* AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
 637 {
 638   // creates an AliESDtrackCuts object and fills it with standard (pre data-taking) values for TPC-only cuts
 639
 640   Printf("AliESDtrackCuts::GetStandardTPCOnlyTrackCuts: Creating track cuts for TPC-only.");
 641
 642   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 643
 644   esdTrackCuts->SetMinNClustersTPC(50);
 645   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
 646   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 647
 648   esdTrackCuts->SetMaxDCAToVertexZ(3.2);
 649   esdTrackCuts->SetMaxDCAToVertexXY(2.4);
 650   esdTrackCuts->SetDCAToVertex2D(kTRUE);
 651
 652   return esdTrackCuts;
 653 }
 654
 655 //____________________________________________________________________
 656 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(Bool_t selPrimaries)
 657 {
 658   // creates an AliESDtrackCuts object and fills it with standard values for ITS-TPC cuts for pp 2009 data
 659
 660   Printf("AliESDtrackCuts::GetStandardITSTPCTrackCuts: Creating track cuts for ITS+TPC.");
 661
 662   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 663
 664   // TPC
 665   esdTrackCuts->SetRequireTPCStandAlone(kTRUE); // to get chi2 and ncls of kTPCin
 666   esdTrackCuts->SetMinNClustersTPC(70);
 667   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
 668   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 669   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 670   // ITS
 671   esdTrackCuts->SetRequireITSRefit(kTRUE);
 672   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 673                                          AliESDtrackCuts::kAny);
 674   if(selPrimaries) {
 675     // 7*(0.0050+0.0060/pt^0.9)
 676     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0350+0.0420/pt^0.9");
 677   }
 678   esdTrackCuts->SetMaxDCAToVertexZ(1.e6);
 679   esdTrackCuts->SetDCAToVertex2D(kFALSE);
 680   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 681   //esdTrackCuts->SetEtaRange(-0.8,+0.8);
 682
 683   return esdTrackCuts;
 684 }
 685
 686 //____________________________________________________________________
 687 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(Bool_t selPrimaries,Int_t clusterCut)
 688 {
 689   // creates an AliESDtrackCuts object and fills it with standard values for ITS-TPC cuts for pp 2010 data
 690   // if clusterCut = 1, the cut on the number of clusters is replaced by
 691   // a cut on the number of crossed rows and on the ration crossed
 692   // rows/findable clusters
 693
 694   Printf("AliESDtrackCuts::GetStandardITSTPCTrackCuts: Creating track cuts for ITS+TPC.");
 695
 696   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 697
 698   // TPC
 699   if(clusterCut == 0)  esdTrackCuts->SetMinNClustersTPC(70);
 700   else if (clusterCut == 1) {
 701     esdTrackCuts->SetMinNCrossedRowsTPC(70);
 702     esdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
 703   }
 704   else {
 705     Printf("Wrong value of the clusterCut parameter (%d), using cut on Nclusters",clusterCut);
 706     esdTrackCuts->SetMinNClustersTPC(70);
 707   }
 708   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
 709   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 710   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 711   // ITS
 712   esdTrackCuts->SetRequireITSRefit(kTRUE);
 713   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 714                                          AliESDtrackCuts::kAny);
 715   if(selPrimaries) {
 716     // 7*(0.0026+0.0050/pt^1.01)
 717     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
 718   }
 719   esdTrackCuts->SetMaxDCAToVertexZ(2);
 720   esdTrackCuts->SetDCAToVertex2D(kFALSE);
 721   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 722
 723   return esdTrackCuts;
 724 }
 725
 726 //____________________________________________________________________
 727
 728
 729
 730 //____________________________________________________________________
 731 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSPureSATrackCuts2009(Bool_t selPrimaries, Bool_t useForPid)
 732 {
 733   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks
 734
 735   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 736   esdTrackCuts->SetRequireITSPureStandAlone(kTRUE);
 737   esdTrackCuts->SetRequireITSRefit(kTRUE);
 738   esdTrackCuts->SetMinNClustersITS(4);
 739   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 740                                          AliESDtrackCuts::kAny);
 741   esdTrackCuts->SetMaxChi2PerClusterITS(1.);
 742
 743   if(selPrimaries) {
 744     // 7*(0.0085+0.0026/pt^1.55)
 745     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0595+0.0182/pt^1.55");
 746   }
 747   if(useForPid){
 748     esdTrackCuts->SetRequireITSPid(kTRUE);
 749   }
 750   return esdTrackCuts;
 751 }
 752
 753 //____________________________________________________________________
 754 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSPureSATrackCuts2010(Bool_t selPrimaries, Bool_t useForPid)
 755 {
 756   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks - pp 2010
 757
 758   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 759   esdTrackCuts->SetRequireITSPureStandAlone(kTRUE);
 760   esdTrackCuts->SetRequireITSRefit(kTRUE);
 761   esdTrackCuts->SetMinNClustersITS(4);
 762   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 763                                          AliESDtrackCuts::kAny);
 764   esdTrackCuts->SetMaxChi2PerClusterITS(2.5);
 765
 766   if(selPrimaries) {
 767     // 7*(0.0033+0.0045/pt^1.3)
 768     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0231+0.0315/pt^1.3");
 769   }
 770   if(useForPid){
 771     esdTrackCuts->SetRequireITSPid(kTRUE);
 772   }
 773   return esdTrackCuts;
 774 }
 775
 776 //____________________________________________________________________
 777 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSSATrackCuts2009(Bool_t selPrimaries, Bool_t useForPid)
 778 {
 779   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks
 780
 781   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 782   esdTrackCuts->SetRequireITSStandAlone(kTRUE);
 783   esdTrackCuts->SetRequireITSPureStandAlone(kFALSE);
 784   esdTrackCuts->SetRequireITSRefit(kTRUE);
 785   esdTrackCuts->SetMinNClustersITS(4);
 786   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 787                                          AliESDtrackCuts::kAny);
 788   esdTrackCuts->SetMaxChi2PerClusterITS(1.);
 789
 790   if(selPrimaries) {
 791     // 7*(0.0085+0.0026/pt^1.55)
 792     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0595+0.0182/pt^1.55");
 793   }
 794   if(useForPid){
 795     esdTrackCuts->SetRequireITSPid(kTRUE);
 796   }
 797   return esdTrackCuts;
 798 }
 799
 800 //____________________________________________________________________
 801 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSSATrackCuts2010(Bool_t selPrimaries, Bool_t useForPid)
 802 {
 803   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks --pp 2010
 804
 805   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 806   esdTrackCuts->SetRequireITSStandAlone(kTRUE);
 807   esdTrackCuts->SetRequireITSPureStandAlone(kFALSE);
 808   esdTrackCuts->SetRequireITSRefit(kTRUE);
 809   esdTrackCuts->SetMinNClustersITS(4);
 810   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 811                                          AliESDtrackCuts::kAny);
 812   esdTrackCuts->SetMaxChi2PerClusterITS(2.5);
 813
 814   if(selPrimaries) {
 815     // 7*(0.0033+0.0045/pt^1.3)
 816     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0231+0.0315/pt^1.3");
 817   }
 818   if(useForPid){
 819     esdTrackCuts->SetRequireITSPid(kTRUE);
 820   }
 821   return esdTrackCuts;
 822 }
 823
 824 //____________________________________________________________________
 825 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSSATrackCutsPbPb2010(Bool_t selPrimaries, Bool_t useForPid)
 826 {
 827   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks -- PbPb 2010
 828
 829   AliESDtrackCuts* esdTrackCuts = GetStandardITSSATrackCuts2010(selPrimaries, useForPid);
 830   esdTrackCuts->SetMaxNOfMissingITSPoints(1);
 831
 832   return esdTrackCuts;
 833 }
 834
 835 //____________________________________________________________________
 836 Int_t AliESDtrackCuts::GetReferenceMultiplicity(AliESDEvent* esd, Bool_t tpcOnly)
 837 {
 838   // Gets reference multiplicity following the standard cuts and a defined fiducial volume
 839   // tpcOnly = kTRUE -> consider TPC-only tracks
 840   //         = kFALSE -> consider global tracks
 841
 842   if (!tpcOnly)
 843   {
 844     Printf("AliESDtrackCuts::GetReferenceMultiplicity: Not implemented for global tracks!");
 845     return -1;
 846   }
 847
 848   static AliESDtrackCuts* esdTrackCuts = 0;
 849   if (!esdTrackCuts)
 850   {
 851     esdTrackCuts = GetStandardTPCOnlyTrackCuts();
 852     esdTrackCuts->SetEtaRange(-0.8, 0.8);
 853     esdTrackCuts->SetPtRange(0.15);
 854   }
 855
 856   Int_t nTracks = esdTrackCuts->CountAcceptedTracks(esd);
 857
 858   return nTracks;
 859 }
 860
 861 //____________________________________________________________________
 862 Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* const esdTrack)
 863 {
 864   // Calculates the number of sigma to the vertex.
 865
 866   Float_t b[2];
 867   Float_t bRes[2];
 868   Float_t bCov[3];
 869   esdTrack->GetImpactParameters(b,bCov);
 870
 871   if (bCov[0]<=0 || bCov[2]<=0) {
 872     AliDebugClass(1, "Estimated b resolution lower or equal zero!");
 873     bCov[0]=0; bCov[2]=0;
 874   }
 875   bRes[0] = TMath::Sqrt(bCov[0]);
 876   bRes[1] = TMath::Sqrt(bCov[2]);
 877
 878   // -----------------------------------
 879   // How to get to a n-sigma cut?
 880   //
 881   // The accumulated statistics from 0 to d is
 882   //
 883   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 884   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 885   //
 886   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
 887   // Can this be expressed in a different way?
 888
 889   if (bRes[0] == 0 || bRes[1] ==0)
 890     return -1;
 891
 892   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 893
 894   // work around precision problem
 895   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 896   // 1e-15 corresponds to nsigma ~ 7.7
 897   if (TMath::Exp(-d * d / 2) < 1e-15)
 898     return 1000;
 899
 900   Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 901   return nSigma;
 902 }
 903
 904 void AliESDtrackCuts::EnableNeededBranches(TTree* tree)
 905 {
 906   // enables the branches needed by AcceptTrack, for a list see comment of AcceptTrack
 907
 908   tree->SetBranchStatus("fTracks.fFlags", 1);
 909   tree->SetBranchStatus("fTracks.fITSncls", 1);
 910   tree->SetBranchStatus("fTracks.fTPCncls", 1);
 911   tree->SetBranchStatus("fTracks.fITSchi2", 1);
 912   tree->SetBranchStatus("fTracks.fTPCchi2", 1);
 913   tree->SetBranchStatus("fTracks.fC*", 1);
 914   tree->SetBranchStatus("fTracks.fD", 1);
 915   tree->SetBranchStatus("fTracks.fZ", 1);
 916   tree->SetBranchStatus("fTracks.fCdd", 1);
 917   tree->SetBranchStatus("fTracks.fCdz", 1);
 918   tree->SetBranchStatus("fTracks.fCzz", 1);
 919   tree->SetBranchStatus("fTracks.fP*", 1);
 920   tree->SetBranchStatus("fTracks.fR*", 1);
 921   tree->SetBranchStatus("fTracks.fKinkIndexes*", 1);
 922 }
 923
 924 //____________________________________________________________________
 925 Bool_t AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack)
 926 {
 927   //
 928   // figure out if the tracks survives all the track cuts defined
 929   //
 930   // the different quality parameter and kinematic values are first
 931   // retrieved from the track. then it is found out what cuts the
 932   // track did not survive and finally the cuts are imposed.
 933
 934   // this function needs the following branches:
 935   // fTracks.fFlags
 936   // fTracks.fITSncls
 937   // fTracks.fTPCncls
 938   // fTracks.fITSchi2
 939   // fTracks.fTPCchi2
 940   // fTracks.fC   //GetExternalCovariance
 941   // fTracks.fD   //GetImpactParameters
 942   // fTracks.fZ   //GetImpactParameters
 943   // fTracks.fCdd //GetImpactParameters
 944   // fTracks.fCdz //GetImpactParameters
 945   // fTracks.fCzz //GetImpactParameters
 946   // fTracks.fP   //GetPxPyPz
 947   // fTracks.fR   //GetMass
 948   // fTracks.fP   //GetMass
 949   // fTracks.fKinkIndexes
 950
 951   UInt_t status = esdTrack->GetStatus();
 952
 953   // getting quality parameters from the ESD track
 954   Int_t nClustersITS = esdTrack->GetITSclusters(0);
 955   Int_t nClustersTPC = -1;
 956   if(fCutRequireTPCStandAlone) {
 957     nClustersTPC = esdTrack->GetTPCNclsIter1();
 958   }
 959   else {
 960     nClustersTPC = esdTrack->GetTPCclusters(0);
 961   }
 962   Float_t nCrossedRowsTPC = esdTrack->GetTPCClusterInfo(2,1);
 963   Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 964   if (esdTrack->GetTPCNclsF()>0) {
 965     ratioCrossedRowsOverFindableClustersTPC = esdTrack->GetTPCClusterInfo(2,1)/esdTrack->GetTPCNclsF();
 966   }
 967
 968   Int_t nClustersTPCShared = esdTrack->GetTPCnclsS();
 969   Float_t fracClustersTPCShared = -1.;
 970
 971   Float_t chi2PerClusterITS = -1;
 972   Float_t chi2PerClusterTPC = -1;
 973   if (nClustersITS!=0)
 974     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
 975   if (nClustersTPC!=0) {
 976     if(fCutRequireTPCStandAlone) {
 977       chi2PerClusterTPC = esdTrack->GetTPCchi2Iter1()/Float_t(nClustersTPC);
 978     } else {
 979       chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
 980     }
 981     fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
 982   }
 983
 984   Double_t extCov[15];
 985   esdTrack->GetExternalCovariance(extCov);
 986
 987   // getting the track to vertex parameters
 988   Float_t nSigmaToVertex = GetSigmaToVertex(esdTrack);
 989
 990   Float_t b[2];
 991   Float_t bCov[3];
 992   esdTrack->GetImpactParameters(b,bCov);
 993   if (bCov[0]<=0 || bCov[2]<=0) {
 994     AliDebug(1, "Estimated b resolution lower or equal zero!");
 995     bCov[0]=0; bCov[2]=0;
 996   }
 997
 998
 999   // set pt-dependent DCA cuts, if requested
1000   SetPtDepDCACuts(esdTrack->Pt());
1001
1002
1003   Float_t dcaToVertexXY = b[0];
1004   Float_t dcaToVertexZ = b[1];
1005
1006   Float_t dcaToVertex = -1;
1007
1008   if (fCutDCAToVertex2D)
1009   {
1010     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1011   }
1012   else
1013     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1014
1015   // getting the kinematic variables of the track
1016   // (assuming the mass is known)
1017   Double_t p[3];
1018   esdTrack->GetPxPyPz(p);
1019
1020   Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
1021   Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
1022   Float_t energy   = TMath::Sqrt(TMath::Power(esdTrack->GetMass(),2) + TMath::Power(momentum,2));
1023
1024   //y-eta related calculations
1025   Float_t eta = -100.;
1026   Float_t y   = -100.;
1027   if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
1028     eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
1029   if((energy != TMath::Abs(p[2]))&&(momentum != 0))
1030     y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
1031
1032   if (extCov[14] < 0)
1033   {
1034     Printf("AliESDtrackCuts::AcceptTrack: WARNING: GetSigma1Pt2() returns negative value for external covariance matrix element fC[14]: %f. Corrupted track information, track will not be accepted!", extCov[14]);
1035     return kFALSE;
1036   }
1037   Float_t relUncertainty1Pt = TMath::Sqrt(extCov[14])*pt;
1038
1039   //########################################################################
1040   // cut the track?
1041
1042   Bool_t cuts[kNCuts];
1043   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1044
1045   // track quality cuts
1046   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1047     cuts[0]=kTRUE;
1048   if (fCutRequireTPCStandAlone && (status&AliESDtrack::kTPCin)==0)
1049     cuts[1]=kTRUE;
1050   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1051     cuts[2]=kTRUE;
1052   if (nClustersTPC<fCutMinNClusterTPC)
1053     cuts[3]=kTRUE;
1054   if (nClustersITS<fCutMinNClusterITS)
1055     cuts[4]=kTRUE;
1056   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1057     cuts[5]=kTRUE;
1058   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1059     cuts[6]=kTRUE;
1060   if (extCov[0]  > fCutMaxC11)
1061     cuts[7]=kTRUE;
1062   if (extCov[2]  > fCutMaxC22)
1063     cuts[8]=kTRUE;
1064   if (extCov[5]  > fCutMaxC33)
1065     cuts[9]=kTRUE;
1066   if (extCov[9]  > fCutMaxC44)
1067     cuts[10]=kTRUE;
1068   if (extCov[14]  > fCutMaxC55)
1069     cuts[11]=kTRUE;
1070   if (nSigmaToVertex > fCutNsigmaToVertex && fCutSigmaToVertexRequired)
1071     cuts[12] = kTRUE;
1072   // if n sigma could not be calculated
1073   if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)
1074     cuts[13]=kTRUE;
1075   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1076     cuts[14]=kTRUE;
1077   // track kinematics cut
1078   if((momentum < fPMin) || (momentum > fPMax))
1079     cuts[15]=kTRUE;
1080   if((pt < fPtMin) || (pt > fPtMax))
1081     cuts[16] = kTRUE;
1082   if((p[0] < fPxMin) || (p[0] > fPxMax))
1083     cuts[17] = kTRUE;
1084   if((p[1] < fPyMin) || (p[1] > fPyMax))
1085     cuts[18] = kTRUE;
1086   if((p[2] < fPzMin) || (p[2] > fPzMax))
1087     cuts[19] = kTRUE;
1088   if((eta < fEtaMin) || (eta > fEtaMax))
1089     cuts[20] = kTRUE;
1090   if((y < fRapMin) || (y > fRapMax))
1091     cuts[21] = kTRUE;
1092   if (fCutDCAToVertex2D && dcaToVertex > 1)
1093     cuts[22] = kTRUE;
1094   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1095     cuts[23] = kTRUE;
1096   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1097     cuts[24] = kTRUE;
1098   if (fCutDCAToVertex2D && fCutMinDCAToVertexXY > 0 && fCutMinDCAToVertexZ > 0 && dcaToVertexXY*dcaToVertexXY/fCutMinDCAToVertexXY/fCutMinDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMinDCAToVertexZ/fCutMinDCAToVertexZ < 1)
1099     cuts[25] = kTRUE;
1100   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) < fCutMinDCAToVertexXY)
1101     cuts[26] = kTRUE;
1102   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) < fCutMinDCAToVertexZ)
1103     cuts[27] = kTRUE;
1104
1105   for (Int_t i = 0; i < 3; i++)
1106     cuts[28+i] = !CheckITSClusterRequirement(fCutClusterRequirementITS[i], esdTrack->HasPointOnITSLayer(i*2), esdTrack->HasPointOnITSLayer(i*2+1));
1107
1108   if(fCutRequireITSStandAlone || fCutRequireITSpureSA){
1109     if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)){
1110       // TPC tracks
1111       cuts[31] = kTRUE;
1112     }else{
1113       // ITS standalone tracks
1114       if(fCutRequireITSStandAlone && !fCutRequireITSpureSA){
1115         if(status & AliESDtrack::kITSpureSA) cuts[31] = kTRUE;
1116       }else if(fCutRequireITSpureSA){
1117         if(!(status & AliESDtrack::kITSpureSA)) cuts[31] = kTRUE;
1118       }
1119     }
1120   }
1121
1122   if (relUncertainty1Pt > fCutMaxRel1PtUncertainty)
1123      cuts[32] = kTRUE;
1124
1125   if (!fCutAcceptSharedTPCClusters && nClustersTPCShared!=0)
1126     cuts[33] = kTRUE;
1127
1128   if (fracClustersTPCShared > fCutMaxFractionSharedTPCClusters)
1129     cuts[34] = kTRUE;
1130
1131   if(fCutRequireITSPid){
1132     UChar_t clumap=esdTrack->GetITSClusterMap();
1133     Int_t nPointsForPid=0;
1134     for(Int_t i=2; i<6; i++){
1135       if(clumap&(1<<i)) ++nPointsForPid;
1136     }
1137     if(nPointsForPid<3) cuts[35] = kTRUE;
1138   }
1139
1140   if (nCrossedRowsTPC<fCutMinNCrossedRowsTPC)
1141     cuts[36]=kTRUE;
1142   if (ratioCrossedRowsOverFindableClustersTPC<fCutMinRatioCrossedRowsOverFindableClustersTPC)
1143     cuts[37]=kTRUE;
1144
1145   Int_t nMissITSpts=0;
1146   Int_t idet,statusLay;
1147   Float_t xloc,zloc;
1148   for(Int_t iLay=0; iLay<6; iLay++){
1149     Bool_t retc=esdTrack->GetITSModuleIndexInfo(iLay,idet,statusLay,xloc,zloc);
1150     if(retc && statusLay==5) ++nMissITSpts;
1151   }
1152   if(nMissITSpts>fCutMaxMissingITSPoints) cuts[38] = kTRUE;
1153
1154   Bool_t cut=kFALSE;
1155   for (Int_t i=0; i<kNCuts; i++)
1156     if (cuts[i]) {cut = kTRUE;}
1157
1158
1159   //########################################################################
1160   // filling histograms
1161   if (fHistogramsOn) {
1162     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
1163     if (cut)
1164       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
1165
1166     for (Int_t i=0; i<kNCuts; i++) {
1167       if (fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i]) < 1)
1168         AliFatal(Form("Inconsistency! Cut %d with name %s not found", i, fgkCutNames[i]));
1169
1170       if (cuts[i])
1171         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
1172
1173       for (Int_t j=i; j<kNCuts; j++) {
1174         if (cuts[i] && cuts[j]) {
1175           Float_t xC = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
1176           Float_t yC = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
1177           fhCutCorrelation->Fill(xC, yC);
1178         }
1179       }
1180     }
1181   }
1182
1183   // now we loop over the filling of the histograms twice: once "before" the cut, once "after"
1184   // the code is not in a function due to too many local variables that would need to be passed
1185
1186   for (Int_t id = 0; id < 2; id++)
1187   {
1188     // id = 0 --> before cut
1189     // id = 1 --> after cut
1190
1191     if (fHistogramsOn)
1192     {
1193       fhNClustersITS[id]->Fill(nClustersITS);
1194       fhNClustersTPC[id]->Fill(nClustersTPC);
1195       fhNSharedClustersTPC[id]->Fill(nClustersTPCShared);
1196       fhNCrossedRowsTPC[id]->Fill(nCrossedRowsTPC);
1197       fhRatioCrossedRowsOverFindableClustersTPC[id]->Fill(ratioCrossedRowsOverFindableClustersTPC);
1198       fhChi2PerClusterITS[id]->Fill(chi2PerClusterITS);
1199       fhChi2PerClusterTPC[id]->Fill(chi2PerClusterTPC);
1200
1201       fhC11[id]->Fill(extCov[0]);
1202       fhC22[id]->Fill(extCov[2]);
1203       fhC33[id]->Fill(extCov[5]);
1204       fhC44[id]->Fill(extCov[9]);
1205       fhC55[id]->Fill(extCov[14]);
1206
1207       fhRel1PtUncertainty[id]->Fill(relUncertainty1Pt);
1208
1209       fhPt[id]->Fill(pt);
1210       fhEta[id]->Fill(eta);
1211
1212       Float_t bRes[2];
1213       bRes[0] = TMath::Sqrt(bCov[0]);
1214       bRes[1] = TMath::Sqrt(bCov[2]);
1215
1216       fhDZ[id]->Fill(b[1]);
1217       fhDXY[id]->Fill(b[0]);
1218       fhDXYDZ[id]->Fill(dcaToVertex);
1219       fhDXYvsDZ[id]->Fill(b[1],b[0]);
1220
1221       if (bRes[0]!=0 && bRes[1]!=0) {
1222         fhDZNormalized[id]->Fill(b[1]/bRes[1]);
1223         fhDXYNormalized[id]->Fill(b[0]/bRes[0]);
1224         fhDXYvsDZNormalized[id]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
1225         fhNSigmaToVertex[id]->Fill(nSigmaToVertex);
1226       }
1227     }
1228
1229     // cut the track
1230     if (cut)
1231       return kFALSE;
1232   }
1233
1234   return kTRUE;
1235 }
1236
1237 //____________________________________________________________________
1238 Bool_t AliESDtrackCuts::CheckITSClusterRequirement(ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2)
1239 {
1240   // checks if the cluster requirement is fullfilled (in this case: return kTRUE)
1241
1242   switch (req)
1243   {
1244         case kOff:        return kTRUE;
1245         case kNone:       return !clusterL1 && !clusterL2;
1246         case kAny:        return clusterL1 || clusterL2;
1247         case kFirst:      return clusterL1;
1248         case kOnlyFirst:  return clusterL1 && !clusterL2;
1249         case kSecond:     return clusterL2;
1250         case kOnlySecond: return clusterL2 && !clusterL1;
1251         case kBoth:       return clusterL1 && clusterL2;
1252   }
1253
1254   return kFALSE;
1255 }
1256
1257 //____________________________________________________________________
1258 AliESDtrack* AliESDtrackCuts::GetTPCOnlyTrack(AliESDEvent* esd, Int_t iTrack)
1259 {
1260
1261   // Utility function to
1262   // create a TPC only track from the given esd track
1263   //
1264   // IMPORTANT: The track has to be deleted by the user
1265   //
1266   // NB. most of the functionality to get a TPC only track from an ESD track is in AliESDtrack, where it should be
1267   // there are only missing propagations here that are needed for old data
1268   // this function will therefore become obsolete
1269   //
1270   // adapted from code provided by CKB
1271
1272   if (!esd->GetPrimaryVertexTPC())
1273     return 0; // No TPC vertex no TPC tracks
1274
1275   if(!esd->GetPrimaryVertexTPC()->GetStatus())
1276     return 0; // TPC Vertex is created by default in AliESDEvent, do not use in this case
1277
1278   AliESDtrack* track = esd->GetTrack(iTrack);
1279   if (!track)
1280     return 0;
1281
1282   AliESDtrack *tpcTrack = new AliESDtrack();
1283
1284   // only true if we have a tpc track
1285   if (!track->FillTPCOnlyTrack(*tpcTrack))
1286   {
1287     delete tpcTrack;
1288     return 0;
1289   }
1290
1291   // propagate to Vertex
1292   // not needed for normal reconstructed ESDs...
1293   // Double_t pTPC[2],covTPC[3];
1294   // tpcTrack->PropagateToDCA(esd->GetPrimaryVertexTPC(), esd->GetMagneticField(), 10000,  pTPC, covTPC);
1295
1296   return tpcTrack;
1297 }
1298
1299 //____________________________________________________________________
1300 TObjArray* AliESDtrackCuts::GetAcceptedTracks(AliESDEvent* esd,Bool_t bTPC)
1301 {
1302   //
1303   // returns an array of all tracks that pass the cuts
1304   // or an array of TPC only tracks (propagated to the TPC vertex during reco)
1305   // tracks that pass the cut
1306   //
1307   // NOTE: List has to be deleted by the user
1308
1309   TObjArray* acceptedTracks = new TObjArray();
1310
1311   // loop over esd tracks
1312   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
1313     if(bTPC){
1314       if(!esd->GetPrimaryVertexTPC())return acceptedTracks; // No TPC vertex no TPC tracks
1315       if(!esd->GetPrimaryVertexTPC()->GetStatus())return acceptedTracks; // No proper TPC vertex, only the default
1316
1317       AliESDtrack *tpcTrack = GetTPCOnlyTrack(esd, iTrack);
1318       if (!tpcTrack)
1319         continue;
1320
1321       if (AcceptTrack(tpcTrack)) {
1322         acceptedTracks->Add(tpcTrack);
1323       }
1324       else
1325         delete tpcTrack;
1326     }
1327     else
1328     {
1329       AliESDtrack* track = esd->GetTrack(iTrack);
1330       if(AcceptTrack(track))
1331         acceptedTracks->Add(track);
1332     }
1333   }
1334   if(bTPC)acceptedTracks->SetOwner(kTRUE);
1335   return acceptedTracks;
1336 }
1337
1338 //____________________________________________________________________
1339 Int_t AliESDtrackCuts::CountAcceptedTracks(AliESDEvent* const esd)
1340 {
1341   //
1342   // returns an the number of tracks that pass the cuts
1343   //
1344
1345   Int_t count = 0;
1346
1347   // loop over esd tracks
1348   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
1349     AliESDtrack* track = esd->GetTrack(iTrack);
1350     if (AcceptTrack(track))
1351       count++;
1352   }
1353
1354   return count;
1355 }
1356
1357 //____________________________________________________________________
1358  void AliESDtrackCuts::DefineHistograms(Int_t color) {
1359    //
1360    // diagnostics histograms are defined
1361    //
1362
1363    fHistogramsOn=kTRUE;
1364
1365    Bool_t oldStatus = TH1::AddDirectoryStatus();
1366    TH1::AddDirectory(kFALSE);
1367
1368    //###################################################################################
1369    // defining histograms
1370
1371    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
1372
1373    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
1374    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
1375
1376    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
1377
1378    for (Int_t i=0; i<kNCuts; i++) {
1379      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
1380      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1381      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1382    }
1383
1384   fhCutStatistics  ->SetLineColor(color);
1385   fhCutCorrelation ->SetLineColor(color);
1386   fhCutStatistics  ->SetLineWidth(2);
1387   fhCutCorrelation ->SetLineWidth(2);
1388
1389   for (Int_t i=0; i<2; i++) {
1390     fhNClustersITS[i]        = new TH1F("nClustersITS"    ,"",8,-0.5,7.5);
1391     fhNClustersTPC[i]        = new TH1F("nClustersTPC"     ,"",165,-0.5,164.5);
1392     fhNSharedClustersTPC[i]  = new TH1F("nSharedClustersTPC"     ,"",165,-0.5,164.5);
1393     fhNCrossedRowsTPC[i]     = new TH1F("nCrossedRowsTPC"  ,"",165,-0.5,164.5);
1394     fhRatioCrossedRowsOverFindableClustersTPC[i]     = new TH1F("ratioCrossedRowsOverFindableClustersTPC"  ,"",60,0,1.5);
1395     fhChi2PerClusterITS[i]   = new TH1F("chi2PerClusterITS","",500,0,10);
1396     fhChi2PerClusterTPC[i]   = new TH1F("chi2PerClusterTPC","",500,0,10);
1397
1398     fhC11[i]                 = new TH1F("covMatrixDiagonal11","",2000,0,20);
1399     fhC22[i]                 = new TH1F("covMatrixDiagonal22","",2000,0,20);
1400     fhC33[i]                 = new TH1F("covMatrixDiagonal33","",1000,0,0.1);
1401     fhC44[i]                 = new TH1F("covMatrixDiagonal44","",1000,0,0.1);
1402     fhC55[i]                 = new TH1F("covMatrixDiagonal55","",1000,0,5);
1403
1404     fhRel1PtUncertainty[i]   = new TH1F("rel1PtUncertainty","",1000,0,5);
1405
1406     fhDXY[i]                 = new TH1F("dXY"    ,"",500,-10,10);
1407     fhDZ[i]                  = new TH1F("dZ"     ,"",500,-10,10);
1408     fhDXYDZ[i]               = new TH1F("dXYDZ"  ,"",500,0,10);
1409     fhDXYvsDZ[i]             = new TH2F("dXYvsDZ","",200,-10,10,200,-10,10);
1410
1411     fhDXYNormalized[i]       = new TH1F("dXYNormalized"    ,"",500,-10,10);
1412     fhDZNormalized[i]        = new TH1F("dZNormalized"     ,"",500,-10,10);
1413     fhDXYvsDZNormalized[i]   = new TH2F("dXYvsDZNormalized","",200,-10,10,200,-10,10);
1414
1415     fhNSigmaToVertex[i]      = new TH1F("nSigmaToVertex","",500,0,10);
1416
1417     fhPt[i]                  = new TH1F("pt"     ,"p_{T} distribution;p_{T} (GeV/c)", 800, 0.0, 10.0);
1418     fhEta[i]                 = new TH1F("eta"     ,"#eta distribution;#eta",40,-2.0,2.0);
1419
1420     fhNClustersITS[i]->SetTitle("n ITS clusters");
1421     fhNClustersTPC[i]->SetTitle("n TPC clusters");
1422     fhNSharedClustersTPC[i]->SetTitle("n TPC shared clusters");
1423     fhChi2PerClusterITS[i]->SetTitle("#Chi^{2} per ITS cluster");
1424     fhChi2PerClusterTPC[i]->SetTitle("#Chi^{2} per TPC cluster");
1425
1426     fhC11[i]->SetTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
1427     fhC22[i]->SetTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
1428     fhC33[i]->SetTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
1429     fhC44[i]->SetTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
1430     fhC55[i]->SetTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
1431
1432     fhRel1PtUncertainty[i]->SetTitle("rel. uncertainty of 1/p_{T}");
1433
1434     fhDXY[i]->SetXTitle("transverse impact parameter (cm)");
1435     fhDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1436     fhDXYDZ[i]->SetTitle("absolute impact parameter;sqrt(dXY**2 + dZ**2)  (cm)");
1437     fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1438     fhDXYvsDZ[i]->SetYTitle("transverse impact parameter (cm)");
1439
1440     fhDXYNormalized[i]->SetTitle("normalized trans impact par (n#sigma)");
1441     fhDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1442     fhDXYvsDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1443     fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par (n#sigma)");
1444     fhNSigmaToVertex[i]->SetTitle("n #sigma to vertex");
1445
1446     fhNClustersITS[i]->SetLineColor(color);   fhNClustersITS[i]->SetLineWidth(2);
1447     fhNClustersTPC[i]->SetLineColor(color);   fhNClustersTPC[i]->SetLineWidth(2);
1448     fhNSharedClustersTPC[i]->SetLineColor(color);  fhNSharedClustersTPC[i]->SetLineWidth(2);
1449     fhChi2PerClusterITS[i]->SetLineColor(color);   fhChi2PerClusterITS[i]->SetLineWidth(2);
1450     fhChi2PerClusterTPC[i]->SetLineColor(color);   fhChi2PerClusterTPC[i]->SetLineWidth(2);
1451
1452     fhC11[i]->SetLineColor(color);   fhC11[i]->SetLineWidth(2);
1453     fhC22[i]->SetLineColor(color);   fhC22[i]->SetLineWidth(2);
1454     fhC33[i]->SetLineColor(color);   fhC33[i]->SetLineWidth(2);
1455     fhC44[i]->SetLineColor(color);   fhC44[i]->SetLineWidth(2);
1456     fhC55[i]->SetLineColor(color);   fhC55[i]->SetLineWidth(2);
1457
1458     fhRel1PtUncertainty[i]->SetLineColor(color); fhRel1PtUncertainty[i]->SetLineWidth(2);
1459
1460     fhDXY[i]->SetLineColor(color);   fhDXY[i]->SetLineWidth(2);
1461     fhDZ[i]->SetLineColor(color);    fhDZ[i]->SetLineWidth(2);
1462     fhDXYDZ[i]->SetLineColor(color); fhDXYDZ[i]->SetLineWidth(2);
1463
1464     fhDXYNormalized[i]->SetLineColor(color);   fhDXYNormalized[i]->SetLineWidth(2);
1465     fhDZNormalized[i]->SetLineColor(color);    fhDZNormalized[i]->SetLineWidth(2);
1466     fhNSigmaToVertex[i]->SetLineColor(color);  fhNSigmaToVertex[i]->SetLineWidth(2);
1467   }
1468
1469   // The number of sigmas to the vertex is per definition gaussian
1470   ffDTheoretical = new TF1("nSigmaToVertexTheoretical","([0]/2.506628274)*exp(-(x**2)/2)",0,50);
1471   ffDTheoretical->SetParameter(0,1);
1472
1473   TH1::AddDirectory(oldStatus);
1474 }
1475
1476 //____________________________________________________________________
1477 Bool_t AliESDtrackCuts::LoadHistograms(const Char_t* dir)
1478 {
1479   //
1480   // loads the histograms from a file
1481   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
1482   //
1483
1484   if (!dir)
1485     dir = GetName();
1486
1487   if (!gDirectory->cd(dir))
1488     return kFALSE;
1489
1490   ffDTheoretical = dynamic_cast<TF1*> (gDirectory->Get("nSigmaToVertexTheory"));
1491
1492   fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
1493   fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
1494
1495   for (Int_t i=0; i<2; i++) {
1496     if (i==0)
1497     {
1498       gDirectory->cd("before_cuts");
1499     }
1500     else
1501       gDirectory->cd("after_cuts");
1502
1503     fhNClustersITS[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersITS"     ));
1504     fhNClustersTPC[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersTPC"     ));
1505     fhNSharedClustersTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSharedClustersTPC"     ));
1506     fhNCrossedRowsTPC[i]   = dynamic_cast<TH1F*> (gDirectory->Get("nCrossedRowsTPC"  ));
1507     fhRatioCrossedRowsOverFindableClustersTPC[i]   = dynamic_cast<TH1F*> (gDirectory->Get("ratioCrossedRowsOverFindableClustersTPC"  ));
1508     fhChi2PerClusterITS[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterITS"));
1509     fhChi2PerClusterTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterTPC"));
1510
1511     fhC11[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal11"));
1512     fhC22[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal22"));
1513     fhC33[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal33"));
1514     fhC44[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal44"));
1515     fhC55[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal55"));
1516
1517     fhRel1PtUncertainty[i] = dynamic_cast<TH1F*> (gDirectory->Get("rel1PtUncertainty"));
1518
1519     fhDXY[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXY"    ));
1520     fhDZ[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZ"     ));
1521     fhDXYDZ[i] =   dynamic_cast<TH1F*> (gDirectory->Get("dXYDZ"));
1522     fhDXYvsDZ[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZ"));
1523
1524     fhDXYNormalized[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXYNormalized"    ));
1525     fhDZNormalized[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZNormalized"     ));
1526     fhDXYvsDZNormalized[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZNormalized"));
1527     fhNSigmaToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSigmaToVertex"));
1528
1529     fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get("pt"));
1530     fhEta[i] = dynamic_cast<TH1F*> (gDirectory->Get("eta"));
1531
1532     gDirectory->cd("../");
1533   }
1534
1535   gDirectory->cd("..");
1536
1537   return kTRUE;
1538 }
1539
1540 //____________________________________________________________________
1541 void AliESDtrackCuts::SaveHistograms(const Char_t* dir) {
1542   //
1543   // saves the histograms in a directory (dir)
1544   //
1545
1546   if (!fHistogramsOn) {
1547     AliDebug(0, "Histograms not on - cannot save histograms!!!");
1548     return;
1549   }
1550
1551   if (!dir)
1552     dir = GetName();
1553
1554   gDirectory->mkdir(dir);
1555   gDirectory->cd(dir);
1556
1557   gDirectory->mkdir("before_cuts");
1558   gDirectory->mkdir("after_cuts");
1559
1560   // a factor of 2 is needed since n sigma is positive
1561   ffDTheoretical->SetParameter(0,2*fhNSigmaToVertex[0]->Integral("width"));
1562   ffDTheoretical->Write("nSigmaToVertexTheory");
1563
1564   fhCutStatistics->Write();
1565   fhCutCorrelation->Write();
1566
1567   for (Int_t i=0; i<2; i++) {
1568     if (i==0)
1569       gDirectory->cd("before_cuts");
1570     else
1571       gDirectory->cd("after_cuts");
1572
1573     fhNClustersITS[i]        ->Write();
1574     fhNClustersTPC[i]        ->Write();
1575     fhNSharedClustersTPC[i]  ->Write();
1576     fhNCrossedRowsTPC[i]     ->Write();
1577     fhRatioCrossedRowsOverFindableClustersTPC[i]     ->Write();
1578     fhChi2PerClusterITS[i]   ->Write();
1579     fhChi2PerClusterTPC[i]   ->Write();
1580
1581     fhC11[i]                 ->Write();
1582     fhC22[i]                 ->Write();
1583     fhC33[i]                 ->Write();
1584     fhC44[i]                 ->Write();
1585     fhC55[i]                 ->Write();
1586
1587     fhRel1PtUncertainty[i]   ->Write();
1588
1589     fhDXY[i]                 ->Write();
1590     fhDZ[i]                  ->Write();
1591     fhDXYDZ[i]               ->Write();
1592     fhDXYvsDZ[i]             ->Write();
1593
1594     fhDXYNormalized[i]       ->Write();
1595     fhDZNormalized[i]        ->Write();
1596     fhDXYvsDZNormalized[i]   ->Write();
1597     fhNSigmaToVertex[i]      ->Write();
1598
1599     fhPt[i]                  ->Write();
1600     fhEta[i]                 ->Write();
1601
1602     gDirectory->cd("../");
1603   }
1604
1605   gDirectory->cd("../");
1606 }
1607
1608 //____________________________________________________________________
1609 void AliESDtrackCuts::DrawHistograms()
1610 {
1611   // draws some histograms
1612
1613   TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "Track Quality Results1", 800, 800);
1614   canvas1->Divide(2, 2);
1615
1616   canvas1->cd(1);
1617   fhNClustersTPC[0]->SetStats(kFALSE);
1618   fhNClustersTPC[0]->Draw();
1619
1620   canvas1->cd(2);
1621   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1622   fhChi2PerClusterTPC[0]->Draw();
1623
1624   canvas1->cd(3);
1625   fhNSigmaToVertex[0]->SetStats(kFALSE);
1626   fhNSigmaToVertex[0]->GetXaxis()->SetRangeUser(0, 10);
1627   fhNSigmaToVertex[0]->Draw();
1628
1629   canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
1630
1631   TCanvas* canvas2 = new TCanvas(Form("%s_2", GetName()), "Track Quality Results2", 1200, 800);
1632   canvas2->Divide(3, 2);
1633
1634   canvas2->cd(1);
1635   fhC11[0]->SetStats(kFALSE);
1636   gPad->SetLogy();
1637   fhC11[0]->Draw();
1638
1639   canvas2->cd(2);
1640   fhC22[0]->SetStats(kFALSE);
1641   gPad->SetLogy();
1642   fhC22[0]->Draw();
1643
1644   canvas2->cd(3);
1645   fhC33[0]->SetStats(kFALSE);
1646   gPad->SetLogy();
1647   fhC33[0]->Draw();
1648
1649   canvas2->cd(4);
1650   fhC44[0]->SetStats(kFALSE);
1651   gPad->SetLogy();
1652   fhC44[0]->Draw();
1653
1654   canvas2->cd(5);
1655   fhC55[0]->SetStats(kFALSE);
1656   gPad->SetLogy();
1657   fhC55[0]->Draw();
1658
1659   canvas2->cd(6);
1660   fhRel1PtUncertainty[0]->SetStats(kFALSE);
1661   gPad->SetLogy();
1662   fhRel1PtUncertainty[0]->Draw();
1663
1664   canvas2->SaveAs(Form("%s_%s.gif", GetName(), canvas2->GetName()));
1665
1666   TCanvas* canvas3 = new TCanvas(Form("%s_3", GetName()), "Track Quality Results3", 1200, 800);
1667   canvas3->Divide(3, 2);
1668
1669   canvas3->cd(1);
1670   fhDXY[0]->SetStats(kFALSE);
1671   gPad->SetLogy();
1672   fhDXY[0]->Draw();
1673
1674   canvas3->cd(2);
1675   fhDZ[0]->SetStats(kFALSE);
1676   gPad->SetLogy();
1677   fhDZ[0]->Draw();
1678
1679   canvas3->cd(3);
1680   fhDXYvsDZ[0]->SetStats(kFALSE);
1681   gPad->SetLogz();
1682   gPad->SetRightMargin(0.15);
1683   fhDXYvsDZ[0]->Draw("COLZ");
1684
1685   canvas3->cd(4);
1686   fhDXYNormalized[0]->SetStats(kFALSE);
1687   gPad->SetLogy();
1688   fhDXYNormalized[0]->Draw();
1689
1690   canvas3->cd(5);
1691   fhDZNormalized[0]->SetStats(kFALSE);
1692   gPad->SetLogy();
1693   fhDZNormalized[0]->Draw();
1694
1695   canvas3->cd(6);
1696   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1697   gPad->SetLogz();
1698   gPad->SetRightMargin(0.15);
1699   fhDXYvsDZNormalized[0]->Draw("COLZ");
1700
1701   canvas3->SaveAs(Form("%s_%s.gif", GetName(), canvas3->GetName()));
1702
1703   TCanvas* canvas4 = new TCanvas(Form("%s_4", GetName()), "Track Quality Results4", 800, 500);
1704   canvas4->Divide(2, 1);
1705
1706   canvas4->cd(1);
1707   fhCutStatistics->SetStats(kFALSE);
1708   fhCutStatistics->LabelsOption("v");
1709   gPad->SetBottomMargin(0.3);
1710   fhCutStatistics->Draw();
1711
1712   canvas4->cd(2);
1713   fhCutCorrelation->SetStats(kFALSE);
1714   fhCutCorrelation->LabelsOption("v");
1715   gPad->SetBottomMargin(0.3);
1716   gPad->SetLeftMargin(0.3);
1717   fhCutCorrelation->Draw("COLZ");
1718
1719   canvas4->SaveAs(Form("%s_%s.gif", GetName(), canvas4->GetName()));
1720
1721   /*canvas->cd(1);
1722   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1723   fhDXYvsDZNormalized[0]->DrawCopy("COLZ");
1724
1725   canvas->cd(2);
1726   fhNClustersTPC[0]->SetStats(kFALSE);
1727   fhNClustersTPC[0]->DrawCopy();
1728
1729   canvas->cd(3);
1730   fhChi2PerClusterITS[0]->SetStats(kFALSE);
1731   fhChi2PerClusterITS[0]->DrawCopy();
1732   fhChi2PerClusterITS[1]->SetLineColor(2);
1733   fhChi2PerClusterITS[1]->DrawCopy("SAME");
1734
1735   canvas->cd(4);
1736   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1737   fhChi2PerClusterTPC[0]->DrawCopy();
1738   fhChi2PerClusterTPC[1]->SetLineColor(2);
1739   fhChi2PerClusterTPC[1]->DrawCopy("SAME");*/
1740 }
1741 //--------------------------------------------------------------------------
1742 void AliESDtrackCuts::SetPtDepDCACuts(Double_t pt) {
1743   //
1744   // set the pt-dependent DCA cuts
1745   //
1746
1747   if(f1CutMaxDCAToVertexXYPtDep) {
1748      fCutMaxDCAToVertexXY=f1CutMaxDCAToVertexXYPtDep->Eval(pt);
1749   }
1750
1751   if(f1CutMaxDCAToVertexZPtDep) {
1752     fCutMaxDCAToVertexZ=f1CutMaxDCAToVertexZPtDep->Eval(pt);
1753   }
1754
1755   if(f1CutMinDCAToVertexXYPtDep) {
1756     fCutMinDCAToVertexXY=f1CutMinDCAToVertexXYPtDep->Eval(pt);
1757   }
1758
1759   if(f1CutMinDCAToVertexZPtDep) {
1760     fCutMinDCAToVertexZ=f1CutMinDCAToVertexZPtDep->Eval(pt);
1761   }
1762
1763
1764   return;
1765 }
1766
1767
1768
1769 //--------------------------------------------------------------------------
1770 Bool_t AliESDtrackCuts::CheckPtDepDCA(TString dist,Bool_t print) const {
1771   //
1772   // Check the correctness of the string syntax
1773   //
1774   Bool_t retval=kTRUE;
1775
1776   if(!dist.Contains("pt")) {
1777     if(print) printf("AliESDtrackCuts::CheckPtDepDCA(): string must contain \"pt\"\n");
1778     retval= kFALSE;
1779   }
1780   return retval;
1781 }
1782
1783  void AliESDtrackCuts::SetMaxDCAToVertexXYPtDep(const char *dist){
1784
1785    if(f1CutMaxDCAToVertexXYPtDep){
1786      delete f1CutMaxDCAToVertexXYPtDep;
1787      // resetiing both
1788      f1CutMaxDCAToVertexXYPtDep = 0;
1789      fCutMaxDCAToVertexXYPtDep = "";
1790    }
1791    if(!CheckPtDepDCA(dist,kTRUE)){
1792      return;
1793    }
1794    fCutMaxDCAToVertexXYPtDep = dist;
1795    TString tmp(dist);
1796    tmp.ReplaceAll("pt","x");
1797    f1CutMaxDCAToVertexXYPtDep = new TFormula("f1CutMaxDCAToVertexXYPtDep",tmp.Data());
1798
1799 }
1800
1801  void AliESDtrackCuts::SetMaxDCAToVertexZPtDep(const char *dist){
1802
1803
1804    if(f1CutMaxDCAToVertexZPtDep){
1805      delete f1CutMaxDCAToVertexZPtDep;
1806      // resetiing both
1807      f1CutMaxDCAToVertexZPtDep = 0;
1808      fCutMaxDCAToVertexZPtDep = "";
1809    }
1810    if(!CheckPtDepDCA(dist,kTRUE))return;
1811
1812    fCutMaxDCAToVertexZPtDep = dist;
1813    TString tmp(dist);
1814    tmp.ReplaceAll("pt","x");
1815    f1CutMaxDCAToVertexZPtDep = new TFormula("f1CutMaxDCAToVertexZPtDep",tmp.Data());
1816
1817
1818 }
1819
1820
1821  void AliESDtrackCuts::SetMinDCAToVertexXYPtDep(const char *dist){
1822
1823
1824    if(f1CutMinDCAToVertexXYPtDep){
1825      delete f1CutMinDCAToVertexXYPtDep;
1826      // resetiing both
1827      f1CutMinDCAToVertexXYPtDep = 0;
1828      fCutMinDCAToVertexXYPtDep = "";
1829    }
1830    if(!CheckPtDepDCA(dist,kTRUE))return;
1831
1832    fCutMinDCAToVertexXYPtDep = dist;
1833    TString tmp(dist);
1834    tmp.ReplaceAll("pt","x");
1835    f1CutMinDCAToVertexXYPtDep = new TFormula("f1CutMinDCAToVertexXYPtDep",tmp.Data());
1836
1837 }
1838
1839
1840  void AliESDtrackCuts::SetMinDCAToVertexZPtDep(const char *dist){
1841
1842
1843
1844    if(f1CutMinDCAToVertexZPtDep){
1845      delete f1CutMinDCAToVertexZPtDep;
1846      // resetiing both
1847      f1CutMinDCAToVertexZPtDep = 0;
1848      fCutMinDCAToVertexZPtDep = "";
1849    }
1850    if(!CheckPtDepDCA(dist,kTRUE))return;
1851    fCutMinDCAToVertexZPtDep = dist;
1852    TString tmp(dist);
1853    tmp.ReplaceAll("pt","x");
1854    f1CutMinDCAToVertexZPtDep = new TFormula("f1CutMinDCAToVertexZPtDep",tmp.Data());
1855 }